OpenVDB  2.0.0
RayIntersector.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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 //
57 
58 
59 #ifndef OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
60 #define OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
61 
62 #include <openvdb/math/Ray.h>
63 #include <openvdb/math/Stencils.h>
64 #include <openvdb/Grid.h>
65 #include <openvdb/Types.h>
66 #include <openvdb/math/Math.h>
67 #include <boost/utility.hpp>
68 #include <boost/type_traits/is_floating_point.hpp>
69 
70 namespace openvdb {
72 namespace OPENVDB_VERSION_NAME {
73 namespace tools {
74 
75 // Helper class that implements hierarchical Digital Differential Analyzers
76 // specialized for ray intersections with level sets
77 template <typename GridT, int NodeLevel> struct LevelSetHDDA;
78 
79 
82 template <typename GridT, int NodeLevel> struct VolumeHDDA;
83 
84 // Helper class that implements the actual search of the zero-crossing
85 // of the level set along the direction of a ray. This particular
86 // implementation uses iterative linear search.
87 template<typename GridT, int Iterations = 0, typename RealT = double>
89 
90 
92 
93 
111 template<typename GridT,
112  typename SearchImplT = LinearSearchImpl<GridT>,
113  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
114  typename RayT = math::Ray<Real> >
116 {
117 public:
118  typedef GridT GridType;
119  typedef RayT RayType;
120  typedef typename RayT::Vec3T Vec3Type;
121  typedef typename GridT::ValueType ValueT;
122  typedef typename GridT::TreeType TreeT;
123 
124  BOOST_STATIC_ASSERT( NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1);
125  BOOST_STATIC_ASSERT(boost::is_floating_point<ValueT>::value);
126 
129  LevelSetRayIntersector(const GridT& grid): mTester(grid)
130  {
131  if (!grid.hasUniformVoxels() ) {
133  "LevelSetRayIntersector only supports uniform voxels!");
134  }
135  if (grid.getGridClass() != GRID_LEVEL_SET) {
137  "LevelSetRayIntersector only supports level sets!"
138  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
139  }
140  }
141 
144  bool intersectsIS(const RayType& iRay) const
145  {
146  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
147  return LevelSetHDDA<TreeT, NodeLevel>::test(mTester);
148  }
149 
154  bool intersectsIS(const RayType& iRay, Vec3Type& xyz) const
155  {
156  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
157  if (!LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
158  mTester.getIndexPos(xyz);
159  return true;
160  }
161 
164  bool intersectsWS(const RayType& wRay) const
165  {
166  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
167  return LevelSetHDDA<TreeT, NodeLevel>::test(mTester);
168  }
169 
174  bool intersectsWS(const RayType& wRay, Vec3Type& world) const
175  {
176  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
177  if (!LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
178  mTester.getWorldPos(world);
179  return true;
180  }
181 
188  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal) const
189  {
190  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
191  if (!LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
192  mTester.getWorldPosAndNml(world, normal);
193  return true;
194  }
195 
196 private:
197 
198  mutable SearchImplT mTester;
199 
200 };// LevelSetRayIntersector
201 
202 
204 
205 
232 template<typename GridT,
233  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
234  typename RayT = math::Ray<Real> >
236 {
237 public:
238  typedef GridT GridType;
239  typedef RayT RayType;
240  typedef typename GridT::TreeType TreeT;
241 
242  BOOST_STATIC_ASSERT( NodeLevel >= 0 && NodeLevel < int(TreeT::DEPTH)-1);
243 
247  VolumeRayIntersector(const GridT& grid): mGrid(&grid), mAccessor(grid.tree())
248  {
249  if (!grid.hasUniformVoxels() ) {
251  "VolumeRayIntersector only supports uniform voxels!");
252  }
253  if ( grid.empty() ) {
254  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
255  }
256  grid.tree().root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
257 
258  mBBox.max().offset(1);//padding so the bbox of a node becomes (origin,origin + node_dim)
259  }
260 
266  inline bool setIndexRay(const RayT& iRay)
267  {
268  mRay = iRay;
269  const bool hit = mRay.clip(mBBox);
270  if (hit) mTmax = mRay.t1();
271  return hit;
272  }
273 
285  inline bool setWorldRay(const RayT& wRay)
286  {
287  return this->setIndexRay(wRay.worldToIndex(*mGrid));
288  }
289 
299  inline int march(Real& t0, Real& t1)
300  {
301  const int n = mRay.test() ? VolumeHDDA<TreeT, NodeLevel>::test(*this) : 0;
302  if (n>0) {
303  t0 = mRay.t0();
304  t1 = mRay.t1();
305  }
306  mRay.setTimes(mRay.t1() + 1e-9, mTmax);
307  return n;
308  }
309 
312  inline Vec3R getIndexPos(Real time) const { return mRay(time); }
313 
316  inline Vec3R getWorldPos(Real time) const { return mGrid->indexToWorld(mRay(time)); }
317 
318 private:
319 
320  inline void setRange(Real t0, Real t1) { mRay.setTimes(t0, t1); }
321 
323  inline const RayT& ray() const { return mRay; }
324 
326  template <typename NodeT>
327  inline bool hasNode(const Coord& ijk)
328  {
329  return mAccessor.template probeConstNode<NodeT>(ijk) != NULL;
330  }
331 
332  bool isValueOn(const Coord& ijk) { return mAccessor.isValueOn(ijk); }
333 
334  template<typename, int> friend struct VolumeHDDA;
335 
336  typedef typename GridT::ConstAccessor AccessorT;
337 
338  const GridT* mGrid;
339  AccessorT mAccessor;
340  RayT mRay;
341  Real mTmax;
342  math::CoordBBox mBBox;
343 
344 };// VolumeRayIntersector
345 
346 
348 
349 
369 template<typename GridT, int Iterations, typename RealT>
370 class LinearSearchImpl
371 {
372 public:
374  typedef typename GridT::ValueType ValueT;
375  typedef typename GridT::ConstAccessor AccessorT;
377  typedef typename StencilT::Vec3Type Vec3T;
378 
380  LinearSearchImpl(const GridT& grid)
381  : mStencil(grid), mThreshold(2*grid.voxelSize()[0])
382  {
383  if ( grid.empty() ) {
384  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
385  }
386  grid.tree().root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
387  }
388 
392  inline bool setIndexRay(const RayT& iRay)
393  {
394  mRay = iRay;
395  return mRay.clip(mBBox);//did it hit the bbox
396  }
397 
401  inline bool setWorldRay(const RayT& wRay)
402  {
403  mRay = wRay.worldToIndex(mStencil.grid());
404  return mRay.clip(mBBox);//did it hit the bbox
405  }
406 
409  inline void getIndexPos(Vec3d& xyz) const { xyz = mRay(mTime); }
410 
413  inline void getWorldPos(Vec3d& xyz) const { xyz = mStencil.grid().indexToWorld(mRay(mTime)); }
414 
418  inline void getWorldPosAndNml(Vec3d& xyz, Vec3d& nml)
419  {
420  this->getIndexPos(xyz);
421  mStencil.moveTo(xyz);
422  nml = mStencil.gradient(xyz);
423  nml.normalize();
424  xyz = mStencil.grid().indexToWorld(xyz);
425  }
426 
427 private:
428 
431  inline void init(RealT t0)
432  {
433  mT[0] = t0;
434  mV[0] = this->interpValue(t0);
435  }
436 
437  inline void setRange(RealT t0, RealT t1) { mRay.setTimes(t0, t1); }
438 
440  inline const RayT& ray() const { return mRay; }
441 
443  template <typename NodeT>
444  inline bool hasNode(const Coord& ijk)
445  {
446  return mStencil.accessor().template probeConstNode<NodeT>(ijk) != NULL;
447  }
448 
454  inline bool operator()(const Coord& ijk, RealT time)
455  {
456  ValueT V;
457  if (mStencil.accessor().probeValue(ijk, V) &&//inside narrow band?
458  math::Abs(V)<mThreshold) {// close to zero-crossing?
459  mT[1] = time;
460  mV[1] = this->interpValue(time);
461  if (math::ZeroCrossing(mV[0], mV[1])) {
462  mTime = this->interpTime();
464  for (int n=0; Iterations>0 && n<Iterations; ++n) {//resolved at compile-time
465  V = this->interpValue(mTime);
466  const int m = math::ZeroCrossing(mV[0], V) ? 1 : 0;
467  mV[m] = V;
468  mT[m] = mTime;
469  mTime = this->interpTime();
470  }
472  return true;
473  }
474  mT[0] = mT[1];
475  mV[0] = mV[1];
476  }
477  return false;
478  }
479 
480  inline RealT interpTime()
481  {
482  assert(math::isApproxLarger(mT[1], mT[0], 1e-6));
483  return mT[0]+(mT[1]-mT[0])*mV[0]/(mV[0]-mV[1]);
484  }
485 
486  inline RealT interpValue(RealT time)
487  {
488  const Vec3R pos = mRay(time);
489  mStencil.moveTo(pos);
490  return mStencil.interpolation(pos);
491  }
492 
493  template<typename, int> friend struct LevelSetHDDA;
494 
495  RayT mRay;
496  StencilT mStencil;
497  RealT mTime;
498  ValueT mV[2];
499  RealT mT[2];
500  ValueT mThreshold;
501  math::CoordBBox mBBox;
502 };// LinearSearchImpl
503 
505 
506 
509 template<typename TreeT, int NodeLevel>
510 struct LevelSetHDDA
511 {
512  typedef typename TreeT::RootNodeType::NodeChainType ChainT;
513  typedef typename boost::mpl::at<ChainT, boost::mpl::int_<NodeLevel> >::type NodeT;
514 
515  template <typename TesterT>
516  static bool test(TesterT& tester)
517  {
519  do {
520  if (tester.template hasNode<NodeT>(dda.voxel())) {
521  tester.setRange(dda.time(), dda.next());
522  if (LevelSetHDDA<TreeT, NodeLevel-1>::test(tester)) return true;
523  }
524  } while(dda.step());
525  return false;
526  }
527 };
528 
531 template<typename TreeT>
532 struct LevelSetHDDA<TreeT, -1>
533 {
534  template <typename TesterT>
535  static bool test(TesterT& tester)
536  {
537  math::DDA<typename TesterT::RayT, 0> dda(tester.ray());
538  tester.init(dda.time());
539  do { if (tester(dda.voxel(), dda.next())) return true; } while(dda.step());
540  return false;
541  }
542 };
543 
544 
546 
547 
550 template <typename TreeT, int NodeLevel>
551 struct VolumeHDDA
552 {
553  typedef typename TreeT::RootNodeType::NodeChainType ChainT;
554  typedef typename boost::mpl::at<ChainT, boost::mpl::int_<NodeLevel> >::type NodeT;
555 
556  template <typename TesterT>
557  static int test(TesterT& tester)
558  {
560  do {
561  if (tester.template hasNode<NodeT>(dda.voxel())) {//child node
562  tester.setRange(dda.time(), dda.next());
563  int hit = VolumeHDDA<TreeT, NodeLevel-1>::test(tester);
564  if (hit > 0) return hit;
565  } else if (tester.isValueOn(dda.voxel())) {//active tile
566  tester.setRange(dda.time(), dda.next());
567  return 1;//active tile
568  }
569  } while (dda.step());
570  return 0;//no hits
571  }
572 };
573 
576 template <typename TreeT>
577 struct VolumeHDDA<TreeT, -1>
578 {
579  template <typename TesterT>
580  static int test(TesterT&) { return 2; }//hit leaf so don't traverse voxels.
581 };
582 
583 } // namespace tools
584 } // namespace OPENVDB_VERSION_NAME
585 } // namespace openvdb
586 
587 #endif // OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
588 
589 // Copyright (c) 2012-2013 DreamWorks Animation LLC
590 // All rights reserved. This software is distributed under the
591 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
Vec3R getIndexPos(Real time) const
Return the floating-point index position along the current index ray at the specified time...
Definition: RayIntersector.h:312
Implements linear iterative search for a zero-crossing of the level set along along the direction of ...
Definition: RayIntersector.h:88
This class provides the public API for intersecting a ray with a generic (e.g. density) volume...
Definition: RayIntersector.h:235
LevelSetRayIntersector(const GridT &grid)
Constructor.
Definition: RayIntersector.h:129
Definition: RayIntersector.h:82
void init(const RayT &ray, RealT startTime, RealT maxTime)
Definition: Ray.h:333
GridT::TreeType TreeT
Definition: RayIntersector.h:240
Helper class that implements Hierarchical Digital Differential Analyzers and is specialized for ray i...
Definition: RayIntersector.h:77
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
bool setWorldRay(const RayT &wRay)
Return false the ray misses the bbox of the grid.
Definition: RayIntersector.h:401
Definition: Ray.h:52
double Real
Definition: Types.h:62
math::Ray< RealT > RayT
Definition: RayIntersector.h:373
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:154
Definition: Exceptions.h:86
bool intersectsIS(const RayType &iRay) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:144
GridT GridType
Definition: RayIntersector.h:118
void getWorldPos(Vec3d &xyz) const
Get the intersection point in world space.
Definition: RayIntersector.h:413
TreeT::RootNodeType::NodeChainType ChainT
Definition: RayIntersector.h:512
GridT GridType
Definition: RayIntersector.h:238
bool setIndexRay(const RayT &iRay)
Return false the ray misses the bbox of the grid.
Definition: RayIntersector.h:392
math::BoxStencil< GridT > StencilT
Definition: RayIntersector.h:376
GridT::ValueType ValueT
Definition: RayIntersector.h:121
RayT RayType
Definition: RayIntersector.h:239
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:188
A Digital Differential Analyzer specialized for OpenVDB grids.
Definition: Ray.h:319
#define OPENVDB_VERSION_NAME
Definition: version.h:45
math::Vec3< Real > Vec3R
Definition: Types.h:74
GridT::TreeType TreeT
Definition: RayIntersector.h:122
void getIndexPos(Vec3d &xyz) const
Get the intersection point in index space.
Definition: RayIntersector.h:409
bool intersectsWS(const RayType &wRay, Vec3Type &world) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:174
bool clip(const Vec3T &center, RealT radius)
Return true if this ray intersects the specified sphere.
Definition: Ray.h:210
static int test(TesterT &tester)
Definition: RayIntersector.h:557
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:96
Vec3< double > Vec3d
Definition: Vec3.h:605
VolumeRayIntersector(const GridT &grid)
Constructor.
Definition: RayIntersector.h:247
RayT RayType
Definition: RayIntersector.h:119
bool setIndexRay(const RayT &iRay)
Return false if the index ray misses the bbox of the grid.
Definition: RayIntersector.h:266
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:229
void getWorldPosAndNml(Vec3d &xyz, Vec3d &nml)
Get the intersection point and normal in world space.
Definition: RayIntersector.h:418
static bool test(TesterT &tester)
Definition: RayIntersector.h:535
TreeT::RootNodeType::NodeChainType ChainT
Definition: RayIntersector.h:553
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:97
RayT::Vec3T Vec3Type
Definition: RayIntersector.h:120
static bool test(TesterT &tester)
Definition: RayIntersector.h:516
Vec3R getWorldPos(Real time) const
Return the floating-point world position along the current index ray at the specified time...
Definition: RayIntersector.h:316
StencilT::Vec3Type Vec3T
Definition: RayIntersector.h:377
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:246
static int test(TesterT &)
Definition: RayIntersector.h:580
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 &lt; tolerance...
Definition: Math.h:344
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:164
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:640
int march(Real &t0, Real &t1)
Return 0 if not hit was detected. A return value of 1 means it hit an active tile, and a return value of 2 means it hit a LeafNode. Only when a hit is detected are t0 and t1 updated with the corresponding entry and exit times along the INDEX ray!
Definition: RayIntersector.h:299
This class provides the public API for intersecting a ray with a narrow-band level set...
Definition: RayIntersector.h:115
GridT::ValueType ValueT
Definition: RayIntersector.h:374
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
GridT::ConstAccessor AccessorT
Definition: RayIntersector.h:375
bool intersectsWS(const RayType &wRay) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:164
LinearSearchImpl(const GridT &grid)
Constructor from a grid.
Definition: RayIntersector.h:380
bool setWorldRay(const RayT &wRay)
Return false if the world ray misses the bbox of the grid.
Definition: RayIntersector.h:285
Definition: Types.h:137
boost::mpl::at< ChainT, boost::mpl::int_< NodeLevel > >::type NodeT
Definition: RayIntersector.h:554
boost::mpl::at< ChainT, boost::mpl::int_< NodeLevel > >::type NodeT
Definition: RayIntersector.h:513