Estimation of smoothness based on [residual] images
FORMAT [fwhm,VRpv] = spm_est_smoothness(VResI,[VM]);
P - filenames of [residual] images
PM - filename of mask image
fwhm - estimated FWHM in all image directions
VRpv - handle of Resels per Voxel image
_______________________________________________________________________
spm_est_smoothness returns a spatial smoothness estimator based on the
varianic_es of the normalized spatial derivatives as described in K.
Worsley, (1996). Inputs are a mask image and a number of [residual]
images. Output is a global estimate of the smoothness expressed as the
FWHM of an equivalent Gaussian point spread function. An estimate of
resels per voxels (see spm_spm) is written as an image file ('RPV.img')
to the current directory.
The mask image specifies voxels, used in smoothness estimation, by
assigning them non-zero values. The dimensions, voxel sizes, orientation
of all images must be the same. The dimensions of the images can be of
dimensions 0, 1, 2 and 3.
Note that 1-dim images (lines) must exist in the 1st dimension and
2-dim images (slices) in the first two dimensions. The estimated fwhm
for any non-existing dimension is infinity.
Ref:
K. Worsley (1996). An unbiased estimator for the roughness of a
multivariate Gaussian random field. Technical Report, Department of
Mathematics and Statistics, McGill University
_______________________________________________________________________
Copyright (C) 2005 Wellcome Department of Imaging Neuroscienic_e0001 function [fwhm,VRpv] = nic_spm_est_smoothness(varargin) 0002 % Estimation of smoothness based on [residual] images 0003 % FORMAT [fwhm,VRpv] = spm_est_smoothness(VResI,[VM]); 0004 % 0005 % P - filenames of [residual] images 0006 % PM - filename of mask image 0007 % 0008 % fwhm - estimated FWHM in all image directions 0009 % VRpv - handle of Resels per Voxel image 0010 %_______________________________________________________________________ 0011 % 0012 % spm_est_smoothness returns a spatial smoothness estimator based on the 0013 % varianic_es of the normalized spatial derivatives as described in K. 0014 % Worsley, (1996). Inputs are a mask image and a number of [residual] 0015 % images. Output is a global estimate of the smoothness expressed as the 0016 % FWHM of an equivalent Gaussian point spread function. An estimate of 0017 % resels per voxels (see spm_spm) is written as an image file ('RPV.img') 0018 % to the current directory. 0019 % 0020 % The mask image specifies voxels, used in smoothness estimation, by 0021 % assigning them non-zero values. The dimensions, voxel sizes, orientation 0022 % of all images must be the same. The dimensions of the images can be of 0023 % dimensions 0, 1, 2 and 3. 0024 % 0025 % Note that 1-dim images (lines) must exist in the 1st dimension and 0026 % 2-dim images (slices) in the first two dimensions. The estimated fwhm 0027 % for any non-existing dimension is infinity. 0028 % 0029 % 0030 % Ref: 0031 % 0032 % K. Worsley (1996). An unbiased estimator for the roughness of a 0033 % multivariate Gaussian random field. Technical Report, Department of 0034 % Mathematics and Statistics, McGill University 0035 % 0036 %_______________________________________________________________________ 0037 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscienic_e 0038 0039 % Stefan Kiebel 0040 % $Id: spm_est_smoothness.m 112 2005-05-04 18:20:52Z john $ 0041 0042 0043 % assign input argumants 0044 %----------------------------------------------------------------------- 0045 if nargin > 0, V = varargin{1}; end 0046 if nargin > 1, VM = varargin{2}; end 0047 if nargin > 2 0048 nic_spm('alert!', 'spm_est_smoothness: Wrong number of arguments'); 0049 return; 0050 end 0051 if nargin < 1 0052 V = nic_spm_select(inf, '^ResI.*\.img$', 'Select residual images'); 0053 end 0054 if nargin < 2 0055 VM = nic_spm_select(1, 'mask.img', 'Select mask image'); 0056 end 0057 0058 % intialise 0059 %----------------------------------------------------------------------- 0060 if ~isstruct(VM) 0061 V = nic_spm_vol(V); 0062 end 0063 if ~isstruct(VM) 0064 VM = nic_spm_vol(VM); 0065 end 0066 0067 %-Intialise RESELS per voxel image 0068 %----------------------------------------------------------------------- 0069 VRpv = struct('fname','RPV.img',... 0070 'dim', VM.dim(1:3),... 0071 'dt', [nic_spm_type('float64') nic_spm_platform('bigend')],... 0072 'mat', VM.mat,... 0073 'pinfo', [1 0 0]',... 0074 'descrip', 'spm_spm: resels per voxel'); 0075 VRpv = nic_spm_create_vol(VRpv); 0076 0077 0078 % dimensionality of image 0079 %----------------------------------------------------------------------- 0080 N = 3 - sum(VM.dim(1:3) == 1); 0081 if N == 0 0082 fwhm = [Inf Inf Inf]; 0083 return 0084 end 0085 0086 % find voxels within mask 0087 %----------------------------------------------------------------------- 0088 [x,y] = ndgrid(1:VM.dim(1), 1:VM.dim(2)); 0089 I = []; Ix = []; Iy = []; Iz = []; 0090 for i = 1:VM.dim(3) 0091 z = i*ones(size(x)); 0092 d = nic_spm_sample_vol(VM, x, y, z, 0); 0093 I = find(d); 0094 Ix = [Ix; x(I)]; 0095 Iy = [Iy; y(I)]; 0096 Iz = [Iz; z(I)]; 0097 end 0098 0099 % compute varianic_e of normalized derivatives in all directions 0100 %----------------------------------------------------------------------- 0101 str = 'Spatial non-sphericity (over scans)'; 0102 fprintf('%-40s: %30s',str,'...estimating derivatives') %-# 0103 nic_spm_progress_bar('Init',100,'smoothness estimation',''); 0104 0105 v = zeros(size(Ix,1),N); 0106 ssq = zeros(size(Ix,1),1); 0107 for i = 1:length(V) % for all residual images 0108 0109 [d, dx, dy, dz] = nic_spm_sample_vol(V(i), Ix, Iy, Iz, 1); 0110 0111 if N >= 1. v(:, 1) = v(:, 1) + dx.^2; end 0112 if N >= 2. v(:, 2) = v(:, 2) + dy.^2; end 0113 if N >= 3, v(:, 3) = v(:, 3) + dz.^2; end 0114 0115 ssq = ssq + d.^2; 0116 0117 nic_spm_progress_bar('Set',100*i/length(V)); 0118 0119 end 0120 nic_spm_progress_bar('Clear') 0121 0122 % normalise derivatives 0123 %----------------------------------------------------------------------- 0124 for i = 1:N 0125 v(:,i) = v(:,i)./ssq; 0126 end 0127 0128 % eliminate zero varianic_e voxels 0129 %----------------------------------------------------------------------- 0130 I = find(any(isnan(v'))); 0131 v(I,:) = []; Ix(I) = []; Iy(I) = []; Iz(I) = []; 0132 0133 0134 % resels per voxel (resel) 0135 % resels = resels/voxel = 1/prod(FWHM) 0136 % FWHM = sqrt(4.ln2/|dv/dx|)) 0137 % fwhm = 1/FWHM 0138 %----------------------------------------------------------------------- 0139 fprintf('\r%-40s: %30s\n',str,'...writing resels/voxel image') %-# 0140 0141 fwhm = sqrt(v./(4*log(2))); 0142 resel = prod(fwhm,2); 0143 for i = 1:VM.dim(3) 0144 d = NaN*ones(VM.dim(1:2)); 0145 I = find(Iz == i); 0146 if ~isempty(I) 0147 d(sub2ind(VM.dim(1:2), Ix(I), Iy(I))) = resel(I); 0148 end 0149 VRpv = nic_spm_write_plane(VRpv, d, i); 0150 end 0151 0152 % global equivalent FWHM {prod(1/FWHM) = (unbiased) RESEL estimator} 0153 %----------------------------------------------------------------------- 0154 fwhm = mean(fwhm ); 0155 RESEL = mean(resel); 0156 fwhm = fwhm*((RESEL/prod(fwhm)).^(1/N)); 0157 FWHM = 1./fwhm; 0158 fwhm = [Inf Inf Inf]; 0159 fwhm(1:length(FWHM)) = FWHM; 0160 0161 0162 0163