How can I rotate a 3d array (nxnxn) by x degrees around x, y, and z axes?

I have a 3d array created with numpy, and I was wondering how I can rotate it by a custom angle, not just the rot90 function that numpy has. Can anyone help?

The 3d matrix represents an image (such as a cube, or some other shape) ie

0: 1 1 1 1 1 1 1 1 1: 1 1 1 1 2: 1 1 1 1 1 1 1 1 

EDIT: Moved solution to answer

4 Answers

Have a look at the scipy.ndimage.interpolation.rotate function.

The reason this is in scipy and not in numpy is that rotating an image 90 degrees is done by just chaning the indices of the array. However, if you want to rotate an image by some arbitrary degrees you have to deal with interpolation, which adds a whole new layer of complexity to the problem. This is because all the pixels in the original image "perfectly lines up with" pixels in the rotated image when you rotate it by a factor of 90 degrees. This is not the case in general when you rotate an image.

After some trial and error I came up with some code for my purposes (0 means empty in the array, another number will mean a filled voxel.

def rotate(self, deg_angle, axis): d = len(self.matrix) h = len(self.matrix[0]) w = len(self.matrix[0][0]) min_new_x = 0 max_new_x = 0 min_new_y = 0 max_new_y = 0 min_new_z = 0 max_new_z = 0 new_coords = [] angle = radians(deg_angle) for z in range(d): for y in range(h): for x in range(w): new_x = None new_y = None new_z = None if axis == "x": new_x = int(round(x)) new_y = int(round(y*cos(angle) - z*sin(angle))) new_z = int(round(y*sin(angle) + z*cos(angle))) elif axis == "y": new_x = int(round(z*sin(angle) + x*cos(angle))) new_y = int(round(y)) new_z = int(round(z*cos(angle) - x*sin(angle))) elif axis == "z": new_x = int(round(x*cos(angle) - y*sin(angle))) new_y = int(round(x*sin(angle) + y*cos(angle))) new_z = int(round(z)) val = self.matrix.item((z, y, x)) new_coords.append((val, new_x, new_y, new_z)) if new_x < min_new_x: min_new_x = new_x if new_x > max_new_x: max_new_x = new_x if new_y < min_new_y: min_new_y = new_y if new_y > max_new_y: max_new_y = new_y if new_z < min_new_z: min_new_z = new_z if new_z > max_new_z: max_new_z = new_z new_x_offset = abs(min_new_x) new_y_offset = abs(min_new_y) new_z_offset = abs(min_new_z) new_width = abs(min_new_x - max_new_x) new_height = abs(min_new_y - max_new_y) new_depth = abs(min_new_z - max_new_z) rotated = np.empty((new_depth + 1, new_height + 1, new_width + 1)) rotated.fill(0) for coord in new_coords: val = coord[0] x = coord[1] y = coord[2] z = coord[3] if rotated[new_z_offset + z][new_y_offset + y][new_x_offset + x] == 0: rotated[new_z_offset + z][new_y_offset + y][new_x_offset + x] = val self.matrix = rotated 

The way I use the above code is:

cube = Rect_Prism(20, 20, 20) # creates a 3d array similar to above example, just bigger cube.rotate(20, "x") cube.rotate(60, "y") 

Rect_Prism creates a MxNxD matrix, but in this case NxNxN.

And result when printing:

 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

You have to create a rotation matrix and multiply your this matrix for your array. Here the information

Wikipedea rotation matrix information

An example for 2d rotating

I think you should consider a "vector" representation for your data, instead of the current "raster" representation.

A vector representation would mean that, instead of each "voxel" being defined by its position in a grid, you would have some sort of list of voxels, with actual 3D coordinates.

So instead of having an MxNxD matrix where each voxel is a "black/white" dot, you could have a Mx3 matrix, where each row is a point, with columns being X, Y and Z.

That way, you would multiply the list by a 3x3 rotation matrix, and get another list of transformed coordinates.

You would remain with the problem of "rendering" these vector points (or lines, better yet) to a raster matrix (either pixels or voxels, but your sample image looks like 3D info has been projected to 2D space). There are a lot of techniques for doing this.