function [vert, face] = read_ply(fn)

% READ_PLY reads triangles, tetrahedrons or hexahedrons from a Stanford *.ply file
%
% Use as
%   [vert, face, prop, face_prop] = read_ply(filename)
%
% Documentation is provided on
%   http://paulbourke.net/dataformats/ply/
%   http://en.wikipedia.org/wiki/PLY_(file_format)
%
% See also WRITE_PLY, WRITE_VTK, READ_VTK

% Copyright (C) 2013-2018, Robert Oostenveld
% 
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(fn, 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the file starts with an ascii header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

line = readline(fid);
if ~strcmp(line, 'ply')
  fclose(fid);
  ft_error('unexpected header line');
end
line = readline(fid);

if ~strncmp(line, 'format', 6)
  fclose(fid);
  ft_error('unexpected header line');
else
  format = line(8:end);
end
line = readline(fid);

if ~strncmp(line, 'element vertex', 14)
  fclose(fid);
  ft_error('unexpected header line');
else
  nvert = str2double(line(16:end));
end
line = readline(fid);

prop = [];
while strncmp(line, 'property', 8)
  tok = tokenize(line);
  prop(end+1).format = tok{2};
  prop(end  ).name   = tok{3};
  line = readline(fid);
end

if ~strncmp(line, 'element face', 12)
  fclose(fid);
  ft_error('unexpected header line');
else
  nface = str2double(line(14:end));
end
line = readline(fid);

if ~strcmp(line, 'property list uchar int vertex_index')    && ...
    ~strcmp(line, 'property list uchar int vertex_indices') && ...
    ~strcmp(line, 'property list uchar uint vertex_index')  && ...
    ~strcmp(line, 'property list uchar uint vertex_indices')
  % the wikipedia documentation specifies vertex_index, but the OPTOCAT files have vertex_indices

  % it would not be very difficult to enhance the reader here with another
  % representation of the faces, i.e. something else than "uchar int"
  fclose(fid);
  ft_error('unexpected header line');
end
line = readline(fid);

while ~strcmp(line, 'end_header')
  line = readline(fid);
end

offset = ftell(fid);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the file continues with a section of data, which can be ascii or binary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch format

  case 'ascii 1.0'

    fid = fopen_or_error(fn, 'rt');
    fseek(fid, offset, 'cof');
    dat = fscanf(fid,'%f',[numel(prop), nvert])';
    for j=1:length(prop)
      vert.(prop(j).name) = dat(:,j);
    end
    face = zeros(nface,0);
    num = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fscanf(fid,'%f',1);
      face(i,1:num(i)) = fscanf(fid,'%f',num(i));
    end
    fclose(fid);

  case 'binary_little_endian 1.0'

    fid = fopen_or_error(fn, 'rb', 'l');
    fseek(fid, offset, 'cof');

    dat = zeros(nvert,length(prop));
    for i=1:nvert
      for j=1:length(prop)
        dat(i,j) = fread(fid, 1, prop(j).format);
      end % for each property
    end % for each vertex
    for j=1:length(prop)
      switch prop(j).format
        % the format can be one of the following: char uchar short ushort int uint float double int8 uint8 int16 uint16 int32 uint32 float32 float64
        case 'char'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'uchar'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'short'
          vert.(prop(j).name) = int16(dat(:,j));
        case 'ushort'
          vert.(prop(j).name) = uint16(dat(:,j));
        otherwise
          vert.(prop(j).name) = dat(:,j);
      end
    end
    clear dat;

    face = zeros(nface,0);
    num  = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fread(fid, 1, 'uint8');
      face(i,1:num(i)) = fread(fid, num(i), 'int32');
    end
    fclose(fid);

  case 'binary_big_endian 1.0'
    % this is exactly the same as the code above, except that the file is opened as big endian

    fid = fopen_or_error(fn, 'rb', 'b');
    fseek(fid, offset, 'cof');

    dat = zeros(nvert,length(prop));
    for i=1:nvert
      for j=1:length(prop)
        dat(i,j) = fread(fid, 1, prop(j).format);
      end % for each property
    end % for each vertex
    for j=1:length(prop)
      switch prop(j).format
        % the format can be one of the following: char uchar short ushort int uint float double int8 uint8 int16 uint16 int32 uint32 float32 float64
        case 'char'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'uchar'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'short'
          vert.(prop(j).name) = int16(dat(:,j));
        case 'ushort'
          vert.(prop(j).name) = uint16(dat(:,j));
        otherwise
          vert.(prop(j).name) = dat(:,j);
      end
    end
    clear dat;

    face = zeros(nface,0);
    num  = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fread(fid, 1, 'uint8');
      face(i,1:num(i)) = fread(fid, num(i), 'int32');
    end
    fclose(fid);

  otherwise
    ft_error('unsupported format');
end % switch
  

% each polygon can have a different number of elements
% mark all invalid entries with nan
if any(num<size(face,2))
  for i=1:nface
    face(i,(num(i)+1):end) = nan;
  end
end % if

if numel(face)>0
  % MATLAB indexes start at 1, inside the file they start at 0
  face = face+1;
end

end % function read_ply


function line = readline(fid)
% read the next line from the ascii header, skip all comment lines
line = fgetl(fid);
while strncmp(line, 'comment', 7)
  line = fgetl(fid);
end
end % function readline

