fast interp on 1 dimension of 3d matrix

1. Assume a Matrix A which is (M X N X P) array.

I want to upsample one dimension of the matrix, for example: upsample the third dimension by factor of 2.

Currently I'm using loops to make it happen:

B = zeros(M, N, 2 * P);

for ii = 1 : M

for jj = 1 : N

B(ii, jj, :) = interp(double(A(ii, jj, :)), 2);

end

end

Is there any way to do so which doesn't include any use of loops?

[Bruno Luong]

Attention to those who don't know how to count:

S(1) = 4;
linspace(1,S(1),2*S(1))

ans = 1×8

1.0000 1.4286 1.8571 2.2857 2.7143 3.1429 3.5714 4.0000

is ugly. Only the two end points are from the original data after interpolation.

The much nicer (correct half step upsampling) is

linspace(1,S(1),2*S(1)-1)

ans = 1×7

1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000

Half of the data points are from the original.

[Stephen23]

@Bruno Luong: I moved your comment because this size was specified by the OP in their question.

See also: https://en.wikipedia.org/wiki/Off-by-one_error#Fencepost_error

[Bruno Luong]

@Stephen23 My comment is directed to your code with linspace, obviously violate the rule of integer upsampling https://en.wikipedia.org/wiki/Upsampling

Expansion: Create a sequence, , comprising the original samples,.

OP does not specify the size, he just allocate (wrongly) the size. He obviously want to use interp with upsampling integer factor

B(ii, jj, :) = interp(double(A(ii, jj, :)), 2)

And obviously induce @Lizhao Du to use the same flawed code.

Answer 1

[Stephen23]

Use interp3 or interpn:

>> A = randi(9,5,7,3)
A(:,:,1) =
   3   3   5   3   8   1   2
   8   7   3   4   9   9   2
   3   5   7   8   6   4   9
   2   9   1   7   3   6   6
   2   5   8   9   7   9   2
A(:,:,2) =
   3   6   9   5   3   5   4
   6   8   1   6   9   5   8
   8   5   5   9   4   2   9
   7   6   5   2   7   6   4
   8   4   2   6   1   6   4
A(:,:,3) =
   5   2   1   9   9   8   1
   9   8   5   5   8   8   2
   9   6   2   1   8   9   4
   5   1   4   9   7   4   2
   1   1   4   8   9   8   8
>> S = size(A);
>> Z = interpn(A,1:S(1),1:S(2),1:0.5:S(3)) % halve the step size
Z(:,:,1) =
   3.0   3.0   5.0   3.0   8.0   1.0   2.0
   8.0   7.0   3.0   4.0   9.0   9.0   2.0
   3.0   5.0   7.0   8.0   6.0   4.0   9.0
   2.0   9.0   1.0   7.0   3.0   6.0   6.0
   2.0   5.0   8.0   9.0   7.0   9.0   2.0
Z(:,:,2) =
   3.0   4.5   7.0   4.0   5.5   3.0   3.0
   7.0   7.5   2.0   5.0   9.0   7.0   5.0
   5.5   5.0   6.0   8.5   5.0   3.0   9.0
   4.5   7.5   3.0   4.5   5.0   6.0   5.0
   5.0   4.5   5.0   7.5   4.0   7.5   3.0
Z(:,:,3) =
   3.0   6.0   9.0   5.0   3.0   5.0   4.0
   6.0   8.0   1.0   6.0   9.0   5.0   8.0
   8.0   5.0   5.0   9.0   4.0   2.0   9.0
   7.0   6.0   5.0   2.0   7.0   6.0   4.0
   8.0   4.0   2.0   6.0   1.0   6.0   4.0
Z(:,:,4) =
   4.0   4.0   5.0   7.0   6.0   6.5   2.5
   7.5   8.0   3.0   5.5   8.5   6.5   5.0
   8.5   5.5   3.5   5.0   6.0   5.5   6.5
   6.0   3.5   4.5   5.5   7.0   5.0   3.0
   4.5   2.5   3.0   7.0   5.0   7.0   6.0
Z(:,:,5) =
   5.0   2.0   1.0   9.0   9.0   8.0   1.0
   9.0   8.0   5.0   5.0   8.0   8.0   2.0
   9.0   6.0   2.0   1.0   8.0   9.0   4.0
   5.0   1.0   4.0   9.0   7.0   4.0   2.0
   1.0   1.0   4.0   8.0   9.0   8.0   8.0

Generate the required indices either using the colon operator (as shown, e.g. to halve the step size) or linspace (e.g. to double the number of points):

Z = interpn(A,1:S(1),1:S(2),linspace(1,S(3),2*S(3))) % double the number of points
Z(:,:,1) =
   3.0   3.0   5.0   3.0   8.0   1.0   2.0
   8.0   7.0   3.0   4.0   9.0   9.0   2.0
   3.0   5.0   7.0   8.0   6.0   4.0   9.0
   2.0   9.0   1.0   7.0   3.0   6.0   6.0
   2.0   5.0   8.0   9.0   7.0   9.0   2.0
Z(:,:,2) =
   3.0   4.2   6.6   3.8   6.0   2.6   2.8
   7.2   7.4   2.2   4.8   9.0   7.4   4.4
   5.0   5.0   6.2   8.4   5.2   3.2   9.0
   4.0   7.8   2.6   5.0   4.6   6.0   5.2
   4.4   4.6   5.6   7.8   4.6   7.8   2.8
Z(:,:,3) =
   3.0   5.4   8.2   4.6   4.0   4.2   3.6
   6.4   7.8   1.4   5.6   9.0   5.8   6.8
   7.0   5.0   5.4   8.8   4.4   2.4   9.0
   6.0   6.6   4.2   3.0   6.2   6.0   4.4
   6.8   4.2   3.2   6.6   2.2   6.6   3.6
Z(:,:,4) =
   3.4   5.2   7.4   5.8   4.2   5.6   3.4
   6.6   8.0   1.8   5.8   8.8   5.6   6.8
   8.2   5.2   4.4   7.4   4.8   3.4   8.0
   6.6   5.0   4.8   3.4   7.0   5.6   3.6
   6.6   3.4   2.4   6.4   2.6   6.4   4.8
Z(:,:,5) =
   4.2   3.6   4.2   7.4   6.6   6.8   2.2
   7.8   8.0   3.4   5.4   8.4   6.8   4.4
   8.6   5.6   3.2   4.2   6.4   6.2   6.0
   5.8   3.0   4.4   6.2   7.0   4.8   2.8
   3.8   2.2   3.2   7.2   5.8   7.2   6.4
Z(:,:,6) =
   5.0   2.0   1.0   9.0   9.0   8.0   1.0
   9.0   8.0   5.0   5.0   8.0   8.0   2.0
   9.0   6.0   2.0   1.0   8.0   9.0   4.0
   5.0   1.0   4.0   9.0   7.0   4.0   2.0
   1.0   1.0   4.0   8.0   9.0   8.0   8.0

[Nitai Fingerhut]

Vq = interp3(V,k) returns the interpolated values on a refined grid formed by repeatedly halving the intervals k times in each dimension. This results in 2^k-1 interpolated points between sample values.

I want to interpolate only along the third dimension and not along all 3 dimensions.

[Stephen23]

"I want to interpolate only along the third dimension and not along all 3 dimensions."

That is exactly what my answer does.

Look at the example: the first and second dimensions do not change size. I used a small example array so that you could count the rows/columns easily. Try counting them. They do not change.

The code you show in your comment has nothing to do with my answer.

[Nitai Fingerhut]

I did try it.

In fact I have a 4-D matrix, and I'm trying to interpolate along the first dimension, so I used:

[x,y,z,w] = size(A);

A = interpn(double(A) , 1:(2*x), 1:y, 1:z, 1:w);

but I got:

"Error using griddedInterpolant/subsref

The input data has inconsistent size."

What I did wrong?

Thanks

[Nitai Fingerhut]

As you can see, I tried to double the smapling points of the first dimension only, but with no success.

[Stephen23]

"What I did wrong?"

You ignored my answer. You invented this:

A = interpn(double(A) , 1:(2*x), 1:y, 1:z, 1:w);
                        ^^^^^^^ where do you see this in my answer?

I showed you how to use linspace to double the number of points. Why not try using linspace, exactly like I showed you in my answer?:

B = interpn(double(A) , linspace(1,x,2*x), 1:y, 1:z, 1:w);

[Nitai Fingerhut]

I tried linsapce too...

like this:

A = interpn(double(A) , 1:(2*x), 1:y, 1:z, 1:w);

still didn't work...

[Nitai Fingerhut]

Attaching screen shot

note that "row" is A

[Stephen23]

"I tried linsapce too..."

"like this:"

A = interpn(double(A) , 1:(2*x), 1:y, 1:z, 1:w);

Your code does not have linspace in it anywhere. Where is the linspace ?

[Stephen23]

@Nitai Fingerhut: Please upload your variable raw in a mat file.

[Nitai Fingerhut]

Look at the screenshot, I used linspace.

raw is 4-D int 16, could you please explaine how to update it to a mat file?

and does it nessecary ?

[Stephen23]

First save that variable in a .mat file:

save('temp.mat','raw')

Then make a new comment and upload the file temp.mat by clicking the paperclip button.

[Nitai Fingerhut]

Apperentally it's too big for uploading...

[Nitai Fingerhut]

Apperentally it's too big for uploading...

tried to zip it but it's still too big...

[Stephen23]

@Nitai Fingerhut: Please show the output of these commands:

whos raw
size(raw)

[Nitai Fingerhut]

Attached

[Stephen23]

@Nitai Fingerhut: what MATLAB version are you using?

[Nitai Fingerhut]

matlab 2018a

[Stephen23]

https://www.mathworks.com/help/matlab/release-notes.html?rntext=interpn&startrelease=R2010a&endrelease=R2018b&rntype=incompatibility&groupby=release&sortby=descending&searchHighlight=interpn

Unfortunately for R2012a and later versions you will need to generate index matrices yourself using ndgrid, which you could do like this:

[W,X,Y,Z] = ndgrid(linspace(1,S(1),2*S(1)), 1:S(2), 1:S(3), 1:S(4));
 Z = interpn(A, W,X,Y,Z);
 

Unfortunately this generates large intermediate arrays, which will take a lot of memory and be quite slow.... but here are two (probably faster) alternatives:

1 . Avoid the large index matrices by calling griddedInterpolant directly (this lets you pick which dimension/s you want to interpolate over):

 I = griddedInterpolant({1:S(1), 1:S(2), 1:S(3), 1:S(4)}, A);
 Z = I({linspace(1,S(1),2*S(1)), 1:S(2), 1:S(3), 1:S(4)});
 

2 . Another solution is to call interp1 without a loop as long as you only want to interpolate over the first dimension, in which case it will simply interpolate the columns of the array:

 Z = interp1(A,linspace(1,S(1),2*S(1)))

That is probably the fastest solution, but it only interpolates along the first dimension. If you want to interpolate along another dimension then you could use permute before and after.

[Nitai Fingerhut]

Last ones worked!

So, for future uses- I assume that best solution is to rearrange matrix so the desired dimension will be first, and then use interp1.

Thanks!

[Bruno Luong]

The difference between

Z = interp1(A,linspace(1,S(1),2*S(1)))

and your original INTERP code is that the INTERP1 is not integer-upsampling.

INTERP do integer upsampling and must extrapolate on the right side with (q-1) points where q is up-sampling factor (2 in your example).

It carried out but padding one 0 and filter the padded data, so the extrapolation signal dies out on the right side and can introduce some oscilation.

The difference, though subtle, might be important or not depending on your ultimate goal.

[Lizhao Du]

Z = interp1(A,linspace(1,S(1),2*S(1)))

It does help!

Thanks to Stephen!

Answer 2

[madhan ravi]

 interp3(A,2)

[Nitai Fingerhut]

Vq = interp3(V,k) returns the interpolated values on a refined grid formed by repeatedly halving the intervals k times in each dimension. This results in 2^k-1 interpolated points between sample values.

I want to interpolate only along the third dimension and not along all 3 dimensions.

Answer 3

[Bruno Luong]

[m,n,p] = size(A);
B = reshape(A,[m*n p]).';
Br = resample(B,2,1); % 2 is upsampling factor
Ar = reshape(Br.',m,n,[]); 

Answer 4

[Matt J]

Using imresizen (Download)

B = imresizen(A,[1,1,2]);

[Nitai Fingerhut]

what kind of interpolation does this function uses?

In addition- I noticed that "imresizen" uses a for loop so it won't improve running time...

[Matt J]

By default, it uses linear interpolation, but you can specify any interpolation method that griddedInterpolant implements, e.g.,

B = imresizen(A,[1,1,2],'pchip');

[Matt J]

In addition- I noticed that "imresizen" uses a for loop so it won't improve running time...

Incorrect. The for loop in imresizen only runs in your case from 1 to 3 (the number of dimensions) and does trivially fast pre-computations.