Sorry about that last partial email. I must have pressed the wrong button....
The two attached matlab files should be what you need, although they both
compute the Jacobian determinant field in a different way. Note that the
fields are scaled differently, because one of them corrects for voxel sizes.
Also, if the images are flipped in the spatial normalization, then the
determinant field from write_dets becomes negative.
Good luck,
-John
| Date: Thu, 14 Oct 1999 16:09:38 +1000
| Mime-version: 1.0
| X-Priority: 3
| Content-transfer-encoding: 7bit
| Subject: nonlinear deformations and detecting structural changes
| From: Greig de Zubicaray <[log in to unmask]>
| To: [log in to unmask]
| X-List: [log in to unmask]
| X-Unsub: To leave, send text 'leave spm' to [log in to unmask]
| X-List-Unsubscribe: <mailto:[log in to unmask]>
| Content-Transfer-Encoding: 7bit
|
| Dear SPMers,
|
| I'm hoping to use SPM99b to compare structural images from two groups, and
| would like to use the method described in Gaser et al. (1999). Neuroimage,
| 10, 107-113, based on the analysis of deformation fields in SPM.
|
| I have the code provided by John Ashburner in
| http://www.mailbase.ac.uk/lists/spm/1999-08/att-bin1a004cV that writes the
| deformation fields to x, y, and z images from the sn3d.mat files. However, I
| need to calculate the jacobian determinant of the vectors in the .sn3d.mat
| files to distinguish between volume loss or gain. Does anyone already have
| the MATLAB code for this that writes the output as a single image file that
| can be analysed through the statistics module in SPM99b? Admittedly, I'm
| being lazy asking this, but any assistance would be appreciated.
|
| regards,
|
| Greig de Zubicaray
|
|
| --
| Centre for Magnetic Resonance
| The University of Queensland
| Brisbane
| QLD 4072
| AUSTRALIA
|
| Tel: +61 (0) 7 3365 4250
| Fax: +61 (0) 7 3365 3833
function dets_from_sn3d(matname, vox,bb)
% Modulate spatially normalized images in order to preserve `counts'.
% FORMAT dets_from_sn3d(matname, vox,bb)
% matname - the `_sn3d.mat' file containing the spatial normalization
% parameters.
%_______________________________________________________________________
% %W% John Ashburner %E%
load(deblank(matname))
if (exist('mgc') ~= 1)
error(['Matrix file ' matname ' is the wrong type.']);
end
if (mgc ~= 960209)
error(['Matrix file ' matname ' is the wrong type.']);
end
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 = round(-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;
VO = struct('fname','det1.img',...
'dim',[dim 16], 'mat',mat,...
'pinfo',[1 0 0]',...
'descrip','Determinant field');
% Assume transverse images, and obtain
% position of pixels in millimeters.
%----------------------------------------------------------------------------
x = (1:VO.dim(1))*VO.mat(1,1) + VO.mat(1,4);
y = (1:VO.dim(2))*VO.mat(2,2) + VO.mat(2,4);
z = (1:VO.dim(3))*VO.mat(3,3) + VO.mat(3,4);
% Convert to voxel space of template.
%----------------------------------------------------------------------------
x = x/Dims(3,1) + Dims(4,1);
y = y/Dims(3,2) + Dims(4,2);
z = z/Dims(3,3) + Dims(4,3);
X = x'*ones(1,VO.dim(2));
Y = ones(VO.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');
% Start progress plot
%----------------------------------------------------------------------------
spm_progress_bar('Init',length(z),'Modulating by determinants of Jacobian','planes completed');
M = MF*Affine*inv(MG);
sgn = sign(det(M(1:3,1:3)));
img = zeros(VO.dim(1:3));
% Cycle over planes
%----------------------------------------------------------------------------
for j=1:length(z)
% Nonlinear deformations
%----------------------------------------------------------------------------
% 2D transforms for each plane
tbx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );
tby = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );
tbz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );
tdx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', Dims(2,1), Dims(2,2) );
tdy = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', Dims(2,1), Dims(2,2) );
tdz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', 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.
%---------------------------------------------
J11 = M(1,1)*j11 + M(1,2)*j21 + M(1,3)*j31;
J12 = M(1,1)*j12 + M(1,2)*j22 + M(1,3)*j32;
J13 = M(1,1)*j13 + M(1,2)*j23 + M(1,3)*j33;
J21 = M(2,1)*j11 + M(2,2)*j21 + M(2,3)*j31;
J22 = M(2,1)*j12 + M(2,2)*j22 + M(2,3)*j32;
J23 = M(2,1)*j13 + M(2,2)*j23 + M(2,3)*j33;
J31 = M(3,1)*j11 + M(3,2)*j21 + M(3,3)*j31;
J32 = M(3,1)*j12 + M(3,2)*j22 + M(3,3)*j32;
J33 = M(3,1)*j13 + M(3,2)*j23 + M(3,3)*j33;
% The determinant of the Jacobian reflects relative volume changes.
%------------------------------------------------------------------
detJ = J11.*(J22.*J33 - J23.*J32) - J21.*(J12.*J33 - J13.*J32) + J31.*(J12.*J23 - J13.*J22);
img(:,:,j) = detJ*sgn;
spm_progress_bar('Set',j);
end
% Write the data
%------------------------------------------------------------------
VO = spm_write_vol(VO, img);
spm_progress_bar('Clear');
return;
%_______________________________________________________________________
function write_dets(V)
if nargin==0,
V = spm_vol(spm_get(3,'*_y?.img','Specify deformation field'));
end;
VO = V(1);
VO.fname = 'dets.img';
VO.pinfo = [1 0 0]';
VO.descrip = 'Determinant field';
VO = spm_create_image(VO);
[prev1,prev2,prev3] = read_slices(V,1);
[cur1, cur2, cur3 ] = read_slices(V,2);
write_plane(VO,1,cur1,cur2,cur3,prev1,prev2,prev3,1);
pprev1 = prev1;
pprev2 = prev2;
pprev3 = prev3;
prev1 = cur1;
prev2 = cur2;
prev3 = cur3;
for p=3:VO.dim(3),
[cur1, cur2, cur3]= read_slices(V,p);
write_plane(VO,p-1,cur1,cur2,cur3,pprev1,pprev2,pprev3,2);
pprev1 = prev1;
pprev2 = prev2;
pprev3 = prev3;
prev1 = cur1;
prev2 = cur2;
prev3 = cur3;
end;
write_plane(VO,VO.dim(3),prev1,prev2,prev3,pprev1,pprev2,pprev3,1);
return;
function write_plane(VO,p,cur1,cur2,cur3,pprev1,pprev2,pprev3,d)
[J12,J11] = gradient(cur1);
[J22,J21] = gradient(cur2);
[J32,J31] = gradient(cur3);
J13 = (cur1-pprev1)/d;
J23 = (cur2-pprev2)/d;
J33 = (cur3-pprev3)/d;
dt = J11.*J22.*J33-J11.*J23.*J32-J21.*J12.*J33+J21.*J13.*J32+J31.*J12.*J23-J31.*J13.*J22;
spm_write_plane(VO,dt,p);
return;
function [cur1, cur2, cur3]= read_slices(V,p)
M = spm_matrix([0 0 p]);
cur1 = spm_slice_vol(V(1),M,V(1).dim(1:2),1);
cur2 = spm_slice_vol(V(2),M,V(1).dim(1:2),1);
cur3 = spm_slice_vol(V(3),M,V(1).dim(1:2),1);
return;
|