 |
 |
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
[log in to unmask]">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:
[log in to unmask]">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