get_wdimg

 calculates and returns denoised image in voxel space
 FORMAT [Vdcon, Vderr, pD, changef] = get_wdimg(pD, Ic, wdstruct, fname) 
   
 Inputs (defaults in [square brackets])
 pD         - phido design object
 Ic         - vector of contrast numbers defining statistic image [GUI]
              if length(Ic)>1 = conjunction. In this case the
              contrasts should be in orthogonalization order
 wdstruct   - structure with info for wavelet denoising.  
              Most fields are passed directly to the thresh_calc
              routine, see that routine for detailed comments.
              Fields are:
              'levels' - levels to calulcate / apply  thresholding [](=all)
              'thlev'  - level at which to apply threshold ['level']             
              'thcalc' - wavelet denoising type ['stein']
              'thapp'  - form of wavelet denoising ['linear']
              'ncalc'  - null hypothesis calculation ['n']
              'alpha'  - alpha for t etc Bonferroni etc correction [0.05]
              't2z'    - if not 0, does T to Z transform on T data [0]
              'write_err' - if not 0, writes std and sort-of t image [1]
 fname      - filename for denoised image [via GUI]

 Returns
 Vdcon      - spm vol struct for denoised image
 Vderr      - spm vol struct for error image
 pD         - possibly modified phido object (contrasts entereed)
 changef    - whether the design object was modified or not

 Based on spm_getSPM from the spm99 distribution
 (spm_getSPM, v2.35 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 00/01/2)

 Matthew Brett, Federico Turkheimer, 9/10/00

 $Id: get_wdimg.m,v 1.15 2005/07/01 22:24:07 matthewbrett Exp $

0001 function [Vdcon, Vderr, pD, changef] = get_wdimg(pD, Ic, wdstruct, fname) 
0002 % calculates and returns denoised image in voxel space
0003 % FORMAT [Vdcon, Vderr, pD, changef] = get_wdimg(pD, Ic, wdstruct, fname)
0004 %
0005 % Inputs (defaults in [square brackets])
0006 % pD         - phido design object
0007 % Ic         - vector of contrast numbers defining statistic image [GUI]
0008 %              if length(Ic)>1 = conjunction. In this case the
0009 %              contrasts should be in orthogonalization order
0010 % wdstruct   - structure with info for wavelet denoising.
0011 %              Most fields are passed directly to the thresh_calc
0012 %              routine, see that routine for detailed comments.
0013 %              Fields are:
0014 %              'levels' - levels to calulcate / apply  thresholding [](=all)
0015 %              'thlev'  - level at which to apply threshold ['level']
0016 %              'thcalc' - wavelet denoising type ['stein']
0017 %              'thapp'  - form of wavelet denoising ['linear']
0018 %              'ncalc'  - null hypothesis calculation ['n']
0019 %              'alpha'  - alpha for t etc Bonferroni etc correction [0.05]
0020 %              't2z'    - if not 0, does T to Z transform on T data [0]
0021 %              'write_err' - if not 0, writes std and sort-of t image [1]
0022 % fname      - filename for denoised image [via GUI]
0023 %
0024 % Returns
0025 % Vdcon      - spm vol struct for denoised image
0026 % Vderr      - spm vol struct for error image
0027 % pD         - possibly modified phido object (contrasts entereed)
0028 % changef    - whether the design object was modified or not
0029 %
0030 % Based on spm_getSPM from the spm99 distribution
0031 % (spm_getSPM, v2.35 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 00/01/2)
0032 %
0033 % Matthew Brett, Federico Turkheimer, 9/10/00
0034 %
0035 % $Id: get_wdimg.m,v 1.15 2005/07/01 22:24:07 matthewbrett Exp $
0036   
0037 if nargin < 2
0038   Ic = [];
0039 end
0040 if length(Ic) > 1
0041   error('Can only get one denoised contrast at a time');
0042 end
0043 if nargin < 3
0044   wdstruct = [];
0045 end
0046 if nargin < 4
0047   fname = '';
0048 end
0049 
0050 % default return values
0051 [Vdcon, Vderr] = deal([]);
0052 changef = 0;
0053 
0054 % Backspace macro
0055 bs30 = repmat(sprintf('\b'),1,30);
0056 
0057 % default denoising - see comments above
0058 def_struct = struct(...
0059     'levels', [], ...
0060     'thlev','level',...
0061     'thcalc', 'stein', ...
0062     'thapp', 'linear', ...
0063     'ncalc', 'n', ...
0064     'varpoolf',0,...
0065     'alpha', 0.05,...
0066     't2z', 1, ...
0067     'write_err', 1);
0068 
0069 wdstruct = mars_struct('ffillsplit', def_struct, wdstruct);
0070 
0071 % check if can write to current directory
0072 if ~swd_writable(pD)
0073   error(['Sorry, cannot write to directory: ' swd(pD)]);
0074 end
0075 
0076 %-Get Stats data from SPM.mat
0077 %-----------------------------------------------------------------------
0078 xX     = design_structure(pD);
0079 
0080 %-Get/Compute mm<->voxel matrices & image dimensions from SPM.mat
0081 %-----------------------------------------------------------------------
0082 VY    = get_images(pD);
0083 M     = VY(1).mat;
0084 iM    = inv(M);
0085 DIM   = VY(1).dim(1:3);
0086 
0087 %-Contrasts
0088 %=======================================================================
0089 if isempty(Ic)
0090   % only single t contrast allowed for the moment
0091   [Ic, pD, changef] = ui_get_contrasts(pD, ...
0092                        'T', ...
0093                        1, ...
0094                        'Select contrast...', ...
0095                        '', ...
0096                        1);
0097 end
0098 
0099 % Check we have only T contrasts
0100 xCon = get_contrasts(pD);
0101 if ~all([xCon(Ic).STAT] == 'T')
0102   error('Sorry, we can only denoise t contrasts');
0103 end
0104 
0105 % get filename if necessary
0106 %=======================================================================
0107 if isempty(fname)
0108   str  = sprintf('con%04d_%s',Ic,xCon(Ic).name);
0109   fname = spm_input('Filename for contrast', 1, 's', ...
0110             ['phiw_' mars_utils('str2fname',str)]);
0111   if isempty(fname), return, end
0112 end
0113 Qdir = spm_str_manip(fname,'Hv');
0114 fname = [spm_str_manip(fname,'stv'),'.img'];
0115 if ~strcmp(Qdir, '.')
0116   fprintf('Note: writing image to current directory');
0117 end
0118 
0119 % write con and statistic images
0120 %=======================================================================
0121 [pD, Ic, changef] = write_contrasts(pD, Ic);
0122 
0123 % get contrasts
0124 if changef, xCon = get_contrasts(pD); end
0125 xC1  = xCon(Ic);
0126 erdf = error_df(pD);
0127 edf   = [xC1.eidf erdf];
0128 
0129 % get contrast image
0130 wave = get_wave(pD);
0131 wvcon = phiw_wvimg(full_vol(pD, xC1.Vcon),struct('noproc',0), wave);
0132 
0133 % get and process error image
0134 VResMS = get_vol_field(pD, 'VResMS');
0135 wverr = phiw_wvimg(VResMS,[], wave);
0136 rmsi = as_matrix(wverr);
0137 rmsi(abs(rmsi)<eps) = NaN;
0138 con_cov = xC1.c'*xX.Bcov*xC1.c;
0139 rmsi = sqrt(rmsi .* con_cov);
0140 wverr = as_matrix(wverr, rmsi);
0141 
0142 % do wavelet denoise/inversion
0143 %=======================================================================
0144 % calculate denoising (which also removes NaNs)
0145 statinf = struct('stat','T','df',edf(2));
0146 if wdstruct.t2z, statinf.stat = 'TZ'; end
0147 fprintf('\t%-32s: %30s','Wavelet image','...calculate denoising')         %-#
0148 [th_obj, dndescrip] = thresh_calc(wvcon, wverr, statinf, wdstruct);
0149 fprintf('%s%30s\n',bs30,'...done')               %-#
0150 
0151 % check there is still some signal
0152 if all_null(th_obj)
0153   warning('No wavelet coefficients survive thresholding')
0154   return
0155 end
0156 
0157 % Apply thresholding to contrast image
0158 wvcond = thresh_apply(wvcon, th_obj, dndescrip);
0159 
0160 % inverse wavelet transform and save denoised image
0161 d_swd = swd(pD);
0162 Vdcon = struct(...
0163     'fname',  fullfile(d_swd,fname),...
0164     'dim',    [1 1 1,16],...
0165     'mat',    wave.ovol.mat,...
0166     'pinfo',  [1,0,0]',...
0167     'descrip',sprintf('Phiwave{%c} - %s',...
0168               xC1.STAT,xC1.name));
0169 Vdcon = write_iwtimg(wvcond,Vdcon);
0170 fprintf('%s%30s\n',bs30,'...done')               %-#
0171 
0172 % Write description into text file
0173 write_descrip(wvcond,Vdcon);
0174 
0175 % Create error map, if we have saved residuals
0176 VResI = get_vol_field(pD, 'VResI');
0177 if ~isempty(VResI) & wdstruct.write_err
0178   % Quite a lot of work to do here
0179   nScan = prod(size(VResI));
0180   oi = oimgi(wave);
0181   odim = diff(oi)+1;
0182   std_img = zeros(odim);
0183   sum_img = std_img;
0184     
0185   % get whole thresholding object as image, find mask
0186   th_img = as_matrix(th_obj);
0187   tmp = isfinite(th_img) & ~(th_img == 0);
0188   in_mask = find(tmp);
0189   out_mask = find(~tmp);
0190   th_img = th_img(in_mask);
0191   clear th_obj tmp;
0192 
0193   for i = 1:nScan
0194     obj = phiw_wvimg(VResI(i), [], wave);
0195     img = as_matrix(obj);
0196     img(out_mask) = 0;
0197     img(in_mask) = img(in_mask) .* th_img;
0198     img = invert(img, wave.wavelet, wave.scales, oi);
0199     sum_img = sum_img + img;
0200     std_img = std_img + img .^2;
0201   end
0202   xX = design_structure(pD);
0203   std_img = (std_img - (sum_img .^2 / nScan)) / xX.trRV;
0204   std_img = sqrt(std_img .* con_cov);
0205   
0206   % save std image
0207   efname = fullfile(d_swd, ['std_' fname]);
0208   Vderr = struct(...
0209       'fname',  efname,...
0210       'dim',    [odim,16],...
0211       'mat',    wave.ovol.mat,...
0212       'pinfo',  [1,0,0]',...
0213       'descrip',sprintf('Phiwave{%c} - error:%s',...
0214             xC1.STAT,xC1.name));
0215   Vderr = spm_write_vol(Vderr, std_img);
0216 
0217   % save sort-of t image
0218   con_img = spm_read_vols(Vdcon);
0219   t_msk = abs(con_img) > eps & abs(std_img) > eps;
0220   t_img = zeros(size(con_img));
0221   t_img(t_msk) = con_img(t_msk) ./ std_img(t_msk);
0222   tfname = fullfile(d_swd, ['t_' fname]);
0223   Vdt = struct(...
0224       'fname',  tfname,...
0225       'dim',    [odim,16],...
0226       'mat',    wave.ovol.mat,...
0227       'pinfo',  [1,0,0]',...
0228       'descrip',sprintf('Phiwave{%c} - t:%s',...
0229             xC1.STAT,xC1.name));
0230   Vdt = spm_write_vol(Vdt, t_img);
0231 
0232   fprintf('%s%30s\n',bs30,'...done');
0233 end
0234 return