|
| I would like to know how to set parameters in spm_dbm.m
| when I try to use it to examine hemispheric asymmetry.
To test asymmetry, you need to use a symmetric template
and weighting image (the default is brainmask.img which
is asymmetric). The basis functions of a discrete cosine
transform consist of symmetric and anti-symmetric bases,
so tests for assymmetry are done entirely on the
anti-symmetric ones.
The tests done by spm_dbm.m are multivariate, so they just
test for a global difference.
|
| My interests are:
| (1) to examine whether the brains of a group show any
| significant hemispheric asymmetry of shape.
|
| and
|
| (2) to examine whether the brains of two groups show
| any significant differences in the tendency of
| hemisperic asymmetry of shape.
|
| The number(s) of subjects: case (1): 24
| case (2): 8 and 16
| (The subjects in the two cases are the same.)
|
| Normalization: affine+nonlinear(11*13*10)
I hate to tell you this, but spm_dbm reparameterises the deformation
fields (after removing a rigid body transform and an isotropic zoom)
by only 8x8x8 basis functions. Therefore, only about 1/3 of your
coefficiants will be used.
Because the tests are multivariate over the whole images, these
coefficients are further reduced by singular value decomposition.
After SVD, the deformation fields are only represented by a
small number of parameters (as the number of parameters per image
needs to be much less than the number of images in the study).
|
| Using the xtract_shapes function in spm_dbm.m, I made
| a 24subs_shape.mat file which contains 24 subjects'
| normalization data.
|
| Then I tried to use run_stats in spm_dbm.m. After choosing
| the file 24subs_shape.mat, I was asked to set
| 'Interesting Group weights',
These are effects of interest. If you have a shape.mat file for
the two groups, then you would enter:
1 0; 0 1
If you had four groups that you wanted to identify a difference
among, then you would enter: 1 0 0 0; 0 1 0 0; 0 0 1 0;0 0 0 1
A more complex example may involve looking for handedness
differences where you have male and female subjects. If the
groups were FR (female right-handed), FL, MR & ML, then you
would enter: 1 0 1 0 ; 0 1 0 1
However, you would include sex as a confounding effect by
including this in the 'Uninteresting Covariates' as:
1 1 0 0; 0 0 1 1
or you may also model sex-by-handedness interractions here by:
1 1 0 0; 0 0 1 1; 1 0 0 1; 0 1 1 0
| 'Interesting Covariates',
These are just covariates of interest. So (for example) to include
age and age squared as covariates of interest when you have a column
vector of ages in the workspace, you may enter something like:
a' ; a'.^2
| 'Uninteresting Group weights',
| 'Uninteresting Covariates', and
See above
| 'Restrict to asymmetries?'.
| I guess I should enter Yes in the last, but I am not sure what
| the first four mean. Should I enter, say, all 1s for the case (1),
| and eight 1s and sixteen -1s for the case (2)?
When you restrict to asymmetries, all the basis function coefficients
pertaining to symmetric warps are discarded.
|
| The function run_stats in spm_dbm.m seems to output a P value
| based on chi-square test. Is it the probability of the test agaist
| the hypothesis that there is no over-all(whole-hemisphere) asymmetry?
It is a test with a null hypothesis that the effects of interest do not
account for a significant amount of variance in the data.
| If so, how can we see local asymmetries? By using the function
| show_difference in spm_dbm.m?
show_differences shows a caricature of the differences among the groups.
It is either the absolute differences, scaled by some amount, or it can
be differences that are characterised by canonical variates analysis.
I suspect that what you really want is some kind of SPM that localises
significant shape differences (as done by Gaser et al). The simplest
way to do this would be to do tests of local volume differences, which
can be derived from the Jacobian determinants of the deformation fields.
Somewhere in the attached file is a routine that will write images of
Jacobian determinants.
Best regards,
-John
function spm_Deformations(varargin)
% 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
if nargin==0,
a1 = spm_input('Write what?',1,'m',...
['Deformations|Jacobian Matrices|Jacobian Determinants|' ...
'Strain Tensors'],...
[1 2 3 4],1);
arg = allowable(a1);
if strcmp(arg,'tensor')
tensorder = spm_input('Write what?',2,'m',...
['Strain Tensors, m=-2 (Almansi)|'...
'Strain Tensors, m=-1 |'...
'Strain Tensors, m= 0 (Hencky) |'...
'Strain Tensors, m= 1 (Biot) |'...
'Strain Tensors, m= 2 (Green) '],...
[-2 -1 0 1 2],3);
end;
else,
arg = lower(varargin{1});
switch arg,
case 'tensor',
if nargin >1,
tensorder = varargin{2};
else,
tensorder = 0;
end;
case {'def','jacmat','jacdet'},
otherwise,
error('Unknown argument');
end;
end;
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;
switch arg,
case 'tensor'
spm_progress_bar('Init',length(sn3d),...
['Computing Tensors (m= ' num2str(tensorder) ')'],...
'volumes completed');
case 'def'
spm_progress_bar('Init',length(sn3d),...
['Computing Deformations'],'volumes completed');
case 'jacmat'
spm_progress_bar('Init',length(sn3d),...
['Computing Jacobian Matrices'],...
'volumes completed');
case 'jacdet'
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;
switch arg,
case 'tensor',
write_tensor(sn3d(i),Vox,bb,tensorder);
case 'def',
spm_write_defs(sn3d(i),Vox,bb);
case 'jacmat',
write_jacobianm(sn3d(i),Vox,bb);
case 'jacdet',
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 = MF*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;
%_______________________________________________________________________
|
