readMesh_ply

 Read mesh data stored in PLY (Stanford triangle) format.

   [V, F] = readMesh_ply(FNAME)

   Example
   readMesh_ply

   References
   Wrapper function for Gabriel Peyré's read_ply that is a wrapper
   function of Pascal Getreuer's plyread.

   See also
   meshes3d, readMesh, readMesh_off, readMesh_stl

0001 function varargout = readMesh_ply(fileName)
0002 % Read mesh data stored in PLY (Stanford triangle) format.
0003 %
0004 %   [V, F] = readMesh_ply(FNAME)
0005 %
0006 %   Example
0007 %   readMesh_ply
0008 %
0009 %   References
0010 %   Wrapper function for Gabriel Peyré's read_ply that is a wrapper
0011 %   function of Pascal Getreuer's plyread.
0012 %
0013 %   See also
0014 %   meshes3d, readMesh, readMesh_off, readMesh_stl
0015 
0016 % ------
0017 % Author: David Legland
0018 % e-mail: david.legland@inrae.fr
0019 % Created: 2018-04-26,    using Matlab 9.4.0.813654 (R2018a)
0020 % Copyright 2018 INRA - Cepia Software Platform.
0021 
0022 
0023 %% open file
0024 f = fopen(fileName, 'r');
0025 if f == -1
0026     error('matGeom:readMesh_ply:FileNotFound', ...
0027         ['Could not find file: ' fileName]);
0028 end
0029 
0030 [vertices, faces] = read_ply(fileName);
0031 varargout = formatMeshOutput(nargout, vertices, faces);
0032 
0033 end
0034 
0035 function [vertex,face,d,c] = read_ply(filename)
0036 % read_ply - read data from PLY file.
0037 %
0038 %   [vertex,face] = read_ply(filename);
0039 %
0040 %   'vertex' is a 'nb.vert x 3' array specifying the position of the vertices.
0041 %   'face' is a 'nb.face x 3' array specifying the connectivity of the mesh.
0042 %
0043 %   Copyright (c) 2003 Gabriel Peyré
0044 
0045 possibleFacePropertyNames = ...
0046     {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};
0047 
0048 [d,c] = plyread(filename);
0049 facePropertyNameIdx = find(matches(possibleFacePropertyNames, fieldnames(d.face)));
0050 assert(length(facePropertyNameIdx) == 1)
0051 vi = d.face.(possibleFacePropertyNames{facePropertyNameIdx});
0052 lengths = cellfun('length',vi);
0053 maxLength = max(lengths);
0054 if all(maxLength == lengths)
0055     face = cell2mat(vi)+1;
0056 else
0057     face = cellfun(@(x) x+1, vi, 'uni',0);
0058 end
0059 vertex = [d.vertex.x, d.vertex.y, d.vertex.z];
0060 
0061 end
0062 
0063 
0064 function [Elements,varargout] = plyread(Path,Str)
0065 %PLYREAD   Read a PLY 3D data file.
0066 %   [DATA,COMMENTS] = PLYREAD(FILENAME) reads a version 1.0 PLY file
0067 %   FILENAME and returns a structure DATA.  The fields in this structure
0068 %   are defined by the PLY header; each element type is a field and each
0069 %   element property is a subfield.  If the file contains any comments,
0070 %   they are returned in a cell string array COMMENTS.
0071 %
0072 %   [TRI,PTS] = PLYREAD(FILENAME,'tri') or
0073 %   [TRI,PTS,DATA,COMMENTS] = PLYREAD(FILENAME,'tri') converts vertex
0074 %   and face data into triangular connectivity and vertex arrays.  The
0075 %   mesh can then be displayed using the TRISURF command.
0076 %
0077 %   Note: This function is slow for large mesh files (+50K faces),
0078 %   especially when reading data with list type properties.
0079 %
0080 %   Example:
0081 %   [Tri,Pts] = PLYREAD('cow.ply','tri');
0082 %   trisurf(Tri,Pts(:,1),Pts(:,2),Pts(:,3));
0083 %   colormap(gray); axis equal;
0084 %
0085 %   See also: PLYWRITE
0086 
0087 % Pascal Getreuer 2004
0088 
0089 [fid,Msg] = fopen(Path,'rt'); % open file in read text mode
0090 
0091 if fid == -1, error(Msg); end
0092 
0093 Buf = fscanf(fid,'%s',1);
0094 if ~strcmp(Buf,'ply')
0095     fclose(fid);
0096     error('Not a PLY file.');
0097 end
0098 
0099 %% read header
0100 ftell(fid);
0101 Format = '';
0102 NumComments = 0;
0103 Comments = {};      % for storing any file comments
0104 NumElements = 0;
0105 NumProperties = 0;
0106 Elements = [];      % structure for holding the element data
0107 ElementCount = [];  % number of each type of element in file
0108 PropertyTypes = []; % corresponding structure recording property types
0109 ElementNames = {};  % list of element names in the order they are stored in the file
0110 PropertyNames = [];    % structure of lists of property names
0111 
0112 while 1
0113     Buf = fgetl(fid); % read one line from file
0114     Token = split(Buf); % split line into tokens
0115     Count = length(Token); % count tokens
0116     
0117     if Count % parse line
0118         switch lower(Token{1})
0119             case 'format' % read data format
0120                 if Count >= 2
0121                     Format = lower(Token{2});
0122                     if Count == 3 && ~strcmp(Token{3},'1.0')
0123                         fclose(fid);
0124                         error('Only PLY format version 1.0 supported.');
0125                     end
0126                 end
0127             case 'comment' % read file comment
0128                 NumComments = NumComments + 1;
0129                 Comments{NumComments} = ''; %#ok<AGROW>
0130                 for i = 2:Count
0131                     Comments{NumComments} = [Comments{NumComments},Token{i},' '];
0132                 end
0133             case 'element' % element name
0134                 if Count >= 3
0135                     if isfield(Elements,Token{2})
0136                         fclose(fid);
0137                         error(['Duplicate element name, ''',Token{2},'''.']);
0138                     end
0139                     
0140                     NumElements = NumElements + 1;
0141                     NumProperties = 0;
0142                     Elements.(Token{2}) = [];
0143                     PropertyTypes.(Token{2}) =[];
0144                     ElementNames{NumElements} = Token{2}; %#ok<AGROW>
0145                     PropertyNames.(Token{2}) = {};
0146                     CurElement = Token{2};
0147                     ElementCount(NumElements) = str2double(Token{3}); %#ok<AGROW>
0148                     
0149                     if isnan(ElementCount(NumElements))
0150                         fclose(fid);
0151                         error(['Bad element definition: ',Buf]);
0152                     end
0153                 else
0154                     error(['Bad element definition: ',Buf]);
0155                 end
0156             case 'property'    % element property
0157                 if ~isempty(CurElement) && Count >= 3
0158                     NumProperties = NumProperties + 1;
0159                     
0160                     if isfield(Elements.(CurElement),Token{Count})
0161                         fclose(fid);
0162                         error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
0163                     end
0164                     
0165                     % add property subfield to Elements
0166                     Elements.(CurElement).(Token{Count}) = [];
0167                     % add property subfield to PropertyTypes and save type
0168                     PropertyTypes.(CurElement).(Token{Count}) = Token(2:Count-1);
0169                     % record property name order
0170                     PropertyNames.(CurElement){NumProperties} = Token{Count};
0171                 else
0172                     fclose(fid);
0173                     if isempty(CurElement)
0174                         error(['Property definition without element definition: ',Buf]);
0175                     else
0176                         error(['Bad property definition: ',Buf]);
0177                     end
0178                 end
0179             case 'end_header' % end of header, break from while loop
0180                 break;
0181         end
0182     end
0183 end
0184 
0185 %% set reading for specified data format
0186 if isempty(Format)
0187     warning('Data format unspecified, assuming ASCII.');
0188     Format = 'ascii';
0189 end
0190 
0191 switch Format
0192     case 'ascii'
0193         Format = 0;
0194     case 'binary_little_endian'
0195         Format = 1;
0196     case 'binary_big_endian'
0197         Format = 2;
0198     otherwise
0199         fclose(fid);
0200         error(['Data format ''',Format,''' not supported.']);
0201 end
0202 
0203 if ~Format
0204     Buf = fscanf(fid,'%f'); % read the rest of the file as ASCII data
0205     BufOff = 1;
0206 else
0207     % reopen the file in read binary mode
0208     fclose(fid);
0209     
0210     if Format == 1
0211         fid = fopen(Path,'r','ieee-le.l64'); % little endian
0212     else
0213         fid = fopen(Path,'r','ieee-be.l64'); % big endian
0214     end
0215     
0216     % find the end of the header again (using ftell on the old handle doesn't give the correct position)
0217     BufSize = 8192;
0218     Buf = [blanks(10),char(fread(fid,BufSize,'uchar')')];
0219     i = [];
0220     tmp = -11;
0221     
0222     while isempty(i)
0223         i = strfind(Buf,['end_header',13,10]);              % look for end_header + CR/LF
0224         i = [i,strfind(Buf,['end_header',10])]; %#ok<AGROW> % look for end_header + LF
0225         
0226         if isempty(i)
0227             tmp = tmp + BufSize;
0228             Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
0229         end
0230     end
0231     
0232     % seek to just after the line feed
0233     fseek(fid,i + tmp + 11 + (Buf(i + 10) == 13),-1);
0234 end
0235 
0236 
0237 %% read element data
0238 % PLY and MATLAB data types (for fread)
0239 PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
0240     'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};
0241 MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};
0242 SizeOf = [1,1,2,2,4,4,4,8];    % size in bytes of each type
0243 
0244 for i = 1:NumElements
0245     % get current element property information
0246     CurPropertyNames = PropertyNames.(ElementNames{i});
0247     CurPropertyTypes = PropertyTypes.(ElementNames{i});
0248     NumProperties = size(CurPropertyNames,2);
0249     
0250     % fprintf('Reading %s...\n',ElementNames{i});
0251     
0252     %% read ASCII data
0253     if ~Format
0254         Type = zeros(1,NumProperties);
0255         for j = 1:NumProperties
0256             Token = CurPropertyTypes.(CurPropertyNames{j});
0257             
0258             if strcmpi(Token{1},'list')
0259                 Type(j) = 1;
0260             end
0261         end
0262         
0263         % parse buffer
0264         if ~any(Type)
0265             % no list types
0266             Data = reshape(Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1),NumProperties,ElementCount(i))';
0267             BufOff = BufOff + ElementCount(i)*NumProperties;
0268         else
0269             ListData = cell(NumProperties,1);
0270             
0271             for k = 1:NumProperties
0272                 ListData{k} = cell(ElementCount(i),1);
0273             end
0274             
0275             % list type
0276             for j = 1:ElementCount(i)
0277                 for k = 1:NumProperties
0278                     if ~Type(k)
0279                         Data(j,k) = Buf(BufOff);
0280                         BufOff = BufOff + 1;
0281                     else
0282                         tmp = Buf(BufOff);
0283                         ListData{k}{j} = Buf(BufOff+(1:tmp))';
0284                         BufOff = BufOff + tmp + 1;
0285                     end
0286                 end
0287             end
0288         end
0289     else
0290         %% read binary data
0291         % translate PLY data type names to MATLAB data type names
0292         ListFlag = 0; % = 1 if there is a list type
0293         SameFlag = 1; % = 1 if all types are the same
0294         
0295         Type = cell(1,NumProperties);
0296         Type2 = Type;
0297         TypeSize = zeros(1,NumProperties);
0298         TypeSize2 = TypeSize;
0299         for j = 1:NumProperties
0300             Token = CurPropertyTypes.(CurPropertyNames{j});
0301             
0302             if ~strcmp(Token{1},'list')    % non-list type
0303                 tmp = rem(find(matches(PlyTypeNames,Token{1}))-1,8)+1;
0304                 
0305                 if ~isempty(tmp)
0306                     TypeSize(j) = SizeOf(tmp);
0307                     Type{j} = MatlabTypeNames{tmp};
0308                     TypeSize2(j) = 0;
0309                     Type2{j} = '';
0310                     
0311                     SameFlag = SameFlag & strcmp(Type{1},Type{j});
0312                 else
0313                     fclose(fid);
0314                     error(['Unknown property data type, ''',Token{1},''', in ', ...
0315                         ElementNames{i},'.',CurPropertyNames{j},'.']);
0316                 end
0317             else % list type
0318                 if length(Token) == 3
0319                     ListFlag = 1;
0320                     SameFlag = 0;
0321                     tmp = rem(find(matches(PlyTypeNames,Token{2}))-1,8)+1;
0322                     tmp2 = rem(find(matches(PlyTypeNames,Token{3}))-1,8)+1;
0323                     
0324                     if ~isempty(tmp) && ~isempty(tmp2)
0325                         TypeSize(j) = SizeOf(tmp);
0326                         Type{j} = MatlabTypeNames{tmp};
0327                         TypeSize2(j) = SizeOf(tmp2);
0328                         Type2{j} = MatlabTypeNames{tmp2};
0329                     else
0330                         fclose(fid);
0331                         error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
0332                             ElementNames{i},'.',CurPropertyNames{j},'.']);
0333                     end
0334                 else
0335                     fclose(fid);
0336                     error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);
0337                 end
0338             end
0339         end
0340         
0341         % read file
0342         if ~ListFlag
0343             if SameFlag
0344                 % no list types, all the same type (fast)
0345                 Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
0346             else
0347                 % no list types, mixed type
0348                 Data = zeros(ElementCount(i),NumProperties);
0349                 
0350                 for j = 1:ElementCount(i)
0351                     for k = 1:NumProperties
0352                         Data(j,k) = fread(fid,1,Type{k});
0353                     end
0354                 end
0355             end
0356         else
0357             ListData = cell(NumProperties,1);
0358             
0359             for k = 1:NumProperties
0360                 ListData{k} = cell(ElementCount(i),1);
0361             end
0362             
0363             if NumProperties == 1
0364                 BufSize = 512;
0365                 SkipNum = 4;
0366                 j = 0;
0367                 
0368                 % list type, one property (fast if lists are usually the same length)
0369                 while j < ElementCount(i)
0370                     Position = ftell(fid);
0371                     % read in BufSize count values, assuming all counts = SkipNum
0372                     [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
0373                     Miss = find(Buf ~= SkipNum); % find first count that is not SkipNum
0374                     fseek(fid,Position + TypeSize(1),-1); % seek back to after first count
0375                     
0376                     if isempty(Miss) % all counts are SkipNum
0377                         Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
0378                         fseek(fid,-TypeSize(1),0); % undo last skip
0379                         
0380                         for k = 1:BufSize
0381                             ListData{1}{j+k} = Buf(k,:);
0382                         end
0383                         
0384                         j = j + BufSize;
0385                         BufSize = floor(1.5*BufSize);
0386                     else
0387                         if Miss(1) > 1 % some counts are SkipNum
0388                             Buf2 = fread(fid,[SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1));
0389                             Buf2 = Buf2';
0390                             
0391                             for k = 1:Miss(1)-1
0392                                 ListData{1}{j+k} = Buf2(k,:);
0393                             end
0394                             
0395                             j = j + k;
0396                         end
0397                         % Alec: check if done and rewind one step
0398                         if j >= ElementCount(i)
0399                             fseek(fid,-1,0);
0400                             break;
0401                         end
0402                         
0403                         % read in the list with the missed count
0404                         SkipNum = Buf(Miss(1));
0405                         j = j + 1;
0406                         ListData{1}{j} = fread(fid,[1,SkipNum],Type2{1});
0407                         BufSize = ceil(0.6*BufSize);
0408                     end
0409                 end
0410             else
0411                 % list type(s), multiple properties (slow)
0412                 Data = zeros(ElementCount(i),NumProperties);
0413                 
0414                 for j = 1:ElementCount(i)
0415                     for k = 1:NumProperties
0416                         if isempty(Type2{k})
0417                             Data(j,k) = fread(fid,1,Type{k});
0418                         else
0419                             tmp = fread(fid,1,Type{k});
0420                             ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
0421                         end
0422                     end
0423                 end
0424             end
0425         end
0426     end
0427     
0428     % put data into Elements structure
0429     for k = 1:NumProperties
0430         if (~Format && ~Type(k)) || (Format && isempty(Type2{k}))
0431             Elements.(ElementNames{i}).(CurPropertyNames{k}) = Data(:,k);
0432         else
0433             Elements.(ElementNames{i}).(CurPropertyNames{k}) = ListData{k};
0434         end
0435     end
0436 end
0437 
0438 clear Data ListData;
0439 fclose(fid);
0440 
0441 if (nargin > 1 && strcmpi(Str,'Tri')) || nargout > 2
0442     % find vertex element field
0443     Name = {'vertex','Vertex','point','Point','pts','Pts'};
0444     Names = [];
0445     
0446     for i = 1:length(Name)
0447         if any(strcmp(ElementNames,Name{i}))
0448             Names = PropertyNames.(Name{i});
0449             Name = Name{i};
0450             break;
0451         end
0452     end
0453     
0454     if any(strcmp(Names,'x')) && any(strcmp(Names,'y')) && any(strcmp(Names,'z'))
0455         varargout{1} = [Elements.(Name).x, Elements.(Name).y, Elements.(Name).z];
0456     else
0457         varargout{1} = zeros(1,3);
0458     end
0459     
0460     varargout{2} = Elements;
0461     varargout{3} = Comments;
0462     Elements = [];
0463     
0464     % find face element field
0465     Name = {'face','Face','poly','Poly','tri','Tri'};
0466     Names = [];
0467     
0468     for i = 1:length(Name)
0469         if any(strcmp(ElementNames,Name{i}))
0470             Names = PropertyNames.(Name{i});
0471             Name = Name{i};
0472             break;
0473         end
0474     end
0475     
0476     if ~isempty(Names)
0477         % find vertex indices property subfield
0478         PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};
0479         
0480         for i = 1:length(PropertyName)
0481             if any(strcmp(Names,PropertyName{i}))
0482                 PropertyName = PropertyName{i};
0483                 break;
0484             end
0485         end
0486         
0487         if ~iscell(PropertyName)
0488             % convert face index lists to triangular connectivity
0489             FaceIndices = varargout{2}.(Name).PropertyName;
0490             N = length(FaceIndices);
0491             Elements = zeros(N*2,3);
0492             Extra = 0;
0493             
0494             for k = 1:N
0495                 Elements(k,:) = FaceIndices{k}(1:3);
0496                 
0497                 for j = 4:length(FaceIndices{k})
0498                     Extra = Extra + 1;
0499                     Elements(N + Extra,:) = [Elements(k,[1,j-1]),FaceIndices{k}(j)];
0500                 end
0501             end
0502             Elements = Elements(1:N+Extra,:) + 1;
0503         end
0504     end
0505 else
0506     varargout{1} = Comments;
0507 end
0508 
0509 end