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