How to transform a 3D image by an affine transformation matrix

Question

I have a question concerning Image Processing: I have a stack of images, which I can compose to a 3D image using Image3D. Additionally I have a 4x4 affine transformation matrix. I would like to transform the 3D image using my transformation matrix.

In 2D this would be possible using ImageTransformation[image, transformationMatrix]

t = AffineTransform[{{{RandomReal[{-0.5, 5}], 
 RandomReal[{-0.5, 5}]}, {RandomReal[{-0.5, 5}], 
 RandomReal[{-0.5, 5}]}}, {RandomReal[{0.2, 0.2}], 
RandomReal[{-0.2, 0.2}]}}]
ImageTransformation[ExampleData[{"TestImage", "Lena"}], t, PlotRange -> All, Padding -> None]

So this works fine for 2D images, but unfortunately ImageTransformation is not applicable to Image3D-objects.

My first attempt was to decompose the matrix into several matrices and transform the image stack in 2 steps:

1. transforming each image of the image stack
2. transpose the image stack and transform with the transformation matrix for the third dimension.

As an example I do the following steps for translation:

image3D = Import["ExampleData/CTengine.tiff", "Image3D"]
slices = Image3DSlices[image3D];
translation3D = TranslationTransform[{0.2, 0.1, 0.3}]
translation3D1 = TranslationTransform[{0.2, 0.1}]
translation3D2 = TranslationTransform[{0, 0.3}]
tempimageslices = ImageTransformation[#, translation3D1] & /@ slices;
transposedtempimageeslices =  Image /@ (Transpose[ImageData[#] & /@ tempimageslices]);
tempimageslices2 =   ImageTransformation[#, translation3D2] & /@ transposedtempimageeslices;
final = Image /@ (Transpose[ImageData[#] & /@ tempimageslices2, {2, 1, 3}]);
Image3D[final]

The Steps in this example are:

1. I cut the image into several slices.

2. I split the 4 x 4 transformation matrix "translation3D" into 2 3 x 3 matrices "translation3D1" and "translation3D2".

3. I transform every image in slices using ImageTransformation and the matrix "translation3D1"

4. I transpose the imageData and create new images in the new dimensions

5. I transform the images in the new dimensions using the matrix "translation3D2"

6. I transpose the images in order to have the original ImageDimensions

For a translation or a scaling, the two transformation matrices i need, are easy to read from the matrix. The problem I have with this method is that I can't compute the decomposition of an affine transformation matrix.

Is there an easy way to decompose the affine transformation or another solution in which the full affine transformation can be used to transform the image stack?

Answer 1

Edit: ImageTransformation and ImagePerspectiveTransformation accept Image3D objects and corresponding transformations as of Mathematica 10.

Waiting for ImageTransformation to accept Image3D objects, one can revert to the classic method:

For each voxel (i,j,k) in the result 3D image,

1. Compute the corresponding coordinates (x,y,z) in the input image using the transformation from the output image to the input image.
2. Interpolate the image values around the (x,y,z) position.

End

A top-level implementation is going to be very slow, but I would try compiling it to improve on the speed.

Another possibility for visualization purposes is to use Raster3D. This is adapted from the reference page of Raster3D, section Properties and Relations:

raster3D = Raster3D[
    RawArray["Byte", ImageData[image3D, "Byte"]], 
    {{-1, -1, -1}, {1, 1, 1}}, {0, 255}, ColorFunction -> "GrayLevelOpacity"];

Graphics3D[GeometricTransformation[raster3D, AffineTransform[RandomReal[1, {3, 3}]]]]

Answer 2

According to this Documentation example you can use GeometricTransformation with Image3D. For instance, this 3D image of an engine:

img = Import["ExampleData/CTengine.tiff", "Image3D"]

can be transformed by first constructing Raster3D:

engine = Raster3D[RawArray["Byte", Import["ExampleData/CTengine.tiff", "Data"]], 
{{-1, -1, -1}, {1, 1, 1}}, {0, 255}, ColorFunction -> "GrayLevelOpacity"]

and then using say a ShearingTransform:

Graphics3D[GeometricTransformation[engine, ShearingTransform[Pi/7, {1, 0, 0}, {0, 0, 1}]]]