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,
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>
235 bool staggered =
false);
257 template<
typename BoolTreeType,
typename BoundaryOp>
262 const BoundaryOp& boundaryOp,
264 bool staggered =
false);
270 template<
typename ValueType>
290 template<
typename LeafType>
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
360 const VIndex invalidIdx = -1;
361 typename VIdxTreeT::Ptr
result(
365 result->voxelizeActiveTiles();
381 template<
typename VectorValueType,
typename SourceTreeType>
385 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
386 using LeafT =
typename SourceTreeType::LeafNodeType;
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>
420 inline typename math::pcg::Vector<VectorValueType>::Ptr
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>
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)
493 VIdxLeafMgrT leafManager(idxTree);
495 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
508 template<
typename BoolTreeType,
typename BoundaryOp>
509 struct ISStaggeredLaplacianOp
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) {
541 LaplacianMatrix::RowEditor
row =
laplacian->getRowEditor(rowNum);
544 if (interior.isValueOn(ijk)) {
549 row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
551 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
553 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
555 row.setValue(rowNum, diagonal);
557 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
559 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
561 row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
567 ValueT modifiedDiagonal = 0.f;
570 if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0),
column)) {
571 row.setValue(column, offDiagonal);
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)) {
578 row.setValue(column, offDiagonal);
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)) {
585 row.setValue(column, offDiagonal);
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)) {
592 row.setValue(column, offDiagonal);
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)) {
599 row.setValue(column, offDiagonal);
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)) {
606 row.setValue(column, offDiagonal);
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;
630 using VectorT =
typename math::pcg::Vector<ValueT>;
633 const VIdxTreeT* idxTree;
634 const BoundaryOp boundaryOp;
640 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
642 typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
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();
660 LaplacianMatrix::RowEditor row =
laplacian->getRowEditor(rowNum);
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)
673 && vectorIdx.isValueOn(neighbor));
675 const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
680 row.setValue(column, 1.
f);
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>
717 const BoundaryOp& boundaryOp,
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,
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
Lightweight, variable-length vector.
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
GLsizei const GLfloat * value
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
#define OPENVDB_USE_VERSION_NAMESPACE
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
**But if you need a result
Information about the state of a conjugate gradient solution.
Base class for interrupters.
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
#define OPENVDB_ASSERT(X)
SharedPtr< Preconditioner > Ptr
Preconditioner using incomplete Cholesky factorization.
State solve(const PositiveDefMatrix &A, const Vector< typename PositiveDefMatrix::ValueType > &b, Vector< typename PositiveDefMatrix::ValueType > &x, Preconditioner< typename PositiveDefMatrix::ValueType > &preconditioner, const State &termination=terminationDefaults< typename PositiveDefMatrix::ValueType >())
Solve Ax = b via the preconditioned conjugate gradient method.
GLsizei GLsizei GLchar * source
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
GLboolean GLboolean GLboolean b
SharedPtr< SparseStencilMatrix > Ptr
GLenum GLsizei GLsizei GLint * values
Implementation of morphological dilation and erosion.
GLenum GLenum GLsizei void * row
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
GLenum GLenum GLsizei void GLsizei void * column
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.