HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PoissonSolver.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2017 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,N)
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 /// // C++11
73 /// auto boundary = [](const openvdb::Coord& ijk, const openvdb::Coord& neighbor,
74 /// double& source, double& diagonal)
75 /// {
76 /// if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) {
77 /// if (neighbor.y() < ijk.y()) source -= 1.0;
78 /// else diagonal -= 1.0;
79 /// }
80 /// };
81 ///
82 /// math::pcg::State state = math::pcg::terminationDefaults<double>();
83 /// util::NullInterrupter interrupter;
84 ///
85 /// DoubleTree::Ptr solution = tools::poisson::solveWithBoundaryConditions(
86 /// source, boundary, state, interrupter);
87 /// @endcode
88 
89 #ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
90 #define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
91 
92 #include <openvdb/Types.h>
95 #include <openvdb/tree/Tree.h>
97 
98 #include "Morphology.h" // for erodeVoxels
99 
100 #include <hboost/scoped_array.hpp>
101 
102 
103 namespace openvdb {
105 namespace OPENVDB_VERSION_NAME {
106 namespace tools {
107 namespace poisson {
108 
109 // This type should be at least as wide as math::pcg::SizeType.
110 typedef Int32 VIndex;
111 
112 /// The type of a matrix used to represent a three-dimensional Laplacian operator
114 
115 
116 //@{
117 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>,
118 /// where @e b is a vector comprising the values of all of the active voxels
119 /// in the input tree.
120 /// @return a new tree, with the same active voxel topology as the input tree,
121 /// whose voxel values are the elements of the solution vector <i>x</i>.
122 /// @details On input, the State object should specify convergence criteria
123 /// (minimum error and maximum number of iterations); on output, it gives
124 /// the actual termination conditions.
125 /// @details The solution is computed using the conjugate gradient method
126 /// with (where possible) incomplete Cholesky preconditioning, falling back
127 /// to Jacobi preconditioning.
128 /// @sa solveWithBoundaryConditions
129 template<typename TreeType>
130 inline typename TreeType::Ptr
131 solve(const TreeType&, math::pcg::State&);
132 
133 template<typename TreeType, typename Interrupter>
134 inline typename TreeType::Ptr
135 solve(const TreeType&, math::pcg::State&, Interrupter&);
136 //@}
137 
138 
139 //@{
140 /// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>
141 /// with user-specified boundary conditions, where @e b is a vector comprising
142 /// the values of all of the active voxels in the input tree or domain mask if provided
143 /// @return a new tree, with the same active voxel topology as the input tree,
144 /// whose voxel values are the elements of the solution vector <i>x</i>.
145 /// @details On input, the State object should specify convergence criteria
146 /// (minimum error and maximum number of iterations); on output, it gives
147 /// the actual termination conditions.
148 /// @details The solution is computed using the conjugate gradient method with
149 /// the specified type of preconditioner (default: incomplete Cholesky),
150 /// falling back to Jacobi preconditioning if necessary.
151 /// @details Each thread gets its own copy of the BoundaryOp, which should be
152 /// a functor of the form
153 /// @code
154 /// struct BoundaryOp {
155 /// typedef LaplacianMatrix::ValueType ValueType;
156 /// void operator()(
157 /// const Coord& ijk, // coordinates of a boundary voxel
158 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
159 /// ValueType& source, // element of b corresponding to ijk
160 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
161 /// ) const;
162 /// };
163 /// @endcode
164 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
165 /// and it must specify a boundary condition for @ijk by modifying one or both of two
166 /// provided values: the entry in the source vector @e b corresponding to @ijk and
167 /// the weighting coefficient for @ijk in the Laplacian operator matrix.
168 ///
169 /// @sa solve
170 template<typename TreeType, typename BoundaryOp, typename Interrupter>
171 inline typename TreeType::Ptr
172 solveWithBoundaryConditions(const TreeType&, const BoundaryOp&, math::pcg::State&, Interrupter&);
173 
174 template<typename PreconditionerType, typename TreeType, typename BoundaryOp, typename Interrupter>
175 inline typename TreeType::Ptr
176 solveWithBoundaryConditionsAndPreconditioner(const TreeType&, const BoundaryOp&,
177  math::pcg::State&, Interrupter&);
178 
179 template<typename PreconditionerType, typename TreeType, typename DomainTreeType, typename BoundaryOp, typename Interrupter>
180 inline typename TreeType::Ptr
181 solveWithBoundaryConditionsAndPreconditioner(const TreeType&, const DomainTreeType&, const BoundaryOp&,
182  math::pcg::State&, Interrupter&);
183 //@}
184 
185 
186 /// @name Low-level functions
187 //@{
188 // The following are low-level routines that can be used to assemble custom solvers.
189 
190 /// @brief Overwrite each active voxel in the given scalar tree
191 /// with a sequential index, starting from zero.
192 template<typename VIndexTreeType>
193 inline void populateIndexTree(VIndexTreeType&);
194 
195 /// @brief Iterate over the active voxels of the input tree and for each one
196 /// assign its index in the iteration sequence to the corresponding voxel
197 /// of an integer-valued output tree.
198 template<typename TreeType>
199 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
200 createIndexTree(const TreeType&);
201 
202 
203 /// @brief Return a vector of the active voxel values of the scalar-valued @a source tree.
204 /// @details The <i>n</i>th element of the vector corresponds to the voxel whose value
205 /// in the @a index tree is @e n.
206 /// @param source a tree with a scalar value type
207 /// @param index a tree of the same configuration as @a source but with
208 /// value type VIndex that maps voxels to elements of the output vector
209 template<typename VectorValueType, typename SourceTreeType>
212  const SourceTreeType& source,
213  const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
214 
215 
216 /// @brief Return a tree with the same active voxel topology as the @a index tree
217 /// but whose voxel values are taken from the the given vector.
218 /// @details The voxel whose value in the @a index tree is @e n gets assigned
219 /// the <i>n</i>th element of the vector.
220 /// @param index a tree with value type VIndex that maps voxels to elements of @a values
221 /// @param values a vector of values with which to populate the active voxels of the output tree
222 /// @param background the value for the inactive voxels of the output tree
223 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
224 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
227  const VIndexTreeType& index,
228  const TreeValueType& background);
229 
230 
231 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) Laplacian operator
232 /// using second-order finite differences.
233 /// @details This construction assumes homogeneous Dirichlet boundary conditions
234 /// (exterior grid points are zero).
235 template<typename BoolTreeType>
238  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
239  const BoolTreeType& interiorMask);
240 
241 
242 /// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) Laplacian operator
243 /// with user-specified boundary conditions using second-order finite differences.
244 /// @details Each thread gets its own copy of @a boundaryOp, which should be
245 /// a functor of the form
246 /// @code
247 /// struct BoundaryOp {
248 /// typedef LaplacianMatrix::ValueType ValueType;
249 /// void operator()(
250 /// const Coord& ijk, // coordinates of a boundary voxel
251 /// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
252 /// ValueType& source, // element of source vector corresponding to ijk
253 /// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
254 /// ) const;
255 /// };
256 /// @endcode
257 /// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
258 /// and it must specify a boundary condition for @ijk by modifying one or both of two
259 /// provided values: an entry in the given @a source vector corresponding to @ijk and
260 /// the weighting coefficient for @ijk in the Laplacian matrix.
261 template<typename BoolTreeType, typename BoundaryOp>
264  const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
265  const BoolTreeType& interiorMask,
266  const BoundaryOp& boundaryOp,
268 
269 //@}
270 
271 
272 ////////////////////////////////////////
273 
274 
275 namespace internal {
276 
277 /// @brief Functor for use with LeafManager::foreach() to populate an array
278 /// with per-leaf active voxel counts
279 template<typename LeafType>
281 {
283  LeafCountOp(VIndex* count_): count(count_) {}
284  void operator()(const LeafType& leaf, size_t leafIdx) const {
285  count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
286  }
287 };
288 
289 
290 /// @brief Functor for use with LeafManager::foreach() to populate
291 /// active leaf voxels with sequential indices
292 template<typename LeafType>
294 {
295  const VIndex* count;
296  LeafIndexOp(const VIndex* count_): count(count_) {}
297  void operator()(LeafType& leaf, size_t leafIdx) const {
298  VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
299  for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
300  it.setValue(idx++);
301  }
302  }
303 };
304 
305 } // namespace internal
306 
307 
308 template<typename VIndexTreeType>
309 inline void
310 populateIndexTree(VIndexTreeType& result)
311 {
312  typedef typename VIndexTreeType::LeafNodeType LeafT;
313  typedef typename tree::LeafManager<VIndexTreeType> LeafMgrT;
314 
315  // Linearize the tree.
316  LeafMgrT leafManager(result);
317  const size_t leafCount = leafManager.leafCount();
318 
319  if(leafCount == 0) return;
320 
321  // Count the number of active voxels in each leaf node.
322  hboost::scoped_array<VIndex> perLeafCount(new VIndex[leafCount]);
323  VIndex* perLeafCountPtr = perLeafCount.get();
324  leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
325 
326  // The starting index for each leaf node is the total number
327  // of active voxels in all preceding leaf nodes.
328  for (size_t i = 1; i < leafCount; ++i) {
329  perLeafCount[i] += perLeafCount[i - 1];
330  }
331 
332  // The last accumulated value should be the total of all active voxels.
333  assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
334 
335  // Parallelize over the leaf nodes of the tree, storing a unique index
336  // in each active voxel.
337  leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
338 }
339 
340 
341 template<typename TreeType>
342 inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
343 createIndexTree(const TreeType& tree)
344 {
345  typedef typename TreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
346 
347  // Construct an output tree with the same active voxel topology as the input tree.
348  const VIndex invalidIdx = -1;
349  typename VIdxTreeT::Ptr result(
350  new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
351 
352  // All active voxels are degrees of freedom, including voxels contained in active tiles.
353  result->voxelizeActiveTiles();
354 
355  populateIndexTree(*result);
356 
357  return result;
358 }
359 
360 
361 ////////////////////////////////////////
362 
363 
364 namespace internal {
365 
366 /// @brief Functor for use with LeafManager::foreach() to populate a vector
367 /// with the values of a tree's active voxels
368 template<typename VectorValueType, typename SourceTreeType>
370 {
371  typedef typename SourceTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
372  typedef typename VIdxTreeT::LeafNodeType VIdxLeafT;
373  typedef typename SourceTreeType::LeafNodeType LeafT;
374  typedef typename SourceTreeType::ValueType TreeValueT;
376 
377  const SourceTreeType* tree;
379 
380  CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
381 
382  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
383  {
384  VectorT& vec = *vector;
385  if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
386  // If a corresponding leaf node exists in the source tree,
387  // copy voxel values from the source node to the output vector.
388  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
389  vec[*it] = leaf->getValue(it.pos());
390  }
391  } else {
392  // If no corresponding leaf exists in the source tree,
393  // fill the vector with a uniform value.
394  const TreeValueT& value = tree->getValue(idxLeaf.origin());
395  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
396  vec[*it] = value;
397  }
398  }
399  }
400 };
401 
402 } // namespace internal
403 
404 
405 template<typename VectorValueType, typename SourceTreeType>
407 createVectorFromTree(const SourceTreeType& tree,
408  const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
409 {
410  typedef typename SourceTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
411  typedef tree::LeafManager<const VIdxTreeT> VIdxLeafMgrT;
412  typedef typename math::pcg::Vector<VectorValueType> VectorT;
413 
414  // Allocate the vector.
415  const size_t numVoxels = idxTree.activeVoxelCount();
416  typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
417 
418  // Parallelize over the leaf nodes of the index tree, filling the output vector
419  // with values from corresponding voxels of the source tree.
420  VIdxLeafMgrT leafManager(idxTree);
421  leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
422 
423  return result;
424 }
425 
426 
427 ////////////////////////////////////////
428 
429 
430 namespace internal {
431 
432 /// @brief Functor for use with LeafManager::foreach() to populate a tree
433 /// with values from a vector
434 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
436 {
437  typedef typename VIndexTreeType::template ValueConverter<TreeValueType>::Type OutTreeT;
438  typedef typename OutTreeT::LeafNodeType OutLeafT;
439  typedef typename VIndexTreeType::LeafNodeType VIdxLeafT;
441 
442  const VectorT* vector;
444 
445  CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
446 
447  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
448  {
449  const VectorT& vec = *vector;
450  OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
451  assert(leaf != NULL);
452  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
453  leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
454  }
455  }
456 };
457 
458 } // namespace internal
459 
460 
461 template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
462 inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
465  const VIndexTreeType& idxTree,
466  const TreeValueType& background)
467 {
468  typedef typename VIndexTreeType::template ValueConverter<TreeValueType>::Type OutTreeT;
469  typedef typename tree::LeafManager<const VIndexTreeType> VIdxLeafMgrT;
470 
471  // Construct an output tree with the same active voxel topology as the index tree.
472  typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
473  OutTreeT& tree = *result;
474 
475  // Parallelize over the leaf nodes of the index tree, populating voxels
476  // of the output tree with values from the input vector.
477  VIdxLeafMgrT leafManager(idxTree);
478  leafManager.foreach(
480 
481  return result;
482 }
483 
484 
485 ////////////////////////////////////////
486 
487 
488 namespace internal {
489 
490 /// Constant boundary condition functor
491 template<typename ValueType>
492 struct DirichletOp {
493  inline void operator()(
494  const Coord&, const Coord&, ValueType&, ValueType& diag) const { diag -= 1; }
495 };
496 
497 
498 /// Functor for use with LeafManager::foreach() to populate a sparse Laplacian matrix
499 template<typename BoolTreeType, typename BoundaryOp>
501 {
502  typedef typename BoolTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
503  typedef typename VIdxTreeT::LeafNodeType VIdxLeafT;
506 
509  const BoolTreeType* interiorMask;
510  const BoundaryOp boundaryOp;
512 
514  const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
515  laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
516 
517  void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
518  {
519  // Local accessors
521  typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
522 
523  Coord ijk;
524  VIndex column;
525  const ValueT diagonal = -6.f, offDiagonal = 1.f;
526 
527  // Loop over active voxels in this leaf.
528  for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
529  assert(it.getValue() > -1);
530  const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
531 
533 
534  ijk = it.getCoord();
535  if (interior.isValueOn(ijk)) {
536  // The current voxel is an interior voxel.
537  // All of its neighbors are in the solution domain.
538 
539  // -x direction
540  row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
541  // -y direction
542  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
543  // -z direction
544  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
545  // diagonal
546  row.setValue(rowNum, diagonal);
547  // +z direction
548  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
549  // +y direction
550  row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
551  // +x direction
552  row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
553 
554  } else {
555  // The current voxel is a boundary voxel.
556  // At least one of its neighbors is outside the solution domain.
557 
558  ValueT modifiedDiagonal = 0.f;
559 
560  // -x direction
561  if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
562  row.setValue(column, offDiagonal);
563  modifiedDiagonal -= 1;
564  } else {
565  boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
566  }
567  // -y direction
568  if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
569  row.setValue(column, offDiagonal);
570  modifiedDiagonal -= 1;
571  } else {
572  boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
573  }
574  // -z direction
575  if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
576  row.setValue(column, offDiagonal);
577  modifiedDiagonal -= 1;
578  } else {
579  boundaryOp(ijk, ijk.offsetBy(0, 0, -1), 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  // +y direction
589  if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
590  row.setValue(column, offDiagonal);
591  modifiedDiagonal -= 1;
592  } else {
593  boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
594  }
595  // +x direction
596  if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
597  row.setValue(column, offDiagonal);
598  modifiedDiagonal -= 1;
599  } else {
600  boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
601  }
602  // diagonal
603  row.setValue(rowNum, modifiedDiagonal);
604  }
605  } // end loop over voxels
606  }
607 };
608 
609 } // namespace internal
610 
611 
612 template<typename BoolTreeType>
614 createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
615  const BoolTreeType& interiorMask)
616 {
617  typedef LaplacianMatrix::ValueType ValueT;
619  static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
621  return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused);
622 }
623 
624 
625 template<typename BoolTreeType, typename BoundaryOp>
628  const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
629  const BoolTreeType& interiorMask,
630  const BoundaryOp& boundaryOp,
632 {
633  typedef typename BoolTreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
634  typedef typename tree::LeafManager<const VIdxTreeT> VIdxLeafMgrT;
635 
636  // The number of active voxels is the number of degrees of freedom.
637  const Index64 numDoF = idxTree.activeVoxelCount();
638 
639  // Construct the matrix.
640  LaplacianMatrix::Ptr laplacianPtr(
641  new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
642  LaplacianMatrix& laplacian = *laplacianPtr;
643 
644  // Populate the matrix using a second-order, 7-point CD stencil.
645  VIdxLeafMgrT idxLeafManager(idxTree);
647  laplacian, idxTree, interiorMask, boundaryOp, source));
648 
649  return laplacianPtr;
650 }
651 
652 
653 ////////////////////////////////////////
654 
655 
656 template<typename TreeType>
657 inline typename TreeType::Ptr
658 solve(const TreeType& inTree, math::pcg::State& state)
659 {
660  util::NullInterrupter interrupter;
661  return solve(inTree, state, interrupter);
662 }
663 
664 
665 template<typename TreeType, typename Interrupter>
666 inline typename TreeType::Ptr
667 solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter)
668 {
670  return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter);
671 }
672 
673 
674 template<typename TreeType, typename BoundaryOp, typename Interrupter>
675 inline typename TreeType::Ptr
676 solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
677  math::pcg::State& state, Interrupter& interrupter)
678 {
680  return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
681  inTree, boundaryOp, state, interrupter);
682 }
683 
684 
685 template<typename PreconditionerType, typename TreeType, typename BoundaryOp, typename Interrupter>
686 inline typename TreeType::Ptr
688  const BoundaryOp& boundaryOp, math::pcg::State& state, Interrupter& interrupter)
689 {
690 
691  return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(inTree /*source*/, inTree /*domain mask*/,
692  boundaryOp, state, interrupter);
693 }
694 
695 template<typename PreconditionerType, typename TreeType, typename DomainTreeType, typename BoundaryOp, typename Interrupter>
696 inline typename TreeType::Ptr
698  const DomainTreeType& domainMask,
699  const BoundaryOp& boundaryOp,
700  math::pcg::State& state, Interrupter& interrupter)
701 {
702 
703  typedef typename TreeType::ValueType TreeValueT;
704  typedef LaplacianMatrix::ValueType VecValueT;
705  typedef typename math::pcg::Vector<VecValueT> VectorT;
706  typedef typename TreeType::template ValueConverter<VIndex>::Type VIdxTreeT;
707  typedef typename TreeType::template ValueConverter<bool>::Type MaskTreeT;
708 
709  // 1. Create a mapping from active voxels of the input tree to elements of a vector.
710  typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
711 
712  // 2. Populate a vector with values from the input tree.
713  typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
714 
715  // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
716  typename MaskTreeT::Ptr interiorMask(
717  new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
718  tools::erodeVoxels(*interiorMask, /*iterations=*/1, tools::NN_FACE);
719 
720  // 4. Create the Laplacian matrix.
722  *idxTree, *interiorMask, boundaryOp, *b);
723 
724  // 5. Solve the Poisson equation.
725  laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
726  b->scale(-1.0);
727  typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
729  new PreconditionerType(*laplacian));
730  if (!precond->isValid()) {
731  precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
732  }
733 
734  state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
735 
736  // 6. Populate the output tree with values from the solution vector.
737  /// @todo if (state.success) ... ?
738  return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
739 }
740 
741 } // namespace poisson
742 } // namespace tools
743 } // namespace OPENVDB_VERSION_NAME
744 } // namespace openvdb
745 
746 #endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
747 
748 ///////////////////////////////////////////////////////////////////////////
749 //
750 // Copyright (c) 2012-2017 DreamWorks Animation LLC
751 //
752 // All rights reserved. This software is distributed under the
753 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
754 //
755 // Redistributions of source code must retain the above copyright
756 // and license notice and the following restrictions and disclaimer.
757 //
758 // * Neither the name of DreamWorks Animation nor the names of
759 // its contributors may be used to endorse or promote products derived
760 // from this software without specific prior written permission.
761 //
762 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
763 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
764 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
765 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
766 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
767 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
768 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
769 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
770 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
771 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
772 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
773 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
774 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
775 //
776 ///////////////////////////////////////////////////////////////////////////
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
GLint GLuint mask
Definition: glcorearb.h:123
TreeType::Ptr solveWithBoundaryConditions(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the value...
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:73
BoolTreeType::template ValueConverter< VIndex >::Type VIdxTreeT
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
VIndexTreeType::template ValueConverter< TreeValueType >::Type OutTreeT
LaplacianMatrix::Ptr createISLaplacianWithBoundaryConditions(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, const BoundaryOp &boundaryOp, typename math::pcg::Vector< LaplacianMatrix::ValueType > &source)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator with user-specified boundary ...
TreeType::Ptr solve(const TreeType &, math::pcg::State &)
Solve ∇2x = b for x, where b is a vector comprising the values of all of the active voxels in the inp...
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.
#define OPENVDB_VERSION_NAME
Definition: version.h:43
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
SourceTreeType::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:115
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1601
png_bytepp row
Definition: png.h:1836
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::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:875
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
SizeType setValue(SizeType column, const ValueType &value)
Set the value of the entry in the specified column.
LaplacianMatrix::Ptr createISLaplacian(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator using second-order finite dif...
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.
TreeType::Ptr solveWithBoundaryConditionsAndPreconditioner(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the value...
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:503
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Functor for use with LeafManager::foreach() to populate active leaf voxels with sequential indices...
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
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...
ISLaplacianOp(LaplacianMatrix &m, const VIdxTreeT &idx, const BoolTreeType &mask, const BoundaryOp &op, VectorT &src)
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
void operator()(const VIdxLeafT &idxLeaf, size_t) const
GLenum src
Definition: glcorearb.h:1792