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

lsd.c File Reference


Detailed Description

LSD module code.

Author:
rafael grompone von gioi (grompone@gmail.com)

Definition in file lsd.c.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include "lsd.h"

Include dependency graph for lsd.c:

Include dependency graph

Go to the source code of this file.

Data Structures

struct  coorlist
 Chained list of coordinates. More...

struct  point
 A point (or pixel). More...

struct  rect
 Rectangle structure: line segment with width. More...

struct  rect_iter
 Rectangle points iterator. More...


Defines

#define M_LN10   2.30258509299404568402
 ln(10)

#define M_PI   3.14159265358979323846
 PI.

#define FALSE   0
#define TRUE   1
#define NOTDEF   -1024.0
 Label for pixels with undefined gradient.

#define M_3_2_PI   4.71238898038
 3/2 pi

#define M_2__PI   6.28318530718
 2 pi

#define NOTUSED   0
 Label for pixels not used in yet.

#define USED   1
 Label for pixels already used in detection.

#define RELATIVE_ERROR_FACTOR   100.0
 Doubles relative error factor.

#define log_gamma(x)   ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
 Computes the natural logarithm of the absolute value of the gamma function of x.

#define TABSIZE   100000
 Size of the table to store already computed inverse values.


Functions

void error (char *msg)
 Fatal error, print a message to standard-error output and exit.

int double_equal (double a, double b)
 Compare doubles by relative error.

double dist (double x1, double y1, double x2, double y2)
 Computes Euclidean distance between point (x1,y1) and point (x2,y2).

void free_ntuple_list (ntuple_list in)
 Free memory used in n-tuple 'in'.

ntuple_list new_ntuple_list (unsigned int dim)
 Create an n-tuple list and allocate memory for one element.

void enlarge_ntuple_list (ntuple_list n_tuple)
 Enlarge the allocated memory of an n-tuple list.

void add_5tuple (ntuple_list out, double v1, double v2, double v3, double v4, double v5)
 Add a 5-tuple to an n-tuple list.

void free_image_char (image_char i)
 Free memory used in image_char 'i'.

image_char new_image_char (unsigned int xsize, unsigned int ysize)
 Create a new image_char of size 'xsize' times 'ysize'.

image_char new_image_char_ini (unsigned int xsize, unsigned int ysize, unsigned char fill_value)
 Create a new image_char of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

void free_image_int (image_int i)
 Free memory used in image_int 'i'.

image_int new_image_int (unsigned int xsize, unsigned int ysize)
 Create a new image_int of size 'xsize' times 'ysize'.

image_int new_image_int_ini (unsigned int xsize, unsigned int ysize, int fill_value)
 Create a new image_int of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

void free_image_double (image_double i)
 Free memory used in image_double 'i'.

image_double new_image_double (unsigned int xsize, unsigned int ysize)
 Create a new image_double of size 'xsize' times 'ysize'.

image_double new_image_double_ini (unsigned int xsize, unsigned int ysize, double fill_value)
 Create a new image_double of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

void gaussian_kernel (ntuple_list kernel, double sigma, double mean)
 Compute a Gaussian kernel of length 'kernel->dim', standard deviation 'sigma', and centered at value 'mean'.

image_double gaussian_sampler (image_double in, double scale, double sigma_scale)
 Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.

image_double ll_angle (image_double in, double threshold, struct coorlist **list_p, void **mem_p, image_double *modgrad, unsigned int n_bins, double max_grad)
 Computes the direction of the level line of 'in' at each point.

int isaligned (int x, int y, image_double angles, double theta, double prec)
 Is point (x,y) aligned to angle theta, up to precision 'prec'?

double angle_diff (double a, double b)
 Absolute value angle difference.

double angle_diff_signed (double a, double b)
 Signed angle difference.

double log_gamma_lanczos (double x)
 Computes the natural logarithm of the absolute value of the gamma function of x using the Lanczos approximation.

double log_gamma_windschitl (double x)
 Computes the natural logarithm of the absolute value of the gamma function of x using Windschitl method.

double nfa (int n, int k, double p, double logNT)
 Computes -log10(NFA).

void rect_copy (struct rect *in, struct rect *out)
 Copy one rectangle structure to another.

double inter_low (double x, double x1, double y1, double x2, double y2)
 Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller of 'y1' and 'y2'.

double inter_hi (double x, double x1, double y1, double x2, double y2)
 Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger of 'y1' and 'y2'.

void ri_del (rect_iter *iter)
 Free memory used by a rectangle iterator.

int ri_end (rect_iter *i)
 Check if the iterator finished the full iteration.

void ri_inc (rect_iter *i)
 Increment a rectangle iterator.

rect_iterri_ini (struct rect *r)
 Create and initialize a rectangle iterator.

double rect_nfa (struct rect *rec, image_double angles, double logNT)
 Compute a rectangle's NFA value.

double get_theta (struct point *reg, int reg_size, double x, double y, image_double modgrad, double reg_angle, double prec)
 Compute region's angle as the principal inertia axis of the region.

void region2rect (struct point *reg, int reg_size, image_double modgrad, double reg_angle, double prec, double p, struct rect *rec)
 Computes a rectangle that covers a region of points.

void region_grow (int x, int y, image_double angles, struct point *reg, int *reg_size, double *reg_angle, image_char used, double prec)
 Build a region of pixels that share the same angle, up to a tolerance 'prec', starting at point (x,y).

double rect_improve (struct rect *rec, image_double angles, double logNT, double eps)
 Try some rectangles variations to improve NFA value.

int reduce_region_radius (struct point *reg, int *reg_size, image_double modgrad, double reg_angle, double prec, double p, struct rect *rec, image_char used, image_double angles, double density_th)
 Reduce the region size, by elimination the points far from the starting point, until that leads to rectangle with the right density of region points or to discard the region if too small.

int refine (struct point *reg, int *reg_size, image_double modgrad, double reg_angle, double prec, double p, struct rect *rec, image_char used, image_double angles, double density_th)
 Refine a rectangle.

ntuple_list LineSegmentDetection (image_double image, double scale, double sigma_scale, double quant, double ang_th, double eps, double density_th, int n_bins, double max_grad, image_int *region)
 LSD Full Interface.

ntuple_list lsd_scale (image_double image, double scale)
 LSD Simple Interface with Scale.

ntuple_list lsd (image_double image)
 LSD Simple Interface.


Define Documentation

#define FALSE   0
 

Definition at line 89 of file lsd.c.

Referenced by isaligned(), reduce_region_radius(), and refine().

#define log_gamma  )     ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
 

Computes the natural logarithm of the absolute value of the gamma function of x.

When x>15 use log_gamma_windschitl(), otherwise use log_gamma_lanczos().

Definition at line 913 of file lsd.c.

Referenced by nfa().

#define M_2__PI   6.28318530718
 

2 pi

Definition at line 103 of file lsd.c.

Referenced by angle_diff(), angle_diff_signed(), and isaligned().

#define M_3_2_PI   4.71238898038
 

3/2 pi

Definition at line 100 of file lsd.c.

Referenced by isaligned().

#define M_LN10   2.30258509299404568402
 

ln(10)

Definition at line 80 of file lsd.c.

Referenced by nfa().

#define M_PI   3.14159265358979323846
 

PI.

Definition at line 85 of file lsd.c.

Referenced by angle_diff(), angle_diff_signed(), get_theta(), LineSegmentDetection(), and rect_improve().

#define NOTDEF   -1024.0
 

Label for pixels with undefined gradient.

Definition at line 97 of file lsd.c.

Referenced by isaligned(), LineSegmentDetection(), and ll_angle().

#define NOTUSED   0
 

Label for pixels not used in yet.

Definition at line 106 of file lsd.c.

Referenced by LineSegmentDetection(), reduce_region_radius(), and refine().

#define RELATIVE_ERROR_FACTOR   100.0
 

Doubles relative error factor.

Definition at line 142 of file lsd.c.

Referenced by double_equal().

#define TABSIZE   100000
 

Size of the table to store already computed inverse values.

Definition at line 918 of file lsd.c.

Referenced by nfa().

#define TRUE   1
 

Definition at line 93 of file lsd.c.

Referenced by double_equal(), reduce_region_radius(), and refine().

#define USED   1
 

Label for pixels already used in detection.

Definition at line 109 of file lsd.c.

Referenced by region_grow().


Function Documentation

void add_5tuple ntuple_list  out,
double  v1,
double  v2,
double  v3,
double  v4,
double  v5
[static]
 

Add a 5-tuple to an n-tuple list.

Definition at line 250 of file lsd.c.

References ntuple_list_s::dim, enlarge_ntuple_list(), error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

double angle_diff double  a,
double  b
[static]
 

Absolute value angle difference.

Definition at line 819 of file lsd.c.

References M_2__PI, and M_PI.

Referenced by get_theta().

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 }

double angle_diff_signed double  a,
double  b
[static]
 

Signed angle difference.

Definition at line 831 of file lsd.c.

References M_2__PI, and M_PI.

Referenced by refine().

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 }

double dist double  x1,
double  y1,
double  x2,
double  y2
[static]
 

Computes Euclidean distance between point (x1,y1) and point (x2,y2).

Definition at line 181 of file lsd.c.

Referenced by reduce_region_radius(), and refine().

00182 {
00183   return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00184 }

int double_equal double  a,
double  b
[static]
 

Compare doubles by relative error.

The resulting rounding error after floating point computations depend on the specific operations done. The same number computed by different algorithms could present different rounding errors. For a useful comparison, an estimation of the relative rounding error should be considered and compared to a factor times EPS. The factor should be related to the cumulated rounding error in the chain of computation. Here, as a simplification, a fixed factor is used.

Definition at line 155 of file lsd.c.

References RELATIVE_ERROR_FACTOR, and TRUE.

Referenced by get_theta(), inter_hi(), inter_low(), and nfa().

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 }

void enlarge_ntuple_list ntuple_list  n_tuple  )  [static]
 

Enlarge the allocated memory of an n-tuple list.

Definition at line 232 of file lsd.c.

References ntuple_list_s::dim, error(), ntuple_list_s::max_size, ntuple_list, and ntuple_list_s::values.

Referenced by add_5tuple(), and gaussian_kernel().

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 }

Here is the call graph for this function:

void error char *  msg  )  [static]
 

Fatal error, print a message to standard-error output and exit.

Definition at line 133 of file lsd.c.

Referenced by add_5tuple(), enlarge_ntuple_list(), free_image_char(), free_image_double(), free_image_int(), free_ntuple_list(), gaussian_kernel(), gaussian_sampler(), get_theta(), inter_hi(), inter_low(), isaligned(), LineSegmentDetection(), ll_angle(), new_image_char(), new_image_char_ini(), new_image_double(), new_image_int(), new_ntuple_list(), nfa(), rect_copy(), rect_nfa(), reduce_region_radius(), refine(), region2rect(), region_grow(), ri_del(), ri_end(), ri_inc(), and ri_ini().

00134 {
00135   fprintf(stderr,"LSD Error: %s\n",msg);
00136   exit(EXIT_FAILURE);
00137 }

void free_image_char image_char  i  ) 
 

Free memory used in image_char 'i'.

Definition at line 280 of file lsd.c.

References image_char_s::data, error(), and image_char.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

void free_image_double image_double  i  ) 
 

Free memory used in image_double 'i'.

Definition at line 387 of file lsd.c.

References image_double_s::data, error(), and image_double.

Referenced by gaussian_sampler(), and LineSegmentDetection().

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 }

Here is the call graph for this function:

void free_image_int image_int  i  ) 
 

Free memory used in image_int 'i'.

Definition at line 336 of file lsd.c.

References image_int_s::data, error(), and image_int.

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 }

Here is the call graph for this function:

void free_ntuple_list ntuple_list  in  ) 
 

Free memory used in n-tuple 'in'.

Definition at line 194 of file lsd.c.

References error(), ntuple_list, and ntuple_list_s::values.

Referenced by gaussian_sampler().

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 }

Here is the call graph for this function:

void gaussian_kernel ntuple_list  kernel,
double  sigma,
double  mean
[static]
 

Compute a Gaussian kernel of length 'kernel->dim', standard deviation 'sigma', and centered at value 'mean'.

For example, if mean=0.5, the Gaussian will be centered in the middle point between values 'kernel->values[0]' and 'kernel->values[1]'.

Definition at line 448 of file lsd.c.

References ntuple_list_s::dim, enlarge_ntuple_list(), error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values.

Referenced by gaussian_sampler().

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 }

Here is the call graph for this function:

image_double gaussian_sampler image_double  in,
double  scale,
double  sigma_scale
[static]
 

Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.

For example, scale=0.8 will give a result at 80% of the original size.

The image is convolved with a Gaussian kernel

\[ G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} \]

before the sub-sampling to prevent aliasing.

The standard deviation sigma given by:

  • sigma = sigma_scale / scale, if scale < 1.0
  • sigma = sigma_scale, if scale >= 1.0

To be able to sub-sample at non-integer steps, some interpolation is needed. In this implementation, the interpolation is done by the Gaussian kernel, so both operations (filtering and sampling) are done at the same time. The Gaussian kernel is computed centered on the coordinates of the required sample. In this way, when applied, it gives directly the result of convolving the image with the kernel and interpolated to that particular position.

A fast algorithm is done using the separability of the Gaussian kernel. Applying the 2D Gaussian kernel is equivalent to applying first a horizontal 1D Gaussian kernel and then a vertical 1D Gaussian kernel (or the other way round). The reason is that

\[ G(x,y) = G(x) * G(y) \]

where

\[ G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}}. \]

The algorithm first apply a combined Gaussian kernel and sampling in the x axis, and then the combined Gaussian kernel and sampling in the y axis.

Definition at line 511 of file lsd.c.

References image_double_s::data, ntuple_list_s::dim, error(), free_image_double(), free_ntuple_list(), gaussian_kernel(), image_double, new_image_double(), new_ntuple_list(), ntuple_list, ntuple_list_s::values, image_double_s::xsize, and image_double_s::ysize.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

double get_theta struct point reg,
int  reg_size,
double  x,
double  y,
image_double  modgrad,
double  reg_angle,
double  prec
[static]
 

Compute region's angle as the principal inertia axis of the region.

The following is the region inertia matrix A:

\[ A = \left(\begin{array}{cc} Ixx & Ixy \\ Ixy & Iyy \\ \end{array}\right) \]

where

Ixx = sum_i G(i).(y_i - cx)^2

Iyy = sum_i G(i).(x_i - cy)^2

Ixy = - sum_i G(i).(x_i - cx).(y_i - cy)

and

  • G(i) is the gradient norm at pixel i, used as pixel's weight.
  • x_i and y_i are the coordinates of pixel i.
  • cx and cy are the coordinates of the center of th region.

lambda1 and lambda2 are the eigenvalues of matrix A, with lambda1 >= lambda2. They are found by solving the characteristic polynomial:

det( lambda I - A) = 0

that gives:

lambda1 = ( Ixx + Iyy + sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2

lambda2 = ( Ixx + Iyy - sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2

To get the line segment direction we want to get the angle the eigenvector assotiated to the smaller eigenvalue. We have to solve a,b in:

a.Ixx + b.Ixy = a.lambda2

a.Ixy + b.Iyy = b.lambda2

We want the angle theta = atan(b/a). It can be computed with any of the two equations:

theta = atan( (lambda2-Ixx) / Ixy )

or

theta = atan( Ixy / (lambda2-Iyy) )

When |Ixx| > |Iyy| we use the first, otherwise the second (just to get better numeric precision).

Definition at line 1456 of file lsd.c.

References angle_diff(), image_double_s::data, double_equal(), error(), image_double, M_PI, point::x, image_double_s::xsize, and point::y.

Referenced by region2rect().

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 }

Here is the call graph for this function:

double inter_hi double  x,
double  x1,
double  y1,
double  x2,
double  y2
[static]
 

Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger of 'y1' and 'y2'.

The following restrictions are required:

  • x1 <= x2
  • x1 <= x
  • x <= x2

Definition at line 1187 of file lsd.c.

References double_equal(), and error().

Referenced by ri_inc().

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 }

Here is the call graph for this function:

double inter_low double  x,
double  x1,
double  y1,
double  x2,
double  y2
[static]
 

Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller of 'y1' and 'y2'.

The following restrictions are required:

  • x1 <= x2
  • x1 <= x
  • x <= x2

Definition at line 1165 of file lsd.c.

References double_equal(), and error().

Referenced by ri_inc().

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 }

Here is the call graph for this function:

int isaligned int  x,
int  y,
image_double  angles,
double  theta,
double  prec
[static]
 

Is point (x,y) aligned to angle theta, up to precision 'prec'?

Definition at line 781 of file lsd.c.

References image_double_s::data, error(), FALSE, image_double, M_2__PI, M_3_2_PI, NOTDEF, image_double_s::xsize, and image_double_s::ysize.

Referenced by rect_nfa(), and region_grow().

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 }

Here is the call graph for this function:

ntuple_list LineSegmentDetection image_double  image,
double  scale,
double  sigma_scale,
double  quant,
double  ang_th,
double  eps,
double  density_th,
int  n_bins,
double  max_grad,
image_int region
 

LSD Full Interface.

Parameters:
image Input image.
scale When different than 1.0, LSD will scale the image by Gaussian filtering. Example: is scale=0.8, the input image will be subsampled to 80% of its size, and then the line segment detector will be applied. Suggested value: 0.8
sigma_scale When scale!=1.0, the sigma of the Gaussian filter is: sigma = sigma_scale / scale, if scale < 1.0 sigma = sigma_scale, if scale >= 1.0 Suggested value: 0.6
quant Bound to the quantization error on the gradient norm. Example: if gray level is quantized to integer steps, the gradient (computed by finite differences) error due to quantization will be bounded by 2.0, as the worst case is when the error are 1 and -1, that gives an error of 2.0. Suggested value: 2.0
ang_th Gradient angle tolerance in the region growing algorithm, in degrees. Suggested value: 22.5
eps Detection threshold, -log10(NFA). The bigger, the more strict the detector is, and will result in less detections. (Note that the 'minus sign' makes that this behavior is opposite to the one of NFA.) The value -log10(NFA) is equivalent but more intuitive than NFA:
  • -1.0 corresponds to 10 mean false alarms
  • 0.0 corresponds to 1 mean false alarm
  • 1.0 corresponds to 0.1 mean false alarms
  • 2.0 corresponds to 0.01 mean false alarms
Suggested value: 0.0
density_th Minimal proportion of region points in a rectangle. Suggested value: 0.7
n_bins Number of bins used in the pseudo-ordering of gradient modulus. Suggested value: 1024
max_grad Gradient modulus in the highest bin. For example, for images with integer gray levels in [0,255], the maximum possible gradient value is 255.0. Suggested value: 255.0
region Optional output: an int image where the pixels used in some line support region are marked. Unused pixels have the value '0' while the used ones have the number of the line segment, numbered 1,2,3,... If desired, a non NULL pointer to an image_int should be used. The resulting image has the size of the image used for the processing, that is, the size of the input image scaled by the given factor 'scale'. Suggested value: NULL
Returns:
A 5-tuple list, where each 5-tuple corresponds to a detected line segment. The five values are:
  • x1,y1,x2,y2,width
for a line segment from (x1,y1) to (x2,y2) and a width 'width'.

Definition at line 1913 of file lsd.c.

References add_5tuple(), image_char_s::data, image_double_s::data, error(), free_image_char(), free_image_double(), gaussian_sampler(), image_char, image_double, image_int, ll_angle(), M_PI, new_image_char_ini(), new_image_int_ini(), new_ntuple_list(), NOTDEF, NOTUSED, ntuple_list, rect_improve(), refine(), region2rect(), region_grow(), rect::width, rect::x1, rect::x2, image_char_s::xsize, image_double_s::xsize, rect::y1, rect::y2, and image_double_s::ysize.

Referenced by lsd_scale().

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 }

Here is the call graph for this function:

image_double ll_angle image_double  in,
double  threshold,
struct coorlist **  list_p,
void **  mem_p,
image_double modgrad,
unsigned int  n_bins,
double  max_grad
[static]
 

Computes the direction of the level line of 'in' at each point.

The result is:

  • an image_double with the angle at each pixel, or NOTDEF if not defined.
  • the image_double 'modgrad' (a pointer is passed as argument) with the gradient magnitude at each point.
  • a list of pixels 'list_p' roughly ordered by decreasing gradient magnitude. (The order is made by classifying points into bins by gradient magnitude. The parameters 'n_bins' and 'max_grad' specify the number of bins and the gradient modulus at the highest bin. The pixels in the list would be in decreasing gradient magnitude, up to a precision of the size of the bins.)
  • a pointer 'mem_p' to the memory used by 'list_p' to be able to free the memory when it is not used anymore.

Definition at line 652 of file lsd.c.

References image_double_s::data, error(), image_double, new_image_double(), NOTDEF, image_double_s::xsize, and image_double_s::ysize.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

double log_gamma_lanczos double  x  )  [static]
 

Computes the natural logarithm of the absolute value of the gamma function of x using the Lanczos approximation.

See http://www.rskey.org/gamma.htm

The formula used is

\[ \Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) } (x+5.5)^{x+0.5} e^{-(x+5.5)} \]

so

\[ \log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right) + (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n) \]

and q0 = 75122.6331530, q1 = 80916.6278952, q2 = 36308.2951477, q3 = 8687.24529705, q4 = 1168.92649479, q5 = 83.8676043424, q6 = 2.50662827511.

Definition at line 868 of file lsd.c.

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 }

double log_gamma_windschitl double  x  )  [static]
 

Computes the natural logarithm of the absolute value of the gamma function of x using Windschitl method.

See http://www.rskey.org/gamma.htm

The formula used is

\[ \Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e} \sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x \]

so

\[ \log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x + 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right). \]

This formula is a good approximation when x > 15.

Definition at line 902 of file lsd.c.

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 }

ntuple_list lsd image_double  image  ) 
 

LSD Simple Interface.

Parameters:
image Input image.
Returns:
a 5-tuple list of detected line segments.

Definition at line 2077 of file lsd.c.

References image_double, lsd_scale(), and ntuple_list.

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 }

Here is the call graph for this function:

ntuple_list lsd_scale image_double  image,
double  scale
 

LSD Simple Interface with Scale.

Parameters:
image Input image.
scale When different than 1.0, LSD will scale the image by Gaussian filtering. Example: is scale=0.8, the input image will be subsampled to 80% of its size, and then the line segment detector will be applied. Suggested value: 0.8
Returns:
a 5-tuple list of detected line segments.

Definition at line 2053 of file lsd.c.

References image_double, LineSegmentDetection(), and ntuple_list.

Referenced by lsd().

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 }

Here is the call graph for this function:

image_char new_image_char unsigned int  xsize,
unsigned int  ysize
 

Create a new image_char of size 'xsize' times 'ysize'.

Definition at line 291 of file lsd.c.

References image_char_s::data, error(), image_char, image_char_s::xsize, and image_char_s::ysize.

Referenced by new_image_char_ini().

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 }

Here is the call graph for this function:

image_char new_image_char_ini unsigned int  xsize,
unsigned int  ysize,
unsigned char  fill_value
 

Create a new image_char of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

Definition at line 316 of file lsd.c.

References image_char_s::data, error(), image_char, and new_image_char().

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

image_double new_image_double unsigned int  xsize,
unsigned int  ysize
 

Create a new image_double of size 'xsize' times 'ysize'.

Definition at line 398 of file lsd.c.

References image_double_s::data, error(), image_double, image_double_s::xsize, and image_double_s::ysize.

Referenced by gaussian_sampler(), ll_angle(), and new_image_double_ini().

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 }

Here is the call graph for this function:

image_double new_image_double_ini unsigned int  xsize,
unsigned int  ysize,
double  fill_value
 

Create a new image_double of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

Definition at line 422 of file lsd.c.

References image_double_s::data, image_double, and new_image_double().

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 }

Here is the call graph for this function:

image_int new_image_int unsigned int  xsize,
unsigned int  ysize
 

Create a new image_int of size 'xsize' times 'ysize'.

Definition at line 347 of file lsd.c.

References image_int_s::data, error(), image_int, image_int_s::xsize, and image_int_s::ysize.

Referenced by new_image_int_ini().

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 }

Here is the call graph for this function:

image_int new_image_int_ini unsigned int  xsize,
unsigned int  ysize,
int  fill_value
 

Create a new image_int of size 'xsize' times 'ysize', initialized to the value 'fill_value'.

Definition at line 371 of file lsd.c.

References image_int_s::data, image_int, and new_image_int().

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

ntuple_list new_ntuple_list unsigned int  dim  ) 
 

Create an n-tuple list and allocate memory for one element.

Parameters:
dim the dimension (n) of the n-tuple.

Definition at line 206 of file lsd.c.

References ntuple_list_s::dim, error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values.

Referenced by gaussian_sampler(), and LineSegmentDetection().

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 }

Here is the call graph for this function:

double nfa int  n,
int  k,
double  p,
double  logNT
[static]
 

Computes -log10(NFA).

NFA stands for Number of False Alarms:

\[ \mathrm{NFA} = NT \cdot B(n,k,p) \]

  • NT - number of tests
  • B(n,k,p) - tail of binomial distribution with parameters n,k and p:

    \[ B(n,k,p) = \sum_{j=k}^n \left(\begin{array}{c}n\\j\end{array}\right) p^{j} (1-p)^{n-j} \]

The value -log10(NFA) is equivalent but more intuitive than NFA:

  • -1 corresponds to 10 mean false alarms
  • 0 corresponds to 1 mean false alarm
  • 1 corresponds to 0.1 mean false alarms
  • 2 corresponds to 0.01 mean false alarms
  • ...

Used this way, the bigger the value, better the detection, and a logarithmic scale is used.

Parameters:
n,k,p binomial parameters.
logNT logarithm of Number of Tests
The computation is based in the gamma function by the following relation:

\[ \left(\begin{array}{c}n\\k\end{array}\right) = \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }. \]

We use efficient algorithms to compute the logarithm of the gamma function.

To make the computation faster, not all the sum is computed, part of the terms are neglected based on a bound to the error obtained (an error of 10% in the result is accepted).

Definition at line 962 of file lsd.c.

References double_equal(), error(), log_gamma, M_LN10, and TABSIZE.

Referenced by rect_nfa().

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 }

Here is the call graph for this function:

void rect_copy struct rect in,
struct rect out
[static]
 

Copy one rectangle structure to another.

Definition at line 1071 of file lsd.c.

References rect::dx, rect::dy, error(), rect::p, rect::prec, rect::theta, rect::width, rect::x, rect::x1, rect::x2, rect::y, rect::y1, and rect::y2.

Referenced by rect_improve().

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 }

Here is the call graph for this function:

double rect_improve struct rect rec,
image_double  angles,
double  logNT,
double  eps
[static]
 

Try some rectangles variations to improve NFA value.

Only if the rectangle is not meaningful (i.e., log_nfa <= eps).

Definition at line 1644 of file lsd.c.

References rect::dx, rect::dy, image_double, M_PI, rect::p, rect::prec, rect_copy(), rect_nfa(), rect::width, rect::x1, rect::x2, rect::y1, and rect::y2.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

double rect_nfa struct rect rec,
image_double  angles,
double  logNT
[static]
 

Compute a rectangle's NFA value.

Definition at line 1370 of file lsd.c.

References error(), image_double, isaligned(), nfa(), rect::p, rect::prec, ri_del(), ri_end(), ri_inc(), ri_ini(), rect::theta, rect_iter::x, image_double_s::xsize, rect_iter::y, and image_double_s::ysize.

Referenced by rect_improve().

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 }

Here is the call graph for this function:

int reduce_region_radius struct point reg,
int *  reg_size,
image_double  modgrad,
double  reg_angle,
double  prec,
double  p,
struct rect rec,
image_char  used,
image_double  angles,
double  density_th
[static]
 

Reduce the region size, by elimination the points far from the starting point, until that leads to rectangle with the right density of region points or to discard the region if too small.

Definition at line 1757 of file lsd.c.

References image_double_s::data, image_char_s::data, dist(), error(), FALSE, image_char, image_double, NOTUSED, region2rect(), TRUE, rect::width, point::x, rect::x1, rect::x2, image_char_s::xsize, point::y, rect::y1, and rect::y2.

Referenced by refine().

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 }

Here is the call graph for this function:

int refine struct point reg,
int *  reg_size,
image_double  modgrad,
double  reg_angle,
double  prec,
double  p,
struct rect rec,
image_char  used,
image_double  angles,
double  density_th
[static]
 

Refine a rectangle.

For that, an estimation of the angle tolerance is performed by the standard deviation of the angle at points near the region's starting point. Then, a new region is grown starting from the same point, but using the estimated angle tolerance. If this fails to produce a rectangle with the right density of region points, 'reduce_region_radius' is called to try to satisfy this condition.

Definition at line 1835 of file lsd.c.

References angle_diff_signed(), image_double_s::data, image_char_s::data, dist(), error(), FALSE, image_char, image_double, NOTUSED, reduce_region_radius(), region2rect(), region_grow(), TRUE, rect::width, point::x, rect::x1, rect::x2, image_char_s::xsize, image_double_s::xsize, point::y, rect::y1, and rect::y2.

Referenced by LineSegmentDetection().

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 }

Here is the call graph for this function:

void region2rect struct point reg,
int  reg_size,
image_double  modgrad,
double  reg_angle,
double  prec,
double  p,
struct rect rec
[static]
 

Computes a rectangle that covers a region of points.

Definition at line 1499 of file lsd.c.

References image_double_s::data, rect::dx, rect::dy, error(), get_theta(), image_double, rect::p, rect::prec, rect::theta, rect::width, rect::x, point::x, rect::x1, rect::x2, image_double_s::xsize, rect::y, point::y, rect::y1, and rect::y2.

Referenced by LineSegmentDetection(), reduce_region_radius(), and refine().

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 }

Here is the call graph for this function:

void region_grow int  x,
int  y,
image_double  angles,
struct point reg,
int *  reg_size,
double *  reg_angle,
image_char  used,
double  prec
[static]
 

Build a region of pixels that share the same angle, up to a tolerance 'prec', starting at point (x,y).

Definition at line 1592 of file lsd.c.

References image_char_s::data, image_double_s::data, error(), image_char, image_double, isaligned(), USED, point::x, image_char_s::xsize, image_double_s::xsize, point::y, image_char_s::ysize, and image_double_s::ysize.

Referenced by LineSegmentDetection(), and refine().

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 }

Here is the call graph for this function:

void ri_del rect_iter iter  )  [static]
 

Free memory used by a rectangle iterator.

Definition at line 1202 of file lsd.c.

References error().

Referenced by rect_nfa().

01203 {
01204   if( iter == NULL ) error("ri_del: NULL iterator.");
01205   free( (void *) iter );
01206 }

Here is the call graph for this function:

int ri_end rect_iter i  )  [static]
 

Check if the iterator finished the full iteration.

See details in rect_iter

Definition at line 1213 of file lsd.c.

References error(), rect_iter::vx, and rect_iter::x.

Referenced by rect_nfa(), and ri_inc().

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 }

Here is the call graph for this function:

void ri_inc rect_iter i  )  [static]
 

Increment a rectangle iterator.

See details in rect_iter

Definition at line 1229 of file lsd.c.

References error(), inter_hi(), inter_low(), ri_end(), rect_iter::vx, rect_iter::vy, rect_iter::x, rect_iter::y, rect_iter::ye, and rect_iter::ys.

Referenced by rect_nfa(), and ri_ini().

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 }

Here is the call graph for this function:

rect_iter* ri_ini struct rect r  )  [static]
 

Create and initialize a rectangle iterator.

See details in rect_iter

Definition at line 1299 of file lsd.c.

References rect::dx, rect::dy, error(), ri_inc(), rect_iter::vx, rect_iter::vy, rect::width, rect_iter::x, rect::x1, rect::x2, rect_iter::y, rect::y1, rect::y2, rect_iter::ye, and rect_iter::ys.

Referenced by rect_nfa().

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 }

Here is the call graph for this function:


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