Re: Arrange points so polygon is not self intersecting

In article
<1184080032.594955.39800@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>,
beeworks@xxxxxxxxxxx wrote:

On Jul 9, 7:01 pm, mek...@xxxxxxxxx wrote:
Hi all,

I was wondering if someone knows some sort of algorithm to find an
order for a set of points so that the polygon of those points is not
self-intersecting. Assume this is a 2-D plane with Euclidean geometry.

Thanks!
Dan

CLOCK Alogrithm:

1) Call your set of points S. Find the convex hull of S.

2) Pick any point P within the convex hull that is not coincident with
one of the points in S.

3) Sweep a ray R emminating from P through 360 degrees. Each time R
sweeps through a point in S, put that point next in your ordered list.

Claim: Line segments connecting adjacent points in the list created by
Step 3, plus the line segment connecting the first and last points in
the list, form a (generally star) polygon with non-intersecting sides.

This works. It is necessary to account for multiple
points on a ray. Also, it is not necessary to find the
convex hull to find one point on the hull.

Find the pivot point, p_0, on the hull by locating the
point(s) with least y-coordinate. Then set p_0 to be
that point with least x-coordinate.

Now sort the points by the angles they make with
respect to p_0. Points on the same ray are sorted by
distance from p_0.

We are almost done. If the ray sorted last has only p_0
and p_{n-1}, where n is the number of points, then the
sorted points form a simple polygon. Otherwise there
are points in the sorted list
p_{n-k}, p_{n-k+1}, p_{n-k+2}, ..., p_{n-1}
that are colinear. Reverse these points, and the sorted
list is a simple polygon. For example consider the
integer coordinate points in the closure of the triangle
{(0, 0), (2, 0), (0, 2)}.

Input points.
0: 0 2
1: 0 1
2: 1 1
3: 0 0
4: 1 0
5: 2 0

0


1 2


3 4 5


After sorting.
3: 0 0
4: 1 0
5: 2 0
2: 1 1
1: 0 1
0: 0 2

Resutlt of the simple polygon finder.
3: 0 0
4: 1 0
5: 2 0
2: 1 1
0: 0 2
1: 0 1
The polygon is simple.

We can sort the points in O(n.log(n)) time, so the
algoritm is O(n.log(n)).

After sorting the points, and before re-ordering the
final ray, we can run a Graham sweep to generate the
hull. This routine is included in the exapmle code. The
Graham sweep time is linear in n, making the hull
algoritm O(n.log(n)).


Implementation Notes

We do not find the angles made by point wrt p_0. Rather
we infer the relative angle between a and b by
computing twice the signed area of triangle{a, b, p_0}.
This is in the routine area2 which computes twice the
area
= a cross b + b cross p_0 + p_0 cross a.

Integer arithmetic is used. This computation can easily
overflow the word size used to store the point
coordinates. Even doubling the number of bits for
storing the result can overflow. For instance with four
bit signed word size and points
{(7, 7), (-7, 7), (-7, -7)}
twice the area is 196 and this is greater than
127, the largest integer representble by eight bit
signed integers. For production code, this needs to be
solved. I would probably link to a multi-precision
library and profile the time cost against the native
multiplication plus whatever other solution to overflow
is contemplated. One other solution would be to range
check the point coordinates, but that is like herding
cats.

If the C language implements the type long long int
you can define LONG_LONG_DEFINED to be 1.

The code uses a variant, qsort_r, of the standard
library routine qsort. qsort_r has an extra input
argument that it passes in an extra argument to the
user supplied compare routine. We use this feature
because in order to compare two points wrt p_0, the
compare routine needs to know the value of p_0. If
qsort_r is not available, it is necessary to make sure
that p_0 is in the scope of the compare function.

Define VERIFY_SIMPLE_POLYGON to be 1, and the code
verifies that the result of the polygon finder is a
simple polygon. The verification code is O(n^2) with a
constant term at least as large as the polygon finder.

Define COMPUTE_HULL to be 1, and the code computes the convex hull.

The compiled program attempts to read in data from
files named on the commmand line. If no files are
specified it reads stdin. It reads white space
separated strings of base 10 digits, two at a time,
taking them to be the x and y coordinates of a point.

Suppose the file data contains

0 2
1 2
2 2
0 1
1 1
2 1
0 0
1 0
2 0

These form a 3 by 3 square. The program result is
Input points.
0: 0 2
1: 1 2
2: 2 2
3: 0 1
4: 1 1
5: 2 1
6: 0 0
7: 1 0
8: 2 0

Convex hull.
0: 0 2
6: 0 0
8: 2 0
2: 2 2

Resutlt of the simple polygon finder.
6: 0 0
7: 1 0
8: 2 0
5: 2 1
4: 1 1
2: 2 2
1: 1 2
0: 0 2
3: 0 1
The polygon is simple.



This code is supplied with no guarantees.
Do not taunt Happy Fun Ball.





/*____________________________BEGIN____________________________*/
/* Connect a set of points in the plane into a simple polygon. */
/* cc <this file> -o <executable name> */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#define PMAX 1000
#define LONG_LONG_DEFINED 1
#define VERIFY_SIMPLE_POLYGON 1
#define COMPUTE_HULL 1

#if LONG_LONG_DEFINED
typedef long long int2;
#else
typedef int int2;
#endif

typedef struct
{
int x;
int y;
int tag;
} point_struct;
typedef point_struct point_t[1];

static point_t polygon[PMAX + 1];

/*
Declarations for generating a simple polygon.
*/
int read_points (FILE *, size_t, point_t *);
void find_simple_polygon (int, point_t *);
void find_lowest (int, point_t *);
void fix_up (int, point_t *);
void print_polygon (int, point_t *);
int2 point_length2 (point_t);
void point_assign (point_t, point_t);
void point_sub (point_t, point_t, point_t);
int2 area2 (point_t, point_t, point_t);
int compare (point_t, point_t, point_t);

#if VERIFY_SIMPLE_POLYGON
/*
Declarations for determining if a polygon is simple .
*/
int is_polygon_simple (int, point_t *);
int segment_intersect (point_t a, point_t b, point_t c, point_t d);
#endif

#if COMPUTE_HULL
/*
Declarations for generating the convex hull.
*/
/* Indices into the sorted points delineating the hull */
static int hull[PMAX+1];
static int top;
int run_Graham (int, point_t *);
#endif COMPUTE_HULL



int main (int argc, char *argv[])
{
int n;
int idx;
int retval = 0;
FILE *fp;

if (argc == 1)
{
n = read_points(fp, PMAX, polygon);
find_simple_polygon(n, polygon);
}
else
{
for (idx = 1 ; idx < argc; ++idx)
{
if ((fp = fopen(argv[idx], "r")) == NULL)
{
fprintf(stderr, "Unable to open %s\n", argv[idx]);
}
else
{
n = read_points(fp, PMAX, polygon);
printf("n = %d.\n", n);
find_simple_polygon(n, polygon);
fclose(fp);
}
}
}
return retval;
} /* main */

void find_simple_polygon(int n, point_t *p)
{
int idx;
int t;

printf("Input points.\n");
print_polygon(n, p);
fputc('\n', stdout);

if(n < 4) return;

/* We need a pivot point on the convex hull. */
find_lowest(n, p);

/* Sort the points on the angle they make with the pivot point, p[0]. */
qsort_r(
&p[1],
n-1,
sizeof(point_t),
p[0],
compare
);

#if COMPUTE_HULL
/*
The code here for finding the convex hull is not used
to find a simple polygon from the input points.
*/
t = run_Graham(n, p);
printf("Convex hull.\n");
for(idx = 0; idx < t; ++idx)
{
for(idx = 0; idx < t; ++idx)
printf("%5d%10d%10d\n", p[hull[idx]]->tag,
p[hull[idx]]->x,
p[hull[idx]]->y);
}
fputc('\n', stdout);
#endif COMPUTE_HULL

/* If the final ray has more than two points, rearrange them properly. */
fix_up(n, p);
printf("Resutlt of the simple polygon finder.\n");
print_polygon(n, p);

#if VERIFY_SIMPLE_POLYGON
/* Verify that the polygon is simple. */
t = is_polygon_simple(n, p);
printf("The polygon is %ssimple.\n\n", (t == 0) ? "not " : "");
#endif VERIFY_SIMPLE_POLYGON
}

/*
Find the leftmost of the points with least y coordinate.
Move it to the initial position in the array.
*/
void find_lowest(int n, point_t *p)
{
int idx;
int m;
point_t t;

m = 0;
for(idx = 1; idx < n; ++idx)
{
if( p[idx]->y < p[m]->y
|| ( p[idx]->y == p[m]->y
&& p[idx]->x < p[m]->x
)
) m = idx;
}
point_assign( t, p[0]);
point_assign(p[0], p[m]);
point_assign(p[m], t);
}

/*
Compare routine for qsort_r. Point a is ranked
lower than b when its angle wrt p is less than the
angle of b wrt p. Angle is not calculated
explicity, but rather the difference of angles is
inferred by calculating the signed area of
triangle{abp}. Points with the same angle are
sorted by distance from the pivot point, p.
*/
int compare(point_t p, point_t a, point_t b)
{
int2 z;
int2 l1;
int2 l2;
point_t r1;
point_t r2;

z = area2(p, a, b);
if (z > 0) return -1;
else if (z < 0) return 1;
else
{
point_sub(r1, a, p);
point_sub(r2, b, p);
l1 = point_length2(r1);
l2 = point_length2(r2);
return (l1 < l2) ? -1 : (l1 > l2) ? 1 : 0;
}
}

int2 area2(point_t a, point_t b, point_t c)
{
return (int2)a->x * b->y - a->y * b->x
+ b->x * c->y - b->y * c->x
+ c->x * a->y - c->y * a->x;
}

/*
If the last ray has more than one non-pivot points
it needs to have the non-pivot points reversed so
that they link up with the pivot point in the proper
order.
*/
void fix_up(int n, point_t *p)
{
int idx;
int count;
int nswap;
point_t t;

for(count = 1; area2(p[0], p[n-1], p[n-1 - count]) == 0
&& count < PMAX ; ++count);
nswap = count >> 1;
for(idx = 0; idx < nswap; ++idx)
{
point_assign(t, p[n-1 - idx]);
point_assign(p[n-1 - idx], p[n - count + idx]);
point_assign(p[n - count + idx], t);
}
}

int read_points(FILE *fp, size_t size, point_t *p)
{
int n;

for (n = 0; n < size; ++n)
{
if(fscanf(fp, "%d %d", &p[n]->x, &p[n]->y) < 2) break;
p[n]->tag = n;
}
return n;
}

int2 point_length2(point_t r)
{
int2 t1 = r->x;
int2 t2 = r->y;

return t1 * t1 + t2 * t2;
}

void point_sub(point_t r, point_t a, point_t b)
{
r->x = a->x - b->x;
r->y = a->y - b->y;
}

void point_assign(point_t to, point_t from)
{
memcpy(to, from, sizeof(point_t));
}

void print_polygon(int n, point_t *p)
{
int idx;
int fw1 = 5;
int fw2 = 10;

for(idx = 0; idx < n; ++idx)
printf("%* d%* d%* d\n", fw1, p[idx]->tag, fw2, p[idx]->x, fw2, p[idx]->y);
}

#if VERIFY_SIMPLE_POLYGON
/*
Code for determining if a polygon is simple .
*/
#define between(a,b,c) \
( (a->x != b->x) \
? (a->x <= c->x && c->x <= b->x) \
|| (a->x >= c->x && c->x >= b->x) \
: (a->y <= c->y && c->y <= b->y) \
|| (a->y >= c->y && c->y >= b->y) \
)

/*
Determine if the line segments ab and cd intersect.
*/
int segment_intersect(point_t a, point_t b, point_t c, point_t d)
{
int2 abc = area2(a, b, c);
int2 abd = area2(a, b, d);
int2 cda = area2(c, d, a);
int2 cdb = area2(c, d, b);

if( abc == 0 && between(a, b, c)
|| abd == 0 && between(a, b, d)
|| cda == 0 && between(c, d, a)
|| cdb == 0 && between(c, d, b)
) return 1;
else
return ((abc > 0) ^ (abd > 0)) && ((cda > 0) ^ (cdb > 0));
}

int is_polygon_simple(int n, point_t *p)
{
int idx;
int k;
int limit;
int retval = 1;

point_assign(p[n], p[0]);

/*
First time through we must not check edge (0,1) against edge (n-1, 0).
*/
limit = n-1;
for(idx = 0; idx < n-1 - 2; ++idx, limit = n)
{
for(k = idx + 2; k < limit; ++k)
{
if(segment_intersect(p[idx], p[idx+1], p[k], p[k+1]) == 1)
{
retval = 0;
goto egress;
}
}
}
egress:
return retval;
}
#endif VERIFY_SIMPLE_POLYGON

#if COMPUTE_HULL
/*
Run the Graham sweep on a set of angularly sorted points.
*/
int run_Graham(int n, point_t *p)
{
int idx;

top = 0;
hull[top++] = n-1;
hull[top++] = 0;

for(idx = 1; idx < n; )
{
if(area2(p[hull[top-2]], p[hull[top-1]], p[idx]) > 0) hull[top++] = idx++;
else --top;
}
return --top;
}
#endif COMPUTE_HULL
/*_____________________________END_____________________________*/

--
Michael Press
.

Relevant Pages

  • Re: Point in Polygon testing including boundaries
    ... very close to) the boundaries of the polygon. ... int pnpoly1(int npol, const double*xp, const double*yp, double x, double y) ... pnpoly = pnpoly1; ...
    (comp.graphics.algorithms)
  • Re: Pointer to pointers of structs.
    ... something in my mallocs is wrong... ... typedef struct polygon ... int main ...
    (comp.lang.c)
  • Pointer to pointers of structs.
    ... something in my mallocs is wrong... ... typedef struct polygon ... int main ...
    (comp.lang.c)
  • Re: Creating a shape from a set of points
    ... I have a bunch of points ... Are you looking for the convex hull? ... Polygon is a particularly convenient ... home dot woh dot rr dot com slash jbmatthews ...
    (comp.lang.java.programmer)
  • Re: modification of convex hull
    ... polygon X which contains all points of P. ... don't know how to reduce properly it's number of verticies. ... Perhaps you can say why you were disappointed by starting with the convex hull and reducing it. ... did you experiment with different criteria for choosing the edge to eliminate? ...
    (comp.graphics.algorithms)