Why such a fuss with ndgrid and meshgrid issues?

Folks,

In the world of artificial intelligence and more intelligent programming, I would like to see more intelligence put into place with how function like griddata, imagesc, pcolor, and surf interpret their inputs. There is an unfortunate degree of confusion that has persisted since the early 2000's at least with how to interpret arrays in terms of Cartesian axes. I am troubled today with trying to understand how griddata is working on input that I setup using ndgrid. I appear to have something wrong, and earth layers are not dipping in the direction that I expect. The documentation says griddata can work with either meshgrid or griddata input, which makes my head spin given my understanding of both meshgrid and ndgrid. What is really troubling is that the Mathworks documentation routinely contains only simple examples where the person that created the documentation made "quicky" non-realistic examples where they simply created input arrays that had exactly the same x and y dimensions, which does not help anyone in my situation who is dealing with realistic inputs that do not have the same dimension lengths. I request that future documentation efforts spend a little extra time making non-symmetric input x and y arrays as examples.

[John D'Errico]

Why are you asking this on Answers? Answers is not official tech support. It is not a direct line to MathWorks. Answers is a volunteer forum. Sometimes people from TMW look in. But there is no asurance.

If you have a request, then send it in to support.

https://www.mathworks.com/support/contact_us.html?s_tid=hc_trail

[Stephen23]

"...the person that created the documentation made "quicky" non-realistic examples where they simply created input arrays that had exactly the same x and y dimensions..."

I agree that square matrices are unhelpful examples for this topic. At the bottom of every documentation page is a field for rating the documentation: give it a rating and write the precise reason why. The documentation authors are unlikely to read random forum threads, but they definitely get the rating feedback. Use it!

[Kristoffer Walker]

I'll try to ask some more clarifying questions. I appreciate the attention this post is getting.

In obliqueslice, it requires input 3d array V. Are the assumed dimensions of V (nX, nY, nZ) or (nY, nX, nZ)? I can see in the documentation that the point and normal vectors are in (x, y, z) format. In an experiment, when I have V with dimensions of (130, 150, 567), when I use a normal of [0 1 0], I get an output slice of dimensions (567, 150), which is unexpected and not an XZ slice. I have tried reading the documentation regarding where the "upper left point" maps to, but I could not understand it.

Furthermore, when I assume V should be (nY, nX, nZ) dimensions, the output slice has features in it that do not align with features in the input V array, which was created using ndgrid and griddata. The creation and 3d visualization of V appears consistent with expectations when I plot squeezed xy slices using imagesc or pcolor, but the obliqueslice is not and I have to do weird things to try and make it align with expectations (like reversing the coordinate vectors when I plot the slice), yet the complete alignment has not been achieved yet. Ultimately I want to use obliqueslice so I can get oblique angles. I have tried "slice" as well, but also had problems with it too.

One solution that I tried was using "slice". I already have an array V that took a long time to grid. I do not want to have to regrid is using a meshgrid format. So I thought I could use reshape to convert my arrays from ndgrid to meshgrid format, because that is what splice seems to require. Below shows success using meshgrid format, but an issue arises when using ndgrid format that I do not understand how to fix without using meshgrid format.

clear all; close all

% Using meshgrid approach for non-square slice example
x = [-2:0.2:2]; y = [-5:0.1:5]; z = [-25:2:25];

[X,Y,Z] = meshgrid(x,y,z);

V = X.*exp(-X.^2-Y.^2-Z.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(X,Y,Z,V,xslice,yslice,zslice)


% Using ndgrid approach for non-square slice example
x = [-2:0.2:2]; y = [-5:0.1:5]; z = [-25:2:25];

[X,Y,Z] = ndgrid(x,y,z);

nX = length(x); nY = length(y); nZ = length(z);
xg = reshape(X,[nY,nX,nZ]);
yg = reshape(Y,[nY,nX,nZ]);
zg = reshape(Z,[nY,nX,nZ]);

V = xg.*exp(-xg.^2-yg.^2-zg.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(xg,yg,zg,V,xslice,yslice,zslice)

[Voss]

Don't use reshape for that; use permute:

% Using meshgrid approach for non-square slice example
x = -2:0.2:2; y = -5:0.1:5; z = -25:2:25;

[X,Y,Z] = meshgrid(x,y,z);

V = X.*exp(-X.^2-Y.^2-Z.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(X,Y,Z,V,xslice,yslice,zslice)

% Using ndgrid approach for non-square slice example
x = -2:0.2:2; y = -5:0.1:5; z = -25:2:25;

[X,Y,Z] = ndgrid(x,y,z);

nX = length(x); nY = length(y); nZ = length(z);
xg = permute(X,[2,1,3]);
yg = permute(Y,[2,1,3]);
zg = permute(Z,[2,1,3]);

V = xg.*exp(-xg.^2-yg.^2-zg.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(xg,yg,zg,V,xslice,yslice,zslice)

[Kristoffer Walker]

Thank you, Voss. That fixed my issue with both splice and obliquesplice.

Answer 1

[Paul]

Hi Kristoffer,

The inputs to griddata that can be in either meshgrid or ndgrid format are the query points in Xq and Yq. I think that their specific format doesn't matter because griddata processes each pair ( Xq(i,j),Yq(i,j) ) where the output vq(i,j) = f(Xq(i,j) , Yq(i,j)) is the interpolation based on Delaunay triangulation of scattered sample points and sample values (in contrast to, for example, interp2 that requires the sample points and values to be in meshgrid format).

What would matter is what's done with Xq, Yq, and vq after the call to gridddata. For example, surf prefers (requires?) the input data to be in meshgrid format, so a call to griddata should use query points in meshgrid format if to be followed by a call to surf with the results.

If having a problem with griddata with query points in ndgrid format, you're likely to get more help by posting a simple code example to show what the issue is.

Edit: As shown by @Star Strider in this Answer, surf is agnostic to whether or not the the input is in ndgrid/meshgrid format.

Edit: Use of surf with meshgrid/ndgrid was also explained by @Voss in this thread that I had forgotten about.

I think my confusion with surf (et. al.?) might be the result of a long-ago behavior change that still hasn't solidified with me.

In surf (2019a), the doc stated "To create a matrix for arbitrary domains, use the meshgrid function." No ambiguity there.

In surf (2019b) the doc states: "You can use the meshgrid function to create X and Y matrices." ("can" not "must") which suggests something changed from 2019a to 2019b, though there is no mention of such in the 2019b release notes.

[Kristoffer Walker]

Hi Paul. Please see my comment to my original post above.

Answer 2

[David Goodmanson]

Hi Kristoffer,

I totally agree that Mathworks should use nonsquare examples in the documentation for meshgrid, ndgrid, surf etc. It just increases awareness. And comparing [x y] = meshgrid(a,b) and [x y] = ndhgrid(a,b),

***** the meshgrid x matrix and ndgrid x matrix are transposes of each other. Same for y.

I think you meant to say that griddata works with either meshgrid or ndgrid input, and that's true. Griddata is agnostic on the topic, assuming of course that either meshgrid or ndgrid are used consistently throughout.

[Kristoffer Walker]

Hi David. Please see my comment to my original question above.

Answer 3

[Star Strider]

I’m not certain what the problem is.

I generally prefer scatteredInterpolant to griddata.

The choce of meshgrid or ndgrid does not seem to be an issue, at least in this example —

xv = linspace(-2, 2, 15);
yv = linspace(-1, 4, 19);
[X1,Y1] = meshgrid(xv, yv);
Z1 = exp(-X1.^2) .* exp(-(Y1-2.5).^2*1.5)*2.5;

figure
surf(X1, Y1, Z1)
colormap(turbo)
axis('equal')
title('Original')

F = scatteredInterpolant(X1(:), Y1(:), Z1(:));

xv2 = linspace(-2, 2, 40);
yv2 = linspace(-1, 4, 30);
[X2,Y2] = ndgrid(xv2,yv2);
Z2 = F(X2,Y2);

figure
surf(X2,Y2,Z2)
colormap(turbo)
axis('equal')
title('scatteredInterpolant — ndgrid')

[X2,Y2] = meshgrid(xv2,yv2);
Z2 = F(X2,Y2);

figure
surf(X2,Y2,Z2)
colormap(turbo)
axis('equal')
title('scatteredInterpolant — meshgrid')

Do either of these deal with the issue you’re describing?

.

[Kristoffer Walker]

Hi Star. Please see my comment to my original post above.

[Star Strider]

We would need a more concrete example, specifically with code and data. It would be best to paste the code as text in a Comment to your original question, or edit your original questtion to include it, and attach the data using the ‘paperclip’ icon. It may be necessary for you first to use the zip function to create a .zip file of the data.

I have essentially no experience with the obliqueslice function, since I rarely do volumetric image processing.

[Kristoffer Walker]

Thanks Star. I will try to recreate the problem with a different data set. It is very hard because I cannot share any data.

[Kristoffer Walker]

One solution that I tried was using "slice". I already have an array V that took a long time to grid. I do not want to have to regrid is using a meshgrid format. So I thought I could use reshape to convert my arrays from ndgrid to meshgrid format, because that is what splice seems to require. Below shows success using meshgrid format, but an issue arises when using ndgrid format that I do not understand how to fix without using meshgrid format.

clear all; close all

% Using meshgrid approach for non-square slice example
x = [-2:0.2:2]; y = [-5:0.1:5]; z = [-25:2:25];

[X,Y,Z] = meshgrid(x,y,z);

V = X.*exp(-X.^2-Y.^2-Z.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(X,Y,Z,V,xslice,yslice,zslice)


% Using ndgrid approach for non-square slice example
x = [-2:0.2:2]; y = [-5:0.1:5]; z = [-25:2:25];

[X,Y,Z] = ndgrid(x,y,z);

nX = length(x); nY = length(y); nZ = length(z);
xg = reshape(X,[nY,nX,nZ]);
yg = reshape(Y,[nY,nX,nZ]);
zg = reshape(Z,[nY,nX,nZ]);

V = xg.*exp(-xg.^2-yg.^2-zg.^2);

xslice = [0.8]; yslice = []; zslice = [];
figure; slice(xg,yg,zg,V,xslice,yslice,zslice)

[BANAN]

[x, y] = meshgrid(-2:0.2:2);
z = x .* exp(-x.^2 - y.^2);  

Answer 4

[William Rose]

@Kristoffer Walker, I definitely agree with you. I use inputs that are non-square when i use surf() and related routines. Thus I use vectors X, Y with different lengths in surf(X,Y,Z), etc. If the inputs are square, I may flip things around by mistake, and I might not know it, because there's no error message.

I realize @John D'Errico is right that this forum is not the Suggestion Box for matlab. Try here.

[Kristoffer Walker]

Hi Bill. Please see my comment to my original post above.

Answer 5

[Kristoffer Walker]

Thank you everyone for your input. I will surely look back on this post in the future for the multiple, helpful contributions.