HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RayIntersector.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 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 /// @file RayIntersector.h
32 ///
33 /// @author Ken Museth
34 ///
35 /// @brief Accelerated intersection of a ray with a narrow-band level
36 /// set or a generic (e.g. density) volume. This will of course be
37 /// useful for respectively surface and volume rendering.
38 ///
39 /// @details This file defines two main classes,
40 /// LevelSetRayIntersector and VolumeRayIntersector, as well as the
41 /// three support classes LevelSetHDDA, VolumeHDDA and LinearSearchImpl.
42 /// The LevelSetRayIntersector is templated on the LinearSearchImpl class
43 /// and calls instances of the LevelSetHDDA class. The reason to split
44 /// level set ray intersection into three classes is twofold. First
45 /// LevelSetRayIntersector defines the public API for client code and
46 /// LinearSearchImpl defines the actual algorithm used for the
47 /// ray level-set intersection. In other words this design will allow
48 /// for the public API to be fixed while the intersection algorithm
49 /// can change without resolving to (slow) virtual methods. Second,
50 /// LevelSetHDDA, which implements a hierarchical Differential Digital
51 /// Analyzer, relies on partial template specialization, so it has to
52 /// be a standalone class (as opposed to a member class of
53 /// LevelSetRayIntersector). The VolumeRayIntersector is conceptually
54 /// much simpler than the LevelSetRayIntersector, and hence it only
55 /// depends on VolumeHDDA that implements the hierarchical
56 /// Differential Digital Analyzer.
57 
58 
59 #ifndef OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
60 #define OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
61 
62 #include <openvdb/math/DDA.h>
63 #include <openvdb/math/Math.h>
64 #include <openvdb/math/Ray.h>
65 #include <openvdb/math/Stencils.h>
66 #include <openvdb/Grid.h>
67 #include <openvdb/Types.h>
68 #include "Morphology.h"
69 #include <iostream>
70 #include <type_traits>
71 
72 
73 namespace openvdb {
75 namespace OPENVDB_VERSION_NAME {
76 namespace tools {
77 
78 // Helper class that implements the actual search of the zero-crossing
79 // of the level set along the direction of a ray. This particular
80 // implementation uses iterative linear search.
81 template<typename GridT, int Iterations = 0, typename RealT = double>
83 
84 
85 ///////////////////////////////////// LevelSetRayIntersector /////////////////////////////////////
86 
87 
88 /// @brief This class provides the public API for intersecting a ray
89 /// with a narrow-band level set.
90 ///
91 /// @details It wraps a SearchImplT with a simple public API and
92 /// performs the actual hierarchical tree node and voxel traversal.
93 ///
94 /// @warning Use the (default) copy-constructor to make sure each
95 /// computational thread has their own instance of this class. This is
96 /// important since the SearchImplT contains a ValueAccessor that is
97 /// not thread-safe. However copying is very efficient.
98 ///
99 /// @see tools/RayTracer.h for examples of intended usage.
100 ///
101 /// @todo Add TrilinearSearchImpl, as an alternative to LinearSearchImpl,
102 /// that performs analytical 3D trilinear intersection tests, i.e., solves
103 /// cubic equations. This is slower but also more accurate than the 1D
104 /// linear interpolation in LinearSearchImpl.
105 template<typename GridT,
106  typename SearchImplT = LinearSearchImpl<GridT>,
107  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
108  typename RayT = math::Ray<Real> >
110 {
111 public:
112  using GridType = GridT;
113  using RayType = RayT;
114  using RealType = typename RayT::RealType;
115  using Vec3Type = typename RayT::Vec3T;
116  using ValueT = typename GridT::ValueType;
117  using TreeT = typename GridT::TreeType;
118 
119  static_assert(NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range");
121  "level set grids must have scalar, floating-point value types");
122 
123  /// @brief Constructor
124  /// @param grid level set grid to intersect rays against.
125  /// @param isoValue optional iso-value for the ray-intersection.
126  LevelSetRayIntersector(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>())
127  : mTester(grid, isoValue)
128  {
129  if (!grid.hasUniformVoxels() ) {
130  OPENVDB_THROW(RuntimeError,
131  "LevelSetRayIntersector only supports uniform voxels!");
132  }
133  if (grid.getGridClass() != GRID_LEVEL_SET) {
134  OPENVDB_THROW(RuntimeError,
135  "LevelSetRayIntersector only supports level sets!"
136  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
137  }
138  }
139 
140  /// @brief Return the iso-value used for ray-intersections
141  const ValueT& getIsoValue() const { return mTester.getIsoValue(); }
142 
143  /// @brief Return @c true if the index-space ray intersects the level set.
144  /// @param iRay ray represented in index space.
145  bool intersectsIS(const RayType& iRay) const
146  {
147  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
149  }
150 
151  /// @brief Return @c true if the index-space ray intersects the level set
152  /// @param iRay ray represented in index space.
153  /// @param iTime if an intersection was found it is assigned the time of the
154  /// intersection along the index ray.
155  bool intersectsIS(const RayType& iRay, RealType &iTime) const
156  {
157  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
158  iTime = mTester.getIndexTime();
160  }
161 
162  /// @brief Return @c true if the index-space ray intersects the level set.
163  /// @param iRay ray represented in index space.
164  /// @param xyz if an intersection was found it is assigned the
165  /// intersection point in index space, otherwise it is unchanged.
166  bool intersectsIS(const RayType& iRay, Vec3Type& xyz) const
167  {
168  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
169  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
170  mTester.getIndexPos(xyz);
171  return true;
172  }
173 
174  /// @brief Return @c true if the index-space ray intersects the level set.
175  /// @param iRay ray represented in index space.
176  /// @param xyz if an intersection was found it is assigned the
177  /// intersection point in index space, otherwise it is unchanged.
178  /// @param iTime if an intersection was found it is assigned the time of the
179  /// intersection along the index ray.
180  bool intersectsIS(const RayType& iRay, Vec3Type& xyz, RealType &iTime) const
181  {
182  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
183  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
184  mTester.getIndexPos(xyz);
185  iTime = mTester.getIndexTime();
186  return true;
187  }
188 
189  /// @brief Return @c true if the world-space ray intersects the level set.
190  /// @param wRay ray represented in world space.
191  bool intersectsWS(const RayType& wRay) const
192  {
193  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
195  }
196 
197  /// @brief Return @c true if the world-space ray intersects the level set.
198  /// @param wRay ray represented in world space.
199  /// @param wTime if an intersection was found it is assigned the time of the
200  /// intersection along the world ray.
201  bool intersectsWS(const RayType& wRay, RealType &wTime) const
202  {
203  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
204  wTime = mTester.getWorldTime();
206  }
207 
208  /// @brief Return @c true if the world-space ray intersects the level set.
209  /// @param wRay ray represented in world space.
210  /// @param world if an intersection was found it is assigned the
211  /// intersection point in world space, otherwise it is unchanged
212  bool intersectsWS(const RayType& wRay, Vec3Type& world) const
213  {
214  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
215  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
216  mTester.getWorldPos(world);
217  return true;
218  }
219 
220  /// @brief Return @c true if the world-space ray intersects the level set.
221  /// @param wRay ray represented in world space.
222  /// @param world if an intersection was found it is assigned the
223  /// intersection point in world space, otherwise it is unchanged.
224  /// @param wTime if an intersection was found it is assigned the time of the
225  /// intersection along the world ray.
226  bool intersectsWS(const RayType& wRay, Vec3Type& world, RealType &wTime) const
227  {
228  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
229  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
230  mTester.getWorldPos(world);
231  wTime = mTester.getWorldTime();
232  return true;
233  }
234 
235  /// @brief Return @c true if the world-space ray intersects the level set.
236  /// @param wRay ray represented in world space.
237  /// @param world if an intersection was found it is assigned the
238  /// intersection point in world space, otherwise it is unchanged.
239  /// @param normal if an intersection was found it is assigned the normal
240  /// of the level set surface in world space, otherwise it is unchanged.
241  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal) const
242  {
243  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
244  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
245  mTester.getWorldPosAndNml(world, normal);
246  return true;
247  }
248 
249  /// @brief Return @c true if the world-space ray intersects the level set.
250  /// @param wRay ray represented in world space.
251  /// @param world if an intersection was found it is assigned the
252  /// intersection point in world space, otherwise it is unchanged.
253  /// @param normal if an intersection was found it is assigned the normal
254  /// of the level set surface in world space, otherwise it is unchanged.
255  /// @param wTime if an intersection was found it is assigned the time of the
256  /// intersection along the world ray.
257  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal, RealType &wTime) const
258  {
259  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
260  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
261  mTester.getWorldPosAndNml(world, normal);
262  wTime = mTester.getWorldTime();
263  return true;
264  }
265 
266 private:
267 
268  mutable SearchImplT mTester;
269 
270 };// LevelSetRayIntersector
271 
272 
273 ////////////////////////////////////// VolumeRayIntersector //////////////////////////////////////
274 
275 
276 /// @brief This class provides the public API for intersecting a ray
277 /// with a generic (e.g. density) volume.
278 /// @details Internally it performs the actual hierarchical tree node traversal.
279 /// @warning Use the (default) copy-constructor to make sure each
280 /// computational thread has their own instance of this class. This is
281 /// important since it contains a ValueAccessor that is
282 /// not thread-safe and a CoordBBox of the active voxels that should
283 /// not be re-computed for each thread. However copying is very efficient.
284 /// @par Example:
285 /// @code
286 /// // Create an instance for the master thread
287 /// VolumeRayIntersector inter(grid);
288 /// // For each additional thread use the copy constructor. This
289 /// // amortizes the overhead of computing the bbox of the active voxels!
290 /// VolumeRayIntersector inter2(inter);
291 /// // Before each ray-traversal set the index ray.
292 /// iter.setIndexRay(ray);
293 /// // or world ray
294 /// iter.setWorldRay(ray);
295 /// // Now you can begin the ray-marching using consecutive calls to VolumeRayIntersector::march
296 /// double t0=0, t1=0;// note the entry and exit times are with respect to the INDEX ray
297 /// while ( inter.march(t0, t1) ) {
298 /// // perform line-integration between t0 and t1
299 /// }}
300 /// @endcode
301 template<typename GridT,
302  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
303  typename RayT = math::Ray<Real> >
305 {
306 public:
307  using GridType = GridT;
308  using RayType = RayT;
309  using RealType = typename RayT::RealType;
310  using RootType = typename GridT::TreeType::RootNodeType;
312 
313  static_assert(NodeLevel >= 0 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range");
314 
315  /// @brief Grid constructor
316  /// @param grid Generic grid to intersect rays against.
317  /// @param dilationCount The number of voxel dilations performed
318  /// on (a boolean copy of) the input grid. This allows the
319  /// intersector to account for the size of interpolation kernels
320  /// in client code.
321  /// @throw RuntimeError if the voxels of the grid are not uniform
322  /// or the grid is empty.
323  VolumeRayIntersector(const GridT& grid, int dilationCount = 0)
324  : mIsMaster(true)
325  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
326  , mGrid(&grid)
327  , mAccessor(*mTree)
328  {
329  if (!grid.hasUniformVoxels() ) {
330  OPENVDB_THROW(RuntimeError,
331  "VolumeRayIntersector only supports uniform voxels!");
332  }
333  if ( grid.empty() ) {
334  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
335  }
336 
337  // Dilate active voxels to better account for the size of interpolation kernels
338  tools::dilateVoxels(*mTree, dilationCount);
339 
340  mTree->root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
341 
342  mBBox.max().offset(1);//padding so the bbox of a node becomes (origin,origin + node_dim)
343  }
344 
345  /// @brief Grid and BBox constructor
346  /// @param grid Generic grid to intersect rays against.
347  /// @param bbox The axis-aligned bounding-box in the index space of the grid.
348  /// @warning It is assumed that bbox = (min, min + dim) where min denotes
349  /// to the smallest grid coordinates and dim are the integer length of the bbox.
350  /// @throw RuntimeError if the voxels of the grid are not uniform
351  /// or the grid is empty.
353  : mIsMaster(true)
354  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
355  , mGrid(&grid)
356  , mAccessor(*mTree)
357  , mBBox(bbox)
358  {
359  if (!grid.hasUniformVoxels() ) {
360  OPENVDB_THROW(RuntimeError,
361  "VolumeRayIntersector only supports uniform voxels!");
362  }
363  if ( grid.empty() ) {
364  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
365  }
366  }
367 
368  /// @brief Shallow copy constructor
369  /// @warning This copy constructor creates shallow copies of data
370  /// members of the instance passed as the argument. For
371  /// performance reasons we are not using shared pointers (their
372  /// mutex-lock impairs multi-threading).
374  : mIsMaster(false)
375  , mTree(other.mTree)//shallow copy
376  , mGrid(other.mGrid)//shallow copy
377  , mAccessor(*mTree)//initialize new (vs deep copy)
378  , mRay(other.mRay)//deep copy
379  , mTmax(other.mTmax)//deep copy
380  , mBBox(other.mBBox)//deep copy
381  {
382  }
383 
384  /// @brief Destructor
385  ~VolumeRayIntersector() { if (mIsMaster) delete mTree; }
386 
387  /// @brief Return @c false if the index ray misses the bbox of the grid.
388  /// @param iRay Ray represented in index space.
389  /// @warning Call this method (or setWorldRay) before the ray
390  /// traversal starts and use the return value to decide if further
391  /// marching is required.
392  inline bool setIndexRay(const RayT& iRay)
393  {
394  mRay = iRay;
395  const bool hit = mRay.clip(mBBox);
396  if (hit) mTmax = mRay.t1();
397  return hit;
398  }
399 
400  /// @brief Return @c false if the world ray misses the bbox of the grid.
401  /// @param wRay Ray represented in world space.
402  /// @warning Call this method (or setIndexRay) before the ray
403  /// traversal starts and use the return value to decide if further
404  /// marching is required.
405  /// @details Since hit times are computed with respect to the ray
406  /// represented in index space of the current grid, it is
407  /// recommended that either the client code uses getIndexPos to
408  /// compute index position from hit times or alternatively keeps
409  /// an instance of the index ray and instead uses setIndexRay to
410  /// initialize the ray.
411  inline bool setWorldRay(const RayT& wRay)
412  {
413  return this->setIndexRay(wRay.worldToIndex(*mGrid));
414  }
415 
416  inline typename RayT::TimeSpan march()
417  {
418  const typename RayT::TimeSpan t = mHDDA.march(mRay, mAccessor);
419  if (t.t1>0) mRay.setTimes(t.t1 + math::Delta<RealType>::value(), mTmax);
420  return t;
421  }
422 
423  /// @brief Return @c true if the ray intersects active values,
424  /// i.e. either active voxels or tiles. Only when a hit is
425  /// detected are t0 and t1 updated with the corresponding entry
426  /// and exit times along the INDEX ray!
427  /// @note Note that t0 and t1 are only resolved at the node level
428  /// (e.g. a LeafNode with active voxels) as opposed to the individual
429  /// active voxels.
430  /// @param t0 If the return value > 0 this is the time of the
431  /// first hit of an active tile or leaf.
432  /// @param t1 If the return value > t0 this is the time of the
433  /// first hit (> t0) of an inactive tile or exit point of the
434  /// BBOX for the leaf nodes.
435  /// @warning t0 and t1 are computed with respect to the ray represented in
436  /// index space of the current grid, not world space!
437  inline bool march(RealType& t0, RealType& t1)
438  {
439  const typename RayT::TimeSpan t = this->march();
440  t.get(t0, t1);
441  return t.valid();
442  }
443 
444  /// @brief Generates a list of hits along the ray.
445  ///
446  /// @param list List of hits represented as time spans.
447  ///
448  /// @note ListType is a list of RayType::TimeSpan and is required to
449  /// have the two methods: clear() and push_back(). Thus, it could
450  /// be std::vector<typename RayType::TimeSpan> or
451  /// std::deque<typename RayType::TimeSpan>.
452  template <typename ListType>
453  inline void hits(ListType& list)
454  {
455  mHDDA.hits(mRay, mAccessor, list);
456  }
457 
458  /// @brief Return the floating-point index position along the
459  /// current index ray at the specified time.
460  inline Vec3R getIndexPos(RealType time) const { return mRay(time); }
461 
462  /// @brief Return the floating-point world position along the
463  /// current index ray at the specified time.
464  inline Vec3R getWorldPos(RealType time) const { return mGrid->indexToWorld(mRay(time)); }
465 
467  {
468  return time*mGrid->transform().baseMap()->applyJacobian(mRay.dir()).length();
469  }
470 
471  /// @brief Return a const reference to the input grid.
472  const GridT& grid() const { return *mGrid; }
473 
474  /// @brief Return a const reference to the (potentially dilated)
475  /// bool tree used to accelerate the ray marching.
476  const TreeT& tree() const { return *mTree; }
477 
478  /// @brief Return a const reference to the BBOX of the grid
479  const math::CoordBBox& bbox() const { return mBBox; }
480 
481  /// @brief Print bbox, statistics, memory usage and other information.
482  /// @param os a stream to which to write textual information
483  /// @param verboseLevel 1: print bbox only; 2: include boolean tree
484  /// statistics; 3: include memory usage
485  void print(std::ostream& os = std::cout, int verboseLevel = 1)
486  {
487  if (verboseLevel>0) {
488  os << "BBox: " << mBBox << std::endl;
489  if (verboseLevel==2) {
490  mTree->print(os, 1);
491  } else if (verboseLevel>2) {
492  mTree->print(os, 2);
493  }
494  }
495  }
496 
497 private:
498  using AccessorT = typename tree::ValueAccessor<const TreeT,/*IsSafe=*/false>;
499 
500  const bool mIsMaster;
501  TreeT* mTree;
502  const GridT* mGrid;
503  AccessorT mAccessor;
504  RayT mRay;
505  RealType mTmax;
506  math::CoordBBox mBBox;
508 
509 };// VolumeRayIntersector
510 
511 
512 //////////////////////////////////////// LinearSearchImpl ////////////////////////////////////////
513 
514 
515 /// @brief Implements linear iterative search for an iso-value of
516 /// the level set along the direction of the ray.
517 ///
518 /// @note Since this class is used internally in
519 /// LevelSetRayIntersector (define above) and LevelSetHDDA (defined below)
520 /// client code should never interact directly with its API. This also
521 /// explains why we are not concerned with the fact that several of
522 /// its methods are unsafe to call unless roots were already detected.
523 ///
524 /// @details It is approximate due to the limited number of iterations
525 /// which can can be defined with a template parameter. However the default value
526 /// has proven surprisingly accurate and fast. In fact more iterations
527 /// are not guaranteed to give significantly better results.
528 ///
529 /// @warning Since the root-searching algorithm is approximate
530 /// (first-order) it is possible to miss intersections if the
531 /// iso-value is too close to the inside or outside of the narrow
532 /// band (typically a distance less than a voxel unit).
533 ///
534 /// @warning Since this class internally stores a ValueAccessor it is NOT thread-safe,
535 /// so make sure to give each thread its own instance. This of course also means that
536 /// the cost of allocating an instance should (if possible) be amortized over
537 /// as many ray intersections as possible.
538 template<typename GridT, int Iterations, typename RealT>
539 class LinearSearchImpl
540 {
541 public:
544  using ValueT = typename GridT::ValueType;
545  using AccessorT = typename GridT::ConstAccessor;
547 
548  /// @brief Constructor from a grid.
549  /// @throw RunTimeError if the grid is empty.
550  /// @throw ValueError if the isoValue is not inside the narrow-band.
551  LinearSearchImpl(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>())
552  : mStencil(grid),
553  mIsoValue(isoValue),
554  mMinValue(isoValue - ValueT(2 * grid.voxelSize()[0])),
555  mMaxValue(isoValue + ValueT(2 * grid.voxelSize()[0]))
556  {
557  if ( grid.empty() ) {
558  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
559  }
560  if (mIsoValue<= -grid.background() ||
561  mIsoValue>= grid.background() ){
562  OPENVDB_THROW(ValueError, "The iso-value must be inside the narrow-band!");
563  }
564  grid.tree().root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
565  }
566 
567  /// @brief Return the iso-value used for ray-intersections
568  const ValueT& getIsoValue() const { return mIsoValue; }
569 
570  /// @brief Return @c false if the ray misses the bbox of the grid.
571  /// @param iRay Ray represented in index space.
572  /// @warning Call this method before the ray traversal starts.
573  inline bool setIndexRay(const RayT& iRay)
574  {
575  mRay = iRay;
576  return mRay.clip(mBBox);//did it hit the bbox
577  }
578 
579  /// @brief Return @c false if the ray misses the bbox of the grid.
580  /// @param wRay Ray represented in world space.
581  /// @warning Call this method before the ray traversal starts.
582  inline bool setWorldRay(const RayT& wRay)
583  {
584  mRay = wRay.worldToIndex(mStencil.grid());
585  return mRay.clip(mBBox);//did it hit the bbox
586  }
587 
588  /// @brief Get the intersection point in index space.
589  /// @param xyz The position in index space of the intersection.
590  inline void getIndexPos(VecT& xyz) const { xyz = mRay(mTime); }
591 
592  /// @brief Get the intersection point in world space.
593  /// @param xyz The position in world space of the intersection.
594  inline void getWorldPos(VecT& xyz) const { xyz = mStencil.grid().indexToWorld(mRay(mTime)); }
595 
596  /// @brief Get the intersection point and normal in world space
597  /// @param xyz The position in world space of the intersection.
598  /// @param nml The surface normal in world space of the intersection.
599  inline void getWorldPosAndNml(VecT& xyz, VecT& nml)
600  {
601  this->getIndexPos(xyz);
602  mStencil.moveTo(xyz);
603  nml = mStencil.gradient(xyz);
604  nml.normalize();
605  xyz = mStencil.grid().indexToWorld(xyz);
606  }
607 
608  /// @brief Return the time of intersection along the index ray.
609  inline RealT getIndexTime() const { return mTime; }
610 
611  /// @brief Return the time of intersection along the world ray.
612  inline RealT getWorldTime() const
613  {
614  return mTime*mStencil.grid().transform().baseMap()->applyJacobian(mRay.dir()).length();
615  }
616 
617 private:
618 
619  /// @brief Initiate the local voxel intersection test.
620  /// @warning Make sure to call this method before the local voxel intersection test.
621  inline void init(RealT t0)
622  {
623  mT[0] = t0;
624  mV[0] = static_cast<ValueT>(this->interpValue(t0));
625  }
626 
627  inline void setRange(RealT t0, RealT t1) { mRay.setTimes(t0, t1); }
628 
629  /// @brief Return a const reference to the ray.
630  inline const RayT& ray() const { return mRay; }
631 
632  /// @brief Return true if a node of the specified type exists at ijk.
633  template <typename NodeT>
634  inline bool hasNode(const Coord& ijk)
635  {
636  return mStencil.accessor().template probeConstNode<NodeT>(ijk) != nullptr;
637  }
638 
639  /// @brief Return @c true if an intersection is detected.
640  /// @param ijk Grid coordinate of the node origin or voxel being tested.
641  /// @param time Time along the index ray being tested.
642  /// @warning Only if an intersection is detected is it safe to
643  /// call getIndexPos, getWorldPos and getWorldPosAndNml!
644  inline bool operator()(const Coord& ijk, RealT time)
645  {
646  ValueT V;
647  if (mStencil.accessor().probeValue(ijk, V) &&//within narrow band
648  V>mMinValue && V<mMaxValue) {// and close to iso-value?
649  mT[1] = time;
650  mV[1] = static_cast<ValueT>(this->interpValue(time));
651  if (math::ZeroCrossing(mV[0], mV[1])) {
652  mTime = this->interpTime();
654  for (int n=0; Iterations>0 && n<Iterations; ++n) {//resolved at compile-time
655  V = static_cast<ValueT>(this->interpValue(mTime));
656  const int m = math::ZeroCrossing(mV[0], V) ? 1 : 0;
657  mV[m] = V;
658  mT[m] = mTime;
659  mTime = this->interpTime();
660  }
662  return true;
663  }
664  mT[0] = mT[1];
665  mV[0] = mV[1];
666  }
667  return false;
668  }
669 
670  inline RealT interpTime()
671  {
672  assert( math::isApproxLarger(mT[1], mT[0], RealT(1e-6) ) );
673  return mT[0]+(mT[1]-mT[0])*mV[0]/(mV[0]-mV[1]);
674  }
675 
676  inline RealT interpValue(RealT time)
677  {
678  const VecT pos = mRay(time);
679  mStencil.moveTo(pos);
680  return mStencil.interpolation(pos) - mIsoValue;
681  }
682 
683  template<typename, int> friend struct math::LevelSetHDDA;
684 
685  RayT mRay;
686  StencilT mStencil;
687  RealT mTime;//time of intersection
688  ValueT mV[2];
689  RealT mT[2];
690  const ValueT mIsoValue, mMinValue, mMaxValue;
691  math::CoordBBox mBBox;
692 };// LinearSearchImpl
693 
694 } // namespace tools
695 } // namespace OPENVDB_VERSION_NAME
696 } // namespace openvdb
697 
698 #endif // OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
699 
700 // Copyright (c) 2012-2018 DreamWorks Animation LLC
701 // All rights reserved. This software is distributed under the
702 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
bool setIndexRay(const RayT &iRay)
Return false if the ray misses the bbox of the grid.
LinearSearchImpl(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor from a grid.
void getWorldPos(VecT &xyz) const
Get the intersection point in world space.
void setTimes(RealT t0=math::Delta< RealT >::value(), RealT t1=std::numeric_limits< RealT >::max())
Definition: Ray.h:103
bool clip(const Vec3T &center, RealT radius)
Return true if this ray intersects the specified sphere.
Definition: Ray.h:244
bool intersectsWS(const RayType &wRay) const
Return true if the world-space ray intersects the level set.
bool ZeroCrossing(const Type &a, const Type &b)
Return true if the interval [a, b] includes zero, i.e., if either a or b is zero or if they have diff...
Definition: Math.h:707
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:377
bool setIndexRay(const RayT &iRay)
Return false if the index ray misses the bbox of the grid.
This class provides the public API for intersecting a ray with a generic (e.g. density) volume...
RealT getWorldTime() const
Return the time of intersection along the world ray.
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
GT_API const UT_StringHolder time
const math::CoordBBox & bbox() const
Return a const reference to the BBOX of the grid.
Ray worldToIndex(const GridType &grid) const
Return a new ray in the index space of the specified grid, assuming the existing ray is represented i...
Definition: Ray.h:198
bool intersectsIS(const RayType &iRay) const
Return true if the index-space ray intersects the level set.
static bool test(TesterT &tester)
Definition: DDA.h:178
RootNodeType & root()
Return this tree's root node.
Definition: Tree.h:307
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
const GridT & grid() const
Return a const reference to the input grid.
tree::Tree< typename RootType::template ValueConverter< bool >::Type > TreeT
bool setWorldRay(const RayT &wRay)
Return false if the ray misses the bbox of the grid.
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:363
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:204
Vec3R getWorldPos(RealType time) const
Return the floating-point world position along the current index ray at the specified time...
VolumeRayIntersector(const GridT &grid, const math::CoordBBox &bbox)
Grid and BBox constructor.
GLdouble n
Definition: glcorearb.h:2007
bool intersectsWS(const RayType &wRay, Vec3Type &world) const
Return true if the world-space ray intersects the level set.
void print(std::ostream &os=std::cout, int verboseLevel=1) const override
Print statistics, memory usage and other information about this tree.
Definition: Tree.h:2239
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
const TreeT & tree() const
Return a const reference to the (potentially dilated) bool tree used to accelerate the ray marching...
Delta for small floating-point offsets.
Definition: Math.h:124
void getIndexPos(VecT &xyz) const
Get the intersection point in index space.
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:264
bool intersectsIS(const RayType &iRay, Vec3Type &xyz, RealType &iTime) const
Return true if the index-space ray intersects the level set.
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal, RealType &wTime) const
Return true if the world-space ray intersects the level set.
bool intersectsIS(const RayType &iRay, RealType &iTime) const
Return true if the index-space ray intersects the level set.
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:73
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:129
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal) const
Return true if the world-space ray intersects the level set.
void hits(ListType &list)
Generates a list of hits along the ray.
Implements linear iterative search for an iso-value of the level set along the direction of the ray...
void print(std::ostream &os=std::cout, int verboseLevel=1)
Print bbox, statistics, memory usage and other information.
ValueType interpolation(const math::Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:335
Helper class that implements Hierarchical Digital Differential Analyzers and is specialized for ray i...
Definition: DDA.h:172
GLsizei const GLfloat * value
Definition: glcorearb.h:823
bool intersectsWS(const RayType &wRay, Vec3Type &world, RealType &wTime) const
Return true if the world-space ray intersects the level set.
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:200
const Vec3T & dir() const
Definition: Ray.h:126
This class provides the public API for intersecting a ray with a narrow-band level set...
bool march(RealType &t0, RealType &t1)
Return true if the ray intersects active values, i.e. either active voxels or tiles. Only when a hit is detected are t0 and t1 updated with the corresponding entry and exit times along the INDEX ray!
RealT getIndexTime() const
Return the time of intersection along the index ray.
Implementation of morphological dilation and erosion.
typename GridT::TreeType::RootNodeType RootType
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
bool intersectsWS(const RayType &wRay, RealType &wTime) const
Return true if the world-space ray intersects the level set.
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:130
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:518
void getWorldPosAndNml(VecT &xyz, VecT &nml)
Get the intersection point and normal in world space.
Vec3R getIndexPos(RealType time) const
Return the floating-point index position along the current index ray at the specified time...
Digital Differential Analyzers specialized for VDB.
VolumeRayIntersector(const VolumeRayIntersector &other)
Shallow copy constructor.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
VolumeRayIntersector(const GridT &grid, int dilationCount=0)
Grid constructor.
LevelSetRayIntersector(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor.
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
bool setWorldRay(const RayT &wRay)
Return false if the world ray misses the bbox of the grid.
A Ray class.
Helper class that implements Hierarchical Digital Differential Analyzers for ray intersections agains...
Definition: DDA.h:215
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
OPENVDB_STATIC_SPECIALIZATION void dilateVoxels(TreeType &tree, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically dilate all leaf-level active voxels in a tree using one of three nearest neighbor conne...
Definition: Morphology.h:858
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:794
bool isApproxLarger(const Type &a, const Type &b, const Type &tolerance)
Return true if a is larger than b to within the given tolerance, i.e., if b - a < tolerance...
Definition: Math.h:386