Thanks a bunch John!
I will try your suggestions, they look very efficient.
I also got a bit further with this today, adapting your
get_orig_coord5.m to preallocate xyz2 for speed, and taking 1 reshape
operation out of the loop over coordinates, as it was always doing the
same thing. It already speed up the process ~17-fold, bringing it to
'workable' processing durations for large sets of DTI fibers. I think I
could also do without the loop altogether, applying deformations in 1 go
to the entire basX,basY,basZ, but for that I have to dive in it a bit
deeper.
I adapted get_orig_coord5.m as follows (see % BAS: comments):
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(d(1),size(Tr,1),xyz(1,:)-1);
basY = spm_dctmtx(d(2),size(Tr,2),xyz(2,:)-1);
basZ = spm_dctmtx(d(3),size(Tr,3),xyz(3,:)-1);
end;
xyz2=zeros(3,size(xyz,2)); %BAS: now preallocate for speed
if affine_only,
xyz2 = Mult(1:3,:)*[xyz ; ones(1,size(xyz,2))];
else,
rTr1=reshape(Tr(:,:,:,1),size(Tr,1)*size(Tr,2),size(Tr,3)); %BAS:
removed from loop below (for i=1:..)
rTr2=reshape(Tr(:,:,:,2),size(Tr,1)*size(Tr,2),size(Tr,3)); %BAS:
removed from loop below (for i=1:..)
rTr3=reshape(Tr(:,:,:,3),size(Tr,1)*size(Tr,2),size(Tr,3)); %BAS:
removed from loop below (for i=1:..)
for i=1:size(xyz,2),
bx = basX(i,:);
by = basY(i,:);
bz = basZ(i,:);
tx = reshape(rTr1*bz', size(Tr,1), size(Tr,2) ); %BAS: rTr1
instead of time consuming reshape
ty = reshape(rTr2*bz', size(Tr,1), size(Tr,2) ); %BAS: rTr2
instead of time consuming reshape
tz = reshape(rTr3*bz', size(Tr,1), size(Tr,2) ); %BAS: rTr3
instead of time consuming reshape
xyz2(:,i) = Mult(1:3,:)*[xyz(:,i) + [bx*tx*by' ; bx*ty*by' ;
bx*tz*by']; 1];
end;
end;
orig_coord = xyz2';
return;
Op 11-03-10 14:17, John Ashburner schreef:
> Perhaps it would be easier just to write the entire deformation field
> out using the Deformations utility. You can then sample voxels in the
> resulting deformation to see where they point to. Note that they point
> to some mm coordinates, so you may need to convert these to voxels (via
> the inverse of the voxel-to-world mapping for the image whose voxels you
> are after).
>
> P=spm_select(1,'y_.*\.nii');
> N=nifti(P);
>
> tx = zeros(size(xyz));
> for i=1:size(xyz,2)
> tx(:,i) = squeeze(N.dat(xyz(1,i),xyz(2,i),xyz(3,i),1,:));
> end
>
>
> Note that the above just works for integer values. You could also do it
> via:
>
> P=spm_select(1,'y_.*\.nii');
> N=nifti(P);
> X=N.dat(:,:,:,1,1);
> Y=N.dat(:,:,:,1,2);
> Z=N.dat(:,:,:,1,3);
>
> tx = [spm_bsplins(X,xyz(1,:),xyz(2,:),xyz(3,:),[1 1 1 0 0 0])
> spm_bsplins(Y,xyz(1,:),xyz(2,:),xyz(3,:),[1 1 1 0 0 0])
> spm_bsplins(Z,xyz(1,:),xyz(2,:),xyz(3,:),[1 1 1 0 0 0])];
>
>
> Converting to voxel indices then requires the voxel-to-world mapping of
> the image:
>
> P=spm_select(1,'nifti');
> M=spm_get_space(P);
> M=inv(M);
> tx_vox = M(1:3,:)*[tx; ones(1,size(tx,2))];
>
>
> Best regards,
> -John
>
>
> we are often using the function spm_get_orig_coord.m in spm5 to
> normalize lists of x,y,z coordinates from native space to MNI
> space (based on normalization parameters obtained using unified
> segmentation). This is of course not the same as normalization
> of image volumes for which one uses spm_normalise.m.
>
> It does what it has to do, so far so good.
>
> However, for non-linear warping the procedure seems to use a
> for-loop going through all coordinates in the list (line 50 and
> below):
>
> if affine_only,
> xyz2 = Mult(1:3,:)*[xyz ; ones(1,size(xyz,2))];
> else,
> for i=1:size(xyz,2),
> bx = basX(i,:);
> by = basY(i,:);
> bz = basZ(i,:);
> tx = reshape(...
>
> reshape(Tr(:,:,:,1),size(Tr,1)*size(Tr,2),size(Tr,3))...
> *bz', size(Tr,1), size(Tr,2) );
> ...
>
> So with long lists of coordinates (in our case, DTI fiber
> coordinates) this is really slow, taking a day or so for 1
> subject.
>
> Is there a way to parallelize this whole procedure, or is there
> a faster way now for doing this? Compiling to mex files does not
> help anymore I read, as matlab now seems to use Just-In-Time
> compilation.
>
> Any suggestion is welcome!
>
>
> --
> John Ashburner<[log in to unmask]>
>
>
--
--------------------------------------------------
Dr. S.F.W. Neggers
Division of Brain Research
Rudolf Magnus Institute for Neuroscience
Utrecht University Medical Center
Visiting : Heidelberglaan 100, 3584 CX Utrecht
Room B.01.1.03
Mail : Huispost B01.206, P.O. Box
3508 GA Utrecht, the Netherlands
Tel : +31 (0)88 7559609
Fax : +31 (0)88 7555443
E-mail : [log in to unmask]
Web : http://www.neuromri.nl/people/bas-neggers
--------------------------------------------------
------------------------------------------------------------------------------
De informatie opgenomen in dit bericht kan vertrouwelijk zijn en is
uitsluitend bestemd voor de geadresseerde. Indien u dit bericht onterecht
ontvangt, wordt u verzocht de inhoud niet te gebruiken en de afzender direct
te informeren door het bericht te retourneren. Het Universitair Medisch
Centrum Utrecht is een publiekrechtelijke rechtspersoon in de zin van de W.H.W.
(Wet Hoger Onderwijs en Wetenschappelijk Onderzoek) en staat geregistreerd bij
de Kamer van Koophandel voor Midden-Nederland onder nr. 30244197.
Denk s.v.p aan het milieu voor u deze e-mail afdrukt.
------------------------------------------------------------------------------
This message may contain confidential information and is intended exclusively
for the addressee. If you receive this message unintentionally, please do not
use the contents but notify the sender immediately by return e-mail. University
Medical Center Utrecht is a legal person by public law and is registered at
the Chamber of Commerce for Midden-Nederland under no. 30244197.
Please consider the environment before printing this e-mail.
|