Determine Minimum Number of Line Segments to Solve

Question

I have a problem where I need to define a polygon with the minimum number of vertices that intersects or contains every pixel in an image that is not transparent (Let N be the number of pixels in the image). My only assumptions are that the image cannot contain transparent pixels within its boundary (holes), and that at least two pixels are non-transparent. As an example, lets say I have the following image:

My idea for an algorithm is thus:

1) Determine the edge pixels.
This is done in O(N) time by iterating through each pixel, and determining whether any neighbors (out of the four pixels left, above, right, and below it) are empty. Store the pixel and which neighbors were non-transparent, keyed by linear index into the array of pixels. Let there be P edge pixels, shown in orange below.

2) Get an adjacency list of the edge pixels.
This is done in O(P) time by selecting one of the edge pixels, and chosing a direction based on empty neighbors. For example, if a pixel has a bottom and right neighbor, then the next pixel will be either one in the upper-right corner, or the pixel immediately to the right. Select the next edge pixel from the dictionary of remaining edge pixels. Append that pixel to the list until the algorithm returns to the starting pixel. There are 27 edge pixels in the example image below (some are an edge pixel more than once).

3) Draw a maze that all edges must lie between.
This is done in O(P) time by iterating through the adjacency list, and adding an edge to all edges on those pixels without a neighbor. In addition, an edge is added to the interior of the shape based on the direction to the next edge pixel. If the pixel represents a peninsula with single pixel width, the inner edge is added to the middle of the edge instead of the pixel vertex. The interior of the maze is shown in red. Note that the maze boundary is a super-set of all the edge pixels found in step 2.

4) Find a polygon with almost minimal edges that does not touch the border of the maze.
This is the part I need help with. Does anyone have a suggestion of how you would go from step #3 to a solution such as the following?

Answer 1

1. See find holes in 2D point set
   • it is very similar problem
   • just invert the map and use midpoint of each grid square as point
   • that will lead to this:
   • 
   • as you want the inner polygon then there are 2 choices
   • shrink shape by 1 cell before applying this (can loose some detail in shape)
   • change the H/V lines so they are 1 cell shorter (half from both sides)
   • that will lead to something like this:
   • 
   • after some changes in code I can use 2x multi-sampling now so result is:
   • 
   • now you can join connected lines with the same slope
   • and apply something like ear clipping on the found corners to get more close to your desired polygon

Here the inverted source code in C++ (there may be some hole comments left):

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
/*  usage:
    int i;
    pntcloud_polygons h;
    pnt2D point[points];

    h.scann_beg(); for (i=0;i<points;i++) { p=point[i]; h.scann_pnt(p.x,p.y); } h.scann_end();
    h.cell_size(2.5);   // or (h.x1-h.x0)*0.01  ... cell size >> avg point distance
    h.holes_beg(); for (i=0;i<points;i++) { p=point[i]; h.holes_pnt(p.x,p.y); } h.holes_end();

*/
//---------------------------------------------------------------------------
class pntcloud_polygons
    {
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
        float 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
    pntcloud_polygons()  { 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; };
    pntcloud_polygons(pntcloud_polygons& a){ *this=a; };
    ~pntcloud_polygons() { _free(); };
    pntcloud_polygons* operator = (const pntcloud_polygons *a) { *this=*a; return this; };
    pntcloud_polygons* operator = (const pntcloud_polygons &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 pntcloud_polygons 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); }
    inline void l2g(double &x,double &y,float i,float 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 pntcloud_polygons::scann_beg()
    {
    x0=0.0; y0=0.0; x1=0.0; y1=0.0; n=0;
    }
//---------------------------------------------------------------------------
void pntcloud_polygons::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 pntcloud_polygons::scann_end()
    {
    }
//---------------------------------------------------------------------------
void pntcloud_polygons::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 pntcloud_polygons::holes_beg()
    {
    int i,j;
    for (i=0;i<xs;i++)
     for (j=0;j<ys;j++)
      map[i][j]=0;
    }
//---------------------------------------------------------------------------
void pntcloud_polygons::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 pntcloud_polygons::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 pntcloud_polygons (map[i][j].cnt!=0) or (map[i][j].cnt>=treshold)
    // and create lin[] list of H,V lines covering pntcloud_polygons
    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 polygon
        for (;i<xs;i++) if (map[i][j]==0) break; i1=i;      // find end of polygon
        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+0.5;
        l.i1=i1-0.5;
        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 polygon
        for (;j<ys;j++) if (map[i][j]==0) break; j1=j  ;    // find end of polygon
        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+0.5;
        l.j1=j1-0.5;
        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 losed 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;
        // 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;
            }
        }
    }
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

• it is the same as the code in holes link above
• just the map[][] conditions inverted to search polygons instead of holes
• and found HV lines are shrinked by half of cell from each side
• _lin coordinates are float now so o need for 4x multisampling
• the best results are with 2x multi-sampling (to avoid polygon width=1)
• so each cell in your map add as 2x2 points in the cell area
• I added also 2 corner points to better align map[][] and your image

Answer 2

I have no background in image processing, but I came across the Ramer–Douglas–Peucker algorithm yesterday and I think it might be helpful.

From my quick scan of the Wikipedia article, it reduces the number of point in a curve, so I would run this algorithm on each line where the points of the line are the end points you have and also set as points the borders of squares that the line crosses.

I marked out in this image two lines you could run the algorithm on and I think it would work.

How to find each line and when to stop - not 100% sure, but I hope this was useful.