1/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com)
2 All rights reserved.
3
4
5 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6
7 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
8
9 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
10
11 3. The names of the contributors may not be used to endorse or promote products derived from this software without specific prior written permission.
12
13 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
14 */
15#include "vhacdICHull.h"
16#include <limits>
17
18#ifdef _MSC_VER
19#pragma warning(disable:4456 4706)
20#endif
21
22
23namespace VHACD {
24const double ICHull::sc_eps = 1.0e-15;
25const int32_t ICHull::sc_dummyIndex = std::numeric_limits<int32_t>::max();
26ICHull::ICHull()
27{
28 m_isFlat = false;
29}
30bool ICHull::AddPoints(const Vec3<double>* points, size_t nPoints)
31{
32 if (!points) {
33 return false;
34 }
35 CircularListElement<TMMVertex>* vertex = NULL;
36 for (size_t i = 0; i < nPoints; i++) {
37 vertex = m_mesh.AddVertex();
38 vertex->GetData().m_pos.X() = points[i].X();
39 vertex->GetData().m_pos.Y() = points[i].Y();
40 vertex->GetData().m_pos.Z() = points[i].Z();
41 vertex->GetData().m_name = static_cast<int32_t>(i);
42 }
43 return true;
44}
45bool ICHull::AddPoint(const Vec3<double>& point, int32_t id)
46{
47 if (AddPoints(&point, 1)) {
48 m_mesh.m_vertices.GetData().m_name = id;
49 return true;
50 }
51 return false;
52}
53
54ICHullError ICHull::Process()
55{
56 uint32_t addedPoints = 0;
57 if (m_mesh.GetNVertices() < 3) {
58 return ICHullErrorNotEnoughPoints;
59 }
60 if (m_mesh.GetNVertices() == 3) {
61 m_isFlat = true;
62 CircularListElement<TMMTriangle>* t1 = m_mesh.AddTriangle();
63 CircularListElement<TMMTriangle>* t2 = m_mesh.AddTriangle();
64 CircularListElement<TMMVertex>* v0 = m_mesh.m_vertices.GetHead();
65 CircularListElement<TMMVertex>* v1 = v0->GetNext();
66 CircularListElement<TMMVertex>* v2 = v1->GetNext();
67 // Compute the normal to the plane
68 Vec3<double> p0 = v0->GetData().m_pos;
69 Vec3<double> p1 = v1->GetData().m_pos;
70 Vec3<double> p2 = v2->GetData().m_pos;
71 m_normal = (p1 - p0) ^ (p2 - p0);
72 m_normal.Normalize();
73 t1->GetData().m_vertices[0] = v0;
74 t1->GetData().m_vertices[1] = v1;
75 t1->GetData().m_vertices[2] = v2;
76 t2->GetData().m_vertices[0] = v1;
77 t2->GetData().m_vertices[1] = v2;
78 t2->GetData().m_vertices[2] = v2;
79 return ICHullErrorOK;
80 }
81 if (m_isFlat) {
82 m_mesh.m_edges.Clear();
83 m_mesh.m_triangles.Clear();
84 m_isFlat = false;
85 }
86 if (m_mesh.GetNTriangles() == 0) // we have to create the first polyhedron
87 {
88 ICHullError res = DoubleTriangle();
89 if (res != ICHullErrorOK) {
90 return res;
91 }
92 else {
93 addedPoints += 3;
94 }
95 }
96 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
97 // go to the first added and not processed vertex
98 while (!(vertices.GetHead()->GetPrev()->GetData().m_tag)) {
99 vertices.Prev();
100 }
101 while (!vertices.GetData().m_tag) // not processed
102 {
103 vertices.GetData().m_tag = true;
104 if (ProcessPoint()) {
105 addedPoints++;
106 CleanUp(addedPoints);
107 vertices.Next();
108 if (!GetMesh().CheckConsistancy()) {
109 size_t nV = m_mesh.GetNVertices();
110 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
111 for (size_t v = 0; v < nV; ++v) {
112 if (vertices.GetData().m_name == sc_dummyIndex) {
113 vertices.Delete();
114 break;
115 }
116 vertices.Next();
117 }
118 return ICHullErrorInconsistent;
119 }
120 }
121 }
122 if (m_isFlat) {
123 SArray<CircularListElement<TMMTriangle>*> trianglesToDuplicate;
124 size_t nT = m_mesh.GetNTriangles();
125 for (size_t f = 0; f < nT; f++) {
126 TMMTriangle& currentTriangle = m_mesh.m_triangles.GetHead()->GetData();
127 if (currentTriangle.m_vertices[0]->GetData().m_name == sc_dummyIndex || currentTriangle.m_vertices[1]->GetData().m_name == sc_dummyIndex || currentTriangle.m_vertices[2]->GetData().m_name == sc_dummyIndex) {
128 m_trianglesToDelete.PushBack(m_mesh.m_triangles.GetHead());
129 for (int32_t k = 0; k < 3; k++) {
130 for (int32_t h = 0; h < 2; h++) {
131 if (currentTriangle.m_edges[k]->GetData().m_triangles[h] == m_mesh.m_triangles.GetHead()) {
132 currentTriangle.m_edges[k]->GetData().m_triangles[h] = 0;
133 break;
134 }
135 }
136 }
137 }
138 else {
139 trianglesToDuplicate.PushBack(m_mesh.m_triangles.GetHead());
140 }
141 m_mesh.m_triangles.Next();
142 }
143 size_t nE = m_mesh.GetNEdges();
144 for (size_t e = 0; e < nE; e++) {
145 TMMEdge& currentEdge = m_mesh.m_edges.GetHead()->GetData();
146 if (currentEdge.m_triangles[0] == 0 && currentEdge.m_triangles[1] == 0) {
147 m_edgesToDelete.PushBack(m_mesh.m_edges.GetHead());
148 }
149 m_mesh.m_edges.Next();
150 }
151 size_t nV = m_mesh.GetNVertices();
152 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
153 for (size_t v = 0; v < nV; ++v) {
154 if (vertices.GetData().m_name == sc_dummyIndex) {
155 vertices.Delete();
156 }
157 else {
158 vertices.GetData().m_tag = false;
159 vertices.Next();
160 }
161 }
162 CleanEdges();
163 CleanTriangles();
164 CircularListElement<TMMTriangle>* newTriangle;
165 for (size_t t = 0; t < trianglesToDuplicate.Size(); t++) {
166 newTriangle = m_mesh.AddTriangle();
167 newTriangle->GetData().m_vertices[0] = trianglesToDuplicate[t]->GetData().m_vertices[1];
168 newTriangle->GetData().m_vertices[1] = trianglesToDuplicate[t]->GetData().m_vertices[0];
169 newTriangle->GetData().m_vertices[2] = trianglesToDuplicate[t]->GetData().m_vertices[2];
170 }
171 }
172 return ICHullErrorOK;
173}
174ICHullError ICHull::Process(const uint32_t nPointsCH,
175 const double minVolume)
176{
177 uint32_t addedPoints = 0;
178 if (nPointsCH < 3 || m_mesh.GetNVertices() < 3) {
179 return ICHullErrorNotEnoughPoints;
180 }
181 if (m_mesh.GetNVertices() == 3) {
182 m_isFlat = true;
183 CircularListElement<TMMTriangle>* t1 = m_mesh.AddTriangle();
184 CircularListElement<TMMTriangle>* t2 = m_mesh.AddTriangle();
185 CircularListElement<TMMVertex>* v0 = m_mesh.m_vertices.GetHead();
186 CircularListElement<TMMVertex>* v1 = v0->GetNext();
187 CircularListElement<TMMVertex>* v2 = v1->GetNext();
188 // Compute the normal to the plane
189 Vec3<double> p0 = v0->GetData().m_pos;
190 Vec3<double> p1 = v1->GetData().m_pos;
191 Vec3<double> p2 = v2->GetData().m_pos;
192 m_normal = (p1 - p0) ^ (p2 - p0);
193 m_normal.Normalize();
194 t1->GetData().m_vertices[0] = v0;
195 t1->GetData().m_vertices[1] = v1;
196 t1->GetData().m_vertices[2] = v2;
197 t2->GetData().m_vertices[0] = v1;
198 t2->GetData().m_vertices[1] = v0;
199 t2->GetData().m_vertices[2] = v2;
200 return ICHullErrorOK;
201 }
202
203 if (m_isFlat) {
204 m_mesh.m_triangles.Clear();
205 m_mesh.m_edges.Clear();
206 m_isFlat = false;
207 }
208
209 if (m_mesh.GetNTriangles() == 0) // we have to create the first polyhedron
210 {
211 ICHullError res = DoubleTriangle();
212 if (res != ICHullErrorOK) {
213 return res;
214 }
215 else {
216 addedPoints += 3;
217 }
218 }
219 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
220 while (!vertices.GetData().m_tag && addedPoints < nPointsCH) // not processed
221 {
222 if (!FindMaxVolumePoint((addedPoints > 4) ? minVolume : 0.0)) {
223 break;
224 }
225 vertices.GetData().m_tag = true;
226 if (ProcessPoint()) {
227 addedPoints++;
228 CleanUp(addedPoints);
229 if (!GetMesh().CheckConsistancy()) {
230 size_t nV = m_mesh.GetNVertices();
231 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
232 for (size_t v = 0; v < nV; ++v) {
233 if (vertices.GetData().m_name == sc_dummyIndex) {
234 vertices.Delete();
235 break;
236 }
237 vertices.Next();
238 }
239 return ICHullErrorInconsistent;
240 }
241 vertices.Next();
242 }
243 }
244 // delete remaining points
245 while (!vertices.GetData().m_tag) {
246 vertices.Delete();
247 }
248 if (m_isFlat) {
249 SArray<CircularListElement<TMMTriangle>*> trianglesToDuplicate;
250 size_t nT = m_mesh.GetNTriangles();
251 for (size_t f = 0; f < nT; f++) {
252 TMMTriangle& currentTriangle = m_mesh.m_triangles.GetHead()->GetData();
253 if (currentTriangle.m_vertices[0]->GetData().m_name == sc_dummyIndex || currentTriangle.m_vertices[1]->GetData().m_name == sc_dummyIndex || currentTriangle.m_vertices[2]->GetData().m_name == sc_dummyIndex) {
254 m_trianglesToDelete.PushBack(m_mesh.m_triangles.GetHead());
255 for (int32_t k = 0; k < 3; k++) {
256 for (int32_t h = 0; h < 2; h++) {
257 if (currentTriangle.m_edges[k]->GetData().m_triangles[h] == m_mesh.m_triangles.GetHead()) {
258 currentTriangle.m_edges[k]->GetData().m_triangles[h] = 0;
259 break;
260 }
261 }
262 }
263 }
264 else {
265 trianglesToDuplicate.PushBack(m_mesh.m_triangles.GetHead());
266 }
267 m_mesh.m_triangles.Next();
268 }
269 size_t nE = m_mesh.GetNEdges();
270 for (size_t e = 0; e < nE; e++) {
271 TMMEdge& currentEdge = m_mesh.m_edges.GetHead()->GetData();
272 if (currentEdge.m_triangles[0] == 0 && currentEdge.m_triangles[1] == 0) {
273 m_edgesToDelete.PushBack(m_mesh.m_edges.GetHead());
274 }
275 m_mesh.m_edges.Next();
276 }
277 size_t nV = m_mesh.GetNVertices();
278 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
279 for (size_t v = 0; v < nV; ++v) {
280 if (vertices.GetData().m_name == sc_dummyIndex) {
281 vertices.Delete();
282 }
283 else {
284 vertices.GetData().m_tag = false;
285 vertices.Next();
286 }
287 }
288 CleanEdges();
289 CleanTriangles();
290 CircularListElement<TMMTriangle>* newTriangle;
291 for (size_t t = 0; t < trianglesToDuplicate.Size(); t++) {
292 newTriangle = m_mesh.AddTriangle();
293 newTriangle->GetData().m_vertices[0] = trianglesToDuplicate[t]->GetData().m_vertices[1];
294 newTriangle->GetData().m_vertices[1] = trianglesToDuplicate[t]->GetData().m_vertices[0];
295 newTriangle->GetData().m_vertices[2] = trianglesToDuplicate[t]->GetData().m_vertices[2];
296 }
297 }
298 return ICHullErrorOK;
299}
300bool ICHull::FindMaxVolumePoint(const double minVolume)
301{
302 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
303 CircularListElement<TMMVertex>* vMaxVolume = 0;
304 CircularListElement<TMMVertex>* vHeadPrev = vertices.GetHead()->GetPrev();
305
306 double maxVolume = minVolume;
307 double volume = 0.0;
308 while (!vertices.GetData().m_tag) // not processed
309 {
310 if (ComputePointVolume(volume, false)) {
311 if (maxVolume < volume) {
312 maxVolume = volume;
313 vMaxVolume = vertices.GetHead();
314 }
315 vertices.Next();
316 }
317 }
318 CircularListElement<TMMVertex>* vHead = vHeadPrev->GetNext();
319 vertices.GetHead() = vHead;
320 if (!vMaxVolume) {
321 return false;
322 }
323 if (vMaxVolume != vHead) {
324 Vec3<double> pos = vHead->GetData().m_pos;
325 int32_t id = vHead->GetData().m_name;
326 vHead->GetData().m_pos = vMaxVolume->GetData().m_pos;
327 vHead->GetData().m_name = vMaxVolume->GetData().m_name;
328 vMaxVolume->GetData().m_pos = pos;
329 vHead->GetData().m_name = id;
330 }
331 return true;
332}
333ICHullError ICHull::DoubleTriangle()
334{
335 // find three non colinear points
336 m_isFlat = false;
337 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
338 CircularListElement<TMMVertex>* v0 = vertices.GetHead();
339 while (Colinear(v0->GetData().m_pos,
340 v0->GetNext()->GetData().m_pos,
341 v0->GetNext()->GetNext()->GetData().m_pos)) {
342 if ((v0 = v0->GetNext()) == vertices.GetHead()) {
343 return ICHullErrorCoplanarPoints;
344 }
345 }
346 CircularListElement<TMMVertex>* v1 = v0->GetNext();
347 CircularListElement<TMMVertex>* v2 = v1->GetNext();
348 // mark points as processed
349 v0->GetData().m_tag = v1->GetData().m_tag = v2->GetData().m_tag = true;
350
351 // create two triangles
352 CircularListElement<TMMTriangle>* f0 = MakeFace(v0, v1, v2, 0);
353 MakeFace(v2, v1, v0, f0);
354
355 // find a fourth non-coplanar point to form tetrahedron
356 CircularListElement<TMMVertex>* v3 = v2->GetNext();
357 vertices.GetHead() = v3;
358
359 double vol = ComputeVolume4(v0->GetData().m_pos, v1->GetData().m_pos, v2->GetData().m_pos, v3->GetData().m_pos);
360 while (fabs(vol) < sc_eps && !v3->GetNext()->GetData().m_tag) {
361 v3 = v3->GetNext();
362 vol = ComputeVolume4(v0->GetData().m_pos, v1->GetData().m_pos, v2->GetData().m_pos, v3->GetData().m_pos);
363 }
364 if (fabs(vol) < sc_eps) {
365 // compute the barycenter
366 Vec3<double> bary(0.0, 0.0, 0.0);
367 CircularListElement<TMMVertex>* vBary = v0;
368 do {
369 bary += vBary->GetData().m_pos;
370 } while ((vBary = vBary->GetNext()) != v0);
371 bary /= static_cast<double>(vertices.GetSize());
372
373 // Compute the normal to the plane
374 Vec3<double> p0 = v0->GetData().m_pos;
375 Vec3<double> p1 = v1->GetData().m_pos;
376 Vec3<double> p2 = v2->GetData().m_pos;
377 m_normal = (p1 - p0) ^ (p2 - p0);
378 m_normal.Normalize();
379 // add dummy vertex placed at (bary + normal)
380 vertices.GetHead() = v2;
381 Vec3<double> newPt = bary + m_normal;
382 AddPoint(newPt, sc_dummyIndex);
383 m_isFlat = true;
384 return ICHullErrorOK;
385 }
386 else if (v3 != vertices.GetHead()) {
387 TMMVertex temp;
388 temp.m_name = v3->GetData().m_name;
389 temp.m_pos = v3->GetData().m_pos;
390 v3->GetData().m_name = vertices.GetHead()->GetData().m_name;
391 v3->GetData().m_pos = vertices.GetHead()->GetData().m_pos;
392 vertices.GetHead()->GetData().m_name = temp.m_name;
393 vertices.GetHead()->GetData().m_pos = temp.m_pos;
394 }
395 return ICHullErrorOK;
396}
397CircularListElement<TMMTriangle>* ICHull::MakeFace(CircularListElement<TMMVertex>* v0,
398 CircularListElement<TMMVertex>* v1,
399 CircularListElement<TMMVertex>* v2,
400 CircularListElement<TMMTriangle>* fold)
401{
402 CircularListElement<TMMEdge>* e0;
403 CircularListElement<TMMEdge>* e1;
404 CircularListElement<TMMEdge>* e2;
405 int32_t index = 0;
406 if (!fold) // if first face to be created
407 {
408 e0 = m_mesh.AddEdge(); // create the three edges
409 e1 = m_mesh.AddEdge();
410 e2 = m_mesh.AddEdge();
411 }
412 else // otherwise re-use existing edges (in reverse order)
413 {
414 e0 = fold->GetData().m_edges[2];
415 e1 = fold->GetData().m_edges[1];
416 e2 = fold->GetData().m_edges[0];
417 index = 1;
418 }
419 e0->GetData().m_vertices[0] = v0;
420 e0->GetData().m_vertices[1] = v1;
421 e1->GetData().m_vertices[0] = v1;
422 e1->GetData().m_vertices[1] = v2;
423 e2->GetData().m_vertices[0] = v2;
424 e2->GetData().m_vertices[1] = v0;
425 // create the new face
426 CircularListElement<TMMTriangle>* f = m_mesh.AddTriangle();
427 f->GetData().m_edges[0] = e0;
428 f->GetData().m_edges[1] = e1;
429 f->GetData().m_edges[2] = e2;
430 f->GetData().m_vertices[0] = v0;
431 f->GetData().m_vertices[1] = v1;
432 f->GetData().m_vertices[2] = v2;
433 // link edges to face f
434 e0->GetData().m_triangles[index] = e1->GetData().m_triangles[index] = e2->GetData().m_triangles[index] = f;
435 return f;
436}
437CircularListElement<TMMTriangle>* ICHull::MakeConeFace(CircularListElement<TMMEdge>* e, CircularListElement<TMMVertex>* p)
438{
439 // create two new edges if they don't already exist
440 CircularListElement<TMMEdge>* newEdges[2];
441 for (int32_t i = 0; i < 2; ++i) {
442 if (!(newEdges[i] = e->GetData().m_vertices[i]->GetData().m_duplicate)) { // if the edge doesn't exits add it and mark the vertex as duplicated
443 newEdges[i] = m_mesh.AddEdge();
444 newEdges[i]->GetData().m_vertices[0] = e->GetData().m_vertices[i];
445 newEdges[i]->GetData().m_vertices[1] = p;
446 e->GetData().m_vertices[i]->GetData().m_duplicate = newEdges[i];
447 }
448 }
449 // make the new face
450 CircularListElement<TMMTriangle>* newFace = m_mesh.AddTriangle();
451 newFace->GetData().m_edges[0] = e;
452 newFace->GetData().m_edges[1] = newEdges[0];
453 newFace->GetData().m_edges[2] = newEdges[1];
454 MakeCCW(newFace, e, p);
455 for (int32_t i = 0; i < 2; ++i) {
456 for (int32_t j = 0; j < 2; ++j) {
457 if (!newEdges[i]->GetData().m_triangles[j]) {
458 newEdges[i]->GetData().m_triangles[j] = newFace;
459 break;
460 }
461 }
462 }
463 return newFace;
464}
465bool ICHull::ComputePointVolume(double& totalVolume, bool markVisibleFaces)
466{
467 // mark visible faces
468 CircularListElement<TMMTriangle>* fHead = m_mesh.GetTriangles().GetHead();
469 CircularListElement<TMMTriangle>* f = fHead;
470 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
471 CircularListElement<TMMVertex>* vertex0 = vertices.GetHead();
472 bool visible = false;
473 Vec3<double> pos0 = Vec3<double>(vertex0->GetData().m_pos.X(),
474 vertex0->GetData().m_pos.Y(),
475 vertex0->GetData().m_pos.Z());
476 double vol = 0.0;
477 totalVolume = 0.0;
478 Vec3<double> ver0, ver1, ver2;
479 do {
480 ver0.X() = f->GetData().m_vertices[0]->GetData().m_pos.X();
481 ver0.Y() = f->GetData().m_vertices[0]->GetData().m_pos.Y();
482 ver0.Z() = f->GetData().m_vertices[0]->GetData().m_pos.Z();
483 ver1.X() = f->GetData().m_vertices[1]->GetData().m_pos.X();
484 ver1.Y() = f->GetData().m_vertices[1]->GetData().m_pos.Y();
485 ver1.Z() = f->GetData().m_vertices[1]->GetData().m_pos.Z();
486 ver2.X() = f->GetData().m_vertices[2]->GetData().m_pos.X();
487 ver2.Y() = f->GetData().m_vertices[2]->GetData().m_pos.Y();
488 ver2.Z() = f->GetData().m_vertices[2]->GetData().m_pos.Z();
489 vol = ComputeVolume4(ver0, ver1, ver2, pos0);
490 if (vol < -sc_eps) {
491 vol = fabs(vol);
492 totalVolume += vol;
493 if (markVisibleFaces) {
494 f->GetData().m_visible = true;
495 m_trianglesToDelete.PushBack(f);
496 }
497 visible = true;
498 }
499 f = f->GetNext();
500 } while (f != fHead);
501
502 if (m_trianglesToDelete.Size() == m_mesh.m_triangles.GetSize()) {
503 for (size_t i = 0; i < m_trianglesToDelete.Size(); i++) {
504 m_trianglesToDelete[i]->GetData().m_visible = false;
505 }
506 visible = false;
507 }
508 // if no faces visible from p then p is inside the hull
509 if (!visible && markVisibleFaces) {
510 vertices.Delete();
511 m_trianglesToDelete.Resize(0);
512 return false;
513 }
514 return true;
515}
516bool ICHull::ProcessPoint()
517{
518 double totalVolume = 0.0;
519 if (!ComputePointVolume(totalVolume, true)) {
520 return false;
521 }
522 // Mark edges in interior of visible region for deletion.
523 // Create a new face based on each border edge
524 CircularListElement<TMMVertex>* v0 = m_mesh.GetVertices().GetHead();
525 CircularListElement<TMMEdge>* eHead = m_mesh.GetEdges().GetHead();
526 CircularListElement<TMMEdge>* e = eHead;
527 CircularListElement<TMMEdge>* tmp = 0;
528 int32_t nvisible = 0;
529 m_edgesToDelete.Resize(0);
530 m_edgesToUpdate.Resize(0);
531 do {
532 tmp = e->GetNext();
533 nvisible = 0;
534 for (int32_t k = 0; k < 2; k++) {
535 if (e->GetData().m_triangles[k]->GetData().m_visible) {
536 nvisible++;
537 }
538 }
539 if (nvisible == 2) {
540 m_edgesToDelete.PushBack(e);
541 }
542 else if (nvisible == 1) {
543 e->GetData().m_newFace = MakeConeFace(e, v0);
544 m_edgesToUpdate.PushBack(e);
545 }
546 e = tmp;
547 } while (e != eHead);
548 return true;
549}
550bool ICHull::MakeCCW(CircularListElement<TMMTriangle>* f,
551 CircularListElement<TMMEdge>* e,
552 CircularListElement<TMMVertex>* v)
553{
554 // the visible face adjacent to e
555 CircularListElement<TMMTriangle>* fv;
556 if (e->GetData().m_triangles[0]->GetData().m_visible) {
557 fv = e->GetData().m_triangles[0];
558 }
559 else {
560 fv = e->GetData().m_triangles[1];
561 }
562
563 // set vertex[0] and vertex[1] to have the same orientation as the corresponding vertices of fv.
564 int32_t i; // index of e->m_vertices[0] in fv
565 CircularListElement<TMMVertex>* v0 = e->GetData().m_vertices[0];
566 CircularListElement<TMMVertex>* v1 = e->GetData().m_vertices[1];
567 for (i = 0; fv->GetData().m_vertices[i] != v0; i++)
568 ;
569
570 if (fv->GetData().m_vertices[(i + 1) % 3] != e->GetData().m_vertices[1]) {
571 f->GetData().m_vertices[0] = v1;
572 f->GetData().m_vertices[1] = v0;
573 }
574 else {
575 f->GetData().m_vertices[0] = v0;
576 f->GetData().m_vertices[1] = v1;
577 // swap edges
578 CircularListElement<TMMEdge>* tmp = f->GetData().m_edges[0];
579 f->GetData().m_edges[0] = f->GetData().m_edges[1];
580 f->GetData().m_edges[1] = tmp;
581 }
582 f->GetData().m_vertices[2] = v;
583 return true;
584}
585bool ICHull::CleanUp(uint32_t& addedPoints)
586{
587 bool r0 = CleanEdges();
588 bool r1 = CleanTriangles();
589 bool r2 = CleanVertices(addedPoints);
590 return r0 && r1 && r2;
591}
592bool ICHull::CleanEdges()
593{
594 // integrate the new faces into the data structure
595 CircularListElement<TMMEdge>* e;
596 const size_t ne_update = m_edgesToUpdate.Size();
597 for (size_t i = 0; i < ne_update; ++i) {
598 e = m_edgesToUpdate[i];
599 if (e->GetData().m_newFace) {
600 if (e->GetData().m_triangles[0]->GetData().m_visible) {
601 e->GetData().m_triangles[0] = e->GetData().m_newFace;
602 }
603 else {
604 e->GetData().m_triangles[1] = e->GetData().m_newFace;
605 }
606 e->GetData().m_newFace = 0;
607 }
608 }
609 // delete edges maked for deletion
610 CircularList<TMMEdge>& edges = m_mesh.GetEdges();
611 const size_t ne_delete = m_edgesToDelete.Size();
612 for (size_t i = 0; i < ne_delete; ++i) {
613 edges.Delete(m_edgesToDelete[i]);
614 }
615 m_edgesToDelete.Resize(0);
616 m_edgesToUpdate.Resize(0);
617 return true;
618}
619bool ICHull::CleanTriangles()
620{
621 CircularList<TMMTriangle>& triangles = m_mesh.GetTriangles();
622 const size_t nt_delete = m_trianglesToDelete.Size();
623 for (size_t i = 0; i < nt_delete; ++i) {
624 triangles.Delete(m_trianglesToDelete[i]);
625 }
626 m_trianglesToDelete.Resize(0);
627 return true;
628}
629bool ICHull::CleanVertices(uint32_t& addedPoints)
630{
631 // mark all vertices incident to some undeleted edge as on the hull
632 CircularList<TMMEdge>& edges = m_mesh.GetEdges();
633 CircularListElement<TMMEdge>* e = edges.GetHead();
634 size_t nE = edges.GetSize();
635 for (size_t i = 0; i < nE; i++) {
636 e->GetData().m_vertices[0]->GetData().m_onHull = true;
637 e->GetData().m_vertices[1]->GetData().m_onHull = true;
638 e = e->GetNext();
639 }
640 // delete all the vertices that have been processed but are not on the hull
641 CircularList<TMMVertex>& vertices = m_mesh.GetVertices();
642 CircularListElement<TMMVertex>* vHead = vertices.GetHead();
643 CircularListElement<TMMVertex>* v = vHead;
644 v = v->GetPrev();
645 do {
646 if (v->GetData().m_tag && !v->GetData().m_onHull) {
647 CircularListElement<TMMVertex>* tmp = v->GetPrev();
648 vertices.Delete(v);
649 v = tmp;
650 addedPoints--;
651 }
652 else {
653 v->GetData().m_duplicate = 0;
654 v->GetData().m_onHull = false;
655 v = v->GetPrev();
656 }
657 } while (v->GetData().m_tag && v != vHead);
658 return true;
659}
660void ICHull::Clear()
661{
662 m_mesh.Clear();
663 m_edgesToDelete.Resize(0);
664 m_edgesToUpdate.Resize(0);
665 m_trianglesToDelete.Resize(0);
666 m_isFlat = false;
667}
668const ICHull& ICHull::operator=(ICHull& rhs)
669{
670 if (&rhs != this) {
671 m_mesh.Copy(rhs.m_mesh);
672 m_edgesToDelete = rhs.m_edgesToDelete;
673 m_edgesToUpdate = rhs.m_edgesToUpdate;
674 m_trianglesToDelete = rhs.m_trianglesToDelete;
675 m_isFlat = rhs.m_isFlat;
676 }
677 return (*this);
678}
679bool ICHull::IsInside(const Vec3<double>& pt0, const double eps)
680{
681 const Vec3<double> pt(pt0.X(), pt0.Y(), pt0.Z());
682 if (m_isFlat) {
683 size_t nT = m_mesh.m_triangles.GetSize();
684 Vec3<double> ver0, ver1, ver2, a, b, c;
685 double u, v;
686 for (size_t t = 0; t < nT; t++) {
687 ver0.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.X();
688 ver0.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Y();
689 ver0.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Z();
690 ver1.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.X();
691 ver1.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Y();
692 ver1.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Z();
693 ver2.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.X();
694 ver2.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Y();
695 ver2.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Z();
696 a = ver1 - ver0;
697 b = ver2 - ver0;
698 c = pt - ver0;
699 u = c * a;
700 v = c * b;
701 if (u >= 0.0 && u <= 1.0 && v >= 0.0 && u + v <= 1.0) {
702 return true;
703 }
704 m_mesh.m_triangles.Next();
705 }
706 return false;
707 }
708 else {
709 size_t nT = m_mesh.m_triangles.GetSize();
710 Vec3<double> ver0, ver1, ver2;
711 double vol;
712 for (size_t t = 0; t < nT; t++) {
713 ver0.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.X();
714 ver0.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Y();
715 ver0.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Z();
716 ver1.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.X();
717 ver1.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Y();
718 ver1.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Z();
719 ver2.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.X();
720 ver2.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Y();
721 ver2.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Z();
722 vol = ComputeVolume4(ver0, ver1, ver2, pt);
723 if (vol < eps) {
724 return false;
725 }
726 m_mesh.m_triangles.Next();
727 }
728 return true;
729 }
730}
731}
732