GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/math/DDA.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 81 89 91.0%
Functions: 26 31 83.9%
Branches: 104 194 53.6%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file DDA.h
5 ///
6 /// @author Ken Museth
7 ///
8 /// @brief Digital Differential Analyzers specialized for VDB.
9
10 #ifndef OPENVDB_MATH_DDA_HAS_BEEN_INCLUDED
11 #define OPENVDB_MATH_DDA_HAS_BEEN_INCLUDED
12
13 #include "Coord.h"
14 #include "Math.h"
15 #include "Vec3.h"
16 #include <openvdb/Types.h>
17 #include <iostream> // for std::ostream
18 #include <limits> // for std::numeric_limits<Type>::max()
19
20 namespace openvdb {
21 OPENVDB_USE_VERSION_NAMESPACE
22 namespace OPENVDB_VERSION_NAME {
23 namespace math {
24
25 /// @brief A Digital Differential Analyzer specialized for OpenVDB grids
26 /// @note Conceptually similar to Bresenham's line algorithm applied
27 /// to a 3D Ray intersecting OpenVDB nodes or voxels. Log2Dim = 0
28 /// corresponds to a voxel and Log2Dim a tree node of size 2^Log2Dim.
29 ///
30 /// @note The Ray template class is expected to have the following
31 /// methods: test(time), t0(), t1(), invDir(), and operator()(time).
32 /// See the example Ray class above for their definition.
33 template<typename RayT, Index Log2Dim = 0>
34 class DDA
35 {
36 public:
37 using RealType = typename RayT::RealType;
38 using RealT = RealType;
39 using Vec3Type = typename RayT::Vec3Type;
40 using Vec3T = Vec3Type;
41
42 /// @brief uninitialized constructor
43 122 DDA() {}
44
45 1585747 DDA(const RayT& ray) { this->init(ray); }
46
47 DDA(const RayT& ray, RealT startTime) { this->init(ray, startTime); }
48
49 DDA(const RayT& ray, RealT startTime, RealT maxTime) { this->init(ray, startTime, maxTime); }
50
51 5336366 inline void init(const RayT& ray, RealT startTime, RealT maxTime)
52 {
53
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2668183 times.
 5336366 assert(startTime <= maxTime);
54 static const int DIM = 1 << Log2Dim;
55 5336366 mT0 = startTime;
56 5336366 mT1 = maxTime;
57 5336366 const Vec3T &pos = ray(mT0), &dir = ray.dir(), &inv = ray.invDir();
58 5336366 mVoxel = Coord::floor(pos) & (~(DIM-1));
59
2/2
✓ Branch 0 taken 8004549 times.
✓ Branch 1 taken 2668183 times.
 21345464 for (int axis = 0; axis < 3; ++axis) {
60
2/2
✓ Branch 0 taken 5336346 times.
✓ Branch 1 taken 2668203 times.
 16009098 if (math::isZero(dir[axis])) {//handles dir = +/- 0
61
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5336346 times.
 10672692 mStep[axis] = 0;//dummy value
62 10672692 mNext[axis] = std::numeric_limits<RealT>::max();//i.e. disabled!
63 10672692 mDelta[axis] = std::numeric_limits<RealT>::max();//dummy value
64
2/2
✓ Branch 0 taken 2668152 times.
✓ Branch 1 taken 51 times.
 5336406 } else if (inv[axis] > 0) {
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2668152 times.
 5336304 mStep[axis] = DIM;
66 5336304 mNext[axis] = mT0 + (mVoxel[axis] + DIM - pos[axis]) * inv[axis];
67 5336304 mDelta[axis] = mStep[axis] * inv[axis];
68 } else {
69
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
 102 mStep[axis] = -DIM;
70 102 mNext[axis] = mT0 + (mVoxel[axis] - pos[axis]) * inv[axis];
71 102 mDelta[axis] = mStep[axis] * inv[axis];
72 }
73 }
74 5336366 }
75
76 2668183 inline void init(const RayT& ray) { this->init(ray, ray.t0(), ray.t1()); }
77
78 inline void init(const RayT& ray, RealT startTime) { this->init(ray, startTime, ray.t1()); }
79
80 /// @brief Increment the voxel index to next intersected voxel or node
81 /// and returns true if the step in time does not exceed maxTime.
82 22264409 inline bool step()
83 {
84
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11774736 times.
 22264409 const int stepAxis = static_cast<int>(math::MinIndex(mNext));
85 22264409 mT0 = mNext[stepAxis];
86 22264409 mNext[stepAxis] += mDelta[stepAxis];
87
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11774736 times.
 22264409 mVoxel[stepAxis] += mStep[stepAxis];
88 22264409 return mT0 <= mT1;
89 }
90
91 /// @brief Return the index coordinates of the next node or voxel
92 /// intersected by the ray. If Log2Dim = 0 the return value is the
93 /// actual signed coordinate of the voxel, else it is the origin
94 /// of the corresponding VDB tree node or tile.
95 /// @note Incurs no computational overhead.
96
7/28
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 1285039 times.
 1285174 inline const Coord& voxel() const { return mVoxel; }
97
98 /// @brief Return the time (parameterized along the Ray) of the
99 /// first hit of a tree node of size 2^Log2Dim.
100 /// @details This value is initialized to startTime or ray.t0()
101 /// depending on the constructor used.
102 /// @note Incurs no computational overhead.
103
4/8
✓ Branch 6 taken 60 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 60 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 60 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 54 times.
✗ Branch 16 not taken.
 1822707 inline RealType time() const { return mT0; }
104
105 /// @brief Return the maximum time (parameterized along the Ray).
106 11 inline RealType maxTime() const { return mT1; }
107
108 /// @brief Return the time (parameterized along the Ray) of the
109 /// second (i.e. next) hit of a tree node of size 2^Log2Dim.
110 /// @note Incurs a (small) computational overhead.
111
31/56
✓ Branch 0 taken 3642762 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 537343 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 537323 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 8 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 515868 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 515872 times.
✓ Branch 13 taken 24 times.
✓ Branch 14 taken 17 times.
✓ Branch 15 taken 17 times.
✓ Branch 16 taken 5 times.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 266257 times.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 266256 times.
✓ Branch 21 taken 21 times.
✓ Branch 22 taken 60 times.
✓ Branch 23 taken 38 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 5 times.
✓ Branch 26 taken 5 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 4 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 5 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✓ Branch 53 taken 1285039 times.
✓ Branch 54 taken 1 times.
✓ Branch 55 taken 1285038 times.
 8852071 inline RealType next() const { return math::Min(mT1, mNext[0], mNext[1], mNext[2]); }
112
113 /// @brief Print information about this DDA for debugging.
114 /// @param os a stream to which to write textual information.
115 void print(std::ostream& os = std::cout) const
116 {
117 os << "Dim=" << (1<<Log2Dim) << " time=" << mT0 << " next()="
118 << this->next() << " voxel=" << mVoxel << " next=" << mNext
119 << " delta=" << mDelta << " step=" << mStep << std::endl;
120 }
121
122 private:
123 RealT mT0, mT1;
124 Coord mVoxel, mStep;
125 Vec3T mDelta, mNext;
126 }; // class DDA
127
128 /// @brief Output streaming of the Ray class.
129 /// @note Primarily intended for debugging.
130 template<typename RayT, Index Log2Dim>
131 inline std::ostream& operator<<(std::ostream& os, const DDA<RayT, Log2Dim>& dda)
132 {
133 os << "Dim=" << (1<<Log2Dim) << " time=" << dda.time()
134 << " next()=" << dda.next() << " voxel=" << dda.voxel();
135 return os;
136 }
137
138 /////////////////////////////////////////// LevelSetHDDA ////////////////////////////////////////////
139
140
141 /// @brief Helper class that implements Hierarchical Digital Differential Analyzers
142 /// and is specialized for ray intersections with level sets
143 template<typename TreeT, int NodeLevel>
144 struct LevelSetHDDA
145 {
146 using ChainT = typename TreeT::RootNodeType::NodeChainType;
147 using NodeT = typename ChainT::template Get<NodeLevel>;
148
149 template <typename TesterT>
150 2096788 static bool test(TesterT& tester)
151 {
152 math::DDA<typename TesterT::RayT, NodeT::TOTAL> dda(tester.ray());
153
2/2
✓ Branch 1 taken 6621679 times.
✓ Branch 2 taken 430796 times.
 14104950 do {
154
2/2
✓ Branch 0 taken 1319460 times.
✓ Branch 1 taken 6350613 times.
 15340146 if (tester.template hasNode<NodeT>(dda.voxel())) {
155 tester.setRange(dda.time(), dda.next());
156
2/2
✓ Branch 1 taken 701862 times.
✓ Branch 2 taken 617598 times.
 2638920 if (LevelSetHDDA<TreeT, NodeLevel-1>::test(tester)) return true;
157 }
158 } while(dda.step());
159 return false;
160 }
161 };
162
163 /// @brief Specialization of Hierarchical Digital Differential Analyzer
164 /// class that intersects a ray against the voxels of a level set
165 template<typename TreeT>
166 struct LevelSetHDDA<TreeT, -1>
167 {
168 template <typename TesterT>
169 537330 static bool test(TesterT& tester)
170 {
171 math::DDA<typename TesterT::RayT, 0> dda(tester.ray());
172 tester.init(dda.time());
173
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 3642810 do { if (tester(dda.voxel(), dda.next())) return true; } while(dda.step());
174 return false;
175 }
176 };
177
178 //////////////////////////////////////////// VolumeHDDA /////////////////////////////////////////////
179
180 /// @brief Helper class that implements Hierarchical Digital Differential Analyzers
181 /// for ray intersections against a generic volume.
182 ///
183 /// @details The template argument ChildNodeLevel specifies the entry
184 /// upper node level used for the hierarchical ray-marching. The final
185 /// lowest level is always the leaf node level, i.e. not the voxel level!
186 template <typename TreeT, typename RayT, int ChildNodeLevel>
187 class VolumeHDDA
188 {
189 public:
190
191 using ChainT = typename TreeT::RootNodeType::NodeChainType;
192 using NodeT = typename ChainT::template Get<ChildNodeLevel>;
193 using TimeSpanT = typename RayT::TimeSpan;
194
195 13 VolumeHDDA() {}
196
197 template <typename AccessorT>
198
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 8 times.
 22 TimeSpanT march(RayT& ray, AccessorT &acc)
199 {
200 TimeSpanT t(-1, -1);
201
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 8 times.
 22 if (ray.valid()) this->march(ray, acc, t);
202 22 return t;
203 }
204
205 /// ListType is a list of RayType::TimeSpan and is required to
206 /// have the two methods: clear() and push_back(). Thus, it could
207 /// be std::vector<typename RayType::TimeSpan> or
208 /// std::deque<typename RayType::TimeSpan>.
209 template <typename AccessorT, typename ListT>
210
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 4 void hits(RayT& ray, AccessorT &acc, ListT& times)
211 {
212 TimeSpanT t(-1,-1);
213 2 times.clear();
214 4 this->hits(ray, acc, times, t);
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if (t.valid()) times.push_back(t);
216 4 }
217
218 private:
219
220 friend class VolumeHDDA<TreeT, RayT, ChildNodeLevel+1>;
221
222 template <typename AccessorT>
223 60 bool march(RayT& ray, AccessorT &acc, TimeSpanT& t)
224 {
225 60 mDDA.init(ray);
226
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 4 times.
 26 do {
227
2/2
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 2 times.
 78 if (acc.template probeConstNode<NodeT>(mDDA.voxel()) != nullptr) {//child node
228 ray.setTimes(mDDA.time(), mDDA.next());
229
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 26 times.
 74 if (mHDDA.march(ray, acc, t)) return true;//terminate
230
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 4 } else if (acc.isValueOn(mDDA.voxel())) {//hit an active tile
231 ✗ if (t.t0<0) t.t0 = mDDA.time();//this is the first hit so set t0
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 } else if (t.t0>=0) {//hit an inactive tile after hitting active values
233 ✗ t.t1 = mDDA.time();//set end of active ray segment
234 ✗ if (t.valid()) return true;//terminate
235 t.set(-1, -1);//reset to an empty and invalid time-span
236 }
237 } while (mDDA.step());
238
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 8 if (t.t0>=0) t.t1 = mDDA.maxTime();
239 return false;
240 }
241
242 /// ListType is a list of RayType::TimeSpan and is required to
243 /// have the two methods: clear() and push_back(). Thus, it could
244 /// be std::vector<typename RayType::TimeSpan> or
245 /// std::deque<typename RayType::TimeSpan>.
246 template <typename AccessorT, typename ListT>
247 8 void hits(RayT& ray, AccessorT &acc, ListT& times, TimeSpanT& t)
248 {
249 8 mDDA.init(ray);
250
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 8 do {
251
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 if (acc.template probeConstNode<NodeT>(mDDA.voxel()) != nullptr) {//child node
252 ray.setTimes(mDDA.time(), mDDA.next());
253 8 mHDDA.hits(ray, acc, times, t);
254 ✗ } else if (acc.isValueOn(mDDA.voxel())) {//hit an active tile
255 ✗ if (t.t0<0) t.t0 = mDDA.time();//this is the first hit so set t0
256 ✗ } else if (t.t0>=0) {//hit an inactive tile after hitting active values
257 ✗ t.t1 = mDDA.time();//set end of active ray segment
258 ✗ if (t.valid()) times.push_back(t);
259 t.set(-1,-1);//reset to an empty and invalid time-span
260 }
261 } while (mDDA.step());
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 8 if (t.t0>=0) t.t1 = mDDA.maxTime();
263 8 }
264
265 math::DDA<RayT, NodeT::TOTAL> mDDA;
266 VolumeHDDA<TreeT, RayT, ChildNodeLevel-1> mHDDA;
267 };
268
269 /// @brief Specialization of Hierarchical Digital Differential Analyzer
270 /// class that intersects against the leafs or tiles of a generic volume.
271 template <typename TreeT, typename RayT>
272 class VolumeHDDA<TreeT, RayT, 0>
273 {
274 public:
275
276 using LeafT = typename TreeT::LeafNodeType;
277 using TimeSpanT = typename RayT::TimeSpan;
278
279 VolumeHDDA() {}
280
281 template <typename AccessorT>
282 TimeSpanT march(RayT& ray, AccessorT &acc)
283 {
284 TimeSpanT t(-1, -1);
285 if (ray.valid()) this->march(ray, acc, t);
286 return t;
287 }
288
289 template <typename AccessorT, typename ListT>
290 void hits(RayT& ray, AccessorT &acc, ListT& times)
291 {
292 TimeSpanT t(-1,-1);
293 times.clear();
294 this->hits(ray, acc, times, t);
295 if (t.valid()) times.push_back(t);
296 }
297
298 private:
299
300 friend class VolumeHDDA<TreeT, RayT, 1>;
301
302 template <typename AccessorT>
303 21 bool march(RayT& ray, AccessorT &acc, TimeSpanT& t)
304 {
305 21 mDDA.init(ray);
306
2/2
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 8 times.
 43 do {
307
4/4
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 36 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 19 times.
 56 if (acc.template probeConstNode<LeafT>(mDDA.voxel()) ||
308 acc.isValueOn(mDDA.voxel())) {//hit a leaf or an active tile
309
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 24 times.
 37 if (t.t0<0) t.t0 = mDDA.time();//this is the first hit
310
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 6 times.
 19 } else if (t.t0>=0) {//hit an inactive tile after hitting active values
311
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 t.t1 = mDDA.time();//set end of active ray segment
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 if (t.valid()) return true;//terminate
313 t.set(-1, -1);//reset to an empty and invalid time-span
314 }
315 } while (mDDA.step());
316
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 if (t.t0>=0) t.t1 = mDDA.maxTime();
317 return false;
318 }
319
320 template <typename AccessorT, typename ListT>
321 4 void hits(RayT& ray, AccessorT &acc, ListT& times, TimeSpanT& t)
322 {
323 4 mDDA.init(ray);
324
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
 24 do {
325
4/4
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 4 times.
 24 if (acc.template probeConstNode<LeafT>(mDDA.voxel()) ||
326 acc.isValueOn(mDDA.voxel())) {//hit a leaf or an active tile
327
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 16 if (t.t0<0) t.t0 = mDDA.time();//this is the first hit
328
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 } else if (t.t0>=0) {//hit an inactive tile after hitting active values
329
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 t.t1 = mDDA.time();//set end of active ray segment
330
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 if (t.valid()) times.push_back(t);
331 t.set(-1, -1);//reset to an empty and invalid time-span
332 }
333 } while (mDDA.step());
334
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if (t.t0>=0) t.t1 = mDDA.maxTime();
335 4 }
336 math::DDA<RayT, LeafT::TOTAL> mDDA;
337 };
338
339 } // namespace math
340 } // namespace OPENVDB_VERSION_NAME
341 } // namespace openvdb
342
343 #endif // OPENVDB_MATH_DDA_HAS_BEEN_INCLUDED
344