function Create_SPAMs(AImages, SNFiles, AtlasFile, Output_dir, StrL);
%
% Syntax: 
% Create_SPAMs(AImages, SNFiles, AtlasFile, Output_dir);
%
% This function constructs statistical probability anatomy maps(SPAMs).
% It can be very useful because of the advantage to use the automatic parcellation,
% to build this maps. In this case you can segment automatically a dataset 
% of MRI images and then construct the SPAMs files for your database 
% saving time from the manual segmentation process, just using the results 
% obtained by the automatic labelling for individual subjects.It does not
% implies that this script can't be used to create SPAMs from a manual
% segmented dataset. It can be done too, just knowing the
% normalisation transformation files, a file with the name of the structures 
% manually labelled and a matrix containing the output SPAMs dimensions,
% voxel size and the center of the coordinate system(some changes into the 
% code must be done).
% For a group of subjects previously automatically labelled, the dimension of 
% the output SPAMs files are the same than the dimensions of the atlas files
% employed to parcelate the individual MRIs.  
% The SPAMs files are created using the same method that the employed
% by Evans et al. 1994.  
%
% Input Parameters:
%   AImages    : Individual Atlased Files
%   SNFiles    : Individual Normalisation Transform Matrix
%   AtlasFile  : Reference Atlas File
%   Output_dir : Output Directory to save the SPAMs Files
%   StrL       : Structures code list. It cotains a list with the structures codes 
%    - 'all'        : Computes the SPAMs for all the structures in the CodeFile.
% i.e - '1,78-92,6' : Computes the SPAMs for the following structures. Structure 
%                  with codes iqual to 1, structures with codes from 78 to 92 and 
%                  structure with code iqual to 6.
%                  It offers the chance to choose 2 or more structures with
%                  codes that aren't one after the other, just separating the 
%                  structures codes with a coma(,).
%                  If you want to select a range of codes, you have to put the first 
%                  structure code, then '-' and next the last structure code, 
%                  i.e, 78-92. This expresion computes the SPAMs for all the
%                  structures with codes from 78 to 92.  
% Note: A correct structure selection is very IMPORTANT!!! to compute the SPAMs. 
% 1.- If you want to compute the SPAMs files for individual structures, you must 
% separate the structures codes with a coma (ie. 34,67,89,90).
% 2.- If you want to compute the SPAMs for a range of structures, you must
%  separate the first and the last code with a '-'(i.e, 67-90).
% 3.- Combining the point 1. and 2 the SPAMs for individual and a range of
%  can be obtained (i.e 5,7-56,6,67-78).  
%
% Related reference:
% 1.- Evans AC, Collins DL, Neelin P, MacDonald D, Kamber M, Marrett TS 
%    (1994). Three-Dimensional Correlative Imaging: Applications in Human 
%    Brain Mapping, in: Functional Neuroimaging: Technical Foundations, 
%    Thatcher RW, Hallett M,Zeffiro T, John ER, Huerta M [eds.], 145-162.    
% 2.- J. Mazziotta, A. Toga, A. Evans, P. Fox, and J. Lancaster.
%     A probabilistic atlas of the human brain: theory and rationale for its 
%     development. the international consortium for brain mapping.
%     NeuroImage, 2(2):89-101, 1995.
%
% See also: Create_MaxProb 
%__________________________________________________________________________
% Authors: Yasser Alemán Gómez & Lester Melie García  
% Neuroimaging Department
% Cuban Neuroscience Center
% Last update: November 15th 2005
% Version $1.0

warning off
fclose('all');
%=====================Checking Input Parameters===========================%
if nargin==0
    [AImages,sts] = spm_select([1 Inf],'image','Selecting Atlased Images','',cd);
    [pth,nam,ext] = fileparts(AImages(1,:));
    [SNFiles,sts] = spm_select([1 Inf],'mat','Selecting Normalization Transform Files','',cd);
    [AtlasFileName,AtlasFilePath] = uigetfile({'*.img'},'Reading Reference Atlas File ...');
    AtlasFile.name = [AtlasFilePath AtlasFileName];
    [CodeFileName,CodeFilePath] = uigetfile('*.cod','Reading Strutures Code File ...');
    AtlasFile.cod = [CodeFilePath CodeFileName];
    Output_dir = uigetdir(cd ,'Selecting Output Directory');
    mkdir(Output_dir, 'SPAMs');
    StrL = input('Please enter the structures codes:  ','s');
else
    if isempty(AImages)
        [AImages,sts] = spm_select([1 Inf],'image','Selecting Gray Matter Images','',cd);
        [pth,nam,ext] = fileparts(AImages(1,:));
    end;
    if isempty(StrL)
        StrL = input('Please enter the structures codes:  ','s');
    end;
    if isempty(AtlasFile)
        [AtlasFileName,AtlasFilePath] = uigetfile({'*.mat;*.img'},'Reading Reference Atlas File ...');
        AtlasFile.name = [AtlasFilePath AtlasFileName];
        [CodeFileName,CodeFilePath] = uigetfile('*.cod','Reading Strutures Code File ...');
        AtlasFile.cod = [CodeFilePath CodeFileName];
    end;
    if isempty(SNFiles)
        [SNFiles,sts] = spm_select([1 Inf],'mat','Selecting Normalization Transform Files','',cd);
    end;
    if isempty(Output_dir)
        Output_dir = uigetdir(cd ,'Selecting Output Directory');
        mkdir(Output_dir, 'SPAMs');
    end
end
%=========================================================================%
%
%=========================Main program=====================================
CodeFile = AtlasFile.cod;
AtlasFile = AtlasFile.name;
V = spm_vol(AImages);
[pth,nam,Aext] = fileparts(AtlasFile);
switch Aext
    case {'.img'}
        Vatlas   = spm_vol(AtlasFile);
        VAtlas_mat = Vatlas.mat;
        AtlasDim =[Vatlas.dim(1:3)];
    case {'.nii'}
        Vatlas   = spm_vol(AtlasFile);
        VAtlas_mat = Vatlas.mat;
        AtlasDim =[Vatlas.dim(1:3)];
    otherwise
        disp('Unknown format file ...');
        return;
end;
fid = fopen(CodeFile);
cont = 0;
while 1
    cont = cont + 1; 
    line = fgetl(fid);
    if ~ischar(line),   break,   end
    ind = strfind(line,'=');
    StructNames{1, cont} = line(ind+1:end);
    StructCodes(1,cont) = str2num(deblank(line(1:ind-1)));
end
fclose(fid);
clear tline cont;
NSubj = length(V);
if strcmp(StrL,'All')
    StrList = StructCodes;
else
  StrList = Sel_Structs(StrL,StructCodes);
end
NStruct = length(StrList);
[pth,nm,ext] = fileparts(V(1).fname);
if spm('ver') =='SPM2'
    dim = [AtlasDim 2];
elseif spm('ver') =='SPM5'
    dim = AtlasDim;
end
Vol = struct('fname','','dim',dim,'mat',VAtlas_mat,'pinfo',[1/255 0 0]',...
    'descrip','SPAM image','dt',[2 0]);
Vol = Vol(ones(NStruct,1));
for j = 1:NStruct
    ind = find(StructCodes == StrList(j));
    Struct = StructCodes(1,ind);
    Vol(j).fname = sprintf('%s%sSPAMs%s%.4d_%s',...
        deblank(Output_dir),filesep,filesep,Struct,deblank(StructNames{1,ind}),ext);
end
Vol = spm_create_vol(Vol);
tic
H = waitbar(0,['Slice  0 of  ' num2str(AtlasDim(1,3))],'Resize','on','Position',[233.25 237.75 273 50.25],'Resize','off');
for k = 1:AtlasDim(1,3)
    waitbar(k/AtlasDim(1,3),H,['Slice  ' num2str(k) ' of ' num2str(AtlasDim(1,3))]);
    B = zeros([AtlasDim(1,1) AtlasDim(1,2) NStruct]);
    for i = 1:NSubj
        load('-mat',SNFiles(i,:));
        M = V(i).mat\VF.mat*Affine/VG.mat*VAtlas_mat;
        A = spm_slice_vol(V(i),M*spm_matrix([0 0 k]),AtlasDim(1:2),0);
        for j = 1:NStruct
            ind = find(StructCodes == StrList(j));
            B(:,:,j) = B(:,:,j) + (A == StructCodes(1,ind))/NSubj;
        end
    end
    for j = 1:NStruct
        spm_write_plane(Vol(j),B(:,:,j),k);
    end
end
close(H);
toc
fclose all;
%========================End of main program==============================%
return;

%=====================Internal functions==================================%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T2 = get_2Dtrans(T3,B,j)
d   = [size(T3) 1 1 1];
tmp = reshape(T3,d(1)*d(2),d(3));
T2  = reshape(tmp*B(j,:)',d(1),d(2));
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function StrList = Sel_Structs(StrL,StructCodes);
StrL = lower(StrL);
if ~strcmp(StrL,'all')
    Del = strread(StrL,'%s','delimiter',',');
    StrList = 0;
    for i = 1:length(Del)
        Temp = strread(Del{i,1},'%s','delimiter','-');
        if length(Temp)==2
            StrList = [StrList eval([Temp{1,1} ':' Temp{2,1} ])];
        else
            StrList = [StrList str2num(cell2mat(Temp))];
        end
    end
    StrList(1) =[];
else
    StrList = 1:length(StructCodes);
end