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

- first addition of example and demo files for interpolation and 2D coordinate computation.

This commit is contained in:
Julia Flötotto 2003-11-28 10:37:33 +00:00
parent 059033ddb8
commit 1a0fb4f0bd
10 changed files with 960 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Packages/Interpolation/demo/Interpolation/README Normal file
View File

@ -0,0 +1,30 @@
Some demos use Geomview
[see the chapter Geomview in the cgal manual - support library:
Geomview 1.8.1 is required. The geomview command must be in the user's $PATH,
otherwise the program will not be able to execute.]
------- demo_interpolation_2 -------------------------------------------------
using Geomview
This demo program dumps the plot of the different interpolation functions
using a very simple data set
shows: the differentiablity (or not) at the data points
the closeness to the gradient at the data points
1) Construction of a Delaunay triangulation of the data points.
2) Interpolation on a grid using a user chosen interpolation function.
3) Dumps the plot of the interpolation functions (.off file)
4) Displays the data points in geomview.
--------------------------------------------------------------
------- demo_numerical_2 -----------------------------------------------
This demo program computes the interpolation of a quadratic function
known on a small set of random sample points. The gradients of the function are given.
- verifies the barycentric property
- computes the maximum and mean error of the different interpolation function
compared to the exact function value.
--------------------------------------------------------------

348
Packages/Interpolation/demo/Interpolation/demo_interpolation_2.C Normal file
View File

@ -0,0 +1,348 @@
// ============================================================================
//
// Copyright (c) 1998-1999 The CGAL Consortium
//
// This software and related documentation is part of an INTERNAL release
// of the Computational Geometry Algorithms Library (CGAL). It is not
// intended for general use.
//
// ----------------------------------------------------------------------------
//
// release :
// release_date :
//
// file : demo_2/
// revision : $Revision$
// author(s) : Julia Floetotto
// coordinator : INRIA Sophia Antipolis (Mariette Yvinec)
//
// ============================================================================
// Geomview doesn't work on M$ at the moment, so we don't compile this
// file.
//**********************
//demo 2D Interpolation over the plane - using 2D natural neighbors
//**********************
#if defined(__BORLANDC__) || defined(_MSC_VER)
#include <iostream>
int main()
{
std::cerr << "Geomview doesn't work on Windows, so this demo doesn't work"
<< std::endl;
return 0;
}
#else
#include <CGAL/basic.h>
#include <fstream>
#include <iostream>
#include <utility>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
//#include <CGAL/double.h>
//#include <CGAL/Random.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/copy_n.h>
#include <CGAL/IO/Geomview_stream.h>
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Delaunay_triangulation_2<K> Dt;
typedef CGAL::Interpolation_traits_2<K> ITraits;
typedef K::FT Coord_type;
typedef K::Point_2 Point_2;
typedef K::Vector_2 Vector_2;
typedef K::Point_3 Point_3;
typedef K::Vector_3 Vector_3;
typedef K::Segment_3 Segment_3;
typedef std::map<Point_2, Coord_type, K::Less_xy_2> Point_value_map ;
typedef std::map<Point_2, Vector_2, K::Less_xy_2 > Point_vector_map;
typedef std::vector<Point_3> Point_vector_3;
typedef std::vector<Point_2> Point_vector_2;
//////////////////////
// VISU GEOMVIEW
//////////////////////
template<class Point_vector>
void visu_points(CGAL::Geomview_stream & os, const Point_vector & points)
{
int n = points.size();
for(int i=0; i<n ; i++)
os << points[i];
}
////////////////////////////////
template< class Map >
struct DataAccess : public std::unary_function< typename Map::key_type,
typename Map::mapped_type> {
typedef typename Map::mapped_type Data_type;
typedef typename Map::key_type Point;
//CONSTRUCTOR:
DataAccess< Map >(const Map& m): map(m){};
//Functor
Data_type operator()(const Point& p) {
typename Map::const_iterator mit = map.find(p);
if(mit!= map.end())
return mit->second;
return Data_type();
};
const Map& map;
};
//////////////////////////////////////////////////////////////////////////
////POINT GENERATION:
void generate_grid_points(Point_vector_2& points, int m, float h)
{
//int n = (m+1)*(m+1);
int n = m*m;
points.clear();
points.reserve(n);
std::cout <<" generate " <<n << " grid points in square of height "
<< h <<std::endl;
// Create n points from a 16 x 16 grid. Note that the double
// arithmetic _is_ sufficient to produce exact integer grid points.
// The distance between neighbors is 34 pixel = 510 / 15.
CGAL::points_on_square_grid_2((double) h, n,
std::back_inserter(points),
CGAL::Creator_uniform_2<double,Point_2>());
}
template < class Point_vector>
bool
dump_off_file_quadrilateral(const Point_vector& points, int m)
{
//open file stream
char fname[11];
CGAL_CLIB_STD::strcpy(fname, "result.off");
//dump reconstruction file:
std::ofstream os(fname,std::ios::out);
if (! os) {
std::cout <<"cannot open output file: result.off. "<< std::endl;
return false;
}
int n= points.size();
std::cout << n << " points " << (m-1)*(m-1) <<std::endl;
//first line: OFF
// 2. line : number of vertices and quad_indices:
os << "OFF" << std::endl;
os << n << " " <<(m-1)*(m-1) << " 0" <<
std::endl;
os << std::endl;
for(int i=0; i<n ; i++)
os << points[i]<< std::endl;
//indices
for(int i=0; i< m-1;i++)
for(int j=0; j< m-1; j++)
os << "4 " << i*m + j << " "
<< j+1 + i*m
<< " " <<(j+1) + (i+1)*m << " " <<j + (i+1)*m
<< std::endl;
return true;
}
//////////////////////
int main(int argc, char* argv[])
{
//parameters:
int m = 78;
double h = 10;
double g = 2;
Coord_type w(4);
Point_vector_2 sample;
//3D function graph: 2D points + function value in z-direction:
Point_vector_3 sample_3;
Dt T;
Point_value_map values;
Point_vector_map gradients;
sample.push_back( Point_2(h,h));
sample.push_back( Point_2( -h,h));
sample.push_back( Point_2( h,-h ));
sample.push_back( Point_2( -h,-h));
sample.push_back( Point_2( 0.0, h/2));
sample.push_back( Point_2( 0.0,-h/2));
sample.push_back( Point_2( h/2, 0.0 ));
sample.push_back( Point_2( -h/2, 0.0));
sample.push_back( Point_2(0.0,0.0));
int s = sample.size();
for(int j=0; j<s ; j++){
T.insert(sample[j]);
values.insert(std::make_pair(sample[j], Coord_type(0)));
gradients.insert(std::make_pair(sample[j], CGAL::NULL_VECTOR));
sample_3.push_back( Point_3(sample[j].x(),sample[j].y(),0));
}
Point_2 p = Point_2(h/3,h/3);
T.insert(p);
values.insert(std::make_pair(p, w));
gradients.insert(std::make_pair(p, Vector_2(-g,-g) ));
p = Point_2(-h/3,-h/3);
T.insert(p);
values.insert(std::make_pair(p,w));
gradients.insert(std::make_pair(p, Vector_2(g,g) ));
p = Point_2(-h/3,h/3);
T.insert(p);
values.insert(std::make_pair(p, w));
gradients.insert(std::make_pair(p, Vector_2(g,-g) ));
p = Point_2(h/3,-h/3);
T.insert(p);
values.insert(std::make_pair(p, w));
gradients.insert(std::make_pair(p, Vector_2(-g,g) ));
sample_3.push_back( Point_3(h/3,h/3,w));
sample_3.push_back( Point_3(-h/3,-h/3,w));
sample_3.push_back( Point_3(-h/3,h/3,w));
sample_3.push_back( Point_3(h/3,-h/3,w));
//Interpolation grid:
Point_vector_2 points;
generate_grid_points(points, m, 0.999* h);
Point_vector_3 points_3;
int method;
std::cout << "Enter the choice of the interpolation function: "<<std::endl;
std::cout << " 0 -- linear interpolation" <<std::endl;
std::cout << " 1 -- simple quadratic interpolant" <<std::endl;
std::cout << " 2 -- Sibson's C1 interpolant "<<std::endl;
std::cout << " 3 -- Sibson's C1 interpolant without square-root"<<std::endl;
std::cout << " 4 -- Farin's C1 interpolant" << std::endl;
std::cin >> method;
//INTERPOLATION:
Coord_type exact_value, res, norm;
std::vector< std::pair< Point_2, Coord_type > > coords;
int n = points.size();
ITraits traits;
std::cout << "Interpolation at "<<n <<" grid points " << std::endl;
for(int i=0;i<n;i++){
norm =
CGAL::natural_neighbor_coordinates_2(T, points[i],
std::back_inserter(coords)).second;
assert(norm>0);
switch(method){
case 0:
res = CGAL::linear_interpolation(coords.begin(),coords.end(),norm,
DataAccess<Point_value_map>(values));
break;
case 1:
res = CGAL::quadratic_interpolation(coords.begin(),coords.end(),
norm, points[i],
DataAccess< Point_value_map>
(values),
DataAccess< Point_vector_map>
(gradients),traits); break;
case 2:
res = CGAL::sibson_c1_interpolation(coords.begin(),coords.end(),
norm, points[i],
DataAccess<Point_value_map>
(values),
DataAccess<Point_vector_map>
(gradients), traits); break;
case 3:
res = CGAL::sibson_c1_interpolation_square(coords.begin(),coords.end(),
norm, points[i],
DataAccess<Point_value_map>
(values),
DataAccess<Point_vector_map>
(gradients), traits); break;
case 4:
res = CGAL::farin_c1_interpolation(coords.begin(),coords.end(),
norm, points[i],
DataAccess<Point_value_map>
(values),
DataAccess<Point_vector_map>
(gradients), traits); break;
default: std::cout <<"No valid choice of interpolant." <<
std::endl; break;
}
points_3.push_back(Point_3(points[i].x(), points[i].y(),res));
coords.clear();
}
/************** end of Coordinate computation **************/
dump_off_file_quadrilateral(points_3, m);
char ch;
//viewer
CGAL::Geomview_stream gv(CGAL::Bbox_3(0,0,0, 2, 2, 2));
gv.set_bg_color(CGAL::Color(0, 200, 200));
gv.clear();
gv.set_line_width(2);
std::cout << "The data points are displayed in blue in the geomview"
<< " application." << std::endl;
gv << CGAL::BLUE;
visu_points(gv,sample_3);
//show the gradients
if(method>0){
std::cout << "The function gradients are displayed by red lines "
<<" in the geomview application." << std::endl;
gv <<CGAL::RED;
gv << Segment_3(Point_3(h/3,h/3,w),Point_3(h/3,h/3,w)+ Vector_3(-g,-g,0));
gv << Segment_3(Point_3(-h/3,h/3,w),Point_3(-h/3,h/3,w)+Vector_3(g,-g,0));
gv << Segment_3(Point_3(-h/3,-h/3,w),Point_3(-h/3,-h/3,w)+Vector_3(g,g,0));
gv << Segment_3(Point_3(h/3,-h/3,w),Point_3(h/3,-h/3,w)+Vector_3(-g,g,0));
}
std::cout << "You should open the file 'result.off' with geomview."
<< std::endl;
std::cout << "'result.off' contains the graph of the interpolation"
<< " function over a grid using ";
switch(method){
case 0: std::cout << " linear interpolation"; break;
case 1: std::cout << " simple quadratic interpolant"; break;
case 2: std::cout << " Sibson's C1 interpolant";break;
case 3: std::cout << " Sibson's C1 interpolant without square-root"; break;
case 4: std::cout << " Farin's C1 interpolant"; break;
}
std::cout << std::endl;
std::cout << "Enter any character to quit." << std::endl;
std::cin >> ch;
return 1;
}
#endif // if defined(__BORLANDC__) || defined(_MSC_VER)

239
Packages/Interpolation/demo/Interpolation/demo_numerical_2.C Normal file
View File

@ -0,0 +1,239 @@
//file: demo/Interpolation/demo_numerical_2.C
// compares the result of several interpolation methods
// Author: Julia
#include <CGAL/basic.h>
#include <cstdio>
#include <utility>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Random.h>
#include <CGAL/squared_distance_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/copy_n.h>
#include <CGAL/Origin.h>
struct K : CGAL::Exact_predicates_exact_constructions_kernel {};
typedef K::FT Coord_type;
typedef K::Vector_2 Vector;
typedef K::Point_2 Point;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef std::vector< std::pair<Point, Coord_type> > Coordinate_vector;
typedef std::map<Point, Coord_type, K::Less_xy_2> Point_value_map;
typedef std::map<Point, Vector, K::Less_xy_2 > Point_vector_map;
////////////////////////////////
// functor class that provides access to the function values/ function
// gradients given a data point (using "map"):
template< class Map >
struct DataAccess : public std::unary_function< typename Map::key_type,
typename Map::mapped_type> {
typedef typename Map::mapped_type Data_type;
typedef typename Map::key_type Point;
//CONSTRUCTOR:
DataAccess< Map >(const Map& m): map(m){};
//Functor
Data_type operator()(const Point& p) {
typename Map::const_iterator mit = map.find(p);
if(mit!= map.end())
return mit->second;
return Data_type();
};
const Map& map;
};
//////////////test function:///////////////////////////////
template < class Iterator >
void test_barycenter(Iterator first, Iterator beyond,
typename std::iterator_traits<Iterator>::value_type::
second_type norm,
const typename std::iterator_traits<Iterator>::
value_type::first_type& p)
{
typedef typename
std::iterator_traits<Iterator>::value_type::first_type Point_2;
Point_2 b(0,0);
for(; first !=beyond; first++)
b = b+ (first->second/norm) * (first->first - CGAL::ORIGIN);
std::cout <<"Distance between test point and barycenter " <<
CGAL::squared_distance(p, b ) << " egality: "<< (p==b)<< std::endl;
}
////////////////////////////////
int main()
{
//number of sample points:
int n = 24;
//number of interpolation points:
int m = 20;
std::vector<Point> points;
points.reserve(n+m);
//put four bounding box points:
points.push_back(Point(-3,-3));
points.push_back(Point(3,-3));
points.push_back(Point(-3,3));
points.push_back(Point(3,3));
// Create n+m-4 points within a disc of radius 2
double r_d = 3;
CGAL::Random_points_in_disc_2<Point> g(r_d );
CGAL::copy_n( g, n+m, std::back_inserter(points));
Delaunay_triangulation T;
Point_value_map values;
Point_vector_map gradients;
//parameters for quadratic function:
Coord_type alpha = Coord_type(1.0),
beta1 = Coord_type(2.0),
beta2 = Coord_type(1.0),
gamma1 = Coord_type(0.2),
gamma2 = Coord_type(0.7),
gamma3 = Coord_type(0.4),
gamma4 = Coord_type(0.6);
for(int j=0; j<n ; j++){
T.insert(points[j]);
//determine function value/gradient:
Coord_type x(points[j].x());
Coord_type y(points[j].y());
Coord_type value = alpha + beta1*x + beta2*y + gamma1*(x*x) +
gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
Vector gradient(beta1+ (gamma2+ gamma3)*y + Coord_type(2)*(gamma1*x),
beta2+ (gamma2+ gamma3)*x + Coord_type(2)*(gamma4*y));
values.insert(std::make_pair(points[j], value));
gradients.insert(std::make_pair(points[j], gradient));
}
//variables for statistics:
Coord_type res, error, l_total = Coord_type(0),
q_total(l_total), f_total(l_total), s_total(l_total),
l_max(l_total),
q_max(l_total), f_max(l_total), s_max(l_total),
total_value(l_total);
int failure(0);
//interpolation + error statistics
for(int i=n;i<n+m;i++){
Coord_type x(points[i].x());
Coord_type y(points[i].y());
Coord_type exact_value = alpha + beta1*x + beta2*y + gamma1*(x*x) +
gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
total_value += exact_value;
//Coordinate_vector:
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm =
CGAL::natural_neighbor_coordinates_2(T, points[i],
std::back_inserter(coords)).second;
if(norm>0){
//test barycenter
test_barycenter( coords.begin(), coords.end(),norm, points[i]);
//linear interpolant:
res = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm,
DataAccess<Point_value_map>(values));
error = CGAL_NTS abs(res - exact_value);
l_total += error;
if (error > l_max) l_max = error;
//Farin interpolant:
res = CGAL::farin_c1_interpolation(coords.begin(),
coords.end(), norm,points[i],
DataAccess<Point_value_map>(values),
DataAccess<Point_vector_map>
(gradients),
Traits());
error = CGAL_NTS abs(res - exact_value);
f_total += error;
if (error > f_max) f_max = error;
//quadratic interpolant:
res = CGAL::quadratic_interpolation(coords.begin(), coords.end(),
norm,points[i],
DataAccess<Point_value_map>
(values),
DataAccess<Point_vector_map>
(gradients),
Traits());
error = CGAL_NTS abs(res - exact_value);
q_total += error;
if (error > q_max) q_max = error;
//Sibson interpolant: version without sqrt:
res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(), norm,
points[i],
DataAccess<Point_value_map>
(values),
DataAccess<Point_vector_map>
(gradients),
Traits());
//error statistics
error = CGAL_NTS abs(res - exact_value);
s_total += error;
if (error > s_max) s_max = error;
}else failure++;
}
/************** end of Interpolation: dump statistics **************/
std::cout << "Result: -----------------------------------" << std::endl;
std::cout << "Interpolation of '" << alpha <<" + "
<< beta1<<" x + "
<< beta2 << " y + " << gamma1 <<" x^2 + " << gamma2+ gamma3
<<" xy + " << gamma4 << " y^2'" << std::endl;
std::cout << "Knowing " << m << " sample points. Interpolation on "
<< n <<" random points. "<< std::endl;
std::cout <<"Average function value "
<< (1.0/n) * CGAL::to_double(total_value)
<< ", nb of failures "<< failure << std::endl;
std::cout << "linear interpolant mean error "
<< CGAL::to_double(l_total)/n << " max "
<< CGAL::to_double(l_max) <<std::endl;
std::cout << "quadratic interpolant mean error "
<< CGAL::to_double(q_total)/n << " max "
<< CGAL::to_double(q_max) << std::endl;
std::cout << "Farin interpolant mean error "
<< CGAL::to_double(f_total)/n << " max "
<< CGAL::to_double(f_max) << std::endl;
std::cout << "Sibson interpolant mean error "
<< CGAL::to_double(s_total)/n << " max "
<< CGAL::to_double(s_max) << std::endl;
return 1;
}
//end of file

56
Packages/Interpolation/demo/Interpolation/makefile Normal file
View File

@ -0,0 +1,56 @@
# Created by the script create_makefile
# This is the makefile for compiling a CGAL application.
#---------------------------------------------------------------------#
# include platform specific settings
#---------------------------------------------------------------------#
# Choose the right include file from the <cgalroot>/make directory.
# CGAL_MAKEFILE = ENTER_YOUR_INCLUDE_MAKEFILE_HERE
include $(CGAL_MAKEFILE)
#---------------------------------------------------------------------#
# compiler flags
#---------------------------------------------------------------------#
CXXFLAGS = \
$(CGAL_CXXFLAGS) \
$(LONG_NAME_PROBLEM_CXXFLAGS) \
$(DEBUG_OPT)
#---------------------------------------------------------------------#
# linker flags
#---------------------------------------------------------------------#
LIBPATH = \
$(CGAL_LIBPATH)
LDFLAGS = \
$(LONG_NAME_PROBLEM_LDFLAGS) \
$(CGAL_GEOWIN_LDFLAGS)
#---------------------------------------------------------------------#
# target entries
#---------------------------------------------------------------------#
all: \
demo_interpolation_2$(EXE_EXT) \
demo_numerical_2$(EXE_EXT)
demo_interpolation_2$(EXE_EXT): demo_interpolation_2$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_interpolation_2 demo_interpolation_2$(OBJ_EXT) $(LDFLAGS)
demo_numerical_2$(EXE_EXT): demo_numerical_2$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_numerical_2 demo_numerical_2$(OBJ_EXT) $(LDFLAGS)
clean: \
demo_interpolation_2.clean \
demo_numerical_2.clean
#---------------------------------------------------------------------#
# suffix rules
#---------------------------------------------------------------------#
.C$(OBJ_EXT):
$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<

59
Packages/Interpolation/demo/Interpolation/makefile_local Normal file
View File

@ -0,0 +1,59 @@
# Created by the script create_makefile
# This is the makefile for compiling a CGAL application.
#---------------------------------------------------------------------#
# include platform specific settings
#---------------------------------------------------------------------#
# Choose the right include file from the <cgalroot>/make directory.
# CGAL_MAKEFILE = ENTER_YOUR_INCLUDE_MAKEFILE_HERE
CGAL_MAKEFILE = /0/prisme_util/CGAL/CGAL-I/make/makefile_i686_Linux-2.2.18_g++-3.3.0
include $(CGAL_MAKEFILE)
#---------------------------------------------------------------------#
# compiler flags
#---------------------------------------------------------------------#
CXXFLAGS = \
-I ../../include \
$(CGAL_CXXFLAGS) \
$(LONG_NAME_PROBLEM_CXXFLAGS) \
$(DEBUG_OPT)
#---------------------------------------------------------------------#
# linker flags
#---------------------------------------------------------------------#
LIBPATH = \
$(CGAL_LIBPATH)
LDFLAGS = \
$(LONG_NAME_PROBLEM_LDFLAGS) \
$(CGAL_GEOWIN_LDFLAGS)
#---------------------------------------------------------------------#
# target entries
#---------------------------------------------------------------------#
all: \
demo_interpolation_2$(EXE_EXT) \
demo_numerical_2$(EXE_EXT)
demo_interpolation_2$(EXE_EXT): demo_interpolation_2$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_interpolation_2 demo_interpolation_2$(OBJ_EXT) $(LDFLAGS)
demo_numerical_2$(EXE_EXT): demo_numerical_2$(OBJ_EXT)
$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_numerical_2 demo_numerical_2$(OBJ_EXT) $(LDFLAGS)
clean: \
demo_interpolation_2.clean \
demo_numerical_2.clean
#---------------------------------------------------------------------#
# suffix rules
#---------------------------------------------------------------------#
.C$(OBJ_EXT):
$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<

17
Packages/Interpolation/examples/Interpolation/README Normal file
View File

@ -0,0 +1,17 @@
To compule and run all these examples type :
./cgal_test
To compute and run only some of them type
./cgal_test name-of_wanted_example
nn_coordinates_2: shows how to compute 2D natural coordinates
given a 2D Delaunay triangulation and a query point.
rn_coordinates_2: shows how to compute 2D regular (natural) coordinates
given a 2D Regular triangulation and a (weighted) query point.
linear_interpolation_2: interpolates a linear function using the
linear_interpolation function and 2D natural neighbor coordinates.
sibson_interpolation_2: interpolates a spherical function using the
sibson_gradient_fitting and sibson_c1_interpolation function with
2D natural neighbor coordinates.

67
Packages/Interpolation/examples/Interpolation/linear_interpolation_2.C Normal file
View File

@ -0,0 +1,67 @@
//
//file: examples/Interpolation/linear_interoplation_2.C
//
#include <CGAL/basic.h>
#include <utility>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_traits_2<K> Traits;
typedef K::FT Coord_type;
//Functor for accessing the function values
template< class Map >
struct DataAccess : public std::unary_function< typename Map::key_type,
typename Map::mapped_type> {
typedef typename Map::mapped_type Data_type;
typedef typename Map::key_type Point;
DataAccess< Map >(const Map& m): map(m){};
Data_type operator()(const Point& p) {
typename Map::const_iterator mit = map.find(p);
if(mit!= map.end())
return mit->second;
return Data_type();
};
const Map& map;
};
int main()
{
Delaunay_triangulation T;
std::map<Point, Coord_type, K::Less_xy_2> values;
typedef DataAccess< std::map<Point, Coord_type, K::Less_xy_2 > >
Value_access;
Coord_type a(0.25), bx(1.3), by(-0.7);
for (int y=0 ; y<3 ; y++)
for (int x=0 ; x<3 ; x++){
K::Point_2 p(x,y);
T.insert(p);
values.insert(std::make_pair(p,a + bx* x+ by*y));
}
//coordiante computation
K::Point_2 p(1.3,0.34);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm =
CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;
Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
norm,Value_access(values));
std::cout << " Tested interpolation on " << p << " interpolation: " << res
<< " exact: " << a + bx* p.x()+ by* p.y()<< std::endl;
return 1;
}

29
Packages/Interpolation/examples/Interpolation/nn_coordinates_2.C Normal file
View File

@ -0,0 +1,29 @@
//
//file: examples/Interpolation/nn_coordinates_2.C
//
#include <CGAL/basic.h>
#include <utility>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
int main()
{
Delaunay_triangulation T;
for (int y=0 ; y<3 ; y++)
for (int x=0 ; x<3 ; x++)
T.insert(K::Point_2(x,y));
//coordinate computation
K::Point_2 p(1.2, 0.7);
std::vector< std::pair< K::Point_2, K::FT > > coords;
K::FT norm =
CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;
return 1;
}

34
Packages/Interpolation/examples/Interpolation/rn_coordinates_2.C Normal file
View File

@ -0,0 +1,34 @@
//
//file: examples/Interpolation/rn_coordinates_2.C
//
#include <CGAL/basic.h>
#include <utility>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/Regular_neighbor_coordinates_traits_2.h>
#include <CGAL/regular_neighbor_coordinates_2.h>
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Regular_neighbor_coordinates_traits_2<K> Gt;
typedef CGAL::Regular_triangulation_2<Gt> Regular_triangulation;
typedef Regular_triangulation::Weighted_point Weighted_point;
int main()
{
Regular_triangulation rt;
for (int y=0 ; y<3 ; y++)
for (int x=0 ; x<3 ; x++)
rt.insert(Weighted_point(K::Point_2(x,y), 0));
//coordinate computation
Weighted_point wp(K::Point_2(1.2, 0.7),2);
std::vector< std::pair< K::Point_2, K::FT > > coords;
K::FT norm =
CGAL::regular_neighbor_coordinates_2(rt, wp,std::back_inserter(coords)).second;
return 1;
}

81
Packages/Interpolation/examples/Interpolation/sibson_interpolation_2.C Normal file
View File

@ -0,0 +1,81 @@
//
//file: examples/Interpolation/sibson_interoplation_2.C
//
#include <CGAL/basic.h>
#include <utility>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/Interpolation_gradient_fitting_traits_2.h>
#include <CGAL/sibson_gradient_fitting.h>
#include <CGAL/interpolation_functions.h>
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
typedef K::FT Coord_type;
typedef K::Point_2 Point;
typedef std::map<Point, Coord_type, K::Less_xy_2> Point_value_map ;
typedef std::map<Point, K::Vector_2 , K::Less_xy_2 > Point_vector_map;
//Functor class for accessing the function values/gradients
template< class Map >
struct DataAccess : public std::unary_function< typename Map::key_type,
typename Map::mapped_type> {
typedef typename Map::mapped_type Data_type;
typedef typename Map::key_type Point;
DataAccess< Map >(const Map& m): map(m){};
Data_type operator()(const Point& p) {
typename Map::const_iterator mit = map.find(p);
if(mit!= map.end())
return mit->second;
return Data_type();
};
const Map& map;
};
int main()
{
Delaunay_triangulation T;
Point_value_map values;
Point_vector_map gradients;
//parameters for spherical function:
Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
for (int y=0 ; y<4 ; y++)
for (int x=0 ; x<4 ; x++){
K::Point_2 p(x,y);
T.insert(p);
values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
}
sibson_gradient_fitting_nn_2(T,std::inserter(gradients,gradients.begin()),
DataAccess<Point_value_map>(values), Traits());
//coordiante computation
K::Point_2 p(1.6,1.4);
std::vector< std::pair< Point, Coord_type > > coords;
Coord_type norm =
CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;
//Sibson interpolant: version without sqrt:
Coord_type res = CGAL::sibson_c1_interpolation_square(coords.begin(),
coords.end(),norm,p,
DataAccess<Point_value_map>(values),
DataAccess<Point_vector_map>(gradients),
Traits());
std::cout << " Tested interpolation on " << p << " interpolation: " << res
<< " exact: "
<< a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
<< std::endl;
return 1;
}