Vectorising a Meshgrid Coordinate Construction

I have three coordinate arrays (file attached) of size [M x 2], where M is the number of elements and the 2 represents different points in a 1D system. To find the 3D meshed coordinates upon which I perform calculations I use code of the form:

for ele = 1:size(X_co,1) %Loop Over Each Element
    [X_cube,Y_cube,Z_cube] = meshgrid(X_co(ele,:),Y_co(ele,:),Z_co(ele,:)); %Construct Meshgrid of size [2x2x2] for an Element
    X_save(:,:,:,element) = X_cube; %Save Coordinates
    Y_save(:,:,:,element) = Y_cube;
    Z_save(:,:,:,element) = Z_cube;    
end

The solution then is of size [2x2x2xM]. My best attempt to vectorise this was:

[X_cube,Y_cube,Z_cube] = meshgrid(X_co,Y_co,Z_co);  
X_save = reshape(X_cube,2,2,2,[]); Y_save = reshape(Y_cube,2,2,2,[]);  Z_save = reshape(Z_cube,2,2,2,[]); 

However this solution is of size [2x2x2x(M^3)] which, when plotting both solutions has many many many repeated points so doesn't actually work for my application.

Can this problem be vectorised, or is a for loop the best way to construct this coordinate array?

[KSSV]

Are you looking for ndgrid?

[jessupj]

ndgrid is the way to go here i think. i think it can be done using meshgrid along with some resize + repmat statements.

[ADSW121365]

https://uk.mathworks.com/matlabcentral/answers/588052-vectorising-a-meshgrid-coordinate-construction#comment_994471

This maybe niave implementation, but:

[XO,YO,ZO] = ndgrid(X_co,Y_co,Z_co)

seems to give the same solution as:

[X_cube,Y_cube,Z_cube] = meshgrid(X_co,Y_co,Z_co);  
X_save = reshape(X_cube,2,2,2,[]); Y_save = reshape(Y_cube,2,2,2,[]);  Z_save = reshape(Z_cube,2,2,2,[]); 

which isn't what I'm looking for. Could you guys elaborate a little more please? (Some example arrays are attached to the main post)

Answer 1

[Bruno Luong]

One way

[~,nx] = size(X_co);
[~,ny] = size(Y_co);
[m,nz] = size(Z_co);
[NX,NY,NZ] = meshgrid(1:nx,1:ny,1:nz);
sz = [size(NX),m];
% if you know nx, ny, nz are 2 you can replace the 5 above commands by
% [NX,NY,NZ] = meshgrid(1:2);
% sz = [2,2,2,m];

X_save = X_co.';
Y_save = Y_co.';
Z_save = Z_co.';
X_save = reshape(X_save(NX,:),sz);
Y_save = reshape(Y_save(NY,:),sz);
Z_save = reshape(Z_save(NZ,:),sz);

Note: I honestly prefer your for-loop if you add a preallocation

[ADSW121365]

I like this answer a lot, but what makes you prefer my for-loop approach?

In the main problem, nx = ny = nz but M maybe significantly larger, this is then embedded in two other for-loops. The idea here was this is the longest for-loop by some significant margin, so vectorising this loop would be the best choice to optimise the code for speed. Is my solution just nicer/more compact for this specific problem, or is there a speed advantage to using the for-loop?

[ADSW121365]

I think I can probably vectorise the other two for loops if I try hard enough and get creative with memory management but I'm not sure the conceptual nightmare trying to set that problem up is worth the extra effort.

[Bruno Luong]

The for-loop code is more readable and I would guess is not slow at all if you properly allocate the output array.

In two ways of indexing (resp for-loop and vectorization)

A(:,:,:,elem) % and 
A(NX,:), 

MATLAB prefers the first since it's access contiguous memory in the first case and scattered in the second case, which is not good for speed on modern CPU.

However the for-loop needs to call MESHGRID many time, which in turns is a penalty compared to the vectorization code.

But really up to you.

[ADSW121365]

So by properly allocate you mean initialise the arrays as zeros before the for loop starts?

I'm suprised to learn the first is actually faster but that does make a lot of sense when I think a bit harder about this. Thanks, I really appreciate the thoughts!

[Bruno Luong]

"So by properly allocate you mean initialise the arrays as zeros before the for loop starts?"

Yes intiialize the 3 zeros arrays with SZ from with my code.

[ADSW121365]

Thanks, I was just checking my indexing wasn't the allocation issue.