interp3 gives error in some cases and accept some cases with the same size of meshgrid !

I have strange problem in my program that could not be solved , and I would like to have explanation for it please

I am trying to generate a specific 3D flow and use it in interp3 function for warping process. (in 3D MRI images)

mag=5;
U=linspace(-mag,mag,size_x);
V=linspace(-mag,mag,size_y);
W=linspace(-mag,mag,size_z);
X=1:size_x; %256
Y=1:size_y; %256
Z=1:size_z; %75

if I generate the flow using ndgrid,

[u,v,w]=ndgrid(U,V,W);
[x,y,z]=ndgrid(X,Y,Z);

I got the correct flow that I want

but I cannot use it because I got this error in interp3 function

Error using interp3 (line 146)
Input grid is not a valid MESHGRID.

however, if I generate the flow using meshgrid with this order of arguments

 [w,u,v]=meshgrid(U,V,W);
 [z,x,y]=meshgrid(X,Y,Z);

I got wrong flow , but interp3 do not give any error !!!

In addition, I tried meshgrid with this order ( logic order of the arguments)

[u,v,w]=meshgrid(U,V,W);
[x,y,z]=meshgrid(X,Y,Z);

that give me wrong flow (but organized flow) and error in interp3 function !!!!

it is strange because in all cases all the meshgrids and ndgrid has the same size ( 256 256 75 ) in x,y,z,u,v, w

I need explanation why it works in some cases and other not ( same sizes)? and how can I generate the correct flow that can be used with interp3 without problem

appreciate your help in advanced

Answer 1

[Kelly Kearney]

The difference between ndgrid and meshgrid is how they treat the first two dimensions.

size_x = 3;
size_y = 4;
size_z = 2;

mag=5;
U=linspace(-mag,mag,size_x);
V=linspace(-mag,mag,size_y);
W=linspace(-mag,mag,size_z);
X=1:size_x; %256
Y=1:size_y; %256
Z=1:size_z; %75

ndgrid varies things in order (vector one along first dimension, vector 2 along the second dimension, etc):

[u1,v1,w1]=ndgrid(U,V,W);
[x1,y1,z1]=ndgrid(X,Y,Z);

>> x1

x1(:,:,1) =

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

x1(:,:,2) =

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

>> y1

y1(:,:,1) =

       1     2     3     4
       1     2     3     4
       1     2     3     4

y1(:,:,2) =

       1     2     3     4
       1     2     3     4
       1     2     3     4

while meshgrid flips the first two dimensions (to sort of match x/y conventions when plotting):

[u2,v2,w2]=meshgrid(U,V,W);
[x2,y2,z2]=meshgrid(X,Y,Z);

>> x2

x2(:,:,1) =

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

x2(:,:,2) =

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

>> y2

y2(:,:,1) =

       1     1     1
       2     2     2
       3     3     3
       4     4     4

y2(:,:,2) =

       1     1     1
       2     2     2
       3     3     3
       4     4     4

Most plotting functions, such as quiver3, will accept either of these, as long as the input matrices all have the same input order. And when dealing with more than two dimensions, I definitely prefer to use ndgrid. However, several of the interpolation functions (namely, griddedInterpolant, which underlies interp3), very annoyingly insist on using the meshgrid-style input only.

You can either use the meshgrid setup (but make sure to keep your input and output variables in the same order!), or use ndgrid for setup and then permute the values when interpolating:

x3 = rand(10,1)*(size_x-1) + 1;
y3 = rand(10,1)*(size_y-1) + 1;
z3 = rand(10,1)*(size_z-1) + 1;

u3 = interp3(permute(x1,[2 1 3]), ...
             permute(y1,[2 1 3]), ...
             permute(z1,[2 1 3]), ...
             permute(u1,[2 1 3]), x3,y3,z3);
v3 = interp3(permute(x1,[2 1 3]), ...
             permute(y1,[2 1 3]), ...
             permute(z1,[2 1 3]), ...
             permute(v1,[2 1 3]), x3,y3,z3);
w3 = interp3(permute(x1,[2 1 3]), ...
             permute(y1,[2 1 3]), ...
             permute(z1,[2 1 3]), ...
             permute(w1,[2 1 3]), x3,y3,z3);

u4 = interp3(x2,y2,z2,u2, x3,y3,z3);
v4 = interp3(x2,y2,z2,v2, x3,y3,z3);
z4 = interp3(x2,y2,z2,w2, x3,y3,z3);

>> isequal(u3,u4)

ans =

     1

[Seereen]

appreciate your answer , I tried your code and it's work without error. I already try permute function to change the order, but I still get the error ! that why I confused. if I use ndgrid as I mention in the question, I got 256 256 75 size, then If I use permute I got the same size and I still get the same error in interp3 function !

[Seereen]

TO make it easy I just read about interpn function that accept ndgrid format :D

[Kelly Kearney]

Keep in mind that it's not the size of the array that interp3 gets upset about, but rather the order of the plaid in the first two dimensions.

But going with interpn is probably a good solution.