Getting started with the Bbox_3_Ray_3 do_intersect

Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h

@@ -289,6 +289,232 @@ public:
     return Base::operator()(s,b);
   }
 
+
+
+  result_type 
+  operator()(const Ray_3 &r, const Bbox_3& b) const
+  {
+    CGAL_BRANCH_PROFILER_3("semi-static failures/attempts/calls to   : Do_intersect_3", tmp);
+
+    Get_approx<Point_3> get_approx; // Identity functor for all points
+    // but lazy points.
+    const Point_3& p = r.source(); 
+    const Point_3& q = r.point(1); 
+
+    double px, py, pz, qx, qy, qz;
+    double bxmin = b.xmin(), bymin = b.ymin(), bzmin = b.zmin(), 
+      bxmax = b.xmax(), bymax = b.ymax(), bzmax = b.zmax();
+
+    if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) &&
+        fit_in_double(get_approx(p).z(), pz) &&
+        fit_in_double(get_approx(q).x(), qx) && fit_in_double(get_approx(q).y(), qy) &&
+        fit_in_double(get_approx(q).z(), qz) )   
+    {
+      // std::cerr << "\n"
+      //           << px << " " <<  py << " " <<  pz << "\n" 
+      //           << qx << " " <<  qy << " " <<  qz << "\n" 
+      //           << bxmin << " " <<  bymin << " " <<  bzmin << "\n" 
+      //           << bxmax << " " <<  bymax << " " <<  bzmax << "\n";
+      CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
+
+     
+      // AF: I copy pasted the code to call the max
+
+      // -----------------------------------
+      // treat x coord
+      // -----------------------------------
+      double dmin, dmax, tmin, tmax;
+      if ( qx >= px )  // this is input and needs no epsilon
+      {
+        // different from segment: if(px > bxmax) return false; // segment on the right of bbox
+        //                         if(qx < bxmin) return false; // segment on the left of bbox
+        if (bxmax < px)  return false; // different from segment: added
+        tmax = bxmax - px;
+        
+        dmax = qx - px;
+        if ( bxmin < px ) // tmin < 0 means px is in the x-range of bbox
+        {
+          tmin = 0;
+          dmin = 1;
+        } else {
+          tmin = bxmin - px;
+          dmin = dmax;
+        }
+      }
+      else
+      {
+        // different from segment: if(qx > bxmax) return false; // segment on the right of bbox
+        //                         if(px < bxmin) return false; // segment on the left of bbox
+        if(px < bxmin) return false; // different from segment: added
+        tmax = px - bxmin;
+
+        dmax = px - qx;
+        if ( px < bxmax ) // tmin < 0 means px is in the x-range of bbox
+        {
+          tmin = 0;
+          dmin = 1;
+        } else {
+          tmin = px - bxmax;
+          dmin = dmax;
+        }
+      }   
+
+      // std::cerr << "t1 ";
+   
+      double m = CGAL::abs(tmin), m2;
+      m2 = CGAL::abs(tmax); if(m2 > m) { m = m2; }
+      m2 = CGAL::abs(dmin); if(m2 > m) { m = m2; }
+
+      if(m < 7e-294) {
+        // underflow in the computation of 'error'
+        return Base::operator()(r,b);
+      }
+
+      const double EPS_1 = 3.55618e-15;
+
+      double error =  EPS_1 * m;
+
+      switch(sign_with_error( tmax - dmax, error)) {
+      case POSITIVE:
+        // std::cerr << "t2 ";
+        tmax = 1;
+        dmax = 1;
+        break;
+      case NEGATIVE:
+        break;
+      default:
+        // ambiguity of the comparison tmax > dmin
+        // let's call the exact predicate
+        // std::cerr << "\ntest2 NEED EXACT\n";
+        return Base::operator()(r,b);
+      }
+
+      // -----------------------------------
+      // treat y coord
+      // -----------------------------------
+      double d_, tmin_, tmax_;
+      if ( qy >= py )   // this is input and needs no epsilon
+      {
+        tmin_ = bymin - py;
+        tmax_ = bymax - py;
+        d_ = qy - py;
+      }
+      else
+      {
+        tmin_ = py - bymax;
+        tmax_ = py - bymin;
+        d_ = py - qy;
+      }
+    
+      m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
+      m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
+      m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
+
+      if(m < 3e-147) {
+        // underflow in the computation of 'error'
+        return Base::operator()(r,b);
+      }
+      
+      error =  EPS_1 * m * m;
+
+      // std::cerr << dmin << " " << tmax_ << " " << d_ << " "
+      //           << tmin << " " << dmax << " " << tmin_ << std::endl;
+
+      if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ 
+        // potential overflow on the computation of 'sign1' and 'sign2'
+        return Base::operator()(r,b);
+      }
+      Sign sign1 = sign_with_error( (d_*tmin) - (dmin*tmax_) , error);
+      Sign sign2 = sign_with_error( (dmax*tmin_) - (d_*tmax) , error);
+
+      if(sign1 == POSITIVE || sign2 == POSITIVE) 
+        return false; // We are *sure* the segment is outside the box, on one
+                      // side or the other.
+      if(sign1 == ZERO || sign2 == ZERO) {
+        // std::cerr << "\ntest3 NEED EXACT\n";
+        return Base::operator()(r,b); // We are *unsure*: one *may be*
+                                      // positive.
+      }
+
+      // std::cerr << "t3 ";
+
+      // Here we are sure the two signs are negative. We can continue with
+      // the rest of the function...
+
+      // epsilons needed
+      switch(sign_with_error((dmin*tmin_) - (d_*tmin) , error)) {
+      case POSITIVE:
+        tmin = tmin_;
+        dmin = d_;
+        // std::cerr << "t4 ";
+        break;
+      case NEGATIVE:
+        break;
+      default: // uncertainty
+        // std::cerr << "\ntest4 NEED EXACT\n";
+        return Base::operator()(r,b);
+      }
+    
+      // epsilons needed
+      switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) {
+      case POSITIVE:
+        tmax = tmax_;
+        dmax = d_;
+        // std::cerr << "t5 ";break;
+      case NEGATIVE:
+        break;
+      default: // uncertainty
+        // std::cerr << "\ntest5 NEED EXACT\n";
+        return Base::operator()(r,b);
+      }
+    
+      // -----------------------------------
+      // treat z coord
+      // -----------------------------------
+      if ( qz >= pz )   // this is input and needs no epsilon
+      {
+        tmin_ = bzmin - pz;
+        tmax_ = bzmax - pz;
+        d_ = qz - pz;
+      }
+      else
+      {
+        tmin_ = pz - bzmax;
+        tmax_ = pz - bzmin;
+        d_ = pz - qz;
+      }
+    
+      m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
+      m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
+      m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
+
+      // m may have changed
+      error =  EPS_1 * m * m;
+    
+      if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */ 
+        // potential overflow on the computation of 'sign1' and 'sign2'
+        return Base::operator()(r,b);
+      }
+      sign1 = sign_with_error( (dmin*tmax_) - (d_*tmin) , error);
+      sign2 = sign_with_error( (d_*tmax) - (dmax*tmin_) , error);
+      if(sign1 == NEGATIVE || sign2 == NEGATIVE) {
+        // std::cerr << "f6";
+        return false; // We are *sure* the segment is outside the box, on one
+                      // side or the other.
+      }
+      if(sign1 == ZERO || sign2 == ZERO) {
+        // std::cerr << "\test6 NEED EXACT\n";
+        return Base::operator()(r,b); // We are *unsure*: one *may be*
+                                      // negative.
+      }
+      // std::cerr << "t6";
+      return true; // We are *sure* the two signs are positive.
+    }
+    return Base::operator()(r,b);
+  }
+
+
+
   // Computes the epsilon for Bbox_3_Segment_3_do_intersect.
   static double compute_epsilon_bbox_segment_3()
   {