FindN8Lines.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00012 /***************************************************************************
00013  *   Copyright (C) 2009 Thomas K Sharpless                                 *
00014  *   tksharpless@gmail.com                                                 *
00015  *                                                                         *
00016  *   This program is free software; you can redistribute it and/or modify  *
00017  *   it under the terms of the GNU General Public License as published by  *
00018  *   the Free Software Foundation; either version 2 of the License, or     *
00019  *   (at your option) any later version.                                   *
00020  *                                                                         *
00021  *   This program is distributed in the hope that it will be useful,       *
00022  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00023  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00024  *   GNU General Public License for more details.                          *
00025  *                                                                         *
00026  *   You should have received a copy of the GNU General Public License     *
00027  *   along with this program.  If not, see <http://www.gnu.org/licenses/>. *
00028  ***************************************************************************/
00029 
00030 #include <assert.h>
00031 
00032 #include "FindN8Lines.h"
00033 #include <vigra/copyimage.hxx>
00034 #include <vigra/pixelneighborhood.hxx>
00035 
00036 using namespace vigra;
00037 using namespace std;
00038 
00039 namespace HuginLines
00040 {
00041 // colors used in lines image
00042 #define N8_bg  255
00043 #define N8_end 1
00044 #define N8_mid 96
00045 
00046 BImage edgeMap2linePts(BImage & input)
00047 {
00048     BImage edge(input.width(),input.height(),N8_bg);
00049     BImage line(input.size());
00050     //copy input image, left border pixel in background color
00051     vigra::copyImage(input.upperLeft()+Diff2D(1,1),input.lowerRight()-Diff2D(1,1),input.accessor(),edge.upperLeft(),edge.accessor());
00052 
00053     int width  = input.width(),
00054         height = input.height();
00055 
00056     BImage::traverser eul, lul;
00057 
00058     /* first pass erases "elbow" points  from edge image,
00059        leaving only diagonal connections.  An elbow point
00060        has one horizontal and one vertical neighbor, and no
00061        more than 4 neighbors in all.
00062 
00063     */
00064     eul = edge.upperLeft() + Diff2D( 1, 1 );
00065     for(int y=1; y<height-1; ++y, ++eul.y )
00066     {
00067         BImage::traverser ix = eul;
00068         for(int x=1; x<width-1; ++x, ++ix.x )
00069         {
00070           // center must be an edge point
00071             if( *ix != N8_bg )
00072             {
00073                 int n = 0;
00074                 int nh = 0, nv = 0, nu = 0, nd = 0;
00075                 if( ix( 1, 0 ) != N8_bg ) ++nh, ++n;
00076                 if( ix( 1, -1 ) != N8_bg ) ++nu, ++n;
00077                 if( ix( 0, -1 ) != N8_bg ) ++nv, ++n;
00078                 if( ix( -1, -1 ) != N8_bg ) ++nd, ++n;
00079                 if( ix( -1, 0 ) != N8_bg ) ++nh, ++n;
00080                 if( ix( -1, 1 ) != N8_bg ) ++nu, ++n;
00081                 if( ix( 0, 1 ) != N8_bg ) ++nv, ++n;
00082                 if( ix( 1, 1 ) != N8_bg ) ++nd, ++n;
00083 
00084                 if( nh == 1 && nv == 1 && n > 1 && n < 5 )
00085                 {
00086                     *ix = N8_bg;    // clear this point
00087                 }
00088             }
00089         }
00090     }
00091 
00092     /* second pass copies center points of 3x3 regions having line-like
00093        configurations, to the temp image, marked with the orientation of
00094        the inferred line.  As a result of pass 1, these points will have
00095        not more than 2 marked neighbors, in an oblique configuration.
00096        Orientation codes 1 to 8 correspond to 22.5 degree increments
00097        from horizontal.  We multiply these by 20 to make the differences
00098        visible on a debug image.
00099 
00100     */
00101     eul = (edge.upperLeft() + Diff2D(1,1));
00102     lul = (line.upperLeft() + Diff2D(1,1));
00103     line.init(N8_bg);
00104 
00105     for(int y=1; y<height-1; ++y, ++eul.y, ++lul.y)
00106     {
00107         BImage::traverser ix = eul;
00108         BImage::traverser ox = lul;
00109         for(int x=1; x<width-1; ++x, ++ix.x, ++ox.x)
00110         {
00111             if( *ix != N8_bg )
00112             {
00113                 int n = 0;
00114                 int nh = 0, nv = 0, nu = 0, nd = 0;
00115 
00116                 if( ix( 1, 0 ) != N8_bg ) ++nh, ++n;
00117                 if( ix( 1, -1 ) != N8_bg ) ++nu, ++n;
00118                 if( ix( 0, -1 ) != N8_bg ) ++nv, ++n;
00119                 if( ix( -1, -1 ) != N8_bg ) ++nd, ++n;
00120                 if( ix( -1, 0 ) != N8_bg ) ++nh, ++n;
00121                 if( ix( -1, 1 ) != N8_bg ) ++nu, ++n;
00122                 if( ix( 0, 1 ) != N8_bg ) ++nv, ++n;
00123                 if( ix( 1, 1 ) != N8_bg ) ++nd, ++n;
00124 
00125               // mark the output pixel
00126                 int code = 0;
00127                 if( n == 1 )
00128                 {
00129                     if( nh ) code = 1;
00130                     else if( nu ) code = 3;
00131                     else if( nv ) code = 5;
00132                     else code = 7;
00133                 }
00134                 else
00135                     if( n == 2 )
00136                     {
00137                         if( nh == 2 ) code = 1;
00138                         else if( nu == 2 ) code = 3;
00139                         else if( nv == 2 ) code = 5;
00140                         else if( nd == 2 ) code = 7;
00141                         else if( nh && nu ) code = 2;
00142                         else if( nh && nd ) code = 8;
00143                         else if( nv && nu ) code = 4;
00144                         else if( nv && nd ) code = 6;
00145                     }
00146                 assert( code < 9 );
00147                 *ox = code ? code : N8_bg;
00148             }
00149         }
00150     }
00151 
00152     /* 3rd pass erases points in the temp image having 2 neighbors whose
00153        orientations differ by more than one unit.  This breaks corners.
00154 
00155     */
00156     lul = (line.upperLeft() + Diff2D(1,1));
00157 
00158     for(int y=1; y<height-1; ++y, ++lul.y)
00159     {
00160         BImage::traverser ix = lul;
00161         for(int x=1; x<width-1; ++x, ++ix.x )
00162         {
00163             int c = *ix;
00164             if( c != N8_bg )
00165             {
00166                 assert( c > 0 && c < 9 );
00167                 int n = 0, omin = 9, omax = 0;
00168                 NeighborhoodCirculator<BImage::traverser, EightNeighborCode>
00169                                circulator(ix),
00170                                end(circulator);
00171                 do {
00172                     int nc = *circulator;
00173                     if( nc != N8_bg )
00174                     {
00175                         assert( nc > 0 && nc < 9 );
00176                         if( nc < omin ) omin = nc;
00177                         if( nc > omax ) omax = nc;
00178                         ++n;
00179                     }
00180                 } while(++circulator != end);
00181                 // mark the pixel
00182                 if( n == 2 )
00183                 {
00184                     int d =  omax - omin;
00185                     assert( d < 8 );
00186                     if( d > 4 )
00187                         d = 7 - d;
00188                     if ( d > 1 )
00189                         *ix = N8_bg;
00190                 }
00191             }
00192         }
00193     }
00194 
00195     /* 4th pass changes marks in temp image from orientation to
00196        end = 1, interior = 2, and erases isolated points.
00197     */
00198     lul = (line.upperLeft() + Diff2D(1,1));
00199 
00200     for(int y=1; y<height-1; ++y, ++lul.y)
00201     {
00202         BImage::traverser ix = lul;
00203         for(int x=1; x<width-1; ++x, ++ix.x )
00204         {
00205             int c = *ix;
00206             if( c != N8_bg )
00207             {
00208                 int n = 0;
00209                 NeighborhoodCirculator<BImage::traverser, EightNeighborCode> circulator(ix);
00210                 NeighborhoodCirculator<BImage::traverser, EightNeighborCode> end(circulator);
00211                 do {
00212                     int nc = *circulator;
00213                     if( nc != N8_bg )
00214                     {
00215                         ++n;
00216                     }
00217                 } while(++circulator != end);
00218                 // mark the pixel
00219                 if( n == 0 ) *ix = N8_bg;
00220                 else if( n == 1 ) *ix = 1;
00221                 else *ix = 2;
00222             }
00223         }
00224     }
00225 
00226     /* 5th pass copies line points to edge image, skipping end points that are
00227        neighbors of end points (this removes 2-point lines).  End points are
00228        marked N8_end, interior points N8_mid.
00229     */
00230     lul = (line.upperLeft() + Diff2D(1,1));
00231     eul = (edge.upperLeft() + Diff2D(1,1)); // rewind edge
00232     edge.init(N8_bg);
00233 
00234     for(int y=1; y<height-1; ++y, ++lul.y, ++eul.y )
00235     {
00236         BImage::traverser ox = eul;
00237         BImage::traverser ix = lul;
00238         for(int x=1; x<width-1; ++x, ++ix.x, ++ox.x )
00239         {
00240             int c = *ix;
00241             if( c != N8_bg )
00242             {
00243                 int n = 0;
00244                 NeighborhoodCirculator<BImage::traverser, EightNeighborCode> circulator(ix);
00245                 NeighborhoodCirculator<BImage::traverser, EightNeighborCode> end(circulator);
00246                 do {
00247                     int nc = *circulator;
00248                     if( nc != N8_bg )
00249                     {
00250                         ++n;
00251                         if( c == 1 && nc == 1 ) c = 0;    //skip
00252                     }
00253                 } while(++circulator != end);
00254 
00255               // mark the pixel
00256                 if( c )
00257                 {
00258                     if( n == 1 )
00259                     {
00260                         *ox = N8_end;
00261                     }
00262                     else
00263                     {
00264                         *ox = N8_mid;
00265                     }
00266                 }
00267             }
00268         }
00269     }
00270     return edge;
00271 }
00272 
00273 // chain-code distance
00274 inline float ccdist( int dx, int dy )
00275 {
00276     dx = abs(dx); dy = abs(dy);
00277     return float(max(dx,dy)) + float(0.41421 * min(dx,dy));
00278 }
00279 
00280 // Euclidian distance
00281 inline float eudist( int dx, int dy )
00282 {
00283     register float x = dx, y = dy;
00284     return sqrt( x*x + y*y );
00285 }
00286 
00287 // chord and arc using chain code distance
00288 static double CtoAcc( vector<Point2D> & pts,
00289                     int start, int count,
00290                     double & C, double & A)
00291 {
00292     int n = (int)pts.size() - start;
00293     if( count > n ) count = n;
00294     if( count < 3 ) return 1;
00295     int xl = pts.at(start).x, yl = pts.at(start).y;
00296     int xr = pts.at(start + count - 1).x,
00297         yr = pts.at(start + count - 1).y;
00298     C = ccdist( xr - xl, yr - yl );
00299     A = 0;
00300     for(int i = start + 1; i < count; i++ )
00301     {
00302         xr = pts.at(i).x; yr = pts.at(i).y;
00303         A += ccdist( xr - xl, yr - yl );
00304         xl = xr; yl = yr;
00305     }
00306     return C/A;
00307 }
00308 
00309 // chord and arc using Euclidian distance
00310 static double CtoAeu( vector<Point2D> & pts,
00311                     int start, int count,
00312                     double & C, double & A)
00313 {
00314     int n = (int)pts.size() - start;
00315     if( count > n ) count = n;
00316     if( count < 3 ) return 1;
00317     int xl = pts.at(start).x, yl = pts.at(start).y;
00318     int xr = pts.at(start + count - 1).x,
00319         yr = pts.at(start + count - 1).y;
00320     C = eudist( xr - xl, yr - yl );
00321     A = 0;
00322     for(int i = start + 1; i < count; i++ )
00323     {
00324         xr = pts.at(i).x; yr = pts.at(i).y;
00325         A += eudist( xr - xl, yr - yl );
00326         xl = xr; yl = yr;
00327     }
00328     return C/A;
00329 }
00330 
00331 /* 3-point scalar curvature
00332    positive clockwise
00333 */
00334 inline float scurv( Point2D & l, Point2D & m, Point2D & r )
00335 {
00336     return float(m.x - l.x) * float(r.y - l.y)
00337          + float(m.y - l.y) * float(l.x - r.x );
00338 }
00339 
00340 /* 3-point float vector curvature
00341    positive inward
00342 */
00343 inline void vcurv ( Point2D & l, Point2D & m, Point2D & r, float & vx, float & vy )
00344 {
00345     vx = 0.5 * (l.x + r.x) - float( m.x );
00346     vy = 0.5 * (l.y + r.y) - float( m.y );
00347 }
00348 
00349 
00350 /* final filter for line segments
00351   input: a point list, a minimum size, the f.l. in pixels and
00352     the projection center point.
00353   output: the point list may contain a segment, of at least
00354     mimimum size, having uniform curvature <= a limit that
00355     depends on f.l., position and orientation.
00356   return value is length of the pruned segment, or
00357       0 if rejectd on length or extreme initial curvature.
00358      -1 if rejected on orientation
00359      -2 if rejected on curvature
00360   Note if a segment has multiple smooth parts, only the longest
00361   one will survive.
00362 */
00363 static int lineFilter( vector< Point2D > & pts,
00364                        int lmin,
00365                        double flpix,
00366                        double xcen, double ycen
00367                       )
00368 {
00369     // cos( min angle from radius to chord )
00370     const double crcmin = 0.20;    // about 11 degrees
00371     // max ratio of actual to theoretical curvature
00372     const double cvrmax = 3.0;
00373 
00374     double r90 = 1.57 * flpix;    // 90 degree radius, pixels
00375 
00376     int n = (int)pts.size();
00377     // reject short segments
00378     if( n < lmin )
00379         return 0;    // too short
00380 
00381     // the full chord
00382     double x0 = pts.at(0).x,
00383            y0 = pts.at(0).y;
00384     double dx = pts.at(n-1).x - x0,
00385            dy = pts.at(n-1).y - y0;
00386     double chord = sqrt( dx*dx + dy*dy );
00387     if( chord < 0.6 * n )
00388         return 0;    // way too bent!
00389     // the chord midpoint (in centered coords)
00390     double ccx = x0 + 0.5 * dx - xcen,
00391            ccy = y0 + 0.5 * dy - ycen;
00392     // the arc center point in centered coords
00393     // aka radius vector
00394     double acx = pts.at(n/2).x - xcen,
00395            acy = pts.at(n/2).y - ycen;
00396     double acd = sqrt( acx*acx + acy*acy );
00397     // the curvature vector (points toward convex side)
00398     double cvx = acx - ccx,
00399            cvy = acy - ccy;
00400     double cvd = sqrt( cvx*cvx + cvy*cvy );
00401     // cos( angle from radius to chord )
00402     double cosrc =  fabs( acx * dy - acy * dx ) / (chord * acd);
00403     if( fabs(cosrc) < crcmin )
00404         return -1;    // too radial!
00405 
00406     /* curvature limit test */
00407     // curvature limit from lens model
00408     double S =  1.2 * chord * cosrc / r90;    // pseudo sine
00409     double C = r90 * (1 - sqrt( 1 - S * S ));    // corresp chord->arc dist
00410     // curvature test
00411     if( cvd / C > cvrmax )
00412         return -2;
00413 
00414     // for testing...
00415     return n;
00416 }
00417 
00418 /*  There are 2 line selection parameters:
00419   -- minimum size (in points )
00420   -- the assumed focal length (in pixels) of the lens
00421   The maximum allowed overall curvature depends on the f.l. and
00422   the position and orientation of the segment.  This calculation
00423   assumes the center of projection is in the center of the image.
00424 
00425   Passing a longer f.l. makes the curvature limits smaller. For
00426   rectilinear lenses, I suggest using a value 3 to 4 times the
00427   actual lens f.l.
00428 
00429   Lines whose orientation is radial, or nearly so, are removed
00430   as well, since they can contribute nothing to lens calibration.
00431 
00432   Note flpix == 0 disables the curvature & orientation tests.
00433 
00434   Besides the global curvature test, there are two local tests
00435   to ensure that the curvature is nearly uniform along the line.
00436   The first is a "corner test" that may split a line into smaller
00437   segments while it is being traced.  The second, part of the
00438   curvature filter, just prunes the ends, if possible, giving
00439   a single segment of nearly uniform curvature.
00440 */
00441 
00442 int linePts2lineList( BImage & img, int minsize, double flpix, Lines& lines)
00443 {
00444     int nadd = 0, nrejL = 0;
00445 
00446     static Diff2D offs[8] = {
00447         Diff2D( 1, 0 ),
00448         Diff2D( 1, -1 ),
00449         Diff2D( 0, -1 ),
00450         Diff2D( -1, -1 ),
00451         Diff2D( -1, 0 ),
00452         Diff2D( -1, 1 ),
00453         Diff2D( 0, 1 ),
00454         Diff2D( 1, 1 )
00455     };
00456 
00457     // corner filter parameters
00458     const int span = 10;
00459     const float maxacd = 1.4f;
00460 
00461     if(minsize < span) minsize = span; // else bend filter fails
00462 
00463     int width  = img.width(),
00464         height = img.height();
00465     double xcen = 0.5 * width,
00466            ycen = 0.5 * height;
00467 
00468 
00469     BImage::traverser ul(img.upperLeft() + Diff2D( 1, 1 ));
00470     for(int y=1; y<height-1; ++y, ++ul.y)
00471     {
00472         BImage::traverser ix = ul;
00473         for(int x=1; x<width-1; ++x, ++ix.x )
00474         {
00475             int cd = *ix;    // point code
00476             if( cd == N8_end )
00477             { // new line...
00478               // trace and erase the line
00479                 BImage::traverser tr = ix;
00480                 Point2D pos( x, y );
00481                 vector<Point2D> pts;
00482                 pts.push_back(pos);
00483                 *tr = N8_bg;
00484                 int dir;
00485                 Diff2D dif;
00486                 do {    // scan the neighborhood
00487                     for( dir = 0; dir < 8; dir++ )
00488                     {
00489                         dif = offs[dir];
00490                         if( tr[dif] != N8_bg ) break;
00491                     }
00492                     assert( dir < 8 );
00493                     tr += dif;        // move to next point
00494                     cd = *tr;        // pick up that pixel
00495                     *tr = N8_bg;    // erase it from image
00496                     pos += dif;        // update position
00497                     pts.push_back(pos);    // add same to pointlist
00498                 } while( cd != N8_end );
00499                 // validate the point list
00500                 if( pts.size() < minsize )
00501                 {
00502                     ++nrejL;
00503                 }
00504                 else
00505                     if( flpix > 0 )
00506                     {
00507                         int ncopy = 0, nskip = 0;
00508                         /* bend test: compute running chord and arc on a fixed span.
00509                             Output segments where their difference is small.
00510                             The parameters are set to pass some hooked lines
00511                             that will be cleaned up later.
00512                         */
00513                         int ip = 0, np = (int)pts.size();
00514 
00515                         vector<Point2D> tmp;
00516                         float ccd[32];    // rolling interpoint chaincode distances
00517                         int isql = 0, isqr = 0;    // left & rgt indices to same
00518                         int xl = pts.at( 0 ).x,
00519                             yl = pts.at( 0 ).y;
00520                         float dx = pts.at(span-1).x - xl,
00521                               dy = pts.at(span-1).y - yl;
00522                         float Chord = sqrt(dx*dx + dy*dy);
00523                         float Arc = 0;
00524                         // fill 1st span's d's, initialize Dsq
00525                         for( isqr = 0; isqr < span-1; isqr++ )
00526                         {
00527                             register int x = pts.at( isqr + 1).x,
00528                                          y = pts.at( isqr + 1 ).y;
00529                             float d = ccdist( x - xl, y - yl );
00530                             ccd[isqr] = d;
00531                             Arc += d;
00532                             xl = x; yl = y;
00533                         }
00534                         // copy 1st span pnts if span is good
00535                         if( Arc - Chord <= maxacd )
00536                         {
00537                             for( int i = 0; i < span; i++ )
00538                             {
00539                                 tmp.push_back(Point2D( pts.at(i).x, pts.at(i).y ));
00540                             }
00541                         }
00542 
00543                         // roll span....
00544                         for( ip = 1; ip < np - span; ip++ )
00545                         {
00546                             int irgt = ip + span - 1;
00547                             // update span
00548                             Arc -= ccd[isql];
00549                             isql = (++isql)&31;
00550                             isqr = (++isqr)&31;
00551                             int x = pts.at( irgt ).x,
00552                                 y = pts.at( irgt ).y;
00553                             float d = ccdist( x - xl, y - yl );
00554                             ccd[isqr] = d;
00555                             Arc += d;
00556                             xl = x; yl = y;
00557                             dx = x - pts.at(ip).x; dy  = y - pts.at(ip).y;
00558                             Chord = sqrt( dx*dx + dy*dy );
00559                             /* if this span is good, copy right pnt
00560                             else flush tmp to lines if long enough
00561                             */
00562                             if( Arc - Chord <= maxacd )
00563                             {
00564                                 ++ncopy;
00565                                 tmp.push_back(Point2D( pts.at(irgt).x, pts.at(irgt).y ));
00566                             }
00567                             else
00568                             {
00569                                 int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
00570                                 SingleLine line;
00571                                 line.line=tmp;
00572                                 if( q > 0 )
00573                                 {
00574                                   ++nadd;
00575                                   line.status=valid_line;
00576                                 }
00577                                 else if( q == 0 ) line.status=bad_length;
00578                                 else if( q == -1 ) line.status=bad_orientation;
00579                                 else if( q == -2 ) line.status=bad_curvature;
00580                                 lines.push_back(line);
00581                                 tmp.clear();
00582                                 ++nskip;
00583                             }
00584                         }
00585                         // if the final span is good, copy & flush.
00586                         int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
00587                         SingleLine line;
00588                         line.line=tmp;
00589                         if( q > 0 )
00590                         {
00591                             ++nadd;
00592                             line.status=valid_line;
00593                         }
00594                         else if( q == 0 ) line.status=bad_length;
00595                         else if( q == -1 ) line.status=bad_orientation;
00596                         else if( q == -2 ) line.status=bad_curvature;
00597                         lines.push_back(line);
00598                     } // end curvature test
00599                     else
00600                     {    // length ok, no curvature test
00601                         SingleLine line;
00602                         line.line=pts;
00603                         line.status=valid_line;
00604                         lines.push_back(line);
00605                         ++nadd;
00606                     }
00607                 }  // end trace line
00608         }    // end x scan
00609     }  // end y scan
00610 
00611     return nadd;
00612 }
00613 
00614 }; //namespace

Generated on 4 Sep 2015 for Hugintrunk by  doxygen 1.4.7