> I'm interested in projecting coordinates from a normalized image back
> to the original image coordinates. Does anyone have any advice about this?
The deformation fields estimated during spatial normalisation are actually a
mapping from the spatially normalised image back to the original
un-normalised one, so in principle, mapping a point in a normalised one back
to the original image co-ordinates is reasonably straightforward.
The way a spatially normalised image is created begins with an empty
normalised volume. Each voxel in this volume is filled in by working out
where it should project to in the un-normalised volume. The intensity of the
un-normalised image is sampled at this point, and inserted into the
normalised image.
The alternative would be to work through the original image, and put each
voxel in the appropriate position of the normalised version. This is more
difficult though, as it may leave holes in the normalised image where no data
has been put.
I have attached a function to this email that will hopefully do what you want.
Best regards,
-John
--
Dr John Ashburner.
Wellcome Department of Cognitive Neurology.
12 Queen Square, London WC1N 3BG, UK.
tel: +44 (0)20 78337491 or +44 (0)20 78373611 x4381
fax: +44 (0)20 78131420
http://www.fil.ion.ucl.ac.uk/~john
mail: [log in to unmask]
____________________________________________________________________________
____________________________________________________________________________
function orig_coord = get_orig_coord(coord, matname,PN,PU)
% Determine corresponding co-ordinate in un-normalised image.
% FORMAT orig_coord = get_orig_coord(coord, matname,PN,PU)
% coord - [x y z] in normalised image (voxel).
% matname - File containing transformation information (_sn3d.mat).
% PN - Name of normalised image
% PU - Name of un-normalised image
% orig_coord - Co-ordinate in un-normalised image (voxel).
%_______________________________________________________________________
% %W% John Ashburner %E%
[Dims,Affine,MF,MG,Tr] = load_params(matname);
VN = spm_vol(PN);
VU = spm_vol(PU);
Mat = MG\VN.mat;
coord = coord(:);
x = Mat(1,:)*[coord' 1]';
y = Mat(2,:)*[coord' 1]';
z = Mat(3,:)*[coord' 1]';
if (prod(size(Tr)) == 0),
affine_only = 1;
basX = 0; tx = 0;
basY = 0; ty = 0;
basZ = 0; tz = 0;
else,
affine_only = 0;
basX = spm_dctmtx(Dims(1,1),size(Tr,1),x-1);
basY = spm_dctmtx(Dims(1,2),size(Tr,2),y-1);
basZ = spm_dctmtx(Dims(1,3),size(Tr,3),z-1);
end;
Mult = VU.mat\MF*Affine;
if affine_only,
x2 = Mult(1,:)*[x y z 1]';
y2 = Mult(2,:)*[x y z 1]';
z2 = Mult(3,:)*[x y z 1]';
else,
tx = reshape(...
reshape(Tr(:,:,:,1),size(Tr,1)*size(Tr,2),size(Tr,3))...
*basZ', size(Tr,1), size(Tr,2) );
ty = reshape(...
reshape(Tr(:,:,:,2),size(Tr,1)*size(Tr,2),size(Tr,3))...
*basZ', size(Tr,1), size(Tr,2) );
tz = reshape(...
reshape(Tr(:,:,:,3),size(Tr,1)*size(Tr,2),size(Tr,3))...
*basZ', size(Tr,1), size(Tr,2) );
x1 = x + basX*tx*basY';
y1 = y + basX*ty*basY';
z1 = z + basX*tz*basY';
x2 = Mult(1,:)*[x1 y1 z1 1]';
y2 = Mult(2,:)*[x1 y1 z1 1]';
z2 = Mult(3,:)*[x1 y1 z1 1]';
end;
orig_coord = [x2 y2 z2];
return;
function [Dims,Affine,MF,MG,Tr] = load_params(matname)
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
Tr = reshape(Transform,[Dims(2,:) 3]);
return;
|