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

Jump to bottom

Better faces #1040

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Better faces #1040

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 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
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?

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) {
Copy link
Collaborator

Choose a reason for hiding this comment

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

Considering h = NextHalfedge(halfedge_[interiorHalfedges.back()].pairedHalfedge), will this condition ever hold?

interiorHalfedges.pop_back();
} else {
interiorHalfedges.push_back(h);
}
const int hNext = NextHalfedge(h);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Why are we also working with the front of the queue?

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