function Create_MaxProb(SPAMsIm, Atlas_FileName, Cod_FileName, Thresh);
%
% Syntax:
% Create_MaxProb(SPAMsIm,Atlas_FileName, Cod_FileName, Thresh);
%
% This function computes the Maximun Probability Atlas based on the
% Statistical Probability Anatomy Maps(SPAMs)files. At each voxel of
% the MaxProb  Atlas, only the label of the most likely structure is
% stored.
%
% Input Parameters:
%  SPAMsIm        : Statistical Probability Anatomy Maps(SPAMs) Files.
%  Atlas_FileName : Maximun Probability Atlas File Name.
%  Cod_FileName   : File with the structure's names and codes(*.cod)
%                   constructed for this Maximun Probability Atlas
%  Thresh         : Threshold for the SPAMs voxel values.
%                   Just the voxels with  higher probability than the specified
%                   threshold are taken into acount to construct the Maximun
%                   Probability Atlas. If the threshold isn't specified then
%                   0 value is taken.
%
% See also: Create_SPAMs
%__________________________________________________
% Authors: Yasser Alemán Gómez & Lester Melie García
% Neuroimaging Department
% Cuban Neuroscience Center
% Last update: November 15th 2005
% Version $1.0

%=====================Checking Input Parameters===========================%
if nargin==0
    [SPAMsIm,sts] = spm_select([1 Inf],'image','Selecting SPAMs Files','',cd);
    [pth,nam,ext] = fileparts(SPAMsIm(1,:));
    [OutputFileName,OutputFilePath] = uiputfile('*.img;*.nii','Saving MaxProb Atlas ...');
    Atlas_FileName = [OutputFilePath OutputFileName];
    Cod_FileName = [Atlas_FileName(1:end-4) '.cod'];
    Thresh = input('Please select a threshold for structure probablity:   ');
else
    if isempty(SPAMsIm)
        [SPAMsIm,sts] = spm_select([1 Inf],'image','Selecting SPMAs Files','',cd);
        [pth,nam,ext] = fileparts(SPAMsIm(1,:));
    end;
    if isempty(Atlas_FileName)
        [OutputFileName,OutputFilePath] = uiputfile('*.img','Saving MaxProb Atlas ...');
        Atlas_FileName = [OutputFilePath OutputFileName];
        Cod_FileName = Atlas_FileName([Atlas_FileName(1:end-4) '.cod']);
    end
    if isempty(Thresh)
        Thresh = input('Please select a threshold for the structure probablity:   ');
    end;
end
%=========================================================================%
%
%=========================Main program=====================================
warning off;
fclose('all');
VSt = spm_vol(SPAMsIm);
NStruct = size(VSt,1);
fid = fopen(Cod_FileName,'w');
[pth,nm,ext] = fileparts(VSt(1).fname);
[pth1,nm1,ext1] = fileparts(Atlas_FileName);
if spm('ver') == 'SPM2'
    dim = [VSt(1).dim(1:3) 2];
    Atlas_FileName = [pth1 filesep nm1 ext1];
    dt = [2 0];
elseif spm('ver') == 'SPM5'
    dim = VSt(1).dim(1:3);
    dt = [2 0];
end
Atlas_FileName = [pth1 filesep nm1 ext];
Vatlas = struct('fname',Atlas_FileName,'dim',dim,'mat',VSt(1).mat,'pinfo',[1 0 0]',...
    'descrip','MaxProbAtlas image','dt',dt);
Np = Vatlas.dim(3);
for i = 1:NStruct
    [pth,name,ext] = fileparts(VSt(i).fname);
    if str2num(name(1:4))
        fprintf(fid,'%s\r',[name(1:4) ' = ' deblank(name(6:end))]);
    else
        fprintf(fid,'%s\r',[num2str(i) ' = ' name]);
    end
end
Vatlas = spm_create_vol(Vatlas);
H = waitbar(0,['Slice  0 of  ' num2str(Np)],'Resize','on','Position',[233.25 237.75 273 50.25],'Resize','off');
for z =1:Np
    waitbar(z/Np,H,['Slice  ' num2str(z) ' of ' num2str(Np)]);
    Temp = zeros(Vatlas.dim(1),Vatlas.dim(2));
    Plane = Temp;
    for i = 1:NStruct
        [pth,name,ext] = fileparts(VSt(i).fname);
        if str2num(name(1:4))
            cod = str2num(name(1:4));
        else
            cod = i;
        end
        Stp = spm_slice_vol(VSt(i),spm_matrix([0 0 z]),VSt(i).dim(1:2),0);
        Res = Temp - Stp;
        ind = find((Res<0)&(Stp>Thresh));
        Temp(ind) = Stp(ind);
        Plane(ind) = cod;

    end
    spm_write_plane(Vatlas,Plane,z);
end
close(H);
fclose all;
%========================End of main program==============================%
return
