Hi Keith,

I'm glad that you can replicate the non-cropped version - I think
that is the biggest battle won.  What you really are asking now is
how to modify the flirt matrix for a cropped volume, since if you
had the appropriate flirt matrix then you would be able to get it
to work now.  This is a much easier problem and can be solved
like this:

Say that the 0,0,0 (voxel) coordinate in your cropped image was
originally the x0,y0,z0 (voxel) coordinate, and just ignore qforms.
You can relate these two coordinate systems by a voxel-wise
transformation that is simply
   Mc = [ 1 0 0 -x0 ; 0 1 0 -y0 ; 0 0 1 -z0 ; 0 0 0 1 ]
which converts voxel coords in the original volume to those in
the cropped volume.  If you've worked out the appropriate
matrix for voxel coordinate transforms (what I called "M" in
my previous email) then you can simply concatenate these
transforms by multiplying the matrices.  That is:
   Mc * M * Mc^{-1}
would map a voxel coordinate in the cropped volume to the
equivalent coordinate in the non-cropped volume, then
transform this into the reference volume and then map it
back to the cropped coordinate (assuming that you are
cropping the reference in the same way - if not, just change
the Mc at the front to be the appropriate cropping for the
reference volume).  You should then be able to apply this
matrix in matlab in the way that you can already do with no
cropping.

Hope that's clear.
All the best,
	Mark

P.S. Yes, the double printing of the last coordinate is a minor
bug - sorry.


On 29 Jul 2009, at 18:44, Keith Schneider wrote:

> Hi Mark,
>
> On Jul 29, 2009, at 3:06 AM, Mark Jenkinson wrote:
>
>> Hi,
>>
>> I don't know why you are using 0 for the last entry in (x,y,z,0) -  
>> it should be (x,y,z,1) and
>> it is this in img2imgcoord and all other calls in FSL.
>
> Okay, you're right, I checked this again.  BTW, img2imgcoord prints  
> the last set of coordinates twice (minor bug?).
>
>> Apart from that, what you have written makes sense and describes  
>> what is done in
>> applyxfm4D or any other FSL tool like FLIRT.  However, because you  
>> are trying
>> to use this in matlab, with a volume that is cropped and another  
>> function (tformarray)
>> then there are several things which may be going wrong.  Firstly,  
>> be aware that all of
>> the FSL coordinate conventions are like nifti - they start at 0 -  
>> whereas in matlab
>> will start at 1.  Also, be careful that you are giving it the  
>> matrix in the right direction
>> (FLIRT stores a mapping *from* input *to* reference space, but your  
>> matlab function
>> may want the opposite).  In addition, be careful whether you need  
>> to give it a mapping
>> between voxel coordinates or between mm coordinates.
>>
>> The qform is not used at all by FSL programs for spatial resampling  
>> *except* to check
>> for radiological/neurological storage (former has negative  
>> determinant, the latter has
>> positive determinant).
>
> The transformation matrix must assume the origin is (0,0,0).  But  
> when I crop the volume, the origin is now somewhere else, and so if  
> I attempt to apply the original transformation to the cropped  
> volume, I will get the wrong answer, right? (because the rotations  
> are no longer around the correct point)
>
> So I've been thinking that I should subtract the offsets (translate)  
> then apply the flirt transformation, then add the offsets back  
> again.  Hence, I should use the matrix you call D (which includes  
> the offset translations) rather than S (which includes only the  
> voxel scaling).
>
>> It isn't easy to convert between conventions for spatial  
>> coordinates and coordinate
>> transformations - I know because I've had to do this in the past  
>> and it is painful.
>> Can you get the correct answers for the actual coordinate locations  
>> (not resampling
>> of images) between the matlab and fsl representations?  For  
>> example, using img2imgcoord
>> and matrix maths in matlab?  If so, then it is probably just a  
>> question of working out
>> what tformarray requires - and I'm not familiar with this so I  
>> cannot help.
>
> Yes, img2imgcoord is computing inv(S) * A * S on zero-based voxel  
> coordinates, I can replicate that, it is just the offsets that are  
> giving trouble, I think.
>
> keith
>