From 9d013b3fdaaecc34913ca8f5801b25edcdcb9b09 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 25 Jul 2024 17:52:34 +0530 Subject: [PATCH 01/24] Removed QuickHull Dependency --- bindings/python/third_party/nanobind | 2 +- src/manifold/CMakeLists.txt | 2 +- src/manifold/src/manifold.cpp | 2 +- src/manifold/src/quickhull.cpp | 503 ++++++++++++++ src/manifold/src/quickhull.hpp | 981 +++++++++++++++++++++++++++ src/third_party/CMakeLists.txt | 7 - src/third_party/quickhull | 1 - 7 files changed, 1487 insertions(+), 11 deletions(-) create mode 100644 src/manifold/src/quickhull.cpp create mode 100644 src/manifold/src/quickhull.hpp delete mode 160000 src/third_party/quickhull diff --git a/bindings/python/third_party/nanobind b/bindings/python/third_party/nanobind index 8d7f1ee06..1a309ba44 160000 --- a/bindings/python/third_party/nanobind +++ b/bindings/python/third_party/nanobind @@ -1 +1 @@ -Subproject commit 8d7f1ee0621c17fa370b704b2100ffa0243d5bfb +Subproject commit 1a309ba444a47e081dc6213d72345a2fbbd20795 diff --git a/src/manifold/CMakeLists.txt b/src/manifold/CMakeLists.txt index d0405bf91..5e50295ed 100644 --- a/src/manifold/CMakeLists.txt +++ b/src/manifold/CMakeLists.txt @@ -22,7 +22,7 @@ target_include_directories(${PROJECT_NAME} PUBLIC $) target_link_libraries(${PROJECT_NAME} PUBLIC utilities sdf - PRIVATE collider polygon ${MANIFOLD_INCLUDE} quickhull + PRIVATE collider polygon ${MANIFOLD_INCLUDE} ) target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 07d4b55d4..34af207b9 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -16,7 +16,7 @@ #include #include -#include "QuickHull.hpp" +#include "quickhull.hpp" #include "boolean3.h" #include "csg_tree.h" #include "impl.h" diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp new file mode 100644 index 000000000..e291ea5e2 --- /dev/null +++ b/src/manifold/src/quickhull.cpp @@ -0,0 +1,503 @@ +#include "quickhull.hpp" +// #include "MathUtils.hpp" +#include +#include +#include +#include +#include +// #include "Structs/Mesh.hpp" + +namespace quickhull { + + template<> + float defaultEps() { + return 0.0001f; + } + + template<> + double defaultEps() { + return 0.0000001; + } + + /* + * Implementation of the algorithm + */ + + template + ConvexHull QuickHull::getConvexHull(const std::vector>& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { + VertexDataSource vertexDataSource(pointCloud); + return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); + } + + template + ConvexHull QuickHull::getConvexHull(const Vector3* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) { + VertexDataSource vertexDataSource(vertexData,vertexCount); + return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); + } + + template + ConvexHull QuickHull::getConvexHull(const T* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) { + VertexDataSource vertexDataSource((const vec3*)vertexData,vertexCount); + return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); + } + + template + HalfEdgeMesh QuickHull::getConvexHullAsMesh(const FloatType* vertexData, size_t vertexCount, bool CCW, FloatType epsilon) { + VertexDataSource vertexDataSource((const vec3*)vertexData,vertexCount); + buildMesh(vertexDataSource, CCW, false, epsilon); + return HalfEdgeMesh(m_mesh, m_vertexData); + } + + template + void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { + // CCW is unused for now + (void)CCW; + // useOriginalIndices is unused for now + (void)useOriginalIndices; + + if (pointCloud.size()==0) { + m_mesh = MeshBuilder(); + return; + } + m_vertexData = pointCloud; + + // Very first: find extreme values and use them to compute the scale of the point cloud. + m_extremeValues = getExtremeValues(); + m_scale = getScale(m_extremeValues); + + // Epsilon we use depends on the scale + m_epsilon = epsilon*m_scale; + m_epsilonSquared = m_epsilon*m_epsilon; + + // Reset diagnostics + m_diagnostics = DiagnosticsData(); + + m_planar = false; // The planar case happens when all the points appear to lie on a two dimensional subspace of R^3. + createConvexHalfEdgeMesh(); + if (m_planar) { + const size_t extraPointIndex = m_planarPointCloudTemp.size()-1; + for (auto& he : m_mesh.m_halfEdges) { + if (he.m_endVertex == extraPointIndex) { + he.m_endVertex = 0; + } + } + m_vertexData = pointCloud; + m_planarPointCloudTemp.clear(); + } + } + + template + ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { + buildMesh(pointCloud,CCW,useOriginalIndices,epsilon); + return ConvexHull(m_mesh,m_vertexData, CCW, useOriginalIndices); + } + + template + void QuickHull::createConvexHalfEdgeMesh() { + m_visibleFaces.clear(); + m_horizonEdges.clear(); + m_possiblyVisibleFaces.clear(); + + // Compute base tetrahedron + setupInitialTetrahedron(); + assert(m_mesh.m_faces.size()==4); + + // Init face stack with those faces that have points assigned to them + m_faceList.clear(); + for (size_t i=0;i < 4;i++) { + auto& f = m_mesh.m_faces[i]; + if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size()>0) { + m_faceList.push_back(i); + f.m_inFaceStack = 1; + } + } + + // Process faces until the face list is empty. + size_t iter = 0; + while (!m_faceList.empty()) { + iter++; + if (iter == std::numeric_limits::max()) { + // Visible face traversal marks visited faces with iteration counter (to mark that the face has been visited on this iteration) and the max value represents unvisited faces. At this point we have to reset iteration counter. This shouldn't be an + // issue on 64 bit machines. + iter = 0; + } + + const size_t topFaceIndex = m_faceList.front(); + m_faceList.pop_front(); + + auto& tf = m_mesh.m_faces[topFaceIndex]; + tf.m_inFaceStack = 0; + + assert(!tf.m_pointsOnPositiveSide || tf.m_pointsOnPositiveSide->size()>0); + if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) { + continue; + } + + // Pick the most distant point to this triangle plane as the point to which we extrude + const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint]; + const size_t activePointIndex = tf.m_mostDistantPoint; + + // Find out the faces that have our active point on their positive side (these are the "visible faces"). The face on top of the stack of course is one of them. At the same time, we create a list of horizon edges. + m_horizonEdges.clear(); + m_possiblyVisibleFaces.clear(); + m_visibleFaces.clear(); + m_possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits::max()); + while (m_possiblyVisibleFaces.size()) { + const auto faceData = m_possiblyVisibleFaces.back(); + m_possiblyVisibleFaces.pop_back(); + auto& pvf = m_mesh.m_faces[faceData.m_faceIndex]; + assert(!pvf.isDisabled()); + + if (pvf.m_visibilityCheckedOnIteration == iter) { + if (pvf.m_isVisibleFaceOnCurrentIteration) { + continue; + } + } + else { + const Plane& P = pvf.m_P; + pvf.m_visibilityCheckedOnIteration = iter; + const T d = P.m_N.dotProduct(activePoint)+P.m_D; + if (d>0) { + pvf.m_isVisibleFaceOnCurrentIteration = 1; + pvf.m_horizonEdgesOnCurrentIteration = 0; + m_visibleFaces.push_back(faceData.m_faceIndex); + for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) { + if (m_mesh.m_halfEdges[heIndex].m_opp != faceData.m_enteredFromHalfEdge) { + m_possiblyVisibleFaces.emplace_back( m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex ); + } + } + continue; + } + assert(faceData.m_faceIndex != topFaceIndex); + } + + // The face is not visible. Therefore, the halfedge we came from is part of the horizon edge. + pvf.m_isVisibleFaceOnCurrentIteration = 0; + m_horizonEdges.push_back(faceData.m_enteredFromHalfEdge); + // Store which half edge is the horizon edge. The other half edges of the face will not be part of the final mesh so their data slots can by recycled. + const auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]); + const std::int8_t ind = (halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 : (halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2); + m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration |= (1<begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex); + tf.m_pointsOnPositiveSide->erase(it); + if (tf.m_pointsOnPositiveSide->size()==0) { + reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide); + } + continue; + } + + // Except for the horizon edges, all half edges of the visible faces can be marked as disabled. Their data slots will be reused. + // The faces will be disabled as well, but we need to remember the points that were on the positive side of them - therefore + // we save pointers to them. + m_newFaceIndices.clear(); + m_newHalfEdgeIndices.clear(); + m_disabledFacePointVectors.clear(); + size_t disableCounter = 0; + for (auto faceIndex : m_visibleFaces) { + auto& disabledFace = m_mesh.m_faces[faceIndex]; + auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace); + for (size_t j=0;j<3;j++) { + if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1<size()); // Because we should not assign point vectors to faces unless needed... + m_disabledFacePointVectors.push_back(std::move(t)); + } + } + if (disableCounter < horizonEdgeCount*2) { + const size_t newHalfEdgesNeeded = horizonEdgeCount*2-disableCounter; + for (size_t i=0;i planeNormal = mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint); + newFace.m_P = Plane(planeNormal,activePoint); + newFace.m_he = AB; + + m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ? i*2-1 : 2*horizonEdgeCount-1]; + m_mesh.m_halfEdges[BC].m_opp = m_newHalfEdgeIndices[((i+1)*2) % (horizonEdgeCount*2)]; + } + + // Assign points that were on the positive side of the disabled faces to the new faces. + for (auto& disabledPoints : m_disabledFacePointVectors) { + assert(disabledPoints); + for (const auto& point : *(disabledPoints)) { + if (point == activePointIndex) { + continue; + } + for (size_t j=0;jsize()>0); + if (!newFace.m_inFaceStack) { + m_faceList.push_back(newFaceIndex); + newFace.m_inFaceStack = 1; + } + } + } + } + + // Cleanup + m_indexVectorPool.clear(); + } + + /* + * Private helper functions + */ + + template + std::array QuickHull::getExtremeValues() { + std::array outIndices{0,0,0,0,0,0}; + T extremeVals[6] = {m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z}; + const size_t vCount = m_vertexData.size(); + for (size_t i=1;i& pos = m_vertexData[i]; + if (pos.x>extremeVals[0]) { + extremeVals[0]=pos.x; + outIndices[0]=i; + } + else if (pos.xextremeVals[2]) { + extremeVals[2]=pos.y; + outIndices[2]=i; + } + else if (pos.yextremeVals[4]) { + extremeVals[4]=pos.z; + outIndices[4]=i; + } + else if (pos.z + bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { + const size_t horizonEdgeCount = horizonEdges.size(); + for (size_t i=0;i + T QuickHull::getScale(const std::array& extremeValues) { + T s = 0; + for (size_t i=0;i<6;i++) { + const T* v = (const T*)(&m_vertexData[extremeValues[i]]); + v += i/2; + auto a = std::abs(*v); + if (a>s) { + s = a; + } + } + return s; + } + + template + void QuickHull::setupInitialTetrahedron() { + const size_t vertexCount = m_vertexData.size(); + + // If we have at most 4 points, just return a degenerate tetrahedron: + if (vertexCount <= 4) { + size_t v[4] = {0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)}; + const Vector3 N = mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]); + const Plane trianglePlane(N,m_vertexData[v[0]]); + if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { + std::swap(v[0],v[1]); + } + return m_mesh.setup(v[0],v[1],v[2],v[3]); + } + + // Find two most distant extreme points. + T maxD = m_epsilonSquared; + std::pair selectedPoints; + for (size_t i=0;i<6;i++) { + for (size_t j=i+1;j<6;j++) { + const T d = m_vertexData[ m_extremeValues[i] ].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] ); + if (d > maxD) { + maxD=d; + selectedPoints={m_extremeValues[i],m_extremeValues[j]}; + } + } + } + if (maxD == m_epsilonSquared) { + // A degenerate case: the point cloud seems to consists of a single point + return m_mesh.setup(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)); + } + assert(selectedPoints.first != selectedPoints.second); + + // Find the most distant point to the line between the two chosen extreme points. + const Ray r(m_vertexData[selectedPoints.first], (m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first])); + maxD = m_epsilonSquared; + size_t maxI=std::numeric_limits::max(); + const size_t vCount = m_vertexData.size(); + for (size_t i=0;i maxD) { + maxD=distToRay; + maxI=i; + } + } + if (maxD == m_epsilonSquared) { + // It appears that the point cloud belongs to a 1 dimensional subspace of R^3: convex hull has no volume => return a thin triangle + // Pick any point other than selectedPoints.first and selectedPoints.second as the third point of the triangle + auto it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) { + return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second]; + }); + const size_t thirdPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it); + it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) { + return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint]; + }); + const size_t fourthPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it); + return m_mesh.setup(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint); + } + + // These three points form the base triangle for our tetrahedron. + assert(selectedPoints.first != maxI && selectedPoints.second != maxI); + std::array baseTriangle{selectedPoints.first, selectedPoints.second, maxI}; + const Vector3 baseTriangleVertices[]={ m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]], m_vertexData[baseTriangle[2]] }; + + // Next step is to find the 4th vertex of the tetrahedron. We naturally choose the point farthest away from the triangle plane. + maxD=m_epsilon; + maxI=0; + const Vector3 N = mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]); + Plane trianglePlane(N,baseTriangleVertices[0]); + for (size_t i=0;imaxD) { + maxD=d; + maxI=i; + } + } + if (maxD == m_epsilon) { + // All the points seem to lie on a 2D subspace of R^3. How to handle this? Well, let's add one extra point to the point cloud so that the convex hull will have volume. + m_planar = true; + const vec3 N1 = mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]); + m_planarPointCloudTemp.clear(); + m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end()); + const vec3 extraPoint = N1 + m_vertexData[0]; + m_planarPointCloudTemp.push_back(extraPoint); + maxI = m_planarPointCloudTemp.size()-1; + m_vertexData = VertexDataSource(m_planarPointCloudTemp); + } + + // Enforce CCW orientation (if user prefers clockwise orientation, swap two vertices in each triangle when final mesh is created) + const Plane triPlane(N,baseTriangleVertices[0]); + if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) { + std::swap(baseTriangle[0],baseTriangle[1]); + } + + // Create a tetrahedron half edge mesh and compute planes defined by each triangle + m_mesh.setup(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI); + for (auto& f : m_mesh.m_faces) { + auto v = m_mesh.getVertexIndicesOfFace(f); + const Vector3& va = m_vertexData[v[0]]; + const Vector3& vb = m_vertexData[v[1]]; + const Vector3& vc = m_vertexData[v[2]]; + const Vector3 N1 = mathutils::getTriangleNormal(va, vb, vc); + const Plane plane(N1,va); + f.m_P = plane; + } + + // Finally we assign a face for each vertex outside the tetrahedron (vertices inside the tetrahedron have no role anymore) + for (size_t i=0;i; + template class QuickHull; +} + diff --git a/src/manifold/src/quickhull.hpp b/src/manifold/src/quickhull.hpp new file mode 100644 index 000000000..f3096ce54 --- /dev/null +++ b/src/manifold/src/quickhull.hpp @@ -0,0 +1,981 @@ +#ifndef QUICKHULL_HPP_ +#define QUICKHULL_HPP_ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Pool.hpp + + +namespace quickhull { + + template + class Pool { + std::vector> m_data; + public: + void clear() { + m_data.clear(); + } + + void reclaim(std::unique_ptr& ptr) { + m_data.push_back(std::move(ptr)); + } + + std::unique_ptr get() { + if (m_data.size()==0) { + return std::unique_ptr(new T()); + } + auto it = m_data.end()-1; + std::unique_ptr r = std::move(*it); + m_data.erase(it); + return r; + } + + }; + + +// Vector3.hpp + + + template + class Vector3 + { + public: + Vector3() = default; + + Vector3(T x, T y, T z) : x(x), y(y), z(z) { + + } + + T x,y,z; + + T dotProduct(const Vector3& other) const { + return x*other.x+y*other.y+z*other.z; + } + + void normalize() { + const T len = getLength(); + x/=len; + y/=len; + z/=len; + } + + Vector3 getNormalized() const { + const T len = getLength(); + return Vector3(x/len,y/len,z/len); + } + + T getLength() const { + return std::sqrt(x*x+y*y+z*z); + } + + Vector3 operator-(const Vector3& other) const { + return Vector3(x-other.x,y-other.y,z-other.z); + } + + Vector3 operator+(const Vector3& other) const { + return Vector3(x+other.x,y+other.y,z+other.z); + } + + Vector3& operator+=(const Vector3& other) { + x+=other.x; + y+=other.y; + z+=other.z; + return *this; + } + Vector3& operator-=(const Vector3& other) { + x-=other.x; + y-=other.y; + z-=other.z; + return *this; + } + Vector3& operator*=(T c) { + x*=c; + y*=c; + z*=c; + return *this; + } + + Vector3& operator/=(T c) { + x/=c; + y/=c; + z/=c; + return *this; + } + + Vector3 operator-() const { + return Vector3(-x,-y,-z); + } + + template + Vector3 operator*(S c) const { + return Vector3(x*c,y*c,z*c); + } + + template + Vector3 operator/(S c) const { + return Vector3(x/c,y/c,z/c); + } + + T getLengthSquared() const { + return x*x + y*y + z*z; + } + + bool operator!=(const Vector3& o) const { + return x != o.x || y != o.y || z != o.z; + } + + // Projection onto another vector + Vector3 projection(const Vector3& o) const { + T C = dotProduct(o)/o.getLengthSquared(); + return o*C; + } + + Vector3 crossProduct (const Vector3& rhs ) { + T a = y * rhs.z - z * rhs.y ; + T b = z * rhs.x - x * rhs.z ; + T c = x * rhs.y - y * rhs.x ; + Vector3 product( a , b , c ) ; + return product ; + } + + T getDistanceTo(const Vector3& other) const { + Vector3 diff = *this - other; + return diff.getLength(); + } + + T getSquaredDistanceTo(const Vector3& other) const { + const T dx = x-other.x; + const T dy = y-other.y; + const T dz = z-other.z; + return dx*dx+dy*dy+dz*dz; + } + + }; + + // Overload also << operator for easy printing of debug data + template + inline std::ostream& operator<<(std::ostream& os, const Vector3& vec) { + os << "(" << vec.x << "," << vec.y << "," << vec.z << ")"; + return os; + } + + template + inline Vector3 operator*(T c, const Vector3& v) { + return Vector3(v.x*c,v.y*c,v.z*c); + } + +// Plane.hpp + + template + class Plane { + public: + Vector3 m_N; + + // Signed distance (if normal is of length 1) to the plane from origin + T m_D; + + // Normal length squared + T m_sqrNLength; + + bool isPointOnPositiveSide(const Vector3& Q) const { + T d = m_N.dotProduct(Q)+m_D; + if (d>=0) return true; + return false; + } + + Plane() = default; + + // Construct a plane using normal N and any point P on the plane + Plane(const Vector3& N, const Vector3& P) : m_N(N), m_D(-N.dotProduct(P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { + + } + }; + +// Ray.hpp + + template + struct Ray { + const Vector3 m_S; + const Vector3 m_V; + const T m_VInvLengthSquared; + + Ray(const Vector3& S,const Vector3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/m_V.getLengthSquared()) { + } + }; + + +// VertexDataSource + + + template + class VertexDataSource { + const Vector3* m_ptr; + size_t m_count; + + public: + VertexDataSource(const Vector3* ptr, size_t count) : m_ptr(ptr), m_count(count) { + + } + + VertexDataSource(const std::vector>& vec) : m_ptr(&vec[0]), m_count(vec.size()) { + + } + + VertexDataSource() : m_ptr(nullptr), m_count(0) { + + } + + VertexDataSource& operator=(const VertexDataSource& other) = default; + + size_t size() const { + return m_count; + } + + const Vector3& operator[](size_t index) const { + return m_ptr[index]; + } + + const Vector3* begin() const { + return m_ptr; + } + + const Vector3* end() const { + return m_ptr + m_count; + } + }; + + + +// Mesh.hpp + + + template + class MeshBuilder { + public: + struct HalfEdge { + size_t m_endVertex; + size_t m_opp; + size_t m_face; + size_t m_next; + + void disable() { + m_endVertex = std::numeric_limits::max(); + } + + bool isDisabled() const { + return m_endVertex == std::numeric_limits::max(); + } + }; + + struct Face { + size_t m_he; + Plane m_P{}; + T m_mostDistantPointDist; + size_t m_mostDistantPoint; + size_t m_visibilityCheckedOnIteration; + std::uint8_t m_isVisibleFaceOnCurrentIteration : 1; + std::uint8_t m_inFaceStack : 1; + std::uint8_t m_horizonEdgesOnCurrentIteration : 3; // Bit for each half edge assigned to this face, each being 0 or 1 depending on whether the edge belongs to horizon edge + std::unique_ptr> m_pointsOnPositiveSide; + + Face() : m_he(std::numeric_limits::max()), + m_mostDistantPointDist(0), + m_mostDistantPoint(0), + m_visibilityCheckedOnIteration(0), + m_isVisibleFaceOnCurrentIteration(0), + m_inFaceStack(0), + m_horizonEdgesOnCurrentIteration(0) + { + + } + + void disable() { + m_he = std::numeric_limits::max(); + } + + bool isDisabled() const { + return m_he == std::numeric_limits::max(); + } + }; + + // Mesh data + std::vector m_faces; + std::vector m_halfEdges; + + // When the mesh is modified and faces and half edges are removed from it, we do not actually remove them from the container vectors. + // Insted, they are marked as disabled which means that the indices can be reused when we need to add new faces and half edges to the mesh. + // We store the free indices in the following vectors. + std::vector m_disabledFaces,m_disabledHalfEdges; + + size_t addFace() { + if (m_disabledFaces.size()) { + size_t index = m_disabledFaces.back(); + auto& f = m_faces[index]; + assert(f.isDisabled()); + assert(!f.m_pointsOnPositiveSide); + f.m_mostDistantPointDist = 0; + m_disabledFaces.pop_back(); + return index; + } + m_faces.emplace_back(); + return m_faces.size()-1; + } + + size_t addHalfEdge() { + if (m_disabledHalfEdges.size()) { + const size_t index = m_disabledHalfEdges.back(); + m_disabledHalfEdges.pop_back(); + return index; + } + m_halfEdges.emplace_back(); + return m_halfEdges.size()-1; + } + + // Mark a face as disabled and return a pointer to the points that were on the positive of it. + std::unique_ptr> disableFace(size_t faceIndex) { + auto& f = m_faces[faceIndex]; + f.disable(); + m_disabledFaces.push_back(faceIndex); + return std::move(f.m_pointsOnPositiveSide); + } + + void disableHalfEdge(size_t heIndex) { + auto& he = m_halfEdges[heIndex]; + he.disable(); + m_disabledHalfEdges.push_back(heIndex); + } + + MeshBuilder() = default; + + // Create a mesh with initial tetrahedron ABCD. Dot product of AB with the normal of triangle ABC should be negative. + void setup(size_t a, size_t b, size_t c, size_t d) { + m_faces.clear(); + m_halfEdges.clear(); + m_disabledFaces.clear(); + m_disabledHalfEdges.clear(); + + m_faces.reserve(4); + m_halfEdges.reserve(12); + + // Create halfedges + HalfEdge AB; + AB.m_endVertex = b; + AB.m_opp = 6; + AB.m_face = 0; + AB.m_next = 1; + m_halfEdges.push_back(AB); + + HalfEdge BC; + BC.m_endVertex = c; + BC.m_opp = 9; + BC.m_face = 0; + BC.m_next = 2; + m_halfEdges.push_back(BC); + + HalfEdge CA; + CA.m_endVertex = a; + CA.m_opp = 3; + CA.m_face = 0; + CA.m_next = 0; + m_halfEdges.push_back(CA); + + HalfEdge AC; + AC.m_endVertex = c; + AC.m_opp = 2; + AC.m_face = 1; + AC.m_next = 4; + m_halfEdges.push_back(AC); + + HalfEdge CD; + CD.m_endVertex = d; + CD.m_opp = 11; + CD.m_face = 1; + CD.m_next = 5; + m_halfEdges.push_back(CD); + + HalfEdge DA; + DA.m_endVertex = a; + DA.m_opp = 7; + DA.m_face = 1; + DA.m_next = 3; + m_halfEdges.push_back(DA); + + HalfEdge BA; + BA.m_endVertex = a; + BA.m_opp = 0; + BA.m_face = 2; + BA.m_next = 7; + m_halfEdges.push_back(BA); + + HalfEdge AD; + AD.m_endVertex = d; + AD.m_opp = 5; + AD.m_face = 2; + AD.m_next = 8; + m_halfEdges.push_back(AD); + + HalfEdge DB; + DB.m_endVertex = b; + DB.m_opp = 10; + DB.m_face = 2; + DB.m_next = 6; + m_halfEdges.push_back(DB); + + HalfEdge CB; + CB.m_endVertex = b; + CB.m_opp = 1; + CB.m_face = 3; + CB.m_next = 10; + m_halfEdges.push_back(CB); + + HalfEdge BD; + BD.m_endVertex = d; + BD.m_opp = 8; + BD.m_face = 3; + BD.m_next = 11; + m_halfEdges.push_back(BD); + + HalfEdge DC; + DC.m_endVertex = c; + DC.m_opp = 4; + DC.m_face = 3; + DC.m_next = 9; + m_halfEdges.push_back(DC); + + // Create faces + Face ABC; + ABC.m_he = 0; + m_faces.push_back(std::move(ABC)); + + Face ACD; + ACD.m_he = 3; + m_faces.push_back(std::move(ACD)); + + Face BAD; + BAD.m_he = 6; + m_faces.push_back(std::move(BAD)); + + Face CBD; + CBD.m_he = 9; + m_faces.push_back(std::move(CBD)); + } + + std::array getVertexIndicesOfFace(const Face& f) const { + std::array v; + const HalfEdge* he = &m_halfEdges[f.m_he]; + v[0] = he->m_endVertex; + he = &m_halfEdges[he->m_next]; + v[1] = he->m_endVertex; + he = &m_halfEdges[he->m_next]; + v[2] = he->m_endVertex; + return v; + } + + std::array getVertexIndicesOfHalfEdge(const HalfEdge& he) const { + return {m_halfEdges[he.m_opp].m_endVertex,he.m_endVertex}; + } + + std::array getHalfEdgeIndicesOfFace(const Face& f) const { + return {f.m_he,m_halfEdges[f.m_he].m_next,m_halfEdges[m_halfEdges[f.m_he].m_next].m_next}; + } + }; + + + + + +// ConvexHull.hpp + + template + class ConvexHull { + std::unique_ptr>> m_optimizedVertexBuffer; + VertexDataSource m_vertices; + std::vector m_indices; + public: + ConvexHull() {} + + // Copy constructor + ConvexHull(const ConvexHull& o) { + m_indices = o.m_indices; + if (o.m_optimizedVertexBuffer) { + m_optimizedVertexBuffer.reset(new std::vector>(*o.m_optimizedVertexBuffer)); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + } + else { + m_vertices = o.m_vertices; + } + } + + ConvexHull& operator=(const ConvexHull& o) { + if (&o == this) { + return *this; + } + m_indices = o.m_indices; + if (o.m_optimizedVertexBuffer) { + m_optimizedVertexBuffer.reset(new std::vector>(*o.m_optimizedVertexBuffer)); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + } + else { + m_vertices = o.m_vertices; + } + return *this; + } + + ConvexHull(ConvexHull&& o) { + m_indices = std::move(o.m_indices); + if (o.m_optimizedVertexBuffer) { + m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); + o.m_vertices = VertexDataSource(); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + } + else { + m_vertices = o.m_vertices; + } + } + + ConvexHull& operator=(ConvexHull&& o) { + if (&o == this) { + return *this; + } + m_indices = std::move(o.m_indices); + if (o.m_optimizedVertexBuffer) { + m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); + o.m_vertices = VertexDataSource(); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + } + else { + m_vertices = o.m_vertices; + } + return *this; + } + + // Construct vertex and index buffers from half edge mesh and pointcloud + ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { + if (!useOriginalIndices) { + m_optimizedVertexBuffer.reset(new std::vector>()); + } + + std::vector faceProcessed(mesh.m_faces.size(),false); + std::vector faceStack; + std::unordered_map vertexIndexMapping; // Map vertex indices from original point cloud to the new mesh vertex indices + for (size_t i = 0;ipush_back(pointCloud[v]); + vertexIndexMapping[v] = m_optimizedVertexBuffer->size()-1; + v = m_optimizedVertexBuffer->size()-1; + } + else { + v = itV->second; + } + } + } + m_indices.push_back(vertices[0]); + m_indices.push_back(vertices[1 + iCCW]); + m_indices.push_back(vertices[2 - iCCW]); + } + } + + if (!useOriginalIndices) { + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + } + else { + m_vertices = pointCloud; + } + } + + std::vector& getIndexBuffer() { + return m_indices; + } + + const std::vector& getIndexBuffer() const { + return m_indices; + } + + VertexDataSource& getVertexBuffer() { + return m_vertices; + } + + const VertexDataSource& getVertexBuffer() const { + return m_vertices; + } + + // Export the mesh to a Waveform OBJ file + void writeWaveformOBJ(const std::string& filename, const std::string& objectName = "quickhull") const + { + std::ofstream objFile; + objFile.open (filename); + objFile << "o " << objectName << "\n"; + for (const auto& v : getVertexBuffer()) { + objFile << "v " << v.x << " " << v.y << " " << v.z << "\n"; + } + const auto& indBuf = getIndexBuffer(); + size_t triangleCount = indBuf.size()/3; + for (size_t i=0;i + class HalfEdgeMesh { + public: + + struct HalfEdge { + IndexType m_endVertex; + IndexType m_opp; + IndexType m_face; + IndexType m_next; + }; + + struct Face { + IndexType m_halfEdgeIndex; // Index of one of the half edges of this face + }; + + std::vector> m_vertices; + std::vector m_faces; + std::vector m_halfEdges; + + HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) + { + std::unordered_map faceMapping; + std::unordered_map halfEdgeMapping; + std::unordered_map vertexMapping; + + size_t i=0; + for (const auto& face : builderObject.m_faces) { + if (!face.isDisabled()) { + m_faces.push_back({static_cast(face.m_he)}); + faceMapping[i] = m_faces.size()-1; + + const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face); + for (const auto heIndex : heIndices) { + const IndexType vertexIndex = builderObject.m_halfEdges[heIndex].m_endVertex; + if (vertexMapping.count(vertexIndex)==0) { + m_vertices.push_back(vertexData[vertexIndex]); + vertexMapping[vertexIndex] = m_vertices.size()-1; + } + } + } + i++; + } + + i=0; + for (const auto& halfEdge : builderObject.m_halfEdges) { + if (!halfEdge.isDisabled()) { + m_halfEdges.push_back({static_cast(halfEdge.m_endVertex),static_cast(halfEdge.m_opp),static_cast(halfEdge.m_face),static_cast(halfEdge.m_next)}); + halfEdgeMapping[i] = m_halfEdges.size()-1; + } + i++; + } + + for (auto& face : m_faces) { + assert(halfEdgeMapping.count(face.m_halfEdgeIndex) == 1); + face.m_halfEdgeIndex = halfEdgeMapping[face.m_halfEdgeIndex]; + } + + for (auto& he : m_halfEdges) { + he.m_face = faceMapping[he.m_face]; + he.m_opp = halfEdgeMapping[he.m_opp]; + he.m_next = halfEdgeMapping[he.m_next]; + he.m_endVertex = vertexMapping[he.m_endVertex]; + } + } + + }; + + + +// MathUtils.hpp + + namespace mathutils { + + template + inline T getSquaredDistanceBetweenPointAndRay(const Vector3& p, const Ray& r) { + const Vector3 s = p-r.m_S; + T t = s.dotProduct(r.m_V); + return s.getLengthSquared() - t*t*r.m_VInvLengthSquared; + } + + // Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance). + template + inline T getSignedDistanceToPlane(const Vector3& v, const Plane& p) { + return p.m_N.dotProduct(v) + p.m_D; + } + + template + inline Vector3 getTriangleNormal(const Vector3& a,const Vector3& b,const Vector3& c) { + // We want to get (a-c).crossProduct(b-c) without constructing temp vectors + T x = a.x - c.x; + T y = a.y - c.y; + T z = a.z - c.z; + T rhsx = b.x - c.x; + T rhsy = b.y - c.y; + T rhsz = b.z - c.z; + T px = y * rhsz - z * rhsy ; + T py = z * rhsx - x * rhsz ; + T pz = x * rhsy - y * rhsx ; + return Vector3(px,py,pz); + } + + + } + + +// QuickHull.hpp + +/* + * Implementation of the 3D QuickHull algorithm by Antti Kuukka + * + * No copyrights. What follows is 100% Public Domain. + * + * + * + * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) + * + * OUTPUT: a ConvexHull object which provides vertex and index buffers of the generated convex hull as a triangle mesh. + * + * + * + * The implementation is thread-safe if each thread is using its own QuickHull object. + * + * + * SUMMARY OF THE ALGORITHM: + * - Create initial simplex (tetrahedron) using extreme points. We have four faces now and they form a convex mesh M. + * - For each point, assign them to the first face for which they are on the positive side of (so each point is assigned to at most + * one face). Points inside the initial tetrahedron are left behind now and no longer affect the calculations. + * - Add all faces that have points assigned to them to Face Stack. + * - Iterate until Face Stack is empty: + * - Pop topmost face F from the stack + * - From the points assigned to F, pick the point P that is farthest away from the plane defined by F. + * - Find all faces of M that have P on their positive side. Let us call these the "visible faces". + * - Because of the way M is constructed, these faces are connected. Solve their horizon edge loop. + * - "Extrude to P": Create new faces by connecting P with the points belonging to the horizon edge. Add the new faces to M and remove the visible + * faces from M. + * - Each point that was assigned to visible faces is now assigned to at most one of the newly created faces. + * - Those new faces that have points assigned to them are added to the top of Face Stack. + * - M is now the convex hull. + * + * TO DO: + * - Implement a proper 2D QuickHull and use that to solve the degenerate 2D case (when all the points lie on the same plane in 3D space). + * */ + + struct DiagnosticsData { + size_t m_failedHorizonEdges; // How many times QuickHull failed to solve the horizon edge. Failures lead to degenerated convex hulls. + + DiagnosticsData() : m_failedHorizonEdges(0) { } + }; + + template + FloatType defaultEps(); + + template + class QuickHull { + using vec3 = Vector3; + + FloatType m_epsilon, m_epsilonSquared, m_scale; + bool m_planar; + std::vector m_planarPointCloudTemp; + VertexDataSource m_vertexData; + MeshBuilder m_mesh; + std::array m_extremeValues; + DiagnosticsData m_diagnostics; + + // Temporary variables used during iteration process + std::vector m_newFaceIndices; + std::vector m_newHalfEdgeIndices; + std::vector< std::unique_ptr> > m_disabledFacePointVectors; + std::vector m_visibleFaces; + std::vector m_horizonEdges; + struct FaceData { + size_t m_faceIndex; + size_t m_enteredFromHalfEdge; // If the face turns out not to be visible, this half edge will be marked as horizon edge + FaceData(size_t fi, size_t he) : m_faceIndex(fi),m_enteredFromHalfEdge(he) {} + }; + std::vector m_possiblyVisibleFaces; + std::deque m_faceList; + + // Create a half edge mesh representing the base tetrahedron from which the QuickHull iteration proceeds. m_extremeValues must be properly set up when this is called. + void setupInitialTetrahedron(); + + // Given a list of half edges, try to rearrange them so that they form a loop. Return true on success. + bool reorderHorizonEdges(std::vector& horizonEdges); + + // Find indices of extreme values (max x, min x, max y, min y, max z, min z) for the given point cloud + std::array getExtremeValues(); + + // Compute scale of the vertex data. + FloatType getScale(const std::array& extremeValues); + + // Each face contains a unique pointer to a vector of indices. However, many - often most - faces do not have any points on the positive + // side of them especially at the the end of the iteration. When a face is removed from the mesh, its associated point vector, if such + // exists, is moved to the index vector pool, and when we need to add new faces with points on the positive side to the mesh, + // we reuse these vectors. This reduces the amount of std::vectors we have to deal with, and impact on performance is remarkable. + Pool> m_indexVectorPool; + inline std::unique_ptr> getIndexVectorFromPool(); + inline void reclaimToIndexVectorPool(std::unique_ptr>& ptr); + + // Associates a point with a face if the point resides on the positive side of 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 m_mesh from which we create the ConvexHull object that getConvexHull function returns + void createConvexHalfEdgeMesh(); + + // Constructs the convex hull into a MeshBuilder object which can be converted to a ConvexHull or Mesh object + void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + + // The public getConvexHull functions will setup a VertexDataSource object and call this + ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + public: + // Computes convex hull for a given point cloud. + // Params: + // pointCloud: a vector of of 3D points + // CCW: whether the output mesh triangles should have CCW orientation + // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive of it (for a point cloud with scale 1) + ConvexHull getConvexHull(const std::vector>& pointCloud, + bool CCW, + bool useOriginalIndices, + FloatType eps = defaultEps()); + + // Computes convex hull for a given point cloud. + // Params: + // vertexData: pointer to the first 3D point of the point cloud + // vertexCount: number of vertices in the point cloud + // CCW: whether the output mesh triangles should have CCW orientation + // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) + ConvexHull getConvexHull(const Vector3* vertexData, + size_t vertexCount, + bool CCW, + bool useOriginalIndices, + FloatType eps = defaultEps()); + + // 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 + // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) + ConvexHull getConvexHull(const FloatType* vertexData, + size_t vertexCount, + bool CCW, + bool useOriginalIndices, + FloatType eps = defaultEps()); + + // 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 a mesh object with half edge structure. + HalfEdgeMesh getConvexHullAsMesh(const FloatType* vertexData, + size_t vertexCount, + bool CCW, + FloatType eps = defaultEps()); + + // Get diagnostics about last generated convex hull + const DiagnosticsData& getDiagnostics() { + return m_diagnostics; + } + }; + + /* + * Inline function definitions + */ + + template + std::unique_ptr> QuickHull::getIndexVectorFromPool() { + auto r = m_indexVectorPool.get(); + r->clear(); + return r; + } + + template + void QuickHull::reclaimToIndexVectorPool(std::unique_ptr>& ptr) { + const size_t oldSize = ptr->size(); + if ((oldSize+1)*128 < ptr->capacity()) { + // Reduce memory usage! Huge vectors are needed at the beginning of iteration when faces have many points on their positive side. Later on, smaller vectors will suffice. + ptr.reset(nullptr); + return; + } + m_indexVectorPool.reclaim(ptr); + } + + template + bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { + const T D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); + if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) { + if (!f.m_pointsOnPositiveSide) { + f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool()); + } + f.m_pointsOnPositiveSide->push_back( pointIndex ); + if (D > f.m_mostDistantPointDist) { + f.m_mostDistantPointDist = D; + f.m_mostDistantPoint = pointIndex; + } + return true; + } + return false; + } + +} + + +#endif /* QUICKHULL_HPP_ */ diff --git a/src/third_party/CMakeLists.txt b/src/third_party/CMakeLists.txt index 7d912ac4f..8b1378917 100644 --- a/src/third_party/CMakeLists.txt +++ b/src/third_party/CMakeLists.txt @@ -1,8 +1 @@ -add_library(quickhull OBJECT quickhull/QuickHull.cpp) -target_include_directories(quickhull PUBLIC - $ - $) - -target_compile_features(quickhull PUBLIC cxx_std_17) -install(TARGETS quickhull EXPORT manifoldTargets) diff --git a/src/third_party/quickhull b/src/third_party/quickhull deleted file mode 160000 index 1ffbc6f88..000000000 --- a/src/third_party/quickhull +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 1ffbc6f884ea1da89e104a5996cf8a726db673d5 From 24b51a05d583c123b53f97d8077ffb9ce7bc8ceb Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Wed, 31 Jul 2024 16:28:46 +0530 Subject: [PATCH 02/24] Replaced Vector3 with glm::vec3 and VertexDataSource (and formatting) --- bindings/python/manifold3d.cpp | 9 +- extras/perf_test_cgal.cpp | 2 +- src/manifold/src/manifold.cpp | 6 +- src/manifold/src/quickhull.cpp | 1066 +++++++++++++++++--------------- src/manifold/src/quickhull.hpp | 247 ++------ 5 files changed, 646 insertions(+), 684 deletions(-) diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index f8cb5de04..4dee244f4 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -257,10 +257,7 @@ NB_MODULE(manifold3d, m) { manifold__translate__v) .def("scale", &Manifold::Scale, nb::arg("v"), manifold__scale__v) .def( - "scale", - [](const Manifold &m, float s) { - m.Scale({s, s, s}); - }, + "scale", [](const Manifold &m, float s) { m.Scale({s, s, s}); }, nb::arg("s"), "Scale this Manifold in space. This operation can be chained. " "Transforms are combined and applied lazily.\n\n" @@ -673,9 +670,7 @@ NB_MODULE(manifold3d, m) { cross_section__scale__scale) .def( "scale", - [](const CrossSection &self, float s) { - self.Scale({s, s}); - }, + [](const CrossSection &self, float s) { self.Scale({s, s}); }, nb::arg("s"), "Scale this CrossSection in space. This operation can be chained. " "Transforms are combined and applied lazily." diff --git a/extras/perf_test_cgal.cpp b/extras/perf_test_cgal.cpp index 562511425..6b48b2849 100644 --- a/extras/perf_test_cgal.cpp +++ b/extras/perf_test_cgal.cpp @@ -40,7 +40,7 @@ typedef CGAL::SM_Vertex_index Vertex; void manifoldToCGALSurfaceMesh(Manifold &manifold, TriangleMesh &cgalMesh) { auto maniMesh = manifold.GetMesh(); - const int n = maniMesh.vertPos.size(); + const size_t n = maniMesh.vertPos.size(); std::vector vertices(n); for (size_t i = 0; i < n; i++) { auto &vert = maniMesh.vertPos[i]; diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 34af207b9..ae97ad835 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -16,11 +16,11 @@ #include #include -#include "quickhull.hpp" #include "boolean3.h" #include "csg_tree.h" #include "impl.h" #include "par.h" +#include "quickhull.hpp" namespace { using namespace manifold; @@ -919,12 +919,12 @@ Manifold Manifold::Hull(const std::vector& pts) { const int numVert = pts.size(); if (numVert < 4) return Manifold(); - std::vector> vertices(numVert); + std::vector vertices(numVert); for (int i = 0; i < numVert; i++) { vertices[i] = {pts[i].x, pts[i].y, pts[i].z}; } - quickhull::QuickHull qh; + quickhull::QuickHull qh; // bools: correct triangle winding, and use original indices auto hull = qh.getConvexHull(vertices, false, true); const auto& triangles = hull.getIndexBuffer(); diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index e291ea5e2..4a6120caf 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -1,503 +1,585 @@ #include "quickhull.hpp" // #include "MathUtils.hpp" -#include +#include #include +#include #include -#include #include // #include "Structs/Mesh.hpp" namespace quickhull { - template<> - float defaultEps() { - return 0.0001f; - } - - template<> - double defaultEps() { - return 0.0000001; - } - - /* - * Implementation of the algorithm - */ - - template - ConvexHull QuickHull::getConvexHull(const std::vector>& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { - VertexDataSource vertexDataSource(pointCloud); - return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); - } - - template - ConvexHull QuickHull::getConvexHull(const Vector3* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) { - VertexDataSource vertexDataSource(vertexData,vertexCount); - return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); - } - - template - ConvexHull QuickHull::getConvexHull(const T* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) { - VertexDataSource vertexDataSource((const vec3*)vertexData,vertexCount); - return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon); - } - - template - HalfEdgeMesh QuickHull::getConvexHullAsMesh(const FloatType* vertexData, size_t vertexCount, bool CCW, FloatType epsilon) { - VertexDataSource vertexDataSource((const vec3*)vertexData,vertexCount); - buildMesh(vertexDataSource, CCW, false, epsilon); - return HalfEdgeMesh(m_mesh, m_vertexData); - } - - template - void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { - // CCW is unused for now - (void)CCW; - // useOriginalIndices is unused for now - (void)useOriginalIndices; - - if (pointCloud.size()==0) { - m_mesh = MeshBuilder(); - return; - } - m_vertexData = pointCloud; - - // Very first: find extreme values and use them to compute the scale of the point cloud. - m_extremeValues = getExtremeValues(); - m_scale = getScale(m_extremeValues); - - // Epsilon we use depends on the scale - m_epsilon = epsilon*m_scale; - m_epsilonSquared = m_epsilon*m_epsilon; - - // Reset diagnostics - m_diagnostics = DiagnosticsData(); - - m_planar = false; // The planar case happens when all the points appear to lie on a two dimensional subspace of R^3. - createConvexHalfEdgeMesh(); - if (m_planar) { - const size_t extraPointIndex = m_planarPointCloudTemp.size()-1; - for (auto& he : m_mesh.m_halfEdges) { - if (he.m_endVertex == extraPointIndex) { - he.m_endVertex = 0; - } - } - m_vertexData = pointCloud; - m_planarPointCloudTemp.clear(); - } - } - - template - ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) { - buildMesh(pointCloud,CCW,useOriginalIndices,epsilon); - return ConvexHull(m_mesh,m_vertexData, CCW, useOriginalIndices); - } - - template - void QuickHull::createConvexHalfEdgeMesh() { - m_visibleFaces.clear(); - m_horizonEdges.clear(); - m_possiblyVisibleFaces.clear(); - - // Compute base tetrahedron - setupInitialTetrahedron(); - assert(m_mesh.m_faces.size()==4); - - // Init face stack with those faces that have points assigned to them - m_faceList.clear(); - for (size_t i=0;i < 4;i++) { - auto& f = m_mesh.m_faces[i]; - if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size()>0) { - m_faceList.push_back(i); - f.m_inFaceStack = 1; - } - } - - // Process faces until the face list is empty. - size_t iter = 0; - while (!m_faceList.empty()) { - iter++; - if (iter == std::numeric_limits::max()) { - // Visible face traversal marks visited faces with iteration counter (to mark that the face has been visited on this iteration) and the max value represents unvisited faces. At this point we have to reset iteration counter. This shouldn't be an - // issue on 64 bit machines. - iter = 0; - } - - const size_t topFaceIndex = m_faceList.front(); - m_faceList.pop_front(); - - auto& tf = m_mesh.m_faces[topFaceIndex]; - tf.m_inFaceStack = 0; - - assert(!tf.m_pointsOnPositiveSide || tf.m_pointsOnPositiveSide->size()>0); - if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) { - continue; - } - - // Pick the most distant point to this triangle plane as the point to which we extrude - const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint]; - const size_t activePointIndex = tf.m_mostDistantPoint; - - // Find out the faces that have our active point on their positive side (these are the "visible faces"). The face on top of the stack of course is one of them. At the same time, we create a list of horizon edges. - m_horizonEdges.clear(); - m_possiblyVisibleFaces.clear(); - m_visibleFaces.clear(); - m_possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits::max()); - while (m_possiblyVisibleFaces.size()) { - const auto faceData = m_possiblyVisibleFaces.back(); - m_possiblyVisibleFaces.pop_back(); - auto& pvf = m_mesh.m_faces[faceData.m_faceIndex]; - assert(!pvf.isDisabled()); - - if (pvf.m_visibilityCheckedOnIteration == iter) { - if (pvf.m_isVisibleFaceOnCurrentIteration) { - continue; - } - } - else { - const Plane& P = pvf.m_P; - pvf.m_visibilityCheckedOnIteration = iter; - const T d = P.m_N.dotProduct(activePoint)+P.m_D; - if (d>0) { - pvf.m_isVisibleFaceOnCurrentIteration = 1; - pvf.m_horizonEdgesOnCurrentIteration = 0; - m_visibleFaces.push_back(faceData.m_faceIndex); - for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) { - if (m_mesh.m_halfEdges[heIndex].m_opp != faceData.m_enteredFromHalfEdge) { - m_possiblyVisibleFaces.emplace_back( m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex ); - } - } - continue; - } - assert(faceData.m_faceIndex != topFaceIndex); - } - - // The face is not visible. Therefore, the halfedge we came from is part of the horizon edge. - pvf.m_isVisibleFaceOnCurrentIteration = 0; - m_horizonEdges.push_back(faceData.m_enteredFromHalfEdge); - // Store which half edge is the horizon edge. The other half edges of the face will not be part of the final mesh so their data slots can by recycled. - const auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]); - const std::int8_t ind = (halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 : (halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2); - m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration |= (1<begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex); - tf.m_pointsOnPositiveSide->erase(it); - if (tf.m_pointsOnPositiveSide->size()==0) { - reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide); - } - continue; - } - - // Except for the horizon edges, all half edges of the visible faces can be marked as disabled. Their data slots will be reused. - // The faces will be disabled as well, but we need to remember the points that were on the positive side of them - therefore - // we save pointers to them. - m_newFaceIndices.clear(); - m_newHalfEdgeIndices.clear(); - m_disabledFacePointVectors.clear(); - size_t disableCounter = 0; - for (auto faceIndex : m_visibleFaces) { - auto& disabledFace = m_mesh.m_faces[faceIndex]; - auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace); - for (size_t j=0;j<3;j++) { - if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1<size()); // Because we should not assign point vectors to faces unless needed... - m_disabledFacePointVectors.push_back(std::move(t)); - } - } - if (disableCounter < horizonEdgeCount*2) { - const size_t newHalfEdgesNeeded = horizonEdgeCount*2-disableCounter; - for (size_t i=0;i planeNormal = mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint); - newFace.m_P = Plane(planeNormal,activePoint); - newFace.m_he = AB; - - m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ? i*2-1 : 2*horizonEdgeCount-1]; - m_mesh.m_halfEdges[BC].m_opp = m_newHalfEdgeIndices[((i+1)*2) % (horizonEdgeCount*2)]; - } - - // Assign points that were on the positive side of the disabled faces to the new faces. - for (auto& disabledPoints : m_disabledFacePointVectors) { - assert(disabledPoints); - for (const auto& point : *(disabledPoints)) { - if (point == activePointIndex) { - continue; - } - for (size_t j=0;jsize()>0); - if (!newFace.m_inFaceStack) { - m_faceList.push_back(newFaceIndex); - newFace.m_inFaceStack = 1; - } - } - } - } - - // Cleanup - m_indexVectorPool.clear(); - } - - /* - * Private helper functions - */ - - template - std::array QuickHull::getExtremeValues() { - std::array outIndices{0,0,0,0,0,0}; - T extremeVals[6] = {m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z}; - const size_t vCount = m_vertexData.size(); - for (size_t i=1;i& pos = m_vertexData[i]; - if (pos.x>extremeVals[0]) { - extremeVals[0]=pos.x; - outIndices[0]=i; - } - else if (pos.xextremeVals[2]) { - extremeVals[2]=pos.y; - outIndices[2]=i; - } - else if (pos.yextremeVals[4]) { - extremeVals[4]=pos.z; - outIndices[4]=i; - } - else if (pos.z - bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { - const size_t horizonEdgeCount = horizonEdges.size(); - for (size_t i=0;i - T QuickHull::getScale(const std::array& extremeValues) { - T s = 0; - for (size_t i=0;i<6;i++) { - const T* v = (const T*)(&m_vertexData[extremeValues[i]]); - v += i/2; - auto a = std::abs(*v); - if (a>s) { - s = a; - } - } - return s; - } - - template - void QuickHull::setupInitialTetrahedron() { - const size_t vertexCount = m_vertexData.size(); - - // If we have at most 4 points, just return a degenerate tetrahedron: - if (vertexCount <= 4) { - size_t v[4] = {0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)}; - const Vector3 N = mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]); - const Plane trianglePlane(N,m_vertexData[v[0]]); - if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { - std::swap(v[0],v[1]); - } - return m_mesh.setup(v[0],v[1],v[2],v[3]); - } - - // Find two most distant extreme points. - T maxD = m_epsilonSquared; - std::pair selectedPoints; - for (size_t i=0;i<6;i++) { - for (size_t j=i+1;j<6;j++) { - const T d = m_vertexData[ m_extremeValues[i] ].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] ); - if (d > maxD) { - maxD=d; - selectedPoints={m_extremeValues[i],m_extremeValues[j]}; - } - } - } - if (maxD == m_epsilonSquared) { - // A degenerate case: the point cloud seems to consists of a single point - return m_mesh.setup(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)); - } - assert(selectedPoints.first != selectedPoints.second); - - // Find the most distant point to the line between the two chosen extreme points. - const Ray r(m_vertexData[selectedPoints.first], (m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first])); - maxD = m_epsilonSquared; - size_t maxI=std::numeric_limits::max(); - const size_t vCount = m_vertexData.size(); - for (size_t i=0;i maxD) { - maxD=distToRay; - maxI=i; - } - } - if (maxD == m_epsilonSquared) { - // It appears that the point cloud belongs to a 1 dimensional subspace of R^3: convex hull has no volume => return a thin triangle - // Pick any point other than selectedPoints.first and selectedPoints.second as the third point of the triangle - auto it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) { - return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second]; - }); - const size_t thirdPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it); - it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) { - return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint]; - }); - const size_t fourthPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it); - return m_mesh.setup(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint); - } - - // These three points form the base triangle for our tetrahedron. - assert(selectedPoints.first != maxI && selectedPoints.second != maxI); - std::array baseTriangle{selectedPoints.first, selectedPoints.second, maxI}; - const Vector3 baseTriangleVertices[]={ m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]], m_vertexData[baseTriangle[2]] }; - - // Next step is to find the 4th vertex of the tetrahedron. We naturally choose the point farthest away from the triangle plane. - maxD=m_epsilon; - maxI=0; - const Vector3 N = mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]); - Plane trianglePlane(N,baseTriangleVertices[0]); - for (size_t i=0;imaxD) { - maxD=d; - maxI=i; - } - } - if (maxD == m_epsilon) { - // All the points seem to lie on a 2D subspace of R^3. How to handle this? Well, let's add one extra point to the point cloud so that the convex hull will have volume. - m_planar = true; - const vec3 N1 = mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]); - m_planarPointCloudTemp.clear(); - m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end()); - const vec3 extraPoint = N1 + m_vertexData[0]; - m_planarPointCloudTemp.push_back(extraPoint); - maxI = m_planarPointCloudTemp.size()-1; - m_vertexData = VertexDataSource(m_planarPointCloudTemp); - } - - // Enforce CCW orientation (if user prefers clockwise orientation, swap two vertices in each triangle when final mesh is created) - const Plane triPlane(N,baseTriangleVertices[0]); - if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) { - std::swap(baseTriangle[0],baseTriangle[1]); - } - - // Create a tetrahedron half edge mesh and compute planes defined by each triangle - m_mesh.setup(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI); - for (auto& f : m_mesh.m_faces) { - auto v = m_mesh.getVertexIndicesOfFace(f); - const Vector3& va = m_vertexData[v[0]]; - const Vector3& vb = m_vertexData[v[1]]; - const Vector3& vc = m_vertexData[v[2]]; - const Vector3 N1 = mathutils::getTriangleNormal(va, vb, vc); - const Plane plane(N1,va); - f.m_P = plane; - } - - // Finally we assign a face for each vertex outside the tetrahedron (vertices inside the tetrahedron have no role anymore) - for (size_t i=0;i; - template class QuickHull; +template <> +float defaultEps() { + return 0.0001f; +} + +template <> +double defaultEps() { + return 0.0000001; +} + +/* + * Implementation of the algorithm + */ + +template +ConvexHull QuickHull::getConvexHull( + const std::vector& pointCloud, bool CCW, bool useOriginalIndices, + T epsilon) { + VertexDataSource vertexDataSource(pointCloud); + return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); +} + +template +ConvexHull QuickHull::getConvexHull(const glm::vec3* vertexData, + size_t vertexCount, bool CCW, + bool useOriginalIndices, T epsilon) { + VertexDataSource vertexDataSource(vertexData, vertexCount); + return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); +} + +template +ConvexHull QuickHull::getConvexHull(const T* vertexData, + size_t vertexCount, bool CCW, + bool useOriginalIndices, T epsilon) { + VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); + return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); +} + +template +HalfEdgeMesh QuickHull::getConvexHullAsMesh( + const FloatType* vertexData, size_t vertexCount, bool CCW, + FloatType epsilon) { + VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); + buildMesh(vertexDataSource, CCW, false, epsilon); + return HalfEdgeMesh(m_mesh, m_vertexData); } +template +void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, + bool useOriginalIndices, T epsilon) { + // CCW is unused for now + (void)CCW; + // useOriginalIndices is unused for now + (void)useOriginalIndices; + + if (pointCloud.size() == 0) { + m_mesh = MeshBuilder(); + return; + } + m_vertexData = pointCloud; + + // Very first: find extreme values and use them to compute the scale of the + // point cloud. + m_extremeValues = getExtremeValues(); + m_scale = getScale(m_extremeValues); + + // Epsilon we use depends on the scale + m_epsilon = epsilon * m_scale; + m_epsilonSquared = m_epsilon * m_epsilon; + + // Reset diagnostics + m_diagnostics = DiagnosticsData(); + + m_planar = false; // The planar case happens when all the points appear to + // lie on a two dimensional subspace of R^3. + createConvexHalfEdgeMesh(); + if (m_planar) { + const size_t extraPointIndex = m_planarPointCloudTemp.size() - 1; + for (auto& he : m_mesh.m_halfEdges) { + if (he.m_endVertex == extraPointIndex) { + he.m_endVertex = 0; + } + } + m_vertexData = pointCloud; + m_planarPointCloudTemp.clear(); + } +} + +template +ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, + bool CCW, bool useOriginalIndices, + T epsilon) { + buildMesh(pointCloud, CCW, useOriginalIndices, epsilon); + return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); +} + +template +void QuickHull::createConvexHalfEdgeMesh() { + m_visibleFaces.clear(); + m_horizonEdges.clear(); + m_possiblyVisibleFaces.clear(); + + // Compute base tetrahedron + setupInitialTetrahedron(); + assert(m_mesh.m_faces.size() == 4); + + // Init face stack with those faces that have points assigned to them + m_faceList.clear(); + for (size_t i = 0; i < 4; i++) { + auto& f = m_mesh.m_faces[i]; + if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size() > 0) { + m_faceList.push_back(i); + f.m_inFaceStack = 1; + } + } + + // Process faces until the face list is empty. + size_t iter = 0; + while (!m_faceList.empty()) { + iter++; + if (iter == std::numeric_limits::max()) { + // Visible face traversal marks visited faces with iteration counter (to + // mark that the face has been visited on this iteration) and the max + // value represents unvisited faces. At this point we have to reset + // iteration counter. This shouldn't be an issue on 64 bit machines. + iter = 0; + } + + const size_t topFaceIndex = m_faceList.front(); + m_faceList.pop_front(); + + auto& tf = m_mesh.m_faces[topFaceIndex]; + tf.m_inFaceStack = 0; + + assert(!tf.m_pointsOnPositiveSide || tf.m_pointsOnPositiveSide->size() > 0); + if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) { + continue; + } + + // Pick the most distant point to this triangle plane as the point to which + // we extrude + const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint]; + const size_t activePointIndex = tf.m_mostDistantPoint; + + // Find out the faces that have our active point on their positive side + // (these are the "visible faces"). The face on top of the stack of course + // is one of them. At the same time, we create a list of horizon edges. + m_horizonEdges.clear(); + m_possiblyVisibleFaces.clear(); + m_visibleFaces.clear(); + m_possiblyVisibleFaces.emplace_back(topFaceIndex, + std::numeric_limits::max()); + while (m_possiblyVisibleFaces.size()) { + const auto faceData = m_possiblyVisibleFaces.back(); + m_possiblyVisibleFaces.pop_back(); + auto& pvf = m_mesh.m_faces[faceData.m_faceIndex]; + assert(!pvf.isDisabled()); + + if (pvf.m_visibilityCheckedOnIteration == iter) { + if (pvf.m_isVisibleFaceOnCurrentIteration) { + continue; + } + } else { + const Plane& P = pvf.m_P; + pvf.m_visibilityCheckedOnIteration = iter; + const T d = glm::dot(P.m_N, activePoint) + P.m_D; + if (d > 0) { + pvf.m_isVisibleFaceOnCurrentIteration = 1; + pvf.m_horizonEdgesOnCurrentIteration = 0; + m_visibleFaces.push_back(faceData.m_faceIndex); + for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) { + if (m_mesh.m_halfEdges[heIndex].m_opp != + faceData.m_enteredFromHalfEdge) { + m_possiblyVisibleFaces.emplace_back( + m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face, + heIndex); + } + } + continue; + } + assert(faceData.m_faceIndex != topFaceIndex); + } + + // The face is not visible. Therefore, the halfedge we came from is part + // of the horizon edge. + pvf.m_isVisibleFaceOnCurrentIteration = 0; + m_horizonEdges.push_back(faceData.m_enteredFromHalfEdge); + // Store which half edge is the horizon edge. The other half edges of the + // face will not be part of the final mesh so their data slots can by + // recycled. + const auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace( + m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge] + .m_face]); + const std::int8_t ind = + (halfEdges[0] == faceData.m_enteredFromHalfEdge) + ? 0 + : (halfEdges[1] == faceData.m_enteredFromHalfEdge ? 1 : 2); + m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face] + .m_horizonEdgesOnCurrentIteration |= (1 << ind); + } + const size_t horizonEdgeCount = m_horizonEdges.size(); + + // Order horizon edges so that they form a loop. This may fail due to + // 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(m_horizonEdges)) { + m_diagnostics.m_failedHorizonEdges++; + std::cerr << "Failed to solve horizon edge." << std::endl; + auto it = std::find(tf.m_pointsOnPositiveSide->begin(), + tf.m_pointsOnPositiveSide->end(), activePointIndex); + tf.m_pointsOnPositiveSide->erase(it); + if (tf.m_pointsOnPositiveSide->size() == 0) { + reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide); + } + continue; + } + + // Except for the horizon edges, all half edges of the visible faces can be + // marked as disabled. Their data slots will be reused. The faces will be + // disabled as well, but we need to remember the points that were on the + // positive side of them - therefore we save pointers to them. + m_newFaceIndices.clear(); + m_newHalfEdgeIndices.clear(); + m_disabledFacePointVectors.clear(); + size_t disableCounter = 0; + for (auto faceIndex : m_visibleFaces) { + auto& disabledFace = m_mesh.m_faces[faceIndex]; + auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace); + for (size_t j = 0; j < 3; j++) { + if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1 << j)) == 0) { + if (disableCounter < horizonEdgeCount * 2) { + // Use on this iteration + m_newHalfEdgeIndices.push_back(halfEdges[j]); + disableCounter++; + } else { + // Mark for reusal on later iteration step + m_mesh.disableHalfEdge(halfEdges[j]); + } + } + } + // Disable the face, but retain pointer to the points that were on the + // positive side of it. We need to assign those points to the new faces we + // create shortly. + auto t = m_mesh.disableFace(faceIndex); + if (t) { + assert(t->size()); // Because we should not assign point vectors to + // faces unless needed... + m_disabledFacePointVectors.push_back(std::move(t)); + } + } + if (disableCounter < horizonEdgeCount * 2) { + const size_t newHalfEdgesNeeded = horizonEdgeCount * 2 - disableCounter; + for (size_t i = 0; i < newHalfEdgesNeeded; i++) { + m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge()); + } + } + + // Create new faces using the edgeloop + for (size_t i = 0; i < horizonEdgeCount; i++) { + const size_t AB = m_horizonEdges[i]; + + auto horizonEdgeVertexIndices = + m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]); + size_t A, B, C; + A = horizonEdgeVertexIndices[0]; + B = horizonEdgeVertexIndices[1]; + C = activePointIndex; + + const size_t newFaceIndex = m_mesh.addFace(); + m_newFaceIndices.push_back(newFaceIndex); + + const size_t CA = m_newHalfEdgeIndices[2 * i + 0]; + const size_t BC = m_newHalfEdgeIndices[2 * i + 1]; + + m_mesh.m_halfEdges[AB].m_next = BC; + m_mesh.m_halfEdges[BC].m_next = CA; + m_mesh.m_halfEdges[CA].m_next = AB; + + m_mesh.m_halfEdges[BC].m_face = newFaceIndex; + m_mesh.m_halfEdges[CA].m_face = newFaceIndex; + m_mesh.m_halfEdges[AB].m_face = newFaceIndex; + + m_mesh.m_halfEdges[CA].m_endVertex = A; + m_mesh.m_halfEdges[BC].m_endVertex = C; + + auto& newFace = m_mesh.m_faces[newFaceIndex]; + + const glm::vec3 planeNormal = mathutils::getTriangleNormal( + m_vertexData[A], m_vertexData[B], activePoint); + newFace.m_P = Plane(planeNormal, activePoint); + newFace.m_he = AB; + + m_mesh.m_halfEdges[CA].m_opp = + m_newHalfEdgeIndices[i > 0 ? i * 2 - 1 : 2 * horizonEdgeCount - 1]; + m_mesh.m_halfEdges[BC].m_opp = + m_newHalfEdgeIndices[((i + 1) * 2) % (horizonEdgeCount * 2)]; + } + + // Assign points that were on the positive side of the disabled faces to the + // new faces. + for (auto& disabledPoints : m_disabledFacePointVectors) { + assert(disabledPoints); + for (const auto& point : *(disabledPoints)) { + if (point == activePointIndex) { + continue; + } + for (size_t j = 0; j < horizonEdgeCount; j++) { + if (addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]], point)) { + break; + } + } + } + // The points are no longer needed: we can move them to the vector pool + // for reuse. + reclaimToIndexVectorPool(disabledPoints); + } + + // Increase face stack size if needed + for (const auto newFaceIndex : m_newFaceIndices) { + auto& newFace = m_mesh.m_faces[newFaceIndex]; + if (newFace.m_pointsOnPositiveSide) { + assert(newFace.m_pointsOnPositiveSide->size() > 0); + if (!newFace.m_inFaceStack) { + m_faceList.push_back(newFaceIndex); + newFace.m_inFaceStack = 1; + } + } + } + } + + // Cleanup + m_indexVectorPool.clear(); +} + +/* + * Private helper functions + */ + +template +std::array QuickHull::getExtremeValues() { + std::array outIndices{0, 0, 0, 0, 0, 0}; + T extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, m_vertexData[0].y, + m_vertexData[0].y, m_vertexData[0].z, m_vertexData[0].z}; + const size_t vCount = m_vertexData.size(); + for (size_t i = 1; i < vCount; i++) { + const glm::vec3& pos = m_vertexData[i]; + if (pos.x > extremeVals[0]) { + extremeVals[0] = pos.x; + outIndices[0] = i; + } else if (pos.x < extremeVals[1]) { + extremeVals[1] = pos.x; + outIndices[1] = i; + } + if (pos.y > extremeVals[2]) { + extremeVals[2] = pos.y; + outIndices[2] = i; + } else if (pos.y < extremeVals[3]) { + extremeVals[3] = pos.y; + outIndices[3] = i; + } + if (pos.z > extremeVals[4]) { + extremeVals[4] = pos.z; + outIndices[4] = i; + } else if (pos.z < extremeVals[5]) { + extremeVals[5] = pos.z; + outIndices[5] = i; + } + } + return outIndices; +} + +template +bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { + const size_t horizonEdgeCount = horizonEdges.size(); + for (size_t i = 0; i < horizonEdgeCount - 1; i++) { + const size_t endVertex = m_mesh.m_halfEdges[horizonEdges[i]].m_endVertex; + bool foundNext = false; + for (size_t j = i + 1; j < horizonEdgeCount; j++) { + const size_t beginVertex = + m_mesh.m_halfEdges[m_mesh.m_halfEdges[horizonEdges[j]].m_opp] + .m_endVertex; + if (beginVertex == endVertex) { + std::swap(horizonEdges[i + 1], horizonEdges[j]); + foundNext = true; + break; + } + } + if (!foundNext) { + return false; + } + } + assert( + m_mesh.m_halfEdges[horizonEdges[horizonEdges.size() - 1]].m_endVertex == + m_mesh.m_halfEdges[m_mesh.m_halfEdges[horizonEdges[0]].m_opp] + .m_endVertex); + return true; +} + +template +T QuickHull::getScale(const std::array& extremeValues) { + T s = 0; + for (size_t i = 0; i < 6; i++) { + const T* v = (const T*)(&m_vertexData[extremeValues[i]]); + v += i / 2; + auto a = std::abs(*v); + if (a > s) { + s = a; + } + } + return s; +} + +template +void QuickHull::setupInitialTetrahedron() { + const size_t vertexCount = m_vertexData.size(); + + // If we have at most 4 points, just return a degenerate tetrahedron: + if (vertexCount <= 4) { + size_t v[4] = {0, std::min((size_t)1, vertexCount - 1), + std::min((size_t)2, vertexCount - 1), + std::min((size_t)3, vertexCount - 1)}; + const glm::vec3 N = mathutils::getTriangleNormal( + m_vertexData[v[0]], m_vertexData[v[1]], m_vertexData[v[2]]); + const Plane trianglePlane(N, m_vertexData[v[0]]); + if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { + std::swap(v[0], v[1]); + } + return m_mesh.setup(v[0], v[1], v[2], v[3]); + } + + // Find two most distant extreme points. + T maxD = m_epsilonSquared; + std::pair selectedPoints; + for (size_t i = 0; i < 6; i++) { + for (size_t j = i + 1; j < 6; j++) { + // I found a function for squaredDistance but i can't seem to include it + // like this for some reason + const T d = mathutils::getSquaredDistance( + m_vertexData[m_extremeValues[i]], m_vertexData[m_extremeValues[j]]); + if (d > maxD) { + maxD = d; + selectedPoints = {m_extremeValues[i], m_extremeValues[j]}; + } + } + } + if (maxD == m_epsilonSquared) { + // A degenerate case: the point cloud seems to consists of a single point + return m_mesh.setup(0, std::min((size_t)1, vertexCount - 1), + std::min((size_t)2, vertexCount - 1), + std::min((size_t)3, vertexCount - 1)); + } + assert(selectedPoints.first != selectedPoints.second); + + // Find the most distant point to the line between the two chosen extreme + // points. + const Ray r(m_vertexData[selectedPoints.first], + (m_vertexData[selectedPoints.second] - + m_vertexData[selectedPoints.first])); + maxD = m_epsilonSquared; + size_t maxI = std::numeric_limits::max(); + const size_t vCount = m_vertexData.size(); + for (size_t i = 0; i < vCount; i++) { + const T distToRay = + mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i], r); + if (distToRay > maxD) { + maxD = distToRay; + maxI = i; + } + } + if (maxD == m_epsilonSquared) { + // It appears that the point cloud belongs to a 1 dimensional subspace of + // R^3: convex hull has no volume => return a thin triangle Pick any point + // other than selectedPoints.first and selectedPoints.second as the third + // point of the triangle + auto it = std::find_if(m_vertexData.begin(), m_vertexData.end(), + [&](const vec3& ve) { + return ve != m_vertexData[selectedPoints.first] && + ve != m_vertexData[selectedPoints.second]; + }); + const size_t thirdPoint = (it == m_vertexData.end()) + ? selectedPoints.first + : std::distance(m_vertexData.begin(), it); + it = std::find_if(m_vertexData.begin(), m_vertexData.end(), + [&](const vec3& ve) { + return ve != m_vertexData[selectedPoints.first] && + ve != m_vertexData[selectedPoints.second] && + ve != m_vertexData[thirdPoint]; + }); + const size_t fourthPoint = (it == m_vertexData.end()) + ? selectedPoints.first + : std::distance(m_vertexData.begin(), it); + return m_mesh.setup(selectedPoints.first, selectedPoints.second, thirdPoint, + fourthPoint); + } + + // These three points form the base triangle for our tetrahedron. + assert(selectedPoints.first != maxI && selectedPoints.second != maxI); + std::array baseTriangle{selectedPoints.first, + selectedPoints.second, maxI}; + const glm::vec3 baseTriangleVertices[] = {m_vertexData[baseTriangle[0]], + m_vertexData[baseTriangle[1]], + m_vertexData[baseTriangle[2]]}; + + // Next step is to find the 4th vertex of the tetrahedron. We naturally choose + // the point farthest away from the triangle plane. + maxD = m_epsilon; + maxI = 0; + const glm::vec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], + baseTriangleVertices[1], + baseTriangleVertices[2]); + Plane trianglePlane(N, baseTriangleVertices[0]); + for (size_t i = 0; i < vCount; i++) { + const T d = std::abs( + mathutils::getSignedDistanceToPlane(m_vertexData[i], trianglePlane)); + if (d > maxD) { + maxD = d; + maxI = i; + } + } + if (maxD == m_epsilon) { + // All the points seem to lie on a 2D subspace of R^3. How to handle this? + // Well, let's add one extra point to the point cloud so that the convex + // hull will have volume. + m_planar = true; + const vec3 N1 = mathutils::getTriangleNormal(baseTriangleVertices[1], + baseTriangleVertices[2], + baseTriangleVertices[0]); + m_planarPointCloudTemp.clear(); + m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(), + m_vertexData.begin(), m_vertexData.end()); + const vec3 extraPoint = N1 + m_vertexData[0]; + m_planarPointCloudTemp.push_back(extraPoint); + maxI = m_planarPointCloudTemp.size() - 1; + m_vertexData = VertexDataSource(m_planarPointCloudTemp); + } + + // Enforce CCW orientation (if user prefers clockwise orientation, swap two + // vertices in each triangle when final mesh is created) + const Plane triPlane(N, baseTriangleVertices[0]); + if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) { + std::swap(baseTriangle[0], baseTriangle[1]); + } + + // Create a tetrahedron half edge mesh and compute planes defined by each + // triangle + m_mesh.setup(baseTriangle[0], baseTriangle[1], baseTriangle[2], maxI); + for (auto& f : m_mesh.m_faces) { + auto v = m_mesh.getVertexIndicesOfFace(f); + const glm::vec3& va = m_vertexData[v[0]]; + const glm::vec3& vb = m_vertexData[v[1]]; + const glm::vec3& vc = m_vertexData[v[2]]; + const glm::vec3 N1 = mathutils::getTriangleNormal(va, vb, vc); + const Plane plane(N1, va); + f.m_P = plane; + } + + // Finally we assign a face for each vertex outside the tetrahedron (vertices + // inside the tetrahedron have no role anymore) + for (size_t i = 0; i < vCount; i++) { + for (auto& face : m_mesh.m_faces) { + if (addPointToFace(face, i)) { + break; + } + } + } +} + +/* + * Explicit template specifications for float and double + */ + +template class QuickHull; +template class QuickHull; +} // namespace quickhull diff --git a/src/manifold/src/quickhull.hpp b/src/manifold/src/quickhull.hpp index f3096ce54..27a6af3d7 100644 --- a/src/manifold/src/quickhull.hpp +++ b/src/manifold/src/quickhull.hpp @@ -12,6 +12,8 @@ #include #include +#include "par.h" + // Pool.hpp @@ -43,142 +45,23 @@ namespace quickhull { // Vector3.hpp - - - template - class Vector3 - { - public: - Vector3() = default; - - Vector3(T x, T y, T z) : x(x), y(y), z(z) { - - } - - T x,y,z; - - T dotProduct(const Vector3& other) const { - return x*other.x+y*other.y+z*other.z; - } - - void normalize() { - const T len = getLength(); - x/=len; - y/=len; - z/=len; - } - - Vector3 getNormalized() const { - const T len = getLength(); - return Vector3(x/len,y/len,z/len); - } - T getLength() const { - return std::sqrt(x*x+y*y+z*z); - } - - Vector3 operator-(const Vector3& other) const { - return Vector3(x-other.x,y-other.y,z-other.z); - } - - Vector3 operator+(const Vector3& other) const { - return Vector3(x+other.x,y+other.y,z+other.z); - } - - Vector3& operator+=(const Vector3& other) { - x+=other.x; - y+=other.y; - z+=other.z; - return *this; - } - Vector3& operator-=(const Vector3& other) { - x-=other.x; - y-=other.y; - z-=other.z; - return *this; - } - Vector3& operator*=(T c) { - x*=c; - y*=c; - z*=c; - return *this; - } - - Vector3& operator/=(T c) { - x/=c; - y/=c; - z/=c; - return *this; - } - Vector3 operator-() const { - return Vector3(-x,-y,-z); - } - template - Vector3 operator*(S c) const { - return Vector3(x*c,y*c,z*c); - } - - template - Vector3 operator/(S c) const { - return Vector3(x/c,y/c,z/c); - } - - T getLengthSquared() const { - return x*x + y*y + z*z; - } - - bool operator!=(const Vector3& o) const { - return x != o.x || y != o.y || z != o.z; - } - - // Projection onto another vector - Vector3 projection(const Vector3& o) const { - T C = dotProduct(o)/o.getLengthSquared(); - return o*C; - } - - Vector3 crossProduct (const Vector3& rhs ) { - T a = y * rhs.z - z * rhs.y ; - T b = z * rhs.x - x * rhs.z ; - T c = x * rhs.y - y * rhs.x ; - Vector3 product( a , b , c ) ; - return product ; - } - - T getDistanceTo(const Vector3& other) const { - Vector3 diff = *this - other; - return diff.getLength(); - } - - T getSquaredDistanceTo(const Vector3& other) const { - const T dx = x-other.x; - const T dy = y-other.y; - const T dz = z-other.z; - return dx*dx+dy*dy+dz*dz; - } + // // Projection onto another vector + // glm::vec3 projection(const glm::vec3& o) const { + // T C = dotProduct(o)/o.getLengthSquared(); + // return o*C; + // } + - }; - - // Overload also << operator for easy printing of debug data - template - inline std::ostream& operator<<(std::ostream& os, const Vector3& vec) { - os << "(" << vec.x << "," << vec.y << "," << vec.z << ")"; - return os; - } - - template - inline Vector3 operator*(T c, const Vector3& v) { - return Vector3(v.x*c,v.y*c,v.z*c); - } // Plane.hpp template class Plane { public: - Vector3 m_N; + glm::vec3 m_N; // Signed distance (if normal is of length 1) to the plane from origin T m_D; @@ -186,8 +69,8 @@ namespace quickhull { // Normal length squared T m_sqrNLength; - bool isPointOnPositiveSide(const Vector3& Q) const { - T d = m_N.dotProduct(Q)+m_D; + bool isPointOnPositiveSide(const glm::vec3& Q) const { + T d = glm::dot(m_N,Q)+m_D; if (d>=0) return true; return false; } @@ -195,7 +78,7 @@ namespace quickhull { Plane() = default; // Construct a plane using normal N and any point P on the plane - Plane(const Vector3& N, const Vector3& P) : m_N(N), m_D(-N.dotProduct(P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { + Plane(const glm::vec3& N, const glm::vec3& P) : m_N(N), m_D(glm::dot(-N,P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { } }; @@ -204,11 +87,11 @@ namespace quickhull { template struct Ray { - const Vector3 m_S; - const Vector3 m_V; + const glm::vec3 m_S; + const glm::vec3 m_V; const T m_VInvLengthSquared; - Ray(const Vector3& S,const Vector3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/m_V.getLengthSquared()) { + Ray(const glm::vec3& S,const glm::vec3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/(m_V.x*m_V.x+m_V.y*m_V.y+m_V.z*m_V.z)) { } }; @@ -216,17 +99,17 @@ namespace quickhull { // VertexDataSource - template + // template class VertexDataSource { - const Vector3* m_ptr; + const glm::vec3* m_ptr; size_t m_count; public: - VertexDataSource(const Vector3* ptr, size_t count) : m_ptr(ptr), m_count(count) { + VertexDataSource(const glm::vec3* ptr, size_t count) : m_ptr(ptr), m_count(count) { } - VertexDataSource(const std::vector>& vec) : m_ptr(&vec[0]), m_count(vec.size()) { + VertexDataSource(const std::vector& vec) : m_ptr(&vec[0]), m_count(vec.size()) { } @@ -240,15 +123,15 @@ namespace quickhull { return m_count; } - const Vector3& operator[](size_t index) const { + const glm::vec3& operator[](size_t index) const { return m_ptr[index]; } - const Vector3* begin() const { + const glm::vec3* begin() const { return m_ptr; } - const Vector3* end() const { + const glm::vec3* end() const { return m_ptr + m_count; } }; @@ -497,8 +380,8 @@ namespace quickhull { template class ConvexHull { - std::unique_ptr>> m_optimizedVertexBuffer; - VertexDataSource m_vertices; + std::unique_ptr> m_optimizedVertexBuffer; + VertexDataSource m_vertices; std::vector m_indices; public: ConvexHull() {} @@ -507,8 +390,8 @@ namespace quickhull { ConvexHull(const ConvexHull& o) { m_indices = o.m_indices; if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector>(*o.m_optimizedVertexBuffer)); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { m_vertices = o.m_vertices; @@ -521,8 +404,8 @@ namespace quickhull { } m_indices = o.m_indices; if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector>(*o.m_optimizedVertexBuffer)); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { m_vertices = o.m_vertices; @@ -534,8 +417,8 @@ namespace quickhull { m_indices = std::move(o.m_indices); if (o.m_optimizedVertexBuffer) { m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); - o.m_vertices = VertexDataSource(); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + o.m_vertices = VertexDataSource(); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { m_vertices = o.m_vertices; @@ -549,8 +432,8 @@ namespace quickhull { m_indices = std::move(o.m_indices); if (o.m_optimizedVertexBuffer) { m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); - o.m_vertices = VertexDataSource(); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + o.m_vertices = VertexDataSource(); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { m_vertices = o.m_vertices; @@ -559,9 +442,9 @@ namespace quickhull { } // Construct vertex and index buffers from half edge mesh and pointcloud - ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { + ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { if (!useOriginalIndices) { - m_optimizedVertexBuffer.reset(new std::vector>()); + m_optimizedVertexBuffer.reset(new std::vector()); } std::vector faceProcessed(mesh.m_faces.size(),false); @@ -619,7 +502,7 @@ namespace quickhull { } if (!useOriginalIndices) { - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); + m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { m_vertices = pointCloud; @@ -634,11 +517,11 @@ namespace quickhull { return m_indices; } - VertexDataSource& getVertexBuffer() { + VertexDataSource& getVertexBuffer() { return m_vertices; } - const VertexDataSource& getVertexBuffer() const { + const VertexDataSource& getVertexBuffer() const { return m_vertices; } @@ -680,11 +563,11 @@ namespace quickhull { IndexType m_halfEdgeIndex; // Index of one of the half edges of this face }; - std::vector> m_vertices; + std::vector m_vertices; std::vector m_faces; std::vector m_halfEdges; - HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) + HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) { std::unordered_map faceMapping; std::unordered_map halfEdgeMapping; @@ -739,31 +622,33 @@ namespace quickhull { namespace mathutils { template - inline T getSquaredDistanceBetweenPointAndRay(const Vector3& p, const Ray& r) { - const Vector3 s = p-r.m_S; - T t = s.dotProduct(r.m_V); - return s.getLengthSquared() - t*t*r.m_VInvLengthSquared; + inline T getSquaredDistanceBetweenPointAndRay(const glm::vec3& p, const Ray& r) { + const glm::vec3 s = p-r.m_S; + T t = glm::dot(s,r.m_V); + return (s.x*s.x+s.y*s.y+s.z*s.z) - t*t*r.m_VInvLengthSquared; + } + + inline float getSquaredDistance(const glm::vec3& p1, const glm::vec3& p2) { + return ((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z)); } - // Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance). template - inline T getSignedDistanceToPlane(const Vector3& v, const Plane& p) { - return p.m_N.dotProduct(v) + p.m_D; + inline T getSignedDistanceToPlane(const glm::vec3& v, const Plane& p) { + return glm::dot(p.m_N,v) + p.m_D; } - template - inline Vector3 getTriangleNormal(const Vector3& a,const Vector3& b,const Vector3& c) { + inline glm::vec3 getTriangleNormal(const glm::vec3& a,const glm::vec3& b,const glm::vec3& c) { // We want to get (a-c).crossProduct(b-c) without constructing temp vectors - T x = a.x - c.x; - T y = a.y - c.y; - T z = a.z - c.z; - T rhsx = b.x - c.x; - T rhsy = b.y - c.y; - T rhsz = b.z - c.z; - T px = y * rhsz - z * rhsy ; - T py = z * rhsx - x * rhsz ; - T pz = x * rhsy - y * rhsx ; - return Vector3(px,py,pz); + float x = a.x - c.x; + float y = a.y - c.y; + float z = a.z - c.z; + float rhsx = b.x - c.x; + float rhsy = b.y - c.y; + float rhsz = b.z - c.z; + float px = y * rhsz - z * rhsy ; + float py = z * rhsx - x * rhsz ; + float pz = x * rhsy - y * rhsx ; + return glm::vec3(px,py,pz); } @@ -819,12 +704,12 @@ namespace quickhull { template class QuickHull { - using vec3 = Vector3; + using vec3 = glm::vec3; FloatType m_epsilon, m_epsilonSquared, m_scale; bool m_planar; std::vector m_planarPointCloudTemp; - VertexDataSource m_vertexData; + VertexDataSource m_vertexData; MeshBuilder m_mesh; std::array m_extremeValues; DiagnosticsData m_diagnostics; @@ -870,10 +755,10 @@ namespace quickhull { void createConvexHalfEdgeMesh(); // Constructs the convex hull into a MeshBuilder object which can be converted to a ConvexHull or Mesh object - void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); // The public getConvexHull functions will setup a VertexDataSource object and call this - ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); public: // Computes convex hull for a given point cloud. // Params: @@ -882,7 +767,7 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const std::vector>& pointCloud, + ConvexHull getConvexHull(const std::vector& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps = defaultEps()); @@ -895,7 +780,7 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const Vector3* vertexData, + ConvexHull getConvexHull(const glm::vec3* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, From 0f5eb39f3f00426271fbf76041fad97e691a8ece Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 1 Aug 2024 01:38:58 +0530 Subject: [PATCH 03/24] Revert nanobind --- bindings/python/third_party/nanobind | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bindings/python/third_party/nanobind b/bindings/python/third_party/nanobind index 1a309ba44..8d7f1ee06 160000 --- a/bindings/python/third_party/nanobind +++ b/bindings/python/third_party/nanobind @@ -1 +1 @@ -Subproject commit 1a309ba444a47e081dc6213d72345a2fbbd20795 +Subproject commit 8d7f1ee0621c17fa370b704b2100ffa0243d5bfb From 219905c0d5037b6d2182acb5c08b6f49ca0df224 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 1 Aug 2024 02:17:15 +0530 Subject: [PATCH 04/24] Removing Templates (since using glm::vec3 so only need float) --- src/manifold/src/manifold.cpp | 2 +- src/manifold/src/quickhull.cpp | 95 +++++++++++----------------- src/manifold/src/quickhull.hpp | 110 ++++++++++++++------------------- 3 files changed, 85 insertions(+), 122 deletions(-) diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index ae97ad835..7d2dddef4 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -924,7 +924,7 @@ Manifold Manifold::Hull(const std::vector& pts) { vertices[i] = {pts[i].x, pts[i].y, pts[i].z}; } - quickhull::QuickHull qh; + quickhull::QuickHull qh; // bools: correct triangle winding, and use original indices auto hull = qh.getConvexHull(vertices, false, true); const auto& triangles = hull.getIndexBuffer(); diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 4a6120caf..aedd13d8b 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -9,63 +9,52 @@ namespace quickhull { -template <> float defaultEps() { return 0.0001f; } -template <> -double defaultEps() { - return 0.0000001; -} - /* * Implementation of the algorithm */ -template -ConvexHull QuickHull::getConvexHull( +ConvexHull QuickHull::getConvexHull( const std::vector& pointCloud, bool CCW, bool useOriginalIndices, - T epsilon) { + float epsilon) { VertexDataSource vertexDataSource(pointCloud); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -template -ConvexHull QuickHull::getConvexHull(const glm::vec3* vertexData, +ConvexHull QuickHull::getConvexHull(const glm::vec3* vertexData, size_t vertexCount, bool CCW, - bool useOriginalIndices, T epsilon) { + bool useOriginalIndices, float epsilon) { VertexDataSource vertexDataSource(vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -template -ConvexHull QuickHull::getConvexHull(const T* vertexData, +ConvexHull QuickHull::getConvexHull(const float* vertexData, size_t vertexCount, bool CCW, - bool useOriginalIndices, T epsilon) { + bool useOriginalIndices, float epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -template -HalfEdgeMesh QuickHull::getConvexHullAsMesh( - const FloatType* vertexData, size_t vertexCount, bool CCW, - FloatType epsilon) { +HalfEdgeMesh QuickHull::getConvexHullAsMesh( + const float* vertexData, size_t vertexCount, bool CCW, + float epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); buildMesh(vertexDataSource, CCW, false, epsilon); - return HalfEdgeMesh(m_mesh, m_vertexData); + return HalfEdgeMesh(m_mesh, m_vertexData); } -template -void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, - bool useOriginalIndices, T epsilon) { +void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, + bool useOriginalIndices, float epsilon) { // CCW is unused for now (void)CCW; // useOriginalIndices is unused for now (void)useOriginalIndices; if (pointCloud.size() == 0) { - m_mesh = MeshBuilder(); + m_mesh = MeshBuilder(); return; } m_vertexData = pointCloud; @@ -97,16 +86,14 @@ void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, } } -template -ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, +ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, - T epsilon) { + float epsilon) { buildMesh(pointCloud, CCW, useOriginalIndices, epsilon); - return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); + return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); } -template -void QuickHull::createConvexHalfEdgeMesh() { +void QuickHull::createConvexHalfEdgeMesh() { m_visibleFaces.clear(); m_horizonEdges.clear(); m_possiblyVisibleFaces.clear(); @@ -172,9 +159,9 @@ void QuickHull::createConvexHalfEdgeMesh() { continue; } } else { - const Plane& P = pvf.m_P; + const Plane& P = pvf.m_P; pvf.m_visibilityCheckedOnIteration = iter; - const T d = glm::dot(P.m_N, activePoint) + P.m_D; + const float d = glm::dot(P.m_N, activePoint) + P.m_D; if (d > 0) { pvf.m_isVisibleFaceOnCurrentIteration = 1; pvf.m_horizonEdgesOnCurrentIteration = 0; @@ -298,7 +285,7 @@ void QuickHull::createConvexHalfEdgeMesh() { const glm::vec3 planeNormal = mathutils::getTriangleNormal( m_vertexData[A], m_vertexData[B], activePoint); - newFace.m_P = Plane(planeNormal, activePoint); + newFace.m_P = Plane(planeNormal, activePoint); newFace.m_he = AB; m_mesh.m_halfEdges[CA].m_opp = @@ -347,10 +334,9 @@ void QuickHull::createConvexHalfEdgeMesh() { * Private helper functions */ -template -std::array QuickHull::getExtremeValues() { +std::array QuickHull::getExtremeValues() { std::array outIndices{0, 0, 0, 0, 0, 0}; - T extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, m_vertexData[0].y, + float extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, m_vertexData[0].y, m_vertexData[0].y, m_vertexData[0].z, m_vertexData[0].z}; const size_t vCount = m_vertexData.size(); for (size_t i = 1; i < vCount; i++) { @@ -380,8 +366,7 @@ std::array QuickHull::getExtremeValues() { return outIndices; } -template -bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { +bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { const size_t horizonEdgeCount = horizonEdges.size(); for (size_t i = 0; i < horizonEdgeCount - 1; i++) { const size_t endVertex = m_mesh.m_halfEdges[horizonEdges[i]].m_endVertex; @@ -407,11 +392,10 @@ bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { return true; } -template -T QuickHull::getScale(const std::array& extremeValues) { - T s = 0; +float QuickHull::getScale(const std::array& extremeValues) { + float s = 0; for (size_t i = 0; i < 6; i++) { - const T* v = (const T*)(&m_vertexData[extremeValues[i]]); + const float* v = (const float*)(&m_vertexData[extremeValues[i]]); v += i / 2; auto a = std::abs(*v); if (a > s) { @@ -421,8 +405,7 @@ T QuickHull::getScale(const std::array& extremeValues) { return s; } -template -void QuickHull::setupInitialTetrahedron() { +void QuickHull::setupInitialTetrahedron() { const size_t vertexCount = m_vertexData.size(); // If we have at most 4 points, just return a degenerate tetrahedron: @@ -432,7 +415,7 @@ void QuickHull::setupInitialTetrahedron() { std::min((size_t)3, vertexCount - 1)}; const glm::vec3 N = mathutils::getTriangleNormal( m_vertexData[v[0]], m_vertexData[v[1]], m_vertexData[v[2]]); - const Plane trianglePlane(N, m_vertexData[v[0]]); + const Plane trianglePlane(N, m_vertexData[v[0]]); if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { std::swap(v[0], v[1]); } @@ -440,13 +423,13 @@ void QuickHull::setupInitialTetrahedron() { } // Find two most distant extreme points. - T maxD = m_epsilonSquared; + float maxD = m_epsilonSquared; std::pair selectedPoints; for (size_t i = 0; i < 6; i++) { for (size_t j = i + 1; j < 6; j++) { // I found a function for squaredDistance but i can't seem to include it // like this for some reason - const T d = mathutils::getSquaredDistance( + const float d = mathutils::getSquaredDistance( m_vertexData[m_extremeValues[i]], m_vertexData[m_extremeValues[j]]); if (d > maxD) { maxD = d; @@ -464,14 +447,14 @@ void QuickHull::setupInitialTetrahedron() { // Find the most distant point to the line between the two chosen extreme // points. - const Ray r(m_vertexData[selectedPoints.first], + const Ray r(m_vertexData[selectedPoints.first], (m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first])); maxD = m_epsilonSquared; size_t maxI = std::numeric_limits::max(); const size_t vCount = m_vertexData.size(); for (size_t i = 0; i < vCount; i++) { - const T distToRay = + const float distToRay = mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i], r); if (distToRay > maxD) { maxD = distToRay; @@ -519,9 +502,9 @@ void QuickHull::setupInitialTetrahedron() { const glm::vec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], baseTriangleVertices[1], baseTriangleVertices[2]); - Plane trianglePlane(N, baseTriangleVertices[0]); + Plane trianglePlane(N, baseTriangleVertices[0]); for (size_t i = 0; i < vCount; i++) { - const T d = std::abs( + const float d = std::abs( mathutils::getSignedDistanceToPlane(m_vertexData[i], trianglePlane)); if (d > maxD) { maxD = d; @@ -547,7 +530,7 @@ void QuickHull::setupInitialTetrahedron() { // Enforce CCW orientation (if user prefers clockwise orientation, swap two // vertices in each triangle when final mesh is created) - const Plane triPlane(N, baseTriangleVertices[0]); + const Plane triPlane(N, baseTriangleVertices[0]); if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) { std::swap(baseTriangle[0], baseTriangle[1]); } @@ -561,7 +544,7 @@ void QuickHull::setupInitialTetrahedron() { const glm::vec3& vb = m_vertexData[v[1]]; const glm::vec3& vc = m_vertexData[v[2]]; const glm::vec3 N1 = mathutils::getTriangleNormal(va, vb, vc); - const Plane plane(N1, va); + const Plane plane(N1, va); f.m_P = plane; } @@ -576,10 +559,4 @@ void QuickHull::setupInitialTetrahedron() { } } -/* - * Explicit template specifications for float and double - */ - -template class QuickHull; -template class QuickHull; } // namespace quickhull diff --git a/src/manifold/src/quickhull.hpp b/src/manifold/src/quickhull.hpp index 27a6af3d7..d32374627 100644 --- a/src/manifold/src/quickhull.hpp +++ b/src/manifold/src/quickhull.hpp @@ -19,24 +19,23 @@ namespace quickhull { - template class Pool { - std::vector> m_data; + std::vector>> m_data; public: void clear() { m_data.clear(); } - void reclaim(std::unique_ptr& ptr) { + void reclaim(std::unique_ptr>& ptr) { m_data.push_back(std::move(ptr)); } - std::unique_ptr get() { + std::unique_ptr> get() { if (m_data.size()==0) { - return std::unique_ptr(new T()); + return std::unique_ptr>(new std::vector()); } auto it = m_data.end()-1; - std::unique_ptr r = std::move(*it); + std::unique_ptr> r = std::move(*it); m_data.erase(it); return r; } @@ -58,19 +57,18 @@ namespace quickhull { // Plane.hpp - template class Plane { public: glm::vec3 m_N; // Signed distance (if normal is of length 1) to the plane from origin - T m_D; + float m_D; // Normal length squared - T m_sqrNLength; + float m_sqrNLength; bool isPointOnPositiveSide(const glm::vec3& Q) const { - T d = glm::dot(m_N,Q)+m_D; + float d = glm::dot(m_N,Q)+m_D; if (d>=0) return true; return false; } @@ -85,11 +83,10 @@ namespace quickhull { // Ray.hpp - template struct Ray { const glm::vec3 m_S; const glm::vec3 m_V; - const T m_VInvLengthSquared; + const float m_VInvLengthSquared; Ray(const glm::vec3& S,const glm::vec3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/(m_V.x*m_V.x+m_V.y*m_V.y+m_V.z*m_V.z)) { } @@ -98,8 +95,7 @@ namespace quickhull { // VertexDataSource - - // template + class VertexDataSource { const glm::vec3* m_ptr; size_t m_count; @@ -141,7 +137,6 @@ namespace quickhull { // Mesh.hpp - template class MeshBuilder { public: struct HalfEdge { @@ -161,8 +156,8 @@ namespace quickhull { struct Face { size_t m_he; - Plane m_P{}; - T m_mostDistantPointDist; + Plane m_P{}; + float m_mostDistantPointDist; size_t m_mostDistantPoint; size_t m_visibilityCheckedOnIteration; std::uint8_t m_isVisibleFaceOnCurrentIteration : 1; @@ -378,7 +373,6 @@ namespace quickhull { // ConvexHull.hpp - template class ConvexHull { std::unique_ptr> m_optimizedVertexBuffer; VertexDataSource m_vertices; @@ -442,7 +436,7 @@ namespace quickhull { } // Construct vertex and index buffers from half edge mesh and pointcloud - ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { + ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { if (!useOriginalIndices) { m_optimizedVertexBuffer.reset(new std::vector()); } @@ -548,40 +542,39 @@ namespace quickhull { // HalfEdgeMesh.hpp - template class HalfEdgeMesh { public: struct HalfEdge { - IndexType m_endVertex; - IndexType m_opp; - IndexType m_face; - IndexType m_next; + size_t m_endVertex; + size_t m_opp; + size_t m_face; + size_t m_next; }; struct Face { - IndexType m_halfEdgeIndex; // Index of one of the half edges of this face + size_t m_halfEdgeIndex; // Index of one of the half edges of this face }; std::vector m_vertices; std::vector m_faces; std::vector m_halfEdges; - HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) + HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) { - std::unordered_map faceMapping; - std::unordered_map halfEdgeMapping; - std::unordered_map vertexMapping; + std::unordered_map faceMapping; + std::unordered_map halfEdgeMapping; + std::unordered_map vertexMapping; size_t i=0; for (const auto& face : builderObject.m_faces) { if (!face.isDisabled()) { - m_faces.push_back({static_cast(face.m_he)}); + m_faces.push_back({static_cast(face.m_he)}); faceMapping[i] = m_faces.size()-1; const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face); for (const auto heIndex : heIndices) { - const IndexType vertexIndex = builderObject.m_halfEdges[heIndex].m_endVertex; + const size_t vertexIndex = builderObject.m_halfEdges[heIndex].m_endVertex; if (vertexMapping.count(vertexIndex)==0) { m_vertices.push_back(vertexData[vertexIndex]); vertexMapping[vertexIndex] = m_vertices.size()-1; @@ -594,7 +587,7 @@ namespace quickhull { i=0; for (const auto& halfEdge : builderObject.m_halfEdges) { if (!halfEdge.isDisabled()) { - m_halfEdges.push_back({static_cast(halfEdge.m_endVertex),static_cast(halfEdge.m_opp),static_cast(halfEdge.m_face),static_cast(halfEdge.m_next)}); + m_halfEdges.push_back({static_cast(halfEdge.m_endVertex),static_cast(halfEdge.m_opp),static_cast(halfEdge.m_face),static_cast(halfEdge.m_next)}); halfEdgeMapping[i] = m_halfEdges.size()-1; } i++; @@ -621,10 +614,9 @@ namespace quickhull { namespace mathutils { - template - inline T getSquaredDistanceBetweenPointAndRay(const glm::vec3& p, const Ray& r) { + inline float getSquaredDistanceBetweenPointAndRay(const glm::vec3& p, const Ray& r) { const glm::vec3 s = p-r.m_S; - T t = glm::dot(s,r.m_V); + float t = glm::dot(s,r.m_V); return (s.x*s.x+s.y*s.y+s.z*s.z) - t*t*r.m_VInvLengthSquared; } @@ -632,8 +624,7 @@ namespace quickhull { return ((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z)); } // Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance). - template - inline T getSignedDistanceToPlane(const glm::vec3& v, const Plane& p) { + inline float getSignedDistanceToPlane(const glm::vec3& v, const Plane& p) { return glm::dot(p.m_N,v) + p.m_D; } @@ -699,18 +690,16 @@ namespace quickhull { DiagnosticsData() : m_failedHorizonEdges(0) { } }; - template - FloatType defaultEps(); + float defaultEps(); - template class QuickHull { using vec3 = glm::vec3; - FloatType m_epsilon, m_epsilonSquared, m_scale; + float m_epsilon, m_epsilonSquared, m_scale; bool m_planar; std::vector m_planarPointCloudTemp; VertexDataSource m_vertexData; - MeshBuilder m_mesh; + MeshBuilder m_mesh; std::array m_extremeValues; DiagnosticsData m_diagnostics; @@ -738,27 +727,27 @@ namespace quickhull { std::array getExtremeValues(); // Compute scale of the vertex data. - FloatType getScale(const std::array& extremeValues); + float getScale(const std::array& extremeValues); // Each face contains a unique pointer to a vector of indices. However, many - often most - faces do not have any points on the positive // side of them especially at the the end of the iteration. When a face is removed from the mesh, its associated point vector, if such // exists, is moved to the index vector pool, and when we need to add new faces with points on the positive side to the mesh, // we reuse these vectors. This reduces the amount of std::vectors we have to deal with, and impact on performance is remarkable. - Pool> m_indexVectorPool; + Pool m_indexVectorPool; inline std::unique_ptr> getIndexVectorFromPool(); inline void reclaimToIndexVectorPool(std::unique_ptr>& ptr); // Associates a point with a face if the point resides on the positive side of the plane. Returns true if the points was on the positive side. - inline bool addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex); + inline bool addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex); // This will update m_mesh from which we create the ConvexHull object that getConvexHull function returns void createConvexHalfEdgeMesh(); // Constructs the convex hull into a MeshBuilder object which can be converted to a ConvexHull or Mesh object - void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, float eps); // The public getConvexHull functions will setup a VertexDataSource object and call this - ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps); + ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, float eps); public: // Computes convex hull for a given point cloud. // Params: @@ -767,10 +756,10 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const std::vector& pointCloud, + ConvexHull getConvexHull(const std::vector& pointCloud, bool CCW, bool useOriginalIndices, - FloatType eps = defaultEps()); + float eps = defaultEps()); // Computes convex hull for a given point cloud. // Params: @@ -780,11 +769,11 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const glm::vec3* vertexData, + ConvexHull getConvexHull(const glm::vec3* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, - FloatType eps = defaultEps()); + float eps = defaultEps()); // 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,... @@ -795,11 +784,11 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const FloatType* vertexData, + ConvexHull getConvexHull(const float* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, - FloatType eps = defaultEps()); + float eps = defaultEps()); // 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,... @@ -810,10 +799,10 @@ namespace quickhull { // 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 a mesh object with half edge structure. - HalfEdgeMesh getConvexHullAsMesh(const FloatType* vertexData, + HalfEdgeMesh getConvexHullAsMesh(const float* vertexData, size_t vertexCount, bool CCW, - FloatType eps = defaultEps()); + float eps = defaultEps()); // Get diagnostics about last generated convex hull const DiagnosticsData& getDiagnostics() { @@ -825,15 +814,13 @@ namespace quickhull { * Inline function definitions */ - template - std::unique_ptr> QuickHull::getIndexVectorFromPool() { + std::unique_ptr> QuickHull::getIndexVectorFromPool() { auto r = m_indexVectorPool.get(); r->clear(); return r; } - template - void QuickHull::reclaimToIndexVectorPool(std::unique_ptr>& ptr) { + void QuickHull::reclaimToIndexVectorPool(std::unique_ptr>& ptr) { const size_t oldSize = ptr->size(); if ((oldSize+1)*128 < ptr->capacity()) { // Reduce memory usage! Huge vectors are needed at the beginning of iteration when faces have many points on their positive side. Later on, smaller vectors will suffice. @@ -843,9 +830,8 @@ namespace quickhull { m_indexVectorPool.reclaim(ptr); } - template - bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { - const T D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); + bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { + const float D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) { if (!f.m_pointsOnPositiveSide) { f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool()); From 511c1539bf7d6e9d6309ffeaadda1c54bcf83885 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 1 Aug 2024 02:17:51 +0530 Subject: [PATCH 05/24] Formatting --- bindings/python/manifold3d.cpp | 9 ++++++-- src/manifold/src/quickhull.cpp | 41 +++++++++++++++++----------------- 2 files changed, 27 insertions(+), 23 deletions(-) diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index 4dee244f4..f8cb5de04 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -257,7 +257,10 @@ NB_MODULE(manifold3d, m) { manifold__translate__v) .def("scale", &Manifold::Scale, nb::arg("v"), manifold__scale__v) .def( - "scale", [](const Manifold &m, float s) { m.Scale({s, s, s}); }, + "scale", + [](const Manifold &m, float s) { + m.Scale({s, s, s}); + }, nb::arg("s"), "Scale this Manifold in space. This operation can be chained. " "Transforms are combined and applied lazily.\n\n" @@ -670,7 +673,9 @@ NB_MODULE(manifold3d, m) { cross_section__scale__scale) .def( "scale", - [](const CrossSection &self, float s) { self.Scale({s, s}); }, + [](const CrossSection &self, float s) { + self.Scale({s, s}); + }, nb::arg("s"), "Scale this CrossSection in space. This operation can be chained. " "Transforms are combined and applied lazily." diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index aedd13d8b..9a695890b 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -9,45 +9,43 @@ namespace quickhull { -float defaultEps() { - return 0.0001f; -} +float defaultEps() { return 0.0001f; } /* * Implementation of the algorithm */ -ConvexHull QuickHull::getConvexHull( - const std::vector& pointCloud, bool CCW, bool useOriginalIndices, - float epsilon) { +ConvexHull QuickHull::getConvexHull(const std::vector& pointCloud, + bool CCW, bool useOriginalIndices, + float epsilon) { VertexDataSource vertexDataSource(pointCloud); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } ConvexHull QuickHull::getConvexHull(const glm::vec3* vertexData, - size_t vertexCount, bool CCW, - bool useOriginalIndices, float epsilon) { + size_t vertexCount, bool CCW, + bool useOriginalIndices, float epsilon) { VertexDataSource vertexDataSource(vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -ConvexHull QuickHull::getConvexHull(const float* vertexData, - size_t vertexCount, bool CCW, - bool useOriginalIndices, float epsilon) { +ConvexHull QuickHull::getConvexHull(const float* vertexData, size_t vertexCount, + bool CCW, bool useOriginalIndices, + float epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -HalfEdgeMesh QuickHull::getConvexHullAsMesh( - const float* vertexData, size_t vertexCount, bool CCW, - float epsilon) { +HalfEdgeMesh QuickHull::getConvexHullAsMesh(const float* vertexData, + size_t vertexCount, bool CCW, + float epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); buildMesh(vertexDataSource, CCW, false, epsilon); return HalfEdgeMesh(m_mesh, m_vertexData); } void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, - bool useOriginalIndices, float epsilon) { + bool useOriginalIndices, float epsilon) { // CCW is unused for now (void)CCW; // useOriginalIndices is unused for now @@ -87,8 +85,8 @@ void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, } ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, - bool CCW, bool useOriginalIndices, - float epsilon) { + bool CCW, bool useOriginalIndices, + float epsilon) { buildMesh(pointCloud, CCW, useOriginalIndices, epsilon); return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); } @@ -336,8 +334,9 @@ void QuickHull::createConvexHalfEdgeMesh() { std::array QuickHull::getExtremeValues() { std::array outIndices{0, 0, 0, 0, 0, 0}; - float extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, m_vertexData[0].y, - m_vertexData[0].y, m_vertexData[0].z, m_vertexData[0].z}; + float extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, + m_vertexData[0].y, m_vertexData[0].y, + m_vertexData[0].z, m_vertexData[0].z}; const size_t vCount = m_vertexData.size(); for (size_t i = 1; i < vCount; i++) { const glm::vec3& pos = m_vertexData[i]; @@ -448,8 +447,8 @@ void QuickHull::setupInitialTetrahedron() { // Find the most distant point to the line between the two chosen extreme // points. const Ray r(m_vertexData[selectedPoints.first], - (m_vertexData[selectedPoints.second] - - m_vertexData[selectedPoints.first])); + (m_vertexData[selectedPoints.second] - + m_vertexData[selectedPoints.first])); maxD = m_epsilonSquared; size_t maxI = std::numeric_limits::max(); const size_t vCount = m_vertexData.size(); From 4c3389fd99ec3e11f2c473096bd0b90b57597fae Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 1 Aug 2024 14:07:36 +0530 Subject: [PATCH 06/24] Converted glm::vec3 to glm::vec3d and float to double --- src/manifold/src/manifold.cpp | 2 +- src/manifold/src/quickhull.cpp | 58 +++++++++--------- src/manifold/src/quickhull.hpp | 108 ++++++++++++++++----------------- 3 files changed, 84 insertions(+), 84 deletions(-) diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 7d2dddef4..a34c0bd5c 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -919,7 +919,7 @@ Manifold Manifold::Hull(const std::vector& pts) { const int numVert = pts.size(); if (numVert < 4) return Manifold(); - std::vector vertices(numVert); + std::vector vertices(numVert); for (int i = 0; i < numVert; i++) { vertices[i] = {pts[i].x, pts[i].y, pts[i].z}; } diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 9a695890b..26f2126d4 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -9,43 +9,43 @@ namespace quickhull { -float defaultEps() { return 0.0001f; } +double defaultEps() { return 0.0000001; } /* * Implementation of the algorithm */ -ConvexHull QuickHull::getConvexHull(const std::vector& pointCloud, +ConvexHull QuickHull::getConvexHull(const std::vector& pointCloud, bool CCW, bool useOriginalIndices, - float epsilon) { + double epsilon) { VertexDataSource vertexDataSource(pointCloud); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -ConvexHull QuickHull::getConvexHull(const glm::vec3* vertexData, +ConvexHull QuickHull::getConvexHull(const glm::dvec3* vertexData, size_t vertexCount, bool CCW, - bool useOriginalIndices, float epsilon) { + bool useOriginalIndices, double epsilon) { VertexDataSource vertexDataSource(vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -ConvexHull QuickHull::getConvexHull(const float* vertexData, size_t vertexCount, +ConvexHull QuickHull::getConvexHull(const double* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, - float epsilon) { + double epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -HalfEdgeMesh QuickHull::getConvexHullAsMesh(const float* vertexData, +HalfEdgeMesh QuickHull::getConvexHullAsMesh(const double* vertexData, size_t vertexCount, bool CCW, - float epsilon) { + double epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); buildMesh(vertexDataSource, CCW, false, epsilon); return HalfEdgeMesh(m_mesh, m_vertexData); } void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, - bool useOriginalIndices, float epsilon) { + bool useOriginalIndices, double epsilon) { // CCW is unused for now (void)CCW; // useOriginalIndices is unused for now @@ -86,7 +86,7 @@ void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, - float epsilon) { + double epsilon) { buildMesh(pointCloud, CCW, useOriginalIndices, epsilon); return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); } @@ -159,7 +159,7 @@ void QuickHull::createConvexHalfEdgeMesh() { } else { const Plane& P = pvf.m_P; pvf.m_visibilityCheckedOnIteration = iter; - const float d = glm::dot(P.m_N, activePoint) + P.m_D; + const double d = glm::dot(P.m_N, activePoint) + P.m_D; if (d > 0) { pvf.m_isVisibleFaceOnCurrentIteration = 1; pvf.m_horizonEdgesOnCurrentIteration = 0; @@ -281,7 +281,7 @@ void QuickHull::createConvexHalfEdgeMesh() { auto& newFace = m_mesh.m_faces[newFaceIndex]; - const glm::vec3 planeNormal = mathutils::getTriangleNormal( + const glm::dvec3 planeNormal = mathutils::getTriangleNormal( m_vertexData[A], m_vertexData[B], activePoint); newFace.m_P = Plane(planeNormal, activePoint); newFace.m_he = AB; @@ -334,12 +334,12 @@ void QuickHull::createConvexHalfEdgeMesh() { std::array QuickHull::getExtremeValues() { std::array outIndices{0, 0, 0, 0, 0, 0}; - float extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, + double extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, m_vertexData[0].y, m_vertexData[0].y, m_vertexData[0].z, m_vertexData[0].z}; const size_t vCount = m_vertexData.size(); for (size_t i = 1; i < vCount; i++) { - const glm::vec3& pos = m_vertexData[i]; + const glm::dvec3& pos = m_vertexData[i]; if (pos.x > extremeVals[0]) { extremeVals[0] = pos.x; outIndices[0] = i; @@ -391,10 +391,10 @@ bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { return true; } -float QuickHull::getScale(const std::array& extremeValues) { - float s = 0; +double QuickHull::getScale(const std::array& extremeValues) { + double s = 0; for (size_t i = 0; i < 6; i++) { - const float* v = (const float*)(&m_vertexData[extremeValues[i]]); + const double* v = (const double*)(&m_vertexData[extremeValues[i]]); v += i / 2; auto a = std::abs(*v); if (a > s) { @@ -412,7 +412,7 @@ void QuickHull::setupInitialTetrahedron() { size_t v[4] = {0, std::min((size_t)1, vertexCount - 1), std::min((size_t)2, vertexCount - 1), std::min((size_t)3, vertexCount - 1)}; - const glm::vec3 N = mathutils::getTriangleNormal( + const glm::dvec3 N = mathutils::getTriangleNormal( m_vertexData[v[0]], m_vertexData[v[1]], m_vertexData[v[2]]); const Plane trianglePlane(N, m_vertexData[v[0]]); if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { @@ -422,13 +422,13 @@ void QuickHull::setupInitialTetrahedron() { } // Find two most distant extreme points. - float maxD = m_epsilonSquared; + double maxD = m_epsilonSquared; std::pair selectedPoints; for (size_t i = 0; i < 6; i++) { for (size_t j = i + 1; j < 6; j++) { // I found a function for squaredDistance but i can't seem to include it // like this for some reason - const float d = mathutils::getSquaredDistance( + const double d = mathutils::getSquaredDistance( m_vertexData[m_extremeValues[i]], m_vertexData[m_extremeValues[j]]); if (d > maxD) { maxD = d; @@ -453,7 +453,7 @@ void QuickHull::setupInitialTetrahedron() { size_t maxI = std::numeric_limits::max(); const size_t vCount = m_vertexData.size(); for (size_t i = 0; i < vCount; i++) { - const float distToRay = + const double distToRay = mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i], r); if (distToRay > maxD) { maxD = distToRay; @@ -490,7 +490,7 @@ void QuickHull::setupInitialTetrahedron() { assert(selectedPoints.first != maxI && selectedPoints.second != maxI); std::array baseTriangle{selectedPoints.first, selectedPoints.second, maxI}; - const glm::vec3 baseTriangleVertices[] = {m_vertexData[baseTriangle[0]], + const glm::dvec3 baseTriangleVertices[] = {m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]], m_vertexData[baseTriangle[2]]}; @@ -498,12 +498,12 @@ void QuickHull::setupInitialTetrahedron() { // the point farthest away from the triangle plane. maxD = m_epsilon; maxI = 0; - const glm::vec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], + const glm::dvec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], baseTriangleVertices[1], baseTriangleVertices[2]); Plane trianglePlane(N, baseTriangleVertices[0]); for (size_t i = 0; i < vCount; i++) { - const float d = std::abs( + const double d = std::abs( mathutils::getSignedDistanceToPlane(m_vertexData[i], trianglePlane)); if (d > maxD) { maxD = d; @@ -539,10 +539,10 @@ void QuickHull::setupInitialTetrahedron() { m_mesh.setup(baseTriangle[0], baseTriangle[1], baseTriangle[2], maxI); for (auto& f : m_mesh.m_faces) { auto v = m_mesh.getVertexIndicesOfFace(f); - const glm::vec3& va = m_vertexData[v[0]]; - const glm::vec3& vb = m_vertexData[v[1]]; - const glm::vec3& vc = m_vertexData[v[2]]; - const glm::vec3 N1 = mathutils::getTriangleNormal(va, vb, vc); + const glm::dvec3& va = m_vertexData[v[0]]; + const glm::dvec3& vb = m_vertexData[v[1]]; + const glm::dvec3& vc = m_vertexData[v[2]]; + const glm::dvec3 N1 = mathutils::getTriangleNormal(va, vb, vc); const Plane plane(N1, va); f.m_P = plane; } diff --git a/src/manifold/src/quickhull.hpp b/src/manifold/src/quickhull.hpp index d32374627..0de13b0e2 100644 --- a/src/manifold/src/quickhull.hpp +++ b/src/manifold/src/quickhull.hpp @@ -48,7 +48,7 @@ namespace quickhull { // // Projection onto another vector - // glm::vec3 projection(const glm::vec3& o) const { + // glm::dvec3 projection(const glm::dvec3& o) const { // T C = dotProduct(o)/o.getLengthSquared(); // return o*C; // } @@ -59,16 +59,16 @@ namespace quickhull { class Plane { public: - glm::vec3 m_N; + glm::dvec3 m_N; // Signed distance (if normal is of length 1) to the plane from origin - float m_D; + double m_D; // Normal length squared - float m_sqrNLength; + double m_sqrNLength; - bool isPointOnPositiveSide(const glm::vec3& Q) const { - float d = glm::dot(m_N,Q)+m_D; + bool isPointOnPositiveSide(const glm::dvec3& Q) const { + double d = glm::dot(m_N,Q)+m_D; if (d>=0) return true; return false; } @@ -76,7 +76,7 @@ namespace quickhull { Plane() = default; // Construct a plane using normal N and any point P on the plane - Plane(const glm::vec3& N, const glm::vec3& P) : m_N(N), m_D(glm::dot(-N,P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { + Plane(const glm::dvec3& N, const glm::dvec3& P) : m_N(N), m_D(glm::dot(-N,P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { } }; @@ -84,11 +84,11 @@ namespace quickhull { // Ray.hpp struct Ray { - const glm::vec3 m_S; - const glm::vec3 m_V; - const float m_VInvLengthSquared; + const glm::dvec3 m_S; + const glm::dvec3 m_V; + const double m_VInvLengthSquared; - Ray(const glm::vec3& S,const glm::vec3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/(m_V.x*m_V.x+m_V.y*m_V.y+m_V.z*m_V.z)) { + Ray(const glm::dvec3& S,const glm::dvec3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/(m_V.x*m_V.x+m_V.y*m_V.y+m_V.z*m_V.z)) { } }; @@ -97,15 +97,15 @@ namespace quickhull { class VertexDataSource { - const glm::vec3* m_ptr; + const glm::dvec3* m_ptr; size_t m_count; public: - VertexDataSource(const glm::vec3* ptr, size_t count) : m_ptr(ptr), m_count(count) { + VertexDataSource(const glm::dvec3* ptr, size_t count) : m_ptr(ptr), m_count(count) { } - VertexDataSource(const std::vector& vec) : m_ptr(&vec[0]), m_count(vec.size()) { + VertexDataSource(const std::vector& vec) : m_ptr(&vec[0]), m_count(vec.size()) { } @@ -119,15 +119,15 @@ namespace quickhull { return m_count; } - const glm::vec3& operator[](size_t index) const { + const glm::dvec3& operator[](size_t index) const { return m_ptr[index]; } - const glm::vec3* begin() const { + const glm::dvec3* begin() const { return m_ptr; } - const glm::vec3* end() const { + const glm::dvec3* end() const { return m_ptr + m_count; } }; @@ -157,7 +157,7 @@ namespace quickhull { struct Face { size_t m_he; Plane m_P{}; - float m_mostDistantPointDist; + double m_mostDistantPointDist; size_t m_mostDistantPoint; size_t m_visibilityCheckedOnIteration; std::uint8_t m_isVisibleFaceOnCurrentIteration : 1; @@ -374,7 +374,7 @@ namespace quickhull { // ConvexHull.hpp class ConvexHull { - std::unique_ptr> m_optimizedVertexBuffer; + std::unique_ptr> m_optimizedVertexBuffer; VertexDataSource m_vertices; std::vector m_indices; public: @@ -384,7 +384,7 @@ namespace quickhull { ConvexHull(const ConvexHull& o) { m_indices = o.m_indices; if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); + m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { @@ -398,7 +398,7 @@ namespace quickhull { } m_indices = o.m_indices; if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); + m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); m_vertices = VertexDataSource(*m_optimizedVertexBuffer); } else { @@ -438,7 +438,7 @@ namespace quickhull { // Construct vertex and index buffers from half edge mesh and pointcloud ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { if (!useOriginalIndices) { - m_optimizedVertexBuffer.reset(new std::vector()); + m_optimizedVertexBuffer.reset(new std::vector()); } std::vector faceProcessed(mesh.m_faces.size(),false); @@ -556,7 +556,7 @@ namespace quickhull { size_t m_halfEdgeIndex; // Index of one of the half edges of this face }; - std::vector m_vertices; + std::vector m_vertices; std::vector m_faces; std::vector m_halfEdges; @@ -614,32 +614,32 @@ namespace quickhull { namespace mathutils { - inline float getSquaredDistanceBetweenPointAndRay(const glm::vec3& p, const Ray& r) { - const glm::vec3 s = p-r.m_S; - float t = glm::dot(s,r.m_V); + inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p, const Ray& r) { + const glm::dvec3 s = p-r.m_S; + double t = glm::dot(s,r.m_V); return (s.x*s.x+s.y*s.y+s.z*s.z) - t*t*r.m_VInvLengthSquared; } - inline float getSquaredDistance(const glm::vec3& p1, const glm::vec3& p2) { + inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { return ((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z)); } // Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance). - inline float getSignedDistanceToPlane(const glm::vec3& v, const Plane& p) { + inline double getSignedDistanceToPlane(const glm::dvec3& v, const Plane& p) { return glm::dot(p.m_N,v) + p.m_D; } - inline glm::vec3 getTriangleNormal(const glm::vec3& a,const glm::vec3& b,const glm::vec3& c) { + inline glm::dvec3 getTriangleNormal(const glm::dvec3& a,const glm::dvec3& b,const glm::dvec3& c) { // We want to get (a-c).crossProduct(b-c) without constructing temp vectors - float x = a.x - c.x; - float y = a.y - c.y; - float z = a.z - c.z; - float rhsx = b.x - c.x; - float rhsy = b.y - c.y; - float rhsz = b.z - c.z; - float px = y * rhsz - z * rhsy ; - float py = z * rhsx - x * rhsz ; - float pz = x * rhsy - y * rhsx ; - return glm::vec3(px,py,pz); + double x = a.x - c.x; + double y = a.y - c.y; + double z = a.z - c.z; + double rhsx = b.x - c.x; + double rhsy = b.y - c.y; + double rhsz = b.z - c.z; + double px = y * rhsz - z * rhsy ; + double py = z * rhsx - x * rhsz ; + double pz = x * rhsy - y * rhsx ; + return glm::dvec3(px,py,pz); } @@ -690,12 +690,12 @@ namespace quickhull { DiagnosticsData() : m_failedHorizonEdges(0) { } }; - float defaultEps(); + double defaultEps(); class QuickHull { - using vec3 = glm::vec3; + using vec3 = glm::dvec3; - float m_epsilon, m_epsilonSquared, m_scale; + double m_epsilon, m_epsilonSquared, m_scale; bool m_planar; std::vector m_planarPointCloudTemp; VertexDataSource m_vertexData; @@ -727,7 +727,7 @@ namespace quickhull { std::array getExtremeValues(); // Compute scale of the vertex data. - float getScale(const std::array& extremeValues); + double getScale(const std::array& extremeValues); // Each face contains a unique pointer to a vector of indices. However, many - often most - faces do not have any points on the positive // side of them especially at the the end of the iteration. When a face is removed from the mesh, its associated point vector, if such @@ -744,10 +744,10 @@ namespace quickhull { void createConvexHalfEdgeMesh(); // Constructs the convex hull into a MeshBuilder object which can be converted to a ConvexHull or Mesh object - void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, float eps); + void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, double eps); // The public getConvexHull functions will setup a VertexDataSource object and call this - ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, float eps); + ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, double eps); public: // Computes convex hull for a given point cloud. // Params: @@ -756,10 +756,10 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const std::vector& pointCloud, + ConvexHull getConvexHull(const std::vector& pointCloud, bool CCW, bool useOriginalIndices, - float eps = defaultEps()); + double eps = defaultEps()); // Computes convex hull for a given point cloud. // Params: @@ -769,11 +769,11 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const glm::vec3* vertexData, + ConvexHull getConvexHull(const glm::dvec3* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, - float eps = defaultEps()); + double eps = defaultEps()); // 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,... @@ -784,11 +784,11 @@ namespace quickhull { // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const float* vertexData, + ConvexHull getConvexHull(const double* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, - float eps = defaultEps()); + double eps = defaultEps()); // 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,... @@ -799,10 +799,10 @@ namespace quickhull { // 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 a mesh object with half edge structure. - HalfEdgeMesh getConvexHullAsMesh(const float* vertexData, + HalfEdgeMesh getConvexHullAsMesh(const double* vertexData, size_t vertexCount, bool CCW, - float eps = defaultEps()); + double eps = defaultEps()); // Get diagnostics about last generated convex hull const DiagnosticsData& getDiagnostics() { @@ -831,7 +831,7 @@ namespace quickhull { } bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { - const float D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); + const double D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) { if (!f.m_pointsOnPositiveSide) { f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool()); From a5f74ca968dc9fdd7e8f65e63369dc3fc40898f5 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 1 Aug 2024 14:11:04 +0530 Subject: [PATCH 07/24] Fixed Formatting and Removed third_party --- bindings/python/CMakeLists.txt | 1 - src/CMakeLists.txt | 1 - src/manifold/src/quickhull.cpp | 18 +++++++++--------- src/third_party/CMakeLists.txt | 1 - 4 files changed, 9 insertions(+), 12 deletions(-) delete mode 100644 src/third_party/CMakeLists.txt diff --git a/bindings/python/CMakeLists.txt b/bindings/python/CMakeLists.txt index 2b5f90f59..579a7c6a7 100644 --- a/bindings/python/CMakeLists.txt +++ b/bindings/python/CMakeLists.txt @@ -14,7 +14,6 @@ project(python) -add_subdirectory(third_party) nanobind_add_module( manifold3d NB_STATIC STABLE_ABI LTO diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 49ec7a3da..5abab30f4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,7 +13,6 @@ # limitations under the License. # do not add third_party installations -add_subdirectory(third_party EXCLUDE_FROM_ALL) add_subdirectory(utilities) add_subdirectory(collider) add_subdirectory(cross_section) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 26f2126d4..ff9976ce4 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -29,9 +29,9 @@ ConvexHull QuickHull::getConvexHull(const glm::dvec3* vertexData, return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } -ConvexHull QuickHull::getConvexHull(const double* vertexData, size_t vertexCount, - bool CCW, bool useOriginalIndices, - double epsilon) { +ConvexHull QuickHull::getConvexHull(const double* vertexData, + size_t vertexCount, bool CCW, + bool useOriginalIndices, double epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); return getConvexHull(vertexDataSource, CCW, useOriginalIndices, epsilon); } @@ -335,8 +335,8 @@ void QuickHull::createConvexHalfEdgeMesh() { std::array QuickHull::getExtremeValues() { std::array outIndices{0, 0, 0, 0, 0, 0}; double extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, - m_vertexData[0].y, m_vertexData[0].y, - m_vertexData[0].z, m_vertexData[0].z}; + m_vertexData[0].y, m_vertexData[0].y, + m_vertexData[0].z, m_vertexData[0].z}; const size_t vCount = m_vertexData.size(); for (size_t i = 1; i < vCount; i++) { const glm::dvec3& pos = m_vertexData[i]; @@ -491,16 +491,16 @@ void QuickHull::setupInitialTetrahedron() { std::array baseTriangle{selectedPoints.first, selectedPoints.second, maxI}; const glm::dvec3 baseTriangleVertices[] = {m_vertexData[baseTriangle[0]], - m_vertexData[baseTriangle[1]], - m_vertexData[baseTriangle[2]]}; + m_vertexData[baseTriangle[1]], + m_vertexData[baseTriangle[2]]}; // Next step is to find the 4th vertex of the tetrahedron. We naturally choose // the point farthest away from the triangle plane. maxD = m_epsilon; maxI = 0; const glm::dvec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], - baseTriangleVertices[1], - baseTriangleVertices[2]); + baseTriangleVertices[1], + baseTriangleVertices[2]); Plane trianglePlane(N, baseTriangleVertices[0]); for (size_t i = 0; i < vCount; i++) { const double d = std::abs( diff --git a/src/third_party/CMakeLists.txt b/src/third_party/CMakeLists.txt deleted file mode 100644 index 8b1378917..000000000 --- a/src/third_party/CMakeLists.txt +++ /dev/null @@ -1 +0,0 @@ - From 302f1497df8e782ad9f372626778abeb04c5fb9b Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 01:17:15 +0530 Subject: [PATCH 08/24] Made requested Changes (m_,comments,removing unnecessary code) --- src/CMakeLists.txt | 1 - src/manifold/src/manifold.cpp | 4 +- src/manifold/src/quickhull.cpp | 416 ++++++++-------- src/manifold/src/quickhull.h | 824 +++++++++++++++++++++++++++++++ src/manifold/src/quickhull.hpp | 852 --------------------------------- 5 files changed, 1032 insertions(+), 1065 deletions(-) create mode 100644 src/manifold/src/quickhull.h delete mode 100644 src/manifold/src/quickhull.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5abab30f4..67e37c3ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,7 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -# do not add third_party installations add_subdirectory(utilities) add_subdirectory(collider) add_subdirectory(cross_section) diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index a34c0bd5c..817bcfffe 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -20,7 +20,7 @@ #include "csg_tree.h" #include "impl.h" #include "par.h" -#include "quickhull.hpp" +#include "quickhull.h" namespace { using namespace manifold; @@ -924,7 +924,7 @@ Manifold Manifold::Hull(const std::vector& pts) { vertices[i] = {pts[i].x, pts[i].y, pts[i].z}; } - quickhull::QuickHull qh; + QuickHull qh; // bools: correct triangle winding, and use original indices auto hull = qh.getConvexHull(vertices, false, true); const auto& triangles = hull.getIndexBuffer(); diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index ff9976ce4..d138d5197 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -1,13 +1,10 @@ -#include "quickhull.hpp" -// #include "MathUtils.hpp" +#include "quickhull.h" + #include #include #include #include #include -// #include "Structs/Mesh.hpp" - -namespace quickhull { double defaultEps() { return 0.0000001; } @@ -41,7 +38,7 @@ HalfEdgeMesh QuickHull::getConvexHullAsMesh(const double* vertexData, double epsilon) { VertexDataSource vertexDataSource((const vec3*)vertexData, vertexCount); buildMesh(vertexDataSource, CCW, false, epsilon); - return HalfEdgeMesh(m_mesh, m_vertexData); + return HalfEdgeMesh(mesh, originalVertexData); } void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, @@ -52,35 +49,35 @@ void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, (void)useOriginalIndices; if (pointCloud.size() == 0) { - m_mesh = MeshBuilder(); + mesh = MeshBuilder(); return; } - m_vertexData = pointCloud; + originalVertexData = pointCloud; // Very first: find extreme values and use them to compute the scale of the // point cloud. - m_extremeValues = getExtremeValues(); - m_scale = getScale(m_extremeValues); + extremeValues = getExtremeValues(); + scale = getScale(extremeValues); // Epsilon we use depends on the scale - m_epsilon = epsilon * m_scale; - m_epsilonSquared = m_epsilon * m_epsilon; + epsilon = epsilon * scale; + epsilonSquared = epsilon * epsilon; // Reset diagnostics - m_diagnostics = DiagnosticsData(); - - m_planar = false; // The planar case happens when all the points appear to - // lie on a two dimensional subspace of R^3. + diagnostics = DiagnosticsData(); + // The planar case happens when all the points appear to lie on a two + // dimensional subspace of R^3. + planar = false; createConvexHalfEdgeMesh(); - if (m_planar) { - const size_t extraPointIndex = m_planarPointCloudTemp.size() - 1; - for (auto& he : m_mesh.m_halfEdges) { - if (he.m_endVertex == extraPointIndex) { - he.m_endVertex = 0; + if (planar) { + const size_t extraPointIndex = planarPointCloudTemp.size() - 1; + for (auto& he : mesh.halfEdges) { + if (he.endVertex == extraPointIndex) { + he.endVertex = 0; } } - m_vertexData = pointCloud; - m_planarPointCloudTemp.clear(); + originalVertexData = pointCloud; + planarPointCloudTemp.clear(); } } @@ -88,31 +85,31 @@ ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, double epsilon) { buildMesh(pointCloud, CCW, useOriginalIndices, epsilon); - return ConvexHull(m_mesh, m_vertexData, CCW, useOriginalIndices); + return ConvexHull(mesh, originalVertexData, CCW, useOriginalIndices); } void QuickHull::createConvexHalfEdgeMesh() { - m_visibleFaces.clear(); - m_horizonEdges.clear(); - m_possiblyVisibleFaces.clear(); + visibleFaces.clear(); + horizonEdges.clear(); + possiblyVisibleFaces.clear(); // Compute base tetrahedron setupInitialTetrahedron(); - assert(m_mesh.m_faces.size() == 4); + assert(mesh.faces.size() == 4); // Init face stack with those faces that have points assigned to them - m_faceList.clear(); + faceList.clear(); for (size_t i = 0; i < 4; i++) { - auto& f = m_mesh.m_faces[i]; - if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size() > 0) { - m_faceList.push_back(i); - f.m_inFaceStack = 1; + auto& f = mesh.faces[i]; + if (f.pointsOnPositiveSide && f.pointsOnPositiveSide->size() > 0) { + faceList.push_back(i); + f.inFaceStack = 1; } } // Process faces until the face list is empty. size_t iter = 0; - while (!m_faceList.empty()) { + while (!faceList.empty()) { iter++; if (iter == std::numeric_limits::max()) { // Visible face traversal marks visited faces with iteration counter (to @@ -122,91 +119,88 @@ void QuickHull::createConvexHalfEdgeMesh() { iter = 0; } - const size_t topFaceIndex = m_faceList.front(); - m_faceList.pop_front(); + const size_t topFaceIndex = faceList.front(); + faceList.pop_front(); - auto& tf = m_mesh.m_faces[topFaceIndex]; - tf.m_inFaceStack = 0; + auto& tf = mesh.faces[topFaceIndex]; + tf.inFaceStack = 0; - assert(!tf.m_pointsOnPositiveSide || tf.m_pointsOnPositiveSide->size() > 0); - if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) { + assert(!tf.pointsOnPositiveSide || tf.pointsOnPositiveSide->size() > 0); + if (!tf.pointsOnPositiveSide || tf.isDisabled()) { continue; } // Pick the most distant point to this triangle plane as the point to which // we extrude - const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint]; - const size_t activePointIndex = tf.m_mostDistantPoint; + const vec3& activePoint = originalVertexData[tf.mostDistantPoint]; + const size_t activePointIndex = tf.mostDistantPoint; // Find out the faces that have our active point on their positive side // (these are the "visible faces"). The face on top of the stack of course // is one of them. At the same time, we create a list of horizon edges. - m_horizonEdges.clear(); - m_possiblyVisibleFaces.clear(); - m_visibleFaces.clear(); - m_possiblyVisibleFaces.emplace_back(topFaceIndex, - std::numeric_limits::max()); - while (m_possiblyVisibleFaces.size()) { - const auto faceData = m_possiblyVisibleFaces.back(); - m_possiblyVisibleFaces.pop_back(); - auto& pvf = m_mesh.m_faces[faceData.m_faceIndex]; + horizonEdges.clear(); + possiblyVisibleFaces.clear(); + visibleFaces.clear(); + possiblyVisibleFaces.emplace_back(topFaceIndex, + std::numeric_limits::max()); + while (possiblyVisibleFaces.size()) { + const auto faceData = possiblyVisibleFaces.back(); + possiblyVisibleFaces.pop_back(); + auto& pvf = mesh.faces[faceData.faceIndex]; assert(!pvf.isDisabled()); - if (pvf.m_visibilityCheckedOnIteration == iter) { - if (pvf.m_isVisibleFaceOnCurrentIteration) { + if (pvf.visibilityCheckedOnIteration == iter) { + if (pvf.isVisibleFaceOnCurrentIteration) { continue; } } else { - const Plane& P = pvf.m_P; - pvf.m_visibilityCheckedOnIteration = iter; - const double d = glm::dot(P.m_N, activePoint) + P.m_D; + const Plane& P = pvf.P; + pvf.visibilityCheckedOnIteration = iter; + const double d = glm::dot(P.N, activePoint) + P.D; if (d > 0) { - pvf.m_isVisibleFaceOnCurrentIteration = 1; - pvf.m_horizonEdgesOnCurrentIteration = 0; - m_visibleFaces.push_back(faceData.m_faceIndex); - for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) { - if (m_mesh.m_halfEdges[heIndex].m_opp != - faceData.m_enteredFromHalfEdge) { - m_possiblyVisibleFaces.emplace_back( - m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face, - heIndex); + pvf.isVisibleFaceOnCurrentIteration = 1; + pvf.horizonEdgesOnCurrentIteration = 0; + visibleFaces.push_back(faceData.faceIndex); + for (auto heIndex : mesh.getHalfEdgeIndicesOfFace(pvf)) { + if (mesh.halfEdges[heIndex].opp != faceData.enteredFromHalfEdge) { + possiblyVisibleFaces.emplace_back( + mesh.halfEdges[mesh.halfEdges[heIndex].opp].face, heIndex); } } continue; } - assert(faceData.m_faceIndex != topFaceIndex); + assert(faceData.faceIndex != topFaceIndex); } // The face is not visible. Therefore, the halfedge we came from is part // of the horizon edge. - pvf.m_isVisibleFaceOnCurrentIteration = 0; - m_horizonEdges.push_back(faceData.m_enteredFromHalfEdge); + pvf.isVisibleFaceOnCurrentIteration = 0; + horizonEdges.push_back(faceData.enteredFromHalfEdge); // Store which half edge is the horizon edge. The other half edges of the // face will not be part of the final mesh so their data slots can by // recycled. - const auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace( - m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge] - .m_face]); + const auto halfEdges = mesh.getHalfEdgeIndicesOfFace( + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face]); const std::int8_t ind = - (halfEdges[0] == faceData.m_enteredFromHalfEdge) + (halfEdges[0] == faceData.enteredFromHalfEdge) ? 0 - : (halfEdges[1] == faceData.m_enteredFromHalfEdge ? 1 : 2); - m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face] - .m_horizonEdgesOnCurrentIteration |= (1 << ind); + : (halfEdges[1] == faceData.enteredFromHalfEdge ? 1 : 2); + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face] + .horizonEdgesOnCurrentIteration |= (1 << ind); } - const size_t horizonEdgeCount = m_horizonEdges.size(); + const size_t horizonEdgeCount = horizonEdges.size(); // Order horizon edges so that they form a loop. This may fail due to // 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(m_horizonEdges)) { - m_diagnostics.m_failedHorizonEdges++; + if (!reorderHorizonEdges(horizonEdges)) { + diagnostics.failedHorizonEdges++; std::cerr << "Failed to solve horizon edge." << std::endl; - auto it = std::find(tf.m_pointsOnPositiveSide->begin(), - tf.m_pointsOnPositiveSide->end(), activePointIndex); - tf.m_pointsOnPositiveSide->erase(it); - if (tf.m_pointsOnPositiveSide->size() == 0) { - reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide); + auto it = std::find(tf.pointsOnPositiveSide->begin(), + tf.pointsOnPositiveSide->end(), activePointIndex); + tf.pointsOnPositiveSide->erase(it); + if (tf.pointsOnPositiveSide->size() == 0) { + reclaimToIndexVectorPool(tf.pointsOnPositiveSide); } continue; } @@ -215,93 +209,93 @@ void QuickHull::createConvexHalfEdgeMesh() { // marked as disabled. Their data slots will be reused. The faces will be // disabled as well, but we need to remember the points that were on the // positive side of them - therefore we save pointers to them. - m_newFaceIndices.clear(); - m_newHalfEdgeIndices.clear(); - m_disabledFacePointVectors.clear(); + newFaceIndices.clear(); + newHalfEdgeIndices.clear(); + disabledFacePointVectors.clear(); size_t disableCounter = 0; - for (auto faceIndex : m_visibleFaces) { - auto& disabledFace = m_mesh.m_faces[faceIndex]; - auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace); + for (auto faceIndex : visibleFaces) { + auto& disabledFace = mesh.faces[faceIndex]; + auto halfEdges = mesh.getHalfEdgeIndicesOfFace(disabledFace); for (size_t j = 0; j < 3; j++) { - if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1 << j)) == 0) { + if ((disabledFace.horizonEdgesOnCurrentIteration & (1 << j)) == 0) { if (disableCounter < horizonEdgeCount * 2) { // Use on this iteration - m_newHalfEdgeIndices.push_back(halfEdges[j]); + newHalfEdgeIndices.push_back(halfEdges[j]); disableCounter++; } else { // Mark for reusal on later iteration step - m_mesh.disableHalfEdge(halfEdges[j]); + mesh.disableHalfEdge(halfEdges[j]); } } } // Disable the face, but retain pointer to the points that were on the // positive side of it. We need to assign those points to the new faces we // create shortly. - auto t = m_mesh.disableFace(faceIndex); + auto t = mesh.disableFace(faceIndex); if (t) { - assert(t->size()); // Because we should not assign point vectors to - // faces unless needed... - m_disabledFacePointVectors.push_back(std::move(t)); + // Because we should not assign point vectors to faces unless needed... + assert(t->size()); + disabledFacePointVectors.push_back(std::move(t)); } } if (disableCounter < horizonEdgeCount * 2) { const size_t newHalfEdgesNeeded = horizonEdgeCount * 2 - disableCounter; for (size_t i = 0; i < newHalfEdgesNeeded; i++) { - m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge()); + newHalfEdgeIndices.push_back(mesh.addHalfEdge()); } } // Create new faces using the edgeloop for (size_t i = 0; i < horizonEdgeCount; i++) { - const size_t AB = m_horizonEdges[i]; + const size_t AB = horizonEdges[i]; auto horizonEdgeVertexIndices = - m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]); + mesh.getVertexIndicesOfHalfEdge(mesh.halfEdges[AB]); size_t A, B, C; A = horizonEdgeVertexIndices[0]; B = horizonEdgeVertexIndices[1]; C = activePointIndex; - const size_t newFaceIndex = m_mesh.addFace(); - m_newFaceIndices.push_back(newFaceIndex); + const size_t newFaceIndex = mesh.addFace(); + newFaceIndices.push_back(newFaceIndex); - const size_t CA = m_newHalfEdgeIndices[2 * i + 0]; - const size_t BC = m_newHalfEdgeIndices[2 * i + 1]; + const size_t CA = newHalfEdgeIndices[2 * i + 0]; + const size_t BC = newHalfEdgeIndices[2 * i + 1]; - m_mesh.m_halfEdges[AB].m_next = BC; - m_mesh.m_halfEdges[BC].m_next = CA; - m_mesh.m_halfEdges[CA].m_next = AB; + mesh.halfEdges[AB].next = BC; + mesh.halfEdges[BC].next = CA; + mesh.halfEdges[CA].next = AB; - m_mesh.m_halfEdges[BC].m_face = newFaceIndex; - m_mesh.m_halfEdges[CA].m_face = newFaceIndex; - m_mesh.m_halfEdges[AB].m_face = newFaceIndex; + mesh.halfEdges[BC].face = newFaceIndex; + mesh.halfEdges[CA].face = newFaceIndex; + mesh.halfEdges[AB].face = newFaceIndex; - m_mesh.m_halfEdges[CA].m_endVertex = A; - m_mesh.m_halfEdges[BC].m_endVertex = C; + mesh.halfEdges[CA].endVertex = A; + mesh.halfEdges[BC].endVertex = C; - auto& newFace = m_mesh.m_faces[newFaceIndex]; + auto& newFace = mesh.faces[newFaceIndex]; const glm::dvec3 planeNormal = mathutils::getTriangleNormal( - m_vertexData[A], m_vertexData[B], activePoint); - newFace.m_P = Plane(planeNormal, activePoint); - newFace.m_he = AB; - - m_mesh.m_halfEdges[CA].m_opp = - m_newHalfEdgeIndices[i > 0 ? i * 2 - 1 : 2 * horizonEdgeCount - 1]; - m_mesh.m_halfEdges[BC].m_opp = - m_newHalfEdgeIndices[((i + 1) * 2) % (horizonEdgeCount * 2)]; + originalVertexData[A], originalVertexData[B], activePoint); + newFace.P = Plane(planeNormal, activePoint); + newFace.he = AB; + + mesh.halfEdges[CA].opp = + newHalfEdgeIndices[i > 0 ? i * 2 - 1 : 2 * horizonEdgeCount - 1]; + mesh.halfEdges[BC].opp = + newHalfEdgeIndices[((i + 1) * 2) % (horizonEdgeCount * 2)]; } // Assign points that were on the positive side of the disabled faces to the // new faces. - for (auto& disabledPoints : m_disabledFacePointVectors) { + for (auto& disabledPoints : disabledFacePointVectors) { assert(disabledPoints); for (const auto& point : *(disabledPoints)) { if (point == activePointIndex) { continue; } for (size_t j = 0; j < horizonEdgeCount; j++) { - if (addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]], point)) { + if (addPointToFace(mesh.faces[newFaceIndices[j]], point)) { break; } } @@ -312,20 +306,20 @@ void QuickHull::createConvexHalfEdgeMesh() { } // Increase face stack size if needed - for (const auto newFaceIndex : m_newFaceIndices) { - auto& newFace = m_mesh.m_faces[newFaceIndex]; - if (newFace.m_pointsOnPositiveSide) { - assert(newFace.m_pointsOnPositiveSide->size() > 0); - if (!newFace.m_inFaceStack) { - m_faceList.push_back(newFaceIndex); - newFace.m_inFaceStack = 1; + for (const auto newFaceIndex : newFaceIndices) { + auto& newFace = mesh.faces[newFaceIndex]; + if (newFace.pointsOnPositiveSide) { + assert(newFace.pointsOnPositiveSide->size() > 0); + if (!newFace.inFaceStack) { + faceList.push_back(newFaceIndex); + newFace.inFaceStack = 1; } } } } // Cleanup - m_indexVectorPool.clear(); + indexVectorPool.clear(); } /* @@ -334,12 +328,12 @@ void QuickHull::createConvexHalfEdgeMesh() { std::array QuickHull::getExtremeValues() { std::array outIndices{0, 0, 0, 0, 0, 0}; - double extremeVals[6] = {m_vertexData[0].x, m_vertexData[0].x, - m_vertexData[0].y, m_vertexData[0].y, - m_vertexData[0].z, m_vertexData[0].z}; - const size_t vCount = m_vertexData.size(); + double extremeVals[6] = {originalVertexData[0].x, originalVertexData[0].x, + originalVertexData[0].y, originalVertexData[0].y, + originalVertexData[0].z, originalVertexData[0].z}; + const size_t vCount = originalVertexData.size(); for (size_t i = 1; i < vCount; i++) { - const glm::dvec3& pos = m_vertexData[i]; + const glm::dvec3& pos = originalVertexData[i]; if (pos.x > extremeVals[0]) { extremeVals[0] = pos.x; outIndices[0] = i; @@ -368,12 +362,11 @@ std::array QuickHull::getExtremeValues() { bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { const size_t horizonEdgeCount = horizonEdges.size(); for (size_t i = 0; i < horizonEdgeCount - 1; i++) { - const size_t endVertex = m_mesh.m_halfEdges[horizonEdges[i]].m_endVertex; + const size_t endVertex = mesh.halfEdges[horizonEdges[i]].endVertex; bool foundNext = false; for (size_t j = i + 1; j < horizonEdgeCount; j++) { const size_t beginVertex = - m_mesh.m_halfEdges[m_mesh.m_halfEdges[horizonEdges[j]].m_opp] - .m_endVertex; + mesh.halfEdges[mesh.halfEdges[horizonEdges[j]].opp].endVertex; if (beginVertex == endVertex) { std::swap(horizonEdges[i + 1], horizonEdges[j]); foundNext = true; @@ -384,17 +377,15 @@ bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { return false; } } - assert( - m_mesh.m_halfEdges[horizonEdges[horizonEdges.size() - 1]].m_endVertex == - m_mesh.m_halfEdges[m_mesh.m_halfEdges[horizonEdges[0]].m_opp] - .m_endVertex); + assert(mesh.halfEdges[horizonEdges[horizonEdges.size() - 1]].endVertex == + mesh.halfEdges[mesh.halfEdges[horizonEdges[0]].opp].endVertex); return true; } double QuickHull::getScale(const std::array& extremeValues) { double s = 0; for (size_t i = 0; i < 6; i++) { - const double* v = (const double*)(&m_vertexData[extremeValues[i]]); + const double* v = (const double*)(&originalVertexData[extremeValues[i]]); v += i / 2; auto a = std::abs(*v); if (a > s) { @@ -405,157 +396,162 @@ double QuickHull::getScale(const std::array& extremeValues) { } void QuickHull::setupInitialTetrahedron() { - const size_t vertexCount = m_vertexData.size(); + const size_t vertexCount = originalVertexData.size(); // If we have at most 4 points, just return a degenerate tetrahedron: if (vertexCount <= 4) { size_t v[4] = {0, std::min((size_t)1, vertexCount - 1), std::min((size_t)2, vertexCount - 1), std::min((size_t)3, vertexCount - 1)}; - const glm::dvec3 N = mathutils::getTriangleNormal( - m_vertexData[v[0]], m_vertexData[v[1]], m_vertexData[v[2]]); - const Plane trianglePlane(N, m_vertexData[v[0]]); - if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) { + const glm::dvec3 N = mathutils::getTriangleNormal(originalVertexData[v[0]], + originalVertexData[v[1]], + originalVertexData[v[2]]); + const Plane trianglePlane(N, originalVertexData[v[0]]); + if (trianglePlane.isPointOnPositiveSide(originalVertexData[v[3]])) { std::swap(v[0], v[1]); } - return m_mesh.setup(v[0], v[1], v[2], v[3]); + return mesh.setup(v[0], v[1], v[2], v[3]); } // Find two most distant extreme points. - double maxD = m_epsilonSquared; + double maxD = epsilonSquared; std::pair selectedPoints; for (size_t i = 0; i < 6; i++) { for (size_t j = i + 1; j < 6; j++) { // I found a function for squaredDistance but i can't seem to include it // like this for some reason - const double d = mathutils::getSquaredDistance( - m_vertexData[m_extremeValues[i]], m_vertexData[m_extremeValues[j]]); + const double d = + mathutils::getSquaredDistance(originalVertexData[extremeValues[i]], + originalVertexData[extremeValues[j]]); if (d > maxD) { maxD = d; - selectedPoints = {m_extremeValues[i], m_extremeValues[j]}; + selectedPoints = {extremeValues[i], extremeValues[j]}; } } } - if (maxD == m_epsilonSquared) { + if (maxD == epsilonSquared) { // A degenerate case: the point cloud seems to consists of a single point - return m_mesh.setup(0, std::min((size_t)1, vertexCount - 1), - std::min((size_t)2, vertexCount - 1), - std::min((size_t)3, vertexCount - 1)); + return mesh.setup(0, std::min((size_t)1, vertexCount - 1), + std::min((size_t)2, vertexCount - 1), + std::min((size_t)3, vertexCount - 1)); } assert(selectedPoints.first != selectedPoints.second); // Find the most distant point to the line between the two chosen extreme // points. - const Ray r(m_vertexData[selectedPoints.first], - (m_vertexData[selectedPoints.second] - - m_vertexData[selectedPoints.first])); - maxD = m_epsilonSquared; + const Ray r(originalVertexData[selectedPoints.first], + (originalVertexData[selectedPoints.second] - + originalVertexData[selectedPoints.first])); + maxD = epsilonSquared; size_t maxI = std::numeric_limits::max(); - const size_t vCount = m_vertexData.size(); + const size_t vCount = originalVertexData.size(); for (size_t i = 0; i < vCount; i++) { - const double distToRay = - mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i], r); + const double distToRay = mathutils::getSquaredDistanceBetweenPointAndRay( + originalVertexData[i], r); if (distToRay > maxD) { maxD = distToRay; maxI = i; } } - if (maxD == m_epsilonSquared) { + if (maxD == epsilonSquared) { // It appears that the point cloud belongs to a 1 dimensional subspace of // R^3: convex hull has no volume => return a thin triangle Pick any point // other than selectedPoints.first and selectedPoints.second as the third // point of the triangle - auto it = std::find_if(m_vertexData.begin(), m_vertexData.end(), - [&](const vec3& ve) { - return ve != m_vertexData[selectedPoints.first] && - ve != m_vertexData[selectedPoints.second]; - }); - const size_t thirdPoint = (it == m_vertexData.end()) - ? selectedPoints.first - : std::distance(m_vertexData.begin(), it); - it = std::find_if(m_vertexData.begin(), m_vertexData.end(), - [&](const vec3& ve) { - return ve != m_vertexData[selectedPoints.first] && - ve != m_vertexData[selectedPoints.second] && - ve != m_vertexData[thirdPoint]; - }); - const size_t fourthPoint = (it == m_vertexData.end()) - ? selectedPoints.first - : std::distance(m_vertexData.begin(), it); - return m_mesh.setup(selectedPoints.first, selectedPoints.second, thirdPoint, - fourthPoint); + auto it = + std::find_if(originalVertexData.begin(), originalVertexData.end(), + [&](const vec3& ve) { + return ve != originalVertexData[selectedPoints.first] && + ve != originalVertexData[selectedPoints.second]; + }); + const size_t thirdPoint = + (it == originalVertexData.end()) + ? selectedPoints.first + : std::distance(originalVertexData.begin(), it); + it = + std::find_if(originalVertexData.begin(), originalVertexData.end(), + [&](const vec3& ve) { + return ve != originalVertexData[selectedPoints.first] && + ve != originalVertexData[selectedPoints.second] && + ve != originalVertexData[thirdPoint]; + }); + const size_t fourthPoint = + (it == originalVertexData.end()) + ? selectedPoints.first + : std::distance(originalVertexData.begin(), it); + return mesh.setup(selectedPoints.first, selectedPoints.second, thirdPoint, + fourthPoint); } // These three points form the base triangle for our tetrahedron. assert(selectedPoints.first != maxI && selectedPoints.second != maxI); std::array baseTriangle{selectedPoints.first, selectedPoints.second, maxI}; - const glm::dvec3 baseTriangleVertices[] = {m_vertexData[baseTriangle[0]], - m_vertexData[baseTriangle[1]], - m_vertexData[baseTriangle[2]]}; + const glm::dvec3 baseTriangleVertices[] = { + originalVertexData[baseTriangle[0]], originalVertexData[baseTriangle[1]], + originalVertexData[baseTriangle[2]]}; // Next step is to find the 4th vertex of the tetrahedron. We naturally choose // the point farthest away from the triangle plane. - maxD = m_epsilon; + maxD = epsilon; maxI = 0; const glm::dvec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], baseTriangleVertices[1], baseTriangleVertices[2]); Plane trianglePlane(N, baseTriangleVertices[0]); for (size_t i = 0; i < vCount; i++) { - const double d = std::abs( - mathutils::getSignedDistanceToPlane(m_vertexData[i], trianglePlane)); + const double d = std::abs(mathutils::getSignedDistanceToPlane( + originalVertexData[i], trianglePlane)); if (d > maxD) { maxD = d; maxI = i; } } - if (maxD == m_epsilon) { + if (maxD == epsilon) { // All the points seem to lie on a 2D subspace of R^3. How to handle this? // Well, let's add one extra point to the point cloud so that the convex // hull will have volume. - m_planar = true; + planar = true; const vec3 N1 = mathutils::getTriangleNormal(baseTriangleVertices[1], baseTriangleVertices[2], baseTriangleVertices[0]); - m_planarPointCloudTemp.clear(); - m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(), - m_vertexData.begin(), m_vertexData.end()); - const vec3 extraPoint = N1 + m_vertexData[0]; - m_planarPointCloudTemp.push_back(extraPoint); - maxI = m_planarPointCloudTemp.size() - 1; - m_vertexData = VertexDataSource(m_planarPointCloudTemp); + planarPointCloudTemp.clear(); + planarPointCloudTemp.insert(planarPointCloudTemp.begin(), + originalVertexData.begin(), + originalVertexData.end()); + const vec3 extraPoint = N1 + originalVertexData[0]; + planarPointCloudTemp.push_back(extraPoint); + maxI = planarPointCloudTemp.size() - 1; + originalVertexData = VertexDataSource(planarPointCloudTemp); } // Enforce CCW orientation (if user prefers clockwise orientation, swap two // vertices in each triangle when final mesh is created) const Plane triPlane(N, baseTriangleVertices[0]); - if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) { + if (triPlane.isPointOnPositiveSide(originalVertexData[maxI])) { std::swap(baseTriangle[0], baseTriangle[1]); } // Create a tetrahedron half edge mesh and compute planes defined by each // triangle - m_mesh.setup(baseTriangle[0], baseTriangle[1], baseTriangle[2], maxI); - for (auto& f : m_mesh.m_faces) { - auto v = m_mesh.getVertexIndicesOfFace(f); - const glm::dvec3& va = m_vertexData[v[0]]; - const glm::dvec3& vb = m_vertexData[v[1]]; - const glm::dvec3& vc = m_vertexData[v[2]]; + 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 = mathutils::getTriangleNormal(va, vb, vc); const Plane plane(N1, va); - f.m_P = plane; + f.P = plane; } // Finally we assign a face for each vertex outside the tetrahedron (vertices // inside the tetrahedron have no role anymore) for (size_t i = 0; i < vCount; i++) { - for (auto& face : m_mesh.m_faces) { + for (auto& face : mesh.faces) { if (addPointToFace(face, i)) { break; } } } } - -} // namespace quickhull diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h new file mode 100644 index 000000000..28254421a --- /dev/null +++ b/src/manifold/src/quickhull.h @@ -0,0 +1,824 @@ +#ifndef QUICKHULL_HPP_ +#define QUICKHULL_HPP_ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "par.h" + +// Pool.hpp + +class Pool { + std::vector>> data; + + public: + void clear() { data.clear(); } + + void reclaim(std::unique_ptr>& ptr) { + data.push_back(std::move(ptr)); + } + + std::unique_ptr> get() { + if (data.size() == 0) { + return std::unique_ptr>(new std::vector()); + } + auto it = data.end() - 1; + std::unique_ptr> r = std::move(*it); + data.erase(it); + return r; + } +}; + +// Plane.hpp + +class Plane { + public: + glm::dvec3 N; + + // Signed distance (if normal is of length 1) to the plane from origin + double D; + + // Normal length squared + double sqrNLength; + + bool isPointOnPositiveSide(const glm::dvec3& Q) const { + double d = glm::dot(N, Q) + D; + if (d >= 0) return true; + return false; + } + + Plane() = default; + + // Construct a plane using normal N and any point P on the plane + Plane(const glm::dvec3& N, const glm::dvec3& P) + : N(N), + D(glm::dot(-N, P)), + sqrNLength(N.x * N.x + N.y * N.y + N.z * N.z) {} +}; + +// Ray.hpp + +struct Ray { + const glm::dvec3 S; + const glm::dvec3 V; + const double VInvLengthSquared; + + Ray(const glm::dvec3& S, const glm::dvec3& V) + : S(S), + V(V), + VInvLengthSquared(1 / (V.x * V.x + V.y * V.y + V.z * V.z)) {} +}; + +// VertexDataSource + +class VertexDataSource { + const glm::dvec3* ptr; + size_t count; + + public: + VertexDataSource(const glm::dvec3* ptr, size_t count) + : ptr(ptr), count(count) {} + + VertexDataSource(const std::vector& vec) + : ptr(&vec[0]), count(vec.size()) {} + + VertexDataSource() : ptr(nullptr), count(0) {} + + VertexDataSource& operator=(const VertexDataSource& other) = default; + + size_t size() const { return count; } + + const glm::dvec3& operator[](size_t index) const { return ptr[index]; } + + const glm::dvec3* begin() const { return ptr; } + + const glm::dvec3* end() const { return ptr + count; } +}; + +// Mesh.hpp + +class MeshBuilder { + public: + struct HalfEdge { + size_t endVertex; + size_t opp; + size_t face; + size_t next; + + void disable() { endVertex = std::numeric_limits::max(); } + + bool isDisabled() const { + return endVertex == std::numeric_limits::max(); + } + }; + + struct Face { + size_t he; + Plane P{}; + double mostDistantPointDist; + size_t mostDistantPoint; + size_t visibilityCheckedOnIteration; + std::uint8_t isVisibleFaceOnCurrentIteration : 1; + std::uint8_t inFaceStack : 1; + std::uint8_t horizonEdgesOnCurrentIteration : 3; + // Bit for each half edge assigned to this face, each being 0 or 1 depending + // on whether the edge belongs to horizon edge + std::unique_ptr> pointsOnPositiveSide; + + Face() + : he(std::numeric_limits::max()), + mostDistantPointDist(0), + mostDistantPoint(0), + visibilityCheckedOnIteration(0), + isVisibleFaceOnCurrentIteration(0), + inFaceStack(0), + horizonEdgesOnCurrentIteration(0) {} + + void disable() { he = std::numeric_limits::max(); } + + bool isDisabled() const { return he == std::numeric_limits::max(); } + }; + + // Mesh data + std::vector faces; + std::vector halfEdges; + + // When the mesh is modified and faces and half edges are removed from it, we + // do not actually remove them from the container vectors. Insted, they are + // marked as disabled which means that the indices can be reused when we need + // to add new faces and half edges to the mesh. We store the free indices in + // the following vectors. + std::vector disabledFaces, disabledHalfEdges; + + size_t addFace() { + if (disabledFaces.size()) { + size_t index = disabledFaces.back(); + auto& f = faces[index]; + assert(f.isDisabled()); + assert(!f.pointsOnPositiveSide); + f.mostDistantPointDist = 0; + disabledFaces.pop_back(); + return index; + } + faces.emplace_back(); + return faces.size() - 1; + } + + size_t addHalfEdge() { + if (disabledHalfEdges.size()) { + const size_t index = disabledHalfEdges.back(); + disabledHalfEdges.pop_back(); + return index; + } + halfEdges.emplace_back(); + return halfEdges.size() - 1; + } + + // Mark a face as disabled and return a pointer to the points that were on the + // positive of it. + std::unique_ptr> disableFace(size_t faceIndex) { + auto& f = faces[faceIndex]; + f.disable(); + disabledFaces.push_back(faceIndex); + return std::move(f.pointsOnPositiveSide); + } + + void disableHalfEdge(size_t heIndex) { + auto& he = halfEdges[heIndex]; + he.disable(); + disabledHalfEdges.push_back(heIndex); + } + + MeshBuilder() = default; + + // Create a mesh with initial tetrahedron ABCD. Dot product of AB with the + // normal of triangle ABC should be negative. + void setup(size_t a, size_t b, size_t c, size_t d) { + faces.clear(); + halfEdges.clear(); + disabledFaces.clear(); + disabledHalfEdges.clear(); + + faces.reserve(4); + halfEdges.reserve(12); + + // Create halfedges + HalfEdge AB; + AB.endVertex = b; + AB.opp = 6; + AB.face = 0; + AB.next = 1; + halfEdges.push_back(AB); + + HalfEdge BC; + BC.endVertex = c; + BC.opp = 9; + BC.face = 0; + BC.next = 2; + halfEdges.push_back(BC); + + HalfEdge CA; + CA.endVertex = a; + CA.opp = 3; + CA.face = 0; + CA.next = 0; + halfEdges.push_back(CA); + + HalfEdge AC; + AC.endVertex = c; + AC.opp = 2; + AC.face = 1; + AC.next = 4; + halfEdges.push_back(AC); + + HalfEdge CD; + CD.endVertex = d; + CD.opp = 11; + CD.face = 1; + CD.next = 5; + halfEdges.push_back(CD); + + HalfEdge DA; + DA.endVertex = a; + DA.opp = 7; + DA.face = 1; + DA.next = 3; + halfEdges.push_back(DA); + + HalfEdge BA; + BA.endVertex = a; + BA.opp = 0; + BA.face = 2; + BA.next = 7; + halfEdges.push_back(BA); + + HalfEdge AD; + AD.endVertex = d; + AD.opp = 5; + AD.face = 2; + AD.next = 8; + halfEdges.push_back(AD); + + HalfEdge DB; + DB.endVertex = b; + DB.opp = 10; + DB.face = 2; + DB.next = 6; + halfEdges.push_back(DB); + + HalfEdge CB; + CB.endVertex = b; + CB.opp = 1; + CB.face = 3; + CB.next = 10; + halfEdges.push_back(CB); + + HalfEdge BD; + BD.endVertex = d; + BD.opp = 8; + BD.face = 3; + BD.next = 11; + halfEdges.push_back(BD); + + HalfEdge DC; + DC.endVertex = c; + DC.opp = 4; + DC.face = 3; + DC.next = 9; + halfEdges.push_back(DC); + + // Create faces + Face ABC; + ABC.he = 0; + faces.push_back(std::move(ABC)); + + Face ACD; + ACD.he = 3; + faces.push_back(std::move(ACD)); + + Face BAD; + BAD.he = 6; + faces.push_back(std::move(BAD)); + + Face CBD; + CBD.he = 9; + faces.push_back(std::move(CBD)); + } + + std::array getVertexIndicesOfFace(const Face& f) const { + std::array v; + const HalfEdge* he = &halfEdges[f.he]; + v[0] = he->endVertex; + he = &halfEdges[he->next]; + v[1] = he->endVertex; + he = &halfEdges[he->next]; + v[2] = he->endVertex; + return v; + } + + std::array getVertexIndicesOfHalfEdge(const HalfEdge& he) const { + return {halfEdges[he.opp].endVertex, he.endVertex}; + } + + std::array getHalfEdgeIndicesOfFace(const Face& f) const { + return {f.he, halfEdges[f.he].next, halfEdges[halfEdges[f.he].next].next}; + } +}; + +// ConvexHull.hpp + +class ConvexHull { + std::unique_ptr> optimizedVertexBuffer; + VertexDataSource vertices; + std::vector indices; + + public: + ConvexHull() {} + + // Copy constructor + ConvexHull(const ConvexHull& o) { + indices = o.indices; + if (o.optimizedVertexBuffer) { + optimizedVertexBuffer.reset( + new std::vector(*o.optimizedVertexBuffer)); + vertices = VertexDataSource(*optimizedVertexBuffer); + } else { + vertices = o.vertices; + } + } + + ConvexHull& operator=(const ConvexHull& o) { + if (&o == this) { + return *this; + } + indices = o.indices; + if (o.optimizedVertexBuffer) { + optimizedVertexBuffer.reset( + new std::vector(*o.optimizedVertexBuffer)); + vertices = VertexDataSource(*optimizedVertexBuffer); + } else { + vertices = o.vertices; + } + return *this; + } + + ConvexHull(ConvexHull&& o) { + indices = std::move(o.indices); + if (o.optimizedVertexBuffer) { + optimizedVertexBuffer = std::move(o.optimizedVertexBuffer); + o.vertices = VertexDataSource(); + vertices = VertexDataSource(*optimizedVertexBuffer); + } else { + vertices = o.vertices; + } + } + + ConvexHull& operator=(ConvexHull&& o) { + if (&o == this) { + return *this; + } + indices = std::move(o.indices); + if (o.optimizedVertexBuffer) { + optimizedVertexBuffer = std::move(o.optimizedVertexBuffer); + o.vertices = VertexDataSource(); + vertices = VertexDataSource(*optimizedVertexBuffer); + } else { + vertices = o.vertices; + } + return *this; + } + + // Construct vertex and index buffers from half edge mesh and pointcloud + ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, + bool CCW, bool useOriginalIndices) { + if (!useOriginalIndices) { + optimizedVertexBuffer.reset(new std::vector()); + } + + std::vector faceProcessed(mesh.faces.size(), false); + std::vector faceStack; + // Map vertex indices from original point cloud to the new mesh vertex + // indices + std::unordered_map vertexIndexMapping; + for (size_t i = 0; i < mesh.faces.size(); i++) { + if (!mesh.faces[i].isDisabled()) { + faceStack.push_back(i); + break; + } + } + if (faceStack.size() == 0) { + return; + } + + const size_t iCCW = CCW ? 1 : 0; + const size_t finalMeshFaceCount = + mesh.faces.size() - mesh.disabledFaces.size(); + indices.reserve(finalMeshFaceCount * 3); + + while (faceStack.size()) { + auto it = faceStack.end() - 1; + size_t top = *it; + assert(!mesh.faces[top].isDisabled()); + faceStack.erase(it); + if (faceProcessed[top]) { + continue; + } else { + faceProcessed[top] = true; + auto halfEdges = mesh.getHalfEdgeIndicesOfFace(mesh.faces[top]); + size_t adjacent[] = { + mesh.halfEdges[mesh.halfEdges[halfEdges[0]].opp].face, + mesh.halfEdges[mesh.halfEdges[halfEdges[1]].opp].face, + mesh.halfEdges[mesh.halfEdges[halfEdges[2]].opp].face}; + for (auto a : adjacent) { + if (!faceProcessed[a] && !mesh.faces[a].isDisabled()) { + faceStack.push_back(a); + } + } + auto vertices = mesh.getVertexIndicesOfFace(mesh.faces[top]); + if (!useOriginalIndices) { + for (auto& v : vertices) { + auto itV = vertexIndexMapping.find(v); + if (itV == vertexIndexMapping.end()) { + optimizedVertexBuffer->push_back(pointCloud[v]); + vertexIndexMapping[v] = optimizedVertexBuffer->size() - 1; + v = optimizedVertexBuffer->size() - 1; + } else { + v = itV->second; + } + } + } + indices.push_back(vertices[0]); + indices.push_back(vertices[1 + iCCW]); + indices.push_back(vertices[2 - iCCW]); + } + } + + if (!useOriginalIndices) { + vertices = VertexDataSource(*optimizedVertexBuffer); + } else { + vertices = pointCloud; + } + } + + std::vector& getIndexBuffer() { return indices; } + + const std::vector& getIndexBuffer() const { return indices; } + + VertexDataSource& getVertexBuffer() { return vertices; } + + const VertexDataSource& getVertexBuffer() const { return vertices; } +}; + +// HalfEdgeMesh.hpp + +class HalfEdgeMesh { + public: + struct HalfEdge { + size_t endVertex; + size_t opp; + size_t face; + size_t next; + }; + + struct Face { + // Index of one of the half edges of this face + size_t halfEdgeIndex; + }; + + std::vector vertices; + std::vector faces; + std::vector halfEdges; + + HalfEdgeMesh(const MeshBuilder& builderObject, + const VertexDataSource& vertexData) { + std::unordered_map faceMapping; + std::unordered_map halfEdgeMapping; + std::unordered_map vertexMapping; + + size_t i = 0; + for (const auto& face : builderObject.faces) { + if (!face.isDisabled()) { + faces.push_back({static_cast(face.he)}); + faceMapping[i] = faces.size() - 1; + + const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face); + for (const auto heIndex : heIndices) { + const size_t vertexIndex = builderObject.halfEdges[heIndex].endVertex; + if (vertexMapping.count(vertexIndex) == 0) { + vertices.push_back(vertexData[vertexIndex]); + vertexMapping[vertexIndex] = vertices.size() - 1; + } + } + } + i++; + } + + i = 0; + for (const auto& halfEdge : builderObject.halfEdges) { + if (!halfEdge.isDisabled()) { + halfEdges.push_back({static_cast(halfEdge.endVertex), + static_cast(halfEdge.opp), + static_cast(halfEdge.face), + static_cast(halfEdge.next)}); + halfEdgeMapping[i] = halfEdges.size() - 1; + } + i++; + } + + for (auto& face : faces) { + assert(halfEdgeMapping.count(face.halfEdgeIndex) == 1); + face.halfEdgeIndex = halfEdgeMapping[face.halfEdgeIndex]; + } + + for (auto& he : halfEdges) { + he.face = faceMapping[he.face]; + he.opp = halfEdgeMapping[he.opp]; + he.next = halfEdgeMapping[he.next]; + he.endVertex = vertexMapping[he.endVertex]; + } + } +}; + +// MathUtils.hpp + +namespace mathutils { + +inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p, + const Ray& r) { + const glm::dvec3 s = p - r.S; + double t = glm::dot(s, r.V); + return (s.x * s.x + s.y * s.y + s.z * s.z) - t * t * r.VInvLengthSquared; +} + +inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { + return ((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + + (p1.z - p2.z) * (p1.z - p2.z)); +} +// Note that the unit of distance returned is relative to plane's normal's +// length (divide by N.getNormalized() if needed to get the "real" distance). +inline double getSignedDistanceToPlane(const glm::dvec3& v, const Plane& p) { + return glm::dot(p.N, v) + p.D; +} + +inline glm::dvec3 getTriangleNormal(const glm::dvec3& a, const glm::dvec3& b, + const glm::dvec3& c) { + // We want to get (a-c).crossProduct(b-c) without constructing temp vectors + double x = a.x - c.x; + double y = a.y - c.y; + double z = a.z - c.z; + double rhsx = b.x - c.x; + double rhsy = b.y - c.y; + double rhsz = b.z - c.z; + double px = y * rhsz - z * rhsy; + double py = z * rhsx - x * rhsz; + double pz = x * rhsy - y * rhsx; + return glm::dvec3(px, py, pz); +} + +} // namespace mathutils + +// QuickHull.hpp + +/* + * Implementation of the 3D QuickHull algorithm by Antti Kuukka + * + * No copyrights. What follows is 100% Public Domain. + * + * + * + * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) + * + * OUTPUT: a ConvexHull object which provides vertex and index buffers of the + *generated convex hull as a triangle mesh. + * + * + * + * The implementation is thread-safe if each thread is using its own QuickHull + *object. + * + * + * SUMMARY OF THE ALGORITHM: + * - Create initial simplex (tetrahedron) using extreme points. We have + *four faces now and they form a convex mesh M. + * - For each point, assign them to the first face for which they are on + *the positive side of (so each point is assigned to at most one face). Points + *inside the initial tetrahedron are left behind now and no longer affect the + *calculations. + * - Add all faces that have points assigned to them to Face Stack. + * - Iterate until Face Stack is empty: + * - Pop topmost face F from the stack + * - From the points assigned to F, pick the point P that is + *farthest away from the plane defined by F. + * - Find all faces of M that have P on their positive side. Let us + *call these the "visible faces". + * - Because of the way M is constructed, these faces are + *connected. Solve their horizon edge loop. + * - "Extrude to P": Create new faces by connecting + *P with the points belonging to the horizon edge. Add the new faces to M and + *remove the visible faces from M. + * - Each point that was assigned to visible faces is now assigned + *to at most one of the newly created faces. + * - Those new faces that have points assigned to them are added to + *the top of Face Stack. + * - M is now the convex hull. + * + * TO DO: + * - Implement a proper 2D QuickHull and use that to solve the degenerate 2D + *case (when all the points lie on the same plane in 3D space). + * */ + +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 { + using vec3 = glm::dvec3; + + double epsilon, epsilonSquared, scale; + bool planar; + std::vector planarPointCloudTemp; + VertexDataSource originalVertexData; + MeshBuilder mesh; + std::array extremeValues; + DiagnosticsData diagnostics; + + // Temporary variables used during iteration process + std::vector newFaceIndices; + std::vector newHalfEdgeIndices; + std::vector>> disabledFacePointVectors; + std::vector visibleFaces; + std::vector horizonEdges; + struct FaceData { + size_t faceIndex; + // If the face turns out not to be visible, this half edge will be marked as + // horizon edge + size_t enteredFromHalfEdge; + FaceData(size_t fi, size_t he) : faceIndex(fi), enteredFromHalfEdge(he) {} + }; + std::vector possiblyVisibleFaces; + std::deque faceList; + + // Create a half edge mesh representing the base tetrahedron from which the + // QuickHull iteration proceeds. extremeValues must be properly set up when + // this is called. + void setupInitialTetrahedron(); + + // Given a list of half edges, try to rearrange them so that they form a loop. + // Return true on success. + bool reorderHorizonEdges(std::vector& horizonEdges); + + // Find indices of extreme values (max x, min x, max y, min y, max z, min z) + // for the given point cloud + std::array getExtremeValues(); + + // Compute scale of the vertex data. + double getScale(const std::array& extremeValues); + + // Each face contains a unique pointer to a vector of indices. However, many - + // often most - faces do not have any points on the positive side of them + // especially at the the end of the iteration. When a face is removed from the + // mesh, its associated point vector, if such exists, is moved to the index + // vector pool, and when we need to add new faces with points on the positive + // side to the mesh, we reuse these vectors. This reduces the amount of + // std::vectors we have to deal with, and impact on performance is remarkable. + Pool indexVectorPool; + inline std::unique_ptr> getIndexVectorFromPool(); + inline void reclaimToIndexVectorPool( + std::unique_ptr>& ptr); + + // Associates a point with a face if the point resides on the positive side of + // 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 + void createConvexHalfEdgeMesh(); + + // Constructs the convex hull into a MeshBuilder object which can be converted + // to a ConvexHull or Mesh object + void buildMesh(const VertexDataSource& pointCloud, bool CCW, + bool useOriginalIndices, double eps); + + // The public getConvexHull functions will setup a VertexDataSource object and + // call this + ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, + bool useOriginalIndices, double eps); + + public: + // Computes convex hull for a given point cloud. + // Params: + // pointCloud: a vector of of 3D points + // CCW: whether the output mesh triangles should have CCW orientation + // useOriginalIndices: should the output mesh use same vertex indices as the + // original point cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices + // that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive of + // it (for a point cloud with scale 1) + ConvexHull getConvexHull(const std::vector& pointCloud, bool CCW, + bool useOriginalIndices, double eps = defaultEps()); + + // Computes convex hull for a given point cloud. + // Params: + // vertexData: pointer to the first 3D point of the point cloud + // vertexCount: number of vertices in the point cloud + // CCW: whether the output mesh triangles should have CCW orientation + // useOriginalIndices: should the output mesh use same vertex indices as the + // original point cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices + // that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive + // side of it (for a point cloud with scale 1) + ConvexHull getConvexHull(const glm::dvec3* vertexData, size_t vertexCount, + bool CCW, bool useOriginalIndices, + double eps = defaultEps()); + + // 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 useOriginalIndices: + // should the output mesh use same vertex indices as the original point + // cloud. If this is false, + // then we generate a new vertex buffer which contains only the vertices + // that are part of the convex hull. + // eps: minimum distance to a plane to consider a point being on positive + // side of it (for a point cloud with scale 1) + ConvexHull getConvexHull(const double* vertexData, size_t vertexCount, + bool CCW, bool useOriginalIndices, + double eps = defaultEps()); + + // 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 a mesh object with half edge structure. + HalfEdgeMesh getConvexHullAsMesh(const double* vertexData, size_t vertexCount, + bool CCW, double eps = defaultEps()); + + // Get diagnostics about last generated convex hull + const DiagnosticsData& getDiagnostics() { return diagnostics; } +}; + +/* + * Inline function definitions + */ + +std::unique_ptr> QuickHull::getIndexVectorFromPool() { + auto r = indexVectorPool.get(); + r->clear(); + return r; +} + +void QuickHull::reclaimToIndexVectorPool( + std::unique_ptr>& ptr) { + const size_t oldSize = ptr->size(); + if ((oldSize + 1) * 128 < ptr->capacity()) { + // Reduce memory usage! Huge vectors are needed at the beginning of + // iteration when faces have many points on their positive side. Later on, + // smaller vectors will suffice. + ptr.reset(nullptr); + return; + } + indexVectorPool.reclaim(ptr); +} + +bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, + size_t pointIndex) { + const double D = + mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); + if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { + if (!f.pointsOnPositiveSide) { + f.pointsOnPositiveSide = std::move(getIndexVectorFromPool()); + } + f.pointsOnPositiveSide->push_back(pointIndex); + if (D > f.mostDistantPointDist) { + f.mostDistantPointDist = D; + f.mostDistantPoint = pointIndex; + } + return true; + } + return false; +} + +#endif /* QUICKHULL_HPP_ */ diff --git a/src/manifold/src/quickhull.hpp b/src/manifold/src/quickhull.hpp deleted file mode 100644 index 0de13b0e2..000000000 --- a/src/manifold/src/quickhull.hpp +++ /dev/null @@ -1,852 +0,0 @@ -#ifndef QUICKHULL_HPP_ -#define QUICKHULL_HPP_ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "par.h" - -// Pool.hpp - - -namespace quickhull { - - class Pool { - std::vector>> m_data; - public: - void clear() { - m_data.clear(); - } - - void reclaim(std::unique_ptr>& ptr) { - m_data.push_back(std::move(ptr)); - } - - std::unique_ptr> get() { - if (m_data.size()==0) { - return std::unique_ptr>(new std::vector()); - } - auto it = m_data.end()-1; - std::unique_ptr> r = std::move(*it); - m_data.erase(it); - return r; - } - - }; - - -// Vector3.hpp - - - - // // Projection onto another vector - // glm::dvec3 projection(const glm::dvec3& o) const { - // T C = dotProduct(o)/o.getLengthSquared(); - // return o*C; - // } - - - -// Plane.hpp - - class Plane { - public: - glm::dvec3 m_N; - - // Signed distance (if normal is of length 1) to the plane from origin - double m_D; - - // Normal length squared - double m_sqrNLength; - - bool isPointOnPositiveSide(const glm::dvec3& Q) const { - double d = glm::dot(m_N,Q)+m_D; - if (d>=0) return true; - return false; - } - - Plane() = default; - - // Construct a plane using normal N and any point P on the plane - Plane(const glm::dvec3& N, const glm::dvec3& P) : m_N(N), m_D(glm::dot(-N,P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) { - - } - }; - -// Ray.hpp - - struct Ray { - const glm::dvec3 m_S; - const glm::dvec3 m_V; - const double m_VInvLengthSquared; - - Ray(const glm::dvec3& S,const glm::dvec3& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/(m_V.x*m_V.x+m_V.y*m_V.y+m_V.z*m_V.z)) { - } - }; - - -// VertexDataSource - - - class VertexDataSource { - const glm::dvec3* m_ptr; - size_t m_count; - - public: - VertexDataSource(const glm::dvec3* ptr, size_t count) : m_ptr(ptr), m_count(count) { - - } - - VertexDataSource(const std::vector& vec) : m_ptr(&vec[0]), m_count(vec.size()) { - - } - - VertexDataSource() : m_ptr(nullptr), m_count(0) { - - } - - VertexDataSource& operator=(const VertexDataSource& other) = default; - - size_t size() const { - return m_count; - } - - const glm::dvec3& operator[](size_t index) const { - return m_ptr[index]; - } - - const glm::dvec3* begin() const { - return m_ptr; - } - - const glm::dvec3* end() const { - return m_ptr + m_count; - } - }; - - - -// Mesh.hpp - - - class MeshBuilder { - public: - struct HalfEdge { - size_t m_endVertex; - size_t m_opp; - size_t m_face; - size_t m_next; - - void disable() { - m_endVertex = std::numeric_limits::max(); - } - - bool isDisabled() const { - return m_endVertex == std::numeric_limits::max(); - } - }; - - struct Face { - size_t m_he; - Plane m_P{}; - double m_mostDistantPointDist; - size_t m_mostDistantPoint; - size_t m_visibilityCheckedOnIteration; - std::uint8_t m_isVisibleFaceOnCurrentIteration : 1; - std::uint8_t m_inFaceStack : 1; - std::uint8_t m_horizonEdgesOnCurrentIteration : 3; // Bit for each half edge assigned to this face, each being 0 or 1 depending on whether the edge belongs to horizon edge - std::unique_ptr> m_pointsOnPositiveSide; - - Face() : m_he(std::numeric_limits::max()), - m_mostDistantPointDist(0), - m_mostDistantPoint(0), - m_visibilityCheckedOnIteration(0), - m_isVisibleFaceOnCurrentIteration(0), - m_inFaceStack(0), - m_horizonEdgesOnCurrentIteration(0) - { - - } - - void disable() { - m_he = std::numeric_limits::max(); - } - - bool isDisabled() const { - return m_he == std::numeric_limits::max(); - } - }; - - // Mesh data - std::vector m_faces; - std::vector m_halfEdges; - - // When the mesh is modified and faces and half edges are removed from it, we do not actually remove them from the container vectors. - // Insted, they are marked as disabled which means that the indices can be reused when we need to add new faces and half edges to the mesh. - // We store the free indices in the following vectors. - std::vector m_disabledFaces,m_disabledHalfEdges; - - size_t addFace() { - if (m_disabledFaces.size()) { - size_t index = m_disabledFaces.back(); - auto& f = m_faces[index]; - assert(f.isDisabled()); - assert(!f.m_pointsOnPositiveSide); - f.m_mostDistantPointDist = 0; - m_disabledFaces.pop_back(); - return index; - } - m_faces.emplace_back(); - return m_faces.size()-1; - } - - size_t addHalfEdge() { - if (m_disabledHalfEdges.size()) { - const size_t index = m_disabledHalfEdges.back(); - m_disabledHalfEdges.pop_back(); - return index; - } - m_halfEdges.emplace_back(); - return m_halfEdges.size()-1; - } - - // Mark a face as disabled and return a pointer to the points that were on the positive of it. - std::unique_ptr> disableFace(size_t faceIndex) { - auto& f = m_faces[faceIndex]; - f.disable(); - m_disabledFaces.push_back(faceIndex); - return std::move(f.m_pointsOnPositiveSide); - } - - void disableHalfEdge(size_t heIndex) { - auto& he = m_halfEdges[heIndex]; - he.disable(); - m_disabledHalfEdges.push_back(heIndex); - } - - MeshBuilder() = default; - - // Create a mesh with initial tetrahedron ABCD. Dot product of AB with the normal of triangle ABC should be negative. - void setup(size_t a, size_t b, size_t c, size_t d) { - m_faces.clear(); - m_halfEdges.clear(); - m_disabledFaces.clear(); - m_disabledHalfEdges.clear(); - - m_faces.reserve(4); - m_halfEdges.reserve(12); - - // Create halfedges - HalfEdge AB; - AB.m_endVertex = b; - AB.m_opp = 6; - AB.m_face = 0; - AB.m_next = 1; - m_halfEdges.push_back(AB); - - HalfEdge BC; - BC.m_endVertex = c; - BC.m_opp = 9; - BC.m_face = 0; - BC.m_next = 2; - m_halfEdges.push_back(BC); - - HalfEdge CA; - CA.m_endVertex = a; - CA.m_opp = 3; - CA.m_face = 0; - CA.m_next = 0; - m_halfEdges.push_back(CA); - - HalfEdge AC; - AC.m_endVertex = c; - AC.m_opp = 2; - AC.m_face = 1; - AC.m_next = 4; - m_halfEdges.push_back(AC); - - HalfEdge CD; - CD.m_endVertex = d; - CD.m_opp = 11; - CD.m_face = 1; - CD.m_next = 5; - m_halfEdges.push_back(CD); - - HalfEdge DA; - DA.m_endVertex = a; - DA.m_opp = 7; - DA.m_face = 1; - DA.m_next = 3; - m_halfEdges.push_back(DA); - - HalfEdge BA; - BA.m_endVertex = a; - BA.m_opp = 0; - BA.m_face = 2; - BA.m_next = 7; - m_halfEdges.push_back(BA); - - HalfEdge AD; - AD.m_endVertex = d; - AD.m_opp = 5; - AD.m_face = 2; - AD.m_next = 8; - m_halfEdges.push_back(AD); - - HalfEdge DB; - DB.m_endVertex = b; - DB.m_opp = 10; - DB.m_face = 2; - DB.m_next = 6; - m_halfEdges.push_back(DB); - - HalfEdge CB; - CB.m_endVertex = b; - CB.m_opp = 1; - CB.m_face = 3; - CB.m_next = 10; - m_halfEdges.push_back(CB); - - HalfEdge BD; - BD.m_endVertex = d; - BD.m_opp = 8; - BD.m_face = 3; - BD.m_next = 11; - m_halfEdges.push_back(BD); - - HalfEdge DC; - DC.m_endVertex = c; - DC.m_opp = 4; - DC.m_face = 3; - DC.m_next = 9; - m_halfEdges.push_back(DC); - - // Create faces - Face ABC; - ABC.m_he = 0; - m_faces.push_back(std::move(ABC)); - - Face ACD; - ACD.m_he = 3; - m_faces.push_back(std::move(ACD)); - - Face BAD; - BAD.m_he = 6; - m_faces.push_back(std::move(BAD)); - - Face CBD; - CBD.m_he = 9; - m_faces.push_back(std::move(CBD)); - } - - std::array getVertexIndicesOfFace(const Face& f) const { - std::array v; - const HalfEdge* he = &m_halfEdges[f.m_he]; - v[0] = he->m_endVertex; - he = &m_halfEdges[he->m_next]; - v[1] = he->m_endVertex; - he = &m_halfEdges[he->m_next]; - v[2] = he->m_endVertex; - return v; - } - - std::array getVertexIndicesOfHalfEdge(const HalfEdge& he) const { - return {m_halfEdges[he.m_opp].m_endVertex,he.m_endVertex}; - } - - std::array getHalfEdgeIndicesOfFace(const Face& f) const { - return {f.m_he,m_halfEdges[f.m_he].m_next,m_halfEdges[m_halfEdges[f.m_he].m_next].m_next}; - } - }; - - - - - -// ConvexHull.hpp - - class ConvexHull { - std::unique_ptr> m_optimizedVertexBuffer; - VertexDataSource m_vertices; - std::vector m_indices; - public: - ConvexHull() {} - - // Copy constructor - ConvexHull(const ConvexHull& o) { - m_indices = o.m_indices; - if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); - } - else { - m_vertices = o.m_vertices; - } - } - - ConvexHull& operator=(const ConvexHull& o) { - if (&o == this) { - return *this; - } - m_indices = o.m_indices; - if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer.reset(new std::vector(*o.m_optimizedVertexBuffer)); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); - } - else { - m_vertices = o.m_vertices; - } - return *this; - } - - ConvexHull(ConvexHull&& o) { - m_indices = std::move(o.m_indices); - if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); - o.m_vertices = VertexDataSource(); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); - } - else { - m_vertices = o.m_vertices; - } - } - - ConvexHull& operator=(ConvexHull&& o) { - if (&o == this) { - return *this; - } - m_indices = std::move(o.m_indices); - if (o.m_optimizedVertexBuffer) { - m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer); - o.m_vertices = VertexDataSource(); - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); - } - else { - m_vertices = o.m_vertices; - } - return *this; - } - - // Construct vertex and index buffers from half edge mesh and pointcloud - ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { - if (!useOriginalIndices) { - m_optimizedVertexBuffer.reset(new std::vector()); - } - - std::vector faceProcessed(mesh.m_faces.size(),false); - std::vector faceStack; - std::unordered_map vertexIndexMapping; // Map vertex indices from original point cloud to the new mesh vertex indices - for (size_t i = 0;ipush_back(pointCloud[v]); - vertexIndexMapping[v] = m_optimizedVertexBuffer->size()-1; - v = m_optimizedVertexBuffer->size()-1; - } - else { - v = itV->second; - } - } - } - m_indices.push_back(vertices[0]); - m_indices.push_back(vertices[1 + iCCW]); - m_indices.push_back(vertices[2 - iCCW]); - } - } - - if (!useOriginalIndices) { - m_vertices = VertexDataSource(*m_optimizedVertexBuffer); - } - else { - m_vertices = pointCloud; - } - } - - std::vector& getIndexBuffer() { - return m_indices; - } - - const std::vector& getIndexBuffer() const { - return m_indices; - } - - VertexDataSource& getVertexBuffer() { - return m_vertices; - } - - const VertexDataSource& getVertexBuffer() const { - return m_vertices; - } - - // Export the mesh to a Waveform OBJ file - void writeWaveformOBJ(const std::string& filename, const std::string& objectName = "quickhull") const - { - std::ofstream objFile; - objFile.open (filename); - objFile << "o " << objectName << "\n"; - for (const auto& v : getVertexBuffer()) { - objFile << "v " << v.x << " " << v.y << " " << v.z << "\n"; - } - const auto& indBuf = getIndexBuffer(); - size_t triangleCount = indBuf.size()/3; - for (size_t i=0;i m_vertices; - std::vector m_faces; - std::vector m_halfEdges; - - HalfEdgeMesh(const MeshBuilder& builderObject, const VertexDataSource& vertexData ) - { - std::unordered_map faceMapping; - std::unordered_map halfEdgeMapping; - std::unordered_map vertexMapping; - - size_t i=0; - for (const auto& face : builderObject.m_faces) { - if (!face.isDisabled()) { - m_faces.push_back({static_cast(face.m_he)}); - faceMapping[i] = m_faces.size()-1; - - const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face); - for (const auto heIndex : heIndices) { - const size_t vertexIndex = builderObject.m_halfEdges[heIndex].m_endVertex; - if (vertexMapping.count(vertexIndex)==0) { - m_vertices.push_back(vertexData[vertexIndex]); - vertexMapping[vertexIndex] = m_vertices.size()-1; - } - } - } - i++; - } - - i=0; - for (const auto& halfEdge : builderObject.m_halfEdges) { - if (!halfEdge.isDisabled()) { - m_halfEdges.push_back({static_cast(halfEdge.m_endVertex),static_cast(halfEdge.m_opp),static_cast(halfEdge.m_face),static_cast(halfEdge.m_next)}); - halfEdgeMapping[i] = m_halfEdges.size()-1; - } - i++; - } - - for (auto& face : m_faces) { - assert(halfEdgeMapping.count(face.m_halfEdgeIndex) == 1); - face.m_halfEdgeIndex = halfEdgeMapping[face.m_halfEdgeIndex]; - } - - for (auto& he : m_halfEdges) { - he.m_face = faceMapping[he.m_face]; - he.m_opp = halfEdgeMapping[he.m_opp]; - he.m_next = halfEdgeMapping[he.m_next]; - he.m_endVertex = vertexMapping[he.m_endVertex]; - } - } - - }; - - - -// MathUtils.hpp - - namespace mathutils { - - inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p, const Ray& r) { - const glm::dvec3 s = p-r.m_S; - double t = glm::dot(s,r.m_V); - return (s.x*s.x+s.y*s.y+s.z*s.z) - t*t*r.m_VInvLengthSquared; - } - - inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { - return ((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z)); - } - // Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance). - inline double getSignedDistanceToPlane(const glm::dvec3& v, const Plane& p) { - return glm::dot(p.m_N,v) + p.m_D; - } - - inline glm::dvec3 getTriangleNormal(const glm::dvec3& a,const glm::dvec3& b,const glm::dvec3& c) { - // We want to get (a-c).crossProduct(b-c) without constructing temp vectors - double x = a.x - c.x; - double y = a.y - c.y; - double z = a.z - c.z; - double rhsx = b.x - c.x; - double rhsy = b.y - c.y; - double rhsz = b.z - c.z; - double px = y * rhsz - z * rhsy ; - double py = z * rhsx - x * rhsz ; - double pz = x * rhsy - y * rhsx ; - return glm::dvec3(px,py,pz); - } - - - } - - -// QuickHull.hpp - -/* - * Implementation of the 3D QuickHull algorithm by Antti Kuukka - * - * No copyrights. What follows is 100% Public Domain. - * - * - * - * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) - * - * OUTPUT: a ConvexHull object which provides vertex and index buffers of the generated convex hull as a triangle mesh. - * - * - * - * The implementation is thread-safe if each thread is using its own QuickHull object. - * - * - * SUMMARY OF THE ALGORITHM: - * - Create initial simplex (tetrahedron) using extreme points. We have four faces now and they form a convex mesh M. - * - For each point, assign them to the first face for which they are on the positive side of (so each point is assigned to at most - * one face). Points inside the initial tetrahedron are left behind now and no longer affect the calculations. - * - Add all faces that have points assigned to them to Face Stack. - * - Iterate until Face Stack is empty: - * - Pop topmost face F from the stack - * - From the points assigned to F, pick the point P that is farthest away from the plane defined by F. - * - Find all faces of M that have P on their positive side. Let us call these the "visible faces". - * - Because of the way M is constructed, these faces are connected. Solve their horizon edge loop. - * - "Extrude to P": Create new faces by connecting P with the points belonging to the horizon edge. Add the new faces to M and remove the visible - * faces from M. - * - Each point that was assigned to visible faces is now assigned to at most one of the newly created faces. - * - Those new faces that have points assigned to them are added to the top of Face Stack. - * - M is now the convex hull. - * - * TO DO: - * - Implement a proper 2D QuickHull and use that to solve the degenerate 2D case (when all the points lie on the same plane in 3D space). - * */ - - struct DiagnosticsData { - size_t m_failedHorizonEdges; // How many times QuickHull failed to solve the horizon edge. Failures lead to degenerated convex hulls. - - DiagnosticsData() : m_failedHorizonEdges(0) { } - }; - - double defaultEps(); - - class QuickHull { - using vec3 = glm::dvec3; - - double m_epsilon, m_epsilonSquared, m_scale; - bool m_planar; - std::vector m_planarPointCloudTemp; - VertexDataSource m_vertexData; - MeshBuilder m_mesh; - std::array m_extremeValues; - DiagnosticsData m_diagnostics; - - // Temporary variables used during iteration process - std::vector m_newFaceIndices; - std::vector m_newHalfEdgeIndices; - std::vector< std::unique_ptr> > m_disabledFacePointVectors; - std::vector m_visibleFaces; - std::vector m_horizonEdges; - struct FaceData { - size_t m_faceIndex; - size_t m_enteredFromHalfEdge; // If the face turns out not to be visible, this half edge will be marked as horizon edge - FaceData(size_t fi, size_t he) : m_faceIndex(fi),m_enteredFromHalfEdge(he) {} - }; - std::vector m_possiblyVisibleFaces; - std::deque m_faceList; - - // Create a half edge mesh representing the base tetrahedron from which the QuickHull iteration proceeds. m_extremeValues must be properly set up when this is called. - void setupInitialTetrahedron(); - - // Given a list of half edges, try to rearrange them so that they form a loop. Return true on success. - bool reorderHorizonEdges(std::vector& horizonEdges); - - // Find indices of extreme values (max x, min x, max y, min y, max z, min z) for the given point cloud - std::array getExtremeValues(); - - // Compute scale of the vertex data. - double getScale(const std::array& extremeValues); - - // Each face contains a unique pointer to a vector of indices. However, many - often most - faces do not have any points on the positive - // side of them especially at the the end of the iteration. When a face is removed from the mesh, its associated point vector, if such - // exists, is moved to the index vector pool, and when we need to add new faces with points on the positive side to the mesh, - // we reuse these vectors. This reduces the amount of std::vectors we have to deal with, and impact on performance is remarkable. - Pool m_indexVectorPool; - inline std::unique_ptr> getIndexVectorFromPool(); - inline void reclaimToIndexVectorPool(std::unique_ptr>& ptr); - - // Associates a point with a face if the point resides on the positive side of 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 m_mesh from which we create the ConvexHull object that getConvexHull function returns - void createConvexHalfEdgeMesh(); - - // Constructs the convex hull into a MeshBuilder object which can be converted to a ConvexHull or Mesh object - void buildMesh(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, double eps); - - // The public getConvexHull functions will setup a VertexDataSource object and call this - ConvexHull getConvexHull(const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices, double eps); - public: - // Computes convex hull for a given point cloud. - // Params: - // pointCloud: a vector of of 3D points - // CCW: whether the output mesh triangles should have CCW orientation - // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, - // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. - // eps: minimum distance to a plane to consider a point being on positive of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const std::vector& pointCloud, - bool CCW, - bool useOriginalIndices, - double eps = defaultEps()); - - // Computes convex hull for a given point cloud. - // Params: - // vertexData: pointer to the first 3D point of the point cloud - // vertexCount: number of vertices in the point cloud - // CCW: whether the output mesh triangles should have CCW orientation - // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, - // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. - // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const glm::dvec3* vertexData, - size_t vertexCount, - bool CCW, - bool useOriginalIndices, - double eps = defaultEps()); - - // 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 - // useOriginalIndices: should the output mesh use same vertex indices as the original point cloud. If this is false, - // then we generate a new vertex buffer which contains only the vertices that are part of the convex hull. - // eps: minimum distance to a plane to consider a point being on positive side of it (for a point cloud with scale 1) - ConvexHull getConvexHull(const double* vertexData, - size_t vertexCount, - bool CCW, - bool useOriginalIndices, - double eps = defaultEps()); - - // 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 a mesh object with half edge structure. - HalfEdgeMesh getConvexHullAsMesh(const double* vertexData, - size_t vertexCount, - bool CCW, - double eps = defaultEps()); - - // Get diagnostics about last generated convex hull - const DiagnosticsData& getDiagnostics() { - return m_diagnostics; - } - }; - - /* - * Inline function definitions - */ - - std::unique_ptr> QuickHull::getIndexVectorFromPool() { - auto r = m_indexVectorPool.get(); - r->clear(); - return r; - } - - void QuickHull::reclaimToIndexVectorPool(std::unique_ptr>& ptr) { - const size_t oldSize = ptr->size(); - if ((oldSize+1)*128 < ptr->capacity()) { - // Reduce memory usage! Huge vectors are needed at the beginning of iteration when faces have many points on their positive side. Later on, smaller vectors will suffice. - ptr.reset(nullptr); - return; - } - m_indexVectorPool.reclaim(ptr); - } - - bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { - const double D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P); - if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) { - if (!f.m_pointsOnPositiveSide) { - f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool()); - } - f.m_pointsOnPositiveSide->push_back( pointIndex ); - if (D > f.m_mostDistantPointDist) { - f.m_mostDistantPointDist = D; - f.m_mostDistantPoint = pointIndex; - } - return true; - } - return false; - } - -} - - -#endif /* QUICKHULL_HPP_ */ From 1effcd7be05dbd746ecb584a62815b8574853625 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 01:19:05 +0530 Subject: [PATCH 09/24] Resolving Conflict --- src/manifold/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/manifold/CMakeLists.txt b/src/manifold/CMakeLists.txt index 5e50295ed..ed5e30e53 100644 --- a/src/manifold/CMakeLists.txt +++ b/src/manifold/CMakeLists.txt @@ -21,8 +21,8 @@ target_include_directories(${PROJECT_NAME} PUBLIC $ $) target_link_libraries(${PROJECT_NAME} - PUBLIC utilities sdf - PRIVATE collider polygon ${MANIFOLD_INCLUDE} + PUBLIC utilities + PRIVATE $ polygon ${MANIFOLD_INCLUDE} ) target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) From 988a3e331fea64cb479ee6c6b4319768d0cbc027 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 01:29:32 +0530 Subject: [PATCH 10/24] Fixed std::move warning (Shouldn't affect performance, ideally) --- src/manifold/src/quickhull.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 28254421a..f231f3e2b 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -809,7 +809,7 @@ bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { if (!f.pointsOnPositiveSide) { - f.pointsOnPositiveSide = std::move(getIndexVectorFromPool()); + f.pointsOnPositiveSide = getIndexVectorFromPool(); } f.pointsOnPositiveSide->push_back(pointIndex); if (D > f.mostDistantPointDist) { From c1f67c734a96e49bda0bbd41eda622bd1e76bf5c Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 11:42:19 +0530 Subject: [PATCH 11/24] Reverted third_party directory in bindings --- bindings/python/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/bindings/python/CMakeLists.txt b/bindings/python/CMakeLists.txt index d8537fcd8..ad78d338f 100644 --- a/bindings/python/CMakeLists.txt +++ b/bindings/python/CMakeLists.txt @@ -14,6 +14,7 @@ project(python) +add_subdirectory(third_party) nanobind_add_module( manifold3d NB_STATIC STABLE_ABI LTO From a46ca54bf0eab61247dca7d91659532518142605 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 11:55:51 +0530 Subject: [PATCH 12/24] Reverting std::move --- src/manifold/src/quickhull.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index f231f3e2b..28254421a 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -809,7 +809,7 @@ bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { if (!f.pointsOnPositiveSide) { - f.pointsOnPositiveSide = getIndexVectorFromPool(); + f.pointsOnPositiveSide = std::move(getIndexVectorFromPool()); } f.pointsOnPositiveSide->push_back(pointIndex); if (D > f.mostDistantPointDist) { From 3e75a324abb1511201251d8421bbbeb1e66eb773 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 12:20:08 +0530 Subject: [PATCH 13/24] Corrected std::move and added glm::dot to calculate distance squared --- src/manifold/src/quickhull.h | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 28254421a..6df5918fc 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -72,9 +72,7 @@ struct Ray { const double VInvLengthSquared; Ray(const glm::dvec3& S, const glm::dvec3& V) - : S(S), - V(V), - VInvLengthSquared(1 / (V.x * V.x + V.y * V.y + V.z * V.z)) {} + : S(S), V(V), VInvLengthSquared(1 / (glm::dot(V, V))) {} }; // VertexDataSource @@ -555,12 +553,11 @@ inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p, const Ray& r) { const glm::dvec3 s = p - r.S; double t = glm::dot(s, r.V); - return (s.x * s.x + s.y * s.y + s.z * s.z) - t * t * r.VInvLengthSquared; + return glm::dot(s, s) - t * t * r.VInvLengthSquared; } inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { - return ((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + - (p1.z - p2.z) * (p1.z - p2.z)); + return (glm::dot(p1 - p2, p1 - p2)); } // Note that the unit of distance returned is relative to plane's normal's // length (divide by N.getNormalized() if needed to get the "real" distance). @@ -809,7 +806,7 @@ bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { if (!f.pointsOnPositiveSide) { - f.pointsOnPositiveSide = std::move(getIndexVectorFromPool()); + f.pointsOnPositiveSide = getIndexVectorFromPool(); } f.pointsOnPositiveSide->push_back(pointIndex); if (D > f.mostDistantPointDist) { From d818f24d72a23754a94f34ab97baaf1b097ada18 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 12:40:51 +0530 Subject: [PATCH 14/24] Reverted m_epsilon due to epsilon and m_epsilon --- src/manifold/src/quickhull.cpp | 8 ++++---- src/manifold/src/quickhull.h | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index d138d5197..d5dfcd08d 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -60,8 +60,8 @@ void QuickHull::buildMesh(const VertexDataSource& pointCloud, bool CCW, scale = getScale(extremeValues); // Epsilon we use depends on the scale - epsilon = epsilon * scale; - epsilonSquared = epsilon * epsilon; + m_epsilon = epsilon * scale; + epsilonSquared = m_epsilon * m_epsilon; // Reset diagnostics diagnostics = DiagnosticsData(); @@ -493,7 +493,7 @@ void QuickHull::setupInitialTetrahedron() { // Next step is to find the 4th vertex of the tetrahedron. We naturally choose // the point farthest away from the triangle plane. - maxD = epsilon; + maxD = m_epsilon; maxI = 0; const glm::dvec3 N = mathutils::getTriangleNormal(baseTriangleVertices[0], baseTriangleVertices[1], @@ -507,7 +507,7 @@ void QuickHull::setupInitialTetrahedron() { maxI = i; } } - if (maxD == epsilon) { + if (maxD == m_epsilon) { // All the points seem to lie on a 2D subspace of R^3. How to handle this? // Well, let's add one extra point to the point cloud so that the convex // hull will have volume. diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 6df5918fc..4bfee22f7 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -645,7 +645,7 @@ double defaultEps(); class QuickHull { using vec3 = glm::dvec3; - double epsilon, epsilonSquared, scale; + double m_epsilon, epsilonSquared, scale; bool planar; std::vector planarPointCloudTemp; VertexDataSource originalVertexData; From 7cda2726d1a0bbbfffed96882a122ad513d75d77 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 16:38:10 +0530 Subject: [PATCH 15/24] Fixed some Names --- src/manifold/src/quickhull.cpp | 58 +++++++++++----------- src/manifold/src/quickhull.h | 90 +++++++++++++++++----------------- 2 files changed, 74 insertions(+), 74 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index d5dfcd08d..55230b19a 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -90,7 +90,7 @@ ConvexHull QuickHull::getConvexHull(const VertexDataSource& pointCloud, void QuickHull::createConvexHalfEdgeMesh() { visibleFaces.clear(); - horizonEdges.clear(); + horizonEdgesData.clear(); possiblyVisibleFaces.clear(); // Compute base tetrahedron @@ -138,7 +138,7 @@ void QuickHull::createConvexHalfEdgeMesh() { // Find out the faces that have our active point on their positive side // (these are the "visible faces"). The face on top of the stack of course // is one of them. At the same time, we create a list of horizon edges. - horizonEdges.clear(); + horizonEdgesData.clear(); possiblyVisibleFaces.clear(); visibleFaces.clear(); possiblyVisibleFaces.emplace_back(topFaceIndex, @@ -146,7 +146,7 @@ void QuickHull::createConvexHalfEdgeMesh() { while (possiblyVisibleFaces.size()) { const auto faceData = possiblyVisibleFaces.back(); possiblyVisibleFaces.pop_back(); - auto& pvf = mesh.faces[faceData.faceIndex]; + auto& pvf = mesh.faces[faceData.m_faceIndex]; assert(!pvf.isDisabled()); if (pvf.visibilityCheckedOnIteration == iter) { @@ -154,46 +154,46 @@ void QuickHull::createConvexHalfEdgeMesh() { continue; } } else { - const Plane& P = pvf.P; + const Plane& P = pvf.m_P; pvf.visibilityCheckedOnIteration = iter; - const double d = glm::dot(P.N, activePoint) + P.D; + const double d = glm::dot(P.m_N, activePoint) + P.m_D; if (d > 0) { pvf.isVisibleFaceOnCurrentIteration = 1; pvf.horizonEdgesOnCurrentIteration = 0; - visibleFaces.push_back(faceData.faceIndex); + visibleFaces.push_back(faceData.m_faceIndex); for (auto heIndex : mesh.getHalfEdgeIndicesOfFace(pvf)) { if (mesh.halfEdges[heIndex].opp != faceData.enteredFromHalfEdge) { possiblyVisibleFaces.emplace_back( - mesh.halfEdges[mesh.halfEdges[heIndex].opp].face, heIndex); + mesh.halfEdges[mesh.halfEdges[heIndex].opp].m_face, heIndex); } } continue; } - assert(faceData.faceIndex != topFaceIndex); + assert(faceData.m_faceIndex != topFaceIndex); } // The face is not visible. Therefore, the halfedge we came from is part // of the horizon edge. pvf.isVisibleFaceOnCurrentIteration = 0; - horizonEdges.push_back(faceData.enteredFromHalfEdge); + horizonEdgesData.push_back(faceData.enteredFromHalfEdge); // Store which half edge is the horizon edge. The other half edges of the // face will not be part of the final mesh so their data slots can by // recycled. - const auto halfEdges = mesh.getHalfEdgeIndicesOfFace( - mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face]); + const auto halfEdgesMesh = mesh.getHalfEdgeIndicesOfFace( + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].m_face]); const std::int8_t ind = - (halfEdges[0] == faceData.enteredFromHalfEdge) + (halfEdgesMesh[0] == faceData.enteredFromHalfEdge) ? 0 - : (halfEdges[1] == faceData.enteredFromHalfEdge ? 1 : 2); - mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face] + : (halfEdgesMesh[1] == faceData.enteredFromHalfEdge ? 1 : 2); + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].m_face] .horizonEdgesOnCurrentIteration |= (1 << ind); } - const size_t horizonEdgeCount = horizonEdges.size(); + const size_t horizonEdgeCount = horizonEdgesData.size(); // Order horizon edges so that they form a loop. This may fail due to // 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(horizonEdges)) { + if (!reorderHorizonEdges(horizonEdgesData)) { diagnostics.failedHorizonEdges++; std::cerr << "Failed to solve horizon edge." << std::endl; auto it = std::find(tf.pointsOnPositiveSide->begin(), @@ -215,16 +215,16 @@ void QuickHull::createConvexHalfEdgeMesh() { size_t disableCounter = 0; for (auto faceIndex : visibleFaces) { auto& disabledFace = mesh.faces[faceIndex]; - auto halfEdges = mesh.getHalfEdgeIndicesOfFace(disabledFace); + auto halfEdgesMesh = mesh.getHalfEdgeIndicesOfFace(disabledFace); for (size_t j = 0; j < 3; j++) { if ((disabledFace.horizonEdgesOnCurrentIteration & (1 << j)) == 0) { if (disableCounter < horizonEdgeCount * 2) { // Use on this iteration - newHalfEdgeIndices.push_back(halfEdges[j]); + newHalfEdgeIndices.push_back(halfEdgesMesh[j]); disableCounter++; } else { // Mark for reusal on later iteration step - mesh.disableHalfEdge(halfEdges[j]); + mesh.disableHalfEdge(halfEdgesMesh[j]); } } } @@ -247,7 +247,7 @@ void QuickHull::createConvexHalfEdgeMesh() { // Create new faces using the edgeloop for (size_t i = 0; i < horizonEdgeCount; i++) { - const size_t AB = horizonEdges[i]; + const size_t AB = horizonEdgesData[i]; auto horizonEdgeVertexIndices = mesh.getVertexIndicesOfHalfEdge(mesh.halfEdges[AB]); @@ -266,9 +266,9 @@ void QuickHull::createConvexHalfEdgeMesh() { mesh.halfEdges[BC].next = CA; mesh.halfEdges[CA].next = AB; - mesh.halfEdges[BC].face = newFaceIndex; - mesh.halfEdges[CA].face = newFaceIndex; - mesh.halfEdges[AB].face = newFaceIndex; + mesh.halfEdges[BC].m_face = newFaceIndex; + mesh.halfEdges[CA].m_face = newFaceIndex; + mesh.halfEdges[AB].m_face = newFaceIndex; mesh.halfEdges[CA].endVertex = A; mesh.halfEdges[BC].endVertex = C; @@ -277,7 +277,7 @@ void QuickHull::createConvexHalfEdgeMesh() { const glm::dvec3 planeNormal = mathutils::getTriangleNormal( originalVertexData[A], originalVertexData[B], activePoint); - newFace.P = Plane(planeNormal, activePoint); + newFace.m_P = Plane(planeNormal, activePoint); newFace.he = AB; mesh.halfEdges[CA].opp = @@ -362,12 +362,12 @@ std::array QuickHull::getExtremeValues() { bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { const size_t horizonEdgeCount = horizonEdges.size(); for (size_t i = 0; i < horizonEdgeCount - 1; i++) { - const size_t endVertex = mesh.halfEdges[horizonEdges[i]].endVertex; + const size_t endVertexCheck = mesh.halfEdges[horizonEdges[i]].endVertex; bool foundNext = false; for (size_t j = i + 1; j < horizonEdgeCount; j++) { const size_t beginVertex = mesh.halfEdges[mesh.halfEdges[horizonEdges[j]].opp].endVertex; - if (beginVertex == endVertex) { + if (beginVertex == endVertexCheck) { std::swap(horizonEdges[i + 1], horizonEdges[j]); foundNext = true; break; @@ -382,10 +382,10 @@ bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { return true; } -double QuickHull::getScale(const std::array& extremeValues) { +double QuickHull::getScale(const std::array& extremeValuesInput) { double s = 0; for (size_t i = 0; i < 6; i++) { - const double* v = (const double*)(&originalVertexData[extremeValues[i]]); + const double* v = (const double*)(&originalVertexData[extremeValuesInput[i]]); v += i / 2; auto a = std::abs(*v); if (a > s) { @@ -542,7 +542,7 @@ void QuickHull::setupInitialTetrahedron() { const glm::dvec3& vc = originalVertexData[v[2]]; const glm::dvec3 N1 = mathutils::getTriangleNormal(va, vb, vc); const Plane plane(N1, va); - f.P = plane; + f.m_P = plane; } // Finally we assign a face for each vertex outside the tetrahedron (vertices diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 4bfee22f7..9c5cc6ba3 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -41,16 +41,16 @@ class Pool { class Plane { public: - glm::dvec3 N; + glm::dvec3 m_N; // Signed distance (if normal is of length 1) to the plane from origin - double D; + double m_D; // Normal length squared double sqrNLength; bool isPointOnPositiveSide(const glm::dvec3& Q) const { - double d = glm::dot(N, Q) + D; + double d = glm::dot(m_N, Q) + m_D; if (d >= 0) return true; return false; } @@ -59,8 +59,8 @@ class Plane { // Construct a plane using normal N and any point P on the plane Plane(const glm::dvec3& N, const glm::dvec3& P) - : N(N), - D(glm::dot(-N, P)), + : m_N(N), + m_D(glm::dot(-N, P)), sqrNLength(N.x * N.x + N.y * N.y + N.z * N.z) {} }; @@ -108,7 +108,7 @@ class MeshBuilder { struct HalfEdge { size_t endVertex; size_t opp; - size_t face; + size_t m_face; size_t next; void disable() { endVertex = std::numeric_limits::max(); } @@ -120,7 +120,7 @@ class MeshBuilder { struct Face { size_t he; - Plane P{}; + Plane m_P{}; double mostDistantPointDist; size_t mostDistantPoint; size_t visibilityCheckedOnIteration; @@ -212,84 +212,84 @@ class MeshBuilder { HalfEdge AB; AB.endVertex = b; AB.opp = 6; - AB.face = 0; + AB.m_face = 0; AB.next = 1; halfEdges.push_back(AB); HalfEdge BC; BC.endVertex = c; BC.opp = 9; - BC.face = 0; + BC.m_face = 0; BC.next = 2; halfEdges.push_back(BC); HalfEdge CA; CA.endVertex = a; CA.opp = 3; - CA.face = 0; + CA.m_face = 0; CA.next = 0; halfEdges.push_back(CA); HalfEdge AC; AC.endVertex = c; AC.opp = 2; - AC.face = 1; + AC.m_face = 1; AC.next = 4; halfEdges.push_back(AC); HalfEdge CD; CD.endVertex = d; CD.opp = 11; - CD.face = 1; + CD.m_face = 1; CD.next = 5; halfEdges.push_back(CD); HalfEdge DA; DA.endVertex = a; DA.opp = 7; - DA.face = 1; + DA.m_face = 1; DA.next = 3; halfEdges.push_back(DA); HalfEdge BA; BA.endVertex = a; BA.opp = 0; - BA.face = 2; + BA.m_face = 2; BA.next = 7; halfEdges.push_back(BA); HalfEdge AD; AD.endVertex = d; AD.opp = 5; - AD.face = 2; + AD.m_face = 2; AD.next = 8; halfEdges.push_back(AD); HalfEdge DB; DB.endVertex = b; DB.opp = 10; - DB.face = 2; + DB.m_face = 2; DB.next = 6; halfEdges.push_back(DB); HalfEdge CB; CB.endVertex = b; CB.opp = 1; - CB.face = 3; + CB.m_face = 3; CB.next = 10; halfEdges.push_back(CB); HalfEdge BD; BD.endVertex = d; BD.opp = 8; - BD.face = 3; + BD.m_face = 3; BD.next = 11; halfEdges.push_back(BD); HalfEdge DC; DC.endVertex = c; DC.opp = 4; - DC.face = 3; + DC.m_face = 3; DC.next = 9; halfEdges.push_back(DC); @@ -395,19 +395,19 @@ class ConvexHull { } // Construct vertex and index buffers from half edge mesh and pointcloud - ConvexHull(const MeshBuilder& mesh, const VertexDataSource& pointCloud, + ConvexHull(const MeshBuilder& mesh_input, const VertexDataSource& pointCloud, bool CCW, bool useOriginalIndices) { if (!useOriginalIndices) { optimizedVertexBuffer.reset(new std::vector()); } - std::vector faceProcessed(mesh.faces.size(), false); + std::vector faceProcessed(mesh_input.faces.size(), false); std::vector faceStack; // Map vertex indices from original point cloud to the new mesh vertex // indices std::unordered_map vertexIndexMapping; - for (size_t i = 0; i < mesh.faces.size(); i++) { - if (!mesh.faces[i].isDisabled()) { + for (size_t i = 0; i < mesh_input.faces.size(); i++) { + if (!mesh_input.faces[i].isDisabled()) { faceStack.push_back(i); break; } @@ -418,31 +418,31 @@ class ConvexHull { const size_t iCCW = CCW ? 1 : 0; const size_t finalMeshFaceCount = - mesh.faces.size() - mesh.disabledFaces.size(); + mesh_input.faces.size() - mesh_input.disabledFaces.size(); indices.reserve(finalMeshFaceCount * 3); while (faceStack.size()) { auto it = faceStack.end() - 1; size_t top = *it; - assert(!mesh.faces[top].isDisabled()); + assert(!mesh_input.faces[top].isDisabled()); faceStack.erase(it); if (faceProcessed[top]) { continue; } else { faceProcessed[top] = true; - auto halfEdges = mesh.getHalfEdgeIndicesOfFace(mesh.faces[top]); + auto halfEdgesMesh = mesh_input.getHalfEdgeIndicesOfFace(mesh_input.faces[top]); size_t adjacent[] = { - mesh.halfEdges[mesh.halfEdges[halfEdges[0]].opp].face, - mesh.halfEdges[mesh.halfEdges[halfEdges[1]].opp].face, - mesh.halfEdges[mesh.halfEdges[halfEdges[2]].opp].face}; + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[0]].opp].m_face, + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[1]].opp].m_face, + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[2]].opp].m_face}; for (auto a : adjacent) { - if (!faceProcessed[a] && !mesh.faces[a].isDisabled()) { + if (!faceProcessed[a] && !mesh_input.faces[a].isDisabled()) { faceStack.push_back(a); } } - auto vertices = mesh.getVertexIndicesOfFace(mesh.faces[top]); + auto MeshVertices = mesh_input.getVertexIndicesOfFace(mesh_input.faces[top]); if (!useOriginalIndices) { - for (auto& v : vertices) { + for (auto& v : MeshVertices) { auto itV = vertexIndexMapping.find(v); if (itV == vertexIndexMapping.end()) { optimizedVertexBuffer->push_back(pointCloud[v]); @@ -453,9 +453,9 @@ class ConvexHull { } } } - indices.push_back(vertices[0]); - indices.push_back(vertices[1 + iCCW]); - indices.push_back(vertices[2 - iCCW]); + indices.push_back(MeshVertices[0]); + indices.push_back(MeshVertices[1 + iCCW]); + indices.push_back(MeshVertices[2 - iCCW]); } } @@ -482,7 +482,7 @@ class HalfEdgeMesh { struct HalfEdge { size_t endVertex; size_t opp; - size_t face; + size_t m_face; size_t next; }; @@ -524,7 +524,7 @@ class HalfEdgeMesh { if (!halfEdge.isDisabled()) { halfEdges.push_back({static_cast(halfEdge.endVertex), static_cast(halfEdge.opp), - static_cast(halfEdge.face), + static_cast(halfEdge.m_face), static_cast(halfEdge.next)}); halfEdgeMapping[i] = halfEdges.size() - 1; } @@ -537,7 +537,7 @@ class HalfEdgeMesh { } for (auto& he : halfEdges) { - he.face = faceMapping[he.face]; + he.m_face = faceMapping[he.m_face]; he.opp = halfEdgeMapping[he.opp]; he.next = halfEdgeMapping[he.next]; he.endVertex = vertexMapping[he.endVertex]; @@ -562,7 +562,7 @@ inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { // Note that the unit of distance returned is relative to plane's normal's // length (divide by N.getNormalized() if needed to get the "real" distance). inline double getSignedDistanceToPlane(const glm::dvec3& v, const Plane& p) { - return glm::dot(p.N, v) + p.D; + return glm::dot(p.m_N, v) + p.m_D; } inline glm::dvec3 getTriangleNormal(const glm::dvec3& a, const glm::dvec3& b, @@ -658,13 +658,13 @@ class QuickHull { std::vector newHalfEdgeIndices; std::vector>> disabledFacePointVectors; std::vector visibleFaces; - std::vector horizonEdges; + std::vector horizonEdgesData; struct FaceData { - size_t faceIndex; + size_t m_faceIndex; // If the face turns out not to be visible, this half edge will be marked as // horizon edge size_t enteredFromHalfEdge; - FaceData(size_t fi, size_t he) : faceIndex(fi), enteredFromHalfEdge(he) {} + FaceData(size_t fi, size_t he) : m_faceIndex(fi), enteredFromHalfEdge(he) {} }; std::vector possiblyVisibleFaces; std::deque faceList; @@ -683,7 +683,7 @@ class QuickHull { std::array getExtremeValues(); // Compute scale of the vertex data. - double getScale(const std::array& extremeValues); + double getScale(const std::array& extremeValuesInput); // Each face contains a unique pointer to a vector of indices. However, many - // often most - faces do not have any points on the positive side of them @@ -803,8 +803,8 @@ void QuickHull::reclaimToIndexVectorPool( bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { const double D = - mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); - if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { + mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.m_P); + if (D > 0 && D * D > epsilonSquared * f.m_P.sqrNLength) { if (!f.pointsOnPositiveSide) { f.pointsOnPositiveSide = getIndexVectorFromPool(); } From c32dc4ea239c5ec91c312e4d017f6159385c2def Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 16:38:32 +0530 Subject: [PATCH 16/24] Formatting --- src/manifold/src/quickhull.cpp | 3 ++- src/manifold/src/quickhull.h | 19 ++++++++++++------- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 55230b19a..8fd0b70ab 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -385,7 +385,8 @@ bool QuickHull::reorderHorizonEdges(std::vector& horizonEdges) { double QuickHull::getScale(const std::array& extremeValuesInput) { double s = 0; for (size_t i = 0; i < 6; i++) { - const double* v = (const double*)(&originalVertexData[extremeValuesInput[i]]); + const double* v = + (const double*)(&originalVertexData[extremeValuesInput[i]]); v += i / 2; auto a = std::abs(*v); if (a > s) { diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 9c5cc6ba3..a477af682 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -430,17 +430,22 @@ class ConvexHull { continue; } else { faceProcessed[top] = true; - auto halfEdgesMesh = mesh_input.getHalfEdgeIndicesOfFace(mesh_input.faces[top]); + auto halfEdgesMesh = + mesh_input.getHalfEdgeIndicesOfFace(mesh_input.faces[top]); size_t adjacent[] = { - mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[0]].opp].m_face, - mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[1]].opp].m_face, - mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[2]].opp].m_face}; + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[0]].opp] + .m_face, + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[1]].opp] + .m_face, + mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[2]].opp] + .m_face}; for (auto a : adjacent) { if (!faceProcessed[a] && !mesh_input.faces[a].isDisabled()) { faceStack.push_back(a); } } - auto MeshVertices = mesh_input.getVertexIndicesOfFace(mesh_input.faces[top]); + auto MeshVertices = + mesh_input.getVertexIndicesOfFace(mesh_input.faces[top]); if (!useOriginalIndices) { for (auto& v : MeshVertices) { auto itV = vertexIndexMapping.find(v); @@ -802,8 +807,8 @@ void QuickHull::reclaimToIndexVectorPool( bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { - const double D = - mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.m_P); + const double D = mathutils::getSignedDistanceToPlane( + originalVertexData[pointIndex], f.m_P); if (D > 0 && D * D > epsilonSquared * f.m_P.sqrNLength) { if (!f.pointsOnPositiveSide) { f.pointsOnPositiveSide = getIndexVectorFromPool(); From 00dd951c03d0b1117c1efe3f86d27f6f6df5dd22 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 16:49:05 +0530 Subject: [PATCH 17/24] Removed all unecessary m_ --- src/manifold/src/quickhull.cpp | 26 +++++++------- src/manifold/src/quickhull.h | 62 +++++++++++++++++----------------- 2 files changed, 44 insertions(+), 44 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 8fd0b70ab..2f13ef195 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -146,7 +146,7 @@ void QuickHull::createConvexHalfEdgeMesh() { while (possiblyVisibleFaces.size()) { const auto faceData = possiblyVisibleFaces.back(); possiblyVisibleFaces.pop_back(); - auto& pvf = mesh.faces[faceData.m_faceIndex]; + auto& pvf = mesh.faces[faceData.faceIndex]; assert(!pvf.isDisabled()); if (pvf.visibilityCheckedOnIteration == iter) { @@ -154,22 +154,22 @@ void QuickHull::createConvexHalfEdgeMesh() { continue; } } else { - const Plane& P = pvf.m_P; + const Plane& P = pvf.P; pvf.visibilityCheckedOnIteration = iter; - const double d = glm::dot(P.m_N, activePoint) + P.m_D; + const double d = glm::dot(P.N, activePoint) + P.D; if (d > 0) { pvf.isVisibleFaceOnCurrentIteration = 1; pvf.horizonEdgesOnCurrentIteration = 0; - visibleFaces.push_back(faceData.m_faceIndex); + visibleFaces.push_back(faceData.faceIndex); for (auto heIndex : mesh.getHalfEdgeIndicesOfFace(pvf)) { if (mesh.halfEdges[heIndex].opp != faceData.enteredFromHalfEdge) { possiblyVisibleFaces.emplace_back( - mesh.halfEdges[mesh.halfEdges[heIndex].opp].m_face, heIndex); + mesh.halfEdges[mesh.halfEdges[heIndex].opp].face, heIndex); } } continue; } - assert(faceData.m_faceIndex != topFaceIndex); + assert(faceData.faceIndex != topFaceIndex); } // The face is not visible. Therefore, the halfedge we came from is part @@ -180,12 +180,12 @@ void QuickHull::createConvexHalfEdgeMesh() { // face will not be part of the final mesh so their data slots can by // recycled. const auto halfEdgesMesh = mesh.getHalfEdgeIndicesOfFace( - mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].m_face]); + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face]); const std::int8_t ind = (halfEdgesMesh[0] == faceData.enteredFromHalfEdge) ? 0 : (halfEdgesMesh[1] == faceData.enteredFromHalfEdge ? 1 : 2); - mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].m_face] + mesh.faces[mesh.halfEdges[faceData.enteredFromHalfEdge].face] .horizonEdgesOnCurrentIteration |= (1 << ind); } const size_t horizonEdgeCount = horizonEdgesData.size(); @@ -266,9 +266,9 @@ void QuickHull::createConvexHalfEdgeMesh() { mesh.halfEdges[BC].next = CA; mesh.halfEdges[CA].next = AB; - mesh.halfEdges[BC].m_face = newFaceIndex; - mesh.halfEdges[CA].m_face = newFaceIndex; - mesh.halfEdges[AB].m_face = newFaceIndex; + mesh.halfEdges[BC].face = newFaceIndex; + mesh.halfEdges[CA].face = newFaceIndex; + mesh.halfEdges[AB].face = newFaceIndex; mesh.halfEdges[CA].endVertex = A; mesh.halfEdges[BC].endVertex = C; @@ -277,7 +277,7 @@ void QuickHull::createConvexHalfEdgeMesh() { const glm::dvec3 planeNormal = mathutils::getTriangleNormal( originalVertexData[A], originalVertexData[B], activePoint); - newFace.m_P = Plane(planeNormal, activePoint); + newFace.P = Plane(planeNormal, activePoint); newFace.he = AB; mesh.halfEdges[CA].opp = @@ -543,7 +543,7 @@ void QuickHull::setupInitialTetrahedron() { const glm::dvec3& vc = originalVertexData[v[2]]; const glm::dvec3 N1 = mathutils::getTriangleNormal(va, vb, vc); const Plane plane(N1, va); - f.m_P = plane; + f.P = plane; } // Finally we assign a face for each vertex outside the tetrahedron (vertices diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index a477af682..c61073511 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -41,16 +41,16 @@ class Pool { class Plane { public: - glm::dvec3 m_N; + glm::dvec3 N; // Signed distance (if normal is of length 1) to the plane from origin - double m_D; + double D; // Normal length squared double sqrNLength; bool isPointOnPositiveSide(const glm::dvec3& Q) const { - double d = glm::dot(m_N, Q) + m_D; + double d = glm::dot(N, Q) + D; if (d >= 0) return true; return false; } @@ -59,8 +59,8 @@ class Plane { // Construct a plane using normal N and any point P on the plane Plane(const glm::dvec3& N, const glm::dvec3& P) - : m_N(N), - m_D(glm::dot(-N, P)), + : N(N), + D(glm::dot(-N, P)), sqrNLength(N.x * N.x + N.y * N.y + N.z * N.z) {} }; @@ -108,7 +108,7 @@ class MeshBuilder { struct HalfEdge { size_t endVertex; size_t opp; - size_t m_face; + size_t face; size_t next; void disable() { endVertex = std::numeric_limits::max(); } @@ -120,7 +120,7 @@ class MeshBuilder { struct Face { size_t he; - Plane m_P{}; + Plane P{}; double mostDistantPointDist; size_t mostDistantPoint; size_t visibilityCheckedOnIteration; @@ -212,84 +212,84 @@ class MeshBuilder { HalfEdge AB; AB.endVertex = b; AB.opp = 6; - AB.m_face = 0; + AB.face = 0; AB.next = 1; halfEdges.push_back(AB); HalfEdge BC; BC.endVertex = c; BC.opp = 9; - BC.m_face = 0; + BC.face = 0; BC.next = 2; halfEdges.push_back(BC); HalfEdge CA; CA.endVertex = a; CA.opp = 3; - CA.m_face = 0; + CA.face = 0; CA.next = 0; halfEdges.push_back(CA); HalfEdge AC; AC.endVertex = c; AC.opp = 2; - AC.m_face = 1; + AC.face = 1; AC.next = 4; halfEdges.push_back(AC); HalfEdge CD; CD.endVertex = d; CD.opp = 11; - CD.m_face = 1; + CD.face = 1; CD.next = 5; halfEdges.push_back(CD); HalfEdge DA; DA.endVertex = a; DA.opp = 7; - DA.m_face = 1; + DA.face = 1; DA.next = 3; halfEdges.push_back(DA); HalfEdge BA; BA.endVertex = a; BA.opp = 0; - BA.m_face = 2; + BA.face = 2; BA.next = 7; halfEdges.push_back(BA); HalfEdge AD; AD.endVertex = d; AD.opp = 5; - AD.m_face = 2; + AD.face = 2; AD.next = 8; halfEdges.push_back(AD); HalfEdge DB; DB.endVertex = b; DB.opp = 10; - DB.m_face = 2; + DB.face = 2; DB.next = 6; halfEdges.push_back(DB); HalfEdge CB; CB.endVertex = b; CB.opp = 1; - CB.m_face = 3; + CB.face = 3; CB.next = 10; halfEdges.push_back(CB); HalfEdge BD; BD.endVertex = d; BD.opp = 8; - BD.m_face = 3; + BD.face = 3; BD.next = 11; halfEdges.push_back(BD); HalfEdge DC; DC.endVertex = c; DC.opp = 4; - DC.m_face = 3; + DC.face = 3; DC.next = 9; halfEdges.push_back(DC); @@ -434,11 +434,11 @@ class ConvexHull { mesh_input.getHalfEdgeIndicesOfFace(mesh_input.faces[top]); size_t adjacent[] = { mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[0]].opp] - .m_face, + .face, mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[1]].opp] - .m_face, + .face, mesh_input.halfEdges[mesh_input.halfEdges[halfEdgesMesh[2]].opp] - .m_face}; + .face}; for (auto a : adjacent) { if (!faceProcessed[a] && !mesh_input.faces[a].isDisabled()) { faceStack.push_back(a); @@ -487,7 +487,7 @@ class HalfEdgeMesh { struct HalfEdge { size_t endVertex; size_t opp; - size_t m_face; + size_t face; size_t next; }; @@ -529,7 +529,7 @@ class HalfEdgeMesh { if (!halfEdge.isDisabled()) { halfEdges.push_back({static_cast(halfEdge.endVertex), static_cast(halfEdge.opp), - static_cast(halfEdge.m_face), + static_cast(halfEdge.face), static_cast(halfEdge.next)}); halfEdgeMapping[i] = halfEdges.size() - 1; } @@ -542,7 +542,7 @@ class HalfEdgeMesh { } for (auto& he : halfEdges) { - he.m_face = faceMapping[he.m_face]; + he.face = faceMapping[he.face]; he.opp = halfEdgeMapping[he.opp]; he.next = halfEdgeMapping[he.next]; he.endVertex = vertexMapping[he.endVertex]; @@ -567,7 +567,7 @@ inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { // Note that the unit of distance returned is relative to plane's normal's // length (divide by N.getNormalized() if needed to get the "real" distance). inline double getSignedDistanceToPlane(const glm::dvec3& v, const Plane& p) { - return glm::dot(p.m_N, v) + p.m_D; + return glm::dot(p.N, v) + p.D; } inline glm::dvec3 getTriangleNormal(const glm::dvec3& a, const glm::dvec3& b, @@ -665,11 +665,11 @@ class QuickHull { std::vector visibleFaces; std::vector horizonEdgesData; struct FaceData { - size_t m_faceIndex; + size_t faceIndex; // If the face turns out not to be visible, this half edge will be marked as // horizon edge size_t enteredFromHalfEdge; - FaceData(size_t fi, size_t he) : m_faceIndex(fi), enteredFromHalfEdge(he) {} + FaceData(size_t fi, size_t he) : faceIndex(fi), enteredFromHalfEdge(he) {} }; std::vector possiblyVisibleFaces; std::deque faceList; @@ -807,9 +807,9 @@ void QuickHull::reclaimToIndexVectorPool( bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, size_t pointIndex) { - const double D = mathutils::getSignedDistanceToPlane( - originalVertexData[pointIndex], f.m_P); - if (D > 0 && D * D > epsilonSquared * f.m_P.sqrNLength) { + const double D = + mathutils::getSignedDistanceToPlane(originalVertexData[pointIndex], f.P); + if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) { if (!f.pointsOnPositiveSide) { f.pointsOnPositiveSide = getIndexVectorFromPool(); } From ca8484fba1be546eddd5a42a98ccd8d639720bcc Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Fri, 2 Aug 2024 20:46:01 +0530 Subject: [PATCH 18/24] Fixed manifold.yaml --- .github/workflows/manifold.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/manifold.yml b/.github/workflows/manifold.yml index fc9a8390c..e49661898 100644 --- a/.github/workflows/manifold.yml +++ b/.github/workflows/manifold.yml @@ -65,7 +65,7 @@ jobs: cd ../ lcov --capture --gcov-tool gcov-${{ matrix.gcc }} --ignore-errors mismatch --directory . --output-file ./code_coverage_test.info lcov --add-tracefile ./code_coverage_init.info --add-tracefile ./code_coverage_test.info --output-file ./code_coverage_total.info - lcov --remove ./code_coverage_total.info '/usr/*' '*/third_party/*' '*/test/*' '*/extras/*' '*/bindings/*' --output-file ./code_coverage.info + lcov --remove ./code_coverage_total.info '/usr/*' '*/test/*' '*/extras/*' '*/bindings/*' --output-file ./code_coverage.info cd ../ - uses: codecov/codecov-action@v4 if: matrix.parallel_backend == 'NONE' From 9c74ffe3b0788fd4b1a87fa631389249cdd7794c Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Sat, 3 Aug 2024 23:03:42 +0530 Subject: [PATCH 19/24] Added headers and pragma once --- src/manifold/src/quickhull.cpp | 19 ++++++ src/manifold/src/quickhull.h | 112 +++++++++++++++++---------------- 2 files changed, 78 insertions(+), 53 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index 2f13ef195..e3bcb7ba1 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -1,3 +1,22 @@ +// Copyright 2024 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +/* + * Implementation of the 3D QuickHull algorithm by Antti Kuukka + * + * No copyrights. What follows is 100% Public Domain. + * */ #include "quickhull.h" #include diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index c61073511..412fea999 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -1,5 +1,61 @@ -#ifndef QUICKHULL_HPP_ -#define QUICKHULL_HPP_ +// Copyright 2024 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +/* + * Implementation of the 3D QuickHull algorithm by Antti Kuukka + * + * No copyrights. What follows is 100% Public Domain. + * + * + * + * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) + * + * OUTPUT: a ConvexHull object which provides vertex and index buffers of the + *generated convex hull as a triangle mesh. + * + * + * + * The implementation is thread-safe if each thread is using its own QuickHull + *object. + * + * + * SUMMARY OF THE ALGORITHM: + * - Create initial simplex (tetrahedron) using extreme points. We have + *four faces now and they form a convex mesh M. + * - For each point, assign them to the first face for which they are on + *the positive side of (so each point is assigned to at most one face). Points + *inside the initial tetrahedron are left behind now and no longer affect the + *calculations. + * - Add all faces that have points assigned to them to Face Stack. + * - Iterate until Face Stack is empty: + * - Pop topmost face F from the stack + * - From the points assigned to F, pick the point P that is + *farthest away from the plane defined by F. + * - Find all faces of M that have P on their positive side. Let us + *call these the "visible faces". + * - Because of the way M is constructed, these faces are + *connected. Solve their horizon edge loop. + * - "Extrude to P": Create new faces by connecting + *P with the points belonging to the horizon edge. Add the new faces to M and + *remove the visible faces from M. + * - Each point that was assigned to visible faces is now assigned + *to at most one of the newly created faces. + * - Those new faces that have points assigned to them are added to + *the top of Face Stack. + * - M is now the convex hull. + * + * */ +#pragma once #include #include #include @@ -562,7 +618,7 @@ inline double getSquaredDistanceBetweenPointAndRay(const glm::dvec3& p, } inline double getSquaredDistance(const glm::dvec3& p1, const glm::dvec3& p2) { - return (glm::dot(p1 - p2, p1 - p2)); + return glm::dot(p1 - p2, p1 - p2); } // Note that the unit of distance returned is relative to plane's normal's // length (divide by N.getNormalized() if needed to get the "real" distance). @@ -589,54 +645,6 @@ inline glm::dvec3 getTriangleNormal(const glm::dvec3& a, const glm::dvec3& b, // QuickHull.hpp -/* - * Implementation of the 3D QuickHull algorithm by Antti Kuukka - * - * No copyrights. What follows is 100% Public Domain. - * - * - * - * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) - * - * OUTPUT: a ConvexHull object which provides vertex and index buffers of the - *generated convex hull as a triangle mesh. - * - * - * - * The implementation is thread-safe if each thread is using its own QuickHull - *object. - * - * - * SUMMARY OF THE ALGORITHM: - * - Create initial simplex (tetrahedron) using extreme points. We have - *four faces now and they form a convex mesh M. - * - For each point, assign them to the first face for which they are on - *the positive side of (so each point is assigned to at most one face). Points - *inside the initial tetrahedron are left behind now and no longer affect the - *calculations. - * - Add all faces that have points assigned to them to Face Stack. - * - Iterate until Face Stack is empty: - * - Pop topmost face F from the stack - * - From the points assigned to F, pick the point P that is - *farthest away from the plane defined by F. - * - Find all faces of M that have P on their positive side. Let us - *call these the "visible faces". - * - Because of the way M is constructed, these faces are - *connected. Solve their horizon edge loop. - * - "Extrude to P": Create new faces by connecting - *P with the points belonging to the horizon edge. Add the new faces to M and - *remove the visible faces from M. - * - Each point that was assigned to visible faces is now assigned - *to at most one of the newly created faces. - * - Those new faces that have points assigned to them are added to - *the top of Face Stack. - * - M is now the convex hull. - * - * TO DO: - * - Implement a proper 2D QuickHull and use that to solve the degenerate 2D - *case (when all the points lie on the same plane in 3D space). - * */ - struct DiagnosticsData { // How many times QuickHull failed to solve the horizon edge. Failures lead // to degenerated convex hulls. @@ -822,5 +830,3 @@ bool QuickHull::addPointToFace(typename MeshBuilder::Face& f, } return false; } - -#endif /* QUICKHULL_HPP_ */ From 5180b15958ac2cfc23b3e85e9507d69bf56577ab Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Mon, 5 Aug 2024 22:29:30 +0530 Subject: [PATCH 20/24] Modified License --- src/manifold/src/quickhull.cpp | 7 ++----- src/manifold/src/quickhull.h | 9 +++------ 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index e3bcb7ba1..ef0d4a3d0 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -11,12 +11,9 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. +// +// Derived from the public domain work of Antti Kuukka at https://github.com/akuukka/quickhull -/* - * Implementation of the 3D QuickHull algorithm by Antti Kuukka - * - * No copyrights. What follows is 100% Public Domain. - * */ #include "quickhull.h" #include diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 412fea999..9fcb6b5a4 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -11,13 +11,10 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. +// +// Derived from the public domain work of Antti Kuukka at https://github.com/akuukka/quickhull + /* - * Implementation of the 3D QuickHull algorithm by Antti Kuukka - * - * No copyrights. What follows is 100% Public Domain. - * - * - * * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) * * OUTPUT: a ConvexHull object which provides vertex and index buffers of the From 08ec86890fa210f0079e1bc2761568e17af1b6e6 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Mon, 5 Aug 2024 22:30:12 +0530 Subject: [PATCH 21/24] Formatting --- src/manifold/src/quickhull.cpp | 5 +++-- src/manifold/src/quickhull.h | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/manifold/src/quickhull.cpp b/src/manifold/src/quickhull.cpp index ef0d4a3d0..ff4b0c5ea 100644 --- a/src/manifold/src/quickhull.cpp +++ b/src/manifold/src/quickhull.cpp @@ -11,8 +11,9 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. -// -// Derived from the public domain work of Antti Kuukka at https://github.com/akuukka/quickhull +// +// Derived from the public domain work of Antti Kuukka at +// https://github.com/akuukka/quickhull #include "quickhull.h" diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 9fcb6b5a4..8bce5fda9 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -11,8 +11,9 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. -// -// Derived from the public domain work of Antti Kuukka at https://github.com/akuukka/quickhull +// +// Derived from the public domain work of Antti Kuukka at +// https://github.com/akuukka/quickhull /* * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh) From 81f31f02ef96a280aa6f158432afbde3d3fc61ee Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 8 Aug 2024 02:27:43 +0530 Subject: [PATCH 22/24] Activated the two failing tests and added the fix for the algorithm --- bindings/python/third_party/nanobind | 1 + src/manifold/src/quickhull.h | 2 +- test/hull_test.cpp | 22 ++++++++++++++++------ 3 files changed, 18 insertions(+), 7 deletions(-) create mode 160000 bindings/python/third_party/nanobind diff --git a/bindings/python/third_party/nanobind b/bindings/python/third_party/nanobind new file mode 160000 index 000000000..8d7f1ee06 --- /dev/null +++ b/bindings/python/third_party/nanobind @@ -0,0 +1 @@ +Subproject commit 8d7f1ee0621c17fa370b704b2100ffa0243d5bfb diff --git a/src/manifold/src/quickhull.h b/src/manifold/src/quickhull.h index 8bce5fda9..1959a84fc 100644 --- a/src/manifold/src/quickhull.h +++ b/src/manifold/src/quickhull.h @@ -636,7 +636,7 @@ inline glm::dvec3 getTriangleNormal(const glm::dvec3& a, const glm::dvec3& b, double px = y * rhsz - z * rhsy; double py = z * rhsx - x * rhsz; double pz = x * rhsy - y * rhsx; - return glm::dvec3(px, py, pz); + return glm::normalize(glm::dvec3(px, py, pz)); } } // namespace mathutils diff --git a/test/hull_test.cpp b/test/hull_test.cpp index 25bac36b6..50f34b26b 100644 --- a/test/hull_test.cpp +++ b/test/hull_test.cpp @@ -23,7 +23,7 @@ using namespace manifold; // Check if the mesh remains convex after adding new faces -bool isMeshConvex(manifold::Manifold hullManifold) { +bool isMeshConvex(manifold::Manifold hullManifold, double epsilon = 0.0000001) { // Get the mesh from the manifold manifold::Mesh mesh = hullManifold.GetMesh(); @@ -51,7 +51,7 @@ bool isMeshConvex(manifold::Manifold hullManifold) { double distance = glm::dot(normal, v - v0); // If any vertex lies on the opposite side of the normal direction - if (distance > 0) { + if (distance > epsilon) { // The manifold is not convex return false; } @@ -134,7 +134,7 @@ TEST(Hull, Sphere) { sphere.GetProperties().volume); } -TEST(Hull, DISABLED_FailingTest1) { +TEST(Hull, FailingTest1) { // 39202.stl const std::vector hullPts = { {-24.983196259f, -43.272167206f, 52.710712433f}, @@ -158,10 +158,15 @@ TEST(Hull, DISABLED_FailingTest1) { {-20.442623138f, -35.407661438f, 8.2749996185f}, {10.229375839f, -14.717799187f, 10.508025169f}}; auto hull = Manifold::Hull(hullPts); - EXPECT_TRUE(isMeshConvex(hull)); +#ifdef MANIFOLD_EXPORT + if (options.exportModels) { + ExportMesh("failing_test1.glb", hull.GetMesh(), {}); + } +#endif + EXPECT_TRUE(isMeshConvex(hull, 8.99628e-06)); } -TEST(Hull, DISABLED_FailingTest2) { +TEST(Hull, FailingTest2) { // 1750623.stl const std::vector hullPts = { {174.17001343f, -12.022000313f, 29.562002182f}, @@ -186,5 +191,10 @@ TEST(Hull, DISABLED_FailingTest2) { {174.17001343f, -19.502000809f, 29.562002182f}, {190.06401062f, -0.81000006199f, -14.250000954f}}; auto hull = Manifold::Hull(hullPts); - EXPECT_TRUE(isMeshConvex(hull)); +#ifdef MANIFOLD_EXPORT + if (options.exportModels) { + ExportMesh("failing_test2.glb", hull.GetMesh(), {}); + } +#endif + EXPECT_TRUE(isMeshConvex(hull, 2.13966e-05)); } \ No newline at end of file From 462c917df8af61ce00cf73257a564bce09efcf33 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 8 Aug 2024 02:29:51 +0530 Subject: [PATCH 23/24] Modified gitmodules to remove third_party --- .gitmodules | 3 --- 1 file changed, 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index e518ba769..e69de29bb 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +0,0 @@ -[submodule "src/third_party/quickhull"] - path = src/third_party/quickhull - url = https://github.com/akuukka/quickhull From 95f97218e44ee6036fc65b4ccec7f4256945c793 Mon Sep 17 00:00:00 2001 From: Kushal-Shah-03 Date: Thu, 8 Aug 2024 02:38:42 +0530 Subject: [PATCH 24/24] Removed nanobind --- bindings/python/third_party/nanobind | 1 - 1 file changed, 1 deletion(-) delete mode 160000 bindings/python/third_party/nanobind diff --git a/bindings/python/third_party/nanobind b/bindings/python/third_party/nanobind deleted file mode 160000 index 8d7f1ee06..000000000 --- a/bindings/python/third_party/nanobind +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 8d7f1ee0621c17fa370b704b2100ffa0243d5bfb