function [mesh, info] = gtVTKMeshReader(filename, varargin)
% GTVTKMESHREADER  Read mesh (polydata) from a VTK file.
%
%   mesh = gtVTKMeshReader(filename, varargin)
%   --------------------------------------------------------------------------
%
%     INPUT:
%       filename = <string>  Path to the input VTK file
%
%     OUTPUT:
%       mesh     = <struct>  Mesh in a structure, composed of following fields:
%                            - vertices: Nv x 3 array -> vertices coordinates
%                            - faces:    Nf x 3 array -> face vertices indices
%                            - edges:    Ne x 2 array -> edge vertices indices
%
%     Version 001 20-10-2012 by YGuilhem yoann.guilhem@esrf.fr

%% Guessing output directory, filename and extension
[fpath, fname, fext] = fileparts(filename);
if isempty(fext)
    fext = '.vtk';
end
% Trick to handle ~ path in unix
if ~isempty(fpath) && isunix && fpath(1) == '~'
    fpath(1) = '';
    fpath = fullfile(getenv('HOME'), fpath);
end
filename = fullfile(fpath, [fname fext]);

fid = fopen(filename, 'r');
if fid == -1
    gtError('gtVTKMeshReader:bad_input_file', 'Can''t open the file.');
end

%% Read File description and version
str = fgets(fid);
if ~strcmp(str(1:22), '# vtk DataFile Version')
    gtError('gtVTKMeshReader:bad_input_file', 'Not a valid VTK file.');
end
info.version = sscanf(str(24:end), '%f');

%% Read header
str = fgets(fid);
info.header  = strtrim(str);

%% Read mode
str = fgets(fid);
info.mode = sscanf(str, '%s', 1);
if strcmp(info.mode, 'BINARY')
    binary = true;
else
    binary = false;
end

%% Read dataset type
str = fgets(fid);
info.type = sscanf(str, 'DATASET %s', 1);
if ~strcmp(info.type, 'POLYDATA')
    gtError('gtVTKMeshReader:wrong_content', ...
        'This vtk file does not contain PolyData dataset!');
end

%% Read number of vertices and coordinate datatype
str = fgets(fid);
nVert = sscanf(str, 'POINTS %d', 1);
str2 = ['POINTS ' num2str(nVert) ' %s'];
coordType = sscanf(str, str2);
if binary
    switch coordType
    case {'float', 'double'}
        coordFMT = coordType;
    otherwise
        gtError('gtVTKMeshReader:wrong_coord_type', ...
            ['cannot handle coordinate with type: ' coordType]);
    end
else
    switch coordType
    case 'float'
        coordFMT = '%f %f %f';
    case 'double'
        coordFMT = '%f %f %f';
    otherwise
        true;
    end
end
info.Nvertices = nVert;

%% Initialize mesh
mesh = GtMesh();

%% Read vertices
nCoord = 3*nVert;
if binary
    coord = fread(fid, nCoord, ['3*' coordFMT '=>' coordFMT], 0 , 'ieee-be');
else
    [coord, sz] = fscanf(fid, '%f %f %f', nCoord);
    if sz ~= nCoord
        warning('gtVTKMeshReader:vertice_warning', ...
            'Problem in reading vertices.');
    end
end
%mesh.vertices = reshape(coord, 3, nVert)';
mesh.setVertices(reshape(coord, 3, nVert)');

%% Read polygons and lines
while ~feof(fid)
    str = fgets(fid);
    if feof(fid), break, end;
    while ~any(strcmp(str(1), {'P', 'L', 'C', 'F'}))
        str = fgets(fid);
        if feof(fid), break, end;
    end
    if strcmp(str(1:5), 'LINES')
        cellName = 'edges';
        nNodesByCell = 2;
        asciiFMT = '%*d %d %d\n';
    elseif strcmp(str(1:8), 'POLYGONS')
        cellName = 'faces';
        nNodesByCell = 3;
        asciiFMT = '%*d %d %d %d\n';
    elseif any([...
            strcmp(str(1:10), 'POINT_DATA') ...
            strcmp(str(1:9),  'CELL_DATA' ) ...
            strcmp(str(1:5),  'FIELD'     )  ])
        info.hasFields = true;
        warning('gtVTKMeshReader:field_warning', ...
            'Cell data/fields are not supported right now, so skipping...');
        break;
    else
        continue;
    end
    numbers = sscanf(str, '%*s %d %d', 2);
    nCells  = numbers(1);
    %nValues = numbers(2);
    nNodes = nNodesByCell * nCells;
    if binary
        nBytes = 4;
        fseek(fid, nBytes, 0);
        binFMT = sprintf('%d*int32=>uint32', nNodesByCell);
        connectivity = fread(fid, nNodes, binFMT, nBytes, 'ieee-be');
        cnt = length(connectivity);
        fseek(fid, -nNodesByCell, 0);
    else
        [connectivity, cnt] = fscanf(fid, asciiFMT, nNodes);
        connectivity = uint32(connectivity);
    end
    if cnt ~= nNodes
        if cnt > nNodes
            warning('gtVTKMeshReader:cell_warning', ...
                ['Too many values read for ' cellName ' connectivity (' ...
                    num2str(cnt) ')!']);
            connectivity = connectivity(1:nNodes);
        elseif cnt < nNodes
            gtError('gtVTKMeshReader:cell_warning', ...
                ['Not enough values read for ' cellName ' connectivity(' ...
                    num2str(cnt) ')!']);
        end
    end
    %mesh.(cellName) = reshape(connectivity+1, nNodesByCell, nCells)';
    mesh.setEntity(cellName, reshape(connectivity+1, nNodesByCell, nCells)');
    info.(['N' cellName]) = nCells;
end

%% Close the file end terminate
fclose(fid);

end % end of function