Code covered by the BSD License  

Highlights from
Resampling volume or image with affine matrix

4.57143

4.6 | 7 ratings Rate this file 35 Downloads (last 30 days) File Size: 5.38 KB File ID: #21080
image thumbnail

Resampling volume or image with affine matrix

by

 

13 Aug 2008 (Updated )

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

| Watch this File

File Information
Description

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

MATLAB release MATLAB 6.5 (R13)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (14)
28 Nov 2013 ted p teng

thanks for sharing Jimmy!

11 Nov 2013 Jimmy Shen

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.

10 Nov 2013 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!

08 Nov 2013 Jimmy Shen

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.

08 Nov 2013 Myles

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

10 Jun 2013 Samuel Hurley

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

06 May 2013 sun

very good ,I like it

16 Feb 2012 Jimmy Shen

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.

16 Feb 2012 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

16 Apr 2009 Jimmy Shen

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

16 Apr 2009 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

25 Oct 2008 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

20 Aug 2008 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);

19 Aug 2008 jichao zhao

That is excellent, thanks for sharing.

Updates
15 Aug 2008

update program comments

19 Aug 2008

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

22 Aug 2008

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

24 Oct 2008

update didn't show up

16 Apr 2009

change title

Contact us