在2D点集中find洞？

我有一套2d points 。 它们是标准笛卡尔网格系统上的X,Y coordinates （在本例中是UTM zone ）。 我需要find那个点的孔，最好有一些设置find这些孔的algorithm的灵敏度的能力。 通常这些点集非常密集，但有些可能不那么密集。

它们是什么，如果有帮助的话，就是农田里的土壤被取样为农业中的人们显然有用的各种性质的地方。 有时在这些点样品有巨大的岩石或沼泽的地方或充满小湖泊和池塘。

从点集他们希望我find大致定义这些洞的凹多边形。

我已经写出了find外部凹形边界多边形的algorithm，并允许它们设置algorithm形成的多边形的粗糙度或平滑度的灵敏度。 之后，他们希望find洞，并把这些洞作为一个凹多边形，我猜在某些情况下可能是凸的，但这将是边缘情况不规范。

问题是我所听到的关于这个问题的唯一论文是天文学家所做的，他们希望在空间中find大量空洞，我从来没有能够find他们的论文中的实际algorithm，除了粗略的概念性描述外，

我已经尝试过谷歌和各种学术论文search等，但迄今为止我还没有发现很多有用的东西。 这让我想知道是否有这样一个问题的名称，我不知道，所以我正在寻找错误的东西或什么？

所以在这个漫长的解释之后，我的问题是：有没有人知道我应该寻找什么来find这个文件，最好用明确的algorithm，或者有人知道一个algorithm，他们可以指向我？

任何能够帮助我解决这个问题的东西都是非常有用的，并且非常值得赞赏，即使只是正确的search短语，将会产生潜在的algorithm或论文。

这将实施的语言将是C＃，但从Mathematica包到MATLAB or ASM, C, C++, Python, Java or MathCAD等的任何示例都可以，只要示例中没有一些调用去FindTheHole等等FindTheHole没有定义或者是专有的实现软件，例如MathCAD例子通常有很多。

下面是实际点集的两个例子，一个是稠密的，一个是稀疏的，我们需要find它们的区域：

Answer 1

那么像这样的一些位图+vector方法呢？

1. 获得点云区域覆盖的边界框

   如果它不是已知的，就这样做。 通过所有点应该是简单的O(N)循环。

2. 创build该地区的map[N][N]

   这是一个容易进行数据密度计算的“位图”。 只需以简单的比例创buildarea(x,y) -> map[i][j]投影。 网格大小N也是输出的精度 ， 必须大于平均点距离！ 因此map[][]每个单元格至less包含一个点（如果不在孔区域中）。

3. 计算map[][]每个单元的数据密度

   简单的方法就是将map[][].cnt map[i][j].cnt++ （点的计数器） zero然后通过简单的O(N)循环进行计算，其中map[i][j].cnt++用于所有points(x,y)

4. 创build未使用区域列表(map[][].cnt==0)或(map[][].cnt<=treshold)

   我通过水平和垂直线来简化

5. 分割输出

   只需将同一个孔的线条（相交的…向量方法）组合起来，也可以通过填充（位图方法）在子弹＃4中完成，

6. 多边形输出

   取同一个洞/组的H，V线的所有边缘点并创build多边形（对它们进行sorting，使它们的连接不会相交）。 有很多关于这个的库，algorithm和源代码。

我的这种方法的源代码：

 void main_compute(int N) { // cell storage for density computation struct _cell { double x0,x1,y0,y1; // bounding area of points inside cell int cnt; // points inside cell _cell(){}; _cell(_cell& a){ *this=a; }; ~_cell(){}; _cell* operator = (const _cell *a) { *this=*a; return this; }; /*_cell* operator = (const _cell &a) { ...copy... return this; };*/ }; // line storage for hole area struct _line { double x0,y0,x1,y1; // line edge points int id; // id of hole for segmentation/polygonize int i0,i1,j0,j1; // index in map[][] _line(){}; _line(_line& a){ *this=a; }; ~_line(){}; _line* operator = (const _line *a) { *this=*a; return this; }; /*_line* operator = (const _line &a) { ...copy... return this; };*/ }; int i,j,k,M=N*N; // M = max N^2 but usualy is much much less so dynamic list will be better double mx,my; // scale to map _cell *m; // cell ptr glview2D::_pnt *p; // point ptr double x0,x1,y0,y1; // used area (bounding box) _cell **map=NULL; // cell grid _line *lin=NULL; // temp line list for hole segmentation int lins=0; // actual usage/size of lin[M] // scan point cloud for bounding box (if it is known then skip it) p=&view.pnt[0]; x0=p->p[0]; x1=x0; y0=p->p[1]; y1=y0; for (i=0;i<view.pnt.num;i++) { p=&view.pnt[i]; if (x0>p->p[0]) x0=p->p[0]; if (x1<p->p[0]) x1=p->p[0]; if (y0>p->p[1]) y0=p->p[1]; if (y1<p->p[1]) y1=p->p[1]; } // compute scale for coordinate to map index conversion mx=double(N)/(x1-x0); // add avoidance of division by zero if empty point cloud !!! my=double(N)/(y1-y0); // dynamic allocation of map[N][N],lin[M] lin=new _line[M]; map=new _cell*[N]; for (i=0;i<N;i++) map[i]=new _cell[N]; // reset map[N][N] for (i=0;i<N;i++) for (j=0;j<N;j++) map[i][j].cnt=0; // compute point cloud density for (k=0;k<view.pnt.num;k++) { p=&view.pnt[k]; i=double((p->p[0]-x0)*mx); if (i<0) i=0; if (i>=N) i=N-1; j=double((p->p[1]-y0)*my); if (j<0) j=0; if (j>=N) j=N-1; m=&map[i][j]; if (!m->cnt) { m->x0=p->p[0]; m->x1=p->p[0]; m->y0=p->p[1]; m->y1=p->p[1]; } if (m->cnt<0x7FFFFFFF) m->cnt++; // avoid overflow if (m->x0>p->p[0]) m->x0=p->p[0]; if (m->x1<p->p[0]) m->x1=p->p[0]; if (m->y0>p->p[1]) m->y0=p->p[1]; if (m->y1<p->p[1]) m->y1=p->p[1]; } // find holes (map[i][j].cnt==0) or (map[i][j].cnt<=treshold) // and create lin[] list of H,V lines covering holes for (j=0;j<N;j++) // search lines { for (i=0;i<N;) { int i0,i1; for (;i<N;i++) if (map[i][j].cnt==0) break; i0=i-1; // find start of hole for (;i<N;i++) if (map[i][j].cnt!=0) break; i1=i; // find end of hole if (i0< 0) continue; // skip bad circumstances (edges or no hole found) if (i1>=N) continue; if (map[i0][j].cnt==0) continue; if (map[i1][j].cnt==0) continue; _line l; l.i0=i0; l.x0=map[i0][j].x1; l.i1=i1; l.x1=map[i1][j].x0; l.j0=j ; l.y0=0.25*(map[i0][j].y0+map[i0][j].y1+map[i1][j].y0+map[i1][j].y1); l.j1=j ; l.y1=l.y0; lin[lins]=l; lins++; } } for (i=0;i<N;i++) // search columns { for (j=0;j<N;) { int j0,j1; for (;j<N;j++) if (map[i][j].cnt==0) break; j0=j-1; // find start of hole for (;j<N;j++) if (map[i][j].cnt!=0) break; j1=j; // find end of hole if (j0< 0) continue; // skip bad circumstances (edges or no hole found) if (j1>=N) continue; if (map[i][j0].cnt==0) continue; if (map[i][j1].cnt==0) continue; _line l; l.i0=i ; l.y0=map[i][j0].y1; l.i1=i ; l.y1=map[i][j1].y0; l.j0=j0; l.x0=0.25*(map[i][j0].x0+map[i][j0].x1+map[i][j1].x0+map[i][j1].x1); l.j1=j1; l.x1=l.x0; lin[lins]=l; lins++; } } // segmentate lin[] ... group lines of the same hole together by lin[].id // segmentation based on vector lines data // you can also segmentate the map[][] directly as bitmap during hole detection for (i=0;i<lins;i++) lin[i].id=i; // all lines are separate for (;;) // join what you can { int e=0,i0,i1; _line *a,*b; for (a=lin,i=0;i<lins;i++,a++) { for (b=a,j=i;j<lins;j++,b++) if (a->id!=b->id) { // do 2D lines a,b intersect ? double xx0,yy0,xx1,yy1; double kx0,ky0,dx0,dy0,t0; double kx1,ky1,dx1,dy1,t1; double x0=a->x0,y0=a->y0; double x1=a->x1,y1=a->y1; double x2=b->x0,y2=b->y0; double x3=b->x1,y3=b->y1; // discart lines with non intersecting bound rectangles double a0,a1,b0,b1; if (x0<x1) { a0=x0; a1=x1; } else { a0=x1; a1=x0; } if (x2<x3) { b0=x2; b1=x3; } else { b0=x3; b1=x2; } if (a1<b0) continue; if (a0>b1) continue; if (y0<y1) { a0=y0; a1=y1; } else { a0=y1; a1=y0; } if (y2<y3) { b0=y2; b1=y3; } else { b0=y3; b1=y2; } if (a1<b0) continue; if (a0>b1) continue; // compute intersection kx0=x0; ky0=y0; dx0=x1-x0; dy0=y1-y0; kx1=x2; ky1=y2; dx1=x3-x2; dy1=y3-y2; t1=divide(dx0*(ky0-ky1)+dy0*(kx1-kx0),(dx0*dy1)-(dx1*dy0)); xx1=kx1+(dx1*t1); yy1=ky1+(dy1*t1); if (fabs(dx0)>=fabs(dy0)) t0=divide(kx1-kx0+(dx1*t1),dx0); else t0=divide(ky1-ky0+(dy1*t1),dy0); xx0=kx0+(dx0*t0); yy0=ky0+(dy0*t0); // check if intersection exists if (fabs(xx1-xx0)>1e-6) continue; if (fabs(yy1-yy0)>1e-6) continue; if ((t0<0.0)||(t0>1.0)) continue; if ((t1<0.0)||(t1>1.0)) continue; // if yes ... intersection point = xx0,yy0 e=1; break; } if (e) break; // join found ... stop searching } if (!e) break; // no join found ... stop segmentation i0=a->id; // joid ids ... rename i1 to i0 i1=b->id; for (a=lin,i=0;i<lins;i++,a++) if (a->id==i1) a->id=i0; } // visualize lin[] for (i=0;i<lins;i++) { glview2D::_lin l; l.p0.p[0]=lin[i].x0; l.p0.p[1]=lin[i].y0; l.p1.p[0]=lin[i].x1; l.p1.p[1]=lin[i].y1; // l.col=0x0000FF00; l.col=(lin[i].id*0x00D00C10A)+0x00800000; // color is any function of ID view.lin.add(l); } // dynamic deallocation of map[N][N],lin[M] for (i=0;i<N;i++) delete[] map[i]; delete[] map; delete[] lin; } //--------------------------------------------------------------------------- 

只要忽略我的glview2D东西（这是我的几何渲染引擎）

• view.pnt[]是你的点的dynamic列表（随机生成）
• view.lin[]是dynamic列表输出H，V行仅用于可视化
• lin[]是你的线条输出

这是输出：

我懒得添加polygonize现在你可以看到分割工作（着色）。 如果你还需要帮助polygonize然后评论我，但我认为这不应该有任何问题。

复杂度估计取决于整体的空洞覆盖率

但是对于大部分代码来说，它是O(N) ，对于空洞search/分割~O((M^2)+(U^2))其中：

• N是点数
• M是地图网格大小
• U是H，V线数依赖于孔…
• M << N, U << M*M

正如你可以看到上面的图像上的3783点30x30网格，我花了近9ms设置

[Edit1]玩vectorpolygonize一点点

对于简单的孔是好的，但对于更复杂的那些还有一些小腿

[编辑2]终于得到了一点时间，所以这里是它：

这是一个更简单的类/多边形search更令人愉快/可pipe理的forms：

 //--------------------------------------------------------------------------- class holes { public: int xs,ys,n; // cell grid x,y - size and points count int **map; // points density map[xs][ys] // i=(x-x0)*g2l; x=x0+(i*l2g); // j=(y-y0)*g2l; y=y0+(j*l2g); double mg2l,ml2g; // scale to/from global/map space (x,y) <-> map[i][j] double x0,x1,y0,y1; // used area (bounding box) struct _line { int id; // id of hole for segmentation/polygonize int i0,i1,j0,j1; // index in map[][] _line(){}; _line(_line& a){ *this=a; }; ~_line(){}; _line* operator = (const _line *a) { *this=*a; return this; }; /*_line* operator = (const _line &a) { ...copy... return this; };*/ }; List<_line> lin; int lin_i0; // start index for perimeter lines (smaller indexes are the H,V lines inside hole) struct _point { int i,j; // index in map[][] int p0,p1; // previous next point int used; _point(){}; _point(_point& a){ *this=a; }; ~_point(){}; _point* operator = (const _point *a) { *this=*a; return this; }; /*_point* operator = (const _point &a) { ...copy... return this; };*/ }; List<_point> pnt; // class init and internal stuff holes() { xs=0; ys=0; n=0; map=NULL; mg2l=1.0; ml2g=1.0; x0=0.0; y0=0.0; x1=0.0; y1=0.0; lin_i0=0; }; holes(holes& a){ *this=a; }; ~holes() { _free(); }; holes* operator = (const holes *a) { *this=*a; return this; }; holes* operator = (const holes &a) { xs=0; ys=0; n=an; map=NULL; mg2l=a.mg2l; x0=a.x0; x1=a.x1; ml2g=a.ml2g; y0=a.y0; y1=a.y1; _alloc(a.xs,a.ys); for (int i=0;i<xs;i++) for (int j=0;j<ys;j++) map[i][j]=a.map[i][j]; return this; } void _free() { if (map) { for (int i=0;i<xs;i++) if (map[i]) delete[] map[i]; delete[] map; } xs=0; ys=0; } void _alloc(int _xs,int _ys) { int i=0; _free(); xs=_xs; ys=_ys; map=new int*[xs]; if (map) for (i=0;i<xs;i++) { map[i]=new int[ys]; if (map[i]==NULL) { i=-1; break; } } else i=-1; if (i<0) _free(); } // scann boundary box interface void scann_beg(); void scann_pnt(double x,double y); void scann_end(); // dynamic allocations void cell_size(double sz); // compute/allocate grid from grid cell size = sz x sz // scann holes interface void holes_beg(); void holes_pnt(double x,double y); void holes_end(); // global(x,y) <- local map[i][j] + half cell offset inline void l2g(double &x,double &y,int i,int j) { x=x0+((double(i)+0.5)*ml2g); y=y0+((double(j)+0.5)*ml2g); } // local map[i][j] <- global(x,y) inline void g2l(int &i,int &j,double x,double y) { i= double((x-x0) *mg2l); j= double((y-y0) *mg2l); } }; //--------------------------------------------------------------------------- void holes::scann_beg() { x0=0.0; y0=0.0; x1=0.0; y1=0.0; n=0; } //--------------------------------------------------------------------------- void holes::scann_pnt(double x,double y) { if (!n) { x0=x; y0=y; x1=x; y1=y; } if (n<0x7FFFFFFF) n++; // avoid overflow if (x0>x) x0=x; if (x1<x) x1=x; if (y0>y) y0=y; if (y1<y) y1=y; } //--------------------------------------------------------------------------- void holes::scann_end() { } //--------------------------------------------------------------------------- void holes::cell_size(double sz) { int x,y; if (sz<1e-6) sz=1e-6; x=ceil((x1-x0)/sz); y=ceil((y1-y0)/sz); _alloc(x,y); ml2g=sz; mg2l=1.0/sz; } //--------------------------------------------------------------------------- void holes::holes_beg() { int i,j; for (i=0;i<xs;i++) for (j=0;j<ys;j++) map[i][j]=0; } //--------------------------------------------------------------------------- void holes::holes_pnt(double x,double y) { int i,j; g2l(i,j,x,y); if ((i>=0)&&(i<xs)) if ((j>=0)&&(j<ys)) if (map[i][j]<0x7FFFFFFF) map[i][j]++; // avoid overflow } //--------------------------------------------------------------------------- void holes::holes_end() { int i,j,e,i0,i1; List<int> ix; // hole lines start/stop indexes for speed up the polygonization _line *a,*b,l; _point *aa,*bb,p; lin.num=0; lin_i0=0;// clear lines ix.num=0; // clear indexes // find holes (map[i][j].cnt==0) or (map[i][j].cnt<=treshold) // and create lin[] list of H,V lines covering holes for (j=0;j<ys;j++) // search lines for (i=0;i<xs;) { int i0,i1; for (;i<xs;i++) if (map[i][j]==0) break; i0=i-1; // find start of hole for (;i<xs;i++) if (map[i][j]!=0) break; i1=i; // find end of hole if (i0< 0) continue; // skip bad circumstances (edges or no hole found) if (i1>=xs) continue; if (map[i0][j]==0) continue; if (map[i1][j]==0) continue; l.i0=i0; l.i1=i1; l.j0=j ; l.j1=j ; l.id=-1; lin.add(l); } for (i=0;i<xs;i++) // search columns for (j=0;j<ys;) { int j0,j1; for (;j<ys;j++) if (map[i][j]==0) break; j0=j-1; // find start of hole for (;j<ys;j++) if (map[i][j]!=0) break; j1=j ; // find end of hole if (j0< 0) continue; // skip bad circumstances (edges or no hole found) if (j1>=ys) continue; if (map[i][j0]==0) continue; if (map[i][j1]==0) continue; l.i0=i ; l.i1=i ; l.j0=j0; l.j1=j1; l.id=-1; lin.add(l); } // segmentate lin[] ... group lines of the same hole together by lin[].id // segmentation based on vector lines data // you can also segmentate the map[][] directly as bitmap during hole detection for (i=0;i<lin.num;i++) lin[i].id=i; // all lines are separate for (;;) // join what you can { for (e=0,a=lin.dat,i=0;i<lin.num;i++,a++) { for (b=a,j=i;j<lin.num;j++,b++) if (a->id!=b->id) { // if a,b not intersecting or neighbouring if (a->i0>b->i1) continue; if (b->i0>a->i1) continue; if (a->j0>b->j1) continue; if (b->j0>a->j1) continue; // if they do mark e for join groups e=1; break; } if (e) break; // join found ... stop searching } if (!e) break; // no join found ... stop segmentation i0=a->id; // joid ids ... rename i1 to i0 i1=b->id; for (a=lin.dat,i=0;i<lin.num;i++,a++) if (a->id==i1) a->id=i0; } // sort lin[] by id for (e=1;e;) for (e=0,a=&lin[0],b=&lin[1],i=1;i<lin.num;i++,a++,b++) if (a->id>b->id) { l=*a; *a=*b; *b=l; e=1; } // re id lin[] and prepare start/stop indexes for (i0=-1,i1=-1,a=&lin[0],i=0;i<lin.num;i++,a++) if (a->id==i1) a->id=i0; else { i0++; i1=a->id; a->id=i0; ix.add(i); } ix.add(lin.num); // polygonize lin_i0=lin.num; for (j=1;j<ix.num;j++) // process hole { i0=ix[j-1]; i1=ix[j]; // create border pnt[] list (unique points only) pnt.num=0; p.used=0; p.p0=-1; p.p1=-1; for (a=&lin[i0],i=i0;i<i1;i++,a++) { pi=a->i0; pj=a->j0; map[pi][pj]=0; for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++) if ((aa->i==pi)&&(aa->j==pj)) { e=-1; break; } if (e>=0) pnt.add(p); pi=a->i1; pj=a->j1; map[pi][pj]=0; for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++) if ((aa->i==pi)&&(aa->j==pj)) { e=-1; break; } if (e>=0) pnt.add(p); } // mark not border points for (aa=&pnt[0],i=0;i<pnt.num;i++,aa++) if (!aa->used) // ignore marked points if ((aa->i>0)&&(aa->i<xs-1)) // ignore map[][] border points if ((aa->j>0)&&(aa->j<ys-1)) { // ignore if any non hole cell around if (map[aa->i-1][aa->j-1]>0) continue; if (map[aa->i-1][aa->j ]>0) continue; if (map[aa->i-1][aa->j+1]>0) continue; if (map[aa->i ][aa->j-1]>0) continue; if (map[aa->i ][aa->j+1]>0) continue; if (map[aa->i+1][aa->j-1]>0) continue; if (map[aa->i+1][aa->j ]>0) continue; if (map[aa->i+1][aa->j+1]>0) continue; aa->used=1; } // delete marked points for (aa=&pnt[0],e=0,i=0;i<pnt.num;i++,aa++) if (!aa->used) { pnt[e]=*aa; e++; } pnt.num=e; // connect neighbouring points distance=1 for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) if (aa->used<2) for (i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++) if (bb->used<2) { i=aa->i-bb->i; if (i<0) i=-i; e =i; i=aa->j-bb->j; if (i<0) i=-i; e+=i; if (e!=1) continue; aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1; bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0; } // try to connect neighbouring points distance=sqrt(2) for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) if (aa->used<2) for (i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++) if (bb->used<2) if ((aa->p0!=i1)&&(aa->p1!=i1)) if ((bb->p0!=i0)&&(bb->p1!=i0)) { if ((aa->used)&&(aa->p0==bb->p0)) continue; // avoid small closed loops i=aa->i-bb->i; if (i<0) i=-i; e =i*i; i=aa->j-bb->j; if (i<0) i=-i; e+=i*i; if (e!=2) continue; aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1; bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0; } // try to connect to closest point int ii,dd; for (i0= 0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) if (aa->used<2) { for (ii=-1,i1=i0+1,bb=&pnt[i1];i1<pnt.num;i1++,bb++) if (bb->used<2) if ((aa->p0!=i1)&&(aa->p1!=i1)) if ((bb->p0!=i0)&&(bb->p1!=i0)) { i=aa->i-bb->i; if (i<0) i=-i; e =i*i; i=aa->j-bb->j; if (i<0) i=-i; e+=i*i; if ((ii<0)||(e<dd)) { ii=i1; dd=e; } } if (ii<0) continue; i1=ii; bb=&pnt[i1]; aa->used++; if (aa->p0<0) aa->p0=i1; else aa->p1=i1; bb->used++; if (bb->p0<0) bb->p0=i0; else bb->p1=i0; } // add connected points to lin[] ... this is hole perimeter !!! // lines are 2 x duplicated so some additional code for sort the order of line swill be good idea l.id=lin[ix[j-1]].id; for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) { l.i0=aa->i; l.j0=aa->j; // [edit3] this avoid duplicating lines if (aa->p0>i0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } if (aa->p1>i0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } //if (aa->p0>=0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } //if (aa->p1>=0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } } } } //--------------------------------------------------------------------------- 

你只需要将我的List<T>模板replace为std::list或其他（我不能共享的模板）。 这是一个dynamic的一维数组T …

• List<int> x; 和int x[];是一样的int x[];
• x.add(); 将空项添加到x
• x.add(a); 将一个项目添加到x
• x.reset()清除数组
• x.allocate(size)预先分配空间以避免重新分配缓慢的运行
• x.num是项目中x [] …使用的大小中的项目数

在原来的代码只有静态数组，所以如果你感到困惑，请检查它。

现在如何使用它：

 h.scann_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.scann_pnt(p[0],p[1]); } h.scann_end(); h.cell_size(2.5); h.holes_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.holes_pnt(p[0],p[1]); } h.holes_end(); 

其中view.pnt[]是input点列表，并在其中： view.pnt[i].p0.p[ 2 ]= { x,y }

输出在h.lin[]和lin_i0 ，其中：

• h.lin[i] i= < 0,lin_i0 )是内部H，V行
• h.lin[i] i= < lin_i0,h.lin.num )是周长

外围线没有sorting和重复两次，所以只要命令他们，并删除重复（太懒了）。 里面的lin[]是id .. id的孔0,1,2,3,...线属于和i,j坐标里面的地图。 所以为了适当的输出到你的世界坐标做这样的事情：

 int i,j; holes h; // holes class double *p; // input point list ptr h.scann_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.scann_pnt(p[0],p[1]); } h.scann_end(); h.cell_size(2.5); h.holes_beg(); for (i=0;i<view.pnt.num;i++) { p=view.pnt[i].p0.p; h.holes_pnt(p[0],p[1]); } h.holes_end(); DWORD coltab[]= { 0x000000FF, 0x0000FF00, 0x00FF0000, 0x0000FFFF, 0x00FFFF00, 0x00FF00FF, 0x00FFFFFF, 0x00000088, 0x00008800, 0x00880000, 0x00008888, 0x00888800, 0x00880088, 0x00888888, }; for (i=0;i<h.lin.num;i++) // draw lin[] { glview2D::_lin a; holes::_line *b=&h.lin[i]; h.l2g(a.p0.p[0],a.p0.p[1],b->i0,b->j0); h.l2g(a.p1.p[0],a.p1.p[1],b->i1,b->j1); if (i<h.lin_i0) // H,V lines inside hole(b->id) .. gray [edit3] was <= which is wrong and miss-color first perimeter line { a.col=0x00808080; } else{ // hole(b->id) perimeter lines ... each hole different collor if ((b->id>=0)&&(b->id<14)) a.col=coltab[b->id]; if (b->id==-1) a.col=0x00FFFFFF; // special debug lines if (b->id==-2) a.col=0x00AA8040; // special debug lines } view.lin.add(a); // here draw your line or add it to your polygon instead } 

• 我的view.lin[]有成员： p0,p1,它们的点是view.pnt[] ， col是colour

当孔太小(diameter < 3 cells)时，我只看到一个问题，否则就OK

重新排列周界线

要做到这一点，而不是：

  /* add connected points to lin[] ... this is hole perimeter !!! // lines are 2 x duplicated so some additional code for sort the order of line swill be good idea l.id=lin[ix[j-1]].id; for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) { l.i0=aa->i; l.j0=aa->j; // [edit3] this avoid duplicating lines if (aa->p0>i0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } if (aa->p1>i0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } //if (aa->p0>=0) { bb=&pnt[aa->p0]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } //if (aa->p1>=0) { bb=&pnt[aa->p1]; l.i1=bb->i; l.j1=bb->j; lin.add(l); } } */ 

做这个：

  // add connected points to lin[] ... this is hole perimeter !!! l.id=lin[ix[j-1]].id; // add index of points instead points int lin_i1=lin.num; for (i0=0,aa=&pnt[i0];i0<pnt.num;i0++,aa++) { l.i0=i0; if (aa->p0>i0) { l.i1=aa->p0; lin.add(l); } if (aa->p1>i0) { l.i1=aa->p1; lin.add(l); } } // reorder perimeter lines for (i0=lin_i1,a=&lin[i0];i0<lin.num-1;i0++,a++) for (i1=i0+1 ,b=&lin[i1];i1<lin.num ;i1++,b++) { if (a->i1==b->i0) { a++; l=*a; *a=*b; *b=l; a--; break; } if (a->i1==b->i1) { a++; l=*a; *a=*b; *b=l; i=a->i0; a->i0=a->i1; a->i1=i; a--; break; } } // convert point indexes to points for (i0=lin_i1,a=&lin[i0];i0<lin.num;i0++,a++) { bb=&pnt[a->i0]; a->i0=bb->i; a->j0=bb->j; bb=&pnt[a->i1]; a->i1=bb->i; a->j1=bb->j; } 

[编辑5]如何在holes::holes_end多边形

作为这个input，你需要所有的H，V行 lin[]分段/分组/按孔sorting的列表和密度映射map[][] 。

1. 通过所有的孔循环

   2. 循环处理的所有H，V行

      创build所有唯一行端点列表pnt[] （不重复）。 所以每行有两个端点，看每个点是否已经在列表中。 如果没有添加它，否则忽略它。

   3. 从列表中删除所有非边界点

      因此，通过在密度map[][]查看4个邻居来移除所有与非孔洞区域无接触的点

   4. 做点连接组件分析

      1. set used=0; p0=-1; p1=-1; used=0; p0=-1; p1=-1; 对于pnt[]列表中的所有点

      2. 连接distance=1点

         used<2来循环所有点pnt[] ，这意味着它们还没有被完全使用，并且对于每个这样的点再次searchpnt[]以用于distance = 1另一个点。 这意味着它是它的4邻居，应该连接，因此将连接信息添加到他们的展台（使用p0或p1索引，这是未使用的(-1) ），并增加两个点的使用。

      3. 尝试连接点distance=sqrt(2)

         除了现在select8邻居的对angular线的距离之外，几乎与＃2相同。 这次也避免了闭环，所以不要连接已连接的点。

      4. 尝试连接最近的点

         与＃2，＃3几乎一样，但select最接近的点，也避免了闭环。

      5. 从pnt[]形成多边形

         所以select列表中的第一个点并将其添加到多边形。 然后将连接的点添加到它（不pipe你开始p0或p1 ）。 然后添加其连接点（不同于之前添加的点到多边形以避免前后循环）。 添加很多点，因为你有一个pnt[]点pnt[] 。

Answer 2

Delauney三angular测量可以帮助。 它具有三angular形中任何三angular形的外接圆都没有input点的属性。 因此，孔边界点将通过覆盖该孔的更大/更宽的三angular形连接。 在你的情况下，三angular测量将会有很多相似尺寸的三angular形，以及一些尺寸较大的覆盖三angular形的孔。 大概过滤大一点的东西就足够了，连接它们就可以find一个洞了。

Answer 3

这是我的爱好者非科学解决scheme：

1 – 以最小的预定义步长（dx，dy）扫描所有的2D区域。 对于每一步coordfind更大的圆，可以适应没有任何点内。 丢弃半径小于预定义大小的所有圆圈。

2 – 现在find所有的碰撞圈组，容易testing距离和半径，存储和分组列表中。 （问，如果你想要更多关于如何分组他们的细节，真的很容易）

3 – 为每组圆圈找出凹面的边界多边形，这与在已经写好的一组点上find凸多边形的algorithm非常相似，而且vector之间的最后一个问题angular度是相关的。

笔记

优化提示：在步骤1之前，可以将所有点存储在一个网格中，以便简化距离计算并将其限制在给定圆半径的近格方格内。

精度：对于较小的扫描步长值和允许的最小圆弧半径，可以获得更高的精度。

没有自己testing，但我相信它的工作原理。 祝你好运！

Answer 4

使用Delaunay三angular剖分findGabriel图可能会更好。 然后，您可以对Gabriel图进行angular度sorting，并进行圆形步行以生成凸多边形的列表。 然后可以按面积对这些多边形进行sorting。 你会对最大的面积感兴趣。

修改angular度sorting图也是更高效的，你可以沿着从A到B的path，顺时针或逆时针（从angular度sorting）看下一个是什么。 字典的冠词可能是有用的，它被定义为“graphics[A] [B] =（顺时针，逆时针）”。 一个使用字典词典方法的示例algorithm（python）。

 pre_graph = gabriel_graph(point) graph = {} for a in pre_graph: graph[a] = {} angle_sorted = sort(pre_graph[a], key=calc_angle_from(a)) for i,b in enumerate(angle_sorted): clockwise = angle_sorted[(i - 1) % len(angle_sorted)] counterclockwise = angle_sorted[(i + 1) % len(angle_sorted)] graph[a][b] = (clockwise, counterclockwise) polygons = [] for A in points: for B in graph[A]: for direction in [0,1]: polygon = [A] next_point = B: while next != A: polygon.append(next) next_point = graph[A][B][direction] if polygon[0] == min(polygon): # This should avoid duplicates polygons.add(polygon) 

结合Ante的build议也是有用的。

Answer 5

我不知道有什么algorithm可以解决我的头脑问题，但是有一件事情可以尝试（这是我的第一个冲动），类似于平滑粒子stream体力学等无网格方法中的密度计算方法。

如果可以计算空间中任意点的密度值（而不仅仅是在给定的采样点），则可以将孔的边界定义为密度函数的水平集曲线。 也就是说，这些洞是密度下降到某个阈值以下（你可能允许用户configuration）的地方。 你可以使用像行军广场的东西来find边界。

如果你想要说明这些密度插值函数是如何工作的，我可以提供（如果你喜欢的话）（尽我所能和知识）。

我不知道这会有多好，但希望它会给你一些方向。

Answer 6

这是一个想法：

• 对于每个点xfind距离d(x,y) （其中y是与x最接近的邻居）。 如上定义f(x)=d(x,y) 。
• 找出f(x)的均值和方差。
• find“离群值” – 它们的f值离平均值非常远的点，至less远远大于\ alpha标准偏差。 (\alpha is a parameter for the algorithm).
• This will detect the 'holes' – all you have to do now is set the outlier of each hole.

Answer 7

It seems that you can address this by means of (binary) mathematical morphology on images.

Create a white image and plot all your points. Then "inflate" them into rectangles which are larger than the normal horizontal and vertical spacing. You do it by means of a so called erosion operation with a rectangular structuring element.

Doing this you will fill the plane, except at places where the points are too sparse.

The unfilled areas that you detect this way are smaller than the actual voids. You will restore to the full size by applying a dilation, using the same structuring element.

Both transforms combined are called an opening.

http://en.wikipedia.org/wiki/Mathematical_morphology