Bilinear interpolation on grid data

[philipp]

Dear Experts,

I guess it is a very simple task, but somehow I can not get my head around.

There is some data regarding a sensor calibration.
Calibration data consists of:

• reference temperature = x-axis
• sensor signal = y axis
• reference pressure = z-axis

Now:
For a given measurement I know:

• temperature
• sensor signal

The question is: How to get the pressure from the calibration grid?

I understand that the term can either be

• bilinear interpolation or
• intersection point of a line and a surface

Notes:

• the gray surface is plotted via mesh()
• white grid corner dots are plotted via scatter3d()
• green dots - connected with black line - are plotted using scatter3d

I have no function that describes the surface.

My confusion starts, when I try to find the 4 grid points that are next to xi, yi.

Thank you,
Philipp

[mottelet]

Hi Philipp, did you have a look to linear_interpn - N dimensional linear interpolation and splin2d - Bicubic spline gridded 2d interpolation ?

[philipp]

Hello Stéphane,

yes I did, but I am confused.

For my understanding these interpolation methods need a function, that descibes the z-component at any given x-y-component… The examples work with: z = f(x,y)

As I mentioned, I do not have such a function.
I have only the grid points.

I guess I do not see the obvious, but do I not need to find first the 4 corner points that are “around” the intersection of surface and line?
If I know these grid points I could go on with linear interpolation from there.

[mottelet]

Philipp, is is the opposite: 2D interpolation aims to provide a function z=f(x,y) (at least continuous) such that z_{ij}=f(x_i,y_j) where (x_i,y_j,z_{ij}) is your original gridded data. Suppose that you have

x=[273 300 330]; // temperature
y=[1 2 3]; // signal
P=[0.068374    0.7263507   0.2320748
   0.5608486   0.1985144   0.2312237
   0.6623569   0.5442573   0.2164633]; 
// arbitrary random pressure P(i,j) for x(i) and y(j)

then you can obtain e.g. the bilinearly interpolated pressure for x=280 and y=1.5 by

--> p = linear_interpn(280,1.5, x, y, P)
 p  = 

   0.3927784

The algorithm in linear_intern finds itself the relevant patch, i.e. [x_i,x_{i+1}] \times [y_j,y_{j+1}] such that (280,1.5) belong to it and yields the value of the interpolant.

[philipp]

Ok, thank you for the explanations.

I could made the code run, but do not get satisfactory results with linear_interpn.
I guess this is, because for each “P” there is in general a different X-Y-combination.

Best an example:

x = [
-50.069 -40.155 -30.329 -20.494 -10.422 -0.166 9.776 19.835 29.853
-50.192 -40.295 -30.513 -20.658 -10.628 -0.309 9.642 19.694 29.776
-50.274 -40.281 -30.569 -20.685 -10.706 -0.393 9.617 19.644 29.776
-50.261 -40.305 -30.577 -20.766 -10.72 -0.442 9.518 19.607 29.789
-50.273 -40.334 -30.625 -20.806 -10.766 -0.472 9.465 19.545 29.902
-50.248 -40.304 -30.648 -20.788 -10.79 -0.575 9.408 19.519 30.007
-50.273 -40.302 -30.62 -20.824 -10.859 -0.687 9.311 19.374 30.072
-50.158 -40.223 -30.644 -20.849 -10.899 -0.763 9.219 19.332 30.195
-49.975 -40.092 -30.588 -20.804 -10.798 -0.78 9.193 19.266 30.444
-49.836 -39.967 -30.528 -20.753 -10.809 -0.807 9.157 19.263 30.683];

y = [
12774976 12805376 12833152 12862048 12895456 12925600 12952896 12978880 13011680
11989504 12029856 12092608 12130080 12170784 12232416 12268448 12304672 12366464
11204992 11247520 11326240 11387040 11467360 11546592 11586304 11672032 11705312
10376896 10523232 10601664 10714496 10740512 10862688 10957664 10972896 11067936
9592768 9745312 9888000 9997920 10037888 10161696 10252736 10293728 10392096
8832032 8973184 9103264 9197088 9349504 9497312 9557568 9648224 9747936
8135936 8230336 8353664 8528384 8749344 8859328 8907520 9102304 9205088
7323136 7513088 7736416 7905504 8033408 8115264 8298144 8439520 8588896
6631776 6873984 7043200 7282496 7367488 7565632 7712992 7913824 7982368
6291008 6485920 6690656 6928512 7104096 7261024 7399424 7545696 7610368];

P = [
-38.831895 -29.184012 -19.584704 -9.8188522 0.4590955 11.06052 21.283058 31.863397 42.888896
-28.907894 -21.936448 -15.355545 -8.0069926 -0.2270507 8.2110005 16.20254 24.609555 34.192521
-18.879508 -14.481853 -10.775898 -6.0350814 -0.7359974 5.4242796 11.30656 17.954163 25.356559
-8.2909265 -7.6997544 -6.5695943 -4.4490732 -1.236976 2.6588686 6.5510735 10.531062 16.963747
1.5599479 -0.4592106 -2.4972049 -2.8172288 -2.0305449 -0.3484227 1.2123694 3.187839 8.110167
11.065405 6.5586554 1.696941 -1.1204313 -2.9055768 -3.4234306 -4.2476204 -3.886095 -0.3321871
19.455227 13.107699 5.4951288 -0.0756867 -3.9560843 -6.6645999 -9.6718031 -10.20337 -7.4629245
29.237962 19.254572 8.2630596 0.6442607 -5.3648821 -10.649392 -15.090675 -18.004928 -15.605209
37.261181 24.495563 11.224753 1.0880098 -6.9271407 -13.929419 -20.412474 -24.532706 -23.575441
41.079614 27.497264 12.413013 1.0735465 -7.9455568 -16.034039 -23.61702 -29.295384 -28.404332];

// goal: get a “P-Value” for point i
xi = 0;
yi = 9 * 10^6;

Version 1:
pi = linear_interpn(xi, yi, x, y, P); // does not work because x, y are not strictly increasing

Version 2:
// build strictly increasing vectors
x = linspace(min(x), max(x), 9); // there are 9 different values per X row
y = linspace(min(y), max(y),10); // there are 10 different values per Y column

pi = linear_interpn(xi, yi, x, y, P); // works, but gives unsatisfactory result

[mottelet]

OK I thought your x and y were on a cartesian grid but this not the case. I have tried cshep2d on your data but it fails because the x and y values vary in a very correlated manner.
In fact each of x,y and P in your data are naturally mapped from a rectangular grid, i.e. x(i,j), y(i,j), P(i,j) where ì=1:10 and j=1:9 so the idea for a given xk,yk pair is to find the non-integer parameter pair (i,j)such that x(i,j)=xk and y(i,j)=yk. To give a sense to these non-integer coordinates one has to interpolate x,y,P independently with the same method, e.g. bilinear method (2D splines could also work well). I propose the following code:

function out=f(par,xij,yij,x,y)
    i = par(1); j = par(2);
    out = [xij-linear_interpn(i,j,1:10,1:9,x,'natural')
           yij-linear_interpn(i,j,1:10,1:9,y,'natural')];
endfunction

// goal: get a “P-Value” for point k
xk = 0;
yk = 9e6;

// determine parameters value for given (xk,yk)
par = fsolve([5;5],list(f,xk,yk,x,y))
// compute P-value for computed parameters (i,j)
pk = linear_interpn(par(1), par(2), 1:10, 1:9, P)

which yields for your trial values xk and yk

 par  = 

   6.7847854
   6.0663243

 pk  = 

  -6.1353431

the code recovers exact values for xk,yk values matching reference points of your grid:

xk = -50.069;
yk = 12774976;
par = fsolve([5;5],list(f,xk,yk,x,y))
pk = linear_interpn(par(1), par(2), 1:10, 1:9, P, 'natural')

 par  = 

   1.0000000
   1.0000000

 pk  = 

  -38.831895

[philipp]

Excellent.
Thank you very much.

[mottelet]

The "natural" argument allows extrapolation, which seems to be necessary for fsolve, and allows the detection of trial (xk,yk) pairs which could be out of the initial grid. In that case you get values of i or j outside of 1:10 or 1:9, respectively.

[philipp]

Hello Stéphane,

thank you again, but I need to come back to this topic with some questions:

1st.: fsolve has an init value (x0)
→ Is it recommended that the entries of x0 should be within min/max values of the x and y values?

2nd: Extrapolation
Is it somehow possibe to define how many grid points are used for extrapolation?
For “natural” I understand, that the nearest n-linear patch is used …so I assume “natural” will use the two nearest grid points (in each direction) for extrapolation.

3rd: May I ask for an example of a bilinear approach?

[mottelet]

@philipp

1st.: fsolve has an init value (x0)
→ Is it recommended that the entries of x0 should be within min/max values of the x and y values?

No, as your data is considered as a parametric surface with parameters (i,j) in in [1,9]\times[1,10], and as the idea is to first find the parameters corresponding to a given (x,y) trial pair, the starting values should be in the parameter space, hence in [1,9]\times[1,10]. That’s why in the example I took [5;5], i.e. something quite close to the middle for starting ‘fsolve’.

@philipp

2nd: Extrapolation
Is it somehow possibe to define how many grid points are used for extrapolation?
For “natural” I understand, that the nearest n-linear patch is used …so I assume “natural” will use the two nearest grid points (in each direction) for extrapolation.

Exactly, that’s what is explained in the documentation of linear_interpn. If you want to catch the influence of more points you should use splin2d and interp2d.

@philipp

3rd: May I ask for an example of a bilinear approach?

In fact, this is actually bilinear interpolation because here linear_interpn is used with two independant coordinates.

S.

[mottelet]

@mottelet

If you want to catch the influence of more points you should use splin2d and interp2d.

Here is what I meant:

function out=f(par,xij,yij,x,y)
    i = par(1); j = par(2);
    out = [xij-interp2d(i,j,1:10,1:9,sx,'natural')
           yij-interp2d(i,j,1:10,1:9,sy,'natural')];
endfunction

sx = splin2d(1:10,1:9, x,'natural');
sy = splin2d(1:10,1:9, y,'natural');
sP = splin2d(1:10,1:9, P,'natural');

// goal: get a “P-Value” for point k
xk = 0;
yk = 9e6;

// determine parameters value for given (xk,yk)
par = fsolve([5;5],list(f,xk,yk,x,y))
// compute P-value for computed parameters (i,j)
pk = interp2d(par(1), par(2), 1:10, 1:9, sP, 'natural')

Compared to the bilinear interpolation, the 2d spline catches more curvature information and should (or should not) better extrapolate. The default type type of spline (not a knot) does not give good extrapolation that’s why ‘natural’ is used (it has a different meaning than the ‘natural’ of interp2d).

But depending of your goal, maybe interpolation is not the good way to go. I you want to extrapolate, you should rather use a global model (e.g. polynomial) with a small number of parameters.

Here are two illustrations of the extrapolation (mesh is in red) depending of the used interpolation type:

// bilinear interpolation
[i,j]=meshgrid(0:0.5:11,0:0.5:10);
X=linear_interpn(i,j,1:10,1:9,x,'natural');
Y=linear_interpn(i,j,1:10,1:9,y,'natural');
P=linear_interpn(i,j,1:10,1:9,P,'natural');
k = 3:$-2;
mesh(X,Y,P,'edgecolor','red')
mesh(X(k,k),Y(k,k),P(k,k))

// 2d spline interpolation
[i,j]=meshgrid(0:0.5:11,0:0.5:10);
X=interp2d(i,j,1:10,1:9,sx,'natural');
Y=interp2d(i,j,1:10,1:9,sy,'natural');
P=interp2d(i,j,1:10,1:9,sP,'natural');
k = 3:$-2;
mesh(X,Y,P,'edgecolor','red')
mesh(X(k,k),Y(k,k),P(k,k))

I don’t like the amplification of oscillations in the lower right corner…

S.

[philipp]

@mottelet

No, as your data is considered as a parametric surface with parameters (i,j) in [1,9] x [1,10], and as the idea is to first find the parameters corresponding to a given (x,y) trial pair, the starting values should be in the parameter space, hence in [1,9] x [1,10]. That’s why in the example I took [5;5], i.e. something quite close to the middle for starting ‘fsolve’.

Mh…so x0 is defined as an index within the actual x-y-matrices?

More like:
As a start point use the [5,5]-th coordinate of the x-y-coordinates
par = fsolve( [5;5], … )
→ fsolve will find and use -10.766 and 10037888 as a starting point,
So [5;5] tells fsolve the index of where to look at within x and y.

And not:
x0 = x(5,5) // = x(45) = -10.766
y0 = y(5,5) // = y(45) = 10037888
par = fsolve( [ -10.766; 10037888], … ) // equals: par = fsolve( [x0; y0], … )
which would give fsolve some initial values

If so, why does the fsolve help page states:
x0: real vector (initial value of function argument).

[philipp]

@mottelet

Compared to the bilinear interpolation, the 2d spline catches more curvature information and should (or should not) better extrapolate. The default type type of spline (not a knot) does not give good extrapolation that’s why ‘natural’ is used (it has a different meaning than the ‘natural’ of interp2d).

But depending of your goal, maybe interpolation is not the good way to go. I you want to extrapolate, you should rather use a global model (e.g. polynomial) with a small number of parameters.

Well yes, I tried to use a polynomial approach first but wasn’t even satisfied with a quadratic approach. As a result my polynom would not exactly meet the grid points and show non-neglitiable errors when extrapolating.

I might go for an linear interpolation when x-y is within the specified range and - if the 2d-spline appraoch extrapolates Ok - I will use the 2d-spline approach for values outside of x-y.

Thanks again for all explanations and code examples.

Philipp

[mottelet]

@philipp

→ fsolve will find and use -10.766 and 10037888 as a starting point,
So [5;5] tells fsolve the index of where to look at within x and y.

No, in the above code fsolve iterates on par=[i,j] not x,y values.

S.

[mottelet]

@philipp

As a result my polynom would not exactly meet the grid points and show non-neglitiable errors when extrapolating.

How do you compute the error when extrapolating since you do not have data to compare with ? A polynomial approximation will always better extrapolate than using exact interpolation.

S.

[philipp]

Ok, some more background.

The whole thing has to do with sensors and sensor calibration …

There is a different approach to calculate “P”, which is currenty said to be “the truth”.
This approach uses a set of known coefficients to calculate P. This approach is thought to be valid even outside of the x-y-range but the coefficients can make it somewhat confusing when dealing with a lot of sensors.

The goal here is to avoid these coefficients by using some kind of interpolation method based on some known points [x, y, P] . Thanks to your support it is possible to find very good results as long as the measurements stay within the calibration range (within x-y-limits).

However measurements (measure x,y → calculate P) may also lay outside of the given x-y-range…so there is need for at least some extrapolation.

The first idea was to fit the calibration data via a polynom of max 3rd grade.
This would result in a single formula for each sensor, hence less confusion.

But fitting the data with a polynom showed at least two drawbacks:

1.) P-values calculated at the grid-points would not exactly meet the given P-value.
(which is understandable, since it is a datafit)

2.) P-values outside the x-y-range showed unreasonable results. This could be seen due to comparison with the “coefficient approach” but also by just observing the behaviour outside of the x-y-range.
Something like what you mentioned: “I don’t like the amplification of oscillations in the lower right corner…”
→ P should not start to oscillate when leaving the x-y-range.

P.

[mottelet]

@philipp

P-values outside the x-y-range showed unreasonable results. This could be seen due to comparison with the “coefficient approach” but also by just observing the behaviour outside of the x-y-range.

So you have a “coefficient approach” method which can be used for validation. I that case why don’t you use it for training the model in the maximum range of (x,y) ? There is always something unfair when you train a model with restricted data and validate with data outside the domain of training data, unless all the data is used to alternatively train and validate the model by using cross-validation.

S.

[philipp]

@mottelet

There is always something unfair when you train a model with restricted data and validate with data outside the domain of training data

Good point, agree to that.

[philipp]

@mottelet

So you have a “coefficient approach” method which can be used for validation. I that case why don’t you use it for training the model in the maximum range of (x,y) ?

“Which can be used for validation”…
This is questionable:
The “coefficient method” is currently “the truth” but it is not certain that “the truth” is correct. Our measurements leave the operating region the sensors are specified for. The sensors won’t get damaged, only the manufacturer specs do not cover the measurement range.
In contrary the “coefficent method” is itself based on a set of calibration points in a much smaller range. From this limited set of calibration points, the manufacturer expands to a wider region using a set of equations.

Thatswhy the sensors got a dedicated calibration which covers a wider range …still measurements where even outside of this extended calibration range.