HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PointAdvect.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 /// @author Ken Museth, D.J. Hill (openvdb port, added staggered grid support)
32 /// @file PointAdvect.h
33 ///
34 /// @brief Class PointAdvect advects points (with position) in a static velocity field
35 
36 #ifndef OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
37 #define OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
38 
39 #include <openvdb/openvdb.h>
40 #include <openvdb/math/Math.h> // min
41 #include <openvdb/Types.h> // Vec3 types and version number
42 #include <openvdb/Grid.h> // grid
44 #include "Interpolation.h" // sampling
45 #include "VelocityFields.h" // VelocityIntegrator
46 
47 #include <hboost/static_assert.hpp>
48 #include <tbb/blocked_range.h> // threading
49 #include <tbb/parallel_for.h> // threading
50 #include <tbb/task.h> // for cancel
51 
52 
53 namespace openvdb {
55 namespace OPENVDB_VERSION_NAME {
56 namespace tools {
57 
58 /// Class that holds a Vec3 grid, to be interpreted as the closest point to a constraint
59 /// surface. Supports a method to allow a point to be projected onto the closest point
60 /// on the constraint surface. Uses Caching.
61 template<typename CptGridT = Vec3fGrid>
63 {
64 public:
65  typedef CptGridT CptGridType;
66  typedef typename CptGridType::ConstAccessor CptAccessor;
67  typedef typename CptGridType::ValueType CptValueType;
68 
70  mCptIterations(0)
71  {
72  }
73  ClosestPointProjector(const CptGridType& cptGrid, int n):
74  mCptGrid(&cptGrid),
75  mCptAccessor(cptGrid.getAccessor()),
76  mCptIterations(n)
77  {
78  }
80  mCptGrid(other.mCptGrid),
81  mCptAccessor(mCptGrid->getAccessor()),
82  mCptIterations(other.mCptIterations)
83  {
84  }
85  void setConstraintIterations(unsigned int cptIterations) { mCptIterations = cptIterations; }
86  unsigned int numIterations() { return mCptIterations; }
87 
88  // point constraint
89  template <typename LocationType>
90  inline void projectToConstraintSurface(LocationType& W) const
91  {
92  /// Entries in the CPT tree are the closest point to the constraint surface.
93  /// The interpolation step in sample introduces error so that the result
94  /// of a single sample may not lie exactly on the surface. The iterations
95  /// in the loop exist to minimize this error.
96  CptValueType result(W[0], W[1],W[2]);
97  for (unsigned int i = 0; i < mCptIterations; ++i) {
98  const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
99  BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
100  }
101  W[0] = result[0];
102  W[1] = result[1];
103  W[2] = result[2];
104  }
105 
106 private:
107  const CptGridType* mCptGrid; // Closest-Point-Transform vector field
108  CptAccessor mCptAccessor;
109  unsigned int mCptIterations;
110 };// end of ClosestPointProjector class
111 
112 ////////////////////////////////////////
113 
114 
115 /// Performs passive or constrained advection of points in a velocity field
116 /// represented by an OpenVDB grid and an optional closest-point-transform (CPT)
117 /// represented in another OpenVDB grid. Note the CPT is assumed to be
118 /// in world coordinates and NOT index coordinates!
119 /// Supports both collocated velocity grids and staggered velocity grids
120 ///
121 /// The @c PointListT template argument refers to any class with the following
122 /// interface (e.g., std::vector<openvdb::Vec3f>):
123 /// @code
124 /// class PointList {
125 /// ...
126 /// public:
127 /// typedef internal_vector3_type value_type; // must support [] component access
128 /// openvdb::Index size() const; // number of points in list
129 /// value_type& operator[](int n); // world space position of nth point
130 /// };
131 /// @endcode
132 ///
133 /// @note All methods (except size) are assumed to be thread-safe and
134 /// the positions are returned as non-const references since the
135 /// advection method needs to modify them!
136 template<typename GridT = Vec3fGrid,
137  typename PointListT = std::vector<typename GridT::ValueType>,
138  bool StaggeredVelocity = false,
139  typename InterrupterType = util::NullInterrupter>
141 {
142 public:
143  typedef GridT GridType;
144  typedef PointListT PointListType;
147 
148  PointAdvect(const GridT& velGrid, InterrupterType* interrupter=NULL) :
149  mVelGrid(&velGrid),
150  mPoints(NULL),
151  mIntegrationOrder(1),
152  mThreaded(true),
153  mInterrupter(interrupter)
154  {
155  }
156  PointAdvect(const PointAdvect &other) :
157  mVelGrid(other.mVelGrid),
158  mPoints(other.mPoints),
159  mDt(other.mDt),
160  mAdvIterations(other.mAdvIterations),
161  mIntegrationOrder(other.mIntegrationOrder),
162  mThreaded(other.mThreaded),
163  mInterrupter(other.mInterrupter)
164  {
165  }
166  virtual ~PointAdvect()
167  {
168  }
169  /// If the order of the integration is set to zero no advection is performed
170  bool earlyOut() const { return (mIntegrationOrder==0);}
171  /// get & set
172  void setThreaded(bool threaded) { mThreaded = threaded; }
173  bool getThreaded() { return mThreaded; }
174  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
175 
176  /// Constrained advection of a list of points over a time = dt * advIterations
177  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
178  {
179  if (this->earlyOut()) return; // nothing to do!
180  mPoints = &points;
181  mDt = dt;
182  mAdvIterations = advIterations;
183 
184  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
185  if (mThreaded) {
186  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *this);
187  } else {
188  (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
189  }
190  if (mInterrupter) mInterrupter->end();
191  }
192 
193  /// Never call this method directly - it is use by TBB and has to be public!
194  void operator() (const tbb::blocked_range<size_t> &range) const
195  {
196  if (mInterrupter && mInterrupter->wasInterrupted()) {
197  tbb::task::self().cancel_group_execution();
198  }
199 
200  VelocityFieldIntegrator velField(*mVelGrid);
201  switch (mIntegrationOrder) {
202  case 1:
203  {
204  for (size_t n = range.begin(); n != range.end(); ++n) {
205  LocationType& X0 = (*mPoints)[n];
206  // loop over number of time steps
207  for (unsigned int i = 0; i < mAdvIterations; ++i) {
208  velField.template rungeKutta<1>(mDt, X0);
209  }
210  }
211  }
212  break;
213  case 2:
214  {
215  for (size_t n = range.begin(); n != range.end(); ++n) {
216  LocationType& X0 = (*mPoints)[n];
217  // loop over number of time steps
218  for (unsigned int i = 0; i < mAdvIterations; ++i) {
219  velField.template rungeKutta<2>(mDt, X0);
220  }
221  }
222  }
223  break;
224  case 3:
225  {
226  for (size_t n = range.begin(); n != range.end(); ++n) {
227  LocationType& X0 = (*mPoints)[n];
228  // loop over number of time steps
229  for (unsigned int i = 0; i < mAdvIterations; ++i) {
230  velField.template rungeKutta<3>(mDt, X0);
231  }
232  }
233  }
234  break;
235  case 4:
236  {
237  for (size_t n = range.begin(); n != range.end(); ++n) {
238  LocationType& X0 = (*mPoints)[n];
239  // loop over number of time steps
240  for (unsigned int i = 0; i < mAdvIterations; ++i) {
241  velField.template rungeKutta<4>(mDt, X0);
242  }
243  }
244  }
245  break;
246  }
247  }
248 
249 private:
250  // the velocity field
251  const GridType* mVelGrid;
252 
253  // vertex list of all the points
254  PointListT* mPoints;
255 
256  // time integration parameters
257  float mDt; // time step
258  unsigned int mAdvIterations; // number of time steps
259  unsigned int mIntegrationOrder;
260 
261  // operational parameters
262  bool mThreaded;
263  InterrupterType* mInterrupter;
264 
265 };//end of PointAdvect class
266 
267 
268 template<typename GridT = Vec3fGrid,
269  typename PointListT = std::vector<typename GridT::ValueType>,
270  bool StaggeredVelocity = false,
271  typename CptGridType = GridT,
272  typename InterrupterType = util::NullInterrupter>
274 {
275 public:
276  typedef GridT GridType;
280  typedef PointListT PointListType;
281 
283  const GridType& cptGrid, int cptn, InterrupterType* interrupter = NULL):
284  mVelGrid(&velGrid),
285  mCptGrid(&cptGrid),
286  mCptIter(cptn),
287  mInterrupter(interrupter)
288  {
289  }
291  mVelGrid(other.mVelGrid),
292  mCptGrid(other.mCptGrid),
293  mCptIter(other.mCptIter),
294  mPoints(other.mPoints),
295  mDt(other.mDt),
296  mAdvIterations(other.mAdvIterations),
297  mIntegrationOrder(other.mIntegrationOrder),
298  mThreaded(other.mThreaded),
299  mInterrupter(other.mInterrupter)
300  {
301  }
303 
304  void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
305  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
306 
307  void setThreaded(bool threaded) { mThreaded = threaded; }
308  bool getThreaded() { return mThreaded; }
309 
310  /// Constrained Advection a list of points over a time = dt * advIterations
311  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
312  {
313  mPoints = &points;
314  mDt = dt;
315 
316  if (mIntegrationOrder==0 && mCptIter == 0) {
317  return; // nothing to do!
318  }
319  (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
320 
321  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
322  const size_t N = mPoints->size();
323 
324  if (mThreaded) {
325  tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
326  } else {
327  (*this)(tbb::blocked_range<size_t>(0, N));
328  }
329  if (mInterrupter) mInterrupter->end();
330  }
331 
332 
333  /// Never call this method directly - it is use by TBB and has to be public!
334  void operator() (const tbb::blocked_range<size_t> &range) const
335  {
336  if (mInterrupter && mInterrupter->wasInterrupted()) {
337  tbb::task::self().cancel_group_execution();
338  }
339 
340  VelocityIntegratorType velField(*mVelGrid);
341  ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
342  switch (mIntegrationOrder) {
343  case 0://pure CPT projection
344  {
345  for (size_t n = range.begin(); n != range.end(); ++n) {
346  LocationType& X0 = (*mPoints)[n];
347  for (unsigned int i = 0; i < mAdvIterations; ++i) {
348  cptField.projectToConstraintSurface(X0);
349  }
350  }
351  }
352  break;
353  case 1://1'th order advection and CPT projection
354  {
355  for (size_t n = range.begin(); n != range.end(); ++n) {
356  LocationType& X0 = (*mPoints)[n];
357  for (unsigned int i = 0; i < mAdvIterations; ++i) {
358  velField.template rungeKutta<1>(mDt, X0);
359  cptField.projectToConstraintSurface(X0);
360  }
361  }
362  }
363  break;
364  case 2://2'nd order advection and CPT projection
365  {
366  for (size_t n = range.begin(); n != range.end(); ++n) {
367  LocationType& X0 = (*mPoints)[n];
368  for (unsigned int i = 0; i < mAdvIterations; ++i) {
369  velField.template rungeKutta<2>(mDt, X0);
370  cptField.projectToConstraintSurface(X0);
371  }
372  }
373  }
374  break;
375 
376  case 3://3'rd order advection and CPT projection
377  {
378  for (size_t n = range.begin(); n != range.end(); ++n) {
379  LocationType& X0 = (*mPoints)[n];
380  for (unsigned int i = 0; i < mAdvIterations; ++i) {
381  velField.template rungeKutta<3>(mDt, X0);
382  cptField.projectToConstraintSurface(X0);
383  }
384  }
385  }
386  break;
387  case 4://4'th order advection and CPT projection
388  {
389  for (size_t n = range.begin(); n != range.end(); ++n) {
390  LocationType& X0 = (*mPoints)[n];
391  for (unsigned int i = 0; i < mAdvIterations; ++i) {
392  velField.template rungeKutta<4>(mDt, X0);
393  cptField.projectToConstraintSurface(X0);
394  }
395  }
396  }
397  break;
398  }
399  }
400 
401 private:
402  const GridType* mVelGrid; // the velocity field
403  const GridType* mCptGrid;
404  int mCptIter;
405  PointListT* mPoints; // vertex list of all the points
406 
407  // time integration parameters
408  float mDt; // time step
409  unsigned int mAdvIterations; // number of time steps
410  unsigned int mIntegrationOrder; // order of Runge-Kutta integration
411  // operational parameters
412  bool mThreaded;
413  InterrupterType* mInterrupter;
414 };// end of ConstrainedPointAdvect class
415 
416 } // namespace tools
417 } // namespace OPENVDB_VERSION_NAME
418 } // namespace openvdb
419 
420 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
421 
422 // Copyright (c) 2012-2017 DreamWorks Animation LLC
423 // All rights reserved. This software is distributed under the
424 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
VelocityIntegrator< GridT, StaggeredVelocity > VelocityFieldIntegrator
Definition: PointAdvect.h:146
GLenum GLint * range
Definition: glcorearb.h:1924
void setThreaded(bool threaded)
get & set
Definition: PointAdvect.h:172
math::Vec3< Real > Vec3R
Definition: Types.h:75
void operator()(const tbb::blocked_range< size_t > &range) const
Never call this method directly - it is use by TBB and has to be public!
Definition: PointAdvect.h:194
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: PointAdvect.h:73
Dummy NOOP interrupter class defining interface.
png_uint_32 i
Definition: png.h:2877
uint64 value_type
Definition: GA_PrimCompat.h:29
png_infop png_charp png_int_32 * X0
Definition: png.h:2417
GLdouble n
Definition: glcorearb.h:2007
bool earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: PointAdvect.h:170
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=NULL)
Definition: PointAdvect.h:148
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Performs Runge-Kutta time integration of variable order in a static velocity field.
void operator()(const tbb::blocked_range< size_t > &range) const
Never call this method directly - it is use by TBB and has to be public!
Definition: PointAdvect.h:334
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
GLint location
Definition: glcorearb.h:804
ClosestPointProjector< CptGridType > ClosestPointProjectorType
Definition: PointAdvect.h:279
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=NULL)
Definition: PointAdvect.h:282
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained Advection a list of points over a time = dt * advIterations.
Definition: PointAdvect.h:311
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: PointAdvect.h:290
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained advection of a list of points over a time = dt * advIterations.
Definition: PointAdvect.h:177
GA_API const UT_StringHolder N
VelocityIntegrator< GridT, StaggeredVelocity > VelocityIntegratorType
Definition: PointAdvect.h:278
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
ClosestPointProjector(const ClosestPointProjector &other)
Definition: PointAdvect.h:79
void setConstraintIterations(unsigned int cptIterations)
Definition: PointAdvect.h:85