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