HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PoissonSolver.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
29 ///////////////////////////////////////////////////////////////////////////
30 //
31 /// @file PoissonSolver.h
32 ///
33 /// @authors D.J. Hill, Peter Cucka
34 ///
35 /// @brief Solve Poisson's equation &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i>
36 /// for <i>x</i>, where @e b is a vector comprising the values of all of the active voxels
37 /// in a grid.
38 ///
39 /// @par Example:
40 /// Solve for the pressure in a cubic tank of liquid, assuming uniform boundary conditions:
41 /// @code
42 /// FloatTree source(/*background=*/0.0f);
43 /// // Activate voxels to indicate that they contain liquid.
44 /// source.fill(CoordBBox(Coord(0, -10, 0), Coord(10, 0, 10)), /*value=*/0.0f);
45 ///
46 /// math::pcg::State state = math::pcg::terminationDefaults<float>();
47 /// FloatTree::Ptr solution = tools::poisson::solve(source, state);
48 /// @endcode
49 ///
50 /// @par Example:
51 /// Solve for the pressure, <i>P</i>, in a cubic tank of liquid that is open at the top.
52 /// Boundary conditions are <i>P</i>&nbsp;=&nbsp;0 at the top,
53 /// &part;<i>P</i>/&part;<i>y</i>&nbsp;=&nbsp;&minus;1 at the bottom
54 /// and &part;<i>P</i>/&part;<i>x</i>&nbsp;=&nbsp;0 at the sides:
55 /// <pre>
56 /// P = 0
57 /// +--------+ (N,0,N)
58 /// /| /|
59 /// (0,0,0) +--------+ |
60 /// | | | | dP/dx = 0
61 /// dP/dx = 0 | +------|-+
62 /// |/ |/
63 /// (0,-N,0) +--------+ (N,-N,0)
64 /// dP/dy = -1
65 /// </pre>
66 /// @code
67 /// const int N = 10;
68 /// DoubleTree source(/*background=*/0.0);
69 /// // Activate voxels to indicate that they contain liquid.
70 /// source.fill(CoordBBox(Coord(0, -N, 0), Coord(N, 0, N)), /*value=*/0.0);
71 ///
72 /// auto boundary = [](const openvdb::Coord& ijk, const openvdb::Coord& neighbor,
73 /// double& source, double& diagonal)
74 /// {
75 /// if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) {
76 /// if (neighbor.y() < ijk.y()) source -= 1.0;
77 /// else diagonal -= 1.0;
78 /// }
79 /// };
80 ///
81 /// math::pcg::State state = math::pcg::terminationDefaults<double>();
82 /// util::NullInterrupter interrupter;
83 ///
84 /// DoubleTree::Ptr solution = tools::poisson::solveWithBoundaryConditions(
85 /// source, boundary, state, interrupter);
86 /// @endcode
87 
88 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
89 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
90 
91 #include <openvdb/Types.h>
94 #include <openvdb/tree/Tree.h>
96 #include "Morphology.h" // for erodeVoxels
97 
98 
99 namespace openvdb {
101 namespace OPENVDB_VERSION_NAME {
102 namespace tools {
103 namespace poisson {
104 
105 // This type should be at least as wide as math::pcg::SizeType.
106 using VIndex = Int32;
107 
108 /// The type of a matrix used to represent a three-dimensional %Laplacian operator
110 
111 
112 //@{
113 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>,
114 /// where @e b is a vector comprising the values of all of the active voxels
115 /// in the input tree.
116 /// @return a new tree, with the same active voxel topology as the input tree,
117 /// whose voxel values are the elements of the solution vector <i>x</i>.
118 /// @details On input, the State object should specify convergence criteria
119 /// (minimum error and maximum number of iterations); on output, it gives
120 /// the actual termination conditions.
121 /// @details The solution is computed using the conjugate gradient method
122 /// with (where possible) incomplete Cholesky preconditioning, falling back
123 /// to Jacobi preconditioning.
124 /// @sa solveWithBoundaryConditions
125 template<typename TreeType>
126 inline typename TreeType::Ptr
127 solve(const TreeType&, math::pcg::State&, bool staggered = false);
128 
129 template<typename TreeType, typename Interrupter>
130 inline typename TreeType::Ptr
131 solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
132 //@}
133 
134 
135 //@{
136 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>
137 /// with user-specified boundary conditions, where @e b is a vector comprising
138 /// the values of all of the active voxels in the input tree or domain mask if provided
139 /// @return a new tree, with the same active voxel topology as the input tree,
140 /// whose voxel values are the elements of the solution vector <i>x</i>.
141 /// @details On input, the State object should specify convergence criteria
142 /// (minimum error and maximum number of iterations); on output, it gives
143 /// the actual termination conditions.
144 /// @details The solution is computed using the conjugate gradient method with
145 /// the specified type of preconditioner (default: incomplete Cholesky),
146 /// falling back to Jacobi preconditioning if necessary.
147 /// @details Each thread gets its own copy of the BoundaryOp, which should be
148 /// a functor of the form
149 /// @code
150 /// struct BoundaryOp {
151 /// using ValueType = LaplacianMatrix::ValueType;
152 /// void operator()(
153 /// const Coord& ijk, // coordinates of a boundary voxel
154 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
155 /// ValueType& source, // element of b corresponding to ijk
156 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
157 /// ) const;
158 /// };
159 /// @endcode
160 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
161 /// and it must specify a boundary condition for @ijk by modifying one or both of two
162 /// provided values: the entry in the source vector @e b corresponding to @ijk and
163 /// the weighting coefficient for @ijk in the Laplacian operator matrix.
164 ///
165 /// @sa solve
166 template<typename TreeType, typename BoundaryOp, typename Interrupter>
167 inline typename TreeType::Ptr
169  const TreeType&,
170  const BoundaryOp&,
172  Interrupter&,
173  bool staggered = false);
174 
175 template<
176  typename PreconditionerType,
177  typename TreeType,
178  typename BoundaryOp,
179  typename Interrupter>
180 inline typename TreeType::Ptr
182  const TreeType&,
183  const BoundaryOp&,
185  Interrupter&,
186  bool staggered = false);
187 
188 template<
189  typename PreconditionerType,
190  typename TreeType,
191  typename DomainTreeType,
192  typename BoundaryOp,
193  typename Interrupter>
194 inline typename TreeType::Ptr
196  const TreeType&,
197  const DomainTreeType&,
198  const BoundaryOp&,
200  Interrupter&,
201  bool staggered = false);
202 //@}
203 
204 
205 /// @name Low-level functions
206 //@{
207 // The following are low-level routines that can be used to assemble custom solvers.
208 
209 /// @brief Overwrite each active voxel in the given scalar tree
210 /// with a sequential index, starting from zero.
211 template<typename VIndexTreeType>
212 inline void populateIndexTree(VIndexTreeType&);
213 
214 /// @brief Iterate over the active voxels of the input tree and for each one
215 /// assign its index in the iteration sequence to the corresponding voxel
216 /// of an integer-valued output tree.
217 template<typename TreeType>
218 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
219 createIndexTree(const TreeType&);
220 
221 
222 /// @brief Return a vector of the active voxel values of the scalar-valued @a source tree.
223 /// @details The <i>n</i>th element of the vector corresponds to the voxel whose value
224 /// in the @a index tree is @e n.
225 /// @param source a tree with a scalar value type
226 /// @param index a tree of the same configuration as @a source but with
227 /// value type VIndex that maps voxels to elements of the output vector
228 template<typename VectorValueType, typename SourceTreeType>
231  const SourceTreeType& source,
232  const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
233 
234 
235 /// @brief Return a tree with the same active voxel topology as the @a index tree
236 /// but whose voxel values are taken from the the given vector.
237 /// @details The voxel whose value in the @a index tree is @e n gets assigned
238 /// the <i>n</i>th element of the vector.
239 /// @param index a tree with value type VIndex that maps voxels to elements of @a values
240 /// @param values a vector of values with which to populate the active voxels of the output tree
241 /// @param background the value for the inactive voxels of the output tree
242 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
243 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
246  const VIndexTreeType& index,
247  const TreeValueType& background);
248 
249 
250 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
251 /// using second-order finite differences.
252 /// @details This construction assumes homogeneous Dirichlet boundary conditions
253 /// (exterior grid points are zero).
254 template<typename BoolTreeType>
257  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
258  const BoolTreeType& interiorMask,
259  bool staggered = false);
260 
261 
262 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
263 /// with user-specified boundary conditions using second-order finite differences.
264 /// @details Each thread gets its own copy of @a boundaryOp, which should be
265 /// a functor of the form
266 /// @code
267 /// struct BoundaryOp {
268 /// using ValueType = LaplacianMatrix::ValueType;
269 /// void operator()(
270 /// const Coord& ijk, // coordinates of a boundary voxel
271 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
272 /// ValueType& source, // element of source vector corresponding to ijk
273 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
274 /// ) const;
275 /// };
276 /// @endcode
277 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
278 /// and it must specify a boundary condition for @ijk by modifying one or both of two
279 /// provided values: an entry in the given @a source vector corresponding to @ijk and
280 /// the weighting coefficient for @ijk in the %Laplacian matrix.
281 template<typename BoolTreeType, typename BoundaryOp>
284  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
285  const BoolTreeType& interiorMask,
286  const BoundaryOp& boundaryOp,
288  bool staggered = false);
289 
290 
291 /// @brief Dirichlet boundary condition functor
292 /// @details This is useful in describing fluid/air interfaces, where the pressure
293 /// of the air is assumed to be zero.
294 template<typename ValueType>
296  inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
297  // Exterior neighbors are empty, so decrement the weighting coefficient
298  // as for interior neighbors but leave the source vector unchanged.
299  diag -= 1;
300  }
301 };
302 
303 //@}
304 
305 
306 ////////////////////////////////////////
307 
308 
309 namespace internal {
310 
311 /// @brief Functor for use with LeafManager::foreach() to populate an array
312 /// with per-leaf active voxel counts
313 template<typename LeafType>
315 {
317  LeafCountOp(VIndex* count_): count(count_) {}
318  void operator()(const LeafType& leaf, size_t leafIdx) const {
319  count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
320  }
321 };
322 
323 
324 /// @brief Functor for use with LeafManager::foreach() to populate
325 /// active leaf voxels with sequential indices
326 template<typename LeafType>
328 {
329  const VIndex* count;
330  LeafIndexOp(const VIndex* count_): count(count_) {}
331  void operator()(LeafType& leaf, size_t leafIdx) const {
332  VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
333  for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
334  it.setValue(idx++);
335  }
336  }
337 };
338 
339 } // namespace internal
340 
341 
342 template<typename VIndexTreeType>
343 inline void
344 populateIndexTree(VIndexTreeType& result)
345 {
346  using LeafT = typename VIndexTreeType::LeafNodeType;
347  using LeafMgrT = typename tree::LeafManager<VIndexTreeType>;
348 
349  // Linearize the tree.
350  LeafMgrT leafManager(result);
351  const size_t leafCount = leafManager.leafCount();
352 
353  if (leafCount == 0) return;
354 
355  // Count the number of active voxels in each leaf node.
356  std::unique_ptr<VIndex[]> perLeafCount(new VIndex[leafCount]);
357  VIndex* perLeafCountPtr = perLeafCount.get();
358  leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
359 
360  // The starting index for each leaf node is the total number
361  // of active voxels in all preceding leaf nodes.
362  for (size_t i = 1; i < leafCount; ++i) {
363  perLeafCount[i] += perLeafCount[i - 1];
364  }
365 
366  // The last accumulated value should be the total of all active voxels.
367  assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
368 
369  // Parallelize over the leaf nodes of the tree, storing a unique index
370  // in each active voxel.
371  leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
372 }
373 
374 
375 template<typename TreeType>
376 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
377 createIndexTree(const TreeType& tree)
378 {
379  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
380 
381  // Construct an output tree with the same active voxel topology as the input tree.
382  const VIndex invalidIdx = -1;
383  typename VIdxTreeT::Ptr result(
384  new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
385 
386  // All active voxels are degrees of freedom, including voxels contained in active tiles.
387  result->voxelizeActiveTiles();
388 
389  populateIndexTree(*result);
390 
391  return result;
392 }
393 
394 
395 ////////////////////////////////////////
396 
397 
398 namespace internal {
399 
400 /// @brief Functor for use with LeafManager::foreach() to populate a vector
401 /// with the values of a tree's active voxels
402 template<typename VectorValueType, typename SourceTreeType>
404 {
405  using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
406  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
407  using LeafT = typename SourceTreeType::LeafNodeType;
408  using TreeValueT = typename SourceTreeType::ValueType;
410 
411  const SourceTreeType* tree;
413 
414  CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
415 
416  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
417  {
418  VectorT& vec = *vector;
419  if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
420  // If a corresponding leaf node exists in the source tree,
421  // copy voxel values from the source node to the output vector.
422  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
423  vec[*it] = leaf->getValue(it.pos());
424  }
425  } else {
426  // If no corresponding leaf exists in the source tree,
427  // fill the vector with a uniform value.
428  const TreeValueT& value = tree->getValue(idxLeaf.origin());
429  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
430  vec[*it] = value;
431  }
432  }
433  }
434 };
435 
436 } // namespace internal
437 
438 
439 template<typename VectorValueType, typename SourceTreeType>
441 createVectorFromTree(const SourceTreeType& tree,
442  const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
443 {
444  using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
445  using VIdxLeafMgrT = tree::LeafManager<const VIdxTreeT>;
447 
448  // Allocate the vector.
449  const size_t numVoxels = idxTree.activeVoxelCount();
450  typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
451 
452  // Parallelize over the leaf nodes of the index tree, filling the output vector
453  // with values from corresponding voxels of the source tree.
454  VIdxLeafMgrT leafManager(idxTree);
455  leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
456 
457  return result;
458 }
459 
460 
461 ////////////////////////////////////////
462 
463 
464 namespace internal {
465 
466 /// @brief Functor for use with LeafManager::foreach() to populate a tree
467 /// with values from a vector
468 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
470 {
471  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
472  using OutLeafT = typename OutTreeT::LeafNodeType;
473  using VIdxLeafT = typename VIndexTreeType::LeafNodeType;
475 
476  const VectorT* vector;
478 
479  CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
480 
481  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
482  {
483  const VectorT& vec = *vector;
484  OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
485  assert(leaf != nullptr);
486  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
487  leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
488  }
489  }
490 };
491 
492 } // namespace internal
493 
494 
495 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
496 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
499  const VIndexTreeType& idxTree,
500  const TreeValueType& background)
501 {
502  using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
503  using VIdxLeafMgrT = typename tree::LeafManager<const VIndexTreeType>;
504 
505  // Construct an output tree with the same active voxel topology as the index tree.
506  typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
507  OutTreeT& tree = *result;
508 
509  // Parallelize over the leaf nodes of the index tree, populating voxels
510  // of the output tree with values from the input vector.
511  VIdxLeafMgrT leafManager(idxTree);
512  leafManager.foreach(
514 
515  return result;
516 }
517 
518 
519 ////////////////////////////////////////
520 
521 
522 namespace internal {
523 
524 /// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
525 template<typename BoolTreeType, typename BoundaryOp>
527 {
528  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
529  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
532 
535  const BoolTreeType* interiorMask;
536  const BoundaryOp boundaryOp;
538 
540  const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
541  laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
542 
543  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
544  {
545  // Local accessors
547  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
548 
549  Coord ijk;
550  VIndex column;
551  const ValueT diagonal = -6.f, offDiagonal = 1.f;
552 
553  // Loop over active voxels in this leaf.
554  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
555  assert(it.getValue() > -1);
556  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
557 
559 
560  ijk = it.getCoord();
561  if (interior.isValueOn(ijk)) {
562  // The current voxel is an interior voxel.
563  // All of its neighbors are in the solution domain.
564 
565  // -x direction
566  row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
567  // -y direction
568  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
569  // -z direction
570  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
571  // diagonal
572  row.setValue(rowNum, diagonal);
573  // +z direction
574  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
575  // +y direction
576  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
577  // +x direction
578  row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
579 
580  } else {
581  // The current voxel is a boundary voxel.
582  // At least one of its neighbors is outside the solution domain.
583 
584  ValueT modifiedDiagonal = 0.f;
585 
586  // -x direction
587  if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
588  row.setValue(column, offDiagonal);
589  modifiedDiagonal -= 1;
590  } else {
591  boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
592  }
593  // -y direction
594  if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
595  row.setValue(column, offDiagonal);
596  modifiedDiagonal -= 1;
597  } else {
598  boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
599  }
600  // -z direction
601  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
602  row.setValue(column, offDiagonal);
603  modifiedDiagonal -= 1;
604  } else {
605  boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
606  }
607  // +z direction
608  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
609  row.setValue(column, offDiagonal);
610  modifiedDiagonal -= 1;
611  } else {
612  boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
613  }
614  // +y direction
615  if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
616  row.setValue(column, offDiagonal);
617  modifiedDiagonal -= 1;
618  } else {
619  boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
620  }
621  // +x direction
622  if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
623  row.setValue(column, offDiagonal);
624  modifiedDiagonal -= 1;
625  } else {
626  boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
627  }
628  // diagonal
629  row.setValue(rowNum, modifiedDiagonal);
630  }
631  } // end loop over voxels
632  }
633 };
634 
635 
636 // Stencil 1 is the correct stencil, but stencil 2 requires
637 // half as many comparisons and produces smoother results at boundaries.
638 //#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
639 #define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
640 
641 /// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
642 template<typename VIdxTreeT, typename BoundaryOp>
644 {
645  using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
648 
650  const VIdxTreeT* idxTree;
651  const BoundaryOp boundaryOp;
653 
654  ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
655  laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
656 
657  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
658  {
659  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
660 
661  const int kNumOffsets = 6;
662  const Coord ijkOffset[kNumOffsets] = {
663 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
664  Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
665 #else
666  Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
667 #endif
668  };
669 
670  // For each active voxel in this leaf...
671  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
672  assert(it.getValue() > -1);
673 
674  const Coord ijk = it.getCoord();
675  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
676 
678 
679  ValueT modifiedDiagonal = 0.f;
680 
681  // For each of the neighbors of the voxel at (i,j,k)...
682  for (int dir = 0; dir < kNumOffsets; ++dir) {
683  const Coord neighbor = ijk + ijkOffset[dir];
684  VIndex column;
685  // For collocated vector grids, the central differencing stencil requires
686  // access to neighbors at a distance of two voxels in each direction
687  // (-x, +x, -y, +y, -z, +z).
688 #if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
689  const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
690  && vectorIdx.isValueOn(neighbor));
691 #else
692  const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
693 #endif
694  if (ijkIsInterior) {
695  // If (i,j,k) is sufficiently far away from the exterior,
696  // set its weight to one and adjust the center weight accordingly.
697  row.setValue(column, 1.f);
698  modifiedDiagonal -= 1.f;
699  } else {
700  // If (i,j,k) is adjacent to or one voxel away from the exterior,
701  // invoke the boundary condition functor.
702  boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
703  }
704  }
705  // Set the (possibly modified) weight for the voxel at (i,j,k).
706  row.setValue(rowNum, modifiedDiagonal);
707  }
708  }
709 };
710 
711 } // namespace internal
712 
713 
714 template<typename BoolTreeType>
716 createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
717  const BoolTreeType& interiorMask, bool staggered)
718 {
719  using ValueT = LaplacianMatrix::ValueType;
721  static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
723  return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
724 }
725 
726 
727 template<typename BoolTreeType, typename BoundaryOp>
730  const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
731  const BoolTreeType& interiorMask,
732  const BoundaryOp& boundaryOp,
734  bool staggered)
735 {
736  using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
737  using VIdxLeafMgrT = typename tree::LeafManager<const VIdxTreeT>;
738 
739  // The number of active voxels is the number of degrees of freedom.
740  const Index64 numDoF = idxTree.activeVoxelCount();
741 
742  // Construct the matrix.
743  LaplacianMatrix::Ptr laplacianPtr(
744  new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
745  LaplacianMatrix& laplacian = *laplacianPtr;
746 
747  // Populate the matrix using a second-order, 7-point CD stencil.
748  VIdxLeafMgrT idxLeafManager(idxTree);
749  if (staggered) {
751  laplacian, idxTree, interiorMask, boundaryOp, source));
752  } else {
753  idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
754  laplacian, idxTree, boundaryOp, source));
755  }
756 
757  return laplacianPtr;
758 }
759 
760 
761 ////////////////////////////////////////
762 
763 
764 template<typename TreeType>
765 inline typename TreeType::Ptr
766 solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
767 {
768  util::NullInterrupter interrupter;
769  return solve(inTree, state, interrupter, staggered);
770 }
771 
772 
773 template<typename TreeType, typename Interrupter>
774 inline typename TreeType::Ptr
775 solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
776 {
778  return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
779 }
780 
781 
782 template<typename TreeType, typename BoundaryOp, typename Interrupter>
783 inline typename TreeType::Ptr
784 solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
785  math::pcg::State& state, Interrupter& interrupter, bool staggered)
786 {
788  return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
789  inTree, boundaryOp, state, interrupter, staggered);
790 }
791 
792 
793 template<
794  typename PreconditionerType,
795  typename TreeType,
796  typename BoundaryOp,
797  typename Interrupter>
798 inline typename TreeType::Ptr
800  const TreeType& inTree,
801  const BoundaryOp& boundaryOp,
802  math::pcg::State& state,
803  Interrupter& interrupter,
804  bool staggered)
805 {
806  return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
807  /*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
808 }
809 
810 template<
811  typename PreconditionerType,
812  typename TreeType,
813  typename DomainTreeType,
814  typename BoundaryOp,
815  typename Interrupter>
816 inline typename TreeType::Ptr
818  const TreeType& inTree,
819  const DomainTreeType& domainMask,
820  const BoundaryOp& boundaryOp,
821  math::pcg::State& state,
822  Interrupter& interrupter,
823  bool staggered)
824 {
825  using TreeValueT = typename TreeType::ValueType;
826  using VecValueT = LaplacianMatrix::ValueType;
827  using VectorT = typename math::pcg::Vector<VecValueT>;
828  using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
829  using MaskTreeT = typename TreeType::template ValueConverter<bool>::Type;
830 
831  // 1. Create a mapping from active voxels of the input tree to elements of a vector.
832  typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
833 
834  // 2. Populate a vector with values from the input tree.
835  typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
836 
837  // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
838  /// @todo Is this really needed?
839  typename MaskTreeT::Ptr interiorMask(
840  new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
841  tools::erodeVoxels(*interiorMask, /*iterations=*/1, tools::NN_FACE);
842 
843  // 4. Create the Laplacian matrix.
845  *idxTree, *interiorMask, boundaryOp, *b, staggered);
846 
847  // 5. Solve the Poisson equation.
848  laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
849  b->scale(-1.0);
850  typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
852  new PreconditionerType(*laplacian));
853  if (!precond->isValid()) {
854  precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
855  }
856 
857  state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
858 
859  // 6. Populate the output tree with values from the solution vector.
860  /// @todo if (state.success) ... ?
861  return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
862 }
863 
864 } // namespace poisson
865 } // namespace tools
866 } // namespace OPENVDB_VERSION_NAME
867 } // namespace openvdb
868 
869 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
870 
871 // Copyright (c) 2012-2018 DreamWorks Animation LLC
872 // All rights reserved. This software is distributed under the
873 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
RowEditor getRowEditor(SizeType row)
Return a read/write view onto the given row of this matrix.
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Lightweight, variable-length vector.
Definition: ConjGradient.h:64
Functor for use with LeafManager::foreach() to populate an array with per-leaf active voxel counts...
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
void operator()(LeafType &leaf, size_t leafIdx) const
const GLdouble * v
Definition: glcorearb.h:836
Functor for use with LeafManager::foreach() to populate a tree with values from a vector...
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
png_infop png_color_16p * background
Definition: png.h:2326
ISStaggeredLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoolTreeType &mask, const BoundaryOp &op, VectorT &src)
GLint GLuint mask
Definition: glcorearb.h:123
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
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 ...
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:73
Dummy NOOP interrupter class defining interface.
png_uint_32 i
Definition: png.h:2877
Functor for use with LeafManager::foreach() to populate a vector with the values of a tree's active v...
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition: ConjGradient.h:66
typename VIndexTreeType::template ValueConverter< TreeValueType >::Type OutTreeT
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
ISLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoundaryOp &op, VectorT &src)
GLfloat f
Definition: glcorearb.h:1925
void operator()(const VIdxLeafT &idxLeaf, size_t) const
Preconditioner using incomplete Cholesky factorization.
Definition: ConjGradient.h:70
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
State solve(const PositiveDefMatrix &A, const Vector< typename PositiveDefMatrix::ValueType > &b, Vector< typename PositiveDefMatrix::ValueType > &x, Preconditioner< typename PositiveDefMatrix::ValueType > &preconditioner, const State &termination=terminationDefaults< typename PositiveDefMatrix::ValueType >())
Solve Ax = b via the preconditioned conjugate gradient method.
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 ...
GLsizei GLsizei GLchar * source
Definition: glcorearb.h:802
typename BoolTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1601
png_bytepp row
Definition: png.h:1836
typename math::pcg::Vector< VectorValueType > VectorT
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
void operator()(const LeafType &leaf, size_t leafIdx) const
GLint GLsizei count
Definition: glcorearb.h:404
void operator()(const VIdxLeafT &idxLeaf, size_t) const
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:138
TreeType::Ptr solveWithBoundaryConditionsAndPreconditioner(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...
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
OPENVDB_STATIC_SPECIALIZATION void erodeVoxels(TreeType &tree, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically erode all leaf-level active voxels in the given tree.
Definition: Morphology.h:878
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.
GLsizei const GLfloat * value
Definition: glcorearb.h:823
UT_VectorT< T > VectorT
Definition: GU_UVUtils.h:33
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...
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero...
GLuint index
Definition: glcorearb.h:785
GLint GLenum GLint x
Definition: glcorearb.h:408
Implementation of morphological dilation and erosion.
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:518
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
typename SourceTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
TreeType::Ptr solve(const TreeType &, math::pcg::State &, 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...
Functor for use with LeafManager::foreach() to populate active leaf voxels with sequential indices...
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...
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...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix.
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
GLenum src
Definition: glcorearb.h:1792