Fraser Smith wrote:
> Hi Ged,
>
> Cheers for the answer on this. It seems the mapping is pretty isomorphic
> from voxel space to mm space.
Ah, I thought you meant how does the mapping work in terms of
code/mechanics of it. If you mean mathematically, then the answer is
that the spm_read_vols locations use the mapping stored in the vol.mat
field, which is a 12 degree of freedom affine transformation, i.e. the
three-vector of voxel indices is transformed to the three-vector of
world/mm coordinates by a matrix multiplication and a vector addition.
This affine transformation is represented as a single 4-by-4 matrix
(with bottom row [0 0 0 1]) that multiplies a homogeneous 4-vector --
i.e. the voxel indices with a 1 tacked on the end. This is purely for
convenience in this case, since a general 15 DoF homogeneous 4-by-4
matrix cannot be used with SPM. [see also line 7 of my code snippet]
Hence, it is an isomorphism, and indeed a diffeomorphism (with
constant Jacobian matrix given by the 3-by-3 linear part of the 4-by-4
affine matrix).
Hope that helps,
Ged.
> On 17 Nov 2006, at 17:50, Ged Ridgway wrote:
>
>> Fraser Smith wrote:
>>> If I read an image into a matrix A, using spm_read_vols, I get an
>>> XbyYbyZ image, plus a set of XYZ locations returned.
>>> My question is, how exactly do the locations correspond to the
>>> individual voxels in the the XbyYbyZ image?
>>
>> Hi Fraser,
>>
>> The rows of the XYZ matrix are x,y,z (world) coordinates in mm, the
>> columns correspond to the voxels in Matlab's linear-indexing
>> convention (i.e. fastest changing dimension down the columns, then
>> along the rows, then through planes slowest changing dimension).
>>
>> The following might clarify this:
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> V = spm_vol('image.img');
>> [img XYZ] = spm_read_vols(V);
>> vox = 12345; % 12345th voxel (linear index)
>> % convert from linear index to subscripts
>> [x y z] = ind2sub(size(img), vox);
>> XYZ(:, vox) % world coords at voxel vox
>> V.mat*[x y z 1]' % same at voxel (x,y,z)
>> % note that vox and (x,y,z) refer to the same voxel:
>> img(x,y,z)
>> img(vox)
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>
>> Personally, I find the V.mat*[x y z 1]' method of getting world
>> coordinates at voxel (x,y,z) easier than using linear indexing into
>> the XYZ variable anyway.
>>
>> Best,
>> Ged.
>
>
|