readMesh_ply

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

   [V, F] = readMesh_ply(FNAME)
   V is a NV-by-3 numeric array containing vertex coordinates,
   F is a NF-by-3 or NF-by-4 array containg vertex indices of each face.

   MESH = readMesh_ply(FNAME)
   Returns mesh vertex and face information into a single structure with
   fields 'vertices' and 'faces'.

   Example
     [v, f] = readMesh_ply('bunny_F1k.ply');
     trisurf(f, v(:,1),v(:,2),v(:,3));
     colormap(gray); axis equal;

   References
   Adapted from Gabriel Peyré's "read_ply" function, that is was a wrapper
   for the "plyread" function written by Pascal Getreuer.

   See also
       meshes3d, readMesh, readMesh_off, readMesh_stl

0001 function varargout = readMesh_ply(fileName)
0002 %READMESH_PLY Read mesh data stored in PLY (Stanford triangle) format.
0003 %
0004 %   [V, F] = readMesh_ply(FNAME)
0005 %   V is a NV-by-3 numeric array containing vertex coordinates,
0006 %   F is a NF-by-3 or NF-by-4 array containg vertex indices of each face.
0007 %
0008 %   MESH = readMesh_ply(FNAME)
0009 %   Returns mesh vertex and face information into a single structure with
0010 %   fields 'vertices' and 'faces'.
0011 %
0012 %   Example
0013 %     [v, f] = readMesh_ply('bunny_F1k.ply');
0014 %     trisurf(f, v(:,1),v(:,2),v(:,3));
0015 %     colormap(gray); axis equal;
0016 %
0017 %   References
0018 %   Adapted from Gabriel Peyré's "read_ply" function, that is was a wrapper
0019 %   for the "plyread" function written by Pascal Getreuer.
0020 %
0021 %   See also
0022 %       meshes3d, readMesh, readMesh_off, readMesh_stl
0023 
0024 % ------
0025 % Author: David Legland
0026 % E-mail: david.legland@inrae.fr
0027 % Created: 2018-04-26, using Matlab 9.4.0.813654 (R2018a)
0028 % Copyright 2018-2024 INRA - Cepia Software Platform
0029 
0030 %% Open file
0031 
0032 % check existence of file
0033 if ~exist(fileName, 'file')
0034     error('matGeom:readMesh_ply:FileNotFound', ...
0035         ['Could not find file: ' fileName]);
0036 end
0037 
0038 % open file in read text mode
0039 [fid, msg] = fopen(fileName, 'rt'); 
0040 
0041 % check file was properly open
0042 if fid == -1
0043     error(msg); 
0044 end
0045 
0046 % check file format marker
0047 buf = fscanf(fid, '%s', 1);
0048 if ~strcmp(buf, 'ply')
0049     fclose(fid);
0050     error('Not a PLY file.');
0051 end
0052 
0053 
0054 %% Read header
0055 
0056 % initialize empty fields
0057 ftell(fid);
0058 dataFormat = '';
0059 nComments = 0;
0060 comments = {};      % for storing any file comments
0061 nElements = 0;
0062 elements = [];      % structure for holding the element data
0063 elementCounts = []; % number of each type of element in file
0064 elementNames = {};  % list of element names in the order they are stored in the file
0065 nProperties = 0;
0066 propertyNames = [];    % structure of lists of property names
0067 propertyTypes = []; % corresponding structure recording property types
0068 
0069 % iterate over the lines of the file until we find the "end_header" line
0070 while true
0071     buf = fgetl(fid); % read one line from file
0072     tokens = split(buf); % split line into tokens
0073     nToks = length(tokens); % count tokens
0074     
0075     if nToks == 0
0076         continue;
0077     end
0078 
0079     % At the beginning of an element, the first token indicates the nature
0080     % of the element
0081     switch lower(tokens{1})
0082         case 'format'
0083             % read data format (ascii or binary)
0084             if nToks >= 2
0085                 dataFormat = lower(tokens{2});
0086                 if nToks == 3 && ~strcmp(tokens{3},'1.0')
0087                     fclose(fid);
0088                     error('Only PLY format version 1.0 supported.');
0089                 end
0090             end
0091 
0092         case 'comment'
0093             % read file comments
0094             nComments = nComments + 1;
0095             comments{nComments} = ''; %#ok<AGROW>
0096             for i = 2:nToks
0097                 comments{nComments} = [comments{nComments}, tokens{i}, ' '];
0098             end
0099 
0100         case 'element'
0101             % element name
0102             if nToks >= 3
0103                 % check element was not already initialized
0104                 if isfield(elements,tokens{2})
0105                     fclose(fid);
0106                     error(['Duplicate element name, ''',tokens{2},'''.']);
0107                 end
0108 
0109                 nElements = nElements + 1;
0110                 nProperties = 0;
0111                 elements.(tokens{2}) = [];
0112                 propertyTypes.(tokens{2}) =[];
0113                 elementNames{nElements} = tokens{2}; %#ok<AGROW>
0114                 propertyNames.(tokens{2}) = {};
0115                 element = tokens{2};
0116                 elementCounts(nElements) = str2double(tokens{3}); %#ok<AGROW>
0117 
0118                 if isnan(elementCounts(nElements))
0119                     fclose(fid);
0120                     error(['Bad element definition: ', buf]);
0121                 end
0122             else
0123                 error(['Bad element definition: ', buf]);
0124             end
0125 
0126         case 'property'
0127             % element property
0128             if ~isempty(element) && nToks >= 3
0129                 nProperties = nProperties + 1;
0130 
0131                 if isfield(elements.(element), tokens{nToks})
0132                     fclose(fid);
0133                     error(['Duplicate property name, ''',element,'.',tokens{2},'''.']);
0134                 end
0135 
0136                 % add property subfield to elements
0137                 elements.(element).(tokens{nToks}) = [];
0138                 % add property subfield to PropertyTypes and save type
0139                 propertyTypes.(element).(tokens{nToks}) = tokens(2:nToks-1);
0140                 % record property name order
0141                 propertyNames.(element){nProperties} = tokens{nToks};
0142             else
0143                 fclose(fid);
0144                 if isempty(element)
0145                     error(['Property definition without element definition: ', buf]);
0146                 else
0147                     error(['Bad property definition: ', buf]);
0148                 end
0149             end
0150 
0151         case 'end_header' 
0152             % end of header, break infinite loop
0153             break;
0154     end
0155 end
0156 
0157 
0158 %% Set reading for specified data format
0159 if isempty(dataFormat)
0160     warning('Data format unspecified, assuming ASCII.');
0161     dataFormat = 'ascii';
0162 end
0163 
0164 switch dataFormat
0165     case 'ascii'
0166         dataFormat = 0;
0167     case 'binary_little_endian'
0168         dataFormat = 1;
0169     case 'binary_big_endian'
0170         dataFormat = 2;
0171     otherwise
0172         fclose(fid);
0173         error('Data format ''%s'' not supported.', dataFormat);
0174 end
0175 
0176 if dataFormat == 0
0177     % read the rest of the file as ASCII data
0178     buf = fscanf(fid,'%f');
0179     offset = 1;
0180 else
0181     % reopen the file in read binary mode
0182     fclose(fid);
0183     
0184     if dataFormat == 1
0185         fid = fopen(fileName, 'r', 'ieee-le.l64'); % little endian
0186     else
0187         fid = fopen(fileName, 'r', 'ieee-be.l64'); % big endian
0188     end
0189     
0190     % find the end of the header again (using ftell on the old handle doesn't give the correct position)
0191     bufferSize = 8192;
0192     buf = [blanks(10),char(fread(fid,bufferSize,'uchar')')];
0193     i = [];
0194     tmp = -11;
0195     
0196     while isempty(i)
0197         % look for end_header + CR/LF
0198         i = strfind(buf, ['end_header', 13, 10]);
0199         % look for end_header + LF
0200         i = [i, strfind(buf, ['end_header', 10])]; %#ok<AGROW>
0201         
0202         if isempty(i)
0203             tmp = tmp + bufferSize;
0204             buf = [buf(bufferSize+1:bufferSize+10), char(fread(fid,bufferSize,'uchar')')];
0205         end
0206     end
0207     
0208     % seek to just after the line feed
0209     fseek(fid, i + tmp + 11 + (buf(i + 10) == 13), -1);
0210 end
0211 
0212 
0213 %% Read element data
0214 
0215 % PLY and MATLAB data types (for fread)
0216 plyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
0217     'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};
0218 matlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};
0219 dataTypeSizes = [1,1,2,2,4,4,4,8];    % size in bytes of each type
0220 
0221 % iterate over element types
0222 for i = 1:nElements
0223     % get current element property information
0224     currPropNames = propertyNames.(elementNames{i});
0225     currPropTypes = propertyTypes.(elementNames{i});
0226     nProperties = size(currPropNames,2);
0227     
0228     % fprintf('Reading %s...\n',ElementNames{i});
0229     
0230     if dataFormat == 0
0231         % read ASCII data
0232         type = zeros(1, nProperties);
0233         for j = 1:nProperties
0234             tokens = currPropTypes.(currPropNames{j});
0235             
0236             if strcmpi(tokens{1},'list')
0237                 type(j) = 1;
0238             end
0239         end
0240         
0241         % parse buffer
0242         if ~any(type)
0243             % no list types
0244             rawData = reshape(buf(offset:offset+elementCounts(i)*nProperties-1),nProperties,elementCounts(i))';
0245             offset = offset + elementCounts(i)*nProperties;
0246         else
0247             allData = cell(nProperties,1);
0248             
0249             for k = 1:nProperties
0250                 allData{k} = cell(elementCounts(i),1);
0251             end
0252             
0253             % list type
0254             for j = 1:elementCounts(i)
0255                 for k = 1:nProperties
0256                     if ~type(k)
0257                         rawData(j,k) = buf(offset);
0258                         offset = offset + 1;
0259                     else
0260                         tmp = buf(offset);
0261                         allData{k}{j} = buf(offset+(1:tmp))';
0262                         offset = offset + tmp + 1;
0263                     end
0264                 end
0265             end
0266         end
0267     else
0268         % read binary data
0269         % translate PLY data type names to MATLAB data type names
0270         listFlag = 0; % = 1 if there is a list type
0271         sameFlag = 1; % = 1 if all types are the same
0272         
0273         type = cell(1,nProperties);
0274         type2 = type;
0275         typeSize = zeros(1,nProperties);
0276         typeSize2 = typeSize;
0277         for j = 1:nProperties
0278             tokens = currPropTypes.(currPropNames{j});
0279             
0280             if ~strcmp(tokens{1}, 'list')    % non-list type
0281                 tmp = rem(find(matches(plyTypeNames,tokens{1}))-1,8)+1;
0282                 
0283                 if ~isempty(tmp)
0284                     typeSize(j) = dataTypeSizes(tmp);
0285                     type{j} = matlabTypeNames{tmp};
0286                     typeSize2(j) = 0;
0287                     type2{j} = '';
0288                     
0289                     sameFlag = sameFlag & strcmp(type{1},type{j});
0290                 else
0291                     fclose(fid);
0292                     error(['Unknown property data type, ''',tokens{1},''', in ', ...
0293                         elementNames{i},'.',currPropNames{j},'.']);
0294                 end
0295             else % list type
0296                 if length(tokens) == 3
0297                     listFlag = 1;
0298                     sameFlag = 0;
0299                     tmp = rem(find(matches(plyTypeNames,tokens{2}))-1,8)+1;
0300                     tmp2 = rem(find(matches(plyTypeNames,tokens{3}))-1,8)+1;
0301                     
0302                     if ~isempty(tmp) && ~isempty(tmp2)
0303                         typeSize(j) = dataTypeSizes(tmp);
0304                         type{j} = matlabTypeNames{tmp};
0305                         typeSize2(j) = dataTypeSizes(tmp2);
0306                         type2{j} = matlabTypeNames{tmp2};
0307                     else
0308                         fclose(fid);
0309                         error(['Unknown property data type, ''list ',tokens{2},' ',tokens{3},''', in ', ...
0310                             elementNames{i},'.',currPropNames{j},'.']);
0311                     end
0312                 else
0313                     fclose(fid);
0314                     error(['Invalid list syntax in ',elementNames{i},'.',currPropNames{j},'.']);
0315                 end
0316             end
0317         end
0318         
0319         % read file
0320         if ~listFlag
0321             if sameFlag
0322                 % no list types, all the same type (fast)
0323                 rawData = fread(fid,[nProperties,elementCounts(i)],type{1})';
0324             else
0325                 % no list types, mixed type
0326                 rawData = zeros(elementCounts(i),nProperties);
0327                 
0328                 for j = 1:elementCounts(i)
0329                     for k = 1:nProperties
0330                         rawData(j,k) = fread(fid,1,type{k});
0331                     end
0332                 end
0333             end
0334         else
0335             allData = cell(nProperties,1);
0336             
0337             for k = 1:nProperties
0338                 allData{k} = cell(elementCounts(i),1);
0339             end
0340             
0341             if nProperties == 1
0342                 bufferSize = 512;
0343                 nSkip = 4;
0344                 j = 0;
0345                 
0346                 % list type, one property (fast if lists are usually the same length)
0347                 while j < elementCounts(i)
0348                     Position = ftell(fid);
0349                     % read in BufSize count values, assuming all counts = nSkip
0350                     [buf,bufferSize] = fread(fid,bufferSize,type{1},nSkip*typeSize2(1));
0351                     miss = find(buf ~= nSkip); % find first count that is not nSkip
0352                     fseek(fid,Position + typeSize(1),-1); % seek back to after first count
0353                     
0354                     if isempty(miss) % all counts are nSkip
0355                         buf = fread(fid,[nSkip,bufferSize],[int2str(nSkip),'*',type2{1}],typeSize(1))';
0356                         fseek(fid,-typeSize(1),0); % undo last skip
0357                         
0358                         for k = 1:bufferSize
0359                             allData{1}{j+k} = buf(k,:);
0360                         end
0361                         
0362                         j = j + bufferSize;
0363                         bufferSize = floor(1.5*bufferSize);
0364                     else
0365                         if miss(1) > 1 % some counts are numSkip
0366                             Buf2 = fread(fid,[nSkip,miss(1)-1],[int2str(nSkip),'*',type2{1}],typeSize(1));
0367                             Buf2 = Buf2';
0368                             
0369                             for k = 1:miss(1)-1
0370                                 allData{1}{j+k} = Buf2(k,:);
0371                             end
0372                             
0373                             j = j + k;
0374                         end
0375 
0376                         % Alec: check if done and rewind one step
0377                         if j >= elementCounts(i)
0378                             fseek(fid,-1,0);
0379                             break;
0380                         end
0381                         
0382                         % read in the list with the missed count
0383                         nSkip = buf(miss(1));
0384                         j = j + 1;
0385                         allData{1}{j} = fread(fid,[1,nSkip],type2{1});
0386                         bufferSize = ceil(0.6*bufferSize);
0387                     end
0388                 end
0389             else
0390                 % list type(s), multiple properties (slow)
0391                 rawData = zeros(elementCounts(i),nProperties);
0392                 
0393                 for j = 1:elementCounts(i)
0394                     for k = 1:nProperties
0395                         if isempty(type2{k})
0396                             rawData(j,k) = fread(fid,1,type{k});
0397                         else
0398                             tmp = fread(fid,1,type{k});
0399                             allData{k}{j} = fread(fid,[1,tmp],type2{k});
0400                         end
0401                     end
0402                 end
0403             end
0404         end
0405     end
0406     
0407     % put data into 'elements' structure
0408     for k = 1:nProperties
0409         if (~dataFormat && ~type(k)) || (dataFormat && isempty(type2{k}))
0410             elements.(elementNames{i}).(currPropNames{k}) = rawData(:,k);
0411         else
0412             elements.(elementNames{i}).(currPropNames{k}) = allData{k};
0413         end
0414     end
0415 end
0416 
0417 clear rawData allData;
0418 fclose(fid);
0419 
0420 
0421 %% Post-processing
0422 
0423 % find the index of the element corresponding to face vertex indices
0424 possibleFacePropertyNames = ...
0425     {'vertex_indices', 'vertex_indexes', 'vertex_index', 'indices', 'indexes'};
0426 facePropertyNameIdx = find(matches(possibleFacePropertyNames, fieldnames(elements.face)));
0427 assert(isscalar(facePropertyNameIdx))
0428 
0429 % retrieve face vertex data
0430 faces = elements.face.(possibleFacePropertyNames{facePropertyNameIdx});
0431 
0432 % convert face array from 0-indexing into 1-indexing,
0433 % and attempt to convert cell array into numeric array,
0434 lengths = cellfun(@length, faces);
0435 maxLength = max(lengths);
0436 if all(maxLength == lengths)
0437     faces = cell2mat(faces)+1;
0438 else
0439     faces = cellfun(@(x) x+1, faces, 'uni', 0);
0440 end
0441 
0442 % retrieve vertex data
0443 vertices = [elements.vertex.x, elements.vertex.y, elements.vertex.z];
0444 
0445 % format output arguments
0446 varargout = formatMeshOutput(nargout, vertices, faces);
0447 if nargout == 1
0448     varargout{1}.comment = comments;
0449 end