[MOAB-dev] bug in MBOrientedBox::intersect_ray

[Brandon Smith]

Hi All,

While tracking some lost particles I found a bug in MBOrientedBox::intersect_ray. It was due to subtracting a squared distance. More specifically, the fact that a squared distance becomes larger if the distance is greater than one, but smaller if the distance is less than one (assuming distance is always positive). With Jason's help I tried to make the tolerances used in the test more consistent. Below is a diff of the fix. The fix reduces the number of lost particles by ~50%.

Are there any changes/suggestions to this diff?

Brandon

Index: MBOrientedBox.hpp
===================================================================
--- MBOrientedBox.hpp   (revision 3533)
+++ MBOrientedBox.hpp   (working copy)
@@ -53,8 +53,8 @@
 
   inline double inner_radius() const; //!< radius of inscribed sphere
   inline double outer_radius() const; //!< radius of circumscribed sphere
-  inline double outer_radius_squared() const; //!< square of radius of circumsphere
-  inline double inner_radius_squared() const; //!< square of radius if inscribed sphere
+  inline double outer_radius_squared(const double reps) const; //!< square of (radius+at least epsilon) of circumsphere
+  inline double inner_radius_squared(const double reps) const; //!< square of (radius-epsilon) of inscribed sphere
   inline double volume() const;               //!< volume of box
   inline MBCartVect dimensions() const;       //!< number of dimensions for which box is not flat
   inline double area() const;                 //!< largest side area
@@ -176,24 +176,30 @@
 #endif
 }
 
-double MBOrientedBox::outer_radius_squared() const
+// Add at least epsilon to the radius, before squaring it.
+double MBOrientedBox::outer_radius_squared(const double reps) const
 {
 #if MB_ORIENTED_BOX_OUTER_RADIUS
-  return radius * radius;
+  return (radius+reps)*(radius+reps);
 #elif MB_ORIENTED_BOX_UNIT_VECTORS
-  return length % length;
+  MBCartVect tmp(length[0]+reps,length[1]+reps,length[2]+reps);
+  return tmp % tmp;
 #else
-  const MBCartVect half_diag = axis[0] + axis[1] + axis[2];
+  MBCartVect half_diag = axis[0] + axis[1] + axis[2];
+  half_diag += MBCartVect(reps,reps,reps);
   return half_diag % half_diag;
 #endif
 }
 
-double MBOrientedBox::inner_radius_squared() const
+// Subtract epsilon from the length of the shortest axis, before squaring it.
+double MBOrientedBox::inner_radius_squared(const double reps) const
 {
 #if MB_ORIENTED_BOX_UNIT_VECTORS
-  return length[0] * length[0];
+  return (length[0]-reps) * (length[0]-reps);
 #else
-  return axis[0] % axis[0];
+  MBCartVect tmp = axis[0];
+  tmp -= MBCartVect(reps,reps,reps);
+  return (tmp % tmp);
 #endif
 }
 

Index: MBOrientedBox.cpp
===================================================================
--- MBOrientedBox.cpp   (revision 3533)
+++ MBOrientedBox.cpp   (working copy)
@@ -477,16 +477,27 @@
   const MBCartVect cx = center - b;
   double dist_s = cx % m;
   double dist_sq = cx % cx - (dist_s*dist_s);
-  double max_diagsq = outer_radius_squared();
+  double max_diagsq = outer_radius_squared(reps);
   
     // if greater than the longest diagonal, we don't hit
-  if (dist_sq > max_diagsq+reps)
+  if (dist_sq > max_diagsq)
     return false;
-  if (len && dist_s - max_diagsq > *len)
-    return false;
-  
+
+  // If the closest possible hit is farther than len, we don't want the hit.
+  // Problem: the previous method was wrong because max_diagsq will be greater
+  // then max_diag is max_diag>1 but less than max_diag if max_diag<1.
+  if (len) {
+    if (max_diagsq>1) {
+      if (dist_s - max_diagsq > *len)
+        return false;
+    } else {
+      if (dist_s - 1          > *len)
+        return false;
+    }
+  }
+
     // if smaller than shortest diagonal, we do hit
-  if (dist_sq < inner_radius_squared() - reps && dist_s >= 0.0)
+  if (dist_sq < inner_radius_squared(reps) && dist_s >= 0.0)
     return true;
     
     // get transpose of axes


      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/moab-dev/attachments/20100211/bde7ad44/attachment.htm>