OpenVDB  12.0.0
PoissonSolver.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 //
4 /// @file PoissonSolver.h
5 ///
6 /// @authors D.J. Hill, Peter Cucka
7 ///
8 /// @brief Solve Poisson's equation &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i>
9 /// for <i>x</i>, where @e b is a vector comprising the values of all of the active voxels
10 /// in a grid.
11 ///
12 /// @par Example:
13 /// Solve for the pressure in a cubic tank of liquid, assuming uniform boundary conditions:
14 /// @code
15 /// FloatTree source(/*background=*/0.0f);
16 /// // Activate voxels to indicate that they contain liquid.
17 /// source.fill(CoordBBox(Coord(0, -10, 0), Coord(10, 0, 10)), /*value=*/0.0f);
18 ///
19 /// math::pcg::State state = math::pcg::terminationDefaults<float>();
20 /// FloatTree::Ptr solution = tools::poisson::solve(source, state);
21 /// @endcode
22 ///
23 /// @par Example:
24 /// Solve for the pressure, <i>P</i>, in a cubic tank of liquid that is open at the top.
25 /// Boundary conditions are <i>P</i>&nbsp;=&nbsp;0 at the top,
26 /// &part;<i>P</i>/&part;<i>y</i>&nbsp;=&nbsp;&minus;1 at the bottom
27 /// and &part;<i>P</i>/&part;<i>x</i>&nbsp;=&nbsp;0 at the sides:
28 /// <pre>
29 /// P = 0
30 /// +--------+ (N,0,N)
31 /// /| /|
32 /// (0,0,0) +--------+ |
33 /// | | | | dP/dx = 0
34 /// dP/dx = 0 | +------|-+
35 /// |/ |/
36 /// (0,-N,0) +--------+ (N,-N,0)
37 /// dP/dy = -1
38 /// </pre>
39 /// @code
40 /// const int N = 10;
41 /// DoubleTree source(/*background=*/0.0);
42 /// // Activate voxels to indicate that they contain liquid.
43 /// source.fill(CoordBBox(Coord(0, -N, 0), Coord(N, 0, N)), /*value=*/0.0);
44 ///
45 /// auto boundary = [](const openvdb::Coord& ijk, const openvdb::Coord& neighbor,
46 /// double& source, double& diagonal)
47 /// {
48 /// if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) {
49 /// if (neighbor.y() < ijk.y()) source -= 1.0;
50 /// else diagonal -= 1.0;
51 /// }
52 /// };
53 ///
54 /// math::pcg::State state = math::pcg::terminationDefaults<double>();
55 /// util::NullInterrupter interrupter;
56 ///
57 /// DoubleTree::Ptr solution = tools::poisson::solveWithBoundaryConditions(
58 /// source, boundary, state, interrupter);
59 /// @endcode
60 
61 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
63 
64 #include <openvdb/Types.h>
67 #include <openvdb/tree/Tree.h>
69 #include <openvdb/util/Assert.h>
70 
71 #include "Morphology.h" // for erodeActiveValues
72 #include <openvdb/openvdb.h>
73 
74 
75 namespace openvdb {
77 namespace OPENVDB_VERSION_NAME {
78 namespace tools {
79 namespace poisson {
80 
81 // This type should be at least as wide as math::pcg::SizeType.
82 using VIndex = Int32;
83 
84 /// The type of a matrix used to represent a three-dimensional %Laplacian operator
86 
87 
88 //@{
89 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>,
90 /// where @e b is a vector comprising the values of all of the active voxels
91 /// in the input tree.
92 /// @return a new tree, with the same active voxel topology as the input tree,
93 /// whose voxel values are the elements of the solution vector <i>x</i>.
94 /// @details On input, the State object should specify convergence criteria
95 /// (minimum error and maximum number of iterations); on output, it gives
96 /// the actual termination conditions.
97 /// @details The solution is computed using the conjugate gradient method
98 /// with (where possible) incomplete Cholesky preconditioning, falling back
99 /// to Jacobi preconditioning.
100 /// @sa solveWithBoundaryConditions
101 template<typename TreeType>
102 typename TreeType::Ptr
103 solve(const TreeType&, math::pcg::State&, bool staggered = false);
104 
105 template<typename TreeType, typename Interrupter>
106 typename TreeType::Ptr
107 solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
108 //@}
109 
110 
111 //@{
112 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>
113 /// with user-specified boundary conditions, where @e b is a vector comprising
114 /// the values of all of the active voxels in the input tree or domain mask if provided
115 /// @return a new tree, with the same active voxel topology as the input tree,
116 /// whose voxel values are the elements of the solution vector <i>x</i>.
117 /// @details On input, the State object should specify convergence criteria
118 /// (minimum error and maximum number of iterations); on output, it gives
119 /// the actual termination conditions.
120 /// @details The solution is computed using the conjugate gradient method with
121 /// the specified type of preconditioner (default: incomplete Cholesky),
122 /// falling back to Jacobi preconditioning if necessary.
123 /// @details Each thread gets its own copy of the BoundaryOp, which should be
124 /// a functor of the form
125 /// @code
126 /// struct BoundaryOp {
127 /// using ValueType = LaplacianMatrix::ValueType;
128 /// void operator()(
129 /// const Coord& ijk, // coordinates of a boundary voxel
130 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
131 /// ValueType& source, // element of b corresponding to ijk
132 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
133 /// ) const;
134 /// };
135 /// @endcode
136 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
137 /// and it must specify a boundary condition for @ijk by modifying one or both of two
138 /// provided values: the entry in the source vector @e b corresponding to @ijk and
139 /// the weighting coefficient for @ijk in the Laplacian operator matrix.
140 ///
141 /// @sa solve
142 template<typename TreeType, typename BoundaryOp, typename Interrupter>
143 typename TreeType::Ptr
145  const TreeType&,
146  const BoundaryOp&,
148  Interrupter&,
149  bool staggered = false);
150 
151 template<
152  typename PreconditionerType,
153  typename TreeType,
154  typename BoundaryOp,
155  typename Interrupter>
156 typename TreeType::Ptr
158  const TreeType&,
159  const BoundaryOp&,
161  Interrupter&,
162  bool staggered = false);
163 
164 template<
165  typename PreconditionerType,
166  typename TreeType,
167  typename DomainTreeType,
168  typename BoundaryOp,
169  typename Interrupter>
170 typename TreeType::Ptr
172  const TreeType&,
173  const DomainTreeType&,
174  const BoundaryOp&,
176  Interrupter&,
177  bool staggered = false);
178 //@}
179 
180 
181 /// @name Low-level functions
182 //@{
183 // The following are low-level routines that can be used to assemble custom solvers.
184 
185 /// @brief Overwrite each active voxel in the given scalar tree
186 /// with a sequential index, starting from zero.
187 template<typename VIndexTreeType>
188 void populateIndexTree(VIndexTreeType&);
189 
190 /// @brief Iterate over the active voxels of the input tree and for each one
191 /// assign its index in the iteration sequence to the corresponding voxel
192 /// of an integer-valued output tree.
193 template<typename TreeType>
194 typename TreeType::template ValueConverter<VIndex>::Type::Ptr
195 createIndexTree(const TreeType&);
196 
197 
198 /// @brief Return a vector of the active voxel values of the scalar-valued @a source tree.
199 /// @details The <i>n</i>th element of the vector corresponds to the voxel whose value
200 /// in the @a index tree is @e n.
201 /// @param source a tree with a scalar value type
202 /// @param index a tree of the same configuration as @a source but with
203 /// value type VIndex that maps voxels to elements of the output vector
204 template<typename VectorValueType, typename SourceTreeType>
207  const SourceTreeType& source,
208  const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
209 
210 
211 /// @brief Return a tree with the same active voxel topology as the @a index tree
212 /// but whose voxel values are taken from the the given vector.
213 /// @details The voxel whose value in the @a index tree is @e n gets assigned
214 /// the <i>n</i>th element of the vector.
215 /// @param index a tree with value type VIndex that maps voxels to elements of @a values
216 /// @param values a vector of values with which to populate the active voxels of the output tree
217 /// @param background the value for the inactive voxels of the output tree
218 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
219 typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
222  const VIndexTreeType& index,
223  const TreeValueType& background);
224 
225 
226 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
227 /// using second-order finite differences.
228 /// @details This construction assumes homogeneous Dirichlet boundary conditions
229 /// (exterior grid points are zero).
230 template<typename BoolTreeType>
233  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
234  const BoolTreeType& interiorMask,
235  bool staggered = false);
236 
237 
238 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
239 /// with user-specified boundary conditions using second-order finite differences.
240 /// @details Each thread gets its own copy of @a boundaryOp, which should be
241 /// a functor of the form
242 /// @code
243 /// struct BoundaryOp {
244 /// using ValueType = LaplacianMatrix::ValueType;
245 /// void operator()(
246 /// const Coord& ijk, // coordinates of a boundary voxel
247 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
248 /// ValueType& source, // element of source vector corresponding to ijk
249 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
250 /// ) const;
251 /// };
252 /// @endcode
253 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
254 /// and it must specify a boundary condition for @ijk by modifying one or both of two
255 /// provided values: an entry in the given @a source vector corresponding to @ijk and
256 /// the weighting coefficient for @ijk in the %Laplacian matrix.
257 template<typename BoolTreeType, typename BoundaryOp>
260  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
261  const BoolTreeType& interiorMask,
262  const BoundaryOp& boundaryOp,
264  bool staggered = false);
265 
266 
267 /// @brief Dirichlet boundary condition functor
268 /// @details This is useful in describing fluid/air interfaces, where the pressure
269 /// of the air is assumed to be zero.
270 template<typename ValueType>
272  inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
273  // Exterior neighbors are empty, so decrement the weighting coefficient
274  // as for interior neighbors but leave the source vector unchanged.
275  diag -= 1;
276  }
277 };
278 
279 //@}
280 
281 
282 ////////////////////////////////////////
283 
284 /// @cond OPENVDB_DOCS_INTERNAL
285 
286 namespace internal {
287 
288 /// @brief Functor for use with LeafManager::foreach() to populate an array
289 /// with per-leaf active voxel counts
290 template<typename LeafType>
291 struct LeafCountOp
292 {
293  VIndex* count;
294  LeafCountOp(VIndex* count_): count(count_) {}
295  void operator()(const LeafType& leaf, size_t leafIdx) const {
296  count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
297  }
298 };
299 
300 
301 /// @brief Functor for use with LeafManager::foreach() to populate
302 /// active leaf voxels with sequential indices
303 template<typename LeafType>
304 struct LeafIndexOp
305 {
306  const VIndex* count;
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) {
311  it.setValue(idx++);
312  }
313  }
314 };
315 
316 } // namespace internal
317 
318 /// @endcond
319 
320 template<typename VIndexTreeType>
321 inline void
322 populateIndexTree(VIndexTreeType& result)
323 {
324  using LeafT = typename VIndexTreeType::LeafNodeType;
325  using LeafMgrT = typename tree::LeafManager<VIndexTreeType>;
326 
327  // Linearize the tree.
328  LeafMgrT leafManager(result);
329  const size_t leafCount = leafManager.leafCount();
330 
331  if (leafCount == 0) return;
332 
333  // Count the number of active voxels in each leaf node.
334  std::unique_ptr<VIndex[]> perLeafCount(new VIndex[leafCount]);
335  VIndex* perLeafCountPtr = perLeafCount.get();
336  leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
337 
338  // The starting index for each leaf node is the total number
339  // of active voxels in all preceding leaf nodes.
340  for (size_t i = 1; i < leafCount; ++i) {
341  perLeafCount[i] += perLeafCount[i - 1];
342  }
343 
344  // The last accumulated value should be the total of all active voxels.
345  OPENVDB_ASSERT(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
346 
347  // Parallelize over the leaf nodes of the tree, storing a unique index
348  // in each active voxel.
349  leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
350 }
351 
352 
353 template<typename TreeType>
354 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
355 createIndexTree(const TreeType& tree)
356 {
357  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
358 
359  // Construct an output tree with the same active voxel topology as the input tree.
360  const VIndex invalidIdx = -1;
361  typename VIdxTreeT::Ptr result(
362  new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
363 
364  // All active voxels are degrees of freedom, including voxels contained in active tiles.
365  result->voxelizeActiveTiles();
366 
367  populateIndexTree(*result);
368 
369  return result;
370 }
371 
372 
373 ////////////////////////////////////////
374 
375 /// @cond OPENVDB_DOCS_INTERNAL
376 
377 namespace internal {
378 
379 /// @brief Functor for use with LeafManager::foreach() to populate a vector
380 /// with the values of a tree's active voxels
381 template<typename VectorValueType, typename SourceTreeType>
382 struct CopyToVecOp
383 {
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;
388  using VectorT = typename math::pcg::Vector<VectorValueType>;
389 
390  const SourceTreeType* tree;
391  VectorT* vector;
392 
393  CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
394 
395  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
396  {
397  VectorT& vec = *vector;
398  if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
399  // If a corresponding leaf node exists in the source tree,
400  // copy voxel values from the source node to the output vector.
401  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
402  vec[*it] = leaf->getValue(it.pos());
403  }
404  } else {
405  // If no corresponding leaf exists in the source tree,
406  // fill the vector with a uniform value.
407  const TreeValueT& value = tree->getValue(idxLeaf.origin());
408  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
409  vec[*it] = value;
410  }
411  }
412  }
413 };
414 
415 } // namespace internal
416 
417 /// @endcond
418 
419 template<typename VectorValueType, typename SourceTreeType>
421 createVectorFromTree(const SourceTreeType& tree,
422  const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
423 {
424  using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
425  using VIdxLeafMgrT = tree::LeafManager<const VIdxTreeT>;
426  using VectorT = typename math::pcg::Vector<VectorValueType>;
427 
428  // Allocate the vector.
429  const size_t numVoxels = idxTree.activeVoxelCount();
430  typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
431 
432  // Parallelize over the leaf nodes of the index tree, filling the output vector
433  // with values from corresponding voxels of the source tree.
434  VIdxLeafMgrT leafManager(idxTree);
435  leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
436 
437  return result;
438 }
439 
440 
441 ////////////////////////////////////////
442 
443 /// @cond OPENVDB_DOCS_INTERNAL
444 
445 namespace internal {
446 
447 /// @brief Functor for use with LeafManager::foreach() to populate a tree
448 /// with values from a vector
449 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
450 struct CopyFromVecOp
451 {
452  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
453  using OutLeafT = typename OutTreeT::LeafNodeType;
454  using VIdxLeafT = typename VIndexTreeType::LeafNodeType;
455  using VectorT = typename math::pcg::Vector<VectorValueType>;
456 
457  const VectorT* vector;
458  OutTreeT* tree;
459 
460  CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
461 
462  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
463  {
464  const VectorT& vec = *vector;
465  OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
466  OPENVDB_ASSERT(leaf != nullptr);
467  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
468  leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
469  }
470  }
471 };
472 
473 } // namespace internal
474 
475 /// @endcond
476 
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)
483 {
484  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
485  using VIdxLeafMgrT = typename tree::LeafManager<const VIndexTreeType>;
486 
487  // Construct an output tree with the same active voxel topology as the index tree.
488  typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
489  OutTreeT& tree = *result;
490 
491  // Parallelize over the leaf nodes of the index tree, populating voxels
492  // of the output tree with values from the input vector.
493  VIdxLeafMgrT leafManager(idxTree);
494  leafManager.foreach(
495  internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
496 
497  return result;
498 }
499 
500 
501 ////////////////////////////////////////
502 
503 /// @cond OPENVDB_DOCS_INTERNAL
504 
505 namespace internal {
506 
507 /// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
508 template<typename BoolTreeType, typename BoundaryOp>
509 struct ISStaggeredLaplacianOp
510 {
511  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
512  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
513  using ValueT = LaplacianMatrix::ValueType;
514  using VectorT = typename math::pcg::Vector<ValueT>;
515 
517  const VIdxTreeT* idxTree;
518  const BoolTreeType* interiorMask;
519  const BoundaryOp boundaryOp;
520  VectorT* source;
521 
522  ISStaggeredLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx,
523  const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
524  laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
525 
526  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
527  {
528  // Local accessors
529  typename tree::ValueAccessor<const BoolTreeType> interior(*interiorMask);
530  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
531 
532  Coord ijk;
533  VIndex column;
534  const ValueT diagonal = -6.f, offDiagonal = 1.f;
535 
536  // Loop over active voxels in this leaf.
537  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
538  OPENVDB_ASSERT(it.getValue() > -1);
539  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
540 
541  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
542 
543  ijk = it.getCoord();
544  if (interior.isValueOn(ijk)) {
545  // The current voxel is an interior voxel.
546  // All of its neighbors are in the solution domain.
547 
548  // -x direction
549  row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
550  // -y direction
551  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
552  // -z direction
553  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
554  // diagonal
555  row.setValue(rowNum, diagonal);
556  // +z direction
557  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
558  // +y direction
559  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
560  // +x direction
561  row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
562 
563  } else {
564  // The current voxel is a boundary voxel.
565  // At least one of its neighbors is outside the solution domain.
566 
567  ValueT modifiedDiagonal = 0.f;
568 
569  // -x direction
570  if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
571  row.setValue(column, offDiagonal);
572  modifiedDiagonal -= 1;
573  } else {
574  boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
575  }
576  // -y direction
577  if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
578  row.setValue(column, offDiagonal);
579  modifiedDiagonal -= 1;
580  } else {
581  boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
582  }
583  // -z direction
584  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
585  row.setValue(column, offDiagonal);
586  modifiedDiagonal -= 1;
587  } else {
588  boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
589  }
590  // +z direction
591  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
592  row.setValue(column, offDiagonal);
593  modifiedDiagonal -= 1;
594  } else {
595  boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
596  }
597  // +y direction
598  if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
599  row.setValue(column, offDiagonal);
600  modifiedDiagonal -= 1;
601  } else {
602  boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
603  }
604  // +x direction
605  if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
606  row.setValue(column, offDiagonal);
607  modifiedDiagonal -= 1;
608  } else {
609  boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
610  }
611  // diagonal
612  row.setValue(rowNum, modifiedDiagonal);
613  }
614  } // end loop over voxels
615  }
616 };
617 
618 
619 // Stencil 1 is the correct stencil, but stencil 2 requires
620 // half as many comparisons and produces smoother results at boundaries.
621 //#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
622 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
623 
624 /// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
625 template<typename VIdxTreeT, typename BoundaryOp>
626 struct ISLaplacianOp
627 {
628  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
629  using ValueT = LaplacianMatrix::ValueType;
630  using VectorT = typename math::pcg::Vector<ValueT>;
631 
633  const VIdxTreeT* idxTree;
634  const BoundaryOp boundaryOp;
635  VectorT* source;
636 
637  ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
638  laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
639 
640  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
641  {
642  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
643 
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)
648 #else
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)
650 #endif
651  };
652 
653  // For each active voxel in this leaf...
654  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
655  OPENVDB_ASSERT(it.getValue() > -1);
656 
657  const Coord ijk = it.getCoord();
658  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
659 
660  LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
661 
662  ValueT modifiedDiagonal = 0.f;
663 
664  // For each of the neighbors of the voxel at (i,j,k)...
665  for (int dir = 0; dir < kNumOffsets; ++dir) {
666  const Coord neighbor = ijk + ijkOffset[dir];
667  VIndex column;
668  // For collocated vector grids, the central differencing stencil requires
669  // access to neighbors at a distance of two voxels in each direction
670  // (-x, +x, -y, +y, -z, +z).
671 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
672  const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
673  && vectorIdx.isValueOn(neighbor));
674 #else
675  const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
676 #endif
677  if (ijkIsInterior) {
678  // If (i,j,k) is sufficiently far away from the exterior,
679  // set its weight to one and adjust the center weight accordingly.
680  row.setValue(column, 1.f);
681  modifiedDiagonal -= 1.f;
682  } else {
683  // If (i,j,k) is adjacent to or one voxel away from the exterior,
684  // invoke the boundary condition functor.
685  boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
686  }
687  }
688  // Set the (possibly modified) weight for the voxel at (i,j,k).
689  row.setValue(rowNum, modifiedDiagonal);
690  }
691  }
692 };
693 
694 } // namespace internal
695 
696 /// @endcond
697 
698 
699 template<typename BoolTreeType>
701 createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
702  const BoolTreeType& interiorMask, bool staggered)
703 {
704  using ValueT = LaplacianMatrix::ValueType;
706  static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
708  return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
709 }
710 
711 
712 template<typename BoolTreeType, typename BoundaryOp>
715  const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
716  const BoolTreeType& interiorMask,
717  const BoundaryOp& boundaryOp,
719  bool staggered)
720 {
721  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
722  using VIdxLeafMgrT = typename tree::LeafManager<const VIdxTreeT>;
723 
724  // The number of active voxels is the number of degrees of freedom.
725  const Index64 numDoF = idxTree.activeVoxelCount();
726 
727  // Construct the matrix.
728  LaplacianMatrix::Ptr laplacianPtr(
729  new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
730  LaplacianMatrix& laplacian = *laplacianPtr;
731 
732  // Populate the matrix using a second-order, 7-point CD stencil.
733  VIdxLeafMgrT idxLeafManager(idxTree);
734  if (staggered) {
735  idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
736  laplacian, idxTree, interiorMask, boundaryOp, source));
737  } else {
738  idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
739  laplacian, idxTree, boundaryOp, source));
740  }
741 
742  return laplacianPtr;
743 }
744 
745 
746 ////////////////////////////////////////
747 
748 
749 template<typename TreeType>
750 typename TreeType::Ptr
751 solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
752 {
753  util::NullInterrupter interrupter;
754  return solve(inTree, state, interrupter, staggered);
755 }
756 
757 
758 template<typename TreeType, typename Interrupter>
759 typename TreeType::Ptr
760 solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
761 {
763  return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
764 }
765 
766 
767 template<typename TreeType, typename BoundaryOp, typename Interrupter>
768 typename TreeType::Ptr
769 solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
770  math::pcg::State& state, Interrupter& interrupter, bool staggered)
771 {
773  return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
774  inTree, boundaryOp, state, interrupter, staggered);
775 }
776 
777 
778 template<
779  typename PreconditionerType,
780  typename TreeType,
781  typename BoundaryOp,
782  typename Interrupter>
783 typename TreeType::Ptr
785  const TreeType& inTree,
786  const BoundaryOp& boundaryOp,
787  math::pcg::State& state,
788  Interrupter& interrupter,
789  bool staggered)
790 {
791  return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
792  /*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
793 }
794 
795 template<
796  typename PreconditionerType,
797  typename TreeType,
798  typename DomainTreeType,
799  typename BoundaryOp,
800  typename Interrupter>
801 typename TreeType::Ptr
803  const TreeType& inTree,
804  const DomainTreeType& domainMask,
805  const BoundaryOp& boundaryOp,
806  math::pcg::State& state,
807  Interrupter& interrupter,
808  bool staggered)
809 {
810  using TreeValueT = typename TreeType::ValueType;
811  using VecValueT = LaplacianMatrix::ValueType;
812  using VectorT = typename math::pcg::Vector<VecValueT>;
813  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
814  using MaskTreeT = typename TreeType::template ValueConverter<bool>::Type;
815 
816  // 1. Create a mapping from active voxels of the input tree to elements of a vector.
817  typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
818 
819  // 2. Populate a vector with values from the input tree.
820  typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
821 
822  // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
823  /// @todo Is this really needed?
824  typename MaskTreeT::Ptr interiorMask(
825  new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
826  tools::erodeActiveValues(*interiorMask, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
827 
828  // 4. Create the Laplacian matrix.
830  *idxTree, *interiorMask, boundaryOp, *b, staggered);
831 
832  // 5. Solve the Poisson equation.
833  laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
834  b->scale(-1.0);
835  typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
837  new PreconditionerType(*laplacian));
838  if (!precond->isValid()) {
839  precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
840  }
841 
842  state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
843 
844  // 6. Populate the output tree with values from the solution vector.
845  /// @todo if (state.success) ... ?
846  return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
847 }
848 
849 
850 ////////////////////////////////////////
851 
852 
853 // Explicit Template Instantiation
854 
855 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
856 
857 #ifdef OPENVDB_INSTANTIATE_POISSONSOLVER
859 #endif
860 
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)
867 #undef _FUNCTION
868 
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)
875 #undef _FUNCTION
876 
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)
883 #undef _FUNCTION
884 
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)
891 #undef _FUNCTION
892 
893 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
894 
895 
896 } // namespace poisson
897 } // namespace tools
898 } // namespace OPENVDB_VERSION_NAME
899 } // namespace openvdb
900 
901 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
SharedPtr< Preconditioner > Ptr
Definition: ConjGradient.h:468
Definition: Morphology.h:60
void erodeActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically erode all active values (i.e. both voxels and tiles) in a tree using one of three neare...
Definition: Morphology.h:1134
SharedPtr< Vector > Ptr
Definition: ConjGradient.h:142
math::pcg::Vector< VectorValueType >::Ptr createVectorFromTree(const SourceTreeType &source, const typename SourceTreeType::template ValueConverter< VIndex >::Type &index)
Return a vector of the active voxel values of the scalar-valued source tree.
Definition: PoissonSolver.h:421
LaplacianMatrix::Ptr createISLaplacianWithBoundaryConditions(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, const BoundaryOp &boundaryOp, typename math::pcg::Vector< LaplacianMatrix::ValueType > &source, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator with user-specified boundary ...
Definition: PoissonSolver.h:714
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
TreeType::Ptr solveWithBoundaryConditionsAndPreconditioner(const TreeType &, const DomainTreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the valu...
Definition: PoissonSolver.h:802
TreeType::Ptr solve(const TreeType &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x, where b is a vector comprising the values of all of the active voxels in the in...
Definition: PoissonSolver.h:760
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
typename Adapter::TreeType OutTreeT
Definition: ValueTransformer.h:595
Lightweight, variable-length vector.
Definition: ConjGradient.h:37
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
LaplacianMatrix::Ptr createISLaplacian(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator using second-order finite dif...
Definition: PoissonSolver.h:701
int32_t Int32
Definition: Types.h:56
Base class for interrupters.
Definition: NullInterrupter.h:25
OutGridT XformOp & op
Definition: ValueTransformer.h:139
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
Definition: PoissonSolver.h:85
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:1018
GridType::template ValueConverter< bool >::Type::Ptr interiorMask(const GridType &grid, const double isovalue=0.0)
Given an input grid of any type, return a new, boolean grid whose active voxel topology matches the i...
Definition: Mask.h:105
Read/write accessor to a row of this matrix.
Definition: ConjGradient.h:419
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:43
Dirichlet boundary condition functor.
Definition: PoissonSolver.h:271
VIndexTreeType::template ValueConverter< TreeValueType >::Type::Ptr createTreeFromVector(const math::pcg::Vector< VectorValueType > &values, const VIndexTreeType &index, const TreeValueType &background)
Return a tree with the same active voxel topology as the index tree but whose voxel values are taken ...
Definition: PoissonSolver.h:479
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
Definition: Morphology.h:82
Index32 SizeType
Definition: ConjGradient.h:33
Int32 VIndex
Definition: PoissonSolver.h:82
TreeType::Ptr solveWithBoundaryConditions(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the valu...
Definition: PoissonSolver.h:769
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
TreeType::template ValueConverter< VIndex >::Type::Ptr createIndexTree(const TreeType &)
Iterate over the active voxels of the input tree and for each one assign its index in the iteration s...
Definition: PoissonSolver.h:355
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero...
Definition: PoissonSolver.h:322
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
Definition: PoissonSolver.h:272
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&#39;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&#39;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