vox2wt_ana

 Does image and mask conversion for wavelet analysis from voxel analysis
 FORMAT [phiwD] = vox2wt_ana(phiwD, wtinfo)
 
 Inputs 
 phiwD       - phido design object
 params          - structure containing options as fields; [defaults]
                   'wavelet'    - phiwave wavelet object to transform
                      images [phiw_lemarie(2)] 
                   'scales'     - scales for wavelet transform [4]
                   'wtprefix'   - prefix for wavelet transformed files
                      ['wv_']
                   'maskthresh' - threshold for mask image [0.05]

 Outputs
 phiwD       - modified phido design object
 
 The function first calculates the masking for the wavelet analysis, and
 stores results in the design masking structure.  It next checks if the
 images in the design are already WT'ed, by looking in their .mat file
 information.  If not, does WT and stores new vol structs in design

 Matthew Brett 9/10/00

 $Id: vox2wt_ana.m,v 1.4 2005/06/05 04:42:22 matthewbrett Exp $

0001 function [phiwD] = vox2wt_ana(phiwD, params)
0002 % Does image and mask conversion for wavelet analysis from voxel analysis
0003 % FORMAT [phiwD] = vox2wt_ana(phiwD, wtinfo)
0004 %
0005 % Inputs
0006 % phiwD       - phido design object
0007 % params          - structure containing options as fields; [defaults]
0008 %                   'wavelet'    - phiwave wavelet object to transform
0009 %                      images [phiw_lemarie(2)]
0010 %                   'scales'     - scales for wavelet transform [4]
0011 %                   'wtprefix'   - prefix for wavelet transformed files
0012 %                      ['wv_']
0013 %                   'maskthresh' - threshold for mask image [0.05]
0014 %
0015 % Outputs
0016 % phiwD       - modified phido design object
0017 %
0018 % The function first calculates the masking for the wavelet analysis, and
0019 % stores results in the design masking structure.  It next checks if the
0020 % images in the design are already WT'ed, by looking in their .mat file
0021 % information.  If not, does WT and stores new vol structs in design
0022 %
0023 % Matthew Brett 9/10/00
0024 %
0025 % $Id: vox2wt_ana.m,v 1.4 2005/06/05 04:42:22 matthewbrett Exp $
0026   
0027 % Default object structure
0028 defparams = struct('wavelet',  phiw_lemarie(2), ...
0029            'scales',   4, ...
0030            'wtprefix', 'wv_', ...
0031            'maskthresh', 0.05);
0032 
0033 if nargin < 2
0034   params = [];
0035 end
0036 params = mars_struct('ffillsplit', defparams, params);
0037 
0038 % Get images
0039 if ~has_images(phiwD), error('Need images in design'); end
0040 VY = get_images(phiwD);
0041 
0042 % Get wt'ed and non wt'ed images from vols
0043 g_wt_opts = struct('reproc', 1, ...
0044            'verbose', verbose(phiwD), ...
0045            'find_similar', 1);
0046 for i = 1:prod(size(VY))
0047   [v1 v2] = sf_get_wted(VY(i), params, g_wt_opts);
0048   [VY_wt(i) VY_o(i)] = deal(v1, v2);
0049 end
0050 VY_wt = reshape(VY_wt, size(VY));
0051 VY_o  = reshape(VY_o, size(VY));
0052 
0053 % Get masking structure
0054 xM = masking_struct(phiwD);
0055 if isempty(xM), error('Need masking structure'), end
0056 
0057 % Check if the mask matches the current wt information. If not, we will need
0058 % the original images to calculate the mask from.
0059 if isfield(xM, 'wave')
0060   redo_mask = ~same_wtinfo(xM.wave, params)
0061 else
0062   redo_mask = 1;
0063 end
0064 if redo_mask
0065   xM = sf_make_mask(xM, VY_o, params);
0066 end
0067 
0068 % return new design
0069 phiwD = set_images(phiwD, VY_wt);
0070 phiwD = masking_struct(phiwD, xM);
0071 
0072 return
0073 
0074 % Subfunctions
0075 % ------------
0076 
0077 function [Vw, Vo] = sf_get_wted(V, params, options)
0078 % Takes vol struct, returned WT'ed, non WT'ed vol structs
0079 % FORMAT [Vw, Vo] = sf_get_wted(V, params, options)
0080 %
0081 % Input
0082 % V        - WT'ed or non WT'ed vol struct
0083 % params   - wavelet transform infomation structure
0084 % options  - structure containing none or more of the following fields as
0085 %            options:
0086 %              'reproc'  - redo WT on already WT'ed images, if parameters
0087 %                 dont match wtinfo
0088 %              'verbose' - display messages
0089 %              'find_similar' - if processing image, look for appropriate
0090 %                 previously transformed image
0091 %
0092 % Images can be:
0093 % Voxel images (not WTed) -> return to-be-WTed object
0094 % WT'ed images -> return WT'ed object
0095 % Voxel images with matching WT'ed image -> return WT'ed object
0096 
0097 f_s = mars_struct('getifthere', options, 'find_similar');  
0098 wv_opts = struct('noproc', 1, 'wtprefix', params.wtprefix, ...
0099          'find_similar', f_s);
0100 
0101 % Flag whether to reprocess if already WT'ed but with wrong parameters
0102 reproc_f = mars_struct('isthere', options, 'reproc');
0103 
0104 if phiw_wvimg('is_wted', V)  % already transformed
0105 
0106   % Are the images WT'ed as we would like?
0107   Vo = phiw_wvimg('orig_vol', V);
0108   wvobj = phiw_wvimg(V, wv_opts);
0109   if same_wtinfo(wvobj, params) % yes
0110     Vw = V;
0111     return
0112   end
0113   
0114   % no, maybe reprocess
0115   if ~reproc_f
0116     error(sprintf('%s is already WTed with different parameters', ...
0117           V.fname))
0118   end
0119   options.find_similar = 0;
0120   [Vw, Vo] = sf_get_wted(Vo, params, options);
0121   return
0122 
0123 end  
0124 
0125 % Make wvimg object, maybe looking for compatible already processed
0126 % image
0127 wvobj = phiw_wvimg(V, wv_opts, params.wavelet, params.scales);
0128 if ~is_wted(wvobj)
0129   wvobj = write_wtimg(wvobj);
0130 end
0131 Vw = wv_vol(wvobj);
0132 Vo = V;
0133 
0134 return
0135 
0136 function xM = sf_make_mask(xM, VY, params)
0137 % Process and recreate masking structure
0138 % FORMAT xM = sf_make_mask(xM, VY, params)
0139 
0140 % masking, first
0141 %-If xM is not a structure then assumme it's a vector of thresholds
0142 %------------------------------------------------------------------
0143 if ~isstruct(xM)
0144   xM = struct(    'T',    [],...
0145         'TH',    xM,...
0146         'I',    0,...
0147         'VM',    {[]},...
0148         'xs',    struct('Masking','analysis threshold'));
0149 end
0150 
0151 % mask in voxel space
0152 mask = zeros(VY(1).dim(1:3));
0153 
0154 % Create masking image
0155 if isempty(xM.TH)
0156   xM.TH = zeros(size(VY))-Inf;
0157   thmf = 0;
0158 else
0159   thmf = any(xM.TH(:) ~= -Inf);
0160 end
0161 [vol X Y Z] = deal(VY(1),1,2,3);
0162 % plane coordinate stuff
0163 plsz = prod(vol.dim(X:Y));
0164 xords  = [1:vol.dim(X)]'*ones(1,vol.dim(Y)); %-plane X coordinates
0165 yords  = ones(vol.dim(X),1)*[1:vol.dim(Y)];  % plane Y coordinates
0166 tPl = [xords(:) yords(:)]';
0167 % implicit 0 masking flag vector
0168 YNaNrep = spm_type(VY(1,1).dim(4),'nanrep');
0169 IM = xM.I & ~YNaNrep & xM.TH(:,1)<0;
0170     
0171 fprintf(' Calculating masking image...');
0172 % cycle over planes
0173 for z = 1:vol.dim(3)
0174   pPl = tPl; % xyz of current in mask points
0175   pPl(Z,:) = ones(1, plsz)*z;
0176   % intensity thresholding
0177   if thmf
0178     for i = 1:prod(size(VY))
0179       if isempty(pPl), break,end
0180       vals = spm_sample_vol(VY(i),pPl(X,:),pPl(Y,:),pPl(Z,:),0);
0181       msk = vals>xM.TH(i);
0182       if IM(i), msk = msk & abs(vals)>eps; end
0183       pPl = pPl(:, msk);
0184     end
0185   end
0186   % explicit masking images
0187   for i = 1:prod(size(xM.VM))
0188     if isempty(pPl), break,end
0189     tM   = inv(xM.VM(i).mat)*vol.mat;        %-Reorientation matrix
0190     tmp  = tM * [pPl;ones(1,size(pPl,2))];    %-Coords in mask image
0191     
0192     %-Load mask image within current mask & update mask
0193     %--------------------------------------------------
0194     vals = spm_sample_vol(xM.VM(i),tmp(X,:),tmp(Y,:),tmp(Z,:),0);
0195     pPl = pPl(:, vals>0);
0196   end
0197   if ~isempty(pPl)
0198     inds = pPl(X,:)+(pPl(Y,:)-1)*vol.dim(X)+(z-1)*plsz;
0199     mask(inds) = 1;
0200   end
0201 end
0202 fprintf('...done\n');
0203 
0204 % write mask image in voxel space
0205 VM = struct('fname',    'voxmask.img',...
0206         'dim',      [vol.dim(1:3),spm_type('uint8')],...
0207         'mat',      vol.mat,...
0208         'pinfo',    [1 0 0]',...
0209         'descrip',    'phiw:resultant analysis mask');
0210 fprintf(' Saving the intial analysis mask as %s ...', VM.fname);
0211 spm_write_vol(VM, mask);
0212 fprintf('...done\n');
0213 
0214 % Transform, do processing for wavelet space mask
0215 wvVM = pr_wvmask(VM, params);
0216 
0217 % save new mask
0218 fprintf(' Saving the wavelet transformed mask as %s...', wvVM.wvol.fname);
0219 wvm = write_wtimg(wvVM);
0220 fprintf('...done\n');
0221   
0222 % make new masking structure
0223 xM = struct(...
0224     'TH',ones(size(VY))*-Inf,...
0225     'T', -Inf,...
0226     'I', 0,...
0227     'VM',wvm.wvol, ...
0228     'xs', struct(...
0229     'Analysis_threshold','None (-Inf)',...
0230     'Explicit_masking', 'Yes'),...
0231      'wave', thin(wvm));
0232 
0233 return