> Could you give me some details on how to make a template from a series
> of monkey PET scans? I found your script for making templates in the
> mailing list, but the scans are apparently normalized to human
> proportions. Can the script be modified to accommodate monkey anatomy?
> Any help you can provide on this would be wonderful.
You'll probably regret asking this, but....
To begin with, you'll need to define some standard space, or use some
existing space. If you have the chance to start your own co-ordinate system
from scratch and have the necessary mathematical and programming expertise
in your lab, then this space should relate to the average brain or head
size. You also need to define which orientation this template should be in.
Doing this properly would require writing some Matlab code. The simplest
case would just involve defining a space that is based on affine
registrations.
-------------- Creating your own co-ordinate system ---------------
Personally, I would go for an iterative approach, involving affine
registering all the images together. The affine transformations would then
be averaged, and this average factored out from the individual
transformations, before writing out the affine transformed images. These
would then be averaged, and the average used as a template for the next
iteration. The iterations would continue until everything stopped changing.
A procedure similar to the one above could give you an intensity averaged,
affine shape averaged image. You then need to decide how it should be
rotated and translated into your standard space. Perhaps you would
translate it so that the AC is at the origin, the PC lies behind it with the
same Z co-ordinate, and the template is as left-right symmetrical as
possible.
An alternative would be to simultaneously register each image to every other
image, incorporating constraints as in:
{Ashburner J & Friston KJ (1997): "Multimodal image coregistration and
partitioning - a unified framework." NeuroImage. 6(3):209-217
Using M2 to denote the mapping from image I1 to image I2, M3 for the mapping
I1 to I3, etc. Affine registering 4 images together requires 12x3 parameters
as in....
I1->I2 M2
I1->I3 M3
I1->I4 M4
I2->I3 M3*inv(M2)
I2->I4 M4*inv(M2)
I3->I4 M4*inv(M3)
If the affine registration is not framed to be symmetric, then the following
should also be included:
I2->I1 inv(M2)
I3->I1 inv(M3)
I4->I1 inv(M4)
I3->I2 M2*inv(M3)
I4->I2 M2*inv(M4)
I4->I3 M3*inv(M4)
Both approaches should give similar results when the objective function
involves minimising a sum of squares difference. This works because the
objective function that is minimised when registering all images with
each other in the most consistent way:
(a-b)^2 + (a-c)^2 + (a-d)^2 + (b-c)^2 + (b-d)^2 + (c-d)^2
is equivalent to (a scaled version) of the objective function for registering
all images to their own mean:
4*( (a-(a+b+c+d)/4)^2 + (b-(a+b+c+d)/4)^2 + (c-(a+b+c+d)/4)^2 +
(d-(a+b+c+d)/4)^2 )
The above equivalence also explains why I prefer to use a template for
spatial normalisation that is based on the average of loads of images,
rather than based on a single individual. It enforces more consistency
among the parameter estimates.
The tricky thing depends on how you define the average transformation that
should be factored out. There is a certain amount of disagreement about how
this should be done, and what exactly the average means. My own approach
would involve minimising the individual deviations from the average by
penalising shears and zooms, but not rigid body displacements.
Each of the above matrices (M2, M3, and M4) can be thought of as encoding a
shape change from (e.g.) I1 to the average, a rigid rotation of the average
shape, and then a final shape change from the average to (e.g.) I2.
i.e.
M2 = Z1*R2*inv(Z2)
M3 = Z1*R3*inv(Z3)
M4 = Z1*R4*inv(Z4)
The above would then be constraints for an optimisation that minimises the
sum of squares of the logs of the singular values of matrices Z1, Z2, Z3 and
Z4. This can be framed in terms of estimating the 6 parameters required to
define Z1 (see Section 4.2.3 of
http://www.fil.ion.ucl.ac.uk/~john/thesis/chapter4.pdf and Section 6.4.1 of
http://www.fil.ion.ucl.ac.uk/~john/thesis/chapter6.pdf).
Estimating the minimum deformation matrices be done with the following piece
of code (unless I have made a mistake somewhere)....
% This bit estimates the parameters
M = {M2,M3,M4};
p = spm_powell(zeros(6,1), eye(6)/100,ones(1,6)*1e-6,'myfun',M);
% This bit builds up the matrices
clear Z R
Z{1} = [expm([p(1:3)' ; p(2) p(4:5)' ; p(3) p(5) p(6)]) zeros(3,1) ; 0 0 0
1]; R{1} = eye(4);
for i=1:length(M),
RIZ = Z{1}\M{i};
Z{i+1} = [inv(sqrtm(RIZ(1:3,1:3)'*RIZ(1:3,1:3))) zeros(3,1) ; 0 0 0 1];
R{i+1} = RIZ*Z{i+1};
end;
The spm_powell optimiser calls the following function, which returns
the sum of squares of the logs of the singular values of the various
matrices. This should be saved to a file called myfun.m.
------- start of myfun.m -------
function fval = myfun(p,M)
L1 = [p(1:3)' ; p(2) p(4:5)' ; p(3) p(5) p(6)];
Z1 = expm(L1);
fval = sum(L1(:).^2);
for i=1:length(M)
RIZ = Z1\M{i}(1:3,1:3);
tmp = RIZ'*RIZ;
fval = fval + sum(log(svd(tmp)).^2)/4;
end;
------- end of myfun.m -------
The affine mapping from e.g. I2 to the average shape would then be:
inv(R{2})*Z{2}
Note that the orientation of this space is the same as that of image I1,
but the brains are scaled to their average size.
For usage hints for the affine registration code in SPM2b, see the email I
sent yesterday. Note also, that in addition to minimising the mean squared
difference between the images, the affine registration in SPM2b can also be
made to simultaneously minimise the shape changes by minimising the sum of
squares of the logs of the singular values of the affine transformation
matrices. The inclusion of this regularisation may obviate the need to do
the mean correction at each iteration.
-------------- Using an existing co-ordinate system ---------------
If you want to adopt an existing co-ordinate system, then you would need to
somehow get your images into this co-ordinate system. This may involve
picking out certain points, and using these to determine how the images
should be transformed into this space. This requires a minimum of 4 points
in order to get a 12 parameter affine transformation, but more points would
be better.
A further possibility is to find an existing template, register your images
to this, and use the mean of these registred images as your template.
e.g. http://www.jiscmail.ac.uk/cgi-bin/wa.exe?A2=ind0201&L=spm&P=R27107&I=-3
Best regards,
-John
--
Dr John Ashburner.
Functional Imaging Lab., 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
|
