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