GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/PointAdvect.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 100 115 87.0%
Functions: 5 7 71.4%
Branches: 83 125 66.4%
Line Branch Exec Source
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
17 #include "openvdb/util/NullInterrupter.h"
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 {
29 OPENVDB_USE_VERSION_NAMESPACE
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>
37
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 class ClosestPointProjector
38 {
39 public:
40 using CptGridType = CptGridT;
41 using CptAccessor = typename CptGridType::ConstAccessor;
42 using CptValueType = typename CptGridType::ValueType;
43
44 ClosestPointProjector():
45 mCptIterations(0)
46 {
47 }
48 5 ClosestPointProjector(const CptGridType& cptGrid, int n):
49 mCptGrid(&cptGrid),
50 mCptAccessor(cptGrid.getAccessor()),
51 10 mCptIterations(n)
52 {
53 }
54 ClosestPointProjector(const ClosestPointProjector &other):
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 20 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 20 CptValueType result(W[0], W[1],W[2]);
72
2/2
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 20 times.
 120 for (unsigned int i = 0; i < mCptIterations; ++i) {
73 100 const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
74 100 BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
75 }
76 20 W[0] = result[0];
77 20 W[1] = result[1];
78 20 W[2] = result[2];
79 20 }
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>
115 class PointAdvect
116 {
117 public:
118 using GridType = GridT;
119 using PointListType = PointListT;
120 using LocationType = typename PointListT::value_type;
121 using VelocityFieldIntegrator = VelocityIntegrator<GridT, StaggeredVelocity>;
122
123 1 PointAdvect(const GridT& velGrid, InterrupterType* interrupter = nullptr):
124 mVelGrid(&velGrid),
125 mPoints(nullptr),
126 mIntegrationOrder(1),
127 mThreaded(true),
128
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mInterrupter(interrupter)
129 {
130 }
131 101 PointAdvect(const PointAdvect &other) :
132 101 mVelGrid(other.mVelGrid),
133 101 mPoints(other.mPoints),
134 101 mDt(other.mDt),
135 101 mAdvIterations(other.mAdvIterations),
136 101 mIntegrationOrder(other.mIntegrationOrder),
137 101 mThreaded(other.mThreaded),
138 5 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 4 bool earlyOut() const { return (mIntegrationOrder==0);}
146 /// get & set
147 void setThreaded(bool threaded) { mThreaded = threaded; }
148 bool getThreaded() { return mThreaded; }
149
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
150
151 /// Constrained advection of a list of points over a time = dt * advIterations
152 4 void advect(PointListT& points, float dt, unsigned int advIterations = 1)
153 {
154
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (this->earlyOut()) return; // nothing to do!
155 4 mPoints = &points;
156 4 mDt = dt;
157 4 mAdvIterations = advIterations;
158
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
160
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (mThreaded) {
161 4 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (mInterrupter) mInterrupter->end();
166 }
167
168 /// Never call this method directly - it is use by TBB and has to be public!
169 520 void operator() (const tbb::blocked_range<size_t> &range) const
170 {
171
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 520 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 520 if (mInterrupter && mInterrupter->wasInterrupted()) {
172 ✗ thread::cancelGroupExecution();
173 }
174
175 520 VelocityFieldIntegrator velField(*mVelGrid);
176
4/5
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 128 times.
✓ Branch 2 taken 136 times.
✓ Branch 3 taken 128 times.
✗ Branch 4 not taken.
 520 switch (mIntegrationOrder) {
177 case 1:
178 {
179
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 128 times.
 2000128 for (size_t n = range.begin(); n != range.end(); ++n) {
180 2000000 LocationType& X0 = (*mPoints)[n];
181 // loop over number of time steps
182
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 2000000 times.
 4000000 for (unsigned int i = 0; i < mAdvIterations; ++i) {
183
1/2
✓ Branch 1 taken 2000000 times.
✗ Branch 2 not taken.
 2000000 velField.template rungeKutta<1>(mDt, X0);
184 }
185 }
186 }
187 break;
188 case 2:
189 {
190
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 128 times.
 2000128 for (size_t n = range.begin(); n != range.end(); ++n) {
191 2000000 LocationType& X0 = (*mPoints)[n];
192 // loop over number of time steps
193
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 2000000 times.
 4000000 for (unsigned int i = 0; i < mAdvIterations; ++i) {
194
1/2
✓ Branch 1 taken 2000000 times.
✗ Branch 2 not taken.
 2000000 velField.template rungeKutta<2>(mDt, X0);
195 }
196 }
197 }
198 break;
199 case 3:
200 {
201
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 136 times.
 2000136 for (size_t n = range.begin(); n != range.end(); ++n) {
202 2000000 LocationType& X0 = (*mPoints)[n];
203 // loop over number of time steps
204
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 2000000 times.
 4000000 for (unsigned int i = 0; i < mAdvIterations; ++i) {
205
1/2
✓ Branch 1 taken 2000000 times.
✗ Branch 2 not taken.
 2000000 velField.template rungeKutta<3>(mDt, X0);
206 }
207 }
208 }
209 break;
210 case 4:
211 {
212
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 128 times.
 2000128 for (size_t n = range.begin(); n != range.end(); ++n) {
213 2000000 LocationType& X0 = (*mPoints)[n];
214 // loop over number of time steps
215
2/2
✓ Branch 0 taken 2000000 times.
✓ Branch 1 taken 2000000 times.
 4000000 for (unsigned int i = 0; i < mAdvIterations; ++i) {
216
1/2
✓ Branch 1 taken 2000000 times.
✗ Branch 2 not taken.
 2000000 velField.template rungeKutta<4>(mDt, X0);
217 }
218 }
219 }
220 break;
221 }
222 520 }
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>
248 class ConstrainedPointAdvect
249 {
250 public:
251 using GridType = GridT;
252 using LocationType = typename PointListT::value_type;
253 using VelocityIntegratorType = VelocityIntegrator<GridT, StaggeredVelocity>;
254 using ClosestPointProjectorType = ClosestPointProjector<CptGridType>;
255 using PointListType = PointListT;
256
257 1 ConstrainedPointAdvect(const GridType& velGrid,
258 const GridType& cptGrid, int cptn, InterrupterType* interrupter = nullptr):
259 mVelGrid(&velGrid),
260 mCptGrid(&cptGrid),
261 mCptIter(cptn),
262 1 mInterrupter(interrupter)
263 {
264 }
265 ✗ ConstrainedPointAdvect(const ConstrainedPointAdvect& other):
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 }
277
1/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 virtual ~ConstrainedPointAdvect(){}
278
279 1 void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
280
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
281
282
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 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 5 void advect(PointListT& points, float dt, unsigned int advIterations = 1)
287 {
288 5 mPoints = &points;
289 5 mDt = dt;
290
291
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 5 if (mIntegrationOrder==0 && mCptIter == 0) {
292 return; // nothing to do!
293 }
294
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
295
296
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
297
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 const size_t N = mPoints->size();
298
299
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (mThreaded) {
300 ✗ tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
301 } else {
302 5 (*this)(tbb::blocked_range<size_t>(0, N));
303 }
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (mInterrupter) mInterrupter->end();
305 }
306
307
308 /// Never call this method directly - it is use by TBB and has to be public!
309 5 void operator() (const tbb::blocked_range<size_t> &range) const
310 {
311
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 5 if (mInterrupter && mInterrupter->wasInterrupted()) {
312 ✗ thread::cancelGroupExecution();
313 }
314
315 5 VelocityIntegratorType velField(*mVelGrid);
316
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
317
5/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 5 switch (mIntegrationOrder) {
318 case 0://pure CPT projection
319 {
320
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (size_t n = range.begin(); n != range.end(); ++n) {
321 4 LocationType& X0 = (*mPoints)[n];
322
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (unsigned int i = 0; i < mAdvIterations; ++i) {
323
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cptField.projectToConstraintSurface(X0);
324 }
325 }
326 }
327 break;
328 case 1://1'th order advection and CPT projection
329 {
330
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (size_t n = range.begin(); n != range.end(); ++n) {
331 4 LocationType& X0 = (*mPoints)[n];
332
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (unsigned int i = 0; i < mAdvIterations; ++i) {
333
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 velField.template rungeKutta<1>(mDt, X0);
334
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cptField.projectToConstraintSurface(X0);
335 }
336 }
337 }
338 break;
339 case 2://2'nd order advection and CPT projection
340 {
341
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (size_t n = range.begin(); n != range.end(); ++n) {
342 4 LocationType& X0 = (*mPoints)[n];
343
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (unsigned int i = 0; i < mAdvIterations; ++i) {
344
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 velField.template rungeKutta<2>(mDt, X0);
345
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cptField.projectToConstraintSurface(X0);
346 }
347 }
348 }
349 break;
350
351 case 3://3'rd order advection and CPT projection
352 {
353
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (size_t n = range.begin(); n != range.end(); ++n) {
354 4 LocationType& X0 = (*mPoints)[n];
355
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (unsigned int i = 0; i < mAdvIterations; ++i) {
356
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 velField.template rungeKutta<3>(mDt, X0);
357
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cptField.projectToConstraintSurface(X0);
358 }
359 }
360 }
361 break;
362 case 4://4'th order advection and CPT projection
363 {
364
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (size_t n = range.begin(); n != range.end(); ++n) {
365 4 LocationType& X0 = (*mPoints)[n];
366
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (unsigned int i = 0; i < mAdvIterations; ++i) {
367
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 velField.template rungeKutta<4>(mDt, X0);
368
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cptField.projectToConstraintSurface(X0);
369 }
370 }
371 }
372 break;
373 }
374 5 }
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
396