Dear list,
Hello. I have attached two programs that were of
interest and use to me and may be of interest to
others. What the first one does is create a binary
mask image for each cluster in an SPM2 result. The
second program asks for a series of mask images and T
map images. It calculates the mean, variance and
number of voxels for the T maps that are in the masks.
The results are saved to a comma separated file.
The create_cluster_masks is not perfect for it does
not label the masks based on where in the brain they
are. They are not even in the order of most to least
significance. They are in the order that the program
spm_clusters puts them in. I believe this is bottom to
top (in terms of z slice number).
I feel that I have documented them well and if any one
makes some modifications or improvements I would like
to hear about the changes.
best regards,
Jason Steffener
%% create_cluster_masks.m
%% Makes a separate mask image for each cluster in an SPM
%% results. This program is run by typing its name at the
%% MatLab command prompt after running the SPM2 results.
%% (Notes: Only tested on SPM2 and be careful because if
%% you have a lot of clusters you will end up with a lot
%% of output images.)
%%
%% What it does:
%% Take all voxels above threshold and displayed in your results
%% (xSPM.XYZ) and find out which cluster each one belongs to. The result of
%% this is listed in the vector A.
%% This just writes out clusters in the order of bottom to top (I think) it
%% would be better to have them labeled based on their significance. That
%% is the next step.
%% But, how to do it?...
%%
%% written by Jason Steffener, 3/30/06
%%
%% xSPM.XYZ is a variable in the MatLab workspace after results are calculated.
%% spm_clusters performs the clustering analysis of results
A=spm_clusters(xSPM.XYZ);
%% how many clusters are there?
num_cluster=max(A);
%% how many significant voxels are there?
num_voxels=length(xSPM.XYZ);
%% load up the mask image
Vmask=spm_vol(fullfile(pwd,'mask.img'));
%% set an output mask to have the same attributes as the mask image.
Vout=Vmask;
%% cycle over each cluster
for i=1:num_cluster
%% Create the image file for each cluster mask
Vout.fname = fullfile(pwd,['cl_' num2str(i) '_mask.img']);
Vout.descrip = ['cluster ' num2str(i) ' mask'];
Vout = spm_create_vol(Vout);
dim = Vout.dim(1:3);
%% create a volume of zeros
img=uint8(zeros(dim(1:3)));
%% fill in the img matrix
for j=1:num_voxels
if A(j)==i
img(xSPM.XYZ(1,j),xSPM.XYZ(2,j),xSPM.XYZ(3,j))=1;
end
end
%% write the mask to a file
Vout=spm_write_vol(Vout, img);
display(['Finished writing ' Vout.fname]);
end
%% This program is run by typing its name at the MatLab
%% command prompt. It asks you for 1 or more mask images
%% (possibly created with 'create_cluster_masks.m')
%% and then for 1 or more spm_Tmap images. It loads each
%% one up and takes the union of them. So it just
%% multiplies the Tmap by the mask (i1.*i2 in spm_imcalc language).
%% It then calculates
%% the mean, variance and number of voxels in the masked
%% T map (other options are can be easily added). The
%% results are saved to a comma separated file. The
%% output includes the input mask names and the input T
%% map names and paths.
%%
%% written by Jason Steffener 4/4/06
%% Select the mask images
Pmask=spm_get(Inf,'*.img','Select cluster mask(s)');
%% Select you T maps to apply these to
Ptmaps=spm_get(Inf,'spmT*.img','Select T map(s)');
%% How many masks selected
Nmasks=size(Pmask,1);
%% How many T maps selected
Nmaps=size(Ptmaps,1);
%% THe program spends a lot of time loading files in, but this is so it
%% does not overflow your memory by loading them all in.
for i=1:Nmasks
%% load each mask image in once
temp_M=spm_read_vols(spm_vol(Pmask(i,:)));
for j=1:Nmaps
%% load each T map in
temp_T=spm_read_vols(spm_vol(Ptmaps(j,:)));
%% take the union of the mask and T map
temp_TM=temp_T.*temp_M;
%% calculate the mean, variance and number of voxels
M(i,j)=mean(temp_TM(find(temp_TM)));
V(i,j)=var(temp_TM(find(temp_TM)));
N(i,j)=length(temp_TM(find(temp_TM)));
end
end
%% write out data
%% create the output file
fid=fopen(['data_' date '.csv'],'w');
%% first write header
%% first line contains the names of each mask image
fprintf(fid,',,');
for i=1:Nmasks
[temp_path, temp_name]=fileparts(Pmask(i,:));
fprintf(fid,temp_name);
%% The extra commas account for the fact that each mask image will have
%% 3 columns, one for mean, variance and number of voxels
%% So if other measures are included the number of commas should be
%% modified.
fprintf(fid,',,,');
end
%% end the first line
fprintf(fid,'\n');
%% write the second line first and second columns
fprintf(fid,'Path,File,');
%% repeated write the following
for i=1:Nmasks
fprintf(fid,'Mean,');
fprintf(fid,'Variance,');
fprintf(fid,'N voxels,');
end
%% end the second line
fprintf(fid,'\n');
for j=1:Nmaps
%% write the T map path name after some manipulation
[temp_path, temp_name]=fileparts(Ptmaps(j,:));
%% this \ to / conversion is doen so that fprintf does not get confused
%% by the \ in the path name.
temp_path(findstr(temp_path,'\'))='/';
fprintf(fid,temp_path);
fprintf(fid,',');
%% write the T map file name
fprintf(fid,temp_name);
fprintf(fid,',');
%% write the measures
for i=1:Nmasks
fprintf(fid,num2str(M(i,j)));
fprintf(fid,',');
fprintf(fid,num2str(V(i,j)));
fprintf(fid,',');
fprintf(fid,num2str(N(i,j)));
fprintf(fid,',');
end
%% end the line
fprintf(fid,'\n');
end
%% close the file
fclose(fid);
display(['Data written to data_' date '.csv']);
|