c# - 如何找到覆盖多边形的圆的最小扇区？

Question

给定圆 C: (O, r) 和多边形 P，如何找到覆盖 P 的 C 的最小扇区？

假设圆的半径足够大，那么主要的问题是找到扇形的初始角度和最终角度。
我试图从圆心向多边形的每个角度绘制光线，并检查光线是否与多边形重叠。但是可能有两条以上的射线只接触多边形。由于 double ，我不能依赖于根据它们的方向角选择独特的光线。所以在接触光线列表中找到最小和最大角度是没有用的。除此之外，我在选择由两个终端角创建的扇区时遇到问题，因为当由 atan2 计算时，初始角度可能大于最终角度。 .

那么找到这样一个部门的正确方法是什么？
编辑:
三个示例多边形(WKT 格式):

POLYGON((52.87404 30.85613, 42.55699 28.46292, 41.54373 24.319989, 53.57623 21.300564, 62.94891 28.46292, 49.39652 27.550071, 52.87404 30.85613))
POLYGON((52.94294 30.920592, 42.55699 28.46292, 43.61965 35.545578, 55.85037 34.862696, 59.12524 36.621547, 47.68664 39.877048, 35.69973 36.198265, 37.30512 29.196711, 31.09762 28.46292, 41.54373 24.319989, 53.57623 21.300564, 62.94891 28.46292, 49.39652 27.550071, 52.94294 30.920592))
POLYGON((52.94294 30.920592, 42.55699 28.46292, 43.61965 35.545578, 52.45594 37.266299, 59.30560 29.196711, 64.12177 33.290489, 58.81733 36.554277, 47.68664 39.877048, 35.69973 36.198265, 37.30512 29.196711, 31.09762 28.46292, 41.54373 24.319989, 53.57623 21.300564, 62.94891 28.46292, 49.39652 27.550071, 52.94294 30.920592))

所有示例的圆心和半径:

O: (45, 30)
r: 25

Answer 1

首先，我们可以将您的数据处理为点云(多边形顶点)p[i]和一些由中心定义的圆 p0和半径 r .如果您的点云完全在圆内，您可以忽略半径。
我们可以利用 atan2然而，为了避免交叉和扇区选择的问题，我们不会像往常一样扩大标准笛卡尔 BBOX 计算的最小/最大边界:

计算 atan2每个点的角度并记住它在数组中 a[]

排序 a[]

在 a[] 中找到后续角度之间的最大距离
不要忘记角度差可以是|Pi|顶部所以如果它更多你需要+/- 2*PI .也处理a[]作为循环缓冲区。

这是我简单的 C++/VCL 尝试:

//---------------------------------------------------------------------------
float p0[]={52.87404,30.856130,42.55699,28.46292,41.54373,24.319989,53.57623,21.300564,62.94891,28.46292,49.39652,27.550071,52.87404,30.85613,};
float p1[]={52.94294,30.920592,42.55699,28.46292,43.61965,35.545578,55.85037,34.862696,59.12524,36.621547,47.68664,39.877048,35.69973,36.198265,37.30512,29.196711,31.09762,28.46292,41.54373,24.319989,53.57623,21.300564,62.94891,28.46292,49.39652,27.550071,52.94294,30.920592,};
float p2[]={52.94294,30.920592,42.55699,28.46292,43.61965,35.545578,52.45594,37.266299,59.30560,29.196711,64.12177,33.290489,58.81733,36.554277,47.68664,39.877048,35.69973,36.198265,37.30512,29.196711,31.09762,28.46292,41.54373,24.319989,53.57623,21.300564,62.94891,28.46292,49.39652,27.550071,52.94294,30.920592,};
float x0=45.0,y0=30.0,R=25.0;
//---------------------------------------------------------------------------
template <class T> void sort_asc_bubble(T *a,int n)
    {
    int i,e; T a0,a1;
    for (e=1;e;n--)                                     // loop until no swap occurs
     for (e=0,a0=a[0],a1=a[1],i=1;i<n;a0=a1,i++,a1=a[i])// proces unsorted part of array
      if (a0>a1)                                        // condition if swap needed
      { a[i-1]=a1; a[i]=a0; a1=a0; e=1; }               // swap and allow to process array again
    }
//---------------------------------------------------------------------------
void get_sector(float x0,float y0,float r,float *p,int n,float &a0,float &a1)
    {
    // x0,y0    circle center
    // r        circle radius
    // p[n]     polyline vertexes
    // a0,a1    output angle range a0<=a1
    int i,j,m=n>>1;
    float x,y,*a;
    a=new float[m];
    // process points and compute angles
    for (j=0,i=0;i<n;j++)
        {
        x=p[i]-x0; i++;
        y=p[i]-y0; i++;
        a[j]=atan2(y,x);
        }
    // sort by angle
    sort_asc_bubble(a,m);
    // get max distance
    a0=a[m-1]; a1=a[0]; x=a1-a0;
    while (x<-M_PI) x+=2.0*M_PI;
    while (x>+M_PI) x-=2.0*M_PI;
    if (x<0.0) x=-x;
    for (j=1;j<m;j++)
        {
        y=a[j]-a[j-1];
        while (y<-M_PI) y+=2.0*M_PI;
        while (y>+M_PI) y-=2.0*M_PI;
        if (y<0.0) y=-y;
        if (y>x){ a0=a[j-1]; a1=a[j]; x=y; }
        }
    }
//---------------------------------------------------------------------------
void TMain::draw()
    {
    int i,n;
    float x,y,r,*p,a0=0.0,a1=0.0;
    float ax,ay,bx,by;
    float zoom=7.0;
    p=p0; n=sizeof(p0)/sizeof(p0[0]);
//  p=p1; n=sizeof(p1)/sizeof(p1[0]);
//  p=p2; n=sizeof(p2)/sizeof(p2[0]);

    get_sector(x0,y0,R,p,n,a0,a1);

    // clear buffer
    bmp->Canvas->Brush->Color=clBlack;
    bmp->Canvas->FillRect(TRect(0,0,xs,ys));

    // circle
    x=x0; y=y0; r=R;
    ax=x+R*cos(a0);
    ay=y+R*sin(a0);
    bx=x+R*cos(a1);
    by=y+R*sin(a1);
    x*=zoom; y*=zoom; r*=zoom;
    ax*=zoom; ay*=zoom;
    bx*=zoom; by*=zoom;
    bmp->Canvas->Pen->Color=clBlue;
    bmp->Canvas->Brush->Color=TColor(0x00101010);
    bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
    bmp->Canvas->Pen->Color=clAqua;
    bmp->Canvas->Brush->Color=TColor(0x00202020);
    bmp->Canvas->Pie(x-r,y-r,x+r,y+r,ax,ay,bx,by);

    // PCL
    r=2.0;
    bmp->Canvas->Pen->Color=clAqua;
    bmp->Canvas->Brush->Color=clAqua;
    for (i=0;i<n;)
        {
        x=p[i]; i++;
        y=p[i]; i++;
        x*=zoom; y*=zoom;
        bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
        }

    // render backbuffer
    Main->Canvas->Draw(0,0,bmp);
    }
//---------------------------------------------------------------------------

您可以忽略 void TMain::draw()函数只是使用示例，这是预览:

但是，由于您有多边形(线)以避免错误结果，您有两个简单的选择:

超过 2 个点的采样线
这样角度间隙应该大于点云中点之间的距离。因此，如果您对具有足够点的线进行采样，结果将是正确的。然而，错误地选择每条线的点数将导致边缘情况下的错误结果。另一方面，实现这只是添加到当前代码中的简单 DDA 插值。

转换为处理角度间隔而不是角度 a[]
所以对于每条线计算角度间隔<a0,a1>具有预定义的多边形缠绕规则(因此 CW 或 CCW 但一致)。而不是数组 a[]您将订购间隔列表，您可以在其中插入新间隔或在重叠时与现有间隔合并。这种方法是安全的，但合并角度间隔并不那么容易。如果输入数据是折线(就像你的一样)，这意味着下一行从前一行端点开始，所以你可以忽略间隔列表，只放大单个，但你仍然需要正确处理放大，这不是微不足道的。

[Edit1] 使用第一种方法更新后的函数如下所示:

void get_sector_pol(float x0,float y0,float r,float *p,int n,float &a0,float &a1)
    {
    // x0,y0    circle center
    // r        circle radius
    // p[n]     point cloud
    // a0,a1    output angle range a0<=a1
    int i,j,k,N=10,m=(n>>1)*N;
    float ax,ay,bx,by,x,y,dx,dy,*a,_N=1.0/N;
    a=new float[m];
    // process points and compute angles
    bx=p[n-2]-x0; i++;
    by=p[n-1]-y0; i++;
    for (j=0,i=0;i<n;)
        {
        ax=bx; ay=by;
        bx=p[i]-x0; i++;
        by=p[i]-y0; i++;
        dx=_N*(bx-ax); x=ax;
        dy=_N*(by-ay); y=ay;
        for (k=0;k<N;k++,x+=dx,y+=dy,j++) a[j]=atan2(y,x);
        }
    // sort by angle
    sort_asc_bubble(a,m);
    // get max distance
    a0=a[m-1]; a1=a[0]; x=a1-a0;
    while (x<-M_PI) x+=2.0*M_PI;
    while (x>+M_PI) x-=2.0*M_PI;
    if (x<0.0) x=-x;
    for (j=1;j<m;j++)
        {
        y=a[j]-a[j-1];
        while (y<-M_PI) y+=2.0*M_PI;
        while (y>+M_PI) y-=2.0*M_PI;
        if (y<0.0) y=-y;
        if (y>x){ a0=a[j-1]; a1=a[j]; x=y; }
        }
    }

如您所见，它几乎相同，只是简单 DDA 添加到第一个循环 win N每行点数。这是仅使用点云方法失败的第二个多边形的预览: