Hi Doug and list, I'm trying to map surface/mesh vertices (e.g. rh.pial) to voxels (cols, rows, slices) in a functional volume in matlab, and just wanted to a) get a check that what I'm doing sounds right and b) share it in case it's useful to anyone else (I searched the mailing list and found most of my answers there, just not all in one post). Of course, mri_surf2vol does all this, probably better, but I want to iterate through many, many times. I've tested it and it seems to work.
Functional data has been registered to anatomical (T1.mgz? orig.mgz? w/e) using bbregister.
Doug Greve 19 Dec 2008 said:
if you have the voxel index (col, row, slice) then you can:
V2R = [ -1.00000 0.00000 0.00000 128.00000 0.00000 0.00000 1.00000 -128.00000 0.00000 -1.00000 0.00000 128.00000 0.00000 0.00000 0.00000 1.00000 ]
crs = [col row slice]';
xyz = V2R * [crs+1; 1];
then find the vertex whose xyz coords are closest to xyz
So: to find CRS from vertex RAS:
crs = round(inv(V2R)*[xyz; 1]) %maybe crs -1 to make a closer match? but then this messes with below
------------------------------------------------------------------------------
From fscoordinates.pdf page 6, plus above...
To go from Func CRS to T1 CRS (& vice versa)
bbmtx is the transform matrix from register.dat
Q = inv(Func_tkrvox2ras)*bbmtx*T1_tkrvox2ras;
Q is now the transform from Func_tkrvox2ras to T1tkrvox2ras
To find T1 CRS from Func CRS: T1_crs = inv(Q)*[Func_crs; 1];
To find Tfunc CRS from T1 CRS:
Func_crs = Q*[T1_crs; 1]; %T1_crs + 1 if you subtracted 1 above
So...to go from vertex RAS to Func CRS:
Func_crs = Q*(round(inv(V2R)*[xyz; 1]);
This seem right?
best, Jim T