function [ Elements, varargout ] = ply_read ( Path, Str )
%% PLY_READ reads a PLY 3D data file.
%
%   [DATA,COMMENTS] = PLY_READ(FILENAME) reads a version 1.0 PLY file
%   FILENAME and returns a structure DATA.  The fields in this structure
%   are defined by the PLY header; each element type is a field and each
%   element property is a subfield.  If the file contains any comments,
%   they are returned in a cell string array COMMENTS.
%
%   [TRI,PTS] = PLY_READ(FILENAME,'tri') or
%   [TRI,PTS,DATA,COMMENTS] = PLY_READ(FILENAME,'tri') converts vertex
%   and face data into triangular connectivity and vertex arrays.  The
%   mesh can then be displayed using the TRISURF command.
%
%   Note: This function is slow for large mesh files (+50K faces),
%   especially when reading data with list type properties.
%
%   Example:
%   [Tri,Pts] = PLY_READ('cow.ply','tri');
%   trisurf(Tri,Pts(:,1),Pts(:,2),Pts(:,3));
%   colormap(gray); axis equal;
%
%  Discussion:
%
%    The original version of this program had a mistake that meant it
%    did not properly triangulate files whose faces were not already triangular.
%    This has been corrected (JVB, 25 February 2007).
%
%    Glenn Ramsey pointed out and corrected a problem that occurred
%    with an uninitialized value of Type2, 27 August 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 August 2012
%
%  Author:
%
%    Pascal Getreuer 2004
%
%  Parameters:
%
%  Local Parameters:
%
%    COMMENTS, any comments from the file.
%
%    ELEMENTCOUNT, the number of each type of element in file.
%
%    ELEMENTS, the element data.
%
%    PROPERTYTYPES, the element property types.
%
%    SIZEOF, size in bytes of each type.
%

%
%  Open the input file in "read text" mode.
%
[ fid, Msg ] = fopen ( Path, 'rt' );

if ( fid == -1 )
  error ( Msg );
end

Buf = fscanf ( fid, '%s', 1 );

if ( ~strcmp ( Buf, 'ply' ) )
  fclose ( fid );
  error('Not a PLY file.');
end
%
%  Read the header.
%
Position = ftell(fid);
Format = '';
NumComments = 0;
Comments = {};
NumElements = 0;
NumProperties = 0;
Elements = [];
ElementCount = [];
PropertyTypes = [];
ElementNames = {};  % list of element names in the order they are stored in the file
PropertyNames = [];  % structure of lists of property names

while ( 1 )
  %
  %  Read a line from the file.
  %
  Buf = fgetl ( fid );
  BufRem = Buf;
  Token = {};
  Count = 0;
  %
  %  Split the line into tokens.
  %
  while ( ~isempty(BufRem) )
    
    [ tmp, BufRem ] = strtok(BufRem);
    %
    %  Count the tokens.
    %
    if ( ~isempty ( tmp ) )
      Count = Count + 1;
      Token{Count} = tmp;
    end
    
  end
  %
  %  Parse the line.
  %
  if ( Count )
    
    switch lower ( Token{1} )
      %
      %  Read the data format.
      %
      case 'format'
        
        if ( 2 <= Count )
          
          Format = lower ( Token{2} );
          
          if ( Count == 3 & ~strcmp ( Token{3}, '1.0' ) )
            fclose ( fid );
            error('Only PLY format version 1.0 supported.');
          end
        end
        %
        %  Read a comment.
        %
      case 'comment'
        
        NumComments = NumComments + 1;
        Comments{NumComments} = '';
        for i = 2 : Count
          Comments{NumComments} = [Comments{NumComments},Token{i},' '];
        end
        %
        %  Read an element name.
        %
      case 'element'
        
        if ( 3 <= Count )
          
          if ( isfield(Elements,Token{2}) )
            fclose ( fid );
            error(['Duplicate element name, ''',Token{2},'''.']);
          end
          
          NumElements = NumElements + 1;
          NumProperties = 0;
          Elements = setfield(Elements,Token{2},[]);
          PropertyTypes = setfield(PropertyTypes,Token{2},[]);
          ElementNames{NumElements} = Token{2};
          PropertyNames = setfield(PropertyNames,Token{2},{});
          CurElement = Token{2};
          ElementCount(NumElements) = str2double(Token{3});
          
          if ( isnan(ElementCount(NumElements)) )
            fclose ( fid );
            error(['Bad element definition: ',Buf]);
          end
          
        else
          
          error(['Bad element definition: ',Buf]);
          
        end
        %
        %  Read an element property.
        %
      case 'property'
        
        if ( ~isempty(CurElement) & Count >= 3 )
          
          NumProperties = NumProperties + 1;
          eval(['tmp=isfield(Elements.',CurElement,',Token{Count});'],...
            'fclose(fid);error([''Error reading property: '',Buf])');
          
          if ( tmp )
            error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
          end
          %
          %  Add property subfield to Elements.
          %
          eval(['Elements.',CurElement,'.',Token{Count},'=[];'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');
          %
          %  Add property subfield to PropertyTypes and save type.
          %
          eval(['PropertyTypes.',CurElement,'.',Token{Count},'={Token{2:Count-1}};'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');
          %
          %  Record property name order.
          %
          eval(['PropertyNames.',CurElement,'{NumProperties}=Token{Count};'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');
          
        else
          
          fclose ( fid );
          
          if ( isempty(CurElement) )
            error(['Property definition without element definition: ',Buf]);
          else
            error(['Bad property definition: ',Buf]);
          end
          
        end
        %
        %  End of header.
        %
      case 'end_header'
        break;
        
    end
  end
end
%
%  Set reading for specified data format.
%
if ( isempty ( Format ) )
  warning('Data format unspecified, assuming ASCII.');
  Format = 'ascii';
end

switch Format
  
  case 'ascii'
    Format = 0;
  case 'binary_little_endian'
    Format = 1;
  case 'binary_big_endian'
    Format = 2;
  otherwise
    fclose ( fid );
    error(['Data format ''',Format,''' not supported.']);
    
end
%
%  Read the rest of the file as ASCII data...
%
if ( ~Format )
  Buf = fscanf ( fid, '%f' );
  BufOff = 1;
else
  %
  %  ...or, close the file, and reopen in "read binary" mode.
  %
  fclose ( fid );
  %
  %  Reopen the binary file as LITTLE_ENDIAN or BIG_ENDIAN.
  %
  if ( Format == 1 )
    fid = fopen ( Path, 'r', 'ieee-le.l64' );
  else
    fid = fopen ( Path, 'r', 'ieee-be.l64' );
  end
  %
  %  Find the end of the header again.
  %  Using ftell on the old handle doesn't give the correct position.
  %
  BufSize = 8192;
  Buf = [ blanks(10), char(fread(fid,BufSize,'uchar')') ];
  i = [];
  tmp = -11;
  
  while ( isempty(i) )
    
    i = findstr(Buf,['end_header',13,10]);   % look for end_header + CR/LF
    i = [i,findstr(Buf,['end_header',10])];  % look for end_header + LF
    
    if ( isempty(i) )
      tmp = tmp + BufSize;
      Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
    end
    
  end
  %
  %  seek to just after the line feed
  %
  fseek ( fid, i + tmp + 11 + (Buf(i + 10) == 13), -1 );
  
end
%
%  Read element data.
%
%  PLY and MATLAB data types (for fread)
%
PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
  'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};

MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};

SizeOf = [1,1,2,2,4,4,4,8];

for i = 1 : NumElements
  %
  %  get current element property information
  %
  eval(['CurPropertyNames=PropertyNames.',ElementNames{i},';']);
  eval(['CurPropertyTypes=PropertyTypes.',ElementNames{i},';']);
  NumProperties = size(CurPropertyNames,2);
  
  %   fprintf('Reading %s...\n',ElementNames{i});
  %
  %  Read ASCII data.
  %
  if ( ~Format )
    
    for j = 1 : NumProperties
      
      Token = getfield(CurPropertyTypes,CurPropertyNames{j});
      
      if ( strcmpi(Token{1},'list') )
        Type(j) = 1;
      else
        Type(j) = 0;
      end
      %
      %  Glenn Ramsey 20120827
      %  Initialise Type2{} to prevent uninitialised value error.
      %
      Type2{j} = '';
      
    end
    %
    %  Parse the buffer.
    %
    if ( ~any(Type) )
      
      % no list types
      
      Data = reshape ( ...
        Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1), ...
        NumProperties, ElementCount(i) )';
      
      BufOff = BufOff + ElementCount(i) * NumProperties;
      
    else
      
      ListData = cell(NumProperties,1);
      
      for k = 1 : NumProperties
        ListData{k} = cell(ElementCount(i),1);
      end
      %
      % list type
      %
      for j = 1 : ElementCount(i)
        for k = 1 : NumProperties
          
          if ( ~Type(k) )
            Data(j,k) = Buf(BufOff);
            BufOff = BufOff + 1;
          else
            tmp = Buf(BufOff);
            ListData{k}{j} = Buf(BufOff+(1:tmp))';
            BufOff = BufOff + tmp + 1;
          end
          
        end
      end
      
    end
    %
    %  Read binary data.
    %
  else
    % translate PLY data type names to MATLAB data type names
    ListFlag = 0;  % = 1 if there is a list type
    SameFlag = 1;     % = 1 if all types are the same
    
    for j = 1 : NumProperties
      
      Token = getfield(CurPropertyTypes,CurPropertyNames{j});
      %
      %  Non-list type.
      %
      if ( ~strcmp(Token{1},'list' ) )
        
        tmp = rem(strmatch(Token{1},PlyTypeNames,'exact')-1,8)+1;
        
        if ( ~isempty(tmp) )
          
          TypeSize(j) = SizeOf(tmp);
          Type{j} = MatlabTypeNames{tmp};
          TypeSize2(j) = 0;
          Type2{j} = '';
          
          SameFlag = SameFlag & strcmp(Type{1},Type{j});
          
        else
          
          fclose(fid);
          error(['Unknown property data type, ''',Token{1},''', in ', ...
            ElementNames{i},'.',CurPropertyNames{j},'.']);
          
        end
        
      else           % list type
        
        if ( length(Token) == 3 )
          
          ListFlag = 1;
          SameFlag = 0;
          tmp = rem(strmatch(Token{2},PlyTypeNames,'exact')-1,8)+1;
          tmp2 = rem(strmatch(Token{3},PlyTypeNames,'exact')-1,8)+1;
          
          if ( ~isempty(tmp) & ~isempty(tmp2) )
            TypeSize(j) = SizeOf(tmp);
            Type{j} = MatlabTypeNames{tmp};
            TypeSize2(j) = SizeOf(tmp2);
            Type2{j} = MatlabTypeNames{tmp2};
          else
            fclose(fid);
            error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
              ElementNames{i},'.',CurPropertyNames{j},'.']);
          end
          
        else
          
          fclose(fid);
          error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);
          
        end
        
      end
      
    end
    
    % read file
    
    if ( ~ListFlag )
      %
      %  No list types, all the same type (fast)
      %
      if ( SameFlag )
        
        Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
        %
        %  No list types, mixed type.
        %
      else
        
        Data = zeros(ElementCount(i),NumProperties);
        
        for j = 1 : ElementCount(i)
          for k = 1 : NumProperties
            Data(j,k) = fread(fid,1,Type{k});
          end
        end
        
      end
      
    else
      
      ListData = cell(NumProperties,1);
      
      for k = 1 : NumProperties
        ListData{k} = cell(ElementCount(i),1);
      end
      
      if ( NumProperties == 1 )
        
        BufSize = 512;
        SkipNum = 4;
        j = 0;
        %
        %  List type, one property (fast if lists are usually the same length)
        %
        while ( j < ElementCount(i) )
          
          BufSize = min(ElementCount(i)-j,BufSize);
          Position = ftell(fid);
          %
          %  Read in BufSize count values, assuming all counts = SkipNum
          %
          [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
          Miss = find(Buf ~= SkipNum);     % find first count that is not SkipNum
          fseek(fid,Position + TypeSize(1),-1);   % seek back to after first count
          
          if ( isempty(Miss) )
            % all counts are SkipNum
            Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
            fseek(fid,-TypeSize(1),0);     % undo last skip
            
            for k = 1:BufSize
              ListData{1}{j+k} = Buf(k,:);
            end
            
            j = j + BufSize;
            BufSize = floor(1.5*BufSize);
            
          else
            %
            %  Some counts are SkipNum.
            %
            if ( 1 < Miss(1) )
              
              Buf2 = fread(fid,[SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
              
              for k = 1:Miss(1)-1
                ListData{1}{j+k} = Buf2(k,:);
              end
              
              j = j + k;
              
            end
            %
            %  Read in the list with the missed count.
            %
            SkipNum = Buf(Miss(1));
            j = j + 1;
            ListData{1}{j} = fread(fid,[1,SkipNum],Type2{1});
            BufSize = ceil(0.6*BufSize);
            
          end
        end
        
      else
        %
        %  List type(s), multiple properties (slow)
        %
        Data = zeros(ElementCount(i),NumProperties);
        
        for j = 1:ElementCount(i)
          for k = 1:NumProperties
            
            if ( isempty(Type2{k}) )
              Data(j,k) = fread(fid,1,Type{k});
            else
              tmp = fread(fid,1,Type{k});
              ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
            end
            
          end
        end
      end
    end
  end
  %
  %  Put data into Elements structure
  %
  if iscell(Type)
    TypeTest = false(size(Type));
  else
    TypeTest = Type;
  end
  for k = 1 : NumProperties
    
    if ( ~exist('ListData','var') || ( ~Format & ~TypeTest(k) ) || (Format & isempty(Type2{k})) )
      eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=Data(:,k);']);
    else
      eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=ListData{k};']);
    end
    
  end
  
end

clear Data
clear ListData;

fclose ( fid );
%
%  Output the data as a triangular mesh pair.
%
if ( ( nargin > 1 && strcmpi(Str,'Tri') ) || nargout > 2 )
  %
  %  Find vertex element field
  %
  Name = {'vertex','Vertex','point','Point','pts','Pts'};
  Names = [];
  
  for i = 1 : length(Name)
    
    if ( any ( strcmp ( ElementNames, Name{i} ) ) )
      Names = getfield(PropertyNames,Name{i});
      Name = Name{i};
      break;
    end
    
  end
  
  if ( any(strcmp(Names,'x')) & any(strcmp(Names,'y')) & any(strcmp(Names,'z')) )
    eval(['varargout{1}=[Elements.',Name,'.x,Elements.',Name,'.y,Elements.',Name,'.z];']);
  else
    varargout{1} = zeros(1,3);
  end
  varargout{1} = varargout{1}';
  
  varargout{2} = Elements;
  varargout{3} = Comments;
  Elements = [];
  %
  % Find face element field
  %
  Name = {'face','Face','poly','Poly','tri','Tri'};
  Names = [];
  
  for i = 1 : length(Name)
    if ( any(strcmp(ElementNames,Name{i})) )
      Names = getfield(PropertyNames,Name{i});
      Name = Name{i};
      break;
    end
  end
  
  if ( ~isempty(Names) )
    % find vertex indices property subfield
    PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};
    
    for i = 1 : length(PropertyName)
      if ( any(strcmp(Names,PropertyName{i})) )
        PropertyName = PropertyName{i};
        break;
      end
    end
    %
    %  Convert face index list to triangular connectivity.
    %
    if ( ~iscell(PropertyName) )
      
      eval(['FaceIndices=varargout{2}.',Name,'.',PropertyName,';']);
      N = length(FaceIndices);
      Elements = zeros(3,N*2);
      Extra = 0;
      
      for k = 1 : N
        
        Elements(1:3,k) = FaceIndices{k}(1:3)';
        %
        %  The original code had an error in the following loop.
        %
        for j = 4 : length(FaceIndices{k})
          Extra = Extra + 1;
          Elements(1,N + Extra) = FaceIndices{k}(1);
          Elements(2,N + Extra) = FaceIndices{k}(j-1);
          Elements(3,N + Extra) = FaceIndices{k}(j);
        end
        
      end
      %
      %  Add 1 to each vertex value; PLY vertices are zero based.
      %
      Elements = Elements(:,1:N+Extra) + 1;
      
    end
  end
  
else
  
  varargout{1} = Comments;
  
end

return
end