61 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED 62 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED 101 template<
typename TreeType>
102 typename TreeType::Ptr
105 template<
typename TreeType,
typename Interrupter>
106 typename TreeType::Ptr
142 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
143 typename TreeType::Ptr
149 bool staggered =
false);
152 typename PreconditionerType,
155 typename Interrupter>
156 typename TreeType::Ptr
162 bool staggered =
false);
165 typename PreconditionerType,
167 typename DomainTreeType,
169 typename Interrupter>
170 typename TreeType::Ptr
173 const DomainTreeType&,
177 bool staggered =
false);
187 template<
typename VIndexTreeType>
193 template<
typename TreeType>
194 typename TreeType::template ValueConverter<VIndex>::Type::Ptr
204 template<
typename VectorValueType,
typename SourceTreeType>
207 const SourceTreeType& source,
208 const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
218 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
219 typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
222 const VIndexTreeType& index,
223 const TreeValueType& background);
230 template<
typename BoolTreeType>
233 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
235 bool staggered =
false);
257 template<
typename BoolTreeType,
typename BoundaryOp>
260 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
262 const BoundaryOp& boundaryOp,
264 bool staggered =
false);
270 template<
typename ValueType>
272 inline void operator()(
const Coord&,
const Coord&, ValueType&, ValueType& diag)
const {
290 template<
typename LeafType>
294 LeafCountOp(
VIndex* count_): count(count_) {}
295 void operator()(
const LeafType& leaf,
size_t leafIdx)
const {
296 count[leafIdx] =
static_cast<VIndex>(leaf.onVoxelCount());
303 template<
typename LeafType>
307 LeafIndexOp(
const VIndex* count_): count(count_) {}
308 void operator()(LeafType& leaf,
size_t leafIdx)
const {
309 VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
310 for (
typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
320 template<
typename VIndexTreeType>
324 using LeafT =
typename VIndexTreeType::LeafNodeType;
328 LeafMgrT leafManager(result);
329 const size_t leafCount = leafManager.leafCount();
331 if (leafCount == 0)
return;
334 std::unique_ptr<VIndex[]> perLeafCount(
new VIndex[leafCount]);
335 VIndex* perLeafCountPtr = perLeafCount.get();
336 leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
340 for (
size_t i = 1; i < leafCount; ++i) {
341 perLeafCount[i] += perLeafCount[i - 1];
349 leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
353 template<
typename TreeType>
354 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
357 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
360 const VIndex invalidIdx = -1;
361 typename VIdxTreeT::Ptr result(
365 result->voxelizeActiveTiles();
381 template<
typename VectorValueType,
typename SourceTreeType>
384 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
385 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
386 using LeafT =
typename SourceTreeType::LeafNodeType;
387 using TreeValueT =
typename SourceTreeType::ValueType;
390 const SourceTreeType* tree;
393 CopyToVecOp(
const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
395 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const 397 VectorT& vec = *vector;
398 if (
const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
401 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
402 vec[*it] = leaf->getValue(it.pos());
407 const TreeValueT& value = tree->getValue(idxLeaf.origin());
408 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
419 template<
typename VectorValueType,
typename SourceTreeType>
422 const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
424 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
429 const size_t numVoxels = idxTree.activeVoxelCount();
430 typename VectorT::Ptr result(
new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
434 VIdxLeafMgrT leafManager(idxTree);
435 leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
449 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
452 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
453 using OutLeafT =
typename OutTreeT::LeafNodeType;
454 using VIdxLeafT =
typename VIndexTreeType::LeafNodeType;
457 const VectorT* vector;
460 CopyFromVecOp(
const VectorT& v,
OutTreeT& t): vector(&v), tree(&t) {}
462 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const 464 const VectorT& vec = *vector;
465 OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
467 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
468 leaf->setValueOnly(it.pos(),
static_cast<TreeValueType
>(vec[*it]));
477 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
478 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
481 const VIndexTreeType& idxTree,
482 const TreeValueType& background)
484 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
493 VIdxLeafMgrT leafManager(idxTree);
495 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
508 template<
typename BoolTreeType,
typename BoundaryOp>
509 struct ISStaggeredLaplacianOp
511 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
512 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
517 const VIdxTreeT* idxTree;
519 const BoundaryOp boundaryOp;
523 const BoolTreeType& mask,
const BoundaryOp&
op, VectorT& src):
526 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const 534 const ValueT diagonal = -6.f, offDiagonal = 1.f;
537 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
567 ValueT modifiedDiagonal = 0.f;
570 if (vectorIdx.
probeValue(ijk.offsetBy(-1, 0, 0), column)) {
572 modifiedDiagonal -= 1;
574 boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
577 if (vectorIdx.
probeValue(ijk.offsetBy(0, -1, 0), column)) {
579 modifiedDiagonal -= 1;
581 boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
584 if (vectorIdx.
probeValue(ijk.offsetBy(0, 0, -1), column)) {
586 modifiedDiagonal -= 1;
588 boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
591 if (vectorIdx.
probeValue(ijk.offsetBy(0, 0, 1), column)) {
593 modifiedDiagonal -= 1;
595 boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
598 if (vectorIdx.
probeValue(ijk.offsetBy(0, 1, 0), column)) {
600 modifiedDiagonal -= 1;
602 boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
605 if (vectorIdx.
probeValue(ijk.offsetBy(1, 0, 0), column)) {
607 modifiedDiagonal -= 1;
609 boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
612 row.
setValue(rowNum, modifiedDiagonal);
622 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2 625 template<
typename VIdxTreeT,
typename BoundaryOp>
628 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
633 const VIdxTreeT* idxTree;
634 const BoundaryOp boundaryOp;
637 ISLaplacianOp(
LaplacianMatrix& m,
const VIdxTreeT& idx,
const BoundaryOp&
op, VectorT& src):
638 laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
640 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const 644 const int kNumOffsets = 6;
645 const Coord ijkOffset[kNumOffsets] = {
646 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1 647 Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
649 Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
654 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
657 const Coord ijk = it.getCoord();
662 ValueT modifiedDiagonal = 0.f;
665 for (
int dir = 0; dir < kNumOffsets; ++dir) {
666 const Coord neighbor = ijk + ijkOffset[dir];
671 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1 672 const bool ijkIsInterior = (vectorIdx.
probeValue(neighbor + ijkOffset[dir], column)
675 const bool ijkIsInterior = vectorIdx.
probeValue(neighbor, column);
681 modifiedDiagonal -= 1.f;
685 boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
689 row.
setValue(rowNum, modifiedDiagonal);
699 template<
typename BoolTreeType>
706 static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
712 template<
typename BoolTreeType,
typename BoundaryOp>
715 const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
717 const BoundaryOp& boundaryOp,
721 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
725 const Index64 numDoF = idxTree.activeVoxelCount();
733 VIdxLeafMgrT idxLeafManager(idxTree);
735 idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
736 laplacian, idxTree, interiorMask, boundaryOp, source));
738 idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
739 laplacian, idxTree, boundaryOp, source));
749 template<
typename TreeType>
750 typename TreeType::Ptr
754 return solve(inTree, state, interrupter, staggered);
758 template<
typename TreeType,
typename Interrupter>
759 typename TreeType::Ptr
767 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
768 typename TreeType::Ptr
773 return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
774 inTree, boundaryOp, state, interrupter, staggered);
779 typename PreconditionerType,
782 typename Interrupter>
783 typename TreeType::Ptr
785 const TreeType& inTree,
786 const BoundaryOp& boundaryOp,
788 Interrupter& interrupter,
791 return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
792 inTree, inTree, boundaryOp, state, interrupter, staggered);
796 typename PreconditionerType,
798 typename DomainTreeType,
800 typename Interrupter>
801 typename TreeType::Ptr
803 const TreeType& inTree,
804 const DomainTreeType& domainMask,
805 const BoundaryOp& boundaryOp,
807 Interrupter& interrupter,
810 using TreeValueT =
typename TreeType::ValueType;
813 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
814 using MaskTreeT =
typename TreeType::template ValueConverter<bool>::Type;
820 typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
830 *idxTree, *interiorMask, boundaryOp, *b, staggered);
833 laplacian->scale(-1.0);
835 typename VectorT::Ptr x(
new VectorT(b->size(), zeroVal<VecValueT>()));
837 new PreconditionerType(*laplacian));
838 if (!precond->isValid()) {
846 return createTreeFromVector<TreeValueT>(*x, *idxTree, zeroVal<TreeValueT>());
855 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 857 #ifdef OPENVDB_INSTANTIATE_POISSONSOLVER 861 #define _FUNCTION(TreeT) \ 862 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \ 863 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \ 864 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \ 865 math::pcg::State&, util::NullInterrupter&, bool) 869 #define _FUNCTION(TreeT) \ 870 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \ 871 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \ 872 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \ 873 math::pcg::State&, util::NullInterrupter&, bool) 877 #define _FUNCTION(TreeT) \ 878 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \ 879 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \ 880 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \ 881 math::pcg::State&, util::NullInterrupter&, bool) 885 #define _FUNCTION(TreeT) \ 886 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \ 887 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \ 888 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \ 889 math::pcg::State&, util::NullInterrupter&, bool) 893 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 901 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:468
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:142
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
Definition: ConjGradient.h:1065
uint64_t Index64
Definition: Types.h:53
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:39
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
Definition: ConjGradient.h:1221
Lightweight, variable-length vector.
Definition: ConjGradient.h:37
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
int32_t Int32
Definition: Types.h:56
Base class for interrupters.
Definition: NullInterrupter.h:25
Read/write accessor to a row of this matrix.
Definition: ConjGradient.h:419
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:43
SharedPtr< SparseStencilMatrix > Ptr
Definition: ConjGradient.h:242
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
Definition: ValueAccessor.h:492
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
Index32 SizeType
Definition: ConjGradient.h:33
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:455
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:162
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
ValueType_ ValueType
Definition: ConjGradient.h:240
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Diagonal preconditioner.
Definition: ConjGradient.h:42
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:480
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218