c#algorithm sorting geometry computational-geometry

Sorting points by polar angle

I have a problem with sorting points by the angle they create with the X axis. The points look like this

Here is the code I have:

public static List<Point> SortPoints(List<Point> points)
        {
            List<Point> result = new List<Point>();
            List<KeyValuePair<Point, double>> valuePairs = new List<KeyValuePair<Point, double>>();
            foreach (var point in points)
            {
                valuePairs.Add(new KeyValuePair<Point, double>(point, Math.Atan2(point.Y, point.X)));
            }
            valuePairs = valuePairs.OrderByDescending(x => x.Value).ToList();
            foreach (var valuepair in valuePairs)
            {
                result.Add(valuepair.Key);
            }
            return result;
        }

It should sort points in the way, they close up. It works for most points, but it doesn't for some of them. It crashes mostly on these fragments:

Is my thinking correct for that kind of problem or do I miss something? I am still new to geometry in programming.

Solution

Ok I got it working sorting by atan2 angle by using 2 centers only (no need for 3 as I can detect the problem zones directly from the first center alone). This is the algorithm (shape must not self intersect and angle around selected centers must be monotonic !!!):

compute BBOX (x0,y0,x1,y1) for your data

this is needed to properly compute correct center locations for atan2 usage
compute angle by atan2 for each point using BBOX center as center

the center should be inside your shape so center of BBOX is the obvious choice. However as Yves Daoust pointed out this will not work for arbitrary concave shapes only for those shapes and centers where the angle is monotonic.
sort your points by this angle
detect problematic zones

simply in problematic zones the consequent points after the sort has almost the same angle so just threshold that.
compute atan2 angle for each problem zone with different center

again center must be inside ... and should be shifted in any of the multiple of 90 degrees angle from original center. I chose shift toward x0 by 20% of shape x size. The bigger the shift the more ang difference the problem zones will get.
sort the problem zones by new angle
reverse problem zone order after sort if needed

the shifted center might cause the angle direction reversal in comparison to original angles. So after sort if you compute distance between point before zone and zone first and last point if the last point of zone is closer it means you need to reverse the zone points order.

Here preview of output:

Here C++/OpenGL/VCL code for this:

//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
#include "data.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
const int n2=sizeof(pnt)/sizeof(pnt[0]);    // size of pnt[]
const int n=n2>>1;                          // point in pnt[]
int idx[n];                                 // index sort
int ix=1055;                                // just debug X cursor on mouse wheel
//---------------------------------------------------------------------------
void compute()
    {
    int i,j,k,i0,i1,e,m;
    float a,x0,y0,x1,y1,x,y;
    float ang[n];   // atan2 angles per point
    int flag[n];    // 0 or problem zone ix per point
    const float deg=M_PI/180.0;
    const float thr=0.02*deg;
    // shuffle input data for debug as its already ordered
    for (i=0;i<n2;)
        {
        j=Random(n2)&0xFFFFFFFE;
        a=pnt[i]; pnt[i]=pnt[j]; pnt[j]=a; i++; j++;
        a=pnt[i]; pnt[i]=pnt[j]; pnt[j]=a; i++; j++;
        }
    // init index sort table
    for (i=0;i<n;i++) idx[i]=i+i;
    // compute BBOX of data
    x0=x1=pnt[0];
    y0=y1=pnt[1];
    for (i=0;i<n2;)
        {
        x=pnt[i]; i++;
        y=pnt[i]; i++;
        if (x0>x) x0=x;
        if (x1<x) x1=x;
        if (y0>y) y0=y;
        if (y1<y) y1=y;
        }
    // compute atan2 for center set to center of BBOX
    x=0.5*(x0+x1);
    y=0.5*(y0+y1);
    for (i=0,j=0;i<n;i++,j+=2)
        {
        a=atan2(pnt[j+1]-y,pnt[j+0]-x);                 // this might need handling edge cases of atan2
        ang[i]=a;
        }
    // index (bubble) sort idx[] asc by ang[]
    for (e=1,j=n;e;j--)                                 // loop until no swap occurs
     for (e=0,i=1;i<j;i++)                              // process unsorted part of array
      if (ang[idx[i-1]>>1]>ang[idx[i]>>1])              // condition if swap needed
      { e=idx[i-1]; idx[i-1]=idx[i]; idx[i]=e; e=1; }   // swap and allow to process array again
    // detect/label problematic zones m = number of zones +1
    for (i=0;i<n;i++) flag[i]=0;
    for (e=0,j=1,i=1;i<n;i++)
     if (fabs(ang[idx[i]>>1]-ang[idx[i-1]>>1])<thr)
      { flag[idx[i]>>1]=j; flag[idx[i-1]>>1]=j; e=1; }
      else if (e){ e=0; j++; }
    if (e) j++; m=j;
    // compute atan2 for center shifted toward x0
    // so it still inside but not too close to (0,0)
    // so there is some ang diference on problematic zones
    x=x0+(0.3*(x1-x0));
    y=0.5*(y0+y1);
    for (i=0,j=0;i<n;i++,j+=2)
     if (flag[i])                                       // only for problematic zones no need to recompute finished parts
        {
        a=atan2(pnt[j+1]-y,pnt[j+0]-x);                 // this might need handling edge cases of atan2
        ang[i]=a;
        }
    // loop through problematic zones
    for (k=0;k<n;)
        {
        for (;k<n;k++) if (flag[idx[k]>>1])                     // zone start
            {
            i0=i1=k;
            for (;k<n;k++) if (flag[idx[k]>>1]) i1=k;           // zone end
              else break;
            // index (bubble) sort idx[] asc by ang[]
            if (i0!=i1)
             for (e=1,j=i1-i0+1;e;j--)                          // loop until no swap occurs
              for (e=0,i=i0+1;i<i0+j;i++)                       // process unsorted part of array
               if (ang[idx[i-1]>>1]>ang[idx[i]>>1])             // condition if swap needed
                { e=idx[i-1]; idx[i-1]=idx[i]; idx[i]=e; e=1; } // swap and allow to process array again
            // different center for atan2 might reverse the ang order
            // so test if start or end of zone is closer to the point before it
            j=i0-1; if (j<0) j=n-1;                             // 45 deg is never at start or end of data so this should never happen
            x=pnt[idx[j]+0]-pnt[idx[i0]+0];
            y=pnt[idx[j]+1]-pnt[idx[i0]+1];
            a=(x*x)+(y*y);
            x=pnt[idx[j]+0]-pnt[idx[i1]+0];
            y=pnt[idx[j]+1]-pnt[idx[i1]+1];
            x=(x*x)+(y*y);
            // reverse if not in correct order
            if (x<a) for (;i0<i1;i0++,i1--)
             { j=idx[i0]; idx[i0]=idx[i1]; idx[i1]=j; }
            }
        }
    }
//---------------------------------------------------------------------------
void gl_draw()
    {
    int i,j;
    float a,da=1.0/float(n-1),x,y,r;
    glClear(GL_COLOR_BUFFER_BIT);
    glDisable(GL_DEPTH_TEST);
    glDisable(GL_TEXTURE_2D);

    // set view to 2D
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glScalef(1.0/120.0,1.0/120.0,1.0);

    // render points from list
    glBegin(GL_POINTS);
//  glBegin(GL_LINE_STRIP);
//  glBegin(GL_TRIANGLE_FAN);
    glColor3f(0.0,0.0,0.0);
    glVertex2f(0.0,0.0);
    for (a=0.0,i=0;i<n;i++,a+=da)
        {
        glColor3f(a,a,a);
        glVertex2fv(pnt+idx[i]);
        }
    glEnd();

    // render debug index (on mouse wheel)
    x=pnt[idx[ix]+0];
    y=pnt[idx[ix]+1];
    r=5.0;
    glBegin(GL_LINES);
    glColor3f(0.0,1.0,0.0);
    glVertex2f(x-r,y-r);
    glVertex2f(x+r,y+r);
    glVertex2f(x-r,y+r);
    glVertex2f(x+r,y-r);
    glEnd();

    glFinish();
    SwapBuffers(hdc);
    Form1->Caption=ix;
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    // Init of program
    gl_init(Handle);    // init OpenGL
    compute();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    // Exit of program
    gl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    // repaint
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    // resize
    gl_resize(ClientWidth,ClientHeight);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormMouseWheel(TObject *Sender, TShiftState Shift, int WheelDelta, TPoint &MousePos, bool &Handled)
    {
    Handled=true;
    int dix=1; if (Shift.Contains(ssShift)) dix=10;
    if (WheelDelta>0){ ix+=dix; if (ix>=n) ix=  0; }
    if (WheelDelta<0){ ix-=dix; if (ix< 0) ix=n-1; }
    gl_draw();
    }
//---------------------------------------------------------------------------

Just ignore the VCL and OpenGL stuff. The only important stuff here is the function compute() which works as described above. The global variables just above it. I used index sort however so the points order is not changed at all instead idx[i] holds the index of i-th point in the input data. I wanted to keep this as simple as I could so no dynamic allocations nor containers or funny stuff... I used your data as input in form:

float pnt[]=
    {
    -98.622,0.4532042,
    -98.622,1.64291,
    -98.612097,3.0877569,
    ...
    -98.618994,-3.2649391,
    -98.6260115,-1.9205891
    };

The gl_simple.h I used for OpenGL can be found here:

complete GL+GLSL+VAO/VBO C++ example

I tested this by using the cursor (green cross on image) and mouse wheel through the whole shape looking for jumping back and forward ... In previous versions I got the flag[] array as global and rendered different colors for the problem zones so I can debug directly them and not looking for whole shape again and again ... Also I used bubble sort ... in case you got huge data use quick sort instead... however this data you provided is computed instantly on my old computer so I didn't bother to add recursion.

For arbitrary concave non self-intersecting shapes you need to use different approach like for example connected component analysis:

for each point compute its 2 nearest neighbor points

this is very slow and should be speed-ed up by using spatial sorting of points to lower the complexity.

in case of very non uniform sampling the two closest points should not lie on or near the same direction!!! So dot between their unit direction should be as far from +1.0 as it can.
select any start point and add it to output
select one of its yet unused neighbors and add it to output
set link between selected and its predecessing point as used.
loop #4 until you get to the starting point again