Merge pull request #5978 from janetournois/Mesh_2-add_predicate_for_lloyd-jtournois

Mesh 2 and Kernel - add predicate oriented_side_2(segment, triangle)

Cartesian_kernel/include/CGAL/Cartesian/function_objects.h

@@ -4380,6 +4380,8 @@ namespace CartesianKernelFunctors {
     typedef typename K::Circle_2       Circle_2;
     typedef typename K::Line_2         Line_2;
     typedef typename K::Triangle_2     Triangle_2;
+    typedef typename K::Segment_2      Segment_2;
+    typedef typename K::FT             FT;
   public:
     typedef typename K::Oriented_side  result_type;
 
@@ -4416,6 +4418,30 @@ namespace CartesianKernelFunctors {
         ? result_type(ON_ORIENTED_BOUNDARY)
         : opposite(ot);
     }
+
+    result_type
+    operator()(const Segment_2& s, const Triangle_2& t) const
+    {
+      typename K::Construct_source_2 source;
+      typename K::Construct_target_2 target;
+      const Point_2& a = source(s);
+      const Point_2& b = target(s);
+      CGAL_assertion(a != b);
+
+      typename K::Construct_vertex_2 vertex;
+      const Point_2& p0 = vertex(t, 0);
+      const Point_2& p1 = vertex(t, 1);
+      const Point_2& p2 = vertex(t, 2);
+      CGAL_assertion(p0 != p1 && p1 != p2 && p2 != p0);
+
+      return circumcenter_oriented_side_of_oriented_segmentC2(
+                  a.x(), a.y(),
+                  b.x(), b.y(),
+                  p0.x(), p0.y(),
+                  p1.x(), p1.y(),
+                  p2.x(), p2.y()
+      );
+    }
   };
 
   template <typename K>

Cartesian_kernel/include/CGAL/predicates/kernel_ftC2.h

@@ -708,6 +708,28 @@ power_side_of_oriented_power_circleC2(const FT &px, const FT &py, const FT &pwt,
   return cmpy * sign_of_determinant(dpy, dpz, dqy, dqz);
 }
 
+template <class FT>
+Oriented_side
+circumcenter_oriented_side_of_oriented_segmentC2(const FT& ax,  const FT& ay,
+                                                 const FT& bx,  const FT& by,
+                                                 const FT& p0x, const FT& p0y,
+                                                 const FT& p1x, const FT& p1y,
+                                                 const FT& p2x, const FT& p2y)
+{
+  const FT dX = bx - ax;
+  const FT dY = by - ay;
+  const FT R0 = p0x * p0x + p0y * p0y;
+  const FT R1 = p1x * p1x + p1y * p1y;
+  const FT R2 = p2x * p2x + p2y * p2y;
+  const FT denominator = (p1x - p0x) * (p2y - p0y) +
+                         (p0x - p2x) * (p1y - p0y);
+  const FT det = 2 * denominator * (ax * dY - ay * dX)
+                   - (R2 - R1) * (p0x * dX + p0y * dY)
+                   - (R0 - R2) * (p1x * dX + p1y * dY)
+                   - (R1 - R0) * (p2x * dX + p2y * dY);
+  return CGAL::sign(det);
+}
+
 } //namespace CGAL
 
 #endif  // CGAL_PREDICATES_KERNEL_FTC2_H

Homogeneous_kernel/include/CGAL/Homogeneous/function_objects.h

@@ -4683,6 +4683,7 @@ namespace HomogeneousKernelFunctors {
     typedef typename K::Circle_2       Circle_2;
     typedef typename K::Line_2         Line_2;
     typedef typename K::Triangle_2     Triangle_2;
+    typedef typename K::Segment_2      Segment_2;
   public:
     typedef typename K::Oriented_side  result_type;
 
@@ -4722,6 +4723,41 @@ namespace HomogeneousKernelFunctors {
         ? ON_ORIENTED_BOUNDARY
         : -ot;
     }
+
+    result_type
+    operator()(const Segment_2& s, const Triangle_2& t) const
+    {
+      typename K::Construct_source_2 source;
+      typename K::Construct_target_2 target;
+      typename K::Construct_vertex_2 vertex;
+      typename K::Construct_circumcenter_2 circumcenter;
+      typename K::Orientation_2 orientation;
+
+      const Point_2& a = source(s);
+      const Point_2& b = target(s);
+      CGAL_assertion(a != b);
+
+      const Point_2& p0 = vertex(t, 0);
+      const Point_2& p1 = vertex(t, 1);
+      const Point_2& p2 = vertex(t, 2);
+      CGAL_assertion(p0 != p1 && p1 != p2 && p2 != p0);
+
+      const Point_2 cc = circumcenter(t);
+
+      CGAL::Orientation o_abc = orientation(a, b, cc);
+      if (o_abc == CGAL::COLLINEAR)
+        return CGAL::ON_ORIENTED_BOUNDARY;
+
+      CGAL::Orientation o_abt = orientation(a, b, p0);
+      if (o_abt == CGAL::COLLINEAR)
+        o_abt = orientation(a, b, p1);
+      if (o_abt == CGAL::COLLINEAR)
+        o_abt = orientation(a, b, p2);
+      CGAL_assertion(o_abt != CGAL::COLLINEAR);
+
+      if (o_abc == o_abt) return CGAL::ON_POSITIVE_SIDE;
+      else                return CGAL::ON_NEGATIVE_SIDE;
+    }
   };
 
 

Kernel_23/doc/Kernel_23/CGAL/Projection_traits_3.h

@@ -19,6 +19,7 @@ provides exact predicates.
 \cgalModels `DelaunayTriangulationTraits_2`
 \cgalModels `ConstrainedTriangulationTraits_2`
 \cgalModels `PolygonTraits_2`
+\cgalModels `ConformingDelaunayTriangulationTraits_2`
 
 \sa `CGAL::Projection_traits_xy_3`
 \sa `CGAL::Projection_traits_xz_3`

Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h

@@ -9258,6 +9258,16 @@ public:
   Oriented_side operator()(const Kernel::Triangle_2&t,
                            const Kernel::Point_2&p);
 
+  /*!
+  * returns \ref CGAL::ON_ORIENTED_BOUNDARY,
+  * \ref CGAL::ON_NEGATIVE_SIDE, or the constant \ref CGAL::ON_POSITIVE_SIDE,
+  * depending on the position of the circumcenter of `t` relative
+  * to the oriented supporting line of `s`. The orientation of the
+  * supporting line is the same as the orientation of `s`.
+  */
+  Oriented_side operator()(const Kernel::Segment_2& s,
+                           const Kernel::Triangle_2& t);
+
   /// @}
 
 }; /* end Kernel::OrientedSide_2 */

Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h

@@ -115,6 +115,41 @@ public:
     }
 };
 
+template <class R, int dim>
+class Oriented_side_projected_3
+{
+public:
+  typedef typename R::Segment_2   Segment_2;
+  typedef typename R::Triangle_2  Triangle_2;
+
+  typedef typename R::Point_3     Point;
+  typedef typename R::Segment_3   Segment_3;
+  typedef typename R::Triangle_3  Triangle_3;
+
+  typename R::FT x(const Point& p) const { return Projector<R, dim>::x(p); }
+  typename R::FT y(const Point& p) const { return Projector<R, dim>::y(p); }
+
+  typename R::Point_2 project(const Point& p) const
+  {
+    return typename R::Point_2(x(p), y(p));
+  }
+  Triangle_2 project(const Triangle_3& t) const
+  {
+    typename R::Construct_vertex_3 v;
+    return Triangle_2(project(v(t, 0)), project(v(t, 1)), project(v(t, 2)));
+  }
+  Segment_2 project(const Segment_3& s) const
+  {
+    typename R::Construct_source_3 source;
+    typename R::Construct_target_3 target;
+    return Segment_2(project(source(s)), project(target(s)));
+  }
+  CGAL::Oriented_side operator()(const Segment_3& s, const Triangle_3& t) const
+  {
+    return typename R::Oriented_side_2()(project(s), project(t));
+  }
+};
+
 template <class R,int dim>
 class Side_of_oriented_circle_projected_3
 {
@@ -854,6 +889,7 @@ public:
   typedef typename Projector<R,dim>::Compare_x_2              Compare_x_2;
   typedef typename Projector<R,dim>::Compare_y_2              Compare_y_2;
   typedef Orientation_projected_3<Rp,dim>                     Orientation_2;
+  typedef Oriented_side_projected_3<Rp,dim>                   Oriented_side_2;
   typedef Angle_projected_3<Rp,dim>                           Angle_2;
   typedef Side_of_oriented_circle_projected_3<Rp,dim>         Side_of_oriented_circle_2;
   typedef Less_signed_distance_to_line_projected_3<Rp,dim>    Less_signed_distance_to_line_2;
@@ -1014,6 +1050,10 @@ public:
   orientation_2_object() const
     { return Orientation_2();}
 
+  Oriented_side_2
+  oriented_side_2_object() const
+    { return Oriented_side_2();}
+
   Side_of_oriented_circle_2
   side_of_oriented_circle_2_object() const
     {return Side_of_oriented_circle_2();}

Mesh_2/doc/Mesh_2/Concepts/ConformingDelaunayTriangulationTraits_2.h

@@ -92,6 +92,17 @@ points `p`, `q`, `r` (`q` being the vertex of the angle).
 */
 typedef unspecified_type Angle_2;
 
+/*!
+Predicate object. Must provide the operator
+`CGAL::Oriented_side operator()(Segment_2 s, Triangle_2 t)` that
+returns \ref ON_ORIENTED_BOUNDARY, \ref ON_NEGATIVE_SIDE,
+or \ref ON_POSITIVE_SIDE,
+depending on the position of the circumcenter of `t` relative
+to the oriented supporting line of `s`. The orientation of the
+supporting line is the same as the orientation of `s`.
+*/
+typedef unspecified_type Oriented_side_2;
+
 /// @}
 
 

Mesh_2/include/CGAL/Constrained_voronoi_diagram_2.h

@@ -242,25 +242,9 @@ private:
   bool segment_hides_circumcenter(const Segment& seg,
                                   const Triangle& tr)
   {
-    Point a = seg.source();
-    Point b = seg.target();
-    double dX = b.x() - a.x();
-    double dY = b.y() - a.y();
-
-    const Point& p0 = tr[0];
-    const Point& p1 = tr[1];
-    const Point& p2 = tr[2];
-    double R0 = p0.x()*p0.x() + p0.y()*p0.y();
-    double R1 = p1.x()*p1.x() + p1.y()*p1.y();
-    double R2 = p2.x()*p2.x() + p2.y()*p2.y();
-    double denominator = (p1.x()-p0.x())*(p2.y()-p0.y()) +
-                                   (p0.x()-p2.x())*(p1.y()-p0.y());
-
-    double det = 2*denominator * (a.x()*dY - a.y()*dX)
-                    - (R2-R1) * (p0.x()*dX + p0.y()*dY)
-                    - (R0-R2) * (p1.x()*dX + p1.y()*dY)
-                    - (R1-R0) * (p2.x()*dX + p2.y()*dY);
-    return (det <= 0);
+    typename Geom_traits::Oriented_side_2 os
+      = m_cdt.geom_traits().oriented_side_2_object();
+    return (os(seg, tr) != CGAL::ON_POSITIVE_SIDE);
   }
 
   // tags with their sights, with respect to the Edge constraint,

Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h

@@ -352,7 +352,7 @@ private:
       sum += CGAL::sqrt(*it);
 
 #ifdef CGAL_MESH_2_OPTIMIZER_VERBOSE
-    sum_moves_ = sum/big_moves_.size();
+    sum_moves_ = sum/FT(big_moves_.size());
 #endif
 
     return ( sum/FT(big_moves_.size()) < convergence_ratio_ );

Mesh_2/include/CGAL/Mesh_2/Mesh_sizing_field.h

@@ -169,7 +169,7 @@ private:
     typename Tr::Edge_circulator end = ec;
 
     FT sum_len(0.);
-    FT nb = 0.;
+    unsigned int nb = 0;
     do
     {
       Edge e = *ec;
@@ -187,7 +187,7 @@ private:
     while(++ec != end);
     // nb == 0 could happen if there is an isolated point.
     if( 0 != nb )
-      return sum_len/nb;
+      return sum_len/FT(nb);
     else
      // Use outside faces to compute size of point
       return 1.;//todo