Hi,
We recommend you use flirt with the -applyxfm option to reapply a flirt
matrix. The conventions of the coordinates used by flirt are not
entirely
trivial, and they do depend a little on the nifti qform/sform
information.
If you are desperate to use your own matlab implementation then there
is a thread on doing exactly this a little while ago (at least 4 months)
in the FSL list.
All the best,
Mark
On 6 May 2010, at 22:19, Dan Golding wrote:
> Hi all,
>
> I need to reapply the affine transformation generated by flirt on
> different data (dti fibres). I think it should be quite straight
> forward. I have been trying to reapply the transformation myself to
> the t1 data that I used flirt on just to make certain I understand
> it properly and I can't seem to get the same output image.
>
> I think it has something to do with the fact that I completely
> ignore all the information in the nifti header. Is there anything in
> there that I have to consider before applying the affine
> transformation?
>
> Here is the matlab script I have been using - perhaps someone can
> see where I have gone wrong?
>
> -----------------------------------------------------------------------------------------------------------------------------------------------------------
> %Load the affine transformation matirx from flirt
> load e2fflirt -ASCII;
> A = e2fflirt;
>
> %load the source image
> source = load_nii('EF1T_t1.nii');
>
> %make an empty output image
> output = source;
> output.img = output.img*0;
>
> Io = output.img;
> Is = source.img;
>
> %find the output image dimensions so we know how many voxels to
> iterate through
> OutDim = size(output.img);
>
> %Iterate through each output voxel
> for i = 1:OutDim(1)
> for j = 1:OutDim(2)
> for k = 1:OutDim(3)
>
> %Perfom the affine transformation on voxel (i,j,k). Note
> that
> %half of the size of each dimension is subtracted, this
> is so
> %that the centre of rotation is in the centre of the image
> %rather than at voxel (0,0,0).
> Y = A*[i - OutDim(1)/2 + 1; j - OutDim(2)/2 + 1; k -
> OutDim(3)/2 + 1; 1];
>
> %This is used for trilinear interpolation
> id = mod(Y(1),1);
> jd = mod(Y(2),1);
> kd = mod(Y(3),1);
>
> %Add half the dimension back so voxel don't half -ve
> %coordinates.
> Y = floor(Y + [OutDim'./2;0]);
>
> try
> %------------------------------------------------
> %Warping output locations to source
>
> %Trilinear -- I've tested this, it definitely works
> C = Is(Y(1):Y(1)+1,Y(2):Y(2)+1,Y(3):Y(3)+1);
> Io(i,j,k) = (1-kd)*((1-jd)*((1-id)*C(1,1,1) +
> id*C(2,1,1)) + jd*((1-id)*C(1,2,1) + id*C(2,2,1))) + kd*((1-jd)*((1-
> id)*C(1,1,2) + id*C(2,1,2)) + jd*((1-id)*C(1,2,2) + id*C(2,2,2)));
>
> %Nearest Neighbour
> %Io(i,j,k) = Is(Y(1), Y(2), Y(3));
>
> catch ERR %in the cases where the vectors point off the
> image, i.e. I,J,K == 0 | 160(200) as the case may be
>
> continue;
> end;
>
> end;
> end;
> end;
>
> %put the image back into a nifti structure. Note that the image is
> warped
> %but the header remains the same as the source image!
> output.img = Io;
> -----------------------------------------------------------------------------------------------------------------------------------------------------------
>
>
> Any suggestions are welcome
>
> Kind Regards
> Dan Golding
>
|