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:
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_s * | ntuple_list |
'list of n-tuple' data type | |
typedef image_char_s * | image_char |
char image data type | |
typedef image_int_s * | image_int |
int image data type | |
typedef image_double_s * | image_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_iter * | ri_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. |
|
Definition at line 115 of file lsd.c. Referenced by isaligned(), reduce_region_radius(), and refine(). |
|
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(). |
|
2 pi
Definition at line 129 of file lsd.c. Referenced by angle_diff(), angle_diff_signed(), and isaligned(). |
|
3/2 pi
Definition at line 126 of file lsd.c. Referenced by isaligned(). |
|
ln(10)
Definition at line 106 of file lsd.c. Referenced by nfa(). |
|
PI.
Definition at line 111 of file lsd.c. Referenced by angle_diff(), angle_diff_signed(), get_theta(), LineSegmentDetection(), and rect_improve(). |
|
Label for pixels with undefined gradient.
Definition at line 123 of file lsd.c. Referenced by isaligned(), LineSegmentDetection(), and ll_angle(). |
|
Label for pixels not used in yet.
Definition at line 132 of file lsd.c. Referenced by LineSegmentDetection(), reduce_region_radius(), and refine(). |
|
Doubles relative error factor.
Definition at line 168 of file lsd.c. Referenced by double_equal(). |
|
Size of the table to store already computed inverse values.
Definition at line 1030 of file lsd.c. Referenced by nfa(). |
|
Definition at line 119 of file lsd.c. Referenced by double_equal(), reduce_region_radius(), and refine(). |
|
Label for pixels already used in detection.
Definition at line 135 of file lsd.c. Referenced by region_grow(). |
|
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(). |
|
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(). |
|
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(). |
|
'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(). |
|
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:
|
Absolute value angle difference.
Definition at line 931 of file lsd.c. Referenced by get_theta().
|
|
Signed angle difference.
Definition at line 943 of file lsd.c. Referenced by refine().
|
|
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 }
|
|
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 } |
|
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:
|
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 }
|
|
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().
|
Here is the call graph for this function:
|
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().
|
Here is the call graph for this function:
|
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().
|
Here is the call graph for this function:
|
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:
|
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
before the sub-sampling to prevent aliasing. The standard deviation sigma given by:
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
where
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:
|
Compute region's angle as the principal inertia axis of the region. The following is the region inertia matrix A:
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
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:
|
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:
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:
|
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:
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:
|
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:
|
LSD Full Interface.
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, ®_size, 02117 ®_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, ®_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:
|
Computes the direction of the level line of 'in' at each point. The result is:
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:
|
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
so
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 } |
|
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
so
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 }
|
|
LSD Simple Interface.
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:
|
LSD Simple Interface with Scale.
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:
|
LSD Simple Interface with Scale and Region output.
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:
|
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:
|
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:
|
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:
|
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:
|
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:
|
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:
|
Create an n-tuple list and allocate memory for one element.
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:
|
Computes -log10(NFA). NFA stands for Number of False Alarms:
The value -log10(NFA) is equivalent but more intuitive than NFA:
Used this way, the bigger the value, better the detection, and a logarithmic scale is used.
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:
|
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:
|
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:
|
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:
|
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:
|
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,®_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:
|
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:
|
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:
|
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:
|
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:
|
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:
|
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: