Document where sqrt is done approximately

Algebraic_foundations/include/CGAL/number_utils.h

@@ -104,6 +104,27 @@ sqrt( const AS& x ) {
 }
 
 
+template <typename NT>
+NT approximate_sqrt(const NT& nt, CGAL::Field_tag)
+{
+  return NT(sqrt(CGAL::to_double(nt)));
+}
+
+template <typename NT>
+NT approximate_sqrt(const NT& nt, CGAL::Field_with_sqrt_tag)
+{
+  return sqrt(nt);
+}
+
+template <typename NT>
+NT approximate_sqrt(const NT& nt)
+{
+  typedef CGAL::Algebraic_structure_traits<NT> AST;
+  typedef typename AST::Algebraic_category Algebraic_category; 
+  return approximate_sqrt(nt, Algebraic_category());
+}
+
+
 
 // AST-Functor adapting functions BINARY
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h

@@ -86,13 +86,17 @@ void sum_normals(const PM& pmesh,
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-*    \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
 * \cgalNamedParamsEnd
 *
 * @return the computed normal. The return type is a 3D vector type. It is
 * either deduced from the `geom_traits` \ref namedparameters if provided,
 * or from the geometric traits class deduced from the point property map
 * of `pmesh`.
+*
+* \warning This function involves a square root computation.
+* If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+* will be done approximately.
 */
 template <typename PolygonMesh, typename NamedParameters>
 #ifdef DOXYGEN_RUNNING
@@ -138,9 +142,13 @@ compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-*    \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
 * \cgalNamedParamsEnd
 *
+* \warning This function involves a square root computation.
+* If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+* will be done approximately.
+*/
 */
 template <typename PolygonMesh
           , typename FaceNormalMap
@@ -171,13 +179,18 @@ compute_face_normals(const PolygonMesh& pmesh
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-*    \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
 * \cgalNamedParamsEnd
 *
 * @return the computed normal. The return type is a 3D vector type. It is
 * either deduced from the `geom_traits` \ref namedparameters if provided,
 * or the geometric traits class deduced from the point property map
 * of `pmesh`.
+*
+* \warning This function involves a square root computation.
+* If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+* will be done approximately.
+*/
 */
 template<typename PolygonMesh, typename NamedParameters>
 #ifdef DOXYGEN_RUNNING
@@ -248,9 +261,14 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-*    \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
 * \cgalNamedParamsEnd
 *
+* \warning This function involves a square root computation.
+* If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+* will be done approximately.
+*/
+*
 */
 template <typename PolygonMesh
           , typename VertexNormalMap
@@ -290,9 +308,13 @@ compute_vertex_normals(const PolygonMesh& pmesh
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-*    \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
 * \cgalNamedParamsEnd
 *
+* \warning This function involves a square root computation.
+* If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+* will be done approximately.
+*/
 */
 template <typename PolygonMesh
           , typename VertexNormalMap

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h

@@ -53,6 +53,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
+  *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
   * \cgalNamedParamsEnd
   *
   * @return the length of `h`. The return type `FT` is a number type. It is
@@ -60,6 +61,10 @@ namespace Polygon_mesh_processing {
   * or the geometric traits class deduced from the point property map
   * of `pmesh`.
   *
+  * \warning This function involves a square root computation.
+  * If `FT` does not have a `sqrt()` operation, the square root computation
+  * will be done approximately.
+  *
   * @sa `face_border_length()`
   */
   template<typename PolygonMesh,
@@ -81,7 +86,7 @@ namespace Polygon_mesh_processing {
                             pmesh,
                             CGAL::vertex_point);
 
-    return CGAL::sqrt(CGAL::squared_distance(get(vpm, source(h, pmesh)),
+    return CGAL::approximate_sqrt(CGAL::squared_distance(get(vpm, source(h, pmesh)),
                                                          get(vpm, target(h, pmesh))));
   }
 
@@ -129,6 +134,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
+*    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
   * \cgalNamedParamsEnd
   *
   * @return the length of the sequence of border edges of `face(h, pmesh)`.
@@ -137,6 +143,10 @@ namespace Polygon_mesh_processing {
   * or the geometric traits class deduced from the point property map
   * of `pmesh`.
   *
+  * \warning This function involves a square root computation.
+  * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+  * will be done approximately.
+  *
   * @sa `edge_length()`
   */
   template<typename PolygonMesh,
@@ -151,7 +161,7 @@ namespace Polygon_mesh_processing {
               , const PolygonMesh& pmesh
               , const NamedParameters& np)
   {
-    double result = 0.;
+    typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT result = 0.;
     BOOST_FOREACH(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor haf,
                   halfedges_around_face(h, pmesh))
     {
@@ -187,7 +197,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-  *  \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+  *  \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
   * \cgalNamedParamsEnd
   *
   *@pre `f != boost::graph_traits<TriangleMesh>::%null_face()`
@@ -227,9 +237,10 @@ namespace Polygon_mesh_processing {
 
     typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type traits;
 
-    return traits.compute_area_3_object()(get(vpm, source(hd, tmesh)),
+    return approximate_sqrt(
+             traits.compute_squared_area_3_object()(get(vpm, source(hd, tmesh)),
                                                     get(vpm, target(hd, tmesh)),
-                                          get(vpm, target(nhd, tmesh)));
+                                                    get(vpm, target(nhd, tmesh))));
   }
 
   template<typename TriangleMesh>
@@ -260,7 +271,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-  *  \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+  *  \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel` \cgalParamEnd
   * \cgalNamedParamsEnd
   *
   * @return sum of face areas of `faces`.
@@ -269,6 +280,10 @@ namespace Polygon_mesh_processing {
   * or the geometric traits class deduced from the point property map
   * of `tmesh`.
   *
+  * \warning This function involves a square root computation.
+  * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+  * will be done approximately.
+  *
   * @sa `face_area()`
   */
   template<typename FaceRange,
@@ -284,7 +299,7 @@ namespace Polygon_mesh_processing {
      , const CGAL_PMP_NP_CLASS& np)
   {
     typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
-    double result = 0.;
+    typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT result = 0.;
     BOOST_FOREACH(face_descriptor f, face_range)
     {
       result += face_area(f, tmesh, np);
@@ -292,11 +307,16 @@ namespace Polygon_mesh_processing {
     return result;
   }
 
-  template<typename FaceRange, typename TriangleMesh>
+  template<typename FaceRange, typename TriangleMesh > 
-  double area(FaceRange face_range, const TriangleMesh& tmesh)
+  /* AF   what is wrong here???
+   typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
+                                                            CGAL::vertex_point_t>::type>::Kernel::FT
+  */
+     double 
+area(FaceRange face_range, const TriangleMesh& tmesh)
   {
-    return area(face_range, tmesh,
+    return to_double(area(face_range, tmesh,
-      CGAL::Polygon_mesh_processing::parameters::all_default());
+                          CGAL::Polygon_mesh_processing::parameters::all_default()));
   }
 
   /**
@@ -312,7 +332,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-  *  \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+  *  \cgalParamBegin{geom_traits}an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
   * \cgalNamedParamsEnd
   *
   * @return the surface area of `tmesh`.
@@ -321,6 +341,10 @@ namespace Polygon_mesh_processing {
   * or the geometric traits class deduced from the point property map
   * of `tmesh`.
   *
+  * \warning This function involves a square root computation.
+  * If `Kernel::FT` does not have a `sqrt()` operation, the square root computation
+  * will be done approximately.
+  *
   * @sa `face_area()`
   */
   template<typename TriangleMesh
@@ -359,7 +383,7 @@ namespace Polygon_mesh_processing {
   *
   * \cgalNamedParamsBegin
   *  \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd
-  *  \cgalParamBegin{geom_traits} a geometric traits class instance \cgalParamEnd
+  *  \cgalParamBegin{geom_traits}an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
   * \cgalNamedParamsEnd
   *
   * @return the volume bounded by `tmesh`.
@@ -392,7 +416,7 @@ namespace Polygon_mesh_processing {
 
     typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
 
-    double volume = 0.;
+    typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT volume = 0.;
     BOOST_FOREACH(face_descriptor f, faces(tmesh))
     {
       volume += CGAL::volume(origin,

Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp

@@ -22,6 +22,7 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 typedef K::Point_3 Point;
+typedef K::FT FT;
 
 typedef CGAL::Polyhedron_3<K>       Polyhedron;
 typedef CGAL::Surface_mesh<Point>   Surface_mesh;
@@ -43,7 +44,7 @@ void test(const Mesh& pmesh)
       break;
     }
   }
-  double border_l = PMP::face_border_length(border_he, pmesh);
+  FT border_l = PMP::face_border_length(border_he, pmesh);
   std::cout << "length of hole border = " << border_l << std::endl;
 
   face_descriptor valid_patch_face;
@@ -54,7 +55,7 @@ void test(const Mesh& pmesh)
       continue;
     else
     {
-      double face_area = PMP::face_area(face(h, pmesh), pmesh);
+      FT face_area = PMP::face_area(face(h, pmesh), pmesh);
       std::cout << "face area = " << face_area << std::endl;
 
       if(++count == 20)
@@ -81,13 +82,13 @@ void test(const Mesh& pmesh)
       break;//back to starting point
   }
 
-  double patch_area = PMP::area(patch, pmesh);
+  FT patch_area = PMP::area(patch, pmesh);
   std::cout << "patch area = " << patch_area << std::endl;
   
-  double mesh_area = PMP::area(pmesh);
+  FT mesh_area = PMP::area(pmesh);
   std::cout << "mesh area = " << mesh_area << std::endl;
   
-  double mesh_area_np = PMP::area(pmesh,
+  FT mesh_area_np = PMP::area(pmesh,
     PMP::parameters::geom_traits(K()));
   std::cout << "mesh area (NP) = " << mesh_area_np << std::endl;
 
@@ -131,7 +132,7 @@ void test_closed_surface_mesh(const char* filename)
 
   test(sm);
 
-  double vol = PMP::volume(sm);
+  FT vol = PMP::volume(sm);
   std::cout << "volume = " << vol << std::endl;
 }