diff --git a/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h index cb9349ba2a5..299e037a558 100644 --- a/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h +++ b/Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h @@ -281,6 +281,7 @@ namespace CGAL { typedef Advancing_front_surface_reconstruction Extract; typedef typename Triangulation_3::Geom_traits Geom_traits; + typedef typename Kernel::FT FT; typedef typename Kernel::FT coord_type; typedef typename Kernel::Point_3 Point; @@ -390,7 +391,23 @@ namespace CGAL { std::list nbe_pool; std::list ist_pool; + public: + Vector construct_vector(const Point& p, const Point& q) const + { + return T.geom_traits().construct_vector_3_object()(p, q); + } + Vector construct_cross_product(const Vector& v, const Vector& w) const + { + return T.geom_traits().construct_cross_product_vector_3_object()(v, w); + } + + FT compute_scalar_product(const Vector& v, const Vector& w) const + { + return T.geom_traits().compute_scalar_product_3_object()(v, w); + } + + private: Intern_successors_type* new_border() { nbe_pool.resize(nbe_pool.size()+1); @@ -692,12 +709,14 @@ namespace CGAL { ++it; }while(collinear(p,q,it->point())); const Point& r = it->point(); - Vector u = q-r; - Vector v = q-p; - Vector w = r-p; - Vector vw = cross_product(v,w); - double len = (std::max)(u*u,(std::max)(v*v,w*w)); - Point s = p + 10* len * (vw/(vw*vw)); + Vector u = construct_vector(r, q); + Vector v = construct_vector(p, q); + Vector w = construct_vector(p, r); + Vector vw = construct_cross_product(v,w); + double len = (std::max)(compute_scalar_product(u,u), + (std::max)(compute_scalar_product(v,v), + compute_scalar_product(w,w))); + Point s = p + 10 * len * (vw/compute_scalar_product(vw,vw)); added_vertex = T.insert(s); } } @@ -1199,7 +1218,7 @@ namespace CGAL { \param index index of the facet in `c` */ - coord_type + FT smallest_radius_delaunay_sphere(const Cell_handle& c, const int& index) const { @@ -1262,16 +1281,16 @@ namespace CGAL { const Point& pp2 = cc->vertex(i2)->point(); const Point& pp3 = cc->vertex(i3)->point(); - Sphere facet_sphere(pp1, pp2, pp3); - if (squared_distance(facet_sphere.center(), pp0) < - facet_sphere.squared_radius()) + Sphere facet_sphere = T.geom_traits().construct_sphere_3_object()(pp1, pp2, pp3); + if (squared_distance(T.geom_traits().construct_center_3_object()(facet_sphere), pp0) < + T.geom_traits().compute_squared_radius_3_object()(facet_sphere)) { #ifdef AFSR_LAZY value = lazy_squared_radius(cc); #else // qualified with CGAL, to avoid a compilation error with clang if(volume(pp0, pp1, pp2, pp3) != 0){ - value = CGAL::squared_radius(pp0, pp1, pp2, pp3); + value = T.geom_traits().compute_squared_radius_3_object()(pp0, pp1, pp2, pp3); } else { typedef Exact_predicates_exact_constructions_kernel EK; Cartesian_converter to_exact; @@ -1293,26 +1312,30 @@ namespace CGAL { cc = lazy_circumcenter(c); cn = lazy_circumcenter(n); #else - cc = CGAL::circumcenter(cp0, cp1, cp2, cp3); - cn = CGAL::circumcenter(np0, np1, np2, np3); + cc = T.geom_traits().construct_circumcenter_3_object()(cp0, cp1, cp2, cp3); + cn = T.geom_traits().construct_circumcenter_3_object()(np0, np1, np2, np3); #endif // computation of the distance of cp1 to the dual segment cc, cn... - Vector V(cc - cn), Vc(cc - cp1), Vn(cp1 - cn); - coord_type ac(V * Vc), an(V * Vn), norm_V(V * V); + Vector V = construct_vector(cn, cc), + Vc = construct_vector(cp1, cc), + Vn = construct_vector(cn, cp1); + coord_type ac = compute_scalar_product(V, Vc), + an = compute_scalar_product(V, Vn), + norm_V = compute_scalar_product(V, V); if ((ac > 0) && (an > 0)) { - value = (Vc*Vc) - ac*ac/norm_V; + value = compute_scalar_product(Vc, Vc) - ac*ac/norm_V; if ((value < 0)||(norm_V > inv_eps_2)){ // qualified with CGAL, to avoid a compilation error with clang - value = CGAL::squared_radius(cp1, cp2, cp3); + value = T.geom_traits().compute_squared_radius_3_object()(cp1, cp2, cp3); } } else { if (ac <= 0) - value = squared_distance(cc, cp1); + value = T.geom_traits().compute_squared_distance_3_object()(cc, cp1); else // (an <= 0) - value = squared_distance(cn, cp1); + value = T.geom_traits().compute_squared_distance_3_object()(cn, cp1); } } } @@ -1354,9 +1377,9 @@ namespace CGAL { const Point& p2 = c->vertex(i2)->point(); const Point& pc = c->vertex(i3)->point(); - Vector P2P1 = p1-p2, P2Pn, PnP1; + Vector P2P1 = construct_vector(p2, p1), P2Pn, PnP1; - Vector v2, v1 = cross_product(pc-p2, P2P1); + Vector v2, v1 = construct_cross_product(construct_vector(p2, pc), P2P1); coord_type norm, norm1 = v1*v1; coord_type norm12 = P2P1*P2P1; @@ -1388,12 +1411,12 @@ namespace CGAL { { const Point& pn = neigh->vertex(n_i3)->point(); - P2Pn = pn-p2; - v2 = cross_product(P2P1,P2Pn); + P2Pn = construct_vector(p2, pn); + v2 = construct_cross_product(P2P1,P2Pn); //pas necessaire de normer pour un bon echantillon: // on peut alors tester v1*v2 >= 0 - norm = sqrt(norm1 * (v2*v2)); + norm = sqrt(norm1 * compute_scalar_product(v2,v2)); pscal = v1*v2; // check if the triangle will produce a sliver on the surface bool sliver_facet = (pscal <= COS_ALPHA_SLIVER*norm); @@ -1407,10 +1430,9 @@ namespace CGAL { // We skip triangles having an internal angle along e // whose cosinus is smaller than -DELTA // that is the angle is larger than arcos(-DELTA) - border_facet = !((P2P1*P2Pn >= - -DELTA*sqrt(norm12*(P2Pn*P2Pn)))&& - (P2P1*PnP1 >= - -DELTA*sqrt(norm12*(PnP1*PnP1)))); + border_facet = + !((P2P1*P2Pn >= -DELTA*sqrt(norm12*compute_scalar_product(P2Pn,P2Pn))) && + (P2P1*PnP1 >= -DELTA*sqrt(norm12*compute_scalar_product(PnP1,PnP1)))); // \todo investigate why we simply do not skip this triangle // but continue looking for a better candidate // if (!border_facet){ @@ -1582,9 +1604,11 @@ namespace CGAL { int n_i3 = (6 - n_ind - n_i1 - n_i2); const Point& pn = neigh->vertex(n_i3)->point(); - Vector v1 = cross_product(pc-p2,p1-p2), - v2 = cross_product(p1-p2,pn-p2); - coord_type norm = sqrt((v1*v1)*(v2*v2)); + Vector v1 = construct_cross_product(construct_vector(p2, pc), + construct_vector(p2, p1)), + v2 = construct_cross_product(construct_vector(p2, p1), + construct_vector(p2, pn)); + coord_type norm = sqrt(compute_scalar_product(v1, v1) * compute_scalar_product(v2, v2)); if (v1*v2 > COS_BETA*norm) return 1; // label bonne pliure sinon: