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

Generated on Wed Jul 16 01:25:38 2014 for Hugintrunk by  doxygen 1.3.9.1