How interp2 deal with edges on bicubic interpolation?

Hi, I'm trying to make my own algorithm for 2d interpolation.

The upper left figure shows the data used for the interpolation. The upper right is the data interpolated with a linear method. The lower left is the result of matlab interp 2 cubic method and the lower right is an algorithm I made for cubic interpolation in 2d.

I wanted to know what interp2 does when you have data on the edge (for example the point where the X is). You need 2 points both sides to interpolate with the cubic method so what I did is interpolate linearly on the edge of the map. I'm on MATLAB 2016.

At the moment, my method uses separable convolution so I interpolate on rows first and then I interpolate on columns.

From what I understand from the code, interp2 extrapolates to be able to interpolate near edges, but I don't know which method is used to do so? Is it simply linear extrapolation?

Answer 1

[Bruno Luong]

The boundary handling is descriibed in the section Boundary Condition of this Cubic interpolation reference

Especially the equation (19) gives the value of the left outside pixel from the inside pixels

c(-1) = 3*f(x0) - 3*f(x1) + f(x2)

Likewise eqt (25) is for the right-side

c(N+1) = 3*f(x_N) - 3*f(x_{N-1}) + f(x_{N-2})

[Matt J]

That does indeed seem to be a closed form formula for quadratic extrapolation,

 format long 

 f=rand(1,3);

 p=polyfit(0:2, f,2); 
 
 extrap=polyval(p,-1), %quadratic extrapolation

extrap = 
   0.151149787053714

 
 extrap0=dot([3,-3,1],f)  %official source

extrap0 = 
   0.151149787053714

[Bruno Luong]

Yes the equation is derived by zeroing the cubic term, it is written in the paper

[Jérémy Talbot-Pâquet]

Thank you that's what I was searching for. The article make it very clear!

[Matt J]

The boundary handling is descriibed in the section Boundary Condition of this ...

But how do we know authoritatively that this article is what Matlab uses? I can't find its citation anywhere.

[Bruno Luong]

In this thread that you have participated Matt.

I have also showed later the tensorial form in 2D here.

[Matt J]

Ah, such nostalgia. Interestingly, though, the cubic option has changed since R2013 as mentioned in the documentation,

https://www.mathworks.com/help/matlab/ref/interp1.html#mw_f7a8cabf-4c61-4fd7-989b-13e037f87e5b

so I wonder if the article is still an accurate reference.

[Bruno Luong]

My test in the 2D thread shows the formula used in the paper seems still to hold.

[Bruno Luong]

Dig a little deeper, I dig out a version R2018a and the

interp1(..., 'cubic') 

is the same as

interp1(..., 'pchip')

for R2021a.

Here is the script to check it that needed to be run on R201a, then on R2018b

x = 1:10;
y = [7     1     3     1     1     8     7     3     9     1]; % generate using y = ceil(9*rand(size(x)))

xq = linspace(min(x),max(x),129);
s=ver('matlab');
if strcmp(s.Version, '9.10') % R2021a
    yq_R2021a_pchip = interp1(x, y, xq, 'pchip');
    yq_R2021a_cubic = interp1(x, y, xq, 'cubic');
    save yq_R2021a.mat yq_R2021a_pchip yq_R2021a_cubic
else % strcmp(s.Version, '9.5') % R2018b
    yq = interp1(x, y, xq, 'cubic');
    close all
    figure
    load('yq_R2021a.mat');
    plot(x, y, 'or', ...
        xq, yq, '-r', ...
        xq, yq_R2021a_pchip, 'b.', ...
        xq, yq_R2021a_cubic, 'g+' ...
        );
    legend('data', 'cubic R2018b', 'pchip R2021a', 'cubic R2021a', 'location', 'northwest');
end

NOTE: The reference paper is applicable for R2021a 'cubic' method for those who still follow.

[Matt J]

But you dug up the reference paper from an Answers post from back in 2013. How can a post from the year 2013 anticipate what 'cubic' method would change to in the year 2021??

[Bruno Luong]

But that why I don't understand how this paper is mentionned in 2013 by the TMW staff as 'cubic'. I guess this is actually 'v5cubic' for all versions.

What I know is this paper gives the right formula for R2021a cubic. That is the fact.

IMO TMW mess up with the doc.

[Josh Meyer]

It's important to note that the reference paper mentioned in the 2013 thread applied to interp2(..,'cubic') and interp1(..,'v5cubic') when it was originally posted. However, at that time interp1(..,'cubic') had different behavior and was the same as interp1(..,'pchip'), so the reference paper didn't apply there. In other words, interp2 and interp1 had similar options named 'cubic' that actually were different methods, and I believe this is the source of your confusion.

As of R2020b the behavior of interp1(..,'cubic') is the same as interp1(..,'v5cubic'), so the reference paper now applies for interp1's cubic option as well. The details of this change are at the bottom of the interp1 page.

[Bruno Luong]

interp(...'pchip') is NOT the same as interp1(..., 'v5cubic'). The later uses convolution formula on the paper. And 'pchip' preserves data monotony and it is a non-linear method.

So 'cubic' in R2013 cannot match both IMO.

Anyway it is not big deal, I'm only interested in what happens now.

Answer 2

[Kelly Kearney]

interp2 does not extrapolate. By default, it sets any points outside the data bounds to NaN; you can alternatively choose a different scalar value to assign to all extrapolated points. Alternatively, you could use griddedInterpolant, which allows you to explicitly assign an algorithm to the extrapolated points:

[x,y,z] = peaks(5);
 
[xi,yi] = meshgrid(linspace(-4,4,50));
zi1 = interp2(x,y,z,xi,yi,'cubic');
F = griddedInterpolant(x',y',z,'cubic','linear');
zi2 = F(xi',yi');
 
subplot(2,1,1);
 
hi = imagesc(xi(1,:),yi(:,1),zi1, 'AlphaData', ~isnan(zi1), 'alphadatamapping', 'scaled');
hold on;
scatter(x(:), y(:), [], z(:), 'filled', 'markeredgecolor', 'w');
title('interp2: cubic')
 
subplot(2,1,2);
hi = imagesc(xi(1,:),yi(:,1),zi2, 'AlphaData', ~isnan(zi2), 'alphadatamapping', 'scaled');
hold on;
scatter(x(:), y(:), [], z(:), 'filled', 'markeredgecolor', 'w');
title('griddedInterpolant: cubic, linear')

[Jérémy Talbot-Pâquet]

Thank you for your answer. Sorry if I wasn't clear, I meant that somewhere in interp2 for the cubic method, the array is expanded outside the domain to be able to interpolate on the edges since you need 2 points both sides. But you're right, interp2 does not return an extrapolated array unless you give an extrapval argument.

I haven't thought of GriddedInterpolant it could help me to expand the array to more than 1 point. Thank you

Answer 3

[Matt J]

Based on a little reverse engineering, I believe I have discovered that the 'cubic' method of interp1 uses quadratic extrapolation to pad 1 element at the beginning and end of the array. I think we can be pretty confident that the obvious generalization to the 2D case is used for interp2.

To verify this, the test below checks that the extrapolation rule leads to the same interpolated values regardless of whether the signal starts at the beginning of the array or in the middle:

v1=[rand(1,5) 0 0 0 0 0],

v1 = 1×10

0.1348 0.1215 0.0821 0.6977 0.1237 0 0 0 0 0

t1=0:numel(v1)-1;
F1=@(x) interp1(t1,v1,x,'cubic');

 p=polyfit(0:2, v1(1:3),2); 
 extrap=polyval(p,-1), %quadratic extrapolation

extrap = 0.1222

v2=circshift(v1,5);  %shift v1 and insert extrapolated value
  v2(5)=extrap,

v2 = 1×10

0 0 0 0 0.1222 0.1348 0.1215 0.0821 0.6977 0.1237

t2=t1-5; 
F2=@(x) interp1(t2,v2,x,'cubic');

t=linspace(0,1,30);
plot(t,F1(t),'-',t,F2(t),'o')
xlabel t 
ylabel 'Interpolated Values'

legend('v1(t1)','v2(t2)')