% read an image files (e.g. an filtered or imcalced  image)
% return coordinates in voxels and in mm

hold=1;

P = spm_get(1, '.img', 'Select image files.');
V= spm_vol(P);

% read volume and return mm and intensities, without masking
[Y,XYZmm] = spm_read_vols(V,0);

% remove the masked out voxels
XYZmm = XYZmm(:,find(Y));


% convert mm coordinates into voxel coordinates
XYZ   = (  ([[1 0 0; 0 1 0; 0 0 1] V.mat(1:3,4)]...
           *[XYZmm;-ones(1,size(XYZmm,2))])'...
              / (V.mat(1:3,1:3))...
        )';

% convert voxel coordinate into mm coordinates
% XYZmm = V.mat(1:3,:) * [XYZ;    ones(1,size(XYZmm,2))];

% offer an option to convert from MNI to Talairach Atlas with mni2tal
% -------------------------------------------------------------------

fprintf('\nWould you like to convert the coordinates from MNI to Talairach');
fprintf('\nspace using mni2tal? ');

convert=input('Type [y] to convert, return otherwise: ','s');

if strcmp('y',convert); XYZmm=mni2tal(XYZmm); end;


% formulate the clusters
%----------------------
AA=spm_clusters(XYZ);
num_clusters=max(AA);


% Headers for text table...
%----------------------------
TabDat.tit = ['Cluster information for', sprintf('\n'),sprintf(V.fname)];
TabDat.hdr = {	'cluster',	'   k';...
		'centroid',	'x,y,z {mm}'}';

TabDat.fmt = {	'%5.0f', ...			%-Cluster
		'%10.0f', '%3.0f', '%3.0f'};				%-XYZ
TabLin = 1;


% Calculate and save centroids and number of voxels in each cluster
% -----------------------------------------------------------------

for cluster = 1:num_clusters

      % find mm coordinate for this voxels
      TMPmm= XYZmm(:,find(AA==cluster));
      
      %  fprintf('%3.0f\t %3.0f\t  %3.0f\t\n',round(TMPmm))

      centroid=round(mean(TMPmm,2));
      num_vox =length( find(AA==cluster));

      %-Update table
      %-------------
      [TabDat.dat{TabLin,1:4}] = deal(num_vox, centroid(1), ...
                                      centroid(2), centroid(3));
      TabLin = TabLin + 1;
 
      clear TMPmm;

end


% print table information and save as file
% -------------------------------------------------------

      fprintf('\n\n %s\n\n',TabDat.tit)

      fprintf('%s\t\t',TabDat.hdr{1,1}), fprintf('%s\n',TabDat.hdr{1,2});
      fprintf('%s\t\t',TabDat.hdr{2,1}), fprintf('%s\n',TabDat.hdr{2,2});

      
      for i = 1:size(TabDat.dat,1)
	for j=1:size(TabDat.dat,2)
		fprintf(TabDat.fmt{j},TabDat.dat{i,j})
		fprintf('\t')
	end
	fprintf('\n')
      end


      % write a .dat text file with list of clusters info

      if strcmp('y',convert);      name_dat=sprintf('%s',V.fname,'_tal');
         else                      name_dat=sprintf('%s',V.fname,'_mni');
      end;

      %if name_dat(:,length(name_dat)-3:length(name_dat) )=='.img'
      %     name_dat=sprintf('%s',...
      %              name_dat(:,1:length(name_dat)-4),'.dat');
      %else 
      name_dat=sprintf('%s',name_dat,'.dat');
      %end;


      fid=fopen(name_dat,'w');

      for i = 1:size(TabDat.dat,1)
	for j=1:size(TabDat.dat,2)
		fprintf(fid, TabDat.fmt{j},TabDat.dat{i,j});
		fprintf(fid,'\t');
	end;
ls        fprintf(fid,'\n');
      end

      fclose(fid);

fprintf('\n\nThe above information has also been save as a text file in\n');
fprintf(name_dat); fprintf('\n\n');

clear
