"Fossies" - the Fresh Open Source Software Archive

Member "darktable-3.6.1/src/iop/ashift_lsd.c" (10 Sep 2021, 80711 Bytes) of package /linux/misc/darktable-3.6.1.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ashift_lsd.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 3.4.1.1_vs_3.6.0.

    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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
  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  */
 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  */
 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 
 1069 __attribute__((constructor)) static void invConstructor()
 1070 {
 1071   if(inv) return;
 1072   inv = malloc(sizeof(double) * TABSIZE);
 1073 }
 1074 
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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  */
 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,&reg_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
 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, &reg_size,
 2172                      &reg_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, &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;