问题描述
给定圆C:(O,r)和多边形P,我如何找到覆盖P的C的最小扇区?
假设圆的半径足够大,所以主要问题是找到扇形的初始和最终角度。
我试图从圆心向多边形的每个角度绘制射线,并检查射线是否与多边形重叠。但是可能有两条以上的光线仅接触多边形。由于双精度,我不能依靠基于它们的方向角选择唯一的射线。因此,在接触光线列表中找到最小和最大角度是没有用的。除此之外,我在选择由两个终端角度创建的扇区中遇到任何问题,因为当由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,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,52.94294 30.920592))
polyGON((52.94294 30.920592,52.45594 37.266299,59.30560 29.196711,64.12177 33.290489,58.81733 36.554277,52.94294 30.920592))
所有示例的圆心和半径:
O: (45,30)
r: 25
解决方法
对于初学者,我们可以将您的数据作为点云(多边形顶点)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,49.39652,27.550071,52.87404,30.85613,};
float p1[]={52.94294,30.920592,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,52.94294,};
float p2[]={52.94294,52.45594,37.266299,59.30560,64.12177,33.290489,58.81733,36.554277,};
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,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,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,ax,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);
}
// render backbuffer
Main->Canvas->Draw(0,bmp);
}
//---------------------------------------------------------------------------
您可以忽略void TMain::draw()
函数,它只是用法的一个示例,这是预览:
但是,为了避免错误的结果,您有多边形(线),您有两个简单的选择:
-
采样点多于2点的行
这样,角度间隙应大于点云中点之间的距离。因此,如果您对具有足够点的线进行采样,则结果将是正确的。但是,在边缘情况下,每条线选择的点数错误会导致错误的结果。另一方面,这只是添加到当前代码中的简单DDA插值。
-
转换为处理角度间隔,而不是角度
a[]
因此,对于每条线,请使用预定义的多边形缠绕规则(因此CW或CCW但一致)计算角度间隔
<a0,a1>
。而不是使用数组a[]
,而是对间隔列表进行了排序,在该列表中,您可以插入新的间隔,或者在重叠时与现有间隔合并。这种方法很安全,但是合并角度间隔并不是那么容易。如果输入数据是折线(如您的折线),则意味着每条下一行都从上一行端点开始,因此您可以忽略间隔列表,而只放大单个行,但是仍然需要正确处理不小的琐事。
[Edit1]使用第一种方法,更新后的功能如下:
void get_sector_pol(float x0,y0 circle center
// r circle radius
// p[n] point cloud
// a0,k,N=10,m=(n>>1)*N;
float ax,by,x,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,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
点。在这里预览第二个多边形,该多边形由于点云方法而失败: