File Exchange

image thumbnail

Resampling volume or image with affine matrix

version 1.4 (5.38 KB) by

Use affine matrix to convert 3D vol or 2D img to orthogonal one with corresponding affine matrix.

16 Downloads

Updated

View License

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 ];

Comments and Ratings (16)

2one

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

merel (view profile)

ted p teng

ted p teng (view profile)

thanks for sharing Jimmy!

Jimmy Shen

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

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

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

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

This is great. I replaced a hacked solution of my own that worked terribly with this in literally less that 5 minutes.

sun

sun (view profile)

very good ,I like it

Jimmy Shen

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

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

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

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

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

gg kiss

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);

jichao zhao

That is excellent, thanks for sharing.

Updates

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

MATLAB Release
MATLAB 6.5 (R13)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video