MATLAB Examples

Theory of convex hulls, Delaunay triangulations and Voronoi diagrams

In this section theoretical details about convex hulls, Delaunay triangulations and Voronoi Diagrams, as well as their implementation in this package are presented.

Contents

Convex hull

help convhull_nd
  Calculate the CONVex HULL of a set of points in N Dimensions.
 
  Description
      #varargout#=convhull_nd(#varargin#)
      gives the facets of the convex hull of given points.
 
  Definitions:
      (1) A set S is convex if for any two points p,q in S, the line
      segment pq lies into S.
      (2) For a subset P of Rd, the convex hull conv(P) is defined as:
          (2.1) The (unique) smallest convex set C in Rd containing P
          (P<=C)
          (2.2) The intersection of all convex sets containing P
          (2.3) The set of all convex combinations of points in P
          (2.4) The union of all simplices with vertices in P
      (3) The convex hull computation means the determination of conv(S)
      for a given finite set of n points (#points#) in Rd.
 
  Properties:
      The convex-hull operator conv() has the characteristic properties of
      a closure operator:
      (1) It is extensive, meaning that the convex hull of every set P is a
      superset of P.
      (2) It is non-decreasing, meaning that, for every two sets X and Y
      with X<=Y, the convex hull of X is a subset of the convex hull of Y.
      (3) It is idempotent, meaning that for every P, the convex hull of
      the convex hull of P is the same as the convex hull of P.
      (4) The Delaunay triangulation of a point set and its dual, the
      Voronoi diagram, are mathematically related to convex hulls: the
      Delaunay triangulation of a point set in R(#d#) can be viewed as the
      projection of a convex hull in R(#d#+1) (Brown, K. Q., 1979).
 
  Algorithm:
      To compute the convex hull (#faces#) of a set of n points (#points#)
      this algorithm uses the following scheme:
      () set #faces#, #cf# and #df# to an arbitrary simplex defined by the
         first #d#+1 points in #points#
      () for each face #k# of the convex hull created (#faces#)
          () define the simplex formed by the face #k# and the remaining
             point not belonging to the hyperplane defined by the face
          () calculate its volume
          () if the volume is negative
              () reverse the order of the last two vertices of the #k#th
                 face in #faces# to change the volume sign
              () change the sign of the plane coefficients #cf# and #df# of
                 the #k#th face in #faces#
          () end if
      () end for
      () find the relative distance of the remaining points from the center
         of the point set. 
      () sort the remaining points from the highest to the lowest relative
         distance from the center of the point set. 
      () set a remaining point vector (#pleft#) equal to the above order
      () while #pleft# is not empty
          () select the first point #i# from #pleft#. Delete it from
             #pleft#.
          () * find which of the faces of the current convex hull are
             visible from this point
          () if there are any visible faces (#visible#)
              () for each visible face #v# from #visible#
                  () find the nonvisible faces #u# connected to (i.e.
                     having one common edge with) the visible face #v#, if
                     any
                  () for each nonvisible face #k# from #u#
                      () add the boundary between the visible face #v# and
                         the nonvisible face #k# connected to the visible
                         face #v# to the horizon (#horizon#)
                  () end for
              () end for
              () delete the visible faces (#visible#) and the corresponding
                 plane coefficients #cf# and #df#
              () for each boundary #j# from #horizon#
                  () add the face defined by horizon and the selected point
                     #i# to the convex hull (#faces#)
                  () calculate the plane coefficients of the above face and
                     store them appropriately in #cf# and #df#
              () end for
              () for each face #k# from the new faces added to the convex
                 hull by the previous for-loop
                  () choose a point non-coplaner to the face
                  () define the simplex formed by the face #k# and the
                     above point
                  () calculate its volume
                  () if the volume is negative
                      () reverse the order of the last two vertices of the
                         #k#th face to change the volume sign
                      () change the sign of the plane coefficients #cf# and
                         #df# of the #k#th face
                  () end if
              () end for
              clear #horizon#
          () end if
          () if the number of iterations has reached the specified number
             #iter2del#
              () ** delete the points which are not visible by all of the
                 faces of the current convex hull
          () end if
      () end while
 
  Input arguments
      #varargin# ({#points# #iter2del#}) contains the following arguments:
          #points# ([#n# x #d#]) is a matrix containing the coordinates of
          the points the convex hull of which is to be calculated. #n# is
          the number of points and #d# is the dimensionality of the
          problem.
          #iter2del# (scalar) is the number of iterations every which
          deletion of the points inside the current convex hull occurs. It
          is recommended that it is specified (i.e. it is assigned a finite
          value) if the number of points becomes relatively large. There is
          an optimum value for this variable in general. Low values lead to
          increased number of internal point deletions which minimizes the
          number of points deleted in each repetition of the deletion
          process, and is generally more computationally costly if the
          number of points is small; the same number of iterations (one for
          each point) would have much lower computational cost. On the
          other hand, if this argument takes large values, more iterations
          take place before each deletion, which would possibly check
          internal points and would not have been done if these points were
          deleted. The final selection of the value of this parameter is
          judged according to the ratio of the execution time of step ()*
          to the execution time of step ()** shown in the algorithm above.
          If this ratio is large, #iter2del# should be decreased, and in
          the opposite case #iter2del# should be increased. If the points
          belonging to the convex hull is a significant portion of the
          initial point set, point deletion is not large, and
          consequently this parameter has little or no effect. It should be
          ignored or set to a large value. The opposite happens if the
          portion of the initial point set belonging to the convex hull
          is small.
 
  Output arguments
      #varargout# ({#faces# #cf# #df# #Vol#}) contains the following
      arguments:
          #faces# ([#nfaces# x #d#]) is a matrix containing the identities
          of the points defining each facet of the convex hull calculated.
          #nfaces# is the number of facets of the convex hull.
          #cf# ([#nfaces# x #d#]) contains the coefficients of the planes
          defined by #faces#. Each row of #cf# contains the coefficients of
          the corresponding row of #faces#. The equation of the plane #i#
          is: #cf#(#i#,:)*[x1;x2;...xd]+#df#(#i#)=0 where #d# is the
          dimension of the problem.
          #df# ([#nfaces# x 1]) contains the constant terms of the planes
          defined by #faces#. Each row of #df# contains the constant term
          of the corresponding row of #faces#. The equation of the plane
          #i# is: #cf#(#i#,:)*[x1;x2;...xd]+#df#(#i#)=0 where #d# is the
          dimension of the problem.
          #Vol# (scalar) is the volume of the convex hull.
 
  Example:
      [X,Y,Z]=meshgrid(0:1);
      points=[X(:),Y(:),Z(:)];
      points=points+0.001*rand(size(points));
      [faces,~,~,V]=convhull_nd(points)
      trisurf(faces,points(:,1),points(:,2),points(:,3))
 
  Parents (calling functions)
      delaunay_nd > convhull_nd
      voronoi_nd  > convhull_nd
 
  Children (called functions)
      convhull_nd > ismembc
      convhull_nd > plane_nd

Delaunay triangulation

help delaunay_nd
  Calculate the DELAUNAY triangulation of a set of points in N Dimensions.
 
  Description
      #DT#=delaunay_nd(#varargin#)
      gives the point identities defining the Delaunay triangulation of
      #points#.
 
  Definitions:
      (1) A triangulation of a finite point set S is a subdivision of the
      #d#-dimensional domain of the convex hull of S, whose bounded faces
      are #d#-1 entities (each defined by #d# points) and whose vertices
      are the points of S.
      (2) For a point set P in the (#d#-dimensional) Euclidean space, a
      Delaunay triangulation is a triangulation DT(P) such that no point in
      P is inside the circum-hypersphere of any simplex in DT(P). There
      exists a unique Delaunay triangulation for P, if P is a set of points
      in general position; that is, there exists no k-flat containing k + 2
      points nor a k-sphere containing k + 3 points, for 1<=k<=#d#-1
      (e.g., for a point set in R3; no three points are on a line, no four
      on a plane, no four are on a circle, and no five on a sphere).
 
  Properties:
      (1) The union of all simplices in the triangulation is the convex
      hull of the points.
      (2) The Delaunay triangulation contains O(#n#^(#d#/2)) simplices.
      (3) In the plane (#d#=2), if there are b vertices on the convex hull,
      then any triangulation of the points has at most 2#n#-2-b triangles,
      plus one exterior face.
      (4) In the plane, each vertex has on average six surrounding
      triangles.
      (5) In the plane, the Delaunay triangulation maximizes the minimum
      angle. Compared to any other triangulation of the points, the
      smallest angle in the Delaunay triangulation is at least as large as
      the smallest angle in any other. However, the Delaunay triangulation
      does not necessarily minimize the maximum angle. The Delaunay
      triangulation also does not necessarily minimize the length of the
      edges.
      (6) A circle circumscribing any Delaunay triangle does not contain
      any other input points in its interior.
      (7) If a circle passing through two of the input points does not
      contain any other of them in its interior, then the segment
      connecting the two points is an edge of a Delaunay triangulation of
      the given points.
      (8) Each simplex of the Delaunay triangulation of a set of points in
      #d#-dimensional space corresponds to a facet of convex hull of the
      projection of the points onto a (#d#+1)-dimensional paraboloid, and
      vice versa (Brown, K. Q., 1979).
      (9) The closest neighbor b to any point p is on an edge bp in the
      Delaunay triangulation.
      (10) The shortest path between two vertices, along Delaunay edges, is
      known to be no longer than $\frac{4\pi}{3\sqrt{3}} \approx 2.418$
      times the Euclidean distance between them. The plane tangent to the
      point (x0,y0,z0) of the paraboloid is used to find the point where
      this tangent plane crosses the (#d#+1)th axis (w). The facets of the
      convex hull which are visible from this point are projected
      appropriately to give the Delaunay triangulation.
 
  Input arguments
      #varargin# ({#points# #iter2del#}) contains the following arguments:
          #points# ([#n# x #d#]) is a matrix containing the coordinates of
          the points the Delaunay triangulation of which is to be
          calculated. #n# is the number of points and #d# is the
          dimensionality of the problem.
          #iter2del# (scalar) is a parameter controlling internal point
          deletion at the convex hull algorithm. Type |help convhull_nd|
          for more details.
 
  Output arguments
      #DT# ([#u# x #d#]) is a matrix containing the identities of the
      points defining each simplex of the Delaunay triangulation
      calculated. #u# is the number of simplices of the Delaunay
      triangulation.
 
  Example:
      [X,Y,Z]=meshgrid(0:1);
      points=[X(:),Y(:),Z(:)];
      points=points+0.001*rand(size(points));
      DT=delaunay_nd(points);
      % Plot the mesh
      tetramesh(DT,points)
      % Calculate the volume of the mesh
      V=0;
      for i=1:size(DT,1)
          V=V+1/6*abs(det([points(DT(i,:),:) ones(4,1)]));
      end
  
  Parents (calling functions)
  
  Children (called functions)
      delaunay_nd > convhull_nd
      delaunay_nd > convhull_nd > ismembc
      delaunay_nd > convhull_nd > plane_nd
 

Voronoi diagram

help voronoi_nd
  Calculate the vertices of the VORONOI diagram in N Dimensions.
 
  Description
      #V#=voronoi_nd(#varargin#)
      gives the coordinates of the vertices of the Voronoi diagram of
      #points#.
 
  Definitions:
      (1) A Voronoi diagram is a partitioning of a #n#-dimensional space
      into regions based on 'closeness' to points in a specific subset of
      the space. That set of points (called seeds, sites, or generators) is
      specified beforehand, and for each seed there is a corresponding
      region consisting of all points closer to that seed than to any
      other. These regions are called Voronoi cells.
      (2) The segments of the Voronoi diagram are all the points that are
      equidistant to the two nearest sites.
      (3) The Voronoi vertices (nodes) are the points equidistant to three
      (or more) sites.
      
  Properties:
      (1) The dual graph for a Voronoi diagram (in the case of a Euclidean
      space with point sites) corresponds to the Delaunay triangulation for
      the same set of points.
      (2) The closest pair of points corresponds to two adjacent cells in
      the Voronoi diagram.
      (3) Assume the setting is the Euclidean plane and a group of
      different points are given. Then two points are adjacent on the
      convex hull if and only if their Voronoi cells share an infinitely
      long side. Similarly for higher dimensions.
 
  Input arguments
      #varargin# ({#points# #iter2del#}) contains the following arguments:
          #points# ([#n# x #d#]) is a matrix containing the coordinates of
          the points the Voronoi diagram of which is to be calculated. #n#
          is the number of points and #d# is the dimensionality of the
          problem.
          #iter2del# (scalar) is a parameter controlling internal point
          deletion at the convex hull algorithm. Type |help convhull_nd|
          for more details.
 
  Output arguments
      #V# ([#u# x #d#]) is a matrix containing the coordinates of the
      Voronoi diagram vertices. #u# is the number of vertices.
 
  Example:
      points=[0 0;2 1;1 2;4 0;0 4;4 4];
      V=voronoi_nd(points);
      % Plot the initial points and the Voronoi vertices
      scatter(points(:,1),points(:,2))
      hold on
      scatter(V(:,1),V(:,2))
  
  Parents (calling functions)
  
  Children (called functions)
      voronoi_nd > convhull_nd
      voronoi_nd > convhull_nd > ismembc
      voronoi_nd > convhull_nd > plane_nd
 

Contact author

(c) 2014 by George Papazafeiropoulos
First Lieutenant, Infrastructure Engineer, Hellenic Air Force
Civil Engineer, M.Sc., Ph.D. candidate, NTUA

Email: gpapazafeiropoulos@yahoo.gr

Website: http://users.ntua.gr/gpapazaf/