OpenVDB  1.1.0
MeshToVolume.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 
31 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/Types.h>
35 #include <openvdb/math/FiniteDifference.h>
36 #include <openvdb/math/Operators.h>
37 #include <openvdb/math/Proximity.h>
38 #include <openvdb/tools/LevelSetUtil.h>
39 #include <openvdb/tools/Morphology.h>
40 #include <openvdb/util/NullInterrupter.h>
41 #include <openvdb/util/Util.h>
42 
43 #include <tbb/blocked_range.h>
44 #include <tbb/parallel_for.h>
45 #include <tbb/parallel_reduce.h>
46 
47 #include <list>
48 #include <deque>
49 #include <limits>
50 
51 namespace openvdb {
53 namespace OPENVDB_VERSION_NAME {
54 namespace tools {
55 
57 enum { GENERATE_PRIM_INDEX_GRID = 0x1 };
58 
59 
60 // MeshToVolume
61 template<typename DistGridT, typename InterruptT = util::NullInterrupter>
63 {
64 public:
67  typedef typename DistGridT::TreeType DistTreeT;
68  typedef typename DistTreeT::ValueType DistValueT;
69  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
71  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
74 
75  MeshToVolume(openvdb::math::Transform::Ptr&, int conversionFlags = 0,
76  InterruptT *interrupter = NULL, int signSweeps = 1);
77 
89  void convertToLevelSet(
90  const std::vector<Vec3s>& pointList,
91  const std::vector<Vec4I>& polygonList,
94 
103  void convertToUnsignedDistanceField(const std::vector<Vec3s>& pointList,
104  const std::vector<Vec4I>& polygonList, DistValueT exBandWidth);
105 
106  void clear();
107 
109  typename DistGridT::Ptr distGridPtr() const { return mDistGrid; }
110 
113  typename IndexGridT::Ptr indexGridPtr() const { return mIndexGrid; }
114 
115 private:
116  // disallow copy by assignment
117  void operator=(const MeshToVolume<DistGridT, InterruptT>&) {}
118 
119  void doConvert(const std::vector<Vec3s>&, const std::vector<Vec4I>&,
120  DistValueT exBandWidth, DistValueT inBandWidth, bool unsignedDistField = false);
121 
122  openvdb::math::Transform::Ptr mTransform;
123  int mConversionFlags, mSignSweeps;
124 
125  typename DistGridT::Ptr mDistGrid;
126  typename IndexGridT::Ptr mIndexGrid;
127  typename StencilGridT::Ptr mIntersectingVoxelsGrid;
128 
129  InterruptT *mInterrupter;
130 };
131 
132 
135 
136 
137 // Internal utility objects and implementation details
138 
139 namespace internal {
140 
141 
142 template<typename ValueType>
143 struct Tolerance
144 {
145  static ValueType epsilon() { return ValueType(1e-7); }
146  static ValueType minNarrowBandWidth() { return ValueType(1.0 + 1e-6); }
147 };
148 
149 
150 template<typename DistTreeT, typename IndexTreeT>
151 inline void
152 combine(DistTreeT& lhsDist, IndexTreeT& lhsIndex, DistTreeT& rhsDist, IndexTreeT& rhsIndex)
153 {
154  typedef typename DistTreeT::ValueType DistValueT;
155  typename tree::ValueAccessor<DistTreeT> lhsDistAccessor(lhsDist);
156  typename tree::ValueAccessor<IndexTreeT> lhsIndexAccessor(lhsIndex);
157  typename tree::ValueAccessor<IndexTreeT> rhsIndexAccessor(rhsIndex);
158  typename DistTreeT::LeafCIter iter = rhsDist.cbeginLeaf();
159 
160  DistValueT rhsValue;
161  Coord ijk;
162 
163  for ( ; iter; ++iter) {
164  typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
165 
166  for ( ; it; ++it) {
167 
168  ijk = it.getCoord();
169  rhsValue = it.getValue();
170  DistValueT& lhsValue = const_cast<DistValueT&>(lhsDistAccessor.getValue(ijk));
171 
172  if (-rhsValue < std::abs(lhsValue)) {
173  lhsValue = rhsValue;
174  lhsIndexAccessor.setValue(ijk, rhsIndexAccessor.getValue(ijk));
175  }
176  }
177  }
178 }
179 
180 
182 
183 
191 template<typename DistTreeT, typename InterruptT = util::NullInterrupter>
193 {
194 public:
197  typedef typename DistTreeT::ValueType DistValueT;
199  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
201  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
204 
205  MeshVoxelizer(const std::vector<Vec3s>& pointList,
206  const std::vector<Vec4I>& polygonList, InterruptT *interrupter = NULL);
207 
209 
210  void runParallel();
211  void runSerial();
212 
214  void operator() (const tbb::blocked_range<size_t> &range);
215  void join(MeshVoxelizer<DistTreeT, InterruptT>& rhs);
216 
217  DistTreeT& sqrDistTree() { return mSqrDistTree; }
218  IndexTreeT& primIndexTree() { return mPrimIndexTree; }
219  StencilTreeT& intersectionTree() { return mIntersectionTree; }
220 
221 private:
222  // disallow copy by assignment
223  void operator=(const MeshVoxelizer<DistTreeT, InterruptT>&) {}
224 
225  inline bool shortEdge(const Vec3d&, const Vec3d&, const Vec3d&) const;
226  bool evalVoxel(const Coord& ijk, const Int32 polyIdx);
227 
228  std::vector<Vec3s> const * const mPointList;
229  std::vector<Vec4I> const * const mPolygonList;
230 
231  DistTreeT mSqrDistTree;
232  DistAccessorT mSqrDistAccessor;
233 
234  IndexTreeT mPrimIndexTree;
235  IndexAccessorT mPrimIndexAccessor;
236 
237  StencilTreeT mIntersectionTree;
238  StencilAccessorT mIntersectionAccessor;
239 
240  // Used internally for acceleration
241  IndexTreeT mLastPrimTree;
242  IndexAccessorT mLastPrimAccessor;
243 
244  InterruptT *mInterrupter;
245 };
246 
247 
248 template<typename DistTreeT, typename InterruptT>
249 void
251 {
252  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList->size()), *this);
253 }
254 
255 template<typename DistTreeT, typename InterruptT>
256 void
258 {
259  (*this)(tbb::blocked_range<size_t>(0, mPolygonList->size()));
260 }
261 
262 template<typename DistTreeT, typename InterruptT>
264  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
265  InterruptT *interrupter)
266  : mPointList(&pointList)
267  , mPolygonList(&polygonList)
268  , mSqrDistTree(std::numeric_limits<DistValueT>::max())
269  , mSqrDistAccessor(mSqrDistTree)
270  , mPrimIndexTree(Int32(util::INVALID_IDX))
271  , mPrimIndexAccessor(mPrimIndexTree)
272  , mIntersectionTree(false)
273  , mIntersectionAccessor(mIntersectionTree)
274  , mLastPrimTree(Int32(util::INVALID_IDX))
275  , mLastPrimAccessor(mLastPrimTree)
276  , mInterrupter(interrupter)
277 {
278 }
279 
280 template<typename DistTreeT, typename InterruptT>
282  MeshVoxelizer<DistTreeT, InterruptT>& rhs, tbb::split)
283  : mPointList(rhs.mPointList)
284  , mPolygonList(rhs.mPolygonList)
285  , mSqrDistTree(std::numeric_limits<DistValueT>::max())
286  , mSqrDistAccessor(mSqrDistTree)
287  , mPrimIndexTree(Int32(util::INVALID_IDX))
288  , mPrimIndexAccessor(mPrimIndexTree)
289  , mIntersectionTree(false)
290  , mIntersectionAccessor(mIntersectionTree)
291  , mLastPrimTree(Int32(util::INVALID_IDX))
292  , mLastPrimAccessor(mLastPrimTree)
293  , mInterrupter(rhs.mInterrupter)
294 {
295 }
296 
297 template<typename DistTreeT, typename InterruptT>
298 inline bool
300  const Vec3d& v0, const Vec3d& v1, const Vec3d& v2) const
301 {
302  double edge_max = std::abs(v1[0] - v0[0]);
303  edge_max = std::max(edge_max, std::abs(v1[1] - v0[1]));
304  edge_max = std::max(edge_max, std::abs(v1[2] - v0[2]));
305  edge_max = std::max(edge_max, std::abs(v0[0] - v2[0]));
306  edge_max = std::max(edge_max, std::abs(v0[1] - v2[1]));
307  edge_max = std::max(edge_max, std::abs(v0[2] - v2[2]));
308  return edge_max < 200.0;
309 }
310 
311 template<typename DistTreeT, typename InterruptT>
312 void
313 MeshVoxelizer<DistTreeT, InterruptT>::operator()(const tbb::blocked_range<size_t> &range)
314 {
315  std::deque<Coord> coordList;
316  StencilTreeT auxTree(false);
317  StencilAccessorT auxAcc(auxTree);
318  Coord ijk, n_ijk;
319 
320  for (size_t n = range.begin(); n < range.end(); ++n) {
321 
322  if (mInterrupter && mInterrupter->wasInterrupted()) {
323  tbb::task::self().cancel_group_execution();
324  break;
325  }
326 
327  const Int32 primIdx = n;
328  const Vec4I verts = (*mPolygonList)[n];
329 
330  Vec3d p0((*mPointList)[verts[0]]);
331  Vec3d p1((*mPointList)[verts[1]]);
332  Vec3d p2((*mPointList)[verts[2]]);
333 
334  if (shortEdge(p0, p1, p2)) {
335  coordList.clear();
336 
337 
338  ijk = util::nearestCoord(p0);
339  evalVoxel(ijk, primIdx);
340  coordList.push_back(ijk);
341 
342 
343  ijk = util::nearestCoord(p1);
344  evalVoxel(ijk, primIdx);
345  coordList.push_back(ijk);
346 
347 
348  ijk = util::nearestCoord(p2);
349  evalVoxel(ijk, primIdx);
350  coordList.push_back(ijk);
351 
352  if (util::INVALID_IDX != verts[3]) {
353  Vec3d p3((*mPointList)[verts[3]]);
354  ijk = util::nearestCoord(p3);
355  evalVoxel(ijk, primIdx);
356  coordList.push_back(ijk);
357  }
358 
359  while (!coordList.empty()) {
360  if (mInterrupter && mInterrupter->wasInterrupted()) {
361  break;
362  }
363 
364  ijk = coordList.back();
365  coordList.pop_back();
366 
367  mIntersectionAccessor.setActiveState(ijk, true);
368 
369  for (Int32 i = 0; i < 26; ++i) {
370  n_ijk = ijk + util::COORD_OFFSETS[i];
371 
372  if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
373  mLastPrimAccessor.setValue(n_ijk, n);
374  if(evalVoxel(n_ijk, n)) coordList.push_back(n_ijk);
375  }
376  }
377  }
378 
379  } else {
380 
381  ijk = util::nearestCoord(p0);
382  evalVoxel(ijk, primIdx);
383 
384  mLastPrimAccessor.setValue(ijk, primIdx);
385  auxAcc.setActiveState(ijk, true);
386 
387  while (!auxTree.empty()) {
388 
389  if (mInterrupter && mInterrupter->wasInterrupted()) {
390  break;
391  }
392 
393  typename StencilTreeT::LeafIter leafIter = auxTree.beginLeaf();
394  for (; leafIter; leafIter.next()) {
395 
396  if (mInterrupter && mInterrupter->wasInterrupted()) {
397  break;
398  }
399 
400  typename StencilTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
401  for (; iter; iter.next()) {
402  ijk = iter.getCoord();
403  iter.setValueOff();
404 
405  mIntersectionAccessor.setActiveState(ijk, true);
406 
407  for (Int32 i = 0; i < 26; ++i) {
408  n_ijk = ijk + util::COORD_OFFSETS[i];
409 
410  if (primIdx != mLastPrimAccessor.getValue(n_ijk)) {
411  mLastPrimAccessor.setValue(n_ijk, n);
412  if (evalVoxel(n_ijk, n)) auxAcc.setActiveState(n_ijk, true);
413  }
414  }
415  }
416  }
417 
418  auxTree.pruneInactive();
419  }
420  }
421  }
422 }
423 
424 template<typename DistTreeT, typename InterruptT>
425 bool
427 {
428  Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
429  const Vec4I& verts = (*mPolygonList)[polyIdx];
430  const std::vector<Vec3s>& points = *mPointList;
431 
432  // Evaluate first triangle
433  const Vec3d a(points[verts[0]]);
434  const Vec3d b(points[verts[1]]);
435  const Vec3d c(points[verts[2]]);
436 
437  Vec3d uvw;
438  double dist = (voxelCenter -
439  closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
440 
441  // Split-up quad into a second triangle and calac distance.
442  if (util::INVALID_IDX != verts[3]) {
443  const Vec3d d(points[verts[3]]);
444 
445  double secondDist = (voxelCenter -
446  closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw)).lengthSqr();
447 
448  if (secondDist < dist) dist = secondDist;
449  }
450 
451  const DistValueT tmp(dist);
452  if (tmp < std::abs(mSqrDistAccessor.getValue(ijk))) {
453  mSqrDistAccessor.setValue(ijk, -tmp);
454  mPrimIndexAccessor.setValue(ijk, polyIdx);
455  }
456 
457  return (dist < 0.86602540378443861);
458 }
459 
460 template<typename DistTreeT, typename InterruptT>
461 void
463 {
464  typename DistTreeT::LeafCIter iter = rhs.mSqrDistTree.cbeginLeaf();
465  DistValueT rhsDist;
466  Coord ijk;
467 
468  for ( ; iter; ++iter) {
469  typename DistTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
470 
471  for ( ; it; ++it) {
472 
473  ijk = it.getCoord();
474  rhsDist = it.getValue();
475  DistValueT lhsDist = mSqrDistAccessor.getValue(ijk);
476 
477  if (-rhsDist < std::abs(lhsDist)) {
478  mSqrDistAccessor.setValue(ijk, rhsDist);
479  mPrimIndexAccessor.setValue(ijk, rhs.mPrimIndexAccessor.getValue(ijk));
480  }
481  }
482  }
483 
484  mIntersectionTree.merge(rhs.mIntersectionTree);
485 }
486 
487 
489 
490 
491 // ContourTracer
494 template<typename DistTreeT, typename InterruptT = util::NullInterrupter>
496 {
497 public:
500  typedef typename DistTreeT::ValueType DistValueT;
502  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
505 
506  ContourTracer(DistTreeT&, const StencilTreeT&, InterruptT *interrupter = NULL);
508 
509  void runParallel();
510  void runSerial();
511 
513  void operator()(const tbb::blocked_range<int> &range) const;
514 
515 private:
516  void operator=(const ContourTracer<DistTreeT, InterruptT>&) {}
517 
518  int sparseScan(int slice) const;
519 
520  DistTreeT& mDistTree;
521  DistAccessorT mDistAccessor;
522 
523  const StencilTreeT& mIntersectionTree;
524  StencilAccessorT mIntersectionAccessor;
525 
526  CoordBBox mBBox;
527 
529  std::vector<Index> mStepSize;
530 
531  InterruptT *mInterrupter;
532 };
533 
534 template<typename DistTreeT, typename InterruptT>
535 void
537 {
538  tbb::parallel_for(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1), *this);
539 }
540 
541 template<typename DistTreeT, typename InterruptT>
542 void
544 {
545  (*this)(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1));
546 }
547 
548 template<typename DistTreeT, typename InterruptT>
550  DistTreeT& distTree, const StencilTreeT& intersectionTree, InterruptT *interrupter)
551  : mDistTree(distTree)
552  , mDistAccessor(mDistTree)
553  , mIntersectionTree(intersectionTree)
554  , mIntersectionAccessor(mIntersectionTree)
555  , mBBox(CoordBBox())
556  , mStepSize(0)
557  , mInterrupter(interrupter)
558 {
559  // Build the step size table for different tree value depths.
560  std::vector<Index> dims;
561  mDistTree.getNodeLog2Dims(dims);
562 
563  mStepSize.resize(dims.size()+1, 1);
564  Index exponent = 0;
565  for (int idx = static_cast<int>(dims.size()) - 1; idx > -1; --idx) {
566  exponent += dims[idx];
567  mStepSize[idx] = 1 << exponent;
568  }
569 
570  mDistTree.evalLeafBoundingBox(mBBox);
571 
572  // Make sure that mBBox coincides with the min and max corners of the internal nodes.
573  const int tileDim = mStepSize[0];
574 
575  for (size_t i = 0; i < 3; ++i) {
576 
577  int n;
578  double diff = std::abs(double(mBBox.min()[i])) / double(tileDim);
579 
580  if (mBBox.min()[i] <= tileDim) {
581  n = int(std::ceil(diff));
582  mBBox.min()[i] = - n * tileDim;
583  } else {
584  n = int(std::floor(diff));
585  mBBox.min()[i] = n * tileDim;
586  }
587 
588  n = int(std::ceil(std::abs(double(mBBox.max()[i] - mBBox.min()[i])) / double(tileDim)));
589  mBBox.max()[i] = mBBox.min()[i] + n * tileDim;
590  }
591 }
592 
593 template<typename DistTreeT, typename InterruptT>
596  : mDistTree(rhs.mDistTree)
597  , mDistAccessor(mDistTree)
598  , mIntersectionTree(rhs.mIntersectionTree)
599  , mIntersectionAccessor(mIntersectionTree)
600  , mBBox(rhs.mBBox)
601  , mStepSize(rhs.mStepSize)
602  , mInterrupter(rhs.mInterrupter)
603 {
604 }
605 
606 template<typename DistTreeT, typename InterruptT>
607 void
608 ContourTracer<DistTreeT, InterruptT>::operator()(const tbb::blocked_range<int> &range) const
609 {
610  // Slice up the volume and trace contours.
611  int iStep = 1;
612  for (int n = range.begin(); n < range.end(); n += iStep) {
613 
614  if (mInterrupter && mInterrupter->wasInterrupted()) {
615  tbb::task::self().cancel_group_execution();
616  break;
617  }
618 
619  iStep = sparseScan(n);
620  }
621 }
622 
623 template<typename DistTreeT, typename InterruptT>
624 int
626 {
627  bool lastVoxelWasOut = true;
628  int last_k;
629 
630  Coord ijk(slice, mBBox.min()[1], mBBox.min()[2]);
631  Coord step(mStepSize[mDistAccessor.getValueDepth(ijk) + 1]);
632  Coord n_ijk;
633 
634  for (ijk[1] = mBBox.min()[1]; ijk[1] <= mBBox.max()[1]; ijk[1] += step[1]) { // j
635 
636  if (mInterrupter && mInterrupter->wasInterrupted()) {
637  break;
638  }
639 
640  step[1] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
641  step[0] = std::min(step[0], step[1]);
642 
643  for (ijk[2] = mBBox.min()[2]; ijk[2] <= mBBox.max()[2]; ijk[2] += step[2]) { // k
644 
645  step[2] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
646  step[1] = std::min(step[1], step[2]);
647  step[0] = std::min(step[0], step[2]);
648 
649  // If the current voxel is set?
650  if (mDistAccessor.isValueOn(ijk)) {
651 
652  // Is this a boundary voxel?
653  if (mIntersectionAccessor.isValueOn(ijk)) {
654 
655  lastVoxelWasOut = false;
656  last_k = ijk[2];
657 
658  } else if (lastVoxelWasOut) {
659 
660  DistValueT& val = const_cast<DistValueT&>(mDistAccessor.getValue(ijk));
661  val = -val; // flip sign
662 
663  } else {
664 
665  DistValueT val;
666  for (Int32 n = 3; n < 6; n += 2) {
667  n_ijk = ijk + util::COORD_OFFSETS[n];
668 
669  if (mDistAccessor.probeValue(n_ijk, val) && val > 0) {
670  lastVoxelWasOut = true;
671  break;
672  }
673  }
674 
675  if (lastVoxelWasOut) {
676 
677  DistValueT& v = const_cast<DistValueT&>(mDistAccessor.getValue(ijk));
678  v = -v; // flip sign
679 
680  const int tmp_k = ijk[2];
681 
682  // backtrace
683  for (--ijk[2]; ijk[2] >= last_k; --ijk[2]) {
684  if (mIntersectionAccessor.isValueOn(ijk)) break;
685  DistValueT& vb = const_cast<DistValueT&>(mDistAccessor.getValue(ijk));
686  if(vb < DistValueT(0.0)) vb = -vb; // flip sign
687  }
688 
689  last_k = tmp_k;
690  ijk[2] = tmp_k;
691 
692  } else {
693  last_k = std::min(ijk[2], last_k);
694  }
695 
696  }
697 
698  } // end isValueOn check
699  } // end k
700  } // end j
701  return step[0];
702 }
703 
704 
706 
707 
708 // IntersectingVoxelSign
712 template<typename DistTreeT>
714 {
715 public:
718  typedef typename DistTreeT::ValueType DistValueT;
720  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
722  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
726 
728  const std::vector<Vec3s>& pointList,
729  const std::vector<Vec4I>& polygonList,
730  DistTreeT& distTree,
731  IndexTreeT& indexTree,
732  StencilTreeT& intersectionTree,
733  StencilArrayT& leafs);
734 
736 
737  void runParallel();
738  void runSerial();
739 
741  void operator()(const tbb::blocked_range<size_t>&) const;
742 
743 private:
744  void operator=(const IntersectingVoxelSign<DistTreeT>&) {}
745 
746  void evalVoxel(const Coord& ijk) const;
747  Vec3d getClosestPointDir(const Coord& ijk) const;
748 
749  std::vector<Vec3s> const * const mPointList;
750  std::vector<Vec4I> const * const mPolygonList;
751 
752  DistTreeT& mDistTree;
753  DistAccessorT mDistAccessor;
754 
755  IndexTreeT& mIndexTree;
756  IndexAccessorT mIndexAccessor;
757 
758  StencilTreeT& mIntersectionTree;
759  StencilAccessorT mIntersectionAccessor;
760  StencilArrayT& mLeafs;
761 };
762 
763 template<typename DistTreeT>
764 void
766 {
767  tbb::parallel_for(mLeafs.getRange(), *this);
768 }
769 
770 template<typename DistTreeT>
771 void
773 {
774  (*this)(mLeafs.getRange());
775 }
776 
777 template<typename DistTreeT>
779  const std::vector<Vec3s>& pointList,
780  const std::vector<Vec4I>& polygonList,
781  DistTreeT& distTree,
782  IndexTreeT& indexTree,
783  StencilTreeT& intersectionTree,
784  StencilArrayT& leafs)
785  : mPointList(&pointList)
786  , mPolygonList(&polygonList)
787  , mDistTree(distTree)
788  , mDistAccessor(mDistTree)
789  , mIndexTree(indexTree)
790  , mIndexAccessor(mIndexTree)
791  , mIntersectionTree(intersectionTree)
792  , mIntersectionAccessor(mIntersectionTree)
793  , mLeafs(leafs)
794 {
795 }
796 
797 template<typename DistTreeT>
800  : mPointList(rhs.mPointList)
801  , mPolygonList(rhs.mPolygonList)
802  , mDistTree(rhs.mDistTree)
803  , mDistAccessor(mDistTree)
804  , mIndexTree(rhs.mIndexTree)
805  , mIndexAccessor(mIndexTree)
806  , mIntersectionTree(rhs.mIntersectionTree)
807  , mIntersectionAccessor(mIntersectionTree)
808  , mLeafs(rhs.mLeafs)
809 {
810 }
811 
812 template<typename DistTreeT>
813 void
815  const tbb::blocked_range<size_t>& range) const
816 {
817  typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
818 
819  for (size_t n = range.begin(); n < range.end(); ++n) {
820  iter = mLeafs.leaf(n).cbeginValueOn();
821  for (; iter; ++iter) {
822  evalVoxel(iter.getCoord());
823  }
824  }
825 }
826 
827 template<typename DistTreeT>
828 void
830 {
831  const DistValueT val = mDistAccessor.getValue(ijk), zeroVal(0.0);
832 
833  if(!(val < zeroVal)) return;
834 
835  Vec3d dir = getClosestPointDir(ijk), n_dir;
836  DistValueT n_val;
837  Coord n_ijk;
838 
839  // Check voxel-face adjacent neighbours.
840  for (Int32 n = 0; n < 26; ++n) {
841  n_ijk = ijk + util::COORD_OFFSETS[n];
842 
843  if (mIntersectionAccessor.isValueOn(n_ijk)) continue;
844  if (!mDistAccessor.probeValue(n_ijk, n_val)) continue;
845  if (n_val < zeroVal) continue;
846 
847  n_dir = getClosestPointDir(n_ijk);
848 
849  if (n_dir.dot(dir) > 0.0 ) {
850  const_cast<IntersectingVoxelSign<DistTreeT> *>(this)->
851  mDistAccessor.setValue(ijk, -val);
852  break;
853  }
854  }
855 }
856 
857 template<typename DistTreeT>
858 Vec3d
859 IntersectingVoxelSign<DistTreeT>::getClosestPointDir(const Coord& ijk) const
860 {
861  Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
862  const Vec4I& prim = (*mPolygonList)[mIndexAccessor.getValue(ijk)];
863  const std::vector<Vec3s>& points = *mPointList;
864 
865  // Evaluate first triangle
866  const Vec3d a(points[prim[0]]);
867  const Vec3d b(points[prim[1]]);
868  const Vec3d c(points[prim[2]]);
869 
870  Vec3d uvw;
871  Vec3d cpt = closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw);
872  Vec3d diff = voxelCenter - cpt;
873 
874  // Evaluate second triangle if quad.
875  if (prim[3] != util::INVALID_IDX) {
876  const Vec3d d(points[prim[3]]);
877 
878  cpt = closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw);
879  Vec3d tmpDiff = voxelCenter - cpt;
880 
881  if (tmpDiff.lengthSqr() < diff.lengthSqr()) {
882  diff = tmpDiff;
883  }
884  }
885 
886  diff.normalize();
887  return diff;
888 }
889 
890 
892 
893 
894 // IntersectingVoxelCleaner
897 template<typename DistTreeT>
899 {
900 public:
903  typedef typename DistTreeT::ValueType DistValueT;
905  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
907  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
911 
912  IntersectingVoxelCleaner(DistTreeT& distTree, IndexTreeT& indexTree,
913  StencilTreeT& intersectionTree, StencilArrayT& leafs);
914 
916 
917  void runParallel();
918  void runSerial();
919 
921  void operator()(const tbb::blocked_range<size_t>&) const;
922 
923 private:
924  void operator=(const IntersectingVoxelCleaner<DistTreeT>&) {}
925 
926  DistTreeT& mDistTree;
927  DistAccessorT mDistAccessor;
928 
929  IndexTreeT& mIndexTree;
930  IndexAccessorT mIndexAccessor;
931 
932  StencilTreeT& mIntersectionTree;
933  StencilAccessorT mIntersectionAccessor;
934  StencilArrayT& mLeafs;
935 };
936 
937 template<typename DistTreeT>
938 void
940 {
941  tbb::parallel_for(mLeafs.getRange(), *this);
942  mIntersectionTree.pruneInactive();
943 }
944 
945 template<typename DistTreeT>
946 void
948 {
949  (*this)(mLeafs.getRange());
950  mIntersectionTree.pruneInactive();
951 }
952 
953 template<typename DistTreeT>
955  DistTreeT& distTree,
956  IndexTreeT& indexTree,
957  StencilTreeT& intersectionTree,
958  StencilArrayT& leafs)
959  : mDistTree(distTree)
960  , mDistAccessor(mDistTree)
961  , mIndexTree(indexTree)
962  , mIndexAccessor(mIndexTree)
963  , mIntersectionTree(intersectionTree)
964  , mIntersectionAccessor(mIntersectionTree)
965  , mLeafs(leafs)
966 {
967 }
968 
969 template<typename DistTreeT>
972  : mDistTree(rhs.mDistTree)
973  , mDistAccessor(mDistTree)
974  , mIndexTree(rhs.mIndexTree)
975  , mIndexAccessor(mIndexTree)
976  , mIntersectionTree(rhs.mIntersectionTree)
977  , mIntersectionAccessor(mIntersectionTree)
978  , mLeafs(rhs.mLeafs)
979 {
980 }
981 
982 template<typename DistTreeT>
983 void
985  const tbb::blocked_range<size_t>& range) const
986 {
987  Coord ijk, m_ijk;
988  bool turnOff;
989  DistValueT value, bg = mDistTree.background();
990 
991  typename StencilTreeT::LeafNodeType::ValueOnCIter iter;
992 
993  for (size_t n = range.begin(); n < range.end(); ++n) {
994 
995  typename StencilTreeT::LeafNodeType& leaf = mLeafs.leaf(n);
996  iter = leaf.cbeginValueOn();
997 
998  for (; iter; ++iter) {
999 
1000  ijk = iter.getCoord();
1001 
1002  turnOff = true;
1003  for (Int32 m = 0; m < 26; ++m) {
1004  m_ijk = ijk + util::COORD_OFFSETS[m];
1005  if (mDistAccessor.probeValue(m_ijk, value)) {
1006  if (value > 0.0) {
1007  turnOff = false;
1008  break;
1009  }
1010  }
1011  }
1012 
1013  if (turnOff) leaf.setValueOff(ijk, bg);
1014  }
1015  }
1016 }
1017 
1018 
1020 
1021 
1022 // ShellVoxelCleaner
1025 template<typename DistTreeT>
1027 {
1028 public:
1031  typedef typename DistTreeT::ValueType DistValueT;
1034  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
1036  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1039 
1040  ShellVoxelCleaner(DistTreeT& distTree, DistArrayT& leafs, IndexTreeT& indexTree,
1041  StencilTreeT& intersectionTree);
1042 
1044 
1045  void runParallel();
1046  void runSerial();
1047 
1049  void operator()(const tbb::blocked_range<size_t>&) const;
1050 
1051 private:
1052  void operator=(const ShellVoxelCleaner<DistTreeT>&) {}
1053 
1054  DistTreeT& mDistTree;
1055  DistArrayT& mLeafs;
1056  DistAccessorT mDistAccessor;
1057 
1058  IndexTreeT& mIndexTree;
1059  IndexAccessorT mIndexAccessor;
1060 
1061  StencilTreeT& mIntersectionTree;
1062  StencilAccessorT mIntersectionAccessor;
1063 };
1064 
1065 template<typename DistTreeT>
1066 void
1068 {
1069  tbb::parallel_for(mLeafs.getRange(), *this);
1070  mDistTree.pruneInactive();
1071  mIndexTree.pruneInactive();
1072 }
1073 
1074 template<typename DistTreeT>
1075 void
1077 {
1078  (*this)(mLeafs.getRange());
1079  mDistTree.pruneInactive();
1080  mIndexTree.pruneInactive();
1081 }
1082 
1083 template<typename DistTreeT>
1085  DistTreeT& distTree,
1086  DistArrayT& leafs,
1087  IndexTreeT& indexTree,
1088  StencilTreeT& intersectionTree)
1089  : mDistTree(distTree)
1090  , mLeafs(leafs)
1091  , mDistAccessor(mDistTree)
1092  , mIndexTree(indexTree)
1093  , mIndexAccessor(mIndexTree)
1094  , mIntersectionTree(intersectionTree)
1095  , mIntersectionAccessor(mIntersectionTree)
1096 {
1097 }
1098 
1099 template<typename DistTreeT>
1101  const ShellVoxelCleaner<DistTreeT> &rhs)
1102  : mDistTree(rhs.mDistTree)
1103  , mLeafs(rhs.mLeafs)
1104  , mDistAccessor(mDistTree)
1105  , mIndexTree(rhs.mIndexTree)
1106  , mIndexAccessor(mIndexTree)
1107  , mIntersectionTree(rhs.mIntersectionTree)
1108  , mIntersectionAccessor(mIntersectionTree)
1109 {
1110 }
1111 
1112 template<typename DistTreeT>
1113 void
1115  const tbb::blocked_range<size_t>& range) const
1116 {
1117 
1118  Coord ijk, m_ijk;
1119  bool turnOff;
1120  DistValueT value;
1121 
1122  const DistValueT distC = -0.86602540378443861;
1123  const DistValueT distBG = mDistTree.background();
1124  const Int32 indexBG = mIntersectionTree.background();
1125 
1126  typename DistTreeT::LeafNodeType::ValueOnCIter iter;
1127  for (size_t n = range.begin(); n < range.end(); ++n) {
1128 
1129  typename DistTreeT::LeafNodeType& leaf = mLeafs.leaf(n);
1130  iter = leaf.cbeginValueOn();
1131  for (; iter; ++iter) {
1132 
1133  value = iter.getValue();
1134  if(value > 0.0) continue;
1135 
1136  ijk = iter.getCoord();
1137  if (mIntersectionAccessor.isValueOn(ijk)) continue;
1138 
1139  turnOff = true;
1140  for (Int32 m = 0; m < 18; ++m) {
1141  m_ijk = ijk + util::COORD_OFFSETS[m];
1142 
1143  if (mIntersectionAccessor.isValueOn(m_ijk)) {
1144  turnOff = false;
1145  break;
1146  }
1147  }
1148 
1149  if (turnOff) {
1150  leaf.setValueOff(ijk, distBG);
1151 
1152  const_cast<ShellVoxelCleaner<DistTreeT> *>(this)->
1153  mIndexAccessor.setValueOff(ijk, indexBG);
1154 
1155  } else {
1156  if (value > distC) leaf.setValue(ijk, distC);
1157  }
1158  }
1159  }
1160 }
1161 
1162 
1164 
1165 
1166 // ExpandNB
1169 template<typename DistTreeT>
1171 {
1172 public:
1175  typedef typename DistTreeT::ValueType DistValueT;
1177  typedef typename DistTreeT::template ValueConverter<Int32>::Type IndexTreeT;
1179  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1183 
1184  ExpandNB(const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
1185  DistTreeT& distTree, IndexTreeT& indexTree, StencilTreeT& maskTree, StencilArrayT& leafs,
1186  DistValueT exteriorBandWidth, DistValueT interiorBandWidth, DistValueT voxelSize);
1187 
1188  ExpandNB(const ExpandNB<DistTreeT>& rhs, tbb::split);
1189 
1191 
1192  void runParallel();
1193  void runSerial();
1194 
1195  void operator()(const tbb::blocked_range<size_t>&) const;
1196 
1197 private:
1198  void operator=(const ExpandNB<DistTreeT>&) {}
1199  double getDist(const Coord&, DistAccessorT&, IndexAccessorT&, Int32& primIndex) const;
1200  double getDistToPrim(const Coord& ijk, const Int32 polyIdx) const;
1201 
1202  std::vector<Vec3s> const * const mPointList;
1203  std::vector<Vec4I> const * const mPolygonList;
1204 
1205  DistTreeT& mDistTree;
1206  IndexTreeT& mIndexTree;
1207  StencilTreeT& mMaskTree;
1208  StencilArrayT& mLeafs;
1209 
1210  const DistValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1211 };
1212 
1213 template<typename DistTreeT>
1214 void
1216 {
1217  tbb::parallel_for(mLeafs.getRange(), *this);
1218  mMaskTree.pruneInactive();
1219 }
1220 
1221 template<typename DistTreeT>
1222 void
1224 {
1225  (*this)(mLeafs.getRange());
1226  mMaskTree.pruneInactive();
1227 }
1228 
1229 template<typename DistTreeT>
1231  const std::vector<Vec3s>& pointList,
1232  const std::vector<Vec4I>& polygonList,
1233  DistTreeT& distTree,
1234  IndexTreeT& indexTree,
1235  StencilTreeT& maskTree,
1236  StencilArrayT& leafs,
1237  DistValueT exteriorBandWidth, DistValueT interiorBandWidth,
1238  DistValueT voxelSize)
1239  : mPointList(&pointList)
1240  , mPolygonList(&polygonList)
1241  , mDistTree(distTree)
1242  , mIndexTree(indexTree)
1243  , mMaskTree(maskTree)
1244  , mLeafs(leafs)
1245  , mExteriorBandWidth(exteriorBandWidth)
1246  , mInteriorBandWidth(interiorBandWidth)
1247  , mVoxelSize(voxelSize)
1248 {
1249 }
1250 
1251 template<typename DistTreeT>
1253  : mPointList(rhs.mPointList)
1254  , mPolygonList(rhs.mPolygonList)
1255  , mDistTree(rhs.mDistTree)
1256  , mIndexTree(rhs.mIndexTree)
1257  , mMaskTree(rhs.mMaskTree)
1258  , mLeafs(rhs.mLeafs)
1259  , mExteriorBandWidth(rhs.mExteriorBandWidth)
1260  , mInteriorBandWidth(rhs.mInteriorBandWidth)
1261  , mVoxelSize(rhs.mVoxelSize)
1262 {
1263 }
1264 
1265 template<typename DistTreeT>
1266 void
1267 ExpandNB<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
1268 {
1269  typedef typename DistTreeT::LeafNodeType DistLeafT;
1270  typedef typename IndexTreeT::LeafNodeType IndexLeafT;
1271  typedef typename StencilTreeT::LeafNodeType StencilLeafT;
1272 
1273  Coord ijk, n_ijk;
1274  Int32 closestPrimIndex = 0;
1275 
1276  DistAccessorT distAcc(mDistTree);
1277  IndexAccessorT indexAcc(mIndexTree);
1278 
1279  for (size_t n = range.begin(); n < range.end(); ++n) {
1280 
1281  StencilLeafT& maskLeaf = mLeafs.leaf(n);
1282 
1283  const Coord origin = maskLeaf.getOrigin();
1284 
1285  DistLeafT* distLeafPt = distAcc.probeLeaf(origin);
1286  IndexLeafT* indexLeafPt = indexAcc.probeLeaf(origin);
1287 
1288  if (distLeafPt != NULL && indexLeafPt != NULL) {
1289 
1290  DistLeafT& distLeaf = *distLeafPt;
1291  IndexLeafT& indexLeaf = *indexLeafPt;
1292 
1293 
1294  typename StencilLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1295  for (; iter; ++iter) {
1296 
1297  const Index pos = iter.pos();
1298 
1299  if (distLeaf.isValueOn(pos)) {
1300  iter.setValueOff();
1301  continue;
1302  }
1303 
1305  getDist(iter.getCoord(), distAcc, indexAcc, closestPrimIndex);
1306 
1307  const bool inside = distLeaf.getValue(pos) < DistValueT(0.0);
1308 
1309  if (!inside && distance < mExteriorBandWidth) {
1310  distLeaf.setValueOn(pos, distance);
1311  indexLeaf.setValueOn(pos, closestPrimIndex);
1312  } else if (inside && distance < mInteriorBandWidth) {
1313  distLeaf.setValueOn(pos, -distance);
1314  indexLeaf.setValueOn(pos, closestPrimIndex);
1315  } else {
1316  iter.setValueOff();
1317  }
1318  }
1319 
1320  } else {
1321  maskLeaf.setValuesOff();
1322  }
1323  }
1324 }
1325 
1326 template<typename DistTreeT>
1327 double
1328 ExpandNB<DistTreeT>::getDist(const Coord& ijk, DistAccessorT& distAcc,
1329  IndexAccessorT& indexAcc, Int32& primIndex) const
1330 {
1331  Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1332  DistValueT nDist, dist = std::numeric_limits<double>::max();
1333 
1334  // Find neighbour with closest face point
1335  Coord n_ijk;
1336  for (Int32 n = 0; n < 18; ++n) {
1337  n_ijk = ijk + util::COORD_OFFSETS[n];
1338  if (distAcc.probeValue(n_ijk, nDist)) {
1339  nDist = std::abs(nDist);
1340  if (nDist < dist) {
1341  dist = nDist;
1342  primIndex = indexAcc.getValue(n_ijk);
1343  }
1344  }
1345  }
1346 
1347  // Calc. this voxels distance to the closest primitive.
1348  return getDistToPrim(ijk, primIndex);
1349 }
1350 
1351 
1352 template<typename DistTreeT>
1353 double
1354 ExpandNB<DistTreeT>::getDistToPrim(const Coord& ijk, const Int32 polyIdx) const
1355 {
1356  Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1357  const Vec4I& verts = (*mPolygonList)[polyIdx];
1358  const std::vector<Vec3s>& points = *mPointList;
1359 
1360  // Evaluate first triangle
1361  const Vec3d a(points[verts[0]]);
1362  const Vec3d b(points[verts[1]]);
1363  const Vec3d c(points[verts[2]]);
1364 
1365  Vec3d uvw;
1366  double dist = (voxelCenter -
1367  closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
1368 
1369  // Split-up quad into a second triangle and calac distance.
1370  if (util::INVALID_IDX != verts[3]) {
1371  const Vec3d d(points[verts[3]]);
1372 
1373  double secondDist = (voxelCenter -
1374  closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw)).lengthSqr();
1375 
1376  if (secondDist < dist) dist = secondDist;
1377  }
1378 
1379  return std::sqrt(dist) * double(mVoxelSize);
1380 }
1381 
1382 
1384 
1385 
1386 // Helper methods
1387 
1394 template<typename DistTreeT>
1395 inline void
1396 surfaceTracer(const Coord &seed, DistTreeT& distTree,
1397  typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree)
1398 {
1399  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1400  typedef typename tree::ValueAccessor<StencilTreeT> StencilAccessorT;
1401  typedef typename tree::ValueAccessor<DistTreeT> DistAccessorT;
1402  typedef typename DistTreeT::ValueType DistValueT;
1403 
1404  StencilAccessorT intrAccessor(intersectionTree);
1405  DistAccessorT distAccessor(distTree);
1406 
1407  std::deque<Coord> coordList;
1408  coordList.push_back(seed);
1409  Coord ijk, n_ijk;
1410 
1411  while (!coordList.empty()) {
1412  ijk = coordList.back();
1413  coordList.pop_back();
1414 
1415  if (!distAccessor.isValueOn(ijk)) continue;
1416 
1417  DistValueT& dist = const_cast<DistValueT&>(distAccessor.getValue(ijk));
1418  if (!(dist < 0.0)) continue;
1419  dist = -dist; // flip sign
1420 
1421  for (int n = 0; n < 6; ++n) {
1422  n_ijk = ijk + util::COORD_OFFSETS[n];
1423 
1424  if (!intrAccessor.isValueOn(n_ijk)) { // Don't cross the interface
1425  if (distAccessor.isValueOn(n_ijk)) { // Is part of the narrow band
1426  if (distAccessor.getValue(n_ijk) < 0.0) { // Marked as outside.
1427  coordList.push_back(n_ijk);
1428  }
1429  }
1430  }
1431 
1432  } // END neighbour voxel loop.
1433  } // END coordList loop.
1434 }
1435 
1436 
1444 template<typename DistTreeT, typename InterruptT>
1445 inline void
1446 propagateSign(DistTreeT& distTree,
1447  typename DistTreeT::template ValueConverter<bool>::Type& intersectionTree,
1448  InterruptT *interrupter = NULL)
1449 {
1450  typedef typename DistTreeT::template ValueConverter<bool>::Type StencilTreeT;
1451  typedef typename tree::ValueAccessor<StencilTreeT> StencilAccessorT;
1452  typedef typename tree::ValueAccessor<DistTreeT> DistAccessorT;
1453  typedef typename DistTreeT::ValueType DistValueT;
1454 
1455  StencilAccessorT intrAccessor(intersectionTree);
1456  DistAccessorT distAccessor(distTree);
1457  Coord ijk, n_ijk;
1458 
1459  typename DistTreeT::LeafIter leafIter = distTree.beginLeaf();
1460  for (; leafIter; leafIter.next()) {
1461 
1462  if (interrupter && interrupter->wasInterrupted()) break;
1463 
1464  typename DistTreeT::LeafNodeType::ValueOnIter iter = leafIter->beginValueOn();
1465  for (; iter; iter.next()) {
1466 
1467  ijk = iter.getCoord();
1468 
1469  // Ignore intersecting voxels.
1470  if (intrAccessor.isValueOn(ijk)) continue;
1471 
1472  if (iter.getValue() < 0.0) {
1473  for (Int32 n = 0; n < 6; ++n) {
1474  n_ijk = ijk + util::COORD_OFFSETS[n];
1475 
1476  if (distAccessor.isValueOn(n_ijk) && distAccessor.getValue(n_ijk) > 0.0) {
1477  surfaceTracer(ijk, distTree, intersectionTree);
1478  break;
1479  }
1480  }
1481  }
1482 
1483  } // END voxel iteration
1484  } // END leaf iteration
1485 }
1486 
1487 
1488 template<typename ValueType>
1490 {
1491  SqrtAndScaleOp(ValueType voxelSize, bool unsignedDist = false)
1492  : mVoxelSize(voxelSize)
1493  , mUnsigned(unsignedDist)
1494  {
1495  }
1496 
1497  template <typename LeafNodeType>
1498  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
1499  {
1500  ValueType w[2];
1501  w[0] = mVoxelSize;
1502  w[1] = -mVoxelSize;
1503 
1504  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1505  for (; iter; ++iter) {
1506  ValueType& val = const_cast<ValueType&>(iter.getValue());
1507  val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
1508  }
1509  }
1510 
1511 private:
1512  ValueType mVoxelSize;
1513  const bool mUnsigned;
1514 };
1515 
1516 
1517 template<typename ValueType>
1519 {
1520  VoxelSignOp(ValueType exBandWidth, ValueType inBandWidth)
1521  : mExBandWidth(exBandWidth)
1522  , mInBandWidth(inBandWidth)
1523  {
1524  }
1525 
1526  template <typename LeafNodeType>
1527  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
1528  {
1529  ValueType bgValues[2];
1530  bgValues[0] = mExBandWidth;
1531  bgValues[1] = -mInBandWidth;
1532 
1533  typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
1534 
1535  for (; iter; ++iter) {
1536  ValueType& val = const_cast<ValueType&>(iter.getValue());
1537  val = bgValues[int(val < ValueType(0.0))];
1538  }
1539  }
1540 
1541 private:
1542  ValueType mExBandWidth, mInBandWidth;
1543 };
1544 
1545 
1546 template<typename ValueType>
1547 struct TrimOp
1548 {
1549  TrimOp(ValueType exBandWidth, ValueType inBandWidth)
1550  : mExBandWidth(exBandWidth)
1551  , mInBandWidth(inBandWidth)
1552  {
1553  }
1554 
1555  template <typename LeafNodeType>
1556  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
1557  {
1558  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1559 
1560  for (; iter; ++iter) {
1561  const ValueType& val = iter.getValue();
1562  const bool inside = val < ValueType(0.0);
1563 
1564  if (inside && !(val > -mInBandWidth)) {
1565  iter.setValue(-mInBandWidth);
1566  iter.setValueOff();
1567  } else if (!inside && !(val < mInBandWidth)) {
1568  iter.setValue(mExBandWidth);
1569  iter.setValueOff();
1570  }
1571  }
1572  }
1573 
1574 private:
1575  ValueType mExBandWidth, mInBandWidth;
1576 };
1577 
1578 
1579 template<typename ValueType>
1580 struct OffsetOp
1581 {
1582  OffsetOp(ValueType offset): mOffset(offset) {}
1583 
1584  void resetOffset(ValueType offset) { mOffset = offset; }
1585 
1586  template <typename LeafNodeType>
1587  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
1588  {
1589  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1590  for (; iter; ++iter) {
1591  ValueType& val = const_cast<ValueType&>(iter.getValue());
1592  val += mOffset;
1593  }
1594  }
1595 
1596 private:
1597  ValueType mOffset;
1598 };
1599 
1600 
1601 template<typename GridType, typename ValueType>
1602 struct RenormOp
1603 {
1605  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
1608 
1609  RenormOp(GridType& grid, LeafManagerType& leafs, ValueType voxelSize, ValueType cfl = 1.0)
1610  : mGrid(grid)
1611  , mLeafs(leafs)
1612  , mVoxelSize(voxelSize)
1613  , mCFL(cfl)
1614  {
1615  }
1616 
1617  void resetCFL(ValueType cfl) { mCFL = cfl; }
1618 
1619  template <typename LeafNodeType>
1620  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1621  {
1622  const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
1623  Stencil stencil(mGrid);
1624 
1625  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
1626 
1627  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1628  for (; iter; ++iter) {
1629  stencil.moveTo(iter);
1630 
1631  const ValueType normSqGradPhi =
1633 
1634  const ValueType phi0 = stencil.getValue();
1635  const ValueType diff = math::Sqrt(normSqGradPhi) * invDx - one;
1636  const ValueType S = phi0 / (math::Sqrt(math::Pow2(phi0) + normSqGradPhi));
1637 
1638  buffer.setValue(iter.pos(), phi0 - dt * S * diff);
1639  }
1640  }
1641 
1642 private:
1643  GridType& mGrid;
1644  LeafManagerType& mLeafs;
1645  ValueType mVoxelSize, mCFL;
1646 };
1647 
1648 
1649 template<typename TreeType, typename ValueType>
1650 struct MinOp
1651 {
1654 
1655  MinOp(LeafManagerType& leafs): mLeafs(leafs) {}
1656 
1657  template <typename LeafNodeType>
1658  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1659  {
1660  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
1661 
1662  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1663  for (; iter; ++iter) {
1664  ValueType& val = const_cast<ValueType&>(iter.getValue());
1665  val = std::min(val, buffer.getValue(iter.pos()));
1666  }
1667  }
1668 
1669 private:
1670  LeafManagerType& mLeafs;
1671 };
1672 
1673 
1674 template<typename TreeType, typename ValueType>
1676 {
1679 
1680  MergeBufferOp(LeafManagerType& leafs, size_t bufferIndex = 1)
1681  : mLeafs(leafs)
1682  , mBufferIndex(bufferIndex)
1683  {
1684  }
1685 
1686  template <typename LeafNodeType>
1687  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1688  {
1689  BufferType& buffer = mLeafs.getBuffer(leafIndex, mBufferIndex);
1690 
1691  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1692  for (; iter; ++iter) {
1693  leaf.setValueOnly(iter.pos(), buffer.getValue(iter.pos()));
1694  }
1695  }
1696 
1697 private:
1698  LeafManagerType& mLeafs;
1699  const size_t mBufferIndex;
1700 };
1701 
1702 } // internal namespace
1703 
1704 
1706 
1707 
1708 // MeshToVolume
1709 
1710 template<typename DistGridT, typename InterruptT>
1712  openvdb::math::Transform::Ptr& transform, int conversionFlags,
1713  InterruptT *interrupter, int signSweeps)
1714  : mTransform(transform)
1715  , mConversionFlags(conversionFlags)
1716  , mSignSweeps(signSweeps)
1717  , mInterrupter(interrupter)
1718 {
1719  clear();
1720  mSignSweeps = std::min(mSignSweeps, 1);
1721 }
1722 
1723 
1724 template<typename DistGridT, typename InterruptT>
1725 void
1727 {
1728  mDistGrid = DistGridT::create(std::numeric_limits<DistValueT>::max());
1729  mIndexGrid = IndexGridT::create(Int32(util::INVALID_IDX));
1730  mIntersectingVoxelsGrid = StencilGridT::create(false);
1731 }
1732 
1733 
1734 template<typename DistGridT, typename InterruptT>
1735 inline void
1737  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
1738  DistValueT exBandWidth, DistValueT inBandWidth)
1739 {
1740  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
1741  exBandWidth = std::max(internal::Tolerance<DistValueT>::minNarrowBandWidth(), exBandWidth);
1742  inBandWidth = std::max(internal::Tolerance<DistValueT>::minNarrowBandWidth(), inBandWidth);
1743  const DistValueT vs = mTransform->voxelSize()[0];
1744  doConvert(pointList, polygonList, vs * exBandWidth, vs * inBandWidth);
1745  mDistGrid->setGridClass(GRID_LEVEL_SET);
1746 }
1747 
1748 
1749 template<typename DistGridT, typename InterruptT>
1750 inline void
1752  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
1753  DistValueT exBandWidth)
1754 {
1755  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
1756  exBandWidth = std::max(internal::Tolerance<DistValueT>::minNarrowBandWidth(), exBandWidth);
1757  const DistValueT vs = mTransform->voxelSize()[0];
1758  doConvert(pointList, polygonList, vs * exBandWidth, 0.0, true);
1759  mDistGrid->setGridClass(GRID_UNKNOWN);
1760 }
1761 
1762 template<typename DistGridT, typename InterruptT>
1763 void
1765  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
1766  DistValueT exBandWidth, DistValueT inBandWidth, bool unsignedDistField)
1767 {
1768  mDistGrid->setTransform(mTransform);
1769  mIndexGrid->setTransform(mTransform);
1770 
1771  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1772 
1773  // Voxelize mesh
1774  {
1776  voxelizer(pointList, polygonList, mInterrupter);
1777 
1778  voxelizer.runParallel();
1779 
1780  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1781 
1782  mDistGrid->tree().merge(voxelizer.sqrDistTree());
1783  mIndexGrid->tree().merge(voxelizer.primIndexTree());
1784  mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
1785  }
1786 
1787  if (!unsignedDistField) {
1788 
1789  // Determine the inside/outside state for the narrow band of voxels.
1790  {
1791  // Slices up the volume and label the exterior contour of each slice in parallel.
1792  internal::ContourTracer<DistTreeT, InterruptT> trace(
1793  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1794 
1795  for (int i = 0; i < mSignSweeps; ++i) {
1796 
1797  if (mInterrupter && mInterrupter->wasInterrupted()) break;
1798  trace.runParallel();
1799 
1800  if (mInterrupter && mInterrupter->wasInterrupted()) break;
1801 
1802  // Propagate sign information between the slices.
1803  internal::propagateSign<DistTreeT, InterruptT>
1804  (mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
1805  }
1806  }
1807 
1808  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1809 
1810  {
1811  tree::LeafManager<StencilTreeT> leafs(mIntersectingVoxelsGrid->tree());
1812 
1813  // Determine the sign of the mesh intersecting voxels.
1814  internal::IntersectingVoxelSign<DistTreeT> sign(pointList, polygonList,
1815  mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1816 
1817  sign.runParallel();
1818 
1819  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1820 
1821  // Remove mesh intersecting voxels that where set by rasterizing
1822  // self-intersecting portions of the mesh.
1823  internal::IntersectingVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1824  mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
1825 
1826  cleaner.runParallel();
1827  }
1828 
1829  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1830 
1831  {
1832  // Remove shell voxels that where set by rasterizing
1833  // self-intersecting portions of the mesh.
1834 
1835  tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1836 
1837  internal::ShellVoxelCleaner<DistTreeT> cleaner(mDistGrid->tree(),
1838  leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
1839 
1840  cleaner.runParallel();
1841  }
1842 
1843  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1844 
1845  } else { // if unsigned dist. field
1846  inBandWidth = DistValueT(0.0);
1847  }
1848 
1849  if (mDistGrid->activeVoxelCount() == 0) return;
1850 
1851  const DistValueT voxelSize(mTransform->voxelSize()[0]);
1852 
1853  // Transform values (world space scaling etc.)
1854  {
1855  typedef internal::SqrtAndScaleOp<DistValueT> XForm;
1856  tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1857  XForm op(voxelSize, unsignedDistField);
1858  LeafTransformer<DistTreeT, XForm> transform(leafs, op);
1859 
1860  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1861 
1862  transform.runParallel();
1863  }
1864 
1865  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1866 
1867  if (!unsignedDistField) {
1868  // Propagate sign information to inactive values.
1869  mDistGrid->tree().signedFloodFill();
1870 
1871  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1872 
1873  // Update the background value (inactive values)
1874  {
1875  tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
1876 
1877  typedef internal::VoxelSignOp<DistValueT> SignXForm;
1878  SignXForm op(exBandWidth, inBandWidth);
1879 
1880  LeafTransformer<DistTreeT, SignXForm> transform(leafs, op);
1881  transform.runParallel();
1882 
1883  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1884 
1885  DistValueT bgValues[2];
1886  bgValues[0] = exBandWidth;
1887  bgValues[1] = -inBandWidth;
1888 
1889  typename DistTreeT::ValueAllIter tileIt(mDistGrid->tree());
1890  tileIt.setMaxDepth(DistTreeT::ValueAllIter::LEAF_DEPTH - 1);
1891 
1892  for ( ; tileIt; ++tileIt) {
1893  DistValueT& val = const_cast<DistValueT&>(tileIt.getValue());
1894  val = bgValues[int(val < DistValueT(0.0))];
1895  }
1896 
1897  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1898 
1899  // fast bg value swap
1900  typename DistTreeT::Ptr newTree(new DistTreeT(/*background=*/exBandWidth));
1901  newTree->merge(mDistGrid->tree());
1902  mDistGrid->setTree(newTree);
1903  }
1904 
1905  // Smooth out bumps caused by self-intersecting and
1906  // overlapping portions of the mesh and renormalize the level set.
1907  {
1908  typedef internal::OffsetOp<DistValueT> OffsetOp;
1909  typedef internal::RenormOp<DistGridT, DistValueT> RenormOp;
1910  typedef internal::MinOp<DistTreeT, DistValueT> MinOp;
1911  typedef internal::MergeBufferOp<DistTreeT, DistValueT> MergeBufferOp;
1912 
1913  tree::LeafManager<DistTreeT> leafs(mDistGrid->tree(), 1);
1914 
1915  const DistValueT offset = 0.8 * voxelSize;
1916 
1917  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1918 
1919  OffsetOp offsetOp(-offset);
1920  LeafTransformer<DistTreeT, OffsetOp> offsetXform(leafs, offsetOp);
1921  offsetXform.runParallel();
1922 
1923  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1924 
1925  RenormOp renormOp(*mDistGrid, leafs, voxelSize);
1926  LeafTransformer<DistTreeT, RenormOp> renormXform(leafs, renormOp);
1927  renormXform.runParallel();
1928 
1929  MinOp minOp(leafs);
1930  LeafTransformer<DistTreeT, MinOp> minXform(leafs, minOp);
1931  minXform.runParallel();
1932 
1933  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1934 
1935  offsetOp.resetOffset(offset - internal::Tolerance<DistValueT>::epsilon());
1936  offsetXform.runParallel();
1937  }
1938 
1939  mIntersectingVoxelsGrid->clear();
1940  }
1941 
1942  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1943 
1944  // Narrow-band dilation
1945  const DistValueT minWidth = voxelSize * 2.0;
1946  if (inBandWidth > minWidth || exBandWidth > minWidth) {
1947 
1948  // Create the initial voxel mask.
1949  StencilTreeT maskTree(false);
1950  tree::ValueAccessor<StencilTreeT> acc(maskTree);
1951  maskTree.topologyUnion(mDistGrid->tree());
1952 
1953  // Preallocate leafs.
1954  {
1955  typedef typename DistTreeT::LeafNodeType DistLeafType;
1956 
1957  std::vector<DistLeafType*> distLeafs;
1958  distLeafs.reserve(mDistGrid->tree().leafCount());
1959 
1960  typename DistTreeT::LeafIter iter = mDistGrid->tree().beginLeaf();
1961  for ( ; iter; ++iter) distLeafs.push_back(iter.getLeaf());
1962 
1963  tree::ValueAccessor<DistTreeT> distAcc(mDistGrid->tree());
1964 
1965  DistValueT leafSize = DistValueT(DistLeafType::DIM - 1) * voxelSize;
1966 
1967  const double inLeafsRatio = double(inBandWidth) / double(leafSize);
1968  size_t inLeafs = std::numeric_limits<size_t>::max();
1969  if (double(inLeafs) > (inLeafsRatio + 1.0)) {
1970  inLeafs = size_t(std::ceil(inLeafsRatio)) + 1;
1971  }
1972  size_t exLeafs = size_t(std::ceil(exBandWidth / leafSize)) + 1;
1973  size_t numLeafs = std::max(inLeafs, exLeafs);
1974 
1975  for (size_t i = 0; i < numLeafs; ++i) {
1976 
1977  if (mInterrupter && mInterrupter->wasInterrupted()) return;
1978 
1979  std::vector<DistLeafType*> newDistLeafs;
1980  newDistLeafs.reserve(2 * distLeafs.size());
1981 
1982  for (size_t n = 0, N = distLeafs.size(); n < N; ++n) {
1983 
1984  Coord ijk = distLeafs[n]->getOrigin();
1985 
1986  const bool inside = distLeafs[n]->getValue(ijk) < DistValueT(0.0);
1987 
1988  if (inside && !(i < inLeafs)) continue;
1989  else if (!inside && !(i < exLeafs)) continue;
1990 
1991  ijk[0] -= 1;
1992  if (distAcc.probeLeaf(ijk) == NULL) {
1993  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
1994  }
1995 
1996  ijk[0] += 1;
1997  ijk[1] -= 1;
1998  if (distAcc.probeLeaf(ijk) == NULL) {
1999  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2000  }
2001 
2002  ijk[1] += 1;
2003  ijk[2] -= 1;
2004  if (distAcc.probeLeaf(ijk) == NULL) {
2005  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2006  }
2007 
2008  ijk[2] += 1;
2009  ijk[0] += DistLeafType::DIM;
2010  if (distAcc.probeLeaf(ijk) == NULL) {
2011  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2012  }
2013 
2014  ijk[0] -= DistLeafType::DIM;
2015  ijk[1] += DistLeafType::DIM;
2016  if (distAcc.probeLeaf(ijk) == NULL) {
2017  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2018  }
2019 
2020  ijk[1] -= DistLeafType::DIM;
2021  ijk[2] += DistLeafType::DIM;
2022  if (distAcc.probeLeaf(ijk) == NULL) {
2023  newDistLeafs.push_back(distAcc.touchLeaf(ijk));
2024  }
2025  }
2026 
2027  if (newDistLeafs.empty()) break;
2028  distLeafs.swap(newDistLeafs);
2029  }
2030  }
2031 
2032  if (mInterrupter && mInterrupter->wasInterrupted()) return;
2033 
2034  mIndexGrid->tree().topologyUnion(mDistGrid->tree());
2035 
2036  while (maskTree.activeVoxelCount() > 0) {
2037 
2038  if (mInterrupter && mInterrupter->wasInterrupted()) break;
2039 
2040  openvdb::tools::dilateVoxels(maskTree);
2041  tree::LeafManager<StencilTreeT> leafs(maskTree);
2042 
2043  internal::ExpandNB<DistTreeT> expand(pointList, polygonList, mDistGrid->tree(),
2044  mIndexGrid->tree(), maskTree, leafs, exBandWidth, inBandWidth, voxelSize);
2045 
2046  expand.runParallel();
2047  }
2048  }
2049 
2050  if (!bool(GENERATE_PRIM_INDEX_GRID & mConversionFlags)) mIndexGrid->clear();
2051 
2052  const DistValueT minTrimWidth = voxelSize * 4.0;
2053  if (inBandWidth < minTrimWidth || exBandWidth < minTrimWidth) {
2054 
2055  // If the narrow band was not expanded, we might need to trim off
2056  // some of the active voxels in order to respect the narrow band limits.
2057  // (The mesh voxelization step generates some extra 'shell' voxels)
2058 
2059  tree::LeafManager<DistTreeT> leafs(mDistGrid->tree());
2060 
2061  typedef internal::TrimOp<DistValueT> TrimOp;
2062 
2063  TrimOp op(exBandWidth, unsignedDistField ? exBandWidth : inBandWidth);
2064  LeafTransformer<DistTreeT, TrimOp> transform(leafs, op);
2065  transform.runParallel();
2066  }
2067 
2068  if (mInterrupter && mInterrupter->wasInterrupted()) return;
2069 
2070  mDistGrid->tree().pruneLevelSet();
2071 }
2072 
2073 } // namespace tools
2074 } // namespace OPENVDB_VERSION_NAME
2075 } // namespace openvdb
2076 
2077 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
2078 
2079 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2080 // All rights reserved. This software is distributed under the
2081 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )