61 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
99 template<
typename TreeType>
100 typename TreeType::Ptr
103 template<
typename TreeType,
typename Interrupter>
104 typename TreeType::Ptr
140 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
141 typename TreeType::Ptr
147 bool staggered =
false);
150 typename PreconditionerType,
153 typename Interrupter>
154 typename TreeType::Ptr
160 bool staggered =
false);
163 typename PreconditionerType,
165 typename DomainTreeType,
167 typename Interrupter>
168 typename TreeType::Ptr
171 const DomainTreeType&,
175 bool staggered =
false);
185 template<
typename VIndexTreeType>
191 template<
typename TreeType>
192 typename TreeType::template ValueConverter<VIndex>::Type::Ptr
202 template<
typename VectorValueType,
typename SourceTreeType>
205 const SourceTreeType&
source,
216 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
217 typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
220 const VIndexTreeType&
index,
221 const TreeValueType& background);
228 template<
typename BoolTreeType>
233 bool staggered =
false);
255 template<
typename BoolTreeType,
typename BoundaryOp>
260 const BoundaryOp& boundaryOp,
262 bool staggered =
false);
268 template<
typename ValueType>
288 template<
typename LeafType>
293 void operator()(
const LeafType& leaf,
size_t leafIdx)
const {
294 count[leafIdx] =
static_cast<VIndex>(leaf.onVoxelCount());
301 template<
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) {
318 template<
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));
351 template<
typename TreeType>
352 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
358 const VIndex invalidIdx = -1;
359 typename VIdxTreeT::Ptr
result(
363 result->voxelizeActiveTiles();
379 template<
typename VectorValueType,
typename SourceTreeType>
383 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
384 using LeafT =
typename SourceTreeType::LeafNodeType;
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) {
417 template<
typename VectorValueType,
typename SourceTreeType>
418 inline typename math::pcg::Vector<VectorValueType>::Ptr
427 const size_t numVoxels = idxTree.activeVoxelCount();
428 typename VectorT::Ptr
result(
new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
432 VIdxLeafMgrT leafManager(idxTree);
433 leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
447 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
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]));
475 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
476 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
479 const VIndexTreeType& idxTree,
480 const TreeValueType& background)
491 VIdxLeafMgrT leafManager(idxTree);
493 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
506 template<
typename BoolTreeType,
typename BoundaryOp>
507 struct ISStaggeredLaplacianOp
510 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
515 const VIdxTreeT* idxTree;
517 const BoundaryOp boundaryOp;
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
623 template<
typename VIdxTreeT,
typename BoundaryOp>
626 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
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):
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();
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);
697 template<
typename BoolTreeType>
704 static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
710 template<
typename BoolTreeType,
typename BoundaryOp>
715 const BoundaryOp& boundaryOp,
723 const Index64 numDoF = idxTree.activeVoxelCount();
731 VIdxLeafMgrT idxLeafManager(idxTree);
733 idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
734 laplacian, idxTree, interiorMask, boundaryOp, source));
736 idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
737 laplacian, idxTree, boundaryOp, source));
747 template<
typename TreeType>
748 typename TreeType::Ptr
752 return solve(inTree, state, interrupter, staggered);
756 template<
typename TreeType,
typename Interrupter>
757 typename TreeType::Ptr
765 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
766 typename TreeType::Ptr
771 return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
772 inTree, boundaryOp, state, interrupter, staggered);
777 typename PreconditionerType,
780 typename Interrupter>
781 typename 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>
799 typename TreeType::Ptr
801 const TreeType& inTree,
802 const DomainTreeType& domainMask,
803 const BoundaryOp& boundaryOp,
805 Interrupter& interrupter,
818 typename VectorT::Ptr
b = createVectorFromTree<VecValueT>(inTree, *idxTree);
828 *idxTree, *interiorMask, boundaryOp, *b, staggered);
831 laplacian->scale(-1.0);
833 typename VectorT::Ptr
x(
new VectorT(b->size(), zeroVal<VecValueT>()));
835 new PreconditionerType(*laplacian));
836 if (!precond->isValid()) {
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)
891 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
899 #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.
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.