How to put 3D array data into mesh grid and find grid cell index of a cube enclosing a given point?

I have three, 3D arrays containing x, y, z positions (coordinates) inside a 3D grid space.

I am given a point (x,y,z) which is somewhere inside the grid space. I want to find indices (i,j,k) of the grid cell ("cube") enclosing that point.

My question is how to find the grid cell indices, and 8 vertices of that cube.

I have "found" them, but I want to confirm with someone, as I have not worked with 3D grids before and am not too confident.

***

The orientation of the entire grid space is such that:

(i)Zero (0) is in the middle of all axes.

(ii)For X-axis, positive is away from the viewer (towards the screen), negative is toward the viewer.

(iii)Y-axis runs from right to the left, right side is negative.

(iv)Z-axis is vertical, upwards is positive.

For a small test case, let's assume that all the six arrays are of dimension 5x5x5.

If I look inside 'X3D' array (the one that has x-coordinate data), it looks like this:

val (:,:,4) and val (:,:,5) are exactly the same.

For 'Y3D' array (the one that has y-coordinate data), it looks like this:

val (:,:,3), val (:,:,4) and val (:,:,5) are exactly the same.

Finally, for 'Z3D':

val (:,:,4) and val (:,:,5) are same values like val (:,:,2) and val (:,:,1) , but with flipped sign.

The following are steps I used to find the grid cell indices of the cube that encloses a given point.

(1)Read data

load X3D.mat;load Y3D.mat; load Z3D.mat; %coordinate data

(2) Calculate length of each axis, and grid cell width. Calculation for just one axis will do, as it is a cube.

%get the first row of X3D array
first_row = X3D(1,:,1) % returns -1 -0.5  0  0.5 1 (trailing zeros removed for clarity)
axis_half_length = abs(first_row(:,1)); % I know that the first item (1st col.) is the "half length", pls. read axis orientation

%axis length (same for all axes)
axis_length = 2 * axis_half_length;

(3) Calculate grid cell width

no_of_grid_points = size(X3D, 1); % size of the 1st dim, it will be 5 in this test case
grid_cell_width = axis_length / (no_of_grid_points - 1);

(4) Define the grid parameters based on the cell width

%Define the grid parameters based on the cell width
x = -axis_half_length:grid_cell_width:axis_half_length;
y = -axis_half_length:grid_cell_width:axis_half_length;
z = -axis_half_length:grid_cell_width:axis_half_length;

(5) Create meshgrid

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

(6) given position coordinates for the point

px = 0.52; py = -0.22; pz = 0.29;

(7) Plotting the grid

% Plotting the grid
figure;
hold on;
% Plot the grid lines in 3D
for xi = x
    for yi = y
        plot3([xi, xi], [yi, yi], [min(z), max(z)], 'k');
    end
    for zi = z
        plot3([xi, xi], [min(y), max(y)], [zi, zi], 'k');
    end
end
for yi = y
    for zi = z
        plot3([min(x), max(x)], [yi, yi], [zi, zi], 'k');
    end
end

(8) Highlight the point

plot3(px, py, pz, 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b');
text(px, py, pz, sprintf('(%0.2f, %0.2f, %0.2f)', px, py, pz), ...
    'FontSize', 12, 'Color', 'blue', 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left');

(9)Annotate the grid points using the coordinates (Fig. 1)

% Annotate all grid points (using coordinates)
for i = 1:length(x)  
    for j = 1:length(y)
        for k = 1:length(z)
            % Calculate coordinates
            x_pos = x(i);
            y_pos = y(j);
            z_pos = z(k);

            %Annotate with coordinates
            text(x_pos, y_pos, z_pos, ...  
            sprintf('(%g, %g, %g)', x_pos, y_pos, z_pos), ...
            'FontSize', 10, 'Color', 'black');
        end
    end
end

Here are the grid points, using coordinates fom the data file.

The blue dot is a point with coordinates (0.52,-0.22,0.29)

Fig. 1

As can be seen above, each axis of the grid space is 2 units long, from -1 to +1. There are 5 grid points (-1, -0.5, 0, +0.5, +1) in each direction.

Each grid cell width will be (2 *1) / (5-1) = 0.5 units.

Apologies for misalignment of the cubes and the axes.

(10) Annotate the entire grid points using the grid cell indices (Fig. 2) (instead of the coordinates)

% Annotate all grid points using grid cell indices
for i = 1:length(x)  
    for j = 1:length(y)
        for k = 1:length(z)
            
            % coordinates
            x_pos = x(i);
            y_pos = y(j);
            z_pos = z(k);
            
            % Calculate grid indices
            grid_i = i - 1; % Shift to origin (0,0,0)
            grid_j = j - 1;
            grid_k = k - 1;
            
            % Annotate with grid indices
            text(x_pos, y_pos, z_pos, ...
                sprintf('{(%d, %d, %d)}', grid_i, grid_j, grid_k), ...
                'FontSize', 8, 'Color', 'black');
        end
    end
end

Here is the same plot (as Fig. 1), but using grid cell indices. Note the origin (0,0,0).

Fig. 2

(11) Calculate grid cell indices for cube enclosing the point (px,py,pz) (Fig. 3, below)

% I assumed the origin is at the bottom right, front
i = floor((px + axis_half_length) / grid_cell_width) + 1;
j = floor((py + axis_half_length) / grid_cell_width) + 1;
k = floor((pz + axis_half_length) / grid_cell_width) + 1;

% Vertices of the enclosing cube
%https://blogs.mathworks.com/cleve/2022/01/03/anatomy-of-a-cube/
%vert = [0 0 0;  1 0 0;  0 1 0;  1 1 0; bottom v1, v2, v3, v4
%        0 0 1;  1 0 1;  0 1 1;  1 1 1]; top v5, v6, v7, v8
vertices = [
    x(i), y(j), z(k);         % Vertex 1: bottom right, front %000
    x(i+1), y(j), z(k);       % Vertex 2: bottom right, back  %100
    x(i), y(j+1), z(k);       % Vertex 3: bottom left, front %010
    x(i+1), y(j+1), z(k);     % Vertex 4: bottom left, back %110
    x(i), y(j), z(k+1);       % Vertex 5: top right, front %001
    x(i+1), y(j), z(k+1);     % Vertex 6: top right, back %101
    x(i), y(j+1), z(k+1);     % Vertex 7: top left, front %011
    x(i+1), y(j+1), z(k+1);   % Vertex 8: top left, back %111
    ];

(12) Draw the cube that contains the point

% Plot the enclosing cube
edges = [
    1, 2; 1, 3; 1, 5;
    2, 4; 2, 6;
    3, 4; 3, 7;
    4, 8;
    5, 6; 5, 7;
    6, 8;
    7, 8
    ];
for e = 1:size(edges, 1)
    plot3(vertices(edges(e, :), 1), vertices(edges(e, :), 2), vertices(edges(e, :), 3), 'r', 'LineWidth', 2);
end


% Annotate the cube's vertices with grid indices
vertex_labels = {
    sprintf('[%d, %d, %d]', i, j, k);
    sprintf('[%d, %d, %d]', i+1, j, k);
    sprintf('[%d, %d, %d]', i, j+1, k);
    sprintf('[%d, %d, %d]', i+1, j+1, k);
    sprintf('[%d, %d, %d]', i, j, k+1);
    sprintf('[%d, %d, %d]', i+1, j, k+1);
    sprintf('[%d, %d, %d]', i, j+1, k+1);
    sprintf('[%d, %d, %d]', i+1, j+1, k+1)
};

for v = 1:size(vertices, 1)
    text(vertices(v, 1), vertices(v, 2), vertices(v, 3), vertex_labels{v}, ...
        'FontSize', 10, 'Color', 'black', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end

Fig. 3

(13) My real data has many more grids, but the axes orientation is the same.

My question is: did I get those indices correct?

I mean is it how one calculates the grid cell indices?

I thought I did because, for example, if I line up the X3D or Y3D values, I would find a correct "cell" to insert the px(0.52), py(-0.22), pz(0.29) values, and they would be cell indices (4, 2, 3), as shown below in an illustration. (Only for 2 directions- X and Y- are shown.)

Fig. 4

Is it that straightforward (as implemented in (11) and as verified in Fig. 4)? Or is there more to it, and I got it completely wrong and/or only partially right?

Thank you for reading this long post.

Answer 1

[Matt J]

It would be much simpler and more efficient to skip the meshgrid altogether and use discretize.

i=discretize(px,x);
j=discretize(py,y);
k=discretize(pz,z);

which can also be vectorized across many sets of (px,py,pz).

Also, if you are doing this for the purposes of interpolation, you should not be doing the index search at all. You should just be using griddedInterpolant or interp3.

[K.]

Thank you for your comments, MJ!

The discretize method is also giving me the same i,j,k, so I can validate my result.

I don't have that many sets of points, so I will leave my code as it is for now.

You are right about the purpose! Yes, I am trying to find the magnetic field values at that point using interpolation. I had prepared the second half of the question also, but I felt it was getting longer, so I decided to keep the index finding part only.

Once I find the indices of the 8 vertices of the cube that contains the point, I can extract field values at those points from another set of 3D arrays (BX3D, BY3D, BZ3D). Then, I want to use a scheme described by Tristan (trilinear interpolation/volume weighting). Now it seems that I might have some other choices as well (your suggestions).

[Matt J]

You are quite welcome. If your question has been addressed, though, please Accept-click the answer.