Computational Geometry Project #4

Write an Implementation of the Hertel-Mehlhorn Algorithm

Final Report 11/10/2003 by Philip J. Porvaznik

(also in same group: Hiroki Niizato, Plamen Stoyanov, Byron Turner, Thai Vo)


Introduction

Polygon partitioning is an important preprocessing step for many geometric algorithms, because most geometric problems are simpler and faster on convex objects than on non-convex ones. We are better off whenever we can partition a non-convex object into a small number of convex pieces, because it is easier to work with the pieces independently than with the original object.

The "mother" of all polygon partitioning problems is triangulation, where the interior of the polygon is completely partitioned into triangles. Triangles are always convex and have only three sides, so any geometric operation performed on a triangle is destined to be as simple as it can be.

Because triangles all have three sides, all triangulations of polygons contain exactly the same number of pieces. Therefore, triangulation cannot be the answer if we seek a small number of convex pieces.

Convex Partitioning

A partition into triangles can be viewed as a special case of a partition into convex polygons. Because there is an optimal-time triangulation algorithm, there is an optimal-time convex partitioning algorithm. There are two goals of partitions into convex pieces: (1) partition a polygon into as few convex pieces as possible, and (2) do so as quickly as possible. Obviously, these goals conflict and one must compromise on either the number of pieces or on the time complexity (O'Rourke, page 58).

Two types of partition of a polygon P may be distinguished: a partition by diagonals or a partition by segments. The distinction is that diagonal endpoints must be vertices, while segment endpoints need only lie on polygon P. The Hertel-Mehlhorn Algorithm starts with a triangulation so its partitioning is one by diagonals not segments.

History of Optimal Partition Algorithms

Finding a convex partition optimal in the number of pieces is much more time consuming than finding a suboptimal one. The first algorithm for finding an optimal convex partition of a polygon with diagonals was due to Greene (1983). His algorithm runs in O(N4) time. This was subsequently improved by Keil (1985) to O(N3 log N) time. Both of these employ "dynamic programming," a particular algorithm technique. The problem is even more difficult if the partition may be formed with arbitrary segments since it might be necessary to employ partition segments that do not touch the polygon boundary. Nevertheless, Chazelle (1980) solved this problem in his thesis with an intricate O(N3) algorithm (also Chazelle/Dobkin 1985).


Hertel-Mehlhorn Algorithm

The Hertel-Mehlhorn heuristic for convex decomposition using diagonals is simple, efficient, and always produces no more than four times the minimum number of convex pieces. It starts with an arbitrary triangulation of the polygon and then deletes a diagonal or chord that leaves only convex pieces.

Hertel and Mehlhorn's algorithm is simply this: 

(1) Start with a triangulation of polygon P;
(2) remove an inessential diagonal;
(3) repeat.

The decision of whether a diagonal deletion will create a non-convex piece can be made locally from the chords and edges surrounding the diagonal, in constant time. A vertex in a polygon is reflex if the internal angle is greater than 180 degrees. We can remove any diagonal that does not create a reflex vertex. In some convex partition of a polygon by diagonals, call a diagonal d essential for vertex v if removal of d creates a piece that is non-convex at v. Clearly if d is essential it must be incident (i.e. connected) to v, and v must be reflex (i.e. non-convex). A diagonal that is not essential for either of its endpoints is called inessential.

There are two Lemmas associated with the H-M Algorithm. These are

Lemma (1): There can be at most two diagonals essential for any reflex vertex. This is proved by observing there can be only one essential diagonal in the half-plane (H+) to the left of the vertex (v v+), and there can be only one essential diagonal in the half-plane (H-) to the left of vertex v_v (O'Rourke, page 60, diagram).

Lemma (2): The H-M Algorithm is never worse than four-times optimal in the number of convex pieces. Since each vertex can be responsible for at most two essential diagonals, the number of essential diagonals can be no more than 2 * r, where r is the number of reflex vertices. Thus the number of convex pieces M produced by the algorithm satisfies 2 * r + 1 >= M (O'Rourke, page 60).

Hertel and Mehlhorn (1983, see paper titled "Fast Triangulation of Simple Polygons" below) found a very clean algorithm that partitions with diagonals quickly and has bounded "badness" in terms of the number of convex pieces. Clearly this algorithm results in a partition of P by diagonals into convex pieces. It can be accomplished in linear time with the use of appropriate data structures. So the only issue is how far from the optimum might it be (O'Rourke, page 60-61).

Programming Implementation

For the triangulation step, I used the O'Rourke triangulation C code (tri.c an O(N2) implementation) which I was able to compile with 0 errors, 0 warnings in MS Visual C++ (selecting "console application"). The only changes I made from the original was to #include <process.h> for the exit( ) statement; commented out the line referring to ba, ca which are unreferenced local variables in his Between( ) function; added a printf prompt when entering vertices (you enter 99 for x to exit); and added a QuitProgram function. The triangulation part does appear to work for any simple polygon, receives input of the vertices of a polygon, creates output of the diagonals, and draws the polygon and diagonals in PostScript format.

For the "removal of inessential diagonals" step, I modified O'Rourke's Triangulate( ) function. I added two arrays PolygonDiagonals[MAX_VERTICES] and PolygonVertices[MAX_VERTICES] with max #define at 50, saved the vertices while reading them in (an array copy of his linked list), saved all the diagonals while incrementing ndiagonals, and determined for each new diagonal whether it was incident to a convex (i.e. can be removed) or reflex (i.e. do not remove) vertex.

Two counters were required for the number of vertices, and number of diagonals:

int nvertices = 0;    // Total number of polygon vertices
int ndiagonals = 0; // Total number of diagonals after triangulation

The struct for the array PolygonDiagonals[MAX_VERTICES] looks like this:

typedef struct

{
   bool exist;           // whether the diagonal exists, useful at end when it has to be removed
   int  vfrom;            // the index to the array PolygonVertices of "from" point of the diagonal
   bool convexfrom;  // whether the "from" endpoint is convex
   int vto;                 // the index to the array PolygonVertices of "to" point of the diagonal
   bool convexto;      // whether the "to" endpoint is convex

}  DIAGONAL_STRUCT;

DIAGONAL_STRUCT PolygonDiagonals[MAX_VERTICES];

The struct for the array PolygonVertices[MAX_VERTICES] looks like this:

typedef struct

{
   int  xv;  // the x vertex of a point of the polygon
   int  yv;  // the y vertex of a point of the polygon

}  POLYGON_STRUCT;

POLYGON_STRUCT PolygonVertices[MAX_VERTICES];

Determine Whether a Convex or Reflex Vertex

To determine whether a diagonal endpoint was incident to a convex or reflex vertex, I used the LeftOn function, called in his InCone function as follows:

// If a is a convex vertex ...
if (LeftOn(a->v, a1->v, a0->v)) etc....

// Else a is reflex ...

Immediately after the line in Triangulate( ) -- PrintDiagonal( v1, v3 ) --  I added the following to determine whether the diagonal endpoints were convex or reflex:

tVertex a1,a0;

// get the next and previous vertex to v1
a1 = v1->next; a0 = v1->prev;

// check diagonal endpoint v1
if (LeftOn(v1->v, a1->v, a0->v)) then we have convex at v1
else non-convex (or reflex) at v1

// get the next and previous vertex to v3
a1 = v3->next; a0 = v3->prev;

// check diagonal endpoint v3
if (LeftOn(v3->v, a1->v, a0->v)) then we have convex at v3
else non-convex (or reflex) at v3

Also necessary was a function called at the end HMAlgor(  ) that goes from 0 to ndiagonals and removes the "inessential" diagonals. If the vertex at both endpoints of the diagonal was determined to be convex (i.e. less than or equal to 180 degrees), I removed that inessential diagonal by setting .exist to FALSE. If the vertex at either endpoint is reflex (i.e. greater than 180 degrees), it becomes essential and is left alone.

The removal of inessential diagonals can be done in "linear" or O(n) time. The goal is to wind up with a smaller number of diagonals forming only interior convex pieces in your original polygon. It turns out that simply determining whether the diagonal endpoints are incident (i.e. connected) to a convex vertex is not sufficient to remove all the "inessential" diagonals; there are some that remain depending on the shape of the polygon and how it was initially triangulated. My implementation could be more sophisticated if it took into account how the diagonals were related to one another.

The printing of the polygon with diagonals removed is simple: go from 0 to ndiagonals, checking if .exist = TRUE, then print that diagonal (vfrom, vto) since it wasn't removed. I duplicated the PostScript output using the two arrays containing the polygon vertex points and the diagonals. So the final output contains two polygons in PostScript format:

(1) the polygon triangulated using O'Rourke's O(N2) triangulation;

(2) the same polygon with (most?) inessential diagonals removed in O(N) time.


Examples of Input/Output and Results

Original Polygon Triangulated Polygon with Diagonals Removed
After H-M Algorithm Run
Vertices of Polygon
A Square in1.txt note: all diagonals removed out1.ps number of vertices = 4
x= 0 y= 0
x= 25 y= 0
x= 25 y= 25
x= 0 y= 25
A Convex Polygon in2.txt note: all diagonals removed out2.ps number of vertices = 8
x= 20 y= 0
x= 30 y= 10
x= 30 y= 20
x= 20 y= 30
x= 10 y= 30
x= 0 y= 20
x= 0 y= 10
x= 10 y= 0
One Reflex Vertex in3.txt note: split into two convex pieces out3.ps number of vertices = 7
x= 20 y= 0
x= 25 y= 10
x= 20 y= 20
x= 15 y= 15
x= 10 y= 20
x= 0 y= 10
x= 10 y= 0
Star Polygon in4.txt note: two diagonals removed (four convex pieces) out4.ps number of vertices = 8
x= 25 y= 0
x= 30 y= 20
x= 50 y= 25
x= 30 y= 30
x= 25 y= 50
x= 20 y= 30
x= 0 y= 25
x= 20 y= 20
"Snake" Polygon in5.txt note: only one diagonal removed (not optimal) out5.ps number of vertices = 12
x= 10 y= 0
x= 20 y= 10
x= 30 y= 0
x= 40 y= 10
x= 50 y= 0
x= 50 y= 10
x= 40 y= 20
x= 30 y= 10
x= 20 y= 20
x= 10 y= 10
x= 0 y= 20
x= 0 y= 10

Source Code and Sample Input/Output Online

triphilvaz.c -- C source of O'Rourke's Triangulation code with Hertel-Mehlhorn Implementation
macrosphilvaz.h -- header for triphilvaz.c (ADD and NEW macros)

triphilvaz.exe -- Windows 98/ME/XP executable (runs as a "console app" or DOS window)

in1.txt  |  in2.txt  |  in3.txt  |  in4.txt  |  in5.txt
(99 is used to terminate the list, x and y are entered on separate lines)

out1.ps  |  out2.ps  |  out3.ps  |  out4.ps  |  out5.ps
(these PostScript files can be read with GhostView)

note: to run this outside of Visual Studio, the best way is to use re-direction

For example, in a DOS Prompt window do

triphilvaz   <in1.txt   >out1.ps

triphilvaz   <in2.txt   >out2.ps

etc...

The < is for input file redirect, and > is for output file redirect.


Books and Online Sources Consulted:

Computational Geometry in C by Joseph O'Rourke (Cambridge Univ Press, 1998 2nd edition), chapter 1 on Polygon Triangulation, and pages 58-61 on the Hertel-Mehlhorn Algorithm and Convex Partitioning

Joseph O'Rourke's Home Page (where the Triangulation code tri.c is available for download)
http://cs.smith.edu/~orourke/

Kurt Mehlhorn's Home Page
http://www.mpi-sb.mpg.de/~mehlhorn/

Original Paper outlining the Hertel-Mehlhorn Algorithm

Stefan Hertel, Kurt Mehlhorn: Fast Triangulation of Simple Polygons. FCT (Fundamentals of Computation Theory) 1983: 207-218

Marek Karpinski (Ed.): Fundamentals of Computation Theory, Proceedings of the 1983 International FCT-Conference, Borgholm, Sweden, August 21-27, 1983. Lecture Notes in Computer Science 158 Springer-Verlag 1983, ISBN 3-540-12689-9

http://www.informatik.uni-trier.de/~ley/db/conf/fct/fct83.html

Lecture Notes in Computer Science
http://www.informatik.uni-trier.de/~ley/db/journals/lncs.html

The Computational Geometry Algorithms Library
http://www.cgal.org/index2.html


Phil Porvaznik, B.S. Computer Science (December 2003)

this page available online at www.bringyou.to/compgeom

(c) 2003 www.philvaz.com