HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RayIntersector.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2017 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 <hboost/utility.hpp>
70 #include <hboost/type_traits/is_floating_point.hpp>
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  typedef GridT GridType;
113  typedef RayT RayType;
114  typedef typename RayT::RealType RealType;
115  typedef typename RayT::Vec3T Vec3Type;
116  typedef typename GridT::ValueType ValueT;
117  typedef typename GridT::TreeType TreeT;
118 
119  HBOOST_STATIC_ASSERT( NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1);
121 
122  /// @brief Constructor
123  /// @param grid level set grid to intersect rays against.
124  /// @param isoValue optional iso-value for the ray-intersection.
125  LevelSetRayIntersector(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>())
126  : mTester(grid, isoValue)
127  {
128  if (!grid.hasUniformVoxels() ) {
129  OPENVDB_THROW(RuntimeError,
130  "LevelSetRayIntersector only supports uniform voxels!");
131  }
132  if (grid.getGridClass() != GRID_LEVEL_SET) {
133  OPENVDB_THROW(RuntimeError,
134  "LevelSetRayIntersector only supports level sets!"
135  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
136  }
137  }
138 
139  /// @brief Return the iso-value used for ray-intersections
140  const ValueT& getIsoValue() const { return mTester.getIsoValue(); }
141 
142  /// @brief Return @c true if the index-space ray intersects the level set.
143  /// @param iRay ray represented in index space.
144  bool intersectsIS(const RayType& iRay) const
145  {
146  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
148  }
149 
150  /// @brief Return @c true if the index-space ray intersects the level set
151  /// @param iRay ray represented in index space.
152  /// @param iTime if an intersection was found it is assigned the time of the
153  /// intersection along the index ray.
154  bool intersectsIS(const RayType& iRay, RealType &iTime) const
155  {
156  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
157  iTime = mTester.getIndexTime();
159  }
160 
161  /// @brief Return @c true if the index-space ray intersects the level set.
162  /// @param iRay ray represented in index space.
163  /// @param xyz if an intersection was found it is assigned the
164  /// intersection point in index space, otherwise it is unchanged.
165  bool intersectsIS(const RayType& iRay, Vec3Type& xyz) const
166  {
167  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
168  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
169  mTester.getIndexPos(xyz);
170  return true;
171  }
172 
173  /// @brief Return @c true if the index-space ray intersects the level set.
174  /// @param iRay ray represented in index space.
175  /// @param xyz if an intersection was found it is assigned the
176  /// intersection point in index space, otherwise it is unchanged.
177  /// @param iTime if an intersection was found it is assigned the time of the
178  /// intersection along the index ray.
179  bool intersectsIS(const RayType& iRay, Vec3Type& xyz, RealType &iTime) const
180  {
181  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
182  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
183  mTester.getIndexPos(xyz);
184  iTime = mTester.getIndexTime();
185  return true;
186  }
187 
188  /// @brief Return @c true if the world-space ray intersects the level set.
189  /// @param wRay ray represented in world space.
190  bool intersectsWS(const RayType& wRay) const
191  {
192  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
194  }
195 
196  /// @brief Return @c true if the world-space ray intersects the level set.
197  /// @param wRay ray represented in world space.
198  /// @param wTime if an intersection was found it is assigned the time of the
199  /// intersection along the world ray.
200  bool intersectsWS(const RayType& wRay, RealType &wTime) const
201  {
202  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
203  wTime = mTester.getWorldTime();
205  }
206 
207  /// @brief Return @c true if the world-space ray intersects the level set.
208  /// @param wRay ray represented in world space.
209  /// @param world if an intersection was found it is assigned the
210  /// intersection point in world space, otherwise it is unchanged
211  bool intersectsWS(const RayType& wRay, Vec3Type& world) const
212  {
213  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
214  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
215  mTester.getWorldPos(world);
216  return true;
217  }
218 
219  /// @brief Return @c true if the world-space ray intersects the level set.
220  /// @param wRay ray represented in world space.
221  /// @param world if an intersection was found it is assigned the
222  /// intersection point in world space, otherwise it is unchanged.
223  /// @param wTime if an intersection was found it is assigned the time of the
224  /// intersection along the world ray.
225  bool intersectsWS(const RayType& wRay, Vec3Type& world, RealType &wTime) const
226  {
227  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
228  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
229  mTester.getWorldPos(world);
230  wTime = mTester.getWorldTime();
231  return true;
232  }
233 
234  /// @brief Return @c true if the world-space ray intersects the level set.
235  /// @param wRay ray represented in world space.
236  /// @param world if an intersection was found it is assigned the
237  /// intersection point in world space, otherwise it is unchanged.
238  /// @param normal if an intersection was found it is assigned the normal
239  /// of the level set surface in world space, otherwise it is unchanged.
240  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal) const
241  {
242  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
243  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
244  mTester.getWorldPosAndNml(world, normal);
245  return true;
246  }
247 
248  /// @brief Return @c true if the world-space ray intersects the level set.
249  /// @param wRay ray represented in world space.
250  /// @param world if an intersection was found it is assigned the
251  /// intersection point in world space, otherwise it is unchanged.
252  /// @param normal if an intersection was found it is assigned the normal
253  /// of the level set surface in world space, otherwise it is unchanged.
254  /// @param wTime if an intersection was found it is assigned the time of the
255  /// intersection along the world ray.
256  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal, RealType &wTime) const
257  {
258  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
259  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
260  mTester.getWorldPosAndNml(world, normal);
261  wTime = mTester.getWorldTime();
262  return true;
263  }
264 
265 private:
266 
267  mutable SearchImplT mTester;
268 
269 };// LevelSetRayIntersector
270 
271 
272 ////////////////////////////////////// VolumeRayIntersector //////////////////////////////////////
273 
274 
275 /// @brief This class provides the public API for intersecting a ray
276 /// with a generic (e.g. density) volume.
277 /// @details Internally it performs the actual hierarchical tree node traversal.
278 /// @warning Use the (default) copy-constructor to make sure each
279 /// computational thread has their own instance of this class. This is
280 /// important since it contains a ValueAccessor that is
281 /// not thread-safe and a CoordBBox of the active voxels that should
282 /// not be re-computed for each thread. However copying is very efficient.
283 /// @par Example:
284 /// @code
285 /// // Create an instance for the master thread
286 /// VolumeRayIntersector inter(grid);
287 /// // For each additional thread use the copy constructor. This
288 /// // amortizes the overhead of computing the bbox of the active voxels!
289 /// VolumeRayIntersector inter2(inter);
290 /// // Before each ray-traversal set the index ray.
291 /// iter.setIndexRay(ray);
292 /// // or world ray
293 /// iter.setWorldRay(ray);
294 /// // Now you can begin the ray-marching using consecutive calls to VolumeRayIntersector::march
295 /// double t0=0, t1=0;// note the entry and exit times are with respect to the INDEX ray
296 /// while ( inter.march(t0, t1) ) {
297 /// // perform line-integration between t0 and t1
298 /// }}
299 /// @endcode
300 template<typename GridT,
301  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
302  typename RayT = math::Ray<Real> >
304 {
305 public:
306  typedef GridT GridType;
307  typedef RayT RayType;
308  typedef typename RayT::RealType RealType;
309  typedef typename GridT::TreeType::RootNodeType RootType;
311 
312  HBOOST_STATIC_ASSERT( NodeLevel >= 0 && NodeLevel < int(TreeT::DEPTH)-1);
313 
314  /// @brief Grid constructor
315  /// @param grid Generic grid to intersect rays against.
316  /// @param dilationCount The number of voxel dilations performed
317  /// on (a boolean copy of) the input grid. This allows the
318  /// intersector to account for the size of interpolation kernels
319  /// in client code.
320  /// @throw RuntimeError if the voxels of the grid are not uniform
321  /// or the grid is empty.
322  VolumeRayIntersector(const GridT& grid, int dilationCount = 0)
323  : mIsMaster(true)
324  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
325  , mGrid(&grid)
326  , mAccessor(*mTree)
327  {
328  if (!grid.hasUniformVoxels() ) {
329  OPENVDB_THROW(RuntimeError,
330  "VolumeRayIntersector only supports uniform voxels!");
331  }
332  if ( grid.empty() ) {
333  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
334  }
335 
336  // Dilate active voxels to better account for the size of interpolation kernels
337  tools::dilateVoxels(*mTree, dilationCount);
338 
339  mTree->root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
340 
341  mBBox.max().offset(1);//padding so the bbox of a node becomes (origin,origin + node_dim)
342  }
343 
344  /// @brief Grid and BBox constructor
345  /// @param grid Generic grid to intersect rays against.
346  /// @param bbox The axis-aligned bounding-box in the index space of the grid.
347  /// @warning It is assumed that bbox = (min, min + dim) where min denotes
348  /// to the smallest grid coordinates and dim are the integer length of the bbox.
349  /// @throw RuntimeError if the voxels of the grid are not uniform
350  /// or the grid is empty.
352  : mIsMaster(true)
353  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
354  , mGrid(&grid)
355  , mAccessor(*mTree)
356  , mBBox(bbox)
357  {
358  if (!grid.hasUniformVoxels() ) {
359  OPENVDB_THROW(RuntimeError,
360  "VolumeRayIntersector only supports uniform voxels!");
361  }
362  if ( grid.empty() ) {
363  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
364  }
365  }
366 
367  /// @brief Shallow copy constructor
368  /// @warning This copy constructor creates shallow copies of data
369  /// members of the instance passed as the argument. For
370  /// performance reasons we are not using shared pointers (their
371  /// mutex-lock impairs multi-threading).
373  : mIsMaster(false)
374  , mTree(other.mTree)//shallow copy
375  , mGrid(other.mGrid)//shallow copy
376  , mAccessor(*mTree)//initialize new (vs deep copy)
377  , mRay(other.mRay)//deep copy
378  , mTmax(other.mTmax)//deep copy
379  , mBBox(other.mBBox)//deep copy
380  {
381  }
382 
383  /// @brief Destructor
384  ~VolumeRayIntersector() { if (mIsMaster) delete mTree; }
385 
386  /// @brief Return @c false if the index ray misses the bbox of the grid.
387  /// @param iRay Ray represented in index space.
388  /// @warning Call this method (or setWorldRay) before the ray
389  /// traversal starts and use the return value to decide if further
390  /// marching is required.
391  inline bool setIndexRay(const RayT& iRay)
392  {
393  mRay = iRay;
394  const bool hit = mRay.clip(mBBox);
395  if (hit) mTmax = mRay.t1();
396  return hit;
397  }
398 
399  /// @brief Return @c false if the world ray misses the bbox of the grid.
400  /// @param wRay Ray represented in world space.
401  /// @warning Call this method (or setIndexRay) before the ray
402  /// traversal starts and use the return value to decide if further
403  /// marching is required.
404  /// @details Since hit times are computed with respect to the ray
405  /// represented in index space of the current grid, it is
406  /// recommended that either the client code uses getIndexPos to
407  /// compute index position from hit times or alternatively keeps
408  /// an instance of the index ray and instead uses setIndexRay to
409  /// initialize the ray.
410  inline bool setWorldRay(const RayT& wRay)
411  {
412  return this->setIndexRay(wRay.worldToIndex(*mGrid));
413  }
414 
415  inline typename RayT::TimeSpan march()
416  {
417  const typename RayT::TimeSpan t = mHDDA.march(mRay, mAccessor);
418  if (t.t1>0) mRay.setTimes(t.t1 + math::Delta<RealType>::value(), mTmax);
419  return t;
420  }
421 
422  /// @brief Return @c true if the ray intersects active values,
423  /// i.e. either active voxels or tiles. Only when a hit is
424  /// detected are t0 and t1 updated with the corresponding entry
425  /// and exit times along the INDEX ray!
426  /// @note Note that t0 and t1 are only resolved at the node level
427  /// (e.g. a LeafNode with active voxels) as opposed to the individual
428  /// active voxels.
429  /// @param t0 If the return value > 0 this is the time of the
430  /// first hit of an active tile or leaf.
431  /// @param t1 If the return value > t0 this is the time of the
432  /// first hit (> t0) of an inactive tile or exit point of the
433  /// BBOX for the leaf nodes.
434  /// @warning t0 and t1 are computed with respect to the ray represented in
435  /// index space of the current grid, not world space!
436  inline bool march(RealType& t0, RealType& t1)
437  {
438  const typename RayT::TimeSpan t = this->march();
439  t.get(t0, t1);
440  return t.valid();
441  }
442 
443  /// @brief Generates a list of hits along the ray.
444  ///
445  /// @param list List of hits represented as time spans.
446  ///
447  /// @note ListType is a list of RayType::TimeSpan and is required to
448  /// have the two methods: clear() and push_back(). Thus, it could
449  /// be std::vector<typename RayType::TimeSpan> or
450  /// std::deque<typename RayType::TimeSpan>.
451  template <typename ListType>
452  inline void hits(ListType& list)
453  {
454  mHDDA.hits(mRay, mAccessor, list);
455  }
456 
457  /// @brief Return the floating-point index position along the
458  /// current index ray at the specified time.
459  inline Vec3R getIndexPos(RealType time) const { return mRay(time); }
460 
461  /// @brief Return the floating-point world position along the
462  /// current index ray at the specified time.
463  inline Vec3R getWorldPos(RealType time) const { return mGrid->indexToWorld(mRay(time)); }
464 
465  inline RealType getWorldTime(RealType time) const
466  {
467  return time*mGrid->transform().baseMap()->applyJacobian(mRay.dir()).length();
468  }
469 
470  /// @brief Return a const reference to the input grid.
471  const GridT& grid() const { return *mGrid; }
472 
473  /// @brief Return a const reference to the (potentially dilated)
474  /// bool tree used to accelerate the ray marching.
475  const TreeT& tree() const { return *mTree; }
476 
477  /// @brief Return a const reference to the BBOX of the grid
478  const math::CoordBBox& bbox() const { return mBBox; }
479 
480  /// @brief Print bbox, statistics, memory usage and other information.
481  /// @param os a stream to which to write textual information
482  /// @param verboseLevel 1: print bbox only; 2: include boolean tree
483  /// statistics; 3: include memory usage
484  void print(std::ostream& os = std::cout, int verboseLevel = 1)
485  {
486  if (verboseLevel>0) {
487  os << "BBox: " << mBBox << std::endl;
488  if (verboseLevel==2) {
489  mTree->print(os, 1);
490  } else if (verboseLevel>2) {
491  mTree->print(os, 2);
492  }
493  }
494  }
495 
496 private:
497 
498  typedef typename tree::ValueAccessor<const TreeT,/*IsSafe=*/false> AccessorT;
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  typedef typename GridT::ValueType ValueT;
545  typedef typename GridT::ConstAccessor AccessorT;
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) != NULL;
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-2017 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:100
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:719
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:376
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.
const math::CoordBBox & bbox() const
Return a const reference to the BBOX of the grid.
HBOOST_STATIC_ASSERT(NodeLevel >=0 &&NodeLevel< int(TreeT::DEPTH)-1)
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.
HBOOST_STATIC_ASSERT(NodeLevel >=-1 &&NodeLevel< int(TreeT::DEPTH)-1)
static bool test(TesterT &tester)
Definition: DDA.h:176
RootNodeType & root()
Return this tree's root node.
Definition: Tree.h:309
const GridT & grid() const
Return a const reference to the input grid.
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:2238
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
#define OPENVDB_VERSION_NAME
Definition: version.h:43
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:132
void getIndexPos(VecT &xyz) const
Get the intersection point in index space.
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:266
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
tree::Tree< typename RootType::template ValueConverter< bool >::Type > TreeT
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:170
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:121
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.
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:503
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
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.
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:213
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:101
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:855
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:398