Finding holes in 2d point sets?

Solution 1:

what about some bitmap+vector approach like this:

1. obtain bounding box of point cloud area coverage

   Do this if it is not already known. It should be simple O(N) cycle through all points.

2. create map[N][N] of the area

   It is a 'bitmap' of the area for easy data density computation. Just create projection from area(x,y) -> map[i][j] for example with simple scale. Grid size N is also the accuracy of the output and must be bigger then average point distance !!! so each cell inside map[][] covers area with at least one point (if not in hole area).

3. compute data density for each cell of map[][]

   Easy as pie just clear map[][].cnt (counter of points) to zero and compute by simple O(N) cycle where do map[i][j].cnt++ for all points(x,y)

4. create list of unused area (map[][].cnt==0) or (map[][].cnt<=treshold)

   I do it by Horizontal and Vertical lines for simplicity

5. segmentate output

   Just group lines of the same hole together (intersecting ones ... vector approach) and also can be done in bullet #4 by flood fill (bitmap approach)

6. polygonize output

   Take all edge points of H,V lines of the same hole/group and create polygon (sort them so their connection does not intersect anything). There are lot of libs,algorithms and source code about this.

My source code for this approach:

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

Just ignore my glview2D stuff (it is my gfx render engine for geometry)

• view.pnt[] is dynamic list of your points (generated by random)
• view.lin[] is dynamic list output H,V lines for visualization only
• lin[] is your lines output

This is output:

I am too lazy to add polygonize for now you can see that segmentation works (coloring). If you need also help with polygonize then comment me but I think that should not be any problem.

Complexity estimation depends on the overall hole coverage

but for most of the code it is O(N) and for hole search/segmentation ~O((M^2)+(U^2)) where:

• N is point count
• M is map grid size
• U is H,V lines count dependent on holes ...
• M << N, U << M*M

as you can see for 3783 points 30x30 grid on the image above it took almost 9ms on my setup

[Edit1] played with vector polygonize a little

for simple holes is fine but for more complicated ones there are some hick-ups yet

[Edit2] finally got a little time for this so here is it:

This is simple class for hole/polygon search in more pleasant/manageable form:

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

You just need to replace my List<T> template with std::list or whatever (that template I cannot share). It is an dynamic 1D array of T ...

• List<int> x; is the same as int x[];
• x.add(); add empty item to x
• x.add(a); add a item to x
• x.reset() clears the array
• x.allocate(size) preallocate space to avoid reallocations on the run which is slow
• x.num is number of items in x[] ... used size in items

in the original code are only static arrays so if you are confused check with it instead.

Now how to use it:

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();

where view.pnt[] is list of input points and inside it: view.pnt[i].p0.p[ 2 ]= { x,y }

Output is in h.lin[] and lin_i0 where:

• h.lin[i] i= < 0,lin_i0 ) are the inside H,V lines
• h.lin[i] i= < lin_i0,h.lin.num ) are the perimeter

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[]

         so pick first point in list and add it to the polygon. then add the connected point to it (does not matter which way you start p0 or p1). Then add its connected point (different then previous added point to polygon to avoid back and forward loops). Add so many points as you have points in a pnt[].

Solution 2:

Delauney triangulation can help. It has property that no input point is inside the circumcircle of any triangle in triangulation. Because of that, hole boundary points will be connected by larger/wider triangles covering that hole. In your cases triangulation will have lot of triangles of similar size, and some triangles of larger size that cover holes. Probably it is enough to filter larger ones and connect them to find a hole.

Solution 3:

This is my enthusiast non-scientific solution:

1 - Scan all the 2D area with minimum predefined step (dx, dy). For each step coord find the bigger circle that could fit without any point inside. Discard all the circles with radius less than a predefined size.

2 - Now find all groups of colliding circles, easy test of distance and radius, store and group in separated lists. (Ask, if you want more details about how to grouping them, is really easy )

3 - Find the concave bounding polygon for each group of circles, very similar to the algorithm to find the convex polygon around a group of points you already wrote, and your last question angles between vectors was related.

Notes

Optimization tips: Before step 1, you can store all points in a grid|matrix so distance calculation are simplified and limited to near grid squares of the given circle radius.

Precision: You gain more precision for smaller values of scan step and minimum circle radius allowed.

Not tested by myself but, I'm sure it works. Good luck!

Solution 4:

You'd probably be better off using the Delaunay triangulation to find the Gabriel graph. You then angle sort the Gabriel graph and do circular walks to generate a list of convex polygons. You can then sort those polygons by area. You'd be interested in the ones with the largest area.

It'll also be more efficient to modify the angle sorted graph that that you can follow the path from A to B, and see what is next either clockwise or counter clockwise (from the angle sort). A dictonary of dictionaries might be useful, that is defined something like "graph[A][B] = (clockwise, counterclockwise)". An example algorithm (python) using the dictionary of dictionaries approach.

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)

It may also be useful to combine with Ante's suggestion as well.

Solution 5:

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