PLY文件格式及其MATLAB读写操作

     PLY是一种电脑档案格式，全名为多边形档案（Polygon File Format）或 斯坦福三角形档案（Stanford Triangle Format）。 



    史丹佛大学的 The Digital Michelangelo Project计划采用PLY格式储存极高分辨率之米开朗基罗的作品"大卫"雕塑。



   该格式主要用以储存立体扫描结果的三维数值，透过多边形片面的集合描述三维物体，与其他格式相较之下这是较为简单的方法。它可以储存的资讯包含颜色、透明度、表面法向量、材质座标与资料可信度，并能对多边形的正反两面设定不同的属性。



   在档案内容的储存上PLY有两种版本，分别是纯文字（ASCII）版本与二元码（binary）版本，其差异在储存时是否以ASCII编码表示元素资讯。

档案格式

（本文并未提供完整的格式描述，以下仅介绍PLY的基本概念与格式）

每个PLY档都包含档头（header），用以设定网格模型的“元素”与“属性”，以及在档头下方接着一连串的元素“数值资料”。一般而言，网格模型的“元素”就是顶点（vertices）、面（faces），另外还可能包含有边（edges）、深度图样本（samples of range maps）与三角带（triangle strips）等元素。无论是纯文字与二元码的PLY档，档头资讯都是以ASCII编码编写，接续其后的数值资料才有编码之分。PLY档案以此行：

?

ply

开头作为PLY格式的识别。接着第二行是版本资讯，目前有三种写法：

?

format ascii 1.0

format binary_little_endian 1.0

format binary_big_endian 1.0

其中ascii, binary_little_endian, binary_big_endian是档案储存的编码方式，而1.0是遵循的标准版本（现阶段仅有PLY 1.0版）。在档头中可使用’comment’作为一行的开头以编写注解，例如：

?

comment Thisisa comment!

描述元素及属性，必须使用’element’及’property’的关键字，一般的格式为element下方接着属性列表，例如:

?

element

property

property

property

‘property’不仅定义了资料的型态，其出现顺序亦定义了资料的顺序。内定的资料形态有两种写法：一种是char uchar short ushort int uint float double,另外一种是具有位元长度的int8 uint8 int16 uint16 int32 uint32 float32 float64。 例如，描述一个包含12个顶点的物体，每个顶点使用3个单精度浮点数 (x,y,z）代表点的座标，使用3个unsigned char代表顶点颜色，颜色顺序为 (B, G, R),则档头的写法为：

?

element vertex 12

propertyfloatx

propertyfloaty

propertyfloatz

property uchar blue

property uchar green

property uchar red

其中vertex是内定的元素类型，接续的6行property描述构成vertex元素的数值字段顺序代表的意义，及其资料形态。

另一个常使用的元素是面。由于一个面是由3个以上的顶点所组成，因此使用一个“顶点列表”即可描述一个面, PLY格式使用一个特殊关键字’property list’定义之。 例如，一个具有10个面的物体，其PLY档头可能包含：

?

element face 10

property list ucharintvertex_indices

‘property list’表示该元素face的特性是由一行的顶点列表来描述。列表开头以uchar型态的数值表示列表的项目数，后面接着资料型态为int的顶点索引值（vertex_indices），顶点索引值从0开始。

最后，标头必须以此行结尾：

?

end_header

档头后接着的是元素资料（端点座标、拓朴连结等）。在ASCII格式中各个端点与面的资讯都是以独立的一行描述，而二元编码格式则连续储存这些资料，加载时须以’element’定义的元素数目以及’property’中设定的资料形态计算各笔字段的长度。

范例

一个典型的PLY档案结构分成三部分：

?

档头 (从ply开始到end_header）

顶点元素列表

面元素列表

其中的顶点元素列表一般以x y z方式排列，形态如档头所定义；而面元素列表是以下列格式表示。

?

<組成面的端點數N> <端點#1的索引> <端點#2的索引> … <端點#N的索引>

例如画出一个有4个顶点，4个面的四面体，档案内容为:

?

ply

format ascii 1.0

comment這是一個正四面體

element vertex 4

propertyfloatx

propertyfloaty

propertyfloatz

element face 4

property list ucharintvertex_index

end_header

0 3 0

2.449 -1.0 -1.414

0 -1 2.828

-2.449 -1.0 -1.414

3 0 1 3

3 0 2 1

3 0 3 2

3 1 2 3

其中1~10行是档头, 11~14行是顶点元素列表, 15~18行则是面元素列表。

其中： 0 3 0是顶点

历史

PLY格式发展于90年代中期，在史丹佛大学图学实验室的Marc Levoy教授指导下，由Greg Turk及其他成员开发出来。PLY格式受Wavefront .obj格式的启发，但改进了Obj格式所缺少的对任意属性及群组的扩充性。因此PLY格式发明了”property”及”element”这两个关键词，来概括“顶点、面、相关资讯、群组”的概念。

注意

ply文件不支持中文格式的文件名字,所以在使用过程中避免使用中文来命名。

使用MATLAB对PLY文件进行读操作

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");
%   [Tri,Pts] = PLY_READ("bunny.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...
",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

**使用MATLAB对PLY文件进行写操作**
实现将矩阵的内容转换到封装的结构体中

function [Data1,Data2]=change(A,B)
Data1.vertex.x = A(:,1);
Data1.vertex.y = A(:,2);
Data1.vertex.z = A(:,3);
Data2.vertex.x = B(:,1);
Data2.vertex.y = B(:,2);
Data2.vertex.z = B(:,3);

function ply_write ( Elements, Path, Format, Str )

%*****************************************************************************80
%
%% PLY_WRITE writes 3D data as a PLY file.
%
%  Discussion:
%
%    PLY_WRITE(DATA,FILENAME) writes the structure DATA as a binary
%    PLY file.  Every field of DATA is interpreted as an element
%    and every subfield as an element property.  Each subfield of
%    property data must either be an array or a cell array of
%    arrays.  All property data in an element must have the same
%    length.
%
%    A common PLY data structure has the following fields:
%      DATA.vertex.x = x coordinates, [Nx1] real array
%      DATA.vertex.y = y coordinates, [Nx1] real array
%      DATA.vertex.z = z coordinates, [Nx1] real array
%
%      DATA.face.vertex_indices = vertex index lists,
%         an {Mx1} cell array where each cell holds a one-
%         dimesional array (of any length) of vertex indices.
%
%    Some other common data fields:
%
%      DATA.vertex.nx = x coordinate of normal, [Nx1] real array
%      DATA.vertex.ny = y coordinate of normal, [Nx1] real array
%      DATA.vertex.nz = z coordinate of normal, [Nx1] real array
%
%      DATA.edge.vertex1 = index to a vertex, [Px1] integer array
%      DATA.edge.vertex2 = second vertex index, [Px1] integer array
%
%    Many other fields and properties can be added.  The PLY format
%    is not limited to the naming in the examples above -- they are
%    only the conventional naming.
%
%    PLY_WRITE(DATA,FILENAME,FORMAT) write the PLY with a specified
%    data format, where FORMAT is
%      "ascii"                  ASCII text data
%      "binary_little_endian"   binary data, little endian
%      "binary_big_endian"      binary data, big endian (default)
%
%    PLY_WRITE(DATA,FILENAME,FORMAT,"double") or
%
%    PLY_WRITE(DATA,FILENAME,"double") write floating-point data as
%    double precision rather than in the default single precision.
%
%  Example:
%
%    % make a cube
%    clear Data;
%    Data.vertex.x = [0;0;0;0;1;1;1;1];
%    Data.vertex.y = [0;0;1;1;0;0;1;1];
%    Data.vertex.z = [0;1;1;0;0;1;1;0];
%    Data.face.vertex_indices = {[0,1,2,3],[7,6,5,4], ...
%         [0,4,5,1],[1,5,6,2],[2,6,7,3],[3,7,4,0]};
%    ply_write(Data,"cube.ply","ascii");
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2016    Seth Billings   bug fix: "ushort" / "uint16" data type was
%                                    not included in the min/max value arrays; now
%                                    excluding "ushort" from the integer data types 
%                                    (rather than fixing the min/max arrays), since
%                                    some programs such as MeshLab do not handle the 
%                                    "ushort" data type properly anyway
%
%    01 July 2016    Seth Billings   bug fixes for max data value of "short" / "int16" 
%                                    data type and for selection of output data type
%    01 March 2007
%
%  Author:
%
%    Pascal Getreuer 2004
%    Seth Billings 2016
%
  if ( nargin < 4 )

   Str = "";

    if ( nargin < 3 )
      Format = "binary_big_endian";
    elseif ( strcmpi ( Format, "double" ) )
      Str = "double";
      Format = "binary_big_endian";
    end

  end

  [ fid, Msg ] = fopen ( Path, "wt" );

  if ( fid == -1 )
    error(Msg);
  end
%
%  Bug Fix: Seth Billings
%  - correcting max integer value for "short" / "int16" from 2^16-1 to 2^15-1
%  - excluding "ushort" from integer types since MeshLab does not
%    read "ushort" data type properly.
%
  PlyTypeNames = {"char","uchar","short","int","uint","float","double", ...
    "char8","uchar8","short16","int32","uint32","float32","double64"};
  FWriteTypeNames = {"schar","uchar","int16","int32","uint32","single","double"};
  MatlabTypeNames = {"int8","uint8","int16","int32","uint32","single","double"};
  PrintfTypeChar = {"%d","%u","%d","%d","%u","%-.6f","%-.14e"};
  IntegerDataMin =  [-128,   0, -2^15,  -2^31,   0];
  IntegerDataMax =  [ 127, 255,  2^15-1, 2^31-1, 2^32-1];
  numTypes = length(MatlabTypeNames);

  %PlyTypeNames = {"char","uchar","short","ushort","int","uint","float","double", ...
  %  "char8","uchar8","short16","ushort16","int32","uint32","float32","double64"};
  %FWriteTypeNames = {"schar","uchar","int16","uint16","int32","uint32","single","double"};
  %MatlabTypeNames = {"int8","uint8","int16","uint16","int32","uint32","single","double"};
  %PrintfTypeChar = {"%d","%u","%d","%u","%d","%u","%-.6f","%-.14e"};  
  %IntegerDataMin = [-128, 0, -2^15, -2^31, 0];
  %IntegerDataMax = [127, 255, 2^16-1, 2^31-1, 2^32-1];
%
%  Write the PLY header.
%

  fprintf(fid,"ply
format %s 1.0
comment created by MATLAB ply_write
",Format);
  ElementNames = fieldnames(Elements);
  NumElements = length(ElementNames);
  Data = cell(NumElements,1);

  for i = 1 : NumElements

    eval(["tmp=isa(Elements.",ElementNames{i},",""struct"");"]);

    if ( tmp )
      eval(["PropertyNames{i}=fieldnames(Elements.",ElementNames{i},");"]);
    else
      PropertyNames{i} = [];
    end

    if ( ~isempty(PropertyNames{i}) )
      eval(["Data{i}{1}=Elements.",ElementNames{i},".",PropertyNames{i}{1},";"]);
      ElementCount(i) = prod(size(Data{i}{1}));
      Type{i} = zeros(length(PropertyNames{i}),1);
    else
      ElementCount(i) = 0;
    end

    fprintf(fid,"element %s %u
",ElementNames{i},ElementCount(i));

    for j = 1 : length(PropertyNames{i})

      eval(["Data{i}{j}=Elements.",ElementNames{i},".",PropertyNames{i}{j},";"]);

      if ( ElementCount(i) ~= prod(size(Data{i}{j})) )
        fclose(fid);
        error("All property data in an element must have the same length.");
      end

      if ( iscell(Data{i}{j}) )
        Type{i}(j) = numTypes + 1;  % Seth Billings
        %Type{i}(j) = 9;
        Data{i}{j} = Data{i}{j}{1};
      end

      for k = 1 : length(MatlabTypeNames)
        if ( isa(Data{i}{j},MatlabTypeNames{k}) )
          Type{i}(j) = Type{i}(j) + k;
          break;
        end
      end

      if ( ~rem(Type{i}(j),numTypes+1) ) % Seth Billings
      %if ( ~rem(Type{i}(j),9) )
        fclose(fid);
        error("Unsupported data structure.");
      end
%
%  Try to convert float data to integer data
%
%  Array data.
%
      if ( Type{i}(j) <= numTypes)  % Seth Billings
      %if ( Type{i}(j) <= 8 )
        if any(strcmp({"single","double"},MatlabTypeNames{Type{i}(j)}))
          if ~any(floor(Data{i}{j}) ~= Data{i}{j})  % data is integer
            MinValue = min(min(Data{i}{j}));
            MaxValue = max(max(Data{i}{j}));

            % choose smallest possible integer data format            
            % Bug Fix: Seth Billings
            tmp = min(find(MinValue >= IntegerDataMin & MaxValue <= IntegerDataMax));
            %tmp = max(min(find(MinValue >= IntegerDataMin)),min(find(MaxValue <= IntegerDataMax)));

            if ~isempty(tmp)
              Type{i}(j) = tmp;
            end
          end
        end
      else        % cell array data
        eval(["Data{i}{j}=Elements.",ElementNames{i},".",PropertyNames{i}{j},";"]);
        tmp = 1;

        for k = 1:prod(size(Data{i}{j}))
          tmp = tmp & all(floor(Data{i}{j}{k}) == Data{i}{j}{k});
        end

        if tmp  % data is integer
          MinValue = inf;
          MaxValue = -inf;

          for k = 1:prod(size(Data{i}{j}))
            MinValue = min(MinValue,min(Data{i}{j}{k}));
            MaxValue = max(MaxValue,max(Data{i}{j}{k}));
          end
%
%  choose smallest possible integer data format
%  Bug Fix: Seth Billings
%
          tmp = min(find(MinValue >= IntegerDataMin & MaxValue <= IntegerDataMax));
          %tmp = max(min(find(MinValue >= IntegerDataMin)),min(find(MaxValue <= IntegerDataMax)));

          if ( ~ isempty ( tmp ) )
            Type{i}(j) = tmp + numTypes + 1;
            %Type{i}(j) = tmp + 9;
          end

        end
      end

      % convert double to single if specified
      if rem(Type{i}(j),numTypes+1) == numTypes & ~strcmpi(Str,"double") % Seth Billings
      %if rem(Type{i}(j),9) == 8 & ~strcmpi(Str,"double")
        Type{i}(j) = Type{i}(j) - 1;
      end

      if Type{i}(j) <= numTypes  % Seth Billings
      %if Type{i}(j) <= 8
        fprintf(fid,"property %s %s
",PlyTypeNames{Type{i}(j)},PropertyNames{i}{j});
      else
        fprintf(fid,"property list uchar %s %s
",PlyTypeNames{Type{i}(j)-(numTypes+1)},PropertyNames{i}{j}); % Seth Billings
        %fprintf(fid,"property list uchar %s %s
",PlyTypeNames{Type{i}(j)-9},PropertyNames{i}{j});
      end
    end
  end

  fprintf(fid,"end_header
");

  switch Format
    case "ascii"
      Format = 0;
    case "binary_little_endian"
      fclose(fid);
      fid = fopen(Path,"a","ieee-le");
      Format = 1;
    case "binary_big_endian"
      fclose(fid);
      fid = fopen(Path,"a","ieee-be");
      Format = 2;
  end

  for i = 1 : NumElements
    if ~isempty(PropertyNames{i})
      if ~Format          % write ASCII data
        for k = 1:ElementCount(i)
          for j = 1:length(PropertyNames{i})
            if Type{i}(j) <= numTypes  % Seth Billings
            %if Type{i}(j) <= 8
              fprintf(fid,[PrintfTypeChar{Type{i}(j)}," "],Data{i}{j}(k));
            else
              fprintf(fid,"%u%s ",length(Data{i}{j}{k}),sprintf([" ",PrintfTypeChar{Type{i}(j)-(numTypes+1)}],Data{i}{j}{k})); % Seth Billings
              %fprintf(fid,"%u%s ",length(Data{i}{j}{k}),sprintf([" ",PrintfTypeChar{Type{i}(j)-9}],Data{i}{j}{k}));
            end
          end

          fprintf(fid,"
");
        end
      else            % write binary data
        if all(Type{i} <= numTypes) & all(Type{i} == Type{i}(1))  % Seth Billings
        %if all(Type{i} <= 8) & all(Type{i} == Type{i}(1))
          % property data without list types (fast)
          tmp = zeros(length(PropertyNames{i}),ElementCount(i));

          for j = 1:length(PropertyNames{i})
            tmp(j,:) = Data{i}{j}(:)";
          end

          fwrite(fid,tmp,FWriteTypeNames{Type{i}(j)});
        elseif all(Type{i} > numTypes) % Seth Billings
        %elseif all(Type{i} > 8)
        % only list types
          Type{i} = Type{i} - (numTypes+1); % Seth Billings
          %Type{i} = Type{i} - 9;

          if length(PropertyNames{i}) == 1
           % only one list property
            tmp = FWriteTypeNames{Type{i}(1)};

            for k = 1:ElementCount(i)
              fwrite(fid,length(Data{i}{1}{k}),"uchar");
              fwrite(fid,Data{i}{1}{k},tmp);
            end
          else
           % multiple list properties
            for k = 1:ElementCount(i)
              for j = 1:length(PropertyNames{i})
                fwrite(fid,length(Data{i}{j}{k}),"uchar");
                fwrite(fid,Data{i}{j}{k},FWriteTypeNames{Type{i}(j)});
              end
            end
          end
        else
        % mixed type
          for k = 1:ElementCount(i)
            for j = 1:length(PropertyNames{i})
              if Type{i}(j) <= numTypes % Seth Billings
              %if Type{i}(j) <= 8
                fwrite(fid,Data{i}{j}(k),FWriteTypeNames{Type{i}(j)});
              else
                fwrite(fid,length(Data{i}{j}{k}),"uchar");
                fwrite(fid,Data{i}{j}{k},FWriteTypeNames{Type{i}(j)-(numTypes+1)}); % Seth Billings
                %fwrite(fid,Data{i}{j}{k},FWriteTypeNames{Type{i}(j)-9});
              end
            end
          end
        end
      end
    end
  end

  fclose ( fid );

  return
end