function vol = Struct_Vol(AtlasFile,CodeFile,voxsize);
%
% Syntax: 
% vol = struct_vol(AtlasFile,CodeFile,voxsize);
%
% This script computes the volume for different structures based on individual 
% Atlas previously obtained by the atlasing process.
%
% Input Parameters:
%    AtlasFile : Individual Atlas 
%    CodeFile  : File with names and codes of different structures (*.cod) 
%    voxsize   : Voxel size (in millimeter [mm]); If vox is not defined 
%                then the voxel size of the first atlas image is used.
%
% Output Parameters:
%    vol : volume of structures in cm^3.
%          Struct with fields:
%             vol(i).name : structure name
%             vol(i).code : structure code
%             vol(i).vol  : volume of structure (in cm^3 or cc)
%             vol(i).nvox : number of voxels
%
% See also: vols_stats
%__________________________________________________________________________
% Authors: Lester Melie García & Yasser Alemán Gómez
% Neuroimaging Department
% Cuban Neuroscience Center
% Last update: November 15th 2005
% Version $1.0

%=====================Checking Input Parameters===========================%
if nargin<1
   [AtlasFileName,AtlasFilePath] = uigetfile({'*.img';'*.mat'},'Reading Individual Atlas File ...');
   AtlasFile = [AtlasFilePath AtlasFileName];
   disp(['Individual Atlas File : ' AtlasFileName]);
end;
[pth,nam,ext] = fileparts(AtlasFile);
if ext=='.img'
    warning off;
    Va = spm_vol(AtlasFile);
    Atlas = spm_read_vols(Va);
    warning on;
else % for *.mat files
    Atlas = load(AtlasFile);
    Atlas = Atlas.Atlas;
end;
if ~exist('CodeFile','var')
   [CodeFileName,CodeFilePath] = uigetfile('*.cod','Reading Strutures Code File ...');
   CodeFile = [CodeFilePath CodeFileName];
   disp(['Structures Code File : ' CodeFileName]);
end;
if isempty(CodeFile)
   [CodeFileName,CodeFilePath] = uigetfile('*.cod','Reading Strutures Code File ...');
   CodeFile = [CodeFilePath CodeFileName];
   disp(['Structures Code File : ' CodeFileName]);
end;
if ~exist('voxsize','var')
   voxsize = abs(diag(Va.mat(1:3,1:3))'); 
end;
% ========================================================================%
%
%=========================Main program====================================%    
voxvol = prod(voxsize);
fid = fopen(CodeFile,'r');
i=0;
while ~feof(fid)
   i = i+1;
   jline = fgetl(fid);
   sep = findstr(jline,'=');
   StructCode = str2num(jline(1:sep-1));
   StructName = jline(sep+1:end);
   ind = find(Atlas==StructCode);
   Nvox = length(ind);
   vol(i).name = StructName;
   vol(i).code = StructCode;
   vol(i).vol  = Nvox*voxvol./1000; %(in cm^3)
   vol(i).nvox = Nvox;
end;    
%========================End of main program==============================%
return;


    