Function which allows you to perform a rigid transformation (translation and rotation) on 3-D volumetric intensity data.

rigid_transform_3d.py

from scipy.ndimage import map_coordinates
import numpy as np

def rigid_transform_3d(input_im, alpha,beta,gamma,xt,yt,zt,order=2):
    """
    Function which performs a rigid-body transformation on volumetric intensity data:
    
    -- example usuage --
    %pylab inline #in jupyter or ipython include this for the visualisation.
    inw = np.zeros((30,35,50))
    inw[10:20,10:20,10:20] = 1.0
    inw[11:19,11:19,11:19] = 0.0 
    outw = rigid_transform_3d(inw, 0.1,0.0,0.0,5,0,0)
    subplot(1,2,1) 
    imshow(inw[15,:,:])
    subplot(1,2,2)
    imshow(outw[15,:,:])
    
    -- inputs --
    
    input_im   input_volume, with dimensions [z, y, x]
    alpha      rotate around the z axis (Radians)
    beta       rotate around the y axis  (Radians)
    gamma      rotate around the x axis(Radians)
    xt         translate x axis (pixels)
    yt         translate y axis (pixels)
    zt         translate z axis (pixels)
    order      the order of the interpolation (3-D nearest neighbour = 0, trilinear 1, tricubic = 2 (default) )
    
    -- outputs -- 
    return     transformed volume, with dimensions [z, y, x]
    """
    
    z,y,x = np.meshgrid(np.arange(0,input_im.shape[0]),np.arange(0,input_im.shape[1]),np.arange(0,input_im.shape[2]))

    #Values for populationing transformation matrix.
    a00 = np.cos(alpha)*np.cos(beta)
    a01 = np.cos(alpha)*np.sin(beta)*np.sin(gamma)-np.sin(alpha)*np.cos(gamma)
    a02 = np.cos(alpha)*np.sin(beta)*np.cos(gamma)+np.sin(alpha)*np.sin(gamma)
    a03 = -xt
    a10 = np.sin(alpha)*np.cos(beta)
    a11 = np.sin(alpha)*np.sin(beta)*np.sin(gamma)+np.cos(alpha)*np.cos(gamma)
    a12 = np.sin(alpha)*np.sin(beta)*np.cos(gamma)-np.cos(alpha)*np.sin(gamma)
    a13 = -yt
    a20 = -np.sin(beta)
    a21 = np.cos(beta)*np.sin(gamma)
    a22 = np.cos(beta)*np.cos(gamma)
    a23 = -zt
    a30 = 0
    a31 = 0
    a32 = 0
    a33 = 1.

    matrix = np.mat([[a00,a01,a02,a03],[a10,a11,a12,a13],[a20,a21,a22,a23],[a30,a31,a32,a33]])
    
    #Do the transformation
    input_m = np.mat([np.array(x).reshape(-1),y.reshape(-1),z.reshape(-1),np.ones(x.shape).reshape(-1)])
    out = (matrix)*(input_m)
    
    #This is the function which samples the data at floating point coordinates (i.e. interpolates between neighbouring points).
    #Play close attention to the order:
    # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html

    outw = map_coordinates(input_im,[out[2],out[1],out[0]],order=order).reshape(input_im.shape[1],input_im.shape[0],input_im.shape[2])
    return np.rollaxis(outw,1)