OpenVDB  4.0.1
PoissonSolver.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2017 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 //
87 
88 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
89 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
90 
91 #include <openvdb/Types.h>
94 #include <openvdb/tree/Tree.h>
96 #include "Morphology.h" // for erodeVoxels
97 #include <boost/scoped_array.hpp>
98 
99 
100 namespace openvdb {
102 namespace OPENVDB_VERSION_NAME {
103 namespace tools {
104 namespace poisson {
105 
106 // This type should be at least as wide as math::pcg::SizeType.
107 typedef Int32 VIndex;
108 
111 
112 
114 template<typename TreeType>
127 inline typename TreeType::Ptr
128 solve(const TreeType&, math::pcg::State&, bool staggered = false);
129 
130 template<typename TreeType, typename Interrupter>
131 inline typename TreeType::Ptr
132 solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
134 
135 
137 template<typename TreeType, typename BoundaryOp, typename Interrupter>
168 inline typename TreeType::Ptr
170  const TreeType&,
171  const BoundaryOp&,
173  Interrupter&,
174  bool staggered = false);
175 
176 template<
177  typename PreconditionerType,
178  typename TreeType,
179  typename BoundaryOp,
180  typename Interrupter>
181 inline typename TreeType::Ptr
183  const TreeType&,
184  const BoundaryOp&,
186  Interrupter&,
187  bool staggered = false);
188 
189 template<
190  typename PreconditionerType,
191  typename TreeType,
192  typename DomainTreeType,
193  typename BoundaryOp,
194  typename Interrupter>
195 inline typename TreeType::Ptr
197  const TreeType&,
198  const DomainTreeType&,
199  const BoundaryOp&,
201  Interrupter&,
202  bool staggered = false);
204 
205 
207 
208 // The following are low-level routines that can be used to assemble custom solvers.
209 
212 template<typename VIndexTreeType>
213 inline void populateIndexTree(VIndexTreeType&);
214 
218 template<typename TreeType>
219 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
220 createIndexTree(const TreeType&);
221 
222 
229 template<typename VectorValueType, typename SourceTreeType>
232  const SourceTreeType& source,
233  const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
234 
235 
243 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
244 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
247  const VIndexTreeType& index,
248  const TreeValueType& background);
249 
250 
255 template<typename BoolTreeType>
258  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
259  const BoolTreeType& interiorMask,
260  bool staggered = false);
261 
262 
282 template<typename BoolTreeType, typename BoundaryOp>
285  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
286  const BoolTreeType& interiorMask,
287  const BoundaryOp& boundaryOp,
289  bool staggered = false);
290 
291 
295 template<typename ValueType>
297  inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
298  // Exterior neighbors are empty, so decrement the weighting coefficient
299  // as for interior neighbors but leave the source vector unchanged.
300  diag -= 1;
301  }
302 };
303 
305 
306 
308 
309 
310 namespace internal {
311 
314 template<typename LeafType>
316 {
317  VIndex* count;
318  LeafCountOp(VIndex* count_): count(count_) {}
319  void operator()(const LeafType& leaf, size_t leafIdx) const {
320  count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
321  }
322 };
323 
324 
327 template<typename LeafType>
329 {
330  const VIndex* count;
331  LeafIndexOp(const VIndex* count_): count(count_) {}
332  void operator()(LeafType& leaf, size_t leafIdx) const {
333  VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
334  for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
335  it.setValue(idx++);
336  }
337  }
338 };
339 
340 } // namespace internal
341 
342 
343 template<typename VIndexTreeType>
344 inline void
345 populateIndexTree(VIndexTreeType& result)
346 {
347  typedef typename VIndexTreeType::LeafNodeType LeafT;
348  typedef typename tree::LeafManager<VIndexTreeType> LeafMgrT;
349 
350  // Linearize the tree.
351  LeafMgrT leafManager(result);
352  const size_t leafCount = leafManager.leafCount();
353 
354  if (leafCount == 0) return;
355 
356  // Count the number of active voxels in each leaf node.
357  boost::scoped_array<VIndex> perLeafCount(new VIndex[leafCount]);
358  VIndex* perLeafCountPtr = perLeafCount.get();
359  leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
360 
361  // The starting index for each leaf node is the total number
362  // of active voxels in all preceding leaf nodes.
363  for (size_t i = 1; i < leafCount; ++i) {
364  perLeafCount[i] += perLeafCount[i - 1];
365  }
366 
367  // The last accumulated value should be the total of all active voxels.
368  assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
369 
370  // Parallelize over the leaf nodes of the tree, storing a unique index
371  // in each active voxel.
372  leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
373 }
374 
375 
376 template<typename TreeType>
377 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
378 createIndexTree(const TreeType& tree)
379 {
380  typedef typename TreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
381 
382  // Construct an output tree with the same active voxel topology as the input tree.
383  const VIndex invalidIdx = -1;
384  typename VIdxTreeT::Ptr result(
385  new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
386 
387  // All active voxels are degrees of freedom, including voxels contained in active tiles.
388  result->voxelizeActiveTiles();
389 
390  populateIndexTree(*result);
391 
392  return result;
393 }
394 
395 
397 
398 
399 namespace internal {
400 
403 template<typename VectorValueType, typename SourceTreeType>
405 {
406  typedef typename SourceTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
407  typedef typename VIdxTreeT::LeafNodeType VIdxLeafT;
408  typedef typename SourceTreeType::LeafNodeType LeafT;
409  typedef typename SourceTreeType::ValueType TreeValueT;
411 
412  const SourceTreeType* tree;
413  VectorT* vector;
414 
415  CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
416 
417  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
418  {
419  VectorT& vec = *vector;
420  if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
421  // If a corresponding leaf node exists in the source tree,
422  // copy voxel values from the source node to the output vector.
423  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
424  vec[*it] = leaf->getValue(it.pos());
425  }
426  } else {
427  // If no corresponding leaf exists in the source tree,
428  // fill the vector with a uniform value.
429  const TreeValueT& value = tree->getValue(idxLeaf.origin());
430  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
431  vec[*it] = value;
432  }
433  }
434  }
435 };
436 
437 } // namespace internal
438 
439 
440 template<typename VectorValueType, typename SourceTreeType>
442 createVectorFromTree(const SourceTreeType& tree,
443  const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
444 {
445  typedef typename SourceTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
446  typedef tree::LeafManager<const VIdxTreeT> VIdxLeafMgrT;
447  typedef typename math::pcg::Vector<VectorValueType> VectorT;
448 
449  // Allocate the vector.
450  const size_t numVoxels = idxTree.activeVoxelCount();
451  typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
452 
453  // Parallelize over the leaf nodes of the index tree, filling the output vector
454  // with values from corresponding voxels of the source tree.
455  VIdxLeafMgrT leafManager(idxTree);
456  leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
457 
458  return result;
459 }
460 
461 
463 
464 
465 namespace internal {
466 
469 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
471 {
472  typedef typename VIndexTreeType::template ValueConverter<TreeValueType>::Type OutTreeT;
473  typedef typename OutTreeT::LeafNodeType OutLeafT;
474  typedef typename VIndexTreeType::LeafNodeType VIdxLeafT;
476 
477  const VectorT* vector;
478  OutTreeT* tree;
479 
480  CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
481 
482  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
483  {
484  const VectorT& vec = *vector;
485  OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
486  assert(leaf != NULL);
487  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
488  leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
489  }
490  }
491 };
492 
493 } // namespace internal
494 
495 
496 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
497 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
500  const VIndexTreeType& idxTree,
501  const TreeValueType& background)
502 {
503  typedef typename VIndexTreeType::template ValueConverter<TreeValueType>::Type OutTreeT;
504  typedef typename tree::LeafManager<const VIndexTreeType> VIdxLeafMgrT;
505 
506  // Construct an output tree with the same active voxel topology as the index tree.
507  typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
508  OutTreeT& tree = *result;
509 
510  // Parallelize over the leaf nodes of the index tree, populating voxels
511  // of the output tree with values from the input vector.
512  VIdxLeafMgrT leafManager(idxTree);
513  leafManager.foreach(
515 
516  return result;
517 }
518 
519 
521 
522 
523 namespace internal {
524 
526 template<typename BoolTreeType, typename BoundaryOp>
528 {
529  typedef typename BoolTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
530  typedef typename VIdxTreeT::LeafNodeType VIdxLeafT;
533 
534  LaplacianMatrix* laplacian;
535  const VIdxTreeT* idxTree;
536  const BoolTreeType* interiorMask;
537  const BoundaryOp boundaryOp;
538  VectorT* source;
539 
540  ISStaggeredLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx,
541  const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
542  laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
543 
544  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
545  {
546  // Local accessors
547  typename tree::ValueAccessor<const BoolTreeType> interior(*interiorMask);
548  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
549 
550  Coord ijk;
551  VIndex column;
552  const ValueT diagonal = -6.f, offDiagonal = 1.f;
553 
554  // Loop over active voxels in this leaf.
555  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
556  assert(it.getValue() > -1);
557  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
558 
559  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
560 
561  ijk = it.getCoord();
562  if (interior.isValueOn(ijk)) {
563  // The current voxel is an interior voxel.
564  // All of its neighbors are in the solution domain.
565 
566  // -x direction
567  row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
568  // -y direction
569  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
570  // -z direction
571  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
572  // diagonal
573  row.setValue(rowNum, diagonal);
574  // +z direction
575  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
576  // +y direction
577  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
578  // +x direction
579  row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
580 
581  } else {
582  // The current voxel is a boundary voxel.
583  // At least one of its neighbors is outside the solution domain.
584 
585  ValueT modifiedDiagonal = 0.f;
586 
587  // -x direction
588  if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
589  row.setValue(column, offDiagonal);
590  modifiedDiagonal -= 1;
591  } else {
592  boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
593  }
594  // -y direction
595  if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
596  row.setValue(column, offDiagonal);
597  modifiedDiagonal -= 1;
598  } else {
599  boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
600  }
601  // -z direction
602  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
603  row.setValue(column, offDiagonal);
604  modifiedDiagonal -= 1;
605  } else {
606  boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
607  }
608  // +z direction
609  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
610  row.setValue(column, offDiagonal);
611  modifiedDiagonal -= 1;
612  } else {
613  boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
614  }
615  // +y direction
616  if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
617  row.setValue(column, offDiagonal);
618  modifiedDiagonal -= 1;
619  } else {
620  boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
621  }
622  // +x direction
623  if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
624  row.setValue(column, offDiagonal);
625  modifiedDiagonal -= 1;
626  } else {
627  boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
628  }
629  // diagonal
630  row.setValue(rowNum, modifiedDiagonal);
631  }
632  } // end loop over voxels
633  }
634 };
635 
636 
637 // Stencil 1 is the correct stencil, but stencil 2 requires
638 // half as many comparisons and produces smoother results at boundaries.
639 //#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
640 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
641 
643 template<typename VIdxTreeT, typename BoundaryOp>
645 {
646  typedef typename VIdxTreeT::LeafNodeType VIdxLeafT;
649 
650  LaplacianMatrix* laplacian;
651  const VIdxTreeT* idxTree;
652  const BoundaryOp boundaryOp;
653  VectorT* source;
654 
655  ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
656  laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
657 
658  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
659  {
660  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
661 
662  const int kNumOffsets = 6;
663  const Coord ijkOffset[kNumOffsets] = {
664 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
665  Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
666 #else
667  Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
668 #endif
669  };
670 
671  // For each active voxel in this leaf...
672  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
673  assert(it.getValue() > -1);
674 
675  const Coord ijk = it.getCoord();
676  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
677 
678  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
679 
680  ValueT modifiedDiagonal = 0.f;
681 
682  // For each of the neighbors of the voxel at (i,j,k)...
683  for (int dir = 0; dir < kNumOffsets; ++dir) {
684  const Coord neighbor = ijk + ijkOffset[dir];
685  VIndex column;
686  // For collocated vector grids, the central differencing stencil requires
687  // access to neighbors at a distance of two voxels in each direction
688  // (-x, +x, -y, +y, -z, +z).
689 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
690  const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
691  && vectorIdx.isValueOn(neighbor));
692 #else
693  const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
694 #endif
695  if (ijkIsInterior) {
696  // If (i,j,k) is sufficiently far away from the exterior,
697  // set its weight to one and adjust the center weight accordingly.
698  row.setValue(column, 1.f);
699  modifiedDiagonal -= 1.f;
700  } else {
701  // If (i,j,k) is adjacent to or one voxel away from the exterior,
702  // invoke the boundary condition functor.
703  boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
704  }
705  }
706  // Set the (possibly modified) weight for the voxel at (i,j,k).
707  row.setValue(rowNum, modifiedDiagonal);
708  }
709  }
710 };
711 
712 } // namespace internal
713 
714 
715 template<typename BoolTreeType>
717 createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
718  const BoolTreeType& interiorMask, bool staggered)
719 {
720  typedef LaplacianMatrix::ValueType ValueT;
722  static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
724  return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
725 }
726 
727 
728 template<typename BoolTreeType, typename BoundaryOp>
731  const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
732  const BoolTreeType& interiorMask,
733  const BoundaryOp& boundaryOp,
735  bool staggered)
736 {
737  typedef typename BoolTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
738  typedef typename tree::LeafManager<const VIdxTreeT> VIdxLeafMgrT;
739 
740  // The number of active voxels is the number of degrees of freedom.
741  const Index64 numDoF = idxTree.activeVoxelCount();
742 
743  // Construct the matrix.
744  LaplacianMatrix::Ptr laplacianPtr(
745  new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
746  LaplacianMatrix& laplacian = *laplacianPtr;
747 
748  // Populate the matrix using a second-order, 7-point CD stencil.
749  VIdxLeafMgrT idxLeafManager(idxTree);
750  if (staggered) {
752  laplacian, idxTree, interiorMask, boundaryOp, source));
753  } else {
754  idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
755  laplacian, idxTree, boundaryOp, source));
756  }
757 
758  return laplacianPtr;
759 }
760 
761 
763 
764 
765 template<typename TreeType>
766 inline typename TreeType::Ptr
767 solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
768 {
769  util::NullInterrupter interrupter;
770  return solve(inTree, state, interrupter, staggered);
771 }
772 
773 
774 template<typename TreeType, typename Interrupter>
775 inline typename TreeType::Ptr
776 solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
777 {
779  return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
780 }
781 
782 
783 template<typename TreeType, typename BoundaryOp, typename Interrupter>
784 inline typename TreeType::Ptr
785 solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
786  math::pcg::State& state, Interrupter& interrupter, bool staggered)
787 {
789  return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
790  inTree, boundaryOp, state, interrupter, staggered);
791 }
792 
793 
794 template<
795  typename PreconditionerType,
796  typename TreeType,
797  typename BoundaryOp,
798  typename Interrupter>
799 inline typename TreeType::Ptr
801  const TreeType& inTree,
802  const BoundaryOp& boundaryOp,
803  math::pcg::State& state,
804  Interrupter& interrupter,
805  bool staggered)
806 {
807  return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
808  /*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
809 }
810 
811 template<
812  typename PreconditionerType,
813  typename TreeType,
814  typename DomainTreeType,
815  typename BoundaryOp,
816  typename Interrupter>
817 inline typename TreeType::Ptr
819  const TreeType& inTree,
820  const DomainTreeType& domainMask,
821  const BoundaryOp& boundaryOp,
822  math::pcg::State& state,
823  Interrupter& interrupter,
824  bool staggered)
825 {
826  typedef typename TreeType::ValueType TreeValueT;
827  typedef LaplacianMatrix::ValueType VecValueT;
828  typedef typename math::pcg::Vector<VecValueT> VectorT;
829  typedef typename TreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
830  typedef typename TreeType::template ValueConverter<bool>::Type MaskTreeT;
831 
832  // 1. Create a mapping from active voxels of the input tree to elements of a vector.
833  typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
834 
835  // 2. Populate a vector with values from the input tree.
836  typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
837 
838  // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
840  typename MaskTreeT::Ptr interiorMask(
841  new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
842  tools::erodeVoxels(*interiorMask, /*iterations=*/1, tools::NN_FACE);
843 
844  // 4. Create the Laplacian matrix.
846  *idxTree, *interiorMask, boundaryOp, *b, staggered);
847 
848  // 5. Solve the Poisson equation.
849  laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
850  b->scale(-1.0);
851  typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
853  new PreconditionerType(*laplacian));
854  if (!precond->isValid()) {
855  precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
856  }
857 
858  state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
859 
860  // 6. Populate the output tree with values from the solution vector.
862  return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
863 }
864 
865 } // namespace poisson
866 } // namespace tools
867 } // namespace OPENVDB_VERSION_NAME
868 } // namespace openvdb
869 
870 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
871 
873 //
874 // Copyright (c) 2012-2017 DreamWorks Animation LLC
875 //
876 // All rights reserved. This software is distributed under the
877 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
878 //
879 // Redistributions of source code must retain the above copyright
880 // and license notice and the following restrictions and disclaimer.
881 //
882 // * Neither the name of DreamWorks Animation nor the names of
883 // its contributors may be used to endorse or promote products derived
884 // from this software without specific prior written permission.
885 //
886 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
887 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
888 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
889 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
890 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
891 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
892 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
893 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
894 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
895 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
896 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
897 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
898 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
899 //
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:503
LeafIndexOp(const VIndex *count_)
Definition: PoissonSolver.h:331
Functor for use with LeafManager::foreach() to populate a tree with values from a vector...
Definition: PoissonSolver.h:470
TreeType::Ptr solve(const TreeType &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x, where b is a vector comprising the values of all of the active voxels in the i...
Definition: PoissonSolver.h:776
Functor for use with LeafManager::foreach() to populate active leaf voxels with sequential indices...
Definition: PoissonSolver.h:328
TreeType::Ptr solveWithBoundaryConditions(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the val...
Definition: PoissonSolver.h:785
SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:495
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero...
Definition: PoissonSolver.h:345
VIndexTreeType::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:474
Dirichlet boundary condition functor.
Definition: PoissonSolver.h:296
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:266
Definition: ValueAccessor.h:219
Lightweight, variable-length vector.
Definition: ConjGradient.h:64
VectorT * source
Definition: PoissonSolver.h:653
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:407
VIndexTreeType::template ValueConverter< TreeValueType >::Type::Ptr createTreeFromVector(const math::pcg::Vector< VectorValueType > &values, const VIndexTreeType &index, const TreeValueType &background)
Return a tree with the same active voxel topology as the index tree but whose voxel values are taken ...
Definition: PoissonSolver.h:498
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:115
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:48
SharedPtr< SparseStencilMatrix > Ptr
Definition: ConjGradient.h:269
LaplacianMatrix::Ptr createISLaplacianWithBoundaryConditions(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, const BoundaryOp &boundaryOp, typename math::pcg::Vector< LaplacianMatrix::ValueType > &source, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator with user-specified boundary...
Definition: PoissonSolver.h:730
CopyFromVecOp(const VectorT &v, OutTreeT &t)
Definition: PoissonSolver.h:480
const VIdxTreeT * idxTree
Definition: PoissonSolver.h:535
VectorT * vector
Definition: PoissonSolver.h:413
VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:530
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:66
math::pcg::Vector< VectorValueType > VectorT
Definition: PoissonSolver.h:475
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
Definition: PoissonSolver.h:110
const VIdxTreeT * idxTree
Definition: PoissonSolver.h:651
math::pcg::Vector< VectorValueType >::Ptr createVectorFromTree(const SourceTreeType &source, const typename SourceTreeType::template ValueConverter< VIndex >::Type &index)
Return a vector of the active voxel values of the scalar-valued source tree.
Definition: PoissonSolver.h:442
VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:646
math::pcg::Vector< ValueT > VectorT
Definition: PoissonSolver.h:648
LaplacianMatrix::ValueType ValueT
Definition: PoissonSolver.h:531
TreeType::Ptr solveWithBoundaryConditionsAndPreconditioner(const TreeType &, const DomainTreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the val...
Definition: PoissonSolver.h:818
VIndex * count
Definition: PoissonSolver.h:317
#define OPENVDB_VERSION_NAME
Definition: version.h:43
LaplacianMatrix * laplacian
Definition: PoissonSolver.h:650
void operator()(LeafType &leaf, size_t leafIdx) const
Definition: PoissonSolver.h:332
SourceTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
Definition: PoissonSolver.h:406
ISLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoundaryOp &op, VectorT &src)
Definition: PoissonSolver.h:655
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:73
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
Diagonal preconditioner.
Definition: ConjGradient.h:69
const BoundaryOp boundaryOp
Definition: PoissonSolver.h:652
Read/write accessor to a row of this matrix.
Definition: ConjGradient.h:446
T & at(SizeType i)
Return the value of this vector&#39;s ith element.
Definition: ConjGradient.h:228
Implementation of morphological dilation and erosion.
LeafCountOp(VIndex *count_)
Definition: PoissonSolver.h:318
Definition: Exceptions.h:39
ValueType_ ValueType
Definition: ConjGradient.h:267
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:417
Functor for use with LeafManager::foreach() to populate an array with per-leaf active voxel counts...
Definition: PoissonSolver.h:315
LaplacianMatrix::Ptr createISLaplacian(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator using second-order finite di...
Definition: PoissonSolver.h:717
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
Definition: PoissonSolver.h:644
SourceTreeType::LeafNodeType LeafT
Definition: PoissonSolver.h:408
Int32 VIndex
Definition: PoissonSolver.h:107
LaplacianMatrix::ValueType ValueT
Definition: PoissonSolver.h:647
ISStaggeredLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoolTreeType &mask, const BoundaryOp &op, VectorT &src)
Definition: PoissonSolver.h:540
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:482
OutTreeT * tree
Definition: PoissonSolver.h:478
const BoolTreeType * interiorMask
Definition: PoissonSolver.h:536
VIndexTreeType::template ValueConverter< TreeValueType >::Type OutTreeT
Definition: PoissonSolver.h:472
uint64_t Index64
Definition: Types.h:56
const BoundaryOp boundaryOp
Definition: PoissonSolver.h:537
BoolTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
Definition: PoissonSolver.h:529
LaplacianMatrix * laplacian
Definition: PoissonSolver.h:534
CopyToVecOp(const SourceTreeType &t, VectorT &v)
Definition: PoissonSolver.h:415
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:544
void operator()(const LeafType &leaf, size_t leafIdx) const
Definition: PoissonSolver.h:319
const VectorT * vector
Definition: PoissonSolver.h:477
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:658
void erodeVoxels(TreeType &tree, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically erode all leaf-level active voxels in the given tree.
Definition: Morphology.h:874
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:70
math::pcg::Vector< ValueT > VectorT
Definition: PoissonSolver.h:532
const SourceTreeType * tree
Definition: PoissonSolver.h:412
Index32 SizeType
Definition: ConjGradient.h:60
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
Definition: PoissonSolver.h:297
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition: ConjGradient.h:1096
TreeType::template ValueConverter< VIndex >::Type::Ptr createIndexTree(const TreeType &)
Iterate over the active voxels of the input tree and for each one assign its index in the iteration s...
Definition: PoissonSolver.h:378
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
int32_t Int32
Definition: Types.h:59
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
const VIndex * count
Definition: PoissonSolver.h:330
OutTreeT::LeafNodeType OutLeafT
Definition: PoissonSolver.h:473
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:1040
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
Definition: PoissonSolver.h:527
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:169
Functor for use with LeafManager::foreach() to populate a vector with the values of a tree&#39;s active v...
Definition: PoissonSolver.h:404
Definition: Morphology.h:87
SourceTreeType::ValueType TreeValueT
Definition: PoissonSolver.h:409
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition: ConjGradient.h:1252
math::pcg::Vector< VectorValueType > VectorT
Definition: PoissonSolver.h:410
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:115