Amir M. Tahmasebi wrote:
> So here is what I am trying to do: Suppose I have two datasets, a
> template data (name it as A) and another not normalized data (name it as
> B). I am interested in finding a particular slice image in B that
> corresponds to a chosen slice in A. So first I do the registration using
> spm_coreg function. Lets say I end up with a transformation matrix T_r
> (registering B to A).

Okay, let's make sure we're clear here, at this point, T_r maps from
A's world coords to B's world coords. So:
  [BX;BY;BZ;1] = T_r * [AX;AY;AZ;1] % world coords (mm)
and:
  [Bx;By;Bz;1] = B.mat \ T_r * A.mat * [Ax;Ay;Az;1] % voxel coords

> Now I locate the desired slice in A, save the
> orientation matrix (T_d_A) and calculate the orientation in B's space:
> T_d_B = inv(T_r) * T_d_A. The way I locate the slice in A is by using
> the orthoviews and GUI's reorient button. I mean I use manual
> reorientation to display the desired slice and then press "reorient"
> button. I do *NOT* change the location of blobs using mouse.

Now, by changing A's reorientation (e.g. setting roll to 0.2), you
change its voxel-world mapping (and this will be permanent if you
click "reorient images" and select it -- I would avoid this). What
happens is that A.mat gets replaced with R * A.mat, where R is
spm_matrix(params_entered). I think my R is your T_d_A, but I have
changed notation in case you meant T_d_A to be R*A.mat or something else.

There is a little known feature here that might help, which is that
with the image showing in Display (or Check Reg) you can say:
   global st
and then have access to the image's current R from the reorientation
parameters:
   R = st.vols{1}.premul

At this point you can check that if you enter arbitrary voxel
coordinates in the display GUI, you can match the world coords that
the GUI shows with:
   R*A.mat*[Ax;Ay;Az;1]
This hopefully answers your question about whether orth_views is doing
something weird (I don't think it is; I personally think spm_slice_vol
is weird...)

Note, having changed the orientation, the image you are looking at is
no longer world-registered to B, so let's call the new image C, and
note that:
  [CX;CY;CZ;1] = R * [AX;AY;AZ;1] % world coords (mm)
               = R * inv(T_r) * [BX;BY;BZ;1]
               = R * inv(T_r) * B.mat * [Bx;By;Bz;1]

It's my understanding that you are now looking at the plane CZ=0, and
you want to find the corresponding plane of B, right?

Since we have a mapping from C's world-space to B's voxel-space, we
should certainly be able to do this...

Here however, is where I would use spm_sample_vol, rather than
spm_slice_vol. As I said, I find the latter very confusing, so I'm
afraid you will have to rely on someone else on the list (e.g. John,
if he's around?) to help you with that. To use spm_sample_vol though,
we just need to build up matrices of coordinates of the slice CZ=0,
transform these to B's voxel-space, and then pass these to
spm_sample_vol. Something like this should work:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  [CX CY] = meshgrid(-100:100,100:-1:-150); CZ = 0 * CX;
  % I chose the X and Y ranges above (in mm) roughly, by looking at
  % the corners in the GUI Display (could *calculate* bounding box)
  Cworld = [CX(:) CY(:) CZ(:) ones(numel(CX), 1)]';
  Bvoxel = (R * inv(T_r) * B.mat) \ Cworld; % (as derived earlier)
  Bx = reshape(Bvoxel(1,:), size(CX));
  By = reshape(Bvoxel(2,:), size(CY));
  Bz = reshape(Bvoxel(3,:), size(CZ));
  figure; imagesc(spm_sample_vol(B, Bx, By, Bz, 1));
  % should look like graphics window, with luck...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

> To do a test, I use the template A itself and T_d_A as inputs

Good idea! And indeed, I have tested my above suggestion like this
(with T_r = eye(4) and B==C,) and it seems to do what I want, so
hopefully it will work for you in your more general case. If not, I do
hope this email provides you with the info you need to make progress
(as it's taken a while to write!).

Best of luck,
Ged.