ViennaLS
Loading...
Searching...
No Matches
lsConvexHull.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <algorithm>
6#include <list>
7#include <unordered_map>
8
9#include <lsGeometries.hpp>
10#include <lsMesh.hpp>
11
12#include <vcLogger.hpp>
13#include <vcSmartPointer.hpp>
14#include <vcVectorType.hpp>
15
16namespace viennals {
17
18using namespace viennacore;
19
23template <class T, int D> class ConvexHull {
24 typedef VectorType<unsigned, D - 1> EdgeType;
25
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;
31
32 // go through all points in the set and find the correct element
33 // to do this find the normal vector and dot product to see which side
34 // a point is on
35 unsigned pivotEdge(EdgeType currentEdge) const {
36 auto &points = pointCloud->points;
37 unsigned nextIndex = 0;
38
39 // start from index which is neither one of the edges
40 while (nextIndex == currentEdge[0] ||
41 (D == 3 && nextIndex == currentEdge[1])) {
42 ++nextIndex;
43 }
44
45 for (unsigned i = 0; i < points.size(); ++i) {
46 if (i == currentEdge[0] || i == nextIndex)
47 continue;
48 if (D == 3 && i == currentEdge[1])
49 continue;
50
51 // surface element normal and distance to point i
52 VectorType<T, D> normal;
53 VectorType<T, D> distance = points[i] - points[currentEdge[0]];
54 // create surface element and check if all points are to its right
55 if constexpr (D == 2) {
56 // in 2D normal is a 90 degree rotation
57 normal = points[nextIndex] - points[currentEdge[0]];
58 if (D == 2) {
59 std::swap(normal[0], normal[1]);
60 normal[0] = -normal[0];
61 }
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);
66 Normalize(normal);
67 }
68
69 auto product = DotProduct(distance, normal);
70 // if dot product is very small, point is very close to plane
71 // we need to check if we already have the correct point
72 // or if it is the next correct one
73 if (std::abs(product) < 1e-9) {
74 // check if suggested triangle intersects with any other
75 // in the same plane
76 VectorType<unsigned, D> triangle;
77 triangle[0] = currentEdge[0];
78 triangle[1] = currentEdge[1];
79 triangle[2] = i;
80 if (doesTriangleClip(triangle)) {
81 continue;
82 }
83
84 EdgeType edges[2];
85 edges[0][0] = currentEdge[0];
86 edges[0][1] = i;
87 edges[1][0] = currentEdge[1];
88 edges[1][1] = i;
89
90 if (!wasEdgeVisited(edges[0]) && !wasEdgeVisited(edges[1])) {
91 nextIndex = i;
92 }
93 }
94 // check if point is to the right of current element
95 else if (product > 0) {
96 nextIndex = i;
97 }
98 }
99 return nextIndex;
100 }
101
102 // check if edge was already visited
103 bool wasEdgeVisited(const EdgeType &edge) const {
104 for (unsigned i = 0; i < visitedEdges.size(); ++i) {
105 // in 2D one match is enough, in 3D both points must be the same
106 if (edge[0] == visitedEdges[i][0]) {
107 if constexpr (D == 3) {
108 if (edge[1] == visitedEdges[i][1]) {
109 return true;
110 }
111 } else {
112 return true;
113 }
114 }
115
116 // in 3D, an edge is also the same if they contain the same nodes
117 if constexpr (D == 3) {
118 if (edge[0] == visitedEdges[i][1] && edge[1] == visitedEdges[i][0]) {
119 return true;
120 }
121 }
122 }
123 return false;
124 }
125
126 // find whether an edge has already been used for a triangle once
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) {
132 // check all possible permutations of elements
133 for (unsigned k = 0; k < D - 1; ++k) {
134 if ((*it)[i] == edge[k]) {
135 containsElem[i] = true;
136 }
137 }
138 }
139 if (containsElem[0] && containsElem[1]) {
140 return it;
141 }
142 }
143 return remainingEdges.end();
144 }
145
146 // return whether two triangles with shared nodes intersect
147 bool intersectSharedNode(const VectorType<unsigned, D> &triangle1,
148 const VectorType<unsigned, D> &triangle2) const {
149 // find which nodes are shared and which nodes are others
150 std::vector<unsigned> sharedNodes;
151 // otherNodes2 contains the nodes of triangle1 which are not shared
152 std::vector<unsigned> otherNodes;
153
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]) {
158 // node numbers of t2
159 sharedNodes.push_back(j);
160 }
161 }
162 // if triangle1[i] was not shared, push to other nodes
163 if (shared == sharedNodes.size()) {
164 otherNodes.push_back(i);
165 }
166 }
167
168 assert(!sharedNodes.empty());
169
170 auto &points = pointCloud->points;
171
172 // for each shared node, calculate the edges of t2 leading away from it.
173 // then generate the average vector and its angle to the edges.
174 // then calculate the dot product of each other node of t1 to the average
175 // vector. If it is smaller than the dot product between the average vector
176 // and an edge, the node must be outside of triangle. If it is inside, one
177 // triangle must clip the other.
178 for (unsigned int &sharedNode : sharedNodes) {
179 VectorType<T, D> averageVector;
180 // save the dot product between average and one edge
181 double centerEdgeDot;
182 {
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);
192 }
193
194 // go over other nodes of triangle1
195 for (unsigned int &otherNode : otherNodes) {
196 VectorType<T, D> vector = Normalize(points[triangle1[otherNode]] -
197 points[triangle2[sharedNode]]);
198
199 if (DotProduct(vector, averageVector) > centerEdgeDot + 1e-9) {
200 return true;
201 }
202 }
203 }
204 // if no other node was ever within two edges, the triangles do not clip
205 return false;
206 }
207
208 // check if triangle defined by two edges clips any other triangle
209 bool doesTriangleClip(const VectorType<unsigned, D> &triangle) const {
210 auto &points = pointCloud->points;
211
212 auto triangleNormal =
213 calculateNormal(points[triangle[1]] - points[triangle[0]],
214 points[triangle[2]] - points[triangle[0]]);
215
216 // unsigned shareEdges = 0;
217 for (unsigned i = 0; i < hullElements.size(); ++i) {
218 unsigned inOneTriangle = 0;
219 for (unsigned j = 0; j < D; ++j) {
220 // check all possible permutations of elements
221 for (unsigned k = 0; k < D; ++k) {
222 if (hullElements[i][j] == triangle[k]) {
223 ++inOneTriangle;
224 }
225 }
226 }
227
228 // if they share at least one node, they might clip, so check
229 if (inOneTriangle > 0) {
230 // check if they are in the same plane
231 auto normal2 = calculateNormal(
232 points[hullElements[i][1]] - points[hullElements[i][0]],
233 points[hullElements[i][2]] - points[hullElements[i][0]]);
234
235 bool skip = false;
236 for (unsigned d = 0; d < D; ++d) {
237 double diff =
238 std::abs(std::abs(triangleNormal[d]) - std::abs(normal2[d]));
239 if (diff > 1e-6)
240 skip = true;
241 }
242
243 if (skip)
244 continue;
245
246 if (intersectSharedNode(triangle, hullElements[i]))
247 return true;
248 if (intersectSharedNode(hullElements[i], triangle))
249 return true;
250 }
251 }
252 return false;
253 }
254
255 // calculate the normal vector of two vectors
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];
262 return newNormal;
263 }
264
265 // // calculate Area of triangle defined by 2 vectors
266 // T calculateArea(const VectorType<T, D> &v1,
267 // const VectorType<T, D> &v2) const {
268 // VectorType<T, D> newNormal = CrossProduct(v1, v2);
269 // return std::sqrt(DotProduct(newNormal, newNormal));
270 // }
271
272public:
273 ConvexHull() = default;
274
275 ConvexHull(SmartPointer<Mesh<T>> passedMesh,
276 SmartPointer<PointCloud<T, D>> passedPointCloud)
277 : mesh(passedMesh), pointCloud(passedPointCloud) {}
278
279 void setMesh(SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
280
281 void setPointCloud(SmartPointer<PointCloud<T, D>> passedPointCloud) {
282 pointCloud = passedPointCloud;
283 }
284
285 void apply() {
286 if (mesh == nullptr) {
287 Logger::getInstance()
288 .addWarning("No mesh was passed to ConvexHull.")
289 .print();
290 return;
291 }
292 if (pointCloud == nullptr) {
293 Logger::getInstance()
294 .addWarning("No point cloud was passed to ConvexHull.")
295 .print();
296 return;
297 }
298
299 mesh->clear();
300 auto &points = pointCloud->points;
301 if (points.size() == 0)
302 return;
303
304 {
305 // find first hull point
306 auto it = std::min_element(points.begin(), points.end());
307 unsigned currentIndex = std::distance(points.begin(), it);
308
309 if (D == 2) {
310 remainingEdges.push_back(EdgeType{currentIndex});
311 } else {
312 // need to check if there is second point at same z coord
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)
317 continue;
318 if (std::abs(points[i][2] - points[currentIndex][2]) < 1e-7) {
319 auto diff = points[currentIndex] - points[i];
320 // choose closest point if points are in z plane
321 if (const double currentHullMetric = DotProduct(diff, diff);
322 currentHullMetric < hullMetric) {
323 hullMetric = currentHullMetric;
324 currentPointIndex = i;
325 }
326 }
327 }
328
329 EdgeType edge;
330 edge[0] = currentIndex;
331 // if there was no point at same z, find through pivot
332 if (currentPointIndex == -1) {
333 VectorType<T, D> newPoint;
334 newPoint = points[currentIndex];
335 newPoint[0] += 1.0;
336 // take newPoint for fake edge to find correct first point
337 points.push_back(newPoint);
338 edge[1] = points.size() - 1;
339 edge[1] = pivotEdge(edge);
340 points.pop_back();
341 } else {
342 edge[1] = currentPointIndex;
343 }
344
345 remainingEdges.push_back(edge);
346 }
347 }
348
349 // until there are no more edges to check
350 while (!remainingEdges.empty()) {
351 // for all points in the cloud
352 EdgeType currentEdge = remainingEdges.back();
353 remainingEdges.pop_back();
354
355 if (wasEdgeVisited(currentEdge))
356 continue;
357
358 // find next point in hull
359 unsigned nextIndex = pivotEdge(currentEdge);
360
361 // set new hull element and save
362 VectorType<unsigned, D> element;
363 for (unsigned i = 0; i < D - 1; ++i) {
364 element[i] = currentEdge[i];
365 }
366 element[D - 1] = nextIndex;
367 hullElements.push_back(element);
368
369 // mark edge as visited
370 visitedEdges.push_back(currentEdge);
371 if (D == 2) {
372 remainingEdges.push_back(EdgeType{nextIndex});
373 } else {
374 EdgeType edge;
375 edge[0] = nextIndex;
376 edge[1] = currentEdge[1];
377 // if edge is already part of another element
378 // and we just added it, we need to remove it
379 // from the list
380 auto it = findEdge(edge);
381 if (it != remainingEdges.end()) {
382 remainingEdges.erase(it);
383 } else {
384 remainingEdges.push_back(edge);
385 }
386
387 edge[0] = currentEdge[0];
388 edge[1] = nextIndex;
389 it = findEdge(edge);
390 if (it != remainingEdges.end()) {
391 remainingEdges.erase(it);
392 } else {
393 remainingEdges.push_back(edge);
394 }
395 }
396 }
397
398 // now make the mesh
399 auto &elements = mesh->template getElements<D>();
400 std::unordered_map<unsigned, unsigned> oldToNewNodes;
401 for (unsigned i = 0; i < hullElements.size(); ++i) {
402 // here add translation of old index to new index for new points
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()));
407
408 // if insertion took place, add the point to the nodes
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);
414 }
415
416 newElement[j] = insertion.first->second;
417 }
418
419 elements.push_back(newElement);
420 }
421 }
422};
423
424// add all template specialisations for this class
426
427} // namespace viennals
void apply()
Definition lsConvexHull.hpp:285
ConvexHull(SmartPointer< Mesh< T > > passedMesh, SmartPointer< PointCloud< T, D > > passedPointCloud)
Definition lsConvexHull.hpp:275
void setMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsConvexHull.hpp:279
void setPointCloud(SmartPointer< PointCloud< T, D > > passedPointCloud)
Definition lsConvexHull.hpp:281
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:21
Class describing a point cloud, which can be used to create geometries from its convex hull mesh.
Definition lsGeometries.hpp:128
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68