Main Page | Data Structures | File List | Data Fields | Globals

lsd.c

Go to the documentation of this file.
00001 /*----------------------------------------------------------------------------
00002 
00003   LSD - Line Segment Detector on digital images
00004 
00005   Copyright 2007-2010 rafael grompone von gioi (grompone@gmail.com)
00006 
00007   This program is free software: you can redistribute it and/or modify
00008   it under the terms of the GNU Affero General Public License as
00009   published by the Free Software Foundation, either version 3 of the
00010   License, or (at your option) any later version.
00011 
00012   This program is distributed in the hope that it will be useful,
00013   but WITHOUT ANY WARRANTY; without even the implied warranty of
00014   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00015   GNU Affero General Public License for more details.
00016 
00017   You should have received a copy of the GNU Affero General Public License
00018   along with this program. If not, see <http://www.gnu.org/licenses/>.
00019 
00020   ----------------------------------------------------------------------------*/
00021 
00022 /*----------------------------------------------------------------------------*/
00023 /** @file lsd.c
00024     LSD module code
00025     @author rafael grompone von gioi (grompone@gmail.com)
00026  */
00027 /*----------------------------------------------------------------------------*/
00028 
00029 /*----------------------------------------------------------------------------*/
00030 /** @mainpage LSD code documentation
00031 
00032     This is an implementation of the Line Segment Detector described
00033     in the paper:
00034 
00035       "LSD: A Fast Line Segment Detector with a False Detection Control"
00036       by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel,
00037       and Gregory Randall, IEEE Transactions on Pattern Analysis and
00038       Machine Intelligence, vol. 32, no. 4, pp. 722-732, April, 2010.
00039 
00040     and in more details in the CMLA Technical Report:
00041 
00042       "LSD: A Line Segment Detector, Technical Report",
00043       by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel,
00044       Gregory Randall, CMLA, ENS Cachan, 2010.
00045 
00046     The version implemented here includes some further improvements
00047     described on the LSD page at www.ipol.im. That same page includes
00048     more information, including this code and an online demo version:
00049 
00050       http://www.ipol.im/pub/algo/gjmr_line_segment_detector
00051 
00052     The module's main function is lsd().
00053 
00054     The source code is contained in two files: lsd.h and lsd.c.
00055 
00056     HISTORY:
00057     - version 1.5 - dic 2010: Changes in 'refine', -W option added,
00058                               and more comments added.
00059     - version 1.4 - jul 2010: lsd_scale interface added and doxygen doc.
00060     - version 1.3 - feb 2010: Multiple bug correction and improved code.
00061     - version 1.2 - dic 2009: First full Ansi C Language version.
00062     - version 1.1 - sep 2009: Systematic subsampling to scale 0.8 and
00063                               correction to partially handle"angle problem".
00064     - version 1.0 - jan 2009: First complete Megawave2 and Ansi C Language
00065                               version.
00066 
00067     @author rafael grompone von gioi (grompone@gmail.com)
00068  */
00069 /*----------------------------------------------------------------------------*/
00070 
00071 #include <stdio.h>
00072 #include <stdlib.h>
00073 #include <math.h>
00074 #include <limits.h>
00075 #include <float.h>
00076 #include "lsd.h"
00077 
00078 /** ln(10) */
00079 #ifndef M_LN10
00080 #define M_LN10 2.30258509299404568402
00081 #endif /* !M_LN10 */
00082 
00083 /** PI */
00084 #ifndef M_PI
00085 #define M_PI   3.14159265358979323846
00086 #endif /* !M_PI */
00087 
00088 #ifndef FALSE
00089 #define FALSE 0
00090 #endif /* !FALSE */
00091 
00092 #ifndef TRUE
00093 #define TRUE 1
00094 #endif /* !TRUE */
00095 
00096 /** Label for pixels with undefined gradient. */
00097 #define NOTDEF -1024.0
00098 
00099 /** 3/2 pi */
00100 #define M_3_2_PI 4.71238898038
00101 
00102 /** 2 pi */
00103 #define M_2__PI  6.28318530718
00104 
00105 /** Label for pixels not used in yet. */
00106 #define NOTUSED 0
00107 
00108 /** Label for pixels already used in detection. */
00109 #define USED    1
00110 
00111 /*----------------------------------------------------------------------------*/
00112 /** Chained list of coordinates.
00113  */
00114 struct coorlist
00115 {
00116   int x,y;
00117   struct coorlist * next;
00118 };
00119 
00120 /*----------------------------------------------------------------------------*/
00121 /** A point (or pixel).
00122  */
00123 struct point {int x,y;};
00124 
00125 
00126 /*----------------------------------------------------------------------------*/
00127 /*------------------------- Miscellaneous functions --------------------------*/
00128 /*----------------------------------------------------------------------------*/
00129 
00130 /*----------------------------------------------------------------------------*/
00131 /** Fatal error, print a message to standard-error output and exit.
00132  */
00133 static void error(char * msg)
00134 {
00135   fprintf(stderr,"LSD Error: %s\n",msg);
00136   exit(EXIT_FAILURE);
00137 }
00138 
00139 /*----------------------------------------------------------------------------*/
00140 /** Doubles relative error factor
00141  */
00142 #define RELATIVE_ERROR_FACTOR 100.0
00143 
00144 /*----------------------------------------------------------------------------*/
00145 /** Compare doubles by relative error.
00146 
00147     The resulting rounding error after floating point computations
00148     depend on the specific operations done. The same number computed by
00149     different algorithms could present different rounding errors. For a
00150     useful comparison, an estimation of the relative rounding error
00151     should be considered and compared to a factor times EPS. The factor
00152     should be related to the cumulated rounding error in the chain of
00153     computation. Here, as a simplification, a fixed factor is used.
00154  */
00155 static int double_equal(double a, double b)
00156 {
00157   double abs_diff,aa,bb,abs_max;
00158 
00159   /* trivial case */
00160   if( a == b ) return TRUE;
00161 
00162   abs_diff = fabs(a-b);
00163   aa = fabs(a);
00164   bb = fabs(b);
00165   abs_max = aa > bb ? aa : bb;
00166 
00167   /* DBL_MIN is the smallest normalized number, thus, the smallest
00168      number whose relative error is bounded by DBL_EPSILON. For
00169      smaller numbers, the same quantization steps as for DBL_MIN
00170      are used. Then, for smaller numbers, a meaningful "relative"
00171      error should be computed by dividing the difference by DBL_MIN. */
00172   if( abs_max < DBL_MIN ) abs_max = DBL_MIN;
00173 
00174   /* equal if relative error <= factor x eps */
00175   return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
00176 }
00177 
00178 /*----------------------------------------------------------------------------*/
00179 /** Computes Euclidean distance between point (x1,y1) and point (x2,y2).
00180  */
00181 static double dist(double x1, double y1, double x2, double y2)
00182 {
00183   return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00184 }
00185 
00186 
00187 /*----------------------------------------------------------------------------*/
00188 /*----------------------- 'list of n-tuple' data type ------------------------*/
00189 /*----------------------------------------------------------------------------*/
00190 
00191 /*----------------------------------------------------------------------------*/
00192 /** Free memory used in n-tuple 'in'.
00193  */
00194 void free_ntuple_list(ntuple_list in)
00195 {
00196   if( in == NULL || in->values == NULL )
00197     error("free_ntuple_list: invalid n-tuple input.");
00198   free( (void *) in->values );
00199   free( (void *) in );
00200 }
00201 
00202 /*----------------------------------------------------------------------------*/
00203 /** Create an n-tuple list and allocate memory for one element.
00204     @param dim the dimension (n) of the n-tuple.
00205  */
00206 ntuple_list new_ntuple_list(unsigned int dim)
00207 {
00208   ntuple_list n_tuple;
00209 
00210   /* check parameters */
00211   if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");
00212 
00213   /* get memory for list structure */
00214   n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );
00215   if( n_tuple == NULL ) error("not enough memory.");
00216 
00217   /* initialize list */
00218   n_tuple->size = 0;
00219   n_tuple->max_size = 1;
00220   n_tuple->dim = dim;
00221 
00222   /* get memory for tuples */
00223   n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) );
00224   if( n_tuple->values == NULL ) error("not enough memory.");
00225 
00226   return n_tuple;
00227 }
00228 
00229 /*----------------------------------------------------------------------------*/
00230 /** Enlarge the allocated memory of an n-tuple list.
00231  */
00232 static void enlarge_ntuple_list(ntuple_list n_tuple)
00233 {
00234   /* check parameters */
00235   if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )
00236     error("enlarge_ntuple_list: invalid n-tuple.");
00237 
00238   /* duplicate number of tuples */
00239   n_tuple->max_size *= 2;
00240 
00241   /* realloc memory */
00242   n_tuple->values = (double *) realloc( (void *) n_tuple->values,
00243                       n_tuple->dim * n_tuple->max_size * sizeof(double) );
00244   if( n_tuple->values == NULL ) error("not enough memory.");
00245 }
00246 
00247 /*----------------------------------------------------------------------------*/
00248 /** Add a 5-tuple to an n-tuple list.
00249  */
00250 static void add_5tuple( ntuple_list out, double v1, double v2,
00251                         double v3, double v4, double v5 )
00252 {
00253   /* check parameters */
00254   if( out == NULL ) error("add_5tuple: invalid n-tuple input.");
00255   if( out->dim != 5 ) error("add_5tuple: the n-tuple must be a 5-tuple.");
00256 
00257   /* if needed, alloc more tuples to 'out' */
00258   if( out->size == out->max_size ) enlarge_ntuple_list(out);
00259   if( out->values == NULL ) error("add_5tuple: invalid n-tuple input.");
00260 
00261   /* add new 5-tuple */
00262   out->values[ out->size * out->dim + 0 ] = v1;
00263   out->values[ out->size * out->dim + 1 ] = v2;
00264   out->values[ out->size * out->dim + 2 ] = v3;
00265   out->values[ out->size * out->dim + 3 ] = v4;
00266   out->values[ out->size * out->dim + 4 ] = v5;
00267 
00268   /* update number of tuples counter */
00269   out->size++;
00270 }
00271 
00272 
00273 /*----------------------------------------------------------------------------*/
00274 /*----------------------------- Image Data Types -----------------------------*/
00275 /*----------------------------------------------------------------------------*/
00276 
00277 /*----------------------------------------------------------------------------*/
00278 /** Free memory used in image_char 'i'.
00279  */
00280 void free_image_char(image_char i)
00281 {
00282   if( i == NULL || i->data == NULL )
00283     error("free_image_char: invalid input image.");
00284   free( (void *) i->data );
00285   free( (void *) i );
00286 }
00287 
00288 /*----------------------------------------------------------------------------*/
00289 /** Create a new image_char of size 'xsize' times 'ysize'.
00290  */
00291 image_char new_image_char(unsigned int xsize, unsigned int ysize)
00292 {
00293   image_char image;
00294 
00295   /* check parameters */
00296   if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");
00297 
00298   /* get memory */
00299   image = (image_char) malloc( sizeof(struct image_char_s) );
00300   if( image == NULL ) error("not enough memory.");
00301   image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),
00302                                           sizeof(unsigned char) );
00303   if( image->data == NULL ) error("not enough memory.");
00304 
00305   /* set image size */
00306   image->xsize = xsize;
00307   image->ysize = ysize;
00308 
00309   return image;
00310 }
00311 
00312 /*----------------------------------------------------------------------------*/
00313 /** Create a new image_char of size 'xsize' times 'ysize',
00314     initialized to the value 'fill_value'.
00315  */
00316 image_char new_image_char_ini( unsigned int xsize, unsigned int ysize,
00317                                unsigned char fill_value )
00318 {
00319   image_char image = new_image_char(xsize,ysize); /* create image */
00320   unsigned int N = xsize*ysize;
00321   unsigned int i;
00322 
00323   /* check parameters */
00324   if( image == NULL || image->data == NULL )
00325     error("new_image_char_ini: invalid image.");
00326 
00327   /* initialize */
00328   for(i=0; i<N; i++) image->data[i] = fill_value;
00329 
00330   return image;
00331 }
00332 
00333 /*----------------------------------------------------------------------------*/
00334 /** Free memory used in image_int 'i'.
00335  */
00336 void free_image_int(image_int i)
00337 {
00338   if( i == NULL || i->data == NULL )
00339     error("free_image_int: invalid input image.");
00340   free( (void *) i->data );
00341   free( (void *) i );
00342 }
00343 
00344 /*----------------------------------------------------------------------------*/
00345 /** Create a new image_int of size 'xsize' times 'ysize'.
00346  */
00347 image_int new_image_int(unsigned int xsize, unsigned int ysize)
00348 {
00349   image_int image;
00350 
00351   /* check parameters */
00352   if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");
00353 
00354   /* get memory */
00355   image = (image_int) malloc( sizeof(struct image_int_s) );
00356   if( image == NULL ) error("not enough memory.");
00357   image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );
00358   if( image->data == NULL ) error("not enough memory.");
00359 
00360   /* set image size */
00361   image->xsize = xsize;
00362   image->ysize = ysize;
00363 
00364   return image;
00365 }
00366 
00367 /*----------------------------------------------------------------------------*/
00368 /** Create a new image_int of size 'xsize' times 'ysize',
00369     initialized to the value 'fill_value'.
00370  */
00371 image_int new_image_int_ini( unsigned int xsize, unsigned int ysize,
00372                              int fill_value )
00373 {
00374   image_int image = new_image_int(xsize,ysize); /* create image */
00375   unsigned int N = xsize*ysize;
00376   unsigned int i;
00377 
00378   /* initialize */
00379   for(i=0; i<N; i++) image->data[i] = fill_value;
00380 
00381   return image;
00382 }
00383 
00384 /*----------------------------------------------------------------------------*/
00385 /** Free memory used in image_double 'i'.
00386  */
00387 void free_image_double(image_double i)
00388 {
00389   if( i == NULL || i->data == NULL )
00390     error("free_image_double: invalid input image.");
00391   free( (void *) i->data );
00392   free( (void *) i );
00393 }
00394 
00395 /*----------------------------------------------------------------------------*/
00396 /** Create a new image_double of size 'xsize' times 'ysize'.
00397  */
00398 image_double new_image_double(unsigned int xsize, unsigned int ysize)
00399 {
00400   image_double image;
00401 
00402   /* check parameters */
00403   if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");
00404 
00405   /* get memory */
00406   image = (image_double) malloc( sizeof(struct image_double_s) );
00407   if( image == NULL ) error("not enough memory.");
00408   image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );
00409   if( image->data == NULL ) error("not enough memory.");
00410 
00411   /* set image size */
00412   image->xsize = xsize;
00413   image->ysize = ysize;
00414 
00415   return image;
00416 }
00417 
00418 /*----------------------------------------------------------------------------*/
00419 /** Create a new image_double of size 'xsize' times 'ysize',
00420     initialized to the value 'fill_value'.
00421  */
00422 image_double new_image_double_ini( unsigned int xsize, unsigned int ysize,
00423                                    double fill_value )
00424 {
00425   image_double image = new_image_double(xsize,ysize); /* create image */
00426   unsigned int N = xsize*ysize;
00427   unsigned int i;
00428 
00429   /* initialize */
00430   for(i=0; i<N; i++) image->data[i] = fill_value;
00431 
00432   return image;
00433 }
00434 
00435 
00436 /*----------------------------------------------------------------------------*/
00437 /*----------------------------- Gaussian filter ------------------------------*/
00438 /*----------------------------------------------------------------------------*/
00439 
00440 /*----------------------------------------------------------------------------*/
00441 /** Compute a Gaussian kernel of length 'kernel->dim',
00442     standard deviation 'sigma', and centered at value 'mean'.
00443 
00444     For example, if mean=0.5, the Gaussian will be centered
00445     in the middle point between values 'kernel->values[0]'
00446     and 'kernel->values[1]'.
00447  */
00448 static void gaussian_kernel(ntuple_list kernel, double sigma, double mean)
00449 {
00450   double sum = 0.0;
00451   double val;
00452   unsigned int i;
00453 
00454   /* check parameters */
00455   if( kernel == NULL || kernel->values == NULL )
00456     error("gaussian_kernel: invalid n-tuple 'kernel'.");
00457   if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");
00458 
00459   /* compute Gaussian kernel */
00460   if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);
00461   kernel->size = 1;
00462   for(i=0;i<kernel->dim;i++)
00463     {
00464       val = ( (double) i - mean ) / sigma;
00465       kernel->values[i] = exp( -0.5 * val * val );
00466       sum += kernel->values[i];
00467     }
00468 
00469   /* normalization */
00470   if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
00471 }
00472 
00473 /*----------------------------------------------------------------------------*/
00474 /** Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.
00475 
00476     For example, scale=0.8 will give a result at 80% of the original size.
00477 
00478     The image is convolved with a Gaussian kernel
00479     @f[
00480         G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}
00481     @f]
00482     before the sub-sampling to prevent aliasing.
00483 
00484     The standard deviation sigma given by:
00485     -  sigma = sigma_scale / scale,   if scale <  1.0
00486     -  sigma = sigma_scale,           if scale >= 1.0
00487 
00488     To be able to sub-sample at non-integer steps, some interpolation
00489     is needed. In this implementation, the interpolation is done by
00490     the Gaussian kernel, so both operations (filtering and sampling)
00491     are done at the same time. The Gaussian kernel is computed
00492     centered on the coordinates of the required sample. In this way,
00493     when applied, it gives directly the result of convolving the image
00494     with the kernel and interpolated to that particular position.
00495 
00496     A fast algorithm is done using the separability of the Gaussian
00497     kernel. Applying the 2D Gaussian kernel is equivalent to applying
00498     first a horizontal 1D Gaussian kernel and then a vertical 1D
00499     Gaussian kernel (or the other way round). The reason is that
00500     @f[
00501         G(x,y) = G(x) * G(y)
00502     @f]
00503     where
00504     @f[
00505         G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}}.
00506     @f]
00507     The algorithm first apply a combined Gaussian kernel and sampling
00508     in the x axis, and then the combined Gaussian kernel and sampling
00509     in the y axis.
00510  */
00511 static image_double gaussian_sampler( image_double in, double scale,
00512                                       double sigma_scale )
00513 {
00514   image_double aux,out;
00515   ntuple_list kernel;
00516   unsigned int N,M,h,n,x,y,i;
00517   int xc,yc,j,double_x_size,double_y_size;
00518   double sigma,xx,yy,sum,prec;
00519 
00520   /* check parameters */
00521   if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00522     error("gaussian_sampler: invalid image.");
00523   if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");
00524   if( sigma_scale <= 0.0 )
00525     error("gaussian_sampler: 'sigma_scale' must be positive.");
00526 
00527   /* get memory for images */
00528   if( in->xsize * scale > (double) UINT_MAX ||
00529       in->ysize * scale > (double) UINT_MAX )
00530     error("gaussian_sampler: the output image size exceeds the handled size.");
00531   N = (unsigned int) floor( in->xsize * scale );
00532   M = (unsigned int) floor( in->ysize * scale );
00533   aux = new_image_double(N,in->ysize);
00534   out = new_image_double(N,M);
00535 
00536   /* sigma, kernel size and memory for the kernel */
00537   sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;
00538   /*
00539      The size of the kernel is selected to guarantee that the
00540      the first discarded term is at least 10^prec times smaller
00541      than the central value. For that, h should be larger than x, with
00542        e^(-x^2/2sigma^2) = 1/10^prec.
00543      Then,
00544        x = sigma * sqrt( 2 * prec * ln(10) ).
00545    */
00546   prec = 3.0;
00547   h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
00548   n = 1+2*h; /* kernel size */
00549   kernel = new_ntuple_list(n);
00550 
00551   /* auxiliary double image size variables */
00552   double_x_size = (int) (2 * in->xsize);
00553   double_y_size = (int) (2 * in->ysize);
00554 
00555   /* First subsampling: x axis */
00556   for(x=0;x<aux->xsize;x++)
00557     {
00558       /*
00559          x   is the coordinate in the new image.
00560          xx  is the corresponding x-value in the original size image.
00561          xc  is the integer value, the pixel coordinate of xx.
00562        */
00563       xx = (double) x / scale;
00564       /* coordinate (0.0,0.0) is in the center of pixel (0,0),
00565          so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */
00566       xc = (int) floor( xx + 0.5 );
00567       gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );
00568       /* the kernel must be computed for each x because the fine
00569          offset xx-xc is different in each case */
00570 
00571       for(y=0;y<aux->ysize;y++)
00572         {
00573           sum = 0.0;
00574           for(i=0;i<kernel->dim;i++)
00575             {
00576               j = xc - h + i;
00577 
00578               /* symmetry boundary condition */
00579               while( j < 0 ) j += double_x_size;
00580               while( j >= double_x_size ) j -= double_x_size;
00581               if( j >= (int) in->xsize ) j = double_x_size-1-j;
00582 
00583               sum += in->data[ j + y * in->xsize ] * kernel->values[i];
00584             }
00585           aux->data[ x + y * aux->xsize ] = sum;
00586         }
00587     }
00588 
00589   /* Second subsampling: y axis */
00590   for(y=0;y<out->ysize;y++)
00591     {
00592       /*
00593          y   is the coordinate in the new image.
00594          yy  is the corresponding x-value in the original size image.
00595          yc  is the integer value, the pixel coordinate of xx.
00596        */
00597       yy = (double) y / scale;
00598       /* coordinate (0.0,0.0) is in the center of pixel (0,0),
00599          so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */
00600       yc = (int) floor( yy + 0.5 );
00601       gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );
00602       /* the kernel must be computed for each y because the fine
00603          offset yy-yc is different in each case */
00604 
00605       for(x=0;x<out->xsize;x++)
00606         {
00607           sum = 0.0;
00608           for(i=0;i<kernel->dim;i++)
00609             {
00610               j = yc - h + i;
00611 
00612               /* symmetry boundary condition */
00613               while( j < 0 ) j += double_y_size;
00614               while( j >= double_y_size ) j -= double_y_size;
00615               if( j >= (int) in->ysize ) j = double_y_size-1-j;
00616 
00617               sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];
00618             }
00619           out->data[ x + y * out->xsize ] = sum;
00620         }
00621     }
00622 
00623   /* free memory */
00624   free_ntuple_list(kernel);
00625   free_image_double(aux);
00626 
00627   return out;
00628 }
00629 
00630 
00631 /*----------------------------------------------------------------------------*/
00632 /*--------------------------------- Gradient ---------------------------------*/
00633 /*----------------------------------------------------------------------------*/
00634 
00635 /*----------------------------------------------------------------------------*/
00636 /** Computes the direction of the level line of 'in' at each point.
00637 
00638     The result is:
00639     - an image_double with the angle at each pixel, or NOTDEF if not defined.
00640     - the image_double 'modgrad' (a pointer is passed as argument)
00641       with the gradient magnitude at each point.
00642     - a list of pixels 'list_p' roughly ordered by decreasing
00643       gradient magnitude. (The order is made by classifying points
00644       into bins by gradient magnitude. The parameters 'n_bins' and
00645       'max_grad' specify the number of bins and the gradient modulus
00646       at the highest bin. The pixels in the list would be in
00647       decreasing gradient magnitude, up to a precision of the size of
00648       the bins.)
00649     - a pointer 'mem_p' to the memory used by 'list_p' to be able to
00650       free the memory when it is not used anymore.
00651  */
00652 static image_double ll_angle( image_double in, double threshold,
00653                               struct coorlist ** list_p, void ** mem_p,
00654                               image_double * modgrad, unsigned int n_bins,
00655                               double max_grad )
00656 {
00657   image_double g;
00658   unsigned int n,p,x,y,adr,i;
00659   double com1,com2,gx,gy,norm,norm2;
00660   /* the rest of the variables are used for pseudo-ordering
00661      the gradient magnitude values */
00662   int list_count = 0;
00663   struct coorlist * list;
00664   struct coorlist ** range_l_s; /* array of pointers to start of bin list */
00665   struct coorlist ** range_l_e; /* array of pointers to end of bin list */
00666   struct coorlist * start;
00667   struct coorlist * end;
00668 
00669   /* check parameters */
00670   if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00671     error("ll_angle: invalid image.");
00672   if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");
00673   if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");
00674   if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");
00675   if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");
00676   if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");
00677   if( max_grad <= 0.0 ) error("ll_angle: 'max_grad' must be positive.");
00678 
00679   /* image size shortcuts */
00680   n = in->ysize;
00681   p = in->xsize;
00682 
00683   /* allocate output image */
00684   g = new_image_double(in->xsize,in->ysize);
00685 
00686   /* get memory for the image of gradient modulus */
00687   *modgrad = new_image_double(in->xsize,in->ysize);
00688 
00689   /* get memory for "ordered" list of pixels */
00690   list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );
00691   *mem_p = (void *) list;
00692   range_l_s = (struct coorlist **) calloc( (size_t) n_bins,
00693                                            sizeof(struct coorlist *) );
00694   range_l_e = (struct coorlist **) calloc( (size_t) n_bins,
00695                                            sizeof(struct coorlist *) );
00696   if( list == NULL || range_l_s == NULL || range_l_e == NULL )
00697     error("not enough memory.");
00698   for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;
00699 
00700   /* 'undefined' on the down and right boundaries */
00701   for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;
00702   for(y=0;y<n;y++) g->data[p*y+p-1]   = NOTDEF;
00703 
00704   /* compute gradient on the remaining pixels */
00705   for(x=0;x<p-1;x++)
00706     for(y=0;y<n-1;y++)
00707       {
00708         adr = y*p+x;
00709 
00710         /*
00711            Norm 2 computation using 2x2 pixel window:
00712              A B
00713              C D
00714            and
00715              com1 = D-A,  com2 = B-C.
00716            Then
00717              gx = B+D - (A+C)   horizontal difference
00718              gy = C+D - (A+B)   vertical difference
00719            com1 and com2 are just to avoid 2 additions.
00720          */
00721         com1 = in->data[adr+p+1] - in->data[adr];
00722         com2 = in->data[adr+1]   - in->data[adr+p];
00723 
00724         gx = com1+com2; /* gradient x component */
00725         gy = com1-com2; /* gradient y component */
00726         norm2 = gx*gx+gy*gy;
00727         norm = sqrt( norm2 / 4.0 ); /* gradient norm */
00728 
00729         (*modgrad)->data[adr] = norm; /* store gradient norm */
00730 
00731         if( norm <= threshold ) /* norm too small, gradient no defined */
00732           g->data[adr] = NOTDEF; /* gradient angle not defined */
00733         else
00734           {
00735             /* gradient angle computation */
00736             g->data[adr] = atan2(gx,-gy);
00737 
00738             /* store the point in the right bin according to its norm */
00739             i = (unsigned int) (norm * (double) n_bins / max_grad);
00740             if( i >= n_bins ) i = n_bins-1;
00741             if( range_l_e[i] == NULL )
00742               range_l_s[i] = range_l_e[i] = list+list_count++;
00743             else
00744               {
00745                 range_l_e[i]->next = list+list_count;
00746                 range_l_e[i] = list+list_count++;
00747               }
00748             range_l_e[i]->x = (int) x;
00749             range_l_e[i]->y = (int) y;
00750             range_l_e[i]->next = NULL;
00751           }
00752       }
00753 
00754   /* Make the list of pixels (almost) ordered by norm value.
00755      It starts by the larger bin, so the list starts by the
00756      pixels with higher gradient value. Pixels would be ordered
00757      by norm value, up to a precision given by max_grad/n_bins.
00758    */
00759   for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);
00760   start = range_l_s[i];
00761   end = range_l_e[i];
00762   if( start != NULL )
00763     for(i--;i>0; i--)
00764       if( range_l_s[i] != NULL )
00765         {
00766           end->next = range_l_s[i];
00767           end = range_l_e[i];
00768         }
00769   *list_p = start;
00770 
00771   /* free memory */
00772   free( (void *) range_l_s );
00773   free( (void *) range_l_e );
00774 
00775   return g;
00776 }
00777 
00778 /*----------------------------------------------------------------------------*/
00779 /** Is point (x,y) aligned to angle theta, up to precision 'prec'?
00780  */
00781 static int isaligned( int x, int y, image_double angles, double theta,
00782                       double prec )
00783 {
00784   double a;
00785 
00786   /* check parameters */
00787   if( angles == NULL || angles->data == NULL )
00788     error("isaligned: invalid image 'angles'.");
00789   if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
00790     error("isaligned: (x,y) out of the image.");
00791   if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");
00792 
00793   /* angle at pixel (x,y) */
00794   a = angles->data[ x + y * angles->xsize ];
00795 
00796   /* pixels whose level-line angle is not defined
00797      are considered as NON-aligned */
00798   if( a == NOTDEF ) return FALSE;  /* there is no need to call the function
00799                                       'double_equal' here because there is
00800                                       no risk of problems related to the
00801                                       comparison doubles, we are only
00802                                       interested in the exact NOTDEF value */
00803 
00804   /* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */
00805   theta -= a;
00806   if( theta < 0.0 ) theta = -theta;
00807   if( theta > M_3_2_PI )
00808     {
00809       theta -= M_2__PI;
00810       if( theta < 0.0 ) theta = -theta;
00811     }
00812 
00813   return theta < prec;
00814 }
00815 
00816 /*----------------------------------------------------------------------------*/
00817 /** Absolute value angle difference.
00818  */
00819 static double angle_diff(double a, double b)
00820 {
00821   a -= b;
00822   while( a <= -M_PI ) a += M_2__PI;
00823   while( a >   M_PI ) a -= M_2__PI;
00824   if( a < 0.0 ) a = -a;
00825   return a;
00826 }
00827 
00828 /*----------------------------------------------------------------------------*/
00829 /** Signed angle difference.
00830  */
00831 static double angle_diff_signed(double a, double b)
00832 {
00833   a -= b;
00834   while( a <= -M_PI ) a += M_2__PI;
00835   while( a >   M_PI ) a -= M_2__PI;
00836   return a;
00837 }
00838 
00839 
00840 /*----------------------------------------------------------------------------*/
00841 /*----------------------------- NFA computation ------------------------------*/
00842 /*----------------------------------------------------------------------------*/
00843 
00844 /*----------------------------------------------------------------------------*/
00845 /** Computes the natural logarithm of the absolute value of
00846     the gamma function of x using the Lanczos approximation.
00847     See http://www.rskey.org/gamma.htm
00848 
00849     The formula used is
00850     @f[
00851       \Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) }
00852                   (x+5.5)^{x+0.5} e^{-(x+5.5)}
00853     @f]
00854     so
00855     @f[
00856       \log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right)
00857                       + (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n)
00858     @f]
00859     and
00860       q0 = 75122.6331530,
00861       q1 = 80916.6278952,
00862       q2 = 36308.2951477,
00863       q3 = 8687.24529705,
00864       q4 = 1168.92649479,
00865       q5 = 83.8676043424,
00866       q6 = 2.50662827511.
00867  */
00868 static double log_gamma_lanczos(double x)
00869 {
00870   static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
00871                          8687.24529705, 1168.92649479, 83.8676043424,
00872                          2.50662827511 };
00873   double a = (x+0.5) * log(x+5.5) - (x+5.5);
00874   double b = 0.0;
00875   int n;
00876 
00877   for(n=0;n<7;n++)
00878     {
00879       a -= log( x + (double) n );
00880       b += q[n] * pow( x, (double) n );
00881     }
00882   return a + log(b);
00883 }
00884 
00885 /*----------------------------------------------------------------------------*/
00886 /** Computes the natural logarithm of the absolute value of
00887     the gamma function of x using Windschitl method.
00888     See http://www.rskey.org/gamma.htm
00889 
00890     The formula used is
00891     @f[
00892         \Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e}
00893                     \sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x
00894     @f]
00895     so
00896     @f[
00897         \log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x
00898                       + 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right).
00899     @f]
00900     This formula is a good approximation when x > 15.
00901  */
00902 static double log_gamma_windschitl(double x)
00903 {
00904   return 0.918938533204673 + (x-0.5)*log(x) - x
00905          + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
00906 }
00907 
00908 /*----------------------------------------------------------------------------*/
00909 /** Computes the natural logarithm of the absolute value of
00910     the gamma function of x. When x>15 use log_gamma_windschitl(),
00911     otherwise use log_gamma_lanczos().
00912  */
00913 #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
00914 
00915 /*----------------------------------------------------------------------------*/
00916 /** Size of the table to store already computed inverse values.
00917  */
00918 #define TABSIZE 100000
00919 
00920 /*----------------------------------------------------------------------------*/
00921 /** Computes -log10(NFA).
00922 
00923     NFA stands for Number of False Alarms:
00924     @f[
00925         \mathrm{NFA} = NT \cdot B(n,k,p)
00926     @f]
00927 
00928     - NT       - number of tests
00929     - B(n,k,p) - tail of binomial distribution with parameters n,k and p:
00930     @f[
00931         B(n,k,p) = \sum_{j=k}^n
00932                    \left(\begin{array}{c}n\\j\end{array}\right)
00933                    p^{j} (1-p)^{n-j}
00934     @f]
00935 
00936     The value -log10(NFA) is equivalent but more intuitive than NFA:
00937     - -1 corresponds to 10 mean false alarms
00938     -  0 corresponds to 1 mean false alarm
00939     -  1 corresponds to 0.1 mean false alarms
00940     -  2 corresponds to 0.01 mean false alarms
00941     -  ...
00942 
00943     Used this way, the bigger the value, better the detection,
00944     and a logarithmic scale is used.
00945 
00946     @param n,k,p binomial parameters.
00947     @param logNT logarithm of Number of Tests
00948 
00949     The computation is based in the gamma function by the following
00950     relation:
00951     @f[
00952         \left(\begin{array}{c}n\\k\end{array}\right)
00953         = \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }.
00954     @f]
00955     We use efficient algorithms to compute the logarithm of
00956     the gamma function.
00957 
00958     To make the computation faster, not all the sum is computed, part
00959     of the terms are neglected based on a bound to the error obtained
00960     (an error of 10% in the result is accepted).
00961  */
00962 static double nfa(int n, int k, double p, double logNT)
00963 {
00964   static double inv[TABSIZE];   /* table to keep computed inverse values */
00965   double tolerance = 0.1;       /* an error of 10% in the result is accepted */
00966   double log1term,term,bin_term,mult_term,bin_tail,err,p_term;
00967   int i;
00968 
00969   /* check parameters */
00970   if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )
00971     error("nfa: wrong n, k or p values.");
00972 
00973   /* trivial cases */
00974   if( n==0 || k==0 ) return -logNT;
00975   if( n==k ) return -logNT - (double) n * log10(p);
00976 
00977   /* probability term */
00978   p_term = p / (1.0-p);
00979 
00980   /* compute the first term of the series */
00981   /*
00982      binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
00983      where bincoef(n,i) are the binomial coefficients.
00984      But
00985        bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
00986      We use this to compute the first term. Actually the log of it.
00987    */
00988   log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )
00989            - log_gamma( (double) (n-k) + 1.0 )
00990            + (double) k * log(p) + (double) (n-k) * log(1.0-p);
00991   term = exp(log1term);
00992 
00993   /* in some cases no more computations are needed */
00994   if( double_equal(term,0.0) )              /* the first term is almost zero */
00995     {
00996       if( (double) k > (double) n * p )     /* at begin or end of the tail?  */
00997         return -log1term / M_LN10 - logNT;  /* end: use just the first term  */
00998       else
00999         return -logNT;                      /* begin: the tail is roughly 1  */
01000     }
01001 
01002   /* compute more terms if needed */
01003   bin_tail = term;
01004   for(i=k+1;i<=n;i++)
01005     {
01006       /*
01007          As
01008            term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
01009          and
01010            bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,
01011          then,
01012            term_i / term_i-1 = (n-i+1)/i * p/(1-p)
01013          and
01014            term_i = term_i-1 * (n-i+1)/i * p/(1-p).
01015          1/i is stored in a table as they are computed,
01016          because divisions are expensive.
01017          p/(1-p) is computed only once and stored in 'p_term'.
01018        */
01019       bin_term = (double) (n-i+1) * ( i<TABSIZE ?
01020                    ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :
01021                    1.0 / (double) i );
01022 
01023       mult_term = bin_term * p_term;
01024       term *= mult_term;
01025       bin_tail += term;
01026       if(bin_term<1.0)
01027         {
01028           /* When bin_term<1 then mult_term_j<mult_term_i for j>i.
01029              Then, the error on the binomial tail when truncated at
01030              the i term can be bounded by a geometric series of form
01031              term_i * sum mult_term_i^j.                            */
01032           err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /
01033                          (1.0-mult_term) - 1.0 );
01034 
01035           /* One wants an error at most of tolerance*final_result, or:
01036              tolerance * abs(-log10(bin_tail)-logNT).
01037              Now, the error that can be accepted on bin_tail is
01038              given by tolerance*final_result divided by the derivative
01039              of -log10(x) when x=bin_tail. that is:
01040              tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
01041              Finally, we truncate the tail if the error is less than:
01042              tolerance * abs(-log10(bin_tail)-logNT) * bin_tail        */
01043           if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;
01044         }
01045     }
01046   return -log10(bin_tail) - logNT;
01047 }
01048 
01049 
01050 /*----------------------------------------------------------------------------*/
01051 /*--------------------------- Rectangle structure ----------------------------*/
01052 /*----------------------------------------------------------------------------*/
01053 
01054 /*----------------------------------------------------------------------------*/
01055 /** Rectangle structure: line segment with width.
01056  */
01057 struct rect
01058 {
01059   double x1,y1,x2,y2;  /* first and second point of the line segment */
01060   double width;        /* rectangle width */
01061   double x,y;          /* center of the rectangle */
01062   double theta;        /* angle */
01063   double dx,dy;        /* vector with the line segment angle */
01064   double prec;         /* tolerance angle */
01065   double p;            /* probability of a point with angle within 'prec' */
01066 };
01067 
01068 /*----------------------------------------------------------------------------*/
01069 /** Copy one rectangle structure to another.
01070  */
01071 static void rect_copy(struct rect * in, struct rect * out)
01072 {
01073   /* check parameters */
01074   if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");
01075 
01076   /* copy values */
01077   out->x1 = in->x1;
01078   out->y1 = in->y1;
01079   out->x2 = in->x2;
01080   out->y2 = in->y2;
01081   out->width = in->width;
01082   out->x = in->x;
01083   out->y = in->y;
01084   out->theta = in->theta;
01085   out->dx = in->dx;
01086   out->dy = in->dy;
01087   out->prec = in->prec;
01088   out->p = in->p;
01089 }
01090 
01091 /*----------------------------------------------------------------------------*/
01092 /** Rectangle points iterator.
01093 
01094     The integer coordinates of pixels inside a rectangle are
01095     iteratively explored. This structure keep track of the process and
01096     functions ri_ini(), ri_inc(), ri_end(), and ri_del() are used in
01097     the process. An example of how to use the iterator is as follows:
01098     \code
01099 
01100       struct rect * rec = XXX; // some rectangle
01101       rect_iter * i;
01102       for( i=ri_ini(rec); !ri_end(i); ri_inc(i) )
01103         {
01104           // your code, using 'i->x' and 'i->y' as coordinates
01105         }
01106       ri_del(i); // delete iterator
01107 
01108     \endcode
01109     The pixels are explored 'column' by 'column', where we call
01110     'column' a set of pixels with the same x value that are inside the
01111     rectangle. The following is an schematic representation of a
01112     rectangle, the 'column' being explored is marked by colons, and
01113     the current pixel being explored is 'x,y'.
01114     \verbatim
01115 
01116               vx[1],vy[1]
01117                  *   *
01118                 *       *
01119                *           *
01120               *               ye
01121              *                :  *
01122         vx[0],vy[0]           :     *
01123                *              :        *
01124                   *          x,y          *
01125                      *        :              *
01126                         *     :            vx[2],vy[2]
01127                            *  :                *
01128         y                     ys              *
01129         ^                        *           *
01130         |                           *       *
01131         |                              *   *
01132         +---> x                      vx[3],vy[3]
01133 
01134     \endverbatim
01135     The first 'column' to be explored is the one with the smaller x
01136     value. Each 'column' is explored starting from the pixel of the
01137     'column' (inside the rectangle) with the smaller y value.
01138 
01139     The four corners of the rectangle are stored in order that rotates
01140     around the corners at the arrays 'vx[]' and 'vy[]'. The first
01141     point is always the one with smaller x value.
01142 
01143     'x' and 'y' are the coordinates of the pixel being explored. 'ys'
01144     and 'ye' are the start and end values of the current column being
01145     explored. So, 'ys' < 'ye'.
01146  */
01147 typedef struct
01148 {
01149   double vx[4];  /* rectangle's corner X coordinates in circular order */
01150   double vy[4];  /* rectangle's corner Y coordinates in circular order */
01151   double ys,ye;  /* start and end Y values of current 'column' */
01152   int x,y;       /* coordinates of currently explored pixel */
01153 } rect_iter;
01154 
01155 /*----------------------------------------------------------------------------*/
01156 /** Interpolate y value corresponding to 'x' value given, in
01157     the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller
01158     of 'y1' and 'y2'.
01159 
01160     The following restrictions are required:
01161     - x1 <= x2
01162     - x1 <= x
01163     - x  <= x2
01164  */
01165 static double inter_low(double x, double x1, double y1, double x2, double y2)
01166 {
01167   /* check parameters */
01168   if( x1 > x2 || x < x1 || x > x2 )
01169     error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01170 
01171   /* interpolation */
01172   if( double_equal(x1,x2) && y1<y2 ) return y1;
01173   if( double_equal(x1,x2) && y1>y2 ) return y2;
01174   return y1 + (x-x1) * (y2-y1) / (x2-x1);
01175 }
01176 
01177 /*----------------------------------------------------------------------------*/
01178 /** Interpolate y value corresponding to 'x' value given, in
01179     the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger
01180     of 'y1' and 'y2'.
01181 
01182     The following restrictions are required:
01183     - x1 <= x2
01184     - x1 <= x
01185     - x  <= x2
01186  */
01187 static double inter_hi(double x, double x1, double y1, double x2, double y2)
01188 {
01189   /* check parameters */
01190   if( x1 > x2 || x < x1 || x > x2 )
01191     error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01192 
01193   /* interpolation */
01194   if( double_equal(x1,x2) && y1<y2 ) return y2;
01195   if( double_equal(x1,x2) && y1>y2 ) return y1;
01196   return y1 + (x-x1) * (y2-y1) / (x2-x1);
01197 }
01198 
01199 /*----------------------------------------------------------------------------*/
01200 /** Free memory used by a rectangle iterator.
01201  */
01202 static void ri_del(rect_iter * iter)
01203 {
01204   if( iter == NULL ) error("ri_del: NULL iterator.");
01205   free( (void *) iter );
01206 }
01207 
01208 /*----------------------------------------------------------------------------*/
01209 /** Check if the iterator finished the full iteration.
01210 
01211     See details in \ref rect_iter
01212  */
01213 static int ri_end(rect_iter * i)
01214 {
01215   /* check input */
01216   if( i == NULL ) error("ri_end: NULL iterator.");
01217 
01218   /* if the current x value is larger than the larger
01219      x value in the rectangle (vx[2]), we know the full
01220      exploration of the rectangle is finished. */
01221   return (double)(i->x) > i->vx[2];
01222 }
01223 
01224 /*----------------------------------------------------------------------------*/
01225 /** Increment a rectangle iterator.
01226 
01227     See details in \ref rect_iter
01228  */
01229 static void ri_inc(rect_iter * i)
01230 {
01231   /* check input */
01232   if( i == NULL ) error("ri_inc: NULL iterator.");
01233 
01234   /* if not at end of exploration,
01235      increase y value for next pixel in the 'column' */
01236   if( !ri_end(i) ) i->y++;
01237 
01238   /* if the end of the current 'column' is reached,
01239      and it is not the end of exploration,
01240      advance to the next 'column' */
01241   while( (double) (i->y) > i->ye && !ri_end(i) )
01242     {
01243       /* increase x, next 'column' */
01244       i->x++;
01245 
01246       /* if end of exploration, return */
01247       if( ri_end(i) ) return;
01248 
01249       /* update lower y limit (start) for the new 'column'.
01250 
01251          We need to interpolate the y value that corresponds to the
01252          lower side of the rectangle. The first thing is to decide if
01253          the corresponding side is
01254 
01255            vx[0],vy[0] to vx[3],vy[3] or
01256            vx[3],vy[3] to vx[2],vy[2]
01257 
01258          Then, the side is interpolated for the x value of the
01259          'column'. But, if the side is vertical (as it could happen if
01260          the rectangle is vertical and we are dealing with the first
01261          or last 'columns') then we pick the lower value of the side
01262          by using 'inter_low'.
01263        */
01264       if( (double) i->x < i->vx[3] )
01265         i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);
01266       else
01267         i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);
01268 
01269       /* update upper y limit (end) for the new 'column'.
01270 
01271          We need to interpolate the y value that corresponds to the
01272          upper side of the rectangle. The first thing is to decide if
01273          the corresponding side is
01274 
01275            vx[0],vy[0] to vx[1],vy[1] or
01276            vx[1],vy[1] to vx[2],vy[2]
01277 
01278          Then, the side is interpolated for the x value of the
01279          'column'. But, if the side is vertical (as it could happen if
01280          the rectangle is vertical and we are dealing with the first
01281          or last 'columns') then we pick the lower value of the side
01282          by using 'inter_low'.
01283        */
01284       if( (double)i->x < i->vx[1] )
01285         i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);
01286       else
01287         i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);
01288 
01289       /* new y */
01290       i->y = (int) ceil(i->ys);
01291     }
01292 }
01293 
01294 /*----------------------------------------------------------------------------*/
01295 /** Create and initialize a rectangle iterator.
01296 
01297     See details in \ref rect_iter
01298  */
01299 static rect_iter * ri_ini(struct rect * r)
01300 {
01301   double vx[4],vy[4];
01302   int n,offset;
01303   rect_iter * i;
01304 
01305   /* check parameters */
01306   if( r == NULL ) error("ri_ini: invalid rectangle.");
01307 
01308   /* get memory */
01309   i = (rect_iter *) malloc(sizeof(rect_iter));
01310   if( i == NULL ) error("ri_ini: Not enough memory.");
01311 
01312   /* build list of rectangle corners ordered
01313      in a circular way around the rectangle */
01314   vx[0] = r->x1 - r->dy * r->width / 2.0;
01315   vy[0] = r->y1 + r->dx * r->width / 2.0;
01316   vx[1] = r->x2 - r->dy * r->width / 2.0;
01317   vy[1] = r->y2 + r->dx * r->width / 2.0;
01318   vx[2] = r->x2 + r->dy * r->width / 2.0;
01319   vy[2] = r->y2 - r->dx * r->width / 2.0;
01320   vx[3] = r->x1 + r->dy * r->width / 2.0;
01321   vy[3] = r->y1 - r->dx * r->width / 2.0;
01322 
01323   /* compute rotation of index of corners needed so that the first
01324      point has the smaller x.
01325 
01326      if one side is vertical, thus two corners have the same smaller x
01327      value, the one with the largest y value is selected as the first.
01328    */
01329   if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;
01330   else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;
01331   else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;
01332   else offset = 3;
01333 
01334   /* apply rotation of index. */
01335   for(n=0; n<4; n++)
01336     {
01337       i->vx[n] = vx[(offset+n)%4];
01338       i->vy[n] = vy[(offset+n)%4];
01339     }
01340 
01341   /* Set a initial condition.
01342 
01343      The values are set to values that will cause 'ri_inc' (that will
01344      be called immediately) to initialize correctly the first 'column'
01345      and compute the limits 'ys' and 'ye'.
01346 
01347      'y' is set to the integer value of vy[0], the starting corner.
01348 
01349      'ys' and 'ye' are set to very small values, so 'ri_inc' will
01350      notice that it needs to start a new 'column'.
01351 
01352      The smaller integer coordinate inside of the rectangle is
01353      'ceil(vx[0])'. The current 'x' value is set to that value minus
01354      one, so 'ri_inc' (that will increase x by one) will advance to
01355      the first 'column'.
01356    */
01357   i->x = (int) ceil(i->vx[0]) - 1;
01358   i->y = (int) ceil(i->vy[0]);
01359   i->ys = i->ye = -DBL_MAX;
01360 
01361   /* advance to the first pixel */
01362   ri_inc(i);
01363 
01364   return i;
01365 }
01366 
01367 /*----------------------------------------------------------------------------*/
01368 /** Compute a rectangle's NFA value.
01369  */
01370 static double rect_nfa(struct rect * rec, image_double angles, double logNT)
01371 {
01372   rect_iter * i;
01373   int pts = 0;
01374   int alg = 0;
01375 
01376   /* check parameters */
01377   if( rec == NULL ) error("rect_nfa: invalid rectangle.");
01378   if( angles == NULL ) error("rect_nfa: invalid 'angles'.");
01379 
01380   /* compute the total number of pixels and of aligned points in 'rec' */
01381   for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */
01382     if( i->x >= 0 && i->y >= 0 &&
01383         i->x < (int) angles->xsize && i->y < (int) angles->ysize )
01384       {
01385         ++pts; /* total number of pixels counter */
01386         if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )
01387           ++alg; /* aligned points counter */
01388       }
01389   ri_del(i); /* delete iterator */
01390 
01391   return nfa(pts,alg,rec->p,logNT); /* compute NFA value */
01392 }
01393 
01394 
01395 /*----------------------------------------------------------------------------*/
01396 /*---------------------------------- Regions ---------------------------------*/
01397 /*----------------------------------------------------------------------------*/
01398 
01399 /*----------------------------------------------------------------------------*/
01400 /** Compute region's angle as the principal inertia axis of the region.
01401 
01402     The following is the region inertia matrix A:
01403     @f[
01404 
01405         A = \left(\begin{array}{cc}
01406                                     Ixx & Ixy \\
01407                                     Ixy & Iyy \\
01408              \end{array}\right)
01409 
01410     @f]
01411     where
01412 
01413       Ixx =   sum_i G(i).(y_i - cx)^2
01414 
01415       Iyy =   sum_i G(i).(x_i - cy)^2
01416 
01417       Ixy = - sum_i G(i).(x_i - cx).(y_i - cy)
01418 
01419     and
01420     - G(i) is the gradient norm at pixel i, used as pixel's weight.
01421     - x_i and y_i are the coordinates of pixel i.
01422     - cx and cy are the coordinates of the center of th region.
01423 
01424     lambda1 and lambda2 are the eigenvalues of matrix A,
01425     with lambda1 >= lambda2. They are found by solving the
01426     characteristic polynomial:
01427 
01428       det( lambda I - A) = 0
01429 
01430     that gives:
01431 
01432       lambda1 = ( Ixx + Iyy + sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2
01433 
01434       lambda2 = ( Ixx + Iyy - sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2
01435 
01436     To get the line segment direction we want to get the angle the
01437     eigenvector assotiated to the smaller eigenvalue. We have to solve
01438     a,b in:
01439 
01440       a.Ixx + b.Ixy = a.lambda2
01441 
01442       a.Ixy + b.Iyy = b.lambda2
01443 
01444     We want the angle theta = atan(b/a). It can be computed with
01445     any of the two equations:
01446 
01447       theta = atan( (lambda2-Ixx) / Ixy )
01448 
01449     or
01450 
01451       theta = atan( Ixy / (lambda2-Iyy) )
01452 
01453     When |Ixx| > |Iyy| we use the first, otherwise the second (just to
01454     get better numeric precision).
01455  */
01456 static double get_theta( struct point * reg, int reg_size, double x, double y,
01457                          image_double modgrad, double reg_angle, double prec )
01458 {
01459   double lambda,theta,weight;
01460   double Ixx = 0.0;
01461   double Iyy = 0.0;
01462   double Ixy = 0.0;
01463   int i;
01464 
01465   /* check parameters */
01466   if( reg == NULL ) error("get_theta: invalid region.");
01467   if( reg_size <= 1 ) error("get_theta: region size <= 1.");
01468   if( modgrad == NULL || modgrad->data == NULL )
01469     error("get_theta: invalid 'modgrad'.");
01470   if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");
01471 
01472   /* compute inertia matrix */
01473   for(i=0; i<reg_size; i++)
01474     {
01475       weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01476       Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;
01477       Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;
01478       Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;
01479     }
01480   if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )
01481     error("get_theta: null inertia matrix.");
01482 
01483   /* compute smallest eigenvalue */
01484   lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );
01485 
01486   /* compute angle */
01487   theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);
01488 
01489   /* The previous procedure don't cares about orientation,
01490      so it could be wrong by 180 degrees. Here is corrected if necessary. */
01491   if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;
01492 
01493   return theta;
01494 }
01495 
01496 /*----------------------------------------------------------------------------*/
01497 /** Computes a rectangle that covers a region of points.
01498  */
01499 static void region2rect( struct point * reg, int reg_size,
01500                          image_double modgrad, double reg_angle,
01501                          double prec, double p, struct rect * rec )
01502 {
01503   double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;
01504   int i;
01505 
01506   /* check parameters */
01507   if( reg == NULL ) error("region2rect: invalid region.");
01508   if( reg_size <= 1 ) error("region2rect: region size <= 1.");
01509   if( modgrad == NULL || modgrad->data == NULL )
01510     error("region2rect: invalid image 'modgrad'.");
01511   if( rec == NULL ) error("region2rect: invalid 'rec'.");
01512 
01513   /* center of the region:
01514 
01515      It is computed as the weighted sum of the coordinates
01516      of all the pixels in the region. The norm of the gradient
01517      is used as the weight of a pixel. The sum is as follows:
01518        cx = \sum_i G(i).x_i
01519        cy = \sum_i G(i).y_i
01520      where G(i) is the norm of the gradient of pixel i
01521      and x_i,y_i are its coordinates.
01522    */
01523   x = y = sum = 0.0;
01524   for(i=0; i<reg_size; i++)
01525     {
01526       weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01527       x += (double) reg[i].x * weight;
01528       y += (double) reg[i].y * weight;
01529       sum += weight;
01530     }
01531   if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");
01532   x /= sum;
01533   y /= sum;
01534 
01535   /* theta */
01536   theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);
01537 
01538   /* length and width:
01539 
01540      'l' and 'w' are computed as the distance from the center of the
01541      region to pixel i, projected along the rectangle axis (dx,dy) and
01542      to the orthogonal axis (-dy,dx), respectively.
01543 
01544      The length of the rectangle goes from l_min to l_max, where l_min
01545      and l_max are the minimum and maximum values of l in the region.
01546      Analogously, the width is selected from w_min to w_max, where
01547      w_min and w_max are the minimum and maximum of w for the pixels
01548      in the region.
01549    */
01550   dx = cos(theta);
01551   dy = sin(theta);
01552   l_min = l_max = w_min = w_max = 0.0;
01553   for(i=0; i<reg_size; i++)
01554     {
01555       l =  ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;
01556       w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;
01557 
01558       if( l > l_max ) l_max = l;
01559       if( l < l_min ) l_min = l;
01560       if( w > w_max ) w_max = w;
01561       if( w < w_min ) w_min = w;
01562     }
01563 
01564   /* store values */
01565   rec->x1 = x + l_min * dx;
01566   rec->y1 = y + l_min * dy;
01567   rec->x2 = x + l_max * dx;
01568   rec->y2 = y + l_max * dy;
01569   rec->width = w_max - w_min;
01570   rec->x = x;
01571   rec->y = y;
01572   rec->theta = theta;
01573   rec->dx = dx;
01574   rec->dy = dy;
01575   rec->prec = prec;
01576   rec->p = p;
01577 
01578   /* we impose a minimal width of one pixel
01579 
01580      A sharp horizontal or vertical step would produce a perfectly
01581      horizontal or vertical region. The width computed would be
01582      zero. But that corresponds to a one pixels width transition in
01583      the image.
01584    */
01585   if( rec->width < 1.0 ) rec->width = 1.0;
01586 }
01587 
01588 /*----------------------------------------------------------------------------*/
01589 /** Build a region of pixels that share the same angle, up to a
01590     tolerance 'prec', starting at point (x,y).
01591  */
01592 static void region_grow( int x, int y, image_double angles, struct point * reg,
01593                          int * reg_size, double * reg_angle, image_char used,
01594                          double prec )
01595 {
01596   double sumdx,sumdy;
01597   int xx,yy,i;
01598 
01599   /* check parameters */
01600   if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
01601     error("region_grow: (x,y) out of the image.");
01602   if( angles == NULL || angles->data == NULL )
01603     error("region_grow: invalid image 'angles'.");
01604   if( reg == NULL ) error("region_grow: invalid 'reg'.");
01605   if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");
01606   if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");
01607   if( used == NULL || used->data == NULL )
01608     error("region_grow: invalid image 'used'.");
01609 
01610   /* first point of the region */
01611   *reg_size = 1;
01612   reg[0].x = x;
01613   reg[0].y = y;
01614   *reg_angle = angles->data[x+y*angles->xsize];  /* region's angle */
01615   sumdx = cos(*reg_angle);
01616   sumdy = sin(*reg_angle);
01617   used->data[x+y*used->xsize] = USED;
01618 
01619   /* try neighbors as new region points */
01620   for(i=0; i<*reg_size; i++)
01621     for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)
01622       for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)
01623         if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&
01624             used->data[xx+yy*used->xsize] != USED &&
01625             isaligned(xx,yy,angles,*reg_angle,prec) )
01626           {
01627             /* add point */
01628             used->data[xx+yy*used->xsize] = USED;
01629             reg[*reg_size].x = xx;
01630             reg[*reg_size].y = yy;
01631             ++(*reg_size);
01632 
01633             /* update region's angle */
01634             sumdx += cos( angles->data[xx+yy*angles->xsize] );
01635             sumdy += sin( angles->data[xx+yy*angles->xsize] );
01636             *reg_angle = atan2(sumdy,sumdx);
01637           }
01638 }
01639 
01640 /*----------------------------------------------------------------------------*/
01641 /** Try some rectangles variations to improve NFA value. Only if the
01642     rectangle is not meaningful (i.e., log_nfa <= eps).
01643  */
01644 static double rect_improve( struct rect * rec, image_double angles,
01645                             double logNT, double eps )
01646 {
01647   struct rect r;
01648   double log_nfa,log_nfa_new;
01649   double delta = 0.5;
01650   double delta_2 = delta / 2.0;
01651   int n;
01652 
01653   log_nfa = rect_nfa(rec,angles,logNT);
01654 
01655   if( log_nfa > eps ) return log_nfa;
01656 
01657   /* try finer precisions */
01658   rect_copy(rec,&r);
01659   for(n=0; n<5; n++)
01660     {
01661       r.p /= 2.0;
01662       r.prec = r.p * M_PI;
01663       log_nfa_new = rect_nfa(&r,angles,logNT);
01664       if( log_nfa_new > log_nfa )
01665         {
01666           log_nfa = log_nfa_new;
01667           rect_copy(&r,rec);
01668         }
01669     }
01670 
01671   if( log_nfa > eps ) return log_nfa;
01672 
01673   /* try to reduce width */
01674   rect_copy(rec,&r);
01675   for(n=0; n<5; n++)
01676     {
01677       if( (r.width - delta) >= 0.5 )
01678         {
01679           r.width -= delta;
01680           log_nfa_new = rect_nfa(&r,angles,logNT);
01681           if( log_nfa_new > log_nfa )
01682             {
01683               rect_copy(&r,rec);
01684               log_nfa = log_nfa_new;
01685             }
01686         }
01687     }
01688 
01689   if( log_nfa > eps ) return log_nfa;
01690 
01691   /* try to reduce one side of the rectangle */
01692   rect_copy(rec,&r);
01693   for(n=0; n<5; n++)
01694     {
01695       if( (r.width - delta) >= 0.5 )
01696         {
01697           r.x1 += -r.dy * delta_2;
01698           r.y1 +=  r.dx * delta_2;
01699           r.x2 += -r.dy * delta_2;
01700           r.y2 +=  r.dx * delta_2;
01701           r.width -= delta;
01702           log_nfa_new = rect_nfa(&r,angles,logNT);
01703           if( log_nfa_new > log_nfa )
01704             {
01705               rect_copy(&r,rec);
01706               log_nfa = log_nfa_new;
01707             }
01708         }
01709     }
01710 
01711   if( log_nfa > eps ) return log_nfa;
01712 
01713   /* try to reduce the other side of the rectangle */
01714   rect_copy(rec,&r);
01715   for(n=0; n<5; n++)
01716     {
01717       if( (r.width - delta) >= 0.5 )
01718         {
01719           r.x1 -= -r.dy * delta_2;
01720           r.y1 -=  r.dx * delta_2;
01721           r.x2 -= -r.dy * delta_2;
01722           r.y2 -=  r.dx * delta_2;
01723           r.width -= delta;
01724           log_nfa_new = rect_nfa(&r,angles,logNT);
01725           if( log_nfa_new > log_nfa )
01726             {
01727               rect_copy(&r,rec);
01728               log_nfa = log_nfa_new;
01729             }
01730         }
01731     }
01732 
01733   if( log_nfa > eps ) return log_nfa;
01734 
01735   /* try even finer precisions */
01736   rect_copy(rec,&r);
01737   for(n=0; n<5; n++)
01738     {
01739       r.p /= 2.0;
01740       r.prec = r.p * M_PI;
01741       log_nfa_new = rect_nfa(&r,angles,logNT);
01742       if( log_nfa_new > log_nfa )
01743         {
01744           log_nfa = log_nfa_new;
01745           rect_copy(&r,rec);
01746         }
01747     }
01748 
01749   return log_nfa;
01750 }
01751 
01752 /*----------------------------------------------------------------------------*/
01753 /** Reduce the region size, by elimination the points far from the
01754     starting point, until that leads to rectangle with the right
01755     density of region points or to discard the region if too small.
01756  */
01757 static int reduce_region_radius( struct point * reg, int * reg_size,
01758                                  image_double modgrad, double reg_angle,
01759                                  double prec, double p, struct rect * rec,
01760                                  image_char used, image_double angles,
01761                                  double density_th )
01762 {
01763   double density,rad1,rad2,rad,xc,yc;
01764   int i;
01765 
01766   /* check parameters */
01767   if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");
01768   if( reg_size == NULL )
01769     error("reduce_region_radius: invalid pointer 'reg_size'.");
01770   if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");
01771   if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");
01772   if( used == NULL || used->data == NULL )
01773     error("reduce_region_radius: invalid image 'used'.");
01774   if( angles == NULL || angles->data == NULL )
01775     error("reduce_region_radius: invalid image 'angles'.");
01776 
01777   /* compute region points density */
01778   density = (double) *reg_size /
01779                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01780 
01781   /* if the density criterion is satisfied there is nothing to do */
01782   if( density >= density_th ) return TRUE;
01783 
01784   /* compute region's radius */
01785   xc = (double) reg[0].x;
01786   yc = (double) reg[0].y;
01787   rad1 = dist( xc, yc, rec->x1, rec->y1 );
01788   rad2 = dist( xc, yc, rec->x2, rec->y2 );
01789   rad = rad1 > rad2 ? rad1 : rad2;
01790 
01791   /* while the density criterion is not satisfied, remove farther pixels */
01792   while( density < density_th )
01793     {
01794       rad *= 0.75; /* reduce region's radius to 75% of its value */
01795 
01796       /* remove points from the region and update 'used' map */
01797       for(i=0; i<*reg_size; i++)
01798         if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )
01799           {
01800             /* point not kept, mark it as NOTUSED */
01801             used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01802             /* remove point from the region */
01803             reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */
01804             reg[i].y = reg[*reg_size-1].y;
01805             --(*reg_size);
01806             --i; /* to avoid skipping one point */
01807           }
01808 
01809       /* reject if the region is too small.
01810          2 is the minimal region size for 'region2rect' to work. */
01811       if( *reg_size < 2 ) return FALSE;
01812 
01813       /* re-compute rectangle */
01814       region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01815 
01816       /* re-compute region points density */
01817       density = (double) *reg_size /
01818                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01819     }
01820 
01821   /* if this point is reached, the density criterion is satisfied */
01822   return TRUE;
01823 }
01824 
01825 /*----------------------------------------------------------------------------*/
01826 /** Refine a rectangle.
01827 
01828     For that, an estimation of the angle tolerance is performed by the
01829     standard deviation of the angle at points near the region's
01830     starting point. Then, a new region is grown starting from the same
01831     point, but using the estimated angle tolerance. If this fails to
01832     produce a rectangle with the right density of region points,
01833     'reduce_region_radius' is called to try to satisfy this condition.
01834  */
01835 static int refine( struct point * reg, int * reg_size, image_double modgrad,
01836                    double reg_angle, double prec, double p, struct rect * rec,
01837                    image_char used, image_double angles, double density_th )
01838 {
01839   double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;
01840   int i,n;
01841 
01842   /* check parameters */
01843   if( reg == NULL ) error("refine: invalid pointer 'reg'.");
01844   if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");
01845   if( prec < 0.0 ) error("refine: 'prec' must be positive.");
01846   if( rec == NULL ) error("refine: invalid pointer 'rec'.");
01847   if( used == NULL || used->data == NULL )
01848     error("refine: invalid image 'used'.");
01849   if( angles == NULL || angles->data == NULL )
01850     error("refine: invalid image 'angles'.");
01851 
01852   /* compute region points density */
01853   density = (double) *reg_size /
01854                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01855 
01856   /* if the density criterion is satisfied there is nothing to do */
01857   if( density >= density_th ) return TRUE;
01858 
01859   /*------ First try: reduce angle tolerance ------*/
01860 
01861   /* compute the new mean angle and tolerance */
01862   xc = (double) reg[0].x;
01863   yc = (double) reg[0].y;
01864   ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];
01865   sum = s_sum = 0.0;
01866   n = 0;
01867   for(i=0; i<*reg_size; i++)
01868     {
01869       used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01870       if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width )
01871         {
01872           angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];
01873           ang_d = angle_diff_signed(angle,ang_c);
01874           sum += ang_d;
01875           s_sum += ang_d * ang_d;
01876           ++n;
01877         }
01878     }
01879   mean_angle = sum / (double) n;
01880   tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n
01881                          + mean_angle*mean_angle ); /* 2 * standard deviation */
01882 
01883   /* find a new region from the same starting point and new angle tolerance */
01884   region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,&reg_angle,used,tau);
01885 
01886   /* if the region is too small, reject */
01887   if( *reg_size < 2 ) return FALSE;
01888 
01889   /* re-compute rectangle */
01890   region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01891 
01892   /* re-compute region points density */
01893   density = (double) *reg_size /
01894                       ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01895 
01896   /*------ Second try: reduce region radius ------*/
01897   if( density < density_th )
01898     return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,
01899                                  rec, used, angles, density_th );
01900 
01901   /* if this point is reached, the density criterion is satisfied */
01902   return TRUE;
01903 }
01904 
01905 
01906 /*----------------------------------------------------------------------------*/
01907 /*-------------------------- Line Segment Detector ---------------------------*/
01908 /*----------------------------------------------------------------------------*/
01909 
01910 /*----------------------------------------------------------------------------*/
01911 /** LSD full interface.
01912  */
01913 ntuple_list LineSegmentDetection( image_double image, double scale,
01914                                   double sigma_scale, double quant,
01915                                   double ang_th, double eps, double density_th,
01916                                   int n_bins, double max_grad,
01917                                   image_int * region )
01918 {
01919   ntuple_list out = new_ntuple_list(5);
01920   image_double scaled_image,angles,modgrad;
01921   image_char used;
01922   struct coorlist * list_p;
01923   void * mem_p;
01924   struct rect rec;
01925   struct point * reg;
01926   int reg_size,min_reg_size,i;
01927   unsigned int xsize,ysize;
01928   double rho,reg_angle,prec,p,log_nfa,logNT;
01929   int ls_count = 0;                   /* line segments are numbered 1,2,3,... */
01930 
01931 
01932   /* check parameters */
01933   if( image==NULL || image->data==NULL || image->xsize==0 || image->ysize==0 )
01934     error("invalid image input.");
01935   if( scale <= 0.0 ) error("'scale' value must be positive.");
01936   if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");
01937   if( quant < 0.0 ) error("'quant' value must be positive.");
01938   if( ang_th <= 0.0 || ang_th >= 180.0 )
01939     error("'ang_th' value must be in the range (0,180).");
01940   if( density_th < 0.0 || density_th > 1.0 )
01941     error("'density_th' value must be in the range [0,1].");
01942   if( n_bins <= 0 ) error("'n_bins' value must be positive.");
01943   if( max_grad <= 0.0 ) error("'max_grad' value must be positive.");
01944 
01945 
01946   /* angle tolerance */
01947   prec = M_PI * ang_th / 180.0;
01948   p = ang_th / 180.0;
01949   rho = quant / sin(prec); /* gradient magnitude threshold */
01950 
01951 
01952   /* scale image (if necessary) and compute angle at each pixel */
01953   if( scale != 1.0 )
01954     {
01955       scaled_image = gaussian_sampler( image, scale, sigma_scale );
01956       angles = ll_angle( scaled_image, rho, &list_p, &mem_p,
01957                          &modgrad, (unsigned int) n_bins, max_grad );
01958       free_image_double(scaled_image);
01959     }
01960   else
01961     angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,
01962                        (unsigned int) n_bins, max_grad );
01963   xsize = angles->xsize;
01964   ysize = angles->ysize;
01965   logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0;
01966   min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region
01967                                              that can give a meaningful event */
01968 
01969 
01970   /* initialize some structures */
01971   if( region != NULL ) /* image to output pixel region number, if asked */
01972     *region = new_image_int_ini(angles->xsize,angles->ysize,0);
01973   used = new_image_char_ini(xsize,ysize,NOTUSED);
01974   reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );
01975   if( reg == NULL ) error("not enough memory!");
01976 
01977 
01978   /* search for line segments */
01979   for(; list_p != NULL; list_p = list_p->next )
01980     if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&
01981         angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )
01982        /* there is no risk of double comparison problems here
01983           because we are only interested in the exact NOTDEF value */
01984       {
01985         /* find the region of connected point and ~equal angle */
01986         region_grow( list_p->x, list_p->y, angles, reg, &reg_size,
01987                      &reg_angle, used, prec );
01988 
01989         /* reject small regions */
01990         if( reg_size < min_reg_size ) continue;
01991 
01992         /* construct rectangular approximation for the region */
01993         region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);
01994 
01995         /* Check if the rectangle exceeds the minimal density of
01996            region points. If not, try to improve the region.
01997            The rectangle will be rejected if the final one does
01998            not fulfill the minimal density condition.
01999            This is an addition to the original LSD algorithm published in
02000            "LSD: A Fast Line Segment Detector with a False Detection Control"
02001            by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall.
02002            The original algorithm is obtained with density_th = 0.0.
02003          */
02004         if( !refine( reg, &reg_size, modgrad, reg_angle,
02005                      prec, p, &rec, used, angles, density_th ) ) continue;
02006 
02007         /* compute NFA value */
02008         log_nfa = rect_improve(&rec,angles,logNT,eps);
02009         if( log_nfa <= eps ) continue;
02010 
02011         /* A New Line Segment was found! */
02012         ++ls_count;  /* increase line segment counter */
02013 
02014         /*
02015            The gradient was computed with a 2x2 mask, its value corresponds to
02016            points with an offset of (0.5,0.5), that should be added to output.
02017            The coordinates origin is at the center of pixel (0,0).
02018          */
02019         rec.x1 += 0.5; rec.y1 += 0.5;
02020         rec.x2 += 0.5; rec.y2 += 0.5;
02021 
02022         /* scale the result values if a subsampling was performed */
02023         if( scale != 1.0 )
02024           {
02025             rec.x1 /= scale; rec.y1 /= scale;
02026             rec.x2 /= scale; rec.y2 /= scale;
02027             rec.width /= scale;
02028           }
02029 
02030         /* add line segment found to output */
02031         add_5tuple(out, rec.x1, rec.y1, rec.x2, rec.y2, rec.width);
02032 
02033         /* add region number to 'region' image if needed */
02034         if( region != NULL )
02035           for(i=0; i<reg_size; i++)
02036             (*region)->data[reg[i].x+reg[i].y*(*region)->xsize] = ls_count;
02037       }
02038 
02039 
02040   /* free memory */
02041   free_image_double(angles);
02042   free_image_double(modgrad);
02043   free_image_char(used);
02044   free( (void *) reg );
02045   free( (void *) mem_p );
02046 
02047   return out;
02048 }
02049 
02050 /*----------------------------------------------------------------------------*/
02051 /** LSD Simple Interface with Scale.
02052  */
02053 ntuple_list lsd_scale(image_double image, double scale)
02054 {
02055   /* LSD parameters */
02056   double sigma_scale = 0.6; /* Sigma for Gaussian filter is computed as
02057                                 sigma = sigma_scale/scale.                    */
02058   double quant = 2.0;       /* Bound to the quantization error on the
02059                                 gradient norm.                                */
02060   double ang_th = 22.5;     /* Gradient angle tolerance in degrees.           */
02061   double eps = 0.0;         /* Detection threshold, -log10(NFA).              */
02062   double density_th = 0.7;  /* Minimal density of region points in rectangle. */
02063   int n_bins = 1024;        /* Number of bins in pseudo-ordering of gradient
02064                                modulus.                                       */
02065   double max_grad = 255.0;  /* Gradient modulus in the highest bin. The
02066                                default value corresponds to the highest
02067                                gradient modulus on images with gray
02068                                levels in [0,255].                             */
02069 
02070   return LineSegmentDetection( image, scale, sigma_scale, quant, ang_th, eps,
02071                                density_th, n_bins, max_grad, NULL );
02072 }
02073 
02074 /*----------------------------------------------------------------------------*/
02075 /** LSD Simple Interface.
02076  */
02077 ntuple_list lsd(image_double image)
02078 {
02079   /* LSD parameters */
02080   double scale = 0.8;       /* Scale the image by Gaussian filter to 'scale'. */
02081 
02082   return lsd_scale(image,scale);
02083 }
02084 /*----------------------------------------------------------------------------*/

Generated on Fri Dec 3 10:18:18 2010 for LSD by doxygen 1.3.4