ply_read

*****************************************************************************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).

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    25 February 2007

  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.

0001 function [ Elements, varargout ] = ply_read ( Path, Str )
0002 
0003 %*****************************************************************************80
0004 %
0005 %% PLY_READ reads a PLY 3D data file.
0006 %
0007 %   [DATA,COMMENTS] = PLY_READ(FILENAME) reads a version 1.0 PLY file
0008 %   FILENAME and returns a structure DATA.  The fields in this structure
0009 %   are defined by the PLY header; each element type is a field and each
0010 %   element property is a subfield.  If the file contains any comments,
0011 %   they are returned in a cell string array COMMENTS.
0012 %
0013 %   [TRI,PTS] = PLY_READ(FILENAME,'tri') or
0014 %   [TRI,PTS,DATA,COMMENTS] = PLY_READ(FILENAME,'tri') converts vertex
0015 %   and face data into triangular connectivity and vertex arrays.  The
0016 %   mesh can then be displayed using the TRISURF command.
0017 %
0018 %   Note: This function is slow for large mesh files (+50K faces),
0019 %   especially when reading data with list type properties.
0020 %
0021 %   Example:
0022 %   [Tri,Pts] = PLY_READ('cow.ply','tri');
0023 %   trisurf(Tri,Pts(:,1),Pts(:,2),Pts(:,3));
0024 %   colormap(gray); axis equal;
0025 %
0026 %  Discussion:
0027 %
0028 %    The original version of this program had a mistake that meant it
0029 %    did not properly triangulate files whose faces were not already triangular.
0030 %    This has been corrected (JVB, 25 February 2007).
0031 %
0032 %  Licensing:
0033 %
0034 %    This code is distributed under the GNU LGPL license.
0035 %
0036 %  Modified:
0037 %
0038 %    25 February 2007
0039 %
0040 %  Author:
0041 %
0042 %    Pascal Getreuer 2004
0043 %
0044 %  Parameters:
0045 %
0046 %  Local Parameters:
0047 %
0048 %    COMMENTS, any comments from the file.
0049 %
0050 %    ELEMENTCOUNT, the number of each type of element in file.
0051 %
0052 %    ELEMENTS, the element data.
0053 %
0054 %    PROPERTYTYPES, the element property types.
0055 %
0056 %    SIZEOF, size in bytes of each type.
0057 %
0058 
0059 %
0060 %  Open the input file in "read text" mode.
0061 %
0062   [ fid, Msg ] = fopen ( Path, 'rt' );
0063 
0064   if ( fid == -1 )
0065     error ( Msg );
0066   end
0067 
0068   Buf = fscanf ( fid, '%s', 1 );
0069 
0070   if ( ~strcmp ( Buf, 'ply' ) )
0071     fclose ( fid );
0072     error('Not a PLY file.');
0073   end
0074 %
0075 %  Read the header.
0076 %
0077   Position = ftell(fid);
0078   Format = '';
0079   NumComments = 0;
0080   Comments = {};
0081   NumElements = 0;
0082   NumProperties = 0;
0083   Elements = [];
0084   ElementCount = [];
0085   PropertyTypes = [];
0086   ElementNames = {};  % list of element names in the order they are stored in the file
0087   PropertyNames = [];  % structure of lists of property names
0088 
0089   while ( 1 )
0090 %
0091 %  Read a line from the file.
0092 %
0093     Buf = fgetl ( fid );
0094     BufRem = Buf;
0095     Token = {};
0096     Count = 0;
0097 %
0098 %  Split the line into tokens.
0099 %
0100     while ( ~isempty(BufRem) )
0101 
0102       [ tmp, BufRem ] = strtok(BufRem);
0103 %
0104 %  Count the tokens.
0105 %
0106       if ( ~isempty ( tmp ) )
0107         Count = Count + 1;
0108         Token{Count} = tmp;
0109       end
0110 
0111     end
0112 %
0113 %  Parse the line.
0114 %
0115     if ( Count )
0116 
0117       switch lower ( Token{1} )
0118 %
0119 %  Read the data format.
0120 %
0121       case 'format'
0122 
0123         if ( 2 <= Count )
0124 
0125           Format = lower ( Token{2} );
0126 
0127           if ( Count == 3 & ~strcmp ( Token{3}, '1.0' ) )
0128             fclose ( fid );
0129             error('Only PLY format version 1.0 supported.');
0130           end
0131         end
0132 %
0133 %  Read a comment.
0134 %
0135       case 'comment'
0136 
0137         NumComments = NumComments + 1;
0138         Comments{NumComments} = '';
0139         for i = 2 : Count
0140           Comments{NumComments} = [Comments{NumComments},Token{i},' '];
0141         end
0142 %
0143 %  Read an element name.
0144 %
0145       case 'element'
0146 
0147         if ( 3 <= Count )
0148 
0149           if ( isfield(Elements,Token{2}) )
0150             fclose ( fid );
0151             error(['Duplicate element name, ''',Token{2},'''.']);
0152           end
0153 
0154           NumElements = NumElements + 1;
0155           NumProperties = 0;
0156           Elements = setfield(Elements,Token{2},[]);
0157           PropertyTypes = setfield(PropertyTypes,Token{2},[]);
0158           ElementNames{NumElements} = Token{2};
0159           PropertyNames = setfield(PropertyNames,Token{2},{});
0160           CurElement = Token{2};
0161           ElementCount(NumElements) = str2double(Token{3});
0162 
0163           if ( isnan(ElementCount(NumElements)) )
0164             fclose ( fid );
0165             error(['Bad element definition: ',Buf]);
0166           end
0167 
0168         else
0169 
0170           error(['Bad element definition: ',Buf]);
0171 
0172         end
0173 %
0174 %  Read an element property.
0175 %
0176       case 'property'
0177 
0178         if ( ~isempty(CurElement) & Count >= 3 )
0179 
0180           NumProperties = NumProperties + 1;
0181           eval(['tmp=isfield(Elements.',CurElement,',Token{Count});'],...
0182             'fclose(fid);error([''Error reading property: '',Buf])');
0183 
0184           if ( tmp )
0185             error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
0186           end
0187 %
0188 %  Add property subfield to Elements.
0189 %
0190           eval(['Elements.',CurElement,'.',Token{Count},'=[];'], ...
0191             'fclose(fid);error([''Error reading property: '',Buf])');
0192 %
0193 %  Add property subfield to PropertyTypes and save type.
0194 %
0195           eval(['PropertyTypes.',CurElement,'.',Token{Count},'={Token{2:Count-1}};'], ...
0196             'fclose(fid);error([''Error reading property: '',Buf])');
0197 %
0198 %  Record property name order.
0199 %
0200           eval(['PropertyNames.',CurElement,'{NumProperties}=Token{Count};'], ...
0201             'fclose(fid);error([''Error reading property: '',Buf])');
0202 
0203          else
0204 
0205            fclose ( fid );
0206 
0207            if ( isempty(CurElement) )
0208              error(['Property definition without element definition: ',Buf]);
0209            else
0210              error(['Bad property definition: ',Buf]);
0211            end
0212 
0213          end
0214 %
0215 %  End of header.
0216 %
0217         case 'end_header'
0218           break;
0219 
0220       end
0221     end
0222   end
0223 %
0224 %  Set reading for specified data format.
0225 %
0226   if ( isempty ( Format ) )
0227     warning('Data format unspecified, assuming ASCII.');
0228     Format = 'ascii';
0229   end
0230 
0231   switch Format
0232 
0233     case 'ascii'
0234       Format = 0;
0235     case 'binary_little_endian'
0236       Format = 1;
0237     case 'binary_big_endian'
0238       Format = 2;
0239     otherwise
0240       fclose ( fid );
0241       error(['Data format ''',Format,''' not supported.']);
0242 
0243   end
0244 %
0245 %  Read the rest of the file as ASCII data...
0246 %
0247   if ( ~Format )
0248     Buf = fscanf ( fid, '%f' );
0249     BufOff = 1;
0250   else
0251 %
0252 %  ...or, close the file, and reopen in "read binary" mode.
0253 %
0254     fclose ( fid );
0255 %
0256 %  Reopen the binary file as LITTLE_ENDIAN or BIG_ENDIAN.
0257 %
0258     if ( Format == 1 )
0259       fid = fopen ( Path, 'r', 'ieee-le.l64' );
0260     else
0261       fid = fopen ( Path, 'r', 'ieee-be.l64' );
0262     end
0263 %
0264 %  Find the end of the header again.
0265 %  Using ftell on the old handle doesn't give the correct position.
0266 %
0267     BufSize = 8192;
0268     Buf = [ blanks(10), char(fread(fid,BufSize,'uchar')') ];
0269     i = [];
0270     tmp = -11;
0271 
0272     while ( isempty(i) )
0273 
0274       i = findstr(Buf,['end_header',13,10]);   % look for end_header + CR/LF
0275       i = [i,findstr(Buf,['end_header',10])];  % look for end_header + LF
0276 
0277       if ( isempty(i) )
0278         tmp = tmp + BufSize;
0279         Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
0280       end
0281 
0282     end
0283 %
0284 %  seek to just after the line feed
0285 %
0286     fseek ( fid, i + tmp + 11 + (Buf(i + 10) == 13), -1 );
0287 
0288   end
0289 %
0290 %  Read element data.
0291 %
0292 %  PLY and MATLAB data types (for fread)
0293 %
0294   PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
0295     'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};
0296 
0297   MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};
0298 
0299   SizeOf = [1,1,2,2,4,4,4,8];
0300 
0301   for i = 1 : NumElements
0302 %
0303 %  get current element property information
0304 %
0305     eval(['CurPropertyNames=PropertyNames.',ElementNames{i},';']);
0306     eval(['CurPropertyTypes=PropertyTypes.',ElementNames{i},';']);
0307     NumProperties = size(CurPropertyNames,2);
0308 
0309 %   fprintf('Reading %s...\n',ElementNames{i});
0310 %
0311 %  Read ASCII data.
0312 %
0313     if ( ~Format )
0314 
0315       for j = 1 : NumProperties
0316 
0317         Token = getfield(CurPropertyTypes,CurPropertyNames{j});
0318 
0319         if ( strcmpi(Token{1},'list') )
0320           Type(j) = 1;
0321         else
0322           Type(j) = 0;
0323         end
0324 
0325       end
0326 %
0327 %  Parse the buffer.
0328 %
0329       if ( ~any(Type) )
0330 
0331          % no list types
0332 
0333         Data = reshape ( ...
0334           Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1), ...
0335           NumProperties, ElementCount(i) )';
0336 
0337         BufOff = BufOff + ElementCount(i) * NumProperties;
0338 
0339       else
0340 
0341         ListData = cell(NumProperties,1);
0342 
0343         for k = 1 : NumProperties
0344           ListData{k} = cell(ElementCount(i),1);
0345         end
0346 %
0347 % list type
0348 %
0349         for j = 1 : ElementCount(i)
0350           for k = 1 : NumProperties
0351 
0352             if ( ~Type(k) )
0353               Data(j,k) = Buf(BufOff);
0354               BufOff = BufOff + 1;
0355             else
0356               tmp = Buf(BufOff);
0357               ListData{k}{j} = Buf(BufOff+(1:tmp))';
0358               BufOff = BufOff + tmp + 1;
0359             end
0360 
0361           end
0362         end
0363 
0364       end
0365 %
0366 %  Read binary data.
0367 %
0368     else
0369 % translate PLY data type names to MATLAB data type names
0370       ListFlag = 0;  % = 1 if there is a list type
0371       SameFlag = 1;     % = 1 if all types are the same
0372 
0373       for j = 1 : NumProperties
0374 
0375         Token = getfield(CurPropertyTypes,CurPropertyNames{j});
0376 %
0377 %  Non-list type.
0378 %
0379         if ( ~strcmp(Token{1},'list' ) )
0380 
0381           tmp = rem(strmatch(Token{1},PlyTypeNames,'exact')-1,8)+1;
0382 
0383           if ( ~isempty(tmp) )
0384 
0385             TypeSize(j) = SizeOf(tmp);
0386             Type{j} = MatlabTypeNames{tmp};
0387             TypeSize2(j) = 0;
0388             Type2{j} = '';
0389 
0390             SameFlag = SameFlag & strcmp(Type{1},Type{j});
0391 
0392           else
0393 
0394             fclose(fid);
0395             error(['Unknown property data type, ''',Token{1},''', in ', ...
0396               ElementNames{i},'.',CurPropertyNames{j},'.']);
0397 
0398           end
0399 
0400         else           % list type
0401 
0402           if ( length(Token) == 3 )
0403 
0404             ListFlag = 1;
0405             SameFlag = 0;
0406             tmp = rem(strmatch(Token{2},PlyTypeNames,'exact')-1,8)+1;
0407             tmp2 = rem(strmatch(Token{3},PlyTypeNames,'exact')-1,8)+1;
0408 
0409             if ( ~isempty(tmp) & ~isempty(tmp2) )
0410               TypeSize(j) = SizeOf(tmp);
0411               Type{j} = MatlabTypeNames{tmp};
0412               TypeSize2(j) = SizeOf(tmp2);
0413               Type2{j} = MatlabTypeNames{tmp2};
0414             else
0415               fclose(fid);
0416               error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
0417                 ElementNames{i},'.',CurPropertyNames{j},'.']);
0418             end
0419 
0420           else
0421 
0422             fclose(fid);
0423             error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);
0424 
0425           end
0426 
0427         end
0428 
0429       end
0430 
0431 % read file
0432 
0433       if ( ~ListFlag )
0434 %
0435 %  No list types, all the same type (fast)
0436 %
0437         if ( SameFlag )
0438 
0439           Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
0440 %
0441 %  No list types, mixed type.
0442 %
0443         else
0444 
0445           Data = zeros(ElementCount(i),NumProperties);
0446 
0447           for j = 1 : ElementCount(i)
0448             for k = 1 : NumProperties
0449               Data(j,k) = fread(fid,1,Type{k});
0450             end
0451           end
0452 
0453         end
0454 
0455       else
0456 
0457         ListData = cell(NumProperties,1);
0458 
0459         for k = 1 : NumProperties
0460           ListData{k} = cell(ElementCount(i),1);
0461         end
0462 
0463         if ( NumProperties == 1 )
0464 
0465           BufSize = 512;
0466           SkipNum = 4;
0467           j = 0;
0468 %
0469 %  List type, one property (fast if lists are usually the same length)
0470 %
0471           while ( j < ElementCount(i) )
0472 
0473             BufSize = min(ElementCount(i)-j,BufSize);
0474             Position = ftell(fid);
0475 %
0476 %  Read in BufSize count values, assuming all counts = SkipNum
0477 %
0478             [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
0479             Miss = find(Buf ~= SkipNum);     % find first count that is not SkipNum
0480             fseek(fid,Position + TypeSize(1),-1);   % seek back to after first count
0481 
0482             if ( isempty(Miss) )
0483 % all counts are SkipNum
0484               Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
0485               fseek(fid,-TypeSize(1),0);     % undo last skip
0486 
0487               for k = 1:BufSize
0488                 ListData{1}{j+k} = Buf(k,:);
0489               end
0490 
0491               j = j + BufSize;
0492               BufSize = floor(1.5*BufSize);
0493 
0494             else
0495 %
0496 %  Some counts are SkipNum.
0497 %
0498               if ( 1 < Miss(1) )
0499 
0500                 Buf2 = fread(fid,[SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
0501 
0502                 for k = 1:Miss(1)-1
0503                   ListData{1}{j+k} = Buf2(k,:);
0504                 end
0505 
0506                 j = j + k;
0507 
0508               end
0509 %
0510 %  Read in the list with the missed count.
0511 %
0512               SkipNum = Buf(Miss(1));
0513               j = j + 1;
0514               ListData{1}{j} = fread(fid,[1,SkipNum],Type2{1});
0515               BufSize = ceil(0.6*BufSize);
0516 
0517             end
0518           end
0519 
0520         else
0521 %
0522 %  List type(s), multiple properties (slow)
0523 %
0524           Data = zeros(ElementCount(i),NumProperties);
0525 
0526           for j = 1:ElementCount(i)
0527             for k = 1:NumProperties
0528 
0529               if ( isempty(Type2{k}) )
0530                 Data(j,k) = fread(fid,1,Type{k});
0531               else
0532                 tmp = fread(fid,1,Type{k});
0533                 ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
0534               end
0535 
0536             end
0537           end
0538         end
0539       end
0540     end
0541 %
0542 %  Put data into Elements structure
0543 %
0544     for k = 1 : NumProperties
0545 
0546       if ( ( ~Format && ~Type(k) ) || (Format && isempty(Type2{k})) )
0547         eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=Data(:,k);']);
0548       else
0549         eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=ListData{k};']);
0550       end
0551 
0552     end
0553 
0554   end
0555 
0556   clear Data
0557   clear ListData;
0558 
0559   fclose ( fid );
0560 %
0561 %  Output the data as a triangular mesh pair.
0562 %
0563   if ( ( nargin > 1 && strcmpi(Str,'Tri') ) || nargout > 2 )
0564 %
0565 %  Find vertex element field
0566 %
0567     Name = {'vertex','Vertex','point','Point','pts','Pts'};
0568     Names = [];
0569 
0570     for i = 1 : length(Name)
0571 
0572       if ( any ( strcmp ( ElementNames, Name{i} ) ) )
0573         Names = getfield(PropertyNames,Name{i});
0574         Name = Name{i};
0575         break;
0576       end
0577 
0578     end
0579 
0580     if ( any(strcmp(Names,'x')) & any(strcmp(Names,'y')) & any(strcmp(Names,'z')) )
0581       eval(['varargout{1}=[Elements.',Name,'.x,Elements.',Name,'.y,Elements.',Name,'.z];']);
0582     else
0583       varargout{1} = zeros(1,3);
0584     end
0585     varargout{1} = varargout{1}';
0586 
0587     varargout{2} = Elements;
0588     varargout{3} = Comments;
0589     Elements = [];
0590 %
0591 % Find face element field
0592 %
0593     Name = {'face','Face','poly','Poly','tri','Tri'};
0594     Names = [];
0595 
0596     for i = 1 : length(Name)
0597       if ( any(strcmp(ElementNames,Name{i})) )
0598         Names = getfield(PropertyNames,Name{i});
0599         Name = Name{i};
0600         break;
0601       end
0602     end
0603 
0604     if ( ~isempty(Names) )
0605       % find vertex indices property subfield
0606       PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};
0607 
0608       for i = 1 : length(PropertyName)
0609         if ( any(strcmp(Names,PropertyName{i})) )
0610           PropertyName = PropertyName{i};
0611           break;
0612         end
0613       end
0614 %
0615 %  Convert face index list to triangular connectivity.
0616 %
0617       if ( ~iscell(PropertyName) )
0618 
0619         eval(['FaceIndices=varargout{2}.',Name,'.',PropertyName,';']);
0620         N = length(FaceIndices);
0621         Elements = zeros(3,N*2);
0622         Extra = 0;
0623 
0624         for k = 1 : N
0625 
0626           Elements(1:3,k) = FaceIndices{k}(1:3)';
0627 %
0628 %  The original code had an error in the following loop.
0629 %
0630           for j = 4 : length(FaceIndices{k})
0631             Extra = Extra + 1;
0632             Elements(1,N + Extra) = FaceIndices{k}(1);
0633             Elements(2,N + Extra) = FaceIndices{k}(j-1);
0634             Elements(3,N + Extra) = FaceIndices{k}(j);
0635           end
0636 
0637         end
0638 %
0639 %  Add 1 to each vertex value; PLY vertices are zero based.
0640 %
0641         Elements = Elements(:,1:N+Extra) + 1;
0642 
0643       end
0644     end
0645 
0646   else
0647 
0648     varargout{1} = Comments;
0649 
0650   end
0651 
0652   return
0653 end