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

Generated on 5 May 2016 for Hugintrunk by  doxygen 1.4.7