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
>
>
--
Guillaume Flandin, PhD
Wellcome Centre for Human Neuroimaging
UCL Queen Square Institute of Neurology
London WC1N 3BG
|