function [Vols] = All_Brain_Vol(GMSegFile, WMSegFile, CSSegFile,Images,TempFile, OutputFile, Info);
%
% Syntax :
% [Brain_vol] = All_Brain_Vol(GMSegFile, WMSegFile, CSSegFile, OutputFile, Info);
%
% This function computes the hole brain volume for a group of subjects.
% It can use the files previously obtained in the segmentation step or
% can segment the individual MR images to compute the brain volume. It choose
% one or the other option depending on the value of the variable "Info".
% If Info == y it uses the segmentation files and if Info == n it
% segments the individual images and computes the volume using the segmentation
% results.
%
% Input Parameters:
%  GMSegFile      : Gray Matter Segmentation File.
%  WMSegFile      : White Matter Segmentation File.
%  CSSegFile      : Cerebral Spinal Fluid Segmentation File.
%  SNFiles        : Normalization Transformation File
%  Images         : Individual Images(if dont have the segmentation files).
%  OutputFile     : Path and name for saving the results.
%  Info           : Info contains the answer to this question. Do you have
% %                 the segmentation files? Depending on the value that info
%                   has, this function computes the brain volume for a group
%                   of subjects using the segmentation files or segment the
%                   individual MR images to compute the volumens.
% Output Parameter:
%   Vols : Volume matrix in which are saved the volume for a group of
%               subjects(Ns). The size of Brain_vol is Ns x 1.
%
% See also: spm_segment  Auto_Labelling
%__________________________________________________
% Authors: Yasser Alemán Gómez & Lester Melie Garcia
% Neuroimaging Department
% Cuban Neuroscience Center
% November 15th 2005
% Version $1.0

%=====================Checking Input Parameters===========================%
if ~exist('Info')
    Info = input('Do you have the segmentation files:   ','s');
    Info = lower(deblank(Info));
end
if isempty(Info)
    Info = input('Do you have the segmentation files:   ','s');
    Info = lower(deblank(Info));
end
if ~exist('OutputFile','var')
    [OutputFileName,OutputFilePath] = uiputfile('*.txt','Saving Brain Volumes ...');
    OutputFile = [OutputFilePath OutputFileName];
end
if isempty(OutputFile)
    [OutputFileName,OutputFilePath] = uiputfile('*.txt','Saving Brain Volumes ...');
    OutputFile = [OutputFilePath OutputFileName];
end

%=========================================================================%
%
%=========================Main program=====================================
dir = which('IBASPM');[pth,nm,ext] = fileparts(dir);
VHemi = [pth filesep 'Hemisph_Mask.img']; 
switch Info
    case 'y'
        %=====================Checking Input Parameters===================%
        if ~exist('GMSegFile','var')
            [GMSegFile,sts] = spm_select([1 Inf],'image','Selecting Gray Matter Images','',cd);
        end;
        if isempty(GMSegFile)
            [GMSegFile,sts] = spm_select([1 Inf],'image','Selecting Gray Matter Images','',cd);
        end;
        if ~exist('WMSegFile','var')
            [WMSegFile,sts] = spm_select([1 Inf],'image','Selecting White Matter Images','',cd);
        end;
        if isempty(WMSegFile)
            [WMSegFile,sts] = spm_select([1 Inf],'image','Selecting White Matter Images','',cd);
        end;
        if ~exist('CSSegFile','var')
            [CSSegFile,sts] = spm_select([1 Inf],'image','Selecting Cerebral Spinal Fluid Images','',cd);
        end
        if isempty(CSSegFile)
            [CSSegFile,sts] = spm_select([1 Inf],'image','Selecting Cerebral Spinal Fluid Images','',cd);
        end;
        %=================================================================%
        VSeg = spm_vol(GMSegFile);
        Ns = length(VSeg);
        VWM = spm_vol(WMSegFile);
        VCS = spm_vol(CSSegFile);
        Vols = zeros(Ns,16);
        H = waitbar(0,['Subjects  ' num2str(Ns)],'Resize','on','Position',[233.25 237.75 273 50.25],'Resize','off');
        for i=1:Ns
            waitbar(i/Ns,H, ['Subject  ' num2str(i) ' of ' num2str(Ns)]);
            disp(['Case ---> ' num2str(i)]);
            GmVn = 0;WmVn = 0;CmVn = 0;
            Brain = zeros(VSeg(i).dim(1),VSeg(i).dim(2),VSeg(i).dim(3));
            for z = 1:VSeg(i).dim(3)
                %% tissues
                GM  = spm_slice_vol(VSeg(i),spm_matrix([0 0 z]),VSeg(i).dim(1:2),0);
                GmVn = GmVn + sum(sum(GM));
                WM  = spm_slice_vol(VWM(i),spm_matrix([0 0 z]),VWM(i).dim(1:2),0);
                WmVn = WmVn + sum(sum(WM));
                CSF = spm_slice_vol(VCS(i),spm_matrix([0 0 z]),VCS(i).dim(1:2),0);
                CmVn = CmVn + sum(sum(CSF));
                %% inds
                None = 1 - (GM + WM + CSF);
                indGM = find((GM > WM) & (GM > CSF) & (GM > None));
                indWM = find((WM > GM) & (WM > CSF) & (WM > None));
                indCSF = find((CSF > GM) & (CSF > WM) & (CSF > None));
                %% brain
                if (~isempty(indGM))&(~isempty(indWM))&(~isempty(indCSF))
                    WM(unique([indGM; indWM; indCSF])) = 1;
                    Brain(:,:,z) = WM; 
                end;
            end
            BmVn = length(find(Brain));
            voxsize = abs(diag(VSeg(i).mat(1:3,1:3))');
            Vols(i,1) = BmVn*prod(voxsize)/1000; %Brain Volume
            Vols(i,2) = GmVn*prod(voxsize)/1000; %Gray Matter Volume
            Vols(i,3) = WmVn*prod(voxsize)/1000; %White Matter Volume
            Vols(i,4) = CmVn*prod(voxsize)/1000; %CSF Volume
            Vols(i,5) = Vols(i,2)./Vols(i,1); %GrayM/Brain
            Vols(i,6) = Vols(i,3)./Vols(i,1); %WhiteM/Brain
            Vols(i,7) = Vols(i,4)./Vols(i,1); %CSF/Brain
            Vols(i,8) = Vols(i,2)./Vols(i,3); %GrayM/WhiteM
            [pth,nm,ext] = fileparts(VSeg(i).fname);
            Names{i,1} = nm;
            save([OutputFile(1:end-4) '.mat'],'Vols');
        end
    case 'n'
        spm_defaults
        warning off;
        global defaults
        dseg                   = defaults.segment;
        dseg.write.wrt_cor     = 0;
        dnrm                   = defaults.normalise;
        dnrm.write.vox         = [1 1 1];
        dnrm.write.bb          = reshape([-90 90 -126 90 -72 108],2,3); %bounding box template
        dnrm.write.preserve    = 0;
        dnrm.estimate.graphics = 0;   % Crashes otherwise
        %=====================Checking Input Parameters===================%
        if ~exist('Images','var')
            [Images,sts] = spm_select([1 Inf],'image','Selecting Individual Images','',cd);
        end;
        if isempty(Images)
            [Images,sts] = spm_select([1 Inf],'image','Selecting Individual Images','',cd);
        end;
        if ~exist('TempFile','var')
            [TempFile] = Choosing_TempFile;
        end;
        if isempty(TempFile)
            [TempFile] = Choosing_TempFile;
        end;
        
        %=================================================================%
        V = spm_vol(Images);
        Ns = length(V);
        Vols = zeros(Ns,8);
        H = waitbar(0,['Subjects  ' num2str(Ns)],'Resize','on','Position',[233.25 237.75 273 50.25],'Resize','off');
        for i=1:Ns
            waitbar(i/Ns,H,['Subject  ' num2str(i) ' of ' num2str(Ns)]);
            disp(['Case ---> ' num2str(i)]);
            disp('Segmenting... ');
            V0 =spm_segment(V(i),spm_vol(TempFile),dseg);
            GmVn = 0;WmVn = 0;CmVn = 0;
            Brain = zeros(V(i).dim(1),V(i).dim(2),V(i).dim(3));
            for pp=1:V(i).dim(3),
                GM = double(V0(1).dat(:,:,pp))/255;
                GmVn = GmVn + sum(sum(GM));
                WM = double(V0(2).dat(:,:,pp))/255;
                WmVn = WmVn + sum(sum(WM));
                CSF = double(V0(3).dat(:,:,pp))/255;
                CmVn = CmVn + sum(sum(CSF));
                None = 1 - (GM + WM + CSF);
                indGM = find((GM > WM) & (GM > CSF) & (GM > None));
                indWM = find((WM > GM) & (WM > CSF) & (WM > None));
                indCSF = find((CSF > GM) & (CSF > WM) & (CSF > None));
                %% brain
                if (~isempty(indGM))&(~isempty(indWM))&(~isempty(indCSF))
                    Temp = zeros(size(WM));
                    Temp(unique([indGM; indWM; indCSF])) = 1;
                    Brain(:,:,pp) = Temp; 
                end;
            end;
            BmVn = length(find(Brain));
            voxsize = abs(diag(V(i).mat(1:3,1:3))');
            Vols(i,1) = BmVn*prod(voxsize)/1000; %Brain Volume
            Vols(i,2) = GmVn*prod(voxsize)/1000; %Gray Matter Volume
            Vols(i,3) = WmVn*prod(voxsize)/1000; %White Matter Volume
            Vols(i,4) = CmVn*prod(voxsize)/1000; %CSF Volume
            Vols(i,5) = Vols(i,2)./Vols(i,1); %GrayM/Brain
            Vols(i,6) = Vols(i,3)./Vols(i,1); %WhiteM/Brain
            Vols(i,7) = Vols(i,4)./Vols(i,1); %CSF/Brain
            Vols(i,8) = Vols(i,2)./Vols(i,3); %GrayM/WhiteM
            [pth,nm,ext] = fileparts(V(i).fname);
            Names{i,1} = nm;
            save([OutputFile(1:end-4) '.mat'],'Vols');
        end

end
close(H);
fid = fopen(OutputFile,'w');
fprintf(fid,'%s ', 'Subject         ');fprintf(fid,'%s','BrainVol(cm^3)         ');
fprintf(fid,'%s','GrayMVol(cm^3)         ');fprintf(fid,'%s','WhitMVol(cm^3)         ');
fprintf(fid,'%s','CSFMVol(cm^3)         ');fprintf(fid,'%s','GrayMVol/BrainVol         ');
fprintf(fid,'%s','WhitMVol/BrainVol         ');fprintf(fid,'%s','CSFMVol/BrainVol         ');
fprintf(fid,'%s','GrayMVol/WhitMVol  ');fprintf(fid,'\r');
for i=1:Ns
    fprintf(fid,'%s ',char(Names(i,1)));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,1));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,2));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,3));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,4));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,5));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,6));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,7));
    fprintf(fid,'%s ','          '); fprintf(fid,'%8.4f',Vols(i,8));
    fprintf(fid,'\r');
end;
fclose(fid);
%========================End of main program==============================%
return;

%=====================Internal functions==================================%
function [TempFile] = Choosing_TempFile;
Type = input('Choosing template(T1, T2, PET, SPECT, EPI, PD or Other):   ','s');
Type = lower(Type);
switch Type
    case 't1'
        if exist(fullfile(spm('Dir'),'templates','T1.nii'))
            TempFile = fullfile(spm('Dir'),'templates','T1.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'T1.mnc');
        end
    case 't2'
        if exist(fullfile(spm('Dir'),'templates','T2.nii'))
            TempFile = fullfile(spm('Dir'),'templates','T2.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'T2.mnc');
        end
    case 'pet'
        if exist(fullfile(spm('Dir'),'templates','PET.nii'))
            TempFile = fullfile(spm('Dir'),'templates','PET.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'PET.mnc');
        end
    case 'spect'
        if exist(fullfile(spm('Dir'),'templates','SPECT.nii'))
            TempFile = fullfile(spm('Dir'),'templates','SPECT.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'SPECT.mnc');
        end
    case 'epi'
        if exist(fullfile(spm('Dir'),'templates','EPI.nii'))
            TempFile = fullfile(spm('Dir'),'templates','EPI.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'EPI.mnc');
        end
    case 'pd'
        if exist(fullfile(spm('Dir'),'templates','PD.nii'))
            TempFile = fullfile(spm('Dir'),'templates','PD.nii');
        else
            TempFile = fullfile(spm('Dir') , 'templates' , 'PD.mnc');
        end
    case 'other'
        [TempFileName,TempFilePath] = uigetfile('*.mat;*.img','Reading Reference Template File ...');
        TempFile = [TempFilePath TempFileName];
    otherwise
        disp('Unknown template file ...');
end;
return;


