61#ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62#define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
99template<
typename TreeType>
100typename TreeType::Ptr
103template<
typename TreeType,
typename Interrupter>
104typename TreeType::Ptr
105solve(
const TreeType&,
math::pcg::State&, Interrupter&,
bool staggered =
false);
140template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
141typename TreeType::Ptr
147 bool staggered =
false);
150 typename PreconditionerType,
153 typename Interrupter>
154typename TreeType::Ptr
160 bool staggered =
false);
163 typename PreconditionerType,
165 typename DomainTreeType,
167 typename Interrupter>
168typename TreeType::Ptr
171 const DomainTreeType&,
175 bool staggered =
false);
185template<
typename VIndexTreeType>
191template<
typename TreeType>
192typename TreeType::template ValueConverter<VIndex>::Type::Ptr
202template<
typename VectorValueType,
typename SourceTreeType>
205 const SourceTreeType& source,
206 const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
216template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
217typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
220 const VIndexTreeType& index,
221 const TreeValueType& background);
228template<
typename BoolTreeType>
231 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
232 const BoolTreeType& interiorMask,
233 bool staggered =
false);
255template<
typename BoolTreeType,
typename BoundaryOp>
258 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
259 const BoolTreeType& interiorMask,
260 const BoundaryOp& boundaryOp,
262 bool staggered =
false);
268template<
typename ValueType>
288template<
typename LeafType>
292 LeafCountOp(VIndex* count_): count(count_) {}
293 void operator()(
const LeafType& leaf,
size_t leafIdx)
const {
294 count[leafIdx] =
static_cast<VIndex
>(leaf.onVoxelCount());
301template<
typename LeafType>
305 LeafIndexOp(
const VIndex* count_): count(count_) {}
306 void operator()(LeafType& leaf,
size_t leafIdx)
const {
307 VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
308 for (
typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
318template<
typename VIndexTreeType>
322 using LeafT =
typename VIndexTreeType::LeafNodeType;
326 LeafMgrT leafManager(result);
327 const size_t leafCount = leafManager.leafCount();
329 if (leafCount == 0)
return;
332 std::unique_ptr<VIndex[]> perLeafCount(
new VIndex[leafCount]);
333 VIndex* perLeafCountPtr = perLeafCount.get();
334 leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
338 for (
size_t i = 1; i < leafCount; ++i) {
339 perLeafCount[i] += perLeafCount[i - 1];
343 assert(
Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
347 leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
351template<
typename TreeType>
352inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
355 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
358 const VIndex invalidIdx = -1;
359 typename VIdxTreeT::Ptr result(
363 result->voxelizeActiveTiles();
379template<
typename VectorValueType,
typename SourceTreeType>
382 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
383 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
384 using LeafT =
typename SourceTreeType::LeafNodeType;
385 using TreeValueT =
typename SourceTreeType::ValueType;
388 const SourceTreeType* tree;
391 CopyToVecOp(
const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
393 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
395 VectorT& vec = *vector;
396 if (
const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
399 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
400 vec[*it] = leaf->getValue(it.pos());
405 const TreeValueT&
value = tree->getValue(idxLeaf.origin());
406 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
417template<
typename VectorValueType,
typename SourceTreeType>
418inline typename math::pcg::Vector<VectorValueType>::Ptr
420 const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
422 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
427 const size_t numVoxels = idxTree.activeVoxelCount();
432 VIdxLeafMgrT leafManager(idxTree);
433 leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
447template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
450 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
451 using OutLeafT =
typename OutTreeT::LeafNodeType;
452 using VIdxLeafT =
typename VIndexTreeType::LeafNodeType;
455 const VectorT* vector;
458 CopyFromVecOp(
const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
460 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
462 const VectorT& vec = *vector;
463 OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
464 assert(leaf !=
nullptr);
465 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
466 leaf->setValueOnly(it.pos(),
static_cast<TreeValueType
>(vec[*it]));
475template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
476inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
479 const VIndexTreeType& idxTree,
480 const TreeValueType& background)
482 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
486 typename OutTreeT::Ptr result(
new OutTreeT(idxTree, background,
TopologyCopy()));
487 OutTreeT& tree = *result;
491 VIdxLeafMgrT leafManager(idxTree);
493 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
506template<
typename BoolTreeType,
typename BoundaryOp>
507struct ISStaggeredLaplacianOp
509 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
510 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
511 using ValueT = LaplacianMatrix::ValueType;
515 const VIdxTreeT* idxTree;
517 const BoundaryOp boundaryOp;
520 ISStaggeredLaplacianOp(LaplacianMatrix& m,
const VIdxTreeT& idx,
521 const BoolTreeType& mask,
const BoundaryOp& op, VectorT& src):
524 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
532 const ValueT diagonal = -6.f, offDiagonal = 1.f;
535 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
536 assert(it.getValue() > -1);
539 LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
542 if (interior.isValueOn(ijk)) {
547 row.setValue(vectorIdx.getValue(ijk.
offsetBy(-1, 0, 0)), offDiagonal);
549 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, -1, 0)), offDiagonal);
551 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 0, -1)), offDiagonal);
553 row.setValue(rowNum, diagonal);
555 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 0, 1)), offDiagonal);
557 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 1, 0)), offDiagonal);
559 row.setValue(vectorIdx.getValue(ijk.
offsetBy(1, 0, 0)), offDiagonal);
565 ValueT modifiedDiagonal = 0.f;
568 if (vectorIdx.probeValue(ijk.
offsetBy(-1, 0, 0), column)) {
569 row.setValue(column, offDiagonal);
570 modifiedDiagonal -= 1;
572 boundaryOp(ijk, ijk.
offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
575 if (vectorIdx.probeValue(ijk.
offsetBy(0, -1, 0), column)) {
576 row.setValue(column, offDiagonal);
577 modifiedDiagonal -= 1;
579 boundaryOp(ijk, ijk.
offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
582 if (vectorIdx.probeValue(ijk.
offsetBy(0, 0, -1), column)) {
583 row.setValue(column, offDiagonal);
584 modifiedDiagonal -= 1;
586 boundaryOp(ijk, ijk.
offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
589 if (vectorIdx.probeValue(ijk.
offsetBy(0, 0, 1), column)) {
590 row.setValue(column, offDiagonal);
591 modifiedDiagonal -= 1;
593 boundaryOp(ijk, ijk.
offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
596 if (vectorIdx.probeValue(ijk.
offsetBy(0, 1, 0), column)) {
597 row.setValue(column, offDiagonal);
598 modifiedDiagonal -= 1;
600 boundaryOp(ijk, ijk.
offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
603 if (vectorIdx.probeValue(ijk.
offsetBy(1, 0, 0), column)) {
604 row.setValue(column, offDiagonal);
605 modifiedDiagonal -= 1;
607 boundaryOp(ijk, ijk.
offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
610 row.setValue(rowNum, modifiedDiagonal);
620#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
623template<
typename VIdxTreeT,
typename BoundaryOp>
626 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
627 using ValueT = LaplacianMatrix::ValueType;
628 using VectorT =
typename math::pcg::Vector<ValueT>;
631 const VIdxTreeT* idxTree;
632 const BoundaryOp boundaryOp;
635 ISLaplacianOp(LaplacianMatrix& m,
const VIdxTreeT& idx,
const BoundaryOp& op, VectorT& src):
636 laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
638 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
640 typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
642 const int kNumOffsets = 6;
643 const Coord ijkOffset[kNumOffsets] = {
644#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
645 Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
647 Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
652 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
653 assert(it.getValue() > -1);
655 const Coord ijk = it.getCoord();
656 const math::pcg::SizeType rowNum =
static_cast<math::pcg::SizeType
>(it.getValue());
658 LaplacianMatrix::RowEditor row =
laplacian->getRowEditor(rowNum);
660 ValueT modifiedDiagonal = 0.f;
663 for (
int dir = 0; dir < kNumOffsets; ++dir) {
664 const Coord neighbor = ijk + ijkOffset[dir];
669#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
670 const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
671 && vectorIdx.isValueOn(neighbor));
673 const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
678 row.setValue(column, 1.f);
679 modifiedDiagonal -= 1.f;
683 boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
687 row.setValue(rowNum, modifiedDiagonal);
697template<
typename BoolTreeType>
698inline LaplacianMatrix::Ptr
700 const BoolTreeType& interiorMask,
bool staggered)
710template<
typename BoolTreeType,
typename BoundaryOp>
711inline LaplacianMatrix::Ptr
713 const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
714 const BoolTreeType& interiorMask,
715 const BoundaryOp& boundaryOp,
719 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
723 const Index64 numDoF = idxTree.activeVoxelCount();
731 VIdxLeafMgrT idxLeafManager(idxTree);
733 idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
736 idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
737 laplacian, idxTree, boundaryOp, source));
747template<
typename TreeType>
748typename TreeType::Ptr
752 return solve(inTree, state, interrupter, staggered);
756template<
typename TreeType,
typename Interrupter>
757typename TreeType::Ptr
758solve(
const TreeType& inTree,
math::pcg::State& state, Interrupter& interrupter,
bool staggered)
765template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
766typename TreeType::Ptr
771 return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
772 inTree, boundaryOp, state, interrupter, staggered);
777 typename PreconditionerType,
780 typename Interrupter>
781typename TreeType::Ptr
783 const TreeType& inTree,
784 const BoundaryOp& boundaryOp,
786 Interrupter& interrupter,
789 return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
790 inTree, inTree, boundaryOp, state, interrupter, staggered);
794 typename PreconditionerType,
796 typename DomainTreeType,
798 typename Interrupter>
799typename TreeType::Ptr
801 const TreeType& inTree,
802 const DomainTreeType& domainMask,
803 const BoundaryOp& boundaryOp,
805 Interrupter& interrupter,
808 using TreeValueT =
typename TreeType::ValueType;
811 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
812 using MaskTreeT =
typename TreeType::template ValueConverter<bool>::Type;
818 typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
824 tools::erodeActiveValues(*
interiorMask, 1, tools::NN_FACE, tools::IGNORE_TILES);
833 typename VectorT::Ptr x(
new VectorT(b->size(), zeroVal<VecValueT>()));
840 state = math::pcg::solve(*
laplacian, *b, *x, *precond, interrupter, state);
844 return createTreeFromVector<TreeValueT>(*x, *idxTree, zeroVal<TreeValueT>());
853#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
855#ifdef OPENVDB_INSTANTIATE_POISSONSOLVER
859#define _FUNCTION(TreeT) \
860 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
861 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
862 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
863 math::pcg::State&, util::NullInterrupter&, bool)
867#define _FUNCTION(TreeT) \
868 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
869 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
870 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
871 math::pcg::State&, util::NullInterrupter&, bool)
875#define _FUNCTION(TreeT) \
876 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
877 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
878 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
879 math::pcg::State&, util::NullInterrupter&, bool)
883#define _FUNCTION(TreeT) \
884 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
885 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
886 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
887 math::pcg::State&, util::NullInterrupter&, bool)
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
ValueT value
Definition: GridBuilder.h:1290
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Implementation of morphological dilation and erosion.
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:91
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:1347
Diagonal preconditioner.
Definition: ConjGradient.h:1281
virtual bool isValid() const
Definition: ConjGradient.h:473
SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:468
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:238
ValueType_ ValueType
Definition: ConjGradient.h:240
SharedPtr< SparseStencilMatrix > Ptr
Definition: ConjGradient.h:242
Lightweight, variable-length vector.
Definition: ConjGradient.h:139
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:142
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Definition: ValueAccessor.h:191
Index32 SizeType
Definition: ConjGradient.h:33
int32_t Int32
Definition: Types.h:56
uint64_t Index64
Definition: Types.h:53
Definition: Exceptions.h:13
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
Base class for interrupters.
Definition: NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157