// 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.
//
// Derived from the public domain work of Antti Kuukka at
// https://github.com/akuukka/quickhull
#include "quickhull.h"
#include <algorithm>
#include <cstddef>
#include <limits>
#include <unordered_map>
#include "impl.h"
namespace manifold {
double defaultEps() { return 0.0000001; }
inline double getSquaredDistanceBetweenPointAndRay(const vec3& p,
const Ray& r) {
const vec3 s = p - r.S;
double t = la::dot(s, r.V);
return la::dot(s, s) - t * t * r.VInvLengthSquared;
}
inline double getSquaredDistance(const vec3& p1, const vec3& p2) {
return la::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).
inline double getSignedDistanceToPlane(const vec3& v, const Plane& p) {
return la::dot(p.N, v) + p.D;
}
inline vec3 getTriangleNormal(const vec3& a, const vec3& b, const vec3& 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 la::normalize(vec3(px, py, pz));
}
size_t MeshBuilder::addFace() {
if (disabledFaces.size()) {
size_t index = disabledFaces.back();
auto& f = faces[index];
DEBUG_ASSERT(f.isDisabled(), logicErr, "f should be disabled");
DEBUG_ASSERT(!f.pointsOnPositiveSide, logicErr,
"f should not be on the positive side");
f.mostDistantPointDist = 0;
disabledFaces.pop_back();
return index;
}
faces.emplace_back();
return faces.size() - 1;
}
size_t MeshBuilder::addHalfedge() {
if (disabledHalfedges.size()) {
const size_t index = disabledHalfedges.back();
disabledHalfedges.pop_back();
return index;
}
halfedges.push_back({});
halfedgeToFace.push_back(0);
halfedgeNext.push_back(0);
return halfedges.size() - 1;
}
void MeshBuilder::setup(int a, int b, int c, int d) {
faces.clear();
halfedges.clear();
halfedgeToFace.clear();
halfedgeNext.clear();
disabledFaces.clear();
disabledHalfedges.clear();
faces.reserve(4);
halfedges.reserve(12);
// Create halfedges
// AB
halfedges.push_back({0, b, 6});
halfedgeToFace.push_back(0);
halfedgeNext.push_back(1);
// BC
halfedges.push_back({0, c, 9});
halfedgeToFace.push_back(0);
halfedgeNext.push_back(2);
// CA
halfedges.push_back({0, a, 3});
halfedgeToFace.push_back(0);
halfedgeNext.push_back(0);
// AC
halfedges.push_back({0, c, 2});
halfedgeToFace.push_back(1);
halfedgeNext.push_back(4);
// CD
halfedges.push_back({0, d, 11});
halfedgeToFace.push_back(1);
halfedgeNext.push_back(5);
// DA
halfedges.push_back({0, a, 7});
halfedgeToFace.push_back(1);
halfedgeNext.push_back(3);
// BA
halfedges.push_back({0, a, 0});
halfedgeToFace.push_back(2);
halfedgeNext.push_back(7);
// AD
halfedges.push_back({0, d, 5});
halfedgeToFace.push_back(2);
halfedgeNext.push_back(8);
// DB
halfedges.push_back({0, b, 10});
halfedgeToFace.push_back(2);
halfedgeNext.push_back(6);
// CB
halfedges.push_back({0, b, 1});
halfedgeToFace.push_back(3);
halfedgeNext.push_back(10);
// BD
halfedges.push_back({0, d, 8});
halfedgeToFace.push_back(3);
halfedgeNext.push_back(11);
// DC
halfedges.push_back({0, c, 4});
halfedgeToFace.push_back(3);
halfedgeNext.push_back(9);
// Create faces
faces.emplace_back(0);
faces.emplace_back(3);
faces.emplace_back(6);
faces.emplace_back(9);
}
std::array<int, 3> MeshBuilder::getVertexIndicesOfFace(const Face& f) const {
std::array<int, 3> v;
size_t index = f.he;
auto* he = &halfedges[index];
v[0] = he->endVert;
index = halfedgeNext[index];
he = &halfedges[index];
v[1] = he->endVert;
index = halfedgeNext[index];
he = &halfedges[index];
v[2] = he->endVert;
return v;
}
HalfEdgeMesh::HalfEdgeMesh(const MeshBuilder& builderObject,
const VecView<vec3>& vertexData) {
std::unordered_map<size_t, size_t> faceMapping;
std::unordered_map<size_t, size_t> halfEdgeMapping;
std::unordered_map<size_t, size_t> vertexMapping;
size_t i = 0;
for (const auto& face : builderObject.faces) {
if (!face.isDisabled()) {
halfEdgeIndexFaces.emplace_back(static_cast<size_t>(face.he));
faceMapping[i] = halfEdgeIndexFaces.size() - 1;
const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face);
for (const auto heIndex : heIndices) {
const auto vertexIndex = builderObject.halfedges[heIndex].endVert;
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.pairedHalfedge != -1) {
halfedges.push_back({halfEdge.endVert, halfEdge.pairedHalfedge,
builderObject.halfedgeToFace[i]});
halfedgeToFace.push_back(builderObject.halfedgeToFace[i]);
halfedgeNext.push_back(builderObject.halfedgeNext[i]);
halfEdgeMapping[i] = halfedges.size() - 1;
}
i++;
}
for (auto& halfEdgeIndexFace : halfEdgeIndexFaces) {
DEBUG_ASSERT(halfEdgeMapping.count(halfEdgeIndexFace) == 1, logicErr,
"invalid halfedge mapping");
halfEdgeIndexFace = halfEdgeMapping[halfEdgeIndexFace];
}
for (size_t i = 0; i < halfedges.size(); i++) {
auto& he = halfedges[i];
halfedgeToFace[i] = faceMapping[halfedgeToFace[i]];
he.pairedHalfedge = halfEdgeMapping[he.pairedHalfedge];
halfedgeNext[i] = halfEdgeMapping[halfedgeNext[i]];
he.endVert = vertexMapping[he.endVert];
}
}
/*
* Implementation of the algorithm
*/
std::pair<Vec<Halfedge>, Vec<vec3>> QuickHull::buildMesh(double epsilon) {
if (originalVertexData.size() == 0) {
return {Vec<Halfedge>(), Vec<vec3>()};
}
// Very first: find extreme values and use them to compute the scale of the
// point cloud.
extremeValues = getExtremeValues();
scale = getScale(extremeValues);
// Epsilon we use depends on the scale
m_epsilon = epsilon * scale;
epsilonSquared = m_epsilon * m_epsilon;
// The planar case happens when all the points appear to lie on a two
// dimensional subspace of R^3.
planar = false;
createConvexHalfedgeMesh();
if (planar) {
const int extraPointIndex = planarPointCloudTemp.size() - 1;
for (auto& he : mesh.halfedges) {
if (he.endVert == extraPointIndex) {
he.endVert = 0;
}
}
planarPointCloudTemp.clear();
}
// reorder halfedges
Vec<Halfedge> halfedges(mesh.halfedges.size());
Vec<int> halfedgeToFace(mesh.halfedges.size());
Vec<int> counts(mesh.halfedges.size(), 0);
Vec<int> mapping(mesh.halfedges.size());
Vec<int> faceMap(mesh.faces.size());
// Some faces are disabled and should not go into the halfedge vector, we can
// update the face indices of the halfedges at the end using index/3
int j = 0;
for_each(
autoPolicy(mesh.halfedges.size()), countAt(0_uz),
countAt(mesh.halfedges.size()), [&](size_t i) {
if (mesh.halfedges[i].pairedHalfedge < 0) return;
if (mesh.faces[mesh.halfedgeToFace[i]].isDisabled()) return;
if (AtomicAdd(counts[mesh.halfedgeToFace[i]], 1) > 0) return;
int currIndex = AtomicAdd(j, 3);
mapping[i] = currIndex;
halfedges[currIndex + 0] = mesh.halfedges[i];
halfedgeToFace[currIndex + 0] = mesh.halfedgeToFace[i];
size_t k = mesh.halfedgeNext[i];
mapping[k] = currIndex + 1;
halfedges[currIndex + 1] = mesh.halfedges[k];
halfedgeToFace[currIndex + 1] = mesh.halfedgeToFace[k];
k = mesh.halfedgeNext[k];
mapping[k] = currIndex + 2;
halfedges[currIndex + 2] = mesh.halfedges[k];
halfedgeToFace[currIndex + 2] = mesh.halfedgeToFace[k];
halfedges[currIndex + 0].startVert = halfedges[currIndex + 2].endVert;
halfedges[currIndex + 1].startVert = halfedges[currIndex + 0].endVert;
halfedges[currIndex + 2].startVert = halfedges[currIndex + 1].endVert;
});
halfedges.resize(j);
halfedgeToFace.resize(j);
// fix pairedHalfedge id
for_each(
autoPolicy(halfedges.size()), halfedges.begin(), halfedges.end(),
[&](Halfedge& he) { he.pairedHalfedge = mapping[he.pairedHalfedge]; });
counts.resize_nofill(originalVertexData.size() + 1);
fill(counts.begin(), counts.end(), 0);
// remove unused vertices
for_each(autoPolicy(halfedges.size() / 3), countAt(0_uz),
countAt(halfedges.size() / 3), [&](size_t i) {
AtomicAdd(counts[halfedges[3 * i].startVert], 1);
AtomicAdd(counts[halfedges[3 * i + 1].startVert], 1);
AtomicAdd(counts[halfedges[3 * i + 2].startVert], 1);
});
auto saturate = [](int c) { return c > 0 ? 1 : 0; };
exclusive_scan(TransformIterator(counts.begin(), saturate),
TransformIterator(counts.end(), saturate), counts.begin(), 0);
Vec<vec3> vertices(counts.back());
for_each(autoPolicy(originalVertexData.size()), countAt(0_uz),
countAt(originalVertexData.size()), [&](size_t i) {
if (counts[i + 1] - counts[i] > 0) {
vertices[counts[i]] = originalVertexData[i];
}
});
for_each(autoPolicy(halfedges.size()), halfedges.begin(), halfedges.end(),
[&](Halfedge& he) {
he.startVert = counts[he.startVert];
he.endVert = counts[he.endVert];
});
return {std::move(halfedges), std::move(vertices)};
}
void QuickHull::createConvexHalfedgeMesh() {
visibleFaces.clear();
horizonEdgesData.clear();
possiblyVisibleFaces.clear();
// Compute base tetrahedron
setupInitialTetrahedron();
DEBUG_ASSERT(mesh.faces.size() == 4, logicErr, "not a tetrahedron");
// Init face stack with those faces that have points assigned to them
faceList.clear();
for (size_t i = 0; i < 4; i++) {
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 (!faceList.empty()) {
iter++;
if (iter == std::numeric_limits<size_t>::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 auto topFaceIndex = faceList.front();
faceList.pop_front();
auto& tf = mesh.faces[topFaceIndex];
tf.inFaceStack = 0;
DEBUG_ASSERT(
!tf.pointsOnPositiveSide || tf.pointsOnPositiveSide->size() > 0,
logicErr, "there should be points on the positive side");
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 = 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.
horizonEdgesData.clear();
possiblyVisibleFaces.clear();
visibleFaces.clear();
possiblyVisibleFaces.push_back({topFaceIndex, -1});
while (possiblyVisibleFaces.size()) {
const auto faceData = possiblyVisibleFaces.back();
possiblyVisibleFaces.pop_back();
auto& pvf = mesh.faces[faceData.faceIndex];
DEBUG_ASSERT(!pvf.isDisabled(), logicErr, "pvf should not be disabled");
if (pvf.visibilityCheckedOnIteration == iter) {
if (pvf.isVisibleFaceOnCurrentIteration) {
continue;
}
} else {
const Plane& P = pvf.P;
pvf.visibilityCheckedOnIteration = iter;
const double d = la::dot(P.N, activePoint) + P.D;
if (d > 0) {
pvf.isVisibleFaceOnCurrentIteration = 1;
pvf.horizonEdgesOnCurrentIteration = 0;
visibleFaces.push_back(faceData.faceIndex);
for (auto heIndex : mesh.getHalfEdgeIndicesOfFace(pvf)) {
if (mesh.halfedges[heIndex].pairedHalfedge !=
faceData.enteredFromHalfedge) {
possiblyVisibleFaces.push_back(
{mesh.halfedgeToFace[mesh.halfedges[heIndex].pairedHalfedge],
heIndex});
}
}
continue;
}
DEBUG_ASSERT(faceData.faceIndex != topFaceIndex, logicErr,
"face index invalid");
}
// The face is not visible. Therefore, the halfedge we came from is part
// of the horizon edge.
pvf.isVisibleFaceOnCurrentIteration = 0;
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 halfEdgesMesh = mesh.getHalfEdgeIndicesOfFace(
mesh.faces[mesh.halfedgeToFace[faceData.enteredFromHalfedge]]);
const std::int8_t ind =
(halfEdgesMesh[0] == faceData.enteredFromHalfedge)
? 0
: (halfEdgesMesh[1] == faceData.enteredFromHalfedge ? 1 : 2);
mesh.faces[mesh.halfedgeToFace[faceData.enteredFromHalfedge]]
.horizonEdgesOnCurrentIteration |= (1 << ind);
}
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(horizonEdgesData)) {
failedHorizonEdges++;
int change_flag = 0;
for (size_t index = 0; index < tf.pointsOnPositiveSide->size(); index++) {
if ((*tf.pointsOnPositiveSide)[index] == activePointIndex) {
change_flag = 1;
} else if (change_flag == 1) {
change_flag = 2;
(*tf.pointsOnPositiveSide)[index - 1] =
(*tf.pointsOnPositiveSide)[index];
}
}
if (change_flag == 1)
tf.pointsOnPositiveSide->resize(tf.pointsOnPositiveSide->size() - 1);
if (tf.pointsOnPositiveSide->size() == 0) {
reclaimToIndexVectorPool(tf.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.
newFaceIndices.clear();
newHalfedgeIndices.clear();
disabledFacePointVectors.clear();
size_t disableCounter = 0;
for (auto faceIndex : visibleFaces) {
auto& disabledFace = mesh.faces[faceIndex];
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(halfEdgesMesh[j]);
disableCounter++;
} else {
// Mark for reusal on later iteration step
mesh.disableHalfedge(halfEdgesMesh[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 = mesh.disableFace(faceIndex);
if (t) {
// Because we should not assign point vectors to faces unless needed...
DEBUG_ASSERT(t->size(), logicErr, "t should not be empty");
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++) {
newHalfedgeIndices.push_back(mesh.addHalfedge());
}
}
// Create new faces using the edgeloop
for (size_t i = 0; i < horizonEdgeCount; i++) {
const size_t AB = horizonEdgesData[i];
auto horizonEdgeVertexIndices =
mesh.getVertexIndicesOfHalfEdge(mesh.halfedges[AB]);
size_t A, B, C;
A = horizonEdgeVertexIndices[0];
B = horizonEdgeVertexIndices[1];
C = activePointIndex;
const size_t newFaceIndex = mesh.addFace();
newFaceIndices.push_back(newFaceIndex);
const size_t CA = newHalfedgeIndices[2 * i + 0];
const size_t BC = newHalfedgeIndices[2 * i + 1];
mesh.halfedgeNext[AB] = BC;
mesh.halfedgeNext[BC] = CA;
mesh.halfedgeNext[CA] = AB;
mesh.halfedgeToFace[BC] = newFaceIndex;
mesh.halfedgeToFace[CA] = newFaceIndex;
mesh.halfedgeToFace[AB] = newFaceIndex;
mesh.halfedges[CA].endVert = A;
mesh.halfedges[BC].endVert = C;
auto& newFace = mesh.faces[newFaceIndex];
const vec3 planeNormal = getTriangleNormal(
originalVertexData[A], originalVertexData[B], activePoint);
newFace.P = Plane(planeNormal, activePoint);
newFace.he = AB;
mesh.halfedges[CA].pairedHalfedge =
newHalfedgeIndices[i > 0 ? i * 2 - 1 : 2 * horizonEdgeCount - 1];
mesh.halfedges[BC].pairedHalfedge =
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 : disabledFacePointVectors) {
DEBUG_ASSERT(disabledPoints != nullptr, logicErr,
"disabledPoints should not be null");
for (const auto& point : *(disabledPoints)) {
if (point == activePointIndex) {
continue;
}
for (size_t j = 0; j < horizonEdgeCount; j++) {
if (addPointToFace(mesh.faces[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 : newFaceIndices) {
auto& newFace = mesh.faces[newFaceIndex];
if (newFace.pointsOnPositiveSide) {
DEBUG_ASSERT(newFace.pointsOnPositiveSide->size() > 0, logicErr,
"there should be points on the positive side");
if (!newFace.inFaceStack) {
faceList.push_back(newFaceIndex);
newFace.inFaceStack = 1;
}
}
}
}
// Cleanup
indexVectorPool.clear();
}
/*
* Private helper functions
*/
std::array<size_t, 6> QuickHull::getExtremeValues() {
std::array<size_t, 6> outIndices{0, 0, 0, 0, 0, 0};
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 vec3& pos = originalVertexData[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;
}
bool QuickHull::reorderHorizonEdges(VecView<size_t>& horizonEdges) {
const size_t horizonEdgeCount = horizonEdges.size();
for (size_t i = 0; i + 1 < horizonEdgeCount; i++) {
const size_t endVertexCheck = mesh.halfedges[horizonEdges[i]].endVert;
bool foundNext = false;
for (size_t j = i + 1; j < horizonEdgeCount; j++) {
const size_t beginVertex =
mesh.halfedges[mesh.halfedges[horizonEdges[j]].pairedHalfedge]
.endVert;
if (beginVertex == endVertexCheck) {
std::swap(horizonEdges[i + 1], horizonEdges[j]);
foundNext = true;
break;
}
}
if (!foundNext) {
return false;
}
}
DEBUG_ASSERT(
mesh.halfedges[horizonEdges[horizonEdges.size() - 1]].endVert ==
mesh.halfedges[mesh.halfedges[horizonEdges[0]].pairedHalfedge]
.endVert,
logicErr, "invalid halfedge");
return true;
}
double QuickHull::getScale(const std::array<size_t, 6>& extremeValuesInput) {
double s = 0;
for (size_t i = 0; i < 6; i++) {
const double* v =
(const double*)(&originalVertexData[extremeValuesInput[i]]);
v += i / 2;
auto a = std::abs(*v);
if (a > s) {
s = a;
}
}
return s;
}
void QuickHull::setupInitialTetrahedron() {
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 vec3 N =
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 mesh.setup(v[0], v[1], v[2], v[3]);
}
// Find two most distant extreme points.
double maxD = epsilonSquared;
std::pair<size_t, size_t> 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 = getSquaredDistance(originalVertexData[extremeValues[i]],
originalVertexData[extremeValues[j]]);
if (d > maxD) {
maxD = d;
selectedPoints = {extremeValues[i], extremeValues[j]};
}
}
}
if (maxD == epsilonSquared) {
// A degenerate case: the point cloud seems to consists of a single point
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));
}
DEBUG_ASSERT(selectedPoints.first != selectedPoints.second, logicErr,
"degenerate selectedPoints");
// Find the most distant point to the line between the two chosen extreme
// points.
const Ray r(originalVertexData[selectedPoints.first],
(originalVertexData[selectedPoints.second] -
originalVertexData[selectedPoints.first]));
maxD = epsilonSquared;
size_t maxI = std::numeric_limits<size_t>::max();
const size_t vCount = originalVertexData.size();
for (size_t i = 0; i < vCount; i++) {
const double distToRay =
getSquaredDistanceBetweenPointAndRay(originalVertexData[i], r);
if (distToRay > maxD) {
maxD = distToRay;
maxI = i;
}
}
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(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.
DEBUG_ASSERT(selectedPoints.first != maxI && selectedPoints.second != maxI,
logicErr, "degenerate selectedPoints");
std::array<size_t, 3> baseTriangle{selectedPoints.first,
selectedPoints.second, maxI};
const vec3 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;
maxI = 0;
const vec3 N =
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(
getSignedDistanceToPlane(originalVertexData[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.
planar = true;
const vec3 N1 =
getTriangleNormal(baseTriangleVertices[1], baseTriangleVertices[2],
baseTriangleVertices[0]);
planarPointCloudTemp = Vec<vec3>(originalVertexData);
const vec3 extraPoint = N1 + originalVertexData[0];
planarPointCloudTemp.push_back(extraPoint);
maxI = planarPointCloudTemp.size() - 1;
originalVertexData = 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(originalVertexData[maxI])) {
std::swap(baseTriangle[0], baseTriangle[1]);
}
// Create a tetrahedron half edge mesh and compute planes defined by each
// triangle
mesh.setup(baseTriangle[0], baseTriangle[1], baseTriangle[2], maxI);
for (auto& f : mesh.faces) {
auto v = mesh.getVertexIndicesOfFace(f);
const vec3 N1 =
getTriangleNormal(originalVertexData[v[0]], originalVertexData[v[1]],
originalVertexData[v[2]]);
const Plane plane(N1, originalVertexData[v[0]]);
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 : mesh.faces) {
if (addPointToFace(face, i)) {
break;
}
}
}
}
std::unique_ptr<Vec<size_t>> QuickHull::getIndexVectorFromPool() {
auto r = indexVectorPool.get();
r->clear();
return r;
}
void QuickHull::reclaimToIndexVectorPool(std::unique_ptr<Vec<size_t>>& 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 =
getSignedDistanceToPlane(originalVertexData[pointIndex], f.P);
if (D > 0 && D * D > epsilonSquared * f.P.sqrNLength) {
if (!f.pointsOnPositiveSide) {
f.pointsOnPositiveSide = getIndexVectorFromPool();
}
f.pointsOnPositiveSide->push_back(pointIndex);
if (D > f.mostDistantPointDist) {
f.mostDistantPointDist = D;
f.mostDistantPoint = pointIndex;
}
return true;
}
return false;
}
// Wrapper to call the QuickHull algorithm with the given vertex data to build
// the Impl
void Manifold::Impl::Hull(VecView<vec3> vertPos) {
size_t numVert = vertPos.size();
if (numVert < 4) {
status_ = Error::InvalidConstruction;
return;
}
QuickHull qh(vertPos);
std::tie(halfedge_, vertPos_) = qh.buildMesh();
CalculateBBox();
SetEpsilon();
InitializeOriginal();
Finish();
MarkCoplanar();
}
} // namespace manifold