Bring geometry processing (CDT + arrangement) improvements / fixes back over from dev-physics

#rb none

[CL 6351809 by Jimmy Andrews in Dev-Editor branch]
This commit is contained in:
Jimmy Andrews
2019-05-07 18:34:48 -04:00
parent cfe69331d0
commit 3f2ea1e00e
4 changed files with 255 additions and 27 deletions
@@ -8,59 +8,245 @@
#include "ThirdParty/GTEngine/Mathematics/GteUIntegerFP32.h"
#include <vector>
bool FArrangement2d::AttemptTriangulate(TArray<FIntVector> &Triangles, TArray<int32> &SkippedEdges)
namespace
{
//#define DEBUG_FILE_DUMPING 1
#ifndef DEBUG_FILE_DUMPING
void DumpGraphForDebug(const FDynamicGraph2d& Graph, const FString& PathBase)
{
}
void DumpGraphForDebugAsOBJ(const FDynamicGraph2d& Graph, const FString& PathBase)
{
}
void DumpTriangulationForDebug(const FDynamicGraph2d& Graph, const TArray<FIntVector>& Triangles, const FString& PathBase)
{
}
#else
#include <fstream>
static int num = 0;
void DumpGraphForDebug(const FDynamicGraph2d& Graph, const FString& PathBase)
{
num++;
FString Path = PathBase + FString::FromInt(num) + ".txt";
std::ofstream f(*Path);
for (int32 VertexIdx = 0; VertexIdx < Graph.MaxVertexID(); VertexIdx++)
{
const FVector2d& Vertex = Graph.GetVertex(VertexIdx);
f << "v " << Vertex.X << " " << Vertex.Y << " " << std::endl;
}
for (const FDynamicGraph::FEdge& Edge : Graph.Edges())
{
f << "e " << Edge.A << " " << Edge.B << std::endl;
}
f.close();
}
void DumpGraphForDebugAsOBJ(const FDynamicGraph2d& Graph, const FString& PathBase)
{
num++;
FString Path = PathBase + FString::FromInt(num) + ".obj";
std::ofstream f(*Path);
for (int32 VertexIdx = 0; VertexIdx < Graph.MaxVertexID(); VertexIdx++)
{
const FVector2d& Vertex = Graph.GetVertex(VertexIdx);
f << "v " << Vertex.X << " " << Vertex.Y << " 0" << std::endl;
}
for (int32 VertexIdx = 0; VertexIdx < Graph.MaxVertexID(); VertexIdx++)
{
const FVector2d& Vertex = Graph.GetVertex(VertexIdx);
f << "v " << Vertex.X << " " << Vertex.Y << " .5" << std::endl;
}
for (const FDynamicGraph::FEdge& Edge : Graph.Edges())
{
f << "f " << Edge.A + 1 << " " << Edge.B + 1 << " " << 1 + Edge.A + Graph.MaxVertexID() << std::endl;
}
f.close();
}
void DumpTriangulationForDebug(const FDynamicGraph2d& Graph, const TArray<FIntVector>& Triangles, const FString& PathBase)
{
num++;
FString Path = PathBase + FString::FromInt(num) + ".obj";
std::ofstream f(*Path);
for (int32 VertexIdx = 0; VertexIdx < Graph.MaxVertexID(); VertexIdx++)
{
const FVector2d& Vertex = Graph.GetVertex(VertexIdx);
f << "v " << Vertex.X << " " << Vertex.Y << " 0" << std::endl;
}
for (const FIntVector& Tri : Triangles)
{
f << "f " << 1 + Tri.X << " " << 1 + Tri.Y << " " << 1 + Tri.Z << std::endl;
}
f.close();
}
#endif
}
bool FArrangement2d::AttemptTriangulate(TArray<FIntVector> &Triangles, TArray<int32> &SkippedEdges, int32 BoundaryEdgeGroupID)
{
Triangles.Empty();
// non-robust version is very non-robust
// Use this declaration for the non-robust version of the algorithm; this is very non-robust, and not recommended!
//gte::ConstrainedDelaunay2<double, double> Delaunay;
// Here are the types you could use for robust math (doesn't make the algorithm 100% robust, but improves it a lot)
// Here are the types you could use for robust math (in theory makes the algorithm overall robust; in practice I am not sure)
gte::ConstrainedDelaunay2<double, gte::BSNumber<gte::UIntegerFP32<263>>> Delaunay; // Value of 263 is from comment in GTEngine/Mathematics/GteDelaunay2.h
//gte::ConstrainedDelaunay2<double, gte::BSNumber<gte::UIntegerAP32>> Delaunay; // Full arbitrary precision (slowest method, creates a std::vector per number)
// Copy vertices to a compact buffer, and create a mapping in and out of the compact space (TODO: maybe extract this for general use? or is it already somewhere else?)
std::vector<gte::Vector2<double>> InputVertices;
TArray<int> InputIndices; InputIndices.SetNumZeroed(Graph.MaxVertexID());
TArray<int> OutputIndices; OutputIndices.SetNumZeroed(Graph.VertexCount());
int TotalOffset = 0;
for (int i = 0; i < Graph.MaxVertexID(); i++)
// If there are unused vertices, define a vertex index remapping so the delaunay code won't have to understand that
bool bNeedsRemap = Graph.MaxVertexID() != Graph.VertexCount();
TArray<int> InputIndices, OutputIndices;
if (bNeedsRemap)
{
if (Graph.IsVertex(i))
InputIndices.SetNumZeroed(Graph.MaxVertexID());
OutputIndices.SetNumZeroed(Graph.VertexCount());
int TotalOffset = 0;
for (int i = 0; i < Graph.MaxVertexID(); i++)
{
OutputIndices[i - TotalOffset] = i;
InputIndices[i] = i - TotalOffset;
FVector2d Vertex = Graph.GetVertex(i);
InputVertices.push_back(gte::Vector2<double> {Vertex.X, Vertex.Y});
if (Graph.IsVertex(i))
{
OutputIndices[i - TotalOffset] = i;
InputIndices[i] = i - TotalOffset;
FVector2d Vertex = Graph.GetVertex(i);
InputVertices.push_back(gte::Vector2<double> {{Vertex.X, Vertex.Y}});
}
else
{
InputIndices[i] = -1;
TotalOffset++;
}
}
else
}
else
{ // no index remapping needed; just copy the vertices
for (int i = 0; i < Graph.MaxVertexID(); i++)
{
InputIndices[i] = -1;
TotalOffset++;
FVector2d Vertex = Graph.GetVertex(i);
InputVertices.push_back(gte::Vector2<double> {{Vertex.X, Vertex.Y}});
}
}
if (!Delaunay(Graph.VertexCount(), &InputVertices[0], 0.001f))
{
return false;
}
std::vector<int> OutEdges; // (needed for Delaunay.Insert call, but otherwise unused)
std::vector<int> OutEdges;
bool InsertConstraintFailure = false;
TSet<TPair<int, int>> BoundarySet; // tracks all the boundary edges as they are added, so we can later eat what's outside them
bool bBoundaryTrackingFailure = false;
for (int EdgeIdx : Graph.EdgeIndices())
{
FDynamicGraph::FEdge Edge = Graph.GetEdge(EdgeIdx);
if (!Delaunay.Insert({{ InputIndices[Edge.A], InputIndices[Edge.B] }}, OutEdges))
if (bNeedsRemap)
{
// Note the failed edge; we will try to proceed anyway, just without this edge
ensureMsgf(false, TEXT("CDT edge insertion failed -- possibly bad data?"));
Edge.A = InputIndices[Edge.A];
Edge.B = InputIndices[Edge.B];
}
if (!Delaunay.Insert({{Edge.A, Edge.B}}, OutEdges))
{
// Note the failed edge; we will try to proceed anyway, just without this edge. Hopefully the CDT is robust and this never happens!
ensureMsgf(false, TEXT("CDT edge insertion failed"));
InsertConstraintFailure = true;
SkippedEdges.Add(EdgeIdx);
// specially note if we failed to add a boundary edge; in this case removing the outside-boundary triangles would likely delete all triangles from the output
// in this case I choose (below) to just add all triangles to the output, ignoring boundaries, as this may create less noticeable artifacts in the typical cases
if (Edge.Group == BoundaryEdgeGroupID)
{
bBoundaryTrackingFailure = true;
}
}
else
{
if (Edge.Group == BoundaryEdgeGroupID)
{
for (size_t i = 0; i + 1 < OutEdges.size(); i++)
{
int MinV = FMath::Min(OutEdges[i], OutEdges[i + 1]);
int MaxV = FMath::Max(OutEdges[i], OutEdges[i + 1]);
BoundarySet.Add(TPair<int, int>(MinV, MaxV));
}
}
}
}
std::vector<int> const& Indices = Delaunay.GetIndices();
for (size_t i = 0, n = Indices.size() / 3; i < n; i++)
const std::vector<int>& Indices = Delaunay.GetIndices();
if (!bBoundaryTrackingFailure && BoundarySet.Num() > 0) {
auto IsBoundaryEdge = [&BoundarySet](int V0, int V1)
{
int MinV = FMath::Min(V0, V1);
int MaxV = FMath::Max(V0, V1);
return BoundarySet.Contains(TPair<int, int>(MinV, MaxV));
};
// eat hull to boundary
const std::vector<int>& Adj = Delaunay.GetAdjacencies();
int TriNum = int(Adj.size() / 3);
TArray<bool> Eaten;
Eaten.SetNumZeroed(TriNum);
TArray<int> ToEatQ;
// seed the feed queue with any triangles that are on hull but w/ an edge that is not on intended boundary
for (int TriIdx = 0; TriIdx < TriNum; TriIdx++)
{
int BaseIdx = TriIdx * 3;
bool bEatIt = false;
for (int SubIdx = 0, NextIdx = 2; !bEatIt && SubIdx < 3; NextIdx = SubIdx++)
{
// eat the triangle if it has a hull edge that is not a boundary edge
bEatIt = Adj[BaseIdx + NextIdx] == -1 && !IsBoundaryEdge(Indices[BaseIdx + SubIdx], Indices[BaseIdx + NextIdx]);
}
if (bEatIt)
{
ToEatQ.Add(TriIdx);
Eaten[TriIdx] = true;
}
}
// eat any triangle that we can get to by crossing a non-boundary edge from an already-been-eaten triangle
while (ToEatQ.Num())
{
int EatTri = ToEatQ.Pop();
int BaseIdx = EatTri * 3;
for (int SubIdx = 0, NextIdx = 2; SubIdx < 3; NextIdx = SubIdx++)
{
int AdjTri = Adj[BaseIdx + NextIdx];
if (AdjTri != -1 && !Eaten[AdjTri] && !IsBoundaryEdge(Indices[BaseIdx + SubIdx], Indices[BaseIdx + NextIdx]))
{
ToEatQ.Add(AdjTri);
Eaten[AdjTri] = true;
}
}
}
// copy uneaten triangles
for (size_t i = 0; i < TriNum; i++)
{
if (!Eaten[i])
{
Triangles.Add(FIntVector(Indices[i * 3], Indices[i * 3 + 1], Indices[i * 3 + 2]));
}
}
ensure(Triangles.Num()); // technically it could be valid to not have any triangles after eating the outside-boundary ones, but it would be unusual and could also be a bug
}
else
{
Triangles.Add(FIntVector(OutputIndices[Indices[i * 3]], OutputIndices[Indices[i * 3 + 1]], OutputIndices[Indices[i * 3 + 2]]));
// copy all triangles over
for (size_t i = 0, n = Indices.size() / 3; i < n; i++)
{
Triangles.Add(FIntVector(Indices[i * 3], Indices[i * 3 + 1], Indices[i * 3 + 2]));
}
}
if (bNeedsRemap)
{
for (FIntVector& Face : Triangles)
{
Face.X = OutputIndices[Face.X];
Face.Y = OutputIndices[Face.Y];
Face.Z = OutputIndices[Face.Z];
}
}
return !InsertConstraintFailure;
}
@@ -183,6 +183,10 @@ template <typename InputType, typename ComputeType>
bool ConstrainedDelaunay2<InputType, ComputeType>::Insert(
std::array<int, 2> const& edge, int v0Triangle, std::vector<int>& outEdge)
{
// TODO: this check helps avoid cases where the algorithm gets stuck in endless recursion walking the mesh (resulting in a stack overflow) ... overflow could still happen for a large mesh, so more work to prevent such an endless loop would be valuable!
// TODO: (and/or: rewrite the insert function to not be recursive!)
GTE_CDT_REQUIRE(outEdge.size() <= this->mComputeVertices.size());
// Create the neighborhood of triangles that share the vertex v0. On
// entry we already know one such triangle (v0Triangle).
int v0 = edge[0], v1 = edge[1];
@@ -563,7 +567,7 @@ ComputeType ConstrainedDelaunay2<InputType, ComputeType>::ComputePSD(int v0,
else
{
dot = DotPerp(V2mV0, V1mV0);
psd = sqrlen10 * dot * dot;
psd = dot * dot;
}
}
@@ -576,6 +580,10 @@ int ConstrainedDelaunay2<InputType, ComputeType>::GetLinkTriangle(int v) const
// Remap in case an edge vertex was specified that is a duplicate.
v = this->mDuplicates[v];
// TODO: this algorithm does not work for point location in a CDT; it's only guaranteed for a (non-C)DT.
// See https://hal.inria.fr/inria-00072509/document for a report on different approaches to walking a triangulation,
// including an example of where this approach (visibility walk) fails.
// For now, I've added a brute-force fallback below. It would be better to choose a correct algorithm to put here though!
int tri = 0;
for (int i = 0; i < this->mNumTriangles; ++i)
{
@@ -607,6 +615,22 @@ int ConstrainedDelaunay2<InputType, ComputeType>::GetLinkTriangle(int v) const
}
}
// fallback to brute force search
// TODO: rm this after fixing the above point location to something that works for all triangulations
for (tri = 0; tri < this->mNumTriangles; ++tri)
{
// Test whether v is a vertex of the triangle.
std::array<int, 3> indices;
GTE_CDT_REQUIRE(this->GetIndices(tri, indices));
for (int j = 0; j < 3; ++j)
{
if (v == indices[j])
{
return tri;
}
}
}
// The vertex must be in the triangulation.
GTE_CDT_FAILURE_RET(-1);
}
@@ -51,9 +51,10 @@ struct FArrangement2d
*
* Triangles: Output triangles (as indices into Graph vertices)
* SkippedEdges: Output indices of edges that the algorithm failed to insert
* BoundaryEdgeGroupID: ID of edges corresponding to a boundary; if we have a closed loop of these boundary edges on output triangulation, will discard triangles outside this
* return: false if triangulation algo knows it failed (NOTE: triangulation may still be invalid when function returns true, as function is not robust)
*/
bool AttemptTriangulate(TArray<FIntVector>& Triangles, TArray<int32>& SkippedEdges);
bool GEOMETRYALGORITHMS_API AttemptTriangulate(TArray<FIntVector>& Triangles, TArray<int32>& SkippedEdges, int32 BoundaryEdgeGroupID);
/**
* Check if current Graph has self-intersections; not optimized, only for debugging
@@ -155,7 +156,7 @@ protected:
/**
* insert edge [A,B] into the arrangement, splitting existing edges as necessary
*/
bool insert_segment(const FVector2d& A, const FVector2d& B, int GID = -1, double Tol = 0)
bool insert_segment(FVector2d A, FVector2d B, int GID = -1, double Tol = 0)
{
// handle degenerate edges
int a_idx = find_existing_vertex(A);
@@ -164,6 +165,15 @@ protected:
{
return false;
}
// snap input vertices
if (a_idx >= 0)
{
A = Graph.GetVertex(a_idx);
}
if (b_idx >= 0)
{
B = Graph.GetVertex(b_idx);
}
// ok find all intersections
TArray<FIntersection> Hits;