2

I have a 3D matrix in Matlab of certain size, however I would need to interpolate it to obtain matrix of larger size.

size(M)
ans= 
  50   108    86

I need that matrix to be interpolated and finally obtain a matrix of size 100x213x140. Any ideas of how to do it using interp3? Is this possible at all?

I've tried

Vq = interp3(M,1:100,1:213,1:140)
Error using griddedInterpolant/subsref
The input data has inconsistent size.

Error in interp3 (line 178)
    Vq = F(Xq,Yq,Zq);

If I use meshgrid:

[X,Y,Z] = meshgrid(1:100, 1:213, 1:140);
Vq =interp3(M,X,Y,Z);

Matlab seems to like it, however two things happen:

  1. size(Vq) ans= 213 100 140
  2. I can see NaN values in Vq

The reason behind is because I need to compare two matrices sampled at different frequency. So, I could either interpolate M to obtain a matrix of size 100x213x140 or "reduce" the size of my other matrix M2 of size 100x213x140 to 50x108x86. I thought the former should be easier and safer...

Manolete
  • 3,431
  • 7
  • 54
  • 92
  • Are you trying to **resize** the matrix so that it looks like the original `M` but twice as big, or do you simply want to interpolate values? – rayryeng Jul 15 '14 at 19:57
  • @rayryeng I want to resize it, but I also need values for those new coordinates, hence the interpolation – Manolete Jul 15 '14 at 20:00
  • I'll write a post for you now. – rayryeng Jul 15 '14 at 20:04
  • I had to write code for you to do so. There is no built-in function in MATLAB. I will modify my post. – rayryeng Jul 15 '14 at 20:16
  • I have written code for you. Check my post. Good luck! BTW, not only is this an interpolation, this is also a resizing too. Very good question. – rayryeng Jul 15 '14 at 20:30
  • 1
    To resize a 3D matrix you could also use this solution: http://stackoverflow.com/questions/12520152/resizing-3d-matrix-image-in-matlab/24650504#24650504 – nightlyop Jul 16 '14 at 07:55
  • @nightlyop Can you tell me what is the difference between the method used by rayryeng here and the one explained in that thread? – Manolete Jul 17 '14 at 11:01
  • @Manolete: I'm not quite sure and mybe @rayryeng can give you a better answer. I think he is doing the interpolation by hand (2d interpolation of every slice). The link i posted is how to do it directly. Maybe @rayryeng was not aware of the function `tformarray`. Also I have to admit that I just used `tformarray` without looking up how it is working exactly. – nightlyop Jul 23 '14 at 15:26
  • @nightlyop - You are correct. I 2D interpolated every slice, then I used `interp1` to do a single interpolation in 3D. I think the link that you have posted is faster and more efficient than what I did. I wasn't aware of `tformarray` and that's a nifty function! Thanks for sharing! – rayryeng Jul 23 '14 at 16:50
  • Manolete - I'm still not sure what's wrong. Perhaps you should try one of the methods in the link @nightlyop suggested. Good luck – rayryeng Jul 23 '14 at 16:53

1 Answers1

4

You almost have it right. You need to define 3D Grid of co-ordinates. Creating single vectors is not the right way to do it. You can certainly use interp3 here. Try doing:

[X,Y,Z] = meshgrid(1:213, 1:100, 1:140);
Vq = interp3(M, X, Y, Z);

Note that I have swapped the row (100) and column (213) limits, as the first parameter progresses horizontally while the second parameter progresses vertically.

Also, by using interp3 in this fashion, we are assuming that the limits of X, Y and Z fall within 1:213, 1:100 and 1:140. Should you provide any values outside of these limits, you will get NaN. There are a couple of ways you can avoid this:

  1. Specify the spline flag at the end to allow for spline extrapolation
  2. If you want to resize the matrix (like if you were to resize an image), then there currently is no built-in method that can resize a 3D matrix this way. You'll have to write this yourself.

If you want to do Step #2, you can do the following.

First, you need to figure out the scale factors for each dimension. Basically this is the ratio of the output size of each dimension with the original input size.

After this, you create a 2D grid that has its limits bounded by the original size of the input matrix, but the size of this grid will be of the size of the output matrix. The scale factor is useful here because this effectively gives us what each value in the grid should be to interpolate. We would create new co-ordinates that went from 1 to the output size of each dimension, in increments of 1/scaleFactor. As an example, if we wanted to double the size of our matrix, this is a factor of 2. If we had X and Y co-ordinates that went from 1 to 3 and 1 to 3 respectively, the original grid would look like this:

X =            Y = 

1  2  3        1  1  1
1  2  3        2  2  2
1  2  3        3  3  3

To double this, this would simply be:

X =                         Y = 

1  1.5  2  2.5  3           1   1   1   1   1
1  1.5  2  2.5  3          1.5 1.5 1.5 1.5 1.5
1  1.5  2  2.5  3           2   2   2   2   2
1  1.5  2  2.5  3          2.5 2.5 2.5 2.5 2.5 
1  1.5  2  2.5  3           3   3   3   3   3

Note that this create an output grid of 5 x 5. To make this doubled as 6 x 6, you can do whatever you want, but for the sake of simplicity, just duplicate the last row and last column, and so:

X =                         Y = 

1  1.5  2  2.5  3  3         1   1   1   1   1   1
1  1.5  2  2.5  3  3        1.5 1.5 1.5 1.5 1.5 1.5
1  1.5  2  2.5  3  3         2   2   2   2   2   2
1  1.5  2  2.5  3  3        2.5 2.5 2.5 2.5 2.5 2.5
1  1.5  2  2.5  3  3         3   3   3   3   3   3
1  1.5  2  2.5  3  3         3   3   3   3   3   3

This defines our grid of 2D columns for resizing. Now it's a matter of resize in 3D. What we can do is interpolate in between the slices. We can easily do this using permute in MATLAB, and I'll show you how to do that later. As such, the basic algorithm is this:

  • Determine the output size of the output matrix you want
  • Determine the scale factors in each dimension
  • Create a 2D grid of interpolated access values for each dimension following the procedure above
  • For each 2D slice in your matrix, use interp2 to resize each slice to the output rows and columns using the above 2D grid.
  • After, use interp1 and permute to resize the third dimension.

Without further ado, here is the code to do this:

%// Specify output size of your matrix here
outputSize = [100 213 140];

%//Figure out size of original matrix
d = size(M);

%//Scaling coefficients
scaleCoeff = outputSize ./ d;

%//Indices of original slices in 3D
z = 1:d(3);

%//Output slice indices in 3D
zi=1:1/scaleCoeff(3):d(3);

%//Create gridded interpolated co-ordinates for 1 slice
[X,Y] = meshgrid(1:1/scaleCoeff(2):d(2), 1:1/scaleCoeff(1):d(1));

%//We simply duplicate the last rows and last columns of the grid if
%//by doing meshgrid, we don't get exactly the output size we want
%//This is due to round off when perform 1/scaleCoeff(2) or
%//1/scaleCoeff(1).  We would be off by 1.
if size(X,1) ~= outputSize(1)
    X(end+1,:) = X(end,:);
    Y(end+1,:) = Y(end,:);
end
if size(X,2) ~= outputSize(2)
    X(:,end+1) = X(:,end);
    Y(:,end+1) = X(:,end);
end

%//For each slice...
M2D = zeros(outputSize(1), outputSize(2), d(3));
for ind = z
    %//Interpolate each slice via interp2
    M2D(:,:,ind) = interp2(M(:,:,ind), X, Y);
end

%//Now interpolate in 3D
MFinal = permute(interp1(z,permute(M2D,[3 1 2]),zi),[2 3 1]);

%//If the number of output slices don't match after we interpolate in 3D, we
%//just duplicate the last slice again
if size(MFinal,3) ~= outputSize(3)
    MFinal(:,:,end+1) = MFinal(:,:,end);
end

MFinal would be your final interpolated / resized 3D matrix. The key method to interpolate in 3D is the permute method. What this will do is that for each value of z, we will generate a 2D slice of values. As such, if we had a slice at z = 1 and one at z = 2, if we wanted to find what the 2D grid of values was at slice z = 1.5, this will generate a 2D slice that creates these interpolated values using information between z = 1 and z = 2. We do the first call of permute to do this, then another permute call to undo our permutation and get the original dimensions back.

rayryeng
  • 102,964
  • 22
  • 184
  • 193
  • That's really excellent! Should I always swap the rows and columns limits when using `interp3`? – Manolete Jul 15 '14 at 19:52
  • 1
    @Manolete - Yes. The first parameter always progresses horizontally while the second parameter progresses vertically. Also, if you want to extrapolate, use the `spline` flag at the end. I'll modify my answer. – rayryeng Jul 15 '14 at 19:53
  • I've check the values of the matrix obtained and it seems that there are two type of values, 0s and NaNs. I have probably done something wrong – Manolete Jul 15 '14 at 20:05
  • @Manolete - Don't worry. I will fix my answer – rayryeng Jul 15 '14 at 20:17
  • @_rayryeng : This is an excellent explanation, thank you very much! I have followed your steps, however in MFinal I can still see some NaNs. Could I just convert those NaNs to 0s? – Manolete Jul 16 '14 at 08:55
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/57416/discussion-between-manolete-and-rayryeng). – Manolete Jul 16 '14 at 13:59
  • @Manolete - That shouldn't happen... I tested this on my own end and it is generating what I have expected for my tests. Can you tell me what the nature of `M` is? There could be a case where when `interp` is being called (1D or 3D), it could be producing `NaN`. Maybe specify the `spline` flag for each `interp` call and see if that fixes it. Either way, you're generating weird behaviour. I made it specifically so that it would not have to extrapolate. – rayryeng Jul 16 '14 at 14:27
  • @_rayryeng : `M` represents the portion of energy in a 3D space, i.e. the energy delivered per voxel. The values range from 0 to 268435461. – Manolete Jul 17 '14 at 09:39
  • @_rayryeng : I've tried using `spline` and now I cannot see any `NaN` value. Having said that, I cannot understand why `min(min(min(MFinal))) = -5.2212e+07` when `min(min(min(M))) = 0`. As you said, there is something weird here, but the matrix `M` is not really a special one it just represents a 3D space, a part of the body, with the energy dose delivered to those voxels. – Manolete Jul 17 '14 at 10:59