function glassbrain3G99(varargin)
%glassbrain3.m
%-------------------------------
%Function to dump voxels in the glass brain into a txt file to be 
%read into the Tal daemon.
%FORMAT: glassbrain3(SPM, VOL, -f)
%SPM and VOL are necessary; they must be passed in from the SPM
%workspace and should be the structures from spm_results of
%the same name.
%The -f flag is optional and indicates that glassbrain is to 
%dump all activated voxels from the SPM to text; otherwise, the
%default is given in vox (100 for now).
%Sue Whitfield & Jeff Cooper
%Sept. 2002


if nargin < 2
    error('Insufficient arguments.  Format: glassbrain(SPM, VOL, (''-f'')).');
end
S = 0;
if nargin > 2
    if (ischar(varargin{3}))
        if (varargin{3} == '-f')
            S = 1;
        end
    end
    if S == 0
        warning('Argument string not understood.  Use glassbrain(SPM,VOL,''-f'') for full (all-voxel) output.');
    end
end

%Variables----------------------------------------------------
SPM = varargin{1};
VOL = varargin{2};
R         = VOL.R;
FWHM      = VOL.FWHM;
v2r       = 1/prod(FWHM);			%-voxels to resels
n         = SPM.n;
STAT      = SPM.STAT;
df        = SPM.df;
u         = SPM.u;
if (S)
    vox = inf;
else
    vox = 100;   %default number of voxels per cluster to include
end
voxlim = 10000; %max number of voxels per file

%-Tolerance for p-value underflow, when computing equivalent Z's
tol = eps*10;




%Assemble voxel data structure--------------------

D(:,1) = SPM.Z';
D(:,2:4) = SPM.XYZmm';  
%Note that indices for SPM.Z and SPM.XYZ reference same voxel - thus they do 
%here as well
D(:,5:7) = mni2tal(SPM.XYZmm');
D(:,8:10) = SPM.XYZ';
D(:,11) = (spm_clusters(SPM.XYZ))';

%D is now a double array of n x 11 size, where n is the number of voxels in
%the SPM.  Its first column is the voxel Z-score; its 11th is the voxel cluster
%number.  The columns in between are coordinates for that voxel; MNI, Tal, and 
%voxel coordinates, respectively, each a three-column set.

D = sortrows(D,11);
%D is now sorted into contiguous cluster rows (which it should have been anyways,
%but it's nice to guarantee it); we can extract each cluster
%and treat it as a sub-array of voxels.


p = max(D(:,11));
%p is the number of clusters.

for i = 1:p
    k = find(D(:,11)==i);
    E(i) = {D(k,:)};
    %Each E(i) now is a double array consisting of all and only the voxels in 
    %cluster i.  This will do us nicely.
    Q = sortrows(E{i},1);
    Q = flipud(Q);
    E{i} = Q;
    r(i) = Q(1,1);
end

%E is now a cell array in which each cell is a sorted list of voxels in a given
%cluster, in descending order of Z-score.  That's nice.  We now need to get
%E in order of cluster rank, which is easy enough.  SPM ranks clusters in order
%of lowest single uncorrected probability in that cluster.  Conveniently, r now 
%indexes those peak Z-scores.

r(2,:) = double(1:p);
r = r';
r = sortrows(r,1);
r = flipud(r);

tic
%File creation and writing---------------------------------------

outname = spm_str_manip(spm_input('Select TAL coord output filename ', 1, 's'),'sdv');
resdir = spm_get(-1,'.','Select directory to save output file.');close;
cd(resdir)

b=size(D,1);
m=b/voxlim;
m=ceil(m);
%this stuff handled by spm_str_manip
%q=findstr(outname,' ');
%for i=1:length(q)
%   outname(q(i))='_';
%end
%q=findstr(outname,':');
%for i=1:length(q)
%   outname(q(i))='_';
%end

c = 1;
outfile=cat(2,outname,'_',num2str(c),'.txt')
fid=fopen(outfile,'w');

% Print out the SPM VOLUME Parameters
% ~~~~~~~~~~~~~~~~~~Appended by Yuqing Gong~~~~Start~~~~~~~~~~
%-Volume, resels and smoothness 
%-----------------------------------------------------------------------
R         = VOL.R;
FWHM      = VOL.FWHM;
v2r       = 1/prod(FWHM);			%-voxels to resels
n         = SPM.n;
STAT      = SPM.STAT;
df        = SPM.df;
u         = SPM.u;
k         = SPM.k;
S         = VOL.S; % Appended

FWHMmm          = FWHM.*VOL.VOX'; 			% FWHM {mm}
Pz              = spm_P(1,0,u,df,STAT,1,n,S);
Pu              = spm_P(1,0,u,df,STAT,R,n,S);
[P Pn Em En EN] = spm_P(1,k,u,df,STAT,R,n,S);

%-Footnote with SPM parameters
%-----------------------------------------------------------------------
fprintf(fid,'Height threshold: %c = %0.2f, p = %0.3f (%0.3f corrected)\n',...
    STAT,u,Pz,Pu);
fprintf(fid,'Extent threshold: k = %0.0f voxels, p = %0.3f (%0.3f corrected)\n',...
    k/v2r,Pn,P);
fprintf(fid,'Expected voxels per cluster, <k> = %0.3f\n',En/v2r);
fprintf(fid,'Expected number of clusters, <c> = %0.2f\n',Em*Pn);
fprintf(fid,'Degrees of freedom = [%0.1f, %0.1f]\n',df);
fprintf(fid,['Smoothness FWHM = %0.1f %0.1f %0.1f {mm} ',...
        ' = %0.1f %0.1f %0.1f {voxels}\n'],FWHMmm,FWHM);
fprintf(fid,'Search volume: S = %0.0f mm^3 = %0.0f voxels = %0.1f resels\n',...
    S*prod(VOL.VOX),S,R(end));
fprintf(fid,['Voxel size: [%0.1f, %0.1f, %0.1f] mm ',...
        ' (1 resel = %0.2f voxels)\n\n\n'],...
    VOL.VOX,prod(FWHM));
    
% ~~~~~~~~~~~~~~~~~~~~Appended by Yuqing Gong~~~~End~~~~~~~~~~

a = 0;
for i = 1:p
    curr = E{r(i,2)};
    class_brain = spm_mni2class2( curr(:,2:4),'f' );
    probs = [];
    for j = 1:size(curr,1)
        if j > vox
            break;
        end
        if j==1
            %compute and print cluster-level stats
            Nr = (size(curr,1)).*v2r; %gives us extent in resels for spm_P
            [Pk Pn] = spm_P(1,Nr,u,df,STAT,R,n);	% [un]corrected {based on k)
            fprintf(fid,'%1.3f\t%6.0f\t%1.3f\t',Pk,size(curr,1),Pn);
        else
            fprintf(fid,' \t  \t  \t');
        end
        probs(j,1) = spm_P(1,0,curr(j,1),df,STAT,1,n); %Pz from spm_list
        probs(j,2) = spm_P(1,0,curr(j,1),df,STAT,R,n); %Pu from spm_list
        if probs(j,1)<tol 
            probs(j,3)=Inf;
        else 
            probs(j,3) = spm_invNcdf(1 - probs(j,1)); %Ze from spm_list
        end
        con = '%1.3f\t%3.2f\t%3.2f\t%1.3f\t%2.0f\t%2.0f\t%2.0f\t';
        fprintf(fid,con,probs(j,2),curr(j,1),probs(j,3),probs(j,1),curr(j,2:4));
        fprintf(fid,'%2.0f\t%2.0f\t%2.0f',curr(j,5:7));
        % ~~~~~~~~~~~~~~~~~~Appended by Yuqing Gong~~~~Start~~~~~~~~~~
       
        for index_cat = 1:5,
            fprintf(   fid,'\t\t%s',cell2mat( class_brain{j,index_cat} )   );
        end
        fprintf(fid,'\n');
        % ~~~~~~~~~~~~~~~~~~~~Appended by Yuqing Gong~~~~End~~~~~~~~~~
        a = a+1;
        if a > voxlim
            fclose(fid);
            c = c+1;
            outfile = cat(2,outname,'_',num2str(c),'.txt')
            fid = fopen(outfile,'w');
            a = 0;
        end
    end
end

status=fclose(fid);

%clear n P q con D E r probs curr tol vox S a b m c
%clear i j k R u df Nr Pk Pn STAT outfile Q p v2R FWHM fid


toc

