24 typedef VectorType<unsigned,
D - 1> EdgeType;
26 SmartPointer<Mesh<T>> mesh =
nullptr;
27 SmartPointer<PointCloud<T, D>> pointCloud =
nullptr;
28 std::vector<EdgeType> visitedEdges;
29 std::vector<VectorType<unsigned, D>> hullElements;
30 std::list<EdgeType> remainingEdges;
35 unsigned pivotEdge(EdgeType currentEdge)
const {
36 auto &points = pointCloud->points;
37 unsigned nextIndex = 0;
40 while (nextIndex == currentEdge[0] ||
41 (
D == 3 && nextIndex == currentEdge[1])) {
45 for (
unsigned i = 0; i < points.size(); ++i) {
46 if (i == currentEdge[0] || i == nextIndex)
48 if (
D == 3 && i == currentEdge[1])
52 VectorType<T, D> normal;
53 VectorType<T, D> distance = points[i] - points[currentEdge[0]];
55 if constexpr (
D == 2) {
57 normal = points[nextIndex] - points[currentEdge[0]];
59 std::swap(normal[0], normal[1]);
60 normal[0] = -normal[0];
62 }
else if constexpr (
D == 3) {
63 auto v1 = points[currentEdge[1]] - points[currentEdge[0]];
64 auto v2 = points[nextIndex] - points[currentEdge[0]];
65 normal = CrossProduct(v1, v2);
69 auto product = DotProduct(distance, normal);
73 if (std::abs(product) < 1e-9) {
76 VectorType<unsigned, D> triangle;
77 triangle[0] = currentEdge[0];
78 triangle[1] = currentEdge[1];
80 if (doesTriangleClip(triangle)) {
85 edges[0][0] = currentEdge[0];
87 edges[1][0] = currentEdge[1];
90 if (!wasEdgeVisited(edges[0]) && !wasEdgeVisited(edges[1])) {
95 else if (product > 0) {
103 bool wasEdgeVisited(
const EdgeType &edge)
const {
104 for (
unsigned i = 0; i < visitedEdges.size(); ++i) {
106 if (edge[0] == visitedEdges[i][0]) {
107 if constexpr (
D == 3) {
108 if (edge[1] == visitedEdges[i][1]) {
117 if constexpr (
D == 3) {
118 if (edge[0] == visitedEdges[i][1] && edge[1] == visitedEdges[i][0]) {
127 typename std::list<EdgeType>::iterator findEdge(EdgeType edge) {
128 for (
typename std::list<EdgeType>::iterator it = remainingEdges.begin();
129 it != remainingEdges.end(); ++it) {
130 bool containsElem[2] = {
false,
false};
131 for (
unsigned i = 0; i <
D - 1; ++i) {
133 for (
unsigned k = 0; k <
D - 1; ++k) {
134 if ((*it)[i] == edge[k]) {
135 containsElem[i] =
true;
139 if (containsElem[0] && containsElem[1]) {
143 return remainingEdges.end();
147 bool intersectSharedNode(
const VectorType<unsigned, D> &triangle1,
148 const VectorType<unsigned, D> &triangle2)
const {
150 std::vector<unsigned> sharedNodes;
152 std::vector<unsigned> otherNodes;
154 for (
unsigned i = 0; i <
D; ++i) {
155 unsigned shared = sharedNodes.size();
156 for (
unsigned j = 0; j <
D; ++j) {
157 if (triangle1[i] == triangle2[j]) {
159 sharedNodes.push_back(j);
163 if (shared == sharedNodes.size()) {
164 otherNodes.push_back(i);
168 assert(!sharedNodes.empty());
170 auto &points = pointCloud->points;
178 for (
unsigned int &sharedNode : sharedNodes) {
179 VectorType<T, D> averageVector;
181 double centerEdgeDot;
183 auto shared = sharedNode;
184 VectorType<T, D> edge1 = Normalize(points[triangle2[(shared + 1) %
D]] -
185 points[triangle2[shared]]);
186 VectorType<T, D> edge2 = Normalize(points[triangle2[(shared + 2) %
D]] -
187 points[triangle2[shared]]);
188 averageVector = edge1 + edge2;
189 averageVector = averageVector /
T(2);
190 Normalize(averageVector);
191 centerEdgeDot = DotProduct(averageVector, edge2);
195 for (
unsigned int &otherNode : otherNodes) {
196 VectorType<T, D> vector = Normalize(points[triangle1[otherNode]] -
197 points[triangle2[sharedNode]]);
199 if (DotProduct(vector, averageVector) > centerEdgeDot + 1e-9) {
209 bool doesTriangleClip(
const VectorType<unsigned, D> &triangle)
const {
210 auto &points = pointCloud->points;
212 auto triangleNormal =
213 calculateNormal(points[triangle[1]] - points[triangle[0]],
214 points[triangle[2]] - points[triangle[0]]);
217 for (
unsigned i = 0; i < hullElements.size(); ++i) {
218 unsigned inOneTriangle = 0;
219 for (
unsigned j = 0; j <
D; ++j) {
221 for (
unsigned k = 0; k <
D; ++k) {
222 if (hullElements[i][j] == triangle[k]) {
229 if (inOneTriangle > 0) {
231 auto normal2 = calculateNormal(
232 points[hullElements[i][1]] - points[hullElements[i][0]],
233 points[hullElements[i][2]] - points[hullElements[i][0]]);
236 for (
unsigned d = 0; d <
D; ++d) {
238 std::abs(std::abs(triangleNormal[d]) - std::abs(normal2[d]));
246 if (intersectSharedNode(triangle, hullElements[i]))
248 if (intersectSharedNode(hullElements[i], triangle))
256 VectorType<T, D> calculateNormal(
const VectorType<T, D> &v1,
257 const VectorType<T, D> &v2)
const {
258 VectorType<T, D> newNormal;
259 newNormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
260 newNormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
261 newNormal[2] = v1[0] * v2[1] - v1[1] * v2[0];
277 : mesh(passedMesh), pointCloud(passedPointCloud) {}
282 pointCloud = passedPointCloud;
286 if (mesh ==
nullptr) {
287 Logger::getInstance()
288 .addWarning(
"No mesh was passed to ConvexHull.")
292 if (pointCloud ==
nullptr) {
293 Logger::getInstance()
294 .addWarning(
"No point cloud was passed to ConvexHull.")
300 auto &points = pointCloud->points;
301 if (points.size() == 0)
306 auto it = std::min_element(points.begin(), points.end());
307 unsigned currentIndex = std::distance(points.begin(), it);
310 remainingEdges.push_back(EdgeType{currentIndex});
313 int currentPointIndex = -1;
314 double hullMetric = std::numeric_limits<double>::max();
315 for (
int i = 0; i < points.size(); ++i) {
316 if (i == currentIndex)
318 if (std::abs(points[i][2] - points[currentIndex][2]) < 1e-7) {
319 auto diff = points[currentIndex] - points[i];
321 if (
const double currentHullMetric = DotProduct(diff, diff);
322 currentHullMetric < hullMetric) {
323 hullMetric = currentHullMetric;
324 currentPointIndex = i;
330 edge[0] = currentIndex;
332 if (currentPointIndex == -1) {
333 VectorType<T, D> newPoint;
334 newPoint = points[currentIndex];
337 points.push_back(newPoint);
338 edge[1] = points.size() - 1;
339 edge[1] = pivotEdge(edge);
342 edge[1] = currentPointIndex;
345 remainingEdges.push_back(edge);
350 while (!remainingEdges.empty()) {
352 EdgeType currentEdge = remainingEdges.back();
353 remainingEdges.pop_back();
355 if (wasEdgeVisited(currentEdge))
359 unsigned nextIndex = pivotEdge(currentEdge);
362 VectorType<unsigned, D> element;
363 for (
unsigned i = 0; i <
D - 1; ++i) {
364 element[i] = currentEdge[i];
366 element[
D - 1] = nextIndex;
367 hullElements.push_back(element);
370 visitedEdges.push_back(currentEdge);
372 remainingEdges.push_back(EdgeType{nextIndex});
376 edge[1] = currentEdge[1];
380 auto it = findEdge(edge);
381 if (it != remainingEdges.end()) {
382 remainingEdges.erase(it);
384 remainingEdges.push_back(edge);
387 edge[0] = currentEdge[0];
390 if (it != remainingEdges.end()) {
391 remainingEdges.erase(it);
393 remainingEdges.push_back(edge);
399 auto &elements = mesh->template getElements<D>();
400 std::unordered_map<unsigned, unsigned> oldToNewNodes;
401 for (
unsigned i = 0; i < hullElements.size(); ++i) {
403 std::array<unsigned, D> newElement;
404 for (
unsigned j = 0; j <
D; ++j) {
405 auto insertion = oldToNewNodes.insert(
406 std::make_pair(hullElements[i][j], oldToNewNodes.size()));
409 if (insertion.second) {
410 Vec3D<T> node{points[hullElements[i][j]][0],
411 points[hullElements[i][j]][1],
412 (
D == 2) ?
T(0.) : points[hullElements[i][j]][2]};
413 mesh->nodes.push_back(node);
416 newElement[j] = insertion.first->second;
419 elements.push_back(newElement);