Use affine matrix to convert 3D vol or 2D img to orthogonal one with corresponding affine matrix.
3D Affine matrix is such a 4x4 matrix:
M = [ [R T]; [0 0 0 1] ];
[x2 y2 z2 1]’ = M * [x1 y1 z1 1]’ ;
where, R is a 3x3 rotation matrix, and T is a 3x1 translation matrix. [x1 y1 z1] is a location in original 3D image volume, and [x2 y2 z2] is a location in transformed 3D image volume.
Although affine geometric transformation is only limited to parallel projection, it applies to most of the common geometric transformations, including rotation, translation, scaling, reflection, shearing, etc. Some of the common 3D Affine matrices are listed at the bottom of this description.
It may seem easy if you just want to apply the 3D affine matrix to each [x y z] coordinates in the 3D image volume. However, it turns out to be more complicated when you realize that the voxel (smallest element of 3D image volume) in the transformed image may no longer be an orthogonal cube. Therefore, a 3D interpolation algorithms must also be included to transform 3D image volume. I have implemented three interpolation methods in this program.
Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);
old_img  original 2D image or 3D volume. We assume x for the 1st dimension, y for the 2nd dimension, and z for the 3rd dimension.
old_M  a 3x3 2D affine matrix for 2D image, or a 4x4 3D affine matrix for 3D volume. We assume x for the 1st dimension, y for the 2nd dimension, and z for the 3rd dimension.
new_elem_size (optional)  size of voxel along x y z direction for a transformed 3D volume, or size of pixel along x y for a transformed 2D image. We assume x for the 1st dimension y for the 2nd dimension, and z for the 3rd dimension. 'new_elem_size' is 1 if it is default or empty.
verbose (optional)  1, 0
1: show transforming progress in percentage
2: progress will not be displayed
'verbose' is 1 if it is default or empty.
bg (optional)  background voxel intensity in any extra corner that is caused by 3D interpolation. 0 in most cases. 'bg' will be the average of two corner voxel intensities in original image volume, if it is default or empty.
method (optional)  1, 2, or 3
1: for Trilinear interpolation
2: for Nearest Neighbor interpolation
3: for Fischer's Bresenham interpolation
'method' is 1 if it is default or empty.
I suggest that you should use Method 1 (Trilinear) unless you have good reasons to choose other methods. Method 2 (Nearest Neighbor) is slightly faster, but will bring larger interpolation error. The Method 3 (Fischer's Bresenham) is only an implementation to test the special algorithm, and you also need to download my 3D Bresenham's line generation program from:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21057
new_img  transformed 3D image volume
new_M  transformed affine matrix
This program is inspired by:
SPM5 Software from Wellcome Trust Centre for Neuroimaging
http://www.fil.ion.ucl.ac.uk/spm/software
Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid Transformations to Volume Data, WSCG2004 Conference.
http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
Although this program (affine3d.m) can also run on MATLAB earlier than version 6.5 (R13), the speed will be hundred times slower. So I suggest that you run this program on MATLAB version 6.5 (R13) and above.
The following script is a demo to show you how to use this program. The result is displayed in the above image.
Example 1 (3D rotation):
load mri.mat; old_img = double(squeeze(D));
old_M = [0.88 0.5 3 90; 0.5 0.88 3 126; 0 0 2 72; 0 0 0 1];
new_img = affine(old_img, old_M, 2);
[x y z] = meshgrid(1:128,1:128,1:27);
sz = size(new_img);
[x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
figure; slice(x, y, z, old_img, 64, 64, 13.5);
shading flat; colormap(map); view(66, 66);
figure; slice(x1, y1, z1, new_img, sz(1)/2, sz(2)/2, sz(3)/2);
shading flat; colormap(map); view(66, 66);
Example 2 (2D interpolation):
load mri.mat; old_img=D(:,:,1,13)';
old_M = [1 0 0; 0 1 0; 0 0 1];
new_img = affine(old_img, old_M, [.2 .4]);
figure; image(old_img); colormap(map);
figure; image(new_img); colormap(map);
Appendix: Common 3D Affine matrices:
Translation (displacement of [dX dY dZ] for all voxels):
old_M = [
1 0 0 dX
0 1 0 dY
0 0 1 dZ
0 0 0 1 ];
Rotation (about X axis):
old_M = [
1 0 0 0
0 cosX sinX 0
0 sinX cosX 0
0 0 0 1 ];
Scaling:
old_M = [
sX 0 0 0
0 sY 0 0
0 0 sZ 0
0 0 0 1 ];
1.4  change title 

1.3  update didn't show up 

Added a "verbose" parameter to provide an opportunity to disable the progress display. 

Updated Description above and listed Common 3D Affine matrices, including: Translation, reflection, yaw, pitch, roll, and scaling. 

update program comments 
2one (view profile)
Great code!
I am using it to transform a 3d image volume from one perspective to another.
The original image volume (old_img) is of size 409 x 389 x 162 (where I know the pix>millimeter conversion for x, y and z).
My transformation matrix,
M =
0.0337 0.0958 0.9616 2.5360
0.3706 0.0107 0.1012 150.2009
0.0201 0.3571 0.2550 231.0513
0 0 0 1.0000
All other input were default.
The output image volume (new_img) qualitatively looks correct but the size is now: 206 x 172 x 189 and I am not sure how to work out the new pix>millimeter conversion. new_elem_size was blank so default (=1).
merel (view profile)
ted p teng (view profile)
thanks for sharing Jimmy!
Jimmy Shen (view profile)
The output volume is determined by the original volume, the affine matrix, and the new voxel_size that you choose. If any of the above inputs change, the output will also change. Otherwise, the output should be consistent.
The difference in dimension between output volume and original one is obvious. Just think about a unit square in 2D plane with dimension of [1 1]. If you rotate it 45 degree, you have to use a dimension of [1.414 1.414] to hold the same square. However, the distance from origin to edges are still 1.
In order to properly crop the transformed volume, you need to keep in mind that we are focusing on the unit of distance from the originator based on the affine matrix (e.g. millimeter), and don't worry about voxel size or its dimension. Since I have not used it for registration, I don't have detail procedure for you.
Giacomo (view profile)
Hey Jimmy,
thanks for your submission, quite helpful. However, I need to perform a crop on the new volume in order to compare it with the old one (I'm using this to register volumes, actually).
Trying to perform some series of tests translating the new volume in the z direction with [0:0.1:1] steps, I get an output volume sometimes with a bigger size than the original, and sometimes not, and I have troubles understanding what is happening/what to do to properly crop the transformed volume.
Can you help me out on this please?
Thanks a lot!
Jimmy Shen (view profile)
The new_img is interpolated from the old_img. Therefore, there is no exact mapping between old point and new point. I agree with your thoughts. However, keep in mind that we are talking about voxels (little cubic, with voxel_size) rather than points (with no size). Otherwise, you don't need this program, and a plot3 is enough to do affine transformation for points.
Myles (view profile)
Hi Jimmy,
I'm currently rotating a 3D block of data (a series of 2D slices), and my new image looks correct based on the matrix I supply. However, I am unable to plot a single older point to the 'correct' new image. I feel as though the new point should simply be an older point in XYZ multiplied by my transformation Matrix (as there is no rotation), however this does not provide me with the correct location in the new image. Do I need to use some additional information stored in new_M output?
Thanks
Samuel Hurley (view profile)
This is great. I replaced a hacked solution of my own that worked terribly with this in literally less that 5 minutes.
sun (view profile)
very good ,I like it
Jimmy Shen (view profile)
I think you mean that the dimension of the output matrix is different from the input matrix, which is correct.
Think about you have a cubic, and you rotate it 45 degree along Z axis. If we the element size the same, the dimension of X & Y should be at least sqrt(2) larger than the original X & Y dimension, since the new X & Y dimension is equal to the old diagonal.
In addition, there will be some extra elements in the new corner, which will be assigned a value of [bg], which is 0 by default.
Please let me know if you still have difficulty to understand this.
ucd puri (view profile)
Dear Jimmy, Many thanks for this submission. I am trying to transform a volume using a Rigid matrix as following
m=[ 0.9956 0.0491 0.0800 10.8112
0.0472 0.9986 0.0256 6.8496
0.0812 0.0217 0.9965 6.9082
0 0 0 1.0000];
he output seems screwed as expected, however, the size of the input matrix has changed drastically. That does not seems quite right. Could you kindly suggest.
Many Thanks
Jimmy Shen (view profile)
Hi Hamid:
The translation means the translation of the volume originator, which is specified by your affine matrix (old_M).
If there is no rotation, scaling, shearing, etc, the volume is untouched, and the new affine matrix (new_M) uses the translation that you specified.
In other words, with new_vol/new_M, you can plot a volume with the originator correctly translated. You can also do so with old_vol/old_M. So in your case, new_vol/new_M should be the same as old_vol/old_M, which is what we expected.
Keep in mind that this program is used to convert the oblique volume or image to orthogonal volume or image.
I think I need to change the title of this program to make it clearer.
Thank you for your feedback,
Jimmy
Hamid Sarve (view profile)
Dear Jimmy,
Thank you for the contribution.
Translation seems buggy, at least it doesn't work at all for me. I've been using a
[ 1 0 0 dx;
0 1 0 dy;
0 0 1 dz;
0 0 0 1;]
matrix.
The input and output volume are pretty much the same.
Regards,
Hamid
wafa BT (view profile)
hi ....
it is verey nice ..
i would like to have help plz ..
Ihave 3d array I insert slices to it and I need to compared 2 slices 4 pixel togather if it equal to the value than I will divided it by 8 and so on
in the demo code the the size of new_img is different from the old you get some errors. To correct you can use:
sz = size(new_img)
[x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
figure; slice(x1, y1, z1, new_img, 64, 64, 13);
shading flat; colormap(map); view(66, 66);
That is excellent, thanks for sharing.