在2d点集找到漏洞？

Question

我有一套2d points.它们X,Y coordinates位于标准的笛卡尔网格系统上(在本例中为a UTM zone).我需要在该点集中找到孔,最好能够设置找到这些孔的算法的灵敏度.通常这些点集非常密集,但有些可能密度低得多.

它们是什么,如果它有帮助的话,那就是田地里的土壤被采样的点,农业人们显然发现它们有用.有时在这些点样品中有巨大的岩石或沼泽地或满是小湖泊和池塘.

从点集中,他们希望我找到粗略定义这些孔的凹多边形.

我已经编写了找到外部凹面边界多边形的算法,并允许它们设置由算法形成的多边形粗糙或平滑的灵敏度.在那之后,他们想找到洞并将那些洞作为凹多边形给予它们,我猜在某些情况下可能是凸的,但这将是边缘情况而不是常态.

问题是我听过的关于这个问题的唯一论文是那些希望在太空中找到大空虚的天文学家完成的论文,我从来没有能够找到他们的一篇论文,其中的实际算法显示为除了粗略的概念描述之外的任何东西

我已经尝试了谷歌和各种学术论文搜索等,但到目前为止我还没有找到很多有用的东西.这让我想知道是否有这个问题的名字,我不知道所以我在寻找错误的东西或什么？

所以经过那个冗长的解释之后,我的问题是:有没有人知道我应该寻找什么,最好用定义明确的算法找到这方面的论文,或者有人知道一个他们可以指出我的算法吗？

任何帮助我解决这个问题的东西都会非常有用并且非常感激,即使只是正确的搜索短语,也会发现潜在的算法或论文.

这将实现的语言将是C#,但是从Mathematica软件包到其他任何东西的例子MATLAB or ASM, C, C++, Python, Java or MathCAD都可以,只要示例中没有一些调用就像是FindTheHole等等.FindTheHole没有定义或者是对于实现软件而言是专有的,例如,MathCAD示例通常具有很多.

下面是两个实际点集的例子,一个是密集的,一个是稀疏的,我需要找到它们中的区域:

Answer 1

怎么样的位图+矢量方法如下:

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

   如果还不知道,请执行此操作.它应该是O(N)所有点的简单循环.

2. 创造map[N][N]该地区

   它是该区域的"位图",便于进行数据密度计算.只需area(x,y) -> map[i][j]通过简单的比例创建投影即可 .网格大小N也是输出的精度,必须大于平均点距离!所以每个单元内部map[][]覆盖至少一个点(如果不在孔区域内).

3. 计算每个细胞的数据密度 map[][]

   像馅饼一样容易清楚map[][].cnt(点数的计数器)zero并通过简单的O(N)循环来计算map[i][j].cnt++ 所有的事情points(x,y)

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

   为简单起见,我通过水平和垂直线进行操作

5. 细分输出

   只需将同一个孔的线组合在一起(相交的......矢量方法),也可以通过泛洪填充(位图方法)在子弹#4中完成

6. 多边形输出

   获取相同孔/组的H,V线的所有边缘点并创建多边形(对它们进行排序,使它们的连接不与任何东西相交).有很多关于这个的库,算法和源代码.

我的这种方法的源代码:

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东西(它是几何的我的gfx渲染引擎)

• view.pnt[] 是你的积分的动态列表(随机生成)
• view.lin[]是动态列表输出H,V行仅用于可视化
• lin[] 是你的线输出

这是输出:

我太懒了添加polygonize现在你可以看到分割工作(着色).如果你还需要帮助polygonize然后评论我,但我认为这不应该是任何问题.

复杂性估计取决于整个孔的覆盖范围

但对于大多数代码来说,它是O(N)用于洞搜索/分段的~O((M^2)+(U^2))地方:

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

正如你可以看到上面图像上的3783点30x30网格几乎9ms在我的设置上

[Edit1]使用矢量多边形玩了一点

对于简单的孔是好的,但对于更复杂的孔,还有一些hick-up

[Edit2]终于有了一点时间,所以这就是它:

这是一个简单的孔/多边形搜索类,以更愉快/可管理的形式:

//---------------------------------------------------------------------------
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=a.n; 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++)
            {
            p.i=a->i0;
            p.j=a->j0;
            map[p.i][p.j]=0;
            for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++)
             if ((aa->i==p.i)&&(aa->j==p.j)) { e=-1; break; }
            if (e>=0) pnt.add(p);
            p.i=a->i1;
            p.j=a->j1;
            map[p.i][p.j]=0;
            for (aa=&pnt[0],e=0;e<pnt.num;e++,aa++)
             if ((aa->i==p.i)&&(aa->j==p.j)) { 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>模板替换我的模板std::list(我无法共享的模板).这是一个动态的1D阵列T......

• List<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[]是输入点和里面列表: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 ) 是外围

The perimeter lines are not ordered and are duplicated twice so just order them and remove duplicates (too lazy for that). Inside lin[] are id .. id of hole 0,1,2,3,... to which the line belongs and i,j coordinates inside map. so for proper output into your world coordinates do something like this:

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
    }

• my view.lin[] has members: p0,p1, which are points as view.pnt[] and col which is color

I saw only one issue with this when holes are too small (diameter < 3 cells) otherwise is OK

[edit4] reordering perimeter lines

to do that just instead of this:

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

do this:

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

[Edit5] How polygonize inside holes::holes_end works

As input for this you need the list of all H,V lines lin[] segmentated/grouped/sorted by hole and the density map map[][].

1. loop through all holes

   2. loop through all H,V lines of processed hole

      Create list of all unique line endpoints pnt[] (no duplicates). So take 2 endpoints for each line and look if each point is already in the list. If not add it there else ignore it.

   3. delete all non border points from list

      So remove all points that have no contact with non-hole area by looking into 4 neighbors in the density map[][]

   4. do connected components analysis on the points

      1. set used=0; p0=-1; p1=-1; for all points in pnt[] list

      2. connect points with distance=1

         loop through all points pnt[] with used<2 which means they are not fully used yet and for each such point search pnt[] again for another such point that has also distance = 1 to it. It means it is its 4-neighbors and should be connected so add the connection info to booth of them (use p0 or p1 index which ever is unused (-1)) and increment usage of both points.

      3. try to connect points with distance=sqrt(2)

         is almost the same as #2 except the distance which now selects diagonals of 8-neighbors. This time also avoid closed loops so do not connect point that is already connected to it.

      4. try to connect closest points

         again is almost the same as #2,#3 but select the closest point instead and also avoid closed loops.

      5. form polygon from pnt[]

         所以选择列表中的第一个点并将其添加到多边形.然后添加连接指向它(不要紧,你开始哪种方式p0或p1).然后添加其连接点(不同于先前添加的点到多边形以避免后退和前进循环).添加如此多的积分,因为你有积分pnt[].