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 | point |
A point (or pixel). More... | |
struct | rect |
Rectangle structure: line segment with width. More... | |
struct | rect_iter |
Rectangle points iterator. More... | |
Defines | |
#define | M_LN10 2.30258509299404568402 |
ln(10) | |
#define | M_PI 3.14159265358979323846 |
PI. | |
#define | FALSE 0 |
#define | TRUE 1 |
#define | NOTDEF -1024.0 |
Label for pixels with undefined gradient. | |
#define | M_3_2_PI 4.71238898038 |
3/2 pi | |
#define | M_2__PI 6.28318530718 |
2 pi | |
#define | NOTUSED 0 |
Label for pixels not used in yet. | |
#define | USED 1 |
Label for pixels already used in detection. | |
#define | RELATIVE_ERROR_FACTOR 100.0 |
Doubles relative error factor. | |
#define | log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x)) |
Computes the natural logarithm of the absolute value of the gamma function of x. | |
#define | TABSIZE 100000 |
Size of the table to store already computed inverse values. | |
Functions | |
void | error (char *msg) |
Fatal error, print a message to standard-error output and exit. | |
int | double_equal (double a, double b) |
Compare doubles by relative error. | |
double | dist (double x1, double y1, double x2, double y2) |
Computes Euclidean distance between point (x1,y1) and point (x2,y2). | |
void | free_ntuple_list (ntuple_list in) |
Free memory used in n-tuple 'in'. | |
ntuple_list | new_ntuple_list (unsigned int dim) |
Create an n-tuple list and allocate memory for one element. | |
void | enlarge_ntuple_list (ntuple_list n_tuple) |
Enlarge the allocated memory of an n-tuple list. | |
void | add_5tuple (ntuple_list out, double v1, double v2, double v3, double v4, double v5) |
Add a 5-tuple to an n-tuple list. | |
void | free_image_char (image_char i) |
Free memory used in image_char 'i'. | |
image_char | new_image_char (unsigned int xsize, unsigned int ysize) |
Create a new image_char of size 'xsize' times 'ysize'. | |
image_char | new_image_char_ini (unsigned int xsize, unsigned int ysize, unsigned char fill_value) |
Create a new image_char of size 'xsize' times 'ysize', initialized to the value 'fill_value'. | |
void | free_image_int (image_int i) |
Free memory used in image_int 'i'. | |
image_int | new_image_int (unsigned int xsize, unsigned int ysize) |
Create a new image_int of size 'xsize' times 'ysize'. | |
image_int | new_image_int_ini (unsigned int xsize, unsigned int ysize, int fill_value) |
Create a new image_int of size 'xsize' times 'ysize', initialized to the value 'fill_value'. | |
void | free_image_double (image_double i) |
Free memory used in image_double 'i'. | |
image_double | new_image_double (unsigned int xsize, unsigned int ysize) |
Create a new image_double of size 'xsize' times 'ysize'. | |
image_double | new_image_double_ini (unsigned int xsize, unsigned int ysize, double fill_value) |
Create a new image_double of size 'xsize' times 'ysize', initialized to the value 'fill_value'. | |
void | gaussian_kernel (ntuple_list kernel, double sigma, double mean) |
Compute a Gaussian kernel of length 'kernel->dim', standard deviation 'sigma', and centered at value 'mean'. | |
image_double | gaussian_sampler (image_double in, double scale, double sigma_scale) |
Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling. | |
image_double | ll_angle (image_double in, double threshold, struct coorlist **list_p, void **mem_p, image_double *modgrad, unsigned int n_bins, double max_grad) |
Computes the direction of the level line of 'in' at each point. | |
int | isaligned (int x, int y, image_double angles, double theta, double prec) |
Is point (x,y) aligned to angle theta, up to precision 'prec'? | |
double | angle_diff (double a, double b) |
Absolute value angle difference. | |
double | angle_diff_signed (double a, double b) |
Signed angle difference. | |
double | log_gamma_lanczos (double x) |
Computes the natural logarithm of the absolute value of the gamma function of x using the Lanczos approximation. | |
double | log_gamma_windschitl (double x) |
Computes the natural logarithm of the absolute value of the gamma function of x using Windschitl method. | |
double | nfa (int n, int k, double p, double logNT) |
Computes -log10(NFA). | |
void | rect_copy (struct rect *in, struct rect *out) |
Copy one rectangle structure to another. | |
double | inter_low (double x, double x1, double y1, double x2, double y2) |
Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smaller of 'y1' and 'y2'. | |
double | inter_hi (double x, double x1, double y1, double x2, double y2) |
Interpolate y value corresponding to 'x' value given, in the line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the larger of 'y1' and 'y2'. | |
void | ri_del (rect_iter *iter) |
Free memory used by a rectangle iterator. | |
int | ri_end (rect_iter *i) |
Check if the iterator finished the full iteration. | |
void | ri_inc (rect_iter *i) |
Increment a rectangle iterator. | |
rect_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 eps) |
Try some rectangles variations to improve NFA value. | |
int | reduce_region_radius (struct point *reg, int *reg_size, image_double modgrad, double reg_angle, double prec, double p, struct rect *rec, image_char used, image_double angles, double density_th) |
Reduce the region size, by elimination the points far from the starting point, until that leads to rectangle with the right density of region points or to discard the region if too small. | |
int | refine (struct point *reg, int *reg_size, image_double modgrad, double reg_angle, double prec, double p, struct rect *rec, image_char used, image_double angles, double density_th) |
Refine a rectangle. | |
ntuple_list | LineSegmentDetection (image_double image, double scale, double sigma_scale, double quant, double ang_th, double eps, double density_th, int n_bins, double max_grad, image_int *region) |
LSD Full Interface. | |
ntuple_list | lsd_scale (image_double image, double scale) |
LSD Simple Interface with Scale. | |
ntuple_list | lsd (image_double image) |
LSD Simple Interface. |
|
Definition at line 89 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 913 of file lsd.c. Referenced by nfa(). |
|
2 pi
Definition at line 103 of file lsd.c. Referenced by angle_diff(), angle_diff_signed(), and isaligned(). |
|
3/2 pi
Definition at line 100 of file lsd.c. Referenced by isaligned(). |
|
ln(10)
Definition at line 80 of file lsd.c. Referenced by nfa(). |
|
PI.
Definition at line 85 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 97 of file lsd.c. Referenced by isaligned(), LineSegmentDetection(), and ll_angle(). |
|
Label for pixels not used in yet.
Definition at line 106 of file lsd.c. Referenced by LineSegmentDetection(), reduce_region_radius(), and refine(). |
|
Doubles relative error factor.
Definition at line 142 of file lsd.c. Referenced by double_equal(). |
|
Size of the table to store already computed inverse values.
Definition at line 918 of file lsd.c. Referenced by nfa(). |
|
Definition at line 93 of file lsd.c. Referenced by double_equal(), reduce_region_radius(), and refine(). |
|
Label for pixels already used in detection.
Definition at line 109 of file lsd.c. Referenced by region_grow(). |
|
Add a 5-tuple to an n-tuple list.
Definition at line 250 of file lsd.c. References ntuple_list_s::dim, enlarge_ntuple_list(), error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values. Referenced by LineSegmentDetection().
00252 { 00253 /* check parameters */ 00254 if( out == NULL ) error("add_5tuple: invalid n-tuple input."); 00255 if( out->dim != 5 ) error("add_5tuple: the n-tuple must be a 5-tuple."); 00256 00257 /* if needed, alloc more tuples to 'out' */ 00258 if( out->size == out->max_size ) enlarge_ntuple_list(out); 00259 if( out->values == NULL ) error("add_5tuple: invalid n-tuple input."); 00260 00261 /* add new 5-tuple */ 00262 out->values[ out->size * out->dim + 0 ] = v1; 00263 out->values[ out->size * out->dim + 1 ] = v2; 00264 out->values[ out->size * out->dim + 2 ] = v3; 00265 out->values[ out->size * out->dim + 3 ] = v4; 00266 out->values[ out->size * out->dim + 4 ] = v5; 00267 00268 /* update number of tuples counter */ 00269 out->size++; 00270 } |
Here is the call graph for this function:
|
Absolute value angle difference.
Definition at line 819 of file lsd.c. Referenced by get_theta().
|
|
Signed angle difference.
Definition at line 831 of file lsd.c. Referenced by refine().
|
|
Computes Euclidean distance between point (x1,y1) and point (x2,y2).
Definition at line 181 of file lsd.c. Referenced by reduce_region_radius(), and refine().
00182 {
00183 return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00184 }
|
|
Compare doubles by relative error. The resulting rounding error after floating point computations depend on the specific operations done. The same number computed by different algorithms could present different rounding errors. For a useful comparison, an estimation of the relative rounding error should be considered and compared to a factor times EPS. The factor should be related to the cumulated rounding error in the chain of computation. Here, as a simplification, a fixed factor is used. Definition at line 155 of file lsd.c. References RELATIVE_ERROR_FACTOR, and TRUE. Referenced by get_theta(), inter_hi(), inter_low(), and nfa().
00156 { 00157 double abs_diff,aa,bb,abs_max; 00158 00159 /* trivial case */ 00160 if( a == b ) return TRUE; 00161 00162 abs_diff = fabs(a-b); 00163 aa = fabs(a); 00164 bb = fabs(b); 00165 abs_max = aa > bb ? aa : bb; 00166 00167 /* DBL_MIN is the smallest normalized number, thus, the smallest 00168 number whose relative error is bounded by DBL_EPSILON. For 00169 smaller numbers, the same quantization steps as for DBL_MIN 00170 are used. Then, for smaller numbers, a meaningful "relative" 00171 error should be computed by dividing the difference by DBL_MIN. */ 00172 if( abs_max < DBL_MIN ) abs_max = DBL_MIN; 00173 00174 /* equal if relative error <= factor x eps */ 00175 return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON); 00176 } |
|
Enlarge the allocated memory of an n-tuple list.
Definition at line 232 of file lsd.c. References ntuple_list_s::dim, error(), ntuple_list_s::max_size, ntuple_list, and ntuple_list_s::values. Referenced by add_5tuple(), and gaussian_kernel().
00233 { 00234 /* check parameters */ 00235 if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 ) 00236 error("enlarge_ntuple_list: invalid n-tuple."); 00237 00238 /* duplicate number of tuples */ 00239 n_tuple->max_size *= 2; 00240 00241 /* realloc memory */ 00242 n_tuple->values = (double *) realloc( (void *) n_tuple->values, 00243 n_tuple->dim * n_tuple->max_size * sizeof(double) ); 00244 if( n_tuple->values == NULL ) error("not enough memory."); 00245 } |
Here is the call graph for this function:
|
Fatal error, print a message to standard-error output and exit.
Definition at line 133 of file lsd.c. Referenced by add_5tuple(), enlarge_ntuple_list(), free_image_char(), free_image_double(), free_image_int(), free_ntuple_list(), gaussian_kernel(), gaussian_sampler(), get_theta(), inter_hi(), inter_low(), isaligned(), LineSegmentDetection(), ll_angle(), new_image_char(), new_image_char_ini(), new_image_double(), new_image_int(), new_ntuple_list(), nfa(), rect_copy(), rect_nfa(), reduce_region_radius(), refine(), region2rect(), region_grow(), ri_del(), ri_end(), ri_inc(), and ri_ini().
00134 {
00135 fprintf(stderr,"LSD Error: %s\n",msg);
00136 exit(EXIT_FAILURE);
00137 }
|
|
Free memory used in image_char 'i'.
Definition at line 280 of file lsd.c. References image_char_s::data, error(), and image_char. Referenced by LineSegmentDetection().
|
Here is the call graph for this function:
|
Free memory used in image_double 'i'.
Definition at line 387 of file lsd.c. References image_double_s::data, error(), and image_double. Referenced by gaussian_sampler(), and LineSegmentDetection().
|
Here is the call graph for this function:
|
Free memory used in image_int 'i'.
Definition at line 336 of file lsd.c. References image_int_s::data, error(), and image_int.
|
Here is the call graph for this function:
|
Free memory used in n-tuple 'in'.
Definition at line 194 of file lsd.c. References error(), ntuple_list, and ntuple_list_s::values. Referenced by gaussian_sampler().
|
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 448 of file lsd.c. References ntuple_list_s::dim, enlarge_ntuple_list(), error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values. Referenced by gaussian_sampler().
00449 { 00450 double sum = 0.0; 00451 double val; 00452 unsigned int i; 00453 00454 /* check parameters */ 00455 if( kernel == NULL || kernel->values == NULL ) 00456 error("gaussian_kernel: invalid n-tuple 'kernel'."); 00457 if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive."); 00458 00459 /* compute Gaussian kernel */ 00460 if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel); 00461 kernel->size = 1; 00462 for(i=0;i<kernel->dim;i++) 00463 { 00464 val = ( (double) i - mean ) / sigma; 00465 kernel->values[i] = exp( -0.5 * val * val ); 00466 sum += kernel->values[i]; 00467 } 00468 00469 /* normalization */ 00470 if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum; 00471 } |
Here is the call graph for this function:
|
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 apply a combined Gaussian kernel and sampling in the x axis, and then the combined Gaussian kernel and sampling in the y axis. Definition at line 511 of file lsd.c. References image_double_s::data, ntuple_list_s::dim, error(), free_image_double(), free_ntuple_list(), gaussian_kernel(), image_double, new_image_double(), new_ntuple_list(), ntuple_list, ntuple_list_s::values, image_double_s::xsize, and image_double_s::ysize. Referenced by LineSegmentDetection().
00513 { 00514 image_double aux,out; 00515 ntuple_list kernel; 00516 unsigned int N,M,h,n,x,y,i; 00517 int xc,yc,j,double_x_size,double_y_size; 00518 double sigma,xx,yy,sum,prec; 00519 00520 /* check parameters */ 00521 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 ) 00522 error("gaussian_sampler: invalid image."); 00523 if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive."); 00524 if( sigma_scale <= 0.0 ) 00525 error("gaussian_sampler: 'sigma_scale' must be positive."); 00526 00527 /* get memory for images */ 00528 if( in->xsize * scale > (double) UINT_MAX || 00529 in->ysize * scale > (double) UINT_MAX ) 00530 error("gaussian_sampler: the output image size exceeds the handled size."); 00531 N = (unsigned int) floor( in->xsize * scale ); 00532 M = (unsigned int) floor( in->ysize * scale ); 00533 aux = new_image_double(N,in->ysize); 00534 out = new_image_double(N,M); 00535 00536 /* sigma, kernel size and memory for the kernel */ 00537 sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale; 00538 /* 00539 The size of the kernel is selected to guarantee that the 00540 the first discarded term is at least 10^prec times smaller 00541 than the central value. For that, h should be larger than x, with 00542 e^(-x^2/2sigma^2) = 1/10^prec. 00543 Then, 00544 x = sigma * sqrt( 2 * prec * ln(10) ). 00545 */ 00546 prec = 3.0; 00547 h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) ); 00548 n = 1+2*h; /* kernel size */ 00549 kernel = new_ntuple_list(n); 00550 00551 /* auxiliary double image size variables */ 00552 double_x_size = (int) (2 * in->xsize); 00553 double_y_size = (int) (2 * in->ysize); 00554 00555 /* First subsampling: x axis */ 00556 for(x=0;x<aux->xsize;x++) 00557 { 00558 /* 00559 x is the coordinate in the new image. 00560 xx is the corresponding x-value in the original size image. 00561 xc is the integer value, the pixel coordinate of xx. 00562 */ 00563 xx = (double) x / scale; 00564 /* coordinate (0.0,0.0) is in the center of pixel (0,0), 00565 so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */ 00566 xc = (int) floor( xx + 0.5 ); 00567 gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc ); 00568 /* the kernel must be computed for each x because the fine 00569 offset xx-xc is different in each case */ 00570 00571 for(y=0;y<aux->ysize;y++) 00572 { 00573 sum = 0.0; 00574 for(i=0;i<kernel->dim;i++) 00575 { 00576 j = xc - h + i; 00577 00578 /* symmetry boundary condition */ 00579 while( j < 0 ) j += double_x_size; 00580 while( j >= double_x_size ) j -= double_x_size; 00581 if( j >= (int) in->xsize ) j = double_x_size-1-j; 00582 00583 sum += in->data[ j + y * in->xsize ] * kernel->values[i]; 00584 } 00585 aux->data[ x + y * aux->xsize ] = sum; 00586 } 00587 } 00588 00589 /* Second subsampling: y axis */ 00590 for(y=0;y<out->ysize;y++) 00591 { 00592 /* 00593 y is the coordinate in the new image. 00594 yy is the corresponding x-value in the original size image. 00595 yc is the integer value, the pixel coordinate of xx. 00596 */ 00597 yy = (double) y / scale; 00598 /* coordinate (0.0,0.0) is in the center of pixel (0,0), 00599 so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */ 00600 yc = (int) floor( yy + 0.5 ); 00601 gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc ); 00602 /* the kernel must be computed for each y because the fine 00603 offset yy-yc is different in each case */ 00604 00605 for(x=0;x<out->xsize;x++) 00606 { 00607 sum = 0.0; 00608 for(i=0;i<kernel->dim;i++) 00609 { 00610 j = yc - h + i; 00611 00612 /* symmetry boundary condition */ 00613 while( j < 0 ) j += double_y_size; 00614 while( j >= double_y_size ) j -= double_y_size; 00615 if( j >= (int) in->ysize ) j = double_y_size-1-j; 00616 00617 sum += aux->data[ x + j * aux->xsize ] * kernel->values[i]; 00618 } 00619 out->data[ x + y * out->xsize ] = sum; 00620 } 00621 } 00622 00623 /* free memory */ 00624 free_ntuple_list(kernel); 00625 free_image_double(aux); 00626 00627 return out; 00628 } |
Here is the call graph for this function:
|
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 assotiated to the smaller eigenvalue. We have to solve a,b in: a.Ixx + b.Ixy = a.lambda2 a.Ixy + b.Iyy = b.lambda2 We want the angle theta = atan(b/a). It can be computed with any of the two equations: theta = atan( (lambda2-Ixx) / Ixy ) or theta = atan( Ixy / (lambda2-Iyy) ) When |Ixx| > |Iyy| we use the first, otherwise the second (just to get better numeric precision). Definition at line 1456 of file lsd.c. References angle_diff(), image_double_s::data, double_equal(), error(), image_double, M_PI, point::x, image_double_s::xsize, and point::y. Referenced by region2rect().
01458 { 01459 double lambda,theta,weight; 01460 double Ixx = 0.0; 01461 double Iyy = 0.0; 01462 double Ixy = 0.0; 01463 int i; 01464 01465 /* check parameters */ 01466 if( reg == NULL ) error("get_theta: invalid region."); 01467 if( reg_size <= 1 ) error("get_theta: region size <= 1."); 01468 if( modgrad == NULL || modgrad->data == NULL ) 01469 error("get_theta: invalid 'modgrad'."); 01470 if( prec < 0.0 ) error("get_theta: 'prec' must be positive."); 01471 01472 /* compute inertia matrix */ 01473 for(i=0; i<reg_size; i++) 01474 { 01475 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ]; 01476 Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight; 01477 Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight; 01478 Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight; 01479 } 01480 if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) ) 01481 error("get_theta: null inertia matrix."); 01482 01483 /* compute smallest eigenvalue */ 01484 lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) ); 01485 01486 /* compute angle */ 01487 theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy); 01488 01489 /* The previous procedure don't cares about orientation, 01490 so it could be wrong by 180 degrees. Here is corrected if necessary. */ 01491 if( angle_diff(theta,reg_angle) > prec ) theta += M_PI; 01492 01493 return theta; 01494 } |
Here is the call graph for this function:
|
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 1187 of file lsd.c. References double_equal(), and error(). Referenced by ri_inc().
01188 { 01189 /* check parameters */ 01190 if( x1 > x2 || x < x1 || x > x2 ) 01191 error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'."); 01192 01193 /* interpolation */ 01194 if( double_equal(x1,x2) && y1<y2 ) return y2; 01195 if( double_equal(x1,x2) && y1>y2 ) return y1; 01196 return y1 + (x-x1) * (y2-y1) / (x2-x1); 01197 } |
Here is the call graph for this function:
|
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 1165 of file lsd.c. References double_equal(), and error(). Referenced by ri_inc().
01166 { 01167 /* check parameters */ 01168 if( x1 > x2 || x < x1 || x > x2 ) 01169 error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'."); 01170 01171 /* interpolation */ 01172 if( double_equal(x1,x2) && y1<y2 ) return y1; 01173 if( double_equal(x1,x2) && y1>y2 ) return y2; 01174 return y1 + (x-x1) * (y2-y1) / (x2-x1); 01175 } |
Here is the call graph for this function:
|
Is point (x,y) aligned to angle theta, up to precision 'prec'?
Definition at line 781 of file lsd.c. References image_double_s::data, error(), FALSE, image_double, M_2__PI, M_3_2_PI, NOTDEF, image_double_s::xsize, and image_double_s::ysize. Referenced by rect_nfa(), and region_grow().
00783 { 00784 double a; 00785 00786 /* check parameters */ 00787 if( angles == NULL || angles->data == NULL ) 00788 error("isaligned: invalid image 'angles'."); 00789 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize ) 00790 error("isaligned: (x,y) out of the image."); 00791 if( prec < 0.0 ) error("isaligned: 'prec' must be positive."); 00792 00793 /* angle at pixel (x,y) */ 00794 a = angles->data[ x + y * angles->xsize ]; 00795 00796 /* pixels whose level-line angle is not defined 00797 are considered as NON-aligned */ 00798 if( a == NOTDEF ) return FALSE; /* there is no need to call the function 00799 'double_equal' here because there is 00800 no risk of problems related to the 00801 comparison doubles, we are only 00802 interested in the exact NOTDEF value */ 00803 00804 /* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */ 00805 theta -= a; 00806 if( theta < 0.0 ) theta = -theta; 00807 if( theta > M_3_2_PI ) 00808 { 00809 theta -= M_2__PI; 00810 if( theta < 0.0 ) theta = -theta; 00811 } 00812 00813 return theta < prec; 00814 } |
Here is the call graph for this function:
|
LSD Full Interface.
Definition at line 1913 of file lsd.c. References add_5tuple(), image_char_s::data, image_double_s::data, error(), free_image_char(), free_image_double(), gaussian_sampler(), image_char, image_double, image_int, ll_angle(), M_PI, new_image_char_ini(), new_image_int_ini(), new_ntuple_list(), NOTDEF, NOTUSED, ntuple_list, rect_improve(), refine(), region2rect(), region_grow(), rect::width, rect::x1, rect::x2, image_char_s::xsize, image_double_s::xsize, rect::y1, rect::y2, and image_double_s::ysize. Referenced by lsd_scale().
01918 { 01919 ntuple_list out = new_ntuple_list(5); 01920 image_double scaled_image,angles,modgrad; 01921 image_char used; 01922 struct coorlist * list_p; 01923 void * mem_p; 01924 struct rect rec; 01925 struct point * reg; 01926 int reg_size,min_reg_size,i; 01927 unsigned int xsize,ysize; 01928 double rho,reg_angle,prec,p,log_nfa,logNT; 01929 int ls_count = 0; /* line segments are numbered 1,2,3,... */ 01930 01931 01932 /* check parameters */ 01933 if( image==NULL || image->data==NULL || image->xsize==0 || image->ysize==0 ) 01934 error("invalid image input."); 01935 if( scale <= 0.0 ) error("'scale' value must be positive."); 01936 if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive."); 01937 if( quant < 0.0 ) error("'quant' value must be positive."); 01938 if( ang_th <= 0.0 || ang_th >= 180.0 ) 01939 error("'ang_th' value must be in the range (0,180)."); 01940 if( density_th < 0.0 || density_th > 1.0 ) 01941 error("'density_th' value must be in the range [0,1]."); 01942 if( n_bins <= 0 ) error("'n_bins' value must be positive."); 01943 if( max_grad <= 0.0 ) error("'max_grad' value must be positive."); 01944 01945 01946 /* angle tolerance */ 01947 prec = M_PI * ang_th / 180.0; 01948 p = ang_th / 180.0; 01949 rho = quant / sin(prec); /* gradient magnitude threshold */ 01950 01951 01952 /* scale image (if necessary) and compute angle at each pixel */ 01953 if( scale != 1.0 ) 01954 { 01955 scaled_image = gaussian_sampler( image, scale, sigma_scale ); 01956 angles = ll_angle( scaled_image, rho, &list_p, &mem_p, 01957 &modgrad, (unsigned int) n_bins, max_grad ); 01958 free_image_double(scaled_image); 01959 } 01960 else 01961 angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad, 01962 (unsigned int) n_bins, max_grad ); 01963 xsize = angles->xsize; 01964 ysize = angles->ysize; 01965 logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0; 01966 min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region 01967 that can give a meaningful event */ 01968 01969 01970 /* initialize some structures */ 01971 if( region != NULL ) /* image to output pixel region number, if asked */ 01972 *region = new_image_int_ini(angles->xsize,angles->ysize,0); 01973 used = new_image_char_ini(xsize,ysize,NOTUSED); 01974 reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) ); 01975 if( reg == NULL ) error("not enough memory!"); 01976 01977 01978 /* search for line segments */ 01979 for(; list_p != NULL; list_p = list_p->next ) 01980 if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED && 01981 angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF ) 01982 /* there is no risk of double comparison problems here 01983 because we are only interested in the exact NOTDEF value */ 01984 { 01985 /* find the region of connected point and ~equal angle */ 01986 region_grow( list_p->x, list_p->y, angles, reg, ®_size, 01987 ®_angle, used, prec ); 01988 01989 /* reject small regions */ 01990 if( reg_size < min_reg_size ) continue; 01991 01992 /* construct rectangular approximation for the region */ 01993 region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec); 01994 01995 /* Check if the rectangle exceeds the minimal density of 01996 region points. If not, try to improve the region. 01997 The rectangle will be rejected if the final one does 01998 not fulfill the minimal density condition. 01999 This is an addition to the original LSD algorithm published in 02000 "LSD: A Fast Line Segment Detector with a False Detection Control" 02001 by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall. 02002 The original algorithm is obtained with density_th = 0.0. 02003 */ 02004 if( !refine( reg, ®_size, modgrad, reg_angle, 02005 prec, p, &rec, used, angles, density_th ) ) continue; 02006 02007 /* compute NFA value */ 02008 log_nfa = rect_improve(&rec,angles,logNT,eps); 02009 if( log_nfa <= eps ) continue; 02010 02011 /* A New Line Segment was found! */ 02012 ++ls_count; /* increase line segment counter */ 02013 02014 /* 02015 The gradient was computed with a 2x2 mask, its value corresponds to 02016 points with an offset of (0.5,0.5), that should be added to output. 02017 The coordinates origin is at the center of pixel (0,0). 02018 */ 02019 rec.x1 += 0.5; rec.y1 += 0.5; 02020 rec.x2 += 0.5; rec.y2 += 0.5; 02021 02022 /* scale the result values if a subsampling was performed */ 02023 if( scale != 1.0 ) 02024 { 02025 rec.x1 /= scale; rec.y1 /= scale; 02026 rec.x2 /= scale; rec.y2 /= scale; 02027 rec.width /= scale; 02028 } 02029 02030 /* add line segment found to output */ 02031 add_5tuple(out, rec.x1, rec.y1, rec.x2, rec.y2, rec.width); 02032 02033 /* add region number to 'region' image if needed */ 02034 if( region != NULL ) 02035 for(i=0; i<reg_size; i++) 02036 (*region)->data[reg[i].x+reg[i].y*(*region)->xsize] = ls_count; 02037 } 02038 02039 02040 /* free memory */ 02041 free_image_double(angles); 02042 free_image_double(modgrad); 02043 free_image_char(used); 02044 free( (void *) reg ); 02045 free( (void *) mem_p ); 02046 02047 return out; 02048 } |
Here is the call graph for this function:
|
Computes the direction of the level line of 'in' at each point. The result is:
Definition at line 652 of file lsd.c. References image_double_s::data, error(), image_double, new_image_double(), NOTDEF, image_double_s::xsize, and image_double_s::ysize. Referenced by LineSegmentDetection().
00656 { 00657 image_double g; 00658 unsigned int n,p,x,y,adr,i; 00659 double com1,com2,gx,gy,norm,norm2; 00660 /* the rest of the variables are used for pseudo-ordering 00661 the gradient magnitude values */ 00662 int list_count = 0; 00663 struct coorlist * list; 00664 struct coorlist ** range_l_s; /* array of pointers to start of bin list */ 00665 struct coorlist ** range_l_e; /* array of pointers to end of bin list */ 00666 struct coorlist * start; 00667 struct coorlist * end; 00668 00669 /* check parameters */ 00670 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 ) 00671 error("ll_angle: invalid image."); 00672 if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive."); 00673 if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'."); 00674 if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'."); 00675 if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'."); 00676 if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive."); 00677 if( max_grad <= 0.0 ) error("ll_angle: 'max_grad' must be positive."); 00678 00679 /* image size shortcuts */ 00680 n = in->ysize; 00681 p = in->xsize; 00682 00683 /* allocate output image */ 00684 g = new_image_double(in->xsize,in->ysize); 00685 00686 /* get memory for the image of gradient modulus */ 00687 *modgrad = new_image_double(in->xsize,in->ysize); 00688 00689 /* get memory for "ordered" list of pixels */ 00690 list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) ); 00691 *mem_p = (void *) list; 00692 range_l_s = (struct coorlist **) calloc( (size_t) n_bins, 00693 sizeof(struct coorlist *) ); 00694 range_l_e = (struct coorlist **) calloc( (size_t) n_bins, 00695 sizeof(struct coorlist *) ); 00696 if( list == NULL || range_l_s == NULL || range_l_e == NULL ) 00697 error("not enough memory."); 00698 for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL; 00699 00700 /* 'undefined' on the down and right boundaries */ 00701 for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF; 00702 for(y=0;y<n;y++) g->data[p*y+p-1] = NOTDEF; 00703 00704 /* compute gradient on the remaining pixels */ 00705 for(x=0;x<p-1;x++) 00706 for(y=0;y<n-1;y++) 00707 { 00708 adr = y*p+x; 00709 00710 /* 00711 Norm 2 computation using 2x2 pixel window: 00712 A B 00713 C D 00714 and 00715 com1 = D-A, com2 = B-C. 00716 Then 00717 gx = B+D - (A+C) horizontal difference 00718 gy = C+D - (A+B) vertical difference 00719 com1 and com2 are just to avoid 2 additions. 00720 */ 00721 com1 = in->data[adr+p+1] - in->data[adr]; 00722 com2 = in->data[adr+1] - in->data[adr+p]; 00723 00724 gx = com1+com2; /* gradient x component */ 00725 gy = com1-com2; /* gradient y component */ 00726 norm2 = gx*gx+gy*gy; 00727 norm = sqrt( norm2 / 4.0 ); /* gradient norm */ 00728 00729 (*modgrad)->data[adr] = norm; /* store gradient norm */ 00730 00731 if( norm <= threshold ) /* norm too small, gradient no defined */ 00732 g->data[adr] = NOTDEF; /* gradient angle not defined */ 00733 else 00734 { 00735 /* gradient angle computation */ 00736 g->data[adr] = atan2(gx,-gy); 00737 00738 /* store the point in the right bin according to its norm */ 00739 i = (unsigned int) (norm * (double) n_bins / max_grad); 00740 if( i >= n_bins ) i = n_bins-1; 00741 if( range_l_e[i] == NULL ) 00742 range_l_s[i] = range_l_e[i] = list+list_count++; 00743 else 00744 { 00745 range_l_e[i]->next = list+list_count; 00746 range_l_e[i] = list+list_count++; 00747 } 00748 range_l_e[i]->x = (int) x; 00749 range_l_e[i]->y = (int) y; 00750 range_l_e[i]->next = NULL; 00751 } 00752 } 00753 00754 /* Make the list of pixels (almost) ordered by norm value. 00755 It starts by the larger bin, so the list starts by the 00756 pixels with higher gradient value. Pixels would be ordered 00757 by norm value, up to a precision given by max_grad/n_bins. 00758 */ 00759 for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--); 00760 start = range_l_s[i]; 00761 end = range_l_e[i]; 00762 if( start != NULL ) 00763 for(i--;i>0; i--) 00764 if( range_l_s[i] != NULL ) 00765 { 00766 end->next = range_l_s[i]; 00767 end = range_l_e[i]; 00768 } 00769 *list_p = start; 00770 00771 /* free memory */ 00772 free( (void *) range_l_s ); 00773 free( (void *) range_l_e ); 00774 00775 return g; 00776 } |
Here is the call graph for this function:
|
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 868 of file lsd.c.
00869 { 00870 static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477, 00871 8687.24529705, 1168.92649479, 83.8676043424, 00872 2.50662827511 }; 00873 double a = (x+0.5) * log(x+5.5) - (x+5.5); 00874 double b = 0.0; 00875 int n; 00876 00877 for(n=0;n<7;n++) 00878 { 00879 a -= log( x + (double) n ); 00880 b += q[n] * pow( x, (double) n ); 00881 } 00882 return a + log(b); 00883 } |
|
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 902 of file lsd.c.
00903 {
00904 return 0.918938533204673 + (x-0.5)*log(x) - x
00905 + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
00906 }
|
|
LSD Simple Interface.
Definition at line 2077 of file lsd.c. References image_double, lsd_scale(), and ntuple_list.
02078 { 02079 /* LSD parameters */ 02080 double scale = 0.8; /* Scale the image by Gaussian filter to 'scale'. */ 02081 02082 return lsd_scale(image,scale); 02083 } |
Here is the call graph for this function:
|
LSD Simple Interface with Scale.
Definition at line 2053 of file lsd.c. References image_double, LineSegmentDetection(), and ntuple_list. Referenced by lsd().
02054 { 02055 /* LSD parameters */ 02056 double sigma_scale = 0.6; /* Sigma for Gaussian filter is computed as 02057 sigma = sigma_scale/scale. */ 02058 double quant = 2.0; /* Bound to the quantization error on the 02059 gradient norm. */ 02060 double ang_th = 22.5; /* Gradient angle tolerance in degrees. */ 02061 double eps = 0.0; /* Detection threshold, -log10(NFA). */ 02062 double density_th = 0.7; /* Minimal density of region points in rectangle. */ 02063 int n_bins = 1024; /* Number of bins in pseudo-ordering of gradient 02064 modulus. */ 02065 double max_grad = 255.0; /* Gradient modulus in the highest bin. The 02066 default value corresponds to the highest 02067 gradient modulus on images with gray 02068 levels in [0,255]. */ 02069 02070 return LineSegmentDetection( image, scale, sigma_scale, quant, ang_th, eps, 02071 density_th, n_bins, max_grad, NULL ); 02072 } |
Here is the call graph for this function:
|
Create a new image_char of size 'xsize' times 'ysize'.
Definition at line 291 of file lsd.c. References image_char_s::data, error(), image_char, image_char_s::xsize, and image_char_s::ysize. Referenced by new_image_char_ini().
00292 { 00293 image_char image; 00294 00295 /* check parameters */ 00296 if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size."); 00297 00298 /* get memory */ 00299 image = (image_char) malloc( sizeof(struct image_char_s) ); 00300 if( image == NULL ) error("not enough memory."); 00301 image->data = (unsigned char *) calloc( (size_t) (xsize*ysize), 00302 sizeof(unsigned char) ); 00303 if( image->data == NULL ) error("not enough memory."); 00304 00305 /* set image size */ 00306 image->xsize = xsize; 00307 image->ysize = ysize; 00308 00309 return image; 00310 } |
Here is the call graph for this function:
|
Create a new image_char of size 'xsize' times 'ysize', initialized to the value 'fill_value'.
Definition at line 316 of file lsd.c. References image_char_s::data, error(), image_char, and new_image_char(). Referenced by LineSegmentDetection().
00318 { 00319 image_char image = new_image_char(xsize,ysize); /* create image */ 00320 unsigned int N = xsize*ysize; 00321 unsigned int i; 00322 00323 /* check parameters */ 00324 if( image == NULL || image->data == NULL ) 00325 error("new_image_char_ini: invalid image."); 00326 00327 /* initialize */ 00328 for(i=0; i<N; i++) image->data[i] = fill_value; 00329 00330 return image; 00331 } |
Here is the call graph for this function:
|
Create a new image_double of size 'xsize' times 'ysize'.
Definition at line 398 of file lsd.c. References image_double_s::data, error(), image_double, image_double_s::xsize, and image_double_s::ysize. Referenced by gaussian_sampler(), ll_angle(), and new_image_double_ini().
00399 { 00400 image_double image; 00401 00402 /* check parameters */ 00403 if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size."); 00404 00405 /* get memory */ 00406 image = (image_double) malloc( sizeof(struct image_double_s) ); 00407 if( image == NULL ) error("not enough memory."); 00408 image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) ); 00409 if( image->data == NULL ) error("not enough memory."); 00410 00411 /* set image size */ 00412 image->xsize = xsize; 00413 image->ysize = ysize; 00414 00415 return image; 00416 } |
Here is the call graph for this function:
|
Create a new image_double of size 'xsize' times 'ysize', initialized to the value 'fill_value'.
Definition at line 422 of file lsd.c. References image_double_s::data, image_double, and new_image_double().
00424 { 00425 image_double image = new_image_double(xsize,ysize); /* create image */ 00426 unsigned int N = xsize*ysize; 00427 unsigned int i; 00428 00429 /* initialize */ 00430 for(i=0; i<N; i++) image->data[i] = fill_value; 00431 00432 return image; 00433 } |
Here is the call graph for this function:
|
Create a new image_int of size 'xsize' times 'ysize'.
Definition at line 347 of file lsd.c. References image_int_s::data, error(), image_int, image_int_s::xsize, and image_int_s::ysize. Referenced by new_image_int_ini().
00348 { 00349 image_int image; 00350 00351 /* check parameters */ 00352 if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size."); 00353 00354 /* get memory */ 00355 image = (image_int) malloc( sizeof(struct image_int_s) ); 00356 if( image == NULL ) error("not enough memory."); 00357 image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) ); 00358 if( image->data == NULL ) error("not enough memory."); 00359 00360 /* set image size */ 00361 image->xsize = xsize; 00362 image->ysize = ysize; 00363 00364 return image; 00365 } |
Here is the call graph for this function:
|
Create a new image_int of size 'xsize' times 'ysize', initialized to the value 'fill_value'.
Definition at line 371 of file lsd.c. References image_int_s::data, image_int, and new_image_int(). Referenced by LineSegmentDetection().
00373 { 00374 image_int image = new_image_int(xsize,ysize); /* create image */ 00375 unsigned int N = xsize*ysize; 00376 unsigned int i; 00377 00378 /* initialize */ 00379 for(i=0; i<N; i++) image->data[i] = fill_value; 00380 00381 return image; 00382 } |
Here is the call graph for this function:
|
Create an n-tuple list and allocate memory for one element.
Definition at line 206 of file lsd.c. References ntuple_list_s::dim, error(), ntuple_list_s::max_size, ntuple_list, ntuple_list_s::size, and ntuple_list_s::values. Referenced by gaussian_sampler(), and LineSegmentDetection().
00207 { 00208 ntuple_list n_tuple; 00209 00210 /* check parameters */ 00211 if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive."); 00212 00213 /* get memory for list structure */ 00214 n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) ); 00215 if( n_tuple == NULL ) error("not enough memory."); 00216 00217 /* initialize list */ 00218 n_tuple->size = 0; 00219 n_tuple->max_size = 1; 00220 n_tuple->dim = dim; 00221 00222 /* get memory for tuples */ 00223 n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) ); 00224 if( n_tuple->values == NULL ) error("not enough memory."); 00225 00226 return n_tuple; 00227 } |
Here is the call graph for this function:
|
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 962 of file lsd.c. References double_equal(), error(), log_gamma, M_LN10, and TABSIZE. Referenced by rect_nfa().
00963 { 00964 static double inv[TABSIZE]; /* table to keep computed inverse values */ 00965 double tolerance = 0.1; /* an error of 10% in the result is accepted */ 00966 double log1term,term,bin_term,mult_term,bin_tail,err,p_term; 00967 int i; 00968 00969 /* check parameters */ 00970 if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 ) 00971 error("nfa: wrong n, k or p values."); 00972 00973 /* trivial cases */ 00974 if( n==0 || k==0 ) return -logNT; 00975 if( n==k ) return -logNT - (double) n * log10(p); 00976 00977 /* probability term */ 00978 p_term = p / (1.0-p); 00979 00980 /* compute the first term of the series */ 00981 /* 00982 binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i} 00983 where bincoef(n,i) are the binomial coefficients. 00984 But 00985 bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ). 00986 We use this to compute the first term. Actually the log of it. 00987 */ 00988 log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 ) 00989 - log_gamma( (double) (n-k) + 1.0 ) 00990 + (double) k * log(p) + (double) (n-k) * log(1.0-p); 00991 term = exp(log1term); 00992 00993 /* in some cases no more computations are needed */ 00994 if( double_equal(term,0.0) ) /* the first term is almost zero */ 00995 { 00996 if( (double) k > (double) n * p ) /* at begin or end of the tail? */ 00997 return -log1term / M_LN10 - logNT; /* end: use just the first term */ 00998 else 00999 return -logNT; /* begin: the tail is roughly 1 */ 01000 } 01001 01002 /* compute more terms if needed */ 01003 bin_tail = term; 01004 for(i=k+1;i<=n;i++) 01005 { 01006 /* 01007 As 01008 term_i = bincoef(n,i) * p^i * (1-p)^(n-i) 01009 and 01010 bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i, 01011 then, 01012 term_i / term_i-1 = (n-i+1)/i * p/(1-p) 01013 and 01014 term_i = term_i-1 * (n-i+1)/i * p/(1-p). 01015 1/i is stored in a table as they are computed, 01016 because divisions are expensive. 01017 p/(1-p) is computed only once and stored in 'p_term'. 01018 */ 01019 bin_term = (double) (n-i+1) * ( i<TABSIZE ? 01020 ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) : 01021 1.0 / (double) i ); 01022 01023 mult_term = bin_term * p_term; 01024 term *= mult_term; 01025 bin_tail += term; 01026 if(bin_term<1.0) 01027 { 01028 /* When bin_term<1 then mult_term_j<mult_term_i for j>i. 01029 Then, the error on the binomial tail when truncated at 01030 the i term can be bounded by a geometric series of form 01031 term_i * sum mult_term_i^j. */ 01032 err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) / 01033 (1.0-mult_term) - 1.0 ); 01034 01035 /* One wants an error at most of tolerance*final_result, or: 01036 tolerance * abs(-log10(bin_tail)-logNT). 01037 Now, the error that can be accepted on bin_tail is 01038 given by tolerance*final_result divided by the derivative 01039 of -log10(x) when x=bin_tail. that is: 01040 tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail) 01041 Finally, we truncate the tail if the error is less than: 01042 tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */ 01043 if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break; 01044 } 01045 } 01046 return -log10(bin_tail) - logNT; 01047 } |
Here is the call graph for this function:
|
Copy one rectangle structure to another.
Definition at line 1071 of file lsd.c. References rect::dx, rect::dy, error(), rect::p, rect::prec, rect::theta, rect::width, rect::x, rect::x1, rect::x2, rect::y, rect::y1, and rect::y2. Referenced by rect_improve().
01072 { 01073 /* check parameters */ 01074 if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'."); 01075 01076 /* copy values */ 01077 out->x1 = in->x1; 01078 out->y1 = in->y1; 01079 out->x2 = in->x2; 01080 out->y2 = in->y2; 01081 out->width = in->width; 01082 out->x = in->x; 01083 out->y = in->y; 01084 out->theta = in->theta; 01085 out->dx = in->dx; 01086 out->dy = in->dy; 01087 out->prec = in->prec; 01088 out->p = in->p; 01089 } |
Here is the call graph for this function:
|
Try some rectangles variations to improve NFA value. Only if the rectangle is not meaningful (i.e., log_nfa <= eps). Definition at line 1644 of file lsd.c. References rect::dx, rect::dy, image_double, M_PI, rect::p, rect::prec, rect_copy(), rect_nfa(), rect::width, rect::x1, rect::x2, rect::y1, and rect::y2. Referenced by LineSegmentDetection().
01646 { 01647 struct rect r; 01648 double log_nfa,log_nfa_new; 01649 double delta = 0.5; 01650 double delta_2 = delta / 2.0; 01651 int n; 01652 01653 log_nfa = rect_nfa(rec,angles,logNT); 01654 01655 if( log_nfa > eps ) return log_nfa; 01656 01657 /* try finer precisions */ 01658 rect_copy(rec,&r); 01659 for(n=0; n<5; n++) 01660 { 01661 r.p /= 2.0; 01662 r.prec = r.p * M_PI; 01663 log_nfa_new = rect_nfa(&r,angles,logNT); 01664 if( log_nfa_new > log_nfa ) 01665 { 01666 log_nfa = log_nfa_new; 01667 rect_copy(&r,rec); 01668 } 01669 } 01670 01671 if( log_nfa > eps ) return log_nfa; 01672 01673 /* try to reduce width */ 01674 rect_copy(rec,&r); 01675 for(n=0; n<5; n++) 01676 { 01677 if( (r.width - delta) >= 0.5 ) 01678 { 01679 r.width -= delta; 01680 log_nfa_new = rect_nfa(&r,angles,logNT); 01681 if( log_nfa_new > log_nfa ) 01682 { 01683 rect_copy(&r,rec); 01684 log_nfa = log_nfa_new; 01685 } 01686 } 01687 } 01688 01689 if( log_nfa > eps ) return log_nfa; 01690 01691 /* try to reduce one side of the rectangle */ 01692 rect_copy(rec,&r); 01693 for(n=0; n<5; n++) 01694 { 01695 if( (r.width - delta) >= 0.5 ) 01696 { 01697 r.x1 += -r.dy * delta_2; 01698 r.y1 += r.dx * delta_2; 01699 r.x2 += -r.dy * delta_2; 01700 r.y2 += r.dx * delta_2; 01701 r.width -= delta; 01702 log_nfa_new = rect_nfa(&r,angles,logNT); 01703 if( log_nfa_new > log_nfa ) 01704 { 01705 rect_copy(&r,rec); 01706 log_nfa = log_nfa_new; 01707 } 01708 } 01709 } 01710 01711 if( log_nfa > eps ) return log_nfa; 01712 01713 /* try to reduce the other side of the rectangle */ 01714 rect_copy(rec,&r); 01715 for(n=0; n<5; n++) 01716 { 01717 if( (r.width - delta) >= 0.5 ) 01718 { 01719 r.x1 -= -r.dy * delta_2; 01720 r.y1 -= r.dx * delta_2; 01721 r.x2 -= -r.dy * delta_2; 01722 r.y2 -= r.dx * delta_2; 01723 r.width -= delta; 01724 log_nfa_new = rect_nfa(&r,angles,logNT); 01725 if( log_nfa_new > log_nfa ) 01726 { 01727 rect_copy(&r,rec); 01728 log_nfa = log_nfa_new; 01729 } 01730 } 01731 } 01732 01733 if( log_nfa > eps ) return log_nfa; 01734 01735 /* try even finer precisions */ 01736 rect_copy(rec,&r); 01737 for(n=0; n<5; n++) 01738 { 01739 r.p /= 2.0; 01740 r.prec = r.p * M_PI; 01741 log_nfa_new = rect_nfa(&r,angles,logNT); 01742 if( log_nfa_new > log_nfa ) 01743 { 01744 log_nfa = log_nfa_new; 01745 rect_copy(&r,rec); 01746 } 01747 } 01748 01749 return log_nfa; 01750 } |
Here is the call graph for this function:
|
Compute a rectangle's NFA value.
Definition at line 1370 of file lsd.c. References error(), image_double, isaligned(), nfa(), rect::p, rect::prec, ri_del(), ri_end(), ri_inc(), ri_ini(), rect::theta, rect_iter::x, image_double_s::xsize, rect_iter::y, and image_double_s::ysize. Referenced by rect_improve().
01371 { 01372 rect_iter * i; 01373 int pts = 0; 01374 int alg = 0; 01375 01376 /* check parameters */ 01377 if( rec == NULL ) error("rect_nfa: invalid rectangle."); 01378 if( angles == NULL ) error("rect_nfa: invalid 'angles'."); 01379 01380 /* compute the total number of pixels and of aligned points in 'rec' */ 01381 for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */ 01382 if( i->x >= 0 && i->y >= 0 && 01383 i->x < (int) angles->xsize && i->y < (int) angles->ysize ) 01384 { 01385 ++pts; /* total number of pixels counter */ 01386 if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) ) 01387 ++alg; /* aligned points counter */ 01388 } 01389 ri_del(i); /* delete iterator */ 01390 01391 return nfa(pts,alg,rec->p,logNT); /* compute NFA value */ 01392 } |
Here is the call graph for this function:
|
Reduce the region size, by elimination the points far from the starting point, until that leads to rectangle with the right density of region points or to discard the region if too small.
Definition at line 1757 of file lsd.c. References image_double_s::data, image_char_s::data, dist(), error(), FALSE, image_char, image_double, NOTUSED, region2rect(), TRUE, rect::width, point::x, rect::x1, rect::x2, image_char_s::xsize, point::y, rect::y1, and rect::y2. Referenced by refine().
01762 { 01763 double density,rad1,rad2,rad,xc,yc; 01764 int i; 01765 01766 /* check parameters */ 01767 if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'."); 01768 if( reg_size == NULL ) 01769 error("reduce_region_radius: invalid pointer 'reg_size'."); 01770 if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive."); 01771 if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'."); 01772 if( used == NULL || used->data == NULL ) 01773 error("reduce_region_radius: invalid image 'used'."); 01774 if( angles == NULL || angles->data == NULL ) 01775 error("reduce_region_radius: invalid image 'angles'."); 01776 01777 /* compute region points density */ 01778 density = (double) *reg_size / 01779 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width ); 01780 01781 /* if the density criterion is satisfied there is nothing to do */ 01782 if( density >= density_th ) return TRUE; 01783 01784 /* compute region's radius */ 01785 xc = (double) reg[0].x; 01786 yc = (double) reg[0].y; 01787 rad1 = dist( xc, yc, rec->x1, rec->y1 ); 01788 rad2 = dist( xc, yc, rec->x2, rec->y2 ); 01789 rad = rad1 > rad2 ? rad1 : rad2; 01790 01791 /* while the density criterion is not satisfied, remove farther pixels */ 01792 while( density < density_th ) 01793 { 01794 rad *= 0.75; /* reduce region's radius to 75% of its value */ 01795 01796 /* remove points from the region and update 'used' map */ 01797 for(i=0; i<*reg_size; i++) 01798 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad ) 01799 { 01800 /* point not kept, mark it as NOTUSED */ 01801 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED; 01802 /* remove point from the region */ 01803 reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */ 01804 reg[i].y = reg[*reg_size-1].y; 01805 --(*reg_size); 01806 --i; /* to avoid skipping one point */ 01807 } 01808 01809 /* reject if the region is too small. 01810 2 is the minimal region size for 'region2rect' to work. */ 01811 if( *reg_size < 2 ) return FALSE; 01812 01813 /* re-compute rectangle */ 01814 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec); 01815 01816 /* re-compute region points density */ 01817 density = (double) *reg_size / 01818 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width ); 01819 } 01820 01821 /* if this point is reached, the density criterion is satisfied */ 01822 return TRUE; 01823 } |
Here is the call graph for this function:
|
Refine a rectangle. For that, an estimation of the angle tolerance is performed by the standard deviation of the angle at points near the region's starting point. Then, a new region is grown starting from the same point, but using the estimated angle tolerance. If this fails to produce a rectangle with the right density of region points, 'reduce_region_radius' is called to try to satisfy this condition. Definition at line 1835 of file lsd.c. References angle_diff_signed(), image_double_s::data, image_char_s::data, dist(), error(), FALSE, image_char, image_double, NOTUSED, reduce_region_radius(), region2rect(), region_grow(), TRUE, rect::width, point::x, rect::x1, rect::x2, image_char_s::xsize, image_double_s::xsize, point::y, rect::y1, and rect::y2. Referenced by LineSegmentDetection().
01838 { 01839 double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum; 01840 int i,n; 01841 01842 /* check parameters */ 01843 if( reg == NULL ) error("refine: invalid pointer 'reg'."); 01844 if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'."); 01845 if( prec < 0.0 ) error("refine: 'prec' must be positive."); 01846 if( rec == NULL ) error("refine: invalid pointer 'rec'."); 01847 if( used == NULL || used->data == NULL ) 01848 error("refine: invalid image 'used'."); 01849 if( angles == NULL || angles->data == NULL ) 01850 error("refine: invalid image 'angles'."); 01851 01852 /* compute region points density */ 01853 density = (double) *reg_size / 01854 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width ); 01855 01856 /* if the density criterion is satisfied there is nothing to do */ 01857 if( density >= density_th ) return TRUE; 01858 01859 /*------ First try: reduce angle tolerance ------*/ 01860 01861 /* compute the new mean angle and tolerance */ 01862 xc = (double) reg[0].x; 01863 yc = (double) reg[0].y; 01864 ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ]; 01865 sum = s_sum = 0.0; 01866 n = 0; 01867 for(i=0; i<*reg_size; i++) 01868 { 01869 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED; 01870 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width ) 01871 { 01872 angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ]; 01873 ang_d = angle_diff_signed(angle,ang_c); 01874 sum += ang_d; 01875 s_sum += ang_d * ang_d; 01876 ++n; 01877 } 01878 } 01879 mean_angle = sum / (double) n; 01880 tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n 01881 + mean_angle*mean_angle ); /* 2 * standard deviation */ 01882 01883 /* find a new region from the same starting point and new angle tolerance */ 01884 region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,®_angle,used,tau); 01885 01886 /* if the region is too small, reject */ 01887 if( *reg_size < 2 ) return FALSE; 01888 01889 /* re-compute rectangle */ 01890 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec); 01891 01892 /* re-compute region points density */ 01893 density = (double) *reg_size / 01894 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width ); 01895 01896 /*------ Second try: reduce region radius ------*/ 01897 if( density < density_th ) 01898 return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p, 01899 rec, used, angles, density_th ); 01900 01901 /* if this point is reached, the density criterion is satisfied */ 01902 return TRUE; 01903 } |
Here is the call graph for this function:
|
Computes a rectangle that covers a region of points.
Definition at line 1499 of file lsd.c. References image_double_s::data, rect::dx, rect::dy, error(), get_theta(), image_double, rect::p, rect::prec, rect::theta, rect::width, rect::x, point::x, rect::x1, rect::x2, image_double_s::xsize, rect::y, point::y, rect::y1, and rect::y2. Referenced by LineSegmentDetection(), reduce_region_radius(), and refine().
01502 { 01503 double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max; 01504 int i; 01505 01506 /* check parameters */ 01507 if( reg == NULL ) error("region2rect: invalid region."); 01508 if( reg_size <= 1 ) error("region2rect: region size <= 1."); 01509 if( modgrad == NULL || modgrad->data == NULL ) 01510 error("region2rect: invalid image 'modgrad'."); 01511 if( rec == NULL ) error("region2rect: invalid 'rec'."); 01512 01513 /* center of the region: 01514 01515 It is computed as the weighted sum of the coordinates 01516 of all the pixels in the region. The norm of the gradient 01517 is used as the weight of a pixel. The sum is as follows: 01518 cx = \sum_i G(i).x_i 01519 cy = \sum_i G(i).y_i 01520 where G(i) is the norm of the gradient of pixel i 01521 and x_i,y_i are its coordinates. 01522 */ 01523 x = y = sum = 0.0; 01524 for(i=0; i<reg_size; i++) 01525 { 01526 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ]; 01527 x += (double) reg[i].x * weight; 01528 y += (double) reg[i].y * weight; 01529 sum += weight; 01530 } 01531 if( sum <= 0.0 ) error("region2rect: weights sum equal to zero."); 01532 x /= sum; 01533 y /= sum; 01534 01535 /* theta */ 01536 theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec); 01537 01538 /* length and width: 01539 01540 'l' and 'w' are computed as the distance from the center of the 01541 region to pixel i, projected along the rectangle axis (dx,dy) and 01542 to the orthogonal axis (-dy,dx), respectively. 01543 01544 The length of the rectangle goes from l_min to l_max, where l_min 01545 and l_max are the minimum and maximum values of l in the region. 01546 Analogously, the width is selected from w_min to w_max, where 01547 w_min and w_max are the minimum and maximum of w for the pixels 01548 in the region. 01549 */ 01550 dx = cos(theta); 01551 dy = sin(theta); 01552 l_min = l_max = w_min = w_max = 0.0; 01553 for(i=0; i<reg_size; i++) 01554 { 01555 l = ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy; 01556 w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx; 01557 01558 if( l > l_max ) l_max = l; 01559 if( l < l_min ) l_min = l; 01560 if( w > w_max ) w_max = w; 01561 if( w < w_min ) w_min = w; 01562 } 01563 01564 /* store values */ 01565 rec->x1 = x + l_min * dx; 01566 rec->y1 = y + l_min * dy; 01567 rec->x2 = x + l_max * dx; 01568 rec->y2 = y + l_max * dy; 01569 rec->width = w_max - w_min; 01570 rec->x = x; 01571 rec->y = y; 01572 rec->theta = theta; 01573 rec->dx = dx; 01574 rec->dy = dy; 01575 rec->prec = prec; 01576 rec->p = p; 01577 01578 /* we impose a minimal width of one pixel 01579 01580 A sharp horizontal or vertical step would produce a perfectly 01581 horizontal or vertical region. The width computed would be 01582 zero. But that corresponds to a one pixels width transition in 01583 the image. 01584 */ 01585 if( rec->width < 1.0 ) rec->width = 1.0; 01586 } |
Here is the call graph for this function:
|
Build a region of pixels that share the same angle, up to a tolerance 'prec', starting at point (x,y).
Definition at line 1592 of file lsd.c. References image_char_s::data, image_double_s::data, error(), image_char, image_double, isaligned(), USED, point::x, image_char_s::xsize, image_double_s::xsize, point::y, image_char_s::ysize, and image_double_s::ysize. Referenced by LineSegmentDetection(), and refine().
01595 { 01596 double sumdx,sumdy; 01597 int xx,yy,i; 01598 01599 /* check parameters */ 01600 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize ) 01601 error("region_grow: (x,y) out of the image."); 01602 if( angles == NULL || angles->data == NULL ) 01603 error("region_grow: invalid image 'angles'."); 01604 if( reg == NULL ) error("region_grow: invalid 'reg'."); 01605 if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'."); 01606 if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'."); 01607 if( used == NULL || used->data == NULL ) 01608 error("region_grow: invalid image 'used'."); 01609 01610 /* first point of the region */ 01611 *reg_size = 1; 01612 reg[0].x = x; 01613 reg[0].y = y; 01614 *reg_angle = angles->data[x+y*angles->xsize]; /* region's angle */ 01615 sumdx = cos(*reg_angle); 01616 sumdy = sin(*reg_angle); 01617 used->data[x+y*used->xsize] = USED; 01618 01619 /* try neighbors as new region points */ 01620 for(i=0; i<*reg_size; i++) 01621 for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++) 01622 for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++) 01623 if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize && 01624 used->data[xx+yy*used->xsize] != USED && 01625 isaligned(xx,yy,angles,*reg_angle,prec) ) 01626 { 01627 /* add point */ 01628 used->data[xx+yy*used->xsize] = USED; 01629 reg[*reg_size].x = xx; 01630 reg[*reg_size].y = yy; 01631 ++(*reg_size); 01632 01633 /* update region's angle */ 01634 sumdx += cos( angles->data[xx+yy*angles->xsize] ); 01635 sumdy += sin( angles->data[xx+yy*angles->xsize] ); 01636 *reg_angle = atan2(sumdy,sumdx); 01637 } 01638 } |
Here is the call graph for this function:
|
Free memory used by a rectangle iterator.
Definition at line 1202 of file lsd.c. References error(). Referenced by rect_nfa().
01203 { 01204 if( iter == NULL ) error("ri_del: NULL iterator."); 01205 free( (void *) iter ); 01206 } |
Here is the call graph for this function:
|
Check if the iterator finished the full iteration. See details in rect_iter Definition at line 1213 of file lsd.c. References error(), rect_iter::vx, and rect_iter::x. Referenced by rect_nfa(), and ri_inc().
01214 { 01215 /* check input */ 01216 if( i == NULL ) error("ri_end: NULL iterator."); 01217 01218 /* if the current x value is larger than the larger 01219 x value in the rectangle (vx[2]), we know the full 01220 exploration of the rectangle is finished. */ 01221 return (double)(i->x) > i->vx[2]; 01222 } |
Here is the call graph for this function:
|
Increment a rectangle iterator. See details in rect_iter Definition at line 1229 of file lsd.c. References error(), inter_hi(), inter_low(), ri_end(), rect_iter::vx, rect_iter::vy, rect_iter::x, rect_iter::y, rect_iter::ye, and rect_iter::ys. Referenced by rect_nfa(), and ri_ini().
01230 { 01231 /* check input */ 01232 if( i == NULL ) error("ri_inc: NULL iterator."); 01233 01234 /* if not at end of exploration, 01235 increase y value for next pixel in the 'column' */ 01236 if( !ri_end(i) ) i->y++; 01237 01238 /* if the end of the current 'column' is reached, 01239 and it is not the end of exploration, 01240 advance to the next 'column' */ 01241 while( (double) (i->y) > i->ye && !ri_end(i) ) 01242 { 01243 /* increase x, next 'column' */ 01244 i->x++; 01245 01246 /* if end of exploration, return */ 01247 if( ri_end(i) ) return; 01248 01249 /* update lower y limit (start) for the new 'column'. 01250 01251 We need to interpolate the y value that corresponds to the 01252 lower side of the rectangle. The first thing is to decide if 01253 the corresponding side is 01254 01255 vx[0],vy[0] to vx[3],vy[3] or 01256 vx[3],vy[3] to vx[2],vy[2] 01257 01258 Then, the side is interpolated for the x value of the 01259 'column'. But, if the side is vertical (as it could happen if 01260 the rectangle is vertical and we are dealing with the first 01261 or last 'columns') then we pick the lower value of the side 01262 by using 'inter_low'. 01263 */ 01264 if( (double) i->x < i->vx[3] ) 01265 i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]); 01266 else 01267 i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]); 01268 01269 /* update upper y limit (end) for the new 'column'. 01270 01271 We need to interpolate the y value that corresponds to the 01272 upper side of the rectangle. The first thing is to decide if 01273 the corresponding side is 01274 01275 vx[0],vy[0] to vx[1],vy[1] or 01276 vx[1],vy[1] to vx[2],vy[2] 01277 01278 Then, the side is interpolated for the x value of the 01279 'column'. But, if the side is vertical (as it could happen if 01280 the rectangle is vertical and we are dealing with the first 01281 or last 'columns') then we pick the lower value of the side 01282 by using 'inter_low'. 01283 */ 01284 if( (double)i->x < i->vx[1] ) 01285 i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]); 01286 else 01287 i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]); 01288 01289 /* new y */ 01290 i->y = (int) ceil(i->ys); 01291 } 01292 } |
Here is the call graph for this function:
|
Create and initialize a rectangle iterator. See details in rect_iter Definition at line 1299 of file lsd.c. References rect::dx, rect::dy, error(), ri_inc(), rect_iter::vx, rect_iter::vy, rect::width, rect_iter::x, rect::x1, rect::x2, rect_iter::y, rect::y1, rect::y2, rect_iter::ye, and rect_iter::ys. Referenced by rect_nfa().
01300 { 01301 double vx[4],vy[4]; 01302 int n,offset; 01303 rect_iter * i; 01304 01305 /* check parameters */ 01306 if( r == NULL ) error("ri_ini: invalid rectangle."); 01307 01308 /* get memory */ 01309 i = (rect_iter *) malloc(sizeof(rect_iter)); 01310 if( i == NULL ) error("ri_ini: Not enough memory."); 01311 01312 /* build list of rectangle corners ordered 01313 in a circular way around the rectangle */ 01314 vx[0] = r->x1 - r->dy * r->width / 2.0; 01315 vy[0] = r->y1 + r->dx * r->width / 2.0; 01316 vx[1] = r->x2 - r->dy * r->width / 2.0; 01317 vy[1] = r->y2 + r->dx * r->width / 2.0; 01318 vx[2] = r->x2 + r->dy * r->width / 2.0; 01319 vy[2] = r->y2 - r->dx * r->width / 2.0; 01320 vx[3] = r->x1 + r->dy * r->width / 2.0; 01321 vy[3] = r->y1 - r->dx * r->width / 2.0; 01322 01323 /* compute rotation of index of corners needed so that the first 01324 point has the smaller x. 01325 01326 if one side is vertical, thus two corners have the same smaller x 01327 value, the one with the largest y value is selected as the first. 01328 */ 01329 if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0; 01330 else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1; 01331 else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2; 01332 else offset = 3; 01333 01334 /* apply rotation of index. */ 01335 for(n=0; n<4; n++) 01336 { 01337 i->vx[n] = vx[(offset+n)%4]; 01338 i->vy[n] = vy[(offset+n)%4]; 01339 } 01340 01341 /* Set a initial condition. 01342 01343 The values are set to values that will cause 'ri_inc' (that will 01344 be called immediately) to initialize correctly the first 'column' 01345 and compute the limits 'ys' and 'ye'. 01346 01347 'y' is set to the integer value of vy[0], the starting corner. 01348 01349 'ys' and 'ye' are set to very small values, so 'ri_inc' will 01350 notice that it needs to start a new 'column'. 01351 01352 The smaller integer coordinate inside of the rectangle is 01353 'ceil(vx[0])'. The current 'x' value is set to that value minus 01354 one, so 'ri_inc' (that will increase x by one) will advance to 01355 the first 'column'. 01356 */ 01357 i->x = (int) ceil(i->vx[0]) - 1; 01358 i->y = (int) ceil(i->vy[0]); 01359 i->ys = i->ye = -DBL_MAX; 01360 01361 /* advance to the first pixel */ 01362 ri_inc(i); 01363 01364 return i; 01365 } |
Here is the call graph for this function: