When you smooth, you essentially compute local weighted averages of
the voxel intensities. For simplicity, lets assume that we are just
computing local averages by summing up the values in a region of
voxels and dividing by the number of voxels.
The smoothing that is part of the normalisation to MNI space computes
these average intensities from the original data, rather than the
warped versions. When the data are warped, some voxels will grow and
others will shrink. This will change the regional averages, with more
weighting towards those voxels that have grows.
To compute a regional average from the original data, the algorithm
sums up the values within the region and divides by the number of
voxels, but does it in a way that properly considers the relative
volumes of warped and unwarped voxels. This gives a different answer
to simply smoothing the warped data.
For those who are interested, I've re-pasted something a bit more
mathematical below.
Best regards,
-John
To get a bit more mathsy, here's a bit of MATLAB code to illustrate
what's going on in 1D. 3D is a lot more complicated.
% Generate a random deformation (y) and its inverse (iy)
N=256;
h=[1 1 1; 1 1 1];
L = toeplitz([2.001 -1 zeros(1,N-3) -1]*4); L(1,1)=1e6; L(N,N)=1e6;
u = L\randn(N,1);
o=ones(N,1);
y=(1:N)'+u/2^8; for k=1:8, y = spm_bsplins(y,y,o,o,h); end; %
Scaling and squaring
iy = (1:N)'-u/2^8; for k=1:8, iy = spm_bsplins(iy,iy,o,o,h); end %
Scaling and squaring
% Generate a smoothing kernel. Note that each value in a smoothed image is
% a dot-product between the image and a smoothing kernel. The same rules apply
% when summing the values in a region of interest.
s= spm_smoothkern(40,(1:N)'-N/2,1);
% Generate a random image
f = toeplitz(spm_smoothkern(5,[0:floor((N-1)/2) floor(-(N-1)/2):-1]
,1))*rand(N,1);
% The deformation can be represented as a matrix multiplication (by G). Here we
% compute G (and its approximate inverse iG).
G=zeros(N);
iG=zeros(N);
tmp = [1 zeros(1,N-1)];
for i=1:N,
G(i,:) = spm_bsplins(tmp,o',(y(i)+(N:-1:1)),o',h);
iG(i,:) = spm_bsplins(tmp,o',(iy(i)+(N:-1:1)),o',h);
end
imagesc(G); hold on; plot(y,(1:N)','r',(1:N)',iy,'g'); hold off; axis xy
% Just to check, plot s(y), as well as G*s
sy=spm_bsplins(s,y,o,o,h);
plot([G*s sy]);
% Note that various approximations mean that iG can not be the exact
inverse of G (crude interpolation, discrete approximation of a
continuous thing etc)
imagesc([G iG; iG*G G*iG])
% What we want in our preprocessed images is a dot-product of the
smoothing kernel
% (or some other definition of a region) with the original data. So,
we can warp the
% smoothing kernel (or ROI) back to the native image and compute the
dot-product in
% native space:
dot(G*s,f)
% or we can use the adjoint (transpose) of G to push the values in f into
% normalised space and then do our dot-product
dot(G'*f,s)
% These are all equivalent to:
f'*G*s
% Or we can use the inverse of the transform (ie iy) to warp our image into
% normalised space with a change of variables and then do the dot-product
dot((gradient(iy).*spm_bsplins(f,iy,o,o,h)),s)
dot(diag(gradient(iy))*iG*f,s)
% Note that iG'*G' should be roughly identity, so f*G*s is roughly the same as
% f'*iG'*G'*G*s. This gives the following dot-product:
dot((G'*G)*(iG*f),s)
% G'*G is roughly diagonal, and very approximately the same as
diag(gradient(iy))
imagesc([diag(gradient(iy)) G'*G])
On 31 October 2011 23:03, YAN Chao-Gan <[log in to unmask]> wrote:
> Dear John and DARTEL experts,
>
> Thanks a lot for the excellent tool of DARTEL!
>
> Recently, I am using DARTEL to normalize my functional data into MNI space.
> I found the following two ways would give different results.
>
> 1. DARTEL (Normalize to MNI space): Gaussian FWHM - 8 8 8. I will get the
> swra* images.
>
> 2. DARTEL (Normalize to MNI space): Gaussian FWHM - 0 0 0. I will get the
> wra* images. Then I perform smoothing by the Smooth module in SPM (kernel: 8
> 8 8), and I will get the swra* images.
>
> However, The swra* images of these two approaches are different. From the
> manual: "Note that you can also specify [0 0 0], but any “modulated’ data
> will show aliasing (see eg Wikipedia), which occurs because of the way the
> warped images are generated. "
>
> My images is unmodulated in this case, right? I wonder why this difference
> occurs. Thanks a lot!
>
> Best,
>
> Chao-Gan
>
|