【问题标题】:How to find the smallest sector of a circle that covers a polygon?如何找到覆盖多边形的圆的最小扇区?
【发布时间】:2021-01-25 15:29:13
【问题描述】:

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

假设圆的半径足够大,那么主要的问题就是找到扇形的初角和终角。

我尝试从圆心向多边形的每个角度绘制光线,并检查光线是否与多边形重叠。但是可能有不止两条光线只接触多边形。由于双精度,我不能依赖于根据方向角选择独特的光线。所以在接触光线列表中找到最小和最大角度是没有用的。除此之外,我在选择由两个终端角度创建的任何一个扇区时遇到问题,因为当由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

【问题讨论】:

  • 喜欢这些照片。您可以发布一些示例数据以供使用吗?具有半径的圆的中心点,以及定义一个或多个这些多边形的点列表。此外,你已经在周末发布了,很多人都在为期三天的假期周末发布,所以在可能的周二之前不要指望有太多点击率......
  • @Idle_Mind,感谢您的客气话。添加示例。

标签: c# algorithm geometry polygon nettopologysuite


【解决方案1】:

对于初学者,我们可以将您的数据处理为点云(多边形顶点)p[i] 以及由中心 p0 和半径 r 定义的一些圆。如果您的点云完全在圆内,则可以忽略半径。

我们可以利用atan2,但是为了避免交叉和扇区选择问题,我们不会像往常一样扩大标准笛卡尔 BBOX 计算的最小/最大界限:

  1. 计算每个点的atan2角度并将其记住在数组a[]

  2. 排序a[]

  3. 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() 函数,它只是使用示例,这是预览:

但是,当您使用多边形(线)来避免错误结果时,您有两个简单的选择:

  1. 超过 2 个点的采样线

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

  2. 转换为处理角度间隔而不是角度a[]

    所以对于每条线计算角度间隔&lt;a0,a1&gt; 与预定义的多边形缠绕规则(所以顺时针或逆时针但一致)。而不是数组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 被添加到第一个循环中,每行赢得N 点。此处预览仅使用点云方法失败的第二个多边形:

【讨论】:

  • 您能否将这个惊人的答案用于 saastn?
猜你喜欢
  • 2012-05-25
  • 1970-01-01
  • 2016-08-17
  • 1970-01-01
  • 2012-11-19
  • 1970-01-01
  • 1970-01-01
  • 2011-10-08
  • 1970-01-01
相关资源
最近更新 更多