直线检测算法LSD程序

/*---------------------------------------------------------------------------- LSD - Line Segment Detector on digital images Copyright 2007,2008,2009,2010 rafael grompone von gioi (grompone@gmail.com) This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details. You should have received a copy of the GNU Affero General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. ----------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------- This is an implementation of the Line Segment Detector described in the paper: "LSD: A Fast Line Segment Detector with a False Detection Control" by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 32, no. 4, pp. 722-732, April, 2010. and in more details in the CMLA Technical Report: "LSD: A Line Segment Detector, Technical Report", by Rafael Grompone von Gioi, Jeremie Jakubowicz, Jean-Michel Morel, Gregory Randall, CMLA, ENS Cachan, 2010. HISTORY: version 1.3 - feb 2010: Multiple bug correction and improved code. version 1.2 - dic 2009: First full Ansi C Language version. version 1.1 - sep 2009: Systematic subsampling to scale 0.8 and correction to partially handle "angle problem". version 1.0 - jan 2009: First complete Megawave2 and Ansi C Language version. ----------------------------------------------------------------------------*/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <float.h> #include "lsd.h" #ifndef M_LN10 #define M_LN10 2.30258509299404568402 #endif /* !M_LN10 */ #ifndef M_PI #define M_PI 3.14159265358979323846 #endif /* !M_PI */ #ifndef FALSE #define FALSE 0 #endif /* !FALSE */ #ifndef TRUE #define TRUE 1 #endif /* !TRUE */ #define NOTDEF -1024.0 #define M_3_2_PI 4.71238898038 #define M_2__PI 6.28318530718 #define NOTUSED 0 #define USED 1 /*----------------------------------------------------------------------------*/ struct coorlist { int x,y; struct coorlist * next; }; /*----------------------------------------------------------------------------*/ struct point {int x,y;}; /*----------------------------------------------------------------------------*/ /*------------------------- Miscellaneous functions --------------------------*/ /*----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /* Fatal error, print a message to standard-error output and exit. */ static void error(char * msg) { fprintf(stderr,"LSD Error: %s\n",msg); exit(EXIT_FAILURE); } /*----------------------------------------------------------------------------*/ /* 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. */ #define RELATIVE_ERROR_FACTOR 100.0 static int double_equal(double a, double b) { double abs_diff,aa,bb,abs_max; if( a == b ) return TRUE; abs_diff = fabs(a-b); aa = fabs(a); bb = fabs(b); abs_max = aa > bb ? aa : bb; /* DBL_MIN is the smallest normalized number, thus, the smallest number whose relative error is bounded by DBL_EPSILON. For smaller numbers, the same quantization steps as for DBL_MIN are used. Then, for smaller numbers, a meaningful "relative" error should be computed by dividing the difference by DBL_MIN. */ if( abs_max < DBL_MIN ) abs_max = DBL_MIN; /* equal if relative error <= factor x eps */ return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON); } /*----------------------------------------------------------------------------*/ /* Computes Euclidean distance between point (x1,y1) and point (x2,y2). */ static double dist(double x1, double y1, double x2, double y2) { return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ); } /*----------------------------------------------------------------------------*/ /*----------------------- 'list of n-tuple' data type ------------------------*/ /*----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /* Free memory used in n-tuple 'in'. */ void free_ntuple_list(ntuple_list in) { if( in == NULL || in->values == NULL ) error("free_ntuple_list: invalid n-tuple input."); free( (void *) in->values ); free( (void *) in ); } /*----------------------------------------------------------------------------*/ /* Create an n-tuple list and allocate memory for one element. The parameter 'dim' is the dimension (n) of the n-tuple. */ ntuple_list new_ntuple_list(unsigned int dim) { ntuple_list n_tuple; if( dim <= 0 ) error("new_ntuple_list: 'dim' must be positive."); n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) ); if( n_tuple == NULL ) error("not enough memory."); n_tuple->size = 0; n_tuple->max_size = 1; n_tuple->dim = dim; n_tuple->values = (double *) malloc( dim*n_tuple->max_size * sizeof(double) ); if( n_tuple->values == NULL ) error("not enough memory."); return n_tuple; } /*----------------------------------------------------------------------------*/ /* Enlarge the allocated memory of an n-tuple list. */ static void enlarge_ntuple_list(ntuple_list n_tuple) { if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size <= 0 ) error("enlarge_ntuple_list: invalid n-tuple."); n_tuple->max_size *= 2; n_tuple->values = (double *) realloc( (void *) n_tuple->values, n_tuple->dim * n_tuple->max_size * sizeof(double) ); if( n_tuple->values == NULL ) error("not enough memory."); } /*----------------------------------------------------------------------------*/ /* Add a 5-tuple to an n-tuple list. */ static void add_5tuple( ntuple_list out, double v1, double v2, double v3, double v4, double v5 ) { if( out == NULL ) error("add_5tuple: invalid n-tuple input."); if( out->dim != 5 ) error("add_5tuple: the n-tuple must be a 5-tuple."); if( out->size == out->max_size ) enlarge_ntuple_list(out); if( out->values == NULL ) error("add_5tuple: invalid n-tuple input."); out->values[ out->size * out->dim + 0 ] = v1; out->values[ out->size * out->dim + 1 ] = v2; out->values[ out->size * out->dim + 2 ] = v3; out->values[ out->size * out->dim + 3 ] = v4; out->values[ out->size * out->dim + 4 ] = v5; out->size++; } /*----------------------------------------------------------------------------*/ /*----------------------------- Image Data Types -----------------------------*/ /*----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /* Free memory used in image_char 'i'. */ void free_image_char(image_char i) { if( i == NULL || i->data == NULL ) error("free_image_char: invalid input image."); free( (void