1/*2 This file is part of darktable,3 Copyright (C) 2016-2025 darktable developers.45 darktable is free software: you can redistribute it and/or modify6 it under the terms of the GNU General Public License as published by7 the Free Software Foundation, either version 3 of the License, or8 (at your option) any later version.910 darktable is distributed in the hope that it will be useful,11 but WITHOUT ANY WARRANTY; without even the implied warranty of12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13 GNU General Public License for more details.1415 You should have received a copy of the GNU General Public License16 along with darktable. If not, see <http://www.gnu.org/licenses/>.17*/181920/* 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 static24 * comment out unused interface functions25 * catch (unlikely) division by zero near line 203526 * rename rad1 and rad2 to radius1 and radius2 in reduce_region_radius()27 * to avoid naming conflict in windows build28 *29 */3031/* TODO: LSD employs intensive sanity checks of input data and allocation32 * results. In case of errors it will terminate the program. On the one side33 * this is not optimal for a module within darktable - it should better34 * report errors upstream where they should be handled properly. On the other35 * hand the kind of error which could be triggered in LSD would make it very unlikely36 * that the darktable process could survive anyhow.37 */3839// clang-format off4041/*==================================================================================42 * begin LSD code version 1.6. downloaded on January 30, 201643 *==================================================================================*/4445/*----------------------------------------------------------------------------4647 LSD - Line Segment Detector on digital images4849 This code is part of the following publication and was subject50 to peer review:5152 "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-lsd55 http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd5657 Copyright (c) 2007-2011 rafael grompone von gioi <grompone@gmail.com>5859 This program is free software: you can redistribute it and/or modify60 it under the terms of the GNU Affero General Public License as61 published by the Free Software Foundation, either version 3 of the62 License, or (at your option) any later version.6364 This program is distributed in the hope that it will be useful,65 but WITHOUT ANY WARRANTY; without even the implied warranty of66 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the67 GNU Affero General Public License for more details.6869 You should have received a copy of the GNU Affero General Public License70 along with this program. If not, see <http://www.gnu.org/licenses/>.7172 ----------------------------------------------------------------------------*/7374/*----------------------------------------------------------------------------*/75/** @file lsd.c76 LSD module code77 @author rafael grompone von gioi <grompone@gmail.com>78 */79/*----------------------------------------------------------------------------*/8081/*----------------------------------------------------------------------------*/82/** @mainpage LSD code documentation8384 This is an implementation of the Line Segment Detector described85 in the paper:8687 "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 and90 Machine Intelligence, vol. 32, no. 4, pp. 722-732, April, 2010.9192 and in more details in the CMLA Technical Report:9394 "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.9798 The version implemented here includes some further improvements99 described in the following publication, of which this code is part:100101 "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-lsd104 http://dx.doi.org/10.5201/ipol.2012.gjmr-lsd105106 The module's main function is lsd().107108 The source code is contained in two files: lsd.h and lsd.c.109110 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 test115 to consider the different precision values116 tested,117 - a minor bug corrected in the gradient sorting118 code,119 - the algorithm now also returns p and log_nfa120 for each detection,121 - a minor bug was corrected in the image scaling,122 - the angle comparison in "isaligned" changed123 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 and133 correction to partially handle "angle problem".134 - version 1.0 - jan 2009: First complete Megawave2 and Ansi C Language135 version.136137 @author rafael grompone von gioi <grompone@gmail.com>138 */139/*----------------------------------------------------------------------------*/140141142143144145146147148149150151152/** Label for pixels with undefined gradient. */153154155/** 3/2 pi */156157158/** 2 pi */159160161/** Label for pixels not used in yet. */162163164/** Label for pixels already used in detection. */165166167/*----------------------------------------------------------------------------*/168/** Chained list of coordinates.169 */170struct coorlist171{172 int x,y;173 struct coorlist * next;174};175176/*----------------------------------------------------------------------------*/177/** A point (or pixel).178 */179struct point {int x,y;};180181182/*----------------------------------------------------------------------------*/183/*------------------------- Miscellaneous functions --------------------------*/184/*----------------------------------------------------------------------------*/185186/*----------------------------------------------------------------------------*/187/** Fatal error, print a message to standard-error output and exit.188 */189static void error(char * msg)190{191 dt_print(DT_DEBUG_ALWAYS,"LSD Error: %s",msg);192 exit(EXIT_FAILURE);193}194195/*----------------------------------------------------------------------------*/196/** Doubles relative error factor197 */198199200/*----------------------------------------------------------------------------*/201/** Compare doubles by relative error.202203 The resulting rounding error after floating point computations204 depend on the specific operations done. The same number computed by205 different algorithms could present different rounding errors. For a206 useful comparison, an estimation of the relative rounding error207 should be considered and compared to a factor times EPS. The factor208 should be related to the cumulated rounding error in the chain of209 computation. Here, as a simplification, a fixed factor is used.210 */211static int double_equal(const double a, const double b)212{213 double abs_diff,aa,bb,abs_max;214215 /* trivial case */216 if( a == b ) return TRUE;217218 abs_diff = fabs(a-b);219 aa = fabs(a);220 bb = fabs(b);221 abs_max = aa > bb ? aa : bb;222223 /* DBL_MIN is the smallest normalized number, thus, the smallest224 number whose relative error is bounded by DBL_EPSILON. For225 smaller numbers, the same quantization steps as for DBL_MIN226 are used. Then, for smaller numbers, a meaningful "relative"227 error should be computed by dividing the difference by DBL_MIN. */228 if( abs_max < DBL_MIN ) abs_max = DBL_MIN;229230 /* equal if relative error <= factor x eps */231 return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);232}233234/*----------------------------------------------------------------------------*/235/** Computes Euclidean distance between point (x1,y1) and point (x2,y2).236 */237static double dist(const double x1, const double y1, const double x2, const double y2)238{239 return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );240}241242243/*----------------------------------------------------------------------------*/244/*----------------------- 'list of n-tuple' data type ------------------------*/245/*----------------------------------------------------------------------------*/246247/*----------------------------------------------------------------------------*/248/** 'list of n-tuple' data type249250 The i-th component of the j-th n-tuple of an n-tuple list 'ntl'251 is accessed with:252253 ntl->values[ i + j * ntl->dim ]254255 The dimension of the n-tuple (n) is:256257 ntl->dim258259 The number of n-tuples in the list is:260261 ntl->size262263 The maximum number of n-tuples that can be stored in the264 list with the allocated memory at a given time is given by:265266 ntl->max_size267 */268typedef struct ntuple_list_s269{270 unsigned int size;271 unsigned int max_size;272 unsigned int dim;273 double * values;274} * ntuple_list;275276/*----------------------------------------------------------------------------*/277/** Free memory used in n-tuple 'in'.278 */279static void free_ntuple_list(const ntuple_list in)280{281 if( in == NULL || in->values == NULL )282 error("free_ntuple_list: invalid n-tuple input.");283 free( (void *) in->values );284 free( (void *) in );285}286287/*----------------------------------------------------------------------------*/288/** Create an n-tuple list and allocate memory for one element.289 @param dim the dimension (n) of the n-tuple.290 */291static ntuple_list new_ntuple_list(const unsigned int dim)292{293 ntuple_list n_tuple;294295 /* check parameters */296 if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");297298 /* get memory for list structure */299 n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );300 if( n_tuple == NULL ) error("not enough memory.");301302 /* initialize list */303 n_tuple->size = 0;304 n_tuple->max_size = 1;305 n_tuple->dim = dim;306307 /* get memory for tuples */308 n_tuple->values = (double *) malloc(sizeof(double) * dim * n_tuple->max_size);309 if( n_tuple->values == NULL ) error("not enough memory.");310311 return n_tuple;312}313314/*----------------------------------------------------------------------------*/315/** Enlarge the allocated memory of an n-tuple list.316 */317static void enlarge_ntuple_list(const ntuple_list n_tuple)318{319 /* check parameters */320 if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )321 error("enlarge_ntuple_list: invalid n-tuple.");322323 /* duplicate number of tuples */324 n_tuple->max_size *= 2;325326 /* realloc memory */327 n_tuple->values = (double *) realloc( (void *) n_tuple->values,328 sizeof(double) * n_tuple->dim * n_tuple->max_size );329 if( n_tuple->values == NULL ) error("not enough memory.");330}331332/*----------------------------------------------------------------------------*/333/** Add a 7-tuple to an n-tuple list.334 */335static void add_7tuple(const ntuple_list out, const double v1, const double v2, const double v3,336 const double v4, const double v5, const double v6, const double v7 )337{338 /* check parameters */339 if( out == NULL ) error("add_7tuple: invalid n-tuple input.");340 if( out->dim != 7 ) error("add_7tuple: the n-tuple must be a 7-tuple.");341342 /* if needed, alloc more tuples to 'out' */343 if( out->size == out->max_size ) enlarge_ntuple_list(out);344 if( out->values == NULL ) error("add_7tuple: invalid n-tuple input.");345346 /* add new 7-tuple */347 out->values[ out->size * out->dim + 0 ] = v1;348 out->values[ out->size * out->dim + 1 ] = v2;349 out->values[ out->size * out->dim + 2 ] = v3;350 out->values[ out->size * out->dim + 3 ] = v4;351 out->values[ out->size * out->dim + 4 ] = v5;352 out->values[ out->size * out->dim + 5 ] = v6;353 out->values[ out->size * out->dim + 6 ] = v7;354355 /* update number of tuples counter */356 out->size++;357}358359360/*----------------------------------------------------------------------------*/361/*----------------------------- Image Data Types -----------------------------*/362/*----------------------------------------------------------------------------*/363364/*----------------------------------------------------------------------------*/365/** char image data type366367 The pixel value at (x,y) is accessed by:368369 image->data[ x + y * image->xsize ]370371 with x and y integer.372 */373typedef struct image_char_s374{375 unsigned char * data;376 unsigned int xsize,ysize;377} * image_char;378379/*----------------------------------------------------------------------------*/380/** Free memory used in image_char 'i'.381 */382static void free_image_char(const image_char i)383{384 if( i == NULL || i->data == NULL )385 error("free_image_char: invalid input image.");386 free( (void *) i->data );387 free( (void *) i );388}389390/*----------------------------------------------------------------------------*/391/** Create a new image_char of size 'xsize' times 'ysize'.392 */393static image_char new_image_char(const unsigned int xsize, const unsigned int ysize)394{395 image_char image;396397 /* check parameters */398 if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");399400 /* get memory */401 image = (image_char) malloc( sizeof(struct image_char_s) );402 if( image == NULL ) error("not enough memory.");403 image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),404 sizeof(unsigned char) );405 if( image->data == NULL ) error("not enough memory.");406407 /* set image size */408 image->xsize = xsize;409 image->ysize = ysize;410411 return image;412}413414/*----------------------------------------------------------------------------*/415/** Create a new image_char of size 'xsize' times 'ysize',416 initialized to the value 'fill_value'.417 */418static image_char new_image_char_ini(const unsigned int xsize, const unsigned int ysize,419 const unsigned char fill_value )420{421 const image_char image = new_image_char(xsize,ysize); /* create image */422 const unsigned int N = xsize*ysize;423 unsigned int i;424425 /* check parameters */426 if( image == NULL || image->data == NULL )427 error("new_image_char_ini: invalid image.");428429 /* initialize */430 for(i=0; i<N; i++) image->data[i] = fill_value;431432 return image;433}434435/*----------------------------------------------------------------------------*/436/** int image data type437438 The pixel value at (x,y) is accessed by:439440 image->data[ x + y * image->xsize ]441442 with x and y integer.443 */444typedef struct image_int_s445{446 int * data;447 unsigned int xsize,ysize;448} * image_int;449450/*----------------------------------------------------------------------------*/451/** Create a new image_int of size 'xsize' times 'ysize'.452 */453static image_int new_image_int(const unsigned int xsize, const unsigned int ysize)454{455 image_int image;456457 /* check parameters */458 if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");459460 /* get memory */461 image = (image_int) malloc( sizeof(struct image_int_s) );462 if( image == NULL ) error("not enough memory.");463 image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );464 if( image->data == NULL ) error("not enough memory.");465466 /* set image size */467 image->xsize = xsize;468 image->ysize = ysize;469470 return image;471}472473/*----------------------------------------------------------------------------*/474/** Create a new image_int of size 'xsize' times 'ysize',475 initialized to the value 'fill_value'.476 */477static image_int new_image_int_ini(const unsigned int xsize, const unsigned int ysize,478 const int fill_value )479{480 const image_int image = new_image_int(xsize,ysize); /* create image */481 const unsigned int N = xsize*ysize;482 unsigned int i;483484 /* initialize */485 for(i=0; i<N; i++) image->data[i] = fill_value;486487 return image;488}489490/*----------------------------------------------------------------------------*/491/** double image data type492493 The pixel value at (x,y) is accessed by:494495 image->data[ x + y * image->xsize ]496497 with x and y integer.498 */499typedef struct image_double_s500{501 double * data;502 unsigned int xsize,ysize;503} * image_double;504505/*----------------------------------------------------------------------------*/506/** Free memory used in image_double 'i'.507 */508static void free_image_double(const image_double i)509{510 if( i == NULL || i->data == NULL )511 error("free_image_double: invalid input image.");512 free( (void *) i->data );513 free( (void *) i );514}515516/*----------------------------------------------------------------------------*/517/** Create a new image_double of size 'xsize' times 'ysize'.518 */519static image_double new_image_double(const unsigned int xsize, const unsigned int ysize)520{521 image_double image;522523 /* check parameters */524 if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");525526 /* get memory */527 image = (image_double) malloc( sizeof(struct image_double_s) );528 if( image == NULL ) error("not enough memory.");529 image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );530 if( image->data == NULL ) error("not enough memory.");531532 /* set image size */533 image->xsize = xsize;534 image->ysize = ysize;535536 return image;537}538539/*----------------------------------------------------------------------------*/540/** Create a new image_double of size 'xsize' times 'ysize'541 with the data pointed by 'data'.542 */543static image_double new_image_double_ptr(const unsigned int xsize,544 const unsigned int ysize, double * data )545{546 image_double image;547548 /* check parameters */549 if( xsize == 0 || ysize == 0 )550 error("new_image_double_ptr: invalid image size.");551 if( data == NULL ) error("new_image_double_ptr: NULL data pointer.");552553 /* get memory */554 image = (image_double) malloc( sizeof(struct image_double_s) );555 if( image == NULL ) error("not enough memory.");556557 /* set image */558 image->xsize = xsize;559 image->ysize = ysize;560 image->data = data;561562 return image;563}564565566/*----------------------------------------------------------------------------*/567/*----------------------------- Gaussian filter ------------------------------*/568/*----------------------------------------------------------------------------*/569570/*----------------------------------------------------------------------------*/571/** Compute a Gaussian kernel of length 'kernel->dim',572 standard deviation 'sigma', and centered at value 'mean'.573574 For example, if mean=0.5, the Gaussian will be centered575 in the middle point between values 'kernel->values[0]'576 and 'kernel->values[1]'.577 */578static void gaussian_kernel(const ntuple_list kernel, const double sigma, const double mean)579{580 double sum = 0.0;581 double val;582 unsigned int i;583584 /* check parameters */585 if( kernel == NULL || kernel->values == NULL )586 error("gaussian_kernel: invalid n-tuple 'kernel'.");587 if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");588589 /* compute Gaussian kernel */590 if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);591 kernel->size = 1;592 for(i=0;i<kernel->dim;i++)593 {594 val = ( (double) i - mean ) / sigma;595 kernel->values[i] = exp( -0.5 * val * val );596 sum += kernel->values[i];597 }598599 /* normalization */600 if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;601}602603/*----------------------------------------------------------------------------*/604/** Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.605606 For example, scale=0.8 will give a result at 80% of the original size.607608 The image is convolved with a Gaussian kernel609 @f[610 G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}611 @f]612 before the sub-sampling to prevent aliasing.613614 The standard deviation sigma given by:615 - sigma = sigma_scale / scale, if scale < 1.0616 - sigma = sigma_scale, if scale >= 1.0617618 To be able to sub-sample at non-integer steps, some interpolation619 is needed. In this implementation, the interpolation is done by620 the Gaussian kernel, so both operations (filtering and sampling)621 are done at the same time. The Gaussian kernel is computed622 centered on the coordinates of the required sample. In this way,623 when applied, it gives directly the result of convolving the image624 with the kernel and interpolated to that particular position.625626 A fast algorithm is done using the separability of the Gaussian627 kernel. Applying the 2D Gaussian kernel is equivalent to applying628 first a horizontal 1D Gaussian kernel and then a vertical 1D629 Gaussian kernel (or the other way round). The reason is that630 @f[631 G(x,y) = G(x) * G(y)632 @f]633 where634 @f[635 G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}}.636 @f]637 The algorithm first applies a combined Gaussian kernel and sampling638 in the x axis, and then the combined Gaussian kernel and sampling639 in the y axis.640 */641static image_double gaussian_sampler(const image_double in, const double scale,642 const double sigma_scale )643{644 image_double aux,out;645 ntuple_list kernel;646 unsigned int N,M,h,n,x,y,i;647 int xc,yc,j,double_x_size,double_y_size;648 double sigma,xx,yy,sum,prec;649650 /* check parameters */651 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )652 error("gaussian_sampler: invalid image.");653 if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");654 if( sigma_scale <= 0.0 )655 error("gaussian_sampler: 'sigma_scale' must be positive.");656657 /* compute new image size and get memory for images */658 if( in->xsize * scale > (double) UINT_MAX ||659 in->ysize * scale > (double) UINT_MAX )660 error("gaussian_sampler: the output image size exceeds the handled size.");661 N = (unsigned int) ceil( in->xsize * scale );662 M = (unsigned int) ceil( in->ysize * scale );663 aux = new_image_double(N,in->ysize);664 out = new_image_double(N,M);665666 /* sigma, kernel size and memory for the kernel */667 sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;668 /*669 The size of the kernel is selected to guarantee that the670 the first discarded term is at least 10^prec times smaller671 than the central value. For that, h should be larger than x, with672 e^(-x^2/2sigma^2) = 1/10^prec.673 Then,674 x = sigma * sqrt( 2 * prec * ln(10) ).675 */676 prec = 3.0;677 h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );678 n = 1+2*h; /* kernel size */679 kernel = new_ntuple_list(n);680681 /* auxiliary double image size variables */682 double_x_size = (int) (2 * in->xsize);683 double_y_size = (int) (2 * in->ysize);684685 /* First subsampling: x axis */686 for(x=0;x<aux->xsize;x++)687 {688 /*689 x is the coordinate in the new image.690 xx is the corresponding x-value in the original size image.691 xc is the integer value, the pixel coordinate of xx.692 */693 xx = (double) x / scale;694 /* coordinate (0.0,0.0) is in the center of pixel (0,0),695 so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */696 xc = (int) floor( xx + 0.5 );697 gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );698 /* the kernel must be computed for each x because the fine699 offset xx-xc is different in each case */700701 for(y=0;y<aux->ysize;y++)702 {703 sum = 0.0;704 for(i=0;i<kernel->dim;i++)705 {706 j = xc - h + i;707708 /* symmetry boundary condition */709 while( j < 0 ) j += double_x_size;710 while( j >= double_x_size ) j -= double_x_size;711 if( j >= (int) in->xsize ) j = double_x_size-1-j;712713 sum += in->data[ j + y * in->xsize ] * kernel->values[i];714 }715 aux->data[ x + y * aux->xsize ] = sum;716 }717 }718719 /* Second subsampling: y axis */720 for(y=0;y<out->ysize;y++)721 {722 /*723 y is the coordinate in the new image.724 yy is the corresponding x-value in the original size image.725 yc is the integer value, the pixel coordinate of xx.726 */727 yy = (double) y / scale;728 /* coordinate (0.0,0.0) is in the center of pixel (0,0),729 so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */730 yc = (int) floor( yy + 0.5 );731 gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );732 /* the kernel must be computed for each y because the fine733 offset yy-yc is different in each case */734735 for(x=0;x<out->xsize;x++)736 {737 sum = 0.0;738 for(i=0;i<kernel->dim;i++)739 {740 j = yc - h + i;741742 /* symmetry boundary condition */743 while( j < 0 ) j += double_y_size;744 while( j >= double_y_size ) j -= double_y_size;745 if( j >= (int) in->ysize ) j = double_y_size-1-j;746747 sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];748 }749 out->data[ x + y * out->xsize ] = sum;750 }751 }752753 /* free memory */754 free_ntuple_list(kernel);755 free_image_double(aux);756757 return out;758}759760761/*----------------------------------------------------------------------------*/762/*--------------------------------- Gradient ---------------------------------*/763/*----------------------------------------------------------------------------*/764765/*----------------------------------------------------------------------------*/766/** Computes the direction of the level line of 'in' at each point.767768 The result is:769 - an image_double with the angle at each pixel, or NOTDEF if not defined.770 - the image_double 'modgrad' (a pointer is passed as argument)771 with the gradient magnitude at each point.772 - a list of pixels 'list_p' roughly ordered by decreasing773 gradient magnitude. (The order is made by classifying points774 into bins by gradient magnitude. The parameters 'n_bins' and775 'max_grad' specify the number of bins and the gradient modulus776 at the highest bin. The pixels in the list would be in777 decreasing gradient magnitude, up to a precision of the size of778 the bins.)779 - a pointer 'mem_p' to the memory used by 'list_p' to be able to780 free the memory when it is not used anymore.781 */782static image_double ll_angle(const image_double in, const double threshold,783 struct coorlist ** list_p, void ** mem_p,784 image_double * modgrad, const unsigned int n_bins )785{786 image_double g;787 unsigned int n,p,x,y,adr,i;788 double com1,com2,gx,gy,norm,norm2;789 /* the rest of the variables are used for pseudo-ordering790 the gradient magnitude values */791 int list_count = 0;792 struct coorlist * list;793 struct coorlist ** range_l_s; /* array of pointers to start of bin list */794 struct coorlist ** range_l_e; /* array of pointers to end of bin list */795 struct coorlist * start;796 struct coorlist * end;797 double max_grad = 0.0;798799 /* check parameters */800 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )801 error("ll_angle: invalid image.");802 if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");803 if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");804 if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");805 if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");806 if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");807808 /* image size shortcuts */809 n = in->ysize;810 p = in->xsize;811812 /* allocate output image */813 g = new_image_double(in->xsize,in->ysize);814815 /* get memory for the image of gradient modulus */816 *modgrad = new_image_double(in->xsize,in->ysize);817818 /* get memory for "ordered" list of pixels */819 list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );820 *mem_p = (void *) list;821 range_l_s = (struct coorlist **) calloc( (size_t) n_bins,822 sizeof(struct coorlist *) );823 range_l_e = (struct coorlist **) calloc( (size_t) n_bins,824 sizeof(struct coorlist *) );825 if( list == NULL || range_l_s == NULL || range_l_e == NULL )826 {827 free(list);828 free(range_l_s);829 free(range_l_e);830 error("not enough memory.");831 }832 for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;833834 /* 'undefined' on the down and right boundaries */835 for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;836 for(y=0;y<n;y++) g->data[p*y+p-1] = NOTDEF;837838 /* compute gradient on the remaining pixels */839 for(x=0;x<p-1;x++)840 for(y=0;y<n-1;y++)841 {842 adr = y*p+x;843844 /*845 Norm 2 computation using 2x2 pixel window:846 A B847 C D848 and849 com1 = D-A, com2 = B-C.850 Then851 gx = B+D - (A+C) horizontal difference852 gy = C+D - (A+B) vertical difference853 com1 and com2 are just to avoid 2 additions.854 */855 com1 = in->data[adr+p+1] - in->data[adr];856 com2 = in->data[adr+1] - in->data[adr+p];857858 gx = com1+com2; /* gradient x component */859 gy = com1-com2; /* gradient y component */860 norm2 = gx*gx+gy*gy;861 norm = sqrt( norm2 / 4.0 ); /* gradient norm */862863 (*modgrad)->data[adr] = norm; /* store gradient norm */864865 if( norm <= threshold ) /* norm too small, gradient no defined */866 g->data[adr] = NOTDEF; /* gradient angle not defined */867 else868 {869 /* gradient angle computation */870 g->data[adr] = atan2(gx,-gy);871872 /* look for the maximum of the gradient */873 if( norm > max_grad ) max_grad = norm;874 }875 }876877 /* compute histogram of gradient values */878 for(x=0;x<p-1;x++)879 for(y=0;y<n-1;y++)880 {881 norm = (*modgrad)->data[y*p+x];882883 /* store the point in the right bin according to its norm */884 i = (unsigned int) (norm * (double) n_bins / max_grad);885 if( i >= n_bins ) i = n_bins-1;886 if( range_l_e[i] == NULL )887 range_l_s[i] = range_l_e[i] = list+list_count++;888 else889 {890 range_l_e[i]->next = list+list_count;891 range_l_e[i] = list+list_count++;892 }893 range_l_e[i]->x = (int) x;894 range_l_e[i]->y = (int) y;895 range_l_e[i]->next = NULL;896 }897898 /* Make the list of pixels (almost) ordered by norm value.899 It starts by the larger bin, so the list starts by the900 pixels with the highest gradient value. Pixels would be ordered901 by norm value, up to a precision given by max_grad/n_bins.902 */903 for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);904 start = range_l_s[i];905 end = range_l_e[i];906 if( start != NULL )907 while(i>0)908 {909 --i;910 if( range_l_s[i] != NULL )911 {912 end->next = range_l_s[i];913 end = range_l_e[i];914 }915 }916 *list_p = start;917918 /* free memory */919 free( (void *) range_l_s );920 free( (void *) range_l_e );921922 return g;923}924925/*----------------------------------------------------------------------------*/926/** Is point (x,y) aligned to angle theta, up to precision 'prec'?927 */928static int isaligned(const int x, const int y, const image_double angles, double theta,929 const double prec )930{931 double a;932933 /* check parameters */934 if( angles == NULL || angles->data == NULL )935 error("isaligned: invalid image 'angles'.");936 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )937 error("isaligned: (x,y) out of the image.");938 if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");939940 /* angle at pixel (x,y) */941 a = angles->data[ x + y * angles->xsize ];942943 /* pixels whose level-line angle is not defined944 are considered as NON-aligned */945 if( a == NOTDEF ) return FALSE; /* there is no need to call the function946 'double_equal' here because there is947 no risk of problems related to the948 comparison doubles, we are only949 interested in the exact NOTDEF value */950951 /* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */952 theta -= a;953 if( theta < 0.0 ) theta = -theta;954 if( theta > M_3_2_PI )955 {956 theta -= M_2__PI;957 if( theta < 0.0 ) theta = -theta;958 }959960 return theta <= prec;961}962963/*----------------------------------------------------------------------------*/964/** Absolute value angle difference.965 */966static double angle_diff(double a, const double b)967{968 a -= b;969 while( a <= -M_PI ) a += M_2__PI;970 while( a > M_PI ) a -= M_2__PI;971 if( a < 0.0 ) a = -a;972 return a;973}974975/*----------------------------------------------------------------------------*/976/** Signed angle difference.977 */978static double angle_diff_signed(double a, const double b)979{980 a -= b;981 while( a <= -M_PI ) a += M_2__PI;982 while( a > M_PI ) a -= M_2__PI;983 return a;984}985986987/*----------------------------------------------------------------------------*/988/*----------------------------- NFA computation ------------------------------*/989/*----------------------------------------------------------------------------*/990991/*----------------------------------------------------------------------------*/992/** Computes the natural logarithm of the absolute value of993 the gamma function of x using the Lanczos approximation.994 See http://www.rskey.org/gamma.htm995996 The formula used is997 @f[998 \Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) }999 (x+5.5)^{x+0.5} e^{-(x+5.5)}1000 @f]1001 so1002 @f[1003 \log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right)1004 + (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n)1005 @f]1006 and1007 q0 = 75122.6331530,1008 q1 = 80916.6278952,1009 q2 = 36308.2951477,1010 q3 = 8687.24529705,1011 q4 = 1168.92649479,1012 q5 = 83.8676043424,1013 q6 = 2.50662827511.1014 */1015static double log_gamma_lanczos(const double x)1016{1017 static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,1018 8687.24529705, 1168.92649479, 83.8676043424,1019 2.50662827511 };1020 double a = (x+0.5) * log(x+5.5) - (x+5.5);1021 double b = 0.0;1022 int n;10231024 for(n=0;n<7;n++)1025 {1026 a -= log( x + (double) n );1027 b += q[n] * pow( x, (double) n );1028 }1029 return a + log(b);1030}10311032/*----------------------------------------------------------------------------*/1033/** Computes the natural logarithm of the absolute value of1034 the gamma function of x using Windschitl method.1035 See http://www.rskey.org/gamma.htm10361037 The formula used is1038 @f[1039 \Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e}1040 \sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x1041 @f]1042 so1043 @f[1044 \log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x1045 + 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right).1046 @f]1047 This formula is a good approximation when x > 15.1048 */1049static double log_gamma_windschitl(const double x)1050{1051 return 0.918938533204673 + (x-0.5)*log(x) - x1052 + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );1053}10541055/*----------------------------------------------------------------------------*/1056/** Computes the natural logarithm of the absolute value of1057 the gamma function of x. When x>15 use log_gamma_windschitl(),1058 otherwise use log_gamma_lanczos().1059 */106010611062/*----------------------------------------------------------------------------*/1063/** Size of the table to store already computed inverse values.1064 */106510661067// clang-format on10681069static double *inv = NULL; /* table to keep computed inverse values */10701071__attribute__((constructor)) static void invConstructor()1072{1073 if(inv) return;1074 inv = malloc(sizeof(double) * TABSIZE);1075}10761077__attribute__((destructor)) static void invDestructor()1078{1079 free(inv);1080 inv = NULL;1081}10821083// clang-format off10841085/*----------------------------------------------------------------------------*/1086/** Computes -log10(NFA).10871088 NFA stands for Number of False Alarms:1089 @f[1090 \mathrm{NFA} = NT \cdot B(n,k,p)1091 @f]10921093 - NT - number of tests1094 - B(n,k,p) - tail of binomial distribution with parameters n,k and p:1095 @f[1096 B(n,k,p) = \sum_{j=k}^n1097 \left(\begin{array}{c}n\\j\end{array}\right)1098 p^{j} (1-p)^{n-j}1099 @f]11001101 The value -log10(NFA) is equivalent but more intuitive than NFA:1102 - -1 corresponds to 10 mean false alarms1103 - 0 corresponds to 1 mean false alarm1104 - 1 corresponds to 0.1 mean false alarms1105 - 2 corresponds to 0.01 mean false alarms1106 - ...11071108 Used this way, the bigger the value, better the detection,1109 and a logarithmic scale is used.11101111 @param n,k,p binomial parameters.1112 @param logNT logarithm of Number of Tests11131114 The computation is based in the gamma function by the following1115 relation:1116 @f[1117 \left(\begin{array}{c}n\\k\end{array}\right)1118 = \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }.1119 @f]1120 We use efficient algorithms to compute the logarithm of1121 the gamma function.11221123 To make the computation faster, not all the sum is computed, part1124 of the terms are neglected based on a bound to the error obtained1125 (an error of 10% in the result is accepted).1126 */1127static double nfa(const int n, const int k, const double p, const double logNT)1128{1129 const double tolerance = 0.1; /* an error of 10% in the result is accepted */1130 double log1term,term,bin_term,mult_term,bin_tail,err,p_term;1131 int i;11321133 /* check parameters */1134 if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )1135 error("nfa: wrong n, k or p values.");11361137 /* trivial cases */1138 if( n==0 || k==0 ) return -logNT;1139 if( n==k ) return -logNT - (double) n * log10(p);11401141 /* probability term */1142 p_term = p / (1.0-p);11431144 /* compute the first term of the series */1145 /*1146 binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}1147 where bincoef(n,i) are the binomial coefficients.1148 But1149 bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).1150 We use this to compute the first term. Actually the log of it.1151 */1152 log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )1153 - log_gamma( (double) (n-k) + 1.0 )1154 + (double) k * log(p) + (double) (n-k) * log(1.0-p);1155 term = exp(log1term);11561157 /* in some cases no more computations are needed */1158 if( double_equal(term,0.0) ) /* the first term is almost zero */1159 {1160 if( (double) k > (double) n * p ) /* at begin or end of the tail? */1161 return -log1term / M_LN10 - logNT; /* end: use just the first term */1162 else1163 return -logNT; /* begin: the tail is roughly 1 */1164 }11651166 /* compute more terms if needed */1167 bin_tail = term;1168 for(i=k+1;i<=n;i++)1169 {1170 /*1171 As1172 term_i = bincoef(n,i) * p^i * (1-p)^(n-i)1173 and1174 bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,1175 then,1176 term_i / term_i-1 = (n-i+1)/i * p/(1-p)1177 and1178 term_i = term_i-1 * (n-i+1)/i * p/(1-p).1179 1/i is stored in a table as they are computed,1180 because divisions are expensive.1181 p/(1-p) is computed only once and stored in 'p_term'.1182 */1183 bin_term = (double) (n-i+1) * ( i<TABSIZE ?1184 ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :1185 1.0 / (double) i );11861187 mult_term = bin_term * p_term;1188 term *= mult_term;1189 bin_tail += term;1190 if(bin_term<1.0)1191 {1192 /* When bin_term<1 then mult_term_j<mult_term_i for j>i.1193 Then, the error on the binomial tail when truncated at1194 the i term can be bounded by a geometric series of form1195 term_i * sum mult_term_i^j. */1196 err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /1197 (1.0-mult_term) - 1.0 );11981199 /* One wants an error at most of tolerance*final_result, or:1200 tolerance * abs(-log10(bin_tail)-logNT).1201 Now, the error that can be accepted on bin_tail is1202 given by tolerance*final_result divided by the derivative1203 of -log10(x) when x=bin_tail. that is:1204 tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)1205 Finally, we truncate the tail if the error is less than:1206 tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */1207 if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;1208 }1209 }1210 return -log10(bin_tail) - logNT;1211}121212131214/*----------------------------------------------------------------------------*/1215/*--------------------------- Rectangle structure ----------------------------*/1216/*----------------------------------------------------------------------------*/12171218/*----------------------------------------------------------------------------*/1219/** Rectangle structure: line segment with width.1220 */1221struct rect1222{1223 double x1,y1,x2,y2; /* first and second point of the line segment */1224 double width; /* rectangle width */1225 double x,y; /* center of the rectangle */1226 double theta; /* angle */1227 double dx,dy; /* (dx,dy) is vector oriented as the line segment */1228 double prec; /* tolerance angle */1229 double p; /* probability of a point with angle within 'prec' */1230};12311232/*----------------------------------------------------------------------------*/1233/** Copy one rectangle structure to another.1234 */1235static void rect_copy(const struct rect * in, struct rect * out)1236{1237 /* check parameters */1238 if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");12391240 /* copy values */1241 out->x1 = in->x1;1242 out->y1 = in->y1;1243 out->x2 = in->x2;1244 out->y2 = in->y2;1245 out->width = in->width;1246 out->x = in->x;1247 out->y = in->y;1248 out->theta = in->theta;1249 out->dx = in->dx;1250 out->dy = in->dy;1251 out->prec = in->prec;1252 out->p = in->p;1253}12541255/*----------------------------------------------------------------------------*/1256/** Rectangle points iterator.12571258 The integer coordinates of pixels inside a rectangle are1259 iteratively explored. This structure keep track of the process and1260 functions ri_ini(), ri_inc(), ri_end(), and ri_del() are used in1261 the process. An example of how to use the iterator is as follows:1262 \code12631264 struct rect * rec = XXX; // some rectangle1265 rect_iter * i;1266 for( i=ri_ini(rec); !ri_end(i); ri_inc(i) )1267 {1268 // your code, using 'i->x' and 'i->y' as coordinates1269 }1270 ri_del(i); // delete iterator12711272 \endcode1273 The pixels are explored 'column' by 'column', where we call1274 'column' a set of pixels with the same x value that are inside the1275 rectangle. The following is an schematic representation of a1276 rectangle, the 'column' being explored is marked by colons, and1277 the current pixel being explored is 'x,y'.1278 \verbatim12791280 vx[1],vy[1]1281 * *1282 * *1283 * *1284 * ye1285 * : *1286 vx[0],vy[0] : *1287 * : *1288 * x,y *1289 * : *1290 * : vx[2],vy[2]1291 * : *1292 y ys *1293 ^ * *1294 | * *1295 | * *1296 +---> x vx[3],vy[3]12971298 \endverbatim1299 The first 'column' to be explored is the one with the smaller x1300 value. Each 'column' is explored starting from the pixel of the1301 'column' (inside the rectangle) with the smallest y value.13021303 The four corners of the rectangle are stored in order that rotates1304 around the corners at the arrays 'vx[]' and 'vy[]'. The first1305 point is always the one with smaller x value.13061307 'x' and 'y' are the coordinates of the pixel being explored. 'ys'1308 and 'ye' are the start and end values of the current column being1309 explored. So, 'ys' < 'ye'.1310 */1311typedef struct1312{1313 double vx[4]; /* rectangle's corner X coordinates in circular order */1314 double vy[4]; /* rectangle's corner Y coordinates in circular order */1315 double ys,ye; /* start and end Y values of current 'column' */1316 int x,y; /* coordinates of currently explored pixel */1317} rect_iter;13181319/*----------------------------------------------------------------------------*/1320/** Interpolate y value corresponding to 'x' value given, in1321 the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller1322 of 'y1' and 'y2'.13231324 The following restrictions are required:1325 - x1 <= x21326 - x1 <= x1327 - x <= x21328 */1329static double inter_low(const double x, const double x1, const double y1, const double x2, const double y2)1330{1331 /* check parameters */1332 if( x1 > x2 || x < x1 || x > x2 )1333 error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");13341335 /* interpolation */1336 if( double_equal(x1,x2) && y1<y2 ) return y1;1337 if( double_equal(x1,x2) && y1>y2 ) return y2;1338 return y1 + (x-x1) * (y2-y1) / (x2-x1);1339}13401341/*----------------------------------------------------------------------------*/1342/** Interpolate y value corresponding to 'x' value given, in1343 the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger1344 of 'y1' and 'y2'.13451346 The following restrictions are required:1347 - x1 <= x21348 - x1 <= x1349 - x <= x21350 */1351static double inter_hi(const double x, const double x1, const double y1, const double x2, const double y2)1352{1353 /* check parameters */1354 if( x1 > x2 || x < x1 || x > x2 )1355 error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");13561357 /* interpolation */1358 if( double_equal(x1,x2) && y1<y2 ) return y2;1359 if( double_equal(x1,x2) && y1>y2 ) return y1;1360 return y1 + (x-x1) * (y2-y1) / (x2-x1);1361}13621363/*----------------------------------------------------------------------------*/1364/** Free memory used by a rectangle iterator.1365 */1366static void ri_del(rect_iter * iter)1367{1368 if( iter == NULL ) error("ri_del: NULL iterator.");1369 free( (void *) iter );1370}13711372/*----------------------------------------------------------------------------*/1373/** Check if the iterator finished the full iteration.13741375 See details in \ref rect_iter1376 */1377static int ri_end(const rect_iter * i)1378{1379 /* check input */1380 if( i == NULL ) error("ri_end: NULL iterator.");13811382 /* if the current x value is larger than the largest1383 x value in the rectangle (vx[2]), we know the full1384 exploration of the rectangle is finished. */1385 return (double)(i->x) > i->vx[2];1386}13871388/*----------------------------------------------------------------------------*/1389/** Increment a rectangle iterator.13901391 See details in \ref rect_iter1392 */1393static void ri_inc(rect_iter * i)1394{1395 /* check input */1396 if( i == NULL ) error("ri_inc: NULL iterator.");13971398 /* if not at end of exploration,1399 increase y value for next pixel in the 'column' */1400 if( !ri_end(i) ) i->y++;14011402 /* if the end of the current 'column' is reached,1403 and it is not the end of exploration,1404 advance to the next 'column' */1405 while( (double) (i->y) > i->ye && !ri_end(i) )1406 {1407 /* increase x, next 'column' */1408 i->x++;14091410 /* if end of exploration, return */1411 if( ri_end(i) ) return;14121413 /* update lower y limit (start) for the new 'column'.14141415 We need to interpolate the y value that corresponds to the1416 lower side of the rectangle. The first thing is to decide if1417 the corresponding side is14181419 vx[0],vy[0] to vx[3],vy[3] or1420 vx[3],vy[3] to vx[2],vy[2]14211422 Then, the side is interpolated for the x value of the1423 'column'. But, if the side is vertical (as it could happen if1424 the rectangle is vertical and we are dealing with the first1425 or last 'columns') then we pick the lower value of the side1426 by using 'inter_low'.1427 */1428 if( (double) i->x < i->vx[3] )1429 i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);1430 else1431 i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);14321433 /* update upper y limit (end) for the new 'column'.14341435 We need to interpolate the y value that corresponds to the1436 upper side of the rectangle. The first thing is to decide if1437 the corresponding side is14381439 vx[0],vy[0] to vx[1],vy[1] or1440 vx[1],vy[1] to vx[2],vy[2]14411442 Then, the side is interpolated for the x value of the1443 'column'. But, if the side is vertical (as it could happen if1444 the rectangle is vertical and we are dealing with the first1445 or last 'columns') then we pick the lower value of the side1446 by using 'inter_low'.1447 */1448 if( (double)i->x < i->vx[1] )1449 i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);1450 else1451 i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);14521453 /* new y */1454 i->y = (int) ceil(i->ys);1455 }1456}14571458/*----------------------------------------------------------------------------*/1459/** Create and initialize a rectangle iterator.14601461 See details in \ref rect_iter1462 */1463static rect_iter * ri_ini(const struct rect * r)1464{1465 double vx[4],vy[4];1466 int n,offset;1467 rect_iter * i;14681469 /* check parameters */1470 if( r == NULL ) error("ri_ini: invalid rectangle.");14711472 /* get memory */1473 i = (rect_iter *) malloc(sizeof(rect_iter));1474 if( i == NULL ) error("ri_ini: Not enough memory.");14751476 /* build list of rectangle corners ordered1477 in a circular way around the rectangle */1478 vx[0] = r->x1 - r->dy * r->width / 2.0;1479 vy[0] = r->y1 + r->dx * r->width / 2.0;1480 vx[1] = r->x2 - r->dy * r->width / 2.0;1481 vy[1] = r->y2 + r->dx * r->width / 2.0;1482 vx[2] = r->x2 + r->dy * r->width / 2.0;1483 vy[2] = r->y2 - r->dx * r->width / 2.0;1484 vx[3] = r->x1 + r->dy * r->width / 2.0;1485 vy[3] = r->y1 - r->dx * r->width / 2.0;14861487 /* compute rotation of index of corners needed so that the first1488 point has the smaller x.14891490 if one side is vertical, thus two corners have the same smaller x1491 value, the one with the largest y value is selected as the first.1492 */1493 if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;1494 else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;1495 else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;1496 else offset = 3;14971498 /* apply rotation of index. */1499 for(n=0; n<4; n++)1500 {1501 i->vx[n] = vx[(offset+n)%4];1502 i->vy[n] = vy[(offset+n)%4];1503 }15041505 /* Set an initial condition.15061507 The values are set to values that will cause 'ri_inc' (that will1508 be called immediately) to initialize correctly the first 'column'1509 and compute the limits 'ys' and 'ye'.15101511 'y' is set to the integer value of vy[0], the starting corner.15121513 'ys' and 'ye' are set to very small values, so 'ri_inc' will1514 notice that it needs to start a new 'column'.15151516 The smallest integer coordinate inside of the rectangle is1517 'ceil(vx[0])'. The current 'x' value is set to that value minus1518 one, so 'ri_inc' (that will increase x by one) will advance to1519 the first 'column'.1520 */1521 i->x = (int) ceil(i->vx[0]) - 1;1522 i->y = (int) ceil(i->vy[0]);1523 i->ys = i->ye = -DBL_MAX;15241525 /* advance to the first pixel */1526 ri_inc(i);15271528 return i;1529}15301531/*----------------------------------------------------------------------------*/1532/** Compute a rectangle's NFA value.1533 */1534static double rect_nfa(const struct rect * rec, const image_double angles, const double logNT)1535{1536 rect_iter * i;1537 int pts = 0;1538 int alg = 0;15391540 /* check parameters */1541 if( rec == NULL ) error("rect_nfa: invalid rectangle.");1542 if( angles == NULL ) error("rect_nfa: invalid 'angles'.");15431544 /* compute the total number of pixels and of aligned points in 'rec' */1545 for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */1546 if( i->x >= 0 && i->y >= 0 &&1547 i->x < (int) angles->xsize && i->y < (int) angles->ysize )1548 {1549 ++pts; /* total number of pixels counter */1550 if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )1551 ++alg; /* aligned points counter */1552 }1553 ri_del(i); /* delete iterator */15541555 return nfa(pts,alg,rec->p,logNT); /* compute NFA value */1556}155715581559/*----------------------------------------------------------------------------*/1560/*---------------------------------- Regions ---------------------------------*/1561/*----------------------------------------------------------------------------*/15621563/*----------------------------------------------------------------------------*/1564/** Compute region's angle as the principal inertia axis of the region.15651566 The following is the region inertia matrix A:1567 @f[15681569 A = \left(\begin{array}{cc}1570 Ixx & Ixy \\1571 Ixy & Iyy \\1572 \end{array}\right)15731574 @f]1575 where15761577 Ixx = sum_i G(i).(y_i - cx)^215781579 Iyy = sum_i G(i).(x_i - cy)^215801581 Ixy = - sum_i G(i).(x_i - cx).(y_i - cy)15821583 and1584 - G(i) is the gradient norm at pixel i, used as pixel's weight.1585 - x_i and y_i are the coordinates of pixel i.1586 - cx and cy are the coordinates of the center of th region.15871588 lambda1 and lambda2 are the eigenvalues of matrix A,1589 with lambda1 >= lambda2. They are found by solving the1590 characteristic polynomial:15911592 det( lambda I - A) = 015931594 that gives:15951596 lambda1 = ( Ixx + Iyy + sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 215971598 lambda2 = ( Ixx + Iyy - sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 215991600 To get the line segment direction we want to get the angle the1601 eigenvector associated to the smallest eigenvalue. We have1602 to solve for a,b in:16031604 a.Ixx + b.Ixy = a.lambda216051606 a.Ixy + b.Iyy = b.lambda216071608 We want the angle theta = atan(b/a). It can be computed with1609 any of the two equations:16101611 theta = atan( (lambda2-Ixx) / Ixy )16121613 or16141615 theta = atan( Ixy / (lambda2-Iyy) )16161617 When |Ixx| > |Iyy| we use the first, otherwise the second (just to1618 get better numeric precision).1619 */1620static double get_theta(const struct point * reg, const int reg_size, const double x, const double y,1621 const image_double modgrad, const double reg_angle, const double prec )1622{1623 double lambda,theta,weight;1624 double Ixx = 0.0;1625 double Iyy = 0.0;1626 double Ixy = 0.0;1627 int i;16281629 /* check parameters */1630 if( reg == NULL ) error("get_theta: invalid region.");1631 if( reg_size <= 1 ) error("get_theta: region size <= 1.");1632 if( modgrad == NULL || modgrad->data == NULL )1633 error("get_theta: invalid 'modgrad'.");1634 if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");16351636 /* compute inertia matrix */1637 for(i=0; i<reg_size; i++)1638 {1639 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];1640 Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;1641 Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;1642 Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;1643 }1644 if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )1645 error("get_theta: null inertia matrix.");16461647 /* compute smallest eigenvalue */1648 lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );16491650 /* compute angle */1651 theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);16521653 /* The previous procedure doesn't cares about orientation,1654 so it could be wrong by 180 degrees. Here is corrected if necessary. */1655 if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;16561657 return theta;1658}16591660/*----------------------------------------------------------------------------*/1661/** Computes a rectangle that covers a region of points.1662 */1663static void region2rect(const struct point * reg, const int reg_size,1664 const image_double modgrad, const double reg_angle,1665 const double prec, const double p, struct rect * rec )1666{1667 double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;1668 int i;16691670 /* check parameters */1671 if( reg == NULL ) error("region2rect: invalid region.");1672 if( reg_size <= 1 ) error("region2rect: region size <= 1.");1673 if( modgrad == NULL || modgrad->data == NULL )1674 error("region2rect: invalid image 'modgrad'.");1675 if( rec == NULL ) error("region2rect: invalid 'rec'.");16761677 /* center of the region:16781679 It is computed as the weighted sum of the coordinates1680 of all the pixels in the region. The norm of the gradient1681 is used as the weight of a pixel. The sum is as follows:1682 cx = \sum_i G(i).x_i1683 cy = \sum_i G(i).y_i1684 where G(i) is the norm of the gradient of pixel i1685 and x_i,y_i are its coordinates.1686 */1687 x = y = sum = 0.0;1688 for(i=0; i<reg_size; i++)1689 {1690 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];1691 x += (double) reg[i].x * weight;1692 y += (double) reg[i].y * weight;1693 sum += weight;1694 }1695 if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");1696 x /= sum;1697 y /= sum;16981699 /* theta */1700 theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);17011702 /* length and width:17031704 'l' and 'w' are computed as the distance from the center of the1705 region to pixel i, projected along the rectangle axis (dx,dy) and1706 to the orthogonal axis (-dy,dx), respectively.17071708 The length of the rectangle goes from l_min to l_max, where l_min1709 and l_max are the minimum and maximum values of l in the region.1710 Analogously, the width is selected from w_min to w_max, where1711 w_min and w_max are the minimum and maximum of w for the pixels1712 in the region.1713 */1714 dx = cos(theta);1715 dy = sin(theta);1716 l_min = l_max = w_min = w_max = 0.0;1717 for(i=0; i<reg_size; i++)1718 {1719 l = ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;1720 w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;17211722 if( l > l_max ) l_max = l;1723 if( l < l_min ) l_min = l;1724 if( w > w_max ) w_max = w;1725 if( w < w_min ) w_min = w;1726 }17271728 /* store values */1729 rec->x1 = x + l_min * dx;1730 rec->y1 = y + l_min * dy;1731 rec->x2 = x + l_max * dx;1732 rec->y2 = y + l_max * dy;1733 rec->width = w_max - w_min;1734 rec->x = x;1735 rec->y = y;1736 rec->theta = theta;1737 rec->dx = dx;1738 rec->dy = dy;1739 rec->prec = prec;1740 rec->p = p;17411742 /* we impose a minimal width of one pixel17431744 A sharp horizontal or vertical step would produce a perfectly1745 horizontal or vertical region. The width computed would be1746 zero. But that corresponds to a one pixels width transition in1747 the image.1748 */1749 if( rec->width < 1.0 ) rec->width = 1.0;1750}17511752/*----------------------------------------------------------------------------*/1753/** Build a region of pixels that share the same angle, up to a1754 tolerance 'prec', starting at point (x,y).1755 */1756static void region_grow(const int x, const int y, const image_double angles, struct point * reg,1757 int * reg_size, double * reg_angle, const image_char used,1758 const double prec )1759{1760 double sumdx,sumdy;1761 int xx,yy,i;17621763 /* check parameters */1764 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )1765 error("region_grow: (x,y) out of the image.");1766 if( angles == NULL || angles->data == NULL )1767 error("region_grow: invalid image 'angles'.");1768 if( reg == NULL ) error("region_grow: invalid 'reg'.");1769 if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");1770 if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");1771 if( used == NULL || used->data == NULL )1772 error("region_grow: invalid image 'used'.");17731774 /* first point of the region */1775 *reg_size = 1;1776 reg[0].x = x;1777 reg[0].y = y;1778 *reg_angle = angles->data[x+y*angles->xsize]; /* region's angle */1779 sumdx = cos(*reg_angle);1780 sumdy = sin(*reg_angle);1781 used->data[x+y*used->xsize] = USED;17821783 /* try neighbors as new region points */1784 for(i=0; i<*reg_size; i++)1785 for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)1786 for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)1787 if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&1788 used->data[xx+yy*used->xsize] != USED &&1789 isaligned(xx,yy,angles,*reg_angle,prec) )1790 {1791 /* add point */1792 used->data[xx+yy*used->xsize] = USED;1793 reg[*reg_size].x = xx;1794 reg[*reg_size].y = yy;1795 ++(*reg_size);17961797 /* update region's angle */1798 sumdx += cos( angles->data[xx+yy*angles->xsize] );1799 sumdy += sin( angles->data[xx+yy*angles->xsize] );1800 *reg_angle = atan2(sumdy,sumdx);1801 }1802}18031804/*----------------------------------------------------------------------------*/1805/** Try some rectangles variations to improve NFA value. Only if the1806 rectangle is not meaningful (i.e., log_nfa <= log_eps).1807 */1808static double rect_improve( struct rect * rec, const image_double angles,1809 const double logNT, const double log_eps )1810{1811 struct rect r;1812 double log_nfa,log_nfa_new;1813 const double delta = 0.5;1814 const double delta_2 = delta / 2.0;1815 int n;18161817 log_nfa = rect_nfa(rec,angles,logNT);18181819 if( log_nfa > log_eps ) return log_nfa;18201821 /* try finer precisions */1822 rect_copy(rec,&r);1823 for(n=0; n<5; n++)1824 {1825 r.p /= 2.0;1826 r.prec = r.p * M_PI;1827 log_nfa_new = rect_nfa(&r,angles,logNT);1828 if( log_nfa_new > log_nfa )1829 {1830 log_nfa = log_nfa_new;1831 rect_copy(&r,rec);1832 }1833 }18341835 if( log_nfa > log_eps ) return log_nfa;18361837 /* try to reduce width */1838 rect_copy(rec,&r);1839 for(n=0; n<5; n++)1840 {1841 if( (r.width - delta) >= 0.5 )1842 {1843 r.width -= delta;1844 log_nfa_new = rect_nfa(&r,angles,logNT);1845 if( log_nfa_new > log_nfa )1846 {1847 rect_copy(&r,rec);1848 log_nfa = log_nfa_new;1849 }1850 }1851 }18521853 if( log_nfa > log_eps ) return log_nfa;18541855 /* try to reduce one side of the rectangle */1856 rect_copy(rec,&r);1857 for(n=0; n<5; n++)1858 {1859 if( (r.width - delta) >= 0.5 )1860 {1861 r.x1 += -r.dy * delta_2;1862 r.y1 += r.dx * delta_2;1863 r.x2 += -r.dy * delta_2;1864 r.y2 += r.dx * delta_2;1865 r.width -= delta;1866 log_nfa_new = rect_nfa(&r,angles,logNT);1867 if( log_nfa_new > log_nfa )1868 {1869 rect_copy(&r,rec);1870 log_nfa = log_nfa_new;1871 }1872 }1873 }18741875 if( log_nfa > log_eps ) return log_nfa;18761877 /* try to reduce the other side of the rectangle */1878 rect_copy(rec,&r);1879 for(n=0; n<5; n++)1880 {1881 if( (r.width - delta) >= 0.5 )1882 {1883 r.x1 -= -r.dy * delta_2;1884 r.y1 -= r.dx * delta_2;1885 r.x2 -= -r.dy * delta_2;1886 r.y2 -= r.dx * delta_2;1887 r.width -= delta;1888 log_nfa_new = rect_nfa(&r,angles,logNT);1889 if( log_nfa_new > log_nfa )1890 {1891 rect_copy(&r,rec);1892 log_nfa = log_nfa_new;1893 }1894 }1895 }18961897 if( log_nfa > log_eps ) return log_nfa;18981899 /* try even finer precisions */1900 rect_copy(rec,&r);1901 for(n=0; n<5; n++)1902 {1903 r.p /= 2.0;1904 r.prec = r.p * M_PI;1905 log_nfa_new = rect_nfa(&r,angles,logNT);1906 if( log_nfa_new > log_nfa )1907 {1908 log_nfa = log_nfa_new;1909 rect_copy(&r,rec);1910 }1911 }19121913 return log_nfa;1914}19151916/*----------------------------------------------------------------------------*/1917/** Reduce the region size, by elimination the points far from the1918 starting point, until that leads to rectangle with the right1919 density of region points or to discard the region if too small.1920 */1921static int reduce_region_radius( struct point * reg, int * reg_size,1922 const image_double modgrad, const double reg_angle,1923 const double prec, const double p, struct rect * rec,1924 const image_char used, const image_double angles,1925 const double density_th )1926{1927 double density,radius1,radius2,rad,xc,yc;1928 int i;19291930 /* check parameters */1931 if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");1932 if( reg_size == NULL )1933 error("reduce_region_radius: invalid pointer 'reg_size'.");1934 if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");1935 if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");1936 if( used == NULL || used->data == NULL )1937 error("reduce_region_radius: invalid image 'used'.");1938 if( angles == NULL || angles->data == NULL )1939 error("reduce_region_radius: invalid image 'angles'.");19401941 /* compute region points density */1942 density = (double) *reg_size /1943 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );19441945 /* if the density criterion is satisfied there is nothing to do */1946 if( density >= density_th ) return TRUE;19471948 /* compute region's radius */1949 xc = (double) reg[0].x;1950 yc = (double) reg[0].y;1951 radius1 = dist( xc, yc, rec->x1, rec->y1 );1952 radius2 = dist( xc, yc, rec->x2, rec->y2 );1953 rad = radius1 > radius2 ? radius1 : radius2;19541955 /* while the density criterion is not satisfied, remove farther pixels */1956 while( density < density_th )1957 {1958 rad *= 0.75; /* reduce region's radius to 75% of its value */19591960 /* remove points from the region and update 'used' map */1961 for(i=0; i<*reg_size; i++)1962 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )1963 {1964 /* point not kept, mark it as NOTUSED */1965 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;1966 /* remove point from the region */1967 reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */1968 reg[i].y = reg[*reg_size-1].y;1969 --(*reg_size);1970 --i; /* to avoid skipping one point */1971 }19721973 /* reject if the region is too small.1974 2 is the minimal region size for 'region2rect' to work. */1975 if( *reg_size < 2 ) return FALSE;19761977 /* re-compute rectangle */1978 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);19791980 /* re-compute region points density */1981 density = (double) *reg_size /1982 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );1983 }19841985 /* if this point is reached, the density criterion is satisfied */1986 return TRUE;1987}19881989/*----------------------------------------------------------------------------*/1990/** Refine a rectangle.19911992 For that, an estimation of the angle tolerance is performed by the1993 standard deviation of the angle at points near the region's1994 starting point. Then, a new region is grown starting from the same1995 point, but using the estimated angle tolerance. If this fails to1996 produce a rectangle with the right density of region points,1997 'reduce_region_radius' is called to try to satisfy this condition.1998 */1999static int refine( struct point * reg, int * reg_size, const image_double modgrad,2000 double reg_angle, const double prec, const double p, struct rect * rec,2001 const image_char used, const image_double angles, const double density_th )2002{2003 double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;2004 int i,n;20052006 /* check parameters */2007 if( reg == NULL ) error("refine: invalid pointer 'reg'.");2008 if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");2009 if( prec < 0.0 ) error("refine: 'prec' must be positive.");2010 if( rec == NULL ) error("refine: invalid pointer 'rec'.");2011 if( used == NULL || used->data == NULL )2012 error("refine: invalid image 'used'.");2013 if( angles == NULL || angles->data == NULL )2014 error("refine: invalid image 'angles'.");20152016 /* compute region points density */2017 density = (double) *reg_size /2018 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );20192020 /* if the density criterion is satisfied there is nothing to do */2021 if( density >= density_th ) return TRUE;20222023 /*------ First try: reduce angle tolerance ------*/20242025 /* compute the new mean angle and tolerance */2026 xc = (double) reg[0].x;2027 yc = (double) reg[0].y;2028 ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];2029 sum = s_sum = 0.0;2030 n = 0;2031 for(i=0; i<*reg_size; i++)2032 {2033 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;2034 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width )2035 {2036 angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];2037 ang_d = angle_diff_signed(angle,ang_c);2038 sum += ang_d;2039 s_sum += ang_d * ang_d;2040 ++n;2041 }2042 }20432044 /* should not happen */2045 if(n == 0) return FALSE;20462047 mean_angle = sum / (double) n;2048 tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n2049 + mean_angle*mean_angle ); /* 2 * standard deviation */20502051 /* find a new region from the same starting point and new angle tolerance */2052 region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,®_angle,used,tau);20532054 /* if the region is too small, reject */2055 if( *reg_size < 2 ) return FALSE;20562057 /* re-compute rectangle */2058 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);20592060 /* re-compute region points density */2061 density = (double) *reg_size /2062 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );20632064 /*------ Second try: reduce region radius ------*/2065 if( density < density_th )2066 return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,2067 rec, used, angles, density_th );20682069 /* if this point is reached, the density criterion is satisfied */2070 return TRUE;2071}207220732074/*----------------------------------------------------------------------------*/2075/*-------------------------- Line Segment Detector ---------------------------*/2076/*----------------------------------------------------------------------------*/20772078/*----------------------------------------------------------------------------*/2079/** LSD full interface.2080 */2081static2082double * LineSegmentDetection( int * n_out,2083 double * img, const int X, const int Y,2084 const double scale, const double sigma_scale, const double quant,2085 const double ang_th, const double log_eps, const double density_th,2086 const int n_bins,2087 int ** reg_img, int * reg_x, int * reg_y )2088{2089 image_double image;2090 const ntuple_list out = new_ntuple_list(7);2091 double * return_value;2092 image_double scaled_image,angles,modgrad;2093 image_char used;2094 image_int region = NULL;2095 struct coorlist * list_p;2096 void * mem_p;2097 struct rect rec;2098 struct point * reg;2099 int reg_size,min_reg_size,i;2100 unsigned int xsize,ysize;2101 double rho,reg_angle,prec,p,log_nfa,logNT;2102 int ls_count = 0; /* line segments are numbered 1,2,3,... */210321042105 /* check parameters */2106 if( img == NULL || X <= 0 || Y <= 0 ) error("invalid image input.");2107 if( scale <= 0.0 ) error("'scale' value must be positive.");2108 if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");2109 if( quant < 0.0 ) error("'quant' value must be positive.");2110 if( ang_th <= 0.0 || ang_th >= 180.0 )2111 error("'ang_th' value must be in the range (0,180).");2112 if( density_th < 0.0 || density_th > 1.0 )2113 error("'density_th' value must be in the range [0,1].");2114 if( n_bins <= 0 ) error("'n_bins' value must be positive.");211521162117 /* angle tolerance */2118 prec = deg2rad(ang_th);2119 p = ang_th / 180.0;2120 rho = quant / sin(prec); /* gradient magnitude threshold */212121222123 /* load and scale image (if necessary) and compute angle at each pixel */2124 image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );2125 if( scale != 1.0 )2126 {2127 scaled_image = gaussian_sampler( image, scale, sigma_scale );2128 angles = ll_angle( scaled_image, rho, &list_p, &mem_p,2129 &modgrad, (unsigned int) n_bins );2130 free_image_double(scaled_image);2131 }2132 else2133 angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,2134 (unsigned int) n_bins );2135 xsize = angles->xsize;2136 ysize = angles->ysize;21372138 /* Number of Tests - NT21392140 The theoretical number of tests is Np.(XY)^(5/2)2141 where X and Y are number of columns and rows of the image.2142 Np corresponds to the number of angle precisions considered.2143 As the procedure 'rect_improve' tests 5 times to halve the2144 angle precision, and 5 more times after improving other factors,2145 11 different precision values are potentially tested. Thus,2146 the number of tests is2147 11 * (X*Y)^(5/2)2148 whose logarithm value is2149 log10(11) + 5/2 * (log10(X) + log10(Y)).2150 */2151 logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.02152 + log10(11.0);2153 min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region2154 that can give a meaningful event */215521562157 /* initialize some structures */2158 if( reg_img != NULL && reg_x != NULL && reg_y != NULL ) /* save region data */2159 region = new_image_int_ini(angles->xsize,angles->ysize,0);2160 used = new_image_char_ini(xsize,ysize,NOTUSED);2161 reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );2162 if( reg == NULL ) error("not enough memory!");216321642165 /* search for line segments */2166 for(; list_p != NULL; list_p = list_p->next )2167 if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&2168 angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )2169 /* there is no risk of double comparison problems here2170 because we are only interested in the exact NOTDEF value */2171 {2172 /* find the region of connected point and ~equal angle */2173 region_grow( list_p->x, list_p->y, angles, reg, ®_size,2174 ®_angle, used, prec );21752176 /* reject small regions */2177 if( reg_size < min_reg_size ) continue;21782179 /* construct rectangular approximation for the region */2180 region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);21812182 /* Check if the rectangle exceeds the minimal density of2183 region points. If not, try to improve the region.2184 The rectangle will be rejected if the final one does2185 not fulfill the minimal density condition.2186 This is an addition to the original LSD algorithm published in2187 "LSD: A Fast Line Segment Detector with a False Detection Control"2188 by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall.2189 The original algorithm is obtained with density_th = 0.0.2190 */2191 if( !refine( reg, ®_size, modgrad, reg_angle,2192 prec, p, &rec, used, angles, density_th ) ) continue;21932194 /* compute NFA value */2195 log_nfa = rect_improve(&rec,angles,logNT,log_eps);2196 if( log_nfa <= log_eps ) continue;21972198 /* A New Line Segment was found! */2199 ++ls_count; /* increase line segment counter */22002201 /*2202 The gradient was computed with a 2x2 mask, its value corresponds to2203 points with an offset of (0.5,0.5), that should be added to output.2204 The coordinates origin is at the center of pixel (0,0).2205 */2206 rec.x1 += 0.5; rec.y1 += 0.5;2207 rec.x2 += 0.5; rec.y2 += 0.5;22082209 /* scale the result values if a subsampling was performed */2210 if( scale != 1.0 )2211 {2212 rec.x1 /= scale; rec.y1 /= scale;2213 rec.x2 /= scale; rec.y2 /= scale;2214 rec.width /= scale;2215 }22162217 /* add line segment found to output */2218 add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,2219 rec.width, rec.p, log_nfa );22202221 /* add region number to 'region' image if needed */2222 if( region != NULL )2223 for(i=0; i<reg_size; i++)2224 region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;2225 }222622272228 /* free memory */2229 free( (void *) image ); /* only the double_image structure should be freed,2230 the data pointer was provided to this functions2231 and should not be destroyed. */2232 free_image_double(angles);2233 free_image_double(modgrad);2234 free_image_char(used);2235 free( (void *) reg );2236 free( (void *) mem_p );22372238 /* return the result */2239 if( reg_img != NULL && reg_x != NULL && reg_y != NULL )2240 {2241 if( region == NULL ) error("'region' should be a valid image.");2242 *reg_img = region->data;2243 if( region->xsize > (unsigned int) INT_MAX ||2244 region->ysize > (unsigned int) INT_MAX )2245 error("region image to big to fit in INT sizes.");2246 *reg_x = (int) (region->xsize);2247 *reg_y = (int) (region->ysize);22482249 /* free the 'region' structure.2250 we cannot use the function 'free_image_int' because we need to keep2251 the memory with the image data to be returned by this function. */2252 free( (void *) region );2253 }2254 if( out->size > (unsigned int) INT_MAX )2255 error("too many detections to fit in an INT.");2256 *n_out = (int) (out->size);22572258 return_value = out->values;2259 free( (void *) out ); /* only the 'ntuple_list' structure must be freed,2260 but the 'values' pointer must be keep to return2261 as a result. */22622263 return return_value;2264}22652266/*==================================================================================2267 * end of LSD code2268 *==================================================================================*/22692270// clang-format on22712272227322742275227622772278// clang-format off2279// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py2280// vim: shiftwidth=2 expandtab tabstop=2 cindent2281// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;2282// clang-format on2283