Hi Mike,
At 02:04 PM 5/11/2005 -0400, Mike Yassa wrote:
>I am using one of John's scripts to extract the global max t value in a
>statistical t image. However, I would also like to get the coordinates
>associated with this value (or the closest coordintes). Is there any way
>I can do this in SPM or using a Matlab script?
I modified the script to record the coordinates as well as the maximum
value. I tested it t-images I have and seemed to work fine. Attached is the
modified script. Try and see if it works on your data.
Thanks!
-Satoru
Satoru Hayasaka ==============================================
Post-Doctoral Fellow, MR Unit, UCSF / VA Medical Center
Email: shayasak_at_itsa_dot_ucsf_dot_edu Phone:(415) 221-4810 x4237
Homepage: http://www.sph.umich.edu/~hayasaka
==============================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select one .img file
P = spm_get(1,'*.img');
% Read header/mat file information about the image
V = spm_vol(P);
% Get the transformation matrix from voxel space to mm space
MAT = V.mat;
% Initialise maximum value to the most negative possible
mx = -Inf;
% Initializing the coordinates for the max location
mxXYZvox = -Inf*ones(3,1);
mxXYZmm = -Inf*ones(3,1);
% A very small number to be used later
eps = 0.00001;
% Getting XY coordinate mesh grid (in voxel space)
[y x] = meshgrid([1:V.dim(2)],[1:V.dim(1)]');
% Size of each plane
PlSize = V.dim(1)*V.dim(2);
% loop over planes (dim(1:3) is the size of the image)
for i=1:V.dim(3),
% The dimensions of the output image. Co-ordinates
% in this image are x = 1..d(1) and y = 1..d(2)
d = V.dim(1:2);
% An affine transformation mapping from x,y to some
% position in the image, such that the co-ordinate
% x,y in the output image has the following
% co-ordinate in the 3D volume:
% pos = M(1:3,:)*[x y 0 1]';
M = spm_matrix([0 0 i]);
% Sample the volume V, according to the transformation
% matrix M, to create a d(1) by d(2) image using
% interpolation method 0 (nearest neighbour).
img = spm_slice_vol(V,M,d,0);
% Maximum value in this plane is:
mxp = max(img(:));
% Finding the location of the max in this plane
mxpLoc = find(img>(mxp-eps));
if length(mxpLoc)>1, mxpLoc = mxpLoc(length(mxpLoc)); end;
% The max coordinates in this plane (voxel space)
mxpXYZvox = [x(rem(mxpLoc-1,PlSize)+1);
y(rem(mxpLoc-1,PlSize)+1); i];
% The max coordinates in this plane (mm space)
mxpXYZmm = MAT*[mxpXYZvox; 1];
mxpXYZmm(4) = [];
% Maximum value encountered so far is:
if mxp>mx
mx = mxp;
mxXYZvox = mxpXYZvox;
mxXYZmm = mxpXYZmm;
end;
end;
% Display the maximum value
fprintf('Maximum value is %g\n', mx);
% Displaying the maximum coordinates (voxel space)
mxXYZvox
% Displaying the maximum coordinates (mm space)
mxXYZmm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|