thresh_calc

 calculate thresholding from wt'ed statistic volumes
 FORMAT [th_obj, dndescrip] = thresh_calc(wtobj,errobj,statinf,dninf)

 Threshold calculation assumes that a t statistic is made by dividing the
 data in the phiw_wvimg object wtobj by that in errobj.  We need both the
 numerator and the denominator because we may want to pool the error over
 levels.

 Inputs 
   wtobj   - top half of t statistic
   errobj  - bottom half (error) of t statistic
   statinf - structure with information on generated statistic
             containing fields
             .stat  - statistic type 'T' or 'Z' or 'TZ' (t to Z
             transform)
             .df    - degrees of freedom for T statistic
   dninf    - structure with information on denoising to be done
             containing fields
             .levels  - levels to threshold; should be a vector of
                length wtobj.scales+1, with a value for each scale: 
                value = 1 - apply thresholding to this level
                value = 0 - suppress all coeffs for this level
                value = -Inf - omit thresholding for this level
             .thlev - level of thresholding; one of:
                'image' - whole image (excluding omitted etc levels)
                'level' - level by level
                'quadrant' - by quadrant
             .ncalc - calculation of number of null hypotheses per thlev
                 block; can be:
                 'n'      - number of non zero coefficients
                 'pplot'  - number calculated using pplot algorithm
             .thcalc - threshold calculation method  
             'visu':    MinMax threshold
          'sure':    SURE thresholding (needs soft thresholding)
          'stein':      linear Stein
          'bonf':    Bonferroni correction
          'fdr':    false discovery rate correction
          'hoch':    Hoch correction
             .thapp - threshold application method
                 'hard'        all cooefficients under thresh = 0
                 'soft'        sign(y)(abs(y)-thresh)_+
                 'linear'      y*thresh
             .alpha - alpha level for p value threshold methods
             .varpoolf - flag, if non zero, apply variance pooling

 Outputs
   th_obj     - wv_img object with thresholding
   dndescrip  - text description of denoising applied

 Matthew Brett 2/6/2001, Federico E. Turkheimer 17/9/2000

 $Id: thresh_calc.m,v 1.6 2005/06/05 04:36:38 matthewbrett Exp $

0001 function [th_obj, dndescrip] = thresh_calc(wtobj,errobj,statinf,dninf)
0002 % calculate thresholding from wt'ed statistic volumes
0003 % FORMAT [th_obj, dndescrip] = thresh_calc(wtobj,errobj,statinf,dninf)
0004 %
0005 % Threshold calculation assumes that a t statistic is made by dividing the
0006 % data in the phiw_wvimg object wtobj by that in errobj.  We need both the
0007 % numerator and the denominator because we may want to pool the error over
0008 % levels.
0009 %
0010 % Inputs
0011 %   wtobj   - top half of t statistic
0012 %   errobj  - bottom half (error) of t statistic
0013 %   statinf - structure with information on generated statistic
0014 %             containing fields
0015 %             .stat  - statistic type 'T' or 'Z' or 'TZ' (t to Z
0016 %             transform)
0017 %             .df    - degrees of freedom for T statistic
0018 %   dninf    - structure with information on denoising to be done
0019 %             containing fields
0020 %             .levels  - levels to threshold; should be a vector of
0021 %                length wtobj.scales+1, with a value for each scale:
0022 %                value = 1 - apply thresholding to this level
0023 %                value = 0 - suppress all coeffs for this level
0024 %                value = -Inf - omit thresholding for this level
0025 %             .thlev - level of thresholding; one of:
0026 %                'image' - whole image (excluding omitted etc levels)
0027 %                'level' - level by level
0028 %                'quadrant' - by quadrant
0029 %             .ncalc - calculation of number of null hypotheses per thlev
0030 %                 block; can be:
0031 %                 'n'      - number of non zero coefficients
0032 %                 'pplot'  - number calculated using pplot algorithm
0033 %             .thcalc - threshold calculation method
0034 %             'visu':    MinMax threshold
0035 %          'sure':    SURE thresholding (needs soft thresholding)
0036 %          'stein':      linear Stein
0037 %          'bonf':    Bonferroni correction
0038 %          'fdr':    false discovery rate correction
0039 %          'hoch':    Hoch correction
0040 %             .thapp - threshold application method
0041 %                 'hard'        all cooefficients under thresh = 0
0042 %                 'soft'        sign(y)(abs(y)-thresh)_+
0043 %                 'linear'      y*thresh
0044 %             .alpha - alpha level for p value threshold methods
0045 %             .varpoolf - flag, if non zero, apply variance pooling
0046 %
0047 % Outputs
0048 %   th_obj     - wv_img object with thresholding
0049 %   dndescrip  - text description of denoising applied
0050 %
0051 % Matthew Brett 2/6/2001, Federico E. Turkheimer 17/9/2000
0052 %
0053 % $Id: thresh_calc.m,v 1.6 2005/06/05 04:36:38 matthewbrett Exp $
0054 
0055 if nargin < 2
0056   error('Need at least top and bottom of t statistic, sorry');
0057 end
0058 if nargin < 3
0059   statinf = [];
0060 end
0061 if nargin < 4
0062   dninf = [];
0063 end
0064 
0065 % load objects
0066 wtobj = doproc(wtobj);
0067 errobj = doproc(errobj);
0068 
0069 % Make object for return
0070 th_obj = wtobj;
0071 
0072 % check stat structure
0073 defstat = struct('stat','Z','df',10000);
0074 statinf = mars_struct('ffillsplit', defstat,  statinf);
0075 
0076 % check denoise structure
0077 defdn = struct('levels',ones(1,wtobj.scales+1),...
0078            'thlev','level',...
0079            'ncalc','n',...
0080            'thcalc','sure',...
0081            'thapp','soft',...
0082            'varpoolf',0,...
0083            'alpha',0.05);
0084 dninf = mars_struct('ffillsplit', defdn, dninf);
0085 if length(dninf.levels) ~= wtobj.scales+1
0086   error('Unexpected length of levels spec')
0087 end
0088 if strcmp(dninf.thcalc,'sure') & ~strcmp(dninf.thapp,'soft')
0089   error('sure calculation requires soft thresholding')
0090 end
0091 % apply level thresholding options if present
0092 if isfield(dninf,'lettop') & ~isempty(dninf.lettop) & dninf.lettop > 0
0093   dninf.levels((end-dninf.lettop+1):end) = -Inf;
0094 end
0095 if isfield(dninf,'killbot') & ~isempty(dninf.killbot) & dninf.killbot > 0
0096   dninf.levels(1:(dninf.killbot)) = 0;
0097 end
0098 
0099 % remove NaNs
0100 wtobj.img(isnan(wtobj.img)) = 0;
0101 
0102 % apply levels matrix
0103 suplevs = find(dninf.levels==0);
0104 if ~isempty(suplevs)
0105   wtobj = subsasgn(wtobj,struct('type','()','subs',{{suplevs}}),0);
0106 end
0107 letlevs = find(dninf.levels==-Inf);
0108 dolevs = find(dninf.levels==1);
0109 
0110 % determine blocks
0111 quads = 1:7;
0112 switch dninf.thlev
0113  case 'image'
0114   l{1} = dolevs;
0115   q{1} = quads;
0116  case 'level'
0117   l = num2cell(dolevs);
0118   q{1} = quads;
0119  case 'quadrant'
0120   l = num2cell(dolevs);
0121   q = num2cell(quads);
0122  otherwise
0123   error('Don''t recognize level specification');
0124 end  
0125 
0126 % number of blocks
0127 lq = length(q);
0128 nblks = length(l) * lq;
0129 % take into account fewer blocks for top level
0130 if lq > 1 & any(dolevs == wtobj.scales+1)
0131   nblks = nblks - lq + 1;
0132 end
0133 
0134 % alpha level Bonferroni corrected for no of blocks
0135 alpha = dninf.alpha / nblks;
0136 
0137 % flag for p value calculation
0138 if strcmp(dninf.ncalc,'pplot') | ...
0139   any(strcmp(dninf.thcalc,{'hoch','fdr'}))
0140   pvalf = 1;
0141 else
0142   pvalf = 0;
0143 end
0144 
0145 % cycle over blocks to do denoising
0146 for li = 1:length(l)
0147   if l{li} > wtobj.scales
0148     % at top level -> only one quadrant
0149     q = {1};
0150   end
0151   for qi = 1:length(q)
0152     s = struct('type','()','subs',{{l{li},q{qi}}});
0153     % get data
0154     dblk = subsref(wtobj, s);
0155     eblk = subsref(errobj, s);
0156     full_thblk = zeros(size(dblk));
0157     
0158     idx = isfinite(dblk) & dblk~= 0 & isfinite(eblk) & eblk ~=0;
0159     if ~any(idx), break,end
0160     eblk = eblk(idx);
0161     
0162     % implement pooled variance estimate
0163     if dninf.varpoolf
0164       ncoeff = length(eblk);
0165       eblk  = sqrt(sum(eblk.^2) / ncoeff);
0166       df = statinf.df * ncoeff;
0167     else
0168       df = statinf.df;
0169     end
0170     
0171     % Process on stat type
0172     stat = statinf.stat;
0173     stblk = dblk(idx) ./ eblk;
0174     if ~strcmp(stat, 'Z') & df > 1000, stat = 'Z'; end    
0175     if strcmp(stat, 'TZ')
0176       % do T to Z transform
0177       stblk = spm_t2z(stblk, df);
0178     end
0179     
0180     % calculate statistic and p values
0181     if pvalf
0182       if strcmp(stat, 'T')
0183     pblk = spm_Tcdf(-abs(stblk),df)*2;
0184       else
0185     pblk = spm_Ncdf(-abs(stblk))*2;
0186       end    
0187       [pblk pidx] = sort(pblk);
0188     end
0189     
0190     % do null coeff calc
0191     switch dninf.ncalc
0192      case 'n'
0193       n = length(stblk);
0194      case 'pplot'
0195       n = pr_pplot(pblk,alpha,1);
0196      otherwise
0197       error('Don''t recognize ncalc method')
0198     end
0199     
0200     % find thresholds
0201     switch dninf.thcalc
0202      case 'sure'
0203       crit = log2(n).^(3/2)./sqrt(n);
0204       ss       =   sum(stblk.^2-1)/n;
0205       if(ss<=crit)
0206     thresh =sqrt(2*log(n));
0207       else
0208     thresh =min([pr_sure_thresh(stblk') sqrt(2*log(n))]); 
0209       end 
0210      case 'stein'                         
0211       thresh = 1-n/sum(stblk.^2);
0212       thresh = thresh *sign(thresh);
0213      case 'visu'
0214       thresh = sqrt(2*log(n));
0215      case 'bonf'
0216       if strcmp(stat,'T')
0217     thresh = -spm_invTcdf(alpha/2/n,df);
0218       else
0219     thresh = -spm_invNcdf(alpha/2/n);
0220       end
0221      case 'hoch'
0222       H = pr_hoch(pblk,alpha,n,1);
0223       thresh = min([Inf,min(abs(stblk(pidx(logical(H)))))]);
0224      case 'fdr'
0225       H = pr_fdr(pblk,alpha,n,1);
0226       thresh = min([Inf,min(abs(stblk(pidx(logical(H)))))]);
0227      otherwise 
0228       error('Don''t recognize threshold calculation')
0229     end
0230     
0231     % apply thresholding to t image
0232     switch dninf.thapp 
0233      case 'linear'
0234       thblk = stblk * thresh;
0235      case 'soft'
0236       thblk = sign(stblk) .* max(0,(abs(stblk)-thresh));
0237      case 'hard'
0238       thblk = stblk;
0239       thblk(abs(thblk) < thresh) = 0;
0240      otherwise
0241       error('Don''t recognize threshold application type')
0242     end
0243     full_thblk(idx) = thblk ./ stblk;
0244     
0245     % restore data to object
0246     th_obj = subsasgn(th_obj,s,full_thblk);
0247     
0248   end
0249 end
0250 
0251 % description of denoising
0252 dndescrip = ['Denoising by ' dninf.thlev ...
0253          '; null hypothesis calc = ' dninf.ncalc ...
0254          '; threshold calc = ' dninf.thcalc ...;
0255          '; threshold app = ' dninf.thapp ];
0256 if any(strcmp(dninf.thcalc,{'fdr','hoch','bonf'}))
0257   dndescrip = strvcat(dndescrip, ['Alpha level: ' num2str(dninf.alpha)]);
0258 end
0259 if ~isempty(suplevs)
0260   strvcat(dndescrip, [num2str(suplevs) ' levels suppressed']);
0261 end
0262 if ~isempty(letlevs)
0263   strvcat(dndescrip, [num2str(letlevs) ' levels not thresholded']);
0264 end
0265 
0266 th_obj.descrip = strvcat(th_obj.descrip,dndescrip);