thresholdByFDR

 thresholdByFDR - threshold a stat map using the False Discovery Rate (FDR)

 Usage: thresholdMap = thresholdByFDR(flag, statMap, FDR, signed, npk)

 Thresholds an input stat map for use in creating an ROI.
 flag is any valid fsfast stat map (currently, only implimented
 for 'sig' and 't', and 'sig' has some bugs).  statMap is the
 stat map in the form numSlices by x by y.  Threshold is a number to
 threshold the stat map, which will be different for different types
 of stat maps.  Finally, signed says whether or not the thresholding 
 should be done in a signed or unsigned way (i.e., with actual values
 vs. absolute values).  Returns the threshold map with 1s in the thresholded
 voxels and 0s otherwise.

 This function is part of froi, available from http://froi.sourceforge.net,
 and is governed by the terms of the Artistic License.

 $Id: thresholdByFDR.m,v 1.2 2003/09/23 15:38:19 nknouf Exp $

0001 function thresholdMap = thresholdByFDR(flag, statMap, FDR, signed, np)
0002 % thresholdByFDR - threshold a stat map using the False Discovery Rate (FDR)
0003 %
0004 % Usage: thresholdMap = thresholdByFDR(flag, statMap, FDR, signed, npk)
0005 %
0006 % Thresholds an input stat map for use in creating an ROI.
0007 % flag is any valid fsfast stat map (currently, only implimented
0008 % for 'sig' and 't', and 'sig' has some bugs).  statMap is the
0009 % stat map in the form numSlices by x by y.  Threshold is a number to
0010 % threshold the stat map, which will be different for different types
0011 % of stat maps.  Finally, signed says whether or not the thresholding
0012 % should be done in a signed or unsigned way (i.e., with actual values
0013 % vs. absolute values).  Returns the threshold map with 1s in the thresholded
0014 % voxels and 0s otherwise.
0015 %
0016 % This function is part of froi, available from http://froi.sourceforge.net,
0017 % and is governed by the terms of the Artistic License.
0018 %
0019 % $Id: thresholdByFDR.m,v 1.2 2003/09/23 15:38:19 nknouf Exp $
0020 
0021 %% create basic structure to hold data
0022 stat.statMap = statMap;
0023 stat.FDR = FDR;
0024 stat.np = np;
0025 
0026 %% statMap is a three element vector (numSlices by x by y)
0027 [stat.slices, stat.xsize, stat.ysize] = size(stat.statMap);
0028 
0029 %% create an empty thresholdMap the same size as the input stat map
0030 thresholdMap = zeros(stat.slices, stat.xsize, stat.ysize);
0031 
0032 %% run different thresholding code depending on the input flag
0033 switch (flag)
0034     case {'-sig'},
0035         thresholdMap = quickFDR(stat);
0036     case {'-t'},
0037         thresholdMap = slowFDR(flag, stat);
0038 end
0039 
0040 %% FDR calculations when we already have a sig map
0041 function thresholdMap = quickFDR(stat)
0042     statmap1d = reshape1d(stat.statMap);
0043 
0044     [pID, pN] = FDR(statmap1d, stat.FDR);
0045 
0046     if (stat.np)
0047         %% using non-parametric calculation
0048         p = pN;
0049         fprintf('p value from FDR calculation (non-parametric): %f\n', p);
0050     else
0051         p = pID;
0052         fprintf('p value from FDR calculation (with independence and/or positive dependence assumption): %f\n', p);
0053     end
0054 
0055     thresholdMap = stat.statMap > (-log(p));
0056 
0057         %% check and see if we have anything in the resulting map
0058         if(isempty(find(thresholdMap==1)))
0059                 oops(250);
0060         end
0061 return;
0062 
0063 %% FDR calculations when we need to convert some other sort of map to
0064 %% a sig map
0065 function thresholdMap = slowFDR(flag, stat);
0066     %% ensure that all values are finite
0067     %%stat.statMap(~isfinite(stat.statMap))) = [];
0068 
0069     statmap1d = reshape1d(stat.statMap);
0070 
0071 return;