Description of write

0001 function [phiwD, connos, changef] = write_contrasts(phiwD, connos, flags) 
0002 % Writes contrast and maybe statistic images
0003 % FORMAT [phiwD, connos, changef] = write_contrasts(phiwD, connos, flags)
0004 %
0005 % Inputs
0006 % phiwD      - phido design object
0007 % connos     - vector of contrasts to write (fetched by GUI if empty)
0008 % flags      - flags containing none or more of fields
0009 %              'no_new' - if not 0 forbids new contrasts being entered
0010 %                 (when connos=[])
0011 %              't_only' - not 0 -> allow t contrasts only (where connos=[])
0012 %              'f_only' - not 0 -> allow f contrasts only (where connos=[])
0013 %              'single' - not 0, allows a Single contrast only (where
0014 %                 connos=[])
0015 %              'con_only' - don't write statistic images (i.e. con images only)
0016 %
0017 % Returns
0018 % phiwD      - modified design object
0019 % connos     - vector of contrast numbers (as input, or from GUI)
0020 % changef    - set to 1 is phiwD object has changed
0021 %
0022 % Based (very) closely on spm_getSPM from the spm99 distribution
0023 % (spm_getSPM, v2.35 Andrew Holmes, Karl Friston & Jean-Baptiste Poline
0024 % 00/01/2)
0025 %
0026 % Matthew Brett 9/10/00
0027 %
0028 % $Id: write_contrasts.m,v 1.10 2005/07/01 20:03:45 matthewbrett Exp $
0029 
0030 def_flags = struct(...
0031     'no_new', 0,...
0032     't_only', 0,...
0033     'f_only', 0,...
0034     'single', 0,...
0035     'con_only', 0);
0036   
0037 if ~is_phiw_estimated(phiwD), error('Need phiwave estimated design'); end
0038 changef = 0;
0039 if nargin < 2
0040   connos = [];
0041 end
0042 if nargin < 3
0043   fnames = '';
0044 end
0045 if nargin < 4
0046   flags = [];
0047 end
0048 rmsi = [];
0049 
0050 % flags decoding
0051 flags = mars_struct('ffillsplit', def_flags, flags);
0052 wOK = flags.no_new;   % forbid new contrasts
0053 if flags.t_only  % only F, only t flags for contrast selection
0054   statstr = 'T';
0055 elseif flags.f_only
0056   statstr = 'F';
0057 else
0058   statstr = 'T|F';
0059 end
0060 if flags.single  % only one contrast at a time
0061   ncons = 1;
0062 else
0063   ncons = Inf;
0064 end
0065 
0066 % Backspace macro
0067 bs30 = repmat(sprintf('\b'),1,30);
0068 
0069 % get needed stuff from design
0070 Vbeta  = get_vol_field(phiwD, 'Vbeta');
0071 VResMS = get_vol_field(phiwD, 'VResMS');
0072 VY   = get_images(phiwD);
0073 xM   = masking_struct(phiwD);
0074 erdf = error_df(phiwD);
0075 M    = VY(1).mat;
0076 DIM  = VY(1).dim(1:3);
0077 xX   = design_structure(phiwD);
0078 Swd  = swd(phiwD);
0079 
0080 %-See if can write to current directory
0081 %-----------------------------------------------------------------------
0082 wOK = swd_writable(phiwD);
0083 if ~wOK
0084   str = {'Can''t write to the results directory:',...
0085      ['        ',Swd],...
0086      ' ','-> results restricted to contrasts already computed'};
0087   spm('alert!',str,mfilename,1);
0088 end
0089 
0090 %-Get contrasts if needed (if multivariate there is only one structure)
0091 %-----------------------------------------------------------------------
0092 nVar    = size(VY,2);
0093 if nVar > 1 
0094   connos == 1;
0095 elseif isempty(connos) 
0096   [connos, phiwD, changef] = ui_get_contrasts(phiwD, statstr, ncons,...
0097              '    Select contrasts...','to compute', wOK);
0098 end
0099 
0100 if isempty(connos)
0101   return
0102 end
0103 xCon = get_contrasts(phiwD);
0104 
0105 % get wave transform parameters
0106 wave = get_wave(phiwD);
0107 if isempty(wave)
0108   error(['No wave info for first image - unlikely to be wavelet' ...
0109      ' analysis'])
0110 end
0111 wvp = getfield(wtinfo(wave), 'wtprefix');
0112 
0113 %-Compute & store contrast parameters, contrast/ESS images, & stat images
0114 %=======================================================================
0115 spm('Pointer','Watch')
0116 spm_progress_bar('Init',100,'computing...')                          %-#
0117 nPar   = size(xX.X,2);
0118 I      = unique(connos);
0119 for ii = 1:length(I)
0120   
0121   i  = I(ii);
0122 
0123   % Check exists
0124   if i > length(xCon)
0125     warning(['Contrast ' num2str(i) ' does not exist']);
0126     break;
0127   end
0128   
0129   %-Canonicalise contrast structure with required fields
0130   %-------------------------------------------------------------------
0131   eidf = mars_struct('getifthere', xCon(i), 'eidf');
0132   if isempty(eidf)
0133     X1o           = spm_FcUtil('X1o',xCon(i),xX.xKXs);
0134     [trMV,trMVMV] = spm_SpUtil('trMV',X1o,xX.V);
0135     xCon(i).eidf = trMV^2/trMVMV;
0136   else
0137     trMV = []; trMVMV = [];
0138   end
0139 
0140   % make sensible file suffix from contrast name
0141   fsuff = mars_utils('str2fname', xCon(i).name);
0142   
0143   %-Write contrast/ESS images?
0144   %-------------------------------------------------------------------
0145   Vcon = full_vol(phiwD, ...
0146           mars_struct('getifthere', xCon(i), 'Vcon'));
0147   if isempty(Vcon.dim)
0148 
0149     % We're going to change the contrast structure
0150     changef = 1;
0151     
0152     %-Bomb out (nicely) if can't write to results directory
0153     %---------------------------------------------------------------
0154     if ~wOK, spm('alert*',{    'Can''t write to the results directory:',...
0155             ['        ',Swd],' ','=> can''t compute new contrasts'},...
0156          mfilename,sqrt(-1));
0157       spm('Pointer','Arrow')
0158       error('can''t write contrast image')
0159     end
0160         
0161     switch(xCon(i).STAT)
0162       
0163      case 'T'       %-Implement contrast as sum of scaled beta images
0164                
0165       fprintf('\t%-32s: %-10s%20s',sprintf('contrast image %2d',i),...
0166           '(spm_add)','...initialising') %-#
0167       
0168       Q     = find(abs(xCon(i).c) > 0);
0169       V     = Vbeta(Q);
0170       for j = 1:length(Q)
0171     V(j).pinfo(1,:) = V(j).pinfo(1,:)*xCon(i).c(Q(j));
0172       end
0173       
0174       %-Prepare handle for contrast image
0175       %-----------------------------------------------------------
0176       Vcon = struct(...
0177       'fname',  fullfile(Swd,sprintf('%scon_%04d_%s.img',wvp,i,fsuff)),...
0178       'dim',    [DIM,16],...
0179       'mat',    M,...
0180       'pinfo',  [1,0,0]',...
0181       'descrip',sprintf('Phiwave contrast - %d: %s',i,xCon(i).name));
0182       
0183       %-Write image
0184       %-----------------------------------------------------------
0185       fprintf('%s%30s',bs30,'...computing')%-#
0186       Vcon            = spm_create_image(Vcon);
0187       Vcon.pinfo(1,1) = spm_add(V, Vcon);
0188       Vcon            = spm_close_vol(Vcon);
0189             
0190       fprintf('%s%30s\n',bs30,sprintf(...
0191       '...written %s',spm_str_manip(Vcon.fname,'t')))%-#
0192      case 'F'  %-Implement ESS as sum of squared weighted beta images
0193       fprintf('\t%-32s: %30s',sprintf('ESS image %2d',i),...
0194           '...computing') %-#
0195       
0196       %-Residual (in parameter space) forming mtx
0197       %-----------------------------------------------------------
0198       h       = spm_FcUtil('Hsqr',xCon(i),xX.xKXs);
0199       
0200       %-Prepare handle for ESS image
0201       %-----------------------------------------------------------
0202       Vcon = struct(...
0203       'fname',  fullfile(Swd,sprintf('%sess_%04d_%s.img',wvp,i,fsuff)),...
0204       'dim',    [DIM 16],...
0205       'mat',    M,...
0206       'pinfo',  [1,0,0]',...
0207       'descrip',sprintf('Phiwave ESS - contrast %d: %s',i,xCon(i).name));
0208       
0209       %-Write image
0210       %-----------------------------------------------------------
0211       fprintf('%s',bs30)                   %-#
0212       Vcon  = spm_create_image(Vcon);
0213       Vcon  = spm_resss(Vbeta, Vcon, h);
0214       Vcon  = spm_close_vol(Vcon);
0215       
0216      otherwise
0217       %---------------------------------------------------------------
0218       error(['unknown STAT "',xCon(i).STAT,'"'])
0219       
0220     end % (switch(xCon...)
0221     
0222     % Put into structure array
0223     xCon(i).Vcon = design_vol(phiwD, Vcon);
0224 
0225     % Write wave info for file
0226     putwave(Vcon.fname,wave);
0227         
0228   end % (if isempty...)
0229 
0230   spm_progress_bar('Set',100*(2*ii-1)/(2*length(I)+2))             %-#
0231   
0232   %-Write statistic image(s)
0233   %-------------------------------------------------------------------
0234   Vspm = full_vol(phiwD, ...
0235           mars_struct('getifthere', xCon(i), 'Vspm'));
0236   if ~flags.con_only & isempty(Vspm.dim)
0237     
0238     % We're going to change the contrast structure
0239     changef = 1;
0240     
0241     % Read Residual mean squared image if necessary
0242     if isempty(rmsi)
0243       fprintf('\t%-32s: %30s','ResMS file...','...done');
0244       rmsi = spm_read_vols(VResMS);
0245       fprintf('%s%30s\n',bs30,'...done')               %-#
0246       rmsi(abs(rmsi)<eps) = NaN;
0247     end
0248 
0249     %-Bomb out (nicely) if can't write to results directory
0250     %---------------------------------------------------------------
0251     if ~wOK, spm('alert*',{    'Can''t write to the results directory:',...
0252             ['        ',Swd],' ','=> can''t compute new contrasts'},...
0253          mfilename,sqrt(-1));
0254       spm('Pointer','Arrow')
0255       error('can''t write Phiwave image')
0256     end
0257     
0258     fprintf('\t%-32s: %30s',sprintf('phiw{%c} image %2d',xCon(i).STAT,i),...
0259         '...computing')  %-#
0260     
0261     switch(xCon(i).STAT)
0262      case 'T'                                  %-Compute {t} image
0263       
0264       Z   = spm_read_vols(Vcon)./...
0265         (sqrt(rmsi * (xCon(i).c'*xX.Bcov*xCon(i).c) ));
0266       
0267       str = sprintf('[%.2g]',erdf);
0268       
0269      case 'F'                                  %-Compute {F} image
0270       
0271       if isempty(trMV)
0272     trMV = spm_SpUtil('trMV',spm_FcUtil('X1o',xCon(i),xX.xKXs),xX.V);
0273       end
0274       Z =(spm_read_vols(Vcon)/trMV)./rmsi;
0275       
0276       str = sprintf('[%.2g,%.2g]',xCon(i).eidf,erdf);
0277             
0278      otherwise
0279       %---------------------------------------------------------------
0280       error(['unknown STAT "',xCon(i).STAT,'"'])
0281     end
0282         
0283     %-Write full statistic image
0284     %---------------------------------------------------------------
0285     fprintf('%s%30s',bs30,'...writing')      %-#
0286     Vspm = struct(...
0287     'fname',  fullfile(Swd,sprintf('%sphiw%c_%04d_%s.img',wvp,xCon(i).STAT,i,fsuff)),...
0288     'dim',    [DIM 16],...
0289     'mat',    M,...
0290     'pinfo',  [1,0,0]',...
0291     'descrip',sprintf('Phiwave{%c_%s} - contrast %d: %s',...
0292               xCon(i).STAT,str,i,xCon(i).name));
0293         
0294     Vspm       = spm_write_vol(Vspm,Z);
0295     
0296     % Write wave part of mat file
0297     putwave(Vspm.fname,wave);
0298     
0299     fprintf('%s%30s\n',bs30,sprintf(...
0300     '...written %s',spm_str_manip(Vspm.fname,'t')))  %-#
0301         
0302   end % (if ~flags.con_only & isemptu...)
0303   
0304   % Put into structure array
0305   xCon(i).Vspm = design_vol(phiwD, Vspm);
0306 
0307   % Write wave info for file
0308   putwave(Vspm.fname,wave);
0309         
0310   spm_progress_bar('Set',100*(2*ii-0)/(2*length(I)+2))             %-#
0311 
0312 end % (for ii = 1:length(I))
0313 
0314 spm_progress_bar('Set',100)                                          %-#
0315 
0316 %- put contrasts back into design
0317 %=======================================================================
0318 phiwD = set_contrasts(phiwD, xCon);
0319 
0320 spm_progress_bar('Clear')                                          %-#
0321 
0322 return
0323 
0324 % thank you, really
write_contrasts

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
