| 1. How can SPM best be used to compare a single subject to a group of
| controls in order to establish the pattern of regional abnormalities? I
| have tried using the two sample t-test, two groups, one scan per subject
| model, with success, but was wondering if anyone had ideas about other
| approaches using the software.
This is probably the best approach, but it may be worth also modelling
confounding effects such as age or possibly nonlinear age effects
(by also including age^2 and age^3).
Depending how many controls you have, you may also wish to try a
non-parametric analysis using SNPM.
|
| 2. Can the SVC analysis be applied to any arbitrary contour defined by the
| user (say, a hand-drawn outline of the frontal lobe), and how would one
| accomplish this?
I'll let someone else answer this one.
|
| 3. Can the SPM software generate a deformation map from the spatial
| normalization, so that one could see the patterns of correction in different
| regions?
There is nothing in the current release of SPM99 that will do this, but I have
appended a piece of Matlab code that should do what you want.
Best regards,
-John
function spm_Deformations
% Deformation field utilities (works on sn3d.mat files).
%
% This toolbox includes a number of utilities for extracting information
% from _sn3d.mat files created during spatial normalisation. Much of
% the information can be used for either deformation-based morphometry
% (DBM), or tensor-based morphometry (TBM).
%
% Deformations:
% Writes out the deformation fields as *_y1.img, *_y2.img and *_y3.img.
% These fields are the voxel to voxel mappings between an image
% normalised with the specified bounding-box and voxel sizes, and the
% original images. The deformations can be used for multivariate
% morphometric methods (DBM), but they will need some initial
% corrections for voxel-sizes and for some measure of global position
% and possibly size (Procrustes shape).
%
% Jacobian Matrices:
% Writes out the Jacobian matrix field as images a11_*.img, a12_*.img,
% a13_*.img, a21_*.img, a22_*.img, a23_*.img, a31_*.img, a32_*.img
% and a33_*.img.
%
% Jacobian Determinants:
% Writes out the volume change at each location in j_*.img. Fields
% involving a flip should have begative Jacobian determinants, but
% these are implicitly made positive to make things easier to
% understand. The determinants can be used for morphometric methods
% (TBM) that characterise volumetric differences.
%
% Strain Tensors:
% These are intended for voxelwise multivariate morphometry (TBM).
% Much of the notation and ideas are from:
% "Non-linear Elastic Deformations", by R. W. Ogden. Dover
% Publications, 1984.
%
% The transformation maps from elements in the template (x1,x2,x3)
% to elements in the original images (y1,y2,y3).
% An affine mapping from x to y is given by:
% y = A*x + c, where A is the deformation gradient (second order tensor
% which is the Jacobian matrix) and c is a constant representing
% translations.
% We wish to represent shape changes within the co-ordinate framework of
% the template images (Lagrangean framework, as opposed to the Eulerian
% framework where deformations are relative to the individual images).
% To do this, matrix A is decomposed (using polar decomposition), such
% that A = R*U. Matrix R is a rigid body transformation matrix, and U
% is a matrix that is purely shears and zooms (no rotations).
% This gives a mapping from template voxels to image voxels of the form
% y = R*U*x + c, which means zoom and shear the positions in the template,
% and then do a rigid body rotation (and also add the translation).
% We are not interested in the translations and rotations - only in the
% local zooms and shears. Therefore all the information we need is in
% the matrix U, which is simply obtained by U=(A'*A)^(1/2).
% Note that the zooms and shears are done while in the orientation of
% the template, before rotations to the orientation of the images are
% introduced.
%
% Strain tensors are defined that model the amount of distortion. If
% there is no strain, then the tensors are all zero. Generically,
% the family of Lagrangean strain tensors are given by:
% (U^m-eye(3))/m when m~=0, and logm(U) when m==0.
%
%
%_______________________________________________________________________
% %W% John Ashburner %E%
global sptl_BB sptl_Vx BCH
SPMid = spm('FnBanner',mfilename,'%A%');
[Finter,Fgraph,CmdLine] = spm('FnUIsetup','Deformations');
spm_help('!ContextHelp',mfilename);
a1 = spm_input('Write what?',1,'m',...
['Deformations|Jacobian Matrices|Jacobian Determinants|' ...
'Strain Tensors, m=-2 (Almansi)|Strain Tensors, m=-1 |Strain Tensors, m= 0 (Hencky)|'...
'Strain Tensors, m= 1 (Biot) |Strain Tensors, m= 2 (Green)'],...
[100 101 102 -2 -1 0 1 2],3,'batch',{},'option');
if isempty(BCH),
files = spm_get(Inf,'*_sn3d.mat','Select *_sn3d.mat files');
else,
files = spm_input('batch', {},'sn3d_files');
end;
bboxes = [ -78 78 -112 76 -50 85
-64 64 -104 68 -28 72
-90 91 -126 91 -72 109
-95 95 -112 76 -50 95];
bbprompt = [' -78:78 -112:76 -50:85 (Default)|'...
' -64:64 -104:68 -28:72 (SPM95) |'...
' -90:91 -126:91 -72:109 (Template)|'...
' -95:95 -112:76 -50:95 '];
voxdims = [ 1 1 1 ; 2 2 2 ; 3 3 3 ; 4 4 4 ; 5 5 5 ; 6 6 6 ; 7 7 7 ; 8 8 8 ; 9 9 9 ; 10 10 10];
voxprompts = ' 1 1 1 | 2 2 2 | 3 3 3 | 4 4 4 | 5 5 5 | 6 6 6 | 7 7 7 | 8 8 8 | 9 9 9 | 10 10 10';
bb = sptl_BB;
if prod(size(bb)) == 6,
tmp = find( sptl_BB(1) == bboxes(:,1) & sptl_BB(2) == bboxes(:,2) & ...
sptl_BB(3) == bboxes(:,3) & sptl_BB(4) == bboxes(:,4) & ...
sptl_BB(5) == bboxes(:,5) & sptl_BB(6) == bboxes(:,6));
if isempty(tmp), tmp = size(bboxes,1)+1; end;
else,
tmp = size(bboxes,1)+2;
bb = reshape(bboxes(1,:),2,3);
end;
ans = spm_input('Bounding Box?','+1','m',...
[ bbprompt '|Customise'], [1:size(bboxes,1) 0],...
tmp,'batch',{},'bounding_box');
if ans>0, bb=reshape(bboxes(ans,:),2,3);
else,
if prod(size(bb)) ~= 6, bb = reshape(bboxes(1,:),2,3); end;
directions = 'XYZ';
for d=1:3,
str = sprintf('%d %d', bb(1,d), bb(2,d));
bb(:,d) = spm_input(['Bounding Box ' directions(d) ],....
'+1', 'e',str, 2,'batch',{},...
sprintf('direction%d',d));
end;
end;
Vox = sptl_Vx;
if prod(size(Vox)) == 3,
tmp = find(voxdims(:,1) == Vox(1) & voxdims(:,2) == Vox(2) & voxdims(:,3) == Vox(3));
if isempty(tmp), tmp = size(voxdims,1)+1; end;
else, tmp = size(voxdims,1)+2; end;
ans = spm_input(...
['Voxel Sizes?'], '+1','m', [ voxprompts ....
'|Customise'], [1:size(voxdims,1) 0],....
tmp,'batch',{},'voxel_sizes');
if ans>0, Vox = voxdims(ans,:);
else,
vxdef = [8 8 8];
if (prod(size(Vox)) ~= 3) vxdef = [5 5 5]; end
Vox = spm_input('Voxel Sizes ','+0', 'e', ...
sprintf('%d %d %d', vxdef(1), vxdef(2), vxdef(3)), 3,...
'batch',{},'voxel_sizes_custom')';
Vox = reshape(Vox,1,3);
end;
for i=1:size(files,1),
sn3d(i) = load(deblank(files(i,:)));
end;
for i=1:size(files,1),
sn3d(i).fname = deblank(files(i,:));
end;
if a1<100,
spm_progress_bar('Init',length(sn3d),['Computing Tensors (m= ' num2str(a1) ')'],'volumes completed');
elseif a1==100,
spm_progress_bar('Init',length(sn3d),['Computing Deformations'],'volumes completed');
elseif a1==101,
spm_progress_bar('Init',length(sn3d),['Computing Jacobian Matrices'],'volumes completed');
elseif a1==102,
spm_progress_bar('Init',length(sn3d),['Computing Jacobian Determinants'],'volumes completed');
end;
spm('Pointer','Watch')
for i=1:length(sn3d),
spm('FigName',['Deformations: working on subj ' num2str(i)],Finter,CmdLine);drawnow;
if a1<100,
write_tensor(sn3d(i),Vox,bb,a1);
elseif a1==100,
spm_write_defs(sn3d(i),Vox,bb);
elseif a1==101,
write_jacobianm(sn3d(i),Vox,bb);
elseif a1==102,
write_det(sn3d(i),Vox,bb);
end;
spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');
spm('FigName','Deformations: done',Finter,CmdLine);
spm('Pointer');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function write_jacobianm(sn3d,vox,bb)
vo = init_vo(sn3d,vox,bb);
for i=1:3,
for j=1:3,
[pth,nm,xt,vr] = fileparts(deblank(sn3d.fname));
if length(nm)>6 & strcmp(nm(end-4:end),'_sn3d'), nm=nm(1:end-5); end;
VO(i,j) = vo;
VO(i,j).fname = fullfile(pth,['a' num2str(i) num2str(j) '_' nm '.img']);
VO(i,j).dim(4) = spm_type('float');
VO(i,j).pinfo = [1 0 0]';
VO(i,j).descrip = ['Jacobian_Matrix(' num2str(i) ',' num2str(j) ')'];
VO(i,j) = spm_create_image(VO(i,j));
end;
end;
for p=1:vo.dim(3),
A = get_jacobian(sn3d, vo, p);
for i=1:3,
for j=1:3,
VO(i,j) = spm_write_plane(VO(i,j),A(:,:,i,j),p);
end;
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function write_det(sn3d,vox,bb)
VO = init_vo(sn3d,vox,bb);
[pth,nm,xt,vr] = fileparts(deblank(sn3d.fname));
if length(nm)>6 & strcmp(nm(end-4:end),'_sn3d'), nm=nm(1:end-5); end;
VO.fname = fullfile(pth,['j_' nm '.img']);
VO.dim(4) = spm_type('float');
VO.pinfo = [1 0 0]';
VO.descrip = ['Jacobian_Determinant'];
VO = spm_create_image(VO);
tmp = sn3d.MF*sn3d.Affine*inv(sn3d.MG);
sgn = sign(det(tmp(1:3,1:3)));
for p=1:VO.dim(3),
A = get_jacobian(sn3d, VO(1,1), p);
dtA = A(:,:,1,1).*(A(:,:,2,2).*A(:,:,3,3) - A(:,:,2,3).*A(:,:,3,2)) ...
- A(:,:,2,1).*(A(:,:,1,2).*A(:,:,3,3) - A(:,:,1,3).*A(:,:,3,2)) ...
+ A(:,:,3,1).*(A(:,:,1,2).*A(:,:,2,3) - A(:,:,1,3).*A(:,:,2,2));
if sgn < 0, dtA = -dtA; end;
VO = spm_write_plane(VO,dtA,p);
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function write_tensor(sn3d,vox,bb,m)
vo = init_vo(sn3d,vox,bb);
ij = [1 1; 2 1; 3 1; 2 2; 3 2; 3 3];
for i=1:6,
[pth,nm,xt,vr] = fileparts(deblank(sn3d.fname));
if length(nm)>6 & strcmp(nm(end-4:end),'_sn3d'), nm=nm(1:end-5); end;
VO(i) = vo;
VO(i).fname = fullfile(pth,['e' num2str(ij(i,1)) num2str(ij(i,2)) 'm' num2str(m) '_' nm '.img']);
VO(i).dim(4) = spm_type('float');
VO(i).pinfo = [1 0 0]';
VO(i).descrip = ['Strain_Tensor(' num2str(ij(i,1)) ',' num2str(ij(i,2)) ') - m=' num2str(m)];
VO(i) = spm_create_image(VO(i));
end;
for p=1:vo.dim(3),
A = get_jacobian(sn3d, vo, p);
E = tensors_from_jacobian(A,m);
for i=1:6,
VO(i) = spm_write_plane(VO(i),E(:,:,i),p);
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function VO = init_vo(sn3d,vox,bb)
if nargin>=2,
x = (bb(1,1):vox(1):bb(2,1))/sn3d.Dims(3,1) + sn3d.Dims(4,1);
y = (bb(1,2):vox(2):bb(2,2))/sn3d.Dims(3,2) + sn3d.Dims(4,2);
z = (bb(1,3):vox(3):bb(2,3))/sn3d.Dims(3,3) + sn3d.Dims(4,3);
dim = [length(x) length(y) length(z)];
origin = -bb(1,:)./vox + 1;
off = -vox.*origin;
mat = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off(3) ; 0 0 0 1];
else,
dim = sn3d.Dims(1,:);
x = 1:dim(1);
y = 1:dim(2);
z = 1:dim(3);
mat = sn3d.MG;
end;
VO = struct('fname','',...
'dim',[dim NaN], 'mat',mat,...
'pinfo',[NaN NaN NaN]',...
'descrip','');
return;
%_______________________________________________________________________
%_______________________________________________________________________
function A = get_jacobian(sn3d, VO, j)
% Each element of the Jacobian matrix field consists of something analagous to the following:
% maple diff( m11*(x1+p11*b1(x1,x2)+p12*b2(x1,x2)) + m12*(x2+p21*b1(x1,x2)+p22*b2(x1,x2)) + m13, x1)
% maple diff( m11*(x1+p11*b1(x1,x2)+p12*b2(x1,x2)) + m12*(x2+p21*b1(x1,x2)+p22*b2(x1,x2)) + m13, x2)
% maple diff( m21*(x1+p11*b1(x1,x2)+p12*b2(x1,x2)) + m22*(x2+p21*b1(x1,x2)+p22*b2(x1,x2)) + m23, x1)
% maple diff( m21*(x1+p11*b1(x1,x2)+p12*b2(x1,x2)) + m22*(x2+p21*b1(x1,x2)+p22*b2(x1,x2)) + m23, x2)
Dims = sn3d.Dims;
MG = sn3d.MG;
MF = sn3d.MF;
Transform = sn3d.Transform;
Affine = sn3d.Affine;
dim = VO.dim(1:3);
mat = VO.mat;
% Assume transverse images, and obtain position of pixels in millimeters,
% and convert to voxel space of template.
%----------------------------------------------------------------------------
x = ((1:dim(1))*mat(1,1) + mat(1,4))/Dims(3,1) + Dims(4,1);
y = ((1:dim(2))*mat(2,2) + mat(2,4))/Dims(3,2) + Dims(4,2);
z = ( j*mat(3,3) + mat(3,4))/Dims(3,3) + Dims(4,3);
X = x'*ones(1,dim(2));
Y = ones(dim(1),1)*y;
bbX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1);
bbY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1);
bbZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1);
dbX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1,'diff');
dbY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1,'diff');
dbZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1,'diff');
M = MF*Affine*inv(MG);
sgn = sign(det(M(1:3,1:3)));
% Nonlinear deformations
%----------------------------------------------------------------------------
% 2D transforms for each plane
tbx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ', Dims(2,1), Dims(2,2) );
tby = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ', Dims(2,1), Dims(2,2) );
tbz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ', Dims(2,1), Dims(2,2) );
tdx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ', Dims(2,1), Dims(2,2) );
tdy = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ', Dims(2,1), Dims(2,2) );
tdz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ', Dims(2,1), Dims(2,2) );
% Jacobian of transformation from template
% to affine registered image.
%---------------------------------------------
j11 = dbX*tbx*bbY' + 1;
j12 = bbX*tbx*dbY';
j13 = bbX*tdx*bbY';
j21 = dbX*tby*bbY';
j22 = bbX*tby*dbY' + 1;
j23 = bbX*tdy*bbY';
j31 = dbX*tbz*bbY';
j32 = bbX*tbz*dbY';
j33 = bbX*tdz*bbY' + 1;
% Combine Jacobian of transformation from
% template to affine registered image, with
% Jacobian of transformation from affine
% registered image to original image.
%---------------------------------------------
A = zeros([size(j11) 3 3]);
A(:,:,1,1) = M(1,1)*j11 + M(1,2)*j21 + M(1,3)*j31;
A(:,:,1,2) = M(1,1)*j12 + M(1,2)*j22 + M(1,3)*j32;
A(:,:,1,3) = M(1,1)*j13 + M(1,2)*j23 + M(1,3)*j33;
A(:,:,2,1) = M(2,1)*j11 + M(2,2)*j21 + M(2,3)*j31;
A(:,:,2,2) = M(2,1)*j12 + M(2,2)*j22 + M(2,3)*j32;
A(:,:,2,3) = M(2,1)*j13 + M(2,2)*j23 + M(2,3)*j33;
A(:,:,3,1) = M(3,1)*j11 + M(3,2)*j21 + M(3,3)*j31;
A(:,:,3,2) = M(3,1)*j12 + M(3,2)*j22 + M(3,3)*j32;
A(:,:,3,3) = M(3,1)*j13 + M(3,2)*j23 + M(3,3)*j33;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function tensor = tensors_from_jacobian(A, m)
d = size(A);
tensor = zeros([d(1:2) 6]);
I = eye(3);
if m==0,
% Hencky
for j=1:d(2),
for i=1:d(1),
J = squeeze(A(i,j,:,:));
T = logm(J'*J)*0.5;
tensor(i,j,:) = T([1 2 3 5 6 9]);
end;
end;
elseif m==2,
%
for j=1:d(2),
for i=1:d(1),
J = squeeze(A(i,j,:,:));
T = 0.5*(J'*J - I);
tensor(i,j,:) = T([1 2 3 5 6 9]);
end;
end;
else,
for j=1:d(2),
for i=1:d(1),
J = squeeze(A(i,j,:,:));
T = ((J'*J)^(m/2)-I)/m;
tensor(i,j,:) = T([1 2 3 5 6 9]);
end;
end;
end;
return;
%_______________________________________________________________________
%_______________________________________________________________________
function spm_write_defs(sn3d, vox,bb)
% Write deformation field.
% FORMAT spm_write_defs(matname, vox,bb)
% sn3d - information from the `_sn3d.mat' file containing the spatial
% normalization parameters.
% The deformations are stored in y1.img, y2.img and y3.img
%_______________________________________________________________________
% %W% John Ashburner %E%
matname = deblank(sn3d.fname);
Dims = sn3d.Dims;
MG = sn3d.MG;
MF = sn3d.MF;
Transform = sn3d.Transform;
Affine = sn3d.Affine;
if nargin>=3,
x = (bb(1,1):vox(1):bb(2,1))/Dims(3,1) + Dims(4,1);
y = (bb(1,2):vox(2):bb(2,2))/Dims(3,2) + Dims(4,2);
z = (bb(1,3):vox(3):bb(2,3))/Dims(3,3) + Dims(4,3);
dim = [length(x) length(y) length(z)];
origin = -bb(1,:)./vox + 1;
off = -vox.*origin;
mat = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off(3) ; 0 0 0 1];
else,
dim = Dims(1,:);
x = 1:dim(1);
y = 1:dim(2);
z = 1:dim(3);
mat = MG;
end;
[pth,nm,xt,vr] = fileparts(deblank(matname));
if length(nm)>6 & strcmp(nm(end-4:end),'_sn3d'), nm=nm(1:end-5); end;
VX = struct('fname',fullfile(pth,['y1_' nm '.img']), 'dim',[dim 16], ...
'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - X');
VY = struct('fname',fullfile(pth,['y2_' nm '.img']), 'dim',[dim 16], ...
'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - Y');
VZ = struct('fname',fullfile(pth,['y3_' nm '.img']), 'dim',[dim 16], ...
'mat',mat, 'pinfo',[1 0 0]', 'descrip','Deformation field - Z');
X = x'*ones(1,VX.dim(2));
Y = ones(VX.dim(1),1)*y;
if (prod(Dims(2,:)) == 0),
affine_only = 1;
basX = 0; tx = 0;
basY = 0; ty = 0;
basZ = 0; tz = 0;
str = 'affine';
else
affine_only = 0;
basX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1);
basY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1);
basZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1);
str = 'nonlinear';
end
spm_create_image(VX);
spm_create_image(VY);
spm_create_image(VZ);
% Cycle over planes
%----------------------------------------------------------------------------
for j=1:length(z)
% Nonlinear deformations
%----------------------------------------------------------------------------
if (~affine_only)
% 2D transforms for each plane
tx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
ty = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
tz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *basZ(j,:)', Dims(2,1), Dims(2,2) );
X1 = X + basX*tx*basY';
Y1 = Y + basX*ty*basY';
Z1 = z(j) + basX*tz*basY';
end
% Sample each volume
%----------------------------------------------------------------------------
Mult = Affine;
if (~affine_only)
X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
else
X2= Mult(1,1)*X + Mult(1,2)*Y + (Mult(1,3)*z(j) + Mult(1,4));
Y2= Mult(2,1)*X + Mult(2,2)*Y + (Mult(2,3)*z(j) + Mult(2,4));
Z2= Mult(3,1)*X + Mult(3,2)*Y + (Mult(3,3)*z(j) + Mult(3,4));
end
spm_write_plane(VX,X2,j);
spm_write_plane(VY,Y2,j);
spm_write_plane(VZ,Z2,j);
end
return;
%_______________________________________________________________________
|
