1 /*************************************************************************
2 * *
3 * Vega FEM Simulation Library Version 3.1 *
4 * *
5 * "mesher" library , Copyright (C) 2016 USC *
6 * All rights reserved. *
7 * *
8 * Code authors: Danyong Zhao, Yijing Li, Jernej Barbic *
9 * http://www.jernejbarbic.com/code *
10 * *
11 * Funding: National Science Foundation *
12 * *
13 * This library is free software; you can redistribute it and/or *
14 * modify it under the terms of the BSD-style license that is *
15 * included with this library in the file LICENSE.txt *
16 * *
17 * This library is distributed in the hope that it will be useful, *
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the file *
20 * LICENSE.TXT for more details. *
21 * *
22 *************************************************************************/
23
24 #include <set>
25 #include <cfloat>
26 #include <limits.h>
27 #include <algorithm>
28 #include <fstream>
29
30 using namespace std;
31
32 #include "matrixIO.h"
33 #include "isosurfaceMesher.h"
34 #include "marchingCubes.h"
35 #include "delaunayMesher.h"
36 #include "objMeshOrientable.h"
37 #include "performanceCounter.h"
38
39 #ifndef M_PI
40 #define M_PI 3.1415926525897932384
41 #endif
42
IsosurfaceMesher(DistanceFieldBase * f)43 IsosurfaceMesher::IsosurfaceMesher(DistanceFieldBase * f)
44 : field(f), detailedSurfaceMesh(NULL), isovalue(0.0), angularBound(30.0), radialBound(1.), numInitialSamples(200), splitIsosurfaceMesh(NULL), fieldDiagonal(0), objMeshOctree(NULL), isoMesh(NULL), globalLoopIndex(0)
45 {
46 delaunay.computeVoronoiEdgeModification(true);
47 checkDelaunay = saveMarchingCubesObj = false;
48 ignoredSplitComponentVtxRatio = 0.01;
49 cosAngularBound = 0.5;
50 splitComponentIndex = oneMeshLoopIndex = 0;
51 marchingCubesMesh = NULL;
52 }
53
IsosurfaceMesher(ObjMesh * d)54 IsosurfaceMesher::IsosurfaceMesher(ObjMesh * d)
55 : field(NULL), detailedSurfaceMesh(d), isovalue(0.0), angularBound(30.0), radialBound(1.), numInitialSamples(200), splitIsosurfaceMesh(NULL), fieldDiagonal(0), objMeshOctree(NULL), isoMesh(NULL), globalLoopIndex(0)
56 {
57 delaunay.computeVoronoiEdgeModification(true);
58 checkDelaunay = saveMarchingCubesObj = false;
59 ignoredSplitComponentVtxRatio = 0.01;
60 cosAngularBound = 0.5;
61 splitComponentIndex = oneMeshLoopIndex = 0;
62 marchingCubesMesh = NULL;
63 }
64
~IsosurfaceMesher()65 IsosurfaceMesher::~IsosurfaceMesher()
66 {
67 delete isoMesh;
68 delete objMeshOctree;
69 delete marchingCubesMesh;
70 }
71
setStepping(double iso,double a,double r,size_t nis,double ep,Query::Type qt)72 void IsosurfaceMesher::setStepping(double iso, double a, double r, size_t nis, double ep, Query::Type qt)
73 {
74 queryType = qt;
75 epsilon = ep;
76 assert(epsilon >= 0 && epsilon <= 1);
77 isovalue = iso;
78 angularBound = a;
79 radialBound = r;
80 numInitialSamples = nis;
81
82 //clamp angle bound
83 cosAngularBound = cos(angularBound * (M_PI / 180.0));
84 if (cosAngularBound < sqrt(3.) / 2.)
85 cosAngularBound = sqrt(3.) / 2.;
86
87 if (detailedSurfaceMesh == NULL)
88 {
89 PerformanceCounter pc;
90 pc.StartCounter();
91 marchingCubesMesh = MarchingCubes::compute(field, isovalue, 0);
92 assert(marchingCubesMesh);
93 pc.StopCounter();
94 cout << "Finished marching cubes: " << pc.GetElapsedTime() << endl;
95 //regroup the marching cubes mesh according to connected components
96 splitIsosurfaceMesh = marchingCubesMesh->splitIntoConnectedComponents(0, 0);
97
98 if (saveMarchingCubesObj == false)
99 {
100 delete marchingCubesMesh;
101 marchingCubesMesh = NULL;
102 }
103 }
104 else
105 splitIsosurfaceMesh = detailedSurfaceMesh->splitIntoConnectedComponents(0, 0);
106 assert(splitIsosurfaceMesh);
107
108 delete isoMesh;
109 isoMesh = new ObjMesh;
110
111 if (field)
112 fieldDiagonal = field->diagonal();
113 else
114 {
115 Vec3d bmin, bmax;
116 detailedSurfaceMesh->getBoundingBox(1.0, &bmin, &bmax);
117 fieldDiagonal = len(bmax - bmin);
118 }
119
120 splitComponent.resize(splitIsosurfaceMesh->getNumGroups(), NULL);
121 for (unsigned int i = 0; i < splitIsosurfaceMesh->getNumGroups(); i++)
122 {
123 // get one connected component and extract an isosurface from it
124 ObjMesh * singleMesh = splitIsosurfaceMesh->extractGroup(i);
125 splitComponent[i] = singleMesh;
126 assert(singleMesh);
127 }
128
129 splitComponentIndex = 0;
130 }
131
endStepping() const132 bool IsosurfaceMesher::endStepping() const
133 {
134 return splitComponentIndex == splitComponent.size();
135 }
136
doOneStep()137 void IsosurfaceMesher::doOneStep()
138 {
139 if (splitComponentIndex == splitComponent.size())
140 return;
141 ObjMesh * objMesh = splitComponent[splitComponentIndex];
142 if (oneMeshLoopIndex == 0)
143 {
144 initializeOneMesh(objMesh);
145 }
146
147 if (addOneTriangle() == false) // return false if it stops
148 {
149 ObjMesh * mesh = buildCurrentSurfaceMesh();
150 if (mesh != NULL)
151 isoMesh->appendMesh(mesh);
152 splitComponentIndex++;
153 oneMeshLoopIndex = 0;
154 }
155 else
156 oneMeshLoopIndex++;
157 }
158
initializeOneMesh(ObjMesh * objMesh)159 bool IsosurfaceMesher::initializeOneMesh(ObjMesh * objMesh)
160 {
161 // create initial samples (run marching cubes, then do low-discrepancy sampling)
162 delaunay.clear();
163
164 int maxNumTrianglesInLeafNode = 5;
165 int maxTreeDepth = 10;
166
167 //octreeCounter.start();
168 delete objMeshOctree;
169 objMeshOctree = new ObjMeshOctree<TriangleBasic>(objMesh, maxNumTrianglesInLeafNode, maxTreeDepth);
170 //octreeCounter.stop();
171 vector<Vec3d> points;
172 generateInitialPointSetByDiscrepancy(objMesh, points, numInitialSamples);
173
174 bool success = delaunay.computeDelaunayTetrahedralization(points, epsilon, queryType);
175
176 if (success == false)
177 {
178 return false;
179 }
180
181 // initialize containers to store and query IsoFaces
182 radiusSet.clear();
183 angleSet.clear();
184 faceIters.clear();
185 return true;
186 }
187
compute(double o,double a,double r,int nis,double ep,Query::Type qt,int maxNumberOfIterations,double maxTimeSeconds)188 bool IsosurfaceMesher::compute(double o, double a, double r, int nis, double ep, Query::Type qt, int maxNumberOfIterations, double maxTimeSeconds)
189 {
190 cout << "Running isosurfaceMesher" << endl;
191 setStepping(o, a, r, nis, ep, qt);
192
193 // marchingCubesMesh->save("MarchingCubes.obj");
194 // delaunayUpdateCounter.reset();
195 // octreeCounter.reset();
196 // intersectionCounter.reset();
197
198 bool timeout = false;
199 if (maxTimeSeconds < 0)
200 maxTimeSeconds = DBL_MAX;
201
202 if (maxNumberOfIterations < 0)
203 maxNumberOfIterations = INT_MAX;
204 for (size_t i = 0; i < splitComponent.size(); i++)
205 {
206 // get one connected component and extract an isosurface from it
207 ObjMesh * singleMesh = splitComponent[i];
208 if (singleMesh->getNumVertices() < splitIsosurfaceMesh->getNumVertices() * ignoredSplitComponentVtxRatio)
209 continue;
210
211 ObjMesh * mesh = computeOneMesh(singleMesh, maxNumberOfIterations, maxTimeSeconds, timeout);
212
213 if (mesh)
214 {
215 isoMesh->appendMesh(mesh);
216 delete mesh;
217 }
218
219 if (timeout || maxNumberOfIterations <= 0 || maxTimeSeconds <= 0.0)
220 break;
221 }
222
223 // clean intermediate vars
224 delete splitIsosurfaceMesh;
225 splitIsosurfaceMesh = NULL;
226 for (size_t i = 0; i < splitComponent.size(); i++)
227 delete splitComponent[i];
228 splitComponent.clear();
229
230 return timeout;
231 }
232
addOneTriangle(double * targetRadius)233 bool IsosurfaceMesher::addOneTriangle(double * targetRadius)
234 {
235 // remove faces deleted in delaunay mesh from radiusSet, angleSet and faceIters
236 const VoronoiEdgeMap & deletedVEdges = delaunay.getDeletedVoronoiEdges();
237 for (VEdgeCIter it = deletedVEdges.begin(); it != deletedVEdges.end(); it++)
238 {
239 const UTriKey & tri = it->first;
240 IterPairMap::iterator it2 = faceIters.find(tri);
241 if (it2 != faceIters.end())
242 { // a deleted face is inside face
243 // cout << "erase isosurface: " << it2->first.v[0] << " " << it2->first.v[1] << " " << it2->first.v[2] << endl;
244 RadiusIter rit = it2->second.first;
245 AngleIter ait = it2->second.second;
246 radiusSet.erase(rit);
247 angleSet.erase(ait);
248 faceIters.erase(it2);
249 }
250 }
251
252 // test whether new faces in delaunay should be added to radiusSet, angleSet and faceIters
253 const VoronoiEdgeMap & addedVEdges = delaunay.getAddedVoronoiEdges();
254
255 // cout << numIter << ": added VE Edges: " << addedVEdges.size() << ": ";
256 for (VEdgeCIter it = addedVEdges.begin(); it != addedVEdges.end(); it++)
257 {
258 const UTriKey & tri = it->first;
259 const VoronoiEdge & vedge = it->second;
260 // build a Voronoi edge from the center of one ball to the center of the other ball, or to infinity
261 // test intersection of octree with this Voronoi edge. If yes, the tet face of this Voronoi edge should be one of the IsoFaces
262 Vec3d start, dir, isopoint(0.0); // isopoint is the intersection point
263 bool isRay = false;
264 if (vedge.isFinite())
265 {
266 start = vedge.start;
267 dir = vedge.end - start;
268 isRay = false;
269 }
270 else
271 {
272 isRay = true;
273 start = vedge.start;
274 dir = vedge.direction;
275 }
276 // intersectionCounter.start();
277
278 bool onIsosurface = false;
279 {
280 Vec3d segmentEnd;
281 Vec3d segmentStart = start;
282
283 double minT = DBL_MAX;
284 Vec3d direction = dir;
285 if (isRay)
286 {
287 direction *= 1.1 * fieldDiagonal * 100.0;
288 segmentEnd = segmentStart + direction; // this line will be long enough to determine intersection
289 }
290 else
291 segmentEnd = segmentStart + direction;
292
293 // Too large value at segmentStart may cause precison problems in the octree. So we swap them if needed.
294 if (len2(segmentStart) > len2(segmentEnd))
295 {
296 Vec3d tmp = segmentStart;
297 segmentStart = segmentEnd;
298 segmentEnd = tmp;
299 }
300 double length2 = len2(segmentEnd - segmentStart);
301 std::vector<int> triangleList;
302 std::vector<Vec3d> intersectionList;
303
304 std::vector<TriangleBasic*> tlist;
305 objMeshOctree->lineSegmentIntersection(tlist, segmentStart, segmentEnd, &intersectionList);
306 for (size_t i = 0; i < tlist.size(); i++)
307 triangleList.push_back(tlist[i]->index());
308
309 for (unsigned i = 0; i < triangleList.size(); i++)
310 {
311 double t2 = len2(intersectionList[i] - segmentStart) / length2;
312 if (t2 < minT)
313 {
314 minT = t2;
315 isopoint = intersectionList[i];
316 onIsosurface = true;
317 }
318 }
319 }
320
321 if (onIsosurface)
322 {
323 IsoFace isoface(tri, delaunay.getVertex(tri[0]), delaunay.getVertex(tri[1]), delaunay.getVertex(tri[2]), isopoint);
324 RadiusIter rit = (radiusSet.insert(isoface)).first;
325 AngleIter ait = (angleSet.insert(isoface)).first;
326 faceIters[tri] = IterPair(rit, ait);
327 }
328 }
329
330 // compute maximum radius among the surface faces
331 if (radiusSet.size() == 0 || angleSet.size() == 0)
332 {
333 // strange, no iso faces now
334 if (targetRadius)
335 *targetRadius = 0.;
336 return false;
337 }
338
339 //cout << radiusSet.size() << endl;
340
341 for (IsoFaceRadiusSet::iterator itr = radiusSet.begin(); itr != radiusSet.end(); itr++)
342 {
343 IsoFaceRadiusSet::iterator::value_type targetRadiusFace = *itr;
344 if (targetRadius)
345 *targetRadius = targetRadiusFace.radius;
346 if (targetRadiusFace.radius <= radialBound) // requirement fulfilled
347 {
348 //printf("Fulfilled\n");
349 break;
350 }
351
352 if (delaunay.addOnePoint(targetRadiusFace.isopoint))
353 return true;
354 //else
355 //cout << targetRadiusFace.isopoint << endl;
356 }
357
358 for (IsoFaceRadiusSet::iterator itr = angleSet.begin(); itr != angleSet.end(); itr++)
359 {
360 IsoFaceRadiusSet::iterator::value_type targetAngleFace = *itr;
361 if (targetRadius)
362 *targetRadius = targetAngleFace.radius;
363 if (targetAngleFace.maxCosAngle <= cosAngularBound) // requirement fulfilled
364 {
365 //printf("Fulfilled\n");
366 break;
367 }
368
369 if (delaunay.addOnePoint(targetAngleFace.isopoint))
370 return true;
371 //else
372 //cout << targetRadiusFace.isopoint << endl;
373 }
374
375 return false;
376 }
377
computeOneMesh(ObjMesh * objMesh,int & maxNumberOfIterations,double & maxTimeSeconds,bool & timeout)378 ObjMesh * IsosurfaceMesher::computeOneMesh(ObjMesh * objMesh, int & maxNumberOfIterations, double & maxTimeSeconds, bool & timeout)
379 {
380 // repeat:
381 // build Delaunay mesh
382 // determine faces intersected by the Voronoi diagram edges
383 // check all the faces for the splitting condition (angularBound, radialBound)
384 // split the worst face
385 timeout = false;
386 if (initializeOneMesh(objMesh) == false)
387 return NULL;
388
389 oneMeshLoopIndex = 0;
390 bool done = false;
391 PerformanceCounter loopTimeCounter; // measure elapsed time in the following loop
392 loopTimeCounter.StartCounter();
393 while (!done)
394 {
395 // return true if it won't stop
396 double targetRadius = 0;
397 if (addOneTriangle(&targetRadius) == false)
398 break;
399
400 oneMeshLoopIndex++;
401 if (oneMeshLoopIndex >= maxNumberOfIterations)
402 break;
403 loopTimeCounter.StopCounter();
404 double timeElapsed = loopTimeCounter.GetElapsedTime();
405 if (timeElapsed >= maxTimeSeconds)
406 {
407 timeout = true;
408 break;
409 }
410
411 if (oneMeshLoopIndex % 1000 == 0)
412 printf("Itr %d: radius=%.5f\n", oneMeshLoopIndex, targetRadius);
413 }
414
415 maxNumberOfIterations -= oneMeshLoopIndex;
416 loopTimeCounter.StopCounter();
417 maxTimeSeconds -= loopTimeCounter.GetElapsedTime();
418 cout << "Time spent in refining isosurface: " << loopTimeCounter.GetElapsedTime() << endl;
419
420 if (checkDelaunay)
421 {
422 cout << "Begin Delaunay check in IsosurfaceMesher: " << endl;
423 PerformanceCounter pc;
424 bool ret = delaunay.checkDelaunay();
425 pc.StopCounter();
426 cout << "Time in check: " << pc.GetElapsedTime() << ", ";
427 cout << "Check: " << ret << endl;
428 }
429 return buildCurrentSurfaceMesh();
430 }
431
buildCurrentSurfaceMesh(bool keepAllDelaunayVtx) const432 ObjMesh * IsosurfaceMesher::buildCurrentSurfaceMesh(bool keepAllDelaunayVtx) const
433 {
434 vector<double> vertices;
435 vector<int> faces;
436 if (keepAllDelaunayVtx)
437 {
438 vertices.resize(3 * delaunay.getNumVertices());
439 for (size_t i = 0; i < delaunay.getNumVertices(); i++)
440 {
441 Vec3d pos = delaunay.getVertex(i);
442 pos.convertToArray(&vertices[3 * i]);
443 }
444 faces.resize(faceIters.size() * 3);
445 int k = 0;
446 for (IterPairMap::const_iterator it = faceIters.begin(); it != faceIters.end(); it++, k++)
447 memcpy(&faces[3 * k], it->first.indices(), sizeof(int) * 3);
448 }
449 else
450 {
451 map<int, int> vertexMap; // vertex delaunay ID -> new surface mesh vtx ID
452 // go through all the triangles
453 for (IterPairMap::const_iterator it = faceIters.begin(); it != faceIters.end(); it++)
454 for (int j = 0; j < 3; j++)
455 {
456 int vtx = it->first[j]; // one vtx on one triangle
457 map<int, int>::iterator itr = vertexMap.find(vtx);
458 if (itr == vertexMap.end())
459 { // this is a new vtx in vertexMap
460 vertexMap[vtx] = vertices.size() / 3;
461 Vec3d pos = delaunay.getVertex(vtx);
462 for (int i = 0; i < 3; i++)
463 vertices.push_back(pos[i]);
464 }
465 faces.push_back(vertexMap[vtx]);
466 }
467 }
468 return new ObjMesh(vertices.size() / 3, &vertices[0], faces.size() / 3, &faces[0]);
469 }
470
471 // choose some points from objMesh into points. These points should be as far away from each other as possible.
472 // The number of the points required is numInitialSamples
generateInitialPointSetByDiscrepancy(const ObjMesh * objMesh,std::vector<Vec3d> & points,unsigned int numInitialSamples)473 void IsosurfaceMesher::generateInitialPointSetByDiscrepancy(const ObjMesh * objMesh, std::vector<Vec3d> & points, unsigned int numInitialSamples)
474 {
475 srand(6);
476 int numVertices = 0;
477 double * vertices = NULL;
478 objMesh->exportGeometry(&numVertices, &vertices);
479 vector<bool> selected(numVertices, false);
480
481 // select the first vertex
482 int p = rand() % numVertices;
483
484 points.clear();
485 // put first vertex into points
486 points.push_back(vertices + 3 * p);
487 selected[p] = true;
488 vector<int> selectedInd;
489 selectedInd.push_back(p);
490
491 while (points.size() < numInitialSamples)
492 {
493 double max_distance = 0;
494 p = -1;
495 for (int i = 0; i < numVertices; i++)
496 {
497 if (selected[i] == false)
498 {
499 double cur_distance = DBL_MAX;
500 for (size_t j = 0; j < selectedInd.size(); j++)
501 {
502 Vec3d d;
503 for (int k = 0; k < 3; k++)
504 d[k] = vertices[3 * selectedInd[j] + k] - vertices[3 * i + k];
505 double distance = len2(d);
506 if (distance < cur_distance)
507 cur_distance = distance;
508 }
509 if (cur_distance > max_distance)
510 {
511 max_distance = cur_distance;
512 p = i;
513 }
514 }
515 }
516
517 if (p >= 0)
518 {
519 //new point is found
520 points.push_back(vertices + 3 * p);
521 selected[p] = true;
522 selectedInd.push_back(p);
523 }
524 else
525 {
526 numInitialSamples = points.size();
527 break;
528 }
529 }
530
531 // points.push_back(Vec3d(0.0));
532 free(vertices);
533 }
534
535 // check intersection of a ray and an octree
536 // v = segmentStart
537 // if isRay == true, send a ray starting at v and points toward direction
538 // else, send a line segment from v to (v + direction)
539 // store the intersection distance (from v along the line to the intersection point) in intersectionDistance
intersectionUsingOctree(const Vec3d & segmentStart,const Vec3d & dir,bool isRay,Vec3d & intersectionPoint)540 bool IsosurfaceMesher::intersectionUsingOctree(const Vec3d & segmentStart, const Vec3d & dir, bool isRay, Vec3d & intersectionPoint)
541 {
542 Vec3d segmentEnd;
543
544 double minT = DBL_MAX;
545
546 Vec3d direction = dir;
547 if (isRay)
548 {
549 direction *= 1.1 * fieldDiagonal * 100.0;
550 segmentEnd = segmentStart + direction; // this line will be long enough to determine intersection
551 }
552 else
553 segmentEnd = segmentStart + direction;
554
555 double length = len(segmentEnd - segmentStart);
556
557 std::vector<int> triangleList;
558 std::vector<Vec3d> intersectionList;
559 vector<TriangleBasic*> tlist;
560 objMeshOctree->lineSegmentIntersection(tlist, segmentStart, segmentEnd, &intersectionList);
561 for (size_t i = 0; i < tlist.size(); i++)
562 triangleList.push_back(tlist[i]->index());
563
564 for (unsigned i = 0; i < triangleList.size(); i++)
565 {
566 double t = len(intersectionList[i] - segmentStart) / length;
567 if (t < minT)
568 {
569 minT = t;
570 intersectionPoint = intersectionList[i];
571 }
572 }
573
574 return minT != DBL_MAX;
575 }
576
enforceManifoldnessAndOrientNormals(ObjMesh * & objMesh)577 bool IsosurfaceMesher::enforceManifoldnessAndOrientNormals(ObjMesh* &objMesh)
578 {
579 bool flag1 = false, flag2 = false;
580 //remove non-manifold faces and edges
581 do
582 {
583 flag1 = false;
584 int removeHangingFaceNum = 0;
585 do
586 {
587 removeHangingFaceNum = objMesh->removeHangingFaces();
588 if (removeHangingFaceNum > 0)
589 flag1 = true;
590 }
591 while (removeHangingFaceNum > 0);
592
593 int removeEdge = 0;
594 flag2 = false;
595 do
596 {
597 removeEdge = objMesh->removeNonManifoldEdges();
598 if (removeEdge > 0)
599 flag2 = true;
600 }
601 while (removeEdge > 0);
602 }
603 while (flag1 || flag2);
604
605 // process each connected component
606 ObjMesh * connectedComponentMesh = objMesh->splitIntoConnectedComponents();
607 delete objMesh;
608 objMesh = new ObjMesh;
609
610 // for each connected component
611 for (unsigned int i = 0; i < connectedComponentMesh->getNumGroups(); i++)
612 {
613 ObjMesh * component = connectedComponentMesh->extractGroup(i);
614
615 // convert to oriented mesh
616 ObjMeshOrientable objMeshOrientable(component, 1, NULL, 0);
617 ObjMesh * singleMesh = objMeshOrientable.GenerateOrientedMesh();
618 delete component;
619
620 vector<int> group;
621 for (unsigned int i = 0; i < singleMesh->getNumGroups(); i++)
622 group.push_back(i);
623 singleMesh->mergeGroups(group);
624
625 // flip all face normals if orientation is incorrect
626 if (singleMesh->computeVolume() < 0)
627 for (unsigned int i = 0; i < singleMesh->getNumGroups(); i++)
628 {
629 ObjMesh::Group * group = (ObjMesh::Group*)singleMesh->getGroupHandle(i);
630 for (unsigned int j = 0; j < group->getNumFaces(); j++)
631 {
632 ObjMesh::Face * face = (ObjMesh::Face*)group->getFaceHandle(j);
633 face->reverseVertices();
634 }
635 }
636
637 objMesh->appendMesh(singleMesh);
638 delete singleMesh;
639 }
640
641 delete connectedComponentMesh;
642 return true;
643 }
644
IsoFace(const UTriKey & t,const Vec3d & v0,const Vec3d & v1,const Vec3d & v2,const Vec3d & c)645 IsosurfaceMesher::IsoFace::IsoFace(const UTriKey & t, const Vec3d & v0, const Vec3d & v1, const Vec3d & v2, const Vec3d & c) : tri(t)
646 {
647 isopoint = c;
648 radius = (len(v0 - c) + len(v1 - c) + len(v2 - c)) / 3;
649 Vec3d base1 = v1 - v0, base2 = v2 - v0;
650 Vec3d base3 = base2 - base1; // = v2 - v1
651 base1.normalize();
652 base2.normalize();
653 base3.normalize();
654 double cosAngle[3] = { dot(base1, base2), dot(base1, base3), dot(base2, base3) };
655 maxCosAngle = max(cosAngle[0], max(cosAngle[1], cosAngle[2]));
656 }
657
operator ()(const IsoFace & a,const IsoFace & b) const658 bool IsosurfaceMesher::IsoFaceRadiusCompare::operator() (const IsoFace & a, const IsoFace & b) const
659 {
660 if (a.radius > b.radius)
661 return true;
662 if (a.radius < b.radius)
663 return false;
664 return a.tri < b.tri;
665 }
666
operator ()(const IsoFace & a,const IsoFace & b) const667 bool IsosurfaceMesher::IsoFaceAngleCompare::operator() (const IsoFace & a, const IsoFace & b) const
668 {
669 if (a.maxCosAngle > b.maxCosAngle)
670 return true;
671 if (a.maxCosAngle < b.maxCosAngle)
672 return false;
673 return a.tri < b.tri;
674 }
675
getMesh(bool enforceManifold)676 ObjMesh * IsosurfaceMesher::getMesh(bool enforceManifold)
677 {
678 // cout << "delaunay update time: " << delaunayUpdateCounter.getElapsedTime() << endl;
679 // cout << "octree time: " << octreeCounter.getElapsedTime() << endl;
680 // cout << "intersection time: " << intersectionCounter.getElapsedTime() << endl;
681
682 if (enforceManifold)
683 {
684 // outputMesh->save("non-manifold.obj");
685 enforceManifoldnessAndOrientNormals(isoMesh);
686 }
687 else if (field)
688 {
689 // orient faces consistently with the distance field gradient
690 // iterate over all groups and faces
691 for (unsigned int groupIndex = 0; groupIndex < isoMesh->getNumGroups(); groupIndex++)
692 {
693 const ObjMesh::Group * groupHandle = isoMesh->getGroupHandle(groupIndex);
694 for (unsigned int faceIndex = 0; faceIndex < groupHandle->getNumFaces(); faceIndex++)
695 {
696 ObjMesh::Face * faceHandle = (ObjMesh::Face*) groupHandle->getFaceHandle(faceIndex);
697 if (faceHandle->getNumVertices() < 3)
698 cout << "Warning: encountered a face (group=" << groupIndex << ",face=" << faceIndex << ") with fewer than 3 vertices." << endl;
699
700 // compute face centroid
701 Vec3d centroid(0, 0, 0);
702 for (unsigned int vertexIndex = 0; vertexIndex < faceHandle->getNumVertices(); vertexIndex++)
703 {
704 const ObjMesh::Vertex * vertexHandle = faceHandle->getVertexHandle(vertexIndex);
705 Vec3d vertexPos = isoMesh->getPosition(vertexHandle->getPositionIndex());
706 centroid += vertexPos;
707
708 }
709 if (faceHandle->getNumVertices() > 0)
710 centroid *= 1.0 / faceHandle->getNumVertices();
711
712 // reverse face if normal is in the opposite direction to the distance field gradient
713 Vec3d normal = isoMesh->computeFaceNormal(*faceHandle);
714 Vec3d gradient = field->gradient(centroid);
715 if (dot(normal, gradient) < 0)
716 faceHandle->reverseVertices();
717 }
718 }
719 }
720 return new ObjMesh(*isoMesh);
721 }
722
723