00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 #include <stdio.h>
00072 #include <stdlib.h>
00073 #include <math.h>
00074 #include <limits.h>
00075 #include <float.h>
00076 #include "lsd.h"
00077
00078
00079 #ifndef M_LN10
00080 #define M_LN10 2.30258509299404568402
00081 #endif
00082
00083
00084 #ifndef M_PI
00085 #define M_PI 3.14159265358979323846
00086 #endif
00087
00088 #ifndef FALSE
00089 #define FALSE 0
00090 #endif
00091
00092 #ifndef TRUE
00093 #define TRUE 1
00094 #endif
00095
00096
00097 #define NOTDEF -1024.0
00098
00099
00100 #define M_3_2_PI 4.71238898038
00101
00102
00103 #define M_2__PI 6.28318530718
00104
00105
00106 #define NOTUSED 0
00107
00108
00109 #define USED 1
00110
00111
00112
00113
00114 struct coorlist
00115 {
00116 int x,y;
00117 struct coorlist * next;
00118 };
00119
00120
00121
00122
00123 struct point {int x,y;};
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 static void error(char * msg)
00134 {
00135 fprintf(stderr,"LSD Error: %s\n",msg);
00136 exit(EXIT_FAILURE);
00137 }
00138
00139
00140
00141
00142 #define RELATIVE_ERROR_FACTOR 100.0
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 static int double_equal(double a, double b)
00156 {
00157 double abs_diff,aa,bb,abs_max;
00158
00159
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
00168
00169
00170
00171
00172 if( abs_max < DBL_MIN ) abs_max = DBL_MIN;
00173
00174
00175 return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
00176 }
00177
00178
00179
00180
00181 static double dist(double x1, double y1, double x2, double y2)
00182 {
00183 return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00184 }
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 void free_ntuple_list(ntuple_list in)
00195 {
00196 if( in == NULL || in->values == NULL )
00197 error("free_ntuple_list: invalid n-tuple input.");
00198 free( (void *) in->values );
00199 free( (void *) in );
00200 }
00201
00202
00203
00204
00205
00206 ntuple_list new_ntuple_list(unsigned int dim)
00207 {
00208 ntuple_list n_tuple;
00209
00210
00211 if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");
00212
00213
00214 n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );
00215 if( n_tuple == NULL ) error("not enough memory.");
00216
00217
00218 n_tuple->size = 0;
00219 n_tuple->max_size = 1;
00220 n_tuple->dim = dim;
00221
00222
00223 n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) );
00224 if( n_tuple->values == NULL ) error("not enough memory.");
00225
00226 return n_tuple;
00227 }
00228
00229
00230
00231
00232 static void enlarge_ntuple_list(ntuple_list n_tuple)
00233 {
00234
00235 if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )
00236 error("enlarge_ntuple_list: invalid n-tuple.");
00237
00238
00239 n_tuple->max_size *= 2;
00240
00241
00242 n_tuple->values = (double *) realloc( (void *) n_tuple->values,
00243 n_tuple->dim * n_tuple->max_size * sizeof(double) );
00244 if( n_tuple->values == NULL ) error("not enough memory.");
00245 }
00246
00247
00248
00249
00250 static void add_5tuple( ntuple_list out, double v1, double v2,
00251 double v3, double v4, double v5 )
00252 {
00253
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
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
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
00269 out->size++;
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void free_image_char(image_char i)
00281 {
00282 if( i == NULL || i->data == NULL )
00283 error("free_image_char: invalid input image.");
00284 free( (void *) i->data );
00285 free( (void *) i );
00286 }
00287
00288
00289
00290
00291 image_char new_image_char(unsigned int xsize, unsigned int ysize)
00292 {
00293 image_char image;
00294
00295
00296 if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");
00297
00298
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
00306 image->xsize = xsize;
00307 image->ysize = ysize;
00308
00309 return image;
00310 }
00311
00312
00313
00314
00315
00316 image_char new_image_char_ini( unsigned int xsize, unsigned int ysize,
00317 unsigned char fill_value )
00318 {
00319 image_char image = new_image_char(xsize,ysize);
00320 unsigned int N = xsize*ysize;
00321 unsigned int i;
00322
00323
00324 if( image == NULL || image->data == NULL )
00325 error("new_image_char_ini: invalid image.");
00326
00327
00328 for(i=0; i<N; i++) image->data[i] = fill_value;
00329
00330 return image;
00331 }
00332
00333
00334
00335
00336 void free_image_int(image_int i)
00337 {
00338 if( i == NULL || i->data == NULL )
00339 error("free_image_int: invalid input image.");
00340 free( (void *) i->data );
00341 free( (void *) i );
00342 }
00343
00344
00345
00346
00347 image_int new_image_int(unsigned int xsize, unsigned int ysize)
00348 {
00349 image_int image;
00350
00351
00352 if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");
00353
00354
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
00361 image->xsize = xsize;
00362 image->ysize = ysize;
00363
00364 return image;
00365 }
00366
00367
00368
00369
00370
00371 image_int new_image_int_ini( unsigned int xsize, unsigned int ysize,
00372 int fill_value )
00373 {
00374 image_int image = new_image_int(xsize,ysize);
00375 unsigned int N = xsize*ysize;
00376 unsigned int i;
00377
00378
00379 for(i=0; i<N; i++) image->data[i] = fill_value;
00380
00381 return image;
00382 }
00383
00384
00385
00386
00387 void free_image_double(image_double i)
00388 {
00389 if( i == NULL || i->data == NULL )
00390 error("free_image_double: invalid input image.");
00391 free( (void *) i->data );
00392 free( (void *) i );
00393 }
00394
00395
00396
00397
00398 image_double new_image_double(unsigned int xsize, unsigned int ysize)
00399 {
00400 image_double image;
00401
00402
00403 if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");
00404
00405
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
00412 image->xsize = xsize;
00413 image->ysize = ysize;
00414
00415 return image;
00416 }
00417
00418
00419
00420
00421
00422 image_double new_image_double_ini( unsigned int xsize, unsigned int ysize,
00423 double fill_value )
00424 {
00425 image_double image = new_image_double(xsize,ysize);
00426 unsigned int N = xsize*ysize;
00427 unsigned int i;
00428
00429
00430 for(i=0; i<N; i++) image->data[i] = fill_value;
00431
00432 return image;
00433 }
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448 static void gaussian_kernel(ntuple_list kernel, double sigma, double mean)
00449 {
00450 double sum = 0.0;
00451 double val;
00452 unsigned int i;
00453
00454
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
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
00470 if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
00471 }
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 static image_double gaussian_sampler( image_double in, double scale,
00512 double sigma_scale )
00513 {
00514 image_double aux,out;
00515 ntuple_list kernel;
00516 unsigned int N,M,h,n,x,y,i;
00517 int xc,yc,j,double_x_size,double_y_size;
00518 double sigma,xx,yy,sum,prec;
00519
00520
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
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
00537 sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;
00538
00539
00540
00541
00542
00543
00544
00545
00546 prec = 3.0;
00547 h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
00548 n = 1+2*h;
00549 kernel = new_ntuple_list(n);
00550
00551
00552 double_x_size = (int) (2 * in->xsize);
00553 double_y_size = (int) (2 * in->ysize);
00554
00555
00556 for(x=0;x<aux->xsize;x++)
00557 {
00558
00559
00560
00561
00562
00563 xx = (double) x / scale;
00564
00565
00566 xc = (int) floor( xx + 0.5 );
00567 gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );
00568
00569
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
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
00590 for(y=0;y<out->ysize;y++)
00591 {
00592
00593
00594
00595
00596
00597 yy = (double) y / scale;
00598
00599
00600 yc = (int) floor( yy + 0.5 );
00601 gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );
00602
00603
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
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
00624 free_ntuple_list(kernel);
00625 free_image_double(aux);
00626
00627 return out;
00628 }
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652 static image_double ll_angle( image_double in, double threshold,
00653 struct coorlist ** list_p, void ** mem_p,
00654 image_double * modgrad, unsigned int n_bins,
00655 double max_grad )
00656 {
00657 image_double g;
00658 unsigned int n,p,x,y,adr,i;
00659 double com1,com2,gx,gy,norm,norm2;
00660
00661
00662 int list_count = 0;
00663 struct coorlist * list;
00664 struct coorlist ** range_l_s;
00665 struct coorlist ** range_l_e;
00666 struct coorlist * start;
00667 struct coorlist * end;
00668
00669
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
00680 n = in->ysize;
00681 p = in->xsize;
00682
00683
00684 g = new_image_double(in->xsize,in->ysize);
00685
00686
00687 *modgrad = new_image_double(in->xsize,in->ysize);
00688
00689
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
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
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
00712
00713
00714
00715
00716
00717
00718
00719
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;
00725 gy = com1-com2;
00726 norm2 = gx*gx+gy*gy;
00727 norm = sqrt( norm2 / 4.0 );
00728
00729 (*modgrad)->data[adr] = norm;
00730
00731 if( norm <= threshold )
00732 g->data[adr] = NOTDEF;
00733 else
00734 {
00735
00736 g->data[adr] = atan2(gx,-gy);
00737
00738
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
00755
00756
00757
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
00772 free( (void *) range_l_s );
00773 free( (void *) range_l_e );
00774
00775 return g;
00776 }
00777
00778
00779
00780
00781 static int isaligned( int x, int y, image_double angles, double theta,
00782 double prec )
00783 {
00784 double a;
00785
00786
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
00794 a = angles->data[ x + y * angles->xsize ];
00795
00796
00797
00798 if( a == NOTDEF ) return FALSE;
00799
00800
00801
00802
00803
00804
00805 theta -= a;
00806 if( theta < 0.0 ) theta = -theta;
00807 if( theta > M_3_2_PI )
00808 {
00809 theta -= M_2__PI;
00810 if( theta < 0.0 ) theta = -theta;
00811 }
00812
00813 return theta < prec;
00814 }
00815
00816
00817
00818
00819 static double angle_diff(double a, double b)
00820 {
00821 a -= b;
00822 while( a <= -M_PI ) a += M_2__PI;
00823 while( a > M_PI ) a -= M_2__PI;
00824 if( a < 0.0 ) a = -a;
00825 return a;
00826 }
00827
00828
00829
00830
00831 static double angle_diff_signed(double a, double b)
00832 {
00833 a -= b;
00834 while( a <= -M_PI ) a += M_2__PI;
00835 while( a > M_PI ) a -= M_2__PI;
00836 return a;
00837 }
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868 static double log_gamma_lanczos(double x)
00869 {
00870 static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
00871 8687.24529705, 1168.92649479, 83.8676043424,
00872 2.50662827511 };
00873 double a = (x+0.5) * log(x+5.5) - (x+5.5);
00874 double b = 0.0;
00875 int n;
00876
00877 for(n=0;n<7;n++)
00878 {
00879 a -= log( x + (double) n );
00880 b += q[n] * pow( x, (double) n );
00881 }
00882 return a + log(b);
00883 }
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902 static double log_gamma_windschitl(double x)
00903 {
00904 return 0.918938533204673 + (x-0.5)*log(x) - x
00905 + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
00906 }
00907
00908
00909
00910
00911
00912
00913 #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
00914
00915
00916
00917
00918 #define TABSIZE 100000
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962 static double nfa(int n, int k, double p, double logNT)
00963 {
00964 static double inv[TABSIZE];
00965 double tolerance = 0.1;
00966 double log1term,term,bin_term,mult_term,bin_tail,err,p_term;
00967 int i;
00968
00969
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
00974 if( n==0 || k==0 ) return -logNT;
00975 if( n==k ) return -logNT - (double) n * log10(p);
00976
00977
00978 p_term = p / (1.0-p);
00979
00980
00981
00982
00983
00984
00985
00986
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
00994 if( double_equal(term,0.0) )
00995 {
00996 if( (double) k > (double) n * p )
00997 return -log1term / M_LN10 - logNT;
00998 else
00999 return -logNT;
01000 }
01001
01002
01003 bin_tail = term;
01004 for(i=k+1;i<=n;i++)
01005 {
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
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
01029
01030
01031
01032 err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /
01033 (1.0-mult_term) - 1.0 );
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043 if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;
01044 }
01045 }
01046 return -log10(bin_tail) - logNT;
01047 }
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057 struct rect
01058 {
01059 double x1,y1,x2,y2;
01060 double width;
01061 double x,y;
01062 double theta;
01063 double dx,dy;
01064 double prec;
01065 double p;
01066 };
01067
01068
01069
01070
01071 static void rect_copy(struct rect * in, struct rect * out)
01072 {
01073
01074 if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");
01075
01076
01077 out->x1 = in->x1;
01078 out->y1 = in->y1;
01079 out->x2 = in->x2;
01080 out->y2 = in->y2;
01081 out->width = in->width;
01082 out->x = in->x;
01083 out->y = in->y;
01084 out->theta = in->theta;
01085 out->dx = in->dx;
01086 out->dy = in->dy;
01087 out->prec = in->prec;
01088 out->p = in->p;
01089 }
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147 typedef struct
01148 {
01149 double vx[4];
01150 double vy[4];
01151 double ys,ye;
01152 int x,y;
01153 } rect_iter;
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165 static double inter_low(double x, double x1, double y1, double x2, double y2)
01166 {
01167
01168 if( x1 > x2 || x < x1 || x > x2 )
01169 error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01170
01171
01172 if( double_equal(x1,x2) && y1<y2 ) return y1;
01173 if( double_equal(x1,x2) && y1>y2 ) return y2;
01174 return y1 + (x-x1) * (y2-y1) / (x2-x1);
01175 }
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187 static double inter_hi(double x, double x1, double y1, double x2, double y2)
01188 {
01189
01190 if( x1 > x2 || x < x1 || x > x2 )
01191 error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01192
01193
01194 if( double_equal(x1,x2) && y1<y2 ) return y2;
01195 if( double_equal(x1,x2) && y1>y2 ) return y1;
01196 return y1 + (x-x1) * (y2-y1) / (x2-x1);
01197 }
01198
01199
01200
01201
01202 static void ri_del(rect_iter * iter)
01203 {
01204 if( iter == NULL ) error("ri_del: NULL iterator.");
01205 free( (void *) iter );
01206 }
01207
01208
01209
01210
01211
01212
01213 static int ri_end(rect_iter * i)
01214 {
01215
01216 if( i == NULL ) error("ri_end: NULL iterator.");
01217
01218
01219
01220
01221 return (double)(i->x) > i->vx[2];
01222 }
01223
01224
01225
01226
01227
01228
01229 static void ri_inc(rect_iter * i)
01230 {
01231
01232 if( i == NULL ) error("ri_inc: NULL iterator.");
01233
01234
01235
01236 if( !ri_end(i) ) i->y++;
01237
01238
01239
01240
01241 while( (double) (i->y) > i->ye && !ri_end(i) )
01242 {
01243
01244 i->x++;
01245
01246
01247 if( ri_end(i) ) return;
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
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
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
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
01290 i->y = (int) ceil(i->ys);
01291 }
01292 }
01293
01294
01295
01296
01297
01298
01299 static rect_iter * ri_ini(struct rect * r)
01300 {
01301 double vx[4],vy[4];
01302 int n,offset;
01303 rect_iter * i;
01304
01305
01306 if( r == NULL ) error("ri_ini: invalid rectangle.");
01307
01308
01309 i = (rect_iter *) malloc(sizeof(rect_iter));
01310 if( i == NULL ) error("ri_ini: Not enough memory.");
01311
01312
01313
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
01324
01325
01326
01327
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
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
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
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
01362 ri_inc(i);
01363
01364 return i;
01365 }
01366
01367
01368
01369
01370 static double rect_nfa(struct rect * rec, image_double angles, double logNT)
01371 {
01372 rect_iter * i;
01373 int pts = 0;
01374 int alg = 0;
01375
01376
01377 if( rec == NULL ) error("rect_nfa: invalid rectangle.");
01378 if( angles == NULL ) error("rect_nfa: invalid 'angles'.");
01379
01380
01381 for(i=ri_ini(rec); !ri_end(i); ri_inc(i))
01382 if( i->x >= 0 && i->y >= 0 &&
01383 i->x < (int) angles->xsize && i->y < (int) angles->ysize )
01384 {
01385 ++pts;
01386 if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )
01387 ++alg;
01388 }
01389 ri_del(i);
01390
01391 return nfa(pts,alg,rec->p,logNT);
01392 }
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456 static double get_theta( struct point * reg, int reg_size, double x, double y,
01457 image_double modgrad, double reg_angle, double prec )
01458 {
01459 double lambda,theta,weight;
01460 double Ixx = 0.0;
01461 double Iyy = 0.0;
01462 double Ixy = 0.0;
01463 int i;
01464
01465
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
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
01484 lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );
01485
01486
01487 theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);
01488
01489
01490
01491 if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;
01492
01493 return theta;
01494 }
01495
01496
01497
01498
01499 static void region2rect( struct point * reg, int reg_size,
01500 image_double modgrad, double reg_angle,
01501 double prec, double p, struct rect * rec )
01502 {
01503 double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;
01504 int i;
01505
01506
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
01514
01515
01516
01517
01518
01519
01520
01521
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
01536 theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
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
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
01579
01580
01581
01582
01583
01584
01585 if( rec->width < 1.0 ) rec->width = 1.0;
01586 }
01587
01588
01589
01590
01591
01592 static void region_grow( int x, int y, image_double angles, struct point * reg,
01593 int * reg_size, double * reg_angle, image_char used,
01594 double prec )
01595 {
01596 double sumdx,sumdy;
01597 int xx,yy,i;
01598
01599
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
01611 *reg_size = 1;
01612 reg[0].x = x;
01613 reg[0].y = y;
01614 *reg_angle = angles->data[x+y*angles->xsize];
01615 sumdx = cos(*reg_angle);
01616 sumdy = sin(*reg_angle);
01617 used->data[x+y*used->xsize] = USED;
01618
01619
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
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
01634 sumdx += cos( angles->data[xx+yy*angles->xsize] );
01635 sumdy += sin( angles->data[xx+yy*angles->xsize] );
01636 *reg_angle = atan2(sumdy,sumdx);
01637 }
01638 }
01639
01640
01641
01642
01643
01644 static double rect_improve( struct rect * rec, image_double angles,
01645 double logNT, double eps )
01646 {
01647 struct rect r;
01648 double log_nfa,log_nfa_new;
01649 double delta = 0.5;
01650 double delta_2 = delta / 2.0;
01651 int n;
01652
01653 log_nfa = rect_nfa(rec,angles,logNT);
01654
01655 if( log_nfa > eps ) return log_nfa;
01656
01657
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
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
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
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
01736 rect_copy(rec,&r);
01737 for(n=0; n<5; n++)
01738 {
01739 r.p /= 2.0;
01740 r.prec = r.p * M_PI;
01741 log_nfa_new = rect_nfa(&r,angles,logNT);
01742 if( log_nfa_new > log_nfa )
01743 {
01744 log_nfa = log_nfa_new;
01745 rect_copy(&r,rec);
01746 }
01747 }
01748
01749 return log_nfa;
01750 }
01751
01752
01753
01754
01755
01756
01757 static int reduce_region_radius( struct point * reg, int * reg_size,
01758 image_double modgrad, double reg_angle,
01759 double prec, double p, struct rect * rec,
01760 image_char used, image_double angles,
01761 double density_th )
01762 {
01763 double density,rad1,rad2,rad,xc,yc;
01764 int i;
01765
01766
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
01778 density = (double) *reg_size /
01779 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01780
01781
01782 if( density >= density_th ) return TRUE;
01783
01784
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
01792 while( density < density_th )
01793 {
01794 rad *= 0.75;
01795
01796
01797 for(i=0; i<*reg_size; i++)
01798 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )
01799 {
01800
01801 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01802
01803 reg[i].x = reg[*reg_size-1].x;
01804 reg[i].y = reg[*reg_size-1].y;
01805 --(*reg_size);
01806 --i;
01807 }
01808
01809
01810
01811 if( *reg_size < 2 ) return FALSE;
01812
01813
01814 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01815
01816
01817 density = (double) *reg_size /
01818 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01819 }
01820
01821
01822 return TRUE;
01823 }
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835 static int refine( struct point * reg, int * reg_size, image_double modgrad,
01836 double reg_angle, double prec, double p, struct rect * rec,
01837 image_char used, image_double angles, double density_th )
01838 {
01839 double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;
01840 int i,n;
01841
01842
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
01853 density = (double) *reg_size /
01854 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01855
01856
01857 if( density >= density_th ) return TRUE;
01858
01859
01860
01861
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 );
01882
01883
01884 region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,®_angle,used,tau);
01885
01886
01887 if( *reg_size < 2 ) return FALSE;
01888
01889
01890 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01891
01892
01893 density = (double) *reg_size /
01894 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01895
01896
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
01902 return TRUE;
01903 }
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913 ntuple_list LineSegmentDetection( image_double image, double scale,
01914 double sigma_scale, double quant,
01915 double ang_th, double eps, double density_th,
01916 int n_bins, double max_grad,
01917 image_int * region )
01918 {
01919 ntuple_list out = new_ntuple_list(5);
01920 image_double scaled_image,angles,modgrad;
01921 image_char used;
01922 struct coorlist * list_p;
01923 void * mem_p;
01924 struct rect rec;
01925 struct point * reg;
01926 int reg_size,min_reg_size,i;
01927 unsigned int xsize,ysize;
01928 double rho,reg_angle,prec,p,log_nfa,logNT;
01929 int ls_count = 0;
01930
01931
01932
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
01947 prec = M_PI * ang_th / 180.0;
01948 p = ang_th / 180.0;
01949 rho = quant / sin(prec);
01950
01951
01952
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));
01967
01968
01969
01970
01971 if( region != NULL )
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
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
01983
01984 {
01985
01986 region_grow( list_p->x, list_p->y, angles, reg, ®_size,
01987 ®_angle, used, prec );
01988
01989
01990 if( reg_size < min_reg_size ) continue;
01991
01992
01993 region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004 if( !refine( reg, ®_size, modgrad, reg_angle,
02005 prec, p, &rec, used, angles, density_th ) ) continue;
02006
02007
02008 log_nfa = rect_improve(&rec,angles,logNT,eps);
02009 if( log_nfa <= eps ) continue;
02010
02011
02012 ++ls_count;
02013
02014
02015
02016
02017
02018
02019 rec.x1 += 0.5; rec.y1 += 0.5;
02020 rec.x2 += 0.5; rec.y2 += 0.5;
02021
02022
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
02031 add_5tuple(out, rec.x1, rec.y1, rec.x2, rec.y2, rec.width);
02032
02033
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
02041 free_image_double(angles);
02042 free_image_double(modgrad);
02043 free_image_char(used);
02044 free( (void *) reg );
02045 free( (void *) mem_p );
02046
02047 return out;
02048 }
02049
02050
02051
02052
02053 ntuple_list lsd_scale(image_double image, double scale)
02054 {
02055
02056 double sigma_scale = 0.6;
02057
02058 double quant = 2.0;
02059
02060 double ang_th = 22.5;
02061 double eps = 0.0;
02062 double density_th = 0.7;
02063 int n_bins = 1024;
02064
02065 double max_grad = 255.0;
02066
02067
02068
02069
02070 return LineSegmentDetection( image, scale, sigma_scale, quant, ang_th, eps,
02071 density_th, n_bins, max_grad, NULL );
02072 }
02073
02074
02075
02076
02077 ntuple_list lsd(image_double image)
02078 {
02079
02080 double scale = 0.8;
02081
02082 return lsd_scale(image,scale);
02083 }
02084