| >I suspect that what you really want is some kind of SPM that localises
| >significant shape differences (as done by Gaser et al).
|
| I once tried a Gaser's-like method as follows:
| Creating a symmetric T1 template by averaging L-R flipped and unflipped
| T1 template images.
|
| Creating a symmetric brainmask using the same way.
|
| Smoothing original T1 images with 8mm FWHM (24 subjects).
|
| Normalizing the smoothed images to the symmetric template using
| affine transformation only, with the symmetric brain mask.
You shouldn't need to smooth first. This is done as the first processing
step of spatial normalisation.
|
| Normalizing the affine-normalized images with 11*13*10 nonlinear
| basis functions using the same template and brainmask.
|
| Creating deformation field images (*_y1-3.img) from the *_sn3d.mat
| files produced by the second normalization, using spm_write_defs.m.
|
| Then, I created L-R flipped images of these deformation fields, and
| intended to perform 3-variate T-square test or MANOVA on the flipped
| and unflipped images. However, I encountered a confusing situation.
|
| All the deformation field images have a linear-like slope in the
| corresponding direction. For example, the values in the *_y1.img were
| low in the left hemisphere and high in the right hemisphere. They seemed
| to increase linearly (not exactly, just approximately) from left to right.
| Since this tendency was observed in all subjects, any statistical test on
| the hemispheric differences (flipped vs. unflipped) produced awfully
| low P values over the whole brain.
|
| Do you have any idea what caused this?
| I am afraid I am missing something very important.
The deformation fields contain information about position as well as shape
information. Therefore, it is important that the positional information
is first removed. You may also wish to factor out the effects of size.
This can be done by finding the rotations, translations and possibly the
isotropic zoom that minimises the procrustes distance between the images,
and then transforming the deformation fields so that this information
is removed.
The attached code may provide some clues about how this can be done, but it
will probably not work "as is". Your deformation fields will need to be
either converted so that they map to co-ordinates in millimeters, or you
need to modify the code so that it converts from voxels to millimeters.
I am not that keen on methods that produce SPMs based on absolute
displacements. I would much rather see analyses done on volume
differences, or other analyses done on some features of the Jacobian
matrices of the deformations. This is basically because the results
are so dependent on the way that the positional information is removed
from the deformations.
Best regards,
-John
function spm_procrustes_ui
% Removes pose (and global size effects) from deformation fields.
% Corrected fields are written prefixed by 'p'.
%_______________________________________________________________________
% Example Reference:
% F. L. Bookstein (1997). "Landmark Methods for Forms Without
% Landmarks: Morphometrics of Group Differences in Outline Shape"
% Medical Image Analysis 1(3):225-243
%_______________________________________________________________________
% %W% John Ashburner %E%
n = spm_input('Number of subjects','+0', 'n', '1', 1)';
for i=1:n,
P{i} = spm_get(3,'*y?_*.img',['Select deformation field ' num2str(i)]);
end;
PW = spm_get(1,'*.img','Weighting image');
flg = spm_input('Remove size?','+1','y/n',[1,0],1);
VW = spm_vol(PW);
spm_progress_bar('Init',n,'Tweaking deformations','volumes completed');
for i=1:n,
doit(P{i},VW,flg);
spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear')
return;
%_______________________________________________________________________
%_______________________________________________________________________
function doit(P,VW,flg)
V = spm_vol(P);
y1 = spm_load_float(V(1));
y2 = spm_load_float(V(2));
y3 = spm_load_float(V(3));
T = procrustes(y1,y2,y3,V(1).mat,VW,flg);
spm_affdef(y1,y2,y3,T);
disp(T)
VO = V(1);
VO.fname = prepend(V(1).fname, 'p');
VO.dim(4) = spm_type('float');
VO.pinfo = [1 0 0]';
VO.descrip = 'Shape field - X';
spm_write_vol(VO,y1);
VO = V(2);
VO.fname = prepend(V(2).fname, 'p');
VO.dim(4) = spm_type('float');
VO.pinfo = [1 0 0]';
VO.descrip = 'Shape field - Y';
spm_write_vol(VO,y2);
VO = V(3);
VO.fname = prepend(V(3).fname, 'p');
VO.dim(4) = spm_type('float');
VO.pinfo = [1 0 0]';
VO.descrip = 'Shape field - Y';
spm_write_vol(VO,y3);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function out = prepend(in, pre)
[pth,nme,ext,ver] = fileparts(in);
out = fullfile(pth,[pre nme ext ver]);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function T = procrustes(Y1,Y2,Y3,M1,V,flg)
% Compute the affine transformation that minimises the
% Procrustes distance between transformed and un-transformed
% points.
%
% FORMAT T = procrustes(Y1,Y2,Y3,M1,V,flg)
% Y1, Y2 and Y3 - spatial mapping
% M1 - voxel to mm mapping of Y1, Y2 and Y3.
% V - a handle for a weighting image.
% flg - a flag to say if zooming should be included.
% 1 means no
% T - the determined affine transformation matrix.
%_______________________________________________________________________
% Example Reference:
% F. L. Bookstein (1997). "Landmark Methods for Forms Without
% Landmarks: Morphometrics of Group Differences in Outline Shape"
% Medical Image Analysis 1(3):225-243
%_______________________________________________________________________
% %W% John Ashburner %E%
[z1,z2] = ndgrid(1:size(Y1,1),1:size(Y1,2));
z1 = z1(:);
z2 = z2(:);
%_______________________________________________________________________
% Compute centres of mass.
sm = 0;
c1 = [0 0 0];
c2 = [0 0 0];
for j=1:size(Y1,3),
M = V.mat\M1*spm_matrix([0 0 j]);
wt = spm_slice_vol(V,M,V.dim(1:2),0); wt = wt(:);
y1 = double(Y1(:,:,j)); y1 = y1(:);
y2 = double(Y2(:,:,j)); y2 = y2(:);
y3 = double(Y3(:,:,j)); y3 = y3(:);
x1 = M1(1,1)*z1 + M1(1,2)*z2 + (M1(1,3)*j + M1(1,4));
x2 = M1(2,1)*z1 + M1(2,2)*z2 + (M1(2,3)*j + M1(2,4));
x3 = M1(3,1)*z1 + M1(3,2)*z2 + (M1(3,3)*j + M1(3,4));
sm = sm + sum(wt);
c1 = c1 + sum([x1.*wt x2.*wt x3.*wt]);
c2 = c2 + sum([y1.*wt y2.*wt y3.*wt]);
end;
c1 = c1/sm;
c2 = c2/sm;
%_______________________________________________________________________
if nargin<6 | flg,
% Compute centroid size.
sm = 0;
mom1 = 0;
mom2 = 0;
for j=1:size(Y1,3),
M = V.mat\M1*spm_matrix([0 0 j]);
wt = spm_slice_vol(V,M,V.dim(1:2),0);
wt = wt(:);
y1 = double(Y1(:,:,j)); y1 = y1(:);
y2 = double(Y2(:,:,j)); y2 = y2(:);
y3 = double(Y3(:,:,j)); y3 = y3(:);
x1 = M1(1,1)*z1 + M1(1,2)*z2 + (M1(1,3)*j + M1(1,4));
x2 = M1(2,1)*z1 + M1(2,2)*z2 + (M1(2,3)*j + M1(2,4));
x3 = M1(3,1)*z1 + M1(3,2)*z2 + (M1(3,3)*j + M1(3,4));
sm = sm + sum(wt);
mom1 = mom1 + sum([(x1-c1(1)).^2.*wt (x2-c1(2)).^2.*wt (x3-c1(3)).^2.*wt]);
mom2 = mom2 + sum([(y1-c2(1)).^2.*wt (y2-c2(2)).^2.*wt (y3-c2(3)).^2.*wt]);
end;
zm = sqrt(mom1/mom2);
else,
zm = 1;
end;
ZM = spm_matrix([0 0 0 0 0 0 zm zm zm 0 0 0]);
%_______________________________________________________________________
% Compute rotations
C = zeros(3);
for j=1:size(Y1,3),
M = V.mat\M1*spm_matrix([0 0 j]);
wt = spm_slice_vol(V,M,V.dim(1:2),0);
wt = wt(:);
y1 = double(Y1(:,:,j)); y1 = y1(:);
y2 = double(Y2(:,:,j)); y2 = y2(:);
y3 = double(Y3(:,:,j)); y3 = y3(:);
x1 = M1(1,1)*z1 + M1(1,2)*z2 + (M1(1,3)*j + M1(1,4));
x2 = M1(2,1)*z1 + M1(2,2)*z2 + (M1(2,3)*j + M1(2,4));
x3 = M1(3,1)*z1 + M1(3,2)*z2 + (M1(3,3)*j + M1(3,4));
wt1 = wt*zm;
C = C + [(x1-c1(1)).*wt (x2-c1(2)).*wt (x3-c1(3)).*wt ]' * ...
[(y1-c2(1)).*wt1 (y2-c2(2)).*wt1 (y3-c2(3)).*wt1];
end;
[u,s,v] = svd(C);
R = eye(4);
R(1:3,1:3) = u*v';
%_______________________________________________________________________
% Combine rotations, translations (and zooms) into a single
% affine transformation matrix.
T = spm_matrix(c1)*R*ZM*spm_matrix(-c2);
return;
%_______________________________________________________________________
function dat = spm_load_float(V)
% Load a volume into a floating point array
% FORMAT dat = spm_load_float(V)
% V - handle from spm_vol
% dat - a 3D floating point array
%_______________________________________________________________________
% %W% John Ashburner %E%
dim = V(1).dim(1:3);
dat = single(0);
dat(dim(1),dim(2),dim(3))=0;
for i=1:V(1).dim(3),
M = spm_matrix([0 0 i]);
dat(:,:,i) = single(spm_slice_vol(V(1),M,dim(1:2),1));
end;
return;
|