sorting points: concave polygon

Question

I have a set of points that I'm trying to sort in ccw order or cw order from their angle. I want the points to be sorted in a way that they could form a polygon with no splits in its region or intersections. This is difficult because in most cases, it would be a concave polygon.

point centroid;
int main( int argc, char** argv )
{
  // I read a set of points into a struct point array: points[n]

  // Find centroid
  double sx = 0; double sy = 0;
  for (int i = 0; i < n; i++)
  {
    sx += points[i].x;
    sy += points[i].y;
  }
  centroid.x = sx/n;
  centroid.y = sy/n;

  // sort points using in polar order using centroid as reference
  std::qsort(&points, n, sizeof(point), polarOrder);
}

// -1 ccw, 1 cw, 0 collinear
int orientation(point a, point b, point c)
{
   double area2 = (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
   if      (area2 < 0) return -1;
   else if (area2 > 0) return +1;
   else                return  0;
}

// compare other points relative to polar angle they make with this point
// (where the polar angle is between 0 and 2pi)
int polarOrder(const void *vp1, const void *vp2)
{
  point *p1 = (point *)vp1;
  point *p2 = (point *)vp2;

  // translation
  double dx1 = p1->x - centroid.x;
  double dy1 = p1->y - centroid.y;

  double dx2 = p2->x - centroid.x;
  double dy2 = p2->y - centroid.y;

  if      (dy1 >= 0 && dy2 < 0) { return -1; }  // p1 above and p2 below
  else if (dy2 >= 0 && dy1 < 0) { return  1; }  // p1 below and p2 above
  else if (dy1 == 0 && dy2 ==0) {               // 3-collinear and horizontal
      if      (dx1 >= 0 && dx2 < 0) { return -1; }
      else if (dx2 >= 0 && dx1 < 0) { return  1; }
      else                          { return  0; }
  }
  else return -orientation(centroid,*p1,*p2);   // both above or below
}

It looks like the points are sorted accurately(pink) until they "cave" in, in which case the algorithm skips over these points then continues.. Can anyone point me into the right direction to sort the points so that they form the polygon I'm looking for? Raw Point Plot - Blue, Pink Points - Sorted

Point List: http://pastebin.com/N0Wdn2sm (You can ignore the 3rd component, since all these points lie on the same plane.)

Answer 1

Your original points don't appear form a convex polygon, so simply ordering them by angle around a fixed centroid will not necessarily result in a clean polygon. This is a non-trivial problem, you may want to research Delaunay triangulation and/or gift wrapping algorithms, although both would have to be modified because your polygon is concave. The answer here is an interesting example of a modified gift wrapping algorithm for concave polygons. There is also a C++ library called PCL that may do what you need.

But...if you really do want to do a polar sort, your sorting functions seem more complex than necessary. I would sort using atan2 first, then optimize it later once you get the result you want if necessary. Here is an example using lambda functions:

#include <algorithm>
#include <math.h>
#include <vector>

int main()
{
    struct point
    {
        double x;
        double y;
    };

    std::vector< point > points;
    point centroid;

    // fill in your data...

    auto sort_predicate = [&centroid] (const point& a, const point& b) -> bool {
        return atan2 (a.x - centroid.x, a.y - centroid.y) <
                atan2 (b.x - centroid.x, b.y - centroid.y);
    };

    std::sort (points.begin(), points.end(), sort_predicate);
}

Answer 2

The code below (sorry it's C rather than C++) sorts correctly as you wish with atan2.

The problem with your code may be that it attempts to use the included angle between the two vectors being compared. This is doomed to fail. The array is not circular. It has a first and a final element. With respect to the centroid, sorting an array requires a total polar order: a range of angles such that each point corresponds to a unique angle regardless of the other point. The angles are the total polar order, and comparing them as scalars provides the sort comparison function.

In this manner, the algorithm you proposed is guaranteed to produce a star-shaped polyline. It may oscillate wildly between different radii (...which your data do! Is this what you meant by "caved in"? If so, it's a feature of your algorithm and data, not an implementation error), and points corresponding to exactly the same angle might produce edges that coincide (lie directly on top of each other), but the edges won't cross.

I believe that your choice of centroid as the polar origin is sufficient to guarantee that connecting the ends of the polyline generated as above will produce a full star-shaped polygon, however, I don't have a proof.

Result plotted with Excel

Note you can guess from the nearly radial edges where the centroid is! This is the "star shape" I referred to above.

To illustrate this is really a star-shaped polygon, here is a zoom in to the confusing lower left corner:

If you want a polygon that is "nicer" in some sense, you will need a fancier (probably much fancier) algorithm, e.g. the Delaunay triangulation-based ones others have referred to.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct point {
  double x, y;
};

void print(FILE *f, struct point *p) {
  fprintf(f, "%f,%f\n", p->x, p->y);
}

// Return polar angle of p with respect to origin o
double to_angle(const struct point *p, const struct point *o) {
  return atan2(p->y - o->y, p->x - o->x);
}

void find_centroid(struct point *c, struct point *pts, int n_pts) {
  double x = 0, y = 0;
  for (int i = 0; i < n_pts; i++) {
    x += pts[i].x;
    y += pts[i].y;
  }
  c->x = x / n_pts;
  c->y = y / n_pts;
}

static struct point centroid[1];

int by_polar_angle(const void *va, const void *vb) {
  double theta_a = to_angle(va, centroid);
  double theta_b = to_angle(vb, centroid);
  return theta_a < theta_b ? -1 : theta_a > theta_b ? 1 : 0;
}

void sort_by_polar_angle(struct point *pts, int n_pts) {
  find_centroid(centroid, pts, n_pts);
  qsort(pts, n_pts, sizeof pts[0], by_polar_angle);
}

int main(void) {
  FILE *f = fopen("data.txt", "r");
  if (!f) return 1;
  struct point pts[10000];
  int n_pts, n_read;
  for (n_pts = 0; 
       (n_read = fscanf(f, "%lf%lf%*f", &pts[n_pts].x, &pts[n_pts].y)) != EOF;
       ++n_pts)
    if (n_read != 2) return 2;
  fclose(f);
  sort_by_polar_angle(pts, n_pts);
  for (int i = 0; i < n_pts; i++)
    print(stdout, pts + i);
  return 0;
}

Answer 3

Well, first and foremost, I see centroid declared as a local variable in main. Yet inside polarOrder you are also accessing some centroid variable.

Judging by the code you posted, that second centroid is a file-scope variable that you never initialized to any specific value. Hence the meaningless results from your comparison function.

The second strange detail in your code is that you do return -orientation(centroid,*p1,*p2) if both points are above or below. Since orientation returns -1 for CCW and +1 for CW, it should be just return orientation(centroid,*p1,*p2). Why did you feel the need to negate the result of orientation?