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
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 #include <stdio.h>
00098 #include <stdlib.h>
00099 #include <math.h>
00100 #include <limits.h>
00101 #include <float.h>
00102 #include "lsd.h"
00103
00104
00105 #ifndef M_LN10
00106 #define M_LN10 2.30258509299404568402
00107 #endif
00108
00109
00110 #ifndef M_PI
00111 #define M_PI 3.14159265358979323846
00112 #endif
00113
00114 #ifndef FALSE
00115 #define FALSE 0
00116 #endif
00117
00118 #ifndef TRUE
00119 #define TRUE 1
00120 #endif
00121
00122
00123 #define NOTDEF -1024.0
00124
00125
00126 #define M_3_2_PI 4.71238898038
00127
00128
00129 #define M_2__PI 6.28318530718
00130
00131
00132 #define NOTUSED 0
00133
00134
00135 #define USED 1
00136
00137
00138
00139
00140 struct coorlist
00141 {
00142 int x,y;
00143 struct coorlist * next;
00144 };
00145
00146
00147
00148
00149 struct point {int x,y;};
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 static void error(char * msg)
00160 {
00161 fprintf(stderr,"LSD Error: %s\n",msg);
00162 exit(EXIT_FAILURE);
00163 }
00164
00165
00166
00167
00168 #define RELATIVE_ERROR_FACTOR 100.0
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181 static int double_equal(double a, double b)
00182 {
00183 double abs_diff,aa,bb,abs_max;
00184
00185
00186 if( a == b ) return TRUE;
00187
00188 abs_diff = fabs(a-b);
00189 aa = fabs(a);
00190 bb = fabs(b);
00191 abs_max = aa > bb ? aa : bb;
00192
00193
00194
00195
00196
00197
00198 if( abs_max < DBL_MIN ) abs_max = DBL_MIN;
00199
00200
00201 return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
00202 }
00203
00204
00205
00206
00207 static double dist(double x1, double y1, double x2, double y2)
00208 {
00209 return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 typedef struct ntuple_list_s
00239 {
00240 unsigned int size;
00241 unsigned int max_size;
00242 unsigned int dim;
00243 double * values;
00244 } * ntuple_list;
00245
00246
00247
00248
00249 static void free_ntuple_list(ntuple_list in)
00250 {
00251 if( in == NULL || in->values == NULL )
00252 error("free_ntuple_list: invalid n-tuple input.");
00253 free( (void *) in->values );
00254 free( (void *) in );
00255 }
00256
00257
00258
00259
00260
00261 static ntuple_list new_ntuple_list(unsigned int dim)
00262 {
00263 ntuple_list n_tuple;
00264
00265
00266 if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");
00267
00268
00269 n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );
00270 if( n_tuple == NULL ) error("not enough memory.");
00271
00272
00273 n_tuple->size = 0;
00274 n_tuple->max_size = 1;
00275 n_tuple->dim = dim;
00276
00277
00278 n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) );
00279 if( n_tuple->values == NULL ) error("not enough memory.");
00280
00281 return n_tuple;
00282 }
00283
00284
00285
00286
00287 static void enlarge_ntuple_list(ntuple_list n_tuple)
00288 {
00289
00290 if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )
00291 error("enlarge_ntuple_list: invalid n-tuple.");
00292
00293
00294 n_tuple->max_size *= 2;
00295
00296
00297 n_tuple->values = (double *) realloc( (void *) n_tuple->values,
00298 n_tuple->dim * n_tuple->max_size * sizeof(double) );
00299 if( n_tuple->values == NULL ) error("not enough memory.");
00300 }
00301
00302
00303
00304
00305 static void add_7tuple( ntuple_list out, double v1, double v2, double v3,
00306 double v4, double v5, double v6, double v7 )
00307 {
00308
00309 if( out == NULL ) error("add_7tuple: invalid n-tuple input.");
00310 if( out->dim != 7 ) error("add_7tuple: the n-tuple must be a 7-tuple.");
00311
00312
00313 if( out->size == out->max_size ) enlarge_ntuple_list(out);
00314 if( out->values == NULL ) error("add_7tuple: invalid n-tuple input.");
00315
00316
00317 out->values[ out->size * out->dim + 0 ] = v1;
00318 out->values[ out->size * out->dim + 1 ] = v2;
00319 out->values[ out->size * out->dim + 2 ] = v3;
00320 out->values[ out->size * out->dim + 3 ] = v4;
00321 out->values[ out->size * out->dim + 4 ] = v5;
00322 out->values[ out->size * out->dim + 5 ] = v6;
00323 out->values[ out->size * out->dim + 6 ] = v7;
00324
00325
00326 out->size++;
00327 }
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 typedef struct image_char_s
00344 {
00345 unsigned char * data;
00346 unsigned int xsize,ysize;
00347 } * image_char;
00348
00349
00350
00351
00352 static void free_image_char(image_char i)
00353 {
00354 if( i == NULL || i->data == NULL )
00355 error("free_image_char: invalid input image.");
00356 free( (void *) i->data );
00357 free( (void *) i );
00358 }
00359
00360
00361
00362
00363 static image_char new_image_char(unsigned int xsize, unsigned int ysize)
00364 {
00365 image_char image;
00366
00367
00368 if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");
00369
00370
00371 image = (image_char) malloc( sizeof(struct image_char_s) );
00372 if( image == NULL ) error("not enough memory.");
00373 image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),
00374 sizeof(unsigned char) );
00375 if( image->data == NULL ) error("not enough memory.");
00376
00377
00378 image->xsize = xsize;
00379 image->ysize = ysize;
00380
00381 return image;
00382 }
00383
00384
00385
00386
00387
00388 static image_char new_image_char_ini( unsigned int xsize, unsigned int ysize,
00389 unsigned char fill_value )
00390 {
00391 image_char image = new_image_char(xsize,ysize);
00392 unsigned int N = xsize*ysize;
00393 unsigned int i;
00394
00395
00396 if( image == NULL || image->data == NULL )
00397 error("new_image_char_ini: invalid image.");
00398
00399
00400 for(i=0; i<N; i++) image->data[i] = fill_value;
00401
00402 return image;
00403 }
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414 typedef struct image_int_s
00415 {
00416 int * data;
00417 unsigned int xsize,ysize;
00418 } * image_int;
00419
00420
00421
00422
00423 static image_int new_image_int(unsigned int xsize, unsigned int ysize)
00424 {
00425 image_int image;
00426
00427
00428 if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");
00429
00430
00431 image = (image_int) malloc( sizeof(struct image_int_s) );
00432 if( image == NULL ) error("not enough memory.");
00433 image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );
00434 if( image->data == NULL ) error("not enough memory.");
00435
00436
00437 image->xsize = xsize;
00438 image->ysize = ysize;
00439
00440 return image;
00441 }
00442
00443
00444
00445
00446
00447 static image_int new_image_int_ini( unsigned int xsize, unsigned int ysize,
00448 int fill_value )
00449 {
00450 image_int image = new_image_int(xsize,ysize);
00451 unsigned int N = xsize*ysize;
00452 unsigned int i;
00453
00454
00455 for(i=0; i<N; i++) image->data[i] = fill_value;
00456
00457 return image;
00458 }
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469 typedef struct image_double_s
00470 {
00471 double * data;
00472 unsigned int xsize,ysize;
00473 } * image_double;
00474
00475
00476
00477
00478 static void free_image_double(image_double i)
00479 {
00480 if( i == NULL || i->data == NULL )
00481 error("free_image_double: invalid input image.");
00482 free( (void *) i->data );
00483 free( (void *) i );
00484 }
00485
00486
00487
00488
00489 static image_double new_image_double(unsigned int xsize, unsigned int ysize)
00490 {
00491 image_double image;
00492
00493
00494 if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");
00495
00496
00497 image = (image_double) malloc( sizeof(struct image_double_s) );
00498 if( image == NULL ) error("not enough memory.");
00499 image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );
00500 if( image->data == NULL ) error("not enough memory.");
00501
00502
00503 image->xsize = xsize;
00504 image->ysize = ysize;
00505
00506 return image;
00507 }
00508
00509
00510
00511
00512
00513 static image_double new_image_double_ptr( unsigned int xsize,
00514 unsigned int ysize, double * data )
00515 {
00516 image_double image;
00517
00518
00519 if( xsize == 0 || ysize == 0 )
00520 error("new_image_double_ptr: invalid image size.");
00521 if( data == NULL ) error("new_image_double_ptr: NULL data pointer.");
00522
00523
00524 image = (image_double) malloc( sizeof(struct image_double_s) );
00525 if( image == NULL ) error("not enough memory.");
00526
00527
00528 image->xsize = xsize;
00529 image->ysize = ysize;
00530 image->data = data;
00531
00532 return image;
00533 }
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548 static void gaussian_kernel(ntuple_list kernel, double sigma, double mean)
00549 {
00550 double sum = 0.0;
00551 double val;
00552 unsigned int i;
00553
00554
00555 if( kernel == NULL || kernel->values == NULL )
00556 error("gaussian_kernel: invalid n-tuple 'kernel'.");
00557 if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");
00558
00559
00560 if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);
00561 kernel->size = 1;
00562 for(i=0;i<kernel->dim;i++)
00563 {
00564 val = ( (double) i - mean ) / sigma;
00565 kernel->values[i] = exp( -0.5 * val * val );
00566 sum += kernel->values[i];
00567 }
00568
00569
00570 if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
00571 }
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 static image_double gaussian_sampler( image_double in, double scale,
00612 double sigma_scale )
00613 {
00614 image_double aux,out;
00615 ntuple_list kernel;
00616 unsigned int N,M,h,n,x,y,i;
00617 int xc,yc,j,double_x_size,double_y_size;
00618 double sigma,xx,yy,sum,prec;
00619
00620
00621 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00622 error("gaussian_sampler: invalid image.");
00623 if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");
00624 if( sigma_scale <= 0.0 )
00625 error("gaussian_sampler: 'sigma_scale' must be positive.");
00626
00627
00628 if( in->xsize * scale > (double) UINT_MAX ||
00629 in->ysize * scale > (double) UINT_MAX )
00630 error("gaussian_sampler: the output image size exceeds the handled size.");
00631 N = (unsigned int) ceil( in->xsize * scale );
00632 M = (unsigned int) ceil( in->ysize * scale );
00633 aux = new_image_double(N,in->ysize);
00634 out = new_image_double(N,M);
00635
00636
00637 sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;
00638
00639
00640
00641
00642
00643
00644
00645
00646 prec = 3.0;
00647 h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );
00648 n = 1+2*h;
00649 kernel = new_ntuple_list(n);
00650
00651
00652 double_x_size = (int) (2 * in->xsize);
00653 double_y_size = (int) (2 * in->ysize);
00654
00655
00656 for(x=0;x<aux->xsize;x++)
00657 {
00658
00659
00660
00661
00662
00663 xx = (double) x / scale;
00664
00665
00666 xc = (int) floor( xx + 0.5 );
00667 gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );
00668
00669
00670
00671 for(y=0;y<aux->ysize;y++)
00672 {
00673 sum = 0.0;
00674 for(i=0;i<kernel->dim;i++)
00675 {
00676 j = xc - h + i;
00677
00678
00679 while( j < 0 ) j += double_x_size;
00680 while( j >= double_x_size ) j -= double_x_size;
00681 if( j >= (int) in->xsize ) j = double_x_size-1-j;
00682
00683 sum += in->data[ j + y * in->xsize ] * kernel->values[i];
00684 }
00685 aux->data[ x + y * aux->xsize ] = sum;
00686 }
00687 }
00688
00689
00690 for(y=0;y<out->ysize;y++)
00691 {
00692
00693
00694
00695
00696
00697 yy = (double) y / scale;
00698
00699
00700 yc = (int) floor( yy + 0.5 );
00701 gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );
00702
00703
00704
00705 for(x=0;x<out->xsize;x++)
00706 {
00707 sum = 0.0;
00708 for(i=0;i<kernel->dim;i++)
00709 {
00710 j = yc - h + i;
00711
00712
00713 while( j < 0 ) j += double_y_size;
00714 while( j >= double_y_size ) j -= double_y_size;
00715 if( j >= (int) in->ysize ) j = double_y_size-1-j;
00716
00717 sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];
00718 }
00719 out->data[ x + y * out->xsize ] = sum;
00720 }
00721 }
00722
00723
00724 free_ntuple_list(kernel);
00725 free_image_double(aux);
00726
00727 return out;
00728 }
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752 static image_double ll_angle( image_double in, double threshold,
00753 struct coorlist ** list_p, void ** mem_p,
00754 image_double * modgrad, unsigned int n_bins )
00755 {
00756 image_double g;
00757 unsigned int n,p,x,y,adr,i;
00758 double com1,com2,gx,gy,norm,norm2;
00759
00760
00761 int list_count = 0;
00762 struct coorlist * list;
00763 struct coorlist ** range_l_s;
00764 struct coorlist ** range_l_e;
00765 struct coorlist * start;
00766 struct coorlist * end;
00767 double max_grad = 0.0;
00768
00769
00770 if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )
00771 error("ll_angle: invalid image.");
00772 if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");
00773 if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");
00774 if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");
00775 if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");
00776 if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");
00777
00778
00779 n = in->ysize;
00780 p = in->xsize;
00781
00782
00783 g = new_image_double(in->xsize,in->ysize);
00784
00785
00786 *modgrad = new_image_double(in->xsize,in->ysize);
00787
00788
00789 list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );
00790 *mem_p = (void *) list;
00791 range_l_s = (struct coorlist **) calloc( (size_t) n_bins,
00792 sizeof(struct coorlist *) );
00793 range_l_e = (struct coorlist **) calloc( (size_t) n_bins,
00794 sizeof(struct coorlist *) );
00795 if( list == NULL || range_l_s == NULL || range_l_e == NULL )
00796 error("not enough memory.");
00797 for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;
00798
00799
00800 for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;
00801 for(y=0;y<n;y++) g->data[p*y+p-1] = NOTDEF;
00802
00803
00804 for(x=0;x<p-1;x++)
00805 for(y=0;y<n-1;y++)
00806 {
00807 adr = y*p+x;
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820 com1 = in->data[adr+p+1] - in->data[adr];
00821 com2 = in->data[adr+1] - in->data[adr+p];
00822
00823 gx = com1+com2;
00824 gy = com1-com2;
00825 norm2 = gx*gx+gy*gy;
00826 norm = sqrt( norm2 / 4.0 );
00827
00828 (*modgrad)->data[adr] = norm;
00829
00830 if( norm <= threshold )
00831 g->data[adr] = NOTDEF;
00832 else
00833 {
00834
00835 g->data[adr] = atan2(gx,-gy);
00836
00837
00838 if( norm > max_grad ) max_grad = norm;
00839 }
00840 }
00841
00842
00843 for(x=0;x<p-1;x++)
00844 for(y=0;y<n-1;y++)
00845 {
00846 norm = (*modgrad)->data[y*p+x];
00847
00848
00849 i = (unsigned int) (norm * (double) n_bins / max_grad);
00850 if( i >= n_bins ) i = n_bins-1;
00851 if( range_l_e[i] == NULL )
00852 range_l_s[i] = range_l_e[i] = list+list_count++;
00853 else
00854 {
00855 range_l_e[i]->next = list+list_count;
00856 range_l_e[i] = list+list_count++;
00857 }
00858 range_l_e[i]->x = (int) x;
00859 range_l_e[i]->y = (int) y;
00860 range_l_e[i]->next = NULL;
00861 }
00862
00863
00864
00865
00866
00867
00868 for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);
00869 start = range_l_s[i];
00870 end = range_l_e[i];
00871 if( start != NULL )
00872 while(i>0)
00873 {
00874 --i;
00875 if( range_l_s[i] != NULL )
00876 {
00877 end->next = range_l_s[i];
00878 end = range_l_e[i];
00879 }
00880 }
00881 *list_p = start;
00882
00883
00884 free( (void *) range_l_s );
00885 free( (void *) range_l_e );
00886
00887 return g;
00888 }
00889
00890
00891
00892
00893 static int isaligned( int x, int y, image_double angles, double theta,
00894 double prec )
00895 {
00896 double a;
00897
00898
00899 if( angles == NULL || angles->data == NULL )
00900 error("isaligned: invalid image 'angles'.");
00901 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
00902 error("isaligned: (x,y) out of the image.");
00903 if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");
00904
00905
00906 a = angles->data[ x + y * angles->xsize ];
00907
00908
00909
00910 if( a == NOTDEF ) return FALSE;
00911
00912
00913
00914
00915
00916
00917 theta -= a;
00918 if( theta < 0.0 ) theta = -theta;
00919 if( theta > M_3_2_PI )
00920 {
00921 theta -= M_2__PI;
00922 if( theta < 0.0 ) theta = -theta;
00923 }
00924
00925 return theta <= prec;
00926 }
00927
00928
00929
00930
00931 static double angle_diff(double a, double b)
00932 {
00933 a -= b;
00934 while( a <= -M_PI ) a += M_2__PI;
00935 while( a > M_PI ) a -= M_2__PI;
00936 if( a < 0.0 ) a = -a;
00937 return a;
00938 }
00939
00940
00941
00942
00943 static double angle_diff_signed(double a, double b)
00944 {
00945 a -= b;
00946 while( a <= -M_PI ) a += M_2__PI;
00947 while( a > M_PI ) a -= M_2__PI;
00948 return a;
00949 }
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980 static double log_gamma_lanczos(double x)
00981 {
00982 static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
00983 8687.24529705, 1168.92649479, 83.8676043424,
00984 2.50662827511 };
00985 double a = (x+0.5) * log(x+5.5) - (x+5.5);
00986 double b = 0.0;
00987 int n;
00988
00989 for(n=0;n<7;n++)
00990 {
00991 a -= log( x + (double) n );
00992 b += q[n] * pow( x, (double) n );
00993 }
00994 return a + log(b);
00995 }
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014 static double log_gamma_windschitl(double x)
01015 {
01016 return 0.918938533204673 + (x-0.5)*log(x) - x
01017 + 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
01018 }
01019
01020
01021
01022
01023
01024
01025 #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
01026
01027
01028
01029
01030 #define TABSIZE 100000
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074 static double nfa(int n, int k, double p, double logNT)
01075 {
01076 static double inv[TABSIZE];
01077 double tolerance = 0.1;
01078 double log1term,term,bin_term,mult_term,bin_tail,err,p_term;
01079 int i;
01080
01081
01082 if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )
01083 error("nfa: wrong n, k or p values.");
01084
01085
01086 if( n==0 || k==0 ) return -logNT;
01087 if( n==k ) return -logNT - (double) n * log10(p);
01088
01089
01090 p_term = p / (1.0-p);
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100 log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )
01101 - log_gamma( (double) (n-k) + 1.0 )
01102 + (double) k * log(p) + (double) (n-k) * log(1.0-p);
01103 term = exp(log1term);
01104
01105
01106 if( double_equal(term,0.0) )
01107 {
01108 if( (double) k > (double) n * p )
01109 return -log1term / M_LN10 - logNT;
01110 else
01111 return -logNT;
01112 }
01113
01114
01115 bin_tail = term;
01116 for(i=k+1;i<=n;i++)
01117 {
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131 bin_term = (double) (n-i+1) * ( i<TABSIZE ?
01132 ( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :
01133 1.0 / (double) i );
01134
01135 mult_term = bin_term * p_term;
01136 term *= mult_term;
01137 bin_tail += term;
01138 if(bin_term<1.0)
01139 {
01140
01141
01142
01143
01144 err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /
01145 (1.0-mult_term) - 1.0 );
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155 if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;
01156 }
01157 }
01158 return -log10(bin_tail) - logNT;
01159 }
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169 struct rect
01170 {
01171 double x1,y1,x2,y2;
01172 double width;
01173 double x,y;
01174 double theta;
01175 double dx,dy;
01176 double prec;
01177 double p;
01178 };
01179
01180
01181
01182
01183 static void rect_copy(struct rect * in, struct rect * out)
01184 {
01185
01186 if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");
01187
01188
01189 out->x1 = in->x1;
01190 out->y1 = in->y1;
01191 out->x2 = in->x2;
01192 out->y2 = in->y2;
01193 out->width = in->width;
01194 out->x = in->x;
01195 out->y = in->y;
01196 out->theta = in->theta;
01197 out->dx = in->dx;
01198 out->dy = in->dy;
01199 out->prec = in->prec;
01200 out->p = in->p;
01201 }
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259 typedef struct
01260 {
01261 double vx[4];
01262 double vy[4];
01263 double ys,ye;
01264 int x,y;
01265 } rect_iter;
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277 static double inter_low(double x, double x1, double y1, double x2, double y2)
01278 {
01279
01280 if( x1 > x2 || x < x1 || x > x2 )
01281 error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01282
01283
01284 if( double_equal(x1,x2) && y1<y2 ) return y1;
01285 if( double_equal(x1,x2) && y1>y2 ) return y2;
01286 return y1 + (x-x1) * (y2-y1) / (x2-x1);
01287 }
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299 static double inter_hi(double x, double x1, double y1, double x2, double y2)
01300 {
01301
01302 if( x1 > x2 || x < x1 || x > x2 )
01303 error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");
01304
01305
01306 if( double_equal(x1,x2) && y1<y2 ) return y2;
01307 if( double_equal(x1,x2) && y1>y2 ) return y1;
01308 return y1 + (x-x1) * (y2-y1) / (x2-x1);
01309 }
01310
01311
01312
01313
01314 static void ri_del(rect_iter * iter)
01315 {
01316 if( iter == NULL ) error("ri_del: NULL iterator.");
01317 free( (void *) iter );
01318 }
01319
01320
01321
01322
01323
01324
01325 static int ri_end(rect_iter * i)
01326 {
01327
01328 if( i == NULL ) error("ri_end: NULL iterator.");
01329
01330
01331
01332
01333 return (double)(i->x) > i->vx[2];
01334 }
01335
01336
01337
01338
01339
01340
01341 static void ri_inc(rect_iter * i)
01342 {
01343
01344 if( i == NULL ) error("ri_inc: NULL iterator.");
01345
01346
01347
01348 if( !ri_end(i) ) i->y++;
01349
01350
01351
01352
01353 while( (double) (i->y) > i->ye && !ri_end(i) )
01354 {
01355
01356 i->x++;
01357
01358
01359 if( ri_end(i) ) return;
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376 if( (double) i->x < i->vx[3] )
01377 i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);
01378 else
01379 i->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396 if( (double)i->x < i->vx[1] )
01397 i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);
01398 else
01399 i->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);
01400
01401
01402 i->y = (int) ceil(i->ys);
01403 }
01404 }
01405
01406
01407
01408
01409
01410
01411 static rect_iter * ri_ini(struct rect * r)
01412 {
01413 double vx[4],vy[4];
01414 int n,offset;
01415 rect_iter * i;
01416
01417
01418 if( r == NULL ) error("ri_ini: invalid rectangle.");
01419
01420
01421 i = (rect_iter *) malloc(sizeof(rect_iter));
01422 if( i == NULL ) error("ri_ini: Not enough memory.");
01423
01424
01425
01426 vx[0] = r->x1 - r->dy * r->width / 2.0;
01427 vy[0] = r->y1 + r->dx * r->width / 2.0;
01428 vx[1] = r->x2 - r->dy * r->width / 2.0;
01429 vy[1] = r->y2 + r->dx * r->width / 2.0;
01430 vx[2] = r->x2 + r->dy * r->width / 2.0;
01431 vy[2] = r->y2 - r->dx * r->width / 2.0;
01432 vx[3] = r->x1 + r->dy * r->width / 2.0;
01433 vy[3] = r->y1 - r->dx * r->width / 2.0;
01434
01435
01436
01437
01438
01439
01440
01441 if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;
01442 else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;
01443 else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;
01444 else offset = 3;
01445
01446
01447 for(n=0; n<4; n++)
01448 {
01449 i->vx[n] = vx[(offset+n)%4];
01450 i->vy[n] = vy[(offset+n)%4];
01451 }
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469 i->x = (int) ceil(i->vx[0]) - 1;
01470 i->y = (int) ceil(i->vy[0]);
01471 i->ys = i->ye = -DBL_MAX;
01472
01473
01474 ri_inc(i);
01475
01476 return i;
01477 }
01478
01479
01480
01481
01482 static double rect_nfa(struct rect * rec, image_double angles, double logNT)
01483 {
01484 rect_iter * i;
01485 int pts = 0;
01486 int alg = 0;
01487
01488
01489 if( rec == NULL ) error("rect_nfa: invalid rectangle.");
01490 if( angles == NULL ) error("rect_nfa: invalid 'angles'.");
01491
01492
01493 for(i=ri_ini(rec); !ri_end(i); ri_inc(i))
01494 if( i->x >= 0 && i->y >= 0 &&
01495 i->x < (int) angles->xsize && i->y < (int) angles->ysize )
01496 {
01497 ++pts;
01498 if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )
01499 ++alg;
01500 }
01501 ri_del(i);
01502
01503 return nfa(pts,alg,rec->p,logNT);
01504 }
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568 static double get_theta( struct point * reg, int reg_size, double x, double y,
01569 image_double modgrad, double reg_angle, double prec )
01570 {
01571 double lambda,theta,weight;
01572 double Ixx = 0.0;
01573 double Iyy = 0.0;
01574 double Ixy = 0.0;
01575 int i;
01576
01577
01578 if( reg == NULL ) error("get_theta: invalid region.");
01579 if( reg_size <= 1 ) error("get_theta: region size <= 1.");
01580 if( modgrad == NULL || modgrad->data == NULL )
01581 error("get_theta: invalid 'modgrad'.");
01582 if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");
01583
01584
01585 for(i=0; i<reg_size; i++)
01586 {
01587 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01588 Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;
01589 Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;
01590 Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;
01591 }
01592 if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )
01593 error("get_theta: null inertia matrix.");
01594
01595
01596 lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );
01597
01598
01599 theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);
01600
01601
01602
01603 if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;
01604
01605 return theta;
01606 }
01607
01608
01609
01610
01611 static void region2rect( struct point * reg, int reg_size,
01612 image_double modgrad, double reg_angle,
01613 double prec, double p, struct rect * rec )
01614 {
01615 double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;
01616 int i;
01617
01618
01619 if( reg == NULL ) error("region2rect: invalid region.");
01620 if( reg_size <= 1 ) error("region2rect: region size <= 1.");
01621 if( modgrad == NULL || modgrad->data == NULL )
01622 error("region2rect: invalid image 'modgrad'.");
01623 if( rec == NULL ) error("region2rect: invalid 'rec'.");
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635 x = y = sum = 0.0;
01636 for(i=0; i<reg_size; i++)
01637 {
01638 weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];
01639 x += (double) reg[i].x * weight;
01640 y += (double) reg[i].y * weight;
01641 sum += weight;
01642 }
01643 if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");
01644 x /= sum;
01645 y /= sum;
01646
01647
01648 theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662 dx = cos(theta);
01663 dy = sin(theta);
01664 l_min = l_max = w_min = w_max = 0.0;
01665 for(i=0; i<reg_size; i++)
01666 {
01667 l = ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;
01668 w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;
01669
01670 if( l > l_max ) l_max = l;
01671 if( l < l_min ) l_min = l;
01672 if( w > w_max ) w_max = w;
01673 if( w < w_min ) w_min = w;
01674 }
01675
01676
01677 rec->x1 = x + l_min * dx;
01678 rec->y1 = y + l_min * dy;
01679 rec->x2 = x + l_max * dx;
01680 rec->y2 = y + l_max * dy;
01681 rec->width = w_max - w_min;
01682 rec->x = x;
01683 rec->y = y;
01684 rec->theta = theta;
01685 rec->dx = dx;
01686 rec->dy = dy;
01687 rec->prec = prec;
01688 rec->p = p;
01689
01690
01691
01692
01693
01694
01695
01696
01697 if( rec->width < 1.0 ) rec->width = 1.0;
01698 }
01699
01700
01701
01702
01703
01704 static void region_grow( int x, int y, image_double angles, struct point * reg,
01705 int * reg_size, double * reg_angle, image_char used,
01706 double prec )
01707 {
01708 double sumdx,sumdy;
01709 int xx,yy,i;
01710
01711
01712 if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )
01713 error("region_grow: (x,y) out of the image.");
01714 if( angles == NULL || angles->data == NULL )
01715 error("region_grow: invalid image 'angles'.");
01716 if( reg == NULL ) error("region_grow: invalid 'reg'.");
01717 if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");
01718 if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");
01719 if( used == NULL || used->data == NULL )
01720 error("region_grow: invalid image 'used'.");
01721
01722
01723 *reg_size = 1;
01724 reg[0].x = x;
01725 reg[0].y = y;
01726 *reg_angle = angles->data[x+y*angles->xsize];
01727 sumdx = cos(*reg_angle);
01728 sumdy = sin(*reg_angle);
01729 used->data[x+y*used->xsize] = USED;
01730
01731
01732 for(i=0; i<*reg_size; i++)
01733 for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)
01734 for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)
01735 if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&
01736 used->data[xx+yy*used->xsize] != USED &&
01737 isaligned(xx,yy,angles,*reg_angle,prec) )
01738 {
01739
01740 used->data[xx+yy*used->xsize] = USED;
01741 reg[*reg_size].x = xx;
01742 reg[*reg_size].y = yy;
01743 ++(*reg_size);
01744
01745
01746 sumdx += cos( angles->data[xx+yy*angles->xsize] );
01747 sumdy += sin( angles->data[xx+yy*angles->xsize] );
01748 *reg_angle = atan2(sumdy,sumdx);
01749 }
01750 }
01751
01752
01753
01754
01755
01756 static double rect_improve( struct rect * rec, image_double angles,
01757 double logNT, double log_eps )
01758 {
01759 struct rect r;
01760 double log_nfa,log_nfa_new;
01761 double delta = 0.5;
01762 double delta_2 = delta / 2.0;
01763 int n;
01764
01765 log_nfa = rect_nfa(rec,angles,logNT);
01766
01767 if( log_nfa > log_eps ) return log_nfa;
01768
01769
01770 rect_copy(rec,&r);
01771 for(n=0; n<5; n++)
01772 {
01773 r.p /= 2.0;
01774 r.prec = r.p * M_PI;
01775 log_nfa_new = rect_nfa(&r,angles,logNT);
01776 if( log_nfa_new > log_nfa )
01777 {
01778 log_nfa = log_nfa_new;
01779 rect_copy(&r,rec);
01780 }
01781 }
01782
01783 if( log_nfa > log_eps ) return log_nfa;
01784
01785
01786 rect_copy(rec,&r);
01787 for(n=0; n<5; n++)
01788 {
01789 if( (r.width - delta) >= 0.5 )
01790 {
01791 r.width -= delta;
01792 log_nfa_new = rect_nfa(&r,angles,logNT);
01793 if( log_nfa_new > log_nfa )
01794 {
01795 rect_copy(&r,rec);
01796 log_nfa = log_nfa_new;
01797 }
01798 }
01799 }
01800
01801 if( log_nfa > log_eps ) return log_nfa;
01802
01803
01804 rect_copy(rec,&r);
01805 for(n=0; n<5; n++)
01806 {
01807 if( (r.width - delta) >= 0.5 )
01808 {
01809 r.x1 += -r.dy * delta_2;
01810 r.y1 += r.dx * delta_2;
01811 r.x2 += -r.dy * delta_2;
01812 r.y2 += r.dx * delta_2;
01813 r.width -= delta;
01814 log_nfa_new = rect_nfa(&r,angles,logNT);
01815 if( log_nfa_new > log_nfa )
01816 {
01817 rect_copy(&r,rec);
01818 log_nfa = log_nfa_new;
01819 }
01820 }
01821 }
01822
01823 if( log_nfa > log_eps ) return log_nfa;
01824
01825
01826 rect_copy(rec,&r);
01827 for(n=0; n<5; n++)
01828 {
01829 if( (r.width - delta) >= 0.5 )
01830 {
01831 r.x1 -= -r.dy * delta_2;
01832 r.y1 -= r.dx * delta_2;
01833 r.x2 -= -r.dy * delta_2;
01834 r.y2 -= r.dx * delta_2;
01835 r.width -= delta;
01836 log_nfa_new = rect_nfa(&r,angles,logNT);
01837 if( log_nfa_new > log_nfa )
01838 {
01839 rect_copy(&r,rec);
01840 log_nfa = log_nfa_new;
01841 }
01842 }
01843 }
01844
01845 if( log_nfa > log_eps ) return log_nfa;
01846
01847
01848 rect_copy(rec,&r);
01849 for(n=0; n<5; n++)
01850 {
01851 r.p /= 2.0;
01852 r.prec = r.p * M_PI;
01853 log_nfa_new = rect_nfa(&r,angles,logNT);
01854 if( log_nfa_new > log_nfa )
01855 {
01856 log_nfa = log_nfa_new;
01857 rect_copy(&r,rec);
01858 }
01859 }
01860
01861 return log_nfa;
01862 }
01863
01864
01865
01866
01867
01868
01869 static int reduce_region_radius( struct point * reg, int * reg_size,
01870 image_double modgrad, double reg_angle,
01871 double prec, double p, struct rect * rec,
01872 image_char used, image_double angles,
01873 double density_th )
01874 {
01875 double density,rad1,rad2,rad,xc,yc;
01876 int i;
01877
01878
01879 if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");
01880 if( reg_size == NULL )
01881 error("reduce_region_radius: invalid pointer 'reg_size'.");
01882 if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");
01883 if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");
01884 if( used == NULL || used->data == NULL )
01885 error("reduce_region_radius: invalid image 'used'.");
01886 if( angles == NULL || angles->data == NULL )
01887 error("reduce_region_radius: invalid image 'angles'.");
01888
01889
01890 density = (double) *reg_size /
01891 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01892
01893
01894 if( density >= density_th ) return TRUE;
01895
01896
01897 xc = (double) reg[0].x;
01898 yc = (double) reg[0].y;
01899 rad1 = dist( xc, yc, rec->x1, rec->y1 );
01900 rad2 = dist( xc, yc, rec->x2, rec->y2 );
01901 rad = rad1 > rad2 ? rad1 : rad2;
01902
01903
01904 while( density < density_th )
01905 {
01906 rad *= 0.75;
01907
01908
01909 for(i=0; i<*reg_size; i++)
01910 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad )
01911 {
01912
01913 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01914
01915 reg[i].x = reg[*reg_size-1].x;
01916 reg[i].y = reg[*reg_size-1].y;
01917 --(*reg_size);
01918 --i;
01919 }
01920
01921
01922
01923 if( *reg_size < 2 ) return FALSE;
01924
01925
01926 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
01927
01928
01929 density = (double) *reg_size /
01930 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01931 }
01932
01933
01934 return TRUE;
01935 }
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947 static int refine( struct point * reg, int * reg_size, image_double modgrad,
01948 double reg_angle, double prec, double p, struct rect * rec,
01949 image_char used, image_double angles, double density_th )
01950 {
01951 double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;
01952 int i,n;
01953
01954
01955 if( reg == NULL ) error("refine: invalid pointer 'reg'.");
01956 if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");
01957 if( prec < 0.0 ) error("refine: 'prec' must be positive.");
01958 if( rec == NULL ) error("refine: invalid pointer 'rec'.");
01959 if( used == NULL || used->data == NULL )
01960 error("refine: invalid image 'used'.");
01961 if( angles == NULL || angles->data == NULL )
01962 error("refine: invalid image 'angles'.");
01963
01964
01965 density = (double) *reg_size /
01966 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
01967
01968
01969 if( density >= density_th ) return TRUE;
01970
01971
01972
01973
01974 xc = (double) reg[0].x;
01975 yc = (double) reg[0].y;
01976 ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];
01977 sum = s_sum = 0.0;
01978 n = 0;
01979 for(i=0; i<*reg_size; i++)
01980 {
01981 used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;
01982 if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width )
01983 {
01984 angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];
01985 ang_d = angle_diff_signed(angle,ang_c);
01986 sum += ang_d;
01987 s_sum += ang_d * ang_d;
01988 ++n;
01989 }
01990 }
01991 mean_angle = sum / (double) n;
01992 tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n
01993 + mean_angle*mean_angle );
01994
01995
01996 region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,®_angle,used,tau);
01997
01998
01999 if( *reg_size < 2 ) return FALSE;
02000
02001
02002 region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);
02003
02004
02005 density = (double) *reg_size /
02006 ( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );
02007
02008
02009 if( density < density_th )
02010 return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,
02011 rec, used, angles, density_th );
02012
02013
02014 return TRUE;
02015 }
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025 double * LineSegmentDetection( int * n_out,
02026 double * img, int X, int Y,
02027 double scale, double sigma_scale, double quant,
02028 double ang_th, double log_eps, double density_th,
02029 int n_bins,
02030 int ** reg_img, int * reg_x, int * reg_y )
02031 {
02032 image_double image;
02033 ntuple_list out = new_ntuple_list(7);
02034 double * return_value;
02035 image_double scaled_image,angles,modgrad;
02036 image_char used;
02037 image_int region = NULL;
02038 struct coorlist * list_p;
02039 void * mem_p;
02040 struct rect rec;
02041 struct point * reg;
02042 int reg_size,min_reg_size,i;
02043 unsigned int xsize,ysize;
02044 double rho,reg_angle,prec,p,log_nfa,logNT;
02045 int ls_count = 0;
02046
02047
02048
02049 if( img == NULL || X <= 0 || Y <= 0 ) error("invalid image input.");
02050 if( scale <= 0.0 ) error("'scale' value must be positive.");
02051 if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");
02052 if( quant < 0.0 ) error("'quant' value must be positive.");
02053 if( ang_th <= 0.0 || ang_th >= 180.0 )
02054 error("'ang_th' value must be in the range (0,180).");
02055 if( density_th < 0.0 || density_th > 1.0 )
02056 error("'density_th' value must be in the range [0,1].");
02057 if( n_bins <= 0 ) error("'n_bins' value must be positive.");
02058
02059
02060
02061 prec = M_PI * ang_th / 180.0;
02062 p = ang_th / 180.0;
02063 rho = quant / sin(prec);
02064
02065
02066
02067 image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );
02068 if( scale != 1.0 )
02069 {
02070 scaled_image = gaussian_sampler( image, scale, sigma_scale );
02071 angles = ll_angle( scaled_image, rho, &list_p, &mem_p,
02072 &modgrad, (unsigned int) n_bins );
02073 free_image_double(scaled_image);
02074 }
02075 else
02076 angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,
02077 (unsigned int) n_bins );
02078 xsize = angles->xsize;
02079 ysize = angles->ysize;
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094 logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0
02095 + log10(11.0);
02096 min_reg_size = (int) (-logNT/log10(p));
02097
02098
02099
02100
02101 if( reg_img != NULL && reg_x != NULL && reg_y != NULL )
02102 region = new_image_int_ini(angles->xsize,angles->ysize,0);
02103 used = new_image_char_ini(xsize,ysize,NOTUSED);
02104 reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );
02105 if( reg == NULL ) error("not enough memory!");
02106
02107
02108
02109 for(; list_p != NULL; list_p = list_p->next )
02110 if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&
02111 angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )
02112
02113
02114 {
02115
02116 region_grow( list_p->x, list_p->y, angles, reg, ®_size,
02117 ®_angle, used, prec );
02118
02119
02120 if( reg_size < min_reg_size ) continue;
02121
02122
02123 region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134 if( !refine( reg, ®_size, modgrad, reg_angle,
02135 prec, p, &rec, used, angles, density_th ) ) continue;
02136
02137
02138 log_nfa = rect_improve(&rec,angles,logNT,log_eps);
02139 if( log_nfa <= log_eps ) continue;
02140
02141
02142 ++ls_count;
02143
02144
02145
02146
02147
02148
02149 rec.x1 += 0.5; rec.y1 += 0.5;
02150 rec.x2 += 0.5; rec.y2 += 0.5;
02151
02152
02153 if( scale != 1.0 )
02154 {
02155 rec.x1 /= scale; rec.y1 /= scale;
02156 rec.x2 /= scale; rec.y2 /= scale;
02157 rec.width /= scale;
02158 }
02159
02160
02161 add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,
02162 rec.width, rec.p, log_nfa );
02163
02164
02165 if( region != NULL )
02166 for(i=0; i<reg_size; i++)
02167 region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;
02168 }
02169
02170
02171
02172 free( (void *) image );
02173
02174
02175 free_image_double(angles);
02176 free_image_double(modgrad);
02177 free_image_char(used);
02178 free( (void *) reg );
02179 free( (void *) mem_p );
02180
02181
02182 if( reg_img != NULL && reg_x != NULL && reg_y != NULL )
02183 {
02184 if( region == NULL ) error("'region' should be a valid image.");
02185 *reg_img = region->data;
02186 if( region->xsize > (unsigned int) INT_MAX ||
02187 region->xsize > (unsigned int) INT_MAX )
02188 error("region image to big to fit in INT sizes.");
02189 *reg_x = (int) (region->xsize);
02190 *reg_y = (int) (region->ysize);
02191
02192
02193
02194
02195 free( (void *) region );
02196 }
02197 if( out->size > (unsigned int) INT_MAX )
02198 error("too many detections to fit in an INT.");
02199 *n_out = (int) (out->size);
02200
02201 return_value = out->values;
02202 free( (void *) out );
02203
02204
02205
02206 return return_value;
02207 }
02208
02209
02210
02211
02212 double * lsd_scale_region( int * n_out,
02213 double * img, int X, int Y, double scale,
02214 int ** reg_img, int * reg_x, int * reg_y )
02215 {
02216
02217 double sigma_scale = 0.6;
02218
02219 double quant = 2.0;
02220
02221 double ang_th = 22.5;
02222 double log_eps = 0.0;
02223 double density_th = 0.7;
02224 int n_bins = 1024;
02225
02226
02227 return LineSegmentDetection( n_out, img, X, Y, scale, sigma_scale, quant,
02228 ang_th, log_eps, density_th, n_bins,
02229 reg_img, reg_x, reg_y );
02230 }
02231
02232
02233
02234
02235 double * lsd_scale(int * n_out, double * img, int X, int Y, double scale)
02236 {
02237 return lsd_scale_region(n_out,img,X,Y,scale,NULL,NULL,NULL);
02238 }
02239
02240
02241
02242
02243 double * lsd(int * n_out, double * img, int X, int Y)
02244 {
02245
02246 double scale = 0.8;
02247
02248 return lsd_scale(n_out,img,X,Y,scale);
02249 }
02250