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 // chord and arc using chain code distance
00285 static double CtoAcc( std::vector<vigra::Point2D> & pts,
00286                     int start, int count,
00287                     double & C, double & A)
00288 {
00289     int n = (int)pts.size() - start;
00290     if( count > n ) count = n;
00291     if( count < 3 ) return 1;
00292     int xl = pts.at(start).x, yl = pts.at(start).y;
00293     int xr = pts.at(start + count - 1).x,
00294         yr = pts.at(start + count - 1).y;
00295     C = ccdist( xr - xl, yr - yl );
00296     A = 0;
00297     for(int i = start + 1; i < count; i++ )
00298     {
00299         xr = pts.at(i).x; yr = pts.at(i).y;
00300         A += ccdist( xr - xl, yr - yl );
00301         xl = xr; yl = yr;
00302     }
00303     return C/A;
00304 }
00305 
00306 // chord and arc using Euclidian distance
00307 static double CtoAeu( std::vector<vigra::Point2D> & pts,
00308                     int start, int count,
00309                     double & C, double & A)
00310 {
00311     int n = (int)pts.size() - start;
00312     if( count > n ) count = n;
00313     if( count < 3 ) return 1;
00314     int xl = pts.at(start).x, yl = pts.at(start).y;
00315     int xr = pts.at(start + count - 1).x,
00316         yr = pts.at(start + count - 1).y;
00317     C = eudist( xr - xl, yr - yl );
00318     A = 0;
00319     for(int i = start + 1; i < count; i++ )
00320     {
00321         xr = pts.at(i).x; yr = pts.at(i).y;
00322         A += eudist( xr - xl, yr - yl );
00323         xl = xr; yl = yr;
00324     }
00325     return C/A;
00326 }
00327 
00328 /* 3-point scalar curvature
00329    positive clockwise
00330 */
00331 inline float scurv(vigra::Point2D & l, vigra::Point2D & m, vigra::Point2D & r)
00332 {
00333     return float(m.x - l.x) * float(r.y - l.y)
00334          + float(m.y - l.y) * float(l.x - r.x );
00335 }
00336 
00337 /* 3-point float vector curvature
00338    positive inward
00339 */
00340 inline void vcurv(vigra::Point2D & l, vigra::Point2D & m, vigra::Point2D & r, float & vx, float & vy)
00341 {
00342     vx = 0.5 * (l.x + r.x) - float( m.x );
00343     vy = 0.5 * (l.y + r.y) - float( m.y );
00344 }
00345 
00346 
00347 /* final filter for line segments
00348   input: a point list, a minimum size, the f.l. in pixels and
00349     the projection center point.
00350   output: the point list may contain a segment, of at least
00351     mimimum size, having uniform curvature <= a limit that
00352     depends on f.l., position and orientation.
00353   return value is length of the pruned segment, or
00354       0 if rejectd on length or extreme initial curvature.
00355      -1 if rejected on orientation
00356      -2 if rejected on curvature
00357   Note if a segment has multiple smooth parts, only the longest
00358   one will survive.
00359 */
00360 static int lineFilter(std::vector< vigra::Point2D > & pts,
00361                        int lmin,
00362                        double flpix,
00363                        double xcen, double ycen
00364                       )
00365 {
00366     // cos( min angle from radius to chord )
00367     const double crcmin = 0.20;    // about 11 degrees
00368     // max ratio of actual to theoretical curvature
00369     const double cvrmax = 3.0;
00370 
00371     double r90 = 1.57 * flpix;    // 90 degree radius, pixels
00372 
00373     int n = (int)pts.size();
00374     // reject short segments
00375     if( n < lmin )
00376         return 0;    // too short
00377 
00378     // the full chord
00379     double x0 = pts.at(0).x,
00380            y0 = pts.at(0).y;
00381     double dx = pts.at(n-1).x - x0,
00382            dy = pts.at(n-1).y - y0;
00383     double chord = sqrt( dx*dx + dy*dy );
00384     if( chord < 0.6 * n )
00385         return 0;    // way too bent!
00386     // the chord midpoint (in centered coords)
00387     double ccx = x0 + 0.5 * dx - xcen,
00388            ccy = y0 + 0.5 * dy - ycen;
00389     // the arc center point in centered coords
00390     // aka radius vector
00391     double acx = pts.at(n/2).x - xcen,
00392            acy = pts.at(n/2).y - ycen;
00393     double acd = sqrt( acx*acx + acy*acy );
00394     // the curvature vector (points toward convex side)
00395     double cvx = acx - ccx,
00396            cvy = acy - ccy;
00397     double cvd = sqrt( cvx*cvx + cvy*cvy );
00398     // cos( angle from radius to chord )
00399     double cosrc =  fabs( acx * dy - acy * dx ) / (chord * acd);
00400     if( fabs(cosrc) < crcmin )
00401         return -1;    // too radial!
00402 
00403     /* curvature limit test */
00404     // curvature limit from lens model
00405     double S =  1.2 * chord * cosrc / r90;    // pseudo sine
00406     double C = r90 * (1 - sqrt( 1 - S * S ));    // corresp chord->arc dist
00407     // curvature test
00408     if( cvd / C > cvrmax )
00409         return -2;
00410 
00411     // for testing...
00412     return n;
00413 }
00414 
00415 /*  There are 2 line selection parameters:
00416   -- minimum size (in points )
00417   -- the assumed focal length (in pixels) of the lens
00418   The maximum allowed overall curvature depends on the f.l. and
00419   the position and orientation of the segment.  This calculation
00420   assumes the center of projection is in the center of the image.
00421 
00422   Passing a longer f.l. makes the curvature limits smaller. For
00423   rectilinear lenses, I suggest using a value 3 to 4 times the
00424   actual lens f.l.
00425 
00426   Lines whose orientation is radial, or nearly so, are removed
00427   as well, since they can contribute nothing to lens calibration.
00428 
00429   Note flpix == 0 disables the curvature & orientation tests.
00430 
00431   Besides the global curvature test, there are two local tests
00432   to ensure that the curvature is nearly uniform along the line.
00433   The first is a "corner test" that may split a line into smaller
00434   segments while it is being traced.  The second, part of the
00435   curvature filter, just prunes the ends, if possible, giving
00436   a single segment of nearly uniform curvature.
00437 */
00438 
00439 int linePts2lineList(vigra::BImage & img, int minsize, double flpix, Lines& lines)
00440 {
00441     int nadd = 0, nrejL = 0;
00442 
00443     static vigra::Diff2D offs[8] = {
00444         vigra::Diff2D(1, 0),
00445         vigra::Diff2D(1, -1),
00446         vigra::Diff2D(0, -1),
00447         vigra::Diff2D(-1, -1),
00448         vigra::Diff2D(-1, 0),
00449         vigra::Diff2D(-1, 1),
00450         vigra::Diff2D(0, 1),
00451         vigra::Diff2D(1, 1)
00452     };
00453 
00454     // corner filter parameters
00455     const int span = 10;
00456     const float maxacd = 1.4f;
00457 
00458     if(minsize < span) minsize = span; // else bend filter fails
00459 
00460     int width  = img.width(),
00461         height = img.height();
00462     double xcen = 0.5 * width,
00463            ycen = 0.5 * height;
00464 
00465 
00466     vigra::BImage::traverser ul(img.upperLeft() + vigra::Diff2D(1, 1));
00467     for(int y=1; y<height-1; ++y, ++ul.y)
00468     {
00469         vigra::BImage::traverser ix = ul;
00470         for(int x=1; x<width-1; ++x, ++ix.x )
00471         {
00472             int cd = *ix;    // point code
00473             if( cd == N8_end )
00474             { // new line...
00475               // trace and erase the line
00476                 vigra::BImage::traverser tr = ix;
00477                 vigra::Point2D pos(x, y);
00478                 std::vector<vigra::Point2D> pts;
00479                 pts.push_back(pos);
00480                 *tr = N8_bg;
00481                 int dir;
00482                 vigra::Diff2D dif;
00483                 do {    // scan the neighborhood
00484                     for( dir = 0; dir < 8; dir++ )
00485                     {
00486                         dif = offs[dir];
00487                         if( tr[dif] != N8_bg ) break;
00488                     }
00489                     assert( dir < 8 );
00490                     tr += dif;        // move to next point
00491                     cd = *tr;        // pick up that pixel
00492                     *tr = N8_bg;    // erase it from image
00493                     pos += dif;        // update position
00494                     pts.push_back(pos);    // add same to pointlist
00495                 } while( cd != N8_end );
00496                 // validate the point list
00497                 if( pts.size() < minsize )
00498                 {
00499                     ++nrejL;
00500                 }
00501                 else
00502                     if( flpix > 0 )
00503                     {
00504                         int ncopy = 0, nskip = 0;
00505                         /* bend test: compute running chord and arc on a fixed span.
00506                             Output segments where their difference is small.
00507                             The parameters are set to pass some hooked lines
00508                             that will be cleaned up later.
00509                         */
00510                         int ip = 0, np = (int)pts.size();
00511 
00512                         std::vector<vigra::Point2D> tmp;
00513                         float ccd[32];    // rolling interpoint chaincode distances
00514                         int isql = 0, isqr = 0;    // left & rgt indices to same
00515                         int xl = pts.at( 0 ).x,
00516                             yl = pts.at( 0 ).y;
00517                         float dx = pts.at(span-1).x - xl,
00518                               dy = pts.at(span-1).y - yl;
00519                         float Chord = sqrt(dx*dx + dy*dy);
00520                         float Arc = 0;
00521                         // fill 1st span's d's, initialize Dsq
00522                         for( isqr = 0; isqr < span-1; isqr++ )
00523                         {
00524                             register int x = pts.at( isqr + 1).x,
00525                                          y = pts.at( isqr + 1 ).y;
00526                             float d = ccdist( x - xl, y - yl );
00527                             ccd[isqr] = d;
00528                             Arc += d;
00529                             xl = x; yl = y;
00530                         }
00531                         // copy 1st span pnts if span is good
00532                         if( Arc - Chord <= maxacd )
00533                         {
00534                             for( int i = 0; i < span; i++ )
00535                             {
00536                                 tmp.push_back(vigra::Point2D(pts.at(i).x, pts.at(i).y));
00537                             }
00538                         }
00539 
00540                         // roll span....
00541                         for( ip = 1; ip < np - span; ip++ )
00542                         {
00543                             int irgt = ip + span - 1;
00544                             // update span
00545                             Arc -= ccd[isql];
00546                             isql = (++isql)&31;
00547                             isqr = (++isqr)&31;
00548                             int x = pts.at( irgt ).x,
00549                                 y = pts.at( irgt ).y;
00550                             float d = ccdist( x - xl, y - yl );
00551                             ccd[isqr] = d;
00552                             Arc += d;
00553                             xl = x; yl = y;
00554                             dx = x - pts.at(ip).x; dy  = y - pts.at(ip).y;
00555                             Chord = sqrt( dx*dx + dy*dy );
00556                             /* if this span is good, copy right pnt
00557                             else flush tmp to lines if long enough
00558                             */
00559                             if( Arc - Chord <= maxacd )
00560                             {
00561                                 ++ncopy;
00562                                 tmp.push_back(vigra::Point2D(pts.at(irgt).x, pts.at(irgt).y));
00563                             }
00564                             else
00565                             {
00566                                 int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
00567                                 SingleLine line;
00568                                 line.line=tmp;
00569                                 if( q > 0 )
00570                                 {
00571                                   ++nadd;
00572                                   line.status=valid_line;
00573                                 }
00574                                 else if( q == 0 ) line.status=bad_length;
00575                                 else if( q == -1 ) line.status=bad_orientation;
00576                                 else if( q == -2 ) line.status=bad_curvature;
00577                                 lines.push_back(line);
00578                                 tmp.clear();
00579                                 ++nskip;
00580                             }
00581                         }
00582                         // if the final span is good, copy & flush.
00583                         int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
00584                         SingleLine line;
00585                         line.line=tmp;
00586                         if( q > 0 )
00587                         {
00588                             ++nadd;
00589                             line.status=valid_line;
00590                         }
00591                         else if( q == 0 ) line.status=bad_length;
00592                         else if( q == -1 ) line.status=bad_orientation;
00593                         else if( q == -2 ) line.status=bad_curvature;
00594                         lines.push_back(line);
00595                     } // end curvature test
00596                     else
00597                     {    // length ok, no curvature test
00598                         SingleLine line;
00599                         line.line=pts;
00600                         line.status=valid_line;
00601                         lines.push_back(line);
00602                         ++nadd;
00603                     }
00604                 }  // end trace line
00605         }    // end x scan
00606     }  // end y scan
00607 
00608     return nadd;
00609 }
00610 
00611 }; //namespace

Generated on 13 Feb 2016 for Hugintrunk by  doxygen 1.4.7