songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Improved Interpolation readability (no real changes) 

-- Removed trailing whitespace
-- Fixed (some) includes
-- Fixed indentation
-- Fixed some remaining french
This commit is contained in:
Mael Rouxel-Labbé 2017-04-18 13:59:18 +02:00
parent 8842cc9f54
commit b39201ab5c
24 changed files with 886 additions and 955 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

131
Interpolation/examples/Interpolation/interpolation_2.cpp
View File

@ -1,4 +1,4 @@
// compares the result of several interpolation methods // Compares the result of several interpolation methods
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@ -53,19 +53,17 @@ int main()
Delaunay_triangulation T; Delaunay_triangulation T;
Point_value_map values; Point_value_map values;
Point_vector_map gradients; Point_vector_map gradients;
//parameters for quadratic function: //parameters for quadratic function:
Coord_type alpha = Coord_type(1.0), Coord_type alpha = Coord_type(1.0),
beta1 = Coord_type(2.0), beta1 = Coord_type(2.0),
beta2 = Coord_type(1.0), beta2 = Coord_type(1.0),
gamma1 = Coord_type(0.3), gamma1 = Coord_type(0.3),
gamma2 = Coord_type(0.0), gamma2 = Coord_type(0.0),
gamma3 = Coord_type(0.0), gamma3 = Coord_type(0.0),
gamma4 = Coord_type(0.3); gamma4 = Coord_type(0.3);
for(int j=0; j<n ; j++){ for(int j=0; j<n ; j++){
T.insert(points[j]); T.insert(points[j]);
@ -75,9 +73,9 @@ int main()
Coord_type y(points[j].y()); Coord_type y(points[j].y());
Coord_type value = alpha + beta1*x + beta2*y + gamma1*(x*x) + Coord_type value = alpha + beta1*x + beta2*y + gamma1*(x*x) +
gamma4*(y*y) + (gamma2+ gamma3) *(x*y); gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
Vector gradient(beta1+ (gamma2+ gamma3)*y + Coord_type(2)*(gamma1*x), Vector gradient(beta1+ (gamma2+ gamma3)*y + Coord_type(2)*(gamma1*x),
beta2+ (gamma2+ gamma3)*x + Coord_type(2)*(gamma4*y)); beta2+ (gamma2+ gamma3)*x + Coord_type(2)*(gamma4*y));
values.insert(std::make_pair(points[j], value)); values.insert(std::make_pair(points[j], value));
gradients.insert(std::make_pair(points[j], gradient)); gradients.insert(std::make_pair(points[j], gradient));
} }
@ -85,11 +83,11 @@ int main()
//variables for statistics: //variables for statistics:
std::pair<Coord_type, bool> res; std::pair<Coord_type, bool> res;
Coord_type error, l_total = Coord_type(0), Coord_type error, l_total = Coord_type(0),
q_total(l_total), f_total(l_total), s_total(l_total), q_total(l_total), f_total(l_total), s_total(l_total),
ssquare_total(l_total), l_max(l_total), ssquare_total(l_total), l_max(l_total),
q_max(l_total), f_max(l_total), s_max(l_total), q_max(l_total), f_max(l_total), s_max(l_total),
ssquare_max(l_total), ssquare_max(l_total),
total_value(l_total), l_value(l_total); total_value(l_total), l_value(l_total);
int failure(0); int failure(0);
//interpolation + error statistics //interpolation + error statistics
@ -98,22 +96,22 @@ int main()
Coord_type y(points[i].y()); Coord_type y(points[i].y());
Coord_type exact_value = alpha + beta1*x + beta2*y + gamma1*(x*x) + 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; total_value += exact_value;
//Coordinate_vector: //Coordinate_vector:
std::vector< std::pair< Point, Coord_type > > coords; std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = Coord_type norm =
CGAL::natural_neighbor_coordinates_2(T, points[i], CGAL::natural_neighbor_coordinates_2(T, points[i],
std::back_inserter(coords)).second; std::back_inserter(coords)).second;
assert(norm>0); assert(norm>0);
//linear interpolant: //linear interpolant:
l_value = CGAL::linear_interpolation(coords.begin(), coords.end(), l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm, norm,
CGAL::Data_access<Point_value_map>(values)); CGAL::Data_access<Point_value_map>(values));
error = CGAL_NTS abs(l_value - exact_value); error = CGAL_NTS abs(l_value - exact_value);
l_total += error; l_total += error;
@ -121,42 +119,41 @@ int main()
//Farin interpolant: //Farin interpolant:
res = CGAL::farin_c1_interpolation(coords.begin(), res = CGAL::farin_c1_interpolation(coords.begin(),
coords.end(), norm,points[i], coords.end(), norm,points[i],
CGAL::Data_access<Point_value_map>(values), CGAL::Data_access<Point_value_map>(values),
CGAL::Data_access<Point_vector_map> CGAL::Data_access<Point_vector_map>
(gradients), (gradients),
Traits()); Traits());
if(res.second){ if(res.second){
error = CGAL_NTS abs(res.first - exact_value); error = CGAL_NTS abs(res.first - exact_value);
f_total += error; f_total += error;
if (error > f_max) f_max = error; if (error > f_max) f_max = error;
}else ++failure; } else ++failure;
//quadratic interpolant: //quadratic interpolant:
res = CGAL::quadratic_interpolation(coords.begin(), coords.end(), res = CGAL::quadratic_interpolation(coords.begin(), coords.end(),
norm,points[i], norm,points[i],
CGAL::Data_access<Point_value_map> CGAL::Data_access<Point_value_map>
(values), (values),
CGAL::Data_access<Point_vector_map> CGAL::Data_access<Point_vector_map>
(gradients), (gradients),
Traits()); Traits());
if(res.second){ if(res.second){
error = CGAL_NTS abs(res.first - exact_value); error = CGAL_NTS abs(res.first - exact_value);
q_total += error; q_total += error;
if (error > q_max) q_max = error; if (error > q_max) q_max = error;
}else ++failure; } else ++failure;
//Sibson interpolant: version without sqrt: //Sibson interpolant: version without sqrt:
res = CGAL::sibson_c1_interpolation_square(coords.begin(), res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(), norm, coords.end(), norm,
points[i], points[i],
CGAL::Data_access<Point_value_map> CGAL::Data_access<Point_value_map>
(values), (values),
CGAL::Data_access<Point_vector_map> CGAL::Data_access<Point_vector_map>
(gradients), (gradients),
Traits()); Traits());
//error statistics //error statistics
if(res.second){ if(res.second){
error = CGAL_NTS abs(res.first - exact_value); error = CGAL_NTS abs(res.first - exact_value);
ssquare_total += error; ssquare_total += error;
@ -165,13 +162,13 @@ int main()
//with sqrt(the traditional): //with sqrt(the traditional):
res = CGAL::sibson_c1_interpolation(coords.begin(), res = CGAL::sibson_c1_interpolation(coords.begin(),
coords.end(), norm, coords.end(), norm,
points[i], points[i],
CGAL::Data_access<Point_value_map> CGAL::Data_access<Point_value_map>
(values), (values),
CGAL::Data_access<Point_vector_map> CGAL::Data_access<Point_vector_map>
(gradients), (gradients),
Traits()); Traits());
//error statistics //error statistics
if(res.second){ if(res.second){
@ -184,30 +181,30 @@ int main()
/************** end of Interpolation: dump statistics **************/ /************** end of Interpolation: dump statistics **************/
std::cout << "Result: -----------------------------------" << std::endl; std::cout << "Result: -----------------------------------" << std::endl;
std::cout << "Interpolation of '" << alpha <<" + " std::cout << "Interpolation of '" << alpha <<" + "
<< beta1<<" x + " << beta1<<" x + "
<< beta2 << " y + " << gamma1 <<" x^2 + " << gamma2+ gamma3 << beta2 << " y + " << gamma1 <<" x^2 + " << gamma2+ gamma3
<<" xy + " << gamma4 << " y^2'" << std::endl; <<" xy + " << gamma4 << " y^2'" << std::endl;
std::cout << "Knowing " << m << " sample points. Interpolation on " std::cout << "Knowing " << m << " sample points. Interpolation on "
<< n <<" random points. "<< std::endl; << n <<" random points. "<< std::endl;
std::cout <<"Average function value " std::cout <<"Average function value "
<< (1.0/n) * CGAL::to_double(total_value) << (1.0/n) * CGAL::to_double(total_value)
<< ", nb of failures "<< failure << std::endl; << ", nb of failures "<< failure << std::endl;
std::cout << "linear interpolant mean error " std::cout << "linear interpolant mean error "
<< CGAL::to_double(l_total)/n << " max " << CGAL::to_double(l_total)/n << " max "
<< CGAL::to_double(l_max) <<std::endl; << CGAL::to_double(l_max) <<std::endl;
std::cout << "quadratic interpolant mean error " std::cout << "quadratic interpolant mean error "
<< CGAL::to_double(q_total)/n << " max " << CGAL::to_double(q_total)/n << " max "
<< CGAL::to_double(q_max) << std::endl; << CGAL::to_double(q_max) << std::endl;
std::cout << "Farin interpolant mean error " std::cout << "Farin interpolant mean error "
<< CGAL::to_double(f_total)/n << " max " << CGAL::to_double(f_total)/n << " max "
<< CGAL::to_double(f_max) << std::endl; << CGAL::to_double(f_max) << std::endl;
std::cout << "Sibson interpolant(classic) mean error " std::cout << "Sibson interpolant(classic) mean error "
<< CGAL::to_double(s_total)/n << " max " << CGAL::to_double(s_total)/n << " max "
<< CGAL::to_double(s_max) << std::endl; << CGAL::to_double(s_max) << std::endl;
std::cout << "Sibson interpolant(square_dist) mean error " std::cout << "Sibson interpolant(square_dist) mean error "
<< CGAL::to_double(ssquare_total)/n << " max " << CGAL::to_double(ssquare_total)/n << " max "
<< CGAL::to_double(ssquare_max) << std::endl; << CGAL::to_double(ssquare_max) << std::endl;
std::cout << "done" << std::endl; std::cout << "done" << std::endl;
return 0; return 0;

18
Interpolation/examples/Interpolation/linear_interpolation_2.cpp
View File

@ -16,29 +16,31 @@ int main()
Delaunay_triangulation T; Delaunay_triangulation T;
std::map<Point, Coord_type, K::Less_xy_2> function_values; std::map<Point, Coord_type, K::Less_xy_2> function_values;
typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > > typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > >
Value_access; Value_access;
Coord_type a(0.25), bx(1.3), by(-0.7); 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++){ for (int x=0 ; x<3 ; x++){
K::Point_2 p(x,y); K::Point_2 p(x,y);
T.insert(p); T.insert(p);
function_values.insert(std::make_pair(p,a + bx* x+ by*y)); function_values.insert(std::make_pair(p,a + bx* x+ by*y));
} }
}
//coordinate computation //coordinate computation
K::Point_2 p(1.3,0.34); K::Point_2 p(1.3,0.34);
std::vector< std::pair< Point, Coord_type > > coords; std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = Coord_type norm =
CGAL::natural_neighbor_coordinates_2 CGAL::natural_neighbor_coordinates_2
(T, p,std::back_inserter(coords)).second; (T, p,std::back_inserter(coords)).second;
Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm, norm,
Value_access(function_values)); Value_access(function_values));
std::cout << " Tested interpolation on " << p << " interpolation: " std::cout << "Tested interpolation on " << p << " interpolation: "
<< res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl; << res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl;
std::cout << "done" << std::endl; std::cout << "done" << std::endl;
return 0; return 0;

19
Interpolation/examples/Interpolation/nn_coordinates_2.cpp
View File

@ -3,9 +3,8 @@
#include <CGAL/natural_neighbor_coordinates_2.h> #include <CGAL/natural_neighbor_coordinates_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef std::vector< std::pair< K::Point_2, K::FT > > typedef std::vector< std::pair< K::Point_2, K::FT > > Point_coordinate_vector;
Point_coordinate_vector;
int main() int main()
{ {
@ -18,17 +17,17 @@ int main()
//coordinate computation //coordinate computation
K::Point_2 p(1.2, 0.7); K::Point_2 p(1.2, 0.7);
Point_coordinate_vector coords; Point_coordinate_vector coords;
CGAL::Triple< CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
K::FT, bool> result = CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords));
CGAL::natural_neighbor_coordinates_2(dt, p,
std::back_inserter(coords));
if(!result.third){ if(!result.third){
std::cout << "The coordinate computation was not successful." std::cout << "The coordinate computation was not successful."
<< std::endl; << std::endl;
std::cout << "The point (" <<p << ") lies outside the convex hull." std::cout << "The point (" <<p << ") lies outside the convex hull."
<< std::endl; << std::endl;
} }
K::FT norm = result.second; K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl; std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " <<norm << std::endl; std::cout << "Normalization factor: " <<norm << std::endl;

58
Interpolation/examples/Interpolation/nn_coordinates_3.cpp
View File

@ -42,39 +42,39 @@ int main()
pp[2]=Point3(0,0,0); //outside data/points3 convex hull pp[2]=Point3(0,0,0); //outside data/points3 convex hull
std::cout << "P2 is outside the convex hull" << std::endl; std::cout << "P2 is outside the convex hull" << std::endl;
for(int ii=0;ii<3;ii++) for(int ii=0; ii<3; ++ii)
{ {
std::vector< std::pair< Vertex_handle,NT> > coor_laplace; std::vector< std::pair< Vertex_handle,NT> > coor_laplace;
std::vector< std::pair< Vertex_handle,NT> > coor_sibson; std::vector< std::pair< Vertex_handle,NT> > coor_sibson;
NT norm_coeff_laplace, norm_coeff_sibson; NT norm_coeff_laplace, norm_coeff_sibson;
std::cout << "Point P"<< ii+1 << " : "<<pp[ii].x() << " " std::cout << "Point P"<< ii+1 << " : "<<pp[ii].x() << " "
<< pp[ii].y() << " " << pp[ii].y() << " "
<< pp[ii].z() << std::endl; << pp[ii].z() << std::endl;
laplace_natural_neighbor_coordinates_3(triangulation,pp[ii], laplace_natural_neighbor_coordinates_3(triangulation,pp[ii],
std::back_inserter(coor_laplace), std::back_inserter(coor_laplace),
norm_coeff_laplace); norm_coeff_laplace);
std::cout << "Linear combination of natural neighbors with Laplace natural coordinates"; std::cout << "Linear combination of natural neighbors with Laplace natural coordinates";
std::cout << " + correctness (ensured only with an exact number type supporting sqrt)" << std::endl; std::cout << " + correctness (ensured only with an exact number type supporting sqrt)" << std::endl;
std::cout << is_correct_natural_neighborhood(triangulation,pp[ii], std::cout << is_correct_natural_neighborhood(triangulation,pp[ii],
coor_laplace.begin(), coor_laplace.begin(),
coor_laplace.end(), coor_laplace.end(),
norm_coeff_laplace) norm_coeff_laplace)
<< std::endl; << std::endl;
sibson_natural_neighbor_coordinates_3(triangulation,pp[ii], sibson_natural_neighbor_coordinates_3(triangulation,pp[ii],
std::back_inserter(coor_sibson), std::back_inserter(coor_sibson),
norm_coeff_sibson); norm_coeff_sibson);
std::cout << "Linear combination of natural neighbors with Sibson natural coordinates" << std::endl; std::cout << "Linear combination of natural neighbors with Sibson natural coordinates" << std::endl;
std::cout << " + correctness (ensured only with an exact number type)" << std::endl; std::cout << " + correctness (ensured only with an exact number type)" << std::endl;
std::cout << is_correct_natural_neighborhood(triangulation,pp[ii], std::cout << is_correct_natural_neighborhood(triangulation,pp[ii],
coor_sibson.begin(), coor_sibson.begin(),
coor_sibson.end(), coor_sibson.end(),
norm_coeff_sibson) norm_coeff_sibson)
<< std::endl; << std::endl;
} }
std::cout << "done" << std::endl; std::cout << "done" << std::endl;
return 0; return 0;

17
Interpolation/examples/Interpolation/rn_coordinates_2.cpp
View File

@ -9,8 +9,7 @@ typedef CGAL::Regular_triangulation_euclidean_traits_2<K> Gt;
typedef CGAL::Regular_triangulation_2<Gt> Regular_triangulation; typedef CGAL::Regular_triangulation_2<Gt> Regular_triangulation;
typedef Regular_triangulation::Bare_point Bare_point; typedef Regular_triangulation::Bare_point Bare_point;
typedef Regular_triangulation::Weighted_point Weighted_point; typedef Regular_triangulation::Weighted_point Weighted_point;
typedef std::vector< std::pair< Weighted_point, K::FT > > typedef std::vector< std::pair< Weighted_point, K::FT > > Point_coordinate_vector;
Point_coordinate_vector;
int main() int main()
{ {
@ -23,17 +22,17 @@ int main()
//coordinate computation //coordinate computation
Weighted_point wp(Bare_point(1.2, 0.7),2); Weighted_point wp(Bare_point(1.2, 0.7),2);
Point_coordinate_vector coords; Point_coordinate_vector coords;
CGAL::Triple< CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
K::FT, bool> result = CGAL::regular_neighbor_coordinates_2(rt, wp, std::back_inserter(coords));
CGAL::regular_neighbor_coordinates_2(rt, wp,
std::back_inserter(coords));
if(!result.third){ if(!result.third){
std::cout << "The coordinate computation was not successful." std::cout << "The coordinate computation was not successful."
<< std::endl; << std::endl;
std::cout << "The point (" <<wp.point() << ") lies outside the convex hull." std::cout << "The point (" <<wp.point() << ") lies outside the convex hull."
<< std::endl; << std::endl;
} }
K::FT norm = result.second; K::FT norm = result.second;
std::cout << "Coordinate computation successful." << std::endl; std::cout << "Coordinate computation successful." << std::endl;
std::cout << "Normalization factor: " <<norm << std::endl; std::cout << "Normalization factor: " <<norm << std::endl;

42
Interpolation/examples/Interpolation/sibson_interpolation_2.cpp
View File

@ -24,44 +24,44 @@ int main()
//parameters for spherical function: //parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2); Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++) for (int y=0 ; y<4 ; y++){
for (int x=0 ; x<4 ; x++){ for (int x=0 ; x<4 ; x++){
K::Point_2 p(x,y); K::Point_2 p(x,y);
T.insert(p); T.insert(p);
function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y))); function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
} }
sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients, }
function_gradients.begin()),
CGAL::Data_access<Point_value_map>
(function_values),
Traits());
sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients,
function_gradients.begin()),
CGAL::Data_access<Point_value_map>(function_values),
Traits());
//coordiante computation //coordiante computation
K::Point_2 p(1.6,1.4); K::Point_2 p(1.6,1.4);
std::vector< std::pair< Point, Coord_type > > coords; std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm = Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter
CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter (coords)).second;
(coords)).second;
//Sibson interpolant: version without sqrt: //Sibson interpolant: version without sqrt:
std::pair<Coord_type, bool> res = std::pair<Coord_type, bool> res =
CGAL::sibson_c1_interpolation_square CGAL::sibson_c1_interpolation_square(
(coords.begin(), coords.begin(),
coords.end(),norm,p, coords.end(),norm,p,
CGAL::Data_access<Point_value_map>(function_values), CGAL::Data_access<Point_value_map>(function_values),
CGAL::Data_access<Point_vector_map>(function_gradients), CGAL::Data_access<Point_vector_map>(function_gradients),
Traits()); Traits());
if(res.second) if(res.second)
std::cout << " Tested interpolation on " << p std::cout << "Tested interpolation on " << p
<< " interpolation: " << res.first << " exact: " << " interpolation: " << res.first << " exact: "
<< a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y()) << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< std::endl; << std::endl;
else else
std::cout << "C^1 Interpolation not successful." << std::endl std::cout << "C^1 Interpolation not successful." << std::endl
<< " not all function_gradients are provided." << std::endl << " not all function_gradients are provided." << std::endl
<< " You may resort to linear interpolation." << std::endl; << " You may resort to linear interpolation." << std::endl;
std::cout << "done" << std::endl; std::cout << "done" << std::endl;
return 0; return 0;

50
Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp
View File

@ -8,53 +8,51 @@
#include <CGAL/surface_neighbor_coordinates_3.h> #include <CGAL/surface_neighbor_coordinates_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K; #include <iostream>
typedef K::FT Coord_type; #include <iterator>
typedef K::Point_3 Point_3; #include <vector>
typedef K::Vector_3 Vector_3;
typedef std::vector< std::pair< Point_3, K::FT > > typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
Point_coordinate_vector; 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 main()
{ {
int n=100; int n=100;
std::vector< Point_3> points; std::vector< Point_3> points;
points.reserve(n); points.reserve(n);
std::cout << "Generate " << n << " random points on a sphere." std::cout << "Generate " << n << " random points on a sphere." << std::endl;
<< std::endl;
CGAL::Random_points_on_sphere_3<Point_3> g(1); CGAL::Random_points_on_sphere_3<Point_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); Point_3 p(1, 0,0);
Vector_3 normal(p-CGAL::ORIGIN); Vector_3 normal(p - CGAL::ORIGIN);
std::cout << "Compute surface neighbor coordinates for " std::cout << "Compute surface neighbor coordinates for " << p << std::endl;
<< p << std::endl;
Point_coordinate_vector coords; Point_coordinate_vector coords;
CGAL::Triple< std::back_insert_iterator<Point_coordinate_vector>, CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
K::FT, bool> result = K::FT, bool> result =
CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(), CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(),
p, normal, p, normal,
std::back_inserter(coords), std::back_inserter(coords),
K()); K());
if(!result.third){ if(!result.third){
//Undersampling: //Undersampling:
std::cout << "The coordinate computation was not successful." std::cout << "The coordinate computation was not successful." << std::endl;
<< std::endl;
return 0; return 0;
} }
K::FT norm = result.second; K::FT norm = result.second;
std::cout << "Testing the barycentric property " << std::endl; std::cout << "Testing the barycentric property " << std::endl;
Point_3 b(0, 0,0); Point_3 b(0, 0, 0);
for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator
it = coords.begin(); it!=coords.end(); ++it) it = coords.begin(); it!=coords.end(); ++it)
b = b + (it->second/norm)* (it->first - CGAL::ORIGIN); b = b + (it->second/norm)* (it->first - CGAL::ORIGIN);
std::cout <<" weighted barycenter: " << b <<std::endl; std::cout << " weighted barycenter: " << b <<std::endl;
std::cout << " squared distance: " << std::cout << " squared distance: " << CGAL::squared_distance(p,b) << std::endl;
CGAL::squared_distance(p,b) <<std::endl;
std::cout << "done" << std::endl; std::cout << "done" << std::endl;
return 0; return 0;

40
Interpolation/include/CGAL/Interpolation_gradient_fitting_traits_2.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#ifndef CGAL_INTERPOLATION_GRADIENT_FITTING_TRAITS_2_H #ifndef CGAL_INTERPOLATION_GRADIENT_FITTING_TRAITS_2_H
@ -23,7 +19,6 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <CGAL/aff_transformation_tags.h> #include <CGAL/aff_transformation_tags.h>
namespace CGAL { namespace CGAL {
@ -44,14 +39,14 @@ public:
Aff_transformation_2 Aff_transformation_2
operator()(const Aff_transformation_2& tr1, 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), return Aff_transformation_2(tr1.m(0,0) + tr2.m(0,0),
tr1.m(0,1) + tr2.m(0,1), tr1.m(0,1) + tr2.m(0,1),
tr1.m(0,2) + tr2.m(0,2), tr1.m(0,2) + tr2.m(0,2),
tr1.m(1,0) + tr2.m(1,0), tr1.m(1,0) + tr2.m(1,0),
tr1.m(1,1) + tr2.m(1,1), tr1.m(1,1) + tr2.m(1,1),
tr1.m(1,2) + tr2.m(1,2)); tr1.m(1,2) + tr2.m(1,2));
} }
}; };
@ -91,8 +86,9 @@ public:
Aff_transformation_2 Aff_transformation_2
operator()(const Vector_2& v) const operator()(const Vector_2& v) const
{ {
return Aff_transformation_2(v.x()*v.x(),v.x()*v.y(),v.x()*v.y(), return Aff_transformation_2(v.x()*v.x(),
v.y()*v.y()); 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 typename Rep::Aff_transformation_2 Aff_transformation_d;
typedef Construct_null_matrix_2<Aff_transformation_d> typedef Construct_null_matrix_2<Aff_transformation_d>
Construct_null_matrix_d; Construct_null_matrix_d;
typedef Construct_scaling_matrix_2<Aff_transformation_d> typedef Construct_scaling_matrix_2<Aff_transformation_d>
Construct_scaling_matrix_d; Construct_scaling_matrix_d;
typedef Construct_sum_matrix_2<Aff_transformation_d> Construct_sum_matrix_d; typedef Construct_sum_matrix_2<Aff_transformation_d> Construct_sum_matrix_d;
typedef Construct_outer_product_2<Rep> Construct_outer_product_d; typedef Construct_outer_product_2<Rep> Construct_outer_product_d;
Construct_outer_product_d Construct_outer_product_d
construct_outer_product_d_object() const construct_outer_product_d_object() const
{return Construct_outer_product_d();} {return Construct_outer_product_d();}
Construct_sum_matrix_d Construct_sum_matrix_d
construct_sum_matrix_d_object() const construct_sum_matrix_d_object() const
{return Construct_sum_matrix_d();} {return Construct_sum_matrix_d();}
Construct_scaling_matrix_d Construct_scaling_matrix_d
construct_scaling_matrix_d_object() const construct_scaling_matrix_d_object() const
{return Construct_scaling_matrix_d();} {return Construct_scaling_matrix_d();}
Construct_null_matrix_d Construct_null_matrix_d
construct_null_matrix_d_object() const construct_null_matrix_d_object() const
{return Construct_null_matrix_d();} {return Construct_null_matrix_d();}
//also in the traits without gradient computation: //also in the traits without gradient computation:
Construct_scaled_vector_d Construct_scaled_vector_d
construct_scaled_vector_d_object()const construct_scaled_vector_d_object()const
{return Construct_scaled_vector_d();} {return Construct_scaled_vector_d();}
Construct_vector_d Construct_vector_d
construct_vector_d_object()const construct_vector_d_object()const
{return Construct_vector_d();} {return Construct_vector_d();}
Compute_squared_distance_d Compute_squared_distance_d
compute_squared_distance_d_object()const compute_squared_distance_d_object()const
{return Compute_squared_distance_d();} {return Compute_squared_distance_d();}
}; };
} //namespace CGAL } //namespace CGAL

5
Interpolation/include/CGAL/Interpolation_traits_2.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#ifndef CGAL_INTERPOLATION_TRAITS_2_H #ifndef CGAL_INTERPOLATION_TRAITS_2_H
@ -23,7 +19,6 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
namespace CGAL { namespace CGAL {
//-----------------------------------------------------------------------// //-----------------------------------------------------------------------//

5
Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2003 INRIA Sophia-Antipolis (France). // Copyright (c) 2003, 2017 INRIA Sophia-Antipolis (France).
// All rights reserved. // All rights reserved.
// //
// This file is part of CGAL (www.cgal.org). // This file is part of CGAL (www.cgal.org).
@ -19,7 +19,6 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <CGAL/Origin.h> #include <CGAL/Origin.h>
#include <CGAL/tags.h> #include <CGAL/tags.h>
#include <CGAL/number_utils_classes.h> #include <CGAL/number_utils_classes.h>
@ -116,7 +115,6 @@ private:
const Vector& normal; const Vector& normal;
}; };
template < typename K > template < typename K >
class Construct_plane_intersected_bisector_3 class Construct_plane_intersected_bisector_3
{ {
@ -142,7 +140,6 @@ private:
const Vector& normal; const Vector& normal;
}; };
template < typename K > template < typename K >
class Compare_first_projection_3 class Compare_first_projection_3
{ {

136
Interpolation/include/CGAL/constructions/constructions_for_voronoi_intersection_cartesian_2_3.h
View File

@ -29,50 +29,50 @@ namespace CGAL {
template < class RT> template < class RT>
void void
plane_centered_circumcenter_translateC3(const RT &ax, const RT &ay, plane_centered_circumcenter_translateC3(const RT &ax, const RT &ay,
const RT &az, const RT &az,
const RT &nx, const RT &ny, const RT &nx, const RT &ny,
const RT &nz, const RT &nz,
const RT &qx, const RT &qy, const RT &qx, const RT &qy,
const RT &qz, const RT &qz,
const RT &rx, const RT &ry, const RT &rx, const RT &ry,
const RT &rz, const RT &rz,
RT &x, RT &y, RT &z) RT &x, RT &y, RT &z)
{ {
RT den = RT(2) * determinant(nx,qx,rx, RT den = RT(2) * determinant(nx,qx,rx,
ny,qy,ry, ny,qy,ry,
nz,qz,rz); nz,qz,rz);
// The 3 points aren't collinear. // The 3 points aren't collinear.
// Hopefully, this is already checked at the upper level. // Hopefully, this is already checked at the upper level.
CGAL_assertion ( den != RT(0) ); CGAL_assertion ( den != RT(0) );
RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) + 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) + 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; RT na = nx*ax + ny*ay + nz*az;
na *= RT(2.0); na *= RT(2.0);
x = determinant(ny,nz,na, x = determinant(ny,nz,na,
qy,qz,q2, qy,qz,q2,
ry,rz,r2)/ den ; ry,rz,r2)/ den ;
y = - determinant(nx,nz,na, y = - determinant(nx,nz,na,
qx,qz,q2, qx,qz,q2,
rx,rz,r2)/ den ; rx,rz,r2)/ den ;
z = determinant(nx,ny,na, z = determinant(nx,ny,na,
qx,qy,q2, qx,qy,q2,
rx,ry,r2)/ den ; rx,ry,r2)/ den ;
} }
template < class RT> template < class RT>
void void
plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az, plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az,
const RT &nx, const RT &ny, const RT &nz, const RT &nx, const RT &ny, const RT &nz,
const RT &px, const RT &py, const RT &pz, const RT &px, const RT &py, const RT &pz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
const RT &rx, const RT &ry, const RT &rz, const RT &rx, const RT &ry, const RT &rz,
RT &x, RT &y, RT &z) RT &x, RT &y, RT &z)
{ {
// resolution of the system (where we note c the center) // 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: //method:
// - tranlation of p to the origin. // - tranlation of p to the origin.
plane_centered_circumcenter_translateC3(ax-px, ay-py, az-pz, plane_centered_circumcenter_translateC3(ax-px, ay-py, az-pz,
nx, ny, nz, nx, ny, nz,
qx-px, qy-py,qz-pz, qx-px, qy-py,qz-pz,
rx-px, ry-py,rz-pz, rx-px, ry-py,rz-pz,
x, y, z); x, y, z);
x+=px; x+=px;
y+=py; y+=py;
z+=pz; z+=pz;
@ -95,13 +95,13 @@ plane_centered_circumcenterC3(const RT &ax, const RT &ay, const RT &az,
template < class RT> template < class RT>
void void
bisector_plane_intersection_translateC3(const RT &ax, const RT &ay, bisector_plane_intersection_translateC3(const RT &ax, const RT &ay,
const RT &az, const RT &az,
const RT &nx, const RT &ny, const RT &nx, const RT &ny,
const RT &nz, const RT &nz,
const RT &qx, const RT &qy, const RT &qx, const RT &qy,
const RT &qz, const RT& den, const RT &qz, const RT& den,
RT &x1, RT &y1, RT &x2, RT RT &x1, RT &y1, RT &x2, RT
&y2, bool& swapped) &y2, bool& swapped)
{ {
// c: a point on l must be the center of a sphere passing // c: a point on l must be the center of a sphere passing
// through p and q, c lies in h. 2 equations: // 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); // where RT den = RT(2.0) * determinant(qx,qy,nx, ny);
RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) 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; RT na = nx*ax + ny*ay + nz*az;
na *= RT(2.0); 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 not: permutation of p1 and p2
if((sign_of_determinant(qx,qy,qz, x1,y1,RT(0),x2 ,y2,RT(1)) if((sign_of_determinant(qx,qy,qz, x1,y1,RT(0),x2 ,y2,RT(1))
* CGAL_NTS sign (-na)) > 0 ) * CGAL_NTS sign (-na)) > 0 )
{ {
RT x3(x1),y3(y1); RT x3(x1),y3(y1);
x1 =x2; x1 =x2;
y1 =y2; y1 =y2;
x2 = x3; x2 = x3;
y2 = y3; y2 = y3;
swapped =true; swapped =true;
} }
} }
template < class RT> template < class RT>
void void
bisector_plane_intersection_permuteC3(const RT &ax, const RT &ay, bisector_plane_intersection_permuteC3(const RT &ax, const RT &ay,
const RT &az, const RT &az,
const RT &nx, const RT &ny, const RT &nx, const RT &ny,
const RT &nz, const RT &nz,
const RT &px, const RT &py, const RT &px, const RT &py,
const RT &pz, const RT &pz,
const RT &qx, const RT &qy, const RT &qx, const RT &qy,
const RT &qz, const RT &qz,
const RT &den, const RT &den,
RT &x1, RT &y1, RT& z1, RT &x1, RT &y1, RT& z1,
RT &x2, RT &y2, RT& z2) RT &x2, RT &y2, RT& z2)
{ {
//translation of p to the origin //translation of p to the origin
bool swapped =false; bool swapped =false;
CGAL_precondition((nx!=RT(0) || ny!=RT(0)) && (qx!=px || qy!=py) 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, bisector_plane_intersection_translateC3(ax-px, ay-py, az-pz,
nx, ny, nz, nx, ny, nz,
qx-px, qy-py,qz-pz,den, qx-px, qy-py,qz-pz,den,
x1, y1,x2,y2,swapped); x1, y1,x2,y2,swapped);
// re-translation of the origin to p: // re-translation of the origin to p:
x1+=px; x1+=px;
y1+=py; y1+=py;
@ -183,11 +183,11 @@ bisector_plane_intersection_permuteC3(const RT &ax, const RT &ay,
template < class RT> template < class RT>
void void
bisector_plane_intersectionC3(const RT &ax, const RT &ay, const RT &az, bisector_plane_intersectionC3(const RT &ax, const RT &ay, const RT &az,
const RT &nx, const RT &ny, const RT &nz, const RT &nx, const RT &ny, const RT &nz,
const RT &px, const RT &py, const RT &pz, const RT &px, const RT &py, const RT &pz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
RT &x1, RT &y1, RT& z1, RT &x1, RT &y1, RT& z1,
RT &x2, RT &y2, RT& z2) RT &x2, RT &y2, RT& z2)
{ {
// constructs the line l = (p1,p2)= ((x1,y1,z1),(x2,y2,z2)) // constructs the line l = (p1,p2)= ((x1,y1,z1),(x2,y2,z2))
// the intersection line between the bisector of (p,q) and // 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 //den==0 <=> projections of (qx,qy) and (nx,ny) are identical
//intersection with z=0/z=1 //intersection with z=0/z=1
bisector_plane_intersection_permuteC3(ax,ay,az,nx,ny,nz,px,py,pz, bisector_plane_intersection_permuteC3(ax,ay,az,nx,ny,nz,px,py,pz,
qx,qy,qz,den, qx,qy,qz,den,
x1,y1,z1,x2,y2,z2); x1,y1,z1,x2,y2,z2);
else{ else{
den = RT(2.0) * determinant(qy-py,qz-pz,ny,nz); den = RT(2.0) * determinant(qy-py,qz-pz,ny,nz);
if ((ny!=0 || nz!=0) && (qy!=py || qz!=pz) && den!=RT(0)) if ((ny!=0 || nz!=0) && (qy!=py || qz!=pz) && den!=RT(0))
//intersection with x=0/x=1 => permutations //intersection with x=0/x=1 => permutations
bisector_plane_intersection_permuteC3(ay,az,ax,ny,nz,nx,py,pz,px, bisector_plane_intersection_permuteC3(ay,az,ax,ny,nz,nx,py,pz,px,
qy,qz,qx,den, qy,qz,qx,den,
y1,z1,x1,y2,z2,x2); y1,z1,x1,y2,z2,x2);
else{ else{
den = RT(2.0) * determinant(qz-pz,qx-px,nz,nx); den = RT(2.0) * determinant(qz-pz,qx-px,nz,nx);
CGAL_assertion((nx!=0 || nz!=0) && (qx!=px || qz!=pz) && den!=RT(0)); CGAL_assertion((nx!=0 || nz!=0) && (qx!=px || qz!=pz) && den!=RT(0));
//intersection with y=0/y=1 => permutations //intersection with y=0/y=1 => permutations
bisector_plane_intersection_permuteC3(az,ax,ay,nz,nx,ny,pz,px,py, bisector_plane_intersection_permuteC3(az,ax,ay,nz,nx,ny,pz,px,py,
qz,qx,qy,den, qz,qx,qy,den,
z1,x1,y1,z2,x2,y2); z1,x1,y1,z2,x2,y2);
} }
} }
} }

193
Interpolation/include/CGAL/interpolation_functions.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#ifndef CGAL_INTERPOLATION_FUNCTIONS_H #ifndef CGAL_INTERPOLATION_FUNCTIONS_H
@ -23,24 +19,25 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/double.h> #include <CGAL/double.h>
#include <CGAL/use.h> #include <CGAL/use.h>
#include <iterator>
#include <utility>
#include <vector> #include <vector>
namespace CGAL { namespace CGAL {
//Functor class for accessing the function values/gradients //Functor class for accessing the function values/gradients
template< class Map > template< class Map >
struct Data_access : public std::unary_function< typename Map::key_type, struct Data_access
std::pair< typename Map::mapped_type, bool> > : 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::mapped_type Data_type;
typedef typename Map::key_type Key_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> std::pair< Data_type, bool>
operator()(const Key_type& p) const { 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()) if(mit!= map.end())
return std::make_pair(mit->second, true); return std::make_pair(mit->second, true);
return std::make_pair(Data_type(), false); return std::make_pair(Data_type(), false);
}; }
const Map& map; const Map& map;
}; };
@ -57,9 +54,10 @@ struct Data_access : public std::unary_function< typename Map::key_type,
template < class ForwardIterator, class Functor> template < class ForwardIterator, class Functor>
typename Functor::result_type::first_type typename Functor::result_type::first_type
linear_interpolation(ForwardIterator first, ForwardIterator beyond, linear_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename const typename
std::iterator_traits<ForwardIterator>::value_type:: std::iterator_traits<ForwardIterator>::value_type::
second_type& norm, Functor function_value) second_type& norm,
Functor function_value)
{ {
CGAL_precondition(norm>0); CGAL_precondition(norm>0);
typedef typename Functor::result_type::first_type Value_type; 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, template < class ForwardIterator, class Functor, class GradFunctor, class Traits>
class Traits>
std::pair< typename Functor::result_type::first_type, bool> std::pair< typename Functor::result_type::first_type, bool>
quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, quadratic_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename const typename
std::iterator_traits<ForwardIterator>:: std::iterator_traits<ForwardIterator>::
value_type::second_type& norm, const typename value_type::second_type& norm,
std::iterator_traits<ForwardIterator>::value_type:: const typename
first_type& p, Functor function_value, std::iterator_traits<ForwardIterator>::value_type::
GradFunctor function_gradient, first_type& p,
const Traits& traits) Functor function_value,
GradFunctor function_gradient,
const Traits& traits)
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
typedef typename Functor::result_type::first_type Value_type; typedef typename Functor::result_type::first_type Value_type;
@ -99,27 +98,26 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond,
if(!grad.second) if(!grad.second)
return std::make_pair(Value_type(0), false); return std::make_pair(Value_type(0), false);
result += (first->second/norm) result += (first->second/norm)
*( f.first + grad.first* *( f.first + grad.first*
traits.construct_scaled_vector_d_object() traits.construct_scaled_vector_d_object()
(traits.construct_vector_d_object()(first->first, p),0.5)); (traits.construct_vector_d_object()(first->first, p),0.5));
} }
return std::make_pair(result, true); return std::make_pair(result, true);
} }
template < class ForwardIterator, class Functor, class GradFunctor, template < class ForwardIterator, class Functor, class GradFunctor, class Traits>
class Traits>
std::pair< typename Functor::result_type::first_type, bool> std::pair< typename Functor::result_type::first_type, bool>
sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
const typename const typename
std::iterator_traits<ForwardIterator>:: std::iterator_traits<ForwardIterator>::
value_type::second_type& value_type::second_type& norm,
norm, const typename const typename
std::iterator_traits<ForwardIterator>::value_type:: std::iterator_traits<ForwardIterator>::value_type::
first_type& p, first_type& p,
Functor function_value, Functor function_value,
GradFunctor function_gradient, GradFunctor function_gradient,
const Traits& traits) const Traits& traits)
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
typedef typename Functor::result_type::first_type Value_type; 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 coeff = first->second/norm;
Coord_type squared_dist = traits. 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); Coord_type dist = CGAL_NTS sqrt(squared_dist);
if(squared_dist ==0){ if(squared_dist ==0){
@ -151,22 +149,22 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
} }
//three different terms to mix linear and gradient //three different terms to mix linear and gradient
//interpolation //interpolation
term1 += coeff/dist; term1 += coeff / dist;
term2 += coeff * squared_dist; term2 += coeff * squared_dist;
term3 += coeff * dist; term3 += coeff * dist;
linear_int += coeff * f.first; linear_int += coeff * f.first;
gradient_int += (coeff/dist) gradient_int += (coeff/dist)
* (f.first + grad.first * * (f.first + grad.first *
traits.construct_vector_d_object()(first->first, p)); traits.construct_vector_d_object()(first->first, p));
} }
term4 = term3/ term1; term4 = term3/ term1;
gradient_int = gradient_int / term1; gradient_int = gradient_int / term1;
return std::make_pair((term4* linear_int + term2 * gradient_int)/ return std::make_pair((term4* linear_int + term2 * gradient_int)/
(term4 + term2), true); (term4 + term2), true);
} }
//this method works with rational number types: //this method works with rational number types:
@ -182,19 +180,18 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond,
// (vh->get_value()+ vh->get_gradient() // (vh->get_value()+ vh->get_gradient()
// *(p - vh->point())); // *(p - vh->point()));
template < class ForwardIterator, class Functor, class GradFunctor, template < class ForwardIterator, class Functor, class GradFunctor, class Traits>
class Traits>
std::pair< typename Functor::result_type::first_type, bool> std::pair< typename Functor::result_type::first_type, bool>
sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond,
beyond, const typename const typename
std::iterator_traits<ForwardIterator>:: std::iterator_traits<ForwardIterator>::
value_type::second_type& norm, value_type::second_type& norm,
const typename const typename
std::iterator_traits<ForwardIterator>:: std::iterator_traits<ForwardIterator>::
value_type::first_type& p, value_type::first_type& p,
Functor function_value, Functor function_value,
GradFunctor function_gradient, GradFunctor function_gradient,
const Traits& traits) const Traits& traits)
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
typedef typename Functor::result_type::first_type Value_type; 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 Functor::result_type f;
typename GradFunctor::result_type grad; typename GradFunctor::result_type grad;
for(; first !=beyond; ++first){ for(; first!=beyond; ++first){
f = function_value(first->first); f = function_value(first->first);
grad = function_gradient(first->first); grad = function_gradient(first->first);
CGAL_assertion(f.second); CGAL_assertion(f.second);
@ -215,7 +212,7 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator
Coord_type coeff = first->second/norm; Coord_type coeff = first->second/norm;
Coord_type squared_dist = traits. 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){ if(squared_dist ==0){
ForwardIterator it = first; ForwardIterator it = first;
@ -225,38 +222,38 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator
} }
//three different terms to mix linear and gradient //three different terms to mix linear and gradient
//interpolation //interpolation
term1 += coeff/squared_dist; term1 += coeff / squared_dist;
term2 += coeff * squared_dist; term2 += coeff * squared_dist;
term3 += coeff; term3 += coeff;
linear_int += coeff * f.first; linear_int += coeff * f.first;
gradient_int += (coeff/squared_dist) gradient_int += (coeff/squared_dist) * (f.first + grad.first *
*(f.first + grad.first* traits.construct_vector_d_object()(first->first, p));
traits.construct_vector_d_object()(first->first, p));
} }
term4 = term3/ term1; term4 = term3/ term1;
gradient_int = gradient_int / term1; gradient_int = gradient_int / term1;
return std::make_pair((term4* linear_int + term2 * gradient_int)/ return std::make_pair((term4 * linear_int + term2 * gradient_int)/
(term4 + term2), true); (term4 + term2), true);
} }
template < class RandomAccessIterator, class Functor, class template < class RandomAccessIterator, class Functor, class
GradFunctor, class Traits> GradFunctor, class Traits>
std::pair< typename Functor::result_type::first_type, bool> std::pair< typename Functor::result_type::first_type, bool>
farin_c1_interpolation(RandomAccessIterator first, farin_c1_interpolation(RandomAccessIterator first,
RandomAccessIterator beyond, RandomAccessIterator beyond,
const typename const typename
std::iterator_traits<RandomAccessIterator>:: std::iterator_traits<RandomAccessIterator>::
value_type::second_type& norm, const typename value_type::second_type& norm,
std::iterator_traits<RandomAccessIterator>:: const typename
value_type::first_type& /*p*/, std::iterator_traits<RandomAccessIterator>::
Functor function_value, GradFunctor value_type::first_type& /*p*/,
function_gradient, Functor function_value, GradFunctor
const Traits& traits) function_gradient,
const Traits& traits)
{ {
CGAL_precondition(norm >0); CGAL_precondition(norm >0);
//the function value is available for all points //the function value is available for all points
@ -283,7 +280,7 @@ farin_c1_interpolation(RandomAccessIterator first,
const Coord_type fac3(3); const Coord_type fac3(3);
std::vector< std::vector<Value_type> > std::vector< std::vector<Value_type> >
ordinates(n,std::vector<Value_type>(n, Value_type(0))); ordinates(n,std::vector<Value_type>(n, Value_type(0)));
for(int i =0; i<n; ++i){ for(int i =0; i<n; ++i){
it = first+i; it = first+i;
@ -301,21 +298,21 @@ farin_c1_interpolation(RandomAccessIterator first,
Value_type res_i(0); Value_type res_i(0);
for(int j =0; j<n; ++j){ for(int j =0; j<n; ++j){
if(i!=j){ if(i!=j){
it2 = first+j; it2 = first+j;
grad = function_gradient(it->first); grad = function_gradient(it->first);
if(!grad.second) if(!grad.second)
//the gradient is not known //the gradient is not known
return std::make_pair(Value_type(0), false); return std::make_pair(Value_type(0), false);
//ordinates[i][j] = (p_j - p_i) * g_i //ordinates[i][j] = (p_j - p_i) * g_i
ordinates[i][j] = grad.first * ordinates[i][j] = grad.first *
traits.construct_vector_d_object()(it->first,it2->first); traits.construct_vector_d_object()(it->first,it2->first);
// a point in the tangent plane: // a point in the tangent plane:
// 3( f(p_i) + (1/3)(p_j - p_i) * g_i) // 3( f(p_i) + (1/3)(p_j - p_i) * g_i)
// => 3*f(p_i) + (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 += (fac3 * ordinates[i][i] + ordinates[i][j])* it2->second;
} }
} }
//res_i already multiplied by three //res_i already multiplied by three
@ -326,20 +323,20 @@ farin_c1_interpolation(RandomAccessIterator first,
for(int i=0; i< n; ++i) for(int i=0; i< n; ++i)
for(int j=i+1; j< n; ++j) for(int j=i+1; j< n; ++j)
for(int k=j+1; k<n; ++k){ for(int k=j+1; k<n; ++k){
// add 6* (u_i*u_j*u_k) * b_ijk // add 6* (u_i*u_j*u_k) * b_ijk
// b_ijk = 1.5 * a - 0.5*c // b_ijk = 1.5 * a - 0.5*c
//where //where
//c : average of the three data control points //c : average of the three data control points
//a : 1.5*a = 1/12 * (ord[i][j] + ord[i][k] + ord[j][i] + //a : 1.5*a = 1/12 * (ord[i][j] + ord[i][k] + ord[j][i] +
// ord[j][k] + ord[k][i]+ ord[k][j]) // ord[j][k] + ord[k][i]+ ord[k][j])
// => 6 * b_ijk = 3*(f_i + f_j + f_k) + 0.5*a // => 6 * b_ijk = 3*(f_i + f_j + f_k) + 0.5*a
result += (Coord_type(2.0)*( ordinates[i][i]+ ordinates[j][j]+ result += (Coord_type(2.0)*( ordinates[i][i]+ ordinates[j][j]+
ordinates[k][k]) ordinates[k][k])
+ Coord_type(0.5)*(ordinates[i][j] + ordinates[i][k] + Coord_type(0.5)*(ordinates[i][j] + ordinates[i][k]
+ ordinates[j][i] + + ordinates[j][i] +
ordinates[j][k] + ordinates[k][i]+ ordinates[j][k] + ordinates[k][i]+
ordinates[k][j])) ordinates[k][j]))
*(first+i)->second *(first+j)->second *(first+k)->second ; *(first+i)->second *(first+j)->second *(first+k)->second ;
} }
return std::make_pair(result/(CGAL_NTS square(norm)*norm), true); return std::make_pair(result/(CGAL_NTS square(norm)*norm), true);

206
Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Frank Da, Julia Floetotto // Author(s) : Frank Da, Julia Floetotto
#ifndef CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H #ifndef CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H
@ -23,13 +19,14 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/Iterator_project.h> #include <CGAL/Iterator_project.h>
#include <CGAL/Polygon_2.h> #include <CGAL/Polygon_2.h>
#include <CGAL/number_utils_classes.h> #include <CGAL/number_utils_classes.h>
#include <CGAL/utility.h> #include <CGAL/utility.h>
#include <list>
#include <utility>
namespace CGAL { namespace CGAL {
// the struct "Project_vertex_output_iterator" // the struct "Project_vertex_output_iterator"
@ -46,8 +43,8 @@ struct Project_vertex_output_iterator
// std::pair<Vertex_handle,T> // std::pair<Vertex_handle,T>
// into an output iterator with value type std::pair<Point, T> // into an output iterator with value type std::pair<Point, T>
// Conditions: OutputIterator has value type std::pair<Vertex_handle, T> // Conditions: OutputIterator has value type std::pair<Vertex_handle, T>
// and Vertex_handle has a function ->point() // and Vertex_handle has a function ->point()
// with return type const Point& // with return type const Point&
OutputIterator _base; OutputIterator _base;
@ -75,29 +72,27 @@ struct Project_vertex_output_iterator
// template <class Dt, class OutputIterator, class Traits> // template <class Dt, class OutputIterator, class Traits>
// Triple< OutputIterator, typename Traits::FT, bool > // Triple< OutputIterator, typename Traits::FT, bool >
// natural_neighbor_coordinates_2(const Dt& dt, // natural_neighbor_coordinates_2(const Dt& dt,
// const typename Traits::Point_2& p, // const typename Traits::Point_2& p,
// OutputIterator out, const Traits& traits, // OutputIterator out, const Traits& traits,
// typename Dt::Face_handle start // typename Dt::Face_handle start
// = typename Dt::Face_handle()) // = typename Dt::Face_handle())
//
//template <class Dt, class OutputIterator, class Traits> //template <class Dt, class OutputIterator, class Traits>
//Triple< OutputIterator, typename Traits::FT, bool > //Triple< OutputIterator, typename Traits::FT, bool >
//natural_neighbor_coordinates_2(const Dt& dt, //natural_neighbor_coordinates_2(const Dt& dt,
// typename Dt::Vertex_handle vh, // typename Dt::Vertex_handle vh,
// OutputIterator out, const Traits& traits) // OutputIterator out, const Traits& traits)
//the following two functions suppose that //the following two functions suppose that
// OutputIterator has value type // OutputIterator has value type
// std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT> // std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
//!!!they are not documented!!! //!!!they are not documented!!!
template <class Dt, class OutputIterator> template < class Dt, class OutputIterator >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_vertex_2(const Dt& dt, natural_neighbor_coordinates_vertex_2(const Dt& dt,
const typename Dt::Geom_traits::Point_2& p, const typename Dt::Geom_traits::Point_2& p,
OutputIterator out, typename Dt::Face_handle start OutputIterator out, typename Dt::Face_handle start
= typename Dt::Face_handle()) = typename Dt::Face_handle())
{ {
typedef typename Dt::Geom_traits Traits; typedef typename Dt::Geom_traits Traits;
typedef typename Traits::FT Coord_type; typedef typename Traits::FT Coord_type;
@ -114,40 +109,36 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt,
int li; int li;
Face_handle fh = dt.locate(p, lt, li, start); Face_handle fh = dt.locate(p, lt, li, start);
if (lt == Dt::OUTSIDE_AFFINE_HULL if (lt == Dt::OUTSIDE_AFFINE_HULL || lt == Dt::OUTSIDE_CONVEX_HULL)
|| lt == Dt::OUTSIDE_CONVEX_HULL)
{ {
return make_triple(out, Coord_type(1), false); return make_triple(out, Coord_type(1), false);
} }
if ((lt == Dt::EDGE && if ((lt == Dt::EDGE &&
(dt.is_infinite(fh) || (dt.is_infinite(fh) || dt.is_infinite(fh->neighbor(li)))))
dt.is_infinite(fh->neighbor(li)))))
{ {
Vertex_handle v1 = fh->vertex(dt.cw(li)); Vertex_handle v1 = fh->vertex(dt.cw(li));
Vertex_handle v2 = fh->vertex(dt.ccw(li)); Vertex_handle v2 = fh->vertex(dt.ccw(li));
Point_2 p1(v1->point()),p2(v2->point()); Point_2 p1(v1->point()), p2(v2->point());
Coord_type coef1(0); Coord_type coef1(0);
Coord_type coef2(0); Coord_type coef2(0);
Equal_x_2 equal_x_2; Equal_x_2 equal_x_2;
if(!equal_x_2(p1,p2)) if(!equal_x_2(p1,p2))
{ {
coef1 = (p.x() - p2.x())/(p1.x() - p2.x()) ; coef1 = (p.x() - p2.x()) / (p1.x() - p2.x());
coef2 = 1-coef1; coef2 = 1 - coef1;
*out++= std::make_pair(v1,coef1); *out++ = std::make_pair(v1,coef1);
*out++= std::make_pair(v2,coef2); *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);
}
}else{ return make_triple(out, coef1+coef2, true);
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) if (lt == Dt::VERTEX)
@ -157,22 +148,21 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt,
} }
std::list<Edge> hole; std::list<Edge> hole;
dt.get_boundary_of_conflicts(p, std::back_inserter(hole), fh, false); dt.get_boundary_of_conflicts(p, std::back_inserter(hole), fh, false);
return natural_neighbor_coordinates_vertex_2 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: //function call if the conflict zone is known:
// OutputIterator has value type // OutputIterator has value type
// std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT> // std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
template <class Dt, class OutputIterator, class EdgeIterator > template < class Dt, class OutputIterator, class EdgeIterator >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_vertex_2(const Dt& dt, natural_neighbor_coordinates_vertex_2(const Dt& dt,
const typename Dt::Geom_traits::Point_2& p, const typename Dt::Geom_traits::Point_2& p,
OutputIterator out, EdgeIterator OutputIterator out, EdgeIterator
hole_begin, EdgeIterator hole_end) hole_begin, EdgeIterator hole_end)
{ {
CGAL_precondition(dt.dimension()==2); CGAL_precondition(dt.dimension()==2);
//precondition: p must lie inside the hole //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::Vertex_handle Vertex_handle;
typedef typename Dt::Face_circulator Face_circulator; typedef typename Dt::Face_circulator Face_circulator;
std::vector<Point_2> vor(3); std::vector<Point_2> vor(3);
Coord_type area_sum(0); Coord_type area_sum(0);
@ -196,42 +185,39 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt,
hit = hole_begin; hit = hole_begin;
while (hit != hole_end) 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; ++fc;
vor[1] = dt.dual(fc); vor[2] = dt.dual(fc);
area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits());
while(!fc->has_vertex(prev)) vor[1] = vor[2];
{
++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.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); return make_triple(out, area_sum, true);
} }
@ -242,23 +228,22 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt,
//=> OutputIterator has value type //=> OutputIterator has value type
// std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT>
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
template <class Dt, class OutputIterator> template < class Dt, class OutputIterator >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_2(const Dt& dt, natural_neighbor_coordinates_2(const Dt& dt,
const typename Dt::Geom_traits::Point_2& p, const typename Dt::Geom_traits::Point_2& p,
OutputIterator out, OutputIterator out,
typename Dt::Face_handle start = typename Dt::Face_handle start =
CGAL_TYPENAME_DEFAULT_ARG Dt::Face_handle() ) CGAL_TYPENAME_DEFAULT_ARG Dt::Face_handle() )
{ {
CGAL_precondition(dt.dimension() == 2); CGAL_precondition(dt.dimension() == 2);
Project_vertex_output_iterator<OutputIterator> op(out); Project_vertex_output_iterator<OutputIterator> op(out);
Triple< Project_vertex_output_iterator<OutputIterator>, Triple<Project_vertex_output_iterator<OutputIterator>,
typename Dt::Geom_traits::FT, bool > result = typename Dt::Geom_traits::FT, bool > result =
natural_neighbor_coordinates_vertex_2 natural_neighbor_coordinates_vertex_2(dt, p, op, start);
(dt, p, op, start);
return make_triple(result.first.base(), result.second, result.third); 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 //OutputIterator has value type
// std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT>
//function call if the conflict zone is known: //function call if the conflict zone is known:
template <class Dt, class OutputIterator, class EdgeIterator > template < class Dt, class OutputIterator, class EdgeIterator >
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_2(const Dt& dt, natural_neighbor_coordinates_2(const Dt& dt,
const typename Dt::Geom_traits::Point_2& p, const typename Dt::Geom_traits::Point_2& p,
OutputIterator out, EdgeIterator OutputIterator out, EdgeIterator
hole_begin, EdgeIterator hole_end) hole_begin, EdgeIterator hole_end)
{ {
CGAL_precondition(dt.dimension() == 2); CGAL_precondition(dt.dimension() == 2);
Project_vertex_output_iterator<OutputIterator> op(out); Project_vertex_output_iterator<OutputIterator> op(out);
Triple< Project_vertex_output_iterator<OutputIterator>, Triple<Project_vertex_output_iterator<OutputIterator>,
typename Dt::Geom_traits::FT, bool > result = typename Dt::Geom_traits::FT, bool > result =
natural_neighbor_coordinates_vertex_2 natural_neighbor_coordinates_vertex_2(dt, p, op, hole_begin,hole_end);
(dt, p, op, hole_begin,hole_end);
return make_triple(result.first.base(), result.second, result.third); return make_triple(result.first.base(), result.second, result.third);
} }
@ -293,8 +277,8 @@ natural_neighbor_coordinates_2(const Dt& dt,
template <class Dt, class OutputIterator> template <class Dt, class OutputIterator>
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
natural_neighbor_coordinates_2(const Dt& dt, natural_neighbor_coordinates_2(const Dt& dt,
typename Dt::Vertex_handle vh, typename Dt::Vertex_handle vh,
OutputIterator out) OutputIterator out)
{ {
//this functions creates a small triangulation of the //this functions creates a small triangulation of the
// incident vertices of this vertex and computes the // incident vertices of this vertex and computes the
@ -305,26 +289,26 @@ natural_neighbor_coordinates_2(const Dt& dt,
Dt t2; Dt t2;
Vertex_circulator vc = dt.incident_vertices(vh), Vertex_circulator vc = dt.incident_vertices(vh),
done(vc); done(vc);
do{ do{
CGAL_assertion(!dt.is_infinite(vc)); CGAL_assertion(!dt.is_infinite(vc));
t2.insert(vc->point()); t2.insert(vc->point());
} }
while(++vc!=done); while(++vc!=done);
return natural_neighbor_coordinates_2(t2, vh->point(), out); return natural_neighbor_coordinates_2(t2, vh->point(), out);
} }
//class providing a function object: //class providing a function object:
//OutputIterator has value type //OutputIterator has value type
// std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT> // std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT>
template <class Dt, class OutputIterator> template < class Dt, class OutputIterator >
class natural_neighbor_coordinates_2_object class natural_neighbor_coordinates_2_object
{ {
public: public:
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
operator()(const Dt& dt, operator()(const Dt& dt,
typename Dt::Vertex_handle vh, typename Dt::Vertex_handle vh,
OutputIterator out) const OutputIterator out) const
{ {
return natural_neighbor_coordinates_2(dt, vh, out); return natural_neighbor_coordinates_2(dt, vh, out);
} }

98
Interpolation/include/CGAL/natural_neighbor_coordinates_3.h
View File

@ -23,16 +23,18 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <set>
#include <vector>
#include <CGAL/tags.h> #include <CGAL/tags.h>
#include <CGAL/iterator.h> #include <CGAL/iterator.h>
#include <CGAL/utility.h> #include <CGAL/utility.h>
#include <CGAL/triangulation_assertions.h> #include <CGAL/triangulation_assertions.h>
#include <CGAL/number_utils.h> #include <CGAL/number_utils.h>
#include <algorithm>
#include <iostream> //TO DO : to remove #include <iostream> //TO DO : to remove
#include <map>
#include <set>
#include <utility>
#include <vector>
namespace CGAL { namespace CGAL {
@ -69,12 +71,13 @@ construct_circumcenter(const typename DT::Facet& f,
template <class Dt, class OutputIterator> template <class Dt, class OutputIterator>
Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT> Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
typename Dt::Geom_traits::FT, // Should provide 0 and 1 typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool > bool >
laplace_natural_neighbor_coordinates_3(const Dt& dt, laplace_natural_neighbor_coordinates_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& Q, const typename Dt::Geom_traits::Point_3& Q,
OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff, OutputIterator nn_out,
const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) 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 Dt::Geom_traits Gt;
typedef typename Gt::Point_3 Point; typedef typename Gt::Point_3 Point;
@ -84,19 +87,22 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
typedef typename Dt::Locate_type Locate_type; typedef typename Dt::Locate_type Locate_type;
typedef typename Gt::FT Coord_type; typedef typename Gt::FT Coord_type;
CGAL_triangulation_precondition (dt.dimension()== 3); 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); Cell_handle c = dt.locate( Q, lt, li, lj, start);
if ( lt == Dt::VERTEX ) if ( lt == Dt::VERTEX )
{ {
*nn_out++= std::make_pair(c->vertex(li),Coord_type(1)); *nn_out++= std::make_pair(c->vertex(li), Coord_type(1));
return make_triple(nn_out,norm_coeff=Coord_type(1),true); return make_triple(nn_out, norm_coeff = Coord_type(1),true);
} }
else if (dt.is_infinite(c)) 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<Cell_handle> cells; std::set<Cell_handle> cells;
// To replace the forbidden access to the "in conflict" flag : // To replace the forbidden access to the "in conflict" flag :
@ -105,8 +111,8 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
typename std::vector<Facet>::iterator bound_it; typename std::vector<Facet>::iterator bound_it;
// Find the cells in conflict with Q // Find the cells in conflict with Q
dt.find_conflicts(Q, c, dt.find_conflicts(Q, c,
std::back_inserter(bound_facets), std::back_inserter(bound_facets),
std::inserter(cells,cells.begin())); std::inserter(cells,cells.begin()));
std::map<Vertex_handle,Coord_type> coordinate; std::map<Vertex_handle,Coord_type> coordinate;
typename std::map<Vertex_handle,Coord_type>::iterator coor_it; typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
@ -165,12 +171,13 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
template <class Dt, class OutputIterator> template <class Dt, class OutputIterator>
Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT> Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
typename Dt::Geom_traits::FT, // Should provide 0 and 1 typename Dt::Geom_traits::FT, // Should provide 0 and 1
bool > bool >
sibson_natural_neighbor_coordinates_3(const Dt& dt, sibson_natural_neighbor_coordinates_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& Q, const typename Dt::Geom_traits::Point_3& Q,
OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff, OutputIterator nn_out,
const typename Dt::Cell_handle start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) 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 Dt::Geom_traits Gt;
typedef typename Gt::Point_3 Point; typedef typename Gt::Point_3 Point;
@ -182,17 +189,20 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt,
CGAL_triangulation_precondition (dt.dimension()== 3); 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); Cell_handle c = dt.locate( Q, lt, li, lj, start);
if ( lt == Dt::VERTEX ) if ( lt == Dt::VERTEX )
{ {
*nn_out++= std::make_pair(c->vertex(li),Coord_type(1)); *nn_out++ = std::make_pair(c->vertex(li),Coord_type(1));
return make_triple(nn_out,norm_coeff=Coord_type(1),true); return make_triple(nn_out,norm_coeff=Coord_type(1),true);
} }
else if (dt.is_infinite(c)) 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<Cell_handle> cells; std::set<Cell_handle> cells;
typename std::set<Cell_handle>::iterator cit; typename std::set<Cell_handle>::iterator cit;
@ -201,8 +211,8 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt,
// Find the cells in conflict with Q // Find the cells in conflict with Q
dt.find_conflicts(Q, c, dt.find_conflicts(Q, c,
Emptyset_iterator(), Emptyset_iterator(),
std::inserter(cells,cells.begin())); std::inserter(cells,cells.begin()));
std::map<Vertex_handle,Coord_type> coordinate; std::map<Vertex_handle,Coord_type> coordinate;
typename std::map<Vertex_handle,Coord_type>::iterator coor_it; typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
@ -299,9 +309,9 @@ sibson_natural_neighbor_coordinates_3(const Dt& dt,
template <typename Dt, typename InputIterator> template <typename Dt, typename InputIterator>
bool is_correct_natural_neighborhood(const Dt& /*dt*/, bool is_correct_natural_neighborhood(const Dt& /*dt*/,
const typename Dt::Geom_traits::Point_3 & Q, const typename Dt::Geom_traits::Point_3& Q,
InputIterator it_begin, InputIterator it_end, InputIterator it_begin, InputIterator it_end,
const typename Dt::Geom_traits::FT & norm_coeff) const typename Dt::Geom_traits::FT& norm_coeff)
{ {
typedef typename Dt::Geom_traits Gt; typedef typename Dt::Geom_traits Gt;
typedef typename Gt::FT Coord_type; typedef typename Gt::FT Coord_type;
@ -310,18 +320,18 @@ bool is_correct_natural_neighborhood(const Dt& /*dt*/,
Coord_type sum_z(0); Coord_type sum_z(0);
InputIterator it; InputIterator it;
for(it = it_begin ; it != it_end ; ++it) for(it = it_begin ; it != it_end ; ++it)
{ {
sum_x += it->second*(it->first->point().x()); sum_x += it->second*(it->first->point().x());
sum_y += it->second*(it->first->point().y()); sum_y += it->second*(it->first->point().y());
sum_z += it->second*(it->first->point().z()); sum_z += it->second*(it->first->point().z());
} }
//!!!! to be replaced by a linear combination of points as soon //!!!! to be replaced by a linear combination of points as soon
// as it is available in the kernel. // as it is available in the kernel.
std::cout << sum_x/norm_coeff << " " std::cout << sum_x/norm_coeff << " "
<< sum_y/norm_coeff << " " << sum_y/norm_coeff << " "
<< sum_z/norm_coeff << std::endl; << sum_z/norm_coeff << std::endl;
return ((sum_x==norm_coeff*Q.x())&&(sum_y==norm_coeff*Q.y()) return ((sum_x == norm_coeff*Q.x()) && (sum_y == norm_coeff*Q.y())
&&(sum_z==norm_coeff*Q.z())); && (sum_z == norm_coeff*Q.z()));
} }
// ====================== Geometric Traits utilities ========================================= // ====================== Geometric Traits utilities =========================================

80
Interpolation/include/CGAL/predicates/predicates_for_voronoi_intersection_cartesian_2_3.h
View File

@ -31,11 +31,11 @@ namespace CGAL {
template < class RT> template < class RT>
Oriented_side Oriented_side
side_of_plane_centered_sphere_translateC3( side_of_plane_centered_sphere_translateC3(
const RT &ax, const RT &ay, const RT &az, const RT &ax, const RT &ay, const RT &az,
const RT &nx, const RT &ny, const RT &nz, const RT &nx, const RT &ny, const RT &nz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
const RT &rx, const RT &ry, const RT &rz, const RT &rx, const RT &ry, const RT &rz,
const RT &tx, const RT &ty, const RT &tz) const RT &tx, const RT &ty, const RT &tz)
{ {
RT q2 = CGAL_NTS square(qx) + CGAL_NTS square(qy) + CGAL_NTS square(qz); 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); 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); na *= RT(2.0);
Sign num = sign_of_determinant(rx, ry, rz, r2, Sign num = sign_of_determinant(rx, ry, rz, r2,
qx, qy, qz, q2, qx, qy, qz, q2,
nx, ny, nz, na, nx, ny, nz, na,
tx, ty, tz, t2); tx, ty, tz, t2);
//denumerator: //denumerator:
Sign den = sign_of_determinant(nx,ny,nz, Sign den = sign_of_determinant(nx,ny,nz,
qx,qy,qz, qx,qy,qz,
rx,ry,rz); rx,ry,rz);
CGAL_assertion(den != ZERO); CGAL_assertion(den != ZERO);
return den * num; return den * num;
@ -60,11 +60,11 @@ side_of_plane_centered_sphere_translateC3(
template < class RT> template < class RT>
Oriented_side Oriented_side
side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, 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 &nx, const RT &ny, const RT &nz,
const RT &px, const RT &py, const RT &pz, const RT &px, const RT &py, const RT &pz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
const RT &rx, const RT &ry, const RT &rz, const RT &rx, const RT &ry, const RT &rz,
const RT &tx, const RT &ty, const RT &tz) const RT &tx, const RT &ty, const RT &tz)
{ {
// resolution of the system (where c denotes the sphere's center) // 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 // - seperate computation of det and norm of the expression
return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz, return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz,
nx, ny, nz, nx, ny, nz,
qx-px, qy-py,qz-pz, qx-px, qy-py,qz-pz,
rx-px, ry-py,rz-pz, rx-px, ry-py,rz-pz,
tx-px, ty-py,tz-pz); tx-px, ty-py,tz-pz);
} }
template < class RT> template < class RT>
Oriented_side Oriented_side
side_of_plane_centered_sphere_translateC3( side_of_plane_centered_sphere_translateC3(
const RT &ax, const RT &ay, const RT &az, const RT &ax, const RT &ay, const RT &az,
const RT &nx, const RT &ny, const RT &nz, const RT &nx, const RT &ny, const RT &nz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
const RT &rx, const RT &ry, const RT &rz) const RT &rx, const RT &ry, const RT &rz)
{ {
//first choice of n_ortho: (ny+nz, -nx, -nx) //first choice of n_ortho: (ny+nz, -nx, -nx)
// if it is // if it is
@ -100,24 +100,24 @@ side_of_plane_centered_sphere_translateC3(
na *= RT(2.0); na *= RT(2.0);
Sign num = sign_of_determinant(qx, qy, qz, q2, Sign num = sign_of_determinant(qx, qy, qz, q2,
ny, -nx, RT(0), RT(0), ny, -nx, RT(0), RT(0),
nx, ny, nz, na, nx, ny, nz, na,
rx, ry, rz, r2); rx, ry, rz, r2);
//denumerator: //denumerator:
Sign den = sign_of_determinant(nx,ny,nz, Sign den = sign_of_determinant(nx,ny,nz,
ny,-nx, RT(0), ny,-nx, RT(0),
qx,qy,qz); qx,qy,qz);
if (den==ZERO) { if (den==ZERO) {
// bad choice: (ny,-nx,0) is coplanar with n,q. // bad choice: (ny,-nx,0) is coplanar with n,q.
// by precondition: q and n may not be collinear // by precondition: q and n may not be collinear
// => the cross product q*n is orthogonal to q, n and not coplanar // => the cross product q*n is orthogonal to q, n and not coplanar
num = sign_of_determinant(qx, qy, qz, q2, num = sign_of_determinant(qx, qy, qz, q2,
ny*qz-nz*qy, nz*qx-nx*qz,nx*qy-ny*qx, RT(0), ny*qz-nz*qy, nz*qx-nx*qz,nx*qy-ny*qx, RT(0),
nx, ny, nz, na, nx, ny, nz, na,
rx, ry, rz, r2); rx, ry, rz, r2);
den = sign_of_determinant(nx,ny,nz, den = sign_of_determinant(nx,ny,nz,
ny*qz-nz*qy, nz*qx - nx*qz,nx*qy-ny*qx, ny*qz-nz*qy, nz*qx - nx*qz,nx*qy-ny*qx,
qx,qy,qz); qx,qy,qz);
} }
CGAL_assertion(den != ZERO); CGAL_assertion(den != ZERO);
return den * num; return den * num;
@ -126,10 +126,10 @@ side_of_plane_centered_sphere_translateC3(
template < class RT> template < class RT>
Oriented_side Oriented_side
side_of_plane_centered_sphereC3(const RT &ax, const RT &ay, const RT &az, 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 &nx, const RT &ny, const RT &nz,
const RT &px, const RT &py, const RT &pz, const RT &px, const RT &py, const RT &pz,
const RT &qx, const RT &qy, const RT &qz, const RT &qx, const RT &qy, const RT &qz,
const RT &rx, const RT &ry, const RT &rz) const RT &rx, const RT &ry, const RT &rz)
{ {
// precondition: no two points p,q,r have the same projection // precondition: no two points p,q,r have the same projection
// <=> (p-q),(p-r), (q-r) may not be collinear to n // <=> (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 // - seperate computation of det and nom of the expression
return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz, return side_of_plane_centered_sphere_translateC3(ax-px, ay-py, az-pz,
nx, ny, nz, nx, ny, nz,
qx-px, qy-py,qz-pz, qx-px, qy-py,qz-pz,
rx-px, ry-py,rz-pz); rx-px, ry-py,rz-pz);
} }
} //namespace CGAL } //namespace CGAL

223
Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#ifndef CGAL_REGULAR_NEIGHBOR_COORDINATES_2_H #ifndef CGAL_REGULAR_NEIGHBOR_COORDINATES_2_H
@ -23,23 +19,22 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/Polygon_2.h>
#include <CGAL/iterator.h> #include <CGAL/iterator.h>
//for definition of class Project_vertex_output_iterator //for definition of class Project_vertex_output_iterator
#include <CGAL/natural_neighbor_coordinates_2.h> #include <CGAL/natural_neighbor_coordinates_2.h>
#include <list>
#include <utility>
#include <vector>
namespace CGAL { namespace CGAL {
// in this functions, the traits class is defined via the regular // In these functions, the traits class is defined via the regular triangulation.
// triangulation // See natural_neighbor_coordinates_2 for a proposal for signatures
// see natural_neighbor_coordinates_2 for a proposal for signatures // that allow to pass the traits class as argument.
// that allow to pass the traits class as argument
// The following two functions assume that
//the following two functions suppose that
// OutputIterator has value type // OutputIterator has value type
// std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT> // std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>
//!!!they are not documented!!! //!!!they are not documented!!!
@ -49,41 +44,40 @@ namespace CGAL {
template <class Rt, class OutputIterator> template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_vertex_2(const Rt& rt, regular_neighbor_coordinates_vertex_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out) OutputIterator out)
{ {
return regular_neighbor_coordinates_vertex_2(rt, p, 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 // OutputIterator has value type
// std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT> // std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator> template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_vertex_2(const Rt& rt, regular_neighbor_coordinates_vertex_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
typename Rt::Face_handle start) typename Rt::Face_handle start)
{ {
return regular_neighbor_coordinates_vertex_2(rt, p, out, 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 // OutputIterator has value type
// std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT> // std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class OutputIteratorVorVertices> template <class Rt, class OutputIterator, class OutputIteratorVorVertices>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_vertex_2(const Rt& rt, regular_neighbor_coordinates_vertex_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
OutputIteratorVorVertices vor_vertices, OutputIteratorVorVertices vor_vertices,
typename Rt::Face_handle start) typename Rt::Face_handle start)
{ {
//out: the result of the coordinate computation // out: the result of the coordinate computation
//vor_vertices: the vertices of the power cell (to avoid // vor_vertices: the vertices of the power cell (to avoid recomputation)
// recomputation)
typedef typename Rt::Geom_traits Traits; typedef typename Rt::Geom_traits Traits;
typedef typename Traits::FT Coord_type; typedef typename Traits::FT Coord_type;
@ -98,19 +92,16 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
int li; int li;
Face_handle fh = rt.locate(p, lt, li, start); Face_handle fh = rt.locate(p, lt, li, start);
//the point must lie inside the convex hull // the point must lie inside the convex hull otherwisereturn false:
// sinon return false: if(lt == Rt::OUTSIDE_AFFINE_HULL || lt == Rt::OUTSIDE_CONVEX_HULL
if(lt == Rt::OUTSIDE_AFFINE_HULL || lt == || (lt == Rt::EDGE
Rt::OUTSIDE_CONVEX_HULL && (rt.is_infinite(fh) || rt.is_infinite(fh->neighbor(li)))))
|| (lt == Rt::EDGE && (rt.is_infinite(fh)
|| rt.is_infinite(fh->neighbor(li)))))
return make_triple(out, Coord_type(1), false); return make_triple(out, Coord_type(1), false);
if (lt == Rt::VERTEX) if (lt == Rt::VERTEX)
{ {
//the point must be in conflict: //the point must be in conflict:
CGAL_precondition(rt.power_test(fh->vertex(li)->point(), p) != CGAL_precondition(rt.power_test(fh->vertex(li)->point(), p) != ON_NEGATIVE_SIDE);
ON_NEGATIVE_SIDE);
if (rt.power_test(fh->vertex(li)->point(), p) ==ON_ORIENTED_BOUNDARY) if (rt.power_test(fh->vertex(li)->point(), p) ==ON_ORIENTED_BOUNDARY)
{ {
*out++= std::make_pair(fh->vertex(li),Coord_type(1)); *out++= std::make_pair(fh->vertex(li),Coord_type(1));
@ -119,52 +110,50 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
} }
std::list<Edge> hole; std::list<Edge> hole;
std::list< Vertex_handle > hidden_vertices; std::list<Vertex_handle> hidden_vertices;
rt.get_boundary_of_conflicts_and_hidden_vertices(p, rt.get_boundary_of_conflicts_and_hidden_vertices(p,
std::back_inserter(hole), std::back_inserter(hole),
std::back_inserter std::back_inserter
(hidden_vertices), (hidden_vertices),
fh); fh);
return regular_neighbor_coordinates_vertex_2 return regular_neighbor_coordinates_vertex_2(rt, p, out, vor_vertices,
(rt, p, out, vor_vertices, hole.begin(),hole.end(), hole.begin(),hole.end(),
hidden_vertices.begin(), hidden_vertices.end()); hidden_vertices.begin(),
hidden_vertices.end());
} }
// OutputIterator has value type // OutputIterator has value type
// std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT> // std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class EdgeIterator, template <class Rt, class OutputIterator, class EdgeIterator,
class VertexIterator > class VertexIterator >
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_vertex_2(const Rt& rt, regular_neighbor_coordinates_vertex_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, EdgeIterator OutputIterator out,
hole_begin, EdgeIterator hole_end, EdgeIterator hole_begin, EdgeIterator hole_end,
VertexIterator hidden_vertices_begin, VertexIterator hidden_vertices_begin,
VertexIterator hidden_vertices_end) VertexIterator hidden_vertices_end)
{ {
return regular_neighbor_coordinates_vertex_2(rt, p, return regular_neighbor_coordinates_vertex_2(rt, p, out, Emptyset_iterator(),
out,Emptyset_iterator(), hole_begin, hole_end,
hole_begin, hole_end, hidden_vertices_begin,
hidden_vertices_begin, hidden_vertices_end);
hidden_vertices_end);
} }
// OutputIterator has value type // OutputIterator has value type
// std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT> // std::pair<Rt::Vertex_handle, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class EdgeIterator, template <class Rt, class OutputIterator, class EdgeIterator,
class VertexIterator , class OutputIteratorVorVertices > class VertexIterator , class OutputIteratorVorVertices >
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_vertex_2(const Rt& rt, regular_neighbor_coordinates_vertex_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
OutputIteratorVorVertices vor_vertices, OutputIteratorVorVertices vor_vertices,
EdgeIterator EdgeIterator hole_begin, EdgeIterator hole_end,
hole_begin, EdgeIterator hole_end, VertexIterator hidden_vertices_begin,
VertexIterator hidden_vertices_begin, VertexIterator hidden_vertices_end)
VertexIterator hidden_vertices_end)
{ {
//precondition: p must lie inside the non-empty hole //precondition: p must lie inside the non-empty hole
// (=^ inside convex hull of neighbors) // (=^ 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: //a first Voronoi vertex of the cell of p:
vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object() vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object()
(current->point(), (current->point(),
hit->first->vertex(rt.ccw(hit->second))->point(), p); hit->first->vertex(rt.ccw(hit->second))->point(), p);
*vor_vertices++= vor[0]; *vor_vertices++= vor[0];
//triangulation of the Voronoi subcell: //triangulation of the Voronoi subcell:
@ -215,17 +204,17 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
vor[1] = rt.dual(fc); vor[1] = rt.dual(fc);
// iteration over all other "old" Voronoi vertices // iteration over all other "old" Voronoi vertices
while(!fc->has_vertex(prev)) while(!fc->has_vertex(prev))
{ {
++fc; ++fc;
vor[2] = rt.dual(fc); vor[2] = rt.dual(fc);
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits()); area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
vor[1] = vor[2]; vor[1] = vor[2];
} }
//the second Voronoi vertex of the cell of p: //the second Voronoi vertex of the cell of p:
vor[2] = vor[2] =
rt.geom_traits().construct_weighted_circumcenter_2_object() rt.geom_traits().construct_weighted_circumcenter_2_object()
(prev->point(),current->point(),p); (prev->point(),current->point(),p);
*vor_vertices++= vor[2]; *vor_vertices++= vor[2];
area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits()); area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
@ -279,11 +268,10 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
template <class Rt, class OutputIterator> template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out) OutputIterator out)
{ {
return regular_neighbor_coordinates_2(rt, p, out, return regular_neighbor_coordinates_2(rt, p, out, typename Rt::Face_handle());
typename Rt::Face_handle());
} }
//OutputIterator has value type //OutputIterator has value type
@ -292,12 +280,11 @@ regular_neighbor_coordinates_2(const Rt& rt,
template <class Rt, class OutputIterator> template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
typename Rt::Face_handle start) typename Rt::Face_handle start)
{ {
return regular_neighbor_coordinates_2(rt, p, out, return regular_neighbor_coordinates_2(rt, p, out, Emptyset_iterator(), start);
Emptyset_iterator(), start);
} }
//OutputIterator has value type //OutputIterator has value type
@ -306,10 +293,10 @@ regular_neighbor_coordinates_2(const Rt& rt,
template <class Rt, class OutputIterator, class OutputIteratorVorVertices> template <class Rt, class OutputIterator, class OutputIteratorVorVertices>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
OutputIteratorVorVertices vor_vertices, OutputIteratorVorVertices vor_vertices,
typename Rt::Face_handle start) typename Rt::Face_handle start)
{ {
//out: the result of the coordinate computation //out: the result of the coordinate computation
//vor_vertices: the vertices of the power cell (to avoid //vor_vertices: the vertices of the power cell (to avoid
@ -319,45 +306,42 @@ regular_neighbor_coordinates_2(const Rt& rt,
CGAL_precondition(rt.dimension() == 2); CGAL_precondition(rt.dimension() == 2);
Triple< Project_vertex_output_iterator<OutputIterator>, Triple< Project_vertex_output_iterator<OutputIterator>,
typename Rt::Geom_traits::FT, bool > result = typename Rt::Geom_traits::FT, bool > result =
regular_neighbor_coordinates_vertex_2 regular_neighbor_coordinates_vertex_2(rt, p, op, vor_vertices, start);
(rt, p, op , vor_vertices, start);
return make_triple(result.first.base(), result.second, result.third); return make_triple(result.first.base(), result.second, result.third);
} }
//OutputIterator has value type //OutputIterator has value type
// std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class EdgeIterator, template <class Rt, class OutputIterator, class EdgeIterator,
class VertexIterator > class VertexIterator >
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, EdgeIterator OutputIterator out, EdgeIterator
hole_begin, EdgeIterator hole_end, hole_begin, EdgeIterator hole_end,
VertexIterator hidden_vertices_begin, VertexIterator hidden_vertices_begin,
VertexIterator hidden_vertices_end) VertexIterator hidden_vertices_end)
{ {
return regular_neighbor_coordinates_2(rt, p, return regular_neighbor_coordinates_2(rt, p, out, Emptyset_iterator(),
out,Emptyset_iterator(), hole_begin, hole_end,
hole_begin, hole_end, hidden_vertices_begin,
hidden_vertices_begin, hidden_vertices_end);
hidden_vertices_end);
} }
//OutputIterator has value type //OutputIterator has value type
// std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT>
template <class Rt, class OutputIterator, class EdgeIterator, template <class Rt, class OutputIterator, class EdgeIterator,
class VertexIterator , class OutputIteratorVorVertices > class VertexIterator , class OutputIteratorVorVertices >
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
const typename Rt::Weighted_point& p, const typename Rt::Weighted_point& p,
OutputIterator out, OutputIterator out,
OutputIteratorVorVertices vor_vertices, OutputIteratorVorVertices vor_vertices,
EdgeIterator hole_begin, EdgeIterator hole_end, EdgeIterator hole_begin, EdgeIterator hole_end,
VertexIterator hidden_vertices_begin, VertexIterator hidden_vertices_begin,
VertexIterator hidden_vertices_end) VertexIterator hidden_vertices_end)
{ {
//precondition: p must lie inside the non-empty hole //precondition: p must lie inside the non-empty hole
// (=^ inside convex hull of neighbors) // (=^ inside convex hull of neighbors)
@ -367,10 +351,11 @@ regular_neighbor_coordinates_2(const Rt& rt,
Project_vertex_output_iterator<OutputIterator> op(out); Project_vertex_output_iterator<OutputIterator> op(out);
Triple< Project_vertex_output_iterator<OutputIterator>, Triple< Project_vertex_output_iterator<OutputIterator>,
typename Rt::Geom_traits::FT, bool > result = typename Rt::Geom_traits::FT, bool > result =
regular_neighbor_coordinates_vertex_2 regular_neighbor_coordinates_vertex_2(rt, p, op , vor_vertices,
(rt, p, op , vor_vertices, hole_begin,hole_end, hole_begin, hole_end,
hidden_vertices_begin, hidden_vertices_end); hidden_vertices_begin,
hidden_vertices_end);
return make_triple(result.first.base(), result.second, result.third); return make_triple(result.first.base(), result.second, result.third);
} }
@ -382,8 +367,8 @@ regular_neighbor_coordinates_2(const Rt& rt,
template <class Rt, class OutputIterator> template <class Rt, class OutputIterator>
Triple< OutputIterator, typename Rt::Geom_traits::FT, bool > Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
regular_neighbor_coordinates_2(const Rt& rt, regular_neighbor_coordinates_2(const Rt& rt,
typename Rt::Vertex_handle vh, typename Rt::Vertex_handle vh,
OutputIterator out) OutputIterator out)
{ {
//this functions creates a small triangulation of the //this functions creates a small triangulation of the
// incident vertices of this vertex and computes the // incident vertices of this vertex and computes the
@ -393,8 +378,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
CGAL_precondition(rt.dimension() == 2); CGAL_precondition(rt.dimension() == 2);
Rt t2; Rt t2;
Vertex_circulator vc = rt.incident_vertices(vh), Vertex_circulator vc = rt.incident_vertices(vh), done(vc);
done(vc);
do{ do{
CGAL_assertion(!rt.is_infinite(vc)); CGAL_assertion(!rt.is_infinite(vc));
t2.insert(vc->point()); t2.insert(vc->point());
@ -404,7 +388,6 @@ regular_neighbor_coordinates_2(const Rt& rt,
return regular_neighbor_coordinates_2(t2, vh->point(), out); return regular_neighbor_coordinates_2(t2, vh->point(), out);
} }
//class providing a function object: //class providing a function object:
//OutputIterator has value type //OutputIterator has value type
// std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT> // std::pair< Rt::Geom_traits::Point_2, Rt::Geom_traits::FT>
@ -414,8 +397,8 @@ class regular_neighbor_coordinates_2_object
public: public:
Triple< OutputIterator, typename Rt::Geom_traits::FT , bool > Triple< OutputIterator, typename Rt::Geom_traits::FT , bool >
operator()(const Rt& rt, operator()(const Rt& rt,
typename Rt::Vertex_handle vh, typename Rt::Vertex_handle vh,
OutputIterator out) const OutputIterator out) const
{ {
return regular_neighbor_coordinates_2(rt, vh, out); return regular_neighbor_coordinates_2(rt, vh, out);
} }

101
Interpolation/include/CGAL/sibson_gradient_fitting.h
View File

@ -23,24 +23,26 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/Origin.h> #include <CGAL/Origin.h>
#include <CGAL/natural_neighbor_coordinates_2.h> #include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/regular_neighbor_coordinates_2.h> #include <CGAL/regular_neighbor_coordinates_2.h>
#include <iterator>
#include <utility>
namespace CGAL { namespace CGAL {
template < class ForwardIterator, class Functor, class Traits> template < class ForwardIterator, class Functor, class Traits>
typename Traits::Vector_d typename Traits::Vector_d
sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond, sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond,
const typename const typename
std::iterator_traits<ForwardIterator>:: std::iterator_traits<ForwardIterator>::
value_type::second_type& value_type::second_type& norm,
norm, const typename const typename
std::iterator_traits<ForwardIterator>::value_type std::iterator_traits<ForwardIterator>::value_type
::first_type& p, Functor function_value, ::first_type& p,
const Traits& traits) Functor function_value,
const Traits& traits)
{ {
CGAL_precondition( first!=beyond && norm!=0); CGAL_precondition( first!=beyond && norm!=0);
typedef typename Traits::Aff_transformation_d Aff_transformation; 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 CGAL_assertion(fn.second); //function value of p is valid
typename Traits::Vector_d pn = typename Traits::Vector_d pn =
traits.construct_vector_d_object()(NULL_VECTOR); traits.construct_vector_d_object()(NULL_VECTOR);
Aff_transformation scaling, m, Aff_transformation scaling, m,
Hn(traits.construct_null_matrix_d_object()()); Hn(traits.construct_null_matrix_d_object()());
for(;first!=beyond; ++first){ for(;first!=beyond; ++first){
Coord_type square_dist = traits.compute_squared_distance_d_object() Coord_type square_dist = traits.compute_squared_distance_d_object()
(first->first, p); (first->first, p);
CGAL_assertion(square_dist != 0); CGAL_assertion(square_dist != 0);
Coord_type scale = first->second/(norm*square_dist); Coord_type scale = first->second / (norm*square_dist);
typename Traits::Vector_d d= typename Traits::Vector_d d =
traits.construct_vector_d_object()(p, first->first); traits.construct_vector_d_object()(p, first->first);
//compute the vector pn: //compute the vector pn:
typename Functor::result_type f = function_value(first->first); typename Functor::result_type f = function_value(first->first);
CGAL_assertion(f.second);//function value of first->first is valid CGAL_assertion(f.second);//function value of first->first is valid
pn = pn + traits.construct_scaled_vector_d_object() pn = pn + traits.construct_scaled_vector_d_object()
(d,scale * (f.first - fn.first)); (d,scale * (f.first - fn.first));
//compute the matrix Hn: //compute the matrix Hn:
m = traits.construct_outer_product_d_object()(d); m = traits.construct_outer_product_d_object()(d);
scaling = traits.construct_scaling_matrix_d_object()(scale); 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); return Hn.inverse().transform(pn);
@ -82,10 +84,10 @@ template < class Triangul, class OutputIterator, class Functor,
class CoordFunctor, class Traits> class CoordFunctor, class Traits>
OutputIterator OutputIterator
sibson_gradient_fitting(const Triangul& tr, sibson_gradient_fitting(const Triangul& tr,
OutputIterator out, OutputIterator out,
Functor function_value, Functor function_value,
CoordFunctor compute_coordinates, CoordFunctor compute_coordinates,
const Traits& traits) const Traits& traits)
{ {
typedef typename Traits::Point_d Point; typedef typename Traits::Point_d Point;
typedef typename Traits::FT Coord_type; typedef typename Traits::FT Coord_type;
@ -93,22 +95,19 @@ sibson_gradient_fitting(const Triangul& tr,
std::vector< std::pair< Point, Coord_type > > coords; std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm; Coord_type norm;
typename Triangul::Finite_vertices_iterator typename Triangul::Finite_vertices_iterator vit = tr.finite_vertices_begin();
vit = tr.finite_vertices_begin();
for(; vit != tr.finite_vertices_end(); ++vit){ for(; vit != tr.finite_vertices_end(); ++vit){
//test if vit is a convex hull vertex: //test if vit is a convex hull vertex, otherwise do nothing
//otherwise do nothing
if (!tr.is_edge(vit, tr.infinite_vertex())) if (!tr.is_edge(vit, tr.infinite_vertex()))
{ {
norm = compute_coordinates(tr, vit, std::back_inserter(coords)).second; norm = compute_coordinates(tr, vit, std::back_inserter(coords)).second;
*out++ = std::make_pair(vit->point(), *out++ = std::make_pair(vit->point(),
sibson_gradient_fitting(coords.begin(), sibson_gradient_fitting(coords.begin(),
coords.end(), coords.end(),
norm, vit->point(), norm, vit->point(),
function_value, function_value,
traits)); traits));
coords.clear(); coords.clear();
} }
} }
return out; return out;
@ -122,33 +121,37 @@ sibson_gradient_fitting(const Triangul& tr,
template < class Dt, class OutputIterator, class Functor, class Traits> template < class Dt, class OutputIterator, class Functor, class Traits>
OutputIterator OutputIterator
sibson_gradient_fitting_nn_2(const Dt& dt, sibson_gradient_fitting_nn_2(const Dt& dt,
OutputIterator out, OutputIterator out,
Functor function_value, Functor function_value,
const Traits& traits) const Traits& traits)
{ {
typedef typename std::back_insert_iterator< std::vector< std::pair< typedef typename std::back_insert_iterator<
typename Traits::Point_d,typename Traits::FT > > > CoordInserter; std::vector<
std::pair< typename Traits::Point_d,
typename Traits::FT > > > CoordInserter;
return sibson_gradient_fitting return sibson_gradient_fitting(dt, out, function_value,
(dt, out, function_value, natural_neighbor_coordinates_2_object< Dt,
natural_neighbor_coordinates_2_object< Dt, CoordInserter >(), CoordInserter >(),
traits); traits);
} }
template < class Rt, class OutputIterator, class Functor, class Traits> template < class Rt, class OutputIterator, class Functor, class Traits>
OutputIterator OutputIterator
sibson_gradient_fitting_rn_2(const Rt& rt, sibson_gradient_fitting_rn_2(const Rt& rt,
OutputIterator out, OutputIterator out,
Functor function_value, Functor function_value,
const Traits& traits) const Traits& traits)
{ {
typedef typename std::back_insert_iterator< std::vector< std::pair< typedef typename std::back_insert_iterator<
typename Traits::Point_d,typename Traits::FT > > > CoordInserter; std::vector<
std::pair< typename Traits::Point_d,
typename Traits::FT > > > CoordInserter;
return sibson_gradient_fitting return sibson_gradient_fitting(rt, out, function_value,
(rt, out, function_value, regular_neighbor_coordinates_2_object< Rt,
regular_neighbor_coordinates_2_object< Rt, CoordInserter >(), CoordInserter >(),
traits); traits);
} }
} //namespace CGAL } //namespace CGAL

157
Interpolation/include/CGAL/surface_neighbor_coordinates_3.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
// ATTENTION : the surface is supposed to be a closed surface // ATTENTION : the surface is supposed to be a closed surface
@ -25,40 +21,43 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/Iterator_project.h> #include <CGAL/Iterator_project.h>
#include <CGAL/Voronoi_intersection_2_traits_3.h> #include <CGAL/Voronoi_intersection_2_traits_3.h>
#include <CGAL/Regular_triangulation_2.h> #include <CGAL/Regular_triangulation_2.h>
#include <CGAL/regular_neighbor_coordinates_2.h> #include <CGAL/regular_neighbor_coordinates_2.h>
#include <algorithm>
#include <functional>
#include <iterator>
#include <list>
#include <utility>
#include <vector>
namespace CGAL { namespace CGAL {
template <class OutputIterator, class InputIterator, class Kernel> template <class OutputIterator, class InputIterator, class Kernel>
inline inline
Triple< OutputIterator, typename Kernel::FT, bool > Triple< OutputIterator, typename Kernel::FT, bool >
surface_neighbor_coordinates_3(InputIterator surface_neighbor_coordinates_3(InputIterator first, InputIterator beyond,
first, InputIterator beyond, const typename Kernel::Point_3& p,
const typename Kernel::Point_3& p, const typename Kernel::Vector_3& normal,
const typename Kernel::Vector_3& normal, OutputIterator out,
OutputIterator out, const Kernel&)
const Kernel&)
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> 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 <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
Triple< OutputIterator, typename ITraits::FT, bool > Triple< OutputIterator, typename ITraits::FT, bool >
surface_neighbor_coordinates_3(InputIterator surface_neighbor_coordinates_3(InputIterator first, InputIterator beyond,
first, InputIterator beyond, const typename ITraits::Point_2& p,
const typename ITraits::Point_2& p, OutputIterator out,
OutputIterator out, const ITraits& traits)
const ITraits& traits)
{ {
//definition of the Voronoi intersection triangulation: //definition of the Voronoi intersection triangulation:
typedef Regular_triangulation_2< ITraits> I_triangulation; typedef Regular_triangulation_2<ITraits> I_triangulation;
//build Voronoi intersection triangulation: //build Voronoi intersection triangulation:
I_triangulation it(traits); I_triangulation it(traits);
@ -99,15 +98,15 @@ surface_neighbor_coordinates_3(InputIterator
template <class OutputIterator, class InputIterator, class Kernel> template <class OutputIterator, class InputIterator, class Kernel>
Quadruple< OutputIterator, typename Kernel::FT, bool, bool > Quadruple< OutputIterator, typename Kernel::FT, bool, bool >
surface_neighbor_coordinates_certified_3(InputIterator surface_neighbor_coordinates_certified_3(InputIterator
first, InputIterator beyond, first, InputIterator beyond,
const typename Kernel::Point_3& p, const typename Kernel::Point_3& p,
const typename Kernel::Vector_3& normal, const typename Kernel::Vector_3& normal,
OutputIterator out, OutputIterator out,
const Kernel& ) const Kernel& )
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> I_gt;
return surface_neighbor_coordinates_certified_3 return surface_neighbor_coordinates_certified_3(first, beyond, p, out,
(first, beyond, p, out, I_gt(p,normal)); I_gt(p,normal));
} }
//this function takes the radius of the sphere centered on p //this function takes the radius of the sphere centered on p
@ -116,34 +115,33 @@ surface_neighbor_coordinates_certified_3(InputIterator
template <class OutputIterator, class InputIterator, class Kernel> template <class OutputIterator, class InputIterator, class Kernel>
inline inline
Quadruple< OutputIterator, typename Kernel::FT, bool, bool > Quadruple< OutputIterator, typename Kernel::FT, bool, bool >
surface_neighbor_coordinates_certified_3( surface_neighbor_coordinates_certified_3(InputIterator first, InputIterator beyond,
InputIterator first, InputIterator beyond, const typename Kernel::Point_3& p,
const typename Kernel::Point_3& p, const typename Kernel::Vector_3& normal,
const typename Kernel::Vector_3& normal, const typename Kernel::FT& radius,
const typename Kernel::FT& radius, OutputIterator out, const Kernel& )
OutputIterator out, const Kernel& )
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> I_gt;
return surface_neighbor_coordinates_certified_3 return surface_neighbor_coordinates_certified_3(first, beyond, p, radius, out,
(first, beyond, p, radius, out, I_gt(p,normal)); I_gt(p,normal));
} }
// FIXME : this should probably be replaced by some kernel functor. // FIXME : this should probably be replaced by some kernel functor.
//struct necessary to sort the points by distance to p: //struct necessary to sort the points by distance to p:
//also used in surface_neighbors_3.h //also used in surface_neighbors_3.h
template <class Traits > template <class Traits>
struct closer_to_point struct closer_to_point
: public std::less<typename Traits::Point_2> : public std::less<typename Traits::Point_2>
{ {
typedef typename Traits::Point_2 Point_2; typedef typename Traits::Point_2 Point_2;
closer_to_point(const Point_2& _p, const Traits& t) 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 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: private:
Point_2 p; Point_2 p;
@ -153,34 +151,31 @@ private:
// Versions with instantiated traits class: // Versions with instantiated traits class:
template <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
Quadruple< OutputIterator, typename ITraits::FT, bool, bool > Quadruple< OutputIterator, typename ITraits::FT, bool, bool >
surface_neighbor_coordinates_certified_3(InputIterator surface_neighbor_coordinates_certified_3(InputIterator first,
first, InputIterator beyond, InputIterator beyond,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
OutputIterator out, OutputIterator out,
const ITraits& traits) const ITraits& traits)
{ {
//find the point in [first,beyond) furthest from p: //find the point in [first,beyond) furthest from p:
InputIterator furthest = std::max_element(first, beyond, InputIterator furthest = std::max_element(first, beyond,
closer_to_point<ITraits>(p, traits)); closer_to_point<ITraits>(p, traits));
return surface_neighbor_coordinates_certified_3 return surface_neighbor_coordinates_certified_3(first, beyond, p,
(first, beyond, p, traits.compute_squared_distance_2_object()(p,*furthest),
traits.compute_squared_distance_2_object()(p,*furthest), out, traits);
out, traits);
} }
//with radius(maximal distance from p to [first,beyond)) as //with radius(maximal distance from p to [first,beyond)) as
// add. parameter: // add. parameter:
template <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
Quadruple< OutputIterator, typename ITraits::FT, bool, bool > Quadruple< OutputIterator, typename ITraits::FT, bool, bool >
surface_neighbor_coordinates_certified_3(InputIterator surface_neighbor_coordinates_certified_3(InputIterator first,
first, InputIterator beyond, InputIterator beyond,
const typename const typename ITraits::Point_2& p,
ITraits::Point_2& p, const typename ITraits::FT& radius,
const typename ITraits::FT& OutputIterator out,
radius, const ITraits& traits)
OutputIterator out,
const ITraits& traits)
{ {
//definition of the Voronoi intersection triangulation: //definition of the Voronoi intersection triangulation:
typedef Regular_triangulation_2< ITraits> I_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 //collect the Voronoi vertices of the cell of p in order to
//determine the furthest distance from p to the boundary of its cell //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 //unfortunately, there is no function call without Face_handle
// "start" because this would cause type conflicts because // "start" because this would cause type conflicts because
@ -262,11 +257,10 @@ template <class Dt, class OutputIterator>
inline inline
Triple< OutputIterator, typename Dt::Geom_traits::FT, bool > Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
surface_neighbor_coordinates_3(const Dt& dt, surface_neighbor_coordinates_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& p, const typename Dt::Geom_traits::Point_3& p,
const typename Dt::Geom_traits::Vector_3& normal, const typename Dt::Geom_traits::Vector_3& normal,
OutputIterator out, OutputIterator out,
typename Dt::Cell_handle start typename Dt::Cell_handle start = typename Dt::Cell_handle())
= typename Dt::Cell_handle())
{ {
typedef Voronoi_intersection_2_traits_3<typename Dt::Geom_traits> I_gt; typedef Voronoi_intersection_2_traits_3<typename Dt::Geom_traits> I_gt;
return surface_neighbor_coordinates_3(dt, p, out, I_gt(p,normal), start); 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 <class Dt, class OutputIterator, class ITraits> template <class Dt, class OutputIterator, class ITraits>
Triple< OutputIterator, typename ITraits::FT, bool > Triple< OutputIterator, typename ITraits::FT, bool >
surface_neighbor_coordinates_3(const Dt& dt, surface_neighbor_coordinates_3(const Dt& dt,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
OutputIterator out, const ITraits& traits, OutputIterator out, const ITraits& traits,
typename Dt::Cell_handle start typename Dt::Cell_handle start = typename Dt::Cell_handle())
= typename Dt::Cell_handle())
{ {
typedef typename ITraits::FT Coord_type; typedef typename ITraits::FT Coord_type;
typedef typename ITraits::Point_2 Point_3; 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::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: //the Vertex_handle is, in fact, an iterator over vertex:
typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point; typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point;
typedef Iterator_project< typedef Iterator_project<typename std::list< Vertex_handle >::iterator,
typename std::list< Vertex_handle >::iterator, Proj_point,
Proj_point, const Point_3&,
const Point_3&, const Point_3*,
const Point_3*, std::ptrdiff_t,
std::ptrdiff_t, std::forward_iterator_tag> Point_iterator;
std::forward_iterator_tag> Point_iterator;
Locate_type lt; Locate_type lt;
int li, lj ; 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 p is located on a vertex: the only neighbor is found
if(lt == Dt::VERTEX){ if(lt == Dt::VERTEX){
*out++= std::make_pair(c->vertex(li)->point(), *out++= std::make_pair(c->vertex(li)->point(), Coord_type(1));
Coord_type(1));
return make_triple(out, Coord_type(1), true); return make_triple(out, Coord_type(1), true);
} }
//the candidate points are the points of dt in conflict with p: //the candidate points are the points of dt in conflict with p:
typename std::list< Vertex_handle > conflict_vertices; typename std::list< Vertex_handle > conflict_vertices;
dt.vertices_on_conflict_zone_boundary(p,c, 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(); for (typename std::list< Vertex_handle >::iterator it = conflict_vertices.begin();
it != conflict_vertices.end();){ it != conflict_vertices.end();){
@ -323,10 +313,9 @@ surface_neighbor_coordinates_3(const Dt& dt,
it++; it++;
} }
} }
return surface_neighbor_coordinates_3 return surface_neighbor_coordinates_3(Point_iterator(conflict_vertices.begin()),
(Point_iterator(conflict_vertices.begin()), Point_iterator(conflict_vertices.end()),
Point_iterator(conflict_vertices.end()), p, out, traits);
p, out, traits);
} }
} //namespace CGAL } //namespace CGAL

124
Interpolation/include/CGAL/surface_neighbors_3.h
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#ifndef CGAL_SURFACE_NEIGHBORS_3_H #ifndef CGAL_SURFACE_NEIGHBORS_3_H
@ -23,8 +19,6 @@
#include <CGAL/license/Interpolation.h> #include <CGAL/license/Interpolation.h>
#include <utility>
#include <CGAL/Voronoi_intersection_2_traits_3.h> #include <CGAL/Voronoi_intersection_2_traits_3.h>
#include <CGAL/Regular_triangulation_2.h> #include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Iterator_project.h> #include <CGAL/Iterator_project.h>
@ -33,6 +27,10 @@
// the function object Project_vertex_iterator_to_point // the function object Project_vertex_iterator_to_point
#include <CGAL/surface_neighbor_coordinates_3.h> #include <CGAL/surface_neighbor_coordinates_3.h>
#include <iterator>
#include <list>
#include <utility>
namespace CGAL { namespace CGAL {
//without Delaunay filtering //without Delaunay filtering
@ -40,9 +38,9 @@ template <class OutputIterator, class InputIterator, class Kernel>
inline inline
OutputIterator OutputIterator
surface_neighbors_3(InputIterator first, InputIterator beyond, surface_neighbors_3(InputIterator first, InputIterator beyond,
const typename Kernel::Point_3& p, const typename Kernel::Point_3& p,
const typename Kernel::Vector_3& normal, const typename Kernel::Vector_3& normal,
OutputIterator out, const Kernel& ) OutputIterator out, const Kernel&)
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> I_gt;
return surface_neighbors_3(first, beyond, p, out, I_gt(p,normal)); return surface_neighbors_3(first, beyond, p, out, I_gt(p,normal));
@ -51,8 +49,8 @@ surface_neighbors_3(InputIterator first, InputIterator beyond,
template <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
OutputIterator OutputIterator
surface_neighbors_3(InputIterator first, InputIterator beyond, surface_neighbors_3(InputIterator first, InputIterator beyond,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
OutputIterator out, const ITraits& traits) OutputIterator out, const ITraits& traits)
{ {
//definition of the Voronoi intersection triangulation: //definition of the Voronoi intersection triangulation:
typedef Regular_triangulation_2< ITraits> I_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); Face_handle fh = it.locate(wp, lt, li);
if(lt == I_triangulation::VERTEX){ if(lt == I_triangulation::VERTEX){
*out++ =p; *out++ = p;
return out; return out;
} }
Vertex_handle vh = it.insert(wp, fh); Vertex_handle vh = it.insert(wp, fh);
typename I_triangulation::Vertex_circulator typename I_triangulation::Vertex_circulator vc(it.incident_vertices(vh)),
vc(it.incident_vertices(vh)), done(vc);
done(vc);
do{ do{
*out++= wp2p(vc->point()); *out++= wp2p(vc->point());
CGAL_assertion(! it.is_infinite(vc)); CGAL_assertion(! it.is_infinite(vc));
} }
while(vc++!=done); while(vc++ != done);
return out; return out;
} }
@ -105,10 +102,10 @@ surface_neighbors_3(InputIterator first, InputIterator beyond,
template <class OutputIterator, class InputIterator, class Kernel> template <class OutputIterator, class InputIterator, class Kernel>
std::pair< OutputIterator, bool > std::pair< OutputIterator, bool >
surface_neighbors_certified_3(InputIterator first, surface_neighbors_certified_3(InputIterator first,
InputIterator beyond, InputIterator beyond,
const typename Kernel::Point_3& p, const typename Kernel::Point_3& p,
const typename Kernel::Vector_3& normal, const typename Kernel::Vector_3& normal,
OutputIterator out, const Kernel& ) OutputIterator out, const Kernel&)
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> I_gt;
return surface_neighbors_certified_3(first, beyond, p, out, I_gt(p,normal)); return surface_neighbors_certified_3(first, beyond, p, out, I_gt(p,normal));
@ -120,45 +117,45 @@ surface_neighbors_certified_3(InputIterator first,
template <class OutputIterator, class InputIterator, class Kernel> template <class OutputIterator, class InputIterator, class Kernel>
std::pair< OutputIterator, bool > std::pair< OutputIterator, bool >
surface_neighbors_certified_3(InputIterator first, surface_neighbors_certified_3(InputIterator first,
InputIterator beyond, InputIterator beyond,
const typename Kernel::Point_3& p, const typename Kernel::Point_3& p,
const typename Kernel::Vector_3& normal, const typename Kernel::Vector_3& normal,
const typename Kernel::FT& radius, const typename Kernel::FT& radius,
OutputIterator out, OutputIterator out,
const Kernel& /*K*/) const Kernel& /*K*/)
{ {
typedef Voronoi_intersection_2_traits_3<Kernel> I_gt; typedef Voronoi_intersection_2_traits_3<Kernel> I_gt;
return surface_neighbors_certified_3(first, beyond, p, radius, return surface_neighbors_certified_3(first, beyond, p, radius,
out, I_gt(p,normal)); out, I_gt(p,normal));
} }
// Versions with instantiated traits class: // Versions with instantiated traits class:
template <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
std::pair< OutputIterator, bool > std::pair< OutputIterator, bool >
surface_neighbors_certified_3(InputIterator first, surface_neighbors_certified_3(InputIterator first,
InputIterator beyond, InputIterator beyond,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
OutputIterator out, OutputIterator out,
const ITraits& traits) const ITraits& traits)
{ {
//find the point in [first,beyond) furthest from p: //find the point in [first,beyond) furthest from p:
InputIterator furthest = std::max_element(first, beyond, InputIterator furthest = std::max_element(first, beyond,
closer_to_point<ITraits>(p, traits)); closer_to_point<ITraits>(p, traits));
return surface_neighbors_certified_3 return surface_neighbors_certified_3
(first, beyond, p, (first, beyond, p,
traits.compute_squared_distance_2_object()(p,*furthest), traits.compute_squared_distance_2_object()(p,*furthest),
out, traits); out, traits);
} }
template <class OutputIterator, class InputIterator, class ITraits> template <class OutputIterator, class InputIterator, class ITraits>
std::pair< OutputIterator, bool > std::pair< OutputIterator, bool >
surface_neighbors_certified_3(InputIterator first, surface_neighbors_certified_3(InputIterator first,
InputIterator beyond, InputIterator beyond,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
const typename ITraits::FT& radius, const typename ITraits::FT& radius,
OutputIterator out, OutputIterator out,
const ITraits& traits) const ITraits& traits)
{ {
//definition of the Voronoi intersection triangulation: //definition of the Voronoi intersection triangulation:
typedef Regular_triangulation_2< ITraits> I_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); Face_handle fh = it.locate(wp, lt, li);
if(lt == I_triangulation::VERTEX){ if(lt == I_triangulation::VERTEX){
*out++ =p; *out++ = p;
return std::make_pair(out,true); return std::make_pair(out,true);
} }
Vertex_handle vh = it.insert(wp, fh); 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 //determine the furthest distance from p to a vertex of its cell
bool valid(false); bool valid(false);
Face_circulator fc(it.incident_faces(vh)), fdone(fc); Face_circulator fc(it.incident_faces(vh)), fdone(fc);
do do{
valid = (!it.is_infinite(fc) && valid = (!it.is_infinite(fc) &&
(4*radius > traits.compute_squared_distance_2_object() (4*radius > traits.compute_squared_distance_2_object()
(p, it.dual(fc)))); (p, it.dual(fc))));
while(!valid && ++fc!=fdone); }while(!valid && ++fc!=fdone);
//get the neighbor points: //get the neighbor points:
Vertex_circulator vc(it.incident_vertices(vh)), vdone(vc); 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); return std::make_pair(out, valid);
} }
//using Delaunay triangulation for candidate point filtering: //using Delaunay triangulation for candidate point filtering:
// => no certification is necessary // => no certification is necessary
template <class Dt, class OutputIterator> template <class Dt, class OutputIterator>
inline inline
OutputIterator OutputIterator
surface_neighbors_3(const Dt& dt, surface_neighbors_3(const Dt& dt,
const typename Dt::Geom_traits::Point_3& p, const typename Dt::Geom_traits::Point_3& p,
const typename Dt::Geom_traits::Vector_3& normal, const typename Dt::Geom_traits::Vector_3& normal,
OutputIterator out, OutputIterator out,
typename Dt::Cell_handle start =typename Dt::Cell_handle()) typename Dt::Cell_handle start =typename Dt::Cell_handle())
{ {
typedef Voronoi_intersection_2_traits_3<typename Dt::Geom_traits> I_gt; typedef Voronoi_intersection_2_traits_3<typename Dt::Geom_traits> I_gt;
return surface_neighbors_3(dt, p, out, I_gt(p,normal),start); return surface_neighbors_3(dt, p, out, I_gt(p,normal),start);
@ -231,10 +227,9 @@ surface_neighbors_3(const Dt& dt,
template <class Dt, class OutputIterator, class ITraits> template <class Dt, class OutputIterator, class ITraits>
OutputIterator OutputIterator
surface_neighbors_3(const Dt& dt, surface_neighbors_3(const Dt& dt,
const typename ITraits::Point_2& p, const typename ITraits::Point_2& p,
OutputIterator out, const ITraits& traits, OutputIterator out, const ITraits& traits,
typename Dt::Cell_handle start typename Dt::Cell_handle start = typename Dt::Cell_handle())
= typename Dt::Cell_handle())
{ {
typedef typename ITraits::Point_2 Point_3; typedef typename ITraits::Point_2 Point_3;
@ -243,14 +238,13 @@ surface_neighbors_3(const Dt& dt,
typedef typename Dt::Locate_type Locate_type; typedef typename Dt::Locate_type Locate_type;
//the Vertex_handle is, in fact, an iterator over vertex: //the Vertex_handle is, in fact, an iterator over vertex:
typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point; typedef Project_vertex_iterator_to_point< Vertex_handle> Proj_point;
typedef Iterator_project< typedef Iterator_project<typename std::list< Vertex_handle >::iterator,
typename std::list< Vertex_handle >::iterator, Proj_point,
Proj_point, const Point_3&,
const Point_3&, const Point_3*,
const Point_3*, std::ptrdiff_t,
std::ptrdiff_t, std::forward_iterator_tag> Point_iterator;
std::forward_iterator_tag> Point_iterator;
Locate_type lt; Locate_type lt;
int li, lj ; int li, lj ;
@ -265,7 +259,7 @@ surface_neighbors_3(const Dt& dt,
//otherwise get vertices in conflict //otherwise get vertices in conflict
typename std::list< Vertex_handle > conflict_vertices; typename std::list< Vertex_handle > conflict_vertices;
dt.vertices_on_conflict_zone_boundary(p,c, 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(); for (typename std::list< Vertex_handle >::iterator it = conflict_vertices.begin();
it != conflict_vertices.end();){ it != conflict_vertices.end();){
@ -278,8 +272,8 @@ surface_neighbors_3(const Dt& dt,
} }
} }
return surface_neighbors_3(Point_iterator(conflict_vertices.begin()), return surface_neighbors_3(Point_iterator(conflict_vertices.begin()),
Point_iterator(conflict_vertices.end()), Point_iterator(conflict_vertices.end()),
p, out, traits); p, out, traits);
} }
} //namespace CGAL } //namespace CGAL

8
Interpolation/test/Interpolation/include/CGAL/_test_surface_neighbors_3.cpp
View File

@ -106,10 +106,10 @@ bool compare_neighbors(ForwardIteratorCoord first_coord,
template < class Triangul, class OutputIterator> template < class Triangul, class OutputIterator>
OutputIterator OutputIterator
test_neighbors(const Triangul& T, const typename test_neighbors(const Triangul& T, const typename
Triangul::Geom_traits::Point_3& p, Triangul::Geom_traits::Point_3& p,
const typename Triangul::Geom_traits::Vector_3& n, const typename Triangul::Geom_traits::Vector_3& n,
const int& version, const int& version,
OutputIterator out) OutputIterator out)
{ {
typedef CGAL::Voronoi_intersection_2_traits_3<typename typedef CGAL::Voronoi_intersection_2_traits_3<typename
Triangul::Geom_traits> I_traits; Triangul::Geom_traits> I_traits;

18
Interpolation/test/Interpolation/test_interpolation_functions_2.cpp
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#include <CGAL/Exact_predicates_exact_constructions_kernel.h> #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
@ -24,6 +20,8 @@
#include <CGAL/_test_interpolation_functions_2.cpp> #include <CGAL/_test_interpolation_functions_2.cpp>
#include <iostream>
typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Dt; typedef CGAL::Delaunay_triangulation_2<K> Dt;
@ -33,16 +31,16 @@ typedef CGAL::Delaunay_triangulation_2<K2> Dt2;
int main() int main()
{ {
std::cout << "Testing interpolation functions with 2D NN neighbors " std::cout << "Testing interpolation functions with 2D NN neighbors "
<< std::endl; << std::endl;
std::cout << " using Exact_predicates_exact_constructions_kernel: " std::cout << " using Exact_predicates_exact_constructions_kernel: "
<< std::endl ; << std::endl ;
_test_interpolation_functions_2_delaunay( Dt(), K::FT(1e-10)); _test_interpolation_functions_2_delaunay(Dt(), K::FT(1e-10));
std::cout << "Testing interpolation functions with 2D NN neighbors " std::cout << "Testing interpolation functions with 2D NN neighbors "
<< std::endl; << std::endl;
std::cout << " using Exact_predicates_inexact_constructions_kernel: " std::cout << " using Exact_predicates_inexact_constructions_kernel: "
<< std::endl ; << std::endl ;
_test_interpolation_functions_2_delaunay( Dt2(), K2::FT(1e-10)); _test_interpolation_functions_2_delaunay(Dt2(), K2::FT(1e-10));
std::cout << "test_interpolation_functions_2 is finished" << std::endl; std::cout << "test_interpolation_functions_2 is finished" << std::endl;
return 0; return 0;

6
Interpolation/test/Interpolation/test_natural_neighbors_2.cpp
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Naceur MESKINI. // Author(s) : Naceur MESKINI.
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@ -23,11 +19,11 @@
#include <CGAL/_test_natural_neighbors_2.cpp> #include <CGAL/_test_natural_neighbors_2.cpp>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Dt; typedef CGAL::Delaunay_triangulation_2<K> Dt;
int main() int main()
{ {
std::cout << "Testing NN_neighbors_2 " << std::endl; std::cout << "Testing NN_neighbors_2 " << std::endl;

5
Interpolation/test/Interpolation/test_regular_neighbors_2.cpp
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#include <CGAL/basic.h> #include <CGAL/basic.h>
@ -26,6 +22,7 @@
#include <CGAL/_test_regular_neighbors_2.cpp> #include <CGAL/_test_regular_neighbors_2.cpp>
#include <iostream>
typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Exact_predicates_exact_constructions_kernel K;

41
Interpolation/test/Interpolation/test_surface_neighbor_3.cpp
View File

@ -12,10 +12,6 @@
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
// //
// $URL$
// $Id$
//
//
// Author(s) : Julia Floetotto // Author(s) : Julia Floetotto
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@ -25,6 +21,7 @@
#include <CGAL/_test_surface_neighbors_3.cpp> #include <CGAL/_test_surface_neighbors_3.cpp>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_3<K> Dt; typedef CGAL::Delaunay_triangulation_3<K> Dt;
@ -35,7 +32,6 @@ typedef CGAL::Delaunay_triangulation_3<K2> Dt2;
// Fast_location with exact pred exact const. kernel: // Fast_location with exact pred exact const. kernel:
typedef CGAL::Delaunay_triangulation_3<K2, CGAL::Fast_location> Dh; typedef CGAL::Delaunay_triangulation_3<K2, CGAL::Fast_location> Dh;
// Aff_transformation: // Aff_transformation:
typedef CGAL::Aff_transformation_3<K2> Transformation; typedef CGAL::Aff_transformation_3<K2> Transformation;
@ -46,17 +42,17 @@ int main()
_test_surface_neighbors_3_sphere( Dt() ); _test_surface_neighbors_3_sphere( Dt() );
std::cout << " done." << std::endl << std::endl; std::cout << " done." << std::endl << std::endl;
std::cout << "Using Exact_predicates_exact_constructions_kernel: " std::cout << "Using Exact_predicates_exact_constructions_kernel: "
<< std::endl; << std::endl;
//AXIS ALIGNED CUBE //AXIS ALIGNED CUBE
Transformation identity(CGAL::IDENTITY); Transformation identity(CGAL::IDENTITY);
std::cout << "Testing surface_neighbors_3 on a cube with axis : " std::cout << "Testing surface_neighbors_3 on a cube with axis : "
<< identity(K2::Vector_3(0,0,1)) << ", " << identity(K2::Vector_3(0,0,1)) << ", "
<< identity(K2::Vector_3(0,1,0))<< ", " << identity(K2::Vector_3(0,1,0))<< ", "
<< identity(K2::Vector_3(1,0,0)) << identity(K2::Vector_3(1,0,0))
<< std::endl; << std::endl;
std::cout << " with grid sample points"; std::cout << " with grid sample points";
_test_surface_neighbors_3_cube(Dt2(),identity, 75, K2::FT(1e-29)); _test_surface_neighbors_3_cube(Dt2(),identity, 75, K2::FT(1e-29));
std::cout << " done." << std::endl; std::cout << " done." << std::endl;
@ -67,22 +63,23 @@ int main()
//ROTATED CUBE //ROTATED CUBE
Transformation rotate(K2::RT(1),K2::RT(0),K2::RT(0),K2::RT(0), 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.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.42261826),K2::RT(0.9063),K2::RT(0));
std::cout << "Testing surface_neighbors_3 on a ROTATED cube "<< std::endl; std::cout << "Testing surface_neighbors_3 on a ROTATED cube "<< std::endl;
std::cout << " with grid sample points"; std::cout << " with grid sample points";
_test_surface_neighbors_3_cube(Dh(),rotate, 75, K2::FT(1e-2), true); _test_surface_neighbors_3_cube(Dh(),rotate, 75, K2::FT(1e-2), true);
std::cout << " done." << std::endl << std::endl; std::cout << " done." << std::endl << std::endl;
// //undersampled rotated cube // //undersampled rotated cube
// Transformation rotate3(K2::RT(0.1),K2::RT(0.4),K2::RT(0.6),K2::RT(0), // 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.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)); // 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::cout << "Testing surface_neighbors_3 on an undersampled ROTATED cube "
// << std::endl; // << std::endl;
// std::cout << " with grid sample points"; // std::cout << " with grid sample points";
// _test_surface_neighbors_3_cube(Dh(),rotate3,75, K2::FT(9), true); // _test_surface_neighbors_3_cube(Dh(),rotate3,75, K2::FT(9), true);
// std::cout << " done." << std::endl << std::endl; // std::cout << " done." << std::endl << std::endl;
return 0; return 0;
} }