📄 affine3d.m
字号:
% Using 3D affine matrix to rotate, translate, scale, reflect and
% shear a 3D image volume. 3D image volume is represented by a 3D
% matrix, and data type can be real integer or floating-point.
%
% Usage: [new_img new_M] = ...
% affine3d(old_img, old_M, [method], [voxel_size], [bg], [verbose]);
%
% old_img - original 3D image volume
%
% old_M - 3D affine matrix (4x4 matrix)
%
% 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.
%
% voxel_size (optional) - size of voxel along x y z direction for
% transformed 3D image volume. 'voxel_size' is
% 1 if it is default or empty. You can increase
% its value to decrease the resampling rate, and
% make the 3D image more coarse.
%
% 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.
%
% 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.
%
% new_img - transformed 3D image volume
%
% new_M - transformed affine matrix
%
% Example:
% 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 = affine3d(old_img, old_M, 1, 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);
%
% 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] = ...
affine3d(old_img, old_M, method, new_voxel_size, bg, verbose)
if ~exist('old_img','var') | ~exist('old_M','var')
error('Usage: [new_img new_M] = affine3d(old_img, old_M, [method], [voxel_size], [bg], [verbose])');
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
if ~exist('new_voxel_size','var') | isempty(new_voxel_size)
new_voxel_size = [1 1 1];
elseif length(new_voxel_size) < 3
new_voxel_size = new_voxel_size(1)*ones(1,3);
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
% 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 voxel_size
%
new_M = diag([new_voxel_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 image, 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 image volume. Instead, we traverse
% the transformed image, and backproject each voxel in transformed
% image volume back into the original image. If the backprojected
% voxel in original image is within its boundary, the intensity of
% that voxel can be used by the cursor location in the transformed
% image.
%
% First, we traverse along Z axis of transformed image 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
% image (new_XYZvox) to voxel in the original image (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
return; % affine3d
%--------------------------------------------------------------------
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -