61 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
98 template<
typename TreeType>
99 inline typename TreeType::Ptr
102 template<
typename TreeType,
typename Interrupter>
103 inline typename TreeType::Ptr
139 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
140 inline typename TreeType::Ptr
146 bool staggered =
false);
149 typename PreconditionerType,
152 typename Interrupter>
153 inline typename TreeType::Ptr
159 bool staggered =
false);
162 typename PreconditionerType,
164 typename DomainTreeType,
166 typename Interrupter>
167 inline typename TreeType::Ptr
170 const DomainTreeType&,
174 bool staggered =
false);
184 template<
typename VIndexTreeType>
190 template<
typename TreeType>
191 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
201 template<
typename VectorValueType,
typename SourceTreeType>
204 const SourceTreeType& source,
205 const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
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);
227 template<
typename BoolTreeType>
230 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
232 bool staggered =
false);
254 template<
typename BoolTreeType,
typename BoundaryOp>
257 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
259 const BoundaryOp& boundaryOp,
261 bool staggered =
false);
267 template<
typename ValueType>
287 template<
typename LeafType>
291 LeafCountOp(
VIndex* count_): count(count_) {}
292 void operator()(
const LeafType& leaf,
size_t leafIdx)
const {
293 count[leafIdx] =
static_cast<VIndex>(leaf.onVoxelCount());
300 template<
typename LeafType>
304 LeafIndexOp(
const VIndex* count_): count(count_) {}
305 void operator()(LeafType& leaf,
size_t leafIdx)
const {
306 VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
307 for (
typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
317 template<
typename VIndexTreeType>
321 using LeafT =
typename VIndexTreeType::LeafNodeType;
325 LeafMgrT leafManager(result);
326 const size_t leafCount = leafManager.leafCount();
328 if (leafCount == 0)
return;
331 std::unique_ptr<VIndex[]> perLeafCount(
new VIndex[leafCount]);
332 VIndex* perLeafCountPtr = perLeafCount.get();
333 leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
337 for (
size_t i = 1; i < leafCount; ++i) {
338 perLeafCount[i] += perLeafCount[i - 1];
342 assert(
Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
346 leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
350 template<
typename TreeType>
351 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
354 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
357 const VIndex invalidIdx = -1;
358 typename VIdxTreeT::Ptr result(
362 result->voxelizeActiveTiles();
378 template<
typename VectorValueType,
typename SourceTreeType>
381 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
382 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
383 using LeafT =
typename SourceTreeType::LeafNodeType;
384 using TreeValueT =
typename SourceTreeType::ValueType;
387 const SourceTreeType* tree;
390 CopyToVecOp(
const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
392 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
394 VectorT& vec = *vector;
395 if (
const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
398 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
399 vec[*it] = leaf->getValue(it.pos());
404 const TreeValueT& value = tree->getValue(idxLeaf.origin());
405 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
416 template<
typename VectorValueType,
typename SourceTreeType>
417 inline typename math::pcg::Vector<VectorValueType>::Ptr
419 const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
421 using VIdxTreeT =
typename SourceTreeType::template ValueConverter<VIndex>::Type;
426 const size_t numVoxels = idxTree.activeVoxelCount();
431 VIdxLeafMgrT leafManager(idxTree);
432 leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
446 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
449 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
450 using OutLeafT =
typename OutTreeT::LeafNodeType;
451 using VIdxLeafT =
typename VIndexTreeType::LeafNodeType;
454 const VectorT* vector;
457 CopyFromVecOp(
const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
459 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
461 const VectorT& vec = *vector;
462 OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
463 assert(leaf !=
nullptr);
464 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
465 leaf->setValueOnly(it.pos(),
static_cast<TreeValueType
>(vec[*it]));
474 template<
typename TreeValueType,
typename VIndexTreeType,
typename VectorValueType>
475 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
478 const VIndexTreeType& idxTree,
479 const TreeValueType& background)
481 using OutTreeT =
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
485 typename OutTreeT::Ptr result(
new OutTreeT(idxTree, background,
TopologyCopy()));
486 OutTreeT& tree = *result;
490 VIdxLeafMgrT leafManager(idxTree);
492 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
505 template<
typename BoolTreeType,
typename BoundaryOp>
506 struct ISStaggeredLaplacianOp
508 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
509 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
510 using ValueT = LaplacianMatrix::ValueType;
514 const VIdxTreeT* idxTree;
516 const BoundaryOp boundaryOp;
520 const BoolTreeType& mask,
const BoundaryOp& op, VectorT& src):
523 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
531 const ValueT diagonal = -6.f, offDiagonal = 1.f;
534 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
535 assert(it.getValue() > -1);
538 LaplacianMatrix::RowEditor row =
laplacian->getRowEditor(rowNum);
541 if (interior.isValueOn(ijk)) {
546 row.setValue(vectorIdx.getValue(ijk.
offsetBy(-1, 0, 0)), offDiagonal);
548 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, -1, 0)), offDiagonal);
550 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 0, -1)), offDiagonal);
552 row.setValue(rowNum, diagonal);
554 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 0, 1)), offDiagonal);
556 row.setValue(vectorIdx.getValue(ijk.
offsetBy(0, 1, 0)), offDiagonal);
558 row.setValue(vectorIdx.getValue(ijk.
offsetBy(1, 0, 0)), offDiagonal);
564 ValueT modifiedDiagonal = 0.f;
567 if (vectorIdx.probeValue(ijk.
offsetBy(-1, 0, 0), column)) {
568 row.setValue(column, offDiagonal);
569 modifiedDiagonal -= 1;
571 boundaryOp(ijk, ijk.
offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
574 if (vectorIdx.probeValue(ijk.
offsetBy(0, -1, 0), column)) {
575 row.setValue(column, offDiagonal);
576 modifiedDiagonal -= 1;
578 boundaryOp(ijk, ijk.
offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
581 if (vectorIdx.probeValue(ijk.
offsetBy(0, 0, -1), column)) {
582 row.setValue(column, offDiagonal);
583 modifiedDiagonal -= 1;
585 boundaryOp(ijk, ijk.
offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
588 if (vectorIdx.probeValue(ijk.
offsetBy(0, 0, 1), column)) {
589 row.setValue(column, offDiagonal);
590 modifiedDiagonal -= 1;
592 boundaryOp(ijk, ijk.
offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
595 if (vectorIdx.probeValue(ijk.
offsetBy(0, 1, 0), column)) {
596 row.setValue(column, offDiagonal);
597 modifiedDiagonal -= 1;
599 boundaryOp(ijk, ijk.
offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
602 if (vectorIdx.probeValue(ijk.
offsetBy(1, 0, 0), column)) {
603 row.setValue(column, offDiagonal);
604 modifiedDiagonal -= 1;
606 boundaryOp(ijk, ijk.
offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
609 row.setValue(rowNum, modifiedDiagonal);
619 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
622 template<
typename VIdxTreeT,
typename BoundaryOp>
625 using VIdxLeafT =
typename VIdxTreeT::LeafNodeType;
626 using ValueT = LaplacianMatrix::ValueType;
627 using VectorT =
typename math::pcg::Vector<ValueT>;
630 const VIdxTreeT* idxTree;
631 const BoundaryOp boundaryOp;
634 ISLaplacianOp(
LaplacianMatrix& m,
const VIdxTreeT& idx,
const BoundaryOp& op, VectorT& src):
635 laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
637 void operator()(
const VIdxLeafT& idxLeaf,
size_t )
const
639 typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
641 const int kNumOffsets = 6;
642 const Coord ijkOffset[kNumOffsets] = {
643 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
644 Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
646 Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
651 for (
typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
652 assert(it.getValue() > -1);
654 const Coord ijk = it.getCoord();
657 LaplacianMatrix::RowEditor row =
laplacian->getRowEditor(rowNum);
659 ValueT modifiedDiagonal = 0.f;
662 for (
int dir = 0; dir < kNumOffsets; ++dir) {
663 const Coord neighbor = ijk + ijkOffset[dir];
668 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
669 const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
670 && vectorIdx.isValueOn(neighbor));
672 const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
677 row.setValue(column, 1.f);
678 modifiedDiagonal -= 1.f;
682 boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
686 row.setValue(rowNum, modifiedDiagonal);
696 template<
typename BoolTreeType>
697 inline LaplacianMatrix::Ptr
709 template<
typename BoolTreeType,
typename BoundaryOp>
710 inline LaplacianMatrix::Ptr
712 const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
714 const BoundaryOp& boundaryOp,
718 using VIdxTreeT =
typename BoolTreeType::template ValueConverter<VIndex>::Type;
722 const Index64 numDoF = idxTree.activeVoxelCount();
730 VIdxLeafMgrT idxLeafManager(idxTree);
732 idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
735 idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
736 laplacian, idxTree, boundaryOp, source));
746 template<
typename TreeType>
747 inline typename TreeType::Ptr
751 return solve(inTree, state, interrupter, staggered);
755 template<
typename TreeType,
typename Interrupter>
756 inline typename TreeType::Ptr
764 template<
typename TreeType,
typename BoundaryOp,
typename Interrupter>
765 inline typename TreeType::Ptr
770 return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
771 inTree, boundaryOp, state, interrupter, staggered);
776 typename PreconditionerType,
779 typename Interrupter>
780 inline typename TreeType::Ptr
782 const TreeType& inTree,
783 const BoundaryOp& boundaryOp,
785 Interrupter& interrupter,
788 return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
789 inTree, inTree, boundaryOp, state, interrupter, staggered);
793 typename PreconditionerType,
795 typename DomainTreeType,
797 typename Interrupter>
798 inline typename TreeType::Ptr
800 const TreeType& inTree,
801 const DomainTreeType& domainMask,
802 const BoundaryOp& boundaryOp,
804 Interrupter& interrupter,
807 using TreeValueT =
typename TreeType::ValueType;
810 using VIdxTreeT =
typename TreeType::template ValueConverter<VIndex>::Type;
811 using MaskTreeT =
typename TreeType::template ValueConverter<bool>::Type;
817 typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
832 typename VectorT::Ptr x(
new VectorT(b->size(), zeroVal<VecValueT>()));
843 return createTreeFromVector<TreeValueT>(*x, *idxTree, zeroVal<TreeValueT>());
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
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:564
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:183
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
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:180