按极角排序点

Question

我在按点与 X 轴形成的角度对点进行排序时遇到问题。这些点看起来像这样

这是我的代码：

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;
        }

它应该在途中对点进行排序，它们靠拢。它适用于大多数点，但不适用于其中一些点。它主要在这些片段上崩溃：

对于这类问题我的想法是正确的还是我遗漏了什么？我对编程中的几何还是个新手。

Answer 1

好的，我仅使用 2 个中心就可以按 atan2 角度进行排序（不需要 3 个，因为我可以仅从第一个中心直接检测到问题区域）。这是算法（形状不能自相交并且围绕 selected 中心的角度必须是单调的！！！）：

为您的数据计算 BBOX (x0,y0,x1,y1)

这是正确计算 atan2 使用的正确中心位置所必需的
使用 BBOX 中心作为中心
为每个点计算角度 atan2
中心应该在你的形状内，所以 BBOX 的中心是显而易见的选择。然而，正如 Yves Daoust 指出的那样，这不适用于任意凹形，仅适用于角度单调的那些形状和中心。
按这个角度排序你的点
检测有问题的区域

只是在有问题的区域中，排序后的结果点几乎具有相同的角度，因此只需设置阈值即可。
计算每个具有不同中心的问题区域的atan2角度

同样，中心必须在...内部，并且应该从原始中心移动 90 度角的任意倍数。我选择向 x0 移动形状 x 大小的 20%。偏移越大，问题区域的角度差异就越大。
按新角度对问题区域进行排序
如果需要，排序后反转问题区域顺序

与原始角度相比，移动中心可能会导致角度方向反转。因此，在排序之后，如果您计算区域之前的点与区域第一个点和最后一个点之间的距离，如果区域的最后一个点更近，则意味着您需要反转区域点顺序。

此处输出预览：

此处为 C++/OpenGL/VCL 代码：

//---------------------------------------------------------------------------
#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();
    }
//---------------------------------------------------------------------------

忽略 VCL 和 OpenGL 的东西。这里唯一重要的东西是函数 compute() ，它的工作原理如上所述。它上面的全局变量。但是我使用了索引排序，所以点顺序根本没有改变，而是 idx[i] 保存输入数据中第 i 个点的索引。我想让它尽可能简单，所以没有动态分配，也没有容器或有趣的东西......我使用你的数据作为表单的输入：

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

我用于OpenGL的gl_simple.h可以在这里找到：

我通过使用光标（图像上的绿色十字）和鼠标滚轮在整个形状中寻找来回跳跃来测试这个......在以前的版本中我得到了 flag[] 数组作为全局和渲染问题区域有不同的颜色，所以我可以直接调试它们，而不是一次又一次地寻找整个形状......我还使用了冒泡排序......如果你有大量数据，请改用快速排序......但是你提供的这个数据在我的旧计算机上立即计算，所以我没有费心去添加递归。

对于任意凹形非自相交形状，您需要使用不同的方法，例如连通分量分析：

为每个点计算它的 2 个最近邻点

这非常慢，应该通过使用点的空间排序来降低复杂性来加快速度。

在采样非常不均匀的情况下，两个最近的点不应位于或靠近同一方向！！！所以它们的单位方向之间的点应该尽可能远离+1.0。
select任意起点并将其添加到输出
select 其尚未使用的邻居之一并将其添加到输出
设置 link 在 selected 和它使用的前导点之间。
循环 #4 直到再次到达起点

Answer 2

C# 翻译

public static List<Point> SortPoints(List<Point> points)
    {
        double[] pnt = new double[points.Count * 2];
        for (int z = 0, ii = 0; z < points.Count; z++, ii += 2)
        {
            pnt[ii] = points[z].X;
            pnt[ii + 1] = points[z].Y;
        }
        int n2 = pnt.Length;
        int n = n2 >> 1;
        int[] idx = new int[n];
        int i, j, k, i0, i1, e, m;
        double a, x0, y0, x1, y1, x, y;
        double[] ang = new double[n];
        int[] flag = new int[n];
        const double deg = Math.PI / 180.0;
        const double thr = 0.02 * deg;
        for (i = 0; i < n; i++) idx[i] = i + i;
        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;
        }
        x = 0.5 * (x0 + x1);
        y = 0.5 * (y0 + y1);
        for (i = 0, j = 0; i < n; i++, j += 2)
        {
            a = Math.Atan2(pnt[j + 1] - y, pnt[j + 0] - x);
            ang[i] = a;
        }
        for (e = 1, j = n; e != 0; j--)                                 // loop until jo swap occurs
            for (e = 0, i = 1; i < j; i++)                              // proces 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; }
        for (i = 0; i < n; i++) flag[i] = 0;
        for (e = 0, j = 1, i = 1; i < n; i++)
            if (Math.Abs(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 != 0) { e = 0; j++; }
        if (e != 0) j++; m = j;
        x = x0 + (0.3 * (x1 - x0));
        y = 0.5 * (y0 + y1);
        for (i = 0, j = 0; i < n; i++, j += 2)
            if (flag[i] != 0)                                      // only for problematic zones no need to recompute finished parts
            {
                a = Math.Atan2(pnt[j + 1] - y, pnt[j + 0] - x);                 // this might need handling edge cases of atan2
                ang[i] = a;
            }
        for (k = 0; k < n;)
        {
            for (; k < n; k++) if (flag[idx[k] >> 1] != 0)                     // zone start
                {
                    i0 = i1 = k;
                    for (; k < n; k++) if (flag[idx[k] >> 1] != 0) i1 = k;//           // zone end
                        else break;
                    // index (bubble) sort idx[] asc by ang[]
                    if (i0 != i1)
                        for (e = 1, j = i1 - i0 + 1; e > 0; j--)                          // loop until jo swap occurs
                            for (e = 0, i = i0 + 1; i < i0 + j; i++)                       // proces 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; }
                }
        }
        List<Point> result = new List<Point>();
        for (int h = 0; h < pnt.Length - 1; h += 2)
        {
            result.Add(new Point(pnt[h], pnt[h + 1]));
        }
        return result;
    }

按极角排序点

Sorting points by polar angle

c#

sorting

algorithm

geometry

computational-geometry