ply_read.m

function [ Elements, varargout ] = ply_read ( Path, Str )

%*****************************************************************************80
%
%% 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
%
    for k = 1 : NumProperties

      if ( ( ~Format & ~Type(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