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