Two-dimensional interpolation
There are two types of two-dimensional interpolation classes, the first is based on a function defined on a two-dimensional grid (though the spacings between grid points need not be equal). The class o2scl::interp2_direct implements bilinear or bicubic interpolation, and is based on D. Zaslavsky's routines at https://github.com/diazona/interp2d (licensed under GPLv3). A slightly slower (but a bit more flexible) alternative is successive use of o2scl::interp_base objects, implemented in o2scl::interp2_seq .
If data is arranged without a grid, then o2scl::interp2_neigh performs nearest-neighbor interpolation. At present, the only way to compute contour lines on data which is not defined on a grid is to use this class or one of the multi-dimensional interpolation classes described below the data on a grid and then use o2scl::contour afterwards.
Multi-dimensional interpolation
Multi-dimensional interpolation for table defined on a grid is possible with o2scl::tensor_grid. See the documentation for o2scl::tensor_grid::interpolate(), o2scl::tensor_grid::interp_linear() and o2scl::tensor_grid::rearrange_and_copy() . Also, if you want to interpolate rank-1
indices to get a vector result, you can use o2scl::tensor_grid::interp_linear_vec() .
If the data is not on a grid, then inverse distance weighted interpolation is performed by o2scl::interpm_idw.
An experimental class for multidimensional-dimensional kriging is also provided in o2scl::interpm_krige .
Interpolation on a rectangular grid
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/interp2_seq.h>
#include <o2scl/test_mgr.h>
using namespace std;
double f(double x, double y) {
return pow(sin(0.1*x+0.3*y),2.0);
}
int main(void) {
cout.setf(ios::scientific);
int i,j;
typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
x[0]=0.0;
x[1]=1.0;
x[2]=2.0;
y[0]=3.0;
y[1]=2.0;
y[2]=1.0;
cout << endl;
cout << "Data: " << endl;
cout << " x | ";
for(i=0;i<3;i++) cout << x[i] << " ";
cout << endl;
cout << " y |" << endl;
cout << "-------------|-";
for(i=0;i<3;i++) cout << "-------------";
cout << endl;
for(i=0;i<3;i++) {
cout << y[i] << " | ";
for(j=0;j<3;j++) {
data(i,j)=f(x[j],y[i]);
cout << data(i,j) << " ";
}
cout << endl;
}
cout << endl;
cout << "x y Calc. Exact" << endl;
double tol=0.05;
double tol2=0.4;
x0=0.5; y0=1.5;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
x0=0.99; y0=1.99;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
x0=1.0; y0=2.0;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
cout << endl;
x0=0.5; y0=1.5;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
x0=0.99; y0=1.99;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
x0=1.0; y0=2.0;
cout << x0 << " " << y0 << " "
<< ti.
eval(x0,y0) <<
" " << f(x0,y0) << endl;
cout << endl;
return 0;
}
This example creates a sample 3 by 3 grid of data with the function
and performs some interpolations and compares them with the exact result.
Data:
x | 0.000000e+00 1.000000e+00 2.000000e+00
y |
-------------|----------------------------------------
3.000000e+00 | 6.136010e-01 7.080734e-01 7.942506e-01
2.000000e+00 | 3.188211e-01 4.150164e-01 5.145998e-01
1.000000e+00 | 8.733219e-02 1.516466e-01 2.298488e-01
x y Calc. Exact
5.000000e-01 1.500000e+00 6.070521e-01 2.298488e-01
9.900000e-01 1.990000e+00 4.187808e-01 4.110774e-01
1.000000e+00 2.000000e+00 4.150164e-01 4.150164e-01
5.000000e-01 1.500000e+00 6.070521e-01 2.298488e-01
9.900000e-01 1.990000e+00 4.187808e-01 4.110774e-01
1.000000e+00 2.000000e+00 4.150164e-01 4.150164e-01
0 tests performed.
All tests passed.
Contour lines
This example generates contour lines of the function
#include <iostream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/test_mgr.h>
#include <o2scl/contour.h>
#include <o2scl/hdf_io.h>
using namespace std;
double fun(double x, double y) {
return 15.0*exp(-pow(x-20.0,2.0)/400.0-pow(y-5.0,2.0)/25.0)
+40.0*exp(-pow(x-70.0,2.0)/500.0-pow(y-2.0,2.0)/4.0);
}
ubmatrix &data, vector<contour_line> &conts);
int main(void) {
cout.setf(ios::scientific);
for(size_t i=0;i<10;i++) {
y[i]=((double)i);
}
for(size_t i=0;i<12;i++) {
x[i]=((double)i)*((double)i);
}
for(size_t j=0;j<12;j++) {
for(size_t k=0;k<10;k++) {
data(j,k)=fun(x[j],y[k]);
}
}
print_data(12,10,x,y,data);
levels[0]=5.0;
levels[1]=10.0;
levels[2]=0.002;
levels[3]=20.0;
levels[4]=0.2;
levels[5]=30.0;
levels[6]=2.0;
vector<contour_line> conts;
file_out("c1",x,y,data,conts);
size_t nc=conts.size();
for(size_t i=0;i<nc;i++) {
cout << "Contour " << i << " at level " << conts[i].level << ":" << endl;
size_t cs=conts[i].x.size();
for(size_t j=0;j<cs;j++) {
cout << "(" << conts[i].x[j] << ", " << conts[i].y[j] << ") "
<< fun(conts[i].x[j],conts[i].y[j]) << endl;
t.
test_rel(fun(conts[i].x[j],conts[i].y[j]),conts[i].level,
1.0,"curve");
}
cout << endl;
}
size_t sx, sy;
vector<contour_line> conts2;
file_out(
"c2",*
x2,*y2,*data2,conts2);
return 0;
}
The figure below shows contour lines in the region
. The data grid is represented by plus signs, and the associated generated contours. The figure clearly shows the peaks at
and
.
Contour example plot
The o2scl::contour class can also use interpolation to attempt to refine the data grid. The new contours after a refinement of a factor of 5 is given in the figure below.
Contours after regrid_data()