NIFTI_20090325/affine.m
87b08fd7
 %  Using 2D or 3D affine matrix to rotate, translate, scale, reflect and
 %  shear a 2D image or 3D volume. 2D image is represented by a 2D matrix,
 %  3D volume is represented by a 3D matrix, and data type can be real 
 %  integer or floating-point.
 %
 %  You may notice that MATLAB has a function called 'imtransform.m' for
 %  2D spatial transformation. However, keep in mind that 'imtransform.m'
 %  assumes y for the 1st dimension, and x for the 2nd dimension. They are
 %  equivalent otherwise.
 %
 %  In addition, if you adjust the 'new_elem_size' parameter, this 'affine.m'
 %  is equivalent to 'interp2.m' for 2D image, and equivalent to 'interp3.m'
 %  for 3D volume.
 %
 %  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.
 %
 %		You can increase its value to decrease the resampling rate,
 %		and make the 2D image or 3D volume more coarse. It works
 %		just like 'interp3'.
 %
 %  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 the interpolation. 0 in most cases. If it is
 %		default or empty, 'bg' will be the average of two corner
 %		voxel intensities in original data.
 %
 %  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.
 %
 %  new_img  -	transformed 2D image or 3D volume
 %
 %  new_M  -	transformed affine matrix
 %
 %  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);
 %
 %  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
 %  
 %  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
 %
 function [new_img, new_M] = affine(old_img, old_M, new_elem_size, verbose, bg, method)
 
    if ~exist('old_img','var') | ~exist('old_M','var')
       error('Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);');
    end
 
    if ndims(old_img) == 3
       if ~isequal(size(old_M),[4 4])
          error('old_M should be a 4x4 affine matrix for 3D volume.');
       end
    elseif ndims(old_img) == 2
       if ~isequal(size(old_M),[3 3])
          error('old_M should be a 3x3 affine matrix for 2D image.');
       end
    else
       error('old_img should be either 2D image or 3D volume.');
    end
 
    if ~exist('new_elem_size','var') | isempty(new_elem_size)
       new_elem_size = [1 1 1];
    elseif length(new_elem_size) < 2
       new_elem_size = new_elem_size(1)*ones(1,3);
    elseif length(new_elem_size) < 3
       new_elem_size = [new_elem_size(:); 1]';
    end
 
    if ~exist('method','var') | isempty(method)
       method = 1;
    elseif ~exist('bresenham_line3d.m','file') & method == 3
       error([char(10) char(10) 'Please download 3D Bresenham''s line generation program from:' char(10) char(10) 'http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21057' char(10) char(10) 'to test Fischer''s Bresenham interpolation method.' char(10) char(10)]);
    end
 
    %  Make compatible to MATLAB earlier than version 7 (R14), which
    %  can only perform arithmetic on double data type
    %
    old_img = double(old_img);
    old_dim = size(old_img);
 
    if ~exist('bg','var') | isempty(bg)
       bg = mean([old_img(1) old_img(end)]);
    end
 
    if ~exist('verbose','var') | isempty(verbose)
       verbose = 1;
    end
 
    if ndims(old_img) == 2
       old_dim(3) = 1;
       old_M = old_M(:, [1 2 3 3]);
       old_M = old_M([1 2 3 3], :);
       old_M(3,:) = [0 0 1 0];
       old_M(:,3) = [0 0 1 0]';
    end
 
    %  Vertices of img in voxel
    %
    XYZvox = [	1		1		1
 		1		1		old_dim(3)
 		1		old_dim(2)	1
 		1		old_dim(2)	old_dim(3)
 		old_dim(1)	1		1
 		old_dim(1)	1		old_dim(3)
 		old_dim(1)	old_dim(2)	1
 		old_dim(1)	old_dim(2)	old_dim(3)   ]';
 
    old_R = old_M(1:3,1:3);
    old_T = old_M(1:3,4);
 
    %  Vertices of img in millimeter
    %
    XYZmm = old_R*(XYZvox-1) + repmat(old_T, [1, 8]);
 
    %  Make scale of new_M according to new_elem_size
    %
    new_M = diag([new_elem_size 1]);
 
    %  Make translation so minimum vertex is moved to [1,1,1]
    %
    new_M(1:3,4) = round( min(XYZmm,[],2) );
 
    %  New dimensions will be the maximum vertices in XYZ direction (dim_vox)
    %  i.e. compute   dim_vox   via   dim_mm = R*(dim_vox-1)+T
    %  where, dim_mm = round(max(XYZmm,[],2));
    %
    new_dim = ceil(new_M(1:3,1:3) \ ( round(max(XYZmm,[],2))-new_M(1:3,4) )+1)';
 
    %  Initialize new_img with new_dim
    %
    new_img = zeros(new_dim(1:3));
 
    %  Mask out any changes from Z axis of transformed volume, since we
    %  will traverse it voxel by voxel below. We will only apply unit
    %  increment of mask_Z(3,4) to simulate the cursor movement
    %
    %  i.e. we will use   mask_Z * new_XYZvox   to replace   new_XYZvox
    %
    mask_Z = diag(ones(1,4));
    mask_Z(3,3) = 0;
 
    %  It will be easier to do the interpolation if we invert the process
    %  by not traversing the original volume. Instead, we traverse the
    %  transformed volume, and backproject each voxel in the transformed 
    %  volume back into the original volume. If the backprojected voxel
    %  in original volume is within its boundary, the intensity of that
    %  voxel can be used by the cursor location in the transformed volume.
    %
    %  First, we traverse along Z axis of transformed volume voxel by voxel
    %
    for z = 1:new_dim(3)
 
       if verbose & ~mod(z,10)
          fprintf('%.2f percent is done.\n', 100*z/new_dim(3));
       end
 
       %  We need to find out the mapping from voxel in the transformed
       %  volume (new_XYZvox) to voxel in the original volume (old_XYZvox)
       %
       %  The following equation works, because they all equal to XYZmm:
       %  new_R*(new_XYZvox-1) + new_T  ==  old_R*(old_XYZvox-1) + old_T
       %
       %  We can use modified new_M1 & old_M1 to substitute new_M & old_M
       %      new_M1 * new_XYZvox       ==       old_M1 * old_XYZvox
       %
       %  where: M1 = M;   M1(:,4) = M(:,4) - sum(M(:,1:3),2);
       %  and:             M(:,4) == [T; 1] == sum(M1,2)
       %
       %  Therefore:   old_XYZvox = old_M1 \ new_M1 * new_XYZvox;
       %
       %  Since we are traverse Z axis, and   new_XYZvox   is replaced
       %  by   mask_Z * new_XYZvox, the above formula can be rewritten
       %  as:    old_XYZvox = old_M1 \ new_M1 * mask_Z * new_XYZvox;
       %
       %  i.e. we find the mapping from new_XYZvox to old_XYZvox:
       %  M = old_M1 \ new_M1 * mask_Z;
       %
       %  First, compute modified old_M1 & new_M1
       %
       old_M1 = old_M;   old_M1(:,4) = old_M(:,4) - sum(old_M(:,1:3),2);
       new_M1 = new_M;   new_M1(:,4) = new_M(:,4) - sum(new_M(:,1:3),2);
 
       %  Then, apply unit increment of mask_Z(3,4) to simulate the
       %  cursor movement
       %
       mask_Z(3,4) = z;
 
       %  Here is the mapping from new_XYZvox to old_XYZvox
       %
       M = old_M1 \ new_M1 * mask_Z;
 
       switch method
       case 1
          new_img(:,:,z) = trilinear(old_img, new_dim, old_dim, M, bg);
       case 2
          new_img(:,:,z) = nearest_neighbor(old_img, new_dim, old_dim, M, bg);
       case 3
          new_img(:,:,z) = bresenham(old_img, new_dim, old_dim, M, bg);
       end
 
    end;			% for z
 
    if ndims(old_img) == 2
       new_M(3,:) = [];
       new_M(:,3) = [];
    end
 
    return;					% affine
 
 
 %--------------------------------------------------------------------
 function img_slice = trilinear(img, dim1, dim2, M, bg)
 
    img_slice = zeros(dim1(1:2));
    TINY = 5e-2;					% tolerance
 
    %  Dimension of transformed 3D volume
    %
    xdim1 = dim1(1);
    ydim1 = dim1(2);
 
    %  Dimension of original 3D volume
    %
    xdim2 = dim2(1);
    ydim2 = dim2(2);
    zdim2 = dim2(3);
 
    %  initialize new_Y accumulation
    %
    Y2X = 0;
    Y2Y = 0;
    Y2Z = 0;
 
    for y = 1:ydim1
 
       %  increment of new_Y accumulation
       %
       Y2X = Y2X + M(1,2);		% new_Y to old_X
       Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
       Y2Z = Y2Z + M(3,2);		% new_Y to old_Z
 
       %  backproject new_Y accumulation and translation to old_XYZ
       %
       old_X = Y2X + M(1,4);
       old_Y = Y2Y + M(2,4);
       old_Z = Y2Z + M(3,4);
 
       for x = 1:xdim1
 
          %  accumulate the increment of new_X, and apply it
          %  to the backprojected old_XYZ
          %
          old_X = M(1,1) + old_X  ;
          old_Y = M(2,1) + old_Y  ;
          old_Z = M(3,1) + old_Z  ;
 
          %  within boundary of original image
          %
          if (	old_X > 1-TINY & old_X < xdim2+TINY & ...
 		old_Y > 1-TINY & old_Y < ydim2+TINY & ...
 		old_Z > 1-TINY & old_Z < zdim2+TINY	)
 
             %  Calculate distance of old_XYZ to its neighbors for
             %  weighted intensity average
             %
             dx = old_X - floor(old_X);
             dy = old_Y - floor(old_Y);
             dz = old_Z - floor(old_Z);
 
             x000 = floor(old_X);
             x100 = x000 + 1;
 
             if floor(old_X) < 1
                x000 = 1;
                x100 = x000;
             elseif floor(old_X) > xdim2-1
                x000 = xdim2;
                x100 = x000;
             end
 
             x010 = x000;
             x001 = x000;
             x011 = x000;
 
             x110 = x100;
             x101 = x100;
             x111 = x100;
 
             y000 = floor(old_Y);
             y010 = y000 + 1;
 
             if floor(old_Y) < 1
                y000 = 1;
                y100 = y000;
             elseif floor(old_Y) > ydim2-1
                y000 = ydim2;
                y010 = y000;
             end
 
             y100 = y000;
             y001 = y000;
             y101 = y000;
 
             y110 = y010;
             y011 = y010;
             y111 = y010;
 
             z000 = floor(old_Z);
             z001 = z000 + 1;
 
             if floor(old_Z) < 1
                z000 = 1;
                z001 = z000;
             elseif floor(old_Z) > zdim2-1
                z000 = zdim2;
                z001 = z000;
             end
 
             z100 = z000;
             z010 = z000;
             z110 = z000;
 
             z101 = z001;
             z011 = z001;
             z111 = z001;
 
             x010 = x000;
             x001 = x000;
             x011 = x000;
 
             x110 = x100;
             x101 = x100;
             x111 = x100;
 
             v000 = double(img(x000, y000, z000));
             v010 = double(img(x010, y010, z010));
             v001 = double(img(x001, y001, z001));
             v011 = double(img(x011, y011, z011));
 
             v100 = double(img(x100, y100, z100));
             v110 = double(img(x110, y110, z110));
             v101 = double(img(x101, y101, z101));
             v111 = double(img(x111, y111, z111));
 
             img_slice(x,y) = v000*(1-dx)*(1-dy)*(1-dz) + ...
                v010*(1-dx)*dy*(1-dz) + ...
                v001*(1-dx)*(1-dy)*dz + ...
                v011*(1-dx)*dy*dz + ...
                v100*dx*(1-dy)*(1-dz) + ...
                v110*dx*dy*(1-dz) + ...
                v101*dx*(1-dy)*dz + ...
                v111*dx*dy*dz;
 
          else
             img_slice(x,y) = bg;
 
          end	% if boundary
 
       end	% for x
    end		% for y
 
    return;					% trilinear
 
 
 %--------------------------------------------------------------------
 function img_slice = nearest_neighbor(img, dim1, dim2, M, bg)
 
    img_slice = zeros(dim1(1:2));
 
    %  Dimension of transformed 3D volume
    %
    xdim1 = dim1(1);
    ydim1 = dim1(2);
 
    %  Dimension of original 3D volume
    %
    xdim2 = dim2(1);
    ydim2 = dim2(2);
    zdim2 = dim2(3);
 
    %  initialize new_Y accumulation
    %
    Y2X = 0;
    Y2Y = 0;
    Y2Z = 0;
 
    for y = 1:ydim1
 
       %  increment of new_Y accumulation
       %
       Y2X = Y2X + M(1,2);		% new_Y to old_X
       Y2Y = Y2Y + M(2,2);		% new_Y to old_Y
       Y2Z = Y2Z + M(3,2);		% new_Y to old_Z
 
       %  backproject new_Y accumulation and translation to old_XYZ
       %
       old_X = Y2X + M(1,4);
       old_Y = Y2Y + M(2,4);
       old_Z = Y2Z + M(3,4);
 
       for x = 1:xdim1
 
          %  accumulate the increment of new_X and apply it
          %  to the backprojected old_XYZ
          %
          old_X = M(1,1) + old_X  ;
          old_Y = M(2,1) + old_Y  ;
          old_Z = M(3,1) + old_Z  ;
 
          xi = round(old_X);
          yi = round(old_Y);
          zi = round(old_Z);
 
          %  within boundary of original image
          %
          if (	xi >= 1 & xi <= xdim2 & ...
 		yi >= 1 & yi <= ydim2 & ...
 		zi >= 1 & zi <= zdim2	)
 
             img_slice(x,y) = img(xi,yi,zi);
 
          else
             img_slice(x,y) = bg;
 
          end	% if boundary
 
       end	% for x
    end		% for y
 
    return;					% nearest_neighbor
 
 
 %--------------------------------------------------------------------
 function img_slice = bresenham(img, dim1, dim2, M, bg)
 
    img_slice = zeros(dim1(1:2));
 
    %  Dimension of transformed 3D volume
    %
    xdim1 = dim1(1);
    ydim1 = dim1(2);
 
    %  Dimension of original 3D volume
    %
    xdim2 = dim2(1);
    ydim2 = dim2(2);
    zdim2 = dim2(3);
 
    for y = 1:ydim1
 
       start_old_XYZ = round(M*[0     y 0 1]');
       end_old_XYZ   = round(M*[xdim1 y 0 1]');
 
       [X Y Z] = bresenham_line3d(start_old_XYZ, end_old_XYZ);
 
       %  line error correction
       %
 %      del = end_old_XYZ - start_old_XYZ;
  %     del_dom = max(del);
   %    idx_dom = find(del==del_dom);
    %   idx_dom = idx_dom(1);
     %  idx_other = [1 2 3];
      % idx_other(idx_dom) = [];
       %del_x1 = del(idx_other(1));
 %      del_x2 = del(idx_other(2));
  %     line_slope = sqrt((del_x1/del_dom)^2 + (del_x2/del_dom)^2 + 1);
   %    line_error = line_slope - 1;
 % line error correction removed because it is too slow
 
       for x = 1:xdim1
 
          %  rescale ratio
          %
          i = round(x * length(X) / xdim1);
 
          if i < 1
             i = 1;
          elseif i > length(X)
             i = length(X);
          end
 
          xi = X(i);
          yi = Y(i);
          zi = Z(i);
 
          %  within boundary of the old XYZ space
          %
          if (	xi >= 1 & xi <= xdim2 & ...
 		yi >= 1 & yi <= ydim2 & ...
 		zi >= 1 & zi <= zdim2	)
 
             img_slice(x,y) = img(xi,yi,zi);
 
 %            if line_error > 1
  %              x = x + 1;
 
 %               if x <= xdim1
  %                 img_slice(x,y) = img(xi,yi,zi);
   %                line_error = line_slope - 1;
    %            end
     %        end		% if line_error
 % line error correction removed because it is too slow
 
          else
             img_slice(x,y) = bg;
 
          end	% if boundary
 
       end	% for x
    end		% for y
 
    return;					% bresenham