Skip to content

Commit

Permalink
Cleaning Hull Algorithm Code (#904)
Browse files Browse the repository at this point in the history 
* Removing redundant code (structs and function calls)

* Removed unecessary comments

* Removed ConvexHull, CCW and getConvexHullAsMesh

* Moved Impl::Hull to quickhull.cpp

* Modified some comments
  • Loading branch information
@Kushal-Shah-03
Kushal-Shah-03 authored Aug 25, 2024
1 parent db4470c commit a0e4d1b
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 119 deletions.
29 changes: 0 additions & 29 deletions src/manifold/src/impl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -478,35 +478,6 @@ Manifold::Impl::Impl(const Mesh& mesh, const MeshRelationD& relation,
Finish();
}

void Manifold::Impl::Hull(const std::vector<glm::vec3>& vertPos) {
size_t numVert = vertPos.size();
if (numVert < 4) {
status_ = Error::InvalidConstruction;
return;
}

Vec<glm::dvec3> pointCloudVec(numVert);
manifold::transform(vertPos.begin(), vertPos.end(), pointCloudVec.begin(),
[](const glm::vec3& v) { return glm::dvec3(v); });
QuickHull qh(pointCloudVec);
ConvexHull hull = qh.getConvexHullAsMesh(pointCloudVec, false);
// TODO : Once double PR lands, replace this with move
vertPos_.resize(hull.vertices.size());
manifold::transform(hull.vertices.begin(), hull.vertices.end(),
vertPos_.begin(),
[](const glm::dvec3& v) { return glm::vec3(v); });
halfedge_ = std::move(hull.halfedges);
meshRelation_.originalID = ReserveIDs(1);
CalculateBBox();
SetPrecision(bBox_.Scale() * kTolerance);
SplitPinchedVerts();
CalculateNormals();
InitializeOriginal();
CreateFaces({});
SimplifyTopology();
Finish();
}

/**
* Create either a unit tetrahedron, cube or octahedron. The cube is in the
* first octant, while the others are symmetric about the origin.
Expand Down
82 changes: 53 additions & 29 deletions src/manifold/src/quickhull.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,12 @@
#include <algorithm>
#include <limits>

#include "impl.h"

namespace manifold {

double defaultEps() { return 0.0000001; }

// MathUtils.hpp
inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p,
const Ray& r) {
const glm::dvec3 s = p - r.S;
Expand Down Expand Up @@ -162,8 +163,8 @@ HalfEdgeMesh::HalfEdgeMesh(const MeshBuilder& builderObject,
size_t i = 0;
for (const auto& face : builderObject.faces) {
if (!face.isDisabled()) {
faces.emplace_back(static_cast<size_t>(face.he));
faceMapping[i] = faces.size() - 1;
halfEdgeIndexFaces.emplace_back(static_cast<size_t>(face.he));
faceMapping[i] = halfEdgeIndexFaces.size() - 1;

const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face);
for (const auto heIndex : heIndices) {
Expand All @@ -188,10 +189,10 @@ HalfEdgeMesh::HalfEdgeMesh(const MeshBuilder& builderObject,
i++;
}

for (auto& face : faces) {
ASSERT(halfEdgeMapping.count(face.halfEdgeIndex) == 1,
for (auto& halfEdgeIndexFace : halfEdgeIndexFaces) {
ASSERT(halfEdgeMapping.count(halfEdgeIndexFace) == 1,
logicErr("invalid halfedge mapping"));
face.halfEdgeIndex = halfEdgeMapping[face.halfEdgeIndex];
halfEdgeIndexFace = halfEdgeMapping[halfEdgeIndexFace];
}

for (size_t i = 0; i < halfedges.size(); i++) {
Expand All @@ -206,13 +207,11 @@ HalfEdgeMesh::HalfEdgeMesh(const MeshBuilder& builderObject,
/*
* Implementation of the algorithm
*/

ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
double epsilon) {
if (pointCloud.size() == 0) {
return ConvexHull();
// TODO : Once double PR lands replace glm::vec3 with glm::dvec3
std::pair<Vec<Halfedge>, Vec<glm::vec3>> QuickHull::buildMesh(double epsilon) {
if (originalVertexData.size() == 0) {
return {Vec<Halfedge>(), Vec<glm::vec3>()};
}
originalVertexData = pointCloud;

// Very first: find extreme values and use them to compute the scale of the
// point cloud.
Expand All @@ -223,8 +222,6 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
m_epsilon = epsilon * scale;
epsilonSquared = m_epsilon * m_epsilon;

// Reset diagnostics
diagnostics = DiagnosticsData();
// The planar case happens when all the points appear to lie on a two
// dimensional subspace of R^3.
planar = false;
Expand All @@ -236,7 +233,6 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
he.endVert = 0;
}
}
originalVertexData = pointCloud;
planarPointCloudTemp.clear();
}

Expand All @@ -246,9 +242,8 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
Vec<int> mapping(mesh.halfedges.size());
Vec<int> faceMap(mesh.faces.size());

// reorder halfedges
// Since all faces are not used now (so we should start from 0 and just
// increment, we can later set the face id as index/3 for each halfedge
// Some faces are disabled and should not go into the halfedge vector, we can
// update the face indices of the halfedges at the end using index/3
int j = 0;
for_each(
autoPolicy(mesh.halfedges.size()), countAt(0_uz),
Expand Down Expand Up @@ -276,7 +271,7 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
for_each(
autoPolicy(halfedges.size()), halfedges.begin(), halfedges.end(),
[&](Halfedge& he) { he.pairedHalfedge = mapping[he.pairedHalfedge]; });
counts.resize(pointCloud.size() + 1);
counts.resize(originalVertexData.size() + 1);
fill(counts.begin(), counts.end(), 0);

// remove unused vertices
Expand All @@ -290,10 +285,10 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
exclusive_scan(TransformIterator(counts.begin(), saturate),
TransformIterator(counts.end(), saturate), counts.begin(), 0);
Vec<glm::dvec3> vertices(counts.back());
for_each(autoPolicy(pointCloud.size()), countAt(0_uz),
countAt(pointCloud.size()), [&](size_t i) {
for_each(autoPolicy(originalVertexData.size()), countAt(0_uz),
countAt(originalVertexData.size()), [&](size_t i) {
if (counts[i + 1] - counts[i] > 0) {
vertices[counts[i]] = pointCloud[i];
vertices[counts[i]] = originalVertexData[i];
}
});
for_each(autoPolicy(halfedges.size()), halfedges.begin(), halfedges.end(),
Expand All @@ -305,7 +300,11 @@ ConvexHull QuickHull::buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
for (size_t index = 0; index < halfedges.size(); index++) {
halfedges[index].face = index / 3;
}
return ConvexHull{std::move(halfedges), std::move(vertices)};
// TODO : Once double PR lands remove this code
Vec<glm::vec3> verticesFloat(vertices.size());
manifold::transform(vertices.begin(), vertices.end(), verticesFloat.begin(),
[](const glm::dvec3& v) { return glm::vec3(v); });
return {std::move(halfedges), std::move(verticesFloat)};
}

void QuickHull::createConvexHalfedgeMesh() {
Expand Down Expand Up @@ -417,7 +416,7 @@ void QuickHull::createConvexHalfedgeMesh() {
// numerical instability in which case we give up trying to solve horizon
// edge for this point and accept a minor degeneration in the convex hull.
if (!reorderHorizonEdges(horizonEdgesData)) {
diagnostics.failedHorizonEdges++;
failedHorizonEdges++;
int change_flag = 0;
for (size_t index = 0; index < tf.pointsOnPositiveSide->size(); index++) {
if ((*tf.pointsOnPositiveSide)[index] == activePointIndex) {
Expand Down Expand Up @@ -772,11 +771,10 @@ void QuickHull::setupInitialTetrahedron() {
mesh.setup(baseTriangle[0], baseTriangle[1], baseTriangle[2], maxI);
for (auto& f : mesh.faces) {
auto v = mesh.getVertexIndicesOfFace(f);
const glm::dvec3& va = originalVertexData[v[0]];
const glm::dvec3& vb = originalVertexData[v[1]];
const glm::dvec3& vc = originalVertexData[v[2]];
const glm::dvec3 N1 = getTriangleNormal(va, vb, vc);
const Plane plane(N1, va);
const glm::dvec3 N1 =
getTriangleNormal(originalVertexData[v[0]], originalVertexData[v[1]],
originalVertexData[v[2]]);
const Plane plane(N1, originalVertexData[v[0]]);
f.P = plane;
}

Expand Down Expand Up @@ -826,4 +824,30 @@ bool QuickHull::addPointToFace(typename MeshBuilder::Face& f,
}
return false;
}

// Wrapper to call the QuickHull algorithm with the given vertex data to build
// the Impl
void Manifold::Impl::Hull(const std::vector<glm::vec3>& vertPos) {
size_t numVert = vertPos.size();
if (numVert < 4) {
status_ = Error::InvalidConstruction;
return;
}

Vec<glm::dvec3> pointCloudVec(numVert);
manifold::transform(vertPos.begin(), vertPos.end(), pointCloudVec.begin(),
[](const glm::vec3& v) { return glm::dvec3(v); });
QuickHull qh(pointCloudVec);
std::tie(halfedge_, vertPos_) = qh.buildMesh();
meshRelation_.originalID = ReserveIDs(1);
CalculateBBox();
SetPrecision(bBox_.Scale() * kTolerance);
SplitPinchedVerts();
CalculateNormals();
InitializeOriginal();
CreateFaces({});
SimplifyTopology();
Finish();
}

} // namespace manifold
74 changes: 13 additions & 61 deletions src/manifold/src/quickhull.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,13 +63,6 @@

namespace manifold {

class ConvexHull {
public:
Vec<Halfedge> halfedges;
Vec<glm::dvec3> vertices;
};

// Pool.hpp
class Pool {
std::vector<std::unique_ptr<Vec<size_t>>> data;

Expand All @@ -91,8 +84,6 @@ class Pool {
}
};

// Plane.hpp

class Plane {
public:
glm::dvec3 N;
Expand All @@ -116,8 +107,6 @@ class Plane {
: N(N), D(glm::dot(-N, P)), sqrNLength(glm::dot(N, N)) {}
};

// Ray.hpp

struct Ray {
const glm::dvec3 S;
const glm::dvec3 V;
Expand All @@ -127,8 +116,6 @@ struct Ray {
: S(S), V(V), VInvLengthSquared(1 / (glm::dot(V, V))) {}
};

// Mesh.hpp

class MeshBuilder {
public:
struct Face {
Expand Down Expand Up @@ -209,36 +196,18 @@ class MeshBuilder {
}
};

// HalfEdgeMesh.hpp

class HalfEdgeMesh {
public:
struct Face {
// Index of one of the half edges of this face
size_t halfEdgeIndex;

Face(size_t halfEdgeIndex) : halfEdgeIndex(halfEdgeIndex) {}
};

Vec<glm::dvec3> vertices;
std::vector<Face> faces;
// Index of one of the half edges of the faces
std::vector<size_t> halfEdgeIndexFaces;
Vec<Halfedge> halfedges;
Vec<int> halfedgeNext;

HalfEdgeMesh(const MeshBuilder& builderObject,
const VecView<glm::dvec3>& vertexData);
};

// QuickHull.hpp

struct DiagnosticsData {
// How many times QuickHull failed to solve the horizon edge. Failures lead
// to degenerated convex hulls.
size_t failedHorizonEdges;

DiagnosticsData() : failedHorizonEdges(0) {}
};

double defaultEps();

class QuickHull {
Expand All @@ -247,8 +216,6 @@ class QuickHull {
// If the face turns out not to be visible, this half edge will be marked as
// horizon edge
int enteredFromHalfedge;
FaceData() {}
FaceData(int fi, int he) : faceIndex(fi), enteredFromHalfedge(he) {}
};
using vec3 = glm::dvec3;

Expand All @@ -258,7 +225,7 @@ class QuickHull {
VecView<vec3> originalVertexData;
MeshBuilder mesh;
std::array<size_t, 6> extremeValues;
DiagnosticsData diagnostics;
size_t failedHorizonEdges = 0;

// Temporary variables used during iteration process
Vec<size_t> newFaceIndices;
Expand Down Expand Up @@ -300,36 +267,21 @@ class QuickHull {
// the plane. Returns true if the points was on the positive side.
inline bool addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex);

// This will update mesh from which we create the ConvexHull object that
// getConvexHull function returns
// This will create HalfedgeMesh from which we create the ConvexHull object
// that buildMesh function returns
void createConvexHalfedgeMesh();

// Constructs the convex hull into halfEdges and NewVerts
ConvexHull buildMesh(const VecView<glm::dvec3>& pointCloud, bool CCW,
double eps);

public:
// This function assumes that the pointCloudVec data resides in memory in the
// following format: x_0,y_0,z_0,x_1,y_1,z_1,...
QuickHull(const Vec<glm::dvec3>& pointCloudVec)
: originalVertexData(VecView(pointCloudVec)) {}

// Computes convex hull for a given point cloud. This function assumes that
// the vertex data resides in memory in the following format:
// x_0,y_0,z_0,x_1,y_1,z_1,... Params:
// vertexData: pointer to the X component of the first point of the point
// cloud. vertexCount: number of vertices in the point cloud CCW: whether
// the output mesh triangles should have CCW orientation eps: minimum
// distance to a plane to consider a point being on positive side of it (for
// a point cloud with scale 1)
// Returns:
// Convex hull of the point cloud as halfEdge vector and vertex vector
ConvexHull getConvexHullAsMesh(const Vec<glm::dvec3>& pointCloud, bool CCW,
double epsilon = defaultEps()) {
Vec<glm::dvec3> pointCloudVec(pointCloud);
QuickHull qh(pointCloudVec);
return qh.buildMesh(pointCloudVec, CCW, epsilon);
}

// Get diagnostics about last generated convex hull
const DiagnosticsData& getDiagnostics() { return diagnostics; }
// Computes convex hull for a given point cloud. Params: eps: minimum distance
// to a plane to consider a point being on positive side of it (for a point
// cloud with scale 1) Returns: Convex hull of the point cloud as halfEdge
// vector and vertex vector
std::pair<Vec<Halfedge>, Vec<glm::vec3>> buildMesh(double eps = defaultEps());
};

} // namespace manifold

0 comments on commit a0e4d1b

Please sign in to comment.