OpenVDB  1.1.0
VolumeToMesh.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 
31 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/tree/ValueAccessor.h>
35 #include <openvdb/tools/ValueTransformer.h>
36 #include <openvdb/util/Util.h>
37 #include <openvdb/math/Operators.h>
38 
39 #include <openvdb/tools/Morphology.h>
40 #include <boost/scoped_array.hpp>
41 #include <boost/intrusive/slist.hpp>
42 
43 #include <tbb/mutex.h>
44 #include <tbb/tick_count.h>
45 #include <tbb/blocked_range.h>
46 
47 #include <vector>
48 
49 #include <boost/scoped_ptr.hpp>
50 
51 namespace openvdb {
53 namespace OPENVDB_VERSION_NAME {
54 namespace tools {
55 
56 
58 
59 
61 enum {
63 };
64 
65 
68 {
69 public:
71  : mNumQuads(0)
72  , mNumTriangles(0)
73  , mQuads(NULL)
74  , mTriangles(NULL)
75  , mQuadFlags(NULL)
76  , mTriangleFlags(NULL)
77  {
78  }
79 
80  void resetQuads(size_t size)
81  {
82  mNumQuads = size;
83  mQuads.reset(new openvdb::Vec4I[mNumQuads]);
84  mQuadFlags.reset(new char[mNumQuads]);
85  }
86 
87  void clearQuads()
88  {
89  mNumQuads = 0;
90  mQuads.reset(NULL);
91  mQuadFlags.reset(NULL);
92  }
93 
94  void resetTriangles(size_t size)
95  {
96  mNumTriangles = size;
97  mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
98  mTriangleFlags.reset(new char[mNumTriangles]);
99  }
100 
101  void clearTriangles()
102  {
103  mNumTriangles = 0;
104  mTriangles.reset(NULL);
105  mTriangleFlags.reset(NULL);
106  }
107 
108  const size_t& numQuads() const { return mNumQuads; }
109  const size_t& numTriangles() const { return mNumTriangles; }
110 
111  // polygon accessor methods
112  openvdb::Vec4I& quad(size_t n) { return mQuads[n]; }
113  const openvdb::Vec4I& quad(size_t n) const { return mQuads[n]; }
114 
115  openvdb::Vec3I& triangle(size_t n) { return mTriangles[n]; }
116  const openvdb::Vec3I& triangle(size_t n) const { return mTriangles[n]; }
117 
118  // polygon flags accessor methods
119  char& quadFlags(size_t n) { return mQuadFlags[n]; }
120  const char& quadFlags(size_t n) const { return mQuadFlags[n]; }
121 
122  char& triangleFlags(size_t n) { return mTriangleFlags[n]; }
123  const char& triangleFlags(size_t n) const { return mTriangleFlags[n]; }
124 
125 private:
126  size_t mNumQuads, mNumTriangles;
127  boost::scoped_array<openvdb::Vec4I> mQuads;
128  boost::scoped_array<openvdb::Vec3I> mTriangles;
129  boost::scoped_array<char> mQuadFlags, mTriangleFlags;
130 };
131 
132 
135 typedef boost::scoped_array<openvdb::Vec3s> PointList;
136 typedef boost::scoped_array<PolygonPool> PolygonPoolList;
138 
139 
141 
142 
145 {
146 public:
147 
150  VolumeToMesh(double isovalue = 0, double adaptivity = 0);
151 
152  PointList& pointList();
153  const size_t& pointListSize() const;
154 
155  PolygonPoolList& polygonPoolList();
156  const PolygonPoolList& polygonPoolList() const;
157  const size_t& polygonPoolListSize() const;
158 
161  template<typename GridT>
162  void operator()(const GridT&);
163 
164 
186  void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0);
187 
188  // Deprecated in versions 1.1
189  template<typename GridT>
190  OPENVDB_DEPRECATED void accumulateAuxiliaryData(const GridT&) { return; }
191 
192 private:
193 
194  PointList mPoints;
195  PolygonPoolList mPolygons;
196 
197  size_t mPointListSize, mPolygonPoolListSize;
198  double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
199 
200  GridBase::ConstPtr mRefGrid;
201  TreeBase::Ptr mRefEdgeTree, mRefTopologyMaskTree;
202  bool mRefDataCached;
203 };
204 
205 
207 
208 
209 // Internal utility methods
210 
211 
212 namespace internal {
213 
214 
215 // Bit-flags
216 enum { INSIDE = 0x1, XEDGE = 0x2, YEDGE = 0x4, ZEDGE = 0x8 };
217 
218 const bool sAmbiguous[256] =
219  {0,0,0,0,0,1,0,0,
220  0,0,1,0,0,0,0,0,
221  0,0,1,0,1,1,1,0,
222  1,0,1,0,1,0,1,0,
223  0,1,0,0,1,1,0,0,
224  1,1,1,0,1,1,0,0,
225  0,0,0,0,1,1,0,0,
226  1,0,1,0,1,1,1,0,
227  0,1,1,1,0,1,0,0,
228  1,1,1,1,0,0,0,0,
229  1,1,1,1,1,1,1,1,
230  1,1,1,1,1,1,1,1,
231  0,1,0,0,0,1,0,0,
232  1,1,1,1,0,1,0,0,
233  0,0,0,0,0,1,0,0,
234  1,1,1,1,1,1,1,0,
235  0,1,1,1,1,1,1,1,
236  0,0,1,0,0,0,0,0,
237  0,0,1,0,1,1,1,1,
238  0,0,1,0,0,0,1,0,
239  1,1,1,1,1,1,1,1,
240  1,1,1,1,1,1,1,1,
241  0,0,0,0,1,1,1,1,
242  0,0,1,0,1,1,1,0,
243  0,1,1,1,0,1,0,1,
244  0,0,1,1,0,0,0,0,
245  0,0,1,1,0,1,1,1,
246  0,0,1,1,0,0,1,0,
247  0,1,0,1,0,1,0,1,
248  0,1,1,1,0,1,0,0,
249  0,0,0,0,0,1,0,0,
250  0,0,1,0,0,0,0,0};
251 
252 
253 template<class AccessorT>
254 inline bool isAmbiguous(const AccessorT& accessor, const Coord& ijk,
255  typename AccessorT::ValueType isovalue, int dim)
256 {
257  unsigned signs = 0;
258  Coord coord = ijk; // i, j, k
259  if (accessor.getValue(coord) < isovalue) signs |= 1u;
260  coord[0] += dim; // i+dim, j, k
261  if (accessor.getValue(coord) < isovalue) signs |= 2u;
262  coord[2] += dim; // i+dim, j, k+dim
263  if (accessor.getValue(coord) < isovalue) signs |= 4u;
264  coord[0] = ijk[0]; // i, j, k+dim
265  if (accessor.getValue(coord) < isovalue) signs |= 8u;
266  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
267  if (accessor.getValue(coord) < isovalue) signs |= 16u;
268  coord[0] += dim; // i+dim, j+dim, k
269  if (accessor.getValue(coord) < isovalue) signs |= 32u;
270  coord[2] += dim; // i+dim, j+dim, k+dim
271  if (accessor.getValue(coord) < isovalue) signs |= 64u;
272  coord[0] = ijk[0]; // i, j+dim, k+dim
273  if (accessor.getValue(coord) < isovalue) signs |= 128u;
274  return sAmbiguous[signs];
275 }
276 
277 
278 template<class AccessorT>
279 inline bool
280 isNonManifold(const AccessorT& accessor, const Coord& ijk,
281  typename AccessorT::ValueType isovalue, const int dim)
282 {
283  int hDim = dim >> 1;
284  bool m, p[8]; // Corner signs
285 
286  Coord coord = ijk; // i, j, k
287  p[0] = accessor.getValue(coord) < isovalue;
288  coord[0] += dim; // i+dim, j, k
289  p[1] = accessor.getValue(coord) < isovalue;
290  coord[2] += dim; // i+dim, j, k+dim
291  p[2] = accessor.getValue(coord) < isovalue;
292  coord[0] = ijk[0]; // i, j, k+dim
293  p[3] = accessor.getValue(coord) < isovalue;
294  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
295  p[4] = accessor.getValue(coord) < isovalue;
296  coord[0] += dim; // i+dim, j+dim, k
297  p[5] = accessor.getValue(coord) < isovalue;
298  coord[2] += dim; // i+dim, j+dim, k+dim
299  p[6] = accessor.getValue(coord) < isovalue;
300  coord[0] = ijk[0]; // i, j+dim, k+dim
301  p[7] = accessor.getValue(coord) < isovalue;
302 
303  // Check if the corner sign configuration is ambiguous
304  unsigned signs = 0;
305  if (p[0]) signs |= 1u;
306  if (p[1]) signs |= 2u;
307  if (p[2]) signs |= 4u;
308  if (p[3]) signs |= 8u;
309  if (p[4]) signs |= 16u;
310  if (p[5]) signs |= 32u;
311  if (p[6]) signs |= 64u;
312  if (p[7]) signs |= 128u;
313  if (sAmbiguous[signs]) return true;
314 
315  // Manifold check
316 
317  // Evaluate edges
318  int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
319  int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
320  int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
321 
322  // edge 1
323  coord.reset(ip, j, k);
324  m = accessor.getValue(coord) < isovalue;
325  if (p[0] != m && p[1] != m) return true;
326 
327  // edge 2
328  coord.reset(ipp, j, kp);
329  m = accessor.getValue(coord) < isovalue;
330  if (p[1] != m && p[2] != m) return true;
331 
332  // edge 3
333  coord.reset(ip, j, kpp);
334  m = accessor.getValue(coord) < isovalue;
335  if (p[2] != m && p[3] != m) return true;
336 
337  // edge 4
338  coord.reset(i, j, kp);
339  m = accessor.getValue(coord) < isovalue;
340  if (p[0] != m && p[3] != m) return true;
341 
342  // edge 5
343  coord.reset(ip, jpp, k);
344  m = accessor.getValue(coord) < isovalue;
345  if (p[4] != m && p[5] != m) return true;
346 
347  // edge 6
348  coord.reset(ipp, jpp, kp);
349  m = accessor.getValue(coord) < isovalue;
350  if (p[5] != m && p[6] != m) return true;
351 
352  // edge 7
353  coord.reset(ip, jpp, kpp);
354  m = accessor.getValue(coord) < isovalue;
355  if (p[6] != m && p[7] != m) return true;
356 
357  // edge 8
358  coord.reset(i, jpp, kp);
359  m = accessor.getValue(coord) < isovalue;
360  if (p[7] != m && p[4] != m) return true;
361 
362  // edge 9
363  coord.reset(i, jp, k);
364  m = accessor.getValue(coord) < isovalue;
365  if (p[0] != m && p[4] != m) return true;
366 
367  // edge 10
368  coord.reset(ipp, jp, k);
369  m = accessor.getValue(coord) < isovalue;
370  if (p[1] != m && p[5] != m) return true;
371 
372  // edge 11
373  coord.reset(ipp, jp, kpp);
374  m = accessor.getValue(coord) < isovalue;
375  if (p[2] != m && p[6] != m) return true;
376 
377 
378  // edge 12
379  coord.reset(i, jp, kpp);
380  m = accessor.getValue(coord) < isovalue;
381  if (p[3] != m && p[7] != m) return true;
382 
383 
384  // Evaluate faces
385 
386  // face 1
387  coord.reset(ip, jp, k);
388  m = accessor.getValue(coord) < isovalue;
389  if (p[0] != m && p[1] != m && p[4] != m && p[5] != m) return true;
390 
391  // face 2
392  coord.reset(ipp, jp, kp);
393  m = accessor.getValue(coord) < isovalue;
394  if (p[1] != m && p[2] != m && p[5] != m && p[6] != m) return true;
395 
396  // face 3
397  coord.reset(ip, jp, kpp);
398  m = accessor.getValue(coord) < isovalue;
399  if (p[2] != m && p[3] != m && p[6] != m && p[7] != m) return true;
400 
401  // face 4
402  coord.reset(i, jp, kp);
403  m = accessor.getValue(coord) < isovalue;
404  if (p[0] != m && p[3] != m && p[4] != m && p[7] != m) return true;
405 
406  // face 5
407  coord.reset(ip, j, kp);
408  m = accessor.getValue(coord) < isovalue;
409  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m) return true;
410 
411  // face 6
412  coord.reset(ip, jpp, kp);
413  m = accessor.getValue(coord) < isovalue;
414  if (p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
415 
416  // test cube center
417  coord.reset(ip, jp, kp);
418  m = accessor.getValue(coord) < isovalue;
419  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
420  p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
421 
422  return false;
423 }
424 
425 
427 
428 
429 template <class LeafType>
430 inline void
431 mergeVoxels(LeafType& leaf, const Coord& start, int dim, int regionId)
432 {
433  Coord ijk, end = start;
434  end[0] += dim;
435  end[1] += dim;
436  end[2] += dim;
437 
438  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
439  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
440  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
441  if(leaf.isValueOn(ijk)) leaf.setValue(ijk, regionId);
442  }
443  }
444  }
445 }
446 
447 
448 // Note that we must use ValueType::value_type or else Visual C++ gets confused
449 // thinking that it is a constructor.
450 template <class LeafType>
451 inline bool
452 isMergable(LeafType& leaf, const Coord& start, int dim,
453  typename LeafType::ValueType::value_type adaptivity)
454 {
455  if (adaptivity < 1e-6) return false;
456 
457  typedef typename LeafType::ValueType VecT;
458  Coord ijk, end = start;
459  end[0] += dim;
460  end[1] += dim;
461  end[2] += dim;
462 
463  std::vector<VecT> norms;
464  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
465  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
466  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
467 
468  if(!leaf.isValueOn(ijk)) continue;
469  norms.push_back(leaf.getValue(ijk));
470  }
471  }
472  }
473 
474  size_t N = norms.size();
475  for (size_t ni = 0; ni < N; ++ni) {
476  VecT n_i = norms[ni];
477  for (size_t nj = 0; nj < N; ++nj) {
478  VecT n_j = norms[nj];
479  if ((1.0 - n_i.dot(n_j)) > adaptivity) return false;
480  }
481  }
482  return true;
483 }
484 
485 
487 
488 
489 template <class TreeT>
491 {
492 public:
493  typedef std::vector<const typename TreeT::LeafNodeType *> ListT;
494 
495  LeafCPtrList(const TreeT& tree)
496  {
497  mLeafNodes.reserve(tree.leafCount());
498  typename TreeT::LeafCIter iter = tree.cbeginLeaf();
499  for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
500  }
501 
502  size_t size() const { return mLeafNodes.size(); }
503 
504  const typename TreeT::LeafNodeType* operator[](size_t n) const
505  { return mLeafNodes[n]; }
506 
507  tbb::blocked_range<size_t> getRange() const
508  { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
509 
510  const ListT& getList() const { return mLeafNodes; }
511 
512 private:
513  ListT mLeafNodes;
514 };
515 
516 
517 template <class TreeT>
519 {
520 public:
521  typedef std::vector<typename TreeT::LeafNodeType *> ListT;
522 
523  LeafPtrList(TreeT& tree)
524  {
525  mLeafNodes.reserve(tree.leafCount());
526  typename TreeT::LeafIter iter = tree.beginLeaf();
527  for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
528  }
529 
530  size_t size() const { return mLeafNodes.size(); }
531 
532  typename TreeT::LeafNodeType* operator[](size_t n) const
533  { return mLeafNodes[n]; }
534 
535  tbb::blocked_range<size_t> getRange() const
536  { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
537 
538  const ListT& getList() const { return mLeafNodes; }
539 
540 private:
541  ListT mLeafNodes;
542 };
543 
544 
546 
547 
548 template<typename DistTreeT>
550 {
551  typedef typename DistTreeT::ValueType DistValueT;
552  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
553  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
554 
556  : mDistTree(NULL)
557  , mEdgeTree(NULL)
558  , mTopologyMaskTree(NULL)
559  , mSeamMaskTree(typename BoolTreeT::Ptr())
560  , mInternalAdaptivity(DistValueT(0.0))
561  {
562  }
563 
564  bool isValid() const
565  {
566  return mDistTree && mEdgeTree && mTopologyMaskTree && mSeamMaskTree;
567  }
568 
569  const DistTreeT* mDistTree;
572  typename BoolTreeT::Ptr mSeamMaskTree;
574 };
575 
576 
578 
579 
580 template <class DistTreeT>
581 class Count
582 {
583 public:
584  Count(const LeafPtrList<DistTreeT>&, std::vector<size_t>&);
585  inline Count(const Count<DistTreeT>&);
586 
587  void runParallel();
588  void runSerial();
589 
590  inline void operator()(const tbb::blocked_range<size_t>&) const;
591 private:
592  const LeafPtrList<DistTreeT>& mLeafNodes;
593  std::vector<size_t>& mLeafRegionCount;
594 };
595 
596 
597 template <class DistTreeT>
599  const LeafPtrList<DistTreeT>& leafs,
600  std::vector<size_t>& leafRegionCount)
601  : mLeafNodes(leafs)
602  , mLeafRegionCount(leafRegionCount)
603 {
604 }
605 
606 
607 template <class DistTreeT>
608 inline
610  : mLeafNodes(rhs.mLeafNodes)
611  , mLeafRegionCount(rhs.mLeafRegionCount)
612 {
613 }
614 
615 
616 template <class DistTreeT>
617 void
619 {
620  tbb::parallel_for(mLeafNodes.getRange(), *this);
621 }
622 
623 
624 template <class DistTreeT>
625 void
627 {
628  (*this)(mLeafNodes.getRange());
629 }
630 
631 
632 template <class DistTreeT>
633 inline void
634 Count<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
635 {
636  for (size_t n = range.begin(); n != range.end(); ++n) {
637  mLeafRegionCount[n] = size_t(mLeafNodes[n]->onVoxelCount());
638  }
639 }
640 
641 
643 
644 
645 template <class DistTreeT>
646 class Merge
647 {
648 public:
649  typedef typename DistTreeT::ValueType DistValueT;
650  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
651  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
652 
653  Merge(
654  const DistTreeT& distTree,
655  LeafPtrList<IntTreeT>& auxLeafs,
656  std::vector<size_t>& leafRegionCount,
657  const DistValueT iso,
658  const DistValueT adaptivity);
659 
660  inline Merge(const Merge<DistTreeT>&);
661 
663 
664  void runParallel();
665  void runSerial();
666 
667  void operator()(const tbb::blocked_range<size_t>&) const;
668 
669 private:
670  const DistTreeT& mDistTree;
671  LeafPtrList<IntTreeT>& mAuxLeafs;
672  std::vector<size_t>& mLeafRegionCount;
673  const DistValueT mIsovalue, mAdaptivity;
674  const ReferenceData<DistTreeT>* mRefData;
675 };
676 
677 
678 template <class DistTreeT>
680  const DistTreeT& distTree,
681  LeafPtrList<IntTreeT>& auxLeafs,
682  std::vector<size_t>& leafRegionCount,
683  const DistValueT iso,
684  const DistValueT adaptivity)
685  : mDistTree(distTree)
686  , mAuxLeafs(auxLeafs)
687  , mLeafRegionCount(leafRegionCount)
688  , mIsovalue(iso)
689  , mAdaptivity(adaptivity)
690  , mRefData(NULL)
691 {
692 }
693 
694 
695 template <class DistTreeT>
696 inline
698  : mDistTree(rhs.mDistTree)
699  , mAuxLeafs(rhs.mAuxLeafs)
700  , mLeafRegionCount(rhs.mLeafRegionCount)
701  , mIsovalue(rhs.mIsovalue)
702  , mAdaptivity(rhs.mAdaptivity)
703  , mRefData(rhs.mRefData)
704 {
705 }
706 
707 
708 template <class DistTreeT>
709 void
711 {
712  tbb::parallel_for(mAuxLeafs.getRange(), *this);
713 }
714 
715 
716 template <class DistTreeT>
717 void
719 {
720  (*this)(mAuxLeafs.getRange());
721 }
722 
723 template <class DistTreeT>
724 void
726 {
727  mRefData = &refData;
728 }
729 
730 template <class DistTreeT>
731 void
732 Merge<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
733 {
734  typedef math::Vec3<DistValueT> Vec3T;
735  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
736  typedef typename IntTreeT::LeafNodeType IntLeafT;
737  typedef typename BoolLeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
738 
739  typedef typename IntLeafT::ValueOnIter IntIterT;
740  typedef typename BoolLeafT::ValueOnCIter BoolCIterT;
741 
742  typedef typename tree::ValueAccessor<BoolTreeT> BoolTreeAccessorT;
743  typedef typename tree::ValueAccessor<const BoolTreeT> BoolTreeCAccessorT;
744 
745  boost::scoped_ptr<BoolTreeAccessorT> seamMaskAcc;
746  boost::scoped_ptr<BoolTreeCAccessorT> topologyMaskAcc;
747  if (mRefData && mRefData->isValid()) {
748  seamMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
749  topologyMaskAcc.reset(new BoolTreeCAccessorT(*mRefData->mTopologyMaskTree));
750  }
751  const bool hasRefData = seamMaskAcc && topologyMaskAcc;
752 
753  const int LeafDim = BoolLeafT::DIM;
754  tree::ValueAccessor<const DistTreeT> distAcc(mDistTree);
755 
756  // Allocate reusable leaf buffers
757  BoolLeafT mask;
758  Vec3LeafT gradientBuffer;
759  Coord ijk, coord, end;
760 
761  for (size_t n = range.begin(); n != range.end(); ++n) {
762 
763  DistValueT adaptivity = mAdaptivity;
764  IntLeafT& auxLeaf = *mAuxLeafs[n];
765 
766  const Coord& origin = auxLeaf.getOrigin();
767  end[0] = origin[0] + LeafDim;
768  end[1] = origin[1] + LeafDim;
769  end[2] = origin[2] + LeafDim;
770 
771  mask.setValuesOff();
772 
773  // Mask off seam line adjacent voxels
774  if (hasRefData) {
775  const BoolLeafT* seamMask = seamMaskAcc->probeConstLeaf(origin);
776  if (seamMask != NULL) {
777  for (BoolCIterT it = seamMask->cbeginValueOn(); it; ++it) {
778  ijk = it.getCoord();
779  coord[0] = ijk[0] - (ijk[0] % 2);
780  coord[1] = ijk[1] - (ijk[1] % 2);
781  coord[2] = ijk[2] - (ijk[2] % 2);
782  mask.setActiveState(coord, true);
783  }
784  }
785  if (topologyMaskAcc->probeConstLeaf(origin) == NULL) {
786  adaptivity = mRefData->mInternalAdaptivity;
787  }
788  }
789 
790  // Mask off ambiguous voxels
791  for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
792  ijk = it.getCoord();
793  coord[0] = ijk[0] - (ijk[0] % 2);
794  coord[1] = ijk[1] - (ijk[1] % 2);
795  coord[2] = ijk[2] - (ijk[2] % 2);
796  if(mask.isValueOn(coord)) continue;
797  mask.setActiveState(coord, isAmbiguous(distAcc, ijk, mIsovalue, 1));
798  }
799 
800  int dim = 2;
801  // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
802  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
803  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
804  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
805  if (isNonManifold(distAcc, ijk, mIsovalue, dim)) {
806  mask.setActiveState(ijk, true);
807  }
808  }
809  }
810  }
811 
812  // Compute the gradient for the remaining voxels
813  gradientBuffer.setValuesOff();
814 
815  for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
816 
817  ijk = it.getCoord();
818  coord[0] = ijk[0] - (ijk[0] % dim);
819  coord[1] = ijk[1] - (ijk[1] % dim);
820  coord[2] = ijk[2] - (ijk[2] % dim);
821  if(mask.isValueOn(coord)) continue;
822 
823  Vec3T norm(math::ISGradient<math::CD_2ND>::result(distAcc, ijk));
824  // Normalize (Vec3's normalize uses isApproxEqual, which uses abs and does more work)
825  DistValueT length = norm.length();
826  if (length > DistValueT(1.0e-7)) {
827  norm *= DistValueT(1.0) / length;
828  }
829  gradientBuffer.setValue(ijk, norm);
830  }
831 
832  int regionId = 1, next_dim = dim << 1;
833 
834  // Process the first adaptivity level.
835  for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
836  coord[0] = ijk[0] - (ijk[0] % next_dim);
837  for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
838  coord[1] = ijk[1] - (ijk[1] % next_dim);
839  for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
840  coord[2] = ijk[2] - (ijk[2] % next_dim);
841  if(mask.isValueOn(ijk) || !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
842  mask.setActiveState(coord, true);
843  continue;
844  }
845  mergeVoxels(auxLeaf, ijk, dim, regionId++);
846  }
847  }
848  }
849 
850 
851  // Process remaining adaptivity levels
852  for (dim = 4; dim < LeafDim; dim = dim << 1) {
853  next_dim = dim << 1;
854  coord[0] = ijk[0] - (ijk[0] % next_dim);
855  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
856  coord[1] = ijk[1] - (ijk[1] % next_dim);
857  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
858  coord[2] = ijk[2] - (ijk[2] % next_dim);
859  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
860 
861  if (mask.isValueOn(ijk) || isNonManifold(distAcc, ijk, mIsovalue, dim) ||
862  !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
863  mask.setActiveState(coord, true);
864  continue;
865  }
866  mergeVoxels(auxLeaf, ijk, dim, regionId++);
867  }
868  }
869  }
870  }
871 
872 
873  if (!(mask.isValueOn(origin) || isNonManifold(distAcc, origin, mIsovalue, LeafDim))
874  && isMergable(gradientBuffer, origin, LeafDim, adaptivity)) {
875  mergeVoxels(auxLeaf, origin, LeafDim, regionId++);
876  }
877 
878 
879  // Count unique regions
880  size_t numVoxels = 0;
881  IntLeafT tmpLeaf(auxLeaf);
882  for (IntIterT it = tmpLeaf.beginValueOn(); it; ++it) {
883  if(it.getValue() == 0) {
884  it.setValueOff();
885  ++numVoxels;
886  }
887  }
888 
889  while (tmpLeaf.onVoxelCount() > 0) {
890  ++numVoxels;
891  IntIterT it = tmpLeaf.beginValueOn();
892  regionId = it.getValue();
893  for (; it; ++it) {
894  if (it.getValue() == regionId) it.setValueOff();
895  }
896  }
897 
898  mLeafRegionCount[n] = numVoxels;
899  }
900 }
901 
902 
904 
905 
906 template <class DistTreeT>
907 class PointGen
908 {
909 public:
910  typedef typename DistTreeT::ValueType DistValueT;
912  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
913 
914  PointGen(
915  const DistTreeT& distTree,
916  const LeafPtrList<IntTreeT>& auxLeafs,
917  std::vector<size_t>& leafIndices,
918  const openvdb::math::Transform& xform,
919  PointList& points,
920  double iso = 0.0);
921 
923 
925 
926  void runParallel();
927  void runSerial();
928 
929  void operator()(const tbb::blocked_range<size_t>&) const;
930 
931 private:
932  const DistTreeT& mDistTree;
933  const LeafPtrList<IntTreeT>& mAuxLeafs;
934  const std::vector<size_t>& mLeafIndices;
935  const openvdb::math::Transform& mTransform;
936  const PointList& mPoints;
937  const double mIsovalue;
938  const ReferenceData<DistTreeT>* mRefData;
939 
940  double root(double v0, double v1) const { return (mIsovalue - v0) / (v1 - v0); }
941  int calcAvgPoint(DistTreeAccessorT&, const Coord&, openvdb::Vec3d&) const;
942 };
943 
944 
945 template <class DistTreeT>
947  const DistTreeT& distTree,
948  const LeafPtrList<IntTreeT>& auxLeafs,
949  std::vector<size_t>& leafIndices,
950  const openvdb::math::Transform& xform,
951  PointList& points,
952  double iso)
953  : mDistTree(distTree)
954  , mAuxLeafs(auxLeafs)
955  , mLeafIndices(leafIndices)
956  , mTransform(xform)
957  , mPoints(points)
958  , mIsovalue(iso)
959  , mRefData(NULL)
960 {
961 }
962 
963 
964 template <class DistTreeT>
966  : mDistTree(rhs.mDistTree)
967  , mAuxLeafs(rhs.mAuxLeafs)
968  , mLeafIndices(rhs.mLeafIndices)
969  , mTransform(rhs.mTransform)
970  , mPoints(rhs.mPoints)
971  , mIsovalue(rhs.mIsovalue)
972  , mRefData(rhs.mRefData)
973 {
974 }
975 
976 
977 template <class DistTreeT>
978 void
980  const ReferenceData<DistTreeT>& refData)
981 {
982  mRefData = &refData;
983 }
984 
985 template <class DistTreeT>
986 void
988 {
989  tbb::parallel_for(mAuxLeafs.getRange(), *this);
990 }
991 
992 
993 template <class DistTreeT>
994 void
996 {
997  (*this)(mAuxLeafs.getRange());
998 }
999 
1000 
1001 template <class DistTreeT>
1002 void
1003 PointGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
1004 {
1005  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1006 
1007  typedef tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
1008 
1009  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1010  typedef typename IntTreeT::LeafNodeType IntLeafT;
1011 
1012 
1013  boost::scoped_ptr<DistTreeAccessorT> refDistAcc;
1014  boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc;
1015 
1016  if (mRefData && mRefData->isValid()) {
1017  refDistAcc.reset(new DistTreeAccessorT(*mRefData->mDistTree));
1018  refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mTopologyMaskTree));
1019  }
1020 
1021 
1022  const bool hasRefData = refDistAcc && refMaskAcc;
1023  typename IntTreeT::LeafNodeType::ValueOnIter auxIter;
1024  DistTreeAccessorT distAcc(mDistTree);
1025 
1026  Coord ijk;
1027  openvdb::Vec3d avg, tmp;
1028 
1029  for (size_t n = range.begin(); n != range.end(); ++n) {
1030 
1031  size_t idx = mLeafIndices[n];
1032  IntLeafT& auxLeaf = *mAuxLeafs[n];
1033 
1034  const BoolLeafT* maskLeaf = NULL;
1035  if (hasRefData) {
1036  maskLeaf = refMaskAcc->probeConstLeaf(auxLeaf.getOrigin());
1037  }
1038 
1039  for (auxIter = auxLeaf.beginValueOn(); auxIter; ++auxIter) {
1040 
1041  if(auxIter.getValue() == 0) {
1042 
1043  auxIter.setValue(idx);
1044  auxIter.setValueOff();
1045  ijk = auxIter.getCoord();
1046 
1047  if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1048  calcAvgPoint(*refDistAcc.get(), ijk, avg);
1049  } else {
1050  calcAvgPoint(distAcc, ijk, avg);
1051  }
1052 
1053  // offset by cell-origin
1054  avg[0] += double(ijk[0]);
1055  avg[1] += double(ijk[1]);
1056  avg[2] += double(ijk[2]);
1057 
1058  avg = mTransform.indexToWorld(avg);
1059 
1060  openvdb::Vec3s& ptn = mPoints[idx];
1061  ptn[0] = float(avg[0]);
1062  ptn[1] = float(avg[1]);
1063  ptn[2] = float(avg[2]);
1064 
1065  ++idx;
1066  }
1067  }
1068 
1069  while(auxLeaf.onVoxelCount() > 0) {
1070 
1071  avg[0] = 0;
1072  avg[1] = 0;
1073  avg[2] = 0;
1074 
1075  auxIter = auxLeaf.beginValueOn();
1076  int regionId = auxIter.getValue(), points = 0;
1077 
1078  for (; auxIter; ++auxIter) {
1079  if(auxIter.getValue() == regionId) {
1080 
1081  auxIter.setValue(idx);
1082  auxIter.setValueOff();
1083  ijk = auxIter.getCoord();
1084 
1085  if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1086  calcAvgPoint(*refDistAcc.get(), ijk, tmp);
1087  } else {
1088  calcAvgPoint(distAcc, ijk, tmp);
1089  }
1090 
1091  // offset by cell-origin
1092  tmp[0] += double(ijk[0]);
1093  tmp[1] += double(ijk[1]);
1094  tmp[2] += double(ijk[2]);
1095 
1096  avg += mTransform.indexToWorld(tmp);
1097  ++points;
1098  }
1099  }
1100 
1101  if (points > 1) {
1102  double w = 1.0 / double(points);
1103  avg[0] *= w;
1104  avg[1] *= w;
1105  avg[2] *= w;
1106  }
1107 
1108  openvdb::Vec3s& ptn = mPoints[idx];
1109  ptn[0] = float(avg[0]);
1110  ptn[1] = float(avg[1]);
1111  ptn[2] = float(avg[2]);
1112  ++idx;
1113  }
1114  }
1115 }
1116 
1117 template <class DistTreeT>
1118 int
1119 PointGen<DistTreeT>::calcAvgPoint(DistTreeAccessorT& acc,
1120  const Coord& ijk, openvdb::Vec3d& avg) const
1121 {
1122  double values[8];
1123  bool signMask[8];
1124  Coord coord;
1125 
1126  // Sample corner values
1127  coord = ijk;
1128  values[0] = double(acc.getValue(coord)); // i, j, k
1129 
1130  coord[0] += 1;
1131  values[1] = double(acc.getValue(coord)); // i+1, j, k
1132 
1133  coord[2] += 1;
1134  values[2] = double(acc.getValue(coord)); // i+i, j, k+1
1135 
1136  coord[0] = ijk[0];
1137  values[3] = double(acc.getValue(coord)); // i, j, k+1
1138 
1139  coord[1] += 1; coord[2] = ijk[2];
1140  values[4] = double(acc.getValue(coord)); // i, j+1, k
1141 
1142  coord[0] += 1;
1143  values[5] = double(acc.getValue(coord)); // i+1, j+1, k
1144 
1145  coord[2] += 1;
1146  values[6] = double(acc.getValue(coord)); // i+1, j+1, k+1
1147 
1148  coord[0] = ijk[0];
1149  values[7] = double(acc.getValue(coord)); // i, j+1, k+1
1150 
1151  // init sign mask
1152  for (int n = 0; n < 8; ++n) signMask[n] = (values[n] < mIsovalue);
1153 
1154  int samples = 0, edgeFlags = 0;
1155  avg[0] = 0.0;
1156  avg[1] = 0.0;
1157  avg[2] = 0.0;
1158 
1159  if (signMask[0] != signMask[1]) { // Edged: 0 - 1
1160  avg[0] += root(values[0], values[1]);
1161  ++samples;
1162  edgeFlags |= XEDGE;
1163  }
1164 
1165  if (signMask[1] != signMask[2]) { // Edged: 1 - 2
1166  avg[0] += 1.0;
1167  avg[2] += root(values[1], values[2]);
1168  ++samples;
1169  edgeFlags |= ZEDGE;
1170  }
1171 
1172  if (signMask[3] != signMask[2]) { // Edged: 3 - 2
1173  avg[0] += root(values[3], values[2]);
1174  avg[2] += 1.0;
1175  ++samples;
1176  edgeFlags |= XEDGE;
1177  }
1178 
1179  if (signMask[0] != signMask[3]) { // Edged: 0 - 3
1180  avg[2] += root(values[0], values[3]);
1181  ++samples;
1182  edgeFlags |= ZEDGE;
1183  }
1184 
1185  if (signMask[4] != signMask[5]) { // Edged: 4 - 5
1186  avg[0] += root(values[4], values[5]);
1187  avg[1] += 1.0;
1188  ++samples;
1189  edgeFlags |= XEDGE;
1190  }
1191 
1192  if (signMask[5] != signMask[6]) { // Edged: 5 - 6
1193  avg[0] += 1.0;
1194  avg[1] += 1.0;
1195  avg[2] += root(values[5], values[6]);
1196  ++samples;
1197  edgeFlags |= ZEDGE;
1198  }
1199 
1200  if (signMask[7] != signMask[6]) { // Edged: 7 - 6
1201  avg[0] += root(values[7], values[6]);
1202  avg[1] += 1.0;
1203  avg[2] += 1.0;
1204  ++samples;
1205  edgeFlags |= XEDGE;
1206  }
1207 
1208  if (signMask[4] != signMask[7]) { // Edged: 4 - 7
1209  avg[1] += 1.0;
1210  avg[2] += root(values[4], values[7]);
1211  ++samples;
1212  edgeFlags |= ZEDGE;
1213  }
1214 
1215  if (signMask[0] != signMask[4]) { // Edged: 0 - 4
1216  avg[1] += root(values[0], values[4]);
1217  ++samples;
1218  edgeFlags |= YEDGE;
1219  }
1220 
1221  if (signMask[1] != signMask[5]) { // Edged: 1 - 5
1222  avg[0] += 1.0;
1223  avg[1] += root(values[1], values[5]);
1224  ++samples;
1225  edgeFlags |= YEDGE;
1226  }
1227 
1228  if (signMask[2] != signMask[6]) { // Edged: 2 - 6
1229  avg[0] += 1.0;
1230  avg[1] += root(values[2], values[6]);
1231  avg[2] += 1.0;
1232  ++samples;
1233  edgeFlags |= YEDGE;
1234  }
1235 
1236  if (signMask[3] != signMask[7]) { // Edged: 3 - 7
1237  avg[1] += root(values[3], values[7]);
1238  avg[2] += 1.0;
1239  ++samples;
1240  edgeFlags |= YEDGE;
1241  }
1242 
1243  if (samples > 1) {
1244  double w = 1.0 / double(samples);
1245  avg[0] *= w;
1246  avg[1] *= w;
1247  avg[2] *= w;
1248  }
1249 
1250  return edgeFlags;
1251 }
1252 
1253 
1255 
1256 
1258 {
1259  QuadMeshOp(): mIdx(0), mPolygonPool(NULL) {}
1260 
1261  void init(const size_t upperBound, PolygonPool& quadPool)
1262  {
1263  mPolygonPool = &quadPool;
1264  mPolygonPool->resetQuads(upperBound);
1265  mIdx = 0;
1266  }
1267 
1268  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
1269  {
1270  if (!reverse) {
1271  mPolygonPool->quad(mIdx) = verts;
1272  } else {
1273  Vec4I& quad = mPolygonPool->quad(mIdx);
1274  quad[0] = verts[3];
1275  quad[1] = verts[2];
1276  quad[2] = verts[1];
1277  quad[3] = verts[0];
1278  }
1279  mPolygonPool->quadFlags(mIdx) = flags;
1280  ++mIdx;
1281  }
1282 
1283  void done() {}
1284 
1285 private:
1286  size_t mIdx;
1287  PolygonPool* mPolygonPool;
1288 };
1289 
1290 
1292 {
1293  AdaptiveMeshOp(): mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL), mTmpPolygonPool() {}
1294 
1295  void init(const size_t upperBound, PolygonPool& polygonPool)
1296  {
1297  mPolygonPool = &polygonPool;
1298 
1299  mTmpPolygonPool.resetQuads(upperBound);
1300  mTmpPolygonPool.resetTriangles(upperBound);
1301 
1302  mQuadIdx = 0;
1303  mTriangleIdx = 0;
1304  }
1305 
1306  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
1307  {
1308  if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
1309  && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
1310  mTmpPolygonPool.quadFlags(mQuadIdx) = flags;
1311  addQuad(verts, reverse);
1312  } else if (
1313  verts[0] == verts[3] &&
1314  verts[1] != verts[2] &&
1315  verts[1] != verts[0] &&
1316  verts[2] != verts[0]) {
1317  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1318  addTriangle(verts[0], verts[1], verts[2], reverse);
1319  } else if (
1320  verts[1] == verts[2] &&
1321  verts[0] != verts[3] &&
1322  verts[0] != verts[1] &&
1323  verts[3] != verts[1]) {
1324  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1325  addTriangle(verts[0], verts[1], verts[3], reverse);
1326  } else if (
1327  verts[0] == verts[1] &&
1328  verts[2] != verts[3] &&
1329  verts[2] != verts[0] &&
1330  verts[3] != verts[0]) {
1331  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1332  addTriangle(verts[0], verts[2], verts[3], reverse);
1333  } else if (
1334  verts[2] == verts[3] &&
1335  verts[0] != verts[1] &&
1336  verts[0] != verts[2] &&
1337  verts[1] != verts[2]) {
1338  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1339  addTriangle(verts[0], verts[1], verts[2], reverse);
1340  }
1341  }
1342 
1343 
1344  void done()
1345  {
1346  mPolygonPool->resetQuads(mQuadIdx);
1347  for (size_t i = 0; i < mQuadIdx; ++i) {
1348  mPolygonPool->quad(i) = mTmpPolygonPool.quad(i);
1349  mPolygonPool->quadFlags(i) = mTmpPolygonPool.quadFlags(i);
1350  }
1351  mTmpPolygonPool.clearQuads();
1352 
1353  mPolygonPool->resetTriangles(mTriangleIdx);
1354  for (size_t i = 0; i < mTriangleIdx; ++i) {
1355  mPolygonPool->triangle(i) = mTmpPolygonPool.triangle(i);
1356  mPolygonPool->triangleFlags(i) = mTmpPolygonPool.triangleFlags(i);
1357  }
1358  mTmpPolygonPool.clearTriangles();
1359  }
1360 
1361 private:
1362 
1363  void addQuad(const Vec4I& verts, bool reverse)
1364  {
1365  if (!reverse) {
1366  mTmpPolygonPool.quad(mQuadIdx) = verts;
1367  } else {
1368  Vec4I& quad = mTmpPolygonPool.quad(mQuadIdx);
1369  quad[0] = verts[3];
1370  quad[1] = verts[2];
1371  quad[2] = verts[1];
1372  quad[3] = verts[0];
1373  }
1374  ++mQuadIdx;
1375  }
1376 
1377  void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
1378  {
1379  Vec3I& prim = mTmpPolygonPool.triangle(mTriangleIdx);
1380 
1381  prim[1] = v1;
1382 
1383  if (!reverse) {
1384  prim[0] = v0;
1385  prim[2] = v2;
1386  } else {
1387  prim[0] = v2;
1388  prim[2] = v0;
1389  }
1390  ++mTriangleIdx;
1391  }
1392 
1393  size_t mQuadIdx, mTriangleIdx;
1394  PolygonPool *mPolygonPool;
1395  PolygonPool mTmpPolygonPool;
1396 };
1397 
1398 
1400 
1401 
1402 template<class DistTreeT, class MeshingOp>
1403 class MeshGen
1404 {
1405 public:
1406  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
1407  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
1408 
1409  MeshGen(const LeafCPtrList<CharTreeT>& edgeLeafs, const IntTreeT& auxTree, PolygonPoolList&);
1411 
1412  void setRefData(const ReferenceData<DistTreeT>&);
1413 
1414  void runParallel();
1415  void runSerial();
1416 
1417  void operator()(const tbb::blocked_range<size_t>&) const;
1418 
1419 private:
1420  const LeafCPtrList<CharTreeT>& mEdgeLeafs;
1421  const IntTreeT& mAuxTree;
1422  const PolygonPoolList& mPolygonPoolList;
1423  size_t mID;
1424  const ReferenceData<DistTreeT>* mRefData;
1425 };
1426 
1427 
1428 template<class DistTreeT, class MeshingOp>
1430  const IntTreeT& auxTree, PolygonPoolList& polygonPoolList)
1431  : mEdgeLeafs(edgeLeafs)
1432  , mAuxTree(auxTree)
1433  , mPolygonPoolList(polygonPoolList)
1434  , mRefData(NULL)
1435 {
1436 }
1437 
1438 
1439 template<class DistTreeT, class MeshingOp>
1441  : mEdgeLeafs(rhs.mEdgeLeafs)
1442  , mAuxTree(rhs.mAuxTree)
1443  , mPolygonPoolList(rhs.mPolygonPoolList)
1444  , mRefData(rhs.mRefData)
1445 {
1446 }
1447 
1448 
1449 template<class DistTreeT, class MeshingOp>
1450 void
1452  const ReferenceData<DistTreeT>& refData)
1453 {
1454  mRefData = &refData;
1455 }
1456 
1457 
1458 template<class DistTreeT, class MeshingOp>
1459 void
1461 {
1462  tbb::parallel_for(mEdgeLeafs.getRange(), *this);
1463 }
1464 
1465 
1466 template<class DistTreeT, class MeshingOp>
1467 void
1469 {
1470  (*this)(mEdgeLeafs.getRange());
1471 }
1472 
1473 
1474 template<class DistTreeT, class MeshingOp>
1475 void
1477  const tbb::blocked_range<size_t>& range) const
1478 {
1479  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1480  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1481  typedef typename CharTreeT::LeafNodeType CharLeafT;
1482 
1483  typedef openvdb::tree::ValueAccessor<const CharTreeT> CharTreeAccessorT;
1484  typedef openvdb::tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
1485  typedef openvdb::tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
1486 
1487  boost::scoped_ptr<CharTreeAccessorT> refEdgeAcc;
1488  boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc;
1489  if (mRefData && mRefData->isValid()) {
1490  refEdgeAcc.reset(new CharTreeAccessorT(*mRefData->mEdgeTree));
1491  refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
1492  }
1493  const bool hasRefData = refEdgeAcc && refMaskAcc;
1494 
1495 
1496  typename CharTreeT::LeafNodeType::ValueOnCIter iter;
1497  IntTreeAccessorT auxAcc(mAuxTree);
1498 
1499  Coord ijk, coord;
1500  char refEdgeFlags, isSemLinePoly;
1501  const char isExteriorPoly[2] = {0, char(POLYFLAG_EXTERIOR)};
1502  openvdb::Vec4I quad;
1503  size_t edgeCount;
1504 
1505  MeshingOp mesher;
1506 
1507  for (size_t n = range.begin(); n != range.end(); ++n) {
1508 
1509  const Coord origin = mEdgeLeafs[n]->getOrigin();
1510 
1511  // Get an upper bound on the number of primitives.
1512  edgeCount = 0;
1513  iter = mEdgeLeafs[n]->cbeginValueOn();
1514  for (; iter; ++iter) {
1515  char edgeFlags = iter.getValue() >> 1;
1516  edgeCount += edgeFlags & 0x1;
1517 
1518  edgeFlags = edgeFlags >> 1;
1519  edgeCount += edgeFlags & 0x1;
1520 
1521  edgeFlags = edgeFlags >> 1;
1522  edgeCount += edgeFlags & 0x1;
1523  }
1524 
1525  mesher.init(edgeCount, mPolygonPoolList[n]);
1526 
1527  const CharLeafT* refEdgeLeaf = NULL;
1528  const BoolLeafT* refMaskLeaf = NULL;
1529 
1530  if (hasRefData) {
1531  refEdgeLeaf = refEdgeAcc->probeConstLeaf(origin);
1532  refMaskLeaf = refMaskAcc->probeConstLeaf(origin);
1533  }
1534 
1535  iter = mEdgeLeafs[n]->cbeginValueOn();
1536  for (; iter; ++iter) {
1537  ijk = iter.getCoord();
1538  const char& edgeFlags = iter.getValue();
1539 
1540  const bool isInside = edgeFlags & INSIDE;
1541 
1542  refEdgeFlags = 0;
1543  isSemLinePoly = 0;
1544  if (hasRefData) {
1545  if(refEdgeLeaf) refEdgeFlags = refEdgeLeaf->getValue(ijk);
1546  if (refMaskLeaf && refMaskLeaf->isValueOn(ijk)) {
1547  isSemLinePoly = char(POLYFLAG_FRACTURE_SEAM);
1548  }
1549  }
1550 
1551 
1552  int v0 = auxAcc.getValue(ijk);
1553 
1554  if (edgeFlags & XEDGE) {
1555 
1556  quad[0] = v0;
1557  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
1558  quad[1] = auxAcc.getValue(coord);
1559  coord[2] -= 1; // i, j-1, k-1
1560  quad[2] = auxAcc.getValue(coord);
1561  coord[1] = ijk[1]; // i, j, k-1
1562  quad[3] = auxAcc.getValue(coord);
1563 
1564  mesher.addPrim(quad, isInside,
1565  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & XEDGE)]));
1566  }
1567 
1568 
1569  if (edgeFlags & YEDGE) {
1570 
1571  quad[0] = v0;
1572  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1; // i, j, k-1
1573  quad[1] = auxAcc.getValue(coord);
1574  coord[0] -= 1; // i-1, j, k-1
1575  quad[2] = auxAcc.getValue(coord);
1576  coord[2] = ijk[2]; // i-1, j, k
1577  quad[3] = auxAcc.getValue(coord);
1578 
1579  mesher.addPrim(quad, isInside,
1580  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & YEDGE)]));
1581  }
1582 
1583  if (edgeFlags & ZEDGE) {
1584 
1585  quad[0] = v0;
1586  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
1587  quad[1] = auxAcc.getValue(coord);
1588  coord[0] -= 1; // i-1, j-1, k
1589  quad[2] = auxAcc.getValue(coord);
1590  coord[1] = ijk[1]; // i, j, k
1591  quad[3] = auxAcc.getValue(coord);
1592 
1593  mesher.addPrim(quad, !isInside,
1594  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & ZEDGE)]));
1595  }
1596  }
1597 
1598  mesher.done();
1599  }
1600 }
1601 
1602 
1604 
1605 
1606 template<class DistTreeT, class AuxDataT = int>
1608 {
1609 public:
1610  typedef openvdb::tree::ValueAccessor<const DistTreeT> SourceAccessorT;
1611  typedef typename DistTreeT::ValueType ValueT;
1612 
1613  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
1614  typedef openvdb::tree::ValueAccessor<CharTreeT> EdgeAccessorT;
1615 
1616  typedef typename DistTreeT::template ValueConverter<AuxDataT>::Type AuxTreeT;
1617  typedef openvdb::tree::ValueAccessor<AuxTreeT> AuxAccessorT;
1618 
1619  AuxiliaryData(const DistTreeT&, const LeafCPtrList<DistTreeT>&,
1620  double iso = 0.0, bool extraCheck = false);
1621  AuxiliaryData(AuxiliaryData&, tbb::split);
1622 
1623  void runParallel();
1624  void runSerial();
1625 
1626  typename CharTreeT::Ptr edgeTree() const { return mEdgeTree; }
1627  typename AuxTreeT::Ptr auxTree() const { return mAuxTree; }
1628 
1629  void operator()(const tbb::blocked_range<size_t>&);
1630 
1631  void join(const AuxiliaryData& rhs)
1632  {
1633  mEdgeTree->merge(*rhs.mEdgeTree);
1634  mAuxTree->merge(*rhs.mAuxTree);
1635  }
1636 
1637 private:
1638  const LeafCPtrList<DistTreeT>& mLeafNodes;
1639  const DistTreeT& mSourceTree;
1640  SourceAccessorT mSourceAccessor;
1641 
1642  typename CharTreeT::Ptr mEdgeTree;
1643  EdgeAccessorT mEdgeAccessor;
1644 
1645  typename AuxTreeT::Ptr mAuxTree;
1646  AuxAccessorT mAuxAccessor;
1647 
1648  const double mIsovalue;
1649  const bool mExtraCheck;
1650 
1651  int edgeCheck(const Coord& ijk, const bool thisInside);
1652 };
1653 
1654 template<class DistTreeT, class AuxDataT>
1656  const LeafCPtrList<DistTreeT>& leafNodes, double iso, bool extraCheck)
1657  : mLeafNodes(leafNodes)
1658  , mSourceTree(tree)
1659  , mSourceAccessor(mSourceTree)
1660  , mEdgeTree(new CharTreeT(0))
1661  , mEdgeAccessor(*mEdgeTree)
1662  , mAuxTree(new AuxTreeT(AuxDataT(0)))
1663  , mAuxAccessor(*mAuxTree)
1664  , mIsovalue(iso)
1665  , mExtraCheck(extraCheck)
1666 {
1667 }
1668 
1669 template<class DistTreeT, class AuxDataT>
1671  : mLeafNodes(rhs.mLeafNodes)
1672  , mSourceTree(rhs.mSourceTree)
1673  , mSourceAccessor(mSourceTree)
1674  , mEdgeTree(new CharTreeT(0))
1675  , mEdgeAccessor(*mEdgeTree)
1676  , mAuxTree(new AuxTreeT(AuxDataT(0)))
1677  , mAuxAccessor(*mAuxTree)
1678  , mIsovalue(rhs.mIsovalue)
1679  , mExtraCheck(rhs.mExtraCheck)
1680 {
1681 }
1682 
1683 
1684 
1685 template<class DistTreeT, typename AuxDataT>
1686 void
1688 {
1689  tbb::parallel_reduce(mLeafNodes.getRange(), *this);
1690 }
1691 
1692 template<class DistTreeT, typename AuxDataT>
1693 void
1695 {
1696  (*this)(mLeafNodes.getRange());
1697 }
1698 
1699 template<class DistTreeT, typename AuxDataT>
1700 void
1701 AuxiliaryData<DistTreeT, AuxDataT>::operator()(const tbb::blocked_range<size_t>& range)
1702 {
1703  typename DistTreeT::LeafNodeType::ValueOnCIter iter;
1704  Coord ijk;
1705  bool thisInside;
1706  int edgeFlags;
1707  ValueT val;
1708 
1709  if (!mExtraCheck) {
1710  for (size_t n = range.begin(); n != range.end(); ++n) {
1711  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
1712  ijk = iter.getCoord();
1713  thisInside = iter.getValue() < mIsovalue;
1714  edgeFlags = edgeCheck(ijk, thisInside);
1715 
1716  if (edgeFlags != 0) {
1717  edgeFlags |= int(thisInside);
1718  mEdgeAccessor.setValue(ijk, char(edgeFlags));
1719  }
1720  }
1721  }
1722  } else {
1723  for (size_t n = range.begin(); n != range.end(); ++n) {
1724  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
1725 
1726  ijk = iter.getCoord();
1727  thisInside = iter.getValue() < mIsovalue;
1728  edgeFlags = edgeCheck(ijk, thisInside);
1729 
1730  if (edgeFlags != 0) {
1731  edgeFlags |= int(thisInside);
1732  mEdgeAccessor.setValue(ijk, char(edgeFlags));
1733  }
1734 
1735  --ijk[0];
1736  if (!mSourceAccessor.probeValue(ijk, val)) {
1737  thisInside = val < mIsovalue;
1738  edgeFlags = edgeCheck(ijk, thisInside);
1739 
1740  if (edgeFlags != 0) {
1741  edgeFlags |= int(thisInside);
1742  mEdgeAccessor.setValue(ijk, char(edgeFlags));
1743  }
1744  }
1745 
1746  ++ijk[0];
1747  --ijk[1];
1748  if (!mSourceAccessor.probeValue(ijk, val)) {
1749  thisInside = val < mIsovalue;
1750  edgeFlags = edgeCheck(ijk, thisInside);
1751 
1752  if (edgeFlags != 0) {
1753  edgeFlags |= int(thisInside);
1754  mEdgeAccessor.setValue(ijk, char(edgeFlags));
1755  }
1756  }
1757 
1758  ++ijk[1];
1759  --ijk[2];
1760  if (!mSourceAccessor.probeValue(ijk, val)) {
1761  thisInside = val < mIsovalue;
1762  edgeFlags = edgeCheck(ijk, thisInside);
1763 
1764  if (edgeFlags != 0) {
1765  edgeFlags |= int(thisInside);
1766  mEdgeAccessor.setValue(ijk, char(edgeFlags));
1767  }
1768  }
1769  }
1770  }
1771  }
1772 }
1773 
1774 template<class DistTreeT, typename AuxDataT>
1775 inline int
1776 AuxiliaryData<DistTreeT, AuxDataT>::edgeCheck(const Coord& ijk, const bool thisInside)
1777 {
1778  int edgeFlags = 0;
1779  Coord n_ijk, coord;
1780 
1781  // Eval upwind x-edge
1782  n_ijk = ijk; ++n_ijk[0];
1783  bool otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
1784  if (otherInside != thisInside) {
1785 
1786  edgeFlags = XEDGE;
1787 
1788  mAuxAccessor.setActiveState(ijk, true);
1789 
1790  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1791  mAuxAccessor.setActiveState(coord, true);
1792 
1793  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
1794  mAuxAccessor.setActiveState(coord, true);
1795 
1796  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1797  mAuxAccessor.setActiveState(coord, true);
1798  }
1799 
1800  // Eval upwind y-edge
1801  n_ijk[0] = ijk[0]; ++n_ijk[1];
1802  otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
1803  if (otherInside != thisInside) {
1804 
1805  edgeFlags |= YEDGE;
1806 
1807  mAuxAccessor.setActiveState(ijk, true);
1808 
1809  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1810  mAuxAccessor.setActiveState(coord, true);
1811 
1812  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
1813  mAuxAccessor.setActiveState(coord, true);
1814 
1815  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1816  mAuxAccessor.setActiveState(coord, true);
1817  }
1818 
1819  // Eval upwind z-edge
1820  n_ijk[1] = ijk[1]; ++n_ijk[2];
1821  otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
1822  if (otherInside != thisInside) {
1823 
1824  edgeFlags |= ZEDGE;
1825 
1826  mAuxAccessor.setActiveState(ijk, true);
1827 
1828  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1829  mAuxAccessor.setActiveState(coord, true);
1830 
1831  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
1832  mAuxAccessor.setActiveState(coord, true);
1833 
1834  coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1835  mAuxAccessor.setActiveState(coord, true);
1836  }
1837  return edgeFlags;
1838 }
1839 
1840 
1842 
1843 
1844 template <class DistTreeT>
1846 {
1847 public:
1848  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1849  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
1852 
1853  SeamMaskGen(LeafPtrList<BoolTreeT>& seamMaskLeafs,
1854  const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree);
1855 
1857 
1858  void runParallel();
1859  void runSerial();
1860 
1861  void operator()(const tbb::blocked_range<size_t>&) const;
1862 
1863 private:
1864  LeafPtrList<BoolTreeT>& mSeamMaskLeafs;
1865  const BoolTreeT& mTopologyMaskTree;
1866  BoolTreeAccessorT mTopologyMaskAcc;
1867  const IntTreeT& mAuxTree;
1868  IntTreeAccessorT mAuxAcc;
1869 };
1870 
1871 
1872 template <class DistTreeT>
1874  const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree)
1875  : mSeamMaskLeafs(seamMaskLeafs)
1876  , mTopologyMaskTree(topologyMaskTree)
1877  , mTopologyMaskAcc(mTopologyMaskTree)
1878  , mAuxTree(auxTree)
1879  , mAuxAcc(mAuxTree)
1880 {
1881 }
1882 
1883 template <class DistTreeT>
1885  : mSeamMaskLeafs(rhs.mSeamMaskLeafs)
1886  , mTopologyMaskTree(rhs.mTopologyMaskTree)
1887  , mTopologyMaskAcc(mTopologyMaskTree)
1888  , mAuxTree(rhs.mAuxTree)
1889  , mAuxAcc(mAuxTree)
1890 {
1891 }
1892 
1893 template <class DistTreeT>
1894 void
1896 {
1897  tbb::parallel_for(mSeamMaskLeafs.getRange(), *this);
1898 }
1899 
1900 template <class DistTreeT>
1901 void
1903 {
1904  (*this)(mSeamMaskLeafs.getRange());
1905 }
1906 
1907 template <class DistTreeT>
1908 void
1909 SeamMaskGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
1910 {
1911  typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
1912  Coord ijk, n_ijk;
1913  for (size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
1914  ValueOnIterT it = mSeamMaskLeafs[leafIdx]->beginValueOn();
1915  for (; it; ++it) {
1916  ijk = it.getCoord();
1917  if (!mTopologyMaskAcc.isValueOn(ijk)) {
1918  it.setValueOff();
1919  } else {
1920  bool turnOff = true;
1921  for (size_t n = 0; n < 6; ++n) {
1922  n_ijk = ijk + util::COORD_OFFSETS[n];
1923  if (!mAuxTree.isValueOn(n_ijk) && mTopologyMaskAcc.isValueOn(n_ijk)) {
1924  turnOff = false;
1925  break;
1926  }
1927  }
1928  if (turnOff) it.setValueOff();
1929  }
1930  }
1931  }
1932 }
1933 
1934 
1936 
1937 
1938 template<class CharAccessorT, typename AuxAccessorT>
1940 {
1941 public:
1942  AuxDataGenerator(CharAccessorT& edgeAcc, AuxAccessorT& auxAcc)
1943  : mEdgeAcc(edgeAcc), mAuxAcc(auxAcc) {}
1944 
1945 
1946  void setXEdge(char edgeFlags, const Coord& ijk)
1947  {
1948  mEdgeAcc.setValue(ijk, edgeFlags | XEDGE);
1949 
1950  mAuxAcc.setActiveState(ijk, true);
1951 
1952  Coord coord = ijk;
1953  coord[1] = ijk[1]-1;
1954  mAuxAcc.setActiveState(coord, true);
1955 
1956  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
1957  mAuxAcc.setActiveState(coord, true);
1958 
1959  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1960  mAuxAcc.setActiveState(coord, true);
1961  }
1962 
1963  void setYEdge(char edgeFlags, const Coord& ijk)
1964  {
1965  mEdgeAcc.setValue(ijk, edgeFlags | YEDGE);
1966 
1967  mAuxAcc.setActiveState(ijk, true);
1968 
1969  Coord coord = ijk;
1970  coord[2] = ijk[2]-1;
1971  mAuxAcc.setActiveState(coord, true);
1972 
1973  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
1974  mAuxAcc.setActiveState(coord, true);
1975 
1976  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
1977  mAuxAcc.setActiveState(coord, true);
1978  }
1979 
1980  void setZEdge(char edgeFlags, const Coord& ijk)
1981  {
1982  mEdgeAcc.setValue(ijk, edgeFlags | ZEDGE);
1983 
1984  mAuxAcc.setActiveState(ijk, true);
1985 
1986  Coord coord = ijk;
1987  coord[1] = ijk[1]-1;
1988  mAuxAcc.setActiveState(coord, true);
1989 
1990  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
1991  mAuxAcc.setActiveState(coord, true);
1992 
1993  coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
1994  mAuxAcc.setActiveState(coord, true);
1995  }
1996 
1997 private:
1998  CharAccessorT& mEdgeAcc;
1999  AuxAccessorT& mAuxAcc;
2000 };
2001 
2002 
2004 
2005 
2006 template<class DistTreeT, class AuxTreeT, class CharTreeT>
2007 inline void
2009  const DistTreeT& distTree, CharTreeT& edgeTree, AuxTreeT& auxTree,
2010  double iso)
2011 {
2012  typedef tree::ValueAccessor<const DistTreeT> DistAccessorT;
2013  typedef tree::ValueAccessor<AuxTreeT> AuxTreeAccessorT;
2014  typedef tree::ValueAccessor<CharTreeT> CharTreeAccessorT;
2015 
2016  typename DistTreeT::ValueType isoValue = typename DistTreeT::ValueType(iso);
2017 
2018  DistAccessorT distAcc(distTree);
2019  CharTreeAccessorT edgeAcc(edgeTree);
2020  AuxTreeAccessorT auxAcc(auxTree);
2021 
2023 
2024  Coord ijk, nijk;
2025  typename DistTreeT::ValueType value;
2026  CoordBBox bbox;
2027  bool processTileFace;
2028 
2029  typename DistTreeT::ValueOnCIter tileIter(distTree);
2030  tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
2031 
2032  for ( ; tileIter; ++tileIter) {
2033  tileIter.getBoundingBox(bbox);
2034 
2035  const bool thisInside = (distAcc.getValue(bbox.min()) < isoValue);
2036  const int thisDepth = distAcc.getValueDepth(bbox.min());
2037 
2038 
2039  // Eval x-edges
2040 
2041  ijk = bbox.max();
2042  nijk = ijk;
2043  ++nijk[0];
2044 
2045  processTileFace = true;
2046  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2047  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2048  }
2049 
2050  if (processTileFace) {
2051  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2052  nijk[1] = ijk[1];
2053  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2054  nijk[2] = ijk[2];
2055  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2056  auxData.setXEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2057  }
2058  }
2059  }
2060  }
2061 
2062  ijk = bbox.min();
2063  --ijk[0];
2064 
2065  processTileFace = true;
2066  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2067  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2068  }
2069 
2070  if (processTileFace) {
2071  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2072  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2073  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2074  auxData.setXEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2075  }
2076  }
2077  }
2078  }
2079 
2080 
2081  // Eval y-edges
2082 
2083  ijk = bbox.max();
2084  nijk = ijk;
2085  ++nijk[1];
2086 
2087  processTileFace = true;
2088  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2089  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2090  }
2091 
2092  if (processTileFace) {
2093  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2094  nijk[0] = ijk[0];
2095  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2096  nijk[2] = ijk[2];
2097 
2098  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2099  auxData.setYEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2100  }
2101  }
2102  }
2103  }
2104 
2105 
2106  ijk = bbox.min();
2107  --ijk[1];
2108 
2109  processTileFace = true;
2110  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2111  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2112  }
2113 
2114  if (processTileFace) {
2115  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2116  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2117 
2118  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2119  auxData.setYEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2120  }
2121  }
2122  }
2123  }
2124 
2125 
2126  // Eval z-edges
2127 
2128  ijk = bbox.max();
2129  nijk = ijk;
2130  ++nijk[2];
2131 
2132  processTileFace = true;
2133  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2134  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2135  }
2136 
2137  if (processTileFace) {
2138  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2139  nijk[0] = ijk[0];
2140  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2141  nijk[1] = ijk[1];
2142 
2143  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2144  auxData.setZEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2145  }
2146  }
2147  }
2148  }
2149 
2150  ijk = bbox.min();
2151  --ijk[2];
2152 
2153  processTileFace = true;
2154  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2155  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2156  }
2157 
2158  if (processTileFace) {
2159  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2160  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2161 
2162  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2163  auxData.setZEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2164  }
2165  }
2166  }
2167  }
2168  }
2169 }
2170 
2171 
2172 } // namespace internal
2173 
2174 
2176 
2177 
2178 inline VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity)
2179  : mPointListSize(0)
2180  , mPolygonPoolListSize(0)
2181  , mIsovalue(isovalue)
2182  , mPrimAdaptivity(adaptivity)
2183  , mSecAdaptivity(0.0)
2184  , mRefGrid(GridBase::ConstPtr())
2185  , mRefEdgeTree(TreeBase::Ptr())
2186  , mRefTopologyMaskTree(TreeBase::Ptr())
2187 {
2188 }
2189 
2190 inline PointList&
2192 {
2193  return mPoints;
2194 }
2195 
2196 inline const size_t&
2198 {
2199  return mPointListSize;
2200 }
2201 
2202 
2203 inline PolygonPoolList&
2205 {
2206  return mPolygons;
2207 }
2208 
2209 
2210 inline const PolygonPoolList&
2212 {
2213  return mPolygons;
2214 }
2215 
2216 
2217 inline const size_t&
2219 {
2220  return mPolygonPoolListSize;
2221 }
2222 
2223 
2224 inline void
2225 VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity)
2226 {
2227  mRefGrid = grid;
2228  mSecAdaptivity = secAdaptivity;
2229  // Clear out old auxiliary data
2230  mRefEdgeTree = TreeBase::Ptr();
2231  mRefTopologyMaskTree = TreeBase::Ptr();
2232 }
2233 
2234 
2235 template<typename GridT>
2236 inline void
2237 VolumeToMesh::operator()(const GridT& distGrid)
2238 {
2239  typedef typename GridT::TreeType DistTreeT;
2240  typedef typename DistTreeT::ValueType DistValueT;
2241 
2242  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
2243  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
2244  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
2245 
2246  const bool noAdaptivity = mPrimAdaptivity < 1e-6 && mSecAdaptivity < 1e-6;
2247 
2248  const openvdb::math::Transform& transform = distGrid.transform();
2249  const DistTreeT& distTree = distGrid.tree();
2250  typename CharTreeT::Ptr edgeTree; // edge flags
2251  typename IntTreeT::Ptr auxTree; // auxiliary data
2252 
2253  const bool nonLevelSetGrid = distGrid.getGridClass() != GRID_LEVEL_SET;
2254 
2255  const bool extraSignCheck = nonLevelSetGrid ||
2256  (std::abs(mIsovalue - double(distGrid.background())) < (1.5 * transform.voxelSize()[0]));
2257 
2258  // Collect auxiliary data
2259  {
2260 
2261  internal::LeafCPtrList<DistTreeT> sourceLeafs(distTree);
2262  internal::AuxiliaryData<DistTreeT> op(distTree, sourceLeafs, mIsovalue, extraSignCheck);
2263  op.runParallel();
2264  edgeTree = op.edgeTree();
2265  auxTree = op.auxTree();
2266 
2267  // Collect auxiliary data from active tiles
2268  if (nonLevelSetGrid) {
2269  internal::tileAuxiliaryData(distTree, *edgeTree, *auxTree, mIsovalue);
2270  }
2271  }
2272 
2273  // Optionally collect auxiliary data from a reference surface.
2275  if (mRefGrid) {
2276  const GridT* refGrid = static_cast<const GridT*>(mRefGrid.get());
2277  if (refGrid && refGrid->activeVoxelCount() > 0) {
2278 
2279  refData.mDistTree = &refGrid->tree();
2280  refData.mInternalAdaptivity = DistValueT(mSecAdaptivity);
2281 
2282  // Cache reference data for reuse.
2283  if (!mRefEdgeTree && !mRefTopologyMaskTree) {
2286  *refData.mDistTree, leafs, mIsovalue, extraSignCheck);
2287  op.runParallel();
2288  mRefEdgeTree = op.edgeTree();
2289  mRefTopologyMaskTree = op.auxTree();
2290  }
2291 
2292  if (mRefEdgeTree && mRefTopologyMaskTree) {
2293  refData.mEdgeTree = static_cast<CharTreeT*>(mRefEdgeTree.get());
2294  refData.mTopologyMaskTree = static_cast<BoolTreeT*>(mRefTopologyMaskTree.get());
2295  refData.mSeamMaskTree = typename BoolTreeT::Ptr(new BoolTreeT(false));
2296  }
2297  }
2298  }
2299 
2300  // Generate the seamline mask
2301  if (refData.mSeamMaskTree) {
2302  refData.mSeamMaskTree->topologyUnion(*auxTree.get());
2303 
2304  internal::LeafPtrList<BoolTreeT> leafs(*refData.mSeamMaskTree.get());
2305  internal::SeamMaskGen<DistTreeT> op(leafs, *refData.mTopologyMaskTree, *auxTree.get());
2306  op.runParallel();
2307 
2308  refData.mSeamMaskTree->pruneInactive();
2309  dilateVoxels(*refData.mSeamMaskTree);
2310  dilateVoxels(*refData.mSeamMaskTree);
2311  dilateVoxels(*refData.mSeamMaskTree);
2312  }
2313 
2314  // Process auxiliary data
2315  {
2316  internal::LeafPtrList<IntTreeT> auxLeafs(*auxTree);
2317  std::vector<size_t> regions(auxLeafs.size(), 0);
2318 
2319  {
2320  if (noAdaptivity) {
2321  internal::Count<IntTreeT> count(auxLeafs, regions);
2322  count.runParallel();
2323  } else {
2324  internal::Merge<DistTreeT> merge(distTree, auxLeafs,
2325  regions, DistValueT(mIsovalue), DistValueT(mPrimAdaptivity));
2326  merge.setRefData(refData);
2327  merge.runParallel();
2328  }
2329 
2330  mPointListSize = 0;
2331  size_t tmp = 0;
2332  for (size_t n = 0, N = regions.size(); n < N; ++n) {
2333  tmp = regions[n];
2334  regions[n] = mPointListSize;
2335  mPointListSize += tmp;
2336  }
2337  }
2338 
2339  // Generate the unique point list
2340  mPoints.reset(new openvdb::Vec3s[mPointListSize]);
2341 
2343  pointGen(distTree, auxLeafs, regions, transform, mPoints, mIsovalue);
2344  pointGen.setRefData(refData);
2345  pointGen.runParallel();
2346  }
2347 
2348  // Generate mesh
2349  {
2350  internal::LeafCPtrList<CharTreeT> edgeLeafs(*edgeTree);
2351  mPolygonPoolListSize = edgeLeafs.size();
2352  mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
2353 
2354  if (noAdaptivity) {
2356  meshGen(edgeLeafs, *auxTree, mPolygons);
2357  meshGen.setRefData(refData);
2358  meshGen.runParallel();
2359  } else {
2361  meshGen(edgeLeafs, *auxTree, mPolygons);
2362  meshGen.setRefData(refData);
2363  meshGen.runParallel();
2364  }
2365  }
2366 }
2367 
2368 } // namespace tools
2369 } // namespace OPENVDB_VERSION_NAME
2370 } // namespace openvdb
2371 
2372 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
2373 
2374 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2375 // All rights reserved. This software is distributed under the
2376 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )