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