Thresholds and P-values of peaks and clusters of random fields in any D.
[PEAK_THRESHOLD, EXTENT_THRESHOLD, PEAK_THRESHOLD_1 EXTENT_THRESHOLD_1] =
STAT_THRESHOLD([ SEARCH_VOLUME [, NUM_VOXELS [, FWHM [, DF [, P_VAL_PEAK
[, CLUSTER_THRESHOLD [, P_VAL_EXTENT [, NCONJ [, NVAR [, EC_FILE
[, TRANSFORM [, NPRINT ]]]]]]]]]]]] );
If P-values are supplied, returns the thresholds for local maxima
or peaks (PEAK_THRESHOLD) (if thresholds are supplied then it
returns P-values) for the following SPMs:
- Z, chi^2, t, F, Hotelling's T^2, Roy's maximum root, maximum
canonical correlation, cross- and auto-correlation;
- scale space Z and chi^2;
- conjunctions of NCONJ of these (except cross- and auto-correlation);
and spatial extent of contiguous voxels above CLUSTER_THRESHOLD
(EXTENT_THRESHOLD) except conjunctions. Note that for the multivariate
statistics spatial extent is measured in resels for the higher
dimensional space formed by adding sphere(s) to the search region.
For peaks (local maxima), two methods are used:
random field theory (Worsley et al. 1996, Human Brain Mapping, 4:58-73),
and a simple Bonferroni correction based on the number of voxels
in the search region. The final threshold is the minimum of the two.
For clusters, the method of Cao, 1999, Advances in Applied Probability,
31:577-593 is used. The cluster size is only accurate for large
CLUSTER_THRESHOLD (say >3) and large resels = SEARCH_VOLUME / FWHM^D.
PEAK_THRESHOLD_1 is the height of a single peak chosen in advance, and
EXTENT_THRESHOLD_1 is the extent of a single cluster chosen in advance
of looking at the data, e.g. the nearest peak or cluster to a pre-chosen
voxel or ROI - see Friston KJ. Testing for anatomically specified
regional effects. Human Brain Mapping. 5(2):133-6, 1997.
SEARCH_VOLUME is the volume of the search region in mm^3. The method for
finding PEAK_THRESHOLD works well for any value of SEARCH_VOLUME, even
SEARCH_VOLUME = 0 (default), which gives the threshold for the image at
a single voxel. The random field theory threshold is based on the
assumption that the search region is a sphere in 3D, which is a very
tight lower bound for any non-spherical region.
For a non-spherical D-dimensional region, set SEARCH_VOLUME to a vector
of the D+1 intrinsic volumes of the search region. In D=3 these are
[Euler characteristic, 2 * caliper diameter, 1/2 surface area, volume].
E.g. for a sphere of radius r in 3D, use [1, 4*r, 2*pi*r^2, 4/3*pi*r^3],
which is equivalent to just SEARCH_VOLUME = 4/3*pi*r^3.
For a 2D search region, use [1, 1/2 perimeter length, area].
For cross- and auto-correlations, where correlation is maximized over all
pairs of voxels from two search regions, SEARCH_VOLUME has two rows
interpreted as above - see Cao, J. & Worsley, K.J. (1999). The geometry
of correlation fields, with an application to functional connectivity of
the brain. Annals of Applied Probability, 9:1021-1057.
NUM_VOXELS is the number of voxels (3D) or pixels (2D) in the search
volume. For cross- and auto-correlations, use two rows. Default is 1.
FWHM is the fwhm in mm of a smoothing kernel applied to the data. Default
is 0.0, i.e. no smoothing, which is roughly the case for raw fMRI data.
For motion corrected fMRI data, use at least 6mm; for PET data, this would
be the scanner fwhm of about 6mm. Using the geometric mean of the three
x,y,z FWHM's is a very good approximation if they are different.
If you have already calculated resels, then set SEARCH_VOLUME=resels
and FWHM=1. Cluster extents will then be in resels, rather than mm^3.
For cross- and auto-correlations, use two rows. For scale space, use two
columns for the min and max fwhm (only for Z and chi^2 SPMs).
DF=[DF1 DF2; DFW1 DFW2 0] is a 2 x 2 matrix of degrees of freedom.
If DF2 is 0, then DF1 is the df of the T statistic image.
If DF1=Inf then it calculates thresholds for the Gaussian image.
If DF2>0 then DF1 and DF2 are the df's of the F statistic image.
DFW1 and DFW2 are the numerator and denominator df of the FWHM image.
They are only used for calculating cluster resel p-values and thresholds
(ignored if NVAR>1 or NCONJ>1 since there are no theoretical results).
The random estimate of the local resels adds variability to the summed
resels of the cluster. The distribution of the local resels summed over a
region depends on the unknown spatial correlation structure of the resels,
so we assume that the resels are constant over the cluster, reasonable if the
clusters are small. The cluster resels are then multiplied by the random resel
distribution at a point. This gives an upper bound to p-values and thresholds.
If DF=[DF1 DF2] (row) then DFW1=DFW2=Inf, i.e. FWHM is fixed.
If DF=[DF1; DFW1] (column) then DF2=0 and DFW2=DFW1, i.e. t statistic.
If DF=DF1 (scalar) then DF2=0 and DFW1=DFW2=Inf, i.e. t stat, fixed FWHM.
If any component of DF >= 1000 then it is set to Inf. Default is Inf.
P_VAL_PEAK is the row vector of desired P-values for peaks.
If the first element is greater than 1, then they are
treated as peak values and P-values are returned. Default is 0.05.
To get P-values for peak values < 1, set the first element > 1,
set the second to the deisired peak value, then discard the first result.
CLUSTER_THRESHOLD is the scalar threshold of the image for clusters.
If it is <= 1, it is taken as a probability p, and the
cluter_thresh is chosen so that P( T > cluster_thresh ) = p.
Default is 0.001, i.e. P( T > cluster_thresh ) = 0.001.
P_VAL_EXTENT is the row vector of desired P-values for spatial extents
of clusters of contiguous voxels above CLUSTER_THRESHOLD.
If the first element is greater than 1, then they are treated as
extent values and P-values are returned. Default is 0.05.
NCONJ is the number of conjunctions. If NCONJ > 1, calculates P-values and
thresholds for peaks (but not clusters) of the minimum of NCONJ independent
SPM's - see Friston, K.J., Holmes, A.P., Price, C.J., Buchel, C.,
Worsley, K.J. (1999). Multi-subject fMRI studies and conjunction analyses.
NeuroImage, 10:385-396. Default is NCONJ = 1 (no conjunctions).
NVAR is the number of variables for multivariate equivalents of T and F
statistics, found by maximizing T^2 or F over all linear combinations of
variables, i.e. Hotelling's T^2 for DF1=1, Roy's maximum root for DF1>1.
For maximum canonical cross- and auto-correlations, use two rows.
See Worsley, K.J., Taylor, J.E., Tomaiuolo, F. & Lerch, J. (2004). Unified
univariate and multivariate random field theory. Neuroimage, 23:S189-195.
Default is 1, i.e. univariate statistics.
EC_FILE: file for EC densities in IRIS Explorer latin module format.
Ignored if empty (default).
EXPR: matlab expression applied to threshold in t, e.g. 't=sqrt(t);'.
NPRINT: maximum number of P-values or thresholds to print. Default 5.
Examples:
T statistic, 1000000mm^3 search volume, 1000000 voxels (1mm voxel volume),
20mm smoothing, 30 degrees of freedom, P=0.05 and 0.01 for peaks, cluster
threshold chosen so that P=0.001 at a single voxel, P=0.05 and 0.01 for extent:
stat_threshold(1000000,1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]);
peak_threshold = 5.0826 5.7853
Cluster_threshold = 3.3852
peak_threshold_1 = 4.8782 5.5955
extent_threshold = 3340.3 6315.5
extent_threshold_1 = 2911.5 5719.6
Check: Suppose we are given peak values of 5.0826, 5.7853 and
spatial extents of 3340.3, 6315.5 above a threshold of 3.3852,
then we should get back our original P-values:
stat_threshold(1000000,1000000,20,30,[5.0826 5.7853],3.3852,[3340.3 6315.5]);
P_val_peak = 0.0500 0.0100
P_val_peak_1 = 0.0318 0.0065
P_val_extent = 0.0500 0.0100
P_val_extent_1 = 0.0379 0.0074
Another way of doing this is to use the fact that T^2 has an F
distribution with 1 and 30 degrees of freedom, and double the P values:
stat_threshold(1000000,1000000,20,[1 30],[0.1 0.02],0.002,[0.1 0.02]);
peak_threshold = 25.8330 33.4702
Cluster_threshold = 11.4596
peak_threshold_1 = 20.7886 27.9741
extent_threshold = 3297.8 6305.2
extent_threshold_1 = 1936.4 4420.3
Note that sqrt([25.8330 33.4702 11.4596])=[5.0826 5.7853 3.3852]
in agreement with the T statistic thresholds, but the spatial extent
thresholds are very close but not exactly the same.
If the shape of the search region is known, then the 'intrinisic
volume' must be supplied. E.g. for a spherical search region with
a volume of 1000000 mm^3, radius r:
r=(1000000/(4/3*pi))^(1/3);
stat_threshold([1 4*r 2*pi*r^2 4/3*pi*r^3], ...
1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]);
gives the same results as our first example. A 100 x 100 x 100 cube
stat_threshold([1 300 30000 1000000], ...
1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]);
gives slightly higher thresholds (5.0969, 5.7976) for peaks, but the cluster
thresholds are the same since they depend (so far) only on the volume and
not the shape. For a 2D circular search region of the same radius r, use
stat_threshold([1 pi*r pi*r^2],1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]);
A volume of 0 returns the usual uncorrected threshold for peaks,
which can also be found by setting NUM_VOXELS=1, FWHM = 0:
stat_threshold(0,1,20,30,[0.05 0.01])
stat_threshold(1,1, 0,30,[0.05 0.01]);
For non-isotropic fields, replace the SEARCH_VOLUME by the vector of resels
(see Worsley et al. 1996, Human Brain Mapping, 4:58-73), and set FWHM=1,
so that EXTENT_THRESHOLD's are measured in resels, not mm^3. If the
resels are estimated from residuals, add an extra row to DF (see above).
For the conjunction of 2 and 3 T fields as above:
stat_threshold(1000000,1000000,20,30,[0.05 0.01],[],[],2)
stat_threshold(1000000,1000000,20,30,[0.05 0.01],[],[],3)
returns lower thresholds of [3.0251 3.3984] and [2.2336 2.5141] resp. Note
that there are as yet no theoretical results for extents so they are NaN.
For Hotelling's T^2 e.g. for linear models of deformations (NVAR=3):
stat_threshold(1000000,1000000,20,[1 30],[0.05 0.01],[],[],1,3)
For Roy's max root, e.g. for effective connectivity using deformations,
stat_threshold(1000000,1000000,20,[3 30],[0.05 0.01],[],[],1,3)
There are no theoretical results for mulitvariate extents so they are NaN.
Check: A Hotelling's T^2 with Inf df is the same as a chi^2 with NVAR df
stat_threshold(1000000,1000000,20,[1 Inf],[],[],[],[],NVAR)
stat_threshold(1000000,1000000,20,[NVAR Inf])*NVAR
should give the same answer!
Cross- and auto-correlations: to threshold the auto-correlation of fmrilm
residuals (100 df) over all pairs of 1000000 voxels in a 1000000mm^3
region with 20mm FWHM smoothing, use df=99 (one less for the correlation):
stat_threshold([1000000; 1000000], [1000000; 1000000], 20, 99, 0.05)
which gives a T statistic threshold of 6.7923 with 99 df. To convert to
a correlation, use 6.7923/sqrt(99+6.7923^2)=0.5638. Note that auto-
correlations are symmetric, so the P-value is in fact 0.05/2=0.025; on the
other hand, we usually want a two sided test of both positive and negative
correlations, so the P-value should be doubled back to 0.05. For cross-
correlations, e.g. between n=321 GM density volumes (1000000mm^3 ball,
FWHM=13mm) and cortical thickness on the average surface (250000mm^2,
FWHM=18mm), use 321-2=319 df for removing the constant and the correlation.
Note that we use P=0.025 to threshold both +ve and -ve correlations.
r=(1000000/(4/3*pi))^(1/3);
stat_threshold([1 4*r 2*pi*r^2 4/3*pi*r^3; 1 0 250000 0], ...
[1000000; 40962], [13; 18], 319, 0.025)
This gives a T statistic threshold of 6.7158 with 319 df. To convert to
a correlation, use 6.7158/sqrt(319+6.7158^2)=0.3520. Note that NVAR can be
greater than one for maximum canonical cross- and auto-correlations
between all pairs of multivariate imaging data. Finally, the spatial
extent for 5D clusters of connected correlations is 5.9228e+006 mm^5.
Scale space: suppose we smooth 1000000 voxels in a 1000000mm^3 region
from 6mm to 30mm in 10 steps ( e.g. fwhm=6*(30/6).^((0:9)/9) ):
stat_threshold(1000000, 1000000*10, [6 30])
gives a scale space threshold of 5.0663, only slightly higher than the
5.0038 you get from using the highest resolution data only, i.e.
stat_threshold(1000000, 1000000, 6)0001 function [ peak_threshold, extent_threshold, ... 0002 peak_threshold_1, extent_threshold_1, t, rho ] = ... 0003 stat_threshold( search_volume, num_voxels, fwhm, df, p_val_peak, ... 0004 cluster_threshold, p_val_extent, nconj, nvar, EC_file, expr, nprint ); 0005 0006 %Thresholds and P-values of peaks and clusters of random fields in any D. 0007 % 0008 % [PEAK_THRESHOLD, EXTENT_THRESHOLD, PEAK_THRESHOLD_1 EXTENT_THRESHOLD_1] = 0009 % STAT_THRESHOLD([ SEARCH_VOLUME [, NUM_VOXELS [, FWHM [, DF [, P_VAL_PEAK 0010 % [, CLUSTER_THRESHOLD [, P_VAL_EXTENT [, NCONJ [, NVAR [, EC_FILE 0011 % [, TRANSFORM [, NPRINT ]]]]]]]]]]]] ); 0012 % 0013 % If P-values are supplied, returns the thresholds for local maxima 0014 % or peaks (PEAK_THRESHOLD) (if thresholds are supplied then it 0015 % returns P-values) for the following SPMs: 0016 % - Z, chi^2, t, F, Hotelling's T^2, Roy's maximum root, maximum 0017 % canonical correlation, cross- and auto-correlation; 0018 % - scale space Z and chi^2; 0019 % - conjunctions of NCONJ of these (except cross- and auto-correlation); 0020 % and spatial extent of contiguous voxels above CLUSTER_THRESHOLD 0021 % (EXTENT_THRESHOLD) except conjunctions. Note that for the multivariate 0022 % statistics spatial extent is measured in resels for the higher 0023 % dimensional space formed by adding sphere(s) to the search region. 0024 % 0025 % For peaks (local maxima), two methods are used: 0026 % random field theory (Worsley et al. 1996, Human Brain Mapping, 4:58-73), 0027 % and a simple Bonferroni correction based on the number of voxels 0028 % in the search region. The final threshold is the minimum of the two. 0029 % 0030 % For clusters, the method of Cao, 1999, Advances in Applied Probability, 0031 % 31:577-593 is used. The cluster size is only accurate for large 0032 % CLUSTER_THRESHOLD (say >3) and large resels = SEARCH_VOLUME / FWHM^D. 0033 % 0034 % PEAK_THRESHOLD_1 is the height of a single peak chosen in advance, and 0035 % EXTENT_THRESHOLD_1 is the extent of a single cluster chosen in advance 0036 % of looking at the data, e.g. the nearest peak or cluster to a pre-chosen 0037 % voxel or ROI - see Friston KJ. Testing for anatomically specified 0038 % regional effects. Human Brain Mapping. 5(2):133-6, 1997. 0039 % 0040 % SEARCH_VOLUME is the volume of the search region in mm^3. The method for 0041 % finding PEAK_THRESHOLD works well for any value of SEARCH_VOLUME, even 0042 % SEARCH_VOLUME = 0 (default), which gives the threshold for the image at 0043 % a single voxel. The random field theory threshold is based on the 0044 % assumption that the search region is a sphere in 3D, which is a very 0045 % tight lower bound for any non-spherical region. 0046 % For a non-spherical D-dimensional region, set SEARCH_VOLUME to a vector 0047 % of the D+1 intrinsic volumes of the search region. In D=3 these are 0048 % [Euler characteristic, 2 * caliper diameter, 1/2 surface area, volume]. 0049 % E.g. for a sphere of radius r in 3D, use [1, 4*r, 2*pi*r^2, 4/3*pi*r^3], 0050 % which is equivalent to just SEARCH_VOLUME = 4/3*pi*r^3. 0051 % For a 2D search region, use [1, 1/2 perimeter length, area]. 0052 % For cross- and auto-correlations, where correlation is maximized over all 0053 % pairs of voxels from two search regions, SEARCH_VOLUME has two rows 0054 % interpreted as above - see Cao, J. & Worsley, K.J. (1999). The geometry 0055 % of correlation fields, with an application to functional connectivity of 0056 % the brain. Annals of Applied Probability, 9:1021-1057. 0057 % 0058 % NUM_VOXELS is the number of voxels (3D) or pixels (2D) in the search 0059 % volume. For cross- and auto-correlations, use two rows. Default is 1. 0060 % 0061 % FWHM is the fwhm in mm of a smoothing kernel applied to the data. Default 0062 % is 0.0, i.e. no smoothing, which is roughly the case for raw fMRI data. 0063 % For motion corrected fMRI data, use at least 6mm; for PET data, this would 0064 % be the scanner fwhm of about 6mm. Using the geometric mean of the three 0065 % x,y,z FWHM's is a very good approximation if they are different. 0066 % If you have already calculated resels, then set SEARCH_VOLUME=resels 0067 % and FWHM=1. Cluster extents will then be in resels, rather than mm^3. 0068 % For cross- and auto-correlations, use two rows. For scale space, use two 0069 % columns for the min and max fwhm (only for Z and chi^2 SPMs). 0070 % 0071 % DF=[DF1 DF2; DFW1 DFW2 0] is a 2 x 2 matrix of degrees of freedom. 0072 % If DF2 is 0, then DF1 is the df of the T statistic image. 0073 % If DF1=Inf then it calculates thresholds for the Gaussian image. 0074 % If DF2>0 then DF1 and DF2 are the df's of the F statistic image. 0075 % DFW1 and DFW2 are the numerator and denominator df of the FWHM image. 0076 % They are only used for calculating cluster resel p-values and thresholds 0077 % (ignored if NVAR>1 or NCONJ>1 since there are no theoretical results). 0078 % The random estimate of the local resels adds variability to the summed 0079 % resels of the cluster. The distribution of the local resels summed over a 0080 % region depends on the unknown spatial correlation structure of the resels, 0081 % so we assume that the resels are constant over the cluster, reasonable if the 0082 % clusters are small. The cluster resels are then multiplied by the random resel 0083 % distribution at a point. This gives an upper bound to p-values and thresholds. 0084 % If DF=[DF1 DF2] (row) then DFW1=DFW2=Inf, i.e. FWHM is fixed. 0085 % If DF=[DF1; DFW1] (column) then DF2=0 and DFW2=DFW1, i.e. t statistic. 0086 % If DF=DF1 (scalar) then DF2=0 and DFW1=DFW2=Inf, i.e. t stat, fixed FWHM. 0087 % If any component of DF >= 1000 then it is set to Inf. Default is Inf. 0088 % 0089 % P_VAL_PEAK is the row vector of desired P-values for peaks. 0090 % If the first element is greater than 1, then they are 0091 % treated as peak values and P-values are returned. Default is 0.05. 0092 % To get P-values for peak values < 1, set the first element > 1, 0093 % set the second to the deisired peak value, then discard the first result. 0094 % 0095 % CLUSTER_THRESHOLD is the scalar threshold of the image for clusters. 0096 % If it is <= 1, it is taken as a probability p, and the 0097 % cluter_thresh is chosen so that P( T > cluster_thresh ) = p. 0098 % Default is 0.001, i.e. P( T > cluster_thresh ) = 0.001. 0099 % 0100 % P_VAL_EXTENT is the row vector of desired P-values for spatial extents 0101 % of clusters of contiguous voxels above CLUSTER_THRESHOLD. 0102 % If the first element is greater than 1, then they are treated as 0103 % extent values and P-values are returned. Default is 0.05. 0104 % 0105 % NCONJ is the number of conjunctions. If NCONJ > 1, calculates P-values and 0106 % thresholds for peaks (but not clusters) of the minimum of NCONJ independent 0107 % SPM's - see Friston, K.J., Holmes, A.P., Price, C.J., Buchel, C., 0108 % Worsley, K.J. (1999). Multi-subject fMRI studies and conjunction analyses. 0109 % NeuroImage, 10:385-396. Default is NCONJ = 1 (no conjunctions). 0110 % 0111 % NVAR is the number of variables for multivariate equivalents of T and F 0112 % statistics, found by maximizing T^2 or F over all linear combinations of 0113 % variables, i.e. Hotelling's T^2 for DF1=1, Roy's maximum root for DF1>1. 0114 % For maximum canonical cross- and auto-correlations, use two rows. 0115 % See Worsley, K.J., Taylor, J.E., Tomaiuolo, F. & Lerch, J. (2004). Unified 0116 % univariate and multivariate random field theory. Neuroimage, 23:S189-195. 0117 % Default is 1, i.e. univariate statistics. 0118 % 0119 % EC_FILE: file for EC densities in IRIS Explorer latin module format. 0120 % Ignored if empty (default). 0121 % 0122 % EXPR: matlab expression applied to threshold in t, e.g. 't=sqrt(t);'. 0123 % 0124 % NPRINT: maximum number of P-values or thresholds to print. Default 5. 0125 % 0126 % Examples: 0127 % 0128 % T statistic, 1000000mm^3 search volume, 1000000 voxels (1mm voxel volume), 0129 % 20mm smoothing, 30 degrees of freedom, P=0.05 and 0.01 for peaks, cluster 0130 % threshold chosen so that P=0.001 at a single voxel, P=0.05 and 0.01 for extent: 0131 % 0132 % stat_threshold(1000000,1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]); 0133 % 0134 % peak_threshold = 5.0826 5.7853 0135 % Cluster_threshold = 3.3852 0136 % peak_threshold_1 = 4.8782 5.5955 0137 % extent_threshold = 3340.3 6315.5 0138 % extent_threshold_1 = 2911.5 5719.6 0139 % 0140 % Check: Suppose we are given peak values of 5.0826, 5.7853 and 0141 % spatial extents of 3340.3, 6315.5 above a threshold of 3.3852, 0142 % then we should get back our original P-values: 0143 % 0144 % stat_threshold(1000000,1000000,20,30,[5.0826 5.7853],3.3852,[3340.3 6315.5]); 0145 % 0146 % P_val_peak = 0.0500 0.0100 0147 % P_val_peak_1 = 0.0318 0.0065 0148 % P_val_extent = 0.0500 0.0100 0149 % P_val_extent_1 = 0.0379 0.0074 0150 % 0151 % Another way of doing this is to use the fact that T^2 has an F 0152 % distribution with 1 and 30 degrees of freedom, and double the P values: 0153 % 0154 % stat_threshold(1000000,1000000,20,[1 30],[0.1 0.02],0.002,[0.1 0.02]); 0155 % 0156 % peak_threshold = 25.8330 33.4702 0157 % Cluster_threshold = 11.4596 0158 % peak_threshold_1 = 20.7886 27.9741 0159 % extent_threshold = 3297.8 6305.2 0160 % extent_threshold_1 = 1936.4 4420.3 0161 % 0162 % Note that sqrt([25.8330 33.4702 11.4596])=[5.0826 5.7853 3.3852] 0163 % in agreement with the T statistic thresholds, but the spatial extent 0164 % thresholds are very close but not exactly the same. 0165 % 0166 % If the shape of the search region is known, then the 'intrinisic 0167 % volume' must be supplied. E.g. for a spherical search region with 0168 % a volume of 1000000 mm^3, radius r: 0169 % 0170 % r=(1000000/(4/3*pi))^(1/3); 0171 % stat_threshold([1 4*r 2*pi*r^2 4/3*pi*r^3], ... 0172 % 1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]); 0173 % 0174 % gives the same results as our first example. A 100 x 100 x 100 cube 0175 % 0176 % stat_threshold([1 300 30000 1000000], ... 0177 % 1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]); 0178 % 0179 % gives slightly higher thresholds (5.0969, 5.7976) for peaks, but the cluster 0180 % thresholds are the same since they depend (so far) only on the volume and 0181 % not the shape. For a 2D circular search region of the same radius r, use 0182 % 0183 % stat_threshold([1 pi*r pi*r^2],1000000,20,30,[0.05 0.01],0.001,[0.05 0.01]); 0184 % 0185 % A volume of 0 returns the usual uncorrected threshold for peaks, 0186 % which can also be found by setting NUM_VOXELS=1, FWHM = 0: 0187 % 0188 % stat_threshold(0,1,20,30,[0.05 0.01]) 0189 % stat_threshold(1,1, 0,30,[0.05 0.01]); 0190 % 0191 % For non-isotropic fields, replace the SEARCH_VOLUME by the vector of resels 0192 % (see Worsley et al. 1996, Human Brain Mapping, 4:58-73), and set FWHM=1, 0193 % so that EXTENT_THRESHOLD's are measured in resels, not mm^3. If the 0194 % resels are estimated from residuals, add an extra row to DF (see above). 0195 % 0196 % For the conjunction of 2 and 3 T fields as above: 0197 % 0198 % stat_threshold(1000000,1000000,20,30,[0.05 0.01],[],[],2) 0199 % stat_threshold(1000000,1000000,20,30,[0.05 0.01],[],[],3) 0200 % 0201 % returns lower thresholds of [3.0251 3.3984] and [2.2336 2.5141] resp. Note 0202 % that there are as yet no theoretical results for extents so they are NaN. 0203 % 0204 % For Hotelling's T^2 e.g. for linear models of deformations (NVAR=3): 0205 % 0206 % stat_threshold(1000000,1000000,20,[1 30],[0.05 0.01],[],[],1,3) 0207 % 0208 % For Roy's max root, e.g. for effective connectivity using deformations, 0209 % 0210 % stat_threshold(1000000,1000000,20,[3 30],[0.05 0.01],[],[],1,3) 0211 % 0212 % There are no theoretical results for mulitvariate extents so they are NaN. 0213 % Check: A Hotelling's T^2 with Inf df is the same as a chi^2 with NVAR df 0214 % 0215 % stat_threshold(1000000,1000000,20,[1 Inf],[],[],[],[],NVAR) 0216 % stat_threshold(1000000,1000000,20,[NVAR Inf])*NVAR 0217 % 0218 % should give the same answer! 0219 % 0220 % Cross- and auto-correlations: to threshold the auto-correlation of fmrilm 0221 % residuals (100 df) over all pairs of 1000000 voxels in a 1000000mm^3 0222 % region with 20mm FWHM smoothing, use df=99 (one less for the correlation): 0223 % 0224 % stat_threshold([1000000; 1000000], [1000000; 1000000], 20, 99, 0.05) 0225 % 0226 % which gives a T statistic threshold of 6.7923 with 99 df. To convert to 0227 % a correlation, use 6.7923/sqrt(99+6.7923^2)=0.5638. Note that auto- 0228 % correlations are symmetric, so the P-value is in fact 0.05/2=0.025; on the 0229 % other hand, we usually want a two sided test of both positive and negative 0230 % correlations, so the P-value should be doubled back to 0.05. For cross- 0231 % correlations, e.g. between n=321 GM density volumes (1000000mm^3 ball, 0232 % FWHM=13mm) and cortical thickness on the average surface (250000mm^2, 0233 % FWHM=18mm), use 321-2=319 df for removing the constant and the correlation. 0234 % Note that we use P=0.025 to threshold both +ve and -ve correlations. 0235 % 0236 % r=(1000000/(4/3*pi))^(1/3); 0237 % stat_threshold([1 4*r 2*pi*r^2 4/3*pi*r^3; 1 0 250000 0], ... 0238 % [1000000; 40962], [13; 18], 319, 0.025) 0239 % 0240 % This gives a T statistic threshold of 6.7158 with 319 df. To convert to 0241 % a correlation, use 6.7158/sqrt(319+6.7158^2)=0.3520. Note that NVAR can be 0242 % greater than one for maximum canonical cross- and auto-correlations 0243 % between all pairs of multivariate imaging data. Finally, the spatial 0244 % extent for 5D clusters of connected correlations is 5.9228e+006 mm^5. 0245 % 0246 % Scale space: suppose we smooth 1000000 voxels in a 1000000mm^3 region 0247 % from 6mm to 30mm in 10 steps ( e.g. fwhm=6*(30/6).^((0:9)/9) ): 0248 % 0249 % stat_threshold(1000000, 1000000*10, [6 30]) 0250 % 0251 % gives a scale space threshold of 5.0663, only slightly higher than the 0252 % 5.0038 you get from using the highest resolution data only, i.e. 0253 % 0254 % stat_threshold(1000000, 1000000, 6) 0255 0256 %############################################################################ 0257 % COPYRIGHT: Copyright 2003 K.J. Worsley 0258 % Department of Mathematics and Statistics, 0259 % McConnell Brain Imaging Center, 0260 % Montreal Neurological Institute, 0261 % McGill University, Montreal, Quebec, Canada. 0262 % keith.worsley@mcgill.ca , www.math.mcgill.ca/keith 0263 % 0264 % Permission to use, copy, modify, and distribute this 0265 % software and its documentation for any purpose and without 0266 % fee is hereby granted, provided that this copyright 0267 % notice appears in all copies. The author and McGill University 0268 % make no representations about the suitability of this 0269 % software for any purpose. It is provided "as is" without 0270 % express or implied warranty. 0271 %############################################################################ 0272 0273 % Defaults: 0274 if nargin<1; search_volume=[]; end 0275 if nargin<2; num_voxels=[]; end 0276 if nargin<3; fwhm=[]; end 0277 if nargin<4; df=[]; end 0278 if nargin<5; p_val_peak=[]; end 0279 if nargin<6; cluster_threshold=[]; end 0280 if nargin<7; p_val_extent=[]; end 0281 if nargin<8; nconj=[]; end 0282 if nargin<9; nvar=[]; end 0283 if nargin<10; EC_file=[]; end 0284 if nargin<11; expr=[]; end 0285 if nargin<12; nprint=[]; end 0286 0287 if isempty(search_volume); search_volume=0; end 0288 if isempty(num_voxels); num_voxels=1; end 0289 if isempty(fwhm); fwhm=0.0; end 0290 if isempty(df); df=Inf; end 0291 if isempty(p_val_peak); p_val_peak=0.05; end 0292 if isempty(cluster_threshold); cluster_threshold=0.001; end 0293 if isempty(p_val_extent); p_val_extent=0.05; end 0294 if isempty(nconj); nconj=1; end 0295 if isempty(nvar); nvar=1; end 0296 if isempty(nprint); nprint=5; end 0297 0298 if size(fwhm,1)==1; fwhm(2,:)=fwhm; end 0299 if size(fwhm,2)==1; scale=1; else scale=fwhm(1,2)/fwhm(1,1); fwhm=fwhm(:,1); end; 0300 isscale=(scale>1); 0301 0302 if length(num_voxels)==1; num_voxels(2,1)=1; end 0303 0304 if size(search_volume,2)==1 0305 radius=(search_volume/(4/3*pi)).^(1/3); 0306 search_volume=[ones(length(radius),1) 4*radius 2*pi*radius.^2 search_volume]; 0307 end 0308 if size(search_volume,1)==1 0309 search_volume=[search_volume; [1 zeros(1,size(search_volume,2)-1)]]; 0310 end 0311 lsv=size(search_volume,2); 0312 fwhm_inv=all(fwhm>0)./(fwhm+any(fwhm<=0)); 0313 resels=search_volume.*repmat(fwhm_inv,1,lsv).^repmat(0:lsv-1,2,1); 0314 invol=resels.*(4*log(2)).^(repmat(0:lsv-1,2,1)/2); 0315 for k=1:2 0316 D(k,1)=max(find(invol(k,:)))-1; 0317 end 0318 0319 % determines which method was used to estimate fwhm (see fmrilm or multistat): 0320 df_limit=4; 0321 0322 if length(df)==1; df=[df 0]; end 0323 if size(df,1)==1; df=[df; Inf Inf; Inf Inf]; end 0324 if size(df,2)==1; df=[df df]; df(1,2)=0; end 0325 if size(df,1)==2; df=[df; df(2,:)]; end 0326 0327 % is_tstat=1 if it is a t statistic 0328 is_tstat=(df(1,2)==0); 0329 if is_tstat 0330 df1=1; 0331 df2=df(1,1); 0332 else 0333 df1=df(1,1); 0334 df2=df(1,2); 0335 end 0336 if df2 >= 1000; df2=Inf; end 0337 df0=df1+df2; 0338 0339 dfw1=df(2:3,1); 0340 dfw2=df(2:3,2); 0341 for k=1:2 0342 if dfw1(k) >= 1000; dfw1(k)=Inf; end 0343 if dfw2(k) >= 1000; dfw2(k)=Inf; end 0344 end 0345 0346 if length(nvar)==1; nvar(2,1)=df1; end 0347 0348 if isscale & (D(2)>1 | nvar(1,1)>1 | df2<Inf) 0349 D 0350 nvar 0351 df2 0352 fprintf('Cannot do scale space.'); 0353 return 0354 end 0355 0356 Dlim=D+[scale>1; 0]; 0357 DD=Dlim+nvar-1; 0358 0359 % Values of the F statistic: 0360 t=((1000:-1:1)'/100).^4; 0361 % Find the upper tail probs cumulating the F density using Simpson's rule: 0362 if df2==Inf 0363 u=df1*t; 0364 b=exp(-u/2-log(2*pi)/2+log(u)/4)*df1^(1/4)*4/100; 0365 else 0366 u=df1*t/df2; 0367 b=exp(-df0/2*log(1+u)+log(u)/4-betaln(1/2,(df0-1)/2))*(df1/df2)^(1/4)*4/100; 0368 end 0369 t=[t; 0]; 0370 b=[b; 0]; 0371 n=length(t); 0372 sb=cumsum(b); 0373 sb1=cumsum(b.*(-1).^(1:n)'); 0374 pt1=sb+sb1/3-b/3; 0375 pt2=sb-sb1/3-b/3; 0376 tau=zeros(n,DD(1)+1,DD(2)+1); 0377 tau(1:2:n,1,1)=pt1(1:2:n); 0378 tau(2:2:n,1,1)=pt2(2:2:n); 0379 tau(n,1,1)=1; 0380 tau(:,1,1)=min(tau(:,1,1),1); 0381 0382 % Find the EC densities: 0383 u=df1*t; 0384 for d=1:max(DD) 0385 for e=0:min(min(DD),d) 0386 s1=0; 0387 cons=-((d+e)/2+1)*log(pi)+gammaln(d)+gammaln(e+1); 0388 for k=0:(d-1+e)/2 0389 [i,j]=ndgrid(0:k,0:k); 0390 if df2==Inf 0391 q1=log(pi)/2-((d+e-1)/2+i+j)*log(2); 0392 else 0393 q1=(df0-1-d-e)*log(2)+gammaln((df0-d)/2+i)+gammaln((df0-e)/2+j) ... 0394 -gammalni(df0-d-e+i+j+k)-((d+e-1)/2-k)*log(df2); 0395 end 0396 q2=cons-gammalni(i+1)-gammalni(j+1)-gammalni(k-i-j+1) ... 0397 -gammalni(d-k-i+j)-gammalni(e-k-j+i+1); 0398 s2=sum(sum(exp(q1+q2))); 0399 if s2>0 0400 s1=s1+(-1)^k*u.^((d+e-1)/2-k)*s2; 0401 end 0402 end 0403 if df2==Inf 0404 s1=s1.*exp(-u/2); 0405 else 0406 s1=s1.*exp(-(df0-2)/2*log(1+u/df2)); 0407 end 0408 if DD(1)>=DD(2) 0409 tau(:,d+1,e+1)=s1; 0410 if d<=min(DD) 0411 tau(:,e+1,d+1)=s1; 0412 end 0413 else 0414 tau(:,e+1,d+1)=s1; 0415 if d<=min(DD) 0416 tau(:,d+1,e+1)=s1; 0417 end 0418 end 0419 end 0420 end 0421 0422 % For multivariate statistics, add a sphere to the search region: 0423 a=zeros(2,max(nvar)); 0424 for k=1:2 0425 j=(nvar(k)-1):-2:0; 0426 a(k,j+1)=exp(j*log(2)+j/2*log(pi) ... 0427 +gammaln((nvar(k)+1)/2)-gammaln((nvar(k)+1-j)/2)-gammaln(j+1)); 0428 end 0429 rho=zeros(n,Dlim(1)+1,Dlim(2)+1); 0430 for k=1:nvar(1) 0431 for l=1:nvar(2) 0432 rho=rho+a(1,k)*a(2,l)*tau(:,(0:Dlim(1))+k,(0:Dlim(2))+l); 0433 end 0434 end 0435 0436 if is_tstat 0437 if nvar==1 0438 t=[sqrt(t(1:(n-1))); -flipdim(sqrt(t),1)]; 0439 rho=[rho(1:(n-1),:,:); flipdim(rho,1)]/2; 0440 for i=0:D(1) 0441 for j=0:D(2) 0442 rho(n-1+(1:n),i+1,j+1)=-(-1)^(i+j)*rho(n-1+(1:n),i+1,j+1); 0443 end 0444 end 0445 rho(n-1+(1:n),1,1)=rho(n-1+(1:n),1,1)+1; 0446 n=2*n-1; 0447 else 0448 t=sqrt(t); 0449 end 0450 end 0451 0452 % For scale space: 0453 if scale>1 0454 kappa=D(1)/2; 0455 tau=zeros(n,D(1)+1); 0456 for d=0:D(1) 0457 s1=0; 0458 for k=0:d/2 0459 s1=s1+(-1)^k/(1-2*k)*exp(gammaln(d+1)-gammaln(k+1)-gammaln(d-2*k+1) ... 0460 +(1/2-k)*log(kappa)-k*log(4*pi))*rho(:,d+2-2*k,1); 0461 end 0462 if d==0 0463 cons=log(scale); 0464 else 0465 cons=(1-1/scale^d)/d; 0466 end 0467 tau(:,d+1)=rho(:,d+1,1)*(1+1/scale^d)/2+s1*cons; 0468 end 0469 rho(:,1:(D(1)+1),1)=tau; 0470 end 0471 0472 if D(2)==0 0473 d=D(1); 0474 if nconj>1 0475 % Conjunctions: 0476 b=gamma(((0:d)+1)/2)/gamma(1/2); 0477 for i=1:d+1 0478 rho(:,i,1)=rho(:,i,1)/b(i); 0479 end 0480 m1=zeros(n,d+1,d+1); 0481 for i=1:d+1 0482 j=i:d+1; 0483 m1(:,i,j)=rho(:,j-i+1,1); 0484 end 0485 for k=2:nconj 0486 for i=1:d+1 0487 for j=1:d+1 0488 m2(:,i,j)=sum(rho(:,1:d+2-i,1).*m1(:,i:d+1,j),2); 0489 end 0490 end 0491 m1=m2; 0492 end 0493 for i=1:d+1 0494 rho(:,i,1)=m1(:,1,i)*b(i); 0495 end 0496 end 0497 0498 if ~isempty(EC_file) 0499 if d<3 0500 rho(:,(d+2):4,1)=zeros(n,4-d-2+1); 0501 end 0502 fid=fopen(EC_file,'w'); 0503 % first 3 are dimension sizes as 4-byte integers: 0504 fwrite(fid,[n max(d+2,5) 1],'int'); 0505 % next 6 are bounding box as 4-byte floats: 0506 fwrite(fid,[0 0 0; 1 1 1],'float'); 0507 % rest are the data as 4-byte floats: 0508 if ~isempty(expr) 0509 eval(expr); 0510 end 0511 fwrite(fid,t,'float'); 0512 fwrite(fid,rho,'float'); 0513 fclose(fid); 0514 end 0515 end 0516 0517 if all(fwhm>0) 0518 pval_rf=zeros(n,1); 0519 for i=1:D(1)+1 0520 for j=1:D(2)+1 0521 pval_rf=pval_rf+invol(1,i)*invol(2,j)*rho(:,i,j); 0522 end 0523 end 0524 else 0525 pval_rf=Inf; 0526 end 0527 0528 % Bonferroni 0529 pt=rho(:,1,1); 0530 pval_bon=abs(prod(num_voxels))*pt; 0531 0532 % Minimum of the two: 0533 pval=min(pval_rf,pval_bon); 0534 0535 tlim=1; 0536 if p_val_peak(1) <= tlim 0537 peak_threshold=minterp1(pval,t,p_val_peak); 0538 if length(p_val_peak)<=nprint 0539 peak_threshold 0540 end 0541 else 0542 % p_val_peak is treated as a peak value: 0543 P_val_peak=interp1(t,pval,p_val_peak); 0544 i=isnan(P_val_peak); 0545 P_val_peak(i)=(is_tstat & (p_val_peak(i)<0)); 0546 peak_threshold=P_val_peak; 0547 if length(p_val_peak)<=nprint 0548 P_val_peak 0549 end 0550 end 0551 0552 if all(fwhm<=0) | any(num_voxels<0) 0553 peak_threshold_1=p_val_peak+NaN; 0554 extent_threshold=p_val_extent+NaN; 0555 extent_threshold_1=extent_threshold; 0556 return 0557 end 0558 0559 % Cluster_threshold: 0560 0561 if cluster_threshold > tlim 0562 tt=cluster_threshold; 0563 else 0564 % cluster_threshold is treated as a probability: 0565 tt=minterp1(pt,t,cluster_threshold); 0566 if nprint>0 0567 Cluster_threshold=tt 0568 end 0569 end 0570 0571 d=sum(D); 0572 rhoD=interp1(t,rho(:,D(1)+1,D(2)+1),tt); 0573 p=interp1(t,pt,tt); 0574 0575 % Pre-selected peak: 0576 0577 pval=rho(:,D(1)+1,D(2)+1)./rhoD; 0578 0579 if p_val_peak(1) <= tlim 0580 peak_threshold_1=minterp1(pval,t, p_val_peak); 0581 if length(p_val_peak)<=nprint 0582 peak_threshold_1 0583 end 0584 else 0585 % p_val_peak is treated as a peak value: 0586 P_val_peak_1=interp1(t,pval,p_val_peak); 0587 i=isnan(P_val_peak_1); 0588 P_val_peak_1(i)=(is_tstat & (p_val_peak(i)<0)); 0589 peak_threshold_1=P_val_peak_1; 0590 if length(p_val_peak)<=nprint 0591 P_val_peak_1 0592 end 0593 end 0594 0595 if d==0 | nconj>1 | nvar(1)>1 | scale>1 0596 extent_threshold=p_val_extent+NaN; 0597 extent_threshold_1=extent_threshold; 0598 if length(p_val_extent)<=nprint 0599 extent_threshold 0600 extent_threshold_1 0601 end 0602 return 0603 end 0604 0605 % Expected number of clusters: 0606 0607 EL=invol(1,D(1)+1)*invol(2,D(2)+1)*rhoD; 0608 cons=gamma(d/2+1)*(4*log(2))^(d/2)/fwhm(1)^D(1)/fwhm(2)^D(2)*rhoD/p; 0609 0610 if df2==Inf & dfw1(1)==Inf & dfw1(2)==Inf 0611 if p_val_extent(1) <= tlim 0612 pS=-log(1-p_val_extent)/EL; 0613 extent_threshold=(-log(pS)).^(d/2)/cons; 0614 pS=-log(1-p_val_extent); 0615 extent_threshold_1=(-log(pS)).^(d/2)/cons; 0616 if length(p_val_extent)<=nprint 0617 extent_threshold 0618 extent_threshold_1 0619 end 0620 else 0621 % p_val_extent is now treated as a spatial extent: 0622 pS=exp(-(p_val_extent*cons).^(2/d)); 0623 P_val_extent=1-exp(-pS*EL); 0624 extent_threshold=P_val_extent; 0625 P_val_extent_1=1-exp(-pS); 0626 extent_threshold_1=P_val_extent_1; 0627 if length(p_val_extent)<=nprint 0628 P_val_extent 0629 P_val_extent_1 0630 end 0631 end 0632 else 0633 % Find dbn of S by taking logs then using fft for convolution: 0634 ny=2^12; 0635 a=d/2; 0636 b2=a*10*max(sqrt(2/(min(df1+df2,min(dfw1)))),1); 0637 if df2<Inf 0638 b1=a*log((1-(1-0.000001)^(2/(df2-d)))*df2/2); 0639 else 0640 b1=a*log(-log(1-0.000001)); 0641 end 0642 dy=(b2-b1)/ny; 0643 b1=round(b1/dy)*dy; 0644 y=((1:ny)'-1)*dy+b1; 0645 numrv=1+(d+(D(1)>0)+(D(2)>0))*(df2<Inf)+... 0646 (D(1)*(dfw1(1)<Inf)+(dfw2(1)<Inf))*(D(1)>0)+...; 0647 (D(2)*(dfw1(2)<Inf)+(dfw2(2)<Inf))*(D(2)>0); 0648 f=zeros(ny,numrv); 0649 mu=zeros(1,numrv); 0650 if df2<Inf 0651 % Density of log(Beta(1,(df2-d)/2)^(d/2)): 0652 yy=exp(y./a)/df2*2; 0653 yy=yy.*(yy<1); 0654 f(:,1)=(1-yy).^((df2-d)/2-1).*((df2-d)/2).*yy/a; 0655 mu(1)=exp(gammaln(a+1)+gammaln((df2-d+2)/2)-gammaln((df2+2)/2)+a*log(df2/2)); 0656 else 0657 % Density of log(exp(1)^(d/2)): 0658 yy=exp(y./a); 0659 f(:,1)=exp(-yy).*yy/a; 0660 mu(1)=exp(gammaln(a+1)); 0661 end 0662 0663 nuv=[]; 0664 aav=[]; 0665 if df2<Inf 0666 nuv=df2+2-(1:d); 0667 aav=[repmat(-1/2,1,d)]; 0668 for k=1:2 0669 if D(k)>0; 0670 nuv=[df1+df2-D(k) nuv]; 0671 aav=[D(k)/2 aav]; 0672 end; 0673 end; 0674 end 0675 0676 for k=1:2 0677 if dfw1(k)<Inf & D(k)>0 0678 if dfw1(k)>df_limit 0679 nuv=[nuv dfw1(k)-dfw1(k)/dfw2(k)-(0:(D(k)-1))]; 0680 else 0681 nuv=[nuv repmat(dfw1(k)-dfw1(k)/dfw2(k),1,D(k))]; 0682 end 0683 aav=[aav repmat(1/2,1,D(k))]; 0684 end 0685 if dfw2(k)<Inf 0686 nuv=[nuv dfw2(k)]; 0687 aav=[aav -D(k)/2]; 0688 end 0689 end 0690 0691 for i=1:(numrv-1) 0692 nu=nuv(i); 0693 aa=aav(i); 0694 yy=y/aa+log(nu); 0695 % Density of log((chi^2_nu/nu)^aa): 0696 f(:,i+1)=exp(nu/2*yy-exp(yy)/2-(nu/2)*log(2)-gammaln(nu/2))/abs(aa); 0697 mu(i+1)=exp(gammaln(nu/2+aa)-gammaln(nu/2)-aa*log(nu/2)); 0698 end 0699 % Check: plot(y,f); sum(f*dy,1) should be 1 0700 0701 omega=2*pi*((1:ny)'-1)/ny/dy; 0702 shift=complex(cos(-b1*omega),sin(-b1*omega))*dy; 0703 prodfft=prod(fft(f),2).*shift.^(numrv-1); 0704 % Density of Y=log(B^(d/2)*U^(d/2)/sqrt(det(Q))): 0705 ff=real(ifft(prodfft)); 0706 % Check: plot(y,ff); sum(ff*dy) should be 1 0707 mu0=prod(mu); 0708 % Check: plot(y,ff.*exp(y)); sum(ff.*exp(y)*dy.*(y<10)) should equal mu0 0709 0710 alpha=p/rhoD/mu0*fwhm(1)^D(1)*fwhm(2)^D(2)/(4*log(2))^(d/2); 0711 0712 % Integrate the density to get the p-value for one cluster: 0713 pS=cumsum(ff(ny:-1:1))*dy; 0714 pS=pS(ny:-1:1); 0715 % The number of clusters is Poisson with mean EL: 0716 pSmax=1-exp(-pS*EL); 0717 0718 if p_val_extent(1) <= tlim 0719 yval=minterp1(-pSmax,y,-p_val_extent); 0720 % Spatial extent is alpha*exp(Y) -dy/2 correction for mid-point rule: 0721 extent_threshold=alpha*exp(yval-dy/2); 0722 % For a single cluster: 0723 yval=minterp1(-pS,y,-p_val_extent); 0724 extent_threshold_1=alpha*exp(yval-dy/2); 0725 if length(p_val_extent)<=nprint 0726 extent_threshold 0727 extent_threshold_1 0728 end 0729 else 0730 % p_val_extent is now treated as a spatial extent: 0731 logpval=log(p_val_extent/alpha+(p_val_extent<=0))+dy/2; 0732 P_val_extent=interp1(y,pSmax,logpval); 0733 extent_threshold=P_val_extent.*(p_val_extent>0)+(p_val_extent<=0); 0734 % For a single cluster: 0735 P_val_extent_1=interp1(y,pS,logpval); 0736 extent_threshold_1=P_val_extent_1.*(p_val_extent>0)+(p_val_extent<=0); 0737 if length(p_val_extent)<=nprint 0738 P_val_extent 0739 P_val_extent_1 0740 end 0741 end 0742 0743 end 0744 0745 return 0746 0747 function x=gammalni(n); 0748 i=find(n>=0); 0749 x=Inf+n; 0750 if ~isempty(i) 0751 x(i)=gammaln(n(i)); 0752 end 0753 return 0754 0755 function iy=minterp1(x,y,ix); 0756 % interpolates only the monotonically increasing values of x at ix 0757 n=length(x); 0758 mx=x(1); 0759 my=y(1); 0760 xx=x(1); 0761 for i=2:n 0762 if x(i)>xx 0763 xx=x(i); 0764 mx=[mx xx]; 0765 my=[my y(i)]; 0766 end 0767 end 0768 iy=interp1(mx,my,ix); 0769 return