Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Better faces #1040

Open
wants to merge 13 commits into
base: master
Choose a base branch
from
+209 −335
Open

Better faces #1040

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
6767b42
refactor
elalish Nov 13, 2024
63fb341
separate propVert from createFaces
elalish Nov 13, 2024
1c68d76
working, but very slow
elalish Nov 13, 2024
544696b
no longer slow
elalish Nov 13, 2024
b2358dd
added test
elalish Nov 14, 2024
c4ce729
forgot include
elalish Nov 14, 2024
87a19b5
Merge branch 'master' of github.com:elalish/manifold into betterFaces
elalish Nov 24, 2024
a88f646
don't special-case short edges
elalish Dec 3, 2024
7746a38
complete collapse pass
elalish Dec 3, 2024
a45ff30
remove deque
elalish Dec 3, 2024
9984fcd
merging master
elalish Dec 4, 2024
3992a46
Merge branch 'master' of github.com:elalish/manifold into betterFaces
elalish Dec 17, 2024
41f4f8d
some progress
elalish Dec 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
263 changes: 95 additions & 168 deletions src/impl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@

#include <algorithm>
#include <atomic>
#include <deque>
#include <map>

#include "./hashtable.h"
Expand Down Expand Up @@ -90,98 +91,6 @@ struct UpdateMeshID {
void operator()(TriRef& ref) { ref.meshID = meshIDold2new[ref.meshID]; }
};

struct CoplanarEdge {
VecView<std::pair<int, int>> face2face;
VecView<double> triArea;
VecView<const Halfedge> halfedge;
VecView<const vec3> vertPos;
VecView<const TriRef> triRef;
VecView<const ivec3> triProp;
const int numProp;
const double epsilon;
const double tolerance;

void operator()(const int edgeIdx) {
const Halfedge edge = halfedge[edgeIdx];
const Halfedge pair = halfedge[edge.pairedHalfedge];
const int edgeFace = edgeIdx / 3;
const int pairFace = edge.pairedHalfedge / 3;

if (triRef[edgeFace].meshID != triRef[pairFace].meshID) return;

const vec3 base = vertPos[edge.startVert];
const int baseNum = edgeIdx - 3 * edgeFace;
const int jointNum = edge.pairedHalfedge - 3 * pairFace;

if (numProp > 0) {
if (triProp[edgeFace][baseNum] != triProp[pairFace][Next3(jointNum)] ||
triProp[edgeFace][Next3(baseNum)] != triProp[pairFace][jointNum])
return;
}

if (!edge.IsForward()) return;

const int edgeNum = baseNum == 0 ? 2 : baseNum - 1;
const int pairNum = jointNum == 0 ? 2 : jointNum - 1;
const vec3 jointVec = vertPos[pair.startVert] - base;
const vec3 edgeVec =
vertPos[halfedge[3 * edgeFace + edgeNum].startVert] - base;
const vec3 pairVec =
vertPos[halfedge[3 * pairFace + pairNum].startVert] - base;

const double length = std::max(la::length(jointVec), la::length(edgeVec));
const double lengthPair =
std::max(la::length(jointVec), la::length(pairVec));
vec3 normal = la::cross(jointVec, edgeVec);
const double area = la::length(normal);
const double areaPair = la::length(la::cross(pairVec, jointVec));

// make sure we only write this once
if (edgeIdx % 3 == 0) triArea[edgeFace] = area;
// Don't link degenerate triangles
if (area < length * epsilon || areaPair < lengthPair * epsilon) return;

const double volume = std::abs(la::dot(normal, pairVec));
// Only operate on coplanar triangles
if (volume > std::max(area, areaPair) * tolerance) return;

face2face[edgeIdx] = std::make_pair(edgeFace, pairFace);
}
};

struct CheckCoplanarity {
VecView<int> comp2tri;
VecView<const Halfedge> halfedge;
VecView<const vec3> vertPos;
std::vector<int>* components;
const double tolerance;

void operator()(int tri) {
const int component = (*components)[tri];
const int referenceTri =
reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
->load(std::memory_order_relaxed);
if (referenceTri < 0 || referenceTri == tri) return;

const vec3 origin = vertPos[halfedge[3 * referenceTri].startVert];
const vec3 normal = la::normalize(
la::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin,
vertPos[halfedge[3 * referenceTri + 2].startVert] - origin));

for (const int i : {0, 1, 2}) {
const vec3 vert = vertPos[halfedge[3 * tri + i].startVert];
// If any component vertex is not coplanar with the component's reference
// triangle, unmark the entire component so that none of its triangles are
// marked coplanar.
if (std::abs(la::dot(normal, vert - origin)) > tolerance) {
reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
->store(-1, std::memory_order_relaxed);
break;
}
}
}
};

int GetLabels(std::vector<int>& components,
const Vec<std::pair<int, int>>& edges, int numNodes) {
UnionFind<> uf(numNodes);
Expand All @@ -192,18 +101,6 @@ int GetLabels(std::vector<int>& components,

return uf.connectedComponents(components);
}

void DedupePropVerts(manifold::Vec<ivec3>& triProp,
const Vec<std::pair<int, int>>& vert2vert) {
ZoneScoped;
std::vector<int> vertLabels;
const int numLabels = GetLabels(vertLabels, vert2vert, vert2vert.size());

std::vector<int> label2vert(numLabels);
for (size_t v = 0; v < vert2vert.size(); ++v) label2vert[vertLabels[v]] = v;
for (auto& prop : triProp)
for (int i : {0, 1, 2}) prop[i] = label2vert[vertLabels[prop[i]]];
}
} // namespace

namespace manifold {
Expand Down Expand Up @@ -318,75 +215,105 @@ void Manifold::Impl::InitializeOriginal(bool keepFaceID) {

void Manifold::Impl::CreateFaces() {
ZoneScoped;
Vec<std::pair<int, int>> face2face(halfedge_.size(), {-1, -1});
Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
Vec<double> triArea(NumTri());

const size_t numProp = NumProp();
if (numProp > 0) {
for_each_n(
autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
[&vert2vert, numProp, this](const int edgeIdx) {
const Halfedge edge = halfedge_[edgeIdx];
const Halfedge pair = halfedge_[edge.pairedHalfedge];
const int edgeFace = edgeIdx / 3;
const int pairFace = edge.pairedHalfedge / 3;

if (meshRelation_.triRef[edgeFace].meshID !=
meshRelation_.triRef[pairFace].meshID)
return;

const int baseNum = edgeIdx - 3 * edgeFace;
const int jointNum = edge.pairedHalfedge - 3 * pairFace;

const int prop0 = meshRelation_.triProperties[edgeFace][baseNum];
const int prop1 =
meshRelation_
.triProperties[pairFace][jointNum == 2 ? 0 : jointNum + 1];
bool propEqual = true;
for (size_t p = 0; p < numProp; ++p) {
if (meshRelation_.properties[numProp * prop0 + p] !=
meshRelation_.properties[numProp * prop1 + p]) {
propEqual = false;
break;
}
}
if (propEqual) {
vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
}
});
DedupePropVerts(meshRelation_.triProperties, vert2vert);
}
const int numTri = NumTri();
struct TriPriority {
double area2;
int tri;
};
Vec<TriPriority> triPriority(numTri);
for_each_n(autoPolicy(numTri), countAt(0), numTri,
[&triPriority, this](int tri) {
meshRelation_.triRef[tri].faceID = -1;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this causing the check in line 238 and 249 to always be false?

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, because they get reassigned on lines 240 and 253.

const vec3 v = vertPos_[halfedge_[3 * tri].startVert];
triPriority[tri] = {
length2(cross(vertPos_[halfedge_[3 * tri].endVert] - v,
vertPos_[halfedge_[3 * tri + 1].endVert] - v)),
tri};
});

for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
CoplanarEdge({face2face, triArea, halfedge_, vertPos_,
meshRelation_.triRef, meshRelation_.triProperties,
meshRelation_.numProp, epsilon_, tolerance_}));

std::vector<int> components;
const int numComponent = GetLabels(components, face2face, NumTri());

Vec<int> comp2tri(numComponent, -1);
for (size_t tri = 0; tri < NumTri(); ++tri) {
const int comp = components[tri];
const int current = comp2tri[comp];
if (current < 0 || triArea[tri] > triArea[current]) {
comp2tri[comp] = tri;
triArea[comp] = triArea[tri];
stable_sort(triPriority.begin(), triPriority.end(),
[](auto a, auto b) { return a.area2 > b.area2; });

for (const auto tp : triPriority) {
if (meshRelation_.triRef[tp.tri].faceID >= 0) continue;

meshRelation_.triRef[tp.tri].faceID = tp.tri;
const vec3 base = vertPos_[halfedge_[3 * tp.tri].startVert];
const vec3 normal = faceNormal_[tp.tri];
std::deque<int> interiorHalfedges = {3 * tp.tri, 3 * tp.tri + 1,
3 * tp.tri + 2};
while (!interiorHalfedges.empty()) {
const int h =
NextHalfedge(halfedge_[interiorHalfedges.back()].pairedHalfedge);
interiorHalfedges.pop_back();
if (meshRelation_.triRef[h / 3].faceID >= 0) continue;

const vec3 v = vertPos_[halfedge_[h].endVert];
if (abs(dot(v - base, normal)) < tolerance_) {
meshRelation_.triRef[h / 3].faceID = tp.tri;
if (!interiorHalfedges.empty() &&
h == halfedge_[interiorHalfedges.back()].pairedHalfedge) {
elalish marked this conversation as resolved.
Show resolved Hide resolved
interiorHalfedges.pop_back();
} else {
interiorHalfedges.push_back(h);
}
const int hNext = NextHalfedge(h);
elalish marked this conversation as resolved.
Show resolved Hide resolved
if (!interiorHalfedges.empty() &&
hNext == halfedge_[interiorHalfedges.front()].pairedHalfedge) {
interiorHalfedges.pop_front();
} else {
interiorHalfedges.push_back(hNext);
}
}
}
}
}

for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), NumTri(),
CheckCoplanarity(
{comp2tri, halfedge_, vertPos_, &components, tolerance_}));
void Manifold::Impl::DedupePropVerts() {
ZoneScoped;
const size_t numProp = NumProp();
if (numProp == 0) return;

Vec<TriRef>& triRef = meshRelation_.triRef;
for (size_t tri = 0; tri < NumTri(); ++tri) {
const int referenceTri = comp2tri[components[tri]];
if (referenceTri >= 0) {
triRef[tri].faceID = referenceTri;
}
}
Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
for_each_n(
autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
[&vert2vert, numProp, this](const int edgeIdx) {
const Halfedge edge = halfedge_[edgeIdx];
const Halfedge pair = halfedge_[edge.pairedHalfedge];
const int edgeFace = edgeIdx / 3;
const int pairFace = edge.pairedHalfedge / 3;

if (meshRelation_.triRef[edgeFace].meshID !=
meshRelation_.triRef[pairFace].meshID)
return;

const int baseNum = edgeIdx - 3 * edgeFace;
const int jointNum = edge.pairedHalfedge - 3 * pairFace;

const int prop0 = meshRelation_.triProperties[edgeFace][baseNum];
const int prop1 =
meshRelation_
.triProperties[pairFace][jointNum == 2 ? 0 : jointNum + 1];
bool propEqual = true;
for (size_t p = 0; p < numProp; ++p) {
if (meshRelation_.properties[numProp * prop0 + p] !=
meshRelation_.properties[numProp * prop1 + p]) {
propEqual = false;
break;
}
}
if (propEqual) {
vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
}
});

std::vector<int> vertLabels;
const int numLabels = GetLabels(vertLabels, vert2vert, vert2vert.size());

std::vector<int> label2vert(numLabels);
for (size_t v = 0; v < vert2vert.size(); ++v) label2vert[vertLabels[v]] = v;
for (auto& prop : meshRelation_.triProperties)
for (int i : {0, 1, 2}) prop[i] = label2vert[vertLabels[prop[i]]];
}

/**
Expand Down
2 changes: 2 additions & 0 deletions src/impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,7 @@ struct Manifold::Impl {
InitializeOriginal();
}

DedupePropVerts();
CreateFaces();

SimplifyTopology();
Expand Down Expand Up @@ -246,6 +247,7 @@ struct Manifold::Impl {
}

void CreateFaces();
void DedupePropVerts();
void RemoveUnreferencedVerts();
void InitializeOriginal(bool keepFaceID = false);
void CreateHalfedges(const Vec<ivec3>& triVerts);
Expand Down
18 changes: 16 additions & 2 deletions test/properties_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@ TEST(Properties, Tolerance) {
Manifold imperfect3(mesh);

EXPECT_EQ(imperfect.NumTri(), 28);
EXPECT_EQ(imperfect2.NumTri(), 16); // TODO: should be 12
EXPECT_EQ(imperfect3.NumTri(), 22); // TODO: should be 12
EXPECT_EQ(imperfect2.NumTri(), 12);
EXPECT_EQ(imperfect3.NumTri(), 12);

EXPECT_NEAR(imperfect.Volume(), imperfect2.Volume(), 0.01);
EXPECT_NEAR(imperfect.SurfaceArea(), imperfect2.SurfaceArea(), 0.02);
Expand All @@ -76,6 +76,20 @@ TEST(Properties, Tolerance) {
#endif
}

TEST(Properties, ToleranceSphere) {
const int n = 100;
Manifold sphere = Manifold::Sphere(1, 4 * n);
EXPECT_EQ(sphere.NumTri(), 8 * n * n);

Manifold sphere2 = sphere.SetTolerance(0.01);
EXPECT_EQ(sphere2.NumTri(), 16786);
EXPECT_NEAR(sphere.Volume(), sphere2.Volume(), 0.02);
EXPECT_NEAR(sphere.SurfaceArea(), sphere2.SurfaceArea(), 0.01);
#ifdef MANIFOLD_EXPORT
if (options.exportModels) ExportMesh("sphere.glb", sphere2.GetMeshGL(), {});
#endif
}

/**
* Curvature is the inverse of the radius of curvature, and signed such that
* positive is convex and negative is concave. There are two orthogonal
Expand Down
Loading