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