Hugintrunk  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FindN8Lines.cpp
Go to the documentation of this file.
1 // -*- c-basic-offset: 4 -*-
12 /***************************************************************************
13  * Copyright (C) 2009 Thomas K Sharpless *
14  * tksharpless@gmail.com *
15  * *
16  * This program is free software; you can redistribute it and/or modify *
17  * it under the terms of the GNU General Public License as published by *
18  * the Free Software Foundation; either version 2 of the License, or *
19  * (at your option) any later version. *
20  * *
21  * This program is distributed in the hope that it will be useful, *
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
24  * GNU General Public License for more details. *
25  * *
26  * You should have received a copy of the GNU General Public License *
27  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
28  ***************************************************************************/
29 
30 #include <assert.h>
31 
32 #include "FindN8Lines.h"
33 #include <vigra/copyimage.hxx>
34 #include <vigra/pixelneighborhood.hxx>
35 
36 namespace HuginLines
37 {
38 // colors used in lines image
39 #define N8_bg 255
40 #define N8_end 1
41 #define N8_mid 96
42 
43 vigra::BImage edgeMap2linePts(vigra::BImage & input)
44 {
45  vigra::BImage edge(input.width(), input.height(), N8_bg);
46  vigra::BImage line(input.size());
47  //copy input image, left border pixel in background color
48  vigra::copyImage(input.upperLeft() + vigra::Diff2D(1, 1), input.lowerRight() - vigra::Diff2D(1, 1), input.accessor(), edge.upperLeft(), edge.accessor());
49 
50  int width = input.width(),
51  height = input.height();
52 
53  vigra::BImage::traverser eul, lul;
54 
55  /* first pass erases "elbow" points from edge image,
56  leaving only diagonal connections. An elbow point
57  has one horizontal and one vertical neighbor, and no
58  more than 4 neighbors in all.
59 
60  */
61  eul = edge.upperLeft() + vigra::Diff2D(1, 1);
62  for(int y=1; y<height-1; ++y, ++eul.y )
63  {
64  vigra::BImage::traverser ix = eul;
65  for(int x=1; x<width-1; ++x, ++ix.x )
66  {
67  // center must be an edge point
68  if( *ix != N8_bg )
69  {
70  int n = 0;
71  int nh = 0, nv = 0, nu = 0, nd = 0;
72  if( ix( 1, 0 ) != N8_bg ) ++nh, ++n;
73  if( ix( 1, -1 ) != N8_bg ) ++nu, ++n;
74  if( ix( 0, -1 ) != N8_bg ) ++nv, ++n;
75  if( ix( -1, -1 ) != N8_bg ) ++nd, ++n;
76  if( ix( -1, 0 ) != N8_bg ) ++nh, ++n;
77  if( ix( -1, 1 ) != N8_bg ) ++nu, ++n;
78  if( ix( 0, 1 ) != N8_bg ) ++nv, ++n;
79  if( ix( 1, 1 ) != N8_bg ) ++nd, ++n;
80 
81  if( nh == 1 && nv == 1 && n > 1 && n < 5 )
82  {
83  *ix = N8_bg; // clear this point
84  }
85  }
86  }
87  }
88 
89  /* second pass copies center points of 3x3 regions having line-like
90  configurations, to the temp image, marked with the orientation of
91  the inferred line. As a result of pass 1, these points will have
92  not more than 2 marked neighbors, in an oblique configuration.
93  Orientation codes 1 to 8 correspond to 22.5 degree increments
94  from horizontal. We multiply these by 20 to make the differences
95  visible on a debug image.
96 
97  */
98  eul = (edge.upperLeft() + vigra::Diff2D(1, 1));
99  lul = (line.upperLeft() + vigra::Diff2D(1, 1));
100  line.init(N8_bg);
101 
102  for(int y=1; y<height-1; ++y, ++eul.y, ++lul.y)
103  {
104  vigra::BImage::traverser ix = eul;
105  vigra::BImage::traverser ox = lul;
106  for(int x=1; x<width-1; ++x, ++ix.x, ++ox.x)
107  {
108  if( *ix != N8_bg )
109  {
110  int n = 0;
111  int nh = 0, nv = 0, nu = 0, nd = 0;
112 
113  if( ix( 1, 0 ) != N8_bg ) ++nh, ++n;
114  if( ix( 1, -1 ) != N8_bg ) ++nu, ++n;
115  if( ix( 0, -1 ) != N8_bg ) ++nv, ++n;
116  if( ix( -1, -1 ) != N8_bg ) ++nd, ++n;
117  if( ix( -1, 0 ) != N8_bg ) ++nh, ++n;
118  if( ix( -1, 1 ) != N8_bg ) ++nu, ++n;
119  if( ix( 0, 1 ) != N8_bg ) ++nv, ++n;
120  if( ix( 1, 1 ) != N8_bg ) ++nd, ++n;
121 
122  // mark the output pixel
123  int code = 0;
124  if( n == 1 )
125  {
126  if( nh ) code = 1;
127  else if( nu ) code = 3;
128  else if( nv ) code = 5;
129  else code = 7;
130  }
131  else
132  if( n == 2 )
133  {
134  if( nh == 2 ) code = 1;
135  else if( nu == 2 ) code = 3;
136  else if( nv == 2 ) code = 5;
137  else if( nd == 2 ) code = 7;
138  else if( nh && nu ) code = 2;
139  else if( nh && nd ) code = 8;
140  else if( nv && nu ) code = 4;
141  else if( nv && nd ) code = 6;
142  }
143  assert( code < 9 );
144  *ox = code ? code : N8_bg;
145  }
146  }
147  }
148 
149  /* 3rd pass erases points in the temp image having 2 neighbors whose
150  orientations differ by more than one unit. This breaks corners.
151 
152  */
153  lul = (line.upperLeft() + vigra::Diff2D(1, 1));
154 
155  for(int y=1; y<height-1; ++y, ++lul.y)
156  {
157  vigra::BImage::traverser ix = lul;
158  for(int x=1; x<width-1; ++x, ++ix.x )
159  {
160  int c = *ix;
161  if( c != N8_bg )
162  {
163  assert( c > 0 && c < 9 );
164  int n = 0, omin = 9, omax = 0;
165  vigra::NeighborhoodCirculator<vigra::BImage::traverser, vigra::EightNeighborCode>
166  circulator(ix),
167  end(circulator);
168  do {
169  int nc = *circulator;
170  if( nc != N8_bg )
171  {
172  assert( nc > 0 && nc < 9 );
173  if( nc < omin ) omin = nc;
174  if( nc > omax ) omax = nc;
175  ++n;
176  }
177  } while(++circulator != end);
178  // mark the pixel
179  if( n == 2 )
180  {
181  int d = omax - omin;
182  assert( d < 8 );
183  if( d > 4 )
184  d = 7 - d;
185  if ( d > 1 )
186  *ix = N8_bg;
187  }
188  }
189  }
190  }
191 
192  /* 4th pass changes marks in temp image from orientation to
193  end = 1, interior = 2, and erases isolated points.
194  */
195  lul = (line.upperLeft() + vigra::Diff2D(1, 1));
196 
197  for(int y=1; y<height-1; ++y, ++lul.y)
198  {
199  vigra::BImage::traverser ix = lul;
200  for(int x=1; x<width-1; ++x, ++ix.x )
201  {
202  int c = *ix;
203  if( c != N8_bg )
204  {
205  int n = 0;
206  vigra::NeighborhoodCirculator<vigra::BImage::traverser, vigra::EightNeighborCode> circulator(ix);
207  vigra::NeighborhoodCirculator<vigra::BImage::traverser, vigra::EightNeighborCode> end(circulator);
208  do {
209  int nc = *circulator;
210  if( nc != N8_bg )
211  {
212  ++n;
213  }
214  } while(++circulator != end);
215  // mark the pixel
216  if( n == 0 ) *ix = N8_bg;
217  else if( n == 1 ) *ix = 1;
218  else *ix = 2;
219  }
220  }
221  }
222 
223  /* 5th pass copies line points to edge image, skipping end points that are
224  neighbors of end points (this removes 2-point lines). End points are
225  marked N8_end, interior points N8_mid.
226  */
227  lul = (line.upperLeft() + vigra::Diff2D(1, 1));
228  eul = (edge.upperLeft() + vigra::Diff2D(1, 1)); // rewind edge
229  edge.init(N8_bg);
230 
231  for(int y=1; y<height-1; ++y, ++lul.y, ++eul.y )
232  {
233  vigra::BImage::traverser ox = eul;
234  vigra::BImage::traverser ix = lul;
235  for(int x=1; x<width-1; ++x, ++ix.x, ++ox.x )
236  {
237  int c = *ix;
238  if( c != N8_bg )
239  {
240  int n = 0;
241  vigra::NeighborhoodCirculator<vigra::BImage::traverser, vigra::EightNeighborCode> circulator(ix);
242  vigra::NeighborhoodCirculator<vigra::BImage::traverser, vigra::EightNeighborCode> end(circulator);
243  do {
244  int nc = *circulator;
245  if( nc != N8_bg )
246  {
247  ++n;
248  if( c == 1 && nc == 1 ) c = 0; //skip
249  }
250  } while(++circulator != end);
251 
252  // mark the pixel
253  if( c )
254  {
255  if( n == 1 )
256  {
257  *ox = N8_end;
258  }
259  else
260  {
261  *ox = N8_mid;
262  }
263  }
264  }
265  }
266  }
267  return edge;
268 }
269 
270 // chain-code distance
271 inline float ccdist( int dx, int dy )
272 {
273  dx = abs(dx); dy = abs(dy);
274  return float(std::max(dx,dy)) + float(0.41421 * std::min(dx,dy));
275 }
276 
277 // Euclidian distance
278 inline float eudist( int dx, int dy )
279 {
280  register float x = dx, y = dy;
281  return sqrt( x*x + y*y );
282 }
283 
284 /* 3-point scalar curvature
285  positive clockwise
286 */
287 inline float scurv(vigra::Point2D & l, vigra::Point2D & m, vigra::Point2D & r)
288 {
289  return float(m.x - l.x) * float(r.y - l.y)
290  + float(m.y - l.y) * float(l.x - r.x );
291 }
292 
293 /* 3-point float vector curvature
294  positive inward
295 */
296 inline void vcurv(vigra::Point2D & l, vigra::Point2D & m, vigra::Point2D & r, float & vx, float & vy)
297 {
298  vx = 0.5 * (l.x + r.x) - float( m.x );
299  vy = 0.5 * (l.y + r.y) - float( m.y );
300 }
301 
302 
303 /* final filter for line segments
304  input: a point list, a minimum size, the f.l. in pixels and
305  the projection center point.
306  output: the point list may contain a segment, of at least
307  mimimum size, having uniform curvature <= a limit that
308  depends on f.l., position and orientation.
309  return value is length of the pruned segment, or
310  0 if rejectd on length or extreme initial curvature.
311  -1 if rejected on orientation
312  -2 if rejected on curvature
313  Note if a segment has multiple smooth parts, only the longest
314  one will survive.
315 */
316 static int lineFilter(std::vector< vigra::Point2D > & pts,
317  int lmin,
318  double flpix,
319  double xcen, double ycen
320  )
321 {
322  // cos( min angle from radius to chord )
323  const double crcmin = 0.20; // about 11 degrees
324  // max ratio of actual to theoretical curvature
325  const double cvrmax = 3.0;
326 
327  double r90 = 1.57 * flpix; // 90 degree radius, pixels
328 
329  int n = (int)pts.size();
330  // reject short segments
331  if( n < lmin )
332  return 0; // too short
333 
334  // the full chord
335  double x0 = pts.at(0).x,
336  y0 = pts.at(0).y;
337  double dx = pts.at(n-1).x - x0,
338  dy = pts.at(n-1).y - y0;
339  double chord = sqrt( dx*dx + dy*dy );
340  if( chord < 0.6 * n )
341  return 0; // way too bent!
342  // the chord midpoint (in centered coords)
343  double ccx = x0 + 0.5 * dx - xcen,
344  ccy = y0 + 0.5 * dy - ycen;
345  // the arc center point in centered coords
346  // aka radius vector
347  double acx = pts.at(n/2).x - xcen,
348  acy = pts.at(n/2).y - ycen;
349  double acd = sqrt( acx*acx + acy*acy );
350  // the curvature vector (points toward convex side)
351  double cvx = acx - ccx,
352  cvy = acy - ccy;
353  double cvd = sqrt( cvx*cvx + cvy*cvy );
354  // cos( angle from radius to chord )
355  double cosrc = fabs( acx * dy - acy * dx ) / (chord * acd);
356  if( fabs(cosrc) < crcmin )
357  return -1; // too radial!
358 
359  /* curvature limit test */
360  // curvature limit from lens model
361  double S = 1.2 * chord * cosrc / r90; // pseudo sine
362  double C = r90 * (1 - sqrt( 1 - S * S )); // corresp chord->arc dist
363  // curvature test
364  if( cvd / C > cvrmax )
365  return -2;
366 
367  // for testing...
368  return n;
369 }
370 
371 /* There are 2 line selection parameters:
372  -- minimum size (in points )
373  -- the assumed focal length (in pixels) of the lens
374  The maximum allowed overall curvature depends on the f.l. and
375  the position and orientation of the segment. This calculation
376  assumes the center of projection is in the center of the image.
377 
378  Passing a longer f.l. makes the curvature limits smaller. For
379  rectilinear lenses, I suggest using a value 3 to 4 times the
380  actual lens f.l.
381 
382  Lines whose orientation is radial, or nearly so, are removed
383  as well, since they can contribute nothing to lens calibration.
384 
385  Note flpix == 0 disables the curvature & orientation tests.
386 
387  Besides the global curvature test, there are two local tests
388  to ensure that the curvature is nearly uniform along the line.
389  The first is a "corner test" that may split a line into smaller
390  segments while it is being traced. The second, part of the
391  curvature filter, just prunes the ends, if possible, giving
392  a single segment of nearly uniform curvature.
393 */
394 
395 int linePts2lineList(vigra::BImage & img, int minsize, double flpix, Lines& lines)
396 {
397  int nadd = 0, nrejL = 0;
398 
399  static vigra::Diff2D offs[8] = {
400  vigra::Diff2D(1, 0),
401  vigra::Diff2D(1, -1),
402  vigra::Diff2D(0, -1),
403  vigra::Diff2D(-1, -1),
404  vigra::Diff2D(-1, 0),
405  vigra::Diff2D(-1, 1),
406  vigra::Diff2D(0, 1),
407  vigra::Diff2D(1, 1)
408  };
409 
410  // corner filter parameters
411  const int span = 10;
412  const float maxacd = 2.0f;
413 
414  if(minsize < span) minsize = span; // else bend filter fails
415 
416  int width = img.width(),
417  height = img.height();
418  double xcen = 0.5 * width,
419  ycen = 0.5 * height;
420 
421 
422  vigra::BImage::traverser ul(img.upperLeft() + vigra::Diff2D(1, 1));
423  for(int y=1; y<height-1; ++y, ++ul.y)
424  {
425  vigra::BImage::traverser ix = ul;
426  for(int x=1; x<width-1; ++x, ++ix.x )
427  {
428  int cd = *ix; // point code
429  if( cd == N8_end )
430  { // new line...
431  // trace and erase the line
432  vigra::BImage::traverser tr = ix;
433  vigra::Point2D pos(x, y);
434  std::vector<vigra::Point2D> pts;
435  pts.push_back(pos);
436  *tr = N8_bg;
437  int dir;
438  vigra::Diff2D dif;
439  do { // scan the neighborhood
440  for( dir = 0; dir < 8; dir++ )
441  {
442  dif = offs[dir];
443  if( tr[dif] != N8_bg ) break;
444  }
445  assert( dir < 8 );
446  tr += dif; // move to next point
447  cd = *tr; // pick up that pixel
448  *tr = N8_bg; // erase it from image
449  pos += dif; // update position
450  pts.push_back(pos); // add same to pointlist
451  } while( cd != N8_end );
452  // validate the point list
453  if( pts.size() < minsize )
454  {
455  ++nrejL;
456  }
457  else
458  if( flpix > 0 )
459  {
460  int ncopy = 0, nskip = 0;
461  /* bend test: compute running chord and arc on a fixed span.
462  Output segments where their difference is small.
463  The parameters are set to pass some hooked lines
464  that will be cleaned up later.
465  */
466  int ip = 0, np = (int)pts.size();
467 
468  std::vector<vigra::Point2D> tmp;
469  float ccd[32]{}; // rolling interpoint chaincode distances
470  int isql = 0, isqr = 0; // left & rgt indices to same
471  int xl = pts.at( 0 ).x,
472  yl = pts.at( 0 ).y;
473  float dx = pts.at(span-1).x - xl,
474  dy = pts.at(span-1).y - yl;
475  float Chord = sqrt(dx*dx + dy*dy);
476  float Arc = 0;
477  // fill 1st span's d's, initialize Dsq
478  for( isqr = 0; isqr < span-1; isqr++ )
479  {
480  register int x = pts.at( isqr + 1).x,
481  y = pts.at( isqr + 1 ).y;
482  float d = ccdist( x - xl, y - yl );
483  ccd[isqr] = d;
484  Arc += d;
485  xl = x; yl = y;
486  }
487  // copy 1st span pnts if span is good
488  if( Arc - Chord <= maxacd )
489  {
490  for( int i = 0; i < span; i++ )
491  {
492  tmp.push_back(vigra::Point2D(pts.at(i).x, pts.at(i).y));
493  }
494  }
495 
496  // roll span....
497  for( ip = 1; ip < np - span; ip++ )
498  {
499  int irgt = ip + span - 1;
500  // update span
501  Arc -= ccd[isql];
502  isql = (isql + 1) & 31;
503  isqr = (isqr + 1) & 31;
504  int x = pts.at( irgt ).x,
505  y = pts.at( irgt ).y;
506  float d = ccdist( x - xl, y - yl );
507  ccd[isqr] = d;
508  Arc += d;
509  xl = x; yl = y;
510  dx = x - pts.at(ip).x; dy = y - pts.at(ip).y;
511  Chord = sqrt( dx*dx + dy*dy );
512  /* if this span is good, copy right pnt
513  else flush tmp to lines if long enough
514  */
515  if( Arc - Chord <= maxacd )
516  {
517  ++ncopy;
518  tmp.push_back(vigra::Point2D(pts.at(irgt).x, pts.at(irgt).y));
519  }
520  else
521  {
522  int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
524  line.line=tmp;
525  if( q > 0 )
526  {
527  ++nadd;
528  line.status=valid_line;
529  }
530  else if( q == 0 ) line.status=bad_length;
531  else if( q == -1 ) line.status=bad_orientation;
532  else if( q == -2 ) line.status=bad_curvature;
533  lines.push_back(line);
534  tmp.clear();
535  ++nskip;
536  }
537  }
538  // if the final span is good, copy & flush.
539  int q = lineFilter( tmp, minsize, flpix, xcen, ycen );
541  line.line=tmp;
542  if( q > 0 )
543  {
544  ++nadd;
545  line.status=valid_line;
546  }
547  else if( q == 0 ) line.status=bad_length;
548  else if( q == -1 ) line.status=bad_orientation;
549  else if( q == -2 ) line.status=bad_curvature;
550  lines.push_back(line);
551  } // end curvature test
552  else
553  { // length ok, no curvature test
555  line.line=pts;
556  line.status=valid_line;
557  lines.push_back(line);
558  ++nadd;
559  }
560  } // end trace line
561  } // end x scan
562  } // end y scan
563 
564  return nadd;
565 }
566 
567 }; //namespace
std::vector< vigra::Point2D > line
Definition: LinesTypes.h:45
a single line extracted from image
Definition: LinesTypes.h:43
static char * line
Definition: svm.cpp:2784
#define N8_end
Definition: FindN8Lines.cpp:40
#define N8_bg
Definition: FindN8Lines.cpp:39
declaration of find lines algorithm
float ccdist(int dx, int dy)
#define N8_mid
Definition: FindN8Lines.cpp:41
static int lineFilter(std::vector< vigra::Point2D > &pts, int lmin, double flpix, double xcen, double ycen)
float scurv(vigra::Point2D &l, vigra::Point2D &m, vigra::Point2D &r)
std::vector< SingleLine > Lines
vector of extracted lines from image
Definition: LinesTypes.h:50
static T max(T x, T y)
Definition: svm.cpp:65
void copyImage(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc)
Definition: openmp_vigra.h:305
void vcurv(vigra::Point2D &l, vigra::Point2D &m, vigra::Point2D &r, float &vx, float &vy)
float eudist(int dx, int dy)
static T min(T x, T y)
Definition: svm.cpp:62
int linePts2lineList(vigra::BImage &img, int minsize, double flpix, Lines &lines)
converts a linePts image to a list of lines
vigra::BImage edgeMap2linePts(vigra::BImage &input)
marks line point
Definition: FindN8Lines.cpp:43