环境：c++、opencv、VS2017

源码目录：

lsd.h
lsd.c
main.cpp

opencv配置：https://blog.csdn.net/q651742112/article/details/79769676

主文件main.cpp

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
extern "C" {
#include "lsd.h"
}#include <opencv2/opencv.hpp>int main(void)
{double * image;double * out;int x, y, i, j, n;int x0, y0, x1, y1;int X;  /* x image size */int Y;  /* y image size */cv::Point p0, p1;cv::Mat img = cv::imread("D:/ExperimentalPapers/MyPaperCode/understanding_of_indoor_scenes/datasets/V1.1/0000000004.jpg");//cv::imshow("org_img", img);cv::Mat gray;cv::cvtColor(img, gray, cv::COLOR_BGR2GRAY);X = img.cols;Y = img.rows;image = (double *)malloc(X * Y * sizeof(double));for (x = 0; x < X; x++)for (y = 0; y < Y; y++)image[x + y * X] = gray.at<uchar>(y, x);/* LSD call */out = lsd(&n, image,  X, Y);for (i = 0; i < n; i++){x0 = int(out[7 * i + 0]);y0 = int(out[7 * i + 1]);x1 = int(out[7 * i + 2]);y1 = int(out[7 * i + 3]);p0 = cv::Point(x0, y0);p1 = cv::Point(x1, y1);cv::line(img, p0, p1, cv::Scalar(0, 0, 255), 3, 4);}cv::imshow("lines_img", img);cv::waitKey(0);system("pause");
}

作者源代文件lsd.h

/*----------------------------------------------------------------------------LSD - Line Segment Detector on digital imagesThis code is part of the following publication and was subjectto peer review:"LSD: a Line Segment Detector" by Rafael Grompone von Gioi,Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsdhttp://dx.doi.org/10.5201/ipol.2012.gjmr-lsdCopyright (c) 2007-2011 rafael grompone von gioi <grompone@gmail.com>This program is free software: you can redistribute it and/or modifyit under the terms of the GNU Affero General Public License aspublished by the Free Software Foundation, either version 3 of theLicense, 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 ofMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See theGNU Affero General Public License for more details.You should have received a copy of the GNU Affero General Public Licensealong with this program. If not, see <http://www.gnu.org/licenses/>.----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** @file lsd.hLSD module header@author rafael grompone von gioi <grompone@gmail.com>*/
/*----------------------------------------------------------------------------*/
#ifndef LSD_HEADER
#define LSD_HEADER/*----------------------------------------------------------------------------*/
/** LSD Full Interface@param n_out       Pointer to an int where LSD will store the number ofline segments detected.@param img         Pointer to input image data. It must be an array ofdoubles of size X x Y, and the pixel at coordinates(x,y) is obtained by img[x+y*X].@param X           X size of the image: the number of columns.@param Y           Y size of the image: the number of rows.@param scale       When different from 1.0, LSD will scale the input imageby 'scale' factor by Gaussian filtering, before detectingline segments.Example: if scale=0.8, the input image will be subsampledto 80% of its size, before the line segment detectoris applied.Suggested value: 0.8@param sigma_scale When scale!=1.0, the sigma of the Gaussian filter is:sigma = sigma_scale / scale,   if scale <  1.0sigma = sigma_scale,           if scale >= 1.0Suggested value: 0.6@param quant       Bound to the quantization error on the gradient norm.Example: if gray levels are quantized to integer steps,the gradient (computed by finite differences) errordue to quantization will be bounded by 2.0, as theworst case is when the error are 1 and -1, thatgives an error of 2.0.Suggested value: 2.0@param ang_th      Gradient angle tolerance in the region growingalgorithm, in degrees.Suggested value: 22.5@param log_eps     Detection threshold, accept if -log10(NFA) > log_eps.The larger the value, the more strict the detector is,and will result in less detections.(Note that the 'minus sign' makes that thisbehavior is opposite to the one of NFA.)The value -log10(NFA) is equivalent but moreintuitive than NFA:- -1.0 gives an average of 10 false detections on noise-  0.0 gives an average of 1 false detections on noise-  1.0 gives an average of 0.1 false detections on nose-  2.0 gives an average of 0.01 false detections on noise.Suggested value: 0.0@param density_th  Minimal proportion of 'supporting' points in a rectangle.Suggested value: 0.7@param n_bins      Number of bins used in the pseudo-ordering of gradientmodulus.Suggested value: 1024@param reg_img     Optional output: if desired, LSD will return anint image where each pixel indicates the line segmentto which it belongs. Unused pixels have the value '0',while the used ones have the number of the line segment,numbered 1,2,3,..., in the same order as in theoutput list. If desired, a non NULL int** pointer mustbe assigned, and LSD will make that the pointer pointto an int array of size reg_x x reg_y, where the pixelvalue at (x,y) is obtained with (*reg_img)[x+y*reg_x].Note that the resulting image has the size of the imageused for the processing, that is, the size of the inputimage scaled by the given factor 'scale'. If scale!=1this size differs from XxY and that is the reason whyits value is given by reg_x and reg_y.Suggested value: NULL@param reg_x       Pointer to an int where LSD will put the X size'reg_img' image, when asked for.Suggested value: NULL@param reg_y       Pointer to an int where LSD will put the Y size'reg_img' image, when asked for.Suggested value: NULL@return            A double array of size 7 x n_out, containing the listof line segments detected. The array contains first7 values of line segment number 1, then the 7 valuesof line segment number 2, and so on, and it finishby the 7 values of line segment number n_out.The seven values are:- x1,y1,x2,y2,width,p,-log10(NFA).for a line segment from coordinates (x1,y1) to (x2,y2),a width 'width', an angle precision of p in (0,1) givenby angle_tolerance/180 degree, and NFA value 'NFA'.If 'out' is the returned pointer, the 7 values ofline segment number 'n+1' are obtained with'out[7*n+0]' to 'out[7*n+6]'.*/
double * LineSegmentDetection( int * n_out,double * img, int X, int Y,double scale, double sigma_scale, double quant,double ang_th, double log_eps, double density_th,int n_bins,int ** reg_img, int * reg_x, int * reg_y );/*----------------------------------------------------------------------------*/
/** LSD Simple Interface with Scale and Region output.@param n_out       Pointer to an int where LSD will store the number ofline segments detected.@param img         Pointer to input image data. It must be an array ofdoubles of size X x Y, and the pixel at coordinates(x,y) is obtained by img[x+y*X].@param X           X size of the image: the number of columns.@param Y           Y size of the image: the number of rows.@param scale       When different from 1.0, LSD will scale the input imageby 'scale' factor by Gaussian filtering, before detectingline segments.Example: if scale=0.8, the input image will be subsampledto 80% of its size, before the line segment detectoris applied.Suggested value: 0.8@param reg_img     Optional output: if desired, LSD will return anint image where each pixel indicates the line segmentto which it belongs. Unused pixels have the value '0',while the used ones have the number of the line segment,numbered 1,2,3,..., in the same order as in theoutput list. If desired, a non NULL int** pointer mustbe assigned, and LSD will make that the pointer pointto an int array of size reg_x x reg_y, where the pixelvalue at (x,y) is obtained with (*reg_img)[x+y*reg_x].Note that the resulting image has the size of the imageused for the processing, that is, the size of the inputimage scaled by the given factor 'scale'. If scale!=1this size differs from XxY and that is the reason whyits value is given by reg_x and reg_y.Suggested value: NULL@param reg_x       Pointer to an int where LSD will put the X size'reg_img' image, when asked for.Suggested value: NULL@param reg_y       Pointer to an int where LSD will put the Y size'reg_img' image, when asked for.Suggested value: NULL@return            A double array of size 7 x n_out, containing the listof line segments detected. The array contains first7 values of line segment number 1, then the 7 valuesof line segment number 2, and so on, and it finishby the 7 values of line segment number n_out.The seven values are:- x1,y1,x2,y2,width,p,-log10(NFA).for a line segment from coordinates (x1,y1) to (x2,y2),a width 'width', an angle precision of p in (0,1) givenby angle_tolerance/180 degree, and NFA value 'NFA'.If 'out' is the returned pointer, the 7 values ofline segment number 'n+1' are obtained with'out[7*n+0]' to 'out[7*n+6]'.*/
double * lsd_scale_region( int * n_out,double * img, int X, int Y, double scale,int ** reg_img, int * reg_x, int * reg_y );/*----------------------------------------------------------------------------*/
/** LSD Simple Interface with Scale@param n_out       Pointer to an int where LSD will store the number ofline segments detected.@param img         Pointer to input image data. It must be an array ofdoubles of size X x Y, and the pixel at coordinates(x,y) is obtained by img[x+y*X].@param X           X size of the image: the number of columns.@param Y           Y size of the image: the number of rows.@param scale       When different from 1.0, LSD will scale the input imageby 'scale' factor by Gaussian filtering, before detectingline segments.Example: if scale=0.8, the input image will be subsampledto 80% of its size, before the line segment detectoris applied.Suggested value: 0.8@return            A double array of size 7 x n_out, containing the listof line segments detected. The array contains first7 values of line segment number 1, then the 7 valuesof line segment number 2, and so on, and it finishby the 7 values of line segment number n_out.The seven values are:- x1,y1,x2,y2,width,p,-log10(NFA).for a line segment from coordinates (x1,y1) to (x2,y2),a width 'width', an angle precision of p in (0,1) givenby angle_tolerance/180 degree, and NFA value 'NFA'.If 'out' is the returned pointer, the 7 values ofline segment number 'n+1' are obtained with'out[7*n+0]' to 'out[7*n+6]'.*/
double * lsd_scale(int * n_out, double * img, int X, int Y, double scale);/*----------------------------------------------------------------------------*/
/** LSD Simple Interface@param n_out       Pointer to an int where LSD will store the number ofline segments detected.@param img         Pointer to input image data. It must be an array ofdoubles of size X x Y, and the pixel at coordinates(x,y) is obtained by img[x+y*X].@param X           X size of the image: the number of columns.@param Y           Y size of the image: the number of rows.@return            A double array of size 7 x n_out, containing the listof line segments detected. The array contains first7 values of line segment number 1, then the 7 valuesof line segment number 2, and so on, and it finishby the 7 values of line segment number n_out.The seven values are:- x1,y1,x2,y2,width,p,-log10(NFA).for a line segment from coordinates (x1,y1) to (x2,y2),a width 'width', an angle precision of p in (0,1) givenby angle_tolerance/180 degree, and NFA value 'NFA'.If 'out' is the returned pointer, the 7 values ofline segment number 'n+1' are obtained with'out[7*n+0]' to 'out[7*n+6]'.*/
double * lsd(int * n_out, double * img, int X, int Y);#endif /* !LSD_HEADER */
/*----------------------------------------------------------------------------*/

lsd.c

/*----------------------------------------------------------------------------LSD - Line Segment Detector on digital imagesThis code is part of the following publication and was subjectto peer review:"LSD: a Line Segment Detector" by Rafael Grompone von Gioi,Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsdhttp://dx.doi.org/10.5201/ipol.2012.gjmr-lsdCopyright (c) 2007-2011 rafael grompone von gioi <grompone@gmail.com>This program is free software: you can redistribute it and/or modifyit under the terms of the GNU Affero General Public License aspublished by the Free Software Foundation, either version 3 of theLicense, 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 ofMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See theGNU Affero General Public License for more details.You should have received a copy of the GNU Affero General Public Licensealong with this program. If not, see <http://www.gnu.org/licenses/>.----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** @file lsd.cLSD module code@author rafael grompone von gioi <grompone@gmail.com>*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** @mainpage LSD code documentationThis is an implementation of the Line Segment Detector describedin 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 andMachine 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.The version implemented here includes some further improvementsdescribed in the following publication, of which this code is part:"LSD: a Line Segment Detector" by Rafael Grompone von Gioi,Jeremie Jakubowicz, Jean-Michel Morel, and Gregory Randall,Image Processing On Line, 2012. DOI:10.5201/ipol.2012.gjmr-lsdhttp://dx.doi.org/10.5201/ipol.2012.gjmr-lsdThe module's main function is lsd().The source code is contained in two files: lsd.h and lsd.c.HISTORY:- version 1.6 - nov 2011:- changes in the interface,- max_grad parameter removed,- the factor 11 was added to the number of testto consider the different precision valuestested,- a minor bug corrected in the gradient sortingcode,- the algorithm now also returns p and log_nfafor each detection,- a minor bug was corrected in the image scaling,- the angle comparison in "isaligned" changedfrom < to <=,- "eps" variable renamed "log_eps",- "lsd_scale_region" interface was added,- minor changes to comments.- version 1.5 - dec 2010: Changes in 'refine', -W option added,and more comments added.- version 1.4 - jul 2010: lsd_scale interface added and doxygen doc.- version 1.3 - feb 2010: Multiple bug correction and improved code.- version 1.2 - dec 2009: First full Ansi C Language version.- version 1.1 - sep 2009: Systematic subsampling to scale 0.8 andcorrection to partially handle "angle problem".- version 1.0 - jan 2009: First complete Megawave2 and Ansi C Languageversion.@author rafael grompone von gioi <grompone@gmail.com>*/
/*----------------------------------------------------------------------------*/#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include "lsd.h"
/** ln(10) */
#ifndef M_LN10
#define M_LN10 2.30258509299404568402
#endif /* !M_LN10 *//** PI */
#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 *//** Label for pixels with undefined gradient. */
#define NOTDEF -1024.0/** 3/2 pi */
#define M_3_2_PI 4.71238898038/** 2 pi */
#define M_2__PI  6.28318530718/** Label for pixels not used in yet. */
#define NOTUSED 0/** Label for pixels already used in detection. */
#define USED    1/*----------------------------------------------------------------------------*/
/** Chained list of coordinates.*/
struct coorlist
{int x,y;struct coorlist * next;
};/*----------------------------------------------------------------------------*/
/** A point (or pixel).*/
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);
}/*----------------------------------------------------------------------------*/
/** Doubles relative error factor*/
#define RELATIVE_ERROR_FACTOR 100.0/*----------------------------------------------------------------------------*/
/** Compare doubles by relative error.The resulting rounding error after floating point computationsdepend on the specific operations done. The same number computed bydifferent algorithms could present different rounding errors. For auseful comparison, an estimation of the relative rounding errorshould be considered and compared to a factor times EPS. The factorshould be related to the cumulated rounding error in the chain ofcomputation. Here, as a simplification, a fixed factor is used.*/
static int double_equal(double a, double b)
{double abs_diff,aa,bb,abs_max;/* trivial case */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 smallestnumber whose relative error is bounded by DBL_EPSILON. Forsmaller numbers, the same quantization steps as for DBL_MINare 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 ------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** 'list of n-tuple' data typeThe i-th component of the j-th n-tuple of an n-tuple list 'ntl'is accessed with:ntl->values[ i + j * ntl->dim ]The dimension of the n-tuple (n) is:ntl->dimThe number of n-tuples in the list is:ntl->sizeThe maximum number of n-tuples that can be stored in thelist with the allocated memory at a given time is given by:ntl->max_size*/
typedef struct ntuple_list_s
{unsigned int size;unsigned int max_size;unsigned int dim;double * values;
} * ntuple_list;/*----------------------------------------------------------------------------*/
/** Free memory used in n-tuple 'in'.*/
static 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.@param dim the dimension (n) of the n-tuple.*/
static ntuple_list new_ntuple_list(unsigned int dim)
{ntuple_list n_tuple;/* check parameters */if( dim == 0 ) error("new_ntuple_list: 'dim' must be positive.");/* get memory for list structure */n_tuple = (ntuple_list) malloc( sizeof(struct ntuple_list_s) );if( n_tuple == NULL ) error("not enough memory.");/* initialize list */n_tuple->size = 0;n_tuple->max_size = 1;n_tuple->dim = dim;/* get memory for tuples */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)
{/* check parameters */if( n_tuple == NULL || n_tuple->values == NULL || n_tuple->max_size == 0 )error("enlarge_ntuple_list: invalid n-tuple.");/* duplicate number of tuples */n_tuple->max_size *= 2;/* realloc memory */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 7-tuple to an n-tuple list.*/
static void add_7tuple( ntuple_list out, double v1, double v2, double v3,double v4, double v5, double v6, double v7 )
{/* check parameters */if( out == NULL ) error("add_7tuple: invalid n-tuple input.");if( out->dim != 7 ) error("add_7tuple: the n-tuple must be a 7-tuple.");/* if needed, alloc more tuples to 'out' */if( out->size == out->max_size ) enlarge_ntuple_list(out);if( out->values == NULL ) error("add_7tuple: invalid n-tuple input.");/* add new 7-tuple */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->values[ out->size * out->dim + 5 ] = v6;out->values[ out->size * out->dim + 6 ] = v7;/* update number of tuples counter */out->size++;
}/*----------------------------------------------------------------------------*/
/*----------------------------- Image Data Types -----------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** char image data typeThe pixel value at (x,y) is accessed by:image->data[ x + y * image->xsize ]with x and y integer.*/
typedef struct image_char_s
{unsigned char * data;unsigned int xsize,ysize;
} * image_char;/*----------------------------------------------------------------------------*/
/** Free memory used in image_char 'i'.*/
static void free_image_char(image_char i)
{if( i == NULL || i->data == NULL )error("free_image_char: invalid input image.");free( (void *) i->data );free( (void *) i );
}/*----------------------------------------------------------------------------*/
/** Create a new image_char of size 'xsize' times 'ysize'.*/
static image_char new_image_char(unsigned int xsize, unsigned int ysize)
{image_char image;/* check parameters */if( xsize == 0 || ysize == 0 ) error("new_image_char: invalid image size.");/* get memory */image = (image_char) malloc( sizeof(struct image_char_s) );if( image == NULL ) error("not enough memory.");image->data = (unsigned char *) calloc( (size_t) (xsize*ysize),sizeof(unsigned char) );if( image->data == NULL ) error("not enough memory.");/* set image size */image->xsize = xsize;image->ysize = ysize;return image;
}/*----------------------------------------------------------------------------*/
/** Create a new image_char of size 'xsize' times 'ysize',initialized to the value 'fill_value'.*/
static image_char new_image_char_ini( unsigned int xsize, unsigned int ysize,unsigned char fill_value )
{image_char image = new_image_char(xsize,ysize); /* create image */unsigned int N = xsize*ysize;unsigned int i;/* check parameters */if( image == NULL || image->data == NULL )error("new_image_char_ini: invalid image.");/* initialize */for(i=0; i<N; i++) image->data[i] = fill_value;return image;
}/*----------------------------------------------------------------------------*/
/** int image data typeThe pixel value at (x,y) is accessed by:image->data[ x + y * image->xsize ]with x and y integer.*/
typedef struct image_int_s
{int * data;unsigned int xsize,ysize;
} * image_int;/*----------------------------------------------------------------------------*/
/** Create a new image_int of size 'xsize' times 'ysize'.*/
static image_int new_image_int(unsigned int xsize, unsigned int ysize)
{image_int image;/* check parameters */if( xsize == 0 || ysize == 0 ) error("new_image_int: invalid image size.");/* get memory */image = (image_int) malloc( sizeof(struct image_int_s) );if( image == NULL ) error("not enough memory.");image->data = (int *) calloc( (size_t) (xsize*ysize), sizeof(int) );if( image->data == NULL ) error("not enough memory.");/* set image size */image->xsize = xsize;image->ysize = ysize;return image;
}/*----------------------------------------------------------------------------*/
/** Create a new image_int of size 'xsize' times 'ysize',initialized to the value 'fill_value'.*/
static image_int new_image_int_ini( unsigned int xsize, unsigned int ysize,int fill_value )
{image_int image = new_image_int(xsize,ysize); /* create image */unsigned int N = xsize*ysize;unsigned int i;/* initialize */for(i=0; i<N; i++) image->data[i] = fill_value;return image;
}/*----------------------------------------------------------------------------*/
/** double image data typeThe pixel value at (x,y) is accessed by:image->data[ x + y * image->xsize ]with x and y integer.*/
typedef struct image_double_s
{double * data;unsigned int xsize,ysize;
} * image_double;/*----------------------------------------------------------------------------*/
/** Free memory used in image_double 'i'.*/
static void free_image_double(image_double i)
{if( i == NULL || i->data == NULL )error("free_image_double: invalid input image.");free( (void *) i->data );free( (void *) i );
}/*----------------------------------------------------------------------------*/
/** Create a new image_double of size 'xsize' times 'ysize'.*/
static image_double new_image_double(unsigned int xsize, unsigned int ysize)
{image_double image;/* check parameters */if( xsize == 0 || ysize == 0 ) error("new_image_double: invalid image size.");/* get memory */image = (image_double) malloc( sizeof(struct image_double_s) );if( image == NULL ) error("not enough memory.");image->data = (double *) calloc( (size_t) (xsize*ysize), sizeof(double) );if( image->data == NULL ) error("not enough memory.");/* set image size */image->xsize = xsize;image->ysize = ysize;return image;
}/*----------------------------------------------------------------------------*/
/** Create a new image_double of size 'xsize' times 'ysize'with the data pointed by 'data'.*/
static image_double new_image_double_ptr( unsigned int xsize,unsigned int ysize, double * data )
{image_double image;/* check parameters */if( xsize == 0 || ysize == 0 )error("new_image_double_ptr: invalid image size.");if( data == NULL ) error("new_image_double_ptr: NULL data pointer.");/* get memory */image = (image_double) malloc( sizeof(struct image_double_s) );if( image == NULL ) error("not enough memory.");/* set image */image->xsize = xsize;image->ysize = ysize;image->data = data;return image;
}/*----------------------------------------------------------------------------*/
/*----------------------------- Gaussian filter ------------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** Compute a Gaussian kernel of length 'kernel->dim',standard deviation 'sigma', and centered at value 'mean'.For example, if mean=0.5, the Gaussian will be centeredin the middle point between values 'kernel->values[0]'and 'kernel->values[1]'.*/
static void gaussian_kernel(ntuple_list kernel, double sigma, double mean)
{double sum = 0.0;double val;unsigned int i;/* check parameters */if( kernel == NULL || kernel->values == NULL )error("gaussian_kernel: invalid n-tuple 'kernel'.");if( sigma <= 0.0 ) error("gaussian_kernel: 'sigma' must be positive.");/* compute Gaussian kernel */if( kernel->max_size < 1 ) enlarge_ntuple_list(kernel);kernel->size = 1;for(i=0;i<kernel->dim;i++){val = ( (double) i - mean ) / sigma;kernel->values[i] = exp( -0.5 * val * val );sum += kernel->values[i];}/* normalization */if( sum >= 0.0 ) for(i=0;i<kernel->dim;i++) kernel->values[i] /= sum;
}/*----------------------------------------------------------------------------*/
/** Scale the input image 'in' by a factor 'scale' by Gaussian sub-sampling.For example, scale=0.8 will give a result at 80% of the original size.The image is convolved with a Gaussian kernel@f[G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}@f]before the sub-sampling to prevent aliasing.The standard deviation sigma given by:-  sigma = sigma_scale / scale,   if scale <  1.0-  sigma = sigma_scale,           if scale >= 1.0To be able to sub-sample at non-integer steps, some interpolationis needed. In this implementation, the interpolation is done bythe Gaussian kernel, so both operations (filtering and sampling)are done at the same time. The Gaussian kernel is computedcentered on the coordinates of the required sample. In this way,when applied, it gives directly the result of convolving the imagewith the kernel and interpolated to that particular position.A fast algorithm is done using the separability of the Gaussiankernel. Applying the 2D Gaussian kernel is equivalent to applyingfirst a horizontal 1D Gaussian kernel and then a vertical 1DGaussian kernel (or the other way round). The reason is that@f[G(x,y) = G(x) * G(y)@f]where@f[G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}}.@f]The algorithm first applies a combined Gaussian kernel and samplingin the x axis, and then the combined Gaussian kernel and samplingin the y axis.*/
static image_double gaussian_sampler( image_double in, double scale,double sigma_scale )
{image_double aux,out;ntuple_list kernel;unsigned int N,M,h,n,x,y,i;int xc,yc,j,double_x_size,double_y_size;double sigma,xx,yy,sum,prec;/* check parameters */if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )error("gaussian_sampler: invalid image.");if( scale <= 0.0 ) error("gaussian_sampler: 'scale' must be positive.");if( sigma_scale <= 0.0 )error("gaussian_sampler: 'sigma_scale' must be positive.");/* compute new image size and get memory for images */if( in->xsize * scale > (double) UINT_MAX ||in->ysize * scale > (double) UINT_MAX )error("gaussian_sampler: the output image size exceeds the handled size.");N = (unsigned int) ceil( in->xsize * scale );M = (unsigned int) ceil( in->ysize * scale );aux = new_image_double(N,in->ysize);out = new_image_double(N,M);/* sigma, kernel size and memory for the kernel */sigma = scale < 1.0 ? sigma_scale / scale : sigma_scale;/*The size of the kernel is selected to guarantee that thethe first discarded term is at least 10^prec times smallerthan the central value. For that, h should be larger than x, withe^(-x^2/2sigma^2) = 1/10^prec.Then,x = sigma * sqrt( 2 * prec * ln(10) ).*/prec = 3.0;h = (unsigned int) ceil( sigma * sqrt( 2.0 * prec * log(10.0) ) );n = 1+2*h; /* kernel size */kernel = new_ntuple_list(n);/* auxiliary double image size variables */double_x_size = (int) (2 * in->xsize);double_y_size = (int) (2 * in->ysize);/* First subsampling: x axis */for(x=0;x<aux->xsize;x++){/*x   is the coordinate in the new image.xx  is the corresponding x-value in the original size image.xc  is the integer value, the pixel coordinate of xx.*/xx = (double) x / scale;/* coordinate (0.0,0.0) is in the center of pixel (0,0),so the pixel with xc=0 get the values of xx from -0.5 to 0.5 */xc = (int) floor( xx + 0.5 );gaussian_kernel( kernel, sigma, (double) h + xx - (double) xc );/* the kernel must be computed for each x because the fineoffset xx-xc is different in each case */for(y=0;y<aux->ysize;y++){sum = 0.0;for(i=0;i<kernel->dim;i++){j = xc - h + i;/* symmetry boundary condition */while( j < 0 ) j += double_x_size;while( j >= double_x_size ) j -= double_x_size;if( j >= (int) in->xsize ) j = double_x_size-1-j;sum += in->data[ j + y * in->xsize ] * kernel->values[i];}aux->data[ x + y * aux->xsize ] = sum;}}/* Second subsampling: y axis */for(y=0;y<out->ysize;y++){/*y   is the coordinate in the new image.yy  is the corresponding x-value in the original size image.yc  is the integer value, the pixel coordinate of xx.*/yy = (double) y / scale;/* coordinate (0.0,0.0) is in the center of pixel (0,0),so the pixel with yc=0 get the values of yy from -0.5 to 0.5 */yc = (int) floor( yy + 0.5 );gaussian_kernel( kernel, sigma, (double) h + yy - (double) yc );/* the kernel must be computed for each y because the fineoffset yy-yc is different in each case */for(x=0;x<out->xsize;x++){sum = 0.0;for(i=0;i<kernel->dim;i++){j = yc - h + i;/* symmetry boundary condition */while( j < 0 ) j += double_y_size;while( j >= double_y_size ) j -= double_y_size;if( j >= (int) in->ysize ) j = double_y_size-1-j;sum += aux->data[ x + j * aux->xsize ] * kernel->values[i];}out->data[ x + y * out->xsize ] = sum;}}/* free memory */free_ntuple_list(kernel);free_image_double(aux);return out;
}/*----------------------------------------------------------------------------*/
/*--------------------------------- Gradient ---------------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** Computes the direction of the level line of 'in' at each point.The result is:- an image_double with the angle at each pixel, or NOTDEF if not defined.- the image_double 'modgrad' (a pointer is passed as argument)with the gradient magnitude at each point.- a list of pixels 'list_p' roughly ordered by decreasinggradient magnitude. (The order is made by classifying pointsinto bins by gradient magnitude. The parameters 'n_bins' and'max_grad' specify the number of bins and the gradient modulusat the highest bin. The pixels in the list would be indecreasing gradient magnitude, up to a precision of the size ofthe bins.)- a pointer 'mem_p' to the memory used by 'list_p' to be able tofree the memory when it is not used anymore.*/
static image_double ll_angle( image_double in, double threshold,struct coorlist ** list_p, void ** mem_p,image_double * modgrad, unsigned int n_bins )
{image_double g;unsigned int n,p,x,y,adr,i;double com1,com2,gx,gy,norm,norm2;/* the rest of the variables are used for pseudo-orderingthe gradient magnitude values */int list_count = 0;struct coorlist * list;struct coorlist ** range_l_s; /* array of pointers to start of bin list */struct coorlist ** range_l_e; /* array of pointers to end of bin list */struct coorlist * start;struct coorlist * end;double max_grad = 0.0;/* check parameters */if( in == NULL || in->data == NULL || in->xsize == 0 || in->ysize == 0 )error("ll_angle: invalid image.");if( threshold < 0.0 ) error("ll_angle: 'threshold' must be positive.");if( list_p == NULL ) error("ll_angle: NULL pointer 'list_p'.");if( mem_p == NULL ) error("ll_angle: NULL pointer 'mem_p'.");if( modgrad == NULL ) error("ll_angle: NULL pointer 'modgrad'.");if( n_bins == 0 ) error("ll_angle: 'n_bins' must be positive.");/* image size shortcuts */n = in->ysize;p = in->xsize;/* allocate output image */g = new_image_double(in->xsize,in->ysize);/* get memory for the image of gradient modulus */*modgrad = new_image_double(in->xsize,in->ysize);/* get memory for "ordered" list of pixels */list = (struct coorlist *) calloc( (size_t) (n*p), sizeof(struct coorlist) );*mem_p = (void *) list;range_l_s = (struct coorlist **) calloc( (size_t) n_bins,sizeof(struct coorlist *) );range_l_e = (struct coorlist **) calloc( (size_t) n_bins,sizeof(struct coorlist *) );if( list == NULL || range_l_s == NULL || range_l_e == NULL )error("not enough memory.");for(i=0;i<n_bins;i++) range_l_s[i] = range_l_e[i] = NULL;/* 'undefined' on the down and right boundaries */for(x=0;x<p;x++) g->data[(n-1)*p+x] = NOTDEF;for(y=0;y<n;y++) g->data[p*y+p-1]   = NOTDEF;/* compute gradient on the remaining pixels */for(x=0;x<p-1;x++)for(y=0;y<n-1;y++){adr = y*p+x;/*Norm 2 computation using 2x2 pixel window:A BC Dandcom1 = D-A,  com2 = B-C.Thengx = B+D - (A+C)   horizontal differencegy = C+D - (A+B)   vertical differencecom1 and com2 are just to avoid 2 additions.*/com1 = in->data[adr+p+1] - in->data[adr];com2 = in->data[adr+1]   - in->data[adr+p];gx = com1+com2; /* gradient x component */gy = com1-com2; /* gradient y component */norm2 = gx*gx+gy*gy;norm = sqrt( norm2 / 4.0 ); /* gradient norm */(*modgrad)->data[adr] = norm; /* store gradient norm */if( norm <= threshold ) /* norm too small, gradient no defined */g->data[adr] = NOTDEF; /* gradient angle not defined */else{/* gradient angle computation */g->data[adr] = atan2(gx,-gy);/* look for the maximum of the gradient */if( norm > max_grad ) max_grad = norm;}}/* compute histogram of gradient values */for(x=0;x<p-1;x++)for(y=0;y<n-1;y++){norm = (*modgrad)->data[y*p+x];/* store the point in the right bin according to its norm */i = (unsigned int) (norm * (double) n_bins / max_grad);if( i >= n_bins ) i = n_bins-1;if( range_l_e[i] == NULL )range_l_s[i] = range_l_e[i] = list+list_count++;else{range_l_e[i]->next = list+list_count;range_l_e[i] = list+list_count++;}range_l_e[i]->x = (int) x;range_l_e[i]->y = (int) y;range_l_e[i]->next = NULL;}/* Make the list of pixels (almost) ordered by norm value.It starts by the larger bin, so the list starts by thepixels with the highest gradient value. Pixels would be orderedby norm value, up to a precision given by max_grad/n_bins.*/for(i=n_bins-1; i>0 && range_l_s[i]==NULL; i--);start = range_l_s[i];end = range_l_e[i];if( start != NULL )while(i>0){--i;if( range_l_s[i] != NULL ){end->next = range_l_s[i];end = range_l_e[i];}}*list_p = start;/* free memory */free( (void *) range_l_s );free( (void *) range_l_e );return g;
}/*----------------------------------------------------------------------------*/
/** Is point (x,y) aligned to angle theta, up to precision 'prec'?*/
static int isaligned( int x, int y, image_double angles, double theta,double prec )
{double a;/* check parameters */if( angles == NULL || angles->data == NULL )error("isaligned: invalid image 'angles'.");if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )error("isaligned: (x,y) out of the image.");if( prec < 0.0 ) error("isaligned: 'prec' must be positive.");/* angle at pixel (x,y) */a = angles->data[ x + y * angles->xsize ];/* pixels whose level-line angle is not definedare considered as NON-aligned */if( a == NOTDEF ) return FALSE;  /* there is no need to call the function'double_equal' here because there isno risk of problems related to thecomparison doubles, we are onlyinterested in the exact NOTDEF value *//* it is assumed that 'theta' and 'a' are in the range [-pi,pi] */theta -= a;if( theta < 0.0 ) theta = -theta;if( theta > M_3_2_PI ){theta -= M_2__PI;if( theta < 0.0 ) theta = -theta;}return theta <= prec;
}/*----------------------------------------------------------------------------*/
/** Absolute value angle difference.*/
static double angle_diff(double a, double b)
{a -= b;while( a <= -M_PI ) a += M_2__PI;while( a >   M_PI ) a -= M_2__PI;if( a < 0.0 ) a = -a;return a;
}/*----------------------------------------------------------------------------*/
/** Signed angle difference.*/
static double angle_diff_signed(double a, double b)
{a -= b;while( a <= -M_PI ) a += M_2__PI;while( a >   M_PI ) a -= M_2__PI;return a;
}/*----------------------------------------------------------------------------*/
/*----------------------------- NFA computation ------------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** Computes the natural logarithm of the absolute value ofthe gamma function of x using the Lanczos approximation.See http://www.rskey.org/gamma.htmThe formula used is@f[\Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) }(x+5.5)^{x+0.5} e^{-(x+5.5)}@f]so@f[\log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right)+ (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n)@f]andq0 = 75122.6331530,q1 = 80916.6278952,q2 = 36308.2951477,q3 = 8687.24529705,q4 = 1168.92649479,q5 = 83.8676043424,q6 = 2.50662827511.*/
static double log_gamma_lanczos(double x)
{static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,8687.24529705, 1168.92649479, 83.8676043424,2.50662827511 };double a = (x+0.5) * log(x+5.5) - (x+5.5);double b = 0.0;int n;for(n=0;n<7;n++){a -= log( x + (double) n );b += q[n] * pow( x, (double) n );}return a + log(b);
}/*----------------------------------------------------------------------------*/
/** Computes the natural logarithm of the absolute value ofthe gamma function of x using Windschitl method.See http://www.rskey.org/gamma.htmThe formula used is@f[\Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e}\sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x@f]so@f[\log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x+ 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right).@f]This formula is a good approximation when x > 15.*/
static double log_gamma_windschitl(double x)
{return 0.918938533204673 + (x-0.5)*log(x) - x+ 0.5*x*log( x*sinh(1/x) + 1/(810.0*pow(x,6.0)) );
}/*----------------------------------------------------------------------------*/
/** Computes the natural logarithm of the absolute value ofthe gamma function of x. When x>15 use log_gamma_windschitl(),otherwise use log_gamma_lanczos().*/
#define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))/*----------------------------------------------------------------------------*/
/** Size of the table to store already computed inverse values.*/
#define TABSIZE 100000/*----------------------------------------------------------------------------*/
/** Computes -log10(NFA).NFA stands for Number of False Alarms:@f[\mathrm{NFA} = NT \cdot B(n,k,p)@f]- NT       - number of tests- B(n,k,p) - tail of binomial distribution with parameters n,k and p:@f[B(n,k,p) = \sum_{j=k}^n\left(\begin{array}{c}n\\j\end{array}\right)p^{j} (1-p)^{n-j}@f]The value -log10(NFA) is equivalent but more intuitive than NFA:- -1 corresponds to 10 mean false alarms-  0 corresponds to 1 mean false alarm-  1 corresponds to 0.1 mean false alarms-  2 corresponds to 0.01 mean false alarms-  ...Used this way, the bigger the value, better the detection,and a logarithmic scale is used.@param n,k,p binomial parameters.@param logNT logarithm of Number of TestsThe computation is based in the gamma function by the followingrelation:@f[\left(\begin{array}{c}n\\k\end{array}\right)= \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }.@f]We use efficient algorithms to compute the logarithm ofthe gamma function.To make the computation faster, not all the sum is computed, partof the terms are neglected based on a bound to the error obtained(an error of 10% in the result is accepted).*/
static double nfa(int n, int k, double p, double logNT)
{static double inv[TABSIZE];   /* table to keep computed inverse values */double tolerance = 0.1;       /* an error of 10% in the result is accepted */double log1term,term,bin_term,mult_term,bin_tail,err,p_term;int i;/* check parameters */if( n<0 || k<0 || k>n || p<=0.0 || p>=1.0 )error("nfa: wrong n, k or p values.");/* trivial cases */if( n==0 || k==0 ) return -logNT;if( n==k ) return -logNT - (double) n * log10(p);/* probability term */p_term = p / (1.0-p);/* compute the first term of the series *//*binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}where bincoef(n,i) are the binomial coefficients.Butbincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).We use this to compute the first term. Actually the log of it.*/log1term = log_gamma( (double) n + 1.0 ) - log_gamma( (double) k + 1.0 )- log_gamma( (double) (n-k) + 1.0 )+ (double) k * log(p) + (double) (n-k) * log(1.0-p);term = exp(log1term);/* in some cases no more computations are needed */if( double_equal(term,0.0) )              /* the first term is almost zero */{if( (double) k > (double) n * p )     /* at begin or end of the tail?  */return -log1term / M_LN10 - logNT;  /* end: use just the first term  */elsereturn -logNT;                      /* begin: the tail is roughly 1  */}/* compute more terms if needed */bin_tail = term;for(i=k+1;i<=n;i++){/*Asterm_i = bincoef(n,i) * p^i * (1-p)^(n-i)andbincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,then,term_i / term_i-1 = (n-i+1)/i * p/(1-p)andterm_i = term_i-1 * (n-i+1)/i * p/(1-p).1/i is stored in a table as they are computed,because divisions are expensive.p/(1-p) is computed only once and stored in 'p_term'.*/bin_term = (double) (n-i+1) * ( i<TABSIZE ?( inv[i]!=0.0 ? inv[i] : ( inv[i] = 1.0 / (double) i ) ) :1.0 / (double) i );mult_term = bin_term * p_term;term *= mult_term;bin_tail += term;if(bin_term<1.0){/* When bin_term<1 then mult_term_j<mult_term_i for j>i.Then, the error on the binomial tail when truncated atthe i term can be bounded by a geometric series of formterm_i * sum mult_term_i^j.                            */err = term * ( ( 1.0 - pow( mult_term, (double) (n-i+1) ) ) /(1.0-mult_term) - 1.0 );/* One wants an error at most of tolerance*final_result, or:tolerance * abs(-log10(bin_tail)-logNT).Now, the error that can be accepted on bin_tail isgiven by tolerance*final_result divided by the derivativeof -log10(x) when x=bin_tail. that is:tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)Finally, we truncate the tail if the error is less than:tolerance * abs(-log10(bin_tail)-logNT) * bin_tail        */if( err < tolerance * fabs(-log10(bin_tail)-logNT) * bin_tail ) break;}}return -log10(bin_tail) - logNT;
}/*----------------------------------------------------------------------------*/
/*--------------------------- Rectangle structure ----------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** Rectangle structure: line segment with width.*/
struct rect
{double x1,y1,x2,y2;  /* first and second point of the line segment */double width;        /* rectangle width */double x,y;          /* center of the rectangle */double theta;        /* angle */double dx,dy;        /* (dx,dy) is vector oriented as the line segment */double prec;         /* tolerance angle */double p;            /* probability of a point with angle within 'prec' */
};/*----------------------------------------------------------------------------*/
/** Copy one rectangle structure to another.*/
static void rect_copy(struct rect * in, struct rect * out)
{/* check parameters */if( in == NULL || out == NULL ) error("rect_copy: invalid 'in' or 'out'.");/* copy values */out->x1 = in->x1;out->y1 = in->y1;out->x2 = in->x2;out->y2 = in->y2;out->width = in->width;out->x = in->x;out->y = in->y;out->theta = in->theta;out->dx = in->dx;out->dy = in->dy;out->prec = in->prec;out->p = in->p;
}/*----------------------------------------------------------------------------*/
/** Rectangle points iterator.The integer coordinates of pixels inside a rectangle areiteratively explored. This structure keep track of the process andfunctions ri_ini(), ri_inc(), ri_end(), and ri_del() are used inthe process. An example of how to use the iterator is as follows:\codestruct rect * rec = XXX; // some rectanglerect_iter * i;for( i=ri_ini(rec); !ri_end(i); ri_inc(i) ){// your code, using 'i->x' and 'i->y' as coordinates}ri_del(i); // delete iterator\endcodeThe pixels are explored 'column' by 'column', where we call'column' a set of pixels with the same x value that are inside therectangle. The following is an schematic representation of arectangle, the 'column' being explored is marked by colons, andthe current pixel being explored is 'x,y'.\verbatimvx[1],vy[1]*   **       **           **               ye*                :  *vx[0],vy[0]           :     **              :        **          x,y          **        :              **     :            vx[2],vy[2]*  :                *y                     ys              *^                        *           *|                           *       *|                              *   *+---> x                      vx[3],vy[3]\endverbatimThe first 'column' to be explored is the one with the smaller xvalue. Each 'column' is explored starting from the pixel of the'column' (inside the rectangle) with the smallest y value.The four corners of the rectangle are stored in order that rotatesaround the corners at the arrays 'vx[]' and 'vy[]'. The firstpoint is always the one with smaller x value.'x' and 'y' are the coordinates of the pixel being explored. 'ys'and 'ye' are the start and end values of the current column beingexplored. So, 'ys' < 'ye'.*/
typedef struct
{double vx[4];  /* rectangle's corner X coordinates in circular order */double vy[4];  /* rectangle's corner Y coordinates in circular order */double ys,ye;  /* start and end Y values of current 'column' */int x,y;       /* coordinates of currently explored pixel */
} rect_iter;/*----------------------------------------------------------------------------*/
/** Interpolate y value corresponding to 'x' value given, inthe line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the smallerof 'y1' and 'y2'.The following restrictions are required:- x1 <= x2- x1 <= x- x  <= x2*/
static double inter_low(double x, double x1, double y1, double x2, double y2)
{/* check parameters */if( x1 > x2 || x < x1 || x > x2 )error("inter_low: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");/* interpolation */if( double_equal(x1,x2) && y1<y2 ) return y1;if( double_equal(x1,x2) && y1>y2 ) return y2;return y1 + (x-x1) * (y2-y1) / (x2-x1);
}/*----------------------------------------------------------------------------*/
/** Interpolate y value corresponding to 'x' value given, inthe line 'x1,y1' to 'x2,y2'; if 'x1=x2' return the largerof 'y1' and 'y2'.The following restrictions are required:- x1 <= x2- x1 <= x- x  <= x2*/
static double inter_hi(double x, double x1, double y1, double x2, double y2)
{/* check parameters */if( x1 > x2 || x < x1 || x > x2 )error("inter_hi: unsuitable input, 'x1>x2' or 'x<x1' or 'x>x2'.");/* interpolation */if( double_equal(x1,x2) && y1<y2 ) return y2;if( double_equal(x1,x2) && y1>y2 ) return y1;return y1 + (x-x1) * (y2-y1) / (x2-x1);
}/*----------------------------------------------------------------------------*/
/** Free memory used by a rectangle iterator.*/
static void ri_del(rect_iter * iter)
{if( iter == NULL ) error("ri_del: NULL iterator.");free( (void *) iter );
}/*----------------------------------------------------------------------------*/
/** Check if the iterator finished the full iteration.See details in \ref rect_iter*/
static int ri_end(rect_iter * i)
{/* check input */if( i == NULL ) error("ri_end: NULL iterator.");/* if the current x value is larger than the largestx value in the rectangle (vx[2]), we know the fullexploration of the rectangle is finished. */return (double)(i->x) > i->vx[2];
}/*----------------------------------------------------------------------------*/
/** Increment a rectangle iterator.See details in \ref rect_iter*/
static void ri_inc(rect_iter * i)
{/* check input */if( i == NULL ) error("ri_inc: NULL iterator.");/* if not at end of exploration,increase y value for next pixel in the 'column' */if( !ri_end(i) ) i->y++;/* if the end of the current 'column' is reached,and it is not the end of exploration,advance to the next 'column' */while( (double) (i->y) > i->ye && !ri_end(i) ){/* increase x, next 'column' */i->x++;/* if end of exploration, return */if( ri_end(i) ) return;/* update lower y limit (start) for the new 'column'.We need to interpolate the y value that corresponds to thelower side of the rectangle. The first thing is to decide ifthe corresponding side isvx[0],vy[0] to vx[3],vy[3] orvx[3],vy[3] to vx[2],vy[2]Then, the side is interpolated for the x value of the'column'. But, if the side is vertical (as it could happen ifthe rectangle is vertical and we are dealing with the firstor last 'columns') then we pick the lower value of the sideby using 'inter_low'.*/if( (double) i->x < i->vx[3] )i->ys = inter_low((double)i->x,i->vx[0],i->vy[0],i->vx[3],i->vy[3]);elsei->ys = inter_low((double)i->x,i->vx[3],i->vy[3],i->vx[2],i->vy[2]);/* update upper y limit (end) for the new 'column'.We need to interpolate the y value that corresponds to theupper side of the rectangle. The first thing is to decide ifthe corresponding side isvx[0],vy[0] to vx[1],vy[1] orvx[1],vy[1] to vx[2],vy[2]Then, the side is interpolated for the x value of the'column'. But, if the side is vertical (as it could happen ifthe rectangle is vertical and we are dealing with the firstor last 'columns') then we pick the lower value of the sideby using 'inter_low'.*/if( (double)i->x < i->vx[1] )i->ye = inter_hi((double)i->x,i->vx[0],i->vy[0],i->vx[1],i->vy[1]);elsei->ye = inter_hi((double)i->x,i->vx[1],i->vy[1],i->vx[2],i->vy[2]);/* new y */i->y = (int) ceil(i->ys);}
}/*----------------------------------------------------------------------------*/
/** Create and initialize a rectangle iterator.See details in \ref rect_iter*/
static rect_iter * ri_ini(struct rect * r)
{double vx[4],vy[4];int n,offset;rect_iter * i;/* check parameters */if( r == NULL ) error("ri_ini: invalid rectangle.");/* get memory */i = (rect_iter *) malloc(sizeof(rect_iter));if( i == NULL ) error("ri_ini: Not enough memory.");/* build list of rectangle corners orderedin a circular way around the rectangle */vx[0] = r->x1 - r->dy * r->width / 2.0;vy[0] = r->y1 + r->dx * r->width / 2.0;vx[1] = r->x2 - r->dy * r->width / 2.0;vy[1] = r->y2 + r->dx * r->width / 2.0;vx[2] = r->x2 + r->dy * r->width / 2.0;vy[2] = r->y2 - r->dx * r->width / 2.0;vx[3] = r->x1 + r->dy * r->width / 2.0;vy[3] = r->y1 - r->dx * r->width / 2.0;/* compute rotation of index of corners needed so that the firstpoint has the smaller x.if one side is vertical, thus two corners have the same smaller xvalue, the one with the largest y value is selected as the first.*/if( r->x1 < r->x2 && r->y1 <= r->y2 ) offset = 0;else if( r->x1 >= r->x2 && r->y1 < r->y2 ) offset = 1;else if( r->x1 > r->x2 && r->y1 >= r->y2 ) offset = 2;else offset = 3;/* apply rotation of index. */for(n=0; n<4; n++){i->vx[n] = vx[(offset+n)%4];i->vy[n] = vy[(offset+n)%4];}/* Set an initial condition.The values are set to values that will cause 'ri_inc' (that willbe called immediately) to initialize correctly the first 'column'and compute the limits 'ys' and 'ye'.'y' is set to the integer value of vy[0], the starting corner.'ys' and 'ye' are set to very small values, so 'ri_inc' willnotice that it needs to start a new 'column'.The smallest integer coordinate inside of the rectangle is'ceil(vx[0])'. The current 'x' value is set to that value minusone, so 'ri_inc' (that will increase x by one) will advance tothe first 'column'.*/i->x = (int) ceil(i->vx[0]) - 1;i->y = (int) ceil(i->vy[0]);i->ys = i->ye = -DBL_MAX;/* advance to the first pixel */ri_inc(i);return i;
}/*----------------------------------------------------------------------------*/
/** Compute a rectangle's NFA value.*/
static double rect_nfa(struct rect * rec, image_double angles, double logNT)
{rect_iter * i;int pts = 0;int alg = 0;/* check parameters */if( rec == NULL ) error("rect_nfa: invalid rectangle.");if( angles == NULL ) error("rect_nfa: invalid 'angles'.");/* compute the total number of pixels and of aligned points in 'rec' */for(i=ri_ini(rec); !ri_end(i); ri_inc(i)) /* rectangle iterator */if( i->x >= 0 && i->y >= 0 &&i->x < (int) angles->xsize && i->y < (int) angles->ysize ){++pts; /* total number of pixels counter */if( isaligned(i->x, i->y, angles, rec->theta, rec->prec) )++alg; /* aligned points counter */}ri_del(i); /* delete iterator */return nfa(pts,alg,rec->p,logNT); /* compute NFA value */
}/*----------------------------------------------------------------------------*/
/*---------------------------------- Regions ---------------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** Compute region's angle as the principal inertia axis of the region.The following is the region inertia matrix A:@f[A = \left(\begin{array}{cc}Ixx & Ixy \\Ixy & Iyy \\\end{array}\right)@f]whereIxx =   sum_i G(i).(y_i - cx)^2Iyy =   sum_i G(i).(x_i - cy)^2Ixy = - sum_i G(i).(x_i - cx).(y_i - cy)and- G(i) is the gradient norm at pixel i, used as pixel's weight.- x_i and y_i are the coordinates of pixel i.- cx and cy are the coordinates of the center of th region.lambda1 and lambda2 are the eigenvalues of matrix A,with lambda1 >= lambda2. They are found by solving thecharacteristic polynomial:det( lambda I - A) = 0that gives:lambda1 = ( Ixx + Iyy + sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2lambda2 = ( Ixx + Iyy - sqrt( (Ixx-Iyy)^2 + 4.0*Ixy*Ixy) ) / 2To get the line segment direction we want to get the angle theeigenvector associated to the smallest eigenvalue. We haveto solve for a,b in:a.Ixx + b.Ixy = a.lambda2a.Ixy + b.Iyy = b.lambda2We want the angle theta = atan(b/a). It can be computed withany of the two equations:theta = atan( (lambda2-Ixx) / Ixy )ortheta = atan( Ixy / (lambda2-Iyy) )When |Ixx| > |Iyy| we use the first, otherwise the second (just toget better numeric precision).*/
static double get_theta( struct point * reg, int reg_size, double x, double y,image_double modgrad, double reg_angle, double prec )
{double lambda,theta,weight;double Ixx = 0.0;double Iyy = 0.0;double Ixy = 0.0;int i;/* check parameters */if( reg == NULL ) error("get_theta: invalid region.");if( reg_size <= 1 ) error("get_theta: region size <= 1.");if( modgrad == NULL || modgrad->data == NULL )error("get_theta: invalid 'modgrad'.");if( prec < 0.0 ) error("get_theta: 'prec' must be positive.");/* compute inertia matrix */for(i=0; i<reg_size; i++){weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];Ixx += ( (double) reg[i].y - y ) * ( (double) reg[i].y - y ) * weight;Iyy += ( (double) reg[i].x - x ) * ( (double) reg[i].x - x ) * weight;Ixy -= ( (double) reg[i].x - x ) * ( (double) reg[i].y - y ) * weight;}if( double_equal(Ixx,0.0) && double_equal(Iyy,0.0) && double_equal(Ixy,0.0) )error("get_theta: null inertia matrix.");/* compute smallest eigenvalue */lambda = 0.5 * ( Ixx + Iyy - sqrt( (Ixx-Iyy)*(Ixx-Iyy) + 4.0*Ixy*Ixy ) );/* compute angle */theta = fabs(Ixx)>fabs(Iyy) ? atan2(lambda-Ixx,Ixy) : atan2(Ixy,lambda-Iyy);/* The previous procedure doesn't cares about orientation,so it could be wrong by 180 degrees. Here is corrected if necessary. */if( angle_diff(theta,reg_angle) > prec ) theta += M_PI;return theta;
}/*----------------------------------------------------------------------------*/
/** Computes a rectangle that covers a region of points.*/
static void region2rect( struct point * reg, int reg_size,image_double modgrad, double reg_angle,double prec, double p, struct rect * rec )
{double x,y,dx,dy,l,w,theta,weight,sum,l_min,l_max,w_min,w_max;int i;/* check parameters */if( reg == NULL ) error("region2rect: invalid region.");if( reg_size <= 1 ) error("region2rect: region size <= 1.");if( modgrad == NULL || modgrad->data == NULL )error("region2rect: invalid image 'modgrad'.");if( rec == NULL ) error("region2rect: invalid 'rec'.");/* center of the region:It is computed as the weighted sum of the coordinatesof all the pixels in the region. The norm of the gradientis used as the weight of a pixel. The sum is as follows:cx = \sum_i G(i).x_icy = \sum_i G(i).y_iwhere G(i) is the norm of the gradient of pixel iand x_i,y_i are its coordinates.*/x = y = sum = 0.0;for(i=0; i<reg_size; i++){weight = modgrad->data[ reg[i].x + reg[i].y * modgrad->xsize ];x += (double) reg[i].x * weight;y += (double) reg[i].y * weight;sum += weight;}if( sum <= 0.0 ) error("region2rect: weights sum equal to zero.");x /= sum;y /= sum;/* theta */theta = get_theta(reg,reg_size,x,y,modgrad,reg_angle,prec);/* length and width:'l' and 'w' are computed as the distance from the center of theregion to pixel i, projected along the rectangle axis (dx,dy) andto the orthogonal axis (-dy,dx), respectively.The length of the rectangle goes from l_min to l_max, where l_minand l_max are the minimum and maximum values of l in the region.Analogously, the width is selected from w_min to w_max, wherew_min and w_max are the minimum and maximum of w for the pixelsin the region.*/dx = cos(theta);dy = sin(theta);l_min = l_max = w_min = w_max = 0.0;for(i=0; i<reg_size; i++){l =  ( (double) reg[i].x - x) * dx + ( (double) reg[i].y - y) * dy;w = -( (double) reg[i].x - x) * dy + ( (double) reg[i].y - y) * dx;if( l > l_max ) l_max = l;if( l < l_min ) l_min = l;if( w > w_max ) w_max = w;if( w < w_min ) w_min = w;}/* store values */rec->x1 = x + l_min * dx;rec->y1 = y + l_min * dy;rec->x2 = x + l_max * dx;rec->y2 = y + l_max * dy;rec->width = w_max - w_min;rec->x = x;rec->y = y;rec->theta = theta;rec->dx = dx;rec->dy = dy;rec->prec = prec;rec->p = p;/* we impose a minimal width of one pixelA sharp horizontal or vertical step would produce a perfectlyhorizontal or vertical region. The width computed would bezero. But that corresponds to a one pixels width transition inthe image.*/if( rec->width < 1.0 ) rec->width = 1.0;
}/*----------------------------------------------------------------------------*/
/** Build a region of pixels that share the same angle, up to atolerance 'prec', starting at point (x,y).*/
static void region_grow( int x, int y, image_double angles, struct point * reg,int * reg_size, double * reg_angle, image_char used,double prec )
{double sumdx,sumdy;int xx,yy,i;/* check parameters */if( x < 0 || y < 0 || x >= (int) angles->xsize || y >= (int) angles->ysize )error("region_grow: (x,y) out of the image.");if( angles == NULL || angles->data == NULL )error("region_grow: invalid image 'angles'.");if( reg == NULL ) error("region_grow: invalid 'reg'.");if( reg_size == NULL ) error("region_grow: invalid pointer 'reg_size'.");if( reg_angle == NULL ) error("region_grow: invalid pointer 'reg_angle'.");if( used == NULL || used->data == NULL )error("region_grow: invalid image 'used'.");/* first point of the region */*reg_size = 1;reg[0].x = x;reg[0].y = y;*reg_angle = angles->data[x+y*angles->xsize];  /* region's angle */sumdx = cos(*reg_angle);sumdy = sin(*reg_angle);used->data[x+y*used->xsize] = USED;/* try neighbors as new region points */for(i=0; i<*reg_size; i++)for(xx=reg[i].x-1; xx<=reg[i].x+1; xx++)for(yy=reg[i].y-1; yy<=reg[i].y+1; yy++)if( xx>=0 && yy>=0 && xx<(int)used->xsize && yy<(int)used->ysize &&used->data[xx+yy*used->xsize] != USED &&isaligned(xx,yy,angles,*reg_angle,prec) ){/* add point */used->data[xx+yy*used->xsize] = USED;reg[*reg_size].x = xx;reg[*reg_size].y = yy;++(*reg_size);/* update region's angle */sumdx += cos( angles->data[xx+yy*angles->xsize] );sumdy += sin( angles->data[xx+yy*angles->xsize] );*reg_angle = atan2(sumdy,sumdx);}
}/*----------------------------------------------------------------------------*/
/** Try some rectangles variations to improve NFA value. Only if therectangle is not meaningful (i.e., log_nfa <= log_eps).*/
static double rect_improve( struct rect * rec, image_double angles,double logNT, double log_eps )
{struct rect r;double log_nfa,log_nfa_new;double delta = 0.5;double delta_2 = delta / 2.0;int n;log_nfa = rect_nfa(rec,angles,logNT);if( log_nfa > log_eps ) return log_nfa;/* try finer precisions */rect_copy(rec,&r);for(n=0; n<5; n++){r.p /= 2.0;r.prec = r.p * M_PI;log_nfa_new = rect_nfa(&r,angles,logNT);if( log_nfa_new > log_nfa ){log_nfa = log_nfa_new;rect_copy(&r,rec);}}if( log_nfa > log_eps ) return log_nfa;/* try to reduce width */rect_copy(rec,&r);for(n=0; n<5; n++){if( (r.width - delta) >= 0.5 ){r.width -= delta;log_nfa_new = rect_nfa(&r,angles,logNT);if( log_nfa_new > log_nfa ){rect_copy(&r,rec);log_nfa = log_nfa_new;}}}if( log_nfa > log_eps ) return log_nfa;/* try to reduce one side of the rectangle */rect_copy(rec,&r);for(n=0; n<5; n++){if( (r.width - delta) >= 0.5 ){r.x1 += -r.dy * delta_2;r.y1 +=  r.dx * delta_2;r.x2 += -r.dy * delta_2;r.y2 +=  r.dx * delta_2;r.width -= delta;log_nfa_new = rect_nfa(&r,angles,logNT);if( log_nfa_new > log_nfa ){rect_copy(&r,rec);log_nfa = log_nfa_new;}}}if( log_nfa > log_eps ) return log_nfa;/* try to reduce the other side of the rectangle */rect_copy(rec,&r);for(n=0; n<5; n++){if( (r.width - delta) >= 0.5 ){r.x1 -= -r.dy * delta_2;r.y1 -=  r.dx * delta_2;r.x2 -= -r.dy * delta_2;r.y2 -=  r.dx * delta_2;r.width -= delta;log_nfa_new = rect_nfa(&r,angles,logNT);if( log_nfa_new > log_nfa ){rect_copy(&r,rec);log_nfa = log_nfa_new;}}}if( log_nfa > log_eps ) return log_nfa;/* try even finer precisions */rect_copy(rec,&r);for(n=0; n<5; n++){r.p /= 2.0;r.prec = r.p * M_PI;log_nfa_new = rect_nfa(&r,angles,logNT);if( log_nfa_new > log_nfa ){log_nfa = log_nfa_new;rect_copy(&r,rec);}}return log_nfa;
}/*----------------------------------------------------------------------------*/
/** Reduce the region size, by elimination the points far from thestarting point, until that leads to rectangle with the rightdensity of region points or to discard the region if too small.*/
static int reduce_region_radius( struct point * reg, int * reg_size,image_double modgrad, double reg_angle,double prec, double p, struct rect * rec,image_char used, image_double angles,double density_th )
{double density,rad1,rad2,rad,xc,yc;int i;/* check parameters */if( reg == NULL ) error("reduce_region_radius: invalid pointer 'reg'.");if( reg_size == NULL )error("reduce_region_radius: invalid pointer 'reg_size'.");if( prec < 0.0 ) error("reduce_region_radius: 'prec' must be positive.");if( rec == NULL ) error("reduce_region_radius: invalid pointer 'rec'.");if( used == NULL || used->data == NULL )error("reduce_region_radius: invalid image 'used'.");if( angles == NULL || angles->data == NULL )error("reduce_region_radius: invalid image 'angles'.");/* compute region points density */density = (double) *reg_size /( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );/* if the density criterion is satisfied there is nothing to do */if( density >= density_th ) return TRUE;/* compute region's radius */xc = (double) reg[0].x;yc = (double) reg[0].y;rad1 = dist( xc, yc, rec->x1, rec->y1 );rad2 = dist( xc, yc, rec->x2, rec->y2 );rad = rad1 > rad2 ? rad1 : rad2;/* while the density criterion is not satisfied, remove farther pixels */while( density < density_th ){rad *= 0.75; /* reduce region's radius to 75% of its value *//* remove points from the region and update 'used' map */for(i=0; i<*reg_size; i++)if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) > rad ){/* point not kept, mark it as NOTUSED */used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;/* remove point from the region */reg[i].x = reg[*reg_size-1].x; /* if i==*reg_size-1 copy itself */reg[i].y = reg[*reg_size-1].y;--(*reg_size);--i; /* to avoid skipping one point */}/* reject if the region is too small.2 is the minimal region size for 'region2rect' to work. */if( *reg_size < 2 ) return FALSE;/* re-compute rectangle */region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);/* re-compute region points density */density = (double) *reg_size /( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );}/* if this point is reached, the density criterion is satisfied */return TRUE;
}/*----------------------------------------------------------------------------*/
/** Refine a rectangle.For that, an estimation of the angle tolerance is performed by thestandard deviation of the angle at points near the region'sstarting point. Then, a new region is grown starting from the samepoint, but using the estimated angle tolerance. If this fails toproduce a rectangle with the right density of region points,'reduce_region_radius' is called to try to satisfy this condition.*/
static int refine( struct point * reg, int * reg_size, image_double modgrad,double reg_angle, double prec, double p, struct rect * rec,image_char used, image_double angles, double density_th )
{double angle,ang_d,mean_angle,tau,density,xc,yc,ang_c,sum,s_sum;int i,n;/* check parameters */if( reg == NULL ) error("refine: invalid pointer 'reg'.");if( reg_size == NULL ) error("refine: invalid pointer 'reg_size'.");if( prec < 0.0 ) error("refine: 'prec' must be positive.");if( rec == NULL ) error("refine: invalid pointer 'rec'.");if( used == NULL || used->data == NULL )error("refine: invalid image 'used'.");if( angles == NULL || angles->data == NULL )error("refine: invalid image 'angles'.");/* compute region points density */density = (double) *reg_size /( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );/* if the density criterion is satisfied there is nothing to do */if( density >= density_th ) return TRUE;/*------ First try: reduce angle tolerance ------*//* compute the new mean angle and tolerance */xc = (double) reg[0].x;yc = (double) reg[0].y;ang_c = angles->data[ reg[0].x + reg[0].y * angles->xsize ];sum = s_sum = 0.0;n = 0;for(i=0; i<*reg_size; i++){used->data[ reg[i].x + reg[i].y * used->xsize ] = NOTUSED;if( dist( xc, yc, (double) reg[i].x, (double) reg[i].y ) < rec->width ){angle = angles->data[ reg[i].x + reg[i].y * angles->xsize ];ang_d = angle_diff_signed(angle,ang_c);sum += ang_d;s_sum += ang_d * ang_d;++n;}}mean_angle = sum / (double) n;tau = 2.0 * sqrt( (s_sum - 2.0 * mean_angle * sum) / (double) n+ mean_angle*mean_angle ); /* 2 * standard deviation *//* find a new region from the same starting point and new angle tolerance */region_grow(reg[0].x,reg[0].y,angles,reg,reg_size,&reg_angle,used,tau);/* if the region is too small, reject */if( *reg_size < 2 ) return FALSE;/* re-compute rectangle */region2rect(reg,*reg_size,modgrad,reg_angle,prec,p,rec);/* re-compute region points density */density = (double) *reg_size /( dist(rec->x1,rec->y1,rec->x2,rec->y2) * rec->width );/*------ Second try: reduce region radius ------*/if( density < density_th )return reduce_region_radius( reg, reg_size, modgrad, reg_angle, prec, p,rec, used, angles, density_th );/* if this point is reached, the density criterion is satisfied */return TRUE;
}/*----------------------------------------------------------------------------*/
/*-------------------------- Line Segment Detector ---------------------------*/
/*----------------------------------------------------------------------------*//*----------------------------------------------------------------------------*/
/** LSD full interface.*/
double * LineSegmentDetection( int * n_out,double * img, int X, int Y,double scale, double sigma_scale, double quant,double ang_th, double log_eps, double density_th,int n_bins,int ** reg_img, int * reg_x, int * reg_y )
{image_double image;ntuple_list out = new_ntuple_list(7);double * return_value;image_double scaled_image,angles,modgrad;image_char used;image_int region = NULL;struct coorlist * list_p;void * mem_p;struct rect rec;struct point * reg;int reg_size,min_reg_size,i;unsigned int xsize,ysize;double rho,reg_angle,prec,p,log_nfa,logNT;int ls_count = 0;                   /* line segments are numbered 1,2,3,... *//* check parameters */if( img == NULL || X <= 0 || Y <= 0 ) error("invalid image input.");if( scale <= 0.0 ) error("'scale' value must be positive.");if( sigma_scale <= 0.0 ) error("'sigma_scale' value must be positive.");if( quant < 0.0 ) error("'quant' value must be positive.");if( ang_th <= 0.0 || ang_th >= 180.0 )error("'ang_th' value must be in the range (0,180).");if( density_th < 0.0 || density_th > 1.0 )error("'density_th' value must be in the range [0,1].");if( n_bins <= 0 ) error("'n_bins' value must be positive.");/* angle tolerance */prec = M_PI * ang_th / 180.0;p = ang_th / 180.0;rho = quant / sin(prec); /* gradient magnitude threshold *//* load and scale image (if necessary) and compute angle at each pixel */image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );if( scale != 1.0 ){scaled_image = gaussian_sampler( image, scale, sigma_scale );angles = ll_angle( scaled_image, rho, &list_p, &mem_p,&modgrad, (unsigned int) n_bins );free_image_double(scaled_image);}elseangles = ll_angle( image, rho, &list_p, &mem_p, &modgrad,(unsigned int) n_bins );xsize = angles->xsize;ysize = angles->ysize;/* Number of Tests - NTThe theoretical number of tests is Np.(XY)^(5/2)where X and Y are number of columns and rows of the image.Np corresponds to the number of angle precisions considered.As the procedure 'rect_improve' tests 5 times to halve theangle precision, and 5 more times after improving other factors,11 different precision values are potentially tested. Thus,the number of tests is11 * (X*Y)^(5/2)whose logarithm value islog10(11) + 5/2 * (log10(X) + log10(Y)).*/logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0+ log10(11.0);min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in regionthat can give a meaningful event *//* initialize some structures */if( reg_img != NULL && reg_x != NULL && reg_y != NULL ) /* save region data */region = new_image_int_ini(angles->xsize,angles->ysize,0);used = new_image_char_ini(xsize,ysize,NOTUSED);reg = (struct point *) calloc( (size_t) (xsize*ysize), sizeof(struct point) );if( reg == NULL ) error("not enough memory!");/* search for line segments */for(; list_p != NULL; list_p = list_p->next )if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )/* there is no risk of double comparison problems herebecause we are only interested in the exact NOTDEF value */{/* find the region of connected point and ~equal angle */region_grow( list_p->x, list_p->y, angles, reg, &reg_size,&reg_angle, used, prec );/* reject small regions */if( reg_size < min_reg_size ) continue;/* construct rectangular approximation for the region */region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);/* Check if the rectangle exceeds the minimal density ofregion points. If not, try to improve the region.The rectangle will be rejected if the final one doesnot fulfill the minimal density condition.This is an addition to the original LSD algorithm published in"LSD: A Fast Line Segment Detector with a False Detection Control"by R. Grompone von Gioi, J. Jakubowicz, J.M. Morel, and G. Randall.The original algorithm is obtained with density_th = 0.0.*/if( !refine( reg, &reg_size, modgrad, reg_angle,prec, p, &rec, used, angles, density_th ) ) continue;/* compute NFA value */log_nfa = rect_improve(&rec,angles,logNT,log_eps);if( log_nfa <= log_eps ) continue;/* A New Line Segment was found! */++ls_count;  /* increase line segment counter *//*The gradient was computed with a 2x2 mask, its value corresponds topoints with an offset of (0.5,0.5), that should be added to output.The coordinates origin is at the center of pixel (0,0).*/rec.x1 += 0.5; rec.y1 += 0.5;rec.x2 += 0.5; rec.y2 += 0.5;/* scale the result values if a subsampling was performed */if( scale != 1.0 ){rec.x1 /= scale; rec.y1 /= scale;rec.x2 /= scale; rec.y2 /= scale;rec.width /= scale;}/* add line segment found to output */add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,rec.width, rec.p, log_nfa );/* add region number to 'region' image if needed */if( region != NULL )for(i=0; i<reg_size; i++)region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;}/* free memory */free( (void *) image );   /* only the double_image structure should be freed,the data pointer was provided to this functionsand should not be destroyed.                 */free_image_double(angles);free_image_double(modgrad);free_image_char(used);free( (void *) reg );free( (void *) mem_p );/* return the result */if( reg_img != NULL && reg_x != NULL && reg_y != NULL ){if( region == NULL ) error("'region' should be a valid image.");*reg_img = region->data;if( region->xsize > (unsigned int) INT_MAX ||region->xsize > (unsigned int) INT_MAX )error("region image to big to fit in INT sizes.");*reg_x = (int) (region->xsize);*reg_y = (int) (region->ysize);/* free the 'region' structure.we cannot use the function 'free_image_int' because we need to keepthe memory with the image data to be returned by this function. */free( (void *) region );}if( out->size > (unsigned int) INT_MAX )error("too many detections to fit in an INT.");*n_out = (int) (out->size);return_value = out->values;free( (void *) out );  /* only the 'ntuple_list' structure must be freed,but the 'values' pointer must be keep to returnas a result. */return return_value;
}/*----------------------------------------------------------------------------*/
/** LSD Simple Interface with Scale and Region output.*/
double * lsd_scale_region( int * n_out,double * img, int X, int Y, double scale,int ** reg_img, int * reg_x, int * reg_y )
{/* LSD parameters */double sigma_scale = 0.6; /* Sigma for Gaussian filter is computed assigma = sigma_scale/scale.                    */double quant = 2.0;       /* Bound to the quantization error on thegradient norm.                                */double ang_th = 22.5;     /* Gradient angle tolerance in degrees.           */double log_eps = 0.0;     /* Detection threshold: -log10(NFA) > log_eps     */double density_th = 0.7;  /* Minimal density of region points in rectangle. */int n_bins = 1024;        /* Number of bins in pseudo-ordering of gradientmodulus.                                       */return LineSegmentDetection( n_out, img, X, Y, scale, sigma_scale, quant,ang_th, log_eps, density_th, n_bins,reg_img, reg_x, reg_y );
}/*----------------------------------------------------------------------------*/
/** LSD Simple Interface with Scale.*/
double * lsd_scale(int * n_out, double * img, int X, int Y, double scale)
{return lsd_scale_region(n_out,img,X,Y,scale,NULL,NULL,NULL);
}/*----------------------------------------------------------------------------*/
/** LSD Simple Interface.*/
double * lsd(int * n_out, double * img, int X, int Y)
{/* LSD parameters */double scale = 0.8;       /* Scale the image by Gaussian filter to 'scale'. */return lsd_scale(n_out,img,X,Y,scale);
}
/*----------------------------------------------------------------------------*/

LSD（线段检测测试文件）相关推荐

1. 目标检测之线段检测---lsd line segment detector

   (1)线段检测应用背景 (2)线段检测原理简介 (3)线段检测实例 a line segment detector (4)hough 变换和 lsd 的区别 --------------------- ...

2. 实验四：文件状态测试--动态检测指定文件的状态信息，当文件的大小发生改变时，给出提示信息，并继续前进检测。当文件的大小的变化次数或持续检查无变化次数达到一定值时，退出检查--操作系统原理和实践

   实验目的 熟悉UNIX的基本SHELL程序设计方法,包括: 命令行参数检测 变量设置 文件状态检测与特定信息读取 程序运行控制 实验内容 编写一个SHELL程序,动态检测指定文件的状态信息,当文件的大 ...

3. 武大+CMU最新开源！全面支持平面/鱼眼/球面相机的实时统一线段检测算法

   点击上方"3D视觉工坊",选择"星标" 干货第一时间送达 来源丨计算机视觉life 作者丨司晗骞 大家好,今天给大家介绍一篇武汉大学和卡内基梅隆大学联合发布的论 ...

4. 文档扫描识别——基于M-LSD线段检测的拍照文档校正

   前言 1.拍照文档扫描识别是办公类App里面最常用到的的一类应用,市面上有很多相关的App,及主要技术点有几个要用到图像处理,有边缘检测校正,文档滤镜,和OCR. 2.关于边缘文档连续检测,有用传统算 ...

5. 基于线条特征的机场检测算法——LSD直线检测算法、平行线组提取和聚类

   遥感图像的机场检测是图像处理在军事以及航空领域一个重要的应用,现有一些机场提取方法利用显著性特征获取机场区域的方法容易使得机场提取不够完整,而且会混入过多的虚警区域,原因在于图像的显著性特征并能用来表 ...

6. 测试用html文件是否存在,ASP如何检测某文件夹是否存在,不存在则自动创建

   直接给大家分享一下脚本之家测试正常可以使用的代码,并且支持多级目录创建 代码一 Function CreateMultiFolder(ByVal CFolder) Dim objFSO, PhCrea ...

7. LSD直线检测和霍夫线变换的学习建议

   最近笔者学习霍夫线变换和LSD直线检测算法,有一些学习建议,希望可以给予大家一些帮助.  学习霍夫变换的感想 每个人理解的霍夫变换或许略有差异,但是最主要的是笛卡尔坐标系跟极坐标系的相互转换. 霍夫变 ...

8. 渗透测试-文件上传/下载/包含

   渗透测试-文件上传/下载/包含 概述 一.常见校验上传文件的方法 客户端校验 服务器端校验 1.校验请求头 content-type字段 2. 正则匹配来判断文件幻数(文件头)内容是否符合要求 3. ...

9. C# OpenCV EmguCV LSD直线检测使用Demo

   点击下方卡片,关注"OpenCV与AI深度学习"公众号! 视觉/图像重磅干货,第一时间送达! LSD直线检测原理大家可以自行百度查询,这里给出EmguCV4.5.4使用Demo如下 ...

10. 用Python实现LSD直线检测

    LSD (Line Segment Detector) 是一种用于检测图像中直线段的算法. 要在 Python 中实现 LSD 直线检测,首先需要安装 OpenCV 库.OpenCV 是一个开源的计算 ...

最新文章

1. 从Ops到NoOps，阿里文娱智能运维的关键：自动化应用容量管理
2. linux按时间排序并查看发现,linux下扫描文件并按时间排序
3. 视觉计算理论简介【转】
4. 光纤通信是如何接入网络的？
5. 启动器和选择器学习-----（5）启动器
6. Laravel新建对象的方法：make resolve 辅助函数app()
7. 用css3实现ps蒙版效果+动画
8. idea命令行运行多个客户端_推荐一款神仙颜值的 Redis 客户端工具，开源啦
9. python post form data_python实现发送form-data数据的方法详解
10. 华为云教你7天玩转电商应用性能调优，课程免费速来报名！
11. PHP整站迁移空间,Discuz! X2.5 整站搬家迁移升级教程
12. 数据结构与算法之递推算法 C++与PHP实现
13. 暗黑2魔电西格玛攻略_听说你想刷爆怪物的狗头？或许可以试试《暗黑破坏神》类单机手游...
14. Java基础Arrays类
15. 【2017 United Kingdom and Ireland Programming Contest (UKIEPC 2017)】Knightsbridge Rises【最大流+路径输出】
16. 计算机共享文件夹拒绝访问权限,设置共享文件夹访问权限 拒绝访问的方法
17. 视频录制软件有哪些？4款录制视频软件，免费下载
18. java for冒号_浅谈对Java双冒号::的理解
19. 零基础怎么学习平面设计*
20. Android 布局圆角方案总结

热门文章

1. 小内存电脑安装linux,在低内存的情况下安装CentOS系统的技巧
2. Reverb详细大解说
3. 科学养生：揭秘世界上最健康的作息时间表
4. 10 06 01 繁杂
5. MATLAB 的 colormap 函数详解
6. 计算机不定时黑屏,宏基acer 4736ZG不定时黑屏，时亮时不亮通病维修
7. 快乐去学习「快乐机器学习」
8. eNB、gNB、en-gNB和ng-eNB的区别
9. Linux系统监视与进程管理
10. [日推荐]『车主码』解决临时停车、请人挪车的小麻烦