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