[MOAB-dev] r4332 - in MOAB/trunk: src test

[bmsmith6 at wisc.edu]

Author: bmsmith
Date: 2010-12-07 10:44:44 -0600 (Tue, 07 Dec 2010)
New Revision: 4332

Modified:
   MOAB/trunk/src/OrientedBox.cpp
   MOAB/trunk/src/OrientedBox.hpp
   MOAB/trunk/test/OBBTest.cpp
Log:
-Add ability to search in reverse direction to box.intersect_ray. This
is useful for particle tracking through overlaps and numerical precision
problems.
-Check distance limits (in both directions) against sides of the box.
Makes tree pruning due to distance limits more efficient.
-Existing use of the function should not change.
-Added a few unit tests for using a distance limit in the reverse
direction.



Modified: MOAB/trunk/src/OrientedBox.cpp
===================================================================
--- MOAB/trunk/src/OrientedBox.cpp	2010-12-06 21:45:43 UTC (rev 4331)
+++ MOAB/trunk/src/OrientedBox.cpp	2010-12-07 16:44:44 UTC (rev 4332)
@@ -467,38 +467,123 @@
 //}
 
 
+/* This is a helper function to check limits on ray length, turning the box-ray 
+ * intersection test into a box-segment intersection test. Use this to test the
+ * limits against one side (plane) of the box. The side of the box (plane) is
+ * normal to an axis.
+ *
+ *   normal_par_pos  Coordinate of particle's position along axis normal to side of box
+ *   normal_par_dir  Coordinate of particle's direction along axis normal to side of box
+ *   half_extent     Distance between center of box and side of box
+ *   nonneg_ray_len  Maximum ray length in positive direction (in front of origin)
+ *   neg_ray_len     Maximum ray length in negative direction (behind origin)
+ *   return          true if intersection with plane occurs within distance limits
+ *
+ * ray equation:   intersection = origin + dist*direction
+ * plane equation: intersection.plane_normal = half_extent
+ * 
+ * Assume plane_normal and direction are unit vectors. Combine equations.
+ *
+ *     (origin + dist*direction).plane_normal = half_extent
+ *     origin.plane_normal + dist*direction.plane_normal = half_extent
+ *     dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
+ *
+ * Although this solves for distance, avoid floating point division.
+ *
+ *     dist*direction.plane_normal = half_extent - origin.plane_normal
+ *
+ * Use inequalities to test dist against ray length limits. Be aware that 
+ * inequalities change due to sign of direction.plane_normal.
+ */ 
+inline bool check_ray_limits(const double  normal_par_pos,
+                             const double  normal_par_dir,
+                             const double  half_extent,
+                             const double* nonneg_ray_len,
+                             const double* neg_ray_len ) {
+
+  const double extent_pos_diff = half_extent - normal_par_pos;
+
+  // limit in positive direction
+  if(nonneg_ray_len) { // should be 0 <= t <= nonneg_ray_len
+    assert(0 <= *nonneg_ray_len);
+    if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
+      if(*nonneg_ray_len*normal_par_dir>=extent_pos_diff && extent_pos_diff>=0) return true;
+    } else if(normal_par_dir<0) {
+      if(*nonneg_ray_len*normal_par_dir<=extent_pos_diff && extent_pos_diff<=0) return true;