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