Extracting coordinates from meshgrid data

Question

I have a cubic grid as shown in the picture below.

I would like to list the vertices of each sub-cube, so I would end up with a nested list of sub-cubes with their corresponding list of vertices.

My initial attempt was to use a generator,

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

dims = [9,9,9]    
spacer = 3
subBoxCoords = np.array([(x, y, z) for x in range(0, dims[0], spacer) for y in range(0, dims[1], spacer) for z in range(0, dims[2], spacer)])
ax.scatter(subBoxCoords[:,0], subBoxCoords[:,1], subBoxCoords[:,2], c='k', marker='o')

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

This does give me the desired shape but coordinates are ordered in a manner that vertices extraction of the sub-boxes is not straight forward. Also I would like to generalize this to boxes of arbitrary dimension so hard coding in intervals is not a solution.

So, then I thought I would use meshgrid,

nx,ny, nz = (3,3,3)
x = np.linspace(0, 10, nx)
y = np.linspace(0, 10, ny)
z = np.linspace(0, 10, nz)
xv, yv, zv = np.meshgrid(x, y, z, indexing='xy')

ax.scatter(xv, yv, zv, c='g', marker='^')

This appears to be a very powerful way to achieve what I want but I am getting confused. Is there a direct way access vertices in the meshgrid in the manner vertex(x,y,z)? Or even a straight forward way to extract sub-cubes?

It seems to me that the solution is tantalizingly close but I just cant grasp it!

Answer 1

meshgrid is probably what you need, but the shape of the array returned by meshgrid is where it gets confusing. Meshgrid returns three coordinate arrays, all the same shape. The shape of each of xv, yv, zv is (len(x), len(y), len(z)). So, to extract the coordinate at the corner (0, 2, 1), you would write xv[0, 2, 1], yv[0, 2, 1], zv[0, 2, 1]

To extract all of the subcubes' corners' coordinates, it helps to observe that, because of the way the arrays returned by meshgrid are ordered sequentially, xv[:-1, :-1, :-1] returns only the x-coordinates of the near-left-bottom corners of each subcube. Likewise, xv[1:, 1:, 1:] returns the far-right-top corners of each subcube. The other six corners are given by the other six combinations of the slices :-1 and 1: (xv[:-1, 1:, :-1] gives the far-left-top corner, for example).

So, iterate through all eight combinations of :-1 and 1: to get eight parallel arrays of three parallel arrays of x, y, z coordinates for the eight corners of all len(x)-1 * len(y-1) * len(z-1) subcubes. (If you need your subcube corner coordinate arrays to be in a particular shape or axis order, or if you want to use a single index to specify the subcube rather than three, use rollaxis, swapaxis and shape as needed.)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import itertools

nx, ny, nz = (3,3,3)
x = np.linspace(0, 10, nx)
y = np.linspace(0, 10, ny)
z = np.linspace(0, 10, nz)
xv, yv, zv = np.meshgrid(x, y, z, indexing='xy')


slices = slice(None, -1), slice(1, None)
cornerSlices = list(itertools.product(slices, slices, slices))
corners = np.array([(xv[s], yv[s], zv[s]) for s in cornerSlices])

# The shape of `corners` is `(len(cornerSlices), 3, len(x-1), len(y-1), len(z-1)`
# The axes of `corners` represent, in the same order: the corner index; the cartesian 
# coordinate axis (the index into [x, y, z]); the x, y, and z indexes of the subcube. 

# Plot the first subcube (subcube 0, 0, 0)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

subcube = corners[:, :, 0, 0, 0]
subcubeX = subcube [:, 0]
subcubeY = subcube [:, 1]
subcubeZ = subcube [:, 2]

ax.scatter(subcubeX , subcubeY , subcubeZ , c='g', marker='^')

There's invariably a way to get the indexes into xv, yv, zv instead of getting the values, since the values are duplicated quite a few times in the corners array. It would involve slicing arrays of indexes into xv, yv, zv instead of slicing the arrays themselves. My head is already spinning after getting this far into the ndarray voodoo, so I'll leave that as an exercise.

Answer 2

If I understand the question correctly -

x = np.linspace(0,10,3)
xs = np.array(np.meshgrid(x,x,x))
points = np.transpose(xs.reshape(3,-1))
points