Home > rest_20090422 > rest_spm5_files > nic_spm_est_smoothness.m

nic_spm_est_smoothness

PURPOSE ^

Estimation of smoothness based on [residual] images

SYNOPSIS ^

function [fwhm,VRpv] = nic_spm_est_smoothness(varargin)

DESCRIPTION ^

 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_e

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

nic_spm_est_smoothness.m

SOURCE CODE ^

0001 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

Generated on Wed 29-Apr-2009 01:06:38 by m2html © 2005