OpenVDB  3.2.0
LevelSetUtil.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2016 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 //
37 
38 
39 #ifndef OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
41 
42 #include "MeshToVolume.h" // for traceExteriorBoundaries
43 #include "SignedFloodFill.h" // for signedFloodFillWithValues
44 
45 #include <openvdb/Types.h>
46 #include <openvdb/Grid.h>
47 
48 #include <tbb/blocked_range.h>
49 #include <tbb/parallel_for.h>
50 #include <tbb/parallel_reduce.h>
51 #include <tbb/parallel_sort.h>
52 
53 #include <limits>
54 
55 namespace openvdb {
57 namespace OPENVDB_VERSION_NAME {
58 namespace tools {
59 
60 // MS Visual C++ requires this extra level of indirection in order to compile
61 // THIS MUST EXIST IN AN UNNAMED NAMESPACE IN ORDER TO COMPILE ON WINDOWS
62 namespace {
63 
64 template<typename GridType>
65 inline typename GridType::ValueType lsutilGridMax()
66 {
68 }
69 
70 template<typename GridType>
71 inline typename GridType::ValueType lsutilGridZero()
72 {
73  return zeroVal<typename GridType::ValueType>();
74 }
75 
76 } // unnamed namespace
77 
78 
80 
81 
96 template<class GridType>
97 inline void
99  GridType& grid,
100  typename GridType::ValueType cutoffDistance = lsutilGridMax<GridType>());
101 
102 
113 template<class GridOrTreeType>
114 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
116  const GridOrTreeType& volume,
117  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>());
118 
119 
140 template<typename GridOrTreeType>
141 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
142 extractEnclosedRegion(const GridOrTreeType& volume,
143  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>(),
144  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type* fillMask = NULL);
145 
146 
151 template<typename GridOrTreeType>
152 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
153 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue);
154 
155 
161 template<typename GridOrTreeType>
162 inline void
163 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
164  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks);
165 
166 
174 template<typename GridOrTreeType>
175 inline void
176 segmentActiveVoxels(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
177 
178 
187 template<typename GridOrTreeType>
188 inline void
189 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
190 
191 
194 
195 // Internal utility objects and implementation details
196 
197 
198 namespace level_set_util_internal {
199 
200 
201 template<typename LeafNodeType>
203 
204  typedef typename LeafNodeType::ValueType ValueType;
206 
208  ValueType isovalue, const LeafNodeType ** nodes, BoolLeafNodeType ** maskNodes)
209  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
210  {
211  }
212 
213  void operator()(const tbb::blocked_range<size_t>& range) const {
214 
215  BoolLeafNodeType * maskNodePt = NULL;
216 
217  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
218 
219  mMaskNodes[n] = NULL;
220  const LeafNodeType& node = *mNodes[n];
221 
222  if (!maskNodePt) {
223  maskNodePt = new BoolLeafNodeType(node.origin(), false);
224  } else {
225  maskNodePt->setOrigin(node.origin());
226  }
227 
228  const ValueType* values = &node.getValue(0);
229  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
230  if (values[i] < mIsovalue) maskNodePt->setValueOn(i, true);
231  }
232 
233  if (maskNodePt->onVoxelCount() > 0) {
234  mMaskNodes[n] = maskNodePt;
235  maskNodePt = NULL;
236  }
237  }
238 
239  if (maskNodePt) delete maskNodePt;
240  }
241 
242  LeafNodeType const * const * const mNodes;
243  BoolLeafNodeType ** const mMaskNodes;
244  ValueType const mIsovalue;
245 }; // MaskInteriorVoxels
246 
247 
248 template<typename TreeType, typename InternalNodeType>
250 
251  typedef typename TreeType::ValueType ValueType;
252 
253  MaskInteriorTiles(ValueType isovalue, const TreeType& tree, InternalNodeType ** maskNodes)
254  : mTree(&tree), mMaskNodes(maskNodes), mIsovalue(isovalue) { }
255 
256  void operator()(const tbb::blocked_range<size_t>& range) const {
258  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
259  typename InternalNodeType::ValueAllIter it = mMaskNodes[n]->beginValueAll();
260  for (; it; ++it) {
261  if (acc.getValue(it.getCoord()) < mIsovalue) {
262  it.setValue(true);
263  it.setValueOn(true);
264  }
265  }
266  }
267  }
268 
269  TreeType const * const mTree;
270  InternalNodeType ** const mMaskNodes;
271  ValueType const mIsovalue;
272 }; // MaskInteriorTiles
273 
274 
275 template<typename TreeType>
276 struct PopulateTree {
277 
278  typedef typename TreeType::ValueType ValueType;
279  typedef typename TreeType::LeafNodeType LeafNodeType;
280 
281  PopulateTree(TreeType& tree, LeafNodeType** leafnodes,
282  const size_t * nodexIndexMap, ValueType background)
283  : mNewTree(background)
284  , mTreePt(&tree)
285  , mNodes(leafnodes)
286  , mNodeIndexMap(nodexIndexMap)
287  {
288  }
289 
290  PopulateTree(PopulateTree& rhs, tbb::split)
291  : mNewTree(rhs.mNewTree.background())
292  , mTreePt(&mNewTree)
293  , mNodes(rhs.mNodes)
294  , mNodeIndexMap(rhs.mNodeIndexMap)
295  {
296  }
297 
298  void operator()(const tbb::blocked_range<size_t>& range) {
299 
300  tree::ValueAccessor<TreeType> acc(*mTreePt);
301 
302  if (mNodeIndexMap) {
303  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
304  for (size_t i = mNodeIndexMap[n], I = mNodeIndexMap[n + 1]; i < I; ++i) {
305  if (mNodes[i] != NULL) acc.addLeaf(mNodes[i]);
306  }
307  }
308  } else {
309  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
310  acc.addLeaf(mNodes[n]);
311  }
312  }
313  }
314 
315  void join(PopulateTree& rhs) { mTreePt->merge(*rhs.mTreePt); }
316 
317 private:
318  TreeType mNewTree;
319  TreeType * const mTreePt;
320  LeafNodeType ** const mNodes;
321  size_t const * const mNodeIndexMap;
322 }; // PopulateTree
323 
324 
326 template<typename LeafNodeType>
328 
329  typedef typename LeafNodeType::ValueType ValueType;
331 
333  ValueType isovalue, const LeafNodeType ** nodes, CharLeafNodeType ** maskNodes)
334  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
335  {
336  }
337 
338  void operator()(const tbb::blocked_range<size_t>& range) const {
339 
340  CharLeafNodeType * maskNodePt = NULL;
341 
342  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
343 
344  mMaskNodes[n] = NULL;
345  const LeafNodeType& node = *mNodes[n];
346 
347  if (!maskNodePt) {
348  maskNodePt = new CharLeafNodeType(node.origin(), 1);
349  } else {
350  maskNodePt->setOrigin(node.origin());
351  }
352 
353  typename LeafNodeType::ValueOnCIter it;
354  for (it = node.cbeginValueOn(); it; ++it) {
355  maskNodePt->setValueOn(it.pos(), ((*it - mIsovalue) < 0.0) ? 0 : 1);
356  }
357 
358  if (maskNodePt->onVoxelCount() > 0) {
359  mMaskNodes[n] = maskNodePt;
360  maskNodePt = NULL;
361  }
362  }
363 
364  if (maskNodePt) delete maskNodePt;
365  }
366 
367  LeafNodeType const * const * const mNodes;
368  CharLeafNodeType ** const mMaskNodes;
369  ValueType const mIsovalue;
370 }; // LabelBoundaryVoxels
371 
372 
373 template<typename LeafNodeType>
375  typedef typename LeafNodeType::ValueType ValueType;
376 
377  FlipRegionSign(LeafNodeType ** nodes) : mNodes(nodes) { }
378 
379  void operator()(const tbb::blocked_range<size_t>& range) const {
380  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
381  ValueType* values = const_cast<ValueType*>(&mNodes[n]->getValue(0));
382  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
383  values[i] = values[i] < 0 ? 1 : -1;
384  }
385  }
386  }
387 
388  LeafNodeType ** const mNodes;
389 }; // FlipRegionSign
390 
391 
392 template<typename LeafNodeType>
394 
395  typedef typename LeafNodeType::ValueType ValueType;
396 
397  FindMinVoxelValue(LeafNodeType const * const * const leafnodes)
398  : minValue(std::numeric_limits<ValueType>::max())
399  , mNodes(leafnodes)
400  {
401  }
402 
404  : minValue(std::numeric_limits<ValueType>::max())
405  , mNodes(rhs.mNodes)
406  {
407  }
408 
409  void operator()(const tbb::blocked_range<size_t>& range) {
410  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
411  const ValueType* data = mNodes[n]->buffer().data();
412  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
413  minValue = std::min(minValue, data[i]);
414  }
415  }
416  }
417 
418  void join(FindMinVoxelValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
419 
420  ValueType minValue;
421 
422  LeafNodeType const * const * const mNodes;
423 }; // FindMinVoxelValue
424 
425 
426 template<typename InternalNodeType>
428 
429  typedef typename InternalNodeType::ValueType ValueType;
430 
431  FindMinTileValue(InternalNodeType const * const * const nodes)
432  : minValue(std::numeric_limits<ValueType>::max())
433  , mNodes(nodes)
434  {
435  }
436 
438  : minValue(std::numeric_limits<ValueType>::max())
439  , mNodes(rhs.mNodes)
440  {
441  }
442 
443  void operator()(const tbb::blocked_range<size_t>& range) {
444  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
445  typename InternalNodeType::ValueAllCIter it = mNodes[n]->beginValueAll();
446  for (; it; ++it) {
447  minValue = std::min(minValue, *it);
448  }
449  }
450  }
451 
452  void join(FindMinTileValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
453 
454  ValueType minValue;
455 
456  InternalNodeType const * const * const mNodes;
457 }; // FindMinTileValue
458 
459 
460 template<typename LeafNodeType>
462 
463  typedef typename LeafNodeType::ValueType ValueType;
464 
465  SDFVoxelsToFogVolume(LeafNodeType ** nodes, ValueType cutoffDistance)
466  : mNodes(nodes), mWeight(ValueType(1.0) / cutoffDistance)
467  {
468  }
469 
470  void operator()(const tbb::blocked_range<size_t>& range) const {
471 
472  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
473 
474  LeafNodeType& node = *mNodes[n];
475  node.setValuesOff();
476 
477  ValueType* values = node.buffer().data();
478  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
479  values[i] = values[i] > ValueType(0.0) ? ValueType(0.0) : values[i] * mWeight;
480  if (values[i] > ValueType(0.0)) node.setValueOn(i);
481  }
482 
483  if (node.onVoxelCount() == 0) {
484  delete mNodes[n];
485  mNodes[n] = NULL;
486  }
487  }
488  }
489 
490  LeafNodeType ** const mNodes;
491  ValueType const mWeight;
492 }; // SDFVoxelsToFogVolume
493 
494 
495 template<typename TreeType, typename InternalNodeType>
497 
498  SDFTilesToFogVolume(const TreeType& tree, InternalNodeType ** nodes)
499  : mTree(&tree), mNodes(nodes) { }
500 
501  void operator()(const tbb::blocked_range<size_t>& range) const {
502 
503  typedef typename TreeType::ValueType ValueType;
505 
506  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
507  typename InternalNodeType::ValueAllIter it = mNodes[n]->beginValueAll();
508  for (; it; ++it) {
509  if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
510  it.setValue(ValueType(1.0));
511  it.setValueOn(true);
512  }
513  }
514  }
515  }
516 
517  TreeType const * const mTree;
518  InternalNodeType ** const mNodes;
519 }; // SDFTilesToFogVolume
520 
521 
522 template<typename TreeType>
524 
525  typedef typename TreeType::ValueType ValueType;
526  typedef typename TreeType::LeafNodeType LeafNodeType;
527  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
528  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
529 
530  FillMaskBoundary(const TreeType& tree, ValueType isovalue, const BoolTreeType& fillMask,
531  const BoolLeafNodeType ** fillNodes, BoolLeafNodeType ** newNodes)
532  : mTree(&tree), mFillMask(&fillMask), mFillNodes(fillNodes), mNewNodes(newNodes), mIsovalue(isovalue)
533  {
534  }
535 
536  void operator()(const tbb::blocked_range<size_t>& range) const {
537 
538  tree::ValueAccessor<const BoolTreeType> maskAcc(*mFillMask);
539  tree::ValueAccessor<const TreeType> distAcc(*mTree);
540 
541  boost::scoped_array<char> valueMask(new char[BoolLeafNodeType::SIZE]);
542 
543  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
544 
545  mNewNodes[n] = NULL;
546  const BoolLeafNodeType& node = *mFillNodes[n];
547  const Coord& origin = node.origin();
548 
549  const bool denseNode = node.isDense();
550 
551  // possible early out if the fill mask is dense
552  if (denseNode) {
553 
554  int denseNeighbors = 0;
555 
556  const BoolLeafNodeType* neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(-1, 0, 0));
557  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
558 
559  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(BoolLeafNodeType::DIM, 0, 0));
560  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
561 
562  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, -1, 0));
563  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
564 
565  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, BoolLeafNodeType::DIM, 0));
566  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
567 
568  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, -1));
569  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
570 
571  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, BoolLeafNodeType::DIM));
572  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
573 
574  if (denseNeighbors == 6) continue;
575  }
576 
577  // rest value mask
578  memset(valueMask.get(), 0, sizeof(char) * BoolLeafNodeType::SIZE);
579 
580  const typename TreeType::LeafNodeType* distNode = distAcc.probeConstLeaf(origin);
581 
582  // check internal voxel neighbors
583 
584  bool earlyTermination = false;
585 
586  if (!denseNode) {
587  if (distNode) {
588  evalInternalNeighborsP(valueMask.get(), node, *distNode);
589  evalInternalNeighborsN(valueMask.get(), node, *distNode);
590  } else if (distAcc.getValue(origin) > mIsovalue) {
591  earlyTermination = evalInternalNeighborsP(valueMask.get(), node);
592  if (!earlyTermination) earlyTermination = evalInternalNeighborsN(valueMask.get(), node);
593  }
594  }
595 
596  // check external voxel neighbors
597 
598  if (!earlyTermination) {
599  evalExternalNeighborsX<true>(valueMask.get(), node, maskAcc, distAcc);
600  evalExternalNeighborsX<false>(valueMask.get(), node, maskAcc, distAcc);
601  evalExternalNeighborsY<true>(valueMask.get(), node, maskAcc, distAcc);
602  evalExternalNeighborsY<false>(valueMask.get(), node, maskAcc, distAcc);
603  evalExternalNeighborsZ<true>(valueMask.get(), node, maskAcc, distAcc);
604  evalExternalNeighborsZ<false>(valueMask.get(), node, maskAcc, distAcc);
605  }
606 
607  // Export marked boundary voxels.
608 
609  int numBoundaryValues = 0;
610  for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
611  numBoundaryValues += valueMask[i] == 1;
612  }
613 
614  if (numBoundaryValues > 0) {
615  mNewNodes[n] = new BoolLeafNodeType(origin, false);
616  for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
617  if (valueMask[i] == 1) mNewNodes[n]->setValueOn(i);
618  }
619  }
620  }
621  }
622 
623 private:
624 
625  // Check internal voxel neighbors in positive {x, y, z} directions.
626 
627  void evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node, const LeafNodeType& distNode) const {
628 
629  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
630  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
631  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
632  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
633  for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
634  const Index pos = yPos + z;
635 
636  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
637 
638  if (!node.isValueOn(pos + 1) && distNode.getValue(pos + 1) > mIsovalue) {
639  valueMask[pos] = 1;
640  }
641  }
642  }
643  }
644 
645  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
646  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
647  for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
648  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
649  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
650  const Index pos = yPos + z;
651 
652  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
653 
654  if (!node.isValueOn(pos + BoolLeafNodeType::DIM) &&
655  distNode.getValue(pos + BoolLeafNodeType::DIM) > mIsovalue) {
656  valueMask[pos] = 1;
657  }
658  }
659  }
660  }
661 
662  for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
663  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
664  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
665  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
666  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
667  const Index pos = yPos + z;
668 
669  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
670 
671  if (!node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
672  distNode.getValue(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) > mIsovalue) {
673  valueMask[pos] = 1;
674  }
675  }
676  }
677  }
678  }
679 
680  bool evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node) const {
681 
682  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
683  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
684  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
685  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
686  for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
687  const Index pos = yPos + z;
688 
689  if (node.isValueOn(pos) && !node.isValueOn(pos + 1)) {
690  valueMask[pos] = 1;
691  return true;
692  }
693  }
694  }
695  }
696 
697  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
698  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
699  for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
700  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
701  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
702  const Index pos = yPos + z;
703 
704  if (node.isValueOn(pos) && !node.isValueOn(pos + BoolLeafNodeType::DIM)) {
705  valueMask[pos] = 1;
706  return true;
707  }
708  }
709  }
710  }
711 
712  for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
713  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
714  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
715  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
716  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
717  const Index pos = yPos + z;
718 
719  if (node.isValueOn(pos) &&
720  !node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
721  valueMask[pos] = 1;
722  return true;
723  }
724  }
725  }
726  }
727 
728  return false;
729  }
730 
731  // Check internal voxel neighbors in negative {x, y, z} directions.
732 
733  void evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node, const LeafNodeType& distNode) const {
734 
735  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
736  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
737  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
738  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
739  for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
740  const Index pos = yPos + z;
741 
742  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
743 
744  if (!node.isValueOn(pos - 1) && distNode.getValue(pos - 1) > mIsovalue) {
745  valueMask[pos] = 1;
746  }
747  }
748  }
749  }
750 
751  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
752  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
753  for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
754  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
755  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
756  const Index pos = yPos + z;
757 
758  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
759 
760  if (!node.isValueOn(pos - BoolLeafNodeType::DIM) &&
761  distNode.getValue(pos - BoolLeafNodeType::DIM) > mIsovalue) {
762  valueMask[pos] = 1;
763  }
764  }
765  }
766  }
767 
768  for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
769  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
770  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
771  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
772  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
773  const Index pos = yPos + z;
774 
775  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
776 
777  if (!node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
778  distNode.getValue(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) > mIsovalue) {
779  valueMask[pos] = 1;
780  }
781  }
782  }
783  }
784  }
785 
786 
787  bool evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node) const {
788 
789  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
790  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
791  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
792  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
793  for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
794  const Index pos = yPos + z;
795 
796  if (node.isValueOn(pos) && !node.isValueOn(pos - 1)) {
797  valueMask[pos] = 1;
798  return true;
799  }
800  }
801  }
802  }
803 
804  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
805  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
806  for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
807  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
808  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
809  const Index pos = yPos + z;
810 
811  if (node.isValueOn(pos) && !node.isValueOn(pos - BoolLeafNodeType::DIM)) {
812  valueMask[pos] = 1;
813  return true;
814  }
815  }
816  }
817  }
818 
819  for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
820  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
821  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
822  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
823  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
824  const Index pos = yPos + z;
825 
826  if (node.isValueOn(pos) &&
827  !node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
828  valueMask[pos] = 1;
829  return true;
830  }
831  }
832  }
833  }
834 
835  return false;
836  }
837 
838 
839  // Check external voxel neighbors
840 
841  // If UpWind is true check the X+ oriented node face, else the X- oriented face.
842  template<bool UpWind>
843  void evalExternalNeighborsX(char* valueMask, const BoolLeafNodeType& node,
845  const tree::ValueAccessor<const TreeType>& distAcc) const {
846 
847  const Coord& origin = node.origin();
848  Coord ijk(0, 0, 0), nijk;
849  int step = -1;
850 
851  if (UpWind) {
852  step = 1;
853  ijk[0] = int(BoolLeafNodeType::DIM) - 1;
854  }
855 
856  const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
857 
858  for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
859  const Index yPos = xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
860 
861  for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
862  const Index pos = yPos + ijk[2];
863 
864  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
865 
866  nijk = origin + ijk.offsetBy(step, 0, 0);
867 
868  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
869  valueMask[pos] = 1;
870  }
871  }
872  }
873  }
874  }
875 
876  // If UpWind is true check the Y+ oriented node face, else the Y- oriented face.
877  template<bool UpWind>
878  void evalExternalNeighborsY(char* valueMask, const BoolLeafNodeType& node,
880  const tree::ValueAccessor<const TreeType>& distAcc) const {
881 
882  const Coord& origin = node.origin();
883  Coord ijk(0, 0, 0), nijk;
884  int step = -1;
885 
886  if (UpWind) {
887  step = 1;
888  ijk[1] = int(BoolLeafNodeType::DIM) - 1;
889  }
890 
891  const Index yPos = ijk[1] << int(BoolLeafNodeType::LOG2DIM);
892 
893  for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
894  const Index xPos = yPos + (ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM)));
895 
896  for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
897  const Index pos = xPos + ijk[2];
898 
899  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
900 
901  nijk = origin + ijk.offsetBy(0, step, 0);
902  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
903  valueMask[pos] = 1;
904  }
905  }
906  }
907  }
908  }
909 
910  // If UpWind is true check the Z+ oriented node face, else the Z- oriented face.
911  template<bool UpWind>
912  void evalExternalNeighborsZ(char* valueMask, const BoolLeafNodeType& node,
914  const tree::ValueAccessor<const TreeType>& distAcc) const {
915 
916  const Coord& origin = node.origin();
917  Coord ijk(0, 0, 0), nijk;
918  int step = -1;
919 
920  if (UpWind) {
921  step = 1;
922  ijk[2] = int(BoolLeafNodeType::DIM) - 1;
923  }
924 
925  for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
926  const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
927 
928  for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
929  const Index pos = ijk[2] + xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
930 
931  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
932 
933  nijk = origin + ijk.offsetBy(0, 0, step);
934  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
935  valueMask[pos] = 1;
936  }
937  }
938  }
939  }
940  }
941 
943 
944  TreeType const * const mTree;
945  BoolTreeType const * const mFillMask;
946  BoolLeafNodeType const * const * const mFillNodes;
947  BoolLeafNodeType ** const mNewNodes;
948  ValueType const mIsovalue;
949 }; // FillMaskBoundary
950 
951 
954 template <class TreeType>
955 inline typename TreeType::template ValueConverter<char>::Type::Ptr
956 computeEnclosedRegionMask(const TreeType& tree, typename TreeType::ValueType isovalue,
957  const typename TreeType::template ValueConverter<bool>::Type* fillMask)
958 {
959  typedef typename TreeType::LeafNodeType LeafNodeType;
960  typedef typename TreeType::RootNodeType RootNodeType;
961  typedef typename RootNodeType::NodeChainType NodeChainType;
962  typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
963 
964  typedef typename TreeType::template ValueConverter<char>::Type CharTreeType;
965  typedef typename CharTreeType::LeafNodeType CharLeafNodeType;
966  typedef typename CharTreeType::RootNodeType CharRootNodeType;
967  typedef typename CharRootNodeType::NodeChainType CharNodeChainType;
968 
969  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
970  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
971 
973 
974  const TreeType* treePt = &tree;
975 
976  size_t numLeafNodes = 0, numInternalNodes = 0;
977 
978  std::vector<const LeafNodeType*> nodes;
979  std::vector<size_t> leafnodeCount;
980 
981  {
982  // compute the prefix sum of the leafnode count in each internal node.
983  std::vector<const InternalNodeType*> internalNodes;
984  treePt->getNodes(internalNodes);
985 
986  numInternalNodes = internalNodes.size();
987 
988  leafnodeCount.push_back(0);
989  for (size_t n = 0; n < numInternalNodes; ++n) {
990  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
991  }
992 
993  numLeafNodes = leafnodeCount.back();
994 
995  // extract all leafnodes
996  nodes.reserve(numLeafNodes);
997 
998  for (size_t n = 0; n < numInternalNodes; ++n) {
999  internalNodes[n]->getNodes(nodes);
1000  }
1001  }
1002 
1003  // create mask leafnodes
1004  boost::scoped_array<CharLeafNodeType*> maskNodes(new CharLeafNodeType*[numLeafNodes]);
1005 
1006  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1007  LabelBoundaryVoxels<LeafNodeType>(isovalue, &nodes[0], maskNodes.get()));
1008 
1009  // create mask grid
1010  typename CharTreeType::Ptr maskTree(new CharTreeType(1));
1011 
1012  PopulateTree<CharTreeType> populate(*maskTree, maskNodes.get(), &leafnodeCount[0], 1);
1013  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1014 
1015  // optionally evaluate the fill mask
1016 
1017  std::vector<CharLeafNodeType*> extraMaskNodes;
1018 
1019  if (fillMask) {
1020 
1021  std::vector<const BoolLeafNodeType*> fillMaskNodes;
1022  fillMask->getNodes(fillMaskNodes);
1023 
1024  boost::scoped_array<BoolLeafNodeType*> boundaryMaskNodes(new BoolLeafNodeType*[fillMaskNodes.size()]);
1025 
1026  tbb::parallel_for(tbb::blocked_range<size_t>(0, fillMaskNodes.size()),
1027  FillMaskBoundary<TreeType>(tree, isovalue, *fillMask, &fillMaskNodes[0], boundaryMaskNodes.get()));
1028 
1029  tree::ValueAccessor<CharTreeType> maskAcc(*maskTree);
1030 
1031  for (size_t n = 0, N = fillMaskNodes.size(); n < N; ++n) {
1032 
1033  if (boundaryMaskNodes[n] == NULL) continue;
1034 
1035  const BoolLeafNodeType& boundaryNode = *boundaryMaskNodes[n];
1036  const Coord& origin = boundaryNode.origin();
1037 
1038  CharLeafNodeType* maskNodePt = maskAcc.probeLeaf(origin);
1039 
1040  if (!maskNodePt) {
1041  maskNodePt = maskAcc.touchLeaf(origin);
1042  extraMaskNodes.push_back(maskNodePt);
1043  }
1044 
1045  char* data = maskNodePt->buffer().data();
1046 
1047  typename BoolLeafNodeType::ValueOnCIter it = boundaryNode.cbeginValueOn();
1048  for (; it; ++it) {
1049  if (data[it.pos()] != 0) data[it.pos()] = -1;
1050  }
1051 
1052  delete boundaryMaskNodes[n];
1053  }
1054  }
1055 
1056  // eliminate enclosed regions
1057  tools::traceExteriorBoundaries(*maskTree);
1058 
1059  // flip voxel sign to negative inside and positive outside.
1060  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1061  FlipRegionSign<CharLeafNodeType>(maskNodes.get()));
1062 
1063  if (!extraMaskNodes.empty()) {
1064  tbb::parallel_for(tbb::blocked_range<size_t>(0, extraMaskNodes.size()),
1065  FlipRegionSign<CharLeafNodeType>(&extraMaskNodes[0]));
1066  }
1067 
1068  // propagate sign information into tile region
1069  tools::signedFloodFill(*maskTree);
1070 
1071  return maskTree;
1072 } // computeEnclosedRegionMask()
1073 
1074 
1075 template <class TreeType>
1076 inline typename TreeType::template ValueConverter<bool>::Type::Ptr
1077 computeInteriorMask(const TreeType& tree, typename TreeType::ValueType iso)
1078 {
1079  typedef typename TreeType::LeafNodeType LeafNodeType;
1080  typedef typename TreeType::RootNodeType RootNodeType;
1081  typedef typename RootNodeType::NodeChainType NodeChainType;
1082  typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
1083 
1084  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
1085  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
1086  typedef typename BoolTreeType::RootNodeType BoolRootNodeType;
1087  typedef typename BoolRootNodeType::NodeChainType BoolNodeChainType;
1088  typedef typename boost::mpl::at<BoolNodeChainType, boost::mpl::int_<1> >::type BoolInternalNodeType;
1089 
1091  size_t numLeafNodes = 0, numInternalNodes = 0;
1092 
1093  std::vector<const LeafNodeType*> nodes;
1094  std::vector<size_t> leafnodeCount;
1095 
1096  {
1097  // compute the prefix sum of the leafnode count in each internal node.
1098  std::vector<const InternalNodeType*> internalNodes;
1099  tree.getNodes(internalNodes);
1100 
1101  numInternalNodes = internalNodes.size();
1102 
1103  leafnodeCount.push_back(0);
1104  for (size_t n = 0; n < numInternalNodes; ++n) {
1105  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
1106  }
1107 
1108  numLeafNodes = leafnodeCount.back();
1109 
1110  // extract all leafnodes
1111  nodes.reserve(numLeafNodes);
1112 
1113  for (size_t n = 0; n < numInternalNodes; ++n) {
1114  internalNodes[n]->getNodes(nodes);
1115  }
1116  }
1117 
1118  // create mask leafnodes
1119  boost::scoped_array<BoolLeafNodeType*> maskNodes(new BoolLeafNodeType*[numLeafNodes]);
1120 
1121  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1122  MaskInteriorVoxels<LeafNodeType>(iso, &nodes[0], maskNodes.get()));
1123 
1124 
1125  // create mask grid
1126  typename BoolTreeType::Ptr maskTree(new BoolTreeType(false));
1127 
1128  PopulateTree<BoolTreeType> populate(*maskTree, maskNodes.get(), &leafnodeCount[0], false);
1129  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1130 
1131 
1132  // evaluate tile values
1133  std::vector<BoolInternalNodeType*> internalMaskNodes;
1134  maskTree->getNodes(internalMaskNodes);
1135 
1136  tbb::parallel_for(tbb::blocked_range<size_t>(0, internalMaskNodes.size()),
1137  MaskInteriorTiles<TreeType, BoolInternalNodeType>(iso, tree, &internalMaskNodes[0]));
1138 
1140 
1141  typename BoolTreeType::ValueAllIter it(*maskTree);
1142  it.setMaxDepth(BoolTreeType::ValueAllIter::LEAF_DEPTH - 2);
1143 
1144  for ( ; it; ++it) {
1145  if (acc.getValue(it.getCoord()) < iso) {
1146  it.setValue(true);
1147  it.setActiveState(true);
1148  }
1149  }
1150 
1151  return maskTree;
1152 } // computeInteriorMask()
1153 
1154 
1155 template<typename InputTreeType>
1157 {
1158  typedef typename InputTreeType::ValueType InputValueType;
1159  typedef typename InputTreeType::LeafNodeType InputLeafNodeType;
1160  typedef typename InputTreeType::template ValueConverter<bool>::Type BoolTreeType;
1161  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
1162 
1164  const InputTreeType& inputTree,
1165  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
1166  BoolTreeType& maskTree,
1167  InputValueType iso)
1168  : mInputAccessor(inputTree)
1169  , mInputNodes(!inputLeafNodes.empty() ? &inputLeafNodes.front() : NULL)
1170  , mMaskTree(false)
1171  , mMaskAccessor(maskTree)
1172  , mIsovalue(iso)
1173  {
1174  }
1175 
1177  : mInputAccessor(rhs.mInputAccessor.tree())
1178  , mInputNodes(rhs.mInputNodes)
1179  , mMaskTree(false)
1180  , mMaskAccessor(mMaskTree)
1181  , mIsovalue(rhs.mIsovalue)
1182  {
1183  }
1184 
1185  void operator()(const tbb::blocked_range<size_t>& range) {
1186 
1187  const InputValueType iso = mIsovalue;
1188  Coord ijk(0, 0, 0);
1189 
1190  BoolLeafNodeType* maskNodePt = NULL;
1191 
1192  for (size_t n = range.begin(); mInputNodes && (n != range.end()); ++n) {
1193 
1194  const InputLeafNodeType& node = *mInputNodes[n];
1195 
1196  if (!maskNodePt) maskNodePt = new BoolLeafNodeType(node.origin(), false);
1197  else maskNodePt->setOrigin(node.origin());
1198 
1199  bool collectedData = false;
1200 
1201  for (typename InputLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
1202 
1203  bool isUnder = *it < iso;
1204 
1205  ijk = it.getCoord();
1206 
1207  ++ijk[2];
1208  bool signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +z edge
1209  --ijk[2];
1210 
1211  if (!signChange) {
1212  --ijk[2];
1213  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -z edge
1214  ++ijk[2];
1215  }
1216 
1217  if (!signChange) {
1218  ++ijk[1];
1219  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +y edge
1220  --ijk[1];
1221  }
1222 
1223  if (!signChange) {
1224  --ijk[1];
1225  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -y edge
1226  ++ijk[1];
1227  }
1228 
1229  if (!signChange) {
1230  ++ijk[0];
1231  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +x edge
1232  --ijk[0];
1233  }
1234 
1235  if (!signChange) {
1236  --ijk[0];
1237  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -x edge
1238  ++ijk[0];
1239  }
1240 
1241  if (signChange) {
1242  collectedData = true;
1243  maskNodePt->setValueOn(it.pos(), true);
1244  }
1245  }
1246 
1247  if (collectedData) {
1248  mMaskAccessor.addLeaf(maskNodePt);
1249  maskNodePt = NULL;
1250  }
1251  }
1252 
1253  if (maskNodePt) delete maskNodePt;
1254  }
1255 
1257  mMaskAccessor.tree().merge(rhs.mMaskAccessor.tree());
1258  }
1259 
1260 private:
1262  InputLeafNodeType const * const * const mInputNodes;
1263 
1264  BoolTreeType mMaskTree;
1265  tree::ValueAccessor<BoolTreeType> mMaskAccessor;
1266 
1267  InputValueType mIsovalue;
1268 }; // MaskIsovalueCrossingVoxels
1269 
1270 
1272 
1273 
1274 template<typename NodeType>
1276 {
1277  typedef boost::shared_ptr<NodeMaskSegment> Ptr;
1278  typedef typename NodeType::NodeMaskType NodeMaskType;
1279 
1280  NodeMaskSegment() : connections(), mask(false), origin(0,0,0), visited(false) {}
1281 
1282  std::vector<NodeMaskSegment*> connections;
1283  NodeMaskType mask;
1284  Coord origin;
1285  bool visited;
1286 }; // struct NodeMaskSegment
1287 
1288 
1289 template<typename NodeType>
1290 inline void
1291 nodeMaskSegmentation(const NodeType& node,
1292  std::vector<typename NodeMaskSegment<NodeType>::Ptr>& segments)
1293 {
1294  typedef typename NodeType::NodeMaskType NodeMaskType;
1295  typedef NodeMaskSegment<NodeType> NodeMaskSegmentType;
1296  typedef typename NodeMaskSegmentType::Ptr NodeMaskSegmentTypePtr;
1297 
1298  NodeMaskType nodeMask(node.getValueMask());
1299  std::deque<Index> indexList;
1300 
1301  while (!nodeMask.isOff()) {
1302 
1303  NodeMaskSegmentTypePtr segment(new NodeMaskSegmentType());
1304  segment->origin = node.origin();
1305 
1306  NodeMaskType& mask = segment->mask;
1307 
1308  indexList.push_back(nodeMask.findFirstOn());
1309  nodeMask.setOff(indexList.back()); // mark as visited
1310  Coord ijk(0, 0, 0);
1311 
1312  while (!indexList.empty()) {
1313 
1314  const Index pos = indexList.back();
1315  indexList.pop_back();
1316 
1317  if (mask.isOn(pos)) continue;
1318  mask.setOn(pos);
1319 
1320  ijk = NodeType::offsetToLocalCoord(pos);
1321 
1322  Index npos = pos - 1;
1323  if (ijk[2] != 0 && nodeMask.isOn(npos)) {
1324  nodeMask.setOff(npos);
1325  indexList.push_back(npos);
1326  }
1327 
1328  npos = pos + 1;
1329  if (ijk[2] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1330  nodeMask.setOff(npos);
1331  indexList.push_back(npos);
1332  }
1333 
1334  npos = pos - NodeType::DIM;
1335  if (ijk[1] != 0 && nodeMask.isOn(npos)) {
1336  nodeMask.setOff(npos);
1337  indexList.push_back(npos);
1338  }
1339 
1340  npos = pos + NodeType::DIM;
1341  if (ijk[1] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1342  nodeMask.setOff(npos);
1343  indexList.push_back(npos);
1344  }
1345 
1346  npos = pos - NodeType::DIM * NodeType::DIM;
1347  if (ijk[0] != 0 && nodeMask.isOn(npos)) {
1348  nodeMask.setOff(npos);
1349  indexList.push_back(npos);
1350  }
1351 
1352  npos = pos + NodeType::DIM * NodeType::DIM;
1353  if (ijk[0] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1354  nodeMask.setOff(npos);
1355  indexList.push_back(npos);
1356  }
1357 
1358  }
1359 
1360  segments.push_back(segment);
1361  }
1362 }
1363 
1364 
1365 template<typename NodeType>
1367 {
1370  typedef typename std::vector<NodeMaskSegmentTypePtr> NodeMaskSegmentVector;
1371 
1372  SegmentNodeMask(std::vector<NodeType*>& nodes, NodeMaskSegmentVector* nodeMaskArray)
1373  : mNodes(!nodes.empty() ? &nodes.front() : NULL)
1374  , mNodeMaskArray(nodeMaskArray)
1375  {
1376  }
1377 
1378  void operator()(const tbb::blocked_range<size_t>& range) const {
1379  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1380  NodeType& node = *mNodes[n];
1381  nodeMaskSegmentation(node, mNodeMaskArray[n]);
1382 
1383  // hack origin data to store array offset
1384  Coord& origin = const_cast<Coord&>(node.origin());
1385  origin[0] = static_cast<int>(n);
1386  }
1387  }
1388 
1389  NodeType * const * const mNodes;
1390  NodeMaskSegmentVector * const mNodeMaskArray;
1391 }; // struct SegmentNodeMask
1392 
1393 
1394 template<typename TreeType, typename NodeType>
1396 {
1397  typedef typename NodeType::NodeMaskType NodeMaskType;
1400  typedef typename std::vector<NodeMaskSegmentTypePtr> NodeMaskSegmentVector;
1401 
1402  ConnectNodeMaskSegments(const TreeType& tree, NodeMaskSegmentVector* nodeMaskArray)
1403  : mTree(&tree)
1404  , mNodeMaskArray(nodeMaskArray)
1405  {
1406  }
1407 
1408  void operator()(const tbb::blocked_range<size_t>& range) const {
1409 
1411 
1412  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1413 
1414  NodeMaskSegmentVector& segments = mNodeMaskArray[n];
1415  if (segments.empty()) continue;
1416 
1417  std::vector<std::set<NodeMaskSegmentType*> > connections(segments.size());
1418 
1419  Coord ijk = segments[0]->origin;
1420 
1421  const NodeType* node = acc.template probeConstNode<NodeType>(ijk);
1422  if (!node) continue;
1423 
1424  // get neighbour nodes
1425 
1426  ijk[2] += NodeType::DIM;
1427  const NodeType* nodeZUp = acc.template probeConstNode<NodeType>(ijk);
1428  ijk[2] -= (NodeType::DIM + NodeType::DIM);
1429  const NodeType* nodeZDown = acc.template probeConstNode<NodeType>(ijk);
1430  ijk[2] += NodeType::DIM;
1431 
1432  ijk[1] += NodeType::DIM;
1433  const NodeType* nodeYUp = acc.template probeConstNode<NodeType>(ijk);
1434  ijk[1] -= (NodeType::DIM + NodeType::DIM);
1435  const NodeType* nodeYDown = acc.template probeConstNode<NodeType>(ijk);
1436  ijk[1] += NodeType::DIM;
1437 
1438  ijk[0] += NodeType::DIM;
1439  const NodeType* nodeXUp = acc.template probeConstNode<NodeType>(ijk);
1440  ijk[0] -= (NodeType::DIM + NodeType::DIM);
1441  const NodeType* nodeXDown = acc.template probeConstNode<NodeType>(ijk);
1442  ijk[0] += NodeType::DIM;
1443 
1444  const Index startPos = node->getValueMask().findFirstOn();
1445  for (Index pos = startPos; pos < NodeMaskType::SIZE; ++pos) {
1446 
1447  if (!node->isValueOn(pos)) continue;
1448 
1449  ijk = NodeType::offsetToLocalCoord(pos);
1450  Index npos = 0;
1451 
1452  if (ijk[2] == 0) {
1453  npos = pos + (NodeType::DIM - 1);
1454  if (nodeZDown && nodeZDown->isValueOn(npos)) {
1455  NodeMaskSegmentType* nsegment =
1456  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZDown)], npos);
1457  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1458  connections[idx].insert(nsegment);
1459  }
1460  } else if (ijk[2] == (NodeType::DIM - 1)) {
1461  npos = pos - (NodeType::DIM - 1);
1462  if (nodeZUp && nodeZUp->isValueOn(npos)) {
1463  NodeMaskSegmentType* nsegment =
1464  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZUp)], npos);
1465  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1466  connections[idx].insert(nsegment);
1467  }
1468  }
1469 
1470  if (ijk[1] == 0) {
1471  npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1472  if (nodeYDown && nodeYDown->isValueOn(npos)) {
1473  NodeMaskSegmentType* nsegment =
1474  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYDown)], npos);
1475  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1476  connections[idx].insert(nsegment);
1477  }
1478  } else if (ijk[1] == (NodeType::DIM - 1)) {
1479  npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1480  if (nodeYUp && nodeYUp->isValueOn(npos)) {
1481  NodeMaskSegmentType* nsegment =
1482  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYUp)], npos);
1483  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1484  connections[idx].insert(nsegment);
1485  }
1486  }
1487 
1488  if (ijk[0] == 0) {
1489  npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1490  if (nodeXDown && nodeXDown->isValueOn(npos)) {
1491  NodeMaskSegmentType* nsegment =
1492  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXDown)], npos);
1493  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1494  connections[idx].insert(nsegment);
1495  }
1496  } else if (ijk[0] == (NodeType::DIM - 1)) {
1497  npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1498  if (nodeXUp && nodeXUp->isValueOn(npos)) {
1499  NodeMaskSegmentType* nsegment =
1500  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXUp)], npos);
1501  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1502  connections[idx].insert(nsegment);
1503  }
1504  }
1505  }
1506 
1507  for (size_t i = 0, I = connections.size(); i < I; ++i) {
1508 
1509  typename std::set<NodeMaskSegmentType*>::iterator
1510  it = connections[i].begin(), end = connections[i].end();
1511 
1512  std::vector<NodeMaskSegmentType*>& segmentConnections = segments[i]->connections;
1513  segmentConnections.reserve(connections.size());
1514  for (; it != end; ++it) {
1515  segmentConnections.push_back(*it);
1516  }
1517  }
1518  } // end range loop
1519  }
1520 
1521 private:
1522 
1523  static inline size_t getNodeOffset(const NodeType& node) {
1524  return static_cast<size_t>(node.origin()[0]);
1525  }
1526 
1527  static inline NodeMaskSegmentType*
1528  findNodeMaskSegment(NodeMaskSegmentVector& segments, Index pos)
1529  {
1530  NodeMaskSegmentType* segment = NULL;
1531 
1532  for (size_t n = 0, N = segments.size(); n < N; ++n) {
1533  if (segments[n]->mask.isOn(pos)) {
1534  segment = segments[n].get();
1535  break;
1536  }
1537  }
1538 
1539  return segment;
1540  }
1541 
1542  static inline Index
1543  findNodeMaskSegmentIndex(NodeMaskSegmentVector& segments, Index pos)
1544  {
1545  for (Index n = 0, N = Index(segments.size()); n < N; ++n) {
1546  if (segments[n]->mask.isOn(pos)) return n;
1547  }
1548  return Index(-1);
1549  }
1550 
1551  TreeType const * const mTree;
1552  NodeMaskSegmentVector * const mNodeMaskArray;
1553 }; // struct ConnectNodeMaskSegments
1554 
1555 
1556 template<typename TreeType>
1558 {
1559  typedef typename TreeType::LeafNodeType LeafNodeType;
1560  typedef typename TreeType::Ptr TreeTypePtr;
1562 
1563  MaskSegmentGroup(const std::vector<NodeMaskSegmentType*>& segments)
1564  : mSegments(!segments.empty() ? &segments.front() : NULL)
1565  , mTree(new TreeType(false))
1566  {
1567  }
1568 
1569  MaskSegmentGroup(const MaskSegmentGroup& rhs, tbb::split)
1570  : mSegments(rhs.mSegments)
1571  , mTree(new TreeType(false))
1572  {
1573  }
1574 
1575  TreeTypePtr& mask() { return mTree; }
1576 
1577  void join(MaskSegmentGroup& rhs) { mTree->merge(*rhs.mTree); }
1578 
1579  void operator()(const tbb::blocked_range<size_t>& range) {
1580 
1581  tree::ValueAccessor<TreeType> acc(*mTree);
1582 
1583  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1584  NodeMaskSegmentType& segment = *mSegments[n];
1585  LeafNodeType* node = acc.touchLeaf(segment.origin);
1586  node->getValueMask() |= segment.mask;
1587  }
1588  }
1589 
1590 private:
1591  NodeMaskSegmentType * const * const mSegments;
1592  TreeTypePtr mTree;
1593 }; // struct MaskSegmentGroup
1594 
1595 
1597 
1598 
1599 template<typename TreeType>
1601 {
1602  typedef typename TreeType::ValueType ValueType;
1603  typedef typename TreeType::LeafNodeType LeafNodeType;
1604  typedef typename LeafNodeType::NodeMaskType NodeMaskType;
1605 
1606  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
1607  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
1608 
1610 
1611  ExpandLeafNodeRegion(const TreeType& distTree, BoolTreeType& maskTree, std::vector<BoolLeafNodeType*>& maskNodes)
1612  : mDistTree(&distTree)
1613  , mMaskTree(&maskTree)
1614  , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : NULL)
1615  , mNewMaskTree(false)
1616  {
1617  }
1618 
1620  : mDistTree(rhs.mDistTree)
1621  , mMaskTree(rhs.mMaskTree)
1622  , mMaskNodes(rhs.mMaskNodes)
1623  , mNewMaskTree(false)
1624  {
1625  }
1626 
1627  BoolTreeType& newMaskTree() { return mNewMaskTree; }
1628 
1629  void join(ExpandLeafNodeRegion& rhs) { mNewMaskTree.merge(rhs.mNewMaskTree); }
1630 
1631  void operator()(const tbb::blocked_range<size_t>& range) {
1632 
1633  typedef LeafNodeType NodeType;
1634 
1635  tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1636  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1637  tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
1638 
1639  NodeMaskType maskZUp, maskZDown, maskYUp, maskYDown, maskXUp, maskXDown;
1640 
1641  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1642 
1643  BoolLeafNodeType& maskNode = *mMaskNodes[n];
1644  if (maskNode.isEmpty()) continue;
1645 
1646  Coord ijk = maskNode.origin(), nijk;
1647 
1648  const LeafNodeType* distNode = distAcc.probeConstLeaf(ijk);
1649  if (!distNode) continue;
1650 
1651  const ValueType *dataZUp = NULL, *dataZDown = NULL,
1652  *dataYUp = NULL, *dataYDown = NULL,
1653  *dataXUp = NULL, *dataXDown = NULL;
1654 
1655  ijk[2] += NodeType::DIM;
1656  getData(ijk, distAcc, maskAcc, maskZUp, dataZUp);
1657  ijk[2] -= (NodeType::DIM + NodeType::DIM);
1658  getData(ijk, distAcc, maskAcc, maskZDown, dataZDown);
1659  ijk[2] += NodeType::DIM;
1660 
1661  ijk[1] += NodeType::DIM;
1662  getData(ijk, distAcc, maskAcc, maskYUp, dataYUp);
1663  ijk[1] -= (NodeType::DIM + NodeType::DIM);
1664  getData(ijk, distAcc, maskAcc, maskYDown, dataYDown);
1665  ijk[1] += NodeType::DIM;
1666 
1667  ijk[0] += NodeType::DIM;
1668  getData(ijk, distAcc, maskAcc, maskXUp, dataXUp);
1669  ijk[0] -= (NodeType::DIM + NodeType::DIM);
1670  getData(ijk, distAcc, maskAcc, maskXDown, dataXDown);
1671  ijk[0] += NodeType::DIM;
1672 
1673  for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
1674 
1675  const Index pos = it.pos();
1676  const ValueType val = std::abs(distNode->getValue(pos));
1677 
1678  ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
1679  nijk = ijk + maskNode.origin();
1680 
1681  if (dataZUp && ijk[2] == (BoolLeafNodeType::DIM - 1)) {
1682  const Index npos = pos - (NodeType::DIM - 1);
1683  if (maskZUp.isOn(npos) && std::abs(dataZUp[npos]) > val) {
1684  newMaskAcc.setValueOn(nijk.offsetBy(0, 0, 1));
1685  }
1686  } else if (dataZDown && ijk[2] == 0) {
1687  const Index npos = pos + (NodeType::DIM - 1);
1688  if (maskZDown.isOn(npos) && std::abs(dataZDown[npos]) > val) {
1689  newMaskAcc.setValueOn(nijk.offsetBy(0, 0, -1));
1690  }
1691  }
1692 
1693  if (dataYUp && ijk[1] == (BoolLeafNodeType::DIM - 1)) {
1694  const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1695  if (maskYUp.isOn(npos) && std::abs(dataYUp[npos]) > val) {
1696  newMaskAcc.setValueOn(nijk.offsetBy(0, 1, 0));
1697  }
1698  } else if (dataYDown && ijk[1] == 0) {
1699  const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1700  if (maskYDown.isOn(npos) && std::abs(dataYDown[npos]) > val) {
1701  newMaskAcc.setValueOn(nijk.offsetBy(0, -1, 0));
1702  }
1703  }
1704 
1705  if (dataXUp && ijk[0] == (BoolLeafNodeType::DIM - 1)) {
1706  const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1707  if (maskXUp.isOn(npos) && std::abs(dataXUp[npos]) > val) {
1708  newMaskAcc.setValueOn(nijk.offsetBy(1, 0, 0));
1709  }
1710  } else if (dataXDown && ijk[0] == 0) {
1711  const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1712  if (maskXDown.isOn(npos) && std::abs(dataXDown[npos]) > val) {
1713  newMaskAcc.setValueOn(nijk.offsetBy(-1, 0, 0));
1714  }
1715  }
1716 
1717  } // end value on loop
1718  } // end range loop
1719  }
1720 
1721 private:
1722 
1723  static inline void
1724  getData(const Coord& ijk, tree::ValueAccessor<const TreeType>& distAcc,
1725  tree::ValueAccessor<const BoolTreeType>& maskAcc, NodeMaskType& mask, const ValueType*& data)
1726  {
1727  const LeafNodeType* node = distAcc.probeConstLeaf(ijk);
1728  if (node) {
1729  data = node->buffer().data();
1730  mask = node->getValueMask();
1731  const BoolLeafNodeType* maskNodePt = maskAcc.probeConstLeaf(ijk);
1732  if (maskNodePt) mask -= maskNodePt->getValueMask();
1733  }
1734  }
1735 
1736  TreeType const * const mDistTree;
1737  BoolTreeType * const mMaskTree;
1738  BoolLeafNodeType ** const mMaskNodes;
1739 
1740  BoolTreeType mNewMaskTree;
1741 }; // struct ExpandLeafNodeRegion
1742 
1743 
1744 template<typename TreeType>
1746 {
1747  typedef typename TreeType::ValueType ValueType;
1748  typedef typename TreeType::LeafNodeType LeafNodeType;
1749  typedef typename LeafNodeType::NodeMaskType NodeMaskType;
1751 
1752  FillLeafNodeVoxels(const TreeType& tree, std::vector<BoolLeafNodeType*>& maskNodes)
1753  : mTree(&tree), mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : NULL)
1754  {
1755  }
1756 
1757  void operator()(const tbb::blocked_range<size_t>& range) const {
1758 
1759  tree::ValueAccessor<const TreeType> distAcc(*mTree);
1760 
1761  std::vector<Index> indexList;
1762  indexList.reserve(NodeMaskType::SIZE);
1763 
1764  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1765 
1766  BoolLeafNodeType& maskNode = *mMaskNodes[n];
1767 
1768  const LeafNodeType * distNode = distAcc.probeConstLeaf(maskNode.origin());
1769  if (!distNode) continue;
1770 
1771  NodeMaskType mask(distNode->getValueMask());
1772  NodeMaskType& narrowbandMask = maskNode.getValueMask();
1773 
1774  for (Index pos = narrowbandMask.findFirstOn(); pos < NodeMaskType::SIZE; ++pos) {
1775  if (narrowbandMask.isOn(pos)) indexList.push_back(pos);
1776  }
1777 
1778  mask -= narrowbandMask; // bitwise difference
1779  narrowbandMask.setOff();
1780 
1781  const ValueType* data = distNode->buffer().data();
1782  Coord ijk(0, 0, 0);
1783 
1784  while (!indexList.empty()) {
1785 
1786  const Index pos = indexList.back();
1787  indexList.pop_back();
1788 
1789  if (narrowbandMask.isOn(pos)) continue;
1790  narrowbandMask.setOn(pos);
1791 
1792  const ValueType dist = std::abs(data[pos]);
1793 
1794  ijk = LeafNodeType::offsetToLocalCoord(pos);
1795 
1796  Index npos = pos - 1;
1797  if (ijk[2] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1798  mask.setOff(npos);
1799  indexList.push_back(npos);
1800  }
1801 
1802  npos = pos + 1;
1803  if (ijk[2] != (LeafNodeType::DIM - 1) && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1804  mask.setOff(npos);
1805  indexList.push_back(npos);
1806  }
1807 
1808  npos = pos - LeafNodeType::DIM;
1809  if (ijk[1] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1810  mask.setOff(npos);
1811  indexList.push_back(npos);
1812  }
1813 
1814  npos = pos + LeafNodeType::DIM;
1815  if (ijk[1] != (LeafNodeType::DIM - 1) && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1816  mask.setOff(npos);
1817  indexList.push_back(npos);
1818  }
1819 
1820  npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1821  if (ijk[0] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1822  mask.setOff(npos);
1823  indexList.push_back(npos);
1824  }
1825 
1826  npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1827  if (ijk[0] != (LeafNodeType::DIM - 1) && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1828  mask.setOff(npos);
1829  indexList.push_back(npos);
1830  }
1831  } // end flood fill loop
1832  } // end range loop
1833  }
1834 
1835  TreeType const * const mTree;
1836  BoolLeafNodeType ** const mMaskNodes;
1837 }; // FillLeafNodeVoxels
1838 
1839 
1840 template<typename TreeType>
1842 {
1843  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
1844  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
1845  typedef typename BoolTreeType::Ptr BoolTreeTypePtr;
1846 
1847  ExpandNarrowbandMask(const TreeType& tree, std::vector<BoolTreeTypePtr>& segments)
1848  : mTree(&tree), mSegments(!segments.empty() ? &segments.front() : NULL)
1849  {
1850  }
1851 
1852  void operator()(const tbb::blocked_range<size_t>& range) const {
1853 
1854  const TreeType& distTree = *mTree;
1855  std::vector<BoolLeafNodeType*> nodes;
1856 
1857  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1858 
1859  BoolTreeType& narrowBandMask = *mSegments[n];
1860 
1861  BoolTreeType candidateMask(narrowBandMask, false, TopologyCopy());
1862 
1863  while (true) {
1864 
1865  nodes.clear();
1866  candidateMask.getNodes(nodes);
1867  if (nodes.empty()) break;
1868 
1869  const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
1870 
1871  tbb::parallel_for(nodeRange, FillLeafNodeVoxels<TreeType>(distTree, nodes));
1872 
1873  narrowBandMask.topologyUnion(candidateMask);
1874 
1875  ExpandLeafNodeRegion<TreeType> op(distTree, narrowBandMask, nodes);
1876  tbb::parallel_reduce(nodeRange, op);
1877 
1878  if (op.newMaskTree().empty()) break;
1879 
1880  candidateMask.clear();
1881  candidateMask.merge(op.newMaskTree());
1882  } // end expand loop
1883  } // end range loop
1884  }
1885 
1886  TreeType const * const mTree;
1887  BoolTreeTypePtr * const mSegments;
1888 }; // ExpandNarrowbandMask
1889 
1890 
1891 template<typename TreeType>
1893 {
1894  typedef typename TreeType::Ptr TreeTypePtr;
1895  typedef typename TreeType::ValueType ValueType;
1896  typedef typename TreeType::LeafNodeType LeafNodeType;
1897  typedef typename TreeType::RootNodeType RootNodeType;
1898  typedef typename RootNodeType::NodeChainType NodeChainType;
1899  typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
1900 
1901  FloodFillSign(const TreeType& tree, std::vector<TreeTypePtr>& segments)
1902  : mTree(&tree)
1903  , mSegments(!segments.empty() ? &segments.front() : NULL)
1904  , mMinValue(ValueType(0.0))
1905  {
1906  ValueType minSDFValue = std::numeric_limits<ValueType>::max();
1907 
1908  {
1909  std::vector<const InternalNodeType*> nodes;
1910  tree.getNodes(nodes);
1911 
1912  if (!nodes.empty()) {
1913  FindMinTileValue<InternalNodeType> minOp(&nodes[0]);
1914  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1915  minSDFValue = std::min(minSDFValue, minOp.minValue);
1916  }
1917  }
1918 
1919  if (minSDFValue > ValueType(0.0)) {
1920  std::vector<const LeafNodeType*> nodes;
1921  tree.getNodes(nodes);
1922  if (!nodes.empty()) {
1923  FindMinVoxelValue<LeafNodeType> minOp(&nodes[0]);
1924  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1925  minSDFValue = std::min(minSDFValue, minOp.minValue);
1926  }
1927  }
1928 
1929  mMinValue = minSDFValue;
1930  }
1931 
1932  void operator()(const tbb::blocked_range<size_t>& range) const {
1933  const ValueType interiorValue = -std::abs(mMinValue);
1934  const ValueType exteriorValue = std::abs(mTree->background());
1935  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1936  tools::signedFloodFillWithValues(*mSegments[n], exteriorValue, interiorValue);
1937  }
1938  }
1939 
1940 private:
1941 
1942  TreeType const * const mTree;
1943  TreeTypePtr * const mSegments;
1944  ValueType mMinValue;
1945 }; // FloodFillSign
1946 
1947 
1948 template<typename TreeType>
1950 {
1951  typedef typename TreeType::Ptr TreeTypePtr;
1952  typedef typename TreeType::ValueType ValueType;
1953  typedef typename TreeType::LeafNodeType LeafNodeType;
1954 
1955  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
1956  typedef typename BoolTreeType::Ptr BoolTreeTypePtr;
1957  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
1958 
1959  MaskedCopy(const TreeType& tree, std::vector<TreeTypePtr>& segments, std::vector<BoolTreeTypePtr>& masks)
1960  : mTree(&tree)
1961  , mSegments(!segments.empty() ? &segments.front() : NULL)
1962  , mMasks(!masks.empty() ? &masks.front() : NULL)
1963  {
1964  }
1965 
1966  void operator()(const tbb::blocked_range<size_t>& range) const {
1967 
1968  std::vector<const BoolLeafNodeType*> nodes;
1969 
1970  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1971 
1972  const BoolTreeType& mask = *mMasks[n];
1973 
1974  nodes.clear();
1975  mask.getNodes(nodes);
1976 
1977  Copy op(*mTree, nodes);
1978  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
1979  mSegments[n] = op.outputTree();
1980  }
1981  }
1982 
1983 private:
1984 
1985  struct Copy {
1986  Copy(const TreeType& inputTree, std::vector<const BoolLeafNodeType*>& maskNodes)
1987  : mInputTree(&inputTree)
1988  , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : NULL)
1989  , mOutputTreePtr(new TreeType(inputTree.background()))
1990  {
1991  }
1992 
1993  Copy(const Copy& rhs, tbb::split)
1994  : mInputTree(rhs.mInputTree)
1995  , mMaskNodes(rhs.mMaskNodes)
1996  , mOutputTreePtr(new TreeType(mInputTree->background()))
1997  {
1998  }
1999 
2000  TreeTypePtr& outputTree() { return mOutputTreePtr; }
2001 
2002  void join(Copy& rhs) { mOutputTreePtr->merge(*rhs.mOutputTreePtr); }
2003 
2004  void operator()(const tbb::blocked_range<size_t>& range) {
2005 
2006  tree::ValueAccessor<const TreeType> inputAcc(*mInputTree);
2007  tree::ValueAccessor<TreeType> outputAcc(*mOutputTreePtr);
2008 
2009  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2010 
2011  const BoolLeafNodeType& maskNode = *mMaskNodes[n];
2012  if (maskNode.isEmpty()) continue;
2013 
2014  const Coord& ijk = maskNode.origin();
2015 
2016  const LeafNodeType* inputNode = inputAcc.probeConstLeaf(ijk);
2017  if (inputNode) {
2018 
2019  LeafNodeType* outputNode = outputAcc.touchLeaf(ijk);
2020 
2021  for (typename BoolLeafNodeType::ValueOnCIter it = maskNode.cbeginValueOn(); it; ++it) {
2022  const Index idx = it.pos();
2023  outputNode->setValueOn(idx, inputNode->getValue(idx));
2024  }
2025  } else {
2026  const int valueDepth = inputAcc.getValueDepth(ijk);
2027  if (valueDepth >= 0) {
2028  outputAcc.addTile(TreeType::RootNodeType::LEVEL - valueDepth,
2029  ijk, inputAcc.getValue(ijk), true);
2030  }
2031  }
2032  }
2033  }
2034 
2035  private:
2036  TreeType const * const mInputTree;
2037  BoolLeafNodeType const * const * const mMaskNodes;
2038  TreeTypePtr mOutputTreePtr;
2039  }; // struct Copy
2040 
2041  TreeType const * const mTree;
2042  TreeTypePtr * const mSegments;
2043  BoolTreeTypePtr * const mMasks;
2044 }; // MaskedCopy
2045 
2046 
2048 
2049 
2050 template<typename VolumePtrType>
2052 {
2053  ComputeActiveVoxelCount(std::vector<VolumePtrType>& segments, size_t *countArray)
2054  : mSegments(!segments.empty() ? &segments.front() : NULL)
2055  , mCountArray(countArray)
2056  {
2057  }
2058 
2059  void operator()(const tbb::blocked_range<size_t>& range) const {
2060  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2061  mCountArray[n] = mSegments[n]->activeVoxelCount();
2062  }
2063  }
2064 
2065  VolumePtrType * const mSegments;
2066  size_t * const mCountArray;
2067 };
2068 
2069 
2071 {
2072  GreaterCount(const size_t *countArray) : mCountArray(countArray) {}
2073 
2074  inline bool operator() (const size_t& lhs, const size_t& rhs) const
2075  {
2076  return (mCountArray[lhs] > mCountArray[rhs]);
2077  }
2078 
2079  size_t const * const mCountArray;
2080 };
2081 
2083 
2084 
2085 template<typename TreeType>
2087 {
2088  typedef typename TreeType::Ptr TreeTypePtr;
2089  typedef typename TreeType::template ValueConverter<bool>::Type::Ptr BoolTreePtrType;
2090 
2091  static BoolTreePtrType constructMask(const TreeType&, BoolTreePtrType& maskTree) { return maskTree; }
2092  static TreeTypePtr construct(const TreeType&, TreeTypePtr& tree) { return tree; }
2093 };
2094 
2095 
2096 template<typename TreeType>
2097 struct GridOrTreeConstructor<Grid<TreeType> >
2098 {
2101  typedef typename TreeType::Ptr TreeTypePtr;
2102 
2103  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2104  typedef typename BoolTreeType::Ptr BoolTreePtrType;
2107 
2108  static BoolGridPtrType constructMask(const GridType& grid, BoolTreePtrType& maskTree) {
2109  BoolGridPtrType maskGrid(BoolGridType::create(maskTree));
2110  maskGrid->setTransform(grid.transform().copy());
2111  return maskGrid;
2112  }
2113 
2114  static GridTypePtr construct(const GridType& grid, TreeTypePtr& maskTree) {
2115  GridTypePtr maskGrid(GridType::create(maskTree));
2116  maskGrid->setTransform(grid.transform().copy());
2117  maskGrid->insertMeta(grid);
2118  return maskGrid;
2119  }
2120 };
2121 
2122 
2123 } // namespace level_set_util_internal
2124 
2125 
2127 
2128 
2129 template <class GridType>
2130 inline void
2131 sdfToFogVolume(GridType& grid, typename GridType::ValueType cutoffDistance)
2132 {
2133  typedef typename GridType::ValueType ValueType;
2134  typedef typename GridType::TreeType TreeType;
2135  typedef typename TreeType::LeafNodeType LeafNodeType;
2136  typedef typename TreeType::RootNodeType RootNodeType;
2137  typedef typename RootNodeType::NodeChainType NodeChainType;
2138  typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
2139 
2141 
2142  TreeType& tree = grid.tree();
2143 
2144  size_t numLeafNodes = 0, numInternalNodes = 0;
2145 
2146  std::vector<LeafNodeType*> nodes;
2147  std::vector<size_t> leafnodeCount;
2148 
2149  {
2150  // Compute the prefix sum of the leafnode count in each internal node.
2151  std::vector<InternalNodeType*> internalNodes;
2152  tree.getNodes(internalNodes);
2153 
2154  numInternalNodes = internalNodes.size();
2155 
2156  leafnodeCount.push_back(0);
2157  for (size_t n = 0; n < numInternalNodes; ++n) {
2158  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
2159  }
2160 
2161  numLeafNodes = leafnodeCount.back();
2162 
2163  // Steal all leafnodes (Removes them from the tree and transfers ownership.)
2164  nodes.reserve(numLeafNodes);
2165 
2166  for (size_t n = 0; n < numInternalNodes; ++n) {
2167  internalNodes[n]->stealNodes(nodes, tree.background(), false);
2168  }
2169 
2170  // Clamp cutoffDistance to min sdf value
2171  ValueType minSDFValue = std::numeric_limits<ValueType>::max();
2172 
2173  {
2175  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, internalNodes.size()), minOp);
2176  minSDFValue = std::min(minSDFValue, minOp.minValue);
2177  }
2178 
2179  if (minSDFValue > ValueType(0.0)) {
2181  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
2182  minSDFValue = std::min(minSDFValue, minOp.minValue);
2183  }
2184 
2185  cutoffDistance = -std::abs(cutoffDistance);
2186  cutoffDistance = minSDFValue > cutoffDistance ? minSDFValue : cutoffDistance;
2187  }
2188 
2189  // Transform voxel values and delete leafnodes that are uniformly zero after the transformation.
2190  // (Positive values are set to zero with inactive state and negative values are remapped
2191  // from zero to one with active state.)
2192  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
2194 
2195  // Populate a new tree with the remaining leafnodes
2196  typename TreeType::Ptr newTree(new TreeType(ValueType(0.0)));
2197 
2198  level_set_util_internal::PopulateTree<TreeType> populate(*newTree, &nodes[0], &leafnodeCount[0], 0);
2199  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
2200 
2201  // Transform tile values (Negative valued tiles are set to 1.0 with active state.)
2202  std::vector<InternalNodeType*> internalNodes;
2203  newTree->getNodes(internalNodes);
2204 
2205  tbb::parallel_for(tbb::blocked_range<size_t>(0, internalNodes.size()),
2207 
2208  {
2210 
2211  typename TreeType::ValueAllIter it(*newTree);
2212  it.setMaxDepth(TreeType::ValueAllIter::LEAF_DEPTH - 2);
2213 
2214  for ( ; it; ++it) {
2215  if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
2216  it.setValue(ValueType(1.0));
2217  it.setActiveState(true);
2218  }
2219  }
2220  }
2221 
2222  // Insert missing root level tiles. (The new tree is constructed from the remaining leafnodes
2223  // and will therefore not contain any root level tiles that may exist in the original tree.)
2224  {
2225  typename TreeType::ValueAllIter it(tree);
2226  it.setMaxDepth(TreeType::ValueAllIter::ROOT_DEPTH);
2227  for ( ; it; ++it) {
2228  if (it.getValue() < ValueType(0.0)) {
2229  newTree->addTile(TreeType::ValueAllIter::ROOT_LEVEL, it.getCoord(), ValueType(1.0), true);
2230  }
2231  }
2232  }
2233 
2234  grid.setTree(newTree);
2235  grid.setGridClass(GRID_FOG_VOLUME);
2236 }
2237 
2238 
2240 
2241 
2242 template <class GridOrTreeType>
2243 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2244 sdfInteriorMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2245 {
2246  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2247  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2248 
2249  typedef typename TreeType::template ValueConverter<bool>::Type::Ptr BoolTreePtrType;
2250  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(tree, isovalue);
2251 
2253 }
2254 
2255 
2256 template<typename GridOrTreeType>
2257 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2258 extractEnclosedRegion(const GridOrTreeType& volume,
2259  typename GridOrTreeType::ValueType isovalue,
2260  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type* fillMask)
2261 {
2262  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2263  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2264 
2265  typedef typename TreeType::template ValueConverter<char>::Type::Ptr CharTreePtrType;
2266  CharTreePtrType regionMask = level_set_util_internal::computeEnclosedRegionMask(tree, isovalue, fillMask);
2267 
2268  typedef typename TreeType::template ValueConverter<bool>::Type::Ptr BoolTreePtrType;
2269  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(*regionMask, 0);
2270 
2272 }
2273 
2274 
2276 
2277 
2278 template<typename GridOrTreeType>
2279 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2280 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2281 {
2282  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2283  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2284 
2285  std::vector<const typename TreeType::LeafNodeType*> nodes;
2286  tree.getNodes(nodes);
2287 
2288  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2289  typename BoolTreeType::Ptr mask(new BoolTreeType(false));
2290 
2291  level_set_util_internal::MaskIsovalueCrossingVoxels<TreeType> op(tree, nodes, *mask, isovalue);
2292  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2293 
2295 }
2296 
2297 
2299 
2300 
2301 template<typename GridOrTreeType>
2302 inline void
2303 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
2304  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks)
2305 {
2306  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2307  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2308  typedef typename BoolTreeType::Ptr BoolTreePtrType;
2309  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
2310 
2312  typedef typename NodeMaskSegmentType::Ptr NodeMaskSegmentPtrType;
2313  typedef typename std::vector<NodeMaskSegmentPtrType> NodeMaskSegmentPtrVector;
2314  typedef typename std::vector<NodeMaskSegmentType*> NodeMaskSegmentRawPtrVector;
2315 
2317 
2318  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2319 
2320  BoolTreeType topologyMask(tree, false, TopologyCopy());
2321 
2322  if (topologyMask.hasActiveTiles()) {
2323  topologyMask.voxelizeActiveTiles();
2324  }
2325 
2326  std::vector<BoolLeafNodeType*> leafnodes;
2327  topologyMask.getNodes(leafnodes);
2328 
2329  if (leafnodes.empty()) return;
2330 
2331  // 1. Split node masks into disjoint segments
2332  // Note: The LeafNode origin coord is modified to record the 'leafnodes' array offset.
2333 
2334  boost::scoped_array<NodeMaskSegmentPtrVector> nodeSegmentArray(new NodeMaskSegmentPtrVector[leafnodes.size()]);
2335 
2336  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2337  level_set_util_internal::SegmentNodeMask<BoolLeafNodeType>(leafnodes, nodeSegmentArray.get()));
2338 
2339 
2340  // 2. Compute segment connectivity
2341 
2342  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2344  topologyMask, nodeSegmentArray.get()));
2345 
2346  topologyMask.clear();
2347 
2348  size_t nodeSegmentCount = 0;
2349  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2350  nodeSegmentCount += nodeSegmentArray[n].size();
2351  }
2352 
2353  // 3. Group connected segments
2354 
2355  std::deque<NodeMaskSegmentRawPtrVector> nodeSegmentGroups;
2356 
2357  NodeMaskSegmentType* nextSegment = nodeSegmentArray[0][0].get();
2358  while (nextSegment) {
2359 
2360  nodeSegmentGroups.push_back(NodeMaskSegmentRawPtrVector());
2361 
2362  std::vector<NodeMaskSegmentType*>& segmentGroup = nodeSegmentGroups.back();
2363  segmentGroup.reserve(nodeSegmentCount);
2364 
2365  std::deque<NodeMaskSegmentType*> segmentQueue;
2366  segmentQueue.push_back(nextSegment);
2367  nextSegment = NULL;
2368 
2369  while (!segmentQueue.empty()) {
2370 
2371  NodeMaskSegmentType* segment = segmentQueue.back();
2372  segmentQueue.pop_back();
2373 
2374  if (segment->visited) continue;
2375  segment->visited = true;
2376 
2377  segmentGroup.push_back(segment);
2378 
2379  // queue connected segments
2380  std::vector<NodeMaskSegmentType*>& connections = segment->connections;
2381  for (size_t n = 0, N = connections.size(); n < N; ++n) {
2382  if (!connections[n]->visited) segmentQueue.push_back(connections[n]);
2383  }
2384  }
2385 
2386  // find first unvisited segment
2387  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2388  NodeMaskSegmentPtrVector& nodeSegments = nodeSegmentArray[n];
2389  for (size_t i = 0, I = nodeSegments.size(); i < I; ++i) {
2390  if (!nodeSegments[i]->visited) nextSegment = nodeSegments[i].get();
2391  }
2392  }
2393  }
2394 
2395  // 4. Mask segment groups
2396 
2397  if (nodeSegmentGroups.size() == 1) {
2398 
2399  BoolTreePtrType mask(new BoolTreeType(tree, false, TopologyCopy()));
2400 
2401  if (mask->hasActiveTiles()) {
2402  mask->voxelizeActiveTiles();
2403  }
2404 
2405  masks.push_back(
2407 
2408  } else if (nodeSegmentGroups.size() > 1) {
2409 
2410  for (size_t n = 0, N = nodeSegmentGroups.size(); n < N; ++n) {
2411 
2412  NodeMaskSegmentRawPtrVector& segmentGroup = nodeSegmentGroups[n];
2413 
2415  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, segmentGroup.size()), op);
2416 
2417  masks.push_back(
2419  }
2420  }
2421 
2422  // 5. Sort segments in descending order based on the active voxel count.
2423 
2424  if (masks.size() > 1) {
2425  const size_t segmentCount = masks.size();
2426 
2427  boost::scoped_array<size_t> segmentOrderArray(new size_t[segmentCount]);
2428  boost::scoped_array<size_t> voxelCountArray(new size_t[segmentCount]);
2429 
2430  for (size_t n = 0; n < segmentCount; ++n) {
2431  segmentOrderArray[n] = n;
2432  }
2433 
2434  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
2436 
2437  size_t *begin = segmentOrderArray.get();
2438  tbb::parallel_sort(begin, begin + masks.size(), level_set_util_internal::GreaterCount(voxelCountArray.get()));
2439 
2440  std::vector<BoolTreePtrType> orderedMasks;
2441  orderedMasks.reserve(masks.size());
2442 
2443  for (size_t n = 0; n < segmentCount; ++n) {
2444  orderedMasks.push_back(masks[segmentOrderArray[n]]);
2445  }
2446 
2447  masks.swap(orderedMasks);
2448  }
2449 
2450 } // extractActiveVoxelSegmentMasks()
2451 
2452 
2453 template<typename GridOrTreeType>
2454 inline void
2455 segmentActiveVoxels(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments)
2456 {
2457  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2458  typedef typename TreeType::Ptr TreePtrType;
2459  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2460  typedef typename BoolTreeType::Ptr BoolTreePtrType;
2461 
2462  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2463 
2464  // 1. Segment active topology mask
2465  std::vector<BoolTreePtrType> maskSegmentArray;
2466  extractActiveVoxelSegmentMasks(inputTree, maskSegmentArray);
2467 
2468  const size_t numSegments = maskSegmentArray.size();
2469 
2470  if (numSegments < 2) {
2471  // single segment early-out
2472  TreePtrType segment(new TreeType(inputTree));
2473  segments.push_back(
2475  return;
2476  }
2477 
2478  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2479 
2480  // 2. Export segments
2481  std::vector<TreePtrType> outputSegmentArray(numSegments);
2482 
2483  tbb::parallel_for(segmentRange,
2484  level_set_util_internal::MaskedCopy<TreeType>(inputTree, outputSegmentArray, maskSegmentArray));
2485 
2486  for (size_t n = 0, N = numSegments; n < N; ++n) {
2487  segments.push_back(
2489  }
2490 }
2491 
2492 
2493 template<typename GridOrTreeType>
2494 inline void
2495 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments)
2496 {
2497  typedef typename TreeAdapter<GridOrTreeType>::TreeType TreeType;
2498  typedef typename TreeType::Ptr TreePtrType;
2499  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2500  typedef typename BoolTreeType::Ptr BoolTreePtrType;
2501 
2502  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2503 
2504  // 1. Mask zero crossing voxels
2505  BoolTreePtrType mask = extractIsosurfaceMask(inputTree, lsutilGridZero<GridOrTreeType>());
2506 
2507  // 2. Segment the zero crossing mask
2508  std::vector<BoolTreePtrType> maskSegmentArray;
2509  extractActiveVoxelSegmentMasks(*mask, maskSegmentArray);
2510 
2511  const size_t numSegments = maskSegmentArray.size();
2512 
2513  if (numSegments < 2) {
2514  // single segment early-out
2515  TreePtrType segment(new TreeType(inputTree));
2516  segments.push_back(
2518  return;
2519  }
2520 
2521  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2522 
2523 
2524  // 3. Expand zero crossing mask to capture sdf narrow band
2525  tbb::parallel_for(segmentRange,
2526  level_set_util_internal::ExpandNarrowbandMask<TreeType>(inputTree, maskSegmentArray));
2527 
2528  // 4. Export sdf segments
2529  std::vector<TreePtrType> outputSegmentArray(numSegments);
2530 
2531  tbb::parallel_for(segmentRange,
2532  level_set_util_internal::MaskedCopy<TreeType>(inputTree, outputSegmentArray, maskSegmentArray));
2533 
2534  tbb::parallel_for(segmentRange,
2535  level_set_util_internal::FloodFillSign<TreeType>(inputTree, outputSegmentArray));
2536 
2537 
2538  for (size_t n = 0, N = numSegments; n < N; ++n) {
2539  segments.push_back(
2541  }
2542 }
2543 
2544 
2545 } // namespace tools
2546 } // namespace OPENVDB_VERSION_NAME
2547 } // namespace openvdb
2548 
2549 #endif // OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
2550 
2551 // Copyright (c) 2012-2016 DreamWorks Animation LLC
2552 // All rights reserved. This software is distributed under the
2553 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
2554 
tree::LeafNode< bool, LeafNodeType::LOG2DIM > BoolLeafNodeType
Definition: LevelSetUtil.h:1750
SDFVoxelsToFogVolume(LeafNodeType **nodes, ValueType cutoffDistance)
Definition: LevelSetUtil.h:465
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
Definition: ValueAccessor.h:374
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1843
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:294
NodeMaskSegment< NodeType > NodeMaskSegmentType
Definition: LevelSetUtil.h:1398
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:526
CharLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:368
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:424
Definition: Types.h:214
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:242
const Coord & origin() const
Return the grid index coordinates of this node&#39;s local origin.
Definition: LeafNode.h:475
BoolLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:1836
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:470
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1957
tree::LeafNode< bool, LeafNodeType::LOG2DIM > BoolLeafNodeType
Definition: LevelSetUtil.h:205
FillLeafNodeVoxels(const TreeType &tree, std::vector< BoolLeafNodeType * > &maskNodes)
Definition: LevelSetUtil.h:1752
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1844
Definition: Types.h:442
SegmentNodeMask(std::vector< NodeType * > &nodes, NodeMaskSegmentVector *nodeMaskArray)
Definition: LevelSetUtil.h:1372
const NodeMaskType & getValueMask() const
Definition: LeafNode.h:1109
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1161
void join(FindMinTileValue &rhs)
Definition: LevelSetUtil.h:452
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:367
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:278
NodeMaskSegment< LeafNodeType > NodeMaskSegmentType
Definition: LevelSetUtil.h:1561
void setValueOn(const Coord &xyz)
Mark the voxel at the given coordinates as active but don&#39;t change its value.
Definition: LeafNode.h:714
LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:375
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1757
MaskInteriorTiles(ValueType isovalue, const TreeType &tree, InternalNodeType **maskNodes)
Definition: LevelSetUtil.h:253
void segmentSDF(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint SDF surfaces into distinct grids or trees.
Definition: LevelSetUtil.h:2495
RootNodeType::NodeChainType NodeChainType
Definition: LevelSetUtil.h:1898
PopulateTree(PopulateTree &rhs, tbb::split)
Definition: LevelSetUtil.h:290
static BoolTreePtrType constructMask(const TreeType &, BoolTreePtrType &maskTree)
Definition: LevelSetUtil.h:2091
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:501
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1378
tree::LeafNode< char, LeafNodeType::LOG2DIM > CharLeafNodeType
Definition: LevelSetUtil.h:330
LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:463
NodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1278
ComputeActiveVoxelCount(std::vector< VolumePtrType > &segments, size_t *countArray)
Definition: LevelSetUtil.h:2053
InternalNodeType::ValueType ValueType
Definition: LevelSetUtil.h:429
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1852
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1748
NodeMaskSegmentType::Ptr NodeMaskSegmentTypePtr
Definition: LevelSetUtil.h:1399
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:527
int getValueDepth(const Coord &xyz) const
Definition: ValueAccessor.h:275
TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1951
LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:395
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
TreeType::template ValueConverter< bool >::Type::Ptr BoolTreePtrType
Definition: LevelSetUtil.h:2089
SDFTilesToFogVolume(const TreeType &tree, InternalNodeType **nodes)
Definition: LevelSetUtil.h:498
std::vector< NodeMaskSegment * > connections
Definition: LevelSetUtil.h:1282
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1631
BoolTreeTypePtr *const mSegments
Definition: LevelSetUtil.h:1887
MaskSegmentGroup(const std::vector< NodeMaskSegmentType * > &segments)
Definition: LevelSetUtil.h:1563
TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1894
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1955
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
Ptr copy() const
Definition: Transform.h:77
ValueType const mIsovalue
Definition: LevelSetUtil.h:369
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:292
math::Transform & transform()
Return a reference to this grid&#39;s transform, which might be shared with other grids.
Definition: Grid.h:335
void join(ExpandLeafNodeRegion &rhs)
Definition: LevelSetUtil.h:1629
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1559
MaskedCopy(const TreeType &tree, std::vector< TreeTypePtr > &segments, std::vector< BoolTreeTypePtr > &masks)
Definition: LevelSetUtil.h:1959
void sdfToFogVolume(GridType &grid, typename GridType::ValueType cutoffDistance=lsutilGridMax< GridType >())
Threaded method to convert a sparse level set/SDF into a sparse fog volume.
Definition: LevelSetUtil.h:2131
BoolTreeType::Ptr BoolTreeTypePtr
Definition: LevelSetUtil.h:1845
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
Index32 Index
Definition: Types.h:58
GreaterCount(const size_t *countArray)
Definition: LevelSetUtil.h:2072
ValueType const mWeight
Definition: LevelSetUtil.h:491
MaskSegmentGroup(const MaskSegmentGroup &rhs, tbb::split)
Definition: LevelSetUtil.h:1569
NodeMaskSegment< NodeType > NodeMaskSegmentType
Definition: LevelSetUtil.h:1368
void extractActiveVoxelSegmentMasks(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::template ValueConverter< bool >::Type::Ptr > &masks)
Return a mask for each connected component of the given grid&#39;s active voxels.
Definition: LevelSetUtil.h:2303
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:338
TreeType::template ValueConverter< char >::Type::Ptr computeEnclosedRegionMask(const TreeType &tree, typename TreeType::ValueType isovalue, const typename TreeType::template ValueConverter< bool >::Type *fillMask)
Constructs a memory light char tree that represents the exterior region with +1 and the interior regi...
Definition: LevelSetUtil.h:956
#define OPENVDB_VERSION_NAME
Definition: version.h:43
void join(PopulateTree &rhs)
Definition: LevelSetUtil.h:315
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:2103
void join(MaskIsovalueCrossingVoxels &rhs)
Definition: LevelSetUtil.h:1256
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1952
ExpandLeafNodeRegion(const TreeType &distTree, BoolTreeType &maskTree, std::vector< BoolLeafNodeType * > &maskNodes)
Definition: LevelSetUtil.h:1611
PopulateTree(TreeType &tree, LeafNodeType **leafnodes, const size_t *nodexIndexMap, ValueType background)
Definition: LevelSetUtil.h:281
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1895
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
InputTreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1160
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:65
ConnectNodeMaskSegments(const TreeType &tree, NodeMaskSegmentVector *nodeMaskArray)
Definition: LevelSetUtil.h:1402
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1953
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1606
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:280
NodeMaskSegmentVector *const mNodeMaskArray
Definition: LevelSetUtil.h:1390
FloodFillSign(const TreeType &tree, std::vector< TreeTypePtr > &segments)
Definition: LevelSetUtil.h:1901
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:429
This adapter allows code that is templated on a Tree type to accept either a Tree type or a Grid type...
Definition: Grid.h:898
InternalNodeType **const mMaskNodes
Definition: LevelSetUtil.h:270
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
TreeType::RootNodeType RootNodeType
Definition: LevelSetUtil.h:1897
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1408
InputTreeType::LeafNodeType InputLeafNodeType
Definition: LevelSetUtil.h:1159
boost::shared_ptr< NodeMaskSegment > Ptr
Definition: LevelSetUtil.h:1277
Definition: Exceptions.h:39
FindMinVoxelValue(LeafNodeType const *const *const leafnodes)
Definition: LevelSetUtil.h:397
void addTile(Index level, const Coord &xyz, const ValueType &value, bool state)
Add a tile at the specified tree level that contains voxel (x, y, z), possibly deleting existing node...
Definition: ValueAccessor.h:382
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:422
void join(MaskSegmentGroup &rhs)
Definition: LevelSetUtil.h:1577
ValueType const mIsovalue
Definition: LevelSetUtil.h:244
FillMaskBoundary(const TreeType &tree, ValueType isovalue, const BoolTreeType &fillMask, const BoolLeafNodeType **fillNodes, BoolLeafNodeType **newNodes)
Definition: LevelSetUtil.h:530
InputTreeType::ValueType InputValueType
Definition: LevelSetUtil.h:1158
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:443
void join(FindMinVoxelValue &rhs)
Definition: LevelSetUtil.h:418
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1603
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:528
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2244
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1966
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition: MeshToVolume.h:2994
FlipRegionSign(LeafNodeType **nodes)
Definition: LevelSetUtil.h:377
VolumePtrType *const mSegments
Definition: LevelSetUtil.h:2065
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:256
LeafNodeT * touchLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, create one, but preserve the values and active states of all voxels.
Definition: ValueAccessor.h:393
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1607
std::vector< NodeMaskSegmentTypePtr > NodeMaskSegmentVector
Definition: LevelSetUtil.h:1400
_TreeType TreeType
Definition: Grid.h:900
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:279
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractEnclosedRegion(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >(), const typename TreeAdapter< GridOrTreeType >::TreeType::template ValueConverter< bool >::Type *fillMask=NULL)
Extracts the interior regions of a signed distance field and topologically enclosed (watertight) regi...
Definition: LevelSetUtil.h:2258
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1185
Index64 onVoxelCount() const
Return the number of voxels marked On.
Definition: LeafNode.h:440
Negative active values are set 0, everything else is set to 1.
Definition: LevelSetUtil.h:327
static BoolGridPtrType constructMask(const GridType &grid, BoolTreePtrType &maskTree)
Definition: LevelSetUtil.h:2108
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:2059
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:409
static GridTypePtr construct(const GridType &grid, TreeTypePtr &maskTree)
Definition: LevelSetUtil.h:2114
static TreeType & tree(TreeType &t)
Definition: Grid.h:915
void nodeMaskSegmentation(const NodeType &node, std::vector< typename NodeMaskSegment< NodeType >::Ptr > &segments)
Definition: LevelSetUtil.h:1291
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:525
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1602
LeafNodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1604
NodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1397
ExpandNarrowbandMask(const TreeType &tree, std::vector< BoolTreeTypePtr > &segments)
Definition: LevelSetUtil.h:1847
TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:2088
void setOff(Index32 n)
Set the nth bit off.
Definition: NodeMasks.h:454
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:298
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1747
TreeType const *const mTree
Definition: LevelSetUtil.h:269
size_t const *const mCountArray
Definition: LevelSetUtil.h:2079
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:379
ValueType const mIsovalue
Definition: LevelSetUtil.h:271
TreeType::template ValueConverter< bool >::Type::Ptr computeInteriorMask(const TreeType &tree, typename TreeType::ValueType iso)
Definition: LevelSetUtil.h:1077
MaskInteriorVoxels(ValueType isovalue, const LeafNodeType **nodes, BoolLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:207
TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1560
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:213
LabelBoundaryVoxels(ValueType isovalue, const LeafNodeType **nodes, CharLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:332
BoolTreeType & newMaskTree()
Definition: LevelSetUtil.h:1627
TreeType const *const mTree
Definition: LevelSetUtil.h:517
ExpandLeafNodeRegion(const ExpandLeafNodeRegion &rhs, tbb::split)
Definition: LevelSetUtil.h:1619
TreeType & tree() const
Return a reference to the tree associated with this accessor.
Definition: ValueAccessor.h:147
static TreeTypePtr construct(const TreeType &, TreeTypePtr &tree)
Definition: LevelSetUtil.h:2092
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
TreeTypePtr & mask()
Definition: LevelSetUtil.h:1575
LeafNodeType **const mNodes
Definition: LevelSetUtil.h:388
InternalNodeType **const mNodes
Definition: LevelSetUtil.h:518
BoolTreeType::Ptr BoolTreeTypePtr
Definition: LevelSetUtil.h:1956
LeafNodeType **const mNodes
Definition: LevelSetUtil.h:490
void segmentActiveVoxels(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint active topology components into distinct grids or trees.
Definition: LevelSetUtil.h:2455
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1932
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1579
boost::mpl::at< NodeChainType, boost::mpl::int_< 1 > >::type InternalNodeType
Definition: LevelSetUtil.h:1899
TreeType const *const mTree
Definition: LevelSetUtil.h:1835
FindMinVoxelValue(FindMinVoxelValue &rhs, tbb::split)
Definition: LevelSetUtil.h:403
NodeMaskSegmentType::Ptr NodeMaskSegmentTypePtr
Definition: LevelSetUtil.h:1369
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:536
InternalNodeType const *const *const mNodes
Definition: LevelSetUtil.h:456
BoolLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:243
void setOrigin(const Coord &origin)
Set the grid index coordinates of this node&#39;s local origin.
Definition: LeafNode.h:472
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractIsosurfaceMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue)
Return a mask of the voxels that intersect the implicit surface with the given isovalue.
Definition: LevelSetUtil.h:2280
TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1896
LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:204
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
MaskIsovalueCrossingVoxels(const InputTreeType &inputTree, const std::vector< const InputLeafNodeType * > &inputLeafNodes, BoolTreeType &maskTree, InputValueType iso)
Definition: LevelSetUtil.h:1163
TreeType const *const mTree
Definition: LevelSetUtil.h:1886
FindMinTileValue(FindMinTileValue &rhs, tbb::split)
Definition: LevelSetUtil.h:437
NodeType *const *const mNodes
Definition: LevelSetUtil.h:1389
LeafNodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1749
MaskIsovalueCrossingVoxels(MaskIsovalueCrossingVoxels &rhs, tbb::split)
Definition: LevelSetUtil.h:1176
LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:329
TreeType::ValueType ValueType
Definition: LevelSetUtil.h:251
std::vector< NodeMaskSegmentTypePtr > NodeMaskSegmentVector
Definition: LevelSetUtil.h:1370
NodeMaskType mask
Definition: LevelSetUtil.h:1283
FindMinTileValue(InternalNodeType const *const *const nodes)
Definition: LevelSetUtil.h:431
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:128