PS:
I just had closer look at the code John kindly provided and compared one
dimension of krn calculated with the result from spm_smoothkern.
Is it correct to assume, that their differences are due to the
continuity issue John mentioned?
Additionally, I compared a convolution of spm_smoothkern(4, -39:39,1)
with d1 = squeeze(volume2smooth(40, 48, : )) ;
This nearly perfectly replicates the result from the SPM-gui-smoothed
volume.
Regards,
Peter
Am 05.02.2020 um 11:13 schrieb Peter Sperrer:
>
> Dear Guillaume, Dear John,
>
> thank you for your kind answers.
>
> I have to admit, until now my mind somehow clung to the idea of kernel
> and voxel values being not defined anywhere outside of the volume, but
> of course this is not an imperative and can be done differently.
> Especially the explanation considering the choice of boundary
> conditions was very helpful.
>
> * Considering the observed difference between Matlab's convolution
> and SPM's convolution:
> The test volume was defined as
>
> volume2smooth = zeros(79,95,79) ; volume2smooth(40,48,40) = 1;
>
> and saved with voxel dimensions of 2*2*2 mm as 'testImageDot.nii'.
> This volume was smoothed in SPM's gui using an isotropic FWHM=8mm
> kernel.
> Then, I tried replicating the smoothing process. Since the voxel
> of interest is far away from the borders, I assumed no significant
> boundary effects should be observed.
>
> As I use voxels for distances below, not mm, I divided the FWHM by 2:
>
> FWHM = 8;
> FWHM = FWHM/2;
>
> s = FWHM/(sqrt(8*log(2)));
> s6 = round(6*s);
>
> w1 = 0.5*sqrt(2/s);
> w2 = -0.5/s;
> w3 = sqrt(s/2/pi);
>
> filePath = '/XXX/testImageDot.nii';
> volume2smooth = niftiread(filePath);
>
> filePathSmoothed = '/XXX/stestImageDot.nii';
> smoothedVolumeSPM = niftiread(filePathSmoothed);
> smoothedVolumeSPM = 100*smoothedVolumeSPM/max(smoothedVolumeSPM,
> [], 'all');
>
> % dimensions of the volume
> xi = 79;
> xii = 95;
> xiii = 79;
>
> x = double(bwdist(logical(volume2smooth), 'euclidean'));
>
> % restrict kernel size to 6*s:
> xx = x(ceil(xi/2)-s6:ceil(xi/2)+s6, ceil(xii/2)-s6:ceil(xii/2)+s6,
> ceil(xiii/2)-s6:ceil(xiii/2)+s6);
>
> krn = 0.5*(erf(w1*(xx+1)).*(xx+1) + erf(w1*(xx-1)).*(xx-1) -
> 2*erf(w1*xx).* xx)...
> + w3*(exp(w2*(xx+1).^2) + exp(w2*(xx-1).^2) - 2*exp(w2*xx.^2));
> krn(krn<0) = 0;
>
> g = 1/(sqrt(2*pi)*s)^3*exp(-x.^2/(2*s^2));
>
> smoothedVolumeKrn = convn(volume2smooth, krn, 'same');
> smoothedVolumeKrn = 100*smoothedVolumeKrn/max(smoothedVolumeKrn,
> [], 'all');
>
> smoothedVolumeG = convn(volume2smooth, g, 'same');
> smoothedVolumeG = 100*smoothedVolumeG/max(smoothedVolumeG, [], 'all');
>
> figure
> hold on
> plot(squeeze(smoothedVolumeG(40, 48, :)));
> plot(squeeze(smoothedVolumeKrn(40, 48, :)));
> plot(squeeze(smoothedVolumeSPM(40, 48, :)));
> legend('Gaussian', 'Erf', 'SPM');
> hold off
>
> When comparing smoothedVolumeKrn and smoothedVolumeSPM I would
> expect them to be (nearly) identical, but smoothedVolumeSPM is
> much closer to SmoothedVolumeG.
> Visual comparison was performed by plotting a histogram of the
> normalized intensities (so all the maxima are the same) along the
> y-coordinate with x=40 and z=48.
> Probably the difference is caused by me missing additional
> normalization of the kernel?
>
> * Considering smoothing along different dimensions:
> Yes, you are right, when smoothing a volume with just ones, the
> visible difference along z compared to x or y is (mostly?)
> attributed to the boundary conditions.
> Still, when smoothing a single non-zero voxel in the center of a
> large volume, I would expect the voxels along x, y and z to be
> identical (as they are for Matlab's convolution). To test this I used
>
> a = squeeze(smoothedVolumeSPM(40, 48, 40:60));
> b = squeeze(smoothedVolumeSPM(40, 48:68, 40))';
>
> c = a - b;
>
> Admittedly, max(c) is in the order of 10^-14 and thus is probably
> unlikely to significantly influence the results; But since the
> volume used double precision I think this difference should not be
> be attributed to numerical errors?
>
> Again, thank you very much, I appreciate your help.
>
> Kind regards,
>
> Peter
>
>
>
>
> Am 04.02.2020 um 16:37 schrieb Flandin, Guillaume:
>> Dear Peter,
>>
>> Your description is accurate and you are looking at the right functions
>> (spm_smooth.m, spm_smoothkern.m, spm_conv_vol.m/c).
>> The (separable) convolution kernel is a Gaussian convolved with a 1st
>> degree B-spline (triangular) basis function, thresholded at 6 sigmas and
>> normalised to sum to 1.
>> For the convolution itself, the boundary conditions correspond to use
>> zero-padding in-plane (x and y) and truncation out-of-plane (z). See e.g.:
>> https://en.wikipedia.org/wiki/Kernel_(image_processing)#Convolution
>> https://en.wikipedia.org/wiki/Gaussian_blur#Implementation
>> Concerning your evaluation, it would be useful if you shared the code
>> you used. Could it be that you are using an image not large enough
>> compared to the size of the kernel and you observe the effect of
>> boundary conditions?
>>
>> Best regards,
>> Guillaume.
>>
>>
>> On 30/01/2020 12:05, Peter Sperrer wrote:
>>> Dear SPM-Experts,
>>>
>>> I am currently trying to understand what exactly SPM(12) does from a
>>> mathematical point of view during smoothing (for default values).
>>> Unfortunately, neither code nor documentation answered all my questions
>>> and I could not find related questions in the mailing list.
>>>
>>> Could you please explain to me, how smoothing is performed or refer me
>>> to a paper or manual with such an explanation?
>>> In detail, the question is how the kernel function is exactly defined
>>> and how the convolution is applied.
>>>
>>> Doing some tests (FWHM = 8 mm, voxel size = 2x2x2 mm) and digging into
>>> SPM's code and documentation I came to the following insights so far:
>>>
>>> * SPM seems to use something similar (but not identical) to a Gaussian
>>> kernel:
>>>
>>> spm_smooth.m calls spm_smoothkern.m
>>> This returns (as default) the following kernel krn, which according
>>> to the comments resembles a Gaussian convolved with a 1st degree
>>> B-spline:
>>>
>>> w1 = 0.5*sqrt(2/s);
>>> w2 = -0.5/s;
>>> w3 = sqrt(s/2/pi);
>>> krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) -
>>> 2*erf(w1*x ).* x)+w3*(exp(w2*(x+1).^2) + exp(w2*(x-1).^2)
>>> - 2*exp(w2*x.^2));
>>> krn(krn<0) = 0;
>>>
>>> with s being the standard deviation of a Gaussian with the given
>>> FWHM and x being in the rounded range of [-6 s, 6s].
>>>
>>> Then, spm_conv_vol.c is called, which according to spm_cov_vol.m
>>> should do some sort of convolution.
>>> Unfortunately, my knowledge of C is rather limited, so I am having
>>> trouble understanding what exactly is going on there.
>>>
>>> * Even with isotropic FWHM values, smoothing in z is not the same than
>>> for x and y.
>>>
>>> This can best be seen when smoothing a volume of constant voxel
>>> values, but is also mentioned in the comments of spm_conv_vol.m:
>>> "The convolution assumes zero padding in x and y with truncated
>>> smoothing in z."
>>> I assume, this also hints the use of something different than a
>>> simple 3D convolution?
>>>
>>> * For a single, non-zero voxel in the center of a volume, the
>>> (normalized) results do not match a convolution of the image with
>>> the above function (krn).
>>>
>>> It would be great, if any of you could help me understand what exactly
>>> happens during the smoothing procedure.
>>>
>>> Thank you very much for your time and your expertise!
>>>
>>> Kind regards,
>>>
>>> Peter Sperrer
>>>
>>>
