OpenVDB  7.0.0
PoissonSolver.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
60 
61 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
63 
64 #include <openvdb/Types.h>
67 #include <openvdb/tree/Tree.h>
69 #include "Morphology.h" // for erodeVoxels
70 
71 
72 namespace openvdb {
74 namespace OPENVDB_VERSION_NAME {
75 namespace tools {
76 namespace poisson {
77 
78 // This type should be at least as wide as math::pcg::SizeType.
79 using VIndex = Int32;
80 
83 
84 
86 template<typename TreeType>
99 inline typename TreeType::Ptr
100 solve(const TreeType&, math::pcg::State&, bool staggered = false);
101 
102 template<typename TreeType, typename Interrupter>
103 inline typename TreeType::Ptr
104 solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
106 
107 
109 template<typename TreeType, typename BoundaryOp, typename Interrupter>
140 inline typename TreeType::Ptr
142  const TreeType&,
143  const BoundaryOp&,
145  Interrupter&,
146  bool staggered = false);
147 
148 template<
149  typename PreconditionerType,
150  typename TreeType,
151  typename BoundaryOp,
152  typename Interrupter>
153 inline typename TreeType::Ptr
155  const TreeType&,
156  const BoundaryOp&,
158  Interrupter&,
159  bool staggered = false);
160 
161 template<
162  typename PreconditionerType,
163  typename TreeType,
164  typename DomainTreeType,
165  typename BoundaryOp,
166  typename Interrupter>
167 inline typename TreeType::Ptr
169  const TreeType&,
170  const DomainTreeType&,
171  const BoundaryOp&,
173  Interrupter&,
174  bool staggered = false);
176 
177 
179 
180 // The following are low-level routines that can be used to assemble custom solvers.
181 
184 template<typename VIndexTreeType>
185 inline void populateIndexTree(VIndexTreeType&);
186 
190 template<typename TreeType>
191 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
192 createIndexTree(const TreeType&);
193 
194 
201 template<typename VectorValueType, typename SourceTreeType>
204  const SourceTreeType& source,
205  const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
206 
207 
215 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
216 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
219  const VIndexTreeType& index,
220  const TreeValueType& background);
221 
222 
227 template<typename BoolTreeType>
230  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
231  const BoolTreeType& interiorMask,
232  bool staggered = false);
233 
234 
254 template<typename BoolTreeType, typename BoundaryOp>
257  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
258  const BoolTreeType& interiorMask,
259  const BoundaryOp& boundaryOp,
261  bool staggered = false);
262 
263 
267 template<typename ValueType>
269  inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
270  // Exterior neighbors are empty, so decrement the weighting coefficient
271  // as for interior neighbors but leave the source vector unchanged.
272  diag -= 1;
273  }
274 };
275 
277 
278 
280 
281 
282 namespace internal {
283 
286 template<typename LeafType>
288 {
290  LeafCountOp(VIndex* count_): count(count_) {}
291  void operator()(const LeafType& leaf, size_t leafIdx) const {
292  count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
293  }
294 };
295 
296 
299 template<typename LeafType>
301 {
302  const VIndex* count;
303  LeafIndexOp(const VIndex* count_): count(count_) {}
304  void operator()(LeafType& leaf, size_t leafIdx) const {
305  VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
306  for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
307  it.setValue(idx++);
308  }
309  }
310 };
311 
312 } // namespace internal
313 
314 
315 template<typename VIndexTreeType>
316 inline void
317 populateIndexTree(VIndexTreeType& result)
318 {
319  using LeafT = typename VIndexTreeType::LeafNodeType;
320  using LeafMgrT = typename tree::LeafManager<VIndexTreeType>;
321 
322  // Linearize the tree.
323  LeafMgrT leafManager(result);
324  const size_t leafCount = leafManager.leafCount();
325 
326  if (leafCount == 0) return;
327 
328  // Count the number of active voxels in each leaf node.
329  std::unique_ptr<VIndex[]> perLeafCount(new VIndex[leafCount]);
330  VIndex* perLeafCountPtr = perLeafCount.get();
331  leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
332 
333  // The starting index for each leaf node is the total number
334  // of active voxels in all preceding leaf nodes.
335  for (size_t i = 1; i < leafCount; ++i) {
336  perLeafCount[i] += perLeafCount[i - 1];
337  }
338 
339  // The last accumulated value should be the total of all active voxels.
340  assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
341 
342  // Parallelize over the leaf nodes of the tree, storing a unique index
343  // in each active voxel.
344  leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
345 }
346 
347 
348 template<typename TreeType>
349 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
350 createIndexTree(const TreeType& tree)
351 {
352  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
353 
354  // Construct an output tree with the same active voxel topology as the input tree.
355  const VIndex invalidIdx = -1;
356  typename VIdxTreeT::Ptr result(
357  new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
358 
359  // All active voxels are degrees of freedom, including voxels contained in active tiles.
360  result->voxelizeActiveTiles();
361 
362  populateIndexTree(*result);
363 
364  return result;
365 }
366 
367 
369 
370 
371 namespace internal {
372 
375 template<typename VectorValueType, typename SourceTreeType>
377 {
378  using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
379  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
380  using LeafT = typename SourceTreeType::LeafNodeType;
381  using TreeValueT = typename SourceTreeType::ValueType;
383 
384  const SourceTreeType* tree;
386 
387  CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
388 
389  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
390  {
391  VectorT& vec = *vector;
392  if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
393  // If a corresponding leaf node exists in the source tree,
394  // copy voxel values from the source node to the output vector.
395  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
396  vec[*it] = leaf->getValue(it.pos());
397  }
398  } else {
399  // If no corresponding leaf exists in the source tree,
400  // fill the vector with a uniform value.
401  const TreeValueT& value = tree->getValue(idxLeaf.origin());
402  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
403  vec[*it] = value;
404  }
405  }
406  }
407 };
408 
409 } // namespace internal
410 
411 
412 template<typename VectorValueType, typename SourceTreeType>
414 createVectorFromTree(const SourceTreeType& tree,
415  const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
416 {
417  using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
418  using VIdxLeafMgrT = tree::LeafManager<const VIdxTreeT>;
419  using VectorT = typename math::pcg::Vector<VectorValueType>;
420 
421  // Allocate the vector.
422  const size_t numVoxels = idxTree.activeVoxelCount();
423  typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
424 
425  // Parallelize over the leaf nodes of the index tree, filling the output vector
426  // with values from corresponding voxels of the source tree.
427  VIdxLeafMgrT leafManager(idxTree);
428  leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
429 
430  return result;
431 }
432 
433 
435 
436 
437 namespace internal {
438 
441 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
443 {
444  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
445  using OutLeafT = typename OutTreeT::LeafNodeType;
446  using VIdxLeafT = typename VIndexTreeType::LeafNodeType;
448 
449  const VectorT* vector;
451 
452  CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
453 
454  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
455  {
456  const VectorT& vec = *vector;
457  OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
458  assert(leaf != nullptr);
459  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
460  leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
461  }
462  }
463 };
464 
465 } // namespace internal
466 
467 
468 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
469 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
472  const VIndexTreeType& idxTree,
473  const TreeValueType& background)
474 {
475  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
476  using VIdxLeafMgrT = typename tree::LeafManager<const VIndexTreeType>;
477 
478  // Construct an output tree with the same active voxel topology as the index tree.
479  typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
480  OutTreeT& tree = *result;
481 
482  // Parallelize over the leaf nodes of the index tree, populating voxels
483  // of the output tree with values from the input vector.
484  VIdxLeafMgrT leafManager(idxTree);
485  leafManager.foreach(
487 
488  return result;
489 }
490 
491 
493 
494 
495 namespace internal {
496 
498 template<typename BoolTreeType, typename BoundaryOp>
500 {
501  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
502  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
505 
508  const BoolTreeType* interiorMask;
509  const BoundaryOp boundaryOp;
511 
513  const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
514  laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
515 
516  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
517  {
518  // Local accessors
519  typename tree::ValueAccessor<const BoolTreeType> interior(*interiorMask);
520  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
521 
522  Coord ijk;
523  VIndex column;
524  const ValueT diagonal = -6.f, offDiagonal = 1.f;
525 
526  // Loop over active voxels in this leaf.
527  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
528  assert(it.getValue() > -1);
529  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
530 
531  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
532 
533  ijk = it.getCoord();
534  if (interior.isValueOn(ijk)) {
535  // The current voxel is an interior voxel.
536  // All of its neighbors are in the solution domain.
537 
538  // -x direction
539  row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
540  // -y direction
541  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
542  // -z direction
543  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
544  // diagonal
545  row.setValue(rowNum, diagonal);
546  // +z direction
547  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
548  // +y direction
549  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
550  // +x direction
551  row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
552 
553  } else {
554  // The current voxel is a boundary voxel.
555  // At least one of its neighbors is outside the solution domain.
556 
557  ValueT modifiedDiagonal = 0.f;
558 
559  // -x direction
560  if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
561  row.setValue(column, offDiagonal);
562  modifiedDiagonal -= 1;
563  } else {
564  boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
565  }
566  // -y direction
567  if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
568  row.setValue(column, offDiagonal);
569  modifiedDiagonal -= 1;
570  } else {
571  boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
572  }
573  // -z direction
574  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
575  row.setValue(column, offDiagonal);
576  modifiedDiagonal -= 1;
577  } else {
578  boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
579  }
580  // +z direction
581  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
582  row.setValue(column, offDiagonal);
583  modifiedDiagonal -= 1;
584  } else {
585  boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
586  }
587  // +y direction
588  if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
589  row.setValue(column, offDiagonal);
590  modifiedDiagonal -= 1;
591  } else {
592  boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
593  }
594  // +x direction
595  if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
596  row.setValue(column, offDiagonal);
597  modifiedDiagonal -= 1;
598  } else {
599  boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
600  }
601  // diagonal
602  row.setValue(rowNum, modifiedDiagonal);
603  }
604  } // end loop over voxels
605  }
606 };
607 
608 
609 // Stencil 1 is the correct stencil, but stencil 2 requires
610 // half as many comparisons and produces smoother results at boundaries.
611 //#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
612 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
613 
615 template<typename VIdxTreeT, typename BoundaryOp>
617 {
618  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
621 
623  const VIdxTreeT* idxTree;
624  const BoundaryOp boundaryOp;
626 
627  ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
628  laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
629 
630  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
631  {
632  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
633 
634  const int kNumOffsets = 6;
635  const Coord ijkOffset[kNumOffsets] = {
636 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
637  Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
638 #else
639  Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
640 #endif
641  };
642 
643  // For each active voxel in this leaf...
644  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
645  assert(it.getValue() > -1);
646 
647  const Coord ijk = it.getCoord();
648  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
649 
650  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
651 
652  ValueT modifiedDiagonal = 0.f;
653 
654  // For each of the neighbors of the voxel at (i,j,k)...
655  for (int dir = 0; dir < kNumOffsets; ++dir) {
656  const Coord neighbor = ijk + ijkOffset[dir];
657  VIndex column;
658  // For collocated vector grids, the central differencing stencil requires
659  // access to neighbors at a distance of two voxels in each direction
660  // (-x, +x, -y, +y, -z, +z).
661 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
662  const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
663  && vectorIdx.isValueOn(neighbor));
664 #else
665  const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
666 #endif
667  if (ijkIsInterior) {
668  // If (i,j,k) is sufficiently far away from the exterior,
669  // set its weight to one and adjust the center weight accordingly.
670  row.setValue(column, 1.f);
671  modifiedDiagonal -= 1.f;
672  } else {
673  // If (i,j,k) is adjacent to or one voxel away from the exterior,
674  // invoke the boundary condition functor.
675  boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
676  }
677  }
678  // Set the (possibly modified) weight for the voxel at (i,j,k).
679  row.setValue(rowNum, modifiedDiagonal);
680  }
681  }
682 };
683 
684 } // namespace internal
685 
686 
687 template<typename BoolTreeType>
689 createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
690  const BoolTreeType& interiorMask, bool staggered)
691 {
692  using ValueT = LaplacianMatrix::ValueType;
694  static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
696  return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
697 }
698 
699 
700 template<typename BoolTreeType, typename BoundaryOp>
703  const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
704  const BoolTreeType& interiorMask,
705  const BoundaryOp& boundaryOp,
707  bool staggered)
708 {
709  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
710  using VIdxLeafMgrT = typename tree::LeafManager<const VIdxTreeT>;
711 
712  // The number of active voxels is the number of degrees of freedom.
713  const Index64 numDoF = idxTree.activeVoxelCount();
714 
715  // Construct the matrix.
716  LaplacianMatrix::Ptr laplacianPtr(
717  new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
718  LaplacianMatrix& laplacian = *laplacianPtr;
719 
720  // Populate the matrix using a second-order, 7-point CD stencil.
721  VIdxLeafMgrT idxLeafManager(idxTree);
722  if (staggered) {
724  laplacian, idxTree, interiorMask, boundaryOp, source));
725  } else {
726  idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
727  laplacian, idxTree, boundaryOp, source));
728  }
729 
730  return laplacianPtr;
731 }
732 
733 
735 
736 
737 template<typename TreeType>
738 inline typename TreeType::Ptr
739 solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
740 {
741  util::NullInterrupter interrupter;
742  return solve(inTree, state, interrupter, staggered);
743 }
744 
745 
746 template<typename TreeType, typename Interrupter>
747 inline typename TreeType::Ptr
748 solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
749 {
751  return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
752 }
753 
754 
755 template<typename TreeType, typename BoundaryOp, typename Interrupter>
756 inline typename TreeType::Ptr
757 solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
758  math::pcg::State& state, Interrupter& interrupter, bool staggered)
759 {
761  return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
762  inTree, boundaryOp, state, interrupter, staggered);
763 }
764 
765 
766 template<
767  typename PreconditionerType,
768  typename TreeType,
769  typename BoundaryOp,
770  typename Interrupter>
771 inline typename TreeType::Ptr
773  const TreeType& inTree,
774  const BoundaryOp& boundaryOp,
775  math::pcg::State& state,
776  Interrupter& interrupter,
777  bool staggered)
778 {
779  return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
780  /*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
781 }
782 
783 template<
784  typename PreconditionerType,
785  typename TreeType,
786  typename DomainTreeType,
787  typename BoundaryOp,
788  typename Interrupter>
789 inline typename TreeType::Ptr
791  const TreeType& inTree,
792  const DomainTreeType& domainMask,
793  const BoundaryOp& boundaryOp,
794  math::pcg::State& state,
795  Interrupter& interrupter,
796  bool staggered)
797 {
798  using TreeValueT = typename TreeType::ValueType;
799  using VecValueT = LaplacianMatrix::ValueType;
800  using VectorT = typename math::pcg::Vector<VecValueT>;
801  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
802  using MaskTreeT = typename TreeType::template ValueConverter<bool>::Type;
803 
804  // 1. Create a mapping from active voxels of the input tree to elements of a vector.
805  typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
806 
807  // 2. Populate a vector with values from the input tree.
808  typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
809 
810  // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
812  typename MaskTreeT::Ptr interiorMask(
813  new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
814  tools::erodeVoxels(*interiorMask, /*iterations=*/1, tools::NN_FACE);
815 
816  // 4. Create the Laplacian matrix.
818  *idxTree, *interiorMask, boundaryOp, *b, staggered);
819 
820  // 5. Solve the Poisson equation.
821  laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
822  b->scale(-1.0);
823  typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
825  new PreconditionerType(*laplacian));
826  if (!precond->isValid()) {
827  precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
828  }
829 
830  state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
831 
832  // 6. Populate the output tree with values from the solution vector.
834  return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
835 }
836 
837 } // namespace poisson
838 } // namespace tools
839 } // namespace OPENVDB_VERSION_NAME
840 } // namespace openvdb
841 
842 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:237
CopyFromVecOp(const VectorT &v, OutTreeT &t)
Definition: PoissonSolver.h:452
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:630
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero...
Definition: PoissonSolver.h:317
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:414
Definition: ValueAccessor.h:193
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
Definition: PoissonSolver.h:82
LeafCountOp(VIndex *count_)
Definition: PoissonSolver.h:290
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:846
Functor for use with LeafManager::foreach() to populate a tree with values from a vector...
Definition: PoissonSolver.h:442
Dirichlet boundary condition functor.
Definition: PoissonSolver.h:268
Diagonal preconditioner.
Definition: ConjGradient.h:42
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
Definition: PoissonSolver.h:499
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
Definition: PoissonSolver.h:269
typename SourceTreeType::ValueType TreeValueT
Definition: PoissonSolver.h:381
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:39
Read/write accessor to a row of this matrix.
Definition: ConjGradient.h:419
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:516
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:230
typename math::pcg::Vector< ValueT > VectorT
Definition: PoissonSolver.h:620
ISLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoundaryOp &op, VectorT &src)
Definition: PoissonSolver.h:627
LaplacianMatrix * laplacian
Definition: PoissonSolver.h:506
Functor for use with LeafManager::foreach() to populate an array with per-leaf active voxel counts...
Definition: PoissonSolver.h:287
typename VIndexTreeType::template ValueConverter< TreeValueType >::Type OutTreeT
Definition: PoissonSolver.h:444
ISStaggeredLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoolTreeType &mask, const BoundaryOp &op, VectorT &src)
Definition: PoissonSolver.h:512
Index32 SizeType
Definition: ConjGradient.h:33
LeafIndexOp(const VIndex *count_)
Definition: PoissonSolver.h:303
Functor for use with LeafManager::foreach() to populate active leaf voxels with sequential indices...
Definition: PoissonSolver.h:300
Int32 VIndex
Definition: PoissonSolver.h:79
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 in...
Definition: PoissonSolver.h:748
VIndex * count
Definition: PoissonSolver.h:289
VectorT * vector
Definition: PoissonSolver.h:385
const BoundaryOp boundaryOp
Definition: PoissonSolver.h:509
typename VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:502
CopyToVecOp(const SourceTreeType &t, VectorT &v)
Definition: PoissonSolver.h:387
const VIdxTreeT * idxTree
Definition: PoissonSolver.h:507
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
typename SourceTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
Definition: PoissonSolver.h:378
typename VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:618
int32_t Int32
Definition: Types.h:33
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 valu...
Definition: PoissonSolver.h:790
LaplacianMatrix::ValueType ValueT
Definition: PoissonSolver.h:503
typename OutTreeT::LeafNodeType OutLeafT
Definition: PoissonSolver.h:445
const SourceTreeType * tree
Definition: PoissonSolver.h:384
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:43
typename VIndexTreeType::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:446
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 dif...
Definition: PoissonSolver.h:689
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:82
typename math::pcg::Vector< VectorValueType > VectorT
Definition: PoissonSolver.h:382
void operator()(const LeafType &leaf, size_t leafIdx) const
Definition: PoissonSolver.h:291
typename math::pcg::Vector< VectorValueType > VectorT
Definition: PoissonSolver.h:447
const VectorT * vector
Definition: PoissonSolver.h:449
ValueType_ ValueType
Definition: ConjGradient.h:240
const VIndex * count
Definition: PoissonSolver.h:302
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:681
const VIdxTreeT * idxTree
Definition: PoissonSolver.h:623
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
Definition: PoissonSolver.h:616
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:1013
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:389
typename math::pcg::Vector< ValueT > VectorT
Definition: PoissonSolver.h:504
typename VIdxTreeT::LeafNodeType VIdxLeafT
Definition: PoissonSolver.h:379
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition: ConjGradient.h:1225
const BoolTreeType * interiorMask
Definition: PoissonSolver.h:508
typename SourceTreeType::LeafNodeType LeafT
Definition: PoissonSolver.h:380
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
SharedPtr< SparseStencilMatrix > Ptr
Definition: ConjGradient.h:242
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition: ConjGradient.h:1069
Lightweight, variable-length vector.
Definition: ConjGradient.h:37
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:702
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:470
GridType::template ValueConverter< bool >::Type::Ptr interiorMask(const GridType &grid, const double isovalue=0.0)
Given an input grid of any type, return a new, boolean grid whose active voxel topology matches the i...
Definition: Mask.h:111
OutTreeT * tree
Definition: PoissonSolver.h:450
uint64_t Index64
Definition: Types.h:30
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:240
SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:468
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
LaplacianMatrix * laplacian
Definition: PoissonSolver.h:622
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
const BoundaryOp boundaryOp
Definition: PoissonSolver.h:624
VectorT * source
Definition: PoissonSolver.h:625
typename BoolTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
Definition: PoissonSolver.h:501
void operator()(LeafType &leaf, size_t leafIdx) const
Definition: PoissonSolver.h:304
Definition: Morphology.h:60
Functor for use with LeafManager::foreach() to populate a vector with the values of a tree&#39;s active v...
Definition: PoissonSolver.h:376
LaplacianMatrix::ValueType ValueT
Definition: PoissonSolver.h:619
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:142
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Definition: PoissonSolver.h:454
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 valu...
Definition: PoissonSolver.h:757
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:350