Add funcion intersection_point for Line/Plane

FastEnvelope/examples/FastEnvelope/fast.cpp

@@ -366,16 +366,12 @@ struct Envelope {
     for (int i = 0; i < cutp.size(); i++){
       const Plane& plane_i = prism[cutp[i]];
 
-      CGAL::cpp11::result_of<eIntersect_3(eLine_3, ePlane_3)>::type
+      boost::optional<ePoint_3> op = CGAL::intersection_point(line, plane_i.eplane);
-        result = CGAL::intersection(line, plane_i.eplane);
+      if(! op){
-      if(! result){
         std::cout <<  "there must be an intersection 2" << std::endl;
       }
 
-      const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+      const ePoint_3& ip = *op;
-      CGAL_assertion(ipp != nullptr);
-
-      const ePoint_3& ip = *ipp;
 
       for(int j = 0; j < cutp.size(); j++) {
         if (i == j){
@@ -412,18 +408,14 @@ struct Envelope {
                                                           const std::vector<unsigned int> &prismindex, const int &jump, int &id) const
   {
     CGAL::Oriented_side ori;
-    CGAL::cpp11::result_of<eIntersect_3(eLine_3, ePlane_3)>::type
+    boost::optional<ePoint_3> op = CGAL::intersection_point(eline, plane.eplane);
-      result = CGAL::intersection(eline, plane.eplane);
 #ifdef TRACE
-    if(! result){
+    if(! op){
       std::cout <<  "there must be an intersection 3" << std::endl;
     }
 #endif
 
-    const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+    const ePoint_3& ip = *op;
-    CGAL_assertion(ipp != nullptr);
-
-    const ePoint_3& ip = *ipp;
     for (int i = 0; i < prismindex.size(); i++){
       if (prismindex[i] == jump){
         continue;
@@ -551,10 +543,16 @@ struct Envelope {
                                */
                                ) const
   {
+
+    Point_3 itp(CGAL::approx(tp).x().inf(), CGAL::approx(tp).y().inf(), CGAL::approx(tp).z().inf());
+    Point_3 itq(CGAL::approx(tq).x().inf(), CGAL::approx(tq).y().inf(), CGAL::approx(tq).z().inf());
+    Point_3 itr(CGAL::approx(tr).x().inf(), CGAL::approx(tr).y().inf(), CGAL::approx(tr).z().inf());
+    Point_3 in = itp + CGAL::cross_product((itp - itq), (itp - itr));
+    ePoint_3 n(in.x(),in.y(),in.z());
     // return 2;
     // todo:  what do we test here with n ?
     // todo : do this not with Epeck
-    ePoint_3 n = tp + CGAL::cross_product((tp - tq), (tp - tr));
+    // ePoint_3 n = tp + CGAL::cross_product((tp - tq), (tp - tr));
 
     if (CGAL::orientation(n, tp, tq, tr) == 0){
         std::cout << "todo degeneration handling" << std::endl;
@@ -583,7 +581,13 @@ struct Envelope {
                     const ePoint_3& tr,
                     const ePoint_3& ip) const
   {
-    ePoint_3 n = tp + CGAL::cross_product((tp - tq), (tp - tr));
+    Point_3 itp(CGAL::approx(tp).x().inf(), CGAL::approx(tp).y().inf(), CGAL::approx(tp).z().inf());
+    Point_3 itq(CGAL::approx(tq).x().inf(), CGAL::approx(tq).y().inf(), CGAL::approx(tq).z().inf());
+    Point_3 itr(CGAL::approx(tr).x().inf(), CGAL::approx(tr).y().inf(), CGAL::approx(tr).z().inf());
+    Point_3 in = itp + CGAL::cross_product((itp - itq), (itp - itr));
+    ePoint_3 n(in.x(),in.y(),in.z());
+
+    //ePoint_3 n = tp + CGAL::cross_product((tp - tq), (tp - tr));
 #if 0
     if (Predicates::orient_3d(n, tp, q, tr) == 0)
       {
@@ -725,17 +729,14 @@ struct Envelope {
           // std::cout << *(seg0[k]) << "       " <<   *(seg1[k]) << std::endl;
 
           eLine_3 eline(*(seg0[k]), *(seg1[k]));
-          CGAL::cpp11::result_of<eIntersect_3(eLine_3, ePlane_3)>::type
+          boost::optional<ePoint_3> op = CGAL::intersection_point(eline, plane_i.eplane);
-            result = CGAL::intersection(eline, plane_i.eplane);
+          if(! op){
-          if(! result){
 #ifdef TRACE
             std::cout <<  "there must be an intersection 6" << std::endl;
 #endif
           }
 
-          const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+          const ePoint_3& ip = *op;
-          CGAL_assertion(ipp != nullptr);
-          const ePoint_3& ip = *ipp;
 
           for (int j = 0; j < cutp.size(); j++){
             if (i == j){
@@ -900,18 +901,14 @@ struct Envelope {
   {
     eLine_3 eline(segpoint0,segpoint1); // todo replace parameter of function
 
-    CGAL::cpp11::result_of<eIntersect_3(eLine_3, ePlane_3)>::type
+     boost::optional<ePoint_3> op = CGAL::intersection_point(eline, eplane);
-      result = CGAL::intersection(eline, eplane);
+    if(! op){
-    if(! result){
  #ifdef TRACE
       std::cout <<  "there must be an intersection 9" << std::endl;
  #endif
     }
 
-    const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+    const ePoint_3& ip = *op;
-    CGAL_assertion(ipp != nullptr);
-
-    const ePoint_3& ip = *ipp;
     int tot, fid, ori;
     for (int i = 0; i < prismindex.size(); i++){
       if (prismindex[i] == jump){
@@ -985,18 +982,14 @@ struct Envelope {
   {
     eLine_3 eline(segpoint0,segpoint1); // todo replace parameter of function
 
-    CGAL::cpp11::result_of<eIntersect_3(eLine_3, ePlane_3)>::type
+    boost::optional<ePoint_3> op = CGAL::intersection_point(eline, eplane);
-      result = CGAL::intersection(eline, eplane);
+    if(! op){
-    if(! result){
 #ifdef TRACE
       std::cout <<  "there must be an intersection 10" << std::endl;
 #endif
     }
 
-    const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+    const ePoint_3& ip = *op;
-    CGAL_assertion(ipp != nullptr);
-
-    const ePoint_3& ip = *ipp;
 
     int tot, ori, fid;
 
@@ -1472,20 +1465,15 @@ struct Envelope {
 
             // AF: We moved the intersection here
             // In case there is no intersection point we continue
-            CGAL::cpp11::result_of<eIntersect_3(ePlane_3, ePlane_3, ePlane_3)>::type
+            boost::optional<ePoint_3>
-              result = CGAL::intersection(etriangle_eplane,
+              op = CGAL::intersection_point(etriangle_eplane,
                                             halfspace[jump1][intersect_face[queue[i]][k]].eplane,
                                             halfspace[jump2][intersect_face[queue[j]][h]].eplane);
-            if(! result){
+            if(! op){
               continue;
             }
 
-            const ePoint_3* ipp = boost::get<ePoint_3>(&*result);
+            const ePoint_3& ip = *op;
-            if(ipp == nullptr){
-              continue;
-            }
-
-            const ePoint_3& ip = *ipp;
 
 
             cut = is_3_triangle_cut(etriangle[0], etriangle[1], etriangle[2], ip);
@@ -1945,10 +1933,10 @@ int main(int argc, char* argv[])
             if(! CGAL::collinear(env_vertices[i], env_vertices[j],env_vertices[k])){
               if(envelope(env_vertices[i],  env_vertices[j], env_vertices[k])){
                 inside_count++;
-                inside << i << " " << j << " "<< k <<std::endl;
+                // inside << i << " " << j << " "<< k <<std::endl;
               } else{
                 outside_count++;
-                outside << i << " " << j << " "<< k <<std::endl;
+                // outside << i << " " << j << " "<< k <<std::endl;
               }
             }
           }

FastEnvelope/examples/FastEnvelope/fastE.cpp

@@ -88,10 +88,10 @@ int main(int argc, char* argv[]) {
           }
         }
       }
+
+  }
   std::cout << inside_count << " " << outside_count << std::endl;
   std::cout << t.time() << " sec." << std::endl;
-  }
-
 
 
   return 0;

Filtered_kernel/include/CGAL/Lazy.h

@@ -882,11 +882,46 @@ struct Lazy_construction_optional
 
     } catch (Uncertain_conversion_exception&) {
       Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);
-      boost::optional<typename EK::Point_3> oep = ec(CGAL::exact(l1),CGAL::exact(l2),CGAL::exact(l3));
+      boost::optional<typename EK::Point_3> oep = ec(CGAL::exact(l1), CGAL::exact(l2), CGAL::exact(l3));
       if(oep == boost::none){
         return boost::none;
       }
-      typedef Lazy_rep_0<typename AK::Point_3, typename EK::Point_3,E2A> LazyPointRep;
+      typedef Lazy_rep_0<typename AK::Point_3, typename EK::Point_3, E2A> LazyPointRep;
+      const typename EK::Point_3 ep = *oep;
+      LazyPointRep *rep = new LazyPointRep(ep);
+      typename LK::Point_3 lp(rep);
+      return boost::make_optional(lp);
+    }
+    // AF can we get here??
+    return boost::none;
+  }
+
+  // for Intersect_point_3 with  Plane_3  Line_3
+  template <typename L1, typename L2>
+  result_type operator()(const L1& l1, const L2& l2) const
+  {
+    Protect_FPU_rounding<Protection> P;
+
+    try {
+      boost::optional<typename AK::Point_3> oap = ac(CGAL::approx(l1),CGAL::approx(l2));
+      if(oap == boost::none){
+        return boost::none;
+      }
+      // Now we have to construct a rep for a lazy point with the three lazy planes.
+      typedef Lazy_rep_optional_n<typename AK::Point_3, typename EK::Point_3, AC, EC, E2A, L1, L2> LazyPointRep;
+      const typename AK::Point_3 ap = *oap;
+      LazyPointRep *rep = new LazyPointRep(ap, ec, l1, l2);
+      typename LK::Point_3 lp(rep);
+      return boost::make_optional(lp);
+
+
+    } catch (Uncertain_conversion_exception&) {
+      Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);
+      boost::optional<typename EK::Point_3> oep = ec(CGAL::exact(l1), CGAL::exact(l2));
+      if(oep == boost::none){
+        return boost::none;
+      }
+      typedef Lazy_rep_0<typename AK::Point_3, typename EK::Point_3, E2A> LazyPointRep;
       const typename EK::Point_3 ep = *oep;
       LazyPointRep *rep = new LazyPointRep(ep);
       typename LK::Point_3 lp(rep);

Intersections_3/include/CGAL/Intersections_3/Line_3_Plane_3.h

@@ -22,6 +22,22 @@
 namespace CGAL {
 CGAL_INTERSECTION_FUNCTION(Plane_3, Line_3, 3)
 CGAL_DO_INTERSECT_FUNCTION(Plane_3, Line_3, 3)
+
+template < class K >
+inline
+boost::optional<typename K::Point_3>
+intersection_point(const Plane_3<K>& plane, const Line_3<K>& line)
+{
+  return K().intersect_point_3_object()(plane, line);
+}
+
+  template < class K >
+inline
+boost::optional<typename K::Point_3>
+  intersection_point(const Line_3<K>& line, const Plane_3<K>& plane)
+{
+  return K().intersect_point_3_object()(plane, line);
+}
 }
 
 #endif // CGAL_INTERSECTIONS_3_LINE_PLANE_3_H

Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h

@@ -79,6 +79,33 @@ do_intersect(const Plane_3<R> &plane1, const Plane_3<R> &plane2,
 namespace Intersections {
 namespace internal {
 
+template <class K>
+boost::optional<typename K::Point_3>
+intersection_point(const typename K::Plane_3 &plane,
+                   const typename K::Line_3 &line,
+                   const K& /*k*/)
+{
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Direction_3 Direction_3;
+  typedef typename K::RT RT;
+
+  const Point_3 &line_pt = line.point();
+  const Direction_3 &line_dir = line.direction();
+
+  RT num = plane.a()*line_pt.hx() + plane.b()*line_pt.hy()
+    + plane.c()*line_pt.hz() + wmult_hw((K*)0, plane.d(), line_pt);
+  RT den = plane.a()*line_dir.dx() + plane.b()*line_dir.dy()
+    + plane.c()*line_dir.dz();
+  if (den == 0) {
+    return boost::none;
+  }
+  return boost::make_optional(Point_3(den*line_pt.hx()-num*line_dir.dx(),
+                                      den*line_pt.hy()-num*line_dir.dy(),
+                                      den*line_pt.hz()-num*line_dir.dz(),
+                                      wmult_hw((K*)0, den, line_pt)));
+}
+
+
 template <class K>
 typename Intersection_traits<K, typename K::Plane_3, typename K::Line_3>::result_type
 intersection(const typename K::Plane_3  &plane,

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -3573,6 +3573,7 @@ namespace CommonKernelFunctors {
   {
   public:
     typedef typename K::Point_3     Point_3;
+    typedef typename K::Line_3      Line_3;
     typedef typename K::Plane_3     Plane_3;
     typedef typename boost::optional<Point_3> result_type;
 
@@ -3581,6 +3582,12 @@ namespace CommonKernelFunctors {
     {
       return Intersections::internal::intersection_point(pl1, pl2, pl3, K() );
     }
+
+    result_type
+    operator()(const Plane_3& plane, const Line_3& line) const
+    {
+      return Intersections::internal::intersection_point(plane, line, K() );
+    }
   };