Higher-dimensional Interpolation

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 and o2scl::interp2_planar performs interpolation using planes. At present, the only way to compute contour lines on data which is not defined on a grid is to use one of these two classes to recast 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() and o2scl::tensor_grid::interp_linear() . 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.

Interpolation on a rectangular grid

/* Example: ex_twod_intp.cpp
-------------------------------------------------------------------
A simple example for two-dimensional interpolation using
the interp_twod class.
*/
#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;
using namespace o2scl;
// A function for filling the data and comparing results
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;
// Create the sample data
ubvector x(3), y(3);
ubmatrix data(3,3);
// Set the grid
x[0]=0.0;
x[1]=1.0;
x[2]=2.0;
y[0]=3.0;
y[1]=2.0;
y[2]=1.0;
// Set and print out the data
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;
// Perform the interpolation
cout << "x y Calc. Exact" << endl;
// Interpolation, x-first
double tol=0.05;
double tol2=0.4;
ti.set_data(3,3,x,y,data,itp_cspline);
double x0, y0, x1, y1;
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;
// Interpolation, y-first
ti.set_data(3,3,x,y,data,itp_cspline);
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;
t.report();
return 0;
}
// End of example

This example creates a sample 3 by 3 grid of data with the function $ \left[ \sin \left( x/10 + 3 y/10 \right) \right]^2 $ 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.

Interpolation of randomly spaced points

For example, with 10 random points in the x-y plane with $ -1<x<1 $ and $ -1<y<1 $, the figure contains several polygonal regions, each of which represents the set of all points in the domain which will be mapped to the same plane in order to to approximate the original function.

ex_planar_plot.png
Planes from interp2_planar class

Contour lines

This example generates contour lines of the function

\[ z = f(x,y) = 15 \exp \left[ - \frac{1}{20^2}\left( x-20 \right)^2 - \frac{1}{5^2}\left(y-5\right)^2\right] + 40 \exp \left[ - \frac{1}{500}\left( x-70 \right)^2 - \frac{1}{2^2}\left(y-2\right)^2\right] \]

/* Example: ex_contour.cpp
-------------------------------------------------------------------
Example for generating contour lines
*/
#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;
using namespace o2scl;
using namespace o2scl_hdf;
// A function defining the three-dimensional surface
// for which we want to compute contour levels
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);
}
// A function for outputting the data to cout
int print_data(int nx, int ny, ubvector &x, ubvector &y,
ubmatrix &data);
// A function for printing the contour information to a file
int file_out(string prefix, ubvector &x, ubvector &y,
ubmatrix &data, vector<contour_line> &conts,
vector<edge_crossings> &xecs,
vector<edge_crossings> &yecs);
int main(void) {
cout.setf(ios::scientific);
contour co;
// Initialize the data
ubvector x(12), y(10);
ubmatrix data(12,10);
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]);
}
}
co.set_data(12,10,x,y,data);
// Print out the data
print_data(12,10,x,y,data);
// Set the contour levels
ubvector levels(7);
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;
co.set_levels(7,levels);
// Compute the contours
vector<contour_line> conts;
co.calc_contours(conts);
vector<edge_crossings> xecs, yecs;
co.get_edges(xecs,yecs);
// Output to a file
file_out("c1",x,y,data,conts,xecs,yecs);
// Print the contours to the screen and test to make sure
// that they match the requested level
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;
}
// Refine the data using cubic spline interpolation
co.regrid_data(5,5);
ubvector *x2;
ubvector *y2;
ubmatrix *data2;
size_t sx, sy;
co.get_data(sx,sy,x2,y2,data2);
// Recompute the contours
vector<contour_line> conts2;
co.calc_contours(conts2);
vector<edge_crossings> xecs2, yecs2;
co.get_edges(xecs2,yecs2);
// Output to a file
file_out("c2",*x2,*y2,*data2,conts2,xecs2,yecs2);
t.report();
return 0;
}
// End of example

The figure below shows contour lines in the region $ x\in(0,121), y\in(0,9) $. The data grid is represented by plus signs, and the associated generated contours. The figure clearly shows the peaks at $ (20,5) $ and $ (70,2) $ .

ex_contour_plot1.png
Contour example plot

The next figure below shows the edge crossings constructed by the contour class which are used to obtain the two $ z=10 $ contour lines. The edge crossings are connected to their respective data grid points with line segments. Endpoints are in red and internal contour points are in blue.

ex_contour_plot3.png
Contour edge crossings

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.

ex_contour_plot2.png
Contours after regrid_data()

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).