Dear list (and SPM coders?),
When trying to understand the logic of the code for the SPM residuals smoothness estimation, I ran into a slight snag (part of which is, I believe, due to the over-smoothing of images by spm_smooth due to wrong kernel weights, see my other post of today).
I created a 4D series of images with the following properties:
v = spm_vol('randn_64vols.nii');
v(1)
fname: 'randn_64vols.nii'
dim: [53 65 47]
mat: [4x4 double]
pinfo: [3x1 double]
dt: [16 0]
n: [1 1]
descrip: 'spm - 3D normalized'
private: [1x1 nifti]
v(1).mat
ans =
3 0 0 -81
0 3 0 -117
0 0 3 -57
0 0 0 1
v(1).pinfo
ans =
1
0
352
Then I set the data to
single(randn(53, 65, 47, 64));
and computed the estimated FWHM (and RPV image) thus (using a mask image that had all voxels included):
% creating a list of filenames for SPM
fs = cell(64, 1);
for c=1:64, fs{c} = sprintf('randn_64vols.nii,%d', c); end
% estimating the FWHM and RPV vol
[fwhm, rpv] = spm_est_smoothness(spm_vol(char(fs)), spm_vol('mask.nii'), [64, 62]);
The resulting FWHM was
1.3418 1.3369 1.3453
which, in my opinion, is already to high (it should be closer to 1), and the sum of all RPV voxels over the number of voxels was
rpvvox = spm_read_vols(rpv);
sum(rpvvox(:)) / numel(rpvvox)
0.4144
Thus estimating that there are only 0.41 resels per voxel (although randn hopefully produces entirely independent time series).
Using this "baseline bias", I then applied different amounts of smoothing, starting at 3mm and increasing up to 10mm, followed by subsequent estimation of the smoothness of the "residual" data using the following code:
% create new filenames list
sfs = strrep(fs, 'rand', 'srand');
% iterate over kernels
for k = 3:10
% remove old RPV image (just in case)
!rm RPV*
% smooth data
for c=1:64, spm_smooth(fs{c}, sfs{c}, [k,k,k]); end
% estimate smoothness
[fwhm, rpv] = spm_est_smoothness(spm_vol(char(sfs)), spm_vol('mask.nii'), [64, 62]);
% display nominal and estimated fwhm
disp([k, fwhm]);
% display nominal fwhm and average RPV
rpvv = spm_read_vols(rpv);
disp([k, sum(rpvv(:)) / numel(rpvv)]);
end
I obtained the following results:
Spatial non-sphericity (over scans) : ...done
3.0000 4.7406 4.7227 4.7715
3.0000 0.0094
Spatial non-sphericity (over scans) : ...done
4.0000 7.2801 7.2525 7.3370
4.0000 0.0026
Spatial non-sphericity (over scans) : ...done
5.0000 10.6155 10.5741 10.7112
5.0000 0.0008
Spatial non-sphericity (over scans) : ...done
6.0000 14.8334 14.7717 14.9839
6.0000 0.0003
Spatial non-sphericity (over scans) : ...done
7.0000 20.0935 20.0021 20.3191
7.0000 0.0001
Spatial non-sphericity (over scans) : ...done
8.0000 26.5266 26.3926 26.8528
8.0000 0.0001
Spatial non-sphericity (over scans) : ...done
9.0000 34.2301 34.0364 34.6877
9.0000 0.0000
Spatial non-sphericity (over scans) : ...done
10.0000 43.2869 43.0122 43.9119
10.0000 0.0000
In fact, using a 10mm kernel, the sum of ALL voxels of the RPV image was 1.98 (estimating a total of 2 RESELs in the entire volume).
Can anybody please explain how this makes any sense?? Thanks!
/jochen
|