mars_new_space

 make a new image space to contain image with rotations etc
 FORMAT [dim2, mat2, vox2] = mars_new_space(dim, mat, vox)
 
 dim        - original dimensions in voxels
 mat        - orignal mat file (4x4 transformation matrix)
 vox        - required ouput voxel size
 
 OUTPUT
 dim2       - new dimensions
 mat2       - new mat file
 vox2       - new voxel dimensions

 $Id$

0001 function [dim2, mat2, vox2] = mars_new_space(dim, mat, vox)
0002 % make a new image space to contain image with rotations etc
0003 % FORMAT [dim2, mat2, vox2] = mars_new_space(dim, mat, vox)
0004 %
0005 % dim        - original dimensions in voxels
0006 % mat        - orignal mat file (4x4 transformation matrix)
0007 % vox        - required ouput voxel size
0008 %
0009 % OUTPUT
0010 % dim2       - new dimensions
0011 % mat2       - new mat file
0012 % vox2       - new voxel dimensions
0013 %
0014 % $Id$
0015 
0016 if nargin < 2
0017   error('Need two input args');
0018 end
0019 dim = dim(:);
0020 dim = dim(1:3)';
0021 if nargin < 3
0022   vox = [];
0023 end
0024 
0025 % size, opposite corners of transformed img in mm
0026 [sz mn_mx] = mmsz(dim, mat);
0027 
0028 % select new voxel size if needed
0029 if isempty(vox) 
0030   % original voxel size
0031   vox = sqrt(sum(mat(1:3,1:3).^2));
0032 
0033   % XYZ max difference for one voxel
0034   vxsz = mmsz([1 1 1], mat);
0035 
0036   % reassign original dimensions to best matching
0037   % of new dimensions
0038   [t o] = sort(vxsz);
0039   vox2(o) = sort(vox);
0040 else
0041   vox2 = vox;
0042 end
0043 
0044 % new dimensions add 1 to allow for half voxel at either
0045 % side of voxel centres, allowing for rounding
0046 dim2 = sz./vox2 + 1;
0047 rdim2 = round(dim2);
0048 tiny = 1e-12;
0049 for d = 1:3
0050  if (dim2(d) - rdim2(d))<tiny
0051    dim2(d) = rdim2(d);
0052  else
0053    dim2 = ceil(dim2(d));
0054  end
0055 end
0056 
0057 % set new voxel sizes in output mat
0058 mat2 = diag([vox2 1]);
0059 
0060 % translations are from left post inf vox co-ord to mm coord
0061 % left post inf corner of new image
0062 mat2(1:3,4) = [mn_mx(1,:) - vox2]'; 
0063 
0064 return
0065 
0066 function [sz, mn_mx]  = mmsz(dim, mat);
0067 % returns size in mm of matrix dim;
0068   
0069 % 8 corners in voxels of original image
0070 i = [1 1 1; eye(3)];
0071 i = logical([i; ~i]);
0072 corners = ones(8,1) * dim;
0073 corners(i) = 1;
0074 
0075 % corners in mm
0076 mm_c = mat * [corners'; ones(1, 8)];
0077 
0078 % min and max of XYZ of corners in mm
0079 mm_c = sort(mm_c(1:3, :)');
0080 mn_mx = mm_c([1 end],:);
0081 
0082 % size is the difference
0083 sz = diff(mn_mx);
0084