 Christoph Budziszewski authored on30/03/2009 17:54:25
 1 1 `new file mode 100644` ... ... `@@ -0,0 +1,554 @@` 1 ```+% Using 2D or 3D affine matrix to rotate, translate, scale, reflect and ``` 2 ```+% shear a 2D image or 3D volume. 2D image is represented by a 2D matrix, ``` 3 ```+% 3D volume is represented by a 3D matrix, and data type can be real ``` 4 ```+% integer or floating-point. ``` 5 ```+% ``` 6 ```+% You may notice that MATLAB has a function called 'imtransform.m' for ``` 7 ```+% 2D spatial transformation. However, keep in mind that 'imtransform.m' ``` 8 ```+% assumes y for the 1st dimension, and x for the 2nd dimension. They are ``` 9 ```+% equivalent otherwise. ``` 10 ```+% ``` 11 ```+% In addition, if you adjust the 'new_elem_size' parameter, this 'affine.m' ``` 12 ```+% is equivalent to 'interp2.m' for 2D image, and equivalent to 'interp3.m' ``` 13 ```+% for 3D volume. ``` 14 ```+% ``` 15 ```+% Usage: [new_img new_M] = ... ``` 16 ```+% affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]); ``` 17 ```+% ``` 18 ```+% old_img - original 2D image or 3D volume. We assume x for the 1st ``` 19 ```+% dimension, y for the 2nd dimension, and z for the 3rd ``` 20 ```+% dimension. ``` 21 ```+% ``` 22 ```+% old_M - a 3x3 2D affine matrix for 2D image, or a 4x4 3D affine ``` 23 ```+% matrix for 3D volume. We assume x for the 1st dimension, ``` 24 ```+% y for the 2nd dimension, and z for the 3rd dimension. ``` 25 ```+% ``` 26 ```+% new_elem_size (optional) - size of voxel along x y z direction for ``` 27 ```+% a transformed 3D volume, or size of pixel along x y for ``` 28 ```+% a transformed 2D image. We assume x for the 1st dimension ``` 29 ```+% y for the 2nd dimension, and z for the 3rd dimension. ``` 30 ```+% 'new_elem_size' is 1 if it is default or empty. ``` 31 ```+% ``` 32 ```+% You can increase its value to decrease the resampling rate, ``` 33 ```+% and make the 2D image or 3D volume more coarse. It works ``` 34 ```+% just like 'interp3'. ``` 35 ```+% ``` 36 ```+% verbose (optional) - 1, 0 ``` 37 ```+% 1: show transforming progress in percentage ``` 38 ```+% 2: progress will not be displayed ``` 39 ```+% 'verbose' is 1 if it is default or empty. ``` 40 ```+% ``` 41 ```+% bg (optional) - background voxel intensity in any extra corner that ``` 42 ```+% is caused by the interpolation. 0 in most cases. If it is ``` 43 ```+% default or empty, 'bg' will be the average of two corner ``` 44 ```+% voxel intensities in original data. ``` 45 ```+% ``` 46 ```+% method (optional) - 1, 2, or 3 ``` 47 ```+% 1: for Trilinear interpolation ``` 48 ```+% 2: for Nearest Neighbor interpolation ``` 49 ```+% 3: for Fischer's Bresenham interpolation ``` 50 ```+% 'method' is 1 if it is default or empty. ``` 51 ```+% ``` 52 ```+% new_img - transformed 2D image or 3D volume ``` 53 ```+% ``` 54 ```+% new_M - transformed affine matrix ``` 55 ```+% ``` 56 ```+% Example 1 (3D rotation): ``` 57 ```+% load mri.mat; old_img = double(squeeze(D)); ``` 58 ```+% old_M = [0.88 0.5 3 -90; -0.5 0.88 3 -126; 0 0 2 -72; 0 0 0 1]; ``` 59 ```+% new_img = affine(old_img, old_M, 2); ``` 60 ```+% [x y z] = meshgrid(1:128,1:128,1:27); ``` 61 ```+% sz = size(new_img); ``` 62 ```+% [x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3)); ``` 63 ```+% figure; slice(x, y, z, old_img, 64, 64, 13.5); ``` 64 ```+% shading flat; colormap(map); view(-66, 66); ``` 65 ```+% figure; slice(x1, y1, z1, new_img, sz(1)/2, sz(2)/2, sz(3)/2); ``` 66 ```+% shading flat; colormap(map); view(-66, 66); ``` 67 ```+% ``` 68 ```+% Example 2 (2D interpolation): ``` 69 ```+% load mri.mat; old_img=D(:,:,1,13)'; ``` 70 ```+% old_M = [1 0 0; 0 1 0; 0 0 1]; ``` 71 ```+% new_img = affine(old_img, old_M, [.2 .4]); ``` 72 ```+% figure; image(old_img); colormap(map); ``` 73 ```+% figure; image(new_img); colormap(map); ``` 74 ```+% ``` 75 ```+% This program is inspired by: ``` 76 ```+% SPM5 Software from Wellcome Trust Centre for Neuroimaging ``` 77 ```+% http://www.fil.ion.ucl.ac.uk/spm/software ``` 78 ```+% Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid ``` 79 ```+% Transformations to Volume Data, WSCG2004 Conference. ``` 80 ```+% http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf ``` 81 ```+% ``` 82 ```+% - Jimmy Shen (jimmy@rotman-baycrest.on.ca) ``` 83 ```+% ``` 84 ```+function [new_img, new_M] = affine(old_img, old_M, new_elem_size, verbose, bg, method) ``` 85 ```+ ``` 86 ```+ if ~exist('old_img','var') | ~exist('old_M','var') ``` 87 ```+ error('Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);'); ``` 88 ```+ end ``` 89 ```+ ``` 90 ```+ if ndims(old_img) == 3 ``` 91 ```+ if ~isequal(size(old_M),[4 4]) ``` 92 ```+ error('old_M should be a 4x4 affine matrix for 3D volume.'); ``` 93 ```+ end ``` 94 ```+ elseif ndims(old_img) == 2 ``` 95 ```+ if ~isequal(size(old_M),[3 3]) ``` 96 ```+ error('old_M should be a 3x3 affine matrix for 2D image.'); ``` 97 ```+ end ``` 98 ```+ else ``` 99 ```+ error('old_img should be either 2D image or 3D volume.'); ``` 100 ```+ end ``` 101 ```+ ``` 102 ```+ if ~exist('new_elem_size','var') | isempty(new_elem_size) ``` 103 ```+ new_elem_size = [1 1 1]; ``` 104 ```+ elseif length(new_elem_size) < 2 ``` 105 ```+ new_elem_size = new_elem_size(1)*ones(1,3); ``` 106 ```+ elseif length(new_elem_size) < 3 ``` 107 ```+ new_elem_size = [new_elem_size(:); 1]'; ``` 108 ```+ end ``` 109 ```+ ``` 110 ```+ if ~exist('method','var') | isempty(method) ``` 111 ```+ method = 1; ``` 112 ```+ elseif ~exist('bresenham_line3d.m','file') & method == 3 ``` 113 ```+ 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)]); ``` 114 ```+ end ``` 115 ```+ ``` 116 ```+ % Make compatible to MATLAB earlier than version 7 (R14), which ``` 117 ```+ % can only perform arithmetic on double data type ``` 118 ```+ % ``` 119 ```+ old_img = double(old_img); ``` 120 ```+ old_dim = size(old_img); ``` 121 ```+ ``` 122 ```+ if ~exist('bg','var') | isempty(bg) ``` 123 ```+ bg = mean([old_img(1) old_img(end)]); ``` 124 ```+ end ``` 125 ```+ ``` 126 ```+ if ~exist('verbose','var') | isempty(verbose) ``` 127 ```+ verbose = 1; ``` 128 ```+ end ``` 129 ```+ ``` 130 ```+ if ndims(old_img) == 2 ``` 131 ```+ old_dim(3) = 1; ``` 132 ```+ old_M = old_M(:, [1 2 3 3]); ``` 133 ```+ old_M = old_M([1 2 3 3], :); ``` 134 ```+ old_M(3,:) = [0 0 1 0]; ``` 135 ```+ old_M(:,3) = [0 0 1 0]'; ``` 136 ```+ end ``` 137 ```+ ``` 138 ```+ % Vertices of img in voxel ``` 139 ```+ % ``` 140 ```+ XYZvox = [ 1 1 1 ``` 141 ```+ 1 1 old_dim(3) ``` 142 ```+ 1 old_dim(2) 1 ``` 143 ```+ 1 old_dim(2) old_dim(3) ``` 144 ```+ old_dim(1) 1 1 ``` 145 ```+ old_dim(1) 1 old_dim(3) ``` 146 ```+ old_dim(1) old_dim(2) 1 ``` 147 ```+ old_dim(1) old_dim(2) old_dim(3) ]'; ``` 148 ```+ ``` 149 ```+ old_R = old_M(1:3,1:3); ``` 150 ```+ old_T = old_M(1:3,4); ``` 151 ```+ ``` 152 ```+ % Vertices of img in millimeter ``` 153 ```+ % ``` 154 ```+ XYZmm = old_R*(XYZvox-1) + repmat(old_T, [1, 8]); ``` 155 ```+ ``` 156 ```+ % Make scale of new_M according to new_elem_size ``` 157 ```+ % ``` 158 ```+ new_M = diag([new_elem_size 1]); ``` 159 ```+ ``` 160 ```+ % Make translation so minimum vertex is moved to [1,1,1] ``` 161 ```+ % ``` 162 ```+ new_M(1:3,4) = round( min(XYZmm,[],2) ); ``` 163 ```+ ``` 164 ```+ % New dimensions will be the maximum vertices in XYZ direction (dim_vox) ``` 165 ```+ % i.e. compute dim_vox via dim_mm = R*(dim_vox-1)+T ``` 166 ```+ % where, dim_mm = round(max(XYZmm,[],2)); ``` 167 ```+ % ``` 168 ```+ new_dim = ceil(new_M(1:3,1:3) \ ( round(max(XYZmm,[],2))-new_M(1:3,4) )+1)'; ``` 169 ```+ ``` 170 ```+ % Initialize new_img with new_dim ``` 171 ```+ % ``` 172 ```+ new_img = zeros(new_dim(1:3)); ``` 173 ```+ ``` 174 ```+ % Mask out any changes from Z axis of transformed volume, since we ``` 175 ```+ % will traverse it voxel by voxel below. We will only apply unit ``` 176 ```+ % increment of mask_Z(3,4) to simulate the cursor movement ``` 177 ```+ % ``` 178 ```+ % i.e. we will use mask_Z * new_XYZvox to replace new_XYZvox ``` 179 ```+ % ``` 180 ```+ mask_Z = diag(ones(1,4)); ``` 181 ```+ mask_Z(3,3) = 0; ``` 182 ```+ ``` 183 ```+ % It will be easier to do the interpolation if we invert the process ``` 184 ```+ % by not traversing the original volume. Instead, we traverse the ``` 185 ```+ % transformed volume, and backproject each voxel in the transformed ``` 186 ```+ % volume back into the original volume. If the backprojected voxel ``` 187 ```+ % in original volume is within its boundary, the intensity of that ``` 188 ```+ % voxel can be used by the cursor location in the transformed volume. ``` 189 ```+ % ``` 190 ```+ % First, we traverse along Z axis of transformed volume voxel by voxel ``` 191 ```+ % ``` 192 ```+ for z = 1:new_dim(3) ``` 193 ```+ ``` 194 ```+ if verbose & ~mod(z,10) ``` 195 ```+ fprintf('%.2f percent is done.\n', 100*z/new_dim(3)); ``` 196 ```+ end ``` 197 ```+ ``` 198 ```+ % We need to find out the mapping from voxel in the transformed ``` 199 ```+ % volume (new_XYZvox) to voxel in the original volume (old_XYZvox) ``` 200 ```+ % ``` 201 ```+ % The following equation works, because they all equal to XYZmm: ``` 202 ```+ % new_R*(new_XYZvox-1) + new_T == old_R*(old_XYZvox-1) + old_T ``` 203 ```+ % ``` 204 ```+ % We can use modified new_M1 & old_M1 to substitute new_M & old_M ``` 205 ```+ % new_M1 * new_XYZvox == old_M1 * old_XYZvox ``` 206 ```+ % ``` 207 ```+ % where: M1 = M; M1(:,4) = M(:,4) - sum(M(:,1:3),2); ``` 208 ```+ % and: M(:,4) == [T; 1] == sum(M1,2) ``` 209 ```+ % ``` 210 ```+ % Therefore: old_XYZvox = old_M1 \ new_M1 * new_XYZvox; ``` 211 ```+ % ``` 212 ```+ % Since we are traverse Z axis, and new_XYZvox is replaced ``` 213 ```+ % by mask_Z * new_XYZvox, the above formula can be rewritten ``` 214 ```+ % as: old_XYZvox = old_M1 \ new_M1 * mask_Z * new_XYZvox; ``` 215 ```+ % ``` 216 ```+ % i.e. we find the mapping from new_XYZvox to old_XYZvox: ``` 217 ```+ % M = old_M1 \ new_M1 * mask_Z; ``` 218 ```+ % ``` 219 ```+ % First, compute modified old_M1 & new_M1 ``` 220 ```+ % ``` 221 ```+ old_M1 = old_M; old_M1(:,4) = old_M(:,4) - sum(old_M(:,1:3),2); ``` 222 ```+ new_M1 = new_M; new_M1(:,4) = new_M(:,4) - sum(new_M(:,1:3),2); ``` 223 ```+ ``` 224 ```+ % Then, apply unit increment of mask_Z(3,4) to simulate the ``` 225 ```+ % cursor movement ``` 226 ```+ % ``` 227 ```+ mask_Z(3,4) = z; ``` 228 ```+ ``` 229 ```+ % Here is the mapping from new_XYZvox to old_XYZvox ``` 230 ```+ % ``` 231 ```+ M = old_M1 \ new_M1 * mask_Z; ``` 232 ```+ ``` 233 ```+ switch method ``` 234 ```+ case 1 ``` 235 ```+ new_img(:,:,z) = trilinear(old_img, new_dim, old_dim, M, bg); ``` 236 ```+ case 2 ``` 237 ```+ new_img(:,:,z) = nearest_neighbor(old_img, new_dim, old_dim, M, bg); ``` 238 ```+ case 3 ``` 239 ```+ new_img(:,:,z) = bresenham(old_img, new_dim, old_dim, M, bg); ``` 240 ```+ end ``` 241 ```+ ``` 242 ```+ end; % for z ``` 243 ```+ ``` 244 ```+ if ndims(old_img) == 2 ``` 245 ```+ new_M(3,:) = []; ``` 246 ```+ new_M(:,3) = []; ``` 247 ```+ end ``` 248 ```+ ``` 249 ```+ return; % affine ``` 250 ```+ ``` 251 ```+ ``` 252 ```+%-------------------------------------------------------------------- ``` 253 ```+function img_slice = trilinear(img, dim1, dim2, M, bg) ``` 254 ```+ ``` 255 ```+ img_slice = zeros(dim1(1:2)); ``` 256 ```+ TINY = 5e-2; % tolerance ``` 257 ```+ ``` 258 ```+ % Dimension of transformed 3D volume ``` 259 ```+ % ``` 260 ```+ xdim1 = dim1(1); ``` 261 ```+ ydim1 = dim1(2); ``` 262 ```+ ``` 263 ```+ % Dimension of original 3D volume ``` 264 ```+ % ``` 265 ```+ xdim2 = dim2(1); ``` 266 ```+ ydim2 = dim2(2); ``` 267 ```+ zdim2 = dim2(3); ``` 268 ```+ ``` 269 ```+ % initialize new_Y accumulation ``` 270 ```+ % ``` 271 ```+ Y2X = 0; ``` 272 ```+ Y2Y = 0; ``` 273 ```+ Y2Z = 0; ``` 274 ```+ ``` 275 ```+ for y = 1:ydim1 ``` 276 ```+ ``` 277 ```+ % increment of new_Y accumulation ``` 278 ```+ % ``` 279 ```+ Y2X = Y2X + M(1,2); % new_Y to old_X ``` 280 ```+ Y2Y = Y2Y + M(2,2); % new_Y to old_Y ``` 281 ```+ Y2Z = Y2Z + M(3,2); % new_Y to old_Z ``` 282 ```+ ``` 283 ```+ % backproject new_Y accumulation and translation to old_XYZ ``` 284 ```+ % ``` 285 ```+ old_X = Y2X + M(1,4); ``` 286 ```+ old_Y = Y2Y + M(2,4); ``` 287 ```+ old_Z = Y2Z + M(3,4); ``` 288 ```+ ``` 289 ```+ for x = 1:xdim1 ``` 290 ```+ ``` 291 ```+ % accumulate the increment of new_X, and apply it ``` 292 ```+ % to the backprojected old_XYZ ``` 293 ```+ % ``` 294 ```+ old_X = M(1,1) + old_X ; ``` 295 ```+ old_Y = M(2,1) + old_Y ; ``` 296 ```+ old_Z = M(3,1) + old_Z ; ``` 297 ```+ ``` 298 ```+ % within boundary of original image ``` 299 ```+ % ``` 300 ```+ if ( old_X > 1-TINY & old_X < xdim2+TINY & ... ``` 301 ```+ old_Y > 1-TINY & old_Y < ydim2+TINY & ... ``` 302 ```+ old_Z > 1-TINY & old_Z < zdim2+TINY ) ``` 303 ```+ ``` 304 ```+ % Calculate distance of old_XYZ to its neighbors for ``` 305 ```+ % weighted intensity average ``` 306 ```+ % ``` 307 ```+ dx = old_X - floor(old_X); ``` 308 ```+ dy = old_Y - floor(old_Y); ``` 309 ```+ dz = old_Z - floor(old_Z); ``` 310 ```+ ``` 311 ```+ x000 = floor(old_X); ``` 312 ```+ x100 = x000 + 1; ``` 313 ```+ ``` 314 ```+ if floor(old_X) < 1 ``` 315 ```+ x000 = 1; ``` 316 ```+ x100 = x000; ``` 317 ```+ elseif floor(old_X) > xdim2-1 ``` 318 ```+ x000 = xdim2; ``` 319 ```+ x100 = x000; ``` 320 ```+ end ``` 321 ```+ ``` 322 ```+ x010 = x000; ``` 323 ```+ x001 = x000; ``` 324 ```+ x011 = x000; ``` 325 ```+ ``` 326 ```+ x110 = x100; ``` 327 ```+ x101 = x100; ``` 328 ```+ x111 = x100; ``` 329 ```+ ``` 330 ```+ y000 = floor(old_Y); ``` 331 ```+ y010 = y000 + 1; ``` 332 ```+ ``` 333 ```+ if floor(old_Y) < 1 ``` 334 ```+ y000 = 1; ``` 335 ```+ y100 = y000; ``` 336 ```+ elseif floor(old_Y) > ydim2-1 ``` 337 ```+ y000 = ydim2; ``` 338 ```+ y010 = y000; ``` 339 ```+ end ``` 340 ```+ ``` 341 ```+ y100 = y000; ``` 342 ```+ y001 = y000; ``` 343 ```+ y101 = y000; ``` 344 ```+ ``` 345 ```+ y110 = y010; ``` 346 ```+ y011 = y010; ``` 347 ```+ y111 = y010; ``` 348 ```+ ``` 349 ```+ z000 = floor(old_Z); ``` 350 ```+ z001 = z000 + 1; ``` 351 ```+ ``` 352 ```+ if floor(old_Z) < 1 ``` 353 ```+ z000 = 1; ``` 354 ```+ z001 = z000; ``` 355 ```+ elseif floor(old_Z) > zdim2-1 ``` 356 ```+ z000 = zdim2; ``` 357 ```+ z001 = z000; ``` 358 ```+ end ``` 359 ```+ ``` 360 ```+ z100 = z000; ``` 361 ```+ z010 = z000; ``` 362 ```+ z110 = z000; ``` 363 ```+ ``` 364 ```+ z101 = z001; ``` 365 ```+ z011 = z001; ``` 366 ```+ z111 = z001; ``` 367 ```+ ``` 368 ```+ x010 = x000; ``` 369 ```+ x001 = x000; ``` 370 ```+ x011 = x000; ``` 371 ```+ ``` 372 ```+ x110 = x100; ``` 373 ```+ x101 = x100; ``` 374 ```+ x111 = x100; ``` 375 ```+ ``` 376 ```+ v000 = double(img(x000, y000, z000)); ``` 377 ```+ v010 = double(img(x010, y010, z010)); ``` 378 ```+ v001 = double(img(x001, y001, z001)); ``` 379 ```+ v011 = double(img(x011, y011, z011)); ``` 380 ```+ ``` 381 ```+ v100 = double(img(x100, y100, z100)); ``` 382 ```+ v110 = double(img(x110, y110, z110)); ``` 383 ```+ v101 = double(img(x101, y101, z101)); ``` 384 ```+ v111 = double(img(x111, y111, z111)); ``` 385 ```+ ``` 386 ```+ img_slice(x,y) = v000*(1-dx)*(1-dy)*(1-dz) + ... ``` 387 ```+ v010*(1-dx)*dy*(1-dz) + ... ``` 388 ```+ v001*(1-dx)*(1-dy)*dz + ... ``` 389 ```+ v011*(1-dx)*dy*dz + ... ``` 390 ```+ v100*dx*(1-dy)*(1-dz) + ... ``` 391 ```+ v110*dx*dy*(1-dz) + ... ``` 392 ```+ v101*dx*(1-dy)*dz + ... ``` 393 ```+ v111*dx*dy*dz; ``` 394 ```+ ``` 395 ```+ else ``` 396 ```+ img_slice(x,y) = bg; ``` 397 ```+ ``` 398 ```+ end % if boundary ``` 399 ```+ ``` 400 ```+ end % for x ``` 401 ```+ end % for y ``` 402 ```+ ``` 403 ```+ return; % trilinear ``` 404 ```+ ``` 405 ```+ ``` 406 ```+%-------------------------------------------------------------------- ``` 407 ```+function img_slice = nearest_neighbor(img, dim1, dim2, M, bg) ``` 408 ```+ ``` 409 ```+ img_slice = zeros(dim1(1:2)); ``` 410 ```+ ``` 411 ```+ % Dimension of transformed 3D volume ``` 412 ```+ % ``` 413 ```+ xdim1 = dim1(1); ``` 414 ```+ ydim1 = dim1(2); ``` 415 ```+ ``` 416 ```+ % Dimension of original 3D volume ``` 417 ```+ % ``` 418 ```+ xdim2 = dim2(1); ``` 419 ```+ ydim2 = dim2(2); ``` 420 ```+ zdim2 = dim2(3); ``` 421 ```+ ``` 422 ```+ % initialize new_Y accumulation ``` 423 ```+ % ``` 424 ```+ Y2X = 0; ``` 425 ```+ Y2Y = 0; ``` 426 ```+ Y2Z = 0; ``` 427 ```+ ``` 428 ```+ for y = 1:ydim1 ``` 429 ```+ ``` 430 ```+ % increment of new_Y accumulation ``` 431 ```+ % ``` 432 ```+ Y2X = Y2X + M(1,2); % new_Y to old_X ``` 433 ```+ Y2Y = Y2Y + M(2,2); % new_Y to old_Y ``` 434 ```+ Y2Z = Y2Z + M(3,2); % new_Y to old_Z ``` 435 ```+ ``` 436 ```+ % backproject new_Y accumulation and translation to old_XYZ ``` 437 ```+ % ``` 438 ```+ old_X = Y2X + M(1,4); ``` 439 ```+ old_Y = Y2Y + M(2,4); ``` 440 ```+ old_Z = Y2Z + M(3,4); ``` 441 ```+ ``` 442 ```+ for x = 1:xdim1 ``` 443 ```+ ``` 444 ```+ % accumulate the increment of new_X and apply it ``` 445 ```+ % to the backprojected old_XYZ ``` 446 ```+ % ``` 447 ```+ old_X = M(1,1) + old_X ; ``` 448 ```+ old_Y = M(2,1) + old_Y ; ``` 449 ```+ old_Z = M(3,1) + old_Z ; ``` 450 ```+ ``` 451 ```+ xi = round(old_X); ``` 452 ```+ yi = round(old_Y); ``` 453 ```+ zi = round(old_Z); ``` 454 ```+ ``` 455 ```+ % within boundary of original image ``` 456 ```+ % ``` 457 ```+ if ( xi >= 1 & xi <= xdim2 & ... ``` 458 ```+ yi >= 1 & yi <= ydim2 & ... ``` 459 ```+ zi >= 1 & zi <= zdim2 ) ``` 460 ```+ ``` 461 ```+ img_slice(x,y) = img(xi,yi,zi); ``` 462 ```+ ``` 463 ```+ else ``` 464 ```+ img_slice(x,y) = bg; ``` 465 ```+ ``` 466 ```+ end % if boundary ``` 467 ```+ ``` 468 ```+ end % for x ``` 469 ```+ end % for y ``` 470 ```+ ``` 471 ```+ return; % nearest_neighbor ``` 472 ```+ ``` 473 ```+ ``` 474 ```+%-------------------------------------------------------------------- ``` 475 ```+function img_slice = bresenham(img, dim1, dim2, M, bg) ``` 476 ```+ ``` 477 ```+ img_slice = zeros(dim1(1:2)); ``` 478 ```+ ``` 479 ```+ % Dimension of transformed 3D volume ``` 480 ```+ % ``` 481 ```+ xdim1 = dim1(1); ``` 482 ```+ ydim1 = dim1(2); ``` 483 ```+ ``` 484 ```+ % Dimension of original 3D volume ``` 485 ```+ % ``` 486 ```+ xdim2 = dim2(1); ``` 487 ```+ ydim2 = dim2(2); ``` 488 ```+ zdim2 = dim2(3); ``` 489 ```+ ``` 490 ```+ for y = 1:ydim1 ``` 491 ```+ ``` 492 ```+ start_old_XYZ = round(M*[0 y 0 1]'); ``` 493 ```+ end_old_XYZ = round(M*[xdim1 y 0 1]'); ``` 494 ```+ ``` 495 ```+ [X Y Z] = bresenham_line3d(start_old_XYZ, end_old_XYZ); ``` 496 ```+ ``` 497 ```+ % line error correction ``` 498 ```+ % ``` 499 ```+% del = end_old_XYZ - start_old_XYZ; ``` 500 ```+ % del_dom = max(del); ``` 501 ```+ % idx_dom = find(del==del_dom); ``` 502 ```+ % idx_dom = idx_dom(1); ``` 503 ```+ % idx_other = [1 2 3]; ``` 504 ```+ % idx_other(idx_dom) = []; ``` 505 ```+ %del_x1 = del(idx_other(1)); ``` 506 ```+% del_x2 = del(idx_other(2)); ``` 507 ```+ % line_slope = sqrt((del_x1/del_dom)^2 + (del_x2/del_dom)^2 + 1); ``` 508 ```+ % line_error = line_slope - 1; ``` 509 ```+% line error correction removed because it is too slow ``` 510 ```+ ``` 511 ```+ for x = 1:xdim1 ``` 512 ```+ ``` 513 ```+ % rescale ratio ``` 514 ```+ % ``` 515 ```+ i = round(x * length(X) / xdim1); ``` 516 ```+ ``` 517 ```+ if i < 1 ``` 518 ```+ i = 1; ``` 519 ```+ elseif i > length(X) ``` 520 ```+ i = length(X); ``` 521 ```+ end ``` 522 ```+ ``` 523 ```+ xi = X(i); ``` 524 ```+ yi = Y(i); ``` 525 ```+ zi = Z(i); ``` 526 ```+ ``` 527 ```+ % within boundary of the old XYZ space ``` 528 ```+ % ``` 529 ```+ if ( xi >= 1 & xi <= xdim2 & ... ``` 530 ```+ yi >= 1 & yi <= ydim2 & ... ``` 531 ```+ zi >= 1 & zi <= zdim2 ) ``` 532 ```+ ``` 533 ```+ img_slice(x,y) = img(xi,yi,zi); ``` 534 ```+ ``` 535 ```+% if line_error > 1 ``` 536 ```+ % x = x + 1; ``` 537 ```+ ``` 538 ```+% if x <= xdim1 ``` 539 ```+ % img_slice(x,y) = img(xi,yi,zi); ``` 540 ```+ % line_error = line_slope - 1; ``` 541 ```+ % end ``` 542 ```+ % end % if line_error ``` 543 ```+% line error correction removed because it is too slow ``` 544 ```+ ``` 545 ```+ else ``` 546 ```+ img_slice(x,y) = bg; ``` 547 ```+ ``` 548 ```+ end % if boundary ``` 549 ```+ ``` 550 ```+ end % for x ``` 551 ```+ end % for y ``` 552 ```+ ``` 553 ```+ return; % bresenham ``` 554 ```+ ```