avw_reslice - Linear interpolation (using interp3) Usage: avw = avw_reslice(avw,dim,pixdim) avw is the Analyze struct returned by avw_read dim is the new volume dimensions (default, dim = [256,256,256]) pixdim is the new voxel resolution (default, pixdim = [1,1,1] %mm)
0001 function avw = avw_reslice(avw,dim,pixdim) 0002 0003 % avw_reslice - Linear interpolation (using interp3) 0004 % 0005 % Usage: avw = avw_reslice(avw,dim,pixdim) 0006 % 0007 % avw is the Analyze struct returned by avw_read 0008 % 0009 % dim is the new volume dimensions (default, dim = [256,256,256]) 0010 % pixdim is the new voxel resolution (default, pixdim = [1,1,1] %mm) 0011 % 0012 0013 % $Revision: 1.1 $ $Date: 2005/08/15 21:56:28 $ 0014 0015 % Copyright (C) 2005 Darren L. Weber 0016 % 0017 % This program is free software; you can redistribute it and/or 0018 % modify it under the terms of the GNU General Public License 0019 % as published by the Free Software Foundation; either version 2 0020 % of the License, or (at your option) any later version. 0021 % 0022 % This program is distributed in the hope that it will be useful, 0023 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0024 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0025 % GNU General Public License for more details. 0026 % 0027 % You should have received a copy of the GNU General Public License 0028 % along with this program; if not, write to the Free Software 0029 % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, 0030 % USA. 0031 0032 % History: 04/2005, Darren.Weber_at_radiology.ucsf.edu 0033 % 0034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0035 0036 0037 ver = '$Revision: 1.1 $ $Date: 2005/08/15 21:56:28 $'; 0038 fprintf('AVW_RESLICE [v %s]\n',ver(11:15)); 0039 0040 if ~exist('dim','var'), dim = [256,256,256]; end 0041 if isempty(dim), dim = [256,256,256]; end 0042 0043 if ~exist('pixdim','var'), pixdim = [1.0,1.0,1.0]; end 0044 if isempty(pixdim), pixdim = [1.0,1.0,1.0]; end 0045 0046 range = dim .* pixdim; 0047 0048 old_dim = double(avw.hdr.dime.dim(2:4)); 0049 old_pixdim = double(avw.hdr.dime.pixdim(2:4)); 0050 old_range = old_dim .* old_pixdim; 0051 0052 fprintf('...reslicing from: dim = [%d, %d, %d], pixdim = [%6.3f, %6.3f, %6.3f]\n',old_dim,old_pixdim); 0053 fprintf('...reslicing to: dim = [%d, %d, %d], pixdim = [%6.3f, %6.3f, %6.3f]\n',dim,pixdim); 0054 0055 tic; 0056 0057 % NOTE: 0058 % if we start at zero, we get dim(x)+1 values, one value for each voxel 'edge' 0059 % but this is incompatible with the interp3 function, so we start from the 0060 % first pixdim value below and we get dim(x) values. 0061 %x = 0:old_pixdim(1):old_dim(1)*old_pixdim(1); 0062 0063 0064 % contruct the old meshgrid 0065 x = old_pixdim(1):old_pixdim(1):old_range(1); 0066 y = old_pixdim(2):old_pixdim(2):old_range(2); 0067 z = old_pixdim(3):old_pixdim(3):old_range(3); 0068 [X,Y,Z] = meshgrid(x,y,z); 0069 0070 % contruct the new meshgrid 0071 xi = pixdim(1):pixdim(1):range(1); 0072 yi = pixdim(2):pixdim(2):range(2); 0073 zi = pixdim(3):pixdim(3):range(3); 0074 0075 0076 % NOTE: 0077 % we use EXTRAPVAL below because the range of the original volume may be 0078 % less than the new volume. When this is the case, the interp3 function 0079 % will not offset the old volume into the center of the new volume. For 0080 % this to happen, we have to shift the xi/yi/zi values by half the 0081 % difference in the volume size. That is: 0082 0083 if old_range(1) < range(1), 0084 xi_shift = ( old_range(1) - range(1) ) / 2; 0085 xi = xi + xi_shift; 0086 end 0087 if old_range(2) < range(2), 0088 yi_shift = ( old_range(2) - range(2) ) / 2; 0089 yi = yi + yi_shift; 0090 end 0091 if old_range(3) < range(3), 0092 zi_shift = ( old_range(3) - range(3) ) / 2; 0093 zi = zi + zi_shift; 0094 end 0095 [XI,YI,ZI] = meshgrid(xi,yi,zi); 0096 0097 % run the interpolation 0098 EXTRAPVAL = 0.0; 0099 VI = interp3(X,Y,Z,avw.img,XI,YI,ZI,'linear',EXTRAPVAL); 0100 0101 avw.hdr.dime.dim(2:4) = int16(dim); 0102 avw.hdr.dime.pixdim(2:4) = single(pixdim); 0103 avw.img = VI; 0104 0105 t = toc; fprintf('done (%5.2f sec)\n',t); 0106 0107 return