Dear All:
Adopted from spm_getSPM.m (spm99), I wrote a m-function (wrtTimg.m,
attached file)
to quickly generate t/z score images without using SPM-GUI
for a bunch of (simple) contrasts.
(IT IS FOR PET ONLY)
It would be great if the spm-experts and/or anybody who is interested
can double check to see if there is anything wrong in it.
Thanks,
Kewei Chen
Samaritan PET Center
Phoenix AZ
function wrtTimg(contrasts,t2z);
%adopted from spm_getSPM to generate t/z-score images for given contrasts
%The saved t/z-score images have names like comp1.img for contrasts(1,:), etc.
%contrasts can be eigher:
% a matrix like [1 -1; -1 1] (one raw one contrast)
%or the name of a file (string) from which the contrasts can be read in
% (this file can contain MATLAB-style comments (% first)
%something like:
% 1 -1 %this is contrast 1
% -1 1 %this is contrast 2
%t2z=1: save as z-score; t2z=0: save as t-score
swd=pwd; %please run this in the result directory
%-----------------------------------------------------------------------
xSDM = load(fullfile(swd,'SPM.mat')); %-Contents of SPM.mat (in a struct)
%-Map beta and ResMS images, canonicalising relative pathnames
% (SPM result images stored with relative pathnames, for robustness.)
%-----------------------------------------------------------------------
xSDM.Vbeta = ...
spm_vol([repmat([swd,filesep],length(xSDM.Vbeta),1),char(xSDM.Vbeta)]);
xSDM.VResMS = spm_vol(fullfile(swd,xSDM.VResMS));
%-Get Stats data from SPM.mat
%-----------------------------------------------------------------------
xX = xSDM.xX; %-Design definition structure
XYZ = xSDM.XYZ; %-XYZ coordinates
S = xSDM.S; %-search Volume {voxels}
R = xSDM.R; %-search Volume {resels}
xSDM = rmfield(xSDM,{'xX','XYZ'}); %-Remove (large) duplicate fields
%-Get/Compute mm<->voxel matrices & image dimensions from SPM.mat
%-----------------------------------------------------------------------
M = xSDM.M;
iM = inv(M);
DIM = xSDM.DIM;
%-Build index from XYZ into corresponding Y.mad locations
%-----------------------------------------------------------------------
QQ = zeros(1,S);
if exist(fullfile(swd,'Yidx.mat'),'file') & exist(fullfile(swd,'Y.mad'),'file')
load(fullfile(swd,'Yidx.mat'))
QQ(Yidx) = 1:length(Yidx);
end
%-Compute & store contrast parameters, contrast/ESS images, & SPM images
%=======================================================================
if isstr(contrasts);
contrasts=load(contrasts);
end;
for i = 1:size(contrasts,1);
c=contrasts(i,:);
disp(['generating t/z-score images for contrast ',int2str(i)]);
disp(c);
while length(c)<size(xX.Bcov,1);
c=[c 0];
end;
%fprintf('\t%-32s: %-10s%20s',sprintf('contrast image %2d',i),...
% '(spm_add)','...initialising') %-#
Q = find(abs(c) > 0);
V = xSDM.Vbeta(Q);
for j = 1:length(Q)
V(j).pinfo(1,:) = V(j).pinfo(1,:)*c(Q(j));
end
%-Prepare handle for contrast image
%-----------------------------------------------------------
xCon.Vcon = struct(...
'fname', fullfile(swd,sprintf('phx_con%d.img',i)),...
'dim', [DIM',16],...
'mat', M,...
'pinfo', [1,0,0]',...
'descrip',sprintf('SPM contrast - %d:',i));
%-Write image
%-----------------------------------------------------------
%fprintf('%s%20s',sprintf('\b')*ones(1,20),'...computing')%-#
xCon.Vcon = spm_create_image(xCon.Vcon);
xCon.Vcon.pinfo(1,1) = spm_add(V,xCon.Vcon);
xCon.Vcon = spm_create_image(xCon.Vcon);
%fprintf('%s%30s\n',sprintf('\b')*ones(1,30),sprintf(...
% '...written %s',spm_str_manip(xCon.Vcon.fname,'t')))%-#
%-Write statistic image(s)
%-------------------------------------------------------------------
%---------------------------------------------------------------
Z = spm_sample_vol(xCon.Vcon, XYZ(1,:),XYZ(2,:),XYZ(3,:),0)./...
(sqrt(spm_sample_vol(xSDM.VResMS, XYZ(1,:),XYZ(2,:),XYZ(3,:),0)*...
(c*xX.Bcov*c') ));
if t2z==1;
Z = spm_t2z(Z,xX.erdf);
end;
str = sprintf('[%.2g]',xX.erdf);
%-Write full statistic image
%---------------------------------------------------------------
%fprintf('%s%30s',sprintf('\b')*ones(1,30),'...writing') %-#
xCon.Vspm = struct(...
'fname', fullfile(swd,sprintf('comp%d.img',i)),...
'dim', [DIM',16],...
'mat', M,...
'pinfo', [1,0,0]',...
'descrip',sprintf('SPM{%c_%s} - contrast %d: ',...
'T',str,i));
tmp = zeros(DIM');
tmp(cumprod([1,DIM(1:2)'])*XYZ -sum(cumprod(DIM(1:2)'))) = Z;
xCon.Vspm = spm_write_vol(xCon.Vspm,tmp);
clear tmp Z
%fprintf('%s%30s\n',sprintf('\b')*ones(1,30),sprintf(...
%'...written %s',spm_str_manip(xCon.Vspm.fname,'t'))) %-#
end % (for i = 1:size(contrasts,1);
return;
%to test
for i=1:136;
kw=getimgsli('comp1.img',i,1);%generated by kw_wrtTimg
spmT=getimgsli('spmT_0002.img',i,1);%generated by spm result-user-interface
if kw==spmT;disp([int2str(i),' ok']);
else;disp([int2str(i),' not ok']);end
subplot(1,3,1);imagek(kw);subplot(1,3,2);imagek(spmT);subplot(1,3,3);
linrfit(kw(:),spmT(:));pause(.1);end
|