1 /*----------------------------------------------------------------------------
2
3 LSD - Line Segment Detector on digital images
4
5 This code is part of the following publication and was subject
6 to peer review:
7
8 "LSD: a Line Segment Detector" by Rafael Grompone von Gioi,
9 Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,
10 Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsd
11 http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd
12
13 Copyright (c) 2007-2011 rafael grompone von gioi <grompone@gmail.com>
14
15 This program is free software: you can redistribute it and/or modify
16 it under the terms of the GNU Affero General Public License as
17 published by the Free Software Foundation, either version 3 of the
18 License, or (at your option) any later version.
19
20 This program is distributed in the hope that it will be useful,
21 but WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 GNU Affero General Public License for more details.
24
25 You should have received a copy of the GNU Affero General Public License
26 along with this program. If not, see <http://www.gnu.org/licenses/>.
27
28 ----------------------------------------------------------------------------*/
29
30 /*----------------------------------------------------------------------------*/
31 /** @file lsd.c
32 LSD module code
33 @author rafael grompone von gioi <grompone@gmail.com>
34 */
35 /*----------------------------------------------------------------------------*/
36
37 /*----------------------------------------------------------------------------*/
38 /** @mainpage LSD code documentation
39
40 This is an implementation of the Line Segment Detector described
41 in the paper:
42
43 "LSD: A Fast Line Segment Detector with a False Detection Control"
44 by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel,
45 and Gregory Randall, IEEE Transactions on Pattern Analysis and
46 Machine Intelligence, vol. 32, no. 4, pp. 722-732, April, 2010.
47
48 and in more details in the CMLA Technical Report:
49
50 "LSD: A Line Segment Detector, Technical Report",
51 by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel,
52 Gregory Randall, CMLA, ENS Cachan, 2010.
53
54 The version implemented here includes some further improvements
55 described in the following publication, of which this code is part:
56
57 "LSD: a Line Segment Detector" by Rafael Grompone von Gioi,
58 Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,
59 Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsd
60 http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd
61
62 The module's main function is lsd().
63
64 The source code is contained in two files: lsd.h and lsd.c.
65
66 HISTORY:
67 - version 1.6 - nov 2011:
68 - changes in the interface,
69 - max_grad parameter removed,
70 - the factor 11 was added to the number of test
71 to consider the different precision values
72 tested,
73 - a minor bug corrected in the gradient sorting
74 code,
75 - the algorithm now also returns p and log_nfa
76 for each detection,
77 - a minor bug was corrected in the image scaling,
78 - the angle comparison in "isaligned" changed
79 from < to <=,
80 - "eps" variable renamed "log_eps",
81 - "lsd_scale_region" interface was added,
82 - minor changes to comments.
83 - version 1.5 - dec 2010: Changes in 'refine', -W option added,
84 and more comments added.
85 - version 1.4 - jul 2010: lsd_scale interface added and doxygen doc.
86 - version 1.3 - feb 2010: Multiple bug correction and improved code.
87 - version 1.2 - dec 2009: First full Ansi C Language version.
88 - version 1.1 - sep 2009: Systematic subsampling to scale 0.8 and
89 correction to partially handle "angle problem".
90 - version 1.0 - jan 2009: First complete Megawave2 and Ansi C Language
91 version.
92
93 @author rafael grompone von gioi <grompone@gmail.com>
94 */
95 /*----------------------------------------------------------------------------*/
96
97 #include <stdio.h>
98 #include <stdlib.h>
99 #include <math.h>
100 #include <limits.h>
101 #include <float.h>
102 #include "lsd.h"
103
104 /** ln(10) */
105 #ifndef M_LN10
106 #define M_LN10 2.30258509299404568402
107 #endif /* !M_LN10 */
108
109 /** PI */
110 #ifndef M_PI
111 #define M_PI 3.14159265358979323846
112 #endif /* !M_PI */
113
114 #ifndef FALSE
115 #define FALSE 0
116 #endif /* !FALSE */
117
118 #ifndef TRUE
119 #define TRUE 1
120 #endif /* !TRUE */
121
122 /** Label for pixels with undefined gradient. */
123 #define NOTDEF -1024.0
124
125 /** 3/2 pi */
126 #define M_3_2_PI 4.71238898038
127
128 /** 2 pi */
129 #define M_2__PI 6.28318530718
130
131 /** Label for pixels not used in yet. */
132 #define NOTUSED 0
133
134 /** Label for pixels already used in detection. */
135 #define USED 1
136
137 /*----------------------------------------------------------------------------*/
138 /** Chained list of coordinates.
139 */
140 struct coorlist
141 {
142 int x,y;
143 struct coorlist * next;
144 };
145
146 /*----------------------------------------------------------------------------*/
147 /** A point (or pixel).
148 */
149 struct point {int x,y;};
150
151
152 /*----------------------------------------------------------------------------*/
153 /*------------------------- Miscellaneous functions --------------------------*/
154 /*----------------------------------------------------------------------------*/
155
156 /*----------------------------------------------------------------------------*/
157 /** Fatal error, print a message to standard-error output and exit.
158 */
error(char * msg)159 static void error(char * msg)
160 {
161 fprintf(stderr,"LSD Error: %s\n",msg);
162 exit(EXIT_FAILURE);
163 }
164
165 /*----------------------------------------------------------------------------*/
166 /** Doubles relative error factor
167 */
168 #define RELATIVE_ERROR_FACTOR 100.0
169
170 /*----------------------------------------------------------------------------*/
171 /** Compare doubles by relative error.
172
173 The resulting rounding error after floating point computations
174 depend on the specific operations done. The same number computed by
175 different algorithms could present different rounding errors. For a
176 useful comparison, an estimation of the relative rounding error
177 should be considered and compared to a factor times EPS. The factor
178 should be related to the cumulated rounding error in the chain of
179 computation. Here, as a simplification, a fixed factor is used.
180 */
double_equal(double a,double b)181 static int double_equal(double a, double b)
182 {
183 double abs_diff,aa,bb,abs_max;
184
185 /* trivial case */
186 if( a == b ) return TRUE;
187
188 abs_diff = fabs(a-b);
189 aa = fabs(a);
190 bb = fabs(b);
191 abs_max = aa > bb ? aa : bb;
192
193 /* DBL_MIN is the smallest normalized number, thus, the smallest
194 number whose relative error is bounded by DBL_EPSILON. For
195 smaller numbers, the same quantization steps as for DBL_MIN
196 are used. Then, for smaller numbers, a meaningful "relative"
197 error should be computed by dividing the difference by DBL_MIN. */
198 if( abs_max < DBL_MIN ) abs_max = DBL_MIN;
199
200 /* equal if relative error <= factor x eps */
201 return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
202 }
203
204 /*----------------------------------------------------------------------------*/
205 /** Computes Euclidean distance between point (x1,y1) and point (x2,y2).
206 */
dist(double x1,double y1,double x2,double y2)207 static double dist(double x1, double y1, double x2, double y2)
208 {
209 return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
210 }
211
212
213 /*----------------------------------------------------------------------------*/
214 /*----------------------- 'list of n-tuple' data type ------------------------*/
215 /*----------------------------------------------------------------------------*/
216
217 /*----------------------------------------------------------------------------*/
218 /** 'list of n-tuple' data type
219
220 The i-th component of the j-th n-tuple of an n-tuple list 'ntl'
221 is accessed with:
222
223 ntl->values[ i + j * ntl->dim ]
224
225 The dimension of the n-tuple (n) is:
226
227 ntl->dim
228
229 The number of n-tuples in the list is:
230
231 ntl->size
232
233 The maximum number of n-tuples that can be stored in the
234 list with the allocated memory at a given time is given by:
235
236 ntl->max_size
237 */
238 typedef struct ntuple_list_s
239 {
240 unsigned int size;
241 unsigned int max_size;
242 unsigned int dim;
243 double * values;
244 } * ntuple_list;
245
246 /*----------------------------------------------------------------------------*/
247 /** Free memory used in n-tuple 'in'.
248 */
free_ntuple_list(ntuple_list in)249 static void free_ntuple_list(ntuple_list in)
250 {
251 if( in == NULL || in->values == NULL )
252 error("free_ntuple_list: invalid n-tuple input.");
253 free( (void *) in->values );
254 free( (void *) in );
255 }
256
257 /*----------------------------------------------------------------------------*/
258 /** Create an n-tuple list and allocate memory for one element.
259 @param dim the dimension (n) of the n-tuple.
260 */
new_ntuple_list(unsigned int dim)261 static ntuple_list new_ntuple_list(unsigned int dim)
262 {
263 ntuple_list n_tuple;
264
265 /* check parameters */
266 if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");
267
268 /* get memory for list structure */
269 n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );
270 if( n_tuple == NULL ) error("not enough memory.");
271
272 /* initialize list */
273 n_tuple->size = 0;
274 n_tuple->max_size = 1;
275 n_tuple->dim = dim;
276
277 /* get memory for tuples */
278 n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) );
279 if( n_tuple->values == NULL ) error("not enough memory.");
280
281 return n_tuple;
282 }
283
284 /*----------------------------------------------------------------------------*/
285 /** Enlarge the allocated memory of an n-tuple list.
286 */
enlarge_ntuple_list(ntuple_list n_tuple)287 static void enlarge_ntuple_list(ntuple_list n_tuple)
288 {
289 /* check parameters */
290 if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )
291 error("enlarge_ntuple_list: invalid n-tuple.");
292
293 /* duplicate number of tuples */
294 n_tuple->max_size *= 2;
295
296 /* realloc memory */
297 n_tuple->values = (double *) realloc( (void *) n_tuple->values,
298 n_tuple->dim * n_tuple->max_size * sizeof(double) );
299 if( n_tuple->values == NULL ) error("not enough memory.");
300 }
301
302 /*----------------------------------------------------------------------------*/
303 /** Add a 7-tuple to an n-tuple list.
304 */
add_7tuple(ntuple_list out,double v1,double v2,double v3,double v4,double v5,double v6,double v7)305 static void add_7tuple( ntuple_list out, double v1, double v2, double v3,
306 double v4, double v5, double v6, double v7 )
307 {
308 /* check parameters */
309 if( out == NULL ) error("add_7tuple: invalid n-tuple input.");
310 if( out->dim != 7 ) error("add_7tuple: the n-tuple must be a 7-tuple.");
311
312 /* if needed, alloc more tuples to 'out' */
313 if( out->size == out->max_size ) enlarge_ntuple_list(out);
314 if( out->values == NULL ) error("add_7tuple: invalid n-tuple input.");
315
316 /* add new 7-tuple */
317 out->values[ out->size * out->dim + 0 ] = v1;
318 out->values[ out->size * out->dim + 1 ] = v2;
319 out->values[ out->size * out->dim + 2 ] = v3;
320 out->values[ out->size * out->dim + 3 ] = v4;
321 out->values[ out->size * out->dim + 4 ] = v5;
322 out->values[ out->size * out->dim + 5 ] = v6;
323 out->values[ out->size * out->dim + 6 ] = v7;
324
325 /* update number of tuples counter */
326 out->size++;
327 }
328
329
330 /*----------------------------------------------------------------------------*/
331 /*----------------------------- Image Data Types -----------------------------*/
332 /*----------------------------------------------------------------------------*/
333
334 /*----------------------------------------------------------------------------*/
335 /** char image data type
336
337 The pixel value at (x,y) is accessed by:
338
339 image->data[ x + y * image->xsize ]
340
341 with x and y integer.
342 */
343 typedef struct image_char_s
344 {
345 unsigned char * data;
346 unsigned int xsize,ysize;
347 } * image_char;
348
349 /*----------------------------------------------------------------------------*/
350 /** Free memory used in image_char 'i'.
351 */
free_image_char(image_char i)352 static void free_image_char(image_char i)
353 {
354 if( i == NULL || i->data == NULL )
355 error("free_image_char: invalid input image.");
356 free( (void *) i->data );
357 free( (void *) i );
358 }
359
360 /*----------------------------------------------------------------------------*/
361 /** Create a new image_char of size 'xsize' times 'ysize'.
362 */
new_image_char(unsigned int xsize,unsigned int ysize)363 static image_char new_image_char(unsigned int xsize, unsigned int ysize)
364 {
365 image_char image;
366
367 /* check parameters */
368 if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");
369
370 /* get memory */
371 image = (image_char) malloc( sizeof(struct image_char_s) );
372 if( image == NULL ) error("not enough memory.");
373 image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),
374 sizeof(unsigned char) );
375 if( image->data == NULL ) error("not enough memory.");
376
377 /* set image size */
378 image->xsize = xsize;
379 image->ysize = ysize;
380
381 return image;
382 }
383
384 /*----------------------------------------------------------------------------*/
385 /** Create a new image_char of size 'xsize' times 'ysize',
386 initialized to the value 'fill_value'.
387 */
new_image_char_ini(unsigned int xsize,unsigned int ysize,unsigned char fill_value)388 static image_char new_image_char_ini( unsigned int xsize, unsigned int ysize,
389 unsigned char fill_value )
390 {
391 image_char image = new_image_char(xsize,ysize); /* create image */
392 unsigned int N = xsize*ysize;
393 unsigned int i;
394
395 /* check parameters */
396 if( image == NULL || image->data == NULL )
397 error("new_image_char_ini: invalid image.");
398
399 /* initialize */
400 for(i=0; i<N; i++) image->data[i] = fill_value;
401
402 return image;
403 }
404
405 /*----------------------------------------------------------------------------*/
406 /** int image data type
407
408 The pixel value at (x,y) is accessed by:
409
410 image->data[ x + y * image->xsize ]
411
412 with x and y integer.
413 */
414 typedef struct image_int_s
415 {
416 int * data;
417 unsigned int xsize,ysize;
418 } * image_int;
419
420 /*----------------------------------------------------------------------------*/
421 /** Create a new image_int of size 'xsize' times 'ysize'.
422 */
new_image_int(unsigned int xsize,unsigned int ysize)423 static image_int new_image_int(unsigned int xsize, unsigned int ysize)
424 {
425 image_int image;
426
427 /* check parameters */
428 if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");
429
430 /* get memory */
431 image = (image_int) malloc( sizeof(struct image_int_s) );
432 if( image == NULL ) error("not enough memory.");
433 image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );
434 if( image->data == NULL ) error("not enough memory.");
435
436 /* set image size */
437 image->xsize = xsize;
438 image->ysize = ysize;
439
440 return image;
441 }
442
443 /*----------------------------------------------------------------------------*/
444 /** Create a new image_int of size 'xsize' times 'ysize',
445 initialized to the value 'fill_value'.
446 */
new_image_int_ini(unsigned int xsize,unsigned int ysize,int fill_value)447 static image_int new_image_int_ini( unsigned int xsize, unsigned int ysize,
448 int fill_value )
449 {
450 image_int image = new_image_int(xsize,ysize); /* create image */
451 unsigned int N = xsize*ysize;
452 unsigned int i;
453
454 /* initialize */
455 for(i=0; i<N; i++) image->data[i] = fill_value;
456
457 return image;
458 }
459
460 /*----------------------------------------------------------------------------*/
461 /** double image data type
462
463 The pixel value at (x,y) is accessed by:
464
465 image->data[ x + y * image->xsize ]
466
467 with x and y integer.
468 */
469 typedef struct image_double_s
470 {
471 double * data;
472 unsigned int xsize,ysize;
473 } * image_double;
474
475 /*----------------------------------------------------------------------------*/
476 /** Free memory used in image_double 'i'.
477 */
free_image_double(image_double i)478 static void free_image_double(image_double i)
479 {
480 if( i == NULL || i->data == NULL )
481 error("free_image_double: invalid input image.");
482 free( (void *) i->data );
483 free( (void *) i );
484 }
485
486 /*----------------------------------------------------------------------------*/
487 /** Create a new image_double of size 'xsize' times 'ysize'.
488 */
new_image_double(unsigned int xsize,unsigned int ysize)489 static image_double new_image_double(unsigned int xsize, unsigned int ysize)
490 {
491 image_double image;
492
493 /* check parameters */
494 if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");
495
496 /* get memory */
497 image = (image_double) malloc( sizeof(struct image_double_s) );
498 if( image == NULL ) error("not enough memory.");
499 image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );
500 if( image->data == NULL ) error("not enough memory.");
501
502 /* set image size */
503 image->xsize = xsize;
504 image->ysize = ysize;
505
506 return image;
507 }
508
509 /*----------------------------------------------------------------------------*/
510 /** Create a new image_double of size 'xsize' times 'ysize'
511 with the data pointed by 'data'.
512 */
new_image_double_ptr(unsigned int xsize,unsigned int ysize,double * data)513 static image_double new_image_double_ptr( unsigned int xsize,
514 unsigned int ysize, double * data )
515 {
516 image_double image;
517
518 /* check parameters */
519 if( xsize == 0 || ysize == 0 )
520 error("new_image_double_ptr: invalid image size.");
521 if( data == NULL ) error("new_image_double_ptr: NULL data pointer.");
522
523 /* get memory */
524 image = (image_double) malloc( sizeof(struct image_double_s) );
525 if( image == NULL ) error("not enough memory.");
526
527 /* set image */
528 image->xsize = xsize;
529 image->ysize = ysize;
530 image->data = data;
531
532 return image;
533 }
534
535
536 /*----------------------------------------------------------------------------*/
537 /*----------------------------- Gaussian filter ------------------------------*/
538 /*----------------------------------------------------------------------------*/
539
540 /*----------------------------------------------------------------------------*/
541 /** Compute a Gaussian kernel of length 'kernel->dim',
542 standard deviation 'sigma', and centered at value 'mean'.
543
544 For example, if mean=0.5, the Gaussian will be centered
545 in the middle point between values 'kernel->values[0]'
546 and 'kernel->values[1]'.
547 */
gaussian_kernel(ntuple_list kernel,double sigma,double mean)548 static void gaussian_kernel(ntuple_list kernel, double sigma, double mean)
549 {
550 double sum = 0.0;
551 double val;
552 unsigned int i;
553
554 /* check parameters */
555 if( kernel == NULL || kernel->values == NULL )
556 error("gaussian_kernel: invalid n-tuple 'kernel'.");
557 if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");
558
559 /* compute Gaussian kernel */
560 if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);
561 kernel->size = 1;
562 for(i=0;i<kernel->dim;i++)
563 {
564 val = ( (double) i - mean ) / sigma;
565 kernel->values[i] = exp( -0.5 * val * val );
566 sum += kernel->values[i];
567 }
568
569 /* normalization */
570 if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
571 }
572
573 /*----------------------------------------------------------------------------*/
574 /** Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.
575
576 For example, scale=0.8 will give a result at 80% of the original size.
577
578 The image is convolved with a Gaussian kernel
579 @f[
580 G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}
581 @f]
582 before the sub-sampling to prevent aliasing.
583
584 The standard deviation sigma given by:
585 - sigma = sigma_scale / scale, if scale < 1.0
586 - sigma = sigma_scale, if scale >= 1.0
587
588 To be able to sub-sample at non-integer steps, some interpolation
589 is needed. In this implementation, the interpolation is done by
590 the Gaussian kernel, so both operations (filtering and sampling)
591 are done at the same time. The Gaussian kernel is computed
592 centered on the coordinates of the required sample. In this way,
593 when applied, it gives directly the result of convolving the image
594 with the kernel and interpolated to that particular position.
595
596 A fast algorithm is done using the separability of the Gaussian
597 kernel. Applying the 2D Gaussian kernel is equivalent to applying
598 first a horizontal 1D Gaussian kernel and then a vertical 1D
599 Gaussian kernel (or the other way round). The reason is that
600 @f[
601 G(x,y) = G(x) * G(y)
602 @f]
603 where
604 @f[
605 G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}}.
606 @f]
607 The algorithm first applies a combined Gaussian kernel and sampling
608 in the x axis, and then the combined Gaussian kernel and sampling
609 in the y axis.
610 */
gaussian_sampler(image_double in,double scale,double sigma_scale)611 static image_double gaussian_sampler( image_double in, double scale,
612 double sigma_scale )
613 {
614 image_double aux,out;
615 ntuple_list kernel;
616 unsigned int N,M,h,n,x,y,i;
617 int xc,yc,j,double_x_size,double_y_size;
618 double sigma,xx,yy,sum,prec;
619
620 /* check parameters */
621 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
622 error("gaussian_sampler: invalid image.");
623 if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");
624 if( sigma_scale <= 0.0 )
625 error("gaussian_sampler: 'sigma_scale' must be positive.");
626
627 /* compute new image size and get memory for images */
628 if( in->xsize * scale > (double) UINT_MAX ||
629 in->ysize * scale > (double) UINT_MAX )
630 error("gaussian_sampler: the output image size exceeds the handled size.");
631 N = (unsigned int) ceil( in->xsize * scale );
632 M = (unsigned int) ceil( in->ysize * scale );
633 aux = new_image_double(N,in->ysize);
634 out = new_image_double(N,M);
635
636 /* sigma, kernel size and memory for the kernel */
637 sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;
638 /*
639 The size of the kernel is selected to guarantee that the
640 the first discarded term is at least 10^prec times smaller
641 than the central value. For that, h should be larger than x, with
642 e^(-x^2/2sigma^2) = 1/10^prec.
643 Then,
644 x = sigma * sqrt( 2 * prec * ln(10) ).
645 */
646 prec = 3.0;
647 h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
648 n = 1+2*h; /* kernel size */
649 kernel = new_ntuple_list(n);
650
651 /* auxiliary double image size variables */
652 double_x_size = (int) (2 * in->xsize);
653 double_y_size = (int) (2 * in->ysize);
654
655 /* First subsampling: x axis */
656 for(x=0;x<aux->xsize;x++)
657 {
658 /*
659 x is the coordinate in the new image.
660 xx is the corresponding x-value in the original size image.
661 xc is the integer value, the pixel coordinate of xx.
662 */
663 xx = (double) x / scale;
664 /* coordinate (0.0,0.0) is in the center of pixel (0,0),
665 so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */
666 xc = (int) floor( xx + 0.5 );
667 gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );
668 /* the kernel must be computed for each x because the fine
669 offset xx-xc is different in each case */
670
671 for(y=0;y<aux->ysize;y++)
672 {
673 sum = 0.0;
674 for(i=0;i<kernel->dim;i++)
675 {
676 j = xc - h + i;
677
678 /* symmetry boundary condition */
679 while( j < 0 ) j += double_x_size;
680 while( j >= double_x_size ) j -= double_x_size;
681 if( j >= (int) in->xsize ) j = double_x_size-1-j;
682
683 sum += in->data[ j + y * in->xsize ] * kernel->values[i];
684 }
685 aux->data[ x + y * aux->xsize ] = sum;
686 }
687 }
688
689 /* Second subsampling: y axis */
690 for(y=0;y<out->ysize;y++)
691 {
692 /*
693 y is the coordinate in the new image.
694 yy is the corresponding x-value in the original size image.
695 yc is the integer value, the pixel coordinate of xx.
696 */
697 yy = (double) y / scale;
698 /* coordinate (0.0,0.0) is in the center of pixel (0,0),
699 so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */
700 yc = (int) floor( yy + 0.5 );
701 gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );
702 /* the kernel must be computed for each y because the fine
703 offset yy-yc is different in each case */
704
705 for(x=0;x<out->xsize;x++)
706 {
707 sum = 0.0;
708 for(i=0;i<kernel->dim;i++)
709 {
710 j = yc - h + i;
711
712 /* symmetry boundary condition */
713 while( j < 0 ) j += double_y_size;
714 while( j >= double_y_size ) j -= double_y_size;
715 if( j >= (int) in->ysize ) j = double_y_size-1-j;
716
717 sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];
718 }
719 out->data[ x + y * out->xsize ] = sum;
720 }
721 }
722
723 /* free memory */
724 free_ntuple_list(kernel);
725 free_image_double(aux);
726
727 return out;
728 }
729
730
731 /*----------------------------------------------------------------------------*/
732 /*--------------------------------- Gradient ---------------------------------*/
733 /*----------------------------------------------------------------------------*/
734
735 /*----------------------------------------------------------------------------*/
736 /** Computes the direction of the level line of 'in' at each point.
737
738 The result is:
739 - an image_double with the angle at each pixel, or NOTDEF if not defined.
740 - the image_double 'modgrad' (a pointer is passed as argument)
741 with the gradient magnitude at each point.
742 - a list of pixels 'list_p' roughly ordered by decreasing
743 gradient magnitude. (The order is made by classifying points
744 into bins by gradient magnitude. The parameters 'n_bins' and
745 'max_grad' specify the number of bins and the gradient modulus
746 at the highest bin. The pixels in the list would be in
747 decreasing gradient magnitude, up to a precision of the size of
748 the bins.)
749 - a pointer 'mem_p' to the memory used by 'list_p' to be able to
750 free the memory when it is not used anymore.
751 */
ll_angle(image_double in,double threshold,struct coorlist ** list_p,void ** mem_p,image_double * modgrad,unsigned int n_bins)752 static image_double ll_angle( image_double in, double threshold,
753 struct coorlist ** list_p, void ** mem_p,
754 image_double * modgrad, unsigned int n_bins )
755 {
756 image_double g;
757 unsigned int n,p,x,y,adr,i;
758 double com1,com2,gx,gy,norm,norm2;
759 /* the rest of the variables are used for pseudo-ordering
760 the gradient magnitude values */
761 int list_count = 0;
762 struct coorlist * list;
763 struct coorlist ** range_l_s; /* array of pointers to start of bin list */
764 struct coorlist ** range_l_e; /* array of pointers to end of bin list */
765 struct coorlist * start;
766 struct coorlist * end;
767 double max_grad = 0.0;
768
769 /* check parameters */
770 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
771 error("ll_angle: invalid image.");
772 if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");
773 if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");
774 if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");
775 if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");
776 if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");
777
778 /* image size shortcuts */
779 n = in->ysize;
780 p = in->xsize;
781
782 /* allocate output image */
783 g = new_image_double(in->xsize,in->ysize);
784
785 /* get memory for the image of gradient modulus */
786 *modgrad = new_image_double(in->xsize,in->ysize);
787
788 /* get memory for "ordered" list of pixels */
789 list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );
790 *mem_p = (void *) list;
791 range_l_s = (struct coorlist **) calloc( (size_t) n_bins,
792 sizeof(struct coorlist *) );
793 range_l_e = (struct coorlist **) calloc( (size_t) n_bins,
794 sizeof(struct coorlist *) );
795 if( list == NULL || range_l_s == NULL || range_l_e == NULL )
796 error("not enough memory.");
797 for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;
798
799 /* 'undefined' on the down and right boundaries */
800 for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;
801 for(y=0;y<n;y++) g->data[p*y+p-1] = NOTDEF;
802
803 /* compute gradient on the remaining pixels */
804 for(x=0;x<p-1;x++)
805 for(y=0;y<n-1;y++)
806 {
807 adr = y*p+x;
808
809 /*
810 Norm 2 computation using 2x2 pixel window:
811 A B
812 C D
813 and
814 com1 = D-A, com2 = B-C.
815 Then
816 gx = B+D - (A+C) horizontal difference
817 gy = C+D - (A+B) vertical difference
818 com1 and com2 are just to avoid 2 additions.
819 */
820 com1 = in->data[adr+p+1] - in->data[adr];
821 com2 = in->data[adr+1] - in->data[adr+p];
822
823 gx = com1+com2; /* gradient x component */
824 gy = com1-com2; /* gradient y component */
825 norm2 = gx*gx+gy*gy;
826 norm = sqrt( norm2 / 4.0 ); /* gradient norm */
827
828 (*modgrad)->data[adr] = norm; /* store gradient norm */
829
830 if( norm <= threshold ) /* norm too small, gradient no defined */
831 g->data[adr] = NOTDEF; /* gradient angle not defined */
832 else
833 {
834 /* gradient angle computation */
835 g->data[adr] = atan2(gx,-gy);
836
837 /* look for the maximum of the gradient */
838 if( norm > max_grad ) max_grad = norm;
839 }
840 }
841
842 /* compute histogram of gradient values */
843 for(x=0;x<p-1;x++)
844 for(y=0;y<n-1;y++)
845 {
846 norm = (*modgrad)->data[y*p+x];
847
848 /* store the point in the right bin according to its norm */
849 i = (unsigned int) (norm * (double) n_bins / max_grad);
850 if( i >= n_bins ) i = n_bins-1;
851 if( range_l_e[i] == NULL )
852 range_l_s[i] = range_l_e[i] = list+list_count++;
853 else
854 {
855 range_l_e[i]->next = list+list_count;
856 range_l_e[i] = list+list_count++;
857 }
858 range_l_e[i]->x = (int) x;
859 range_l_e[i]->y = (int) y;
860 range_l_e[i]->next = NULL;
861 }
862
863 /* Make the list of pixels (almost) ordered by norm value.
864 It starts by the larger bin, so the list starts by the
865 pixels with the highest gradient value. Pixels would be ordered
866 by norm value, up to a precision given by max_grad/n_bins.
867 */
868 for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);
869 start = range_l_s[i];
870 end = range_l_e[i];
871 if( start != NULL )
872 while(i>0)
873 {
874 --i;
875 if( range_l_s[i] != NULL )
876 {
877 end->next = range_l_s[i];
878 end = range_l_e[i];
879 }
880 }
881 *list_p = start;
882
883 /* free memory */
884 free( (void *) range_l_s );
885 free( (void *) range_l_e );
886
887 return g;
888 }
889
890 /*----------------------------------------------------------------------------*/
891 /** Is point (x,y) aligned to angle theta, up to precision 'prec'?
892 */
isaligned(int x,int y,image_double angles,double theta,double prec)893 static int isaligned( int x, int y, image_double angles, double theta,
894 double prec )
895 {
896 double a;
897
898 /* check parameters */
899 if( angles == NULL || angles->data == NULL )
900 error("isaligned: invalid image 'angles'.");
901 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
902 error("isaligned: (x,y) out of the image.");
903 if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");
904
905 /* angle at pixel (x,y) */
906 a = angles->data[ x + y * angles->xsize ];
907
908 /* pixels whose level-line angle is not defined
909 are considered as NON-aligned */
910 if( a == NOTDEF ) return FALSE; /* there is no need to call the function
911 'double_equal' here because there is
912 no risk of problems related to the
913 comparison doubles, we are only
914 interested in the exact NOTDEF value */
915
916 /* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */
917 theta -= a;
918 if( theta < 0.0 ) theta = -theta;
919 if( theta > M_3_2_PI )
920 {
921 theta -= M_2__PI;
922 if( theta < 0.0 ) theta = -theta;
923 }
924
925 return theta <= prec;
926 }
927
928 /*----------------------------------------------------------------------------*/
929 /** Absolute value angle difference.
930 */
angle_diff(double a,double b)931 static double angle_diff(double a, double b)
932 {
933 a -= b;
934 while( a <= -M_PI ) a += M_2__PI;
935 while( a > M_PI ) a -= M_2__PI;
936 if( a < 0.0 ) a = -a;
937 return a;
938 }
939
940 /*----------------------------------------------------------------------------*/
941 /** Signed angle difference.
942 */
angle_diff_signed(double a,double b)943 static double angle_diff_signed(double a, double b)
944 {
945 a -= b;
946 while( a <= -M_PI ) a += M_2__PI;
947 while( a > M_PI ) a -= M_2__PI;
948 return a;
949 }
950
951
952 /*----------------------------------------------------------------------------*/
953 /*----------------------------- NFA computation ------------------------------*/
954 /*----------------------------------------------------------------------------*/
955
956 /*----------------------------------------------------------------------------*/
957 /** Computes the natural logarithm of the absolute value of
958 the gamma function of x using the Lanczos approximation.
959 See http://www.rskey.org/gamma.htm
960
961 The formula used is
962 @f[
963 \Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) }
964 (x+5.5)^{x+0.5} e^{-(x+5.5)}
965 @f]
966 so
967 @f[
968 \log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right)
969 + (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n)
970 @f]
971 and
972 q0 = 75122.6331530,
973 q1 = 80916.6278952,
974 q2 = 36308.2951477,
975 q3 = 8687.24529705,
976 q4 = 1168.92649479,
977 q5 = 83.8676043424,
978 q6 = 2.50662827511.
979 */
log_gamma_lanczos(double x)980 static double log_gamma_lanczos(double x)
981 {
982 static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
983 8687.24529705, 1168.92649479, 83.8676043424,
984 2.50662827511 };
985 double a = (x+0.5) * log(x+5.5) - (x+5.5);
986 double b = 0.0;
987 int n;
988
989 for(n=0;n<7;n++)
990 {
991 a -= log( x + (double) n );
992 b += q[n] * pow( x, (double) n );
993 }
994 return a + log(b);
995 }
996
997 /*----------------------------------------------------------------------------*/
998 /** Computes the natural logarithm of the absolute value of
999 the gamma function of x using Windschitl method.
1000 See http://www.rskey.org/gamma.htm
1001
1002 The formula used is
1003 @f[
1004 \Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e}
1005 \sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x
1006 @f]
1007 so
1008 @f[
1009 \log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x
1010 + 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right).
1011 @f]
1012 This formula is a good approximation when x > 15.
1013 */
log_gamma_windschitl(double x)1014 static double log_gamma_windschitl(double x)
1015 {
1016 return 0.918938533204673 + (x-0.5)*log(x) - x
1017 + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
1018 }
1019
1020 /*----------------------------------------------------------------------------*/
1021 /** Computes the natural logarithm of the absolute value of
1022 the gamma function of x. When x>15 use log_gamma_windschitl(),
1023 otherwise use log_gamma_lanczos().
1024 */
1025 #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
1026
1027 /*----------------------------------------------------------------------------*/
1028 /** Size of the table to store already computed inverse values.
1029 */
1030 #define TABSIZE 100000
1031
1032 /*----------------------------------------------------------------------------*/
1033 /** Computes -log10(NFA).
1034
1035 NFA stands for Number of False Alarms:
1036 @f[
1037 \mathrm{NFA} = NT \cdot B(n,k,p)
1038 @f]
1039
1040 - NT - number of tests
1041 - B(n,k,p) - tail of binomial distribution with parameters n,k and p:
1042 @f[
1043 B(n,k,p) = \sum_{j=k}^n
1044 \left(\begin{array}{c}n\\j\end{array}\right)
1045 p^{j} (1-p)^{n-j}
1046 @f]
1047
1048 The value -log10(NFA) is equivalent but more intuitive than NFA:
1049 - -1 corresponds to 10 mean false alarms
1050 - 0 corresponds to 1 mean false alarm
1051 - 1 corresponds to 0.1 mean false alarms
1052 - 2 corresponds to 0.01 mean false alarms
1053 - ...
1054
1055 Used this way, the bigger the value, better the detection,
1056 and a logarithmic scale is used.
1057
1058 @param n,k,p binomial parameters.
1059 @param logNT logarithm of Number of Tests
1060
1061 The computation is based in the gamma function by the following
1062 relation:
1063 @f[
1064 \left(\begin{array}{c}n\\k\end{array}\right)
1065 = \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }.
1066 @f]
1067 We use efficient algorithms to compute the logarithm of
1068 the gamma function.
1069
1070 To make the computation faster, not all the sum is computed, part
1071 of the terms are neglected based on a bound to the error obtained
1072 (an error of 10% in the result is accepted).
1073 */
nfa(int n,int k,double p,double logNT)1074 static double nfa(int n, int k, double p, double logNT)
1075 {
1076 static double inv[TABSIZE]; /* table to keep computed inverse values */
1077 double tolerance = 0.1; /* an error of 10% in the result is accepted */
1078 double log1term,term,bin_term,mult_term,bin_tail,err,p_term;
1079 int i;
1080
1081 /* check parameters */
1082 if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )
1083 error("nfa: wrong n, k or p values.");
1084
1085 /* trivial cases */
1086 if( n==0 || k==0 ) return -logNT;
1087 if( n==k ) return -logNT - (double) n * log10(p);
1088
1089 /* probability term */
1090 p_term = p / (1.0-p);
1091
1092 /* compute the first term of the series */
1093 /*
1094 binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
1095 where bincoef(n,i) are the binomial coefficients.
1096 But
1097 bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
1098 We use this to compute the first term. Actually the log of it.
1099 */
1100 log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )
1101 - log_gamma( (double) (n-k) + 1.0 )
1102 + (double) k * log(p) + (double) (n-k) * log(1.0-p);
1103 term = exp(log1term);
1104
1105 /* in some cases no more computations are needed */
1106 if( double_equal(term,0.0) ) /* the first term is almost zero */
1107 {
1108 if( (double) k > (double) n * p ) /* at begin or end of the tail? */
1109 return -log1term / M_LN10 - logNT; /* end: use just the first term */
1110 else
1111 return -logNT; /* begin: the tail is roughly 1 */
1112 }
1113
1114 /* compute more terms if needed */
1115 bin_tail = term;
1116 for(i=k+1;i<=n;i++)
1117 {
1118 /*
1119 As
1120 term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
1121 and
1122 bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,
1123 then,
1124 term_i / term_i-1 = (n-i+1)/i * p/(1-p)
1125 and
1126 term_i = term_i-1 * (n-i+1)/i * p/(1-p).
1127 1/i is stored in a table as they are computed,
1128 because divisions are expensive.
1129 p/(1-p) is computed only once and stored in 'p_term'.
1130 */
1131 bin_term = (double) (n-i+1) * ( i<TABSIZE ?
1132 ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :
1133 1.0 / (double) i );
1134
1135 mult_term = bin_term * p_term;
1136 term *= mult_term;
1137 bin_tail += term;
1138 if(bin_term<1.0)
1139 {
1140 /* When bin_term<1 then mult_term_j<mult_term_i for j>i.
1141 Then, the error on the binomial tail when truncated at
1142 the i term can be bounded by a geometric series of form
1143 term_i * sum mult_term_i^j. */
1144 err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /
1145 (1.0-mult_term) - 1.0 );
1146
1147 /* One wants an error at most of tolerance*final_result, or:
1148 tolerance * abs(-log10(bin_tail)-logNT).
1149 Now, the error that can be accepted on bin_tail is
1150 given by tolerance*final_result divided by the derivative
1151 of -log10(x) when x=bin_tail. that is:
1152 tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
1153 Finally, we truncate the tail if the error is less than:
1154 tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */
1155 if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;
1156 }
1157 }
1158 return -log10(bin_tail) - logNT;
1159 }
1160
1161
1162 /*----------------------------------------------------------------------------*/
1163 /*--------------------------- Rectangle structure ----------------------------*/
1164 /*----------------------------------------------------------------------------*/
1165
1166 /*----------------------------------------------------------------------------*/
1167 /** Rectangle structure: line segment with width.
1168 */
1169 struct rect
1170 {
1171 double x1,y1,x2,y2; /* first and second point of the line segment */
1172 double width; /* rectangle width */
1173 double x,y; /* center of the rectangle */
1174 double theta; /* angle */
1175 double dx,dy; /* (dx,dy) is vector oriented as the line segment */
1176 double prec; /* tolerance angle */
1177 double p; /* probability of a point with angle within 'prec' */
1178 };
1179
1180 /*----------------------------------------------------------------------------*/
1181 /** Copy one rectangle structure to another.
1182 */
rect_copy(struct rect * in,struct rect * out)1183 static void rect_copy(struct rect * in, struct rect * out)
1184 {
1185 /* check parameters */
1186 if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");
1187
1188 /* copy values */
1189 out->x1 = in->x1;
1190 out->y1 = in->y1;
1191 out->x2 = in->x2;
1192 out->y2 = in->y2;
1193 out->width = in->width;
1194 out->x = in->x;
1195 out->y = in->y;
1196 out->theta = in->theta;
1197 out->dx = in->dx;
1198 out->dy = in->dy;
1199 out->prec = in->prec;
1200 out->p = in->p;
1201 }
1202
1203 /*----------------------------------------------------------------------------*/
1204 /** Rectangle points iterator.
1205
1206 The integer coordinates of pixels inside a rectangle are
1207 iteratively explored. This structure keep track of the process and
1208 functions ri_ini(), ri_inc(), ri_end(), and ri_del() are used in
1209 the process. An example of how to use the iterator is as follows:
1210 \code
1211
1212 struct rect * rec = XXX; // some rectangle
1213 rect_iter * i;
1214 for( i=ri_ini(rec); !ri_end(i); ri_inc(i) )
1215 {
1216 // your code, using 'i->x' and 'i->y' as coordinates
1217 }
1218 ri_del(i); // delete iterator
1219
1220 \endcode
1221 The pixels are explored 'column' by 'column', where we call
1222 'column' a set of pixels with the same x value that are inside the
1223 rectangle. The following is an schematic representation of a
1224 rectangle, the 'column' being explored is marked by colons, and
1225 the current pixel being explored is 'x,y'.
1226 \verbatim
1227
1228 vx[1],vy[1]
1229 * *
1230 * *
1231 * *
1232 * ye
1233 * : *
1234 vx[0],vy[0] : *
1235 * : *
1236 * x,y *
1237 * : *
1238 * : vx[2],vy[2]
1239 * : *
1240 y ys *
1241 ^ * *
1242 | * *
1243 | * *
1244 +---> x vx[3],vy[3]
1245
1246 \endverbatim
1247 The first 'column' to be explored is the one with the smaller x
1248 value. Each 'column' is explored starting from the pixel of the
1249 'column' (inside the rectangle) with the smallest y value.
1250
1251 The four corners of the rectangle are stored in order that rotates
1252 around the corners at the arrays 'vx[]' and 'vy[]'. The first
1253 point is always the one with smaller x value.
1254
1255 'x' and 'y' are the coordinates of the pixel being explored. 'ys'
1256 and 'ye' are the start and end values of the current column being
1257 explored. So, 'ys' < 'ye'.
1258 */
1259 typedef struct
1260 {
1261 double vx[4]; /* rectangle's corner X coordinates in circular order */
1262 double vy[4]; /* rectangle's corner Y coordinates in circular order */
1263 double ys,ye; /* start and end Y values of current 'column' */
1264 int x,y; /* coordinates of currently explored pixel */
1265 } rect_iter;
1266
1267 /*----------------------------------------------------------------------------*/
1268 /** Interpolate y value corresponding to 'x' value given, in
1269 the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller
1270 of 'y1' and 'y2'.
1271
1272 The following restrictions are required:
1273 - x1 <= x2
1274 - x1 <= x
1275 - x <= x2
1276 */
inter_low(double x,double x1,double y1,double x2,double y2)1277 static double inter_low(double x, double x1, double y1, double x2, double y2)
1278 {
1279 /* check parameters */
1280 if( x1 > x2 || x < x1 || x > x2 )
1281 error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
1282
1283 /* interpolation */
1284 if( double_equal(x1,x2) && y1<y2 ) return y1;
1285 if( double_equal(x1,x2) && y1>y2 ) return y2;
1286 return y1 + (x-x1) * (y2-y1) / (x2-x1);
1287 }
1288
1289 /*----------------------------------------------------------------------------*/
1290 /** Interpolate y value corresponding to 'x' value given, in
1291 the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger
1292 of 'y1' and 'y2'.
1293
1294 The following restrictions are required:
1295 - x1 <= x2
1296 - x1 <= x
1297 - x <= x2
1298 */
inter_hi(double x,double x1,double y1,double x2,double y2)1299 static double inter_hi(double x, double x1, double y1, double x2, double y2)
1300 {
1301 /* check parameters */
1302 if( x1 > x2 || x < x1 || x > x2 )
1303 error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
1304
1305 /* interpolation */
1306 if( double_equal(x1,x2) && y1<y2 ) return y2;
1307 if( double_equal(x1,x2) && y1>y2 ) return y1;
1308 return y1 + (x-x1) * (y2-y1) / (x2-x1);
1309 }
1310
1311 /*----------------------------------------------------------------------------*/
1312 /** Free memory used by a rectangle iterator.
1313 */
ri_del(rect_iter * iter)1314 static void ri_del(rect_iter * iter)
1315 {
1316 if( iter == NULL ) error("ri_del: NULL iterator.");
1317 free( (void *) iter );
1318 }
1319
1320 /*----------------------------------------------------------------------------*/
1321 /** Check if the iterator finished the full iteration.
1322
1323 See details in \ref rect_iter
1324 */
ri_end(rect_iter * i)1325 static int ri_end(rect_iter * i)
1326 {
1327 /* check input */
1328 if( i == NULL ) error("ri_end: NULL iterator.");
1329
1330 /* if the current x value is larger than the largest
1331 x value in the rectangle (vx[2]), we know the full
1332 exploration of the rectangle is finished. */
1333 return (double)(i->x) > i->vx[2];
1334 }
1335
1336 /*----------------------------------------------------------------------------*/
1337 /** Increment a rectangle iterator.
1338
1339 See details in \ref rect_iter
1340 */
ri_inc(rect_iter * i)1341 static void ri_inc(rect_iter * i)
1342 {
1343 /* check input */
1344 if( i == NULL ) error("ri_inc: NULL iterator.");
1345
1346 /* if not at end of exploration,
1347 increase y value for next pixel in the 'column' */
1348 if( !ri_end(i) ) i->y++;
1349
1350 /* if the end of the current 'column' is reached,
1351 and it is not the end of exploration,
1352 advance to the next 'column' */
1353 while( (double) (i->y) > i->ye && !ri_end(i) )
1354 {
1355 /* increase x, next 'column' */
1356 i->x++;
1357
1358 /* if end of exploration, return */
1359 if( ri_end(i) ) return;
1360
1361 /* update lower y limit (start) for the new 'column'.
1362
1363 We need to interpolate the y value that corresponds to the
1364 lower side of the rectangle. The first thing is to decide if
1365 the corresponding side is
1366
1367 vx[0],vy[0] to vx[3],vy[3] or
1368 vx[3],vy[3] to vx[2],vy[2]
1369
1370 Then, the side is interpolated for the x value of the
1371 'column'. But, if the side is vertical (as it could happen if
1372 the rectangle is vertical and we are dealing with the first
1373 or last 'columns') then we pick the lower value of the side
1374 by using 'inter_low'.
1375 */
1376 if( (double) i->x < i->vx[3] )
1377 i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);
1378 else
1379 i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);
1380
1381 /* update upper y limit (end) for the new 'column'.
1382
1383 We need to interpolate the y value that corresponds to the
1384 upper side of the rectangle. The first thing is to decide if
1385 the corresponding side is
1386
1387 vx[0],vy[0] to vx[1],vy[1] or
1388 vx[1],vy[1] to vx[2],vy[2]
1389
1390 Then, the side is interpolated for the x value of the
1391 'column'. But, if the side is vertical (as it could happen if
1392 the rectangle is vertical and we are dealing with the first
1393 or last 'columns') then we pick the lower value of the side
1394 by using 'inter_low'.
1395 */
1396 if( (double)i->x < i->vx[1] )
1397 i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);
1398 else
1399 i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);
1400
1401 /* new y */
1402 i->y = (int) ceil(i->ys);
1403 }
1404 }
1405
1406 /*----------------------------------------------------------------------------*/
1407 /** Create and initialize a rectangle iterator.
1408
1409 See details in \ref rect_iter
1410 */
ri_ini(struct rect * r)1411 static rect_iter * ri_ini(struct rect * r)
1412 {
1413 double vx[4],vy[4];
1414 int n,offset;
1415 rect_iter * i;
1416
1417 /* check parameters */
1418 if( r == NULL ) error("ri_ini: invalid rectangle.");
1419
1420 /* get memory */
1421 i = (rect_iter *) malloc(sizeof(rect_iter));
1422 if( i == NULL ) error("ri_ini: Not enough memory.");
1423
1424 /* build list of rectangle corners ordered
1425 in a circular way around the rectangle */
1426 vx[0] = r->x1 - r->dy * r->width / 2.0;
1427 vy[0] = r->y1 + r->dx * r->width / 2.0;
1428 vx[1] = r->x2 - r->dy * r->width / 2.0;
1429 vy[1] = r->y2 + r->dx * r->width / 2.0;
1430 vx[2] = r->x2 + r->dy * r->width / 2.0;
1431 vy[2] = r->y2 - r->dx * r->width / 2.0;
1432 vx[3] = r->x1 + r->dy * r->width / 2.0;
1433 vy[3] = r->y1 - r->dx * r->width / 2.0;
1434
1435 /* compute rotation of index of corners needed so that the first
1436 point has the smaller x.
1437
1438 if one side is vertical, thus two corners have the same smaller x
1439 value, the one with the largest y value is selected as the first.
1440 */
1441 if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;
1442 else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;
1443 else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;
1444 else offset = 3;
1445
1446 /* apply rotation of index. */
1447 for(n=0; n<4; n++)
1448 {
1449 i->vx[n] = vx[(offset+n)%4];
1450 i->vy[n] = vy[(offset+n)%4];
1451 }
1452
1453 /* Set an initial condition.
1454
1455 The values are set to values that will cause 'ri_inc' (that will
1456 be called immediately) to initialize correctly the first 'column'
1457 and compute the limits 'ys' and 'ye'.
1458
1459 'y' is set to the integer value of vy[0], the starting corner.
1460
1461 'ys' and 'ye' are set to very small values, so 'ri_inc' will
1462 notice that it needs to start a new 'column'.
1463
1464 The smallest integer coordinate inside of the rectangle is
1465 'ceil(vx[0])'. The current 'x' value is set to that value minus
1466 one, so 'ri_inc' (that will increase x by one) will advance to
1467 the first 'column'.
1468 */
1469 i->x = (int) ceil(i->vx[0]) - 1;
1470 i->y = (int) ceil(i->vy[0]);
1471 i->ys = i->ye = -DBL_MAX;
1472
1473 /* advance to the first pixel */
1474 ri_inc(i);
1475
1476 return i;
1477 }
1478
1479 /*----------------------------------------------------------------------------*/
1480 /** Compute a rectangle's NFA value.
1481 */
rect_nfa(struct rect * rec,image_double angles,double logNT)1482 static double rect_nfa(struct rect * rec, image_double angles, double logNT)
1483 {
1484 rect_iter * i;
1485 int pts = 0;
1486 int alg = 0;
1487
1488 /* check parameters */
1489 if( rec == NULL ) error("rect_nfa: invalid rectangle.");
1490 if( angles == NULL ) error("rect_nfa: invalid 'angles'.");
1491
1492 /* compute the total number of pixels and of aligned points in 'rec' */
1493 for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */
1494 if( i->x >= 0 && i->y >= 0 &&
1495 i->x < (int) angles->xsize && i->y < (int) angles->ysize )
1496 {
1497 ++pts; /* total number of pixels counter */
1498 if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )
1499 ++alg; /* aligned points counter */
1500 }
1501 ri_del(i); /* delete iterator */
1502
1503 return nfa(pts,alg,rec->p,logNT); /* compute NFA value */
1504 }
1505
1506
1507 /*----------------------------------------------------------------------------*/
1508 /*---------------------------------- Regions ---------------------------------*/
1509 /*----------------------------------------------------------------------------*/
1510
1511 /*----------------------------------------------------------------------------*/
1512 /** Compute region's angle as the principal inertia axis of the region.
1513
1514 The following is the region inertia matrix A:
1515 @f[
1516
1517 A = \left(\begin{array}{cc}
1518 Ixx & Ixy \\
1519 Ixy & Iyy \\
1520 \end{array}\right)
1521
1522 @f]
1523 where
1524
1525 Ixx = sum_i G(i).(y_i - cx)^2
1526
1527 Iyy = sum_i G(i).(x_i - cy)^2
1528
1529 Ixy = - sum_i G(i).(x_i - cx).(y_i - cy)
1530
1531 and
1532 - G(i) is the gradient norm at pixel i, used as pixel's weight.
1533 - x_i and y_i are the coordinates of pixel i.
1534 - cx and cy are the coordinates of the center of th region.
1535
1536 lambda1 and lambda2 are the eigenvalues of matrix A,
1537 with lambda1 >= lambda2. They are found by solving the
1538 characteristic polynomial:
1539
1540 det( lambda I - A) = 0
1541
1542 that gives:
1543
1544 lambda1 = ( Ixx + Iyy + sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2
1545
1546 lambda2 = ( Ixx + Iyy - sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2
1547
1548 To get the line segment direction we want to get the angle the
1549 eigenvector associated to the smallest eigenvalue. We have
1550 to solve for a,b in:
1551
1552 a.Ixx + b.Ixy = a.lambda2
1553
1554 a.Ixy + b.Iyy = b.lambda2
1555
1556 We want the angle theta = atan(b/a). It can be computed with
1557 any of the two equations:
1558
1559 theta = atan( (lambda2-Ixx) / Ixy )
1560
1561 or
1562
1563 theta = atan( Ixy / (lambda2-Iyy) )
1564
1565 When |Ixx| > |Iyy| we use the first, otherwise the second (just to
1566 get better numeric precision).
1567 */
get_theta(struct point * reg,int reg_size,double x,double y,image_double modgrad,double reg_angle,double prec)1568 static double get_theta( struct point * reg, int reg_size, double x, double y,
1569 image_double modgrad, double reg_angle, double prec )
1570 {
1571 double lambda,theta,weight;
1572 double Ixx = 0.0;
1573 double Iyy = 0.0;
1574 double Ixy = 0.0;
1575 int i;
1576
1577 /* check parameters */
1578 if( reg == NULL ) error("get_theta: invalid region.");
1579 if( reg_size <= 1 ) error("get_theta: region size <= 1.");
1580 if( modgrad == NULL || modgrad->data == NULL )
1581 error("get_theta: invalid 'modgrad'.");
1582 if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");
1583
1584 /* compute inertia matrix */
1585 for(i=0; i<reg_size; i++)
1586 {
1587 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
1588 Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;
1589 Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;
1590 Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;
1591 }
1592 if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )
1593 error("get_theta: null inertia matrix.");
1594
1595 /* compute smallest eigenvalue */
1596 lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );
1597
1598 /* compute angle */
1599 theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);
1600
1601 /* The previous procedure doesn't cares about orientation,
1602 so it could be wrong by 180 degrees. Here is corrected if necessary. */
1603 if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;
1604
1605 return theta;
1606 }
1607
1608 /*----------------------------------------------------------------------------*/
1609 /** Computes a rectangle that covers a region of points.
1610 */
region2rect(struct point * reg,int reg_size,image_double modgrad,double reg_angle,double prec,double p,struct rect * rec)1611 static void region2rect( struct point * reg, int reg_size,
1612 image_double modgrad, double reg_angle,
1613 double prec, double p, struct rect * rec )
1614 {
1615 double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;
1616 int i;
1617
1618 /* check parameters */
1619 if( reg == NULL ) error("region2rect: invalid region.");
1620 if( reg_size <= 1 ) error("region2rect: region size <= 1.");
1621 if( modgrad == NULL || modgrad->data == NULL )
1622 error("region2rect: invalid image 'modgrad'.");
1623 if( rec == NULL ) error("region2rect: invalid 'rec'.");
1624
1625 /* center of the region:
1626
1627 It is computed as the weighted sum of the coordinates
1628 of all the pixels in the region. The norm of the gradient
1629 is used as the weight of a pixel. The sum is as follows:
1630 cx = \sum_i G(i).x_i
1631 cy = \sum_i G(i).y_i
1632 where G(i) is the norm of the gradient of pixel i
1633 and x_i,y_i are its coordinates.
1634 */
1635 x = y = sum = 0.0;
1636 for(i=0; i<reg_size; i++)
1637 {
1638 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
1639 x += (double) reg[i].x * weight;
1640 y += (double) reg[i].y * weight;
1641 sum += weight;
1642 }
1643 if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");
1644 x /= sum;
1645 y /= sum;
1646
1647 /* theta */
1648 theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);
1649
1650 /* length and width:
1651
1652 'l' and 'w' are computed as the distance from the center of the
1653 region to pixel i, projected along the rectangle axis (dx,dy) and
1654 to the orthogonal axis (-dy,dx), respectively.
1655
1656 The length of the rectangle goes from l_min to l_max, where l_min
1657 and l_max are the minimum and maximum values of l in the region.
1658 Analogously, the width is selected from w_min to w_max, where
1659 w_min and w_max are the minimum and maximum of w for the pixels
1660 in the region.
1661 */
1662 dx = cos(theta);
1663 dy = sin(theta);
1664 l_min = l_max = w_min = w_max = 0.0;
1665 for(i=0; i<reg_size; i++)
1666 {
1667 l = ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;
1668 w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;
1669
1670 if( l > l_max ) l_max = l;
1671 if( l < l_min ) l_min = l;
1672 if( w > w_max ) w_max = w;
1673 if( w < w_min ) w_min = w;
1674 }
1675
1676 /* store values */
1677 rec->x1 = x + l_min * dx;
1678 rec->y1 = y + l_min * dy;
1679 rec->x2 = x + l_max * dx;
1680 rec->y2 = y + l_max * dy;
1681 rec->width = w_max - w_min;
1682 rec->x = x;
1683 rec->y = y;
1684 rec->theta = theta;
1685 rec->dx = dx;
1686 rec->dy = dy;
1687 rec->prec = prec;
1688 rec->p = p;
1689
1690 /* we impose a minimal width of one pixel
1691
1692 A sharp horizontal or vertical step would produce a perfectly
1693 horizontal or vertical region. The width computed would be
1694 zero. But that corresponds to a one pixels width transition in
1695 the image.
1696 */
1697 if( rec->width < 1.0 ) rec->width = 1.0;
1698 }
1699
1700 /*----------------------------------------------------------------------------*/
1701 /** Build a region of pixels that share the same angle, up to a
1702 tolerance 'prec', starting at point (x,y).
1703 */
region_grow(int x,int y,image_double angles,struct point * reg,int * reg_size,double * reg_angle,image_char used,double prec)1704 static void region_grow( int x, int y, image_double angles, struct point * reg,
1705 int * reg_size, double * reg_angle, image_char used,
1706 double prec )
1707 {
1708 double sumdx,sumdy;
1709 int xx,yy,i;
1710
1711 /* check parameters */
1712 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
1713 error("region_grow: (x,y) out of the image.");
1714 if( angles == NULL || angles->data == NULL )
1715 error("region_grow: invalid image 'angles'.");
1716 if( reg == NULL ) error("region_grow: invalid 'reg'.");
1717 if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");
1718 if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");
1719 if( used == NULL || used->data == NULL )
1720 error("region_grow: invalid image 'used'.");
1721
1722 /* first point of the region */
1723 *reg_size = 1;
1724 reg[0].x = x;
1725 reg[0].y = y;
1726 *reg_angle = angles->data[x+y*angles->xsize]; /* region's angle */
1727 sumdx = cos(*reg_angle);
1728 sumdy = sin(*reg_angle);
1729 used->data[x+y*used->xsize] = USED;
1730
1731 /* try neighbors as new region points */
1732 for(i=0; i<*reg_size; i++)
1733 for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)
1734 for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)
1735 if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&
1736 used->data[xx+yy*used->xsize] != USED &&
1737 isaligned(xx,yy,angles,*reg_angle,prec) )
1738 {
1739 /* add point */
1740 used->data[xx+yy*used->xsize] = USED;
1741 reg[*reg_size].x = xx;
1742 reg[*reg_size].y = yy;
1743 ++(*reg_size);
1744
1745 /* update region's angle */
1746 sumdx += cos( angles->data[xx+yy*angles->xsize] );
1747 sumdy += sin( angles->data[xx+yy*angles->xsize] );
1748 *reg_angle = atan2(sumdy,sumdx);
1749 }
1750 }
1751
1752 /*----------------------------------------------------------------------------*/
1753 /** Try some rectangles variations to improve NFA value. Only if the
1754 rectangle is not meaningful (i.e., log_nfa <= log_eps).
1755 */
rect_improve(struct rect * rec,image_double angles,double logNT,double log_eps)1756 static double rect_improve( struct rect * rec, image_double angles,
1757 double logNT, double log_eps )
1758 {
1759 struct rect r;
1760 double log_nfa,log_nfa_new;
1761 double delta = 0.5;
1762 double delta_2 = delta / 2.0;
1763 int n;
1764
1765 log_nfa = rect_nfa(rec,angles,logNT);
1766
1767 if( log_nfa > log_eps ) return log_nfa;
1768
1769 /* try finer precisions */
1770 rect_copy(rec,&r);
1771 for(n=0; n<5; n++)
1772 {
1773 r.p /= 2.0;
1774 r.prec = r.p * M_PI;
1775 log_nfa_new = rect_nfa(&r,angles,logNT);
1776 if( log_nfa_new > log_nfa )
1777 {
1778 log_nfa = log_nfa_new;
1779 rect_copy(&r,rec);
1780 }
1781 }
1782
1783 if( log_nfa > log_eps ) return log_nfa;
1784
1785 /* try to reduce width */
1786 rect_copy(rec,&r);
1787 for(n=0; n<5; n++)
1788 {
1789 if( (r.width - delta) >= 0.5 )
1790 {
1791 r.width -= delta;
1792 log_nfa_new = rect_nfa(&r,angles,logNT);
1793 if( log_nfa_new > log_nfa )
1794 {
1795 rect_copy(&r,rec);
1796 log_nfa = log_nfa_new;
1797 }
1798 }
1799 }
1800
1801 if( log_nfa > log_eps ) return log_nfa;
1802
1803 /* try to reduce one side of the rectangle */
1804 rect_copy(rec,&r);
1805 for(n=0; n<5; n++)
1806 {
1807 if( (r.width - delta) >= 0.5 )
1808 {
1809 r.x1 += -r.dy * delta_2;
1810 r.y1 += r.dx * delta_2;
1811 r.x2 += -r.dy * delta_2;
1812 r.y2 += r.dx * delta_2;
1813 r.width -= delta;
1814 log_nfa_new = rect_nfa(&r,angles,logNT);
1815 if( log_nfa_new > log_nfa )
1816 {
1817 rect_copy(&r,rec);
1818 log_nfa = log_nfa_new;
1819 }
1820 }
1821 }
1822
1823 if( log_nfa > log_eps ) return log_nfa;
1824
1825 /* try to reduce the other side of the rectangle */
1826 rect_copy(rec,&r);
1827 for(n=0; n<5; n++)
1828 {
1829 if( (r.width - delta) >= 0.5 )
1830 {
1831 r.x1 -= -r.dy * delta_2;
1832 r.y1 -= r.dx * delta_2;
1833 r.x2 -= -r.dy * delta_2;
1834 r.y2 -= r.dx * delta_2;
1835 r.width -= delta;
1836 log_nfa_new = rect_nfa(&r,angles,logNT);
1837 if( log_nfa_new > log_nfa )
1838 {
1839 rect_copy(&r,rec);
1840 log_nfa = log_nfa_new;
1841 }
1842 }
1843 }
1844
1845 if( log_nfa > log_eps ) return log_nfa;
1846
1847 /* try even finer precisions */
1848 rect_copy(rec,&r);
1849 for(n=0; n<5; n++)
1850 {
1851 r.p /= 2.0;
1852 r.prec = r.p * M_PI;
1853 log_nfa_new = rect_nfa(&r,angles,logNT);
1854 if( log_nfa_new > log_nfa )
1855 {
1856 log_nfa = log_nfa_new;
1857 rect_copy(&r,rec);
1858 }
1859 }
1860
1861 return log_nfa;
1862 }
1863
1864 /*----------------------------------------------------------------------------*/
1865 /** Reduce the region size, by elimination the points far from the
1866 starting point, until that leads to rectangle with the right
1867 density of region points or to discard the region if too small.
1868 */
reduce_region_radius(struct point * reg,int * reg_size,image_double modgrad,double reg_angle,double prec,double p,struct rect * rec,image_char used,image_double angles,double density_th)1869 static int reduce_region_radius( struct point * reg, int * reg_size,
1870 image_double modgrad, double reg_angle,
1871 double prec, double p, struct rect * rec,
1872 image_char used, image_double angles,
1873 double density_th )
1874 {
1875 double density,rad1,rad2,rad,xc,yc;
1876 int i;
1877
1878 /* check parameters */
1879 if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");
1880 if( reg_size == NULL )
1881 error("reduce_region_radius: invalid pointer 'reg_size'.");
1882 if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");
1883 if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");
1884 if( used == NULL || used->data == NULL )
1885 error("reduce_region_radius: invalid image 'used'.");
1886 if( angles == NULL || angles->data == NULL )
1887 error("reduce_region_radius: invalid image 'angles'.");
1888
1889 /* compute region points density */
1890 density = (double) *reg_size /
1891 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
1892
1893 /* if the density criterion is satisfied there is nothing to do */
1894 if( density >= density_th ) return TRUE;
1895
1896 /* compute region's radius */
1897 xc = (double) reg[0].x;
1898 yc = (double) reg[0].y;
1899 rad1 = dist( xc, yc, rec->x1, rec->y1 );
1900 rad2 = dist( xc, yc, rec->x2, rec->y2 );
1901 rad = rad1 > rad2 ? rad1 : rad2;
1902
1903 /* while the density criterion is not satisfied, remove farther pixels */
1904 while( density < density_th )
1905 {
1906 rad *= 0.75; /* reduce region's radius to 75% of its value */
1907
1908 /* remove points from the region and update 'used' map */
1909 for(i=0; i<*reg_size; i++)
1910 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )
1911 {
1912 /* point not kept, mark it as NOTUSED */
1913 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
1914 /* remove point from the region */
1915 reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */
1916 reg[i].y = reg[*reg_size-1].y;
1917 --(*reg_size);
1918 --i; /* to avoid skipping one point */
1919 }
1920
1921 /* reject if the region is too small.
1922 2 is the minimal region size for 'region2rect' to work. */
1923 if( *reg_size < 2 ) return FALSE;
1924
1925 /* re-compute rectangle */
1926 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
1927
1928 /* re-compute region points density */
1929 density = (double) *reg_size /
1930 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
1931 }
1932
1933 /* if this point is reached, the density criterion is satisfied */
1934 return TRUE;
1935 }
1936
1937 /*----------------------------------------------------------------------------*/
1938 /** Refine a rectangle.
1939
1940 For that, an estimation of the angle tolerance is performed by the
1941 standard deviation of the angle at points near the region's
1942 starting point. Then, a new region is grown starting from the same
1943 point, but using the estimated angle tolerance. If this fails to
1944 produce a rectangle with the right density of region points,
1945 'reduce_region_radius' is called to try to satisfy this condition.
1946 */
refine(struct point * reg,int * reg_size,image_double modgrad,double reg_angle,double prec,double p,struct rect * rec,image_char used,image_double angles,double density_th)1947 static int refine( struct point * reg, int * reg_size, image_double modgrad,
1948 double reg_angle, double prec, double p, struct rect * rec,
1949 image_char used, image_double angles, double density_th )
1950 {
1951 double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;
1952 int i,n;
1953
1954 /* check parameters */
1955 if( reg == NULL ) error("refine: invalid pointer 'reg'.");
1956 if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");
1957 if( prec < 0.0 ) error("refine: 'prec' must be positive.");
1958 if( rec == NULL ) error("refine: invalid pointer 'rec'.");
1959 if( used == NULL || used->data == NULL )
1960 error("refine: invalid image 'used'.");
1961 if( angles == NULL || angles->data == NULL )
1962 error("refine: invalid image 'angles'.");
1963
1964 /* compute region points density */
1965 density = (double) *reg_size /
1966 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
1967
1968 /* if the density criterion is satisfied there is nothing to do */
1969 if( density >= density_th ) return TRUE;
1970
1971 /*------ First try: reduce angle tolerance ------*/
1972
1973 /* compute the new mean angle and tolerance */
1974 xc = (double) reg[0].x;
1975 yc = (double) reg[0].y;
1976 ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];
1977 sum = s_sum = 0.0;
1978 n = 0;
1979 for(i=0; i<*reg_size; i++)
1980 {
1981 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
1982 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width )
1983 {
1984 angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];
1985 ang_d = angle_diff_signed(angle,ang_c);
1986 sum += ang_d;
1987 s_sum += ang_d * ang_d;
1988 ++n;
1989 }
1990 }
1991 mean_angle = sum / (double) n;
1992 tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n
1993 + mean_angle*mean_angle ); /* 2 * standard deviation */
1994
1995 /* find a new region from the same starting point and new angle tolerance */
1996 region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,®_angle,used,tau);
1997
1998 /* if the region is too small, reject */
1999 if( *reg_size < 2 ) return FALSE;
2000
2001 /* re-compute rectangle */
2002 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
2003
2004 /* re-compute region points density */
2005 density = (double) *reg_size /
2006 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
2007
2008 /*------ Second try: reduce region radius ------*/
2009 if( density < density_th )
2010 return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,
2011 rec, used, angles, density_th );
2012
2013 /* if this point is reached, the density criterion is satisfied */
2014 return TRUE;
2015 }
2016
2017
2018 /*----------------------------------------------------------------------------*/
2019 /*-------------------------- Line Segment Detector ---------------------------*/
2020 /*----------------------------------------------------------------------------*/
2021
2022 /*----------------------------------------------------------------------------*/
2023 /** LSD full interface.
2024 */
LineSegmentDetection(int * n_out,double * img,int X,int Y,double scale,double sigma_scale,double quant,double ang_th,double log_eps,double density_th,int n_bins,int ** reg_img,int * reg_x,int * reg_y)2025 double * LineSegmentDetection( int * n_out,
2026 double * img, int X, int Y,
2027 double scale, double sigma_scale, double quant,
2028 double ang_th, double log_eps, double density_th,
2029 int n_bins,
2030 int ** reg_img, int * reg_x, int * reg_y )
2031 {
2032 image_double image;
2033 ntuple_list out = new_ntuple_list(7);
2034 double * return_value;
2035 image_double scaled_image,angles,modgrad;
2036 image_char used;
2037 image_int region = NULL;
2038 struct coorlist * list_p;
2039 void * mem_p;
2040 struct rect rec;
2041 struct point * reg;
2042 int reg_size,min_reg_size,i;
2043 unsigned int xsize,ysize;
2044 double rho,reg_angle,prec,p,log_nfa,logNT;
2045 int ls_count = 0; /* line segments are numbered 1,2,3,... */
2046
2047
2048 /* check parameters */
2049 if( img == NULL || X <= 0 || Y <= 0 ) error("invalid image input.");
2050 if( scale <= 0.0 ) error("'scale' value must be positive.");
2051 if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");
2052 if( quant < 0.0 ) error("'quant' value must be positive.");
2053 if( ang_th <= 0.0 || ang_th >= 180.0 )
2054 error("'ang_th' value must be in the range (0,180).");
2055 if( density_th < 0.0 || density_th > 1.0 )
2056 error("'density_th' value must be in the range [0,1].");
2057 if( n_bins <= 0 ) error("'n_bins' value must be positive.");
2058
2059
2060 /* angle tolerance */
2061 prec = M_PI * ang_th / 180.0;
2062 p = ang_th / 180.0;
2063 rho = quant / sin(prec); /* gradient magnitude threshold */
2064
2065
2066 /* load and scale image (if necessary) and compute angle at each pixel */
2067 image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );
2068 if( scale != 1.0 )
2069 {
2070 scaled_image = gaussian_sampler( image, scale, sigma_scale );
2071 angles = ll_angle( scaled_image, rho, &list_p, &mem_p,
2072 &modgrad, (unsigned int) n_bins );
2073 free_image_double(scaled_image);
2074 }
2075 else
2076 angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,
2077 (unsigned int) n_bins );
2078 xsize = angles->xsize;
2079 ysize = angles->ysize;
2080
2081 /* Number of Tests - NT
2082
2083 The theoretical number of tests is Np.(XY)^(5/2)
2084 where X and Y are number of columns and rows of the image.
2085 Np corresponds to the number of angle precisions considered.
2086 As the procedure 'rect_improve' tests 5 times to halve the
2087 angle precision, and 5 more times after improving other factors,
2088 11 different precision values are potentially tested. Thus,
2089 the number of tests is
2090 11 * (X*Y)^(5/2)
2091 whose logarithm value is
2092 log10(11) + 5/2 * (log10(X) + log10(Y)).
2093 */
2094 logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0
2095 + log10(11.0);
2096 min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region
2097 that can give a meaningful event */
2098
2099
2100 /* initialize some structures */
2101 if( reg_img != NULL && reg_x != NULL && reg_y != NULL ) /* save region data */
2102 region = new_image_int_ini(angles->xsize,angles->ysize,0);
2103 used = new_image_char_ini(xsize,ysize,NOTUSED);
2104 reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );
2105 if( reg == NULL ) error("not enough memory!");
2106
2107
2108 /* search for line segments */
2109 for(; list_p != NULL; list_p = list_p->next )
2110 if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&
2111 angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )
2112 /* there is no risk of double comparison problems here
2113 because we are only interested in the exact NOTDEF value */
2114 {
2115 /* find the region of connected point and ~equal angle */
2116 region_grow( list_p->x, list_p->y, angles, reg, ®_size,
2117 ®_angle, used, prec );
2118
2119 /* reject small regions */
2120 if( reg_size < min_reg_size ) continue;
2121
2122 /* construct rectangular approximation for the region */
2123 region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);
2124
2125 /* Check if the rectangle exceeds the minimal density of
2126 region points. If not, try to improve the region.
2127 The rectangle will be rejected if the final one does
2128 not fulfill the minimal density condition.
2129 This is an addition to the original LSD algorithm published in
2130 "LSD: A Fast Line Segment Detector with a False Detection Control"
2131 by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall.
2132 The original algorithm is obtained with density_th = 0.0.
2133 */
2134 if( !refine( reg, ®_size, modgrad, reg_angle,
2135 prec, p, &rec, used, angles, density_th ) ) continue;
2136
2137 /* compute NFA value */
2138 log_nfa = rect_improve(&rec,angles,logNT,log_eps);
2139 if( log_nfa <= log_eps ) continue;
2140
2141 /* A New Line Segment was found! */
2142 ++ls_count; /* increase line segment counter */
2143
2144 /*
2145 The gradient was computed with a 2x2 mask, its value corresponds to
2146 points with an offset of (0.5,0.5), that should be added to output.
2147 The coordinates origin is at the center of pixel (0,0).
2148 */
2149 rec.x1 += 0.5; rec.y1 += 0.5;
2150 rec.x2 += 0.5; rec.y2 += 0.5;
2151
2152 /* scale the result values if a subsampling was performed */
2153 if( scale != 1.0 )
2154 {
2155 rec.x1 /= scale; rec.y1 /= scale;
2156 rec.x2 /= scale; rec.y2 /= scale;
2157 rec.width /= scale;
2158 }
2159
2160 /* add line segment found to output */
2161 add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,
2162 rec.width, rec.p, log_nfa );
2163
2164 /* add region number to 'region' image if needed */
2165 if( region != NULL )
2166 for(i=0; i<reg_size; i++)
2167 region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;
2168 }
2169
2170
2171 /* free memory */
2172 free( (void *) image ); /* only the double_image structure should be freed,
2173 the data pointer was provided to this functions
2174 and should not be destroyed. */
2175 free_image_double(angles);
2176 free_image_double(modgrad);
2177 free_image_char(used);
2178 free( (void *) reg );
2179 free( (void *) mem_p );
2180
2181 /* return the result */
2182 if( reg_img != NULL && reg_x != NULL && reg_y != NULL )
2183 {
2184 if( region == NULL ) error("'region' should be a valid image.");
2185 *reg_img = region->data;
2186 if( region->xsize > (unsigned int) INT_MAX ||
2187 region->xsize > (unsigned int) INT_MAX )
2188 error("region image to big to fit in INT sizes.");
2189 *reg_x = (int) (region->xsize);
2190 *reg_y = (int) (region->ysize);
2191
2192 /* free the 'region' structure.
2193 we cannot use the function 'free_image_int' because we need to keep
2194 the memory with the image data to be returned by this function. */
2195 free( (void *) region );
2196 }
2197 if( out->size > (unsigned int) INT_MAX )
2198 error("too many detections to fit in an INT.");
2199 *n_out = (int) (out->size);
2200
2201 return_value = out->values;
2202 free( (void *) out ); /* only the 'ntuple_list' structure must be freed,
2203 but the 'values' pointer must be keep to return
2204 as a result. */
2205
2206 return return_value;
2207 }
2208
2209 /*----------------------------------------------------------------------------*/
2210 /** LSD Simple Interface with Scale and Region output.
2211 */
lsd_scale_region(int * n_out,double * img,int X,int Y,double scale,int ** reg_img,int * reg_x,int * reg_y)2212 double * lsd_scale_region( int * n_out,
2213 double * img, int X, int Y, double scale,
2214 int ** reg_img, int * reg_x, int * reg_y )
2215 {
2216 /* LSD parameters */
2217 double sigma_scale = 0.6; /* Sigma for Gaussian filter is computed as
2218 sigma = sigma_scale/scale. */
2219 double quant = 2.0; /* Bound to the quantization error on the
2220 gradient norm. */
2221 double ang_th = 22.5; /* Gradient angle tolerance in degrees. */
2222 double log_eps = 0.0; /* Detection threshold: -log10(NFA) > log_eps */
2223 double density_th = 0.7; /* Minimal density of region points in rectangle. */
2224 int n_bins = 1024; /* Number of bins in pseudo-ordering of gradient
2225 modulus. */
2226
2227 return LineSegmentDetection( n_out, img, X, Y, scale, sigma_scale, quant,
2228 ang_th, log_eps, density_th, n_bins,
2229 reg_img, reg_x, reg_y );
2230 }
2231
2232 /*----------------------------------------------------------------------------*/
2233 /** LSD Simple Interface with Scale.
2234 */
lsd_scale(int * n_out,double * img,int X,int Y,double scale)2235 double * lsd_scale(int * n_out, double * img, int X, int Y, double scale)
2236 {
2237 return lsd_scale_region(n_out,img,X,Y,scale,NULL,NULL,NULL);
2238 }
2239
2240 /*----------------------------------------------------------------------------*/
2241 /** LSD Simple Interface.
2242 */
lsd(int * n_out,double * img,int X,int Y)2243 double * lsd(int * n_out, double * img, int X, int Y)
2244 {
2245 /* LSD parameters */
2246 double scale = 0.8; /* Scale the image by Gaussian filter to 'scale'. */
2247
2248 return lsd_scale(n_out,img,X,Y,scale);
2249 }
2250 /*----------------------------------------------------------------------------*/
2251