diff --git a/Interpolation/examples/Interpolation/interpolation_2.cpp b/Interpolation/examples/Interpolation/interpolation_2.cpp index 6ac71fb497f..1a0501fa79b 100644 --- a/Interpolation/examples/Interpolation/interpolation_2.cpp +++ b/Interpolation/examples/Interpolation/interpolation_2.cpp @@ -1,4 +1,4 @@ -// compares the result of several interpolation methods +// Compares the result of several interpolation methods #include @@ -53,19 +53,17 @@ int main() Delaunay_triangulation T; - Point_value_map values; Point_vector_map gradients; - //parameters for quadratic function: Coord_type alpha = Coord_type(1.0), - beta1 = Coord_type(2.0), - beta2 = Coord_type(1.0), - gamma1 = Coord_type(0.3), - gamma2 = Coord_type(0.0), - gamma3 = Coord_type(0.0), - gamma4 = Coord_type(0.3); + beta1 = Coord_type(2.0), + beta2 = Coord_type(1.0), + gamma1 = Coord_type(0.3), + gamma2 = Coord_type(0.0), + gamma3 = Coord_type(0.0), + gamma4 = Coord_type(0.3); for(int j=0; j res; Coord_type error, l_total = Coord_type(0), - q_total(l_total), f_total(l_total), s_total(l_total), - ssquare_total(l_total), l_max(l_total), - q_max(l_total), f_max(l_total), s_max(l_total), - ssquare_max(l_total), - total_value(l_total), l_value(l_total); + q_total(l_total), f_total(l_total), s_total(l_total), + ssquare_total(l_total), l_max(l_total), + q_max(l_total), f_max(l_total), s_max(l_total), + ssquare_max(l_total), + total_value(l_total), l_value(l_total); int failure(0); //interpolation + error statistics @@ -98,22 +96,22 @@ int main() Coord_type y(points[i].y()); Coord_type exact_value = alpha + beta1*x + beta2*y + gamma1*(x*x) + - gamma4*(y*y) + (gamma2+ gamma3) *(x*y); + gamma4*(y*y) + (gamma2+ gamma3) *(x*y); total_value += exact_value; //Coordinate_vector: std::vector< std::pair< Point, Coord_type > > coords; Coord_type norm = - CGAL::natural_neighbor_coordinates_2(T, points[i], - std::back_inserter(coords)).second; + CGAL::natural_neighbor_coordinates_2(T, points[i], + std::back_inserter(coords)).second; assert(norm>0); //linear interpolant: - l_value = CGAL::linear_interpolation(coords.begin(), coords.end(), - norm, - CGAL::Data_access(values)); + l_value = CGAL::linear_interpolation(coords.begin(), coords.end(), + norm, + CGAL::Data_access(values)); error = CGAL_NTS abs(l_value - exact_value); l_total += error; @@ -121,42 +119,41 @@ int main() //Farin interpolant: res = CGAL::farin_c1_interpolation(coords.begin(), - coords.end(), norm,points[i], - CGAL::Data_access(values), - CGAL::Data_access - (gradients), - Traits()); + coords.end(), norm,points[i], + CGAL::Data_access(values), + CGAL::Data_access + (gradients), + Traits()); if(res.second){ error = CGAL_NTS abs(res.first - exact_value); f_total += error; if (error > f_max) f_max = error; - }else ++failure; - + } else ++failure; //quadratic interpolant: res = CGAL::quadratic_interpolation(coords.begin(), coords.end(), - norm,points[i], - CGAL::Data_access - (values), - CGAL::Data_access - (gradients), - Traits()); + norm,points[i], + CGAL::Data_access + (values), + CGAL::Data_access + (gradients), + Traits()); if(res.second){ error = CGAL_NTS abs(res.first - exact_value); q_total += error; if (error > q_max) q_max = error; - }else ++failure; + } else ++failure; //Sibson interpolant: version without sqrt: res = CGAL::sibson_c1_interpolation_square(coords.begin(), - coords.end(), norm, - points[i], - CGAL::Data_access - (values), - CGAL::Data_access - (gradients), - Traits()); - //error statistics + coords.end(), norm, + points[i], + CGAL::Data_access + (values), + CGAL::Data_access + (gradients), + Traits()); + //error statistics if(res.second){ error = CGAL_NTS abs(res.first - exact_value); ssquare_total += error; @@ -165,13 +162,13 @@ int main() //with sqrt(the traditional): res = CGAL::sibson_c1_interpolation(coords.begin(), - coords.end(), norm, - points[i], - CGAL::Data_access - (values), - CGAL::Data_access - (gradients), - Traits()); + coords.end(), norm, + points[i], + CGAL::Data_access + (values), + CGAL::Data_access + (gradients), + Traits()); //error statistics if(res.second){ @@ -184,30 +181,30 @@ int main() /************** end of Interpolation: dump statistics **************/ std::cout << "Result: -----------------------------------" << std::endl; std::cout << "Interpolation of '" << alpha <<" + " - << beta1<<" x + " - << beta2 << " y + " << gamma1 <<" x^2 + " << gamma2+ gamma3 - <<" xy + " << gamma4 << " y^2'" << std::endl; + << beta1<<" x + " + << beta2 << " y + " << gamma1 <<" x^2 + " << gamma2+ gamma3 + <<" xy + " << gamma4 << " y^2'" << std::endl; std::cout << "Knowing " << m << " sample points. Interpolation on " - << n <<" random points. "<< std::endl; + << n <<" random points. "<< std::endl; std::cout <<"Average function value " - << (1.0/n) * CGAL::to_double(total_value) - << ", nb of failures "<< failure << std::endl; + << (1.0/n) * CGAL::to_double(total_value) + << ", nb of failures "<< failure << std::endl; std::cout << "linear interpolant mean error " - << CGAL::to_double(l_total)/n << " max " - << CGAL::to_double(l_max) < function_values; typedef CGAL::Data_access< std::map > - Value_access; + Value_access; Coord_type a(0.25), bx(1.3), by(-0.7); - for (int y=0 ; y<3 ; y++) + for (int y=0 ; y<3 ; y++){ for (int x=0 ; x<3 ; x++){ K::Point_2 p(x,y); T.insert(p); function_values.insert(std::make_pair(p,a + bx* x+ by*y)); } + } + //coordinate computation K::Point_2 p(1.3,0.34); std::vector< std::pair< Point, Coord_type > > coords; Coord_type norm = - CGAL::natural_neighbor_coordinates_2 - (T, p,std::back_inserter(coords)).second; + CGAL::natural_neighbor_coordinates_2 + (T, p,std::back_inserter(coords)).second; Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), - norm, - Value_access(function_values)); + norm, + Value_access(function_values)); - std::cout << " Tested interpolation on " << p << " interpolation: " - << res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl; + std::cout << "Tested interpolation on " << p << " interpolation: " + << res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl; std::cout << "done" << std::endl; return 0; diff --git a/Interpolation/examples/Interpolation/nn_coordinates_2.cpp b/Interpolation/examples/Interpolation/nn_coordinates_2.cpp index 710237d8fb6..96e3d6567be 100644 --- a/Interpolation/examples/Interpolation/nn_coordinates_2.cpp +++ b/Interpolation/examples/Interpolation/nn_coordinates_2.cpp @@ -3,9 +3,8 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Delaunay_triangulation_2 Delaunay_triangulation; -typedef std::vector< std::pair< K::Point_2, K::FT > > - Point_coordinate_vector; +typedef CGAL::Delaunay_triangulation_2 Delaunay_triangulation; +typedef std::vector< std::pair< K::Point_2, K::FT > > Point_coordinate_vector; int main() { @@ -18,17 +17,17 @@ int main() //coordinate computation K::Point_2 p(1.2, 0.7); Point_coordinate_vector coords; - CGAL::Triple< - std::back_insert_iterator, - K::FT, bool> result = - CGAL::natural_neighbor_coordinates_2(dt, p, - std::back_inserter(coords)); + CGAL::Triple, + K::FT, bool> result = + CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords)); + if(!result.third){ std::cout << "The coordinate computation was not successful." - << std::endl; + << std::endl; std::cout << "The point (" <

> coor_laplace; - std::vector< std::pair< Vertex_handle,NT> > coor_sibson; - NT norm_coeff_laplace, norm_coeff_sibson; + for(int ii=0; ii<3; ++ii) + { + std::vector< std::pair< Vertex_handle,NT> > coor_laplace; + std::vector< std::pair< Vertex_handle,NT> > coor_sibson; + NT norm_coeff_laplace, norm_coeff_sibson; - std::cout << "Point P"<< ii+1 << " : "< Gt; typedef CGAL::Regular_triangulation_2 Regular_triangulation; typedef Regular_triangulation::Bare_point Bare_point; typedef Regular_triangulation::Weighted_point Weighted_point; -typedef std::vector< std::pair< Weighted_point, K::FT > > - Point_coordinate_vector; +typedef std::vector< std::pair< Weighted_point, K::FT > > Point_coordinate_vector; int main() { @@ -23,17 +22,17 @@ int main() //coordinate computation Weighted_point wp(Bare_point(1.2, 0.7),2); Point_coordinate_vector coords; - CGAL::Triple< - std::back_insert_iterator, - K::FT, bool> result = - CGAL::regular_neighbor_coordinates_2(rt, wp, - std::back_inserter(coords)); + CGAL::Triple, + K::FT, bool> result = + CGAL::regular_neighbor_coordinates_2(rt, wp, std::back_inserter(coords)); + if(!result.third){ std::cout << "The coordinate computation was not successful." - << std::endl; + << std::endl; std::cout << "The point (" < - (function_values), - Traits()); + } + sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients, + function_gradients.begin()), + CGAL::Data_access(function_values), + Traits()); //coordiante computation K::Point_2 p(1.6,1.4); std::vector< std::pair< Point, Coord_type > > coords; - Coord_type norm = - CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter - (coords)).second; + Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter + (coords)).second; //Sibson interpolant: version without sqrt: std::pair res = - CGAL::sibson_c1_interpolation_square - (coords.begin(), - coords.end(),norm,p, - CGAL::Data_access(function_values), - CGAL::Data_access(function_gradients), - Traits()); + CGAL::sibson_c1_interpolation_square( + coords.begin(), + coords.end(),norm,p, + CGAL::Data_access(function_values), + CGAL::Data_access(function_gradients), + Traits()); + if(res.second) - std::cout << " Tested interpolation on " << p - << " interpolation: " << res.first << " exact: " - << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y()) - << std::endl; + std::cout << "Tested interpolation on " << p + << " interpolation: " << res.first << " exact: " + << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y()) + << std::endl; else std::cout << "C^1 Interpolation not successful." << std::endl - << " not all function_gradients are provided." << std::endl - << " You may resort to linear interpolation." << std::endl; + << " not all function_gradients are provided." << std::endl + << " You may resort to linear interpolation." << std::endl; std::cout << "done" << std::endl; return 0; diff --git a/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp b/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp index 4c04c6cd7ef..0b533a66431 100644 --- a/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp +++ b/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp @@ -8,53 +8,51 @@ #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::FT Coord_type; -typedef K::Point_3 Point_3; -typedef K::Vector_3 Vector_3; -typedef std::vector< std::pair< Point_3, K::FT > > - Point_coordinate_vector; +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::FT Coord_type; +typedef K::Point_3 Point_3; +typedef K::Vector_3 Vector_3; +typedef std::vector< std::pair< Point_3, K::FT > > Point_coordinate_vector; int main() { - int n=100; std::vector< Point_3> points; points.reserve(n); - std::cout << "Generate " << n << " random points on a sphere." - << std::endl; + std::cout << "Generate " << n << " random points on a sphere." << std::endl; CGAL::Random_points_on_sphere_3 g(1); - CGAL::cpp11::copy_n( g, n, std::back_inserter(points)); + CGAL::cpp11::copy_n(g, n, std::back_inserter(points)); Point_3 p(1, 0,0); - Vector_3 normal(p-CGAL::ORIGIN); - std::cout << "Compute surface neighbor coordinates for " - << p << std::endl; + Vector_3 normal(p - CGAL::ORIGIN); + std::cout << "Compute surface neighbor coordinates for " << p << std::endl; Point_coordinate_vector coords; - CGAL::Triple< std::back_insert_iterator, - K::FT, bool> result = + CGAL::Triple, + K::FT, bool> result = CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(), - p, normal, - std::back_inserter(coords), - K()); + p, normal, + std::back_inserter(coords), + K()); if(!result.third){ //Undersampling: - std::cout << "The coordinate computation was not successful." - << std::endl; + std::cout << "The coordinate computation was not successful." << std::endl; return 0; } K::FT norm = result.second; std::cout << "Testing the barycentric property " << std::endl; - Point_3 b(0, 0,0); - for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator - it = coords.begin(); it!=coords.end(); ++it) + Point_3 b(0, 0, 0); + for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator + it = coords.begin(); it!=coords.end(); ++it) b = b + (it->second/norm)* (it->first - CGAL::ORIGIN); - std::cout <<" weighted barycenter: " << b < - #include namespace CGAL { @@ -44,14 +39,14 @@ public: Aff_transformation_2 operator()(const Aff_transformation_2& tr1, - const Aff_transformation_2& tr2) const + const Aff_transformation_2& tr2) const { return Aff_transformation_2(tr1.m(0,0) + tr2.m(0,0), - tr1.m(0,1) + tr2.m(0,1), - tr1.m(0,2) + tr2.m(0,2), - tr1.m(1,0) + tr2.m(1,0), - tr1.m(1,1) + tr2.m(1,1), - tr1.m(1,2) + tr2.m(1,2)); + tr1.m(0,1) + tr2.m(0,1), + tr1.m(0,2) + tr2.m(0,2), + tr1.m(1,0) + tr2.m(1,0), + tr1.m(1,1) + tr2.m(1,1), + tr1.m(1,2) + tr2.m(1,2)); } }; @@ -91,8 +86,9 @@ public: Aff_transformation_2 operator()(const Vector_2& v) const { - return Aff_transformation_2(v.x()*v.x(),v.x()*v.y(),v.x()*v.y(), - v.y()*v.y()); + return Aff_transformation_2(v.x()*v.x(), + v.x()*v.y(), v.x()*v.y(), + v.y()*v.y()); } }; @@ -116,41 +112,41 @@ public: typedef typename Rep::Aff_transformation_2 Aff_transformation_d; typedef Construct_null_matrix_2 - Construct_null_matrix_d; + Construct_null_matrix_d; typedef Construct_scaling_matrix_2 - Construct_scaling_matrix_d; + Construct_scaling_matrix_d; typedef Construct_sum_matrix_2 Construct_sum_matrix_d; typedef Construct_outer_product_2 Construct_outer_product_d; Construct_outer_product_d construct_outer_product_d_object() const - {return Construct_outer_product_d();} + {return Construct_outer_product_d();} Construct_sum_matrix_d construct_sum_matrix_d_object() const - {return Construct_sum_matrix_d();} + {return Construct_sum_matrix_d();} Construct_scaling_matrix_d construct_scaling_matrix_d_object() const - {return Construct_scaling_matrix_d();} + {return Construct_scaling_matrix_d();} Construct_null_matrix_d construct_null_matrix_d_object() const - {return Construct_null_matrix_d();} + {return Construct_null_matrix_d();} //also in the traits without gradient computation: Construct_scaled_vector_d construct_scaled_vector_d_object()const - {return Construct_scaled_vector_d();} + {return Construct_scaled_vector_d();} Construct_vector_d construct_vector_d_object()const - {return Construct_vector_d();} + {return Construct_vector_d();} Compute_squared_distance_d compute_squared_distance_d_object()const - {return Compute_squared_distance_d();} + {return Compute_squared_distance_d();} }; } //namespace CGAL diff --git a/Interpolation/include/CGAL/Interpolation_traits_2.h b/Interpolation/include/CGAL/Interpolation_traits_2.h index a5f8b80e0a3..51e90c0e169 100644 --- a/Interpolation/include/CGAL/Interpolation_traits_2.h +++ b/Interpolation/include/CGAL/Interpolation_traits_2.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #ifndef CGAL_INTERPOLATION_TRAITS_2_H @@ -23,7 +19,6 @@ #include - namespace CGAL { //-----------------------------------------------------------------------// diff --git a/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h b/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h index bde832500f5..e892b6b2e39 100644 --- a/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h +++ b/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2003 INRIA Sophia-Antipolis (France). +// Copyright (c) 2003, 2017 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -19,7 +19,6 @@ #include - #include #include #include @@ -116,7 +115,6 @@ private: const Vector& normal; }; - template < typename K > class Construct_plane_intersected_bisector_3 { @@ -142,7 +140,6 @@ private: const Vector& normal; }; - template < typename K > class Compare_first_projection_3 { diff --git a/Interpolation/include/CGAL/constructions/constructions_for_voronoi_intersection_cartesian_2_3.h b/Interpolation/include/CGAL/constructions/constructions_for_voronoi_intersection_cartesian_2_3.h index 3c8bd99843b..f1cf4ec750b 100644 --- a/Interpolation/include/CGAL/constructions/constructions_for_voronoi_intersection_cartesian_2_3.h +++ b/Interpolation/include/CGAL/constructions/constructions_for_voronoi_intersection_cartesian_2_3.h @@ -29,50 +29,50 @@ namespace CGAL { template < class RT> void plane_centered_circumcenter_translateC3(const RT &ax, const RT &ay, - const RT &az, - const RT &nx, const RT &ny, - const RT &nz, - const RT &qx, const RT &qy, - const RT &qz, - const RT &rx, const RT &ry, - const RT &rz, - RT &x, RT &y, RT &z) + const RT &az, + const RT &nx, const RT &ny, + const RT &nz, + const RT &qx, const RT &qy, + const RT &qz, + const RT &rx, const RT &ry, + const RT &rz, + RT &x, RT &y, RT &z) { RT den = RT(2) * determinant(nx,qx,rx, - ny,qy,ry, - nz,qz,rz); + ny,qy,ry, + nz,qz,rz); // The 3 points aren't collinear. // Hopefully, this is already checked at the upper level. CGAL_assertion ( den != RT(0) ); RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) + - CGAL_NTS square(qz); + CGAL_NTS square(qz); RT r2 = CGAL_NTS square(rx) + CGAL_NTS square(ry) + - CGAL_NTS square(rz); + CGAL_NTS square(rz); RT na = nx*ax + ny*ay + nz*az; na *= RT(2.0); x = determinant(ny,nz,na, - qy,qz,q2, - ry,rz,r2)/ den ; + qy,qz,q2, + ry,rz,r2)/ den ; y = - determinant(nx,nz,na, - qx,qz,q2, - rx,rz,r2)/ den ; + qx,qz,q2, + rx,rz,r2)/ den ; z = determinant(nx,ny,na, - qx,qy,q2, - rx,ry,r2)/ den ; + qx,qy,q2, + rx,ry,r2)/ den ; } template < class RT> void plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &px, const RT &py, const RT &pz, - const RT &qx, const RT &qy, const RT &qz, - const RT &rx, const RT &ry, const RT &rz, - RT &x, RT &y, RT &z) + const RT &nx, const RT &ny, const RT &nz, + const RT &px, const RT &py, const RT &pz, + const RT &qx, const RT &qy, const RT &qz, + const RT &rx, const RT &ry, const RT &rz, + RT &x, RT &y, RT &z) { // resolution of the system (where we note c the center) // @@ -83,10 +83,10 @@ plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az, //method: // - tranlation of p to the origin. plane_centered_circumcenter_translateC3(ax-px, ay-py, az-pz, - nx, ny, nz, - qx-px, qy-py,qz-pz, - rx-px, ry-py,rz-pz, - x, y, z); + nx, ny, nz, + qx-px, qy-py,qz-pz, + rx-px, ry-py,rz-pz, + x, y, z); x+=px; y+=py; z+=pz; @@ -95,13 +95,13 @@ plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az, template < class RT> void bisector_plane_intersection_translateC3(const RT &ax, const RT &ay, - const RT &az, - const RT &nx, const RT &ny, - const RT &nz, - const RT &qx, const RT &qy, - const RT &qz, const RT& den, - RT &x1, RT &y1, RT &x2, RT - &y2, bool& swapped) + const RT &az, + const RT &nx, const RT &ny, + const RT &nz, + const RT &qx, const RT &qy, + const RT &qz, const RT& den, + RT &x1, RT &y1, RT &x2, RT + &y2, bool& swapped) { // c: a point on l must be the center of a sphere passing // through p and q, c lies in h. 2 equations: @@ -113,7 +113,7 @@ bisector_plane_intersection_translateC3(const RT &ax, const RT &ay, // where RT den = RT(2.0) * determinant(qx,qy,nx, ny); RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) - + CGAL_NTS square(qz); + + CGAL_NTS square(qz); RT na = nx*ax + ny*ay + nz*az; na *= RT(2.0); @@ -134,39 +134,39 @@ bisector_plane_intersection_translateC3(const RT &ax, const RT &ay, // if not: permutation of p1 and p2 if((sign_of_determinant(qx,qy,qz, x1,y1,RT(0),x2 ,y2,RT(1)) * CGAL_NTS sign (-na)) > 0 ) - { - RT x3(x1),y3(y1); - x1 =x2; - y1 =y2; - x2 = x3; - y2 = y3; - swapped =true; - } + { + RT x3(x1),y3(y1); + x1 =x2; + y1 =y2; + x2 = x3; + y2 = y3; + swapped =true; + } } template < class RT> void bisector_plane_intersection_permuteC3(const RT &ax, const RT &ay, - const RT &az, - const RT &nx, const RT &ny, - const RT &nz, - const RT &px, const RT &py, - const RT &pz, - const RT &qx, const RT &qy, - const RT &qz, - const RT &den, - RT &x1, RT &y1, RT& z1, - RT &x2, RT &y2, RT& z2) + const RT &az, + const RT &nx, const RT &ny, + const RT &nz, + const RT &px, const RT &py, + const RT &pz, + const RT &qx, const RT &qy, + const RT &qz, + const RT &den, + RT &x1, RT &y1, RT& z1, + RT &x2, RT &y2, RT& z2) { //translation of p to the origin bool swapped =false; CGAL_precondition((nx!=RT(0) || ny!=RT(0)) && (qx!=px || qy!=py) - &&den!=RT(0)); + &&den!=RT(0)); bisector_plane_intersection_translateC3(ax-px, ay-py, az-pz, - nx, ny, nz, - qx-px, qy-py,qz-pz,den, - x1, y1,x2,y2,swapped); + nx, ny, nz, + qx-px, qy-py,qz-pz,den, + x1, y1,x2,y2,swapped); // re-translation of the origin to p: x1+=px; y1+=py; @@ -183,11 +183,11 @@ bisector_plane_intersection_permuteC3(const RT &ax, const RT &ay, template < class RT> void bisector_plane_intersectionC3(const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &px, const RT &py, const RT &pz, - const RT &qx, const RT &qy, const RT &qz, - RT &x1, RT &y1, RT& z1, - RT &x2, RT &y2, RT& z2) + const RT &nx, const RT &ny, const RT &nz, + const RT &px, const RT &py, const RT &pz, + const RT &qx, const RT &qy, const RT &qz, + RT &x1, RT &y1, RT& z1, + RT &x2, RT &y2, RT& z2) { // constructs the line l = (p1,p2)= ((x1,y1,z1),(x2,y2,z2)) // the intersection line between the bisector of (p,q) and @@ -208,22 +208,22 @@ bisector_plane_intersectionC3(const RT &ax, const RT &ay, const RT &az, //den==0 <=> projections of (qx,qy) and (nx,ny) are identical //intersection with z=0/z=1 bisector_plane_intersection_permuteC3(ax,ay,az,nx,ny,nz,px,py,pz, - qx,qy,qz,den, - x1,y1,z1,x2,y2,z2); + qx,qy,qz,den, + x1,y1,z1,x2,y2,z2); else{ den = RT(2.0) * determinant(qy-py,qz-pz,ny,nz); if ((ny!=0 || nz!=0) && (qy!=py || qz!=pz) && den!=RT(0)) //intersection with x=0/x=1 => permutations bisector_plane_intersection_permuteC3(ay,az,ax,ny,nz,nx,py,pz,px, - qy,qz,qx,den, - y1,z1,x1,y2,z2,x2); + qy,qz,qx,den, + y1,z1,x1,y2,z2,x2); else{ den = RT(2.0) * determinant(qz-pz,qx-px,nz,nx); CGAL_assertion((nx!=0 || nz!=0) && (qx!=px || qz!=pz) && den!=RT(0)); //intersection with y=0/y=1 => permutations bisector_plane_intersection_permuteC3(az,ax,ay,nz,nx,ny,pz,px,py, - qz,qx,qy,den, - z1,x1,y1,z2,x2,y2); + qz,qx,qy,den, + z1,x1,y1,z2,x2,y2); } } } diff --git a/Interpolation/include/CGAL/interpolation_functions.h b/Interpolation/include/CGAL/interpolation_functions.h index 278b5269029..bf3ab51ca90 100644 --- a/Interpolation/include/CGAL/interpolation_functions.h +++ b/Interpolation/include/CGAL/interpolation_functions.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #ifndef CGAL_INTERPOLATION_FUNCTIONS_H @@ -23,24 +19,25 @@ #include - -#include #include #include +#include +#include #include namespace CGAL { //Functor class for accessing the function values/gradients template< class Map > -struct Data_access : public std::unary_function< typename Map::key_type, - std::pair< typename Map::mapped_type, bool> > +struct Data_access + : public std::unary_function< typename Map::key_type, + std::pair< typename Map::mapped_type, bool> > { typedef typename Map::mapped_type Data_type; typedef typename Map::key_type Key_type; - Data_access< Map >(const Map& m): map(m){}; + Data_access< Map >(const Map& m): map(m){} std::pair< Data_type, bool> operator()(const Key_type& p) const { @@ -48,7 +45,7 @@ struct Data_access : public std::unary_function< typename Map::key_type, if(mit!= map.end()) return std::make_pair(mit->second, true); return std::make_pair(Data_type(), false); - }; + } const Map& map; }; @@ -57,9 +54,10 @@ struct Data_access : public std::unary_function< typename Map::key_type, template < class ForwardIterator, class Functor> typename Functor::result_type::first_type linear_interpolation(ForwardIterator first, ForwardIterator beyond, - const typename - std::iterator_traits::value_type:: - second_type& norm, Functor function_value) + const typename + std::iterator_traits::value_type:: + second_type& norm, + Functor function_value) { CGAL_precondition(norm>0); typedef typename Functor::result_type::first_type Value_type; @@ -74,17 +72,18 @@ linear_interpolation(ForwardIterator first, ForwardIterator beyond, } -template < class ForwardIterator, class Functor, class GradFunctor, - class Traits> +template < class ForwardIterator, class Functor, class GradFunctor, class Traits> std::pair< typename Functor::result_type::first_type, bool> quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, - const typename - std::iterator_traits:: - value_type::second_type& norm, const typename - std::iterator_traits::value_type:: - first_type& p, Functor function_value, - GradFunctor function_gradient, - const Traits& traits) + const typename + std::iterator_traits:: + value_type::second_type& norm, + const typename + std::iterator_traits::value_type:: + first_type& p, + Functor function_value, + GradFunctor function_gradient, + const Traits& traits) { CGAL_precondition(norm >0); typedef typename Functor::result_type::first_type Value_type; @@ -99,27 +98,26 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, if(!grad.second) return std::make_pair(Value_type(0), false); result += (first->second/norm) - *( f.first + grad.first* - traits.construct_scaled_vector_d_object() - (traits.construct_vector_d_object()(first->first, p),0.5)); + *( f.first + grad.first* + traits.construct_scaled_vector_d_object() + (traits.construct_vector_d_object()(first->first, p),0.5)); } return std::make_pair(result, true); } -template < class ForwardIterator, class Functor, class GradFunctor, - class Traits> +template < class ForwardIterator, class Functor, class GradFunctor, class Traits> std::pair< typename Functor::result_type::first_type, bool> sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, - const typename - std::iterator_traits:: - value_type::second_type& - norm, const typename - std::iterator_traits::value_type:: - first_type& p, - Functor function_value, - GradFunctor function_gradient, - const Traits& traits) + const typename + std::iterator_traits:: + value_type::second_type& norm, + const typename + std::iterator_traits::value_type:: + first_type& p, + Functor function_value, + GradFunctor function_gradient, + const Traits& traits) { CGAL_precondition(norm >0); typedef typename Functor::result_type::first_type Value_type; @@ -140,7 +138,7 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, Coord_type coeff = first->second/norm; Coord_type squared_dist = traits. - compute_squared_distance_d_object()(first->first, p); + compute_squared_distance_d_object()(first->first, p); Coord_type dist = CGAL_NTS sqrt(squared_dist); if(squared_dist ==0){ @@ -151,22 +149,22 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, } //three different terms to mix linear and gradient //interpolation - term1 += coeff/dist; + term1 += coeff / dist; term2 += coeff * squared_dist; term3 += coeff * dist; linear_int += coeff * f.first; gradient_int += (coeff/dist) - * (f.first + grad.first * - traits.construct_vector_d_object()(first->first, p)); + * (f.first + grad.first * + traits.construct_vector_d_object()(first->first, p)); } term4 = term3/ term1; gradient_int = gradient_int / term1; return std::make_pair((term4* linear_int + term2 * gradient_int)/ - (term4 + term2), true); + (term4 + term2), true); } //this method works with rational number types: @@ -182,19 +180,18 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, // (vh->get_value()+ vh->get_gradient() // *(p - vh->point())); -template < class ForwardIterator, class Functor, class GradFunctor, - class Traits> +template < class ForwardIterator, class Functor, class GradFunctor, class Traits> std::pair< typename Functor::result_type::first_type, bool> -sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator - beyond, const typename - std::iterator_traits:: - value_type::second_type& norm, - const typename - std::iterator_traits:: - value_type::first_type& p, - Functor function_value, - GradFunctor function_gradient, - const Traits& traits) +sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond, + const typename + std::iterator_traits:: + value_type::second_type& norm, + const typename + std::iterator_traits:: + value_type::first_type& p, + Functor function_value, + GradFunctor function_gradient, + const Traits& traits) { CGAL_precondition(norm >0); typedef typename Functor::result_type::first_type Value_type; @@ -205,7 +202,7 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator typename Functor::result_type f; typename GradFunctor::result_type grad; - for(; first !=beyond; ++first){ + for(; first!=beyond; ++first){ f = function_value(first->first); grad = function_gradient(first->first); CGAL_assertion(f.second); @@ -215,7 +212,7 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator Coord_type coeff = first->second/norm; Coord_type squared_dist = traits. - compute_squared_distance_d_object()(first->first, p); + compute_squared_distance_d_object()(first->first, p); if(squared_dist ==0){ ForwardIterator it = first; @@ -225,38 +222,38 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator } //three different terms to mix linear and gradient //interpolation - term1 += coeff/squared_dist; + term1 += coeff / squared_dist; term2 += coeff * squared_dist; term3 += coeff; linear_int += coeff * f.first; - gradient_int += (coeff/squared_dist) - *(f.first + grad.first* - traits.construct_vector_d_object()(first->first, p)); + gradient_int += (coeff/squared_dist) * (f.first + grad.first * + traits.construct_vector_d_object()(first->first, p)); } term4 = term3/ term1; gradient_int = gradient_int / term1; - return std::make_pair((term4* linear_int + term2 * gradient_int)/ - (term4 + term2), true); + return std::make_pair((term4 * linear_int + term2 * gradient_int)/ + (term4 + term2), true); } template < class RandomAccessIterator, class Functor, class -GradFunctor, class Traits> + GradFunctor, class Traits> std::pair< typename Functor::result_type::first_type, bool> farin_c1_interpolation(RandomAccessIterator first, - RandomAccessIterator beyond, - const typename - std::iterator_traits:: - value_type::second_type& norm, const typename - std::iterator_traits:: - value_type::first_type& /*p*/, - Functor function_value, GradFunctor - function_gradient, - const Traits& traits) + RandomAccessIterator beyond, + const typename + std::iterator_traits:: + value_type::second_type& norm, + const typename + std::iterator_traits:: + value_type::first_type& /*p*/, + Functor function_value, GradFunctor + function_gradient, + const Traits& traits) { CGAL_precondition(norm >0); //the function value is available for all points @@ -283,7 +280,7 @@ farin_c1_interpolation(RandomAccessIterator first, const Coord_type fac3(3); std::vector< std::vector > - ordinates(n,std::vector(n, Value_type(0))); + ordinates(n,std::vector(n, Value_type(0))); for(int i =0; ifirst); - if(!grad.second) - //the gradient is not known - return std::make_pair(Value_type(0), false); + grad = function_gradient(it->first); + if(!grad.second) + //the gradient is not known + return std::make_pair(Value_type(0), false); - //ordinates[i][j] = (p_j - p_i) * g_i - ordinates[i][j] = grad.first * - traits.construct_vector_d_object()(it->first,it2->first); + //ordinates[i][j] = (p_j - p_i) * g_i + ordinates[i][j] = grad.first * + traits.construct_vector_d_object()(it->first,it2->first); - // a point in the tangent plane: - // 3( f(p_i) + (1/3)(p_j - p_i) * g_i) - // => 3*f(p_i) + (p_j - p_i) * g_i - res_i += (fac3 * ordinates[i][i] + ordinates[i][j])* it2->second; + // a point in the tangent plane: + // 3( f(p_i) + (1/3)(p_j - p_i) * g_i) + // => 3*f(p_i) + (p_j - p_i) * g_i + res_i += (fac3 * ordinates[i][i] + ordinates[i][j])* it2->second; } } //res_i already multiplied by three @@ -326,20 +323,20 @@ farin_c1_interpolation(RandomAccessIterator first, for(int i=0; i< n; ++i) for(int j=i+1; j< n; ++j) for(int k=j+1; k 6 * b_ijk = 3*(f_i + f_j + f_k) + 0.5*a - result += (Coord_type(2.0)*( ordinates[i][i]+ ordinates[j][j]+ - ordinates[k][k]) - + Coord_type(0.5)*(ordinates[i][j] + ordinates[i][k] - + ordinates[j][i] + - ordinates[j][k] + ordinates[k][i]+ - ordinates[k][j])) - *(first+i)->second *(first+j)->second *(first+k)->second ; + // add 6* (u_i*u_j*u_k) * b_ijk + // b_ijk = 1.5 * a - 0.5*c + //where + //c : average of the three data control points + //a : 1.5*a = 1/12 * (ord[i][j] + ord[i][k] + ord[j][i] + + // ord[j][k] + ord[k][i]+ ord[k][j]) + // => 6 * b_ijk = 3*(f_i + f_j + f_k) + 0.5*a + result += (Coord_type(2.0)*( ordinates[i][i]+ ordinates[j][j]+ + ordinates[k][k]) + + Coord_type(0.5)*(ordinates[i][j] + ordinates[i][k] + + ordinates[j][i] + + ordinates[j][k] + ordinates[k][i]+ + ordinates[k][j])) + *(first+i)->second *(first+j)->second *(first+k)->second ; } return std::make_pair(result/(CGAL_NTS square(norm)*norm), true); diff --git a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h index 2cba5abf2c7..6e724fa643f 100644 --- a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h +++ b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Frank Da, Julia Floetotto #ifndef CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H @@ -23,13 +19,14 @@ #include - -#include #include #include #include #include +#include +#include + namespace CGAL { // the struct "Project_vertex_output_iterator" @@ -46,8 +43,8 @@ struct Project_vertex_output_iterator // std::pair // into an output iterator with value type std::pair // Conditions: OutputIterator has value type std::pair - // and Vertex_handle has a function ->point() - // with return type const Point& + // and Vertex_handle has a function ->point() + // with return type const Point& OutputIterator _base; @@ -75,29 +72,27 @@ struct Project_vertex_output_iterator // template // Triple< OutputIterator, typename Traits::FT, bool > // natural_neighbor_coordinates_2(const Dt& dt, -// const typename Traits::Point_2& p, -// OutputIterator out, const Traits& traits, -// typename Dt::Face_handle start -// = typename Dt::Face_handle()) -// +// const typename Traits::Point_2& p, +// OutputIterator out, const Traits& traits, +// typename Dt::Face_handle start +// = typename Dt::Face_handle()) + //template //Triple< OutputIterator, typename Traits::FT, bool > //natural_neighbor_coordinates_2(const Dt& dt, -// typename Dt::Vertex_handle vh, -// OutputIterator out, const Traits& traits) - +// typename Dt::Vertex_handle vh, +// OutputIterator out, const Traits& traits) //the following two functions suppose that // OutputIterator has value type -// std::pair +// std::pair //!!!they are not documented!!! -template +template < class Dt, class OutputIterator > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > natural_neighbor_coordinates_vertex_2(const Dt& dt, - const typename Dt::Geom_traits::Point_2& p, - OutputIterator out, typename Dt::Face_handle start - = typename Dt::Face_handle()) - + const typename Dt::Geom_traits::Point_2& p, + OutputIterator out, typename Dt::Face_handle start + = typename Dt::Face_handle()) { typedef typename Dt::Geom_traits Traits; typedef typename Traits::FT Coord_type; @@ -109,47 +104,43 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt, typedef typename Traits::Equal_x_2 Equal_x_2; CGAL_precondition(dt.dimension() == 2); - + Locate_type lt; int li; Face_handle fh = dt.locate(p, lt, li, start); - if (lt == Dt::OUTSIDE_AFFINE_HULL - || lt == Dt::OUTSIDE_CONVEX_HULL) + if (lt == Dt::OUTSIDE_AFFINE_HULL || lt == Dt::OUTSIDE_CONVEX_HULL) { return make_triple(out, Coord_type(1), false); } - + if ((lt == Dt::EDGE && - (dt.is_infinite(fh) || - dt.is_infinite(fh->neighbor(li))))) - { - Vertex_handle v1 = fh->vertex(dt.cw(li)); - Vertex_handle v2 = fh->vertex(dt.ccw(li)); + (dt.is_infinite(fh) || dt.is_infinite(fh->neighbor(li))))) + { + Vertex_handle v1 = fh->vertex(dt.cw(li)); + Vertex_handle v2 = fh->vertex(dt.ccw(li)); - Point_2 p1(v1->point()),p2(v2->point()); - - Coord_type coef1(0); - Coord_type coef2(0); - Equal_x_2 equal_x_2; - if(!equal_x_2(p1,p2)) - { - coef1 = (p.x() - p2.x())/(p1.x() - p2.x()) ; - coef2 = 1-coef1; - *out++= std::make_pair(v1,coef1); - *out++= std::make_pair(v2,coef2); + Point_2 p1(v1->point()), p2(v2->point()); - }else{ - coef1 = (p.y() - p2.y())/(p1.y() - p2.y()) ; - coef2 = 1-coef1; - *out++= std::make_pair(v1,coef1); - *out++= std::make_pair(v2,coef2); - - } - - return make_triple(out, coef1+coef2, true); + Coord_type coef1(0); + Coord_type coef2(0); + Equal_x_2 equal_x_2; + if(!equal_x_2(p1,p2)) + { + coef1 = (p.x() - p2.x()) / (p1.x() - p2.x()); + coef2 = 1 - coef1; + *out++ = std::make_pair(v1,coef1); + *out++ = std::make_pair(v2,coef2); + } else { + coef1 = (p.y() - p2.y()) / (p1.y() - p2.y()); + coef2 = 1-coef1; + *out++ = std::make_pair(v1,coef1); + *out++ = std::make_pair(v2,coef2); + } + + return make_triple(out, coef1+coef2, true); } - + if (lt == Dt::VERTEX) { *out++= std::make_pair(fh->vertex(li), Coord_type(1)); @@ -157,22 +148,21 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt, } std::list hole; - dt.get_boundary_of_conflicts(p, std::back_inserter(hole), fh, false); return natural_neighbor_coordinates_vertex_2 - (dt, p, out, hole.begin(), hole.end()); + (dt, p, out, hole.begin(), hole.end()); } //function call if the conflict zone is known: // OutputIterator has value type // std::pair -template +template < class Dt, class OutputIterator, class EdgeIterator > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > natural_neighbor_coordinates_vertex_2(const Dt& dt, - const typename Dt::Geom_traits::Point_2& p, - OutputIterator out, EdgeIterator - hole_begin, EdgeIterator hole_end) + const typename Dt::Geom_traits::Point_2& p, + OutputIterator out, EdgeIterator + hole_begin, EdgeIterator hole_end) { CGAL_precondition(dt.dimension()==2); //precondition: p must lie inside the hole @@ -184,7 +174,6 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt, typedef typename Dt::Vertex_handle Vertex_handle; typedef typename Dt::Face_circulator Face_circulator; - std::vector vor(3); Coord_type area_sum(0); @@ -196,42 +185,39 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt, hit = hole_begin; while (hit != hole_end) + { + Coord_type area(0); + Vertex_handle current = hit->first->vertex(dt.cw(hit->second)); + + vor[0] = dt.geom_traits().construct_circumcenter_2_object()( + current->point(), + hit->first->vertex(dt.ccw(hit->second))->point(), + p); + + Face_circulator fc = dt.incident_faces(current, hit->first); + ++fc; + vor[1] = dt.dual(fc); + + while(!fc->has_vertex(prev)) { - Coord_type area(0); - Vertex_handle current = hit->first->vertex(dt.cw(hit->second)); - - vor[0] = dt.geom_traits().construct_circumcenter_2_object() - (current->point(), - hit->first->vertex(dt.ccw(hit->second))->point(), - p); - - Face_circulator fc = dt.incident_faces(current, hit->first); ++fc; - vor[1] = dt.dual(fc); - - while(!fc->has_vertex(prev)) - { - ++fc; - vor[2] = dt.dual(fc); - - area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits()); - - vor[1] = vor[2]; - }; - vor[2] = - dt.geom_traits().construct_circumcenter_2_object()(prev->point(), - current->point(),p); - area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits()); - - - *out++= std::make_pair(current,area); - area_sum += area; - - - //update prev and hit: - prev= current; - ++hit; + vor[2] = dt.dual(fc); + area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits()); + vor[1] = vor[2]; } + + vor[2] = dt.geom_traits().construct_circumcenter_2_object()(prev->point(), + current->point(), + p); + + area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits()); + *out++ = std::make_pair(current,area); + area_sum += area; + + //update prev and hit: + prev = current; + ++hit; + } return make_triple(out, area_sum, true); } @@ -242,23 +228,22 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt, //=> OutputIterator has value type // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> ///////////////////////////////////////////////////////////// -template +template < class Dt, class OutputIterator > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > natural_neighbor_coordinates_2(const Dt& dt, - const typename Dt::Geom_traits::Point_2& p, - OutputIterator out, - typename Dt::Face_handle start = - CGAL_TYPENAME_DEFAULT_ARG Dt::Face_handle() ) + const typename Dt::Geom_traits::Point_2& p, + OutputIterator out, + typename Dt::Face_handle start = + CGAL_TYPENAME_DEFAULT_ARG Dt::Face_handle() ) { CGAL_precondition(dt.dimension() == 2); - + Project_vertex_output_iterator op(out); - Triple< Project_vertex_output_iterator, - typename Dt::Geom_traits::FT, bool > result = - natural_neighbor_coordinates_vertex_2 - (dt, p, op, start); + Triple, + typename Dt::Geom_traits::FT, bool > result = + natural_neighbor_coordinates_vertex_2(dt, p, op, start); return make_triple(result.first.base(), result.second, result.third); } @@ -266,21 +251,20 @@ natural_neighbor_coordinates_2(const Dt& dt, //OutputIterator has value type // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> //function call if the conflict zone is known: -template +template < class Dt, class OutputIterator, class EdgeIterator > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > natural_neighbor_coordinates_2(const Dt& dt, - const typename Dt::Geom_traits::Point_2& p, - OutputIterator out, EdgeIterator - hole_begin, EdgeIterator hole_end) + const typename Dt::Geom_traits::Point_2& p, + OutputIterator out, EdgeIterator + hole_begin, EdgeIterator hole_end) { CGAL_precondition(dt.dimension() == 2); - + Project_vertex_output_iterator op(out); - Triple< Project_vertex_output_iterator, - typename Dt::Geom_traits::FT, bool > result = - natural_neighbor_coordinates_vertex_2 - (dt, p, op, hole_begin,hole_end); + Triple, + typename Dt::Geom_traits::FT, bool > result = + natural_neighbor_coordinates_vertex_2(dt, p, op, hole_begin,hole_end); return make_triple(result.first.base(), result.second, result.third); } @@ -293,8 +277,8 @@ natural_neighbor_coordinates_2(const Dt& dt, template Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > natural_neighbor_coordinates_2(const Dt& dt, - typename Dt::Vertex_handle vh, - OutputIterator out) + typename Dt::Vertex_handle vh, + OutputIterator out) { //this functions creates a small triangulation of the // incident vertices of this vertex and computes the @@ -302,29 +286,29 @@ natural_neighbor_coordinates_2(const Dt& dt, typedef typename Dt::Vertex_circulator Vertex_circulator; CGAL_precondition(dt.dimension() == 2); - + Dt t2; Vertex_circulator vc = dt.incident_vertices(vh), - done(vc); + done(vc); do{ CGAL_assertion(!dt.is_infinite(vc)); t2.insert(vc->point()); } while(++vc!=done); return natural_neighbor_coordinates_2(t2, vh->point(), out); - } +} //class providing a function object: //OutputIterator has value type // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> -template +template < class Dt, class OutputIterator > class natural_neighbor_coordinates_2_object { public: Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > operator()(const Dt& dt, - typename Dt::Vertex_handle vh, - OutputIterator out) const + typename Dt::Vertex_handle vh, + OutputIterator out) const { return natural_neighbor_coordinates_2(dt, vh, out); } diff --git a/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h b/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h index 0f8f6fc1ab9..34f2d44c32f 100644 --- a/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h +++ b/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h @@ -23,16 +23,18 @@ #include - -#include -#include #include #include #include #include #include +#include #include //TO DO : to remove +#include +#include +#include +#include namespace CGAL { @@ -57,24 +59,25 @@ construct_circumcenter(const typename DT::Facet& f, // ====================== Natural Neighbors Querries ========================== // === Definitions -// Given a 3D point Q and a 3D Delaunay triangulation dt, +// Given a 3D point Q and a 3D Delaunay triangulation dt, // the next two functions calculate the natural neighbors and coordinates of Q with regard of dt -// +// // OutputIterator has value type // std::pair -// Result : +// Result : // - An OutputIterator providing natural neighbors P_i of Q with unnormalized coordinates a_i associated to them // - The normalizing coefficient (sum over i of the a_i) // - A boolean specifying whether the calculation has succeeded or not template Triple< OutputIterator, // iterator with value type std::pair - typename Dt::Geom_traits::FT, // Should provide 0 and 1 - bool > + typename Dt::Geom_traits::FT, // Should provide 0 and 1 + bool > laplace_natural_neighbor_coordinates_3(const Dt& dt, - const typename Dt::Geom_traits::Point_3& Q, - OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff, - const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) + const typename Dt::Geom_traits::Point_3& Q, + OutputIterator nn_out, + typename Dt::Geom_traits::FT& norm_coeff, + const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) { typedef typename Dt::Geom_traits Gt; typedef typename Gt::Point_3 Point; @@ -84,29 +87,32 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt, typedef typename Dt::Locate_type Locate_type; typedef typename Gt::FT Coord_type; - CGAL_triangulation_precondition (dt.dimension()== 3); - - Locate_type lt; int li, lj; + CGAL_triangulation_precondition (dt.dimension() == 3); + Locate_type lt; + int li, lj; Cell_handle c = dt.locate( Q, lt, li, lj, start); if ( lt == Dt::VERTEX ) - { - *nn_out++= std::make_pair(c->vertex(li),Coord_type(1)); - return make_triple(nn_out,norm_coeff=Coord_type(1),true); - } + { + *nn_out++= std::make_pair(c->vertex(li), Coord_type(1)); + return make_triple(nn_out, norm_coeff = Coord_type(1),true); + } else if (dt.is_infinite(c)) - return make_triple(nn_out, Coord_type(1), false);//point outside the convex-hull + { + //point outside the convex-hull + return make_triple(nn_out, Coord_type(1), false); + } std::set cells; - // To replace the forbidden access to the "in conflict" flag : + // To replace the forbidden access to the "in conflict" flag : // std::find operations on this set std::vector bound_facets; bound_facets.reserve(32); typename std::vector::iterator bound_it; // Find the cells in conflict with Q - dt.find_conflicts(Q, c, - std::back_inserter(bound_facets), - std::inserter(cells,cells.begin())); + dt.find_conflicts(Q, c, + std::back_inserter(bound_facets), + std::inserter(cells,cells.begin())); std::map coordinate; typename std::map::iterator coor_it; @@ -165,12 +171,13 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt, template Triple< OutputIterator, // iterator with value type std::pair - typename Dt::Geom_traits::FT, // Should provide 0 and 1 - bool > + typename Dt::Geom_traits::FT, // Should provide 0 and 1 + bool > sibson_natural_neighbor_coordinates_3(const Dt& dt, - const typename Dt::Geom_traits::Point_3& Q, - OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff, - const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) + const typename Dt::Geom_traits::Point_3& Q, + OutputIterator nn_out, + typename Dt::Geom_traits::FT& norm_coeff, + const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) { typedef typename Dt::Geom_traits Gt; typedef typename Gt::Point_3 Point; @@ -182,27 +189,30 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt, CGAL_triangulation_precondition (dt.dimension()== 3); - Locate_type lt; int li, lj; - + Locate_type lt; + int li, lj; Cell_handle c = dt.locate( Q, lt, li, lj, start); if ( lt == Dt::VERTEX ) - { - *nn_out++= std::make_pair(c->vertex(li),Coord_type(1)); - return make_triple(nn_out,norm_coeff=Coord_type(1),true); - } + { + *nn_out++ = std::make_pair(c->vertex(li),Coord_type(1)); + return make_triple(nn_out,norm_coeff=Coord_type(1),true); + } else if (dt.is_infinite(c)) - return make_triple(nn_out, Coord_type(1), false);//point outside the convex-hull + { + //point outside the convex-hull + return make_triple(nn_out, Coord_type(1), false); + } - std::set cells; + std::set cells; typename std::set::iterator cit; - // To replace the forbidden access to the "in conflict" flag : + // To replace the forbidden access to the "in conflict" flag : // std::find operations on this set // Find the cells in conflict with Q - dt.find_conflicts(Q, c, - Emptyset_iterator(), - std::inserter(cells,cells.begin())); + dt.find_conflicts(Q, c, + Emptyset_iterator(), + std::inserter(cells,cells.begin())); std::map coordinate; typename std::map::iterator coor_it; @@ -296,34 +306,34 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt, } return make_triple(nn_out,norm_coeff,true); } - -template + +template bool is_correct_natural_neighborhood(const Dt& /*dt*/, - const typename Dt::Geom_traits::Point_3 & Q, - InputIterator it_begin, InputIterator it_end, - const typename Dt::Geom_traits::FT & norm_coeff) -{ + const typename Dt::Geom_traits::Point_3& Q, + InputIterator it_begin, InputIterator it_end, + const typename Dt::Geom_traits::FT& norm_coeff) +{ typedef typename Dt::Geom_traits Gt; typedef typename Gt::FT Coord_type; Coord_type sum_x(0); Coord_type sum_y(0); Coord_type sum_z(0); - InputIterator it; - for(it = it_begin ; it != it_end ; ++it) - { - sum_x += it->second*(it->first->point().x()); - sum_y += it->second*(it->first->point().y()); - sum_z += it->second*(it->first->point().z()); - } + InputIterator it; + for(it = it_begin ; it != it_end ; ++it) + { + sum_x += it->second*(it->first->point().x()); + sum_y += it->second*(it->first->point().y()); + sum_z += it->second*(it->first->point().z()); + } //!!!! to be replaced by a linear combination of points as soon // as it is available in the kernel. - std::cout << sum_x/norm_coeff << " " - << sum_y/norm_coeff << " " - << sum_z/norm_coeff << std::endl; - return ((sum_x==norm_coeff*Q.x())&&(sum_y==norm_coeff*Q.y()) - &&(sum_z==norm_coeff*Q.z())); + std::cout << sum_x/norm_coeff << " " + << sum_y/norm_coeff << " " + << sum_z/norm_coeff << std::endl; + return ((sum_x == norm_coeff*Q.x()) && (sum_y == norm_coeff*Q.y()) + && (sum_z == norm_coeff*Q.z())); } - + // ====================== Geometric Traits utilities ========================================= // === Definitions diff --git a/Interpolation/include/CGAL/predicates/predicates_for_voronoi_intersection_cartesian_2_3.h b/Interpolation/include/CGAL/predicates/predicates_for_voronoi_intersection_cartesian_2_3.h index e74bff77181..ce39bae2fed 100644 --- a/Interpolation/include/CGAL/predicates/predicates_for_voronoi_intersection_cartesian_2_3.h +++ b/Interpolation/include/CGAL/predicates/predicates_for_voronoi_intersection_cartesian_2_3.h @@ -31,11 +31,11 @@ namespace CGAL { template < class RT> Oriented_side side_of_plane_centered_sphere_translateC3( - const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &qx, const RT &qy, const RT &qz, - const RT &rx, const RT &ry, const RT &rz, - const RT &tx, const RT &ty, const RT &tz) + const RT &ax, const RT &ay, const RT &az, + const RT &nx, const RT &ny, const RT &nz, + const RT &qx, const RT &qy, const RT &qz, + const RT &rx, const RT &ry, const RT &rz, + const RT &tx, const RT &ty, const RT &tz) { RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) + CGAL_NTS square(qz); RT r2 = CGAL_NTS square(rx) + CGAL_NTS square(ry) + CGAL_NTS square(rz); @@ -44,14 +44,14 @@ side_of_plane_centered_sphere_translateC3( na *= RT(2.0); Sign num = sign_of_determinant(rx, ry, rz, r2, - qx, qy, qz, q2, - nx, ny, nz, na, - tx, ty, tz, t2); + qx, qy, qz, q2, + nx, ny, nz, na, + tx, ty, tz, t2); //denumerator: Sign den = sign_of_determinant(nx,ny,nz, - qx,qy,qz, - rx,ry,rz); + qx,qy,qz, + rx,ry,rz); CGAL_assertion(den != ZERO); return den * num; @@ -60,11 +60,11 @@ side_of_plane_centered_sphere_translateC3( template < class RT> Oriented_side side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &px, const RT &py, const RT &pz, - const RT &qx, const RT &qy, const RT &qz, - const RT &rx, const RT &ry, const RT &rz, - const RT &tx, const RT &ty, const RT &tz) + const RT &nx, const RT &ny, const RT &nz, + const RT &px, const RT &py, const RT &pz, + const RT &qx, const RT &qy, const RT &qz, + const RT &rx, const RT &ry, const RT &rz, + const RT &tx, const RT &ty, const RT &tz) { // resolution of the system (where c denotes the sphere's center) // @@ -78,19 +78,19 @@ side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, // - seperate computation of det and norm of the expression return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz, - nx, ny, nz, - qx-px, qy-py,qz-pz, - rx-px, ry-py,rz-pz, - tx-px, ty-py,tz-pz); + nx, ny, nz, + qx-px, qy-py,qz-pz, + rx-px, ry-py,rz-pz, + tx-px, ty-py,tz-pz); } template < class RT> Oriented_side side_of_plane_centered_sphere_translateC3( - const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &qx, const RT &qy, const RT &qz, - const RT &rx, const RT &ry, const RT &rz) + const RT &ax, const RT &ay, const RT &az, + const RT &nx, const RT &ny, const RT &nz, + const RT &qx, const RT &qy, const RT &qz, + const RT &rx, const RT &ry, const RT &rz) { //first choice of n_ortho: (ny+nz, -nx, -nx) // if it is @@ -100,24 +100,24 @@ side_of_plane_centered_sphere_translateC3( na *= RT(2.0); Sign num = sign_of_determinant(qx, qy, qz, q2, - ny, -nx, RT(0), RT(0), - nx, ny, nz, na, - rx, ry, rz, r2); + ny, -nx, RT(0), RT(0), + nx, ny, nz, na, + rx, ry, rz, r2); //denumerator: Sign den = sign_of_determinant(nx,ny,nz, - ny,-nx, RT(0), - qx,qy,qz); + ny,-nx, RT(0), + qx,qy,qz); if (den==ZERO) { // bad choice: (ny,-nx,0) is coplanar with n,q. // by precondition: q and n may not be collinear // => the cross product q*n is orthogonal to q, n and not coplanar num = sign_of_determinant(qx, qy, qz, q2, - ny*qz-nz*qy, nz*qx-nx*qz,nx*qy-ny*qx, RT(0), - nx, ny, nz, na, - rx, ry, rz, r2); + ny*qz-nz*qy, nz*qx-nx*qz,nx*qy-ny*qx, RT(0), + nx, ny, nz, na, + rx, ry, rz, r2); den = sign_of_determinant(nx,ny,nz, - ny*qz-nz*qy, nz*qx - nx*qz,nx*qy-ny*qx, - qx,qy,qz); + ny*qz-nz*qy, nz*qx - nx*qz,nx*qy-ny*qx, + qx,qy,qz); } CGAL_assertion(den != ZERO); return den * num; @@ -126,10 +126,10 @@ side_of_plane_centered_sphere_translateC3( template < class RT> Oriented_side side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, - const RT &nx, const RT &ny, const RT &nz, - const RT &px, const RT &py, const RT &pz, - const RT &qx, const RT &qy, const RT &qz, - const RT &rx, const RT &ry, const RT &rz) + const RT &nx, const RT &ny, const RT &nz, + const RT &px, const RT &py, const RT &pz, + const RT &qx, const RT &qy, const RT &qz, + const RT &rx, const RT &ry, const RT &rz) { // precondition: no two points p,q,r have the same projection // <=> (p-q),(p-r), (q-r) may not be collinear to n @@ -147,9 +147,9 @@ side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, // - seperate computation of det and nom of the expression return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz, - nx, ny, nz, - qx-px, qy-py,qz-pz, - rx-px, ry-py,rz-pz); + nx, ny, nz, + qx-px, qy-py,qz-pz, + rx-px, ry-py,rz-pz); } } //namespace CGAL diff --git a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h index e20b463172e..1930d8f163d 100644 --- a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h +++ b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #ifndef CGAL_REGULAR_NEIGHBOR_COORDINATES_2_H @@ -23,23 +19,22 @@ #include - -#include -#include #include //for definition of class Project_vertex_output_iterator #include +#include +#include +#include + namespace CGAL { -// in this functions, the traits class is defined via the regular -// triangulation -// see natural_neighbor_coordinates_2 for a proposal for signatures -// that allow to pass the traits class as argument +// In these functions, the traits class is defined via the regular triangulation. +// See natural_neighbor_coordinates_2 for a proposal for signatures +// that allow to pass the traits class as argument. - -//the following two functions suppose that +// The following two functions assume that // OutputIterator has value type // std::pair //!!!they are not documented!!! @@ -49,41 +44,40 @@ namespace CGAL { template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_vertex_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out) + const typename Rt::Weighted_point& p, + OutputIterator out) { return regular_neighbor_coordinates_vertex_2(rt, p, out, - typename Rt::Face_handle()); + typename Rt::Face_handle()); } -//Face_handle start is known: +// Face_handle start is known: // OutputIterator has value type // std::pair template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_vertex_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - typename Rt::Face_handle start) + const typename Rt::Weighted_point& p, + OutputIterator out, + typename Rt::Face_handle start) { return regular_neighbor_coordinates_vertex_2(rt, p, out, - Emptyset_iterator(), start); + Emptyset_iterator(), start); } -//the Voronoi vertices of the power cell are known: +// The Voronoi vertices of the power cell are known: // OutputIterator has value type // std::pair template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_vertex_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - OutputIteratorVorVertices vor_vertices, - typename Rt::Face_handle start) + const typename Rt::Weighted_point& p, + OutputIterator out, + OutputIteratorVorVertices vor_vertices, + typename Rt::Face_handle start) { - //out: the result of the coordinate computation - //vor_vertices: the vertices of the power cell (to avoid - // recomputation) + // out: the result of the coordinate computation + // vor_vertices: the vertices of the power cell (to avoid recomputation) typedef typename Rt::Geom_traits Traits; typedef typename Traits::FT Coord_type; @@ -93,24 +87,21 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt, typedef typename Rt::Locate_type Locate_type; CGAL_precondition(rt.dimension() == 2); - + Locate_type lt; int li; Face_handle fh = rt.locate(p, lt, li, start); - //the point must lie inside the convex hull - // sinon return false: - if(lt == Rt::OUTSIDE_AFFINE_HULL || lt == - Rt::OUTSIDE_CONVEX_HULL - || (lt == Rt::EDGE && (rt.is_infinite(fh) - || rt.is_infinite(fh->neighbor(li))))) + // the point must lie inside the convex hull otherwisereturn false: + if(lt == Rt::OUTSIDE_AFFINE_HULL || lt == Rt::OUTSIDE_CONVEX_HULL + || (lt == Rt::EDGE + && (rt.is_infinite(fh) || rt.is_infinite(fh->neighbor(li))))) return make_triple(out, Coord_type(1), false); if (lt == Rt::VERTEX) { //the point must be in conflict: - CGAL_precondition(rt.power_test(fh->vertex(li)->point(), p) != - ON_NEGATIVE_SIDE); + CGAL_precondition(rt.power_test(fh->vertex(li)->point(), p) != ON_NEGATIVE_SIDE); if (rt.power_test(fh->vertex(li)->point(), p) ==ON_ORIENTED_BOUNDARY) { *out++= std::make_pair(fh->vertex(li),Coord_type(1)); @@ -119,52 +110,50 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt, } std::list hole; - std::list< Vertex_handle > hidden_vertices; + std::list hidden_vertices; rt.get_boundary_of_conflicts_and_hidden_vertices(p, - std::back_inserter(hole), - std::back_inserter - (hidden_vertices), - fh); - return regular_neighbor_coordinates_vertex_2 - (rt, p, out, vor_vertices, hole.begin(),hole.end(), - hidden_vertices.begin(), hidden_vertices.end()); + std::back_inserter(hole), + std::back_inserter + (hidden_vertices), + fh); + return regular_neighbor_coordinates_vertex_2(rt, p, out, vor_vertices, + hole.begin(),hole.end(), + hidden_vertices.begin(), + hidden_vertices.end()); } // OutputIterator has value type // std::pair template + class VertexIterator > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_vertex_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, EdgeIterator - hole_begin, EdgeIterator hole_end, - VertexIterator hidden_vertices_begin, - VertexIterator hidden_vertices_end) + const typename Rt::Weighted_point& p, + OutputIterator out, + EdgeIterator hole_begin, EdgeIterator hole_end, + VertexIterator hidden_vertices_begin, + VertexIterator hidden_vertices_end) { - return regular_neighbor_coordinates_vertex_2(rt, p, - out,Emptyset_iterator(), - hole_begin, hole_end, - hidden_vertices_begin, - hidden_vertices_end); + return regular_neighbor_coordinates_vertex_2(rt, p, out, Emptyset_iterator(), + hole_begin, hole_end, + hidden_vertices_begin, + hidden_vertices_end); } - // OutputIterator has value type // std::pair template + class VertexIterator , class OutputIteratorVorVertices > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_vertex_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - OutputIteratorVorVertices vor_vertices, - EdgeIterator - hole_begin, EdgeIterator hole_end, - VertexIterator hidden_vertices_begin, - VertexIterator hidden_vertices_end) + const typename Rt::Weighted_point& p, + OutputIterator out, + OutputIteratorVorVertices vor_vertices, + EdgeIterator hole_begin, EdgeIterator hole_end, + VertexIterator hidden_vertices_begin, + VertexIterator hidden_vertices_end) { //precondition: p must lie inside the non-empty hole // (=^ inside convex hull of neighbors) @@ -204,8 +193,8 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt, //a first Voronoi vertex of the cell of p: vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object() - (current->point(), - hit->first->vertex(rt.ccw(hit->second))->point(), p); + (current->point(), + hit->first->vertex(rt.ccw(hit->second))->point(), p); *vor_vertices++= vor[0]; //triangulation of the Voronoi subcell: @@ -215,17 +204,17 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt, vor[1] = rt.dual(fc); // iteration over all other "old" Voronoi vertices while(!fc->has_vertex(prev)) - { - ++fc; - vor[2] = rt.dual(fc); - - area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits()); - vor[1] = vor[2]; - } + { + ++fc; + vor[2] = rt.dual(fc); + + area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits()); + vor[1] = vor[2]; + } //the second Voronoi vertex of the cell of p: vor[2] = - rt.geom_traits().construct_weighted_circumcenter_2_object() - (prev->point(),current->point(),p); + rt.geom_traits().construct_weighted_circumcenter_2_object() + (prev->point(),current->point(),p); *vor_vertices++= vor[2]; area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits()); @@ -279,11 +268,10 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt, template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out) + const typename Rt::Weighted_point& p, + OutputIterator out) { - return regular_neighbor_coordinates_2(rt, p, out, - typename Rt::Face_handle()); + return regular_neighbor_coordinates_2(rt, p, out, typename Rt::Face_handle()); } //OutputIterator has value type @@ -292,12 +280,11 @@ regular_neighbor_coordinates_2(const Rt& rt, template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - typename Rt::Face_handle start) + const typename Rt::Weighted_point& p, + OutputIterator out, + typename Rt::Face_handle start) { - return regular_neighbor_coordinates_2(rt, p, out, - Emptyset_iterator(), start); + return regular_neighbor_coordinates_2(rt, p, out, Emptyset_iterator(), start); } //OutputIterator has value type @@ -306,10 +293,10 @@ regular_neighbor_coordinates_2(const Rt& rt, template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - OutputIteratorVorVertices vor_vertices, - typename Rt::Face_handle start) + const typename Rt::Weighted_point& p, + OutputIterator out, + OutputIteratorVorVertices vor_vertices, + typename Rt::Face_handle start) { //out: the result of the coordinate computation //vor_vertices: the vertices of the power cell (to avoid @@ -317,47 +304,44 @@ regular_neighbor_coordinates_2(const Rt& rt, Project_vertex_output_iterator op(out); CGAL_precondition(rt.dimension() == 2); - + Triple< Project_vertex_output_iterator, - typename Rt::Geom_traits::FT, bool > result = - regular_neighbor_coordinates_vertex_2 - (rt, p, op , vor_vertices, start); + typename Rt::Geom_traits::FT, bool > result = + regular_neighbor_coordinates_vertex_2(rt, p, op, vor_vertices, start); + return make_triple(result.first.base(), result.second, result.third); } - //OutputIterator has value type // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, EdgeIterator - hole_begin, EdgeIterator hole_end, - VertexIterator hidden_vertices_begin, - VertexIterator hidden_vertices_end) + const typename Rt::Weighted_point& p, + OutputIterator out, EdgeIterator + hole_begin, EdgeIterator hole_end, + VertexIterator hidden_vertices_begin, + VertexIterator hidden_vertices_end) { - return regular_neighbor_coordinates_2(rt, p, - out,Emptyset_iterator(), - hole_begin, hole_end, - hidden_vertices_begin, - hidden_vertices_end); + return regular_neighbor_coordinates_2(rt, p, out, Emptyset_iterator(), + hole_begin, hole_end, + hidden_vertices_begin, + hidden_vertices_end); } - //OutputIterator has value type // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - const typename Rt::Weighted_point& p, - OutputIterator out, - OutputIteratorVorVertices vor_vertices, - EdgeIterator hole_begin, EdgeIterator hole_end, - VertexIterator hidden_vertices_begin, - VertexIterator hidden_vertices_end) + const typename Rt::Weighted_point& p, + OutputIterator out, + OutputIteratorVorVertices vor_vertices, + EdgeIterator hole_begin, EdgeIterator hole_end, + VertexIterator hidden_vertices_begin, + VertexIterator hidden_vertices_end) { //precondition: p must lie inside the non-empty hole // (=^ inside convex hull of neighbors) @@ -367,10 +351,11 @@ regular_neighbor_coordinates_2(const Rt& rt, Project_vertex_output_iterator op(out); Triple< Project_vertex_output_iterator, - typename Rt::Geom_traits::FT, bool > result = - regular_neighbor_coordinates_vertex_2 - (rt, p, op , vor_vertices, hole_begin,hole_end, - hidden_vertices_begin, hidden_vertices_end); + typename Rt::Geom_traits::FT, bool > result = + regular_neighbor_coordinates_vertex_2(rt, p, op , vor_vertices, + hole_begin, hole_end, + hidden_vertices_begin, + hidden_vertices_end); return make_triple(result.first.base(), result.second, result.third); } @@ -382,8 +367,8 @@ regular_neighbor_coordinates_2(const Rt& rt, template Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > regular_neighbor_coordinates_2(const Rt& rt, - typename Rt::Vertex_handle vh, - OutputIterator out) + typename Rt::Vertex_handle vh, + OutputIterator out) { //this functions creates a small triangulation of the // incident vertices of this vertex and computes the @@ -391,10 +376,9 @@ regular_neighbor_coordinates_2(const Rt& rt, typedef typename Rt::Vertex_circulator Vertex_circulator; CGAL_precondition(rt.dimension() == 2); - + Rt t2; - Vertex_circulator vc = rt.incident_vertices(vh), - done(vc); + Vertex_circulator vc = rt.incident_vertices(vh), done(vc); do{ CGAL_assertion(!rt.is_infinite(vc)); t2.insert(vc->point()); @@ -404,7 +388,6 @@ regular_neighbor_coordinates_2(const Rt& rt, return regular_neighbor_coordinates_2(t2, vh->point(), out); } - //class providing a function object: //OutputIterator has value type // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> @@ -414,8 +397,8 @@ class regular_neighbor_coordinates_2_object public: Triple< OutputIterator, typename Rt::Geom_traits::FT , bool > operator()(const Rt& rt, - typename Rt::Vertex_handle vh, - OutputIterator out) const + typename Rt::Vertex_handle vh, + OutputIterator out) const { return regular_neighbor_coordinates_2(rt, vh, out); } diff --git a/Interpolation/include/CGAL/sibson_gradient_fitting.h b/Interpolation/include/CGAL/sibson_gradient_fitting.h index 266e0bd99e8..fbb49f08cf6 100644 --- a/Interpolation/include/CGAL/sibson_gradient_fitting.h +++ b/Interpolation/include/CGAL/sibson_gradient_fitting.h @@ -14,7 +14,7 @@ // // $URL$ // $Id$ -// +// // // Author(s) : Julia Floetotto @@ -23,24 +23,26 @@ #include - -#include #include #include #include +#include +#include + namespace CGAL { template < class ForwardIterator, class Functor, class Traits> typename Traits::Vector_d sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond, - const typename - std::iterator_traits:: - value_type::second_type& - norm, const typename - std::iterator_traits::value_type - ::first_type& p, Functor function_value, - const Traits& traits) + const typename + std::iterator_traits:: + value_type::second_type& norm, + const typename + std::iterator_traits::value_type + ::first_type& p, + Functor function_value, + const Traits& traits) { CGAL_precondition( first!=beyond && norm!=0); typedef typename Traits::Aff_transformation_d Aff_transformation; @@ -50,29 +52,29 @@ sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond, CGAL_assertion(fn.second); //function value of p is valid typename Traits::Vector_d pn = - traits.construct_vector_d_object()(NULL_VECTOR); + traits.construct_vector_d_object()(NULL_VECTOR); Aff_transformation scaling, m, - Hn(traits.construct_null_matrix_d_object()()); + Hn(traits.construct_null_matrix_d_object()()); for(;first!=beyond; ++first){ Coord_type square_dist = traits.compute_squared_distance_d_object() - (first->first, p); + (first->first, p); CGAL_assertion(square_dist != 0); - Coord_type scale = first->second/(norm*square_dist); - typename Traits::Vector_d d= - traits.construct_vector_d_object()(p, first->first); + Coord_type scale = first->second / (norm*square_dist); + typename Traits::Vector_d d = + traits.construct_vector_d_object()(p, first->first); //compute the vector pn: typename Functor::result_type f = function_value(first->first); CGAL_assertion(f.second);//function value of first->first is valid pn = pn + traits.construct_scaled_vector_d_object() - (d,scale * (f.first - fn.first)); + (d,scale * (f.first - fn.first)); //compute the matrix Hn: m = traits.construct_outer_product_d_object()(d); scaling = traits.construct_scaling_matrix_d_object()(scale); - Hn = traits.construct_sum_matrix_d_object()(Hn, scaling * m); + Hn = traits.construct_sum_matrix_d_object()(Hn, scaling * m); } return Hn.inverse().transform(pn); @@ -82,10 +84,10 @@ template < class Triangul, class OutputIterator, class Functor, class CoordFunctor, class Traits> OutputIterator sibson_gradient_fitting(const Triangul& tr, - OutputIterator out, - Functor function_value, - CoordFunctor compute_coordinates, - const Traits& traits) + OutputIterator out, + Functor function_value, + CoordFunctor compute_coordinates, + const Traits& traits) { typedef typename Traits::Point_d Point; typedef typename Traits::FT Coord_type; @@ -93,22 +95,19 @@ sibson_gradient_fitting(const Triangul& tr, std::vector< std::pair< Point, Coord_type > > coords; Coord_type norm; - typename Triangul::Finite_vertices_iterator - vit = tr.finite_vertices_begin(); + typename Triangul::Finite_vertices_iterator vit = tr.finite_vertices_begin(); for(; vit != tr.finite_vertices_end(); ++vit){ - //test if vit is a convex hull vertex: - //otherwise do nothing + //test if vit is a convex hull vertex, otherwise do nothing if (!tr.is_edge(vit, tr.infinite_vertex())) { norm = compute_coordinates(tr, vit, std::back_inserter(coords)).second; *out++ = std::make_pair(vit->point(), - sibson_gradient_fitting(coords.begin(), - coords.end(), - norm, vit->point(), - function_value, - traits)); + sibson_gradient_fitting(coords.begin(), + coords.end(), + norm, vit->point(), + function_value, + traits)); coords.clear(); - } } return out; @@ -122,33 +121,37 @@ sibson_gradient_fitting(const Triangul& tr, template < class Dt, class OutputIterator, class Functor, class Traits> OutputIterator sibson_gradient_fitting_nn_2(const Dt& dt, - OutputIterator out, - Functor function_value, - const Traits& traits) + OutputIterator out, + Functor function_value, + const Traits& traits) { - typedef typename std::back_insert_iterator< std::vector< std::pair< - typename Traits::Point_d,typename Traits::FT > > > CoordInserter; + typedef typename std::back_insert_iterator< + std::vector< + std::pair< typename Traits::Point_d, + typename Traits::FT > > > CoordInserter; - return sibson_gradient_fitting - (dt, out, function_value, - natural_neighbor_coordinates_2_object< Dt, CoordInserter >(), - traits); + return sibson_gradient_fitting(dt, out, function_value, + natural_neighbor_coordinates_2_object< Dt, + CoordInserter >(), + traits); } template < class Rt, class OutputIterator, class Functor, class Traits> OutputIterator sibson_gradient_fitting_rn_2(const Rt& rt, - OutputIterator out, - Functor function_value, - const Traits& traits) + OutputIterator out, + Functor function_value, + const Traits& traits) { - typedef typename std::back_insert_iterator< std::vector< std::pair< - typename Traits::Point_d,typename Traits::FT > > > CoordInserter; + typedef typename std::back_insert_iterator< + std::vector< + std::pair< typename Traits::Point_d, + typename Traits::FT > > > CoordInserter; - return sibson_gradient_fitting - (rt, out, function_value, - regular_neighbor_coordinates_2_object< Rt, CoordInserter >(), - traits); + return sibson_gradient_fitting(rt, out, function_value, + regular_neighbor_coordinates_2_object< Rt, + CoordInserter >(), + traits); } } //namespace CGAL diff --git a/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h b/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h index 64a0296f6e4..f439e0deef7 100644 --- a/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h +++ b/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto // ATTENTION : the surface is supposed to be a closed surface @@ -25,40 +21,43 @@ #include - -#include #include #include #include #include +#include +#include +#include +#include +#include +#include + namespace CGAL { template inline Triple< OutputIterator, typename Kernel::FT, bool > -surface_neighbor_coordinates_3(InputIterator - first, InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - OutputIterator out, - const Kernel&) +surface_neighbor_coordinates_3(InputIterator first, InputIterator beyond, + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + OutputIterator out, + const Kernel&) { typedef Voronoi_intersection_2_traits_3 I_gt; - return surface_neighbor_coordinates_3(first, beyond, p, out, I_gt(p,normal)); + return surface_neighbor_coordinates_3(first, beyond, p, out, I_gt(p, normal)); } template Triple< OutputIterator, typename ITraits::FT, bool > -surface_neighbor_coordinates_3(InputIterator - first, InputIterator beyond, - const typename ITraits::Point_2& p, - OutputIterator out, - const ITraits& traits) +surface_neighbor_coordinates_3(InputIterator first, InputIterator beyond, + const typename ITraits::Point_2& p, + OutputIterator out, + const ITraits& traits) { //definition of the Voronoi intersection triangulation: - typedef Regular_triangulation_2< ITraits> I_triangulation; + typedef Regular_triangulation_2 I_triangulation; //build Voronoi intersection triangulation: I_triangulation it(traits); @@ -99,15 +98,15 @@ surface_neighbor_coordinates_3(InputIterator template Quadruple< OutputIterator, typename Kernel::FT, bool, bool > surface_neighbor_coordinates_certified_3(InputIterator - first, InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - OutputIterator out, - const Kernel& ) + first, InputIterator beyond, + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + OutputIterator out, + const Kernel& ) { typedef Voronoi_intersection_2_traits_3 I_gt; - return surface_neighbor_coordinates_certified_3 - (first, beyond, p, out, I_gt(p,normal)); + return surface_neighbor_coordinates_certified_3(first, beyond, p, out, + I_gt(p,normal)); } //this function takes the radius of the sphere centered on p @@ -116,34 +115,33 @@ surface_neighbor_coordinates_certified_3(InputIterator template inline Quadruple< OutputIterator, typename Kernel::FT, bool, bool > -surface_neighbor_coordinates_certified_3( - InputIterator first, InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - const typename Kernel::FT& radius, - OutputIterator out, const Kernel& ) +surface_neighbor_coordinates_certified_3(InputIterator first, InputIterator beyond, + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + const typename Kernel::FT& radius, + OutputIterator out, const Kernel& ) { typedef Voronoi_intersection_2_traits_3 I_gt; - return surface_neighbor_coordinates_certified_3 - (first, beyond, p, radius, out, I_gt(p,normal)); + return surface_neighbor_coordinates_certified_3(first, beyond, p, radius, out, + I_gt(p,normal)); } // FIXME : this should probably be replaced by some kernel functor. //struct necessary to sort the points by distance to p: //also used in surface_neighbors_3.h -template +template struct closer_to_point : public std::less { typedef typename Traits::Point_2 Point_2; closer_to_point(const Point_2& _p, const Traits& t) - : p(_p), traits(t) {} + : p(_p), traits(t) { } bool operator()(const Point_2& q, const Point_2& r) const - { - return traits.less_distance_to_point_2_object()(p,q,r); - } + { + return traits.less_distance_to_point_2_object()(p,q,r); + } private: Point_2 p; @@ -153,34 +151,31 @@ private: // Versions with instantiated traits class: template Quadruple< OutputIterator, typename ITraits::FT, bool, bool > -surface_neighbor_coordinates_certified_3(InputIterator - first, InputIterator beyond, - const typename ITraits::Point_2& p, - OutputIterator out, - const ITraits& traits) +surface_neighbor_coordinates_certified_3(InputIterator first, + InputIterator beyond, + const typename ITraits::Point_2& p, + OutputIterator out, + const ITraits& traits) { //find the point in [first,beyond) furthest from p: InputIterator furthest = std::max_element(first, beyond, - closer_to_point(p, traits)); + closer_to_point(p, traits)); - return surface_neighbor_coordinates_certified_3 - (first, beyond, p, - traits.compute_squared_distance_2_object()(p,*furthest), - out, traits); + return surface_neighbor_coordinates_certified_3(first, beyond, p, + traits.compute_squared_distance_2_object()(p,*furthest), + out, traits); } //with radius(maximal distance from p to [first,beyond)) as // add. parameter: template Quadruple< OutputIterator, typename ITraits::FT, bool, bool > -surface_neighbor_coordinates_certified_3(InputIterator - first, InputIterator beyond, - const typename - ITraits::Point_2& p, - const typename ITraits::FT& - radius, - OutputIterator out, - const ITraits& traits) +surface_neighbor_coordinates_certified_3(InputIterator first, + InputIterator beyond, + const typename ITraits::Point_2& p, + const typename ITraits::FT& radius, + OutputIterator out, + const ITraits& traits) { //definition of the Voronoi intersection triangulation: typedef Regular_triangulation_2< ITraits> I_triangulation; @@ -199,7 +194,7 @@ surface_neighbor_coordinates_certified_3(InputIterator //collect the Voronoi vertices of the cell of p in order to //determine the furthest distance from p to the boundary of its cell - std::vector< typename ITraits::Point_2 > vor_vertices; + std::vector< typename ITraits::Point_2 > vor_vertices; //unfortunately, there is no function call without Face_handle // "start" because this would cause type conflicts because @@ -262,11 +257,10 @@ template inline Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > surface_neighbor_coordinates_3(const Dt& dt, - const typename Dt::Geom_traits::Point_3& p, - const typename Dt::Geom_traits::Vector_3& normal, - OutputIterator out, - typename Dt::Cell_handle start - = typename Dt::Cell_handle()) + const typename Dt::Geom_traits::Point_3& p, + const typename Dt::Geom_traits::Vector_3& normal, + OutputIterator out, + typename Dt::Cell_handle start = typename Dt::Cell_handle()) { typedef Voronoi_intersection_2_traits_3 I_gt; return surface_neighbor_coordinates_3(dt, p, out, I_gt(p,normal), start); @@ -275,27 +269,25 @@ surface_neighbor_coordinates_3(const Dt& dt, template Triple< OutputIterator, typename ITraits::FT, bool > surface_neighbor_coordinates_3(const Dt& dt, - const typename ITraits::Point_2& p, - OutputIterator out, const ITraits& traits, - typename Dt::Cell_handle start - = typename Dt::Cell_handle()) + const typename ITraits::Point_2& p, + OutputIterator out, const ITraits& traits, + typename Dt::Cell_handle start = typename Dt::Cell_handle()) { typedef typename ITraits::FT Coord_type; typedef typename ITraits::Point_2 Point_3; - typedef typename Dt::Cell_handle Cell_handle; + typedef typename Dt::Cell_handle Cell_handle; typedef typename Dt::Vertex_handle Vertex_handle; - typedef typename Dt::Locate_type Locate_type; + typedef typename Dt::Locate_type Locate_type; //the Vertex_handle is, in fact, an iterator over vertex: typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point; - typedef Iterator_project< - typename std::list< Vertex_handle >::iterator, - Proj_point, - const Point_3&, - const Point_3*, - std::ptrdiff_t, - std::forward_iterator_tag> Point_iterator; + typedef Iterator_project::iterator, + Proj_point, + const Point_3&, + const Point_3*, + std::ptrdiff_t, + std::forward_iterator_tag> Point_iterator; Locate_type lt; int li, lj ; @@ -303,15 +295,13 @@ surface_neighbor_coordinates_3(const Dt& dt, //if p is located on a vertex: the only neighbor is found if(lt == Dt::VERTEX){ - *out++= std::make_pair(c->vertex(li)->point(), - Coord_type(1)); + *out++= std::make_pair(c->vertex(li)->point(), Coord_type(1)); return make_triple(out, Coord_type(1), true); } //the candidate points are the points of dt in conflict with p: typename std::list< Vertex_handle > conflict_vertices; - dt.vertices_on_conflict_zone_boundary(p,c, - std::back_inserter(conflict_vertices)); + dt.vertices_on_conflict_zone_boundary(p, c, std::back_inserter(conflict_vertices)); for (typename std::list< Vertex_handle >::iterator it = conflict_vertices.begin(); it != conflict_vertices.end();){ @@ -323,10 +313,9 @@ surface_neighbor_coordinates_3(const Dt& dt, it++; } } - return surface_neighbor_coordinates_3 - (Point_iterator(conflict_vertices.begin()), - Point_iterator(conflict_vertices.end()), - p, out, traits); + return surface_neighbor_coordinates_3(Point_iterator(conflict_vertices.begin()), + Point_iterator(conflict_vertices.end()), + p, out, traits); } } //namespace CGAL diff --git a/Interpolation/include/CGAL/surface_neighbors_3.h b/Interpolation/include/CGAL/surface_neighbors_3.h index cf748ee25c5..b6a2814834b 100644 --- a/Interpolation/include/CGAL/surface_neighbors_3.h +++ b/Interpolation/include/CGAL/surface_neighbors_3.h @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #ifndef CGAL_SURFACE_NEIGHBORS_3_H @@ -23,8 +19,6 @@ #include - -#include #include #include #include @@ -33,6 +27,10 @@ // the function object Project_vertex_iterator_to_point #include +#include +#include +#include + namespace CGAL { //without Delaunay filtering @@ -40,9 +38,9 @@ template inline OutputIterator surface_neighbors_3(InputIterator first, InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - OutputIterator out, const Kernel& ) + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + OutputIterator out, const Kernel&) { typedef Voronoi_intersection_2_traits_3 I_gt; return surface_neighbors_3(first, beyond, p, out, I_gt(p,normal)); @@ -51,8 +49,8 @@ surface_neighbors_3(InputIterator first, InputIterator beyond, template OutputIterator surface_neighbors_3(InputIterator first, InputIterator beyond, - const typename ITraits::Point_2& p, - OutputIterator out, const ITraits& traits) + const typename ITraits::Point_2& p, + OutputIterator out, const ITraits& traits) { //definition of the Voronoi intersection triangulation: typedef Regular_triangulation_2< ITraits> I_triangulation; @@ -79,20 +77,19 @@ surface_neighbors_3(InputIterator first, InputIterator beyond, Face_handle fh = it.locate(wp, lt, li); if(lt == I_triangulation::VERTEX){ - *out++ =p; + *out++ = p; return out; } Vertex_handle vh = it.insert(wp, fh); - typename I_triangulation::Vertex_circulator - vc(it.incident_vertices(vh)), - done(vc); + typename I_triangulation::Vertex_circulator vc(it.incident_vertices(vh)), + done(vc); do{ *out++= wp2p(vc->point()); CGAL_assertion(! it.is_infinite(vc)); } - while(vc++!=done); + while(vc++ != done); return out; } @@ -105,10 +102,10 @@ surface_neighbors_3(InputIterator first, InputIterator beyond, template std::pair< OutputIterator, bool > surface_neighbors_certified_3(InputIterator first, - InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - OutputIterator out, const Kernel& ) + InputIterator beyond, + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + OutputIterator out, const Kernel&) { typedef Voronoi_intersection_2_traits_3 I_gt; return surface_neighbors_certified_3(first, beyond, p, out, I_gt(p,normal)); @@ -120,45 +117,45 @@ surface_neighbors_certified_3(InputIterator first, template std::pair< OutputIterator, bool > surface_neighbors_certified_3(InputIterator first, - InputIterator beyond, - const typename Kernel::Point_3& p, - const typename Kernel::Vector_3& normal, - const typename Kernel::FT& radius, - OutputIterator out, - const Kernel& /*K*/) + InputIterator beyond, + const typename Kernel::Point_3& p, + const typename Kernel::Vector_3& normal, + const typename Kernel::FT& radius, + OutputIterator out, + const Kernel& /*K*/) { typedef Voronoi_intersection_2_traits_3 I_gt; return surface_neighbors_certified_3(first, beyond, p, radius, - out, I_gt(p,normal)); + out, I_gt(p,normal)); } // Versions with instantiated traits class: template std::pair< OutputIterator, bool > surface_neighbors_certified_3(InputIterator first, - InputIterator beyond, - const typename ITraits::Point_2& p, - OutputIterator out, - const ITraits& traits) + InputIterator beyond, + const typename ITraits::Point_2& p, + OutputIterator out, + const ITraits& traits) { //find the point in [first,beyond) furthest from p: InputIterator furthest = std::max_element(first, beyond, - closer_to_point(p, traits)); + closer_to_point(p, traits)); return surface_neighbors_certified_3 - (first, beyond, p, - traits.compute_squared_distance_2_object()(p,*furthest), - out, traits); + (first, beyond, p, + traits.compute_squared_distance_2_object()(p,*furthest), + out, traits); } template std::pair< OutputIterator, bool > surface_neighbors_certified_3(InputIterator first, - InputIterator beyond, - const typename ITraits::Point_2& p, - const typename ITraits::FT& radius, - OutputIterator out, - const ITraits& traits) + InputIterator beyond, + const typename ITraits::Point_2& p, + const typename ITraits::FT& radius, + OutputIterator out, + const ITraits& traits) { //definition of the Voronoi intersection triangulation: typedef Regular_triangulation_2< ITraits> I_triangulation; @@ -188,7 +185,7 @@ surface_neighbors_certified_3(InputIterator first, Face_handle fh = it.locate(wp, lt, li); if(lt == I_triangulation::VERTEX){ - *out++ =p; + *out++ = p; return std::make_pair(out,true); } Vertex_handle vh = it.insert(wp, fh); @@ -197,11 +194,11 @@ surface_neighbors_certified_3(InputIterator first, //determine the furthest distance from p to a vertex of its cell bool valid(false); Face_circulator fc(it.incident_faces(vh)), fdone(fc); - do + do{ valid = (!it.is_infinite(fc) && - (4*radius > traits.compute_squared_distance_2_object() - (p, it.dual(fc)))); - while(!valid && ++fc!=fdone); + (4*radius > traits.compute_squared_distance_2_object() + (p, it.dual(fc)))); + }while(!valid && ++fc!=fdone); //get the neighbor points: Vertex_circulator vc(it.incident_vertices(vh)), vdone(vc); @@ -212,17 +209,16 @@ surface_neighbors_certified_3(InputIterator first, return std::make_pair(out, valid); } - //using Delaunay triangulation for candidate point filtering: // => no certification is necessary template inline OutputIterator surface_neighbors_3(const Dt& dt, - const typename Dt::Geom_traits::Point_3& p, - const typename Dt::Geom_traits::Vector_3& normal, - OutputIterator out, - typename Dt::Cell_handle start =typename Dt::Cell_handle()) + const typename Dt::Geom_traits::Point_3& p, + const typename Dt::Geom_traits::Vector_3& normal, + OutputIterator out, + typename Dt::Cell_handle start =typename Dt::Cell_handle()) { typedef Voronoi_intersection_2_traits_3 I_gt; return surface_neighbors_3(dt, p, out, I_gt(p,normal),start); @@ -231,26 +227,24 @@ surface_neighbors_3(const Dt& dt, template OutputIterator surface_neighbors_3(const Dt& dt, - const typename ITraits::Point_2& p, - OutputIterator out, const ITraits& traits, - typename Dt::Cell_handle start - = typename Dt::Cell_handle()) + const typename ITraits::Point_2& p, + OutputIterator out, const ITraits& traits, + typename Dt::Cell_handle start = typename Dt::Cell_handle()) { typedef typename ITraits::Point_2 Point_3; typedef typename Dt::Cell_handle Cell_handle; typedef typename Dt::Vertex_handle Vertex_handle; typedef typename Dt::Locate_type Locate_type; - + //the Vertex_handle is, in fact, an iterator over vertex: - typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point; - typedef Iterator_project< - typename std::list< Vertex_handle >::iterator, - Proj_point, - const Point_3&, - const Point_3*, - std::ptrdiff_t, - std::forward_iterator_tag> Point_iterator; + typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point; + typedef Iterator_project::iterator, + Proj_point, + const Point_3&, + const Point_3*, + std::ptrdiff_t, + std::forward_iterator_tag> Point_iterator; Locate_type lt; int li, lj ; @@ -265,7 +259,7 @@ surface_neighbors_3(const Dt& dt, //otherwise get vertices in conflict typename std::list< Vertex_handle > conflict_vertices; dt.vertices_on_conflict_zone_boundary(p,c, - std::back_inserter(conflict_vertices)); + std::back_inserter(conflict_vertices)); for (typename std::list< Vertex_handle >::iterator it = conflict_vertices.begin(); it != conflict_vertices.end();){ @@ -278,8 +272,8 @@ surface_neighbors_3(const Dt& dt, } } return surface_neighbors_3(Point_iterator(conflict_vertices.begin()), - Point_iterator(conflict_vertices.end()), - p, out, traits); + Point_iterator(conflict_vertices.end()), + p, out, traits); } } //namespace CGAL diff --git a/Interpolation/test/Interpolation/include/CGAL/_test_surface_neighbors_3.cpp b/Interpolation/test/Interpolation/include/CGAL/_test_surface_neighbors_3.cpp index 442bd35b35a..cfea253c43d 100644 --- a/Interpolation/test/Interpolation/include/CGAL/_test_surface_neighbors_3.cpp +++ b/Interpolation/test/Interpolation/include/CGAL/_test_surface_neighbors_3.cpp @@ -106,10 +106,10 @@ bool compare_neighbors(ForwardIteratorCoord first_coord, template < class Triangul, class OutputIterator> OutputIterator test_neighbors(const Triangul& T, const typename - Triangul::Geom_traits::Point_3& p, - const typename Triangul::Geom_traits::Vector_3& n, - const int& version, - OutputIterator out) + Triangul::Geom_traits::Point_3& p, + const typename Triangul::Geom_traits::Vector_3& n, + const int& version, + OutputIterator out) { typedef CGAL::Voronoi_intersection_2_traits_3 I_traits; diff --git a/Interpolation/test/Interpolation/test_interpolation_functions_2.cpp b/Interpolation/test/Interpolation/test_interpolation_functions_2.cpp index c24571207da..37d7b8fb359 100644 --- a/Interpolation/test/Interpolation/test_interpolation_functions_2.cpp +++ b/Interpolation/test/Interpolation/test_interpolation_functions_2.cpp @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #include @@ -24,6 +20,8 @@ #include +#include + typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Delaunay_triangulation_2 Dt; @@ -33,16 +31,16 @@ typedef CGAL::Delaunay_triangulation_2 Dt2; int main() { std::cout << "Testing interpolation functions with 2D NN neighbors " - << std::endl; + << std::endl; std::cout << " using Exact_predicates_exact_constructions_kernel: " - << std::endl ; - _test_interpolation_functions_2_delaunay( Dt(), K::FT(1e-10)); - + << std::endl ; + _test_interpolation_functions_2_delaunay(Dt(), K::FT(1e-10)); + std::cout << "Testing interpolation functions with 2D NN neighbors " - << std::endl; + << std::endl; std::cout << " using Exact_predicates_inexact_constructions_kernel: " - << std::endl ; - _test_interpolation_functions_2_delaunay( Dt2(), K2::FT(1e-10)); + << std::endl ; + _test_interpolation_functions_2_delaunay(Dt2(), K2::FT(1e-10)); std::cout << "test_interpolation_functions_2 is finished" << std::endl; return 0; diff --git a/Interpolation/test/Interpolation/test_natural_neighbors_2.cpp b/Interpolation/test/Interpolation/test_natural_neighbors_2.cpp index 9b341cc99fa..768c5340541 100644 --- a/Interpolation/test/Interpolation/test_natural_neighbors_2.cpp +++ b/Interpolation/test/Interpolation/test_natural_neighbors_2.cpp @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Naceur MESKINI. #include @@ -23,11 +19,11 @@ #include +#include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Delaunay_triangulation_2 Dt; - int main() { std::cout << "Testing NN_neighbors_2 " << std::endl; diff --git a/Interpolation/test/Interpolation/test_regular_neighbors_2.cpp b/Interpolation/test/Interpolation/test_regular_neighbors_2.cpp index f82a9268c86..b4d8db5d539 100644 --- a/Interpolation/test/Interpolation/test_regular_neighbors_2.cpp +++ b/Interpolation/test/Interpolation/test_regular_neighbors_2.cpp @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #include @@ -26,6 +22,7 @@ #include +#include typedef CGAL::Exact_predicates_exact_constructions_kernel K; diff --git a/Interpolation/test/Interpolation/test_surface_neighbor_3.cpp b/Interpolation/test/Interpolation/test_surface_neighbor_3.cpp index c249a1481d7..efc3e18244a 100644 --- a/Interpolation/test/Interpolation/test_surface_neighbor_3.cpp +++ b/Interpolation/test/Interpolation/test_surface_neighbor_3.cpp @@ -12,10 +12,6 @@ // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // -// $URL$ -// $Id$ -// -// // Author(s) : Julia Floetotto #include @@ -25,6 +21,7 @@ #include +#include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Delaunay_triangulation_3 Dt; @@ -35,7 +32,6 @@ typedef CGAL::Delaunay_triangulation_3 Dt2; // Fast_location with exact pred exact const. kernel: typedef CGAL::Delaunay_triangulation_3 Dh; - // Aff_transformation: typedef CGAL::Aff_transformation_3 Transformation; @@ -46,17 +42,17 @@ int main() _test_surface_neighbors_3_sphere( Dt() ); std::cout << " done." << std::endl << std::endl; - std::cout << "Using Exact_predicates_exact_constructions_kernel: " - << std::endl; + << std::endl; //AXIS ALIGNED CUBE Transformation identity(CGAL::IDENTITY); std::cout << "Testing surface_neighbors_3 on a cube with axis : " - << identity(K2::Vector_3(0,0,1)) << ", " - << identity(K2::Vector_3(0,1,0))<< ", " - << identity(K2::Vector_3(1,0,0)) - << std::endl; + << identity(K2::Vector_3(0,0,1)) << ", " + << identity(K2::Vector_3(0,1,0))<< ", " + << identity(K2::Vector_3(1,0,0)) + << std::endl; + std::cout << " with grid sample points"; _test_surface_neighbors_3_cube(Dt2(),identity, 75, K2::FT(1e-29)); std::cout << " done." << std::endl; @@ -67,22 +63,23 @@ int main() //ROTATED CUBE Transformation rotate(K2::RT(1),K2::RT(0),K2::RT(0),K2::RT(0), - K2::RT(0),K2::RT(0.9063),K2::RT(-0.42261826),K2::RT(0), - K2::RT(0),K2::RT(0.42261826),K2::RT(0.9063),K2::RT(0)); + K2::RT(0),K2::RT(0.9063),K2::RT(-0.42261826),K2::RT(0), + K2::RT(0),K2::RT(0.42261826),K2::RT(0.9063),K2::RT(0)); + std::cout << "Testing surface_neighbors_3 on a ROTATED cube "<< std::endl; std::cout << " with grid sample points"; _test_surface_neighbors_3_cube(Dh(),rotate, 75, K2::FT(1e-2), true); std::cout << " done." << std::endl << std::endl; -// //undersampled rotated cube -// Transformation rotate3(K2::RT(0.1),K2::RT(0.4),K2::RT(0.6),K2::RT(0), -// K2::RT(0.3),K2::RT(0.5),K2::RT(0.1),K2::RT(0), -// K2::RT(0.5),K2::RT(0.9),K2::RT(0.8),K2::RT(0)); -// std::cout << "Testing surface_neighbors_3 on an undersampled ROTATED cube " -// << std::endl; -// std::cout << " with grid sample points"; -// _test_surface_neighbors_3_cube(Dh(),rotate3,75, K2::FT(9), true); -// std::cout << " done." << std::endl << std::endl; + // //undersampled rotated cube + // Transformation rotate3(K2::RT(0.1),K2::RT(0.4),K2::RT(0.6),K2::RT(0), + // K2::RT(0.3),K2::RT(0.5),K2::RT(0.1),K2::RT(0), + // K2::RT(0.5),K2::RT(0.9),K2::RT(0.8),K2::RT(0)); + // std::cout << "Testing surface_neighbors_3 on an undersampled ROTATED cube " + // << std::endl; + // std::cout << " with grid sample points"; + // _test_surface_neighbors_3_cube(Dh(),rotate3,75, K2::FT(9), true); + // std::cout << " done." << std::endl << std::endl; return 0; }