Hi Jochen,
I can't explain your findings, but I know that if the following is
pasted into MATLAB, the kernel from spm_smoothkrn gives better
results.
The way that the design matrix is computed for fMRI is based on
similar principles. This is why it includes the concept of
"microtime" in order to deal more precisely with event onsets and
offsets.
All the best,
-John
fwhm = 2;
h = 1; % Results not so good if h=0
% because of the way half-way points are
% handled by SPMs NN interpolation
x = (1:50)';
x1 = (1:0.1:50)';
f = randn(size(x));
o = ones(size(x1));
f1 = spm_sample_vol(f,x1,o,o,h);
subplot(2,2,1); plot(x,f,'o',x1,f1,'.'); title('Original data');
gk = spm_smoothkern(fwhm ,[-6 :6 ]',h);
gk1 = spm_smoothkern(fwhm*10,[-60:60]',h);
s = ( fwhm/sqrt(8*log(2)))^2; g = 1/sqrt(2*pi*s)*exp(-(-6 :6 )'.^2/(2*s));
s = (10*fwhm/sqrt(8*log(2)))^2; g1 = 1/sqrt(2*pi*s)*exp(-(-60:60)'.^2/(2*s));
subplot(2,2,3); plot((-6:6)',[g gk],'o',(-60:60)'/10,[g1 gk1]*10,'.');
title('Smoothing Kernels');
sf = convn(f ,gk ,'same');
sf1 = convn(f1,gk1,'same');
subplot(2,2,2); plot(x,sf,'o',x1,sf1,'.'); title('spm smoothkern');
sf = convn(f,g,'same');
sf1 = convn(f1,g1,'same');
subplot(2,2,4); plot(x,sf,'o',x1,sf1,'.'); title('Simple Gaussian');
On 2 November 2011 19:00, Jochen Weber <[log in to unmask]> wrote:
> Hey John,
>
> Thanks a lot for your quick reply!!
>
> I've started to read on this, but when testing this the following way:
>
> - creating one single image (3mm resolution, 53x65x47 array size) of randn
> data and then
> - applying smoothing (6mm kernel) to the actual image or an image first
> interpolated to 1.5mm
> - and then comparing the voxel values at the same positions (all voxels of
> image 1 and odd index voxels of image 2)
>
> The differences of the values coming out of SPM's smoothing are, by average,
> larger than the differences using a "straight, discrete Gaussian" kernel
> (although, admittedly, spread over a smaller value range).
>
> I'm not entirely sure which way is the better way to smooth discrete data...
> If the smoothness estimation would work properly, I could even use that to
> try and find an answer...
>
> /jochen
>
>
> On Wed, Nov 2, 2011 at 1:29 PM, John Ashburner <[log in to unmask]> wrote:
>>
>> Hi Jochen,
>> I can't answer the second question, but I hope the first one is
>> answered OK. Note also that the integral under a Gaussian should be
>> 1, but if it is only sampled every voxel, then this is not the case.
>>
>> x=[-6:6];
>> sig2 = 0.1;
>> y=1/sqrt(2*pi*sig2)*exp(-0.5/sig2*x.^2);
>> sum(y)
>>
>> The problem with using a straight Gaussian is that stuff gets aliased out.
>>
>> All the best,
>> -John
>>
>>
>> On 2 November 2011 17:12, Jochen Weber <[log in to unmask]>
>> wrote:
>> > Hey John,
>> >
>> > I just posted two emails to the SPM mailing list and would appreciate if
>> > I
>> > got a reply from someone from the FIL. In case you don't mind my asking,
>> > could you please forward this to the person (or persons) responsible for
>> > spm_smoothkern.m and spm_est_smoothness.m ? Thanks a lot!! :)
>> >
>> > Problem 1:
>> >
>> > I was wondering whether this discrepancy has been noted before and
>> > discussed
>> > elsewhere (a brief search didn't yield any more results at least).
>> >
>> > After noticing the difference in smoothed images, I'd like to follow up
>> > on
>> > Robert's email with the following observation:
>> >
>> > When I use a call to spm_smoothkern (incl. the last SPM8 update) to
>> > print
>> > the kernel weights to the command line, like this
>> >
>> > spm_smoothkern(2, -6:6, 1)'
>> >
>> > I get the following weights:
>> >
>> > 0.0000
>> > 0.0000
>> > 0.0000
>> > 0.0026
>> > 0.0447
>> > 0.2417
>> > 0.4221
>> > 0.2417
>> > 0.0447
>> > 0.0026
>> > 0.0000
>> > 0.0000
>> > 0.0000
>> >
>> > Now, my understanding of the FWHM concept is that at full width (that is
>> > to
>> > say the distance between the two points with equal distance from the max
>> > and
>> > equal value), the value should be half the maximum. As is clear, 0.2417
>> > is
>> > (significantly) larger than half of 0.4221 (= 0.2111). According to what
>> > I
>> > can tell, this is more like a kernel of FWHM 2.2.
>> >
>> > In essence, I believe that the current implementation "over-smoothes"
>> > images. To further see to what extent this occurs, I computed the ratios
>> > between the values at FWHM point to maximum for different kernels using
>> > the
>> > following code:
>> >
>> > kernels = 1:0.2:4;
>> > kernel_ratios = zeros(size(kernels));
>> > for kc = 1:numel(kernels)
>> > k = spm_smoothkern(kernels(kc), [-kernels(kc)/2, 0], 1);
>> > kernel_ratios(kc) = k(1) / k(2);
>> > end
>> >
>> > When plotting the results with
>> >
>> > plot(kernels, kernel_ratios)
>> >
>> > I can observe that the smaller the desired kernel's FWHM, the larger the
>> > discrepancy.
>> >
>> > Now the question is: Is this *desired* or is this some kind of flaw in
>> > the
>> > new code of SPM8??
>> >
>> > Looking for answers...
>> >
>> >
>> >
>> > Problem 2:
>> >
>> > 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 Weber
>> > Research Assistant
>> > Social Cognitive Affective Neuroscience Unit
>> > In the Department of Psychology
>> > Columbia University
>> > In the City of New York
>> >
>> > 1190 Amsterdam Avenue
>> > 369 Schermerhorn Hall
>> > New York City, NY-10027
>> >
>> > em: [log in to unmask]
>> > ph: +1-212-854-6962
>> > mo: +1-773-234-8079
>> >
>> >
>> >
>
>
>
> --
> Jochen Weber
> Research Assistant
> Social Cognitive Affective Neuroscience Unit
> In the Department of Psychology
> Columbia University
> In the City of New York
>
> 1190 Amsterdam Avenue
> 369 Schermerhorn Hall
> New York City, NY-10027
>
> em: [log in to unmask]
> ph: +1-212-854-6962
> mo: +1-773-234-8079
>
>
>
|
