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  image_char_s
 char image data type More...

struct  image_double_s
 double image data type More...

struct  image_int_s
 int image data type More...

struct  ntuple_list_s
 'list of n-tuple' data type 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.


Typedefs

typedef ntuple_list_sntuple_list
 'list of n-tuple' data type

typedef image_char_simage_char
 char image data type

typedef image_int_simage_int
 int image data type

typedef image_double_simage_double
 double image data type


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_7tuple (ntuple_list out, double v1, double v2, double v3, double v4, double v5, double v6, double v7)
 Add a 7-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'.

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_ptr (unsigned int xsize, unsigned int ysize, double *data)
 Create a new image_double of size 'xsize' times 'ysize' with the data pointed by 'data'.

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)
 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 log_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.

double * LineSegmentDetection (int *n_out, double *img, int X, int Y, double scale, double sigma_scale, double quant, double ang_th, double log_eps, double density_th, int n_bins, int **reg_img, int *reg_x, int *reg_y)
 LSD Full Interface.

double * lsd_scale_region (int *n_out, double *img, int X, int Y, double scale, int **reg_img, int *reg_x, int *reg_y)
 LSD Simple Interface with Scale and Region output.

double * lsd_scale (int *n_out, double *img, int X, int Y, double scale)
 LSD Simple Interface with Scale.

double * lsd (int *n_out, double *img, int X, int Y)
 LSD Simple Interface.


Define Documentation

#define FALSE   0
 

Definition at line 115 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 1025 of file lsd.c.

Referenced by nfa().

#define M_2__PI   6.28318530718
 

2 pi

Definition at line 129 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 126 of file lsd.c.

Referenced by isaligned().

#define M_LN10   2.30258509299404568402
 

ln(10)

Definition at line 106 of file lsd.c.

Referenced by nfa().

#define M_PI   3.14159265358979323846
 

PI.

Definition at line 111 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 123 of file lsd.c.

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

#define NOTUSED   0
 

Label for pixels not used in yet.

Definition at line 132 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 168 of file lsd.c.

Referenced by double_equal().

#define TABSIZE   100000
 

Size of the table to store already computed inverse values.

Definition at line 1030 of file lsd.c.

Referenced by nfa().

#define TRUE   1
 

Definition at line 119 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 135 of file lsd.c.

Referenced by region_grow().


Typedef Documentation

typedef struct image_char_s * image_char
 

char image data type

The pixel value at (x,y) is accessed by:

image->data[ x + y * image->xsize ]

with x and y integer.

Referenced by free_image_char(), LineSegmentDetection(), new_image_char(), new_image_char_ini(), reduce_region_radius(), refine(), and region_grow().

typedef struct image_double_s * image_double
 

double image data type

The pixel value at (x,y) is accessed by:

image->data[ x + y * image->xsize ]

with x and y integer.

Referenced by free_image_double(), gaussian_sampler(), get_theta(), isaligned(), LineSegmentDetection(), ll_angle(), new_image_double(), new_image_double_ptr(), rect_improve(), rect_nfa(), reduce_region_radius(), refine(), region2rect(), and region_grow().

typedef struct image_int_s * image_int
 

int image data type

The pixel value at (x,y) is accessed by:

image->data[ x + y * image->xsize ]

with x and y integer.

Referenced by LineSegmentDetection(), new_image_int(), and new_image_int_ini().

typedef struct ntuple_list_s * ntuple_list
 

'list of n-tuple' data type

The i-th component of the j-th n-tuple of an n-tuple list 'ntl' is accessed with:

ntl->values[ i + j * ntl->dim ]

The dimension of the n-tuple (n) is:

ntl->dim

The number of n-tuples in the list is:

ntl->size

The maximum number of n-tuples that can be stored in the list with the allocated memory at a given time is given by:

ntl->max_size

Referenced by add_7tuple(), enlarge_ntuple_list(), free_ntuple_list(), gaussian_kernel(), gaussian_sampler(), LineSegmentDetection(), and new_ntuple_list().


Function Documentation

void add_7tuple ntuple_list  out,
double  v1,
double  v2,
double  v3,
double  v4,
double  v5,
double  v6,
double  v7
[static]
 

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

Definition at line 305 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().

00307 {
00308   /* check parameters */
00309   if( out == NULL ) error("add_7tuple: invalid n-tuple input.");
00310   if( out->dim != 7 ) error("add_7tuple: the n-tuple must be a 7-tuple.");
00311 
00312   /* if needed, alloc more tuples to 'out' */
00313   if( out->size == out->max_size ) enlarge_ntuple_list(out);
00314   if( out->values == NULL ) error("add_7tuple: invalid n-tuple input.");
00315 
00316   /* add new 7-tuple */
00317   out->values[ out->size * out->dim + 0 ] = v1;
00318   out->values[ out->size * out->dim + 1 ] = v2;
00319   out->values[ out->size * out->dim + 2 ] = v3;
00320   out->values[ out->size * out->dim + 3 ] = v4;
00321   out->values[ out->size * out->dim + 4 ] = v5;
00322   out->values[ out->size * out->dim + 5 ] = v6;
00323   out->values[ out->size * out->dim + 6 ] = v7;
00324 
00325   /* update number of tuples counter */
00326   out->size++;
00327 }

Here is the call graph for this function:

double angle_diff double  a,
double  b
[static]
 

Absolute value angle difference.

Definition at line 931 of file lsd.c.

References M_2__PI, and M_PI.

Referenced by get_theta().

00932 {
00933   a -= b;
00934   while( a <= -M_PI ) a += M_2__PI;
00935   while( a >   M_PI ) a -= M_2__PI;
00936   if( a < 0.0 ) a = -a;
00937   return a;
00938 }

double angle_diff_signed double  a,
double  b
[static]
 

Signed angle difference.

Definition at line 943 of file lsd.c.

References M_2__PI, and M_PI.

Referenced by refine().

00944 {
00945   a -= b;
00946   while( a <= -M_PI ) a += M_2__PI;
00947   while( a >   M_PI ) a -= M_2__PI;
00948   return a;
00949 }

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 207 of file lsd.c.

Referenced by reduce_region_radius(), and refine().

00208 {
00209   return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00210 }

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 181 of file lsd.c.

References RELATIVE_ERROR_FACTOR, and TRUE.

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

00182 {
00183   double abs_diff,aa,bb,abs_max;
00184 
00185   /* trivial case */
00186   if( a == b ) return TRUE;
00187 
00188   abs_diff = fabs(a-b);
00189   aa = fabs(a);
00190   bb = fabs(b);
00191   abs_max = aa > bb ? aa : bb;
00192 
00193   /* DBL_MIN is the smallest normalized number, thus, the smallest
00194      number whose relative error is bounded by DBL_EPSILON. For
00195      smaller numbers, the same quantization steps as for DBL_MIN
00196      are used. Then, for smaller numbers, a meaningful "relative"
00197      error should be computed by dividing the difference by DBL_MIN. */
00198   if( abs_max < DBL_MIN ) abs_max = DBL_MIN;
00199 
00200   /* equal if relative error <= factor x eps */
00201   return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
00202 }

void enlarge_ntuple_list ntuple_list  n_tuple  )  [static]
 

Enlarge the allocated memory of an n-tuple list.

Definition at line 287 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_7tuple(), and gaussian_kernel().

00288 {
00289   /* check parameters */
00290   if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )
00291     error("enlarge_ntuple_list: invalid n-tuple.");
00292 
00293   /* duplicate number of tuples */
00294   n_tuple->max_size *= 2;
00295 
00296   /* realloc memory */
00297   n_tuple->values = (double *) realloc( (void *) n_tuple->values,
00298                       n_tuple->dim * n_tuple->max_size * sizeof(double) );
00299   if( n_tuple->values == NULL ) error("not enough memory.");
00300 }

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 159 of file lsd.c.

Referenced by add_7tuple(), enlarge_ntuple_list(), free_image_char(), free_image_double(), 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_double_ptr(), 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().

00160 {
00161   fprintf(stderr,"LSD Error: %s\n",msg);
00162   exit(EXIT_FAILURE);
00163 }

void free_image_char image_char  i  )  [static]
 

Free memory used in image_char 'i'.

Definition at line 352 of file lsd.c.

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

Referenced by LineSegmentDetection().

00353 {
00354   if( i == NULL || i->data == NULL )
00355     error("free_image_char: invalid input image.");
00356   free( (void *) i->data );
00357   free( (void *) i );
00358 }

Here is the call graph for this function:

void free_image_double image_double  i  )  [static]
 

Free memory used in image_double 'i'.

Definition at line 478 of file lsd.c.

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

Referenced by gaussian_sampler(), and LineSegmentDetection().

00479 {
00480   if( i == NULL || i->data == NULL )
00481     error("free_image_double: invalid input image.");
00482   free( (void *) i->data );
00483   free( (void *) i );
00484 }

Here is the call graph for this function:

void free_ntuple_list ntuple_list  in  )  [static]
 

Free memory used in n-tuple 'in'.

Definition at line 249 of file lsd.c.

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

Referenced by gaussian_sampler().

00250 {
00251   if( in == NULL || in->values == NULL )
00252     error("free_ntuple_list: invalid n-tuple input.");
00253   free( (void *) in->values );
00254   free( (void *) in );
00255 }

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 548 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().

00549 {
00550   double sum = 0.0;
00551   double val;
00552   unsigned int i;
00553 
00554   /* check parameters */
00555   if( kernel == NULL || kernel->values == NULL )
00556     error("gaussian_kernel: invalid n-tuple 'kernel'.");
00557   if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");
00558 
00559   /* compute Gaussian kernel */
00560   if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);
00561   kernel->size = 1;
00562   for(i=0;i<kernel->dim;i++)
00563     {
00564       val = ( (double) i - mean ) / sigma;
00565       kernel->values[i] = exp( -0.5 * val * val );
00566       sum += kernel->values[i];
00567     }
00568 
00569   /* normalization */
00570   if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
00571 }

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 applies 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 611 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().

00613 {
00614   image_double aux,out;
00615   ntuple_list kernel;
00616   unsigned int N,M,h,n,x,y,i;
00617   int xc,yc,j,double_x_size,double_y_size;
00618   double sigma,xx,yy,sum,prec;
00619 
00620   /* check parameters */
00621   if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00622     error("gaussian_sampler: invalid image.");
00623   if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");
00624   if( sigma_scale <= 0.0 )
00625     error("gaussian_sampler: 'sigma_scale' must be positive.");
00626 
00627   /* compute new image size and get memory for images */
00628   if( in->xsize * scale > (double) UINT_MAX ||
00629       in->ysize * scale > (double) UINT_MAX )
00630     error("gaussian_sampler: the output image size exceeds the handled size.");
00631   N = (unsigned int) ceil( in->xsize * scale );
00632   M = (unsigned int) ceil( in->ysize * scale );
00633   aux = new_image_double(N,in->ysize);
00634   out = new_image_double(N,M);
00635 
00636   /* sigma, kernel size and memory for the kernel */
00637   sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;
00638   /*
00639      The size of the kernel is selected to guarantee that the
00640      the first discarded term is at least 10^prec times smaller
00641      than the central value. For that, h should be larger than x, with
00642        e^(-x^2/2sigma^2) = 1/10^prec.
00643      Then,
00644        x = sigma * sqrt( 2 * prec * ln(10) ).
00645    */
00646   prec = 3.0;
00647   h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
00648   n = 1+2*h; /* kernel size */
00649   kernel = new_ntuple_list(n);
00650 
00651   /* auxiliary double image size variables */
00652   double_x_size = (int) (2 * in->xsize);
00653   double_y_size = (int) (2 * in->ysize);
00654 
00655   /* First subsampling: x axis */
00656   for(x=0;x<aux->xsize;x++)
00657     {
00658       /*
00659          x   is the coordinate in the new image.
00660          xx  is the corresponding x-value in the original size image.
00661          xc  is the integer value, the pixel coordinate of xx.
00662        */
00663       xx = (double) x / scale;
00664       /* coordinate (0.0,0.0) is in the center of pixel (0,0),
00665          so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */
00666       xc = (int) floor( xx + 0.5 );
00667       gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );
00668       /* the kernel must be computed for each x because the fine
00669          offset xx-xc is different in each case */
00670 
00671       for(y=0;y<aux->ysize;y++)
00672         {
00673           sum = 0.0;
00674           for(i=0;i<kernel->dim;i++)
00675             {
00676               j = xc - h + i;
00677 
00678               /* symmetry boundary condition */
00679               while( j < 0 ) j += double_x_size;
00680               while( j >= double_x_size ) j -= double_x_size;
00681               if( j >= (int) in->xsize ) j = double_x_size-1-j;
00682 
00683               sum += in->data[ j + y * in->xsize ] * kernel->values[i];
00684             }
00685           aux->data[ x + y * aux->xsize ] = sum;
00686         }
00687     }
00688 
00689   /* Second subsampling: y axis */
00690   for(y=0;y<out->ysize;y++)
00691     {
00692       /*
00693          y   is the coordinate in the new image.
00694          yy  is the corresponding x-value in the original size image.
00695          yc  is the integer value, the pixel coordinate of xx.
00696        */
00697       yy = (double) y / scale;
00698       /* coordinate (0.0,0.0) is in the center of pixel (0,0),
00699          so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */
00700       yc = (int) floor( yy + 0.5 );
00701       gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );
00702       /* the kernel must be computed for each y because the fine
00703          offset yy-yc is different in each case */
00704 
00705       for(x=0;x<out->xsize;x++)
00706         {
00707           sum = 0.0;
00708           for(i=0;i<kernel->dim;i++)
00709             {
00710               j = yc - h + i;
00711 
00712               /* symmetry boundary condition */
00713               while( j < 0 ) j += double_y_size;
00714               while( j >= double_y_size ) j -= double_y_size;
00715               if( j >= (int) in->ysize ) j = double_y_size-1-j;
00716 
00717               sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];
00718             }
00719           out->data[ x + y * out->xsize ] = sum;
00720         }
00721     }
00722 
00723   /* free memory */
00724   free_ntuple_list(kernel);
00725   free_image_double(aux);
00726 
00727   return out;
00728 }

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 associated to the smallest eigenvalue. We have to solve for 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 1568 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().

01570 {
01571   double lambda,theta,weight;
01572   double Ixx = 0.0;
01573   double Iyy = 0.0;
01574   double Ixy = 0.0;
01575   int i;
01576 
01577   /* check parameters */
01578   if( reg == NULL ) error("get_theta: invalid region.");
01579   if( reg_size <= 1 ) error("get_theta: region size <= 1.");
01580   if( modgrad == NULL || modgrad->data == NULL )
01581     error("get_theta: invalid 'modgrad'.");
01582   if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");
01583 
01584   /* compute inertia matrix */
01585   for(i=0; i<reg_size; i++)
01586     {
01587       weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01588       Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;
01589       Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;
01590       Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;
01591     }
01592   if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )
01593     error("get_theta: null inertia matrix.");
01594 
01595   /* compute smallest eigenvalue */
01596   lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );
01597 
01598   /* compute angle */
01599   theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);
01600 
01601   /* The previous procedure doesn't cares about orientation,
01602      so it could be wrong by 180 degrees. Here is corrected if necessary. */
01603   if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;
01604 
01605   return theta;
01606 }

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 1299 of file lsd.c.

References double_equal(), and error().

Referenced by ri_inc().

01300 {
01301   /* check parameters */
01302   if( x1 > x2 || x < x1 || x > x2 )
01303     error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01304 
01305   /* interpolation */
01306   if( double_equal(x1,x2) && y1<y2 ) return y2;
01307   if( double_equal(x1,x2) && y1>y2 ) return y1;
01308   return y1 + (x-x1) * (y2-y1) / (x2-x1);
01309 }

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 1277 of file lsd.c.

References double_equal(), and error().

Referenced by ri_inc().

01278 {
01279   /* check parameters */
01280   if( x1 > x2 || x < x1 || x > x2 )
01281     error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01282 
01283   /* interpolation */
01284   if( double_equal(x1,x2) && y1<y2 ) return y1;
01285   if( double_equal(x1,x2) && y1>y2 ) return y2;
01286   return y1 + (x-x1) * (y2-y1) / (x2-x1);
01287 }

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 893 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().

00895 {
00896   double a;
00897 
00898   /* check parameters */
00899   if( angles == NULL || angles->data == NULL )
00900     error("isaligned: invalid image 'angles'.");
00901   if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
00902     error("isaligned: (x,y) out of the image.");
00903   if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");
00904 
00905   /* angle at pixel (x,y) */
00906   a = angles->data[ x + y * angles->xsize ];
00907 
00908   /* pixels whose level-line angle is not defined
00909      are considered as NON-aligned */
00910   if( a == NOTDEF ) return FALSE;  /* there is no need to call the function
00911                                       'double_equal' here because there is
00912                                       no risk of problems related to the
00913                                       comparison doubles, we are only
00914                                       interested in the exact NOTDEF value */
00915 
00916   /* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */
00917   theta -= a;
00918   if( theta < 0.0 ) theta = -theta;
00919   if( theta > M_3_2_PI )
00920     {
00921       theta -= M_2__PI;
00922       if( theta < 0.0 ) theta = -theta;
00923     }
00924 
00925   return theta <= prec;
00926 }

Here is the call graph for this function:

double* LineSegmentDetection int *  n_out,
double *  img,
int  X,
int  Y,
double  scale,
double  sigma_scale,
double  quant,
double  ang_th,
double  log_eps,
double  density_th,
int  n_bins,
int **  reg_img,
int *  reg_x,
int *  reg_y
 

LSD Full Interface.

Parameters:
n_out Pointer to an int where LSD will store the number of line segments detected.
img Pointer to input image data. It must be an array of doubles of size X x Y, and the pixel at coordinates (x,y) is obtained by img[x+y*X].
X X size of the image: the number of columns.
Y Y size of the image: the number of rows.
scale When different from 1.0, LSD will scale the input image by 'scale' factor by Gaussian filtering, before detecting line segments. Example: if scale=0.8, the input image will be subsampled to 80% of its size, before the line segment detector is 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 levels are 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
log_eps Detection threshold, accept if -log10(NFA) > log_eps. The larger the value, 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 gives an average of 10 false detections on noise
  • 0.0 gives an average of 1 false detections on noise
  • 1.0 gives an average of 0.1 false detections on nose
  • 2.0 gives an average of 0.01 false detections on noise
Suggested value: 0.0
density_th Minimal proportion of 'supporting' points in a rectangle. Suggested value: 0.7
n_bins Number of bins used in the pseudo-ordering of gradient modulus. Suggested value: 1024
reg_img Optional output: if desired, LSD will return an int image where each pixel indicates the line segment to which it belongs. Unused pixels have the value '0', while the used ones have the number of the line segment, numbered 1,2,3,..., in the same order as in the output list. If desired, a non NULL int** pointer must be assigned, and LSD will make that the pointer point to an int array of size reg_x x reg_y, where the pixel value at (x,y) is obtained with (*reg_img)[x+y*reg_x]. Note that 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'. If scale!=1 this size differs from XxY and that is the reason why its value is given by reg_x and reg_y. Suggested value: NULL
reg_x Pointer to an int where LSD will put the X size 'reg_img' image, when asked for. Suggested value: NULL
reg_y Pointer to an int where LSD will put the Y size 'reg_img' image, when asked for. Suggested value: NULL
Returns:
A double array of size 7 x n_out, containing the list of line segments detected. The array contains first 7 values of line segment number 1, then the 7 values of line segment number 2, and so on, and it finish by the 7 values of line segment number n_out. The seven values are:
  • x1,y1,x2,y2,width,p,-log10(NFA)
for a line segment from coordinates (x1,y1) to (x2,y2), a width 'width', an angle precision of p in (0,1) given by angle_tolerance/180 degree, and NFA value 'NFA'. If 'out' is the returned pointer, the 7 values of line segment number 'n+1' are obtained with 'out[7*n+0]' to 'out[7*n+6]'.

Definition at line 2025 of file lsd.c.

References add_7tuple(), image_int_s::data, image_double_s::data, image_char_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_double_ptr(), new_image_int_ini(), new_ntuple_list(), NOTDEF, NOTUSED, ntuple_list, rect::p, rect_improve(), refine(), region2rect(), region_grow(), ntuple_list_s::size, ntuple_list_s::values, rect::width, rect::x1, rect::x2, image_int_s::xsize, image_char_s::xsize, image_double_s::xsize, rect::y1, rect::y2, image_int_s::ysize, and image_double_s::ysize.

Referenced by lsd_scale_region().

02031 {
02032   image_double image;
02033   ntuple_list out = new_ntuple_list(7);
02034   double * return_value;
02035   image_double scaled_image,angles,modgrad;
02036   image_char used;
02037   image_int region = NULL;
02038   struct coorlist * list_p;
02039   void * mem_p;
02040   struct rect rec;
02041   struct point * reg;
02042   int reg_size,min_reg_size,i;
02043   unsigned int xsize,ysize;
02044   double rho,reg_angle,prec,p,log_nfa,logNT;
02045   int ls_count = 0;                   /* line segments are numbered 1,2,3,... */
02046 
02047 
02048   /* check parameters */
02049   if( img == NULL || X <= 0 || Y <= 0 ) error("invalid image input.");
02050   if( scale <= 0.0 ) error("'scale' value must be positive.");
02051   if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");
02052   if( quant < 0.0 ) error("'quant' value must be positive.");
02053   if( ang_th <= 0.0 || ang_th >= 180.0 )
02054     error("'ang_th' value must be in the range (0,180).");
02055   if( density_th < 0.0 || density_th > 1.0 )
02056     error("'density_th' value must be in the range [0,1].");
02057   if( n_bins <= 0 ) error("'n_bins' value must be positive.");
02058 
02059 
02060   /* angle tolerance */
02061   prec = M_PI * ang_th / 180.0;
02062   p = ang_th / 180.0;
02063   rho = quant / sin(prec); /* gradient magnitude threshold */
02064 
02065 
02066   /* load and scale image (if necessary) and compute angle at each pixel */
02067   image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );
02068   if( scale != 1.0 )
02069     {
02070       scaled_image = gaussian_sampler( image, scale, sigma_scale );
02071       angles = ll_angle( scaled_image, rho, &list_p, &mem_p,
02072                          &modgrad, (unsigned int) n_bins );
02073       free_image_double(scaled_image);
02074     }
02075   else
02076     angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,
02077                        (unsigned int) n_bins );
02078   xsize = angles->xsize;
02079   ysize = angles->ysize;
02080 
02081   /* Number of Tests - NT
02082 
02083      The theoretical number of tests is Np.(XY)^(5/2)
02084      where X and Y are number of columns and rows of the image.
02085      Np corresponds to the number of angle precisions considered.
02086      As the procedure 'rect_improve' tests 5 times to halve the
02087      angle precision, and 5 more times after improving other factors,
02088      11 different precision values are potentially tested. Thus,
02089      the number of tests is
02090        11 * (X*Y)^(5/2)
02091      whose logarithm value is
02092        log10(11) + 5/2 * (log10(X) + log10(Y)).
02093   */
02094   logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0
02095           + log10(11.0);
02096   min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region
02097                                              that can give a meaningful event */
02098 
02099 
02100   /* initialize some structures */
02101   if( reg_img != NULL && reg_x != NULL && reg_y != NULL ) /* save region data */
02102     region = new_image_int_ini(angles->xsize,angles->ysize,0);
02103   used = new_image_char_ini(xsize,ysize,NOTUSED);
02104   reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );
02105   if( reg == NULL ) error("not enough memory!");
02106 
02107 
02108   /* search for line segments */
02109   for(; list_p != NULL; list_p = list_p->next )
02110     if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&
02111         angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )
02112        /* there is no risk of double comparison problems here
02113           because we are only interested in the exact NOTDEF value */
02114       {
02115         /* find the region of connected point and ~equal angle */
02116         region_grow( list_p->x, list_p->y, angles, reg, &reg_size,
02117                      &reg_angle, used, prec );
02118 
02119         /* reject small regions */
02120         if( reg_size < min_reg_size ) continue;
02121 
02122         /* construct rectangular approximation for the region */
02123         region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);
02124 
02125         /* Check if the rectangle exceeds the minimal density of
02126            region points. If not, try to improve the region.
02127            The rectangle will be rejected if the final one does
02128            not fulfill the minimal density condition.
02129            This is an addition to the original LSD algorithm published in
02130            "LSD: A Fast Line Segment Detector with a False Detection Control"
02131            by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall.
02132            The original algorithm is obtained with density_th = 0.0.
02133          */
02134         if( !refine( reg, &reg_size, modgrad, reg_angle,
02135                      prec, p, &rec, used, angles, density_th ) ) continue;
02136 
02137         /* compute NFA value */
02138         log_nfa = rect_improve(&rec,angles,logNT,log_eps);
02139         if( log_nfa <= log_eps ) continue;
02140 
02141         /* A New Line Segment was found! */
02142         ++ls_count;  /* increase line segment counter */
02143 
02144         /*
02145            The gradient was computed with a 2x2 mask, its value corresponds to
02146            points with an offset of (0.5,0.5), that should be added to output.
02147            The coordinates origin is at the center of pixel (0,0).
02148          */
02149         rec.x1 += 0.5; rec.y1 += 0.5;
02150         rec.x2 += 0.5; rec.y2 += 0.5;
02151 
02152         /* scale the result values if a subsampling was performed */
02153         if( scale != 1.0 )
02154           {
02155             rec.x1 /= scale; rec.y1 /= scale;
02156             rec.x2 /= scale; rec.y2 /= scale;
02157             rec.width /= scale;
02158           }
02159 
02160         /* add line segment found to output */
02161         add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,
02162                          rec.width, rec.p, log_nfa );
02163 
02164         /* add region number to 'region' image if needed */
02165         if( region != NULL )
02166           for(i=0; i<reg_size; i++)
02167             region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;
02168       }
02169 
02170 
02171   /* free memory */
02172   free( (void *) image );   /* only the double_image structure should be freed,
02173                                the data pointer was provided to this functions
02174                                and should not be destroyed.                 */
02175   free_image_double(angles);
02176   free_image_double(modgrad);
02177   free_image_char(used);
02178   free( (void *) reg );
02179   free( (void *) mem_p );
02180 
02181   /* return the result */
02182   if( reg_img != NULL && reg_x != NULL && reg_y != NULL )
02183     {
02184       if( region == NULL ) error("'region' should be a valid image.");
02185       *reg_img = region->data;
02186       if( region->xsize > (unsigned int) INT_MAX ||
02187           region->xsize > (unsigned int) INT_MAX )
02188         error("region image to big to fit in INT sizes.");
02189       *reg_x = (int) (region->xsize);
02190       *reg_y = (int) (region->ysize);
02191 
02192       /* free the 'region' structure.
02193          we cannot use the function 'free_image_int' because we need to keep
02194          the memory with the image data to be returned by this function. */
02195       free( (void *) region );
02196     }
02197   if( out->size > (unsigned int) INT_MAX )
02198     error("too many detections to fit in an INT.");
02199   *n_out = (int) (out->size);
02200 
02201   return_value = out->values;
02202   free( (void *) out );  /* only the 'ntuple_list' structure must be freed,
02203                             but the 'values' pointer must be keep to return
02204                             as a result. */
02205 
02206   return return_value;
02207 }

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
[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 752 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().

00755 {
00756   image_double g;
00757   unsigned int n,p,x,y,adr,i;
00758   double com1,com2,gx,gy,norm,norm2;
00759   /* the rest of the variables are used for pseudo-ordering
00760      the gradient magnitude values */
00761   int list_count = 0;
00762   struct coorlist * list;
00763   struct coorlist ** range_l_s; /* array of pointers to start of bin list */
00764   struct coorlist ** range_l_e; /* array of pointers to end of bin list */
00765   struct coorlist * start;
00766   struct coorlist * end;
00767   double max_grad = 0.0;
00768 
00769   /* check parameters */
00770   if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00771     error("ll_angle: invalid image.");
00772   if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");
00773   if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");
00774   if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");
00775   if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");
00776   if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");
00777 
00778   /* image size shortcuts */
00779   n = in->ysize;
00780   p = in->xsize;
00781 
00782   /* allocate output image */
00783   g = new_image_double(in->xsize,in->ysize);
00784 
00785   /* get memory for the image of gradient modulus */
00786   *modgrad = new_image_double(in->xsize,in->ysize);
00787 
00788   /* get memory for "ordered" list of pixels */
00789   list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );
00790   *mem_p = (void *) list;
00791   range_l_s = (struct coorlist **) calloc( (size_t) n_bins,
00792                                            sizeof(struct coorlist *) );
00793   range_l_e = (struct coorlist **) calloc( (size_t) n_bins,
00794                                            sizeof(struct coorlist *) );
00795   if( list == NULL || range_l_s == NULL || range_l_e == NULL )
00796     error("not enough memory.");
00797   for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;
00798 
00799   /* 'undefined' on the down and right boundaries */
00800   for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;
00801   for(y=0;y<n;y++) g->data[p*y+p-1]   = NOTDEF;
00802 
00803   /* compute gradient on the remaining pixels */
00804   for(x=0;x<p-1;x++)
00805     for(y=0;y<n-1;y++)
00806       {
00807         adr = y*p+x;
00808 
00809         /*
00810            Norm 2 computation using 2x2 pixel window:
00811              A B
00812              C D
00813            and
00814              com1 = D-A,  com2 = B-C.
00815            Then
00816              gx = B+D - (A+C)   horizontal difference
00817              gy = C+D - (A+B)   vertical difference
00818            com1 and com2 are just to avoid 2 additions.
00819          */
00820         com1 = in->data[adr+p+1] - in->data[adr];
00821         com2 = in->data[adr+1]   - in->data[adr+p];
00822 
00823         gx = com1+com2; /* gradient x component */
00824         gy = com1-com2; /* gradient y component */
00825         norm2 = gx*gx+gy*gy;
00826         norm = sqrt( norm2 / 4.0 ); /* gradient norm */
00827 
00828         (*modgrad)->data[adr] = norm; /* store gradient norm */
00829 
00830         if( norm <= threshold ) /* norm too small, gradient no defined */
00831           g->data[adr] = NOTDEF; /* gradient angle not defined */
00832         else
00833           {
00834             /* gradient angle computation */
00835             g->data[adr] = atan2(gx,-gy);
00836 
00837             /* look for the maximum of the gradient */
00838             if( norm > max_grad ) max_grad = norm;
00839           }
00840       }
00841 
00842   /* compute histogram of gradient values */
00843   for(x=0;x<p-1;x++)
00844     for(y=0;y<n-1;y++)
00845       {
00846         norm = (*modgrad)->data[y*p+x];
00847 
00848         /* store the point in the right bin according to its norm */
00849         i = (unsigned int) (norm * (double) n_bins / max_grad);
00850         if( i >= n_bins ) i = n_bins-1;
00851         if( range_l_e[i] == NULL )
00852           range_l_s[i] = range_l_e[i] = list+list_count++;
00853         else
00854           {
00855             range_l_e[i]->next = list+list_count;
00856             range_l_e[i] = list+list_count++;
00857           }
00858         range_l_e[i]->x = (int) x;
00859         range_l_e[i]->y = (int) y;
00860         range_l_e[i]->next = NULL;
00861       }
00862 
00863   /* Make the list of pixels (almost) ordered by norm value.
00864      It starts by the larger bin, so the list starts by the
00865      pixels with the highest gradient value. Pixels would be ordered
00866      by norm value, up to a precision given by max_grad/n_bins.
00867    */
00868   for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);
00869   start = range_l_s[i];
00870   end = range_l_e[i];
00871   if( start != NULL )
00872     while(i>0)
00873       {
00874         --i;
00875         if( range_l_s[i] != NULL )
00876           {
00877             end->next = range_l_s[i];
00878             end = range_l_e[i];
00879           }
00880       }
00881   *list_p = start;
00882 
00883   /* free memory */
00884   free( (void *) range_l_s );
00885   free( (void *) range_l_e );
00886 
00887   return g;
00888 }

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 980 of file lsd.c.

00981 {
00982   static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
00983                          8687.24529705, 1168.92649479, 83.8676043424,
00984                          2.50662827511 };
00985   double a = (x+0.5) * log(x+5.5) - (x+5.5);
00986   double b = 0.0;
00987   int n;
00988 
00989   for(n=0;n<7;n++)
00990     {
00991       a -= log( x + (double) n );
00992       b += q[n] * pow( x, (double) n );
00993     }
00994   return a + log(b);
00995 }

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 1014 of file lsd.c.

01015 {
01016   return 0.918938533204673 + (x-0.5)*log(x) - x
01017          + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
01018 }

double* lsd int *  n_out,
double *  img,
int  X,
int  Y
 

LSD Simple Interface.

Parameters:
n_out Pointer to an int where LSD will store the number of line segments detected.
img Pointer to input image data. It must be an array of doubles of size X x Y, and the pixel at coordinates (x,y) is obtained by img[x+y*X].
X X size of the image: the number of columns.
Y Y size of the image: the number of rows.
Returns:
A double array of size 7 x n_out, containing the list of line segments detected. The array contains first 7 values of line segment number 1, then the 7 values of line segment number 2, and so on, and it finish by the 7 values of line segment number n_out. The seven values are:
  • x1,y1,x2,y2,width,p,-log10(NFA)
for a line segment from coordinates (x1,y1) to (x2,y2), a width 'width', an angle precision of p in (0,1) given by angle_tolerance/180 degree, and NFA value 'NFA'. If 'out' is the returned pointer, the 7 values of line segment number 'n+1' are obtained with 'out[7*n+0]' to 'out[7*n+6]'.

Definition at line 2243 of file lsd.c.

References lsd_scale().

02244 {
02245   /* LSD parameters */
02246   double scale = 0.8;       /* Scale the image by Gaussian filter to 'scale'. */
02247 
02248   return lsd_scale(n_out,img,X,Y,scale);
02249 }

Here is the call graph for this function:

double* lsd_scale int *  n_out,
double *  img,
int  X,
int  Y,
double  scale
 

LSD Simple Interface with Scale.

Parameters:
n_out Pointer to an int where LSD will store the number of line segments detected.
img Pointer to input image data. It must be an array of doubles of size X x Y, and the pixel at coordinates (x,y) is obtained by img[x+y*X].
X X size of the image: the number of columns.
Y Y size of the image: the number of rows.
scale When different from 1.0, LSD will scale the input image by 'scale' factor by Gaussian filtering, before detecting line segments. Example: if scale=0.8, the input image will be subsampled to 80% of its size, before the line segment detector is applied. Suggested value: 0.8
Returns:
A double array of size 7 x n_out, containing the list of line segments detected. The array contains first 7 values of line segment number 1, then the 7 values of line segment number 2, and so on, and it finish by the 7 values of line segment number n_out. The seven values are:
  • x1,y1,x2,y2,width,p,-log10(NFA)
for a line segment from coordinates (x1,y1) to (x2,y2), a width 'width', an angle precision of p in (0,1) given by angle_tolerance/180 degree, and NFA value 'NFA'. If 'out' is the returned pointer, the 7 values of line segment number 'n+1' are obtained with 'out[7*n+0]' to 'out[7*n+6]'.

Definition at line 2235 of file lsd.c.

References lsd_scale_region().

Referenced by lsd().

02236 {
02237   return lsd_scale_region(n_out,img,X,Y,scale,NULL,NULL,NULL);
02238 }

Here is the call graph for this function:

double* lsd_scale_region int *  n_out,
double *  img,
int  X,
int  Y,
double  scale,
int **  reg_img,
int *  reg_x,
int *  reg_y
 

LSD Simple Interface with Scale and Region output.

Parameters:
n_out Pointer to an int where LSD will store the number of line segments detected.
img Pointer to input image data. It must be an array of doubles of size X x Y, and the pixel at coordinates (x,y) is obtained by img[x+y*X].
X X size of the image: the number of columns.
Y Y size of the image: the number of rows.
scale When different from 1.0, LSD will scale the input image by 'scale' factor by Gaussian filtering, before detecting line segments. Example: if scale=0.8, the input image will be subsampled to 80% of its size, before the line segment detector is applied. Suggested value: 0.8
reg_img Optional output: if desired, LSD will return an int image where each pixel indicates the line segment to which it belongs. Unused pixels have the value '0', while the used ones have the number of the line segment, numbered 1,2,3,..., in the same order as in the output list. If desired, a non NULL int** pointer must be assigned, and LSD will make that the pointer point to an int array of size reg_x x reg_y, where the pixel value at (x,y) is obtained with (*reg_img)[x+y*reg_x]. Note that 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'. If scale!=1 this size differs from XxY and that is the reason why its value is given by reg_x and reg_y. Suggested value: NULL
reg_x Pointer to an int where LSD will put the X size 'reg_img' image, when asked for. Suggested value: NULL
reg_y Pointer to an int where LSD will put the Y size 'reg_img' image, when asked for. Suggested value: NULL
Returns:
A double array of size 7 x n_out, containing the list of line segments detected. The array contains first 7 values of line segment number 1, then the 7 values of line segment number 2, and so on, and it finish by the 7 values of line segment number n_out. The seven values are:
  • x1,y1,x2,y2,width,p,-log10(NFA)
for a line segment from coordinates (x1,y1) to (x2,y2), a width 'width', an angle precision of p in (0,1) given by angle_tolerance/180 degree, and NFA value 'NFA'. If 'out' is the returned pointer, the 7 values of line segment number 'n+1' are obtained with 'out[7*n+0]' to 'out[7*n+6]'.

Definition at line 2212 of file lsd.c.

References LineSegmentDetection().

Referenced by lsd_scale().

02215 {
02216   /* LSD parameters */
02217   double sigma_scale = 0.6; /* Sigma for Gaussian filter is computed as
02218                                 sigma = sigma_scale/scale.                    */
02219   double quant = 2.0;       /* Bound to the quantization error on the
02220                                 gradient norm.                                */
02221   double ang_th = 22.5;     /* Gradient angle tolerance in degrees.           */
02222   double log_eps = 0.0;     /* Detection threshold: -log10(NFA) > log_eps     */
02223   double density_th = 0.7;  /* Minimal density of region points in rectangle. */
02224   int n_bins = 1024;        /* Number of bins in pseudo-ordering of gradient
02225                                modulus.                                       */
02226 
02227   return LineSegmentDetection( n_out, img, X, Y, scale, sigma_scale, quant,
02228                                ang_th, log_eps, density_th, n_bins,
02229                                reg_img, reg_x, reg_y );
02230 }

Here is the call graph for this function:

image_char new_image_char unsigned int  xsize,
unsigned int  ysize
[static]
 

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

Definition at line 363 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().

00364 {
00365   image_char image;
00366 
00367   /* check parameters */
00368   if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");
00369 
00370   /* get memory */
00371   image = (image_char) malloc( sizeof(struct image_char_s) );
00372   if( image == NULL ) error("not enough memory.");
00373   image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),
00374                                           sizeof(unsigned char) );
00375   if( image->data == NULL ) error("not enough memory.");
00376 
00377   /* set image size */
00378   image->xsize = xsize;
00379   image->ysize = ysize;
00380 
00381   return image;
00382 }

Here is the call graph for this function:

image_char new_image_char_ini unsigned int  xsize,
unsigned int  ysize,
unsigned char  fill_value
[static]
 

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

Definition at line 388 of file lsd.c.

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

Referenced by LineSegmentDetection().

00390 {
00391   image_char image = new_image_char(xsize,ysize); /* create image */
00392   unsigned int N = xsize*ysize;
00393   unsigned int i;
00394 
00395   /* check parameters */
00396   if( image == NULL || image->data == NULL )
00397     error("new_image_char_ini: invalid image.");
00398 
00399   /* initialize */
00400   for(i=0; i<N; i++) image->data[i] = fill_value;
00401 
00402   return image;
00403 }

Here is the call graph for this function:

image_double new_image_double unsigned int  xsize,
unsigned int  ysize
[static]
 

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

Definition at line 489 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(), and ll_angle().

00490 {
00491   image_double image;
00492 
00493   /* check parameters */
00494   if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");
00495 
00496   /* get memory */
00497   image = (image_double) malloc( sizeof(struct image_double_s) );
00498   if( image == NULL ) error("not enough memory.");
00499   image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );
00500   if( image->data == NULL ) error("not enough memory.");
00501 
00502   /* set image size */
00503   image->xsize = xsize;
00504   image->ysize = ysize;
00505 
00506   return image;
00507 }

Here is the call graph for this function:

image_double new_image_double_ptr unsigned int  xsize,
unsigned int  ysize,
double *  data
[static]
 

Create a new image_double of size 'xsize' times 'ysize' with the data pointed by 'data'.

Definition at line 513 of file lsd.c.

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

Referenced by LineSegmentDetection().

00515 {
00516   image_double image;
00517 
00518   /* check parameters */
00519   if( xsize == 0 || ysize == 0 )
00520     error("new_image_double_ptr: invalid image size.");
00521   if( data == NULL ) error("new_image_double_ptr: NULL data pointer.");
00522 
00523   /* get memory */
00524   image = (image_double) malloc( sizeof(struct image_double_s) );
00525   if( image == NULL ) error("not enough memory.");
00526 
00527   /* set image */
00528   image->xsize = xsize;
00529   image->ysize = ysize;
00530   image->data = data;
00531 
00532   return image;
00533 }

Here is the call graph for this function:

image_int new_image_int unsigned int  xsize,
unsigned int  ysize
[static]
 

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

Definition at line 423 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().

00424 {
00425   image_int image;
00426 
00427   /* check parameters */
00428   if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");
00429 
00430   /* get memory */
00431   image = (image_int) malloc( sizeof(struct image_int_s) );
00432   if( image == NULL ) error("not enough memory.");
00433   image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );
00434   if( image->data == NULL ) error("not enough memory.");
00435 
00436   /* set image size */
00437   image->xsize = xsize;
00438   image->ysize = ysize;
00439 
00440   return image;
00441 }

Here is the call graph for this function:

image_int new_image_int_ini unsigned int  xsize,
unsigned int  ysize,
int  fill_value
[static]
 

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

Definition at line 447 of file lsd.c.

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

Referenced by LineSegmentDetection().

00449 {
00450   image_int image = new_image_int(xsize,ysize); /* create image */
00451   unsigned int N = xsize*ysize;
00452   unsigned int i;
00453 
00454   /* initialize */
00455   for(i=0; i<N; i++) image->data[i] = fill_value;
00456 
00457   return image;
00458 }

Here is the call graph for this function:

ntuple_list new_ntuple_list unsigned int  dim  )  [static]
 

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

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

Definition at line 261 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().

00262 {
00263   ntuple_list n_tuple;
00264 
00265   /* check parameters */
00266   if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");
00267 
00268   /* get memory for list structure */
00269   n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );
00270   if( n_tuple == NULL ) error("not enough memory.");
00271 
00272   /* initialize list */
00273   n_tuple->size = 0;
00274   n_tuple->max_size = 1;
00275   n_tuple->dim = dim;
00276 
00277   /* get memory for tuples */
00278   n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) );
00279   if( n_tuple->values == NULL ) error("not enough memory.");
00280 
00281   return n_tuple;
00282 }

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 1074 of file lsd.c.

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

Referenced by rect_nfa().

01075 {
01076   static double inv[TABSIZE];   /* table to keep computed inverse values */
01077   double tolerance = 0.1;       /* an error of 10% in the result is accepted */
01078   double log1term,term,bin_term,mult_term,bin_tail,err,p_term;
01079   int i;
01080 
01081   /* check parameters */
01082   if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )
01083     error("nfa: wrong n, k or p values.");
01084 
01085   /* trivial cases */
01086   if( n==0 || k==0 ) return -logNT;
01087   if( n==k ) return -logNT - (double) n * log10(p);
01088 
01089   /* probability term */
01090   p_term = p / (1.0-p);
01091 
01092   /* compute the first term of the series */
01093   /*
01094      binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
01095      where bincoef(n,i) are the binomial coefficients.
01096      But
01097        bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
01098      We use this to compute the first term. Actually the log of it.
01099    */
01100   log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )
01101            - log_gamma( (double) (n-k) + 1.0 )
01102            + (double) k * log(p) + (double) (n-k) * log(1.0-p);
01103   term = exp(log1term);
01104 
01105   /* in some cases no more computations are needed */
01106   if( double_equal(term,0.0) )              /* the first term is almost zero */
01107     {
01108       if( (double) k > (double) n * p )     /* at begin or end of the tail?  */
01109         return -log1term / M_LN10 - logNT;  /* end: use just the first term  */
01110       else
01111         return -logNT;                      /* begin: the tail is roughly 1  */
01112     }
01113 
01114   /* compute more terms if needed */
01115   bin_tail = term;
01116   for(i=k+1;i<=n;i++)
01117     {
01118       /*
01119          As
01120            term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
01121          and
01122            bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,
01123          then,
01124            term_i / term_i-1 = (n-i+1)/i * p/(1-p)
01125          and
01126            term_i = term_i-1 * (n-i+1)/i * p/(1-p).
01127          1/i is stored in a table as they are computed,
01128          because divisions are expensive.
01129          p/(1-p) is computed only once and stored in 'p_term'.
01130        */
01131       bin_term = (double) (n-i+1) * ( i<TABSIZE ?
01132                    ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :
01133                    1.0 / (double) i );
01134 
01135       mult_term = bin_term * p_term;
01136       term *= mult_term;
01137       bin_tail += term;
01138       if(bin_term<1.0)
01139         {
01140           /* When bin_term<1 then mult_term_j<mult_term_i for j>i.
01141              Then, the error on the binomial tail when truncated at
01142              the i term can be bounded by a geometric series of form
01143              term_i * sum mult_term_i^j.                            */
01144           err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /
01145                          (1.0-mult_term) - 1.0 );
01146 
01147           /* One wants an error at most of tolerance*final_result, or:
01148              tolerance * abs(-log10(bin_tail)-logNT).
01149              Now, the error that can be accepted on bin_tail is
01150              given by tolerance*final_result divided by the derivative
01151              of -log10(x) when x=bin_tail. that is:
01152              tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
01153              Finally, we truncate the tail if the error is less than:
01154              tolerance * abs(-log10(bin_tail)-logNT) * bin_tail        */
01155           if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;
01156         }
01157     }
01158   return -log10(bin_tail) - logNT;
01159 }

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 1183 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().

01184 {
01185   /* check parameters */
01186   if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");
01187 
01188   /* copy values */
01189   out->x1 = in->x1;
01190   out->y1 = in->y1;
01191   out->x2 = in->x2;
01192   out->y2 = in->y2;
01193   out->width = in->width;
01194   out->x = in->x;
01195   out->y = in->y;
01196   out->theta = in->theta;
01197   out->dx = in->dx;
01198   out->dy = in->dy;
01199   out->prec = in->prec;
01200   out->p = in->p;
01201 }

Here is the call graph for this function:

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

Try some rectangles variations to improve NFA value.

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

Definition at line 1756 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().

01758 {
01759   struct rect r;
01760   double log_nfa,log_nfa_new;
01761   double delta = 0.5;
01762   double delta_2 = delta / 2.0;
01763   int n;
01764 
01765   log_nfa = rect_nfa(rec,angles,logNT);
01766 
01767   if( log_nfa > log_eps ) return log_nfa;
01768 
01769   /* try finer precisions */
01770   rect_copy(rec,&r);
01771   for(n=0; n<5; n++)
01772     {
01773       r.p /= 2.0;
01774       r.prec = r.p * M_PI;
01775       log_nfa_new = rect_nfa(&r,angles,logNT);
01776       if( log_nfa_new > log_nfa )
01777         {
01778           log_nfa = log_nfa_new;
01779           rect_copy(&r,rec);
01780         }
01781     }
01782 
01783   if( log_nfa > log_eps ) return log_nfa;
01784 
01785   /* try to reduce width */
01786   rect_copy(rec,&r);
01787   for(n=0; n<5; n++)
01788     {
01789       if( (r.width - delta) >= 0.5 )
01790         {
01791           r.width -= delta;
01792           log_nfa_new = rect_nfa(&r,angles,logNT);
01793           if( log_nfa_new > log_nfa )
01794             {
01795               rect_copy(&r,rec);
01796               log_nfa = log_nfa_new;
01797             }
01798         }
01799     }
01800 
01801   if( log_nfa > log_eps ) return log_nfa;
01802 
01803   /* try to reduce one side of the rectangle */
01804   rect_copy(rec,&r);
01805   for(n=0; n<5; n++)
01806     {
01807       if( (r.width - delta) >= 0.5 )
01808         {
01809           r.x1 += -r.dy * delta_2;
01810           r.y1 +=  r.dx * delta_2;
01811           r.x2 += -r.dy * delta_2;
01812           r.y2 +=  r.dx * delta_2;
01813           r.width -= delta;
01814           log_nfa_new = rect_nfa(&r,angles,logNT);
01815           if( log_nfa_new > log_nfa )
01816             {
01817               rect_copy(&r,rec);
01818               log_nfa = log_nfa_new;
01819             }
01820         }
01821     }
01822 
01823   if( log_nfa > log_eps ) return log_nfa;
01824 
01825   /* try to reduce the other side of the rectangle */
01826   rect_copy(rec,&r);
01827   for(n=0; n<5; n++)
01828     {
01829       if( (r.width - delta) >= 0.5 )
01830         {
01831           r.x1 -= -r.dy * delta_2;
01832           r.y1 -=  r.dx * delta_2;
01833           r.x2 -= -r.dy * delta_2;
01834           r.y2 -=  r.dx * delta_2;
01835           r.width -= delta;
01836           log_nfa_new = rect_nfa(&r,angles,logNT);
01837           if( log_nfa_new > log_nfa )
01838             {
01839               rect_copy(&r,rec);
01840               log_nfa = log_nfa_new;
01841             }
01842         }
01843     }
01844 
01845   if( log_nfa > log_eps ) return log_nfa;
01846 
01847   /* try even finer precisions */
01848   rect_copy(rec,&r);
01849   for(n=0; n<5; n++)
01850     {
01851       r.p /= 2.0;
01852       r.prec = r.p * M_PI;
01853       log_nfa_new = rect_nfa(&r,angles,logNT);
01854       if( log_nfa_new > log_nfa )
01855         {
01856           log_nfa = log_nfa_new;
01857           rect_copy(&r,rec);
01858         }
01859     }
01860 
01861   return log_nfa;
01862 }

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 1482 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().

01483 {
01484   rect_iter * i;
01485   int pts = 0;
01486   int alg = 0;
01487 
01488   /* check parameters */
01489   if( rec == NULL ) error("rect_nfa: invalid rectangle.");
01490   if( angles == NULL ) error("rect_nfa: invalid 'angles'.");
01491 
01492   /* compute the total number of pixels and of aligned points in 'rec' */
01493   for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */
01494     if( i->x >= 0 && i->y >= 0 &&
01495         i->x < (int) angles->xsize && i->y < (int) angles->ysize )
01496       {
01497         ++pts; /* total number of pixels counter */
01498         if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )
01499           ++alg; /* aligned points counter */
01500       }
01501   ri_del(i); /* delete iterator */
01502 
01503   return nfa(pts,alg,rec->p,logNT); /* compute NFA value */
01504 }

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 1869 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().

01874 {
01875   double density,rad1,rad2,rad,xc,yc;
01876   int i;
01877 
01878   /* check parameters */
01879   if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");
01880   if( reg_size == NULL )
01881     error("reduce_region_radius: invalid pointer 'reg_size'.");
01882   if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");
01883   if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");
01884   if( used == NULL || used->data == NULL )
01885     error("reduce_region_radius: invalid image 'used'.");
01886   if( angles == NULL || angles->data == NULL )
01887     error("reduce_region_radius: invalid image 'angles'.");
01888 
01889   /* compute region points density */
01890   density = (double) *reg_size /
01891                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01892 
01893   /* if the density criterion is satisfied there is nothing to do */
01894   if( density >= density_th ) return TRUE;
01895 
01896   /* compute region's radius */
01897   xc = (double) reg[0].x;
01898   yc = (double) reg[0].y;
01899   rad1 = dist( xc, yc, rec->x1, rec->y1 );
01900   rad2 = dist( xc, yc, rec->x2, rec->y2 );
01901   rad = rad1 > rad2 ? rad1 : rad2;
01902 
01903   /* while the density criterion is not satisfied, remove farther pixels */
01904   while( density < density_th )
01905     {
01906       rad *= 0.75; /* reduce region's radius to 75% of its value */
01907 
01908       /* remove points from the region and update 'used' map */
01909       for(i=0; i<*reg_size; i++)
01910         if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )
01911           {
01912             /* point not kept, mark it as NOTUSED */
01913             used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01914             /* remove point from the region */
01915             reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */
01916             reg[i].y = reg[*reg_size-1].y;
01917             --(*reg_size);
01918             --i; /* to avoid skipping one point */
01919           }
01920 
01921       /* reject if the region is too small.
01922          2 is the minimal region size for 'region2rect' to work. */
01923       if( *reg_size < 2 ) return FALSE;
01924 
01925       /* re-compute rectangle */
01926       region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01927 
01928       /* re-compute region points density */
01929       density = (double) *reg_size /
01930                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01931     }
01932 
01933   /* if this point is reached, the density criterion is satisfied */
01934   return TRUE;
01935 }

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 1947 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().

01950 {
01951   double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;
01952   int i,n;
01953 
01954   /* check parameters */
01955   if( reg == NULL ) error("refine: invalid pointer 'reg'.");
01956   if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");
01957   if( prec < 0.0 ) error("refine: 'prec' must be positive.");
01958   if( rec == NULL ) error("refine: invalid pointer 'rec'.");
01959   if( used == NULL || used->data == NULL )
01960     error("refine: invalid image 'used'.");
01961   if( angles == NULL || angles->data == NULL )
01962     error("refine: invalid image 'angles'.");
01963 
01964   /* compute region points density */
01965   density = (double) *reg_size /
01966                          ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01967 
01968   /* if the density criterion is satisfied there is nothing to do */
01969   if( density >= density_th ) return TRUE;
01970 
01971   /*------ First try: reduce angle tolerance ------*/
01972 
01973   /* compute the new mean angle and tolerance */
01974   xc = (double) reg[0].x;
01975   yc = (double) reg[0].y;
01976   ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];
01977   sum = s_sum = 0.0;
01978   n = 0;
01979   for(i=0; i<*reg_size; i++)
01980     {
01981       used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01982       if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width )
01983         {
01984           angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];
01985           ang_d = angle_diff_signed(angle,ang_c);
01986           sum += ang_d;
01987           s_sum += ang_d * ang_d;
01988           ++n;
01989         }
01990     }
01991   mean_angle = sum / (double) n;
01992   tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n
01993                          + mean_angle*mean_angle ); /* 2 * standard deviation */
01994 
01995   /* find a new region from the same starting point and new angle tolerance */
01996   region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,&reg_angle,used,tau);
01997 
01998   /* if the region is too small, reject */
01999   if( *reg_size < 2 ) return FALSE;
02000 
02001   /* re-compute rectangle */
02002   region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
02003 
02004   /* re-compute region points density */
02005   density = (double) *reg_size /
02006                       ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
02007 
02008   /*------ Second try: reduce region radius ------*/
02009   if( density < density_th )
02010     return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,
02011                                  rec, used, angles, density_th );
02012 
02013   /* if this point is reached, the density criterion is satisfied */
02014   return TRUE;
02015 }

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 1611 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().

01614 {
01615   double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;
01616   int i;
01617 
01618   /* check parameters */
01619   if( reg == NULL ) error("region2rect: invalid region.");
01620   if( reg_size <= 1 ) error("region2rect: region size <= 1.");
01621   if( modgrad == NULL || modgrad->data == NULL )
01622     error("region2rect: invalid image 'modgrad'.");
01623   if( rec == NULL ) error("region2rect: invalid 'rec'.");
01624 
01625   /* center of the region:
01626 
01627      It is computed as the weighted sum of the coordinates
01628      of all the pixels in the region. The norm of the gradient
01629      is used as the weight of a pixel. The sum is as follows:
01630        cx = \sum_i G(i).x_i
01631        cy = \sum_i G(i).y_i
01632      where G(i) is the norm of the gradient of pixel i
01633      and x_i,y_i are its coordinates.
01634    */
01635   x = y = sum = 0.0;
01636   for(i=0; i<reg_size; i++)
01637     {
01638       weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01639       x += (double) reg[i].x * weight;
01640       y += (double) reg[i].y * weight;
01641       sum += weight;
01642     }
01643   if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");
01644   x /= sum;
01645   y /= sum;
01646 
01647   /* theta */
01648   theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);
01649 
01650   /* length and width:
01651 
01652      'l' and 'w' are computed as the distance from the center of the
01653      region to pixel i, projected along the rectangle axis (dx,dy) and
01654      to the orthogonal axis (-dy,dx), respectively.
01655 
01656      The length of the rectangle goes from l_min to l_max, where l_min
01657      and l_max are the minimum and maximum values of l in the region.
01658      Analogously, the width is selected from w_min to w_max, where
01659      w_min and w_max are the minimum and maximum of w for the pixels
01660      in the region.
01661    */
01662   dx = cos(theta);
01663   dy = sin(theta);
01664   l_min = l_max = w_min = w_max = 0.0;
01665   for(i=0; i<reg_size; i++)
01666     {
01667       l =  ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;
01668       w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;
01669 
01670       if( l > l_max ) l_max = l;
01671       if( l < l_min ) l_min = l;
01672       if( w > w_max ) w_max = w;
01673       if( w < w_min ) w_min = w;
01674     }
01675 
01676   /* store values */
01677   rec->x1 = x + l_min * dx;
01678   rec->y1 = y + l_min * dy;
01679   rec->x2 = x + l_max * dx;
01680   rec->y2 = y + l_max * dy;
01681   rec->width = w_max - w_min;
01682   rec->x = x;
01683   rec->y = y;
01684   rec->theta = theta;
01685   rec->dx = dx;
01686   rec->dy = dy;
01687   rec->prec = prec;
01688   rec->p = p;
01689 
01690   /* we impose a minimal width of one pixel
01691 
01692      A sharp horizontal or vertical step would produce a perfectly
01693      horizontal or vertical region. The width computed would be
01694      zero. But that corresponds to a one pixels width transition in
01695      the image.
01696    */
01697   if( rec->width < 1.0 ) rec->width = 1.0;
01698 }

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 1704 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().

01707 {
01708   double sumdx,sumdy;
01709   int xx,yy,i;
01710 
01711   /* check parameters */
01712   if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
01713     error("region_grow: (x,y) out of the image.");
01714   if( angles == NULL || angles->data == NULL )
01715     error("region_grow: invalid image 'angles'.");
01716   if( reg == NULL ) error("region_grow: invalid 'reg'.");
01717   if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");
01718   if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");
01719   if( used == NULL || used->data == NULL )
01720     error("region_grow: invalid image 'used'.");
01721 
01722   /* first point of the region */
01723   *reg_size = 1;
01724   reg[0].x = x;
01725   reg[0].y = y;
01726   *reg_angle = angles->data[x+y*angles->xsize];  /* region's angle */
01727   sumdx = cos(*reg_angle);
01728   sumdy = sin(*reg_angle);
01729   used->data[x+y*used->xsize] = USED;
01730 
01731   /* try neighbors as new region points */
01732   for(i=0; i<*reg_size; i++)
01733     for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)
01734       for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)
01735         if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&
01736             used->data[xx+yy*used->xsize] != USED &&
01737             isaligned(xx,yy,angles,*reg_angle,prec) )
01738           {
01739             /* add point */
01740             used->data[xx+yy*used->xsize] = USED;
01741             reg[*reg_size].x = xx;
01742             reg[*reg_size].y = yy;
01743             ++(*reg_size);
01744 
01745             /* update region's angle */
01746             sumdx += cos( angles->data[xx+yy*angles->xsize] );
01747             sumdy += sin( angles->data[xx+yy*angles->xsize] );
01748             *reg_angle = atan2(sumdy,sumdx);
01749           }
01750 }

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 1314 of file lsd.c.

References error().

Referenced by rect_nfa().

01315 {
01316   if( iter == NULL ) error("ri_del: NULL iterator.");
01317   free( (void *) iter );
01318 }

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 1325 of file lsd.c.

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

Referenced by rect_nfa(), and ri_inc().

01326 {
01327   /* check input */
01328   if( i == NULL ) error("ri_end: NULL iterator.");
01329 
01330   /* if the current x value is larger than the largest
01331      x value in the rectangle (vx[2]), we know the full
01332      exploration of the rectangle is finished. */
01333   return (double)(i->x) > i->vx[2];
01334 }

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 1341 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().

01342 {
01343   /* check input */
01344   if( i == NULL ) error("ri_inc: NULL iterator.");
01345 
01346   /* if not at end of exploration,
01347      increase y value for next pixel in the 'column' */
01348   if( !ri_end(i) ) i->y++;
01349 
01350   /* if the end of the current 'column' is reached,
01351      and it is not the end of exploration,
01352      advance to the next 'column' */
01353   while( (double) (i->y) > i->ye && !ri_end(i) )
01354     {
01355       /* increase x, next 'column' */
01356       i->x++;
01357 
01358       /* if end of exploration, return */
01359       if( ri_end(i) ) return;
01360 
01361       /* update lower y limit (start) for the new 'column'.
01362 
01363          We need to interpolate the y value that corresponds to the
01364          lower side of the rectangle. The first thing is to decide if
01365          the corresponding side is
01366 
01367            vx[0],vy[0] to vx[3],vy[3] or
01368            vx[3],vy[3] to vx[2],vy[2]
01369 
01370          Then, the side is interpolated for the x value of the
01371          'column'. But, if the side is vertical (as it could happen if
01372          the rectangle is vertical and we are dealing with the first
01373          or last 'columns') then we pick the lower value of the side
01374          by using 'inter_low'.
01375        */
01376       if( (double) i->x < i->vx[3] )
01377         i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);
01378       else
01379         i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);
01380 
01381       /* update upper y limit (end) for the new 'column'.
01382 
01383          We need to interpolate the y value that corresponds to the
01384          upper side of the rectangle. The first thing is to decide if
01385          the corresponding side is
01386 
01387            vx[0],vy[0] to vx[1],vy[1] or
01388            vx[1],vy[1] to vx[2],vy[2]
01389 
01390          Then, the side is interpolated for the x value of the
01391          'column'. But, if the side is vertical (as it could happen if
01392          the rectangle is vertical and we are dealing with the first
01393          or last 'columns') then we pick the lower value of the side
01394          by using 'inter_low'.
01395        */
01396       if( (double)i->x < i->vx[1] )
01397         i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);
01398       else
01399         i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);
01400 
01401       /* new y */
01402       i->y = (int) ceil(i->ys);
01403     }
01404 }

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 1411 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().

01412 {
01413   double vx[4],vy[4];
01414   int n,offset;
01415   rect_iter * i;
01416 
01417   /* check parameters */
01418   if( r == NULL ) error("ri_ini: invalid rectangle.");
01419 
01420   /* get memory */
01421   i = (rect_iter *) malloc(sizeof(rect_iter));
01422   if( i == NULL ) error("ri_ini: Not enough memory.");
01423 
01424   /* build list of rectangle corners ordered
01425      in a circular way around the rectangle */
01426   vx[0] = r->x1 - r->dy * r->width / 2.0;
01427   vy[0] = r->y1 + r->dx * r->width / 2.0;
01428   vx[1] = r->x2 - r->dy * r->width / 2.0;
01429   vy[1] = r->y2 + r->dx * r->width / 2.0;
01430   vx[2] = r->x2 + r->dy * r->width / 2.0;
01431   vy[2] = r->y2 - r->dx * r->width / 2.0;
01432   vx[3] = r->x1 + r->dy * r->width / 2.0;
01433   vy[3] = r->y1 - r->dx * r->width / 2.0;
01434 
01435   /* compute rotation of index of corners needed so that the first
01436      point has the smaller x.
01437 
01438      if one side is vertical, thus two corners have the same smaller x
01439      value, the one with the largest y value is selected as the first.
01440    */
01441   if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;
01442   else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;
01443   else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;
01444   else offset = 3;
01445 
01446   /* apply rotation of index. */
01447   for(n=0; n<4; n++)
01448     {
01449       i->vx[n] = vx[(offset+n)%4];
01450       i->vy[n] = vy[(offset+n)%4];
01451     }
01452 
01453   /* Set an initial condition.
01454 
01455      The values are set to values that will cause 'ri_inc' (that will
01456      be called immediately) to initialize correctly the first 'column'
01457      and compute the limits 'ys' and 'ye'.
01458 
01459      'y' is set to the integer value of vy[0], the starting corner.
01460 
01461      'ys' and 'ye' are set to very small values, so 'ri_inc' will
01462      notice that it needs to start a new 'column'.
01463 
01464      The smallest integer coordinate inside of the rectangle is
01465      'ceil(vx[0])'. The current 'x' value is set to that value minus
01466      one, so 'ri_inc' (that will increase x by one) will advance to
01467      the first 'column'.
01468    */
01469   i->x = (int) ceil(i->vx[0]) - 1;
01470   i->y = (int) ceil(i->vy[0]);
01471   i->ys = i->ye = -DBL_MAX;
01472 
01473   /* advance to the first pixel */
01474   ri_inc(i);
01475 
01476   return i;
01477 }

Here is the call graph for this function:


Generated on Fri Nov 11 11:11:11 2011 for LSD by doxygen 1.3.4