From 9a52f4626bcfa6247aab8f346e2d294d7609936e Mon Sep 17 00:00:00 2001 From: biagio montesano Date: Thu, 21 Aug 2014 17:43:25 +0200 Subject: [PATCH] Created nested classes. Removed some useless headers. --- .../opencv2/line_descriptor/array32.hpp | 101 -- .../opencv2/line_descriptor/bucket_group.hpp | 76 - .../opencv2/line_descriptor/descriptor.hpp | 1029 ++++++++++-- .../line_descriptor/ed_line_detector.hpp | 476 ------ .../line_descriptor/line_structure.hpp | 123 -- .../opencv2/line_descriptor/mihasher.hpp | 130 -- .../line_descriptor/sparse_hashtable.hpp | 89 - .../samples/radius_matching.cpp | 20 +- modules/line_descriptor/src/array32.cpp | 189 --- .../line_descriptor/src/binary_descriptor.cpp | 1414 +++++++++++++++- .../src/binary_descriptor_matcher.cpp | 510 +++++- .../line_descriptor => src}/bitarray.hpp | 0 .../line_descriptor => src}/bitops.hpp | 8 + modules/line_descriptor/src/bucket_group.cpp | 100 -- .../line_descriptor/src/ed_line_detector.cpp | 1432 ----------------- modules/line_descriptor/src/mihasher.cpp | 270 ---- modules/line_descriptor/src/precomp.hpp | 4 + .../line_descriptor/src/sparse_hashtable.cpp | 85 - .../opencv2/line_descriptor => src}/types.hpp | 0 19 files changed, 2802 insertions(+), 3254 deletions(-) delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/array32.hpp delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/bucket_group.hpp delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/ed_line_detector.hpp delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/line_structure.hpp delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/mihasher.hpp delete mode 100644 modules/line_descriptor/include/opencv2/line_descriptor/sparse_hashtable.hpp delete mode 100644 modules/line_descriptor/src/array32.cpp rename modules/line_descriptor/{include/opencv2/line_descriptor => src}/bitarray.hpp (100%) rename modules/line_descriptor/{include/opencv2/line_descriptor => src}/bitops.hpp (98%) delete mode 100644 modules/line_descriptor/src/bucket_group.cpp delete mode 100644 modules/line_descriptor/src/ed_line_detector.cpp delete mode 100644 modules/line_descriptor/src/mihasher.cpp delete mode 100644 modules/line_descriptor/src/sparse_hashtable.cpp rename modules/line_descriptor/{include/opencv2/line_descriptor => src}/types.hpp (100%) diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/array32.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/array32.hpp deleted file mode 100644 index 9a90d13db..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/array32.hpp +++ /dev/null @@ -1,101 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#ifndef __OPENCV_ARRAY32_HPP -#define __OPENCV_ARRAY32_HPP - -#ifdef _WIN32 -#pragma warning( disable : 4267 ) -#endif - -#include "types.hpp" - -class Array32 -{ - - private: - static double ARRAY_RESIZE_FACTOR; - static double ARRAY_RESIZE_ADD_FACTOR; - - public: - /* set ARRAY_RESIZE_FACTOR */ - static void setArrayResizeFactor( double arf ); - - /* constructor */ - Array32(); - - /* destructor */ - ~Array32(); - - /* cleaning function used in destructor */ - void cleanup(); - - /* push data */ - void push( UINT32 data ); - - /* insert data at given index */ - void insert( UINT32 index, UINT32 data ); - - /* return data */ - UINT32* data(); - - /* return data size */ - UINT32 size(); - - /* return capacity */ - UINT32 capacity(); - - /* definition of operator = */ - void operator=( const Array32& ); - - /* print data */ - void print(); - - /* initializer */ - void init( int size ); - - /* data */ - UINT32 *arr; - -}; - -#endif diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/bucket_group.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/bucket_group.hpp deleted file mode 100644 index 7299de062..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/bucket_group.hpp +++ /dev/null @@ -1,76 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#ifndef __OPENCV_BUCKET_GROUP_HPP -#define __OPENCV_BUCKET_GROUP_HPP - -#ifdef _WIN32 -#pragma warning( disable : 4267 ) -#endif - -#include "types.hpp" -#include "array32.hpp" -#include "bitarray.hpp" - -class BucketGroup -{ - - public: - /* constructor */ - BucketGroup(); - - /* destructor */ - ~BucketGroup(); - - /* insert data into the bucket */ - void insert( int subindex, UINT32 data ); - - /* perform a query to the bucket */ - UINT32* query( int subindex, int *size ); - - /* data fields */ - UINT32 empty; - Array32 *group; - -}; - -#endif diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp index 75a57b5e3..63164478a 100644 --- a/modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp +++ b/modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp @@ -42,23 +42,36 @@ #ifndef __OPENCV_DESCRIPTOR_HPP__ #define __OPENCV_DESCRIPTOR_HPP__ -#include "line_structure.hpp" -#include "array32.hpp" -#include "bitarray.hpp" -#include "bitops.hpp" -#include "bucket_group.hpp" -#include "mihasher.hpp" -#include "sparse_hashtable.hpp" -#include "types.hpp" -#include "ed_line_detector.hpp" #include +#include +#include +#include +#include +#include + +#include "opencv2/core/utility.hpp" +#include "opencv2/core/private.hpp" +#include +#include +#include +#include "opencv2/core.hpp" + +/* define data types */ +typedef uint64_t UINT64; +typedef uint32_t UINT32; +typedef uint16_t UINT16; +typedef uint8_t UINT8; + +/* define constants */ +#define UINT64_1 ((UINT64)0x01) +#define UINT32_1 ((UINT32)0x01) namespace cv { namespace line_descriptor { -//CV_EXPORTS bool initModule_line_descriptor(); +CV_EXPORTS bool initModule_line_descriptor(); struct CV_EXPORTS KeyLine { @@ -228,176 +241,906 @@ class CV_EXPORTS BinaryDescriptor : public Algorithm AlgorithmInfo* info() const; private: - /* compute Gaussian pyramids */ - void computeGaussianPyramid( const Mat& image, const int numOctaves ); + /* struct to represent lines extracted from an octave */ + struct OctaveLine + { + unsigned int octaveCount; //the octave which this line is detected + unsigned int lineIDInOctave; //the line ID in that octave image + unsigned int lineIDInScaleLineVec; //the line ID in Scale line vector + float lineLength; //the length of line in original image scale + }; + + // A 2D line (normal equation parameters). + struct SingleLine + { + //note: rho and theta are based on coordinate origin, i.e. the top-left corner of image + double rho; //unit: pixel length + double theta; //unit: rad + double linePointX; // = rho * cos(theta); + double linePointY; // = rho * sin(theta); + //for EndPoints, the coordinate origin is the top-left corner of image. + double startPointX; + double startPointY; + double endPointX; + double endPointY; + //direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis. + double direction; + //mean gradient magnitude + double gradientMagnitude; + //mean gray value of pixels in dark side of line + double darkSideGrayValue; + //mean gray value of pixels in light side of line + double lightSideGrayValue; + //the length of line + double lineLength; + //the width of line; + double width; + //number of pixels + int numOfPixels; + //the decriptor of line + std::vector descriptor; + }; + + // Specifies a vector of lines. + typedef std::vector Lines_list; + + struct OctaveSingleLine + { + /*endPoints, the coordinate origin is the top-left corner of the original image. + *startPointX = sPointInOctaveX * (factor)^octaveCount; */ + float startPointX; + float startPointY; + float endPointX; + float endPointY; + //endPoints, the coordinate origin is the top-left corner of the octave image. + float sPointInOctaveX; + float sPointInOctaveY; + float ePointInOctaveX; + float ePointInOctaveY; + //direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis. + float direction; + //the summation of gradient magnitudes of pixels on lines + float salience; + //the length of line + float lineLength; + //number of pixels + unsigned int numOfPixels; + //the octave which this line is detected + unsigned int octaveCount; + //the decriptor of line + std::vector descriptor; + }; + + struct Pixel + { + unsigned int x; //X coordinate + unsigned int y; //Y coordinate + }; + struct EdgeChains + { + std::vector xCors; //all the x coordinates of edge points + std::vector yCors; //all the y coordinates of edge points + std::vector sId; //the start index of each edge in the coordinate arrays + unsigned int numOfEdges; //the number of edges whose length are larger than minLineLen; numOfEdges < sId.size; + }; + + struct LineChains + { + std::vector xCors; //all the x coordinates of line points + std::vector yCors; //all the y coordinates of line points + std::vector sId; //the start index of each line in the coordinate arrays + unsigned int numOfLines; //the number of lines whose length are larger than minLineLen; numOfLines < sId.size; + }; + + typedef std::list PixelChain; //each edge is a pixel chain + + struct EDLineParam + { + int ksize; + float sigma; + float gradientThreshold; + float anchorThreshold; + int scanIntervals; + int minLineLen; + double lineFitErrThreshold; + }; - /* compute Sobel's derivatives */ - void computeSobel( const Mat& image, const int numOctaves ); + #define RELATIVE_ERROR_FACTOR 100.0 + #define MLN10 2.30258509299404568402 + #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x)) + + /* This class is used to detect lines from input image. + * First, edges are extracted from input image following the method presented in Cihan Topal and + * Cuneyt Akinlar's paper:"Edge Drawing: A Heuristic Approach to Robust Real-Time Edge Detection", 2010. + * Then, lines are extracted from the edge image following the method presented in Cuneyt Akinlar and + * Cihan Topal's paper:"EDLines: A real-time line segment detector with a false detection control", 2011 + * PS: The linking step of edge detection has a little bit difference with the Edge drawing algorithm + * described in the paper. The edge chain doesn't stop when the pixel direction is changed. + */ + class EDLineDetector + { + public: + EDLineDetector(); + EDLineDetector( EDLineParam param ); + ~EDLineDetector(); + + /*extract edges from image + *image: In, gray image; + *edges: Out, store the edges, each edge is a pixel chain + *return -1: error happen + */ + int EdgeDrawing( cv::Mat &image, EdgeChains &edgeChains ); + + /*extract lines from image + *image: In, gray image; + *lines: Out, store the extracted lines, + *return -1: error happen + */ + int EDline( cv::Mat &image, LineChains &lines ); + + /* extract line from image, and store them */ + int EDline( cv::Mat &image ); + + cv::Mat dxImg_; //store the dxImg; + + cv::Mat dyImg_; //store the dyImg; + + cv::Mat gImgWO_; //store the gradient image without threshold; + + LineChains lines_; //store the detected line chains; + + //store the line Equation coefficients, vec3=[w1,w2,w3] for line w1*x + w2*y + w3=0; + std::vector > lineEquations_; + + //store the line endpoints, [x1,y1,x2,y3] + std::vector > lineEndpoints_; + + //store the line direction + std::vector lineDirection_; + + //store the line salience, which is the summation of gradients of pixels on line + std::vector lineSalience_; + + // image sizes + unsigned int imageWidth; + unsigned int imageHeight; + + /*The threshold of line fit error; + *If lineFitErr is large than this threshold, then + *the pixel chain is not accepted as a single line segment.*/ + double lineFitErrThreshold_; + + /*the threshold of pixel gradient magnitude. + *Only those pixel whose gradient magnitude are larger than this threshold will be + *taken as possible edge points. Default value is 36*/ + short gradienThreshold_; + + /*If the pixel's gradient value is bigger than both of its neighbors by a + *certain threshold (ANCHOR_THRESHOLD), the pixel is marked to be an anchor. + *Default value is 8*/ + unsigned char anchorThreshold_; - /* conversion of an LBD descriptor to its binary representation */ - unsigned char binaryConversion( float* f1, float* f2 ); + /*anchor testing can be performed at different scan intervals, i.e., + *every row/column, every second row/column etc. + *Default value is 2*/ + unsigned int scanIntervals_; - /* compute LBD descriptors using EDLine extractor */ - int computeLBD( ScaleLines &keyLines, bool useDetectionData = false ); + int minLineLen_; //minimal acceptable line length - /* gathers lines in groups using EDLine extractor. - Each group contains the same line, detected in different octaves */ - int OctaveKeyLines( cv::Mat& image, ScaleLines &keyLines ); + private: + void InitEDLine_(); - /* the local gaussian coefficient applied to the orthogonal line direction within each band */ - std::vector gaussCoefL_; + /*For an input edge chain, find the best fit line, the default chain length is minLineLen_ + *xCors: In, pointer to the X coordinates of pixel chain; + *yCors: In, pointer to the Y coordinates of pixel chain; + *offsetS:In, start index of this chain in vector; + *lineEquation: Out, [a,b] which are the coefficient of lines y=ax+b(horizontal) or x=ay+b(vertical); + *return: line fit error; -1:error happens; + */ + double LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, std::vector &lineEquation ); + + /*For an input pixel chain, find the best fit line. Only do the update based on new points. + *For A*x=v, Least square estimation of x = Inv(A^T * A) * (A^T * v); + *If some new observations are added, i.e, [A; A'] * x = [v; v'], + *then x' = Inv(A^T * A + (A')^T * A') * (A^T * v + (A')^T * v'); + *xCors: In, pointer to the X coordinates of pixel chain; + *yCors: In, pointer to the Y coordinates of pixel chain; + *offsetS:In, start index of this chain in vector; + *newOffsetS: In, start index of extended part; + *offsetE:In, end index of this chain in vector; + *lineEquation: Out, [a,b] which are the coefficient of lines y=ax+b(horizontal) or x=ay+b(vertical); + *return: line fit error; -1:error happens; + */ + double LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int newOffsetS, unsigned int offsetE, + std::vector &lineEquation ); + + /* Validate line based on the Helmholtz principle, which basically states that + * for a structure to be perceptually meaningful, the expectation of this structure + * by chance must be very low. + */ + bool LineValidation_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int offsetE, std::vector &lineEquation, + float &direction ); + + bool bValidate_; //flag to decide whether line will be validated + + int ksize_; //the size of Gaussian kernel: ksize X ksize, default value is 5. + + float sigma_; //the sigma of Gaussian kernal, default value is 1.0. + + /*For example, there two edges in the image: + *edge1 = [(7,4), (8,5), (9,6),| (10,7)|, (11, 8), (12,9)] and + *edge2 = [(14,9), (15,10), (16,11), (17,12),| (18, 13)|, (19,14)] ; then we store them as following: + *pFirstPartEdgeX_ = [10, 11, 12, 18, 19];//store the first part of each edge[from middle to end] + *pFirstPartEdgeY_ = [7, 8, 9, 13, 14]; + *pFirstPartEdgeS_ = [0,3,5];// the index of start point of first part of each edge + *pSecondPartEdgeX_ = [10, 9, 8, 7, 18, 17, 16, 15, 14];//store the second part of each edge[from middle to front] + *pSecondPartEdgeY_ = [7, 6, 5, 4, 13, 12, 11, 10, 9];//anchor points(10, 7) and (18, 13) are stored again + *pSecondPartEdgeS_ = [0, 4, 9];// the index of start point of second part of each edge + *This type of storage order is because of the order of edge detection process. + *For each edge, start from one anchor point, first go right, then go left or first go down, then go up*/ + + //store the X coordinates of the first part of the pixels for chains + unsigned int *pFirstPartEdgeX_; + + //store the Y coordinates of the first part of the pixels for chains + unsigned int *pFirstPartEdgeY_; + + //store the start index of every edge chain in the first part arrays + unsigned int *pFirstPartEdgeS_; + + //store the X coordinates of the second part of the pixels for chains + unsigned int *pSecondPartEdgeX_; + + //store the Y coordinates of the second part of the pixels for chains + unsigned int *pSecondPartEdgeY_; + + //store the start index of every edge chain in the second part arrays + unsigned int *pSecondPartEdgeS_; + + //store the X coordinates of anchors + unsigned int *pAnchorX_; + + //store the Y coordinates of anchors + unsigned int *pAnchorY_; + + //edges + cv::Mat edgeImage_; + + cv::Mat gImg_; //store the gradient image; + + cv::Mat dirImg_; //store the direction image + + double logNT_; + + cv::Mat_ ATA; //the previous matrix of A^T * A; + + cv::Mat_ ATV; //the previous vector of A^T * V; + + cv::Mat_ fitMatT; //the matrix used in line fit function; + + cv::Mat_ fitVec; //the vector used in line fit function; + + cv::Mat_ tempMatLineFit; //the matrix used in line fit function; + + cv::Mat_ tempVecLineFit; //the vector used in line fit function; + + /** 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 accumulated rounding error in the chain of + computation. 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 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 the natural logarithm of the absolute value of + the gamma function of x using the Lanczos approximation. + See http://www.rskey.org/gamma.htm + The formula used is + @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] + and + q0 = 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 of + the gamma function of x using Windschitl method. + See http://www.rskey.org/gamma.htm + The 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 -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 Tests + The computation is based in the gamma function by the following + relation: + @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 of + the gamma function. + To make the computation faster, not all the sum is computed, part + of the terms are neglected based on a bound to the error obtained + (an error of 10% in the result is accepted). + */ + static double nfa( int n, int k, double p, double logNT ) + { + 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 ) + { + std::cout << "nfa: wrong n, k or p values." << std::endl; + exit( 0 ); + } + /* 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. + But + bincoef(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 / MLN10 - logNT; /* end: use just the first term */ + else + return -logNT; /* begin: the tail is roughly 1 */ +} + +/* compute more terms if needed */ +bin_tail = term; +for ( i = k + 1; i <= n; i++ ) +{ + /* As + term_i = bincoef(n,i) * p^i * (1-p)^(n-i) + and + bincoef(n,i)/bincoef(n,i-1) = n-i+1 / i, + then, + term_i / term_i-1 = (n-i+1)/i * p/(1-p) + and + term_i = term_i-1 * (n-i+1)/i * p/(1-p). + p/(1-p) is computed only once and stored in 'p_term'. + */ + bin_term = (double) ( n - i + 1 ) / (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_ji. + Then, the error on the binomial tail when truncated at + the i term can be bounded by a geometric series of form + term_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 is + given by tolerance*final_result divided by the derivative + of -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; +} +}; - /* the global gaussian coefficient applied to each row within line support region */ - std::vector gaussCoefG_; + // Specifies a vector of lines. +typedef std::vector LinesVec; - /* descriptor parameters */ - Params params; +// each element in ScaleLines is a vector of lines +// which corresponds the same line detected in different octave images. +typedef std::vector ScaleLines; - /* vector of sizes of downsampled and blurred images */ - std::vector images_sizes; +/* compute Gaussian pyramids */ +void computeGaussianPyramid( const Mat& image, const int numOctaves ); - /*For each octave of image, we define an EDLineDetector, because we can get gradient images (dxImg, dyImg, gImg) - *from the EDLineDetector class without extra computation cost. Another reason is that, if we use - *a single EDLineDetector to detect lines in different octave of images, then we need to allocate and release - *memory for gradient images (dxImg, dyImg, gImg) repeatedly for their varying size*/ - std::vector > edLineVec_; +/* compute Sobel's derivatives */ +void computeSobel( const Mat& image, const int numOctaves ); - /* Sobel's derivatives */ - std::vector dxImg_vector, dyImg_vector; +/* conversion of an LBD descriptor to its binary representation */ +unsigned char binaryConversion( float* f1, float* f2 ); - /* Gaussian pyramid */ - std::vector octaveImages; +/* compute LBD descriptors using EDLine extractor */ +int computeLBD( ScaleLines &keyLines, bool useDetectionData = false ); + +/* gathers lines in groups using EDLine extractor. + Each group contains the same line, detected in different octaves */ +int OctaveKeyLines( cv::Mat& image, ScaleLines &keyLines ); + +/* the local gaussian coefficient applied to the orthogonal line direction within each band */ +std::vector gaussCoefL_; + +/* the global gaussian coefficient applied to each row within line support region */ +std::vector gaussCoefG_; + +/* descriptor parameters */ +Params params; + +/* vector of sizes of downsampled and blurred images */ +std::vector images_sizes; + +/*For each octave of image, we define an EDLineDetector, because we can get gradient images (dxImg, dyImg, gImg) + *from the EDLineDetector class without extra computation cost. Another reason is that, if we use + *a single EDLineDetector to detect lines in different octave of images, then we need to allocate and release + *memory for gradient images (dxImg, dyImg, gImg) repeatedly for their varying size*/ +std::vector > edLineVec_; + +/* Sobel's derivatives */ +std::vector dxImg_vector, dyImg_vector; + +/* Gaussian pyramid */ +std::vector octaveImages; }; class CV_EXPORTS LSDDetector : public Algorithm { - public: +public: - /* constructor */ - /*CV_WRAP*/ - LSDDetector() - { - } - ; +/* constructor */ +/*CV_WRAP*/ +LSDDetector() +{ +} +; - /* constructor with smart pointer */ - /*CV_WRAP*/ - static Ptr createLSDDetector(); +/* constructor with smart pointer */ +/*CV_WRAP*/ +static Ptr createLSDDetector(); - /* requires line detection (only one image) */ - /*CV_WRAP*/ - void detect( const Mat& image, CV_OUT std::vector& keypoints, int scale, int numOctaves, const Mat& mask = Mat() ); +/* requires line detection (only one image) */ +/*CV_WRAP*/ +void detect( const Mat& image, CV_OUT std::vector& keypoints, int scale, int numOctaves, const Mat& mask = Mat() ); - /* requires line detection (more than one image) */ - /*CV_WRAP*/ - void detect( const std::vector& images, std::vector >& keylines, int scale, int numOctaves, - const std::vector& masks = std::vector() ) const; +/* requires line detection (more than one image) */ +/*CV_WRAP*/ +void detect( const std::vector& images, std::vector >& keylines, int scale, int numOctaves, +const std::vector& masks = std::vector() ) const; - private: - /* compute Gaussian pyramid of input image */ - void computeGaussianPyramid( const Mat& image, int numOctaves, int scale ); +private: +/* compute Gaussian pyramid of input image */ +void computeGaussianPyramid( const Mat& image, int numOctaves, int scale ); - /* implementation of line detection */ - void detectImpl( const Mat& imageSrc, std::vector& keylines, int numOctaves, int scale, const Mat& mask ) const; +/* implementation of line detection */ +void detectImpl( const Mat& imageSrc, std::vector& keylines, int numOctaves, int scale, const Mat& mask ) const; - /* matrices for Gaussian pyramids */ - std::vector gaussianPyrs; +/* matrices for Gaussian pyramids */ +std::vector gaussianPyrs; - protected: - /* function inherited from Algorithm */ - AlgorithmInfo* info() const; +protected: +/* function inherited from Algorithm */ +AlgorithmInfo* info() const; }; class CV_EXPORTS BinaryDescriptorMatcher : public Algorithm { - public: - /* for every input descriptor, - find the best matching one (for a pair of images) */ - /*CV_WRAP*/ - void match( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector& matches, const Mat& mask = Mat() ) const; +public: +/* for every input descriptor, + find the best matching one (for a pair of images) */ +/*CV_WRAP*/ +void match( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector& matches, const Mat& mask = Mat() ) const; + +/* for every input descriptor, + find the best matching one (from one image to a set) */ +/*CV_WRAP*/ +void match( const Mat& queryDescriptors, std::vector& matches, const std::vector& masks = std::vector() ); + +/* for every input descriptor, + find the best k matching descriptors (for a pair of images) */ +/*CV_WRAP*/ +void knnMatch( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector >& matches, int k, const Mat& mask = Mat(), +bool compactResult = false ) const; + +/* for every input descriptor, + find the best k matching descriptors (from one image to a set) */ +/*CV_WRAP*/ +void knnMatch( const Mat& queryDescriptors, std::vector >& matches, int k, const std::vector& masks = std::vector(), +bool compactResult = false ); + +/* for every input descriptor, find all the ones falling in a + certain matching radius (for a pair of images) */ +/*CV_WRAP*/ +void radiusMatch( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector >& matches, float maxDistance, +const Mat& mask = Mat(), bool compactResult = false ) const; + +/* for every input descriptor, find all the ones falling in a + certain matching radius (from one image to a set) */ +/*CV_WRAP*/ +void radiusMatch( const Mat& queryDescriptors, std::vector >& matches, float maxDistance, const std::vector& masks = +std::vector(), +bool compactResult = false ); + +/* store new descriptors to be inserted in dataset */ +/*CV_WRAP*/ +void add( const std::vector& descriptors ); + +/* store new descriptors into dataset */ +/*CV_WRAP*/ +void train(); + +/* constructor with smart pointer */ +/*CV_WRAP*/ +static Ptr createBinaryDescriptorMatcher(); + +/* clear dataset and internal data */ +/*CV_WRAP*/ +void clear(); + +/* constructor */ +/*CV_WRAP*/ +BinaryDescriptorMatcher(); + +/* destructor */ +~BinaryDescriptorMatcher() +{ +} - /* for every input descriptor, - find the best matching one (from one image to a set) */ - /*CV_WRAP*/ - void match( const Mat& queryDescriptors, std::vector& matches, const std::vector& masks = std::vector() ); +protected: +/* function inherited from Algorithm */ +AlgorithmInfo* info() const; - /* for every input descriptor, - find the best k matching descriptors (for a pair of images) */ - /*CV_WRAP*/ - void knnMatch( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector >& matches, int k, const Mat& mask = Mat(), - bool compactResult = false ) const; +private: +class Array32 +{ - /* for every input descriptor, - find the best k matching descriptors (from one image to a set) */ - /*CV_WRAP*/ - void knnMatch( const Mat& queryDescriptors, std::vector >& matches, int k, const std::vector& masks = std::vector(), - bool compactResult = false ); +public: +/* set ARRAY_RESIZE_FACTOR */ +static void setArrayResizeFactor( double arf ); - /* for every input descriptor, find all the ones falling in a - certain matching radius (for a pair of images) */ - /*CV_WRAP*/ - void radiusMatch( const Mat& queryDescriptors, const Mat& trainDescriptors, std::vector >& matches, float maxDistance, - const Mat& mask = Mat(), bool compactResult = false ) const; +/* constructor */ +Array32(); - /* for every input descriptor, find all the ones falling in a - certain matching radius (from one image to a set) */ - /*CV_WRAP*/ - void radiusMatch( const Mat& queryDescriptors, std::vector >& matches, float maxDistance, const std::vector& masks = - std::vector(), - bool compactResult = false ); +/* destructor */ +~Array32(); - /* store new descriptors to be inserted in dataset */ - /*CV_WRAP*/ - void add( const std::vector& descriptors ); +/* cleaning function used in destructor */ +void cleanup(); - /* store new descriptors into dataset */ - /*CV_WRAP*/ - void train(); +/* push data */ +void push( UINT32 data ); - /* constructor with smart pointer */ - /*CV_WRAP*/ - static Ptr createBinaryDescriptorMatcher(); +/* insert data at given index */ +void insert( UINT32 index, UINT32 data ); - /* clear dataset and internal data */ - /*CV_WRAP*/ - void clear(); +/* return data */ +UINT32* data(); - /* constructor */ - /*CV_WRAP*/ - BinaryDescriptorMatcher(); +/* return data size */ +UINT32 size(); - /* destructor */ - ~BinaryDescriptorMatcher() - { - } +/* return capacity */ +UINT32 capacity(); - protected: - /* function inherited from Algorithm */ - AlgorithmInfo* info() const; +/* definition of operator = */ +void operator=( const Array32& ); - private: - /* retrieve Hamming distances */ - void checkKDistances( UINT32 * numres, int k, std::vector& k_distances, int row, int string_length ) const; +/* print data */ +void print(); + +/* initializer */ +void init( int size ); - /* matrix to store new descriptors */ - Mat descriptorsMat; +/* data */ +UINT32 *arr; - /* map storing where each bunch of descriptors benins in DS */ - std::map indexesMap; +}; + +class BucketGroup +{ - /* internal MiHaser representing dataset */ - Mihasher* dataset; +public: +/* constructor */ +BucketGroup(); - /* index from which next added descriptors' bunch must begin */ - int nextAddedIndex; +/* destructor */ +~BucketGroup(); - /* number of images whose descriptors are stored in DS */ - int numImages; +/* insert data into the bucket */ +void insert( int subindex, UINT32 data ); - /* number of descriptors in dataset */ - int descrInDS; +/* perform a query to the bucket */ +UINT32* query( int subindex, int *size ); + +/* data fields */ +UINT32 empty; +Array32 *group; + +}; + +class SparseHashtable +{ + +private: + +/* Maximum bits per key before folding the table */ +static const int MAX_B; + +/* Bins (each bin is an Array object for duplicates of the same key) */ +BucketGroup *table; + +public: + +/* constructor */ +SparseHashtable(); + +/* destructor */ +~SparseHashtable(); + +/* initializer */ +int init( int _b ); + +/* insert data */ +void insert( UINT64 index, UINT32 data ); + +/* query data */ +UINT32* query( UINT64 index, int* size ); + +/* Bits per index */ +int b; + +/* Number of bins */ +UINT64 size; + +}; + +/* class defining a sequence of bits */ +class bitarray +{ + +public: +/* pointer to bits sequence and sequence's length */ +UINT32 *arr; +UINT32 length; + +/* constructor setting default values */ +bitarray() +{ +arr = NULL; +length = 0; +} + +/* constructor setting sequence's length */ +bitarray( UINT64 _bits ) +{ +init( _bits ); +} + +/* initializer of private fields */ +void init( UINT64 _bits ) +{ +length = (UINT32) ceil( _bits / 32.00 ); +arr = new UINT32[length]; +erase(); +} + +/* destructor */ +~bitarray() +{ +if( arr ) +delete[] arr; +} + +inline void flip( UINT64 index ) +{ +arr[index >> 5] ^= ( (UINT32) 0x01 ) << ( index % 32 ); +} + +inline void set( UINT64 index ) +{ +arr[index >> 5] |= ( (UINT32) 0x01 ) << ( index % 32 ); +} + +inline UINT8 get( UINT64 index ) +{ +return ( arr[index >> 5] & ( ( (UINT32) 0x01 ) << ( index % 32 ) ) ) != 0; +} + +/* reserve menory for an UINT32 */ +inline void erase() +{ +memset( arr, 0, sizeof(UINT32) * length ); +} + +}; + +class Mihasher +{ + +public: +/* Bits per code */ +int B; + +/* B/8 */ +int B_over_8; + +/* Bits per chunk (must be less than 64) */ +int b; + +/* Number of chunks */ +int m; + +/* Number of chunks with b bits (have 1 bit more than others) */ +int mplus; + +/* Maximum hamming search radius (we use B/2 by default) */ +int D; + +/* Maximum hamming search radius per substring */ +int d; + +/* Maximum results to return */ +int K; + +/* Number of codes */ +UINT64 N; + +/* Table of original full-length codes */ +cv::Mat codes; + +/* Counter for eliminating duplicate results (it is not thread safe) */ +bitarray *counter; + +/* Array of m hashtables */ +SparseHashtable *H; + +/* Volume of a b-bit Hamming ball with radius s (for s = 0 to d) */ +UINT32 *xornum; + +/* Used within generation of binary codes at a certain Hamming distance */ +int power[100]; + +/* constructor */ +Mihasher(); + +/* desctructor */ +~Mihasher(); + +/* constructor 2 */ +Mihasher( int B, int m ); + +/* K setter */ +void setK( int K ); + +/* populate tables */ +void populate( cv::Mat & codes, UINT32 N, int dim1codes ); + +/* execute a batch query */ +void batchquery( UINT32 * results, UINT32 *numres/*, qstat *stats*/, const cv::Mat & q, UINT32 numq, int dim1queries ); + +private: + +/* execute a single query */ +void query( UINT32 * results, UINT32* numres/*, qstat *stats*/, UINT8 *q, UINT64 * chunks, UINT32 * res ); +}; + +/* retrieve Hamming distances */ +void checkKDistances( UINT32 * numres, int k, std::vector& k_distances, int row, int string_length ) const; + +/* matrix to store new descriptors */ +Mat descriptorsMat; + +/* map storing where each bunch of descriptors benins in DS */ +std::map indexesMap; + +/* internal MiHaser representing dataset */ +Mihasher* dataset; + +/* index from which next added descriptors' bunch must begin */ +int nextAddedIndex; + +/* number of images whose descriptors are stored in DS */ +int numImages; + +/* number of descriptors in dataset */ +int descrInDS; }; @@ -408,17 +1151,17 @@ class CV_EXPORTS BinaryDescriptorMatcher : public Algorithm /* struct for drawing options */ struct CV_EXPORTS DrawLinesMatchesFlags { - enum - { - DEFAULT = 0, // Output image matrix will be created (Mat::create), - // i.e. existing memory of output image may be reused. - // Two source images, matches, and single keylines - // will be drawn. - DRAW_OVER_OUTIMG = 1, // Output image matrix will not be - // created (using Mat::create). Matches will be drawn - // on existing content of output image. - NOT_DRAW_SINGLE_LINES = 2 // Single keylines will not be drawn. - }; +enum +{ +DEFAULT = 0, // Output image matrix will be created (Mat::create), + // i.e. existing memory of output image may be reused. + // Two source images, matches, and single keylines + // will be drawn. +DRAW_OVER_OUTIMG = 1,// Output image matrix will not be +// created (using Mat::create). Matches will be drawn +// on existing content of output image. +NOT_DRAW_SINGLE_LINES = 2// Single keylines will not be drawn. +}; }; /* draw matches between two images */ diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/ed_line_detector.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/ed_line_detector.hpp deleted file mode 100644 index 7505e4b13..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/ed_line_detector.hpp +++ /dev/null @@ -1,476 +0,0 @@ -/*IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - - By downloading, copying, installing or using the software you agree to this license. - If you do not agree to this license, do not download, install, - copy or use the software. - - - License Agreement - For Open Source Computer Vision Library - - Copyright (C) 2011-2012, Lilian Zhang, all rights reserved. - Copyright (C) 2013, Manuele Tamburrano, Stefano Fabri, all rights reserved. - Copyright (C) 2014, Biagio Montesano, all rights reserved. - Third party copyrights are property of their respective owners. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * The name of the copyright holders may not be used to endorse or promote products - derived from this software without specific prior written permission. - - This software is provided by the copyright holders and contributors "as is" and - any express or implied warranties, including, but not limited to, the implied - warranties of merchantability and fitness for a particular purpose are disclaimed. - In no event shall the Intel Corporation or contributors be liable for any direct, - indirect, incidental, special, exemplary, or consequential damages - (including, but not limited to, procurement of substitute goods or services; - loss of use, data, or profits; or business interruption) however caused - and on any theory of liability, whether in contract, strict liability, - or tort (including negligence or otherwise) arising in any way out of - the use of this software, even if advised of the possibility of such damage. - */ - -#ifndef __OPENCV_ED_LINE_DETECTOR_HH_ -#define __OPENCV_ED_LINE_DETECTOR_HH_ - -#ifdef _WIN32 -#pragma warning( disable : 4267 ) -#endif - -#include -#include -#include -#include -#include -#include -#include -#include - -struct Pixel -{ - unsigned int x; //X coordinate - unsigned int y; //Y coordinate -}; -struct EdgeChains -{ - std::vector xCors; //all the x coordinates of edge points - std::vector yCors; //all the y coordinates of edge points - std::vector sId; //the start index of each edge in the coordinate arrays - unsigned int numOfEdges; //the number of edges whose length are larger than minLineLen; numOfEdges < sId.size; -}; - -struct LineChains -{ - std::vector xCors; //all the x coordinates of line points - std::vector yCors; //all the y coordinates of line points - std::vector sId; //the start index of each line in the coordinate arrays - unsigned int numOfLines; //the number of lines whose length are larger than minLineLen; numOfLines < sId.size; -}; - -typedef std::list PixelChain; //each edge is a pixel chain - -struct EDLineParam -{ - int ksize; - float sigma; - float gradientThreshold; - float anchorThreshold; - int scanIntervals; - int minLineLen; - double lineFitErrThreshold; -}; - -#define RELATIVE_ERROR_FACTOR 100.0 -#define MLN10 2.30258509299404568402 -#define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x)) - -/* This class is used to detect lines from input image. - * First, edges are extracted from input image following the method presented in Cihan Topal and - * Cuneyt Akinlar's paper:"Edge Drawing: A Heuristic Approach to Robust Real-Time Edge Detection", 2010. - * Then, lines are extracted from the edge image following the method presented in Cuneyt Akinlar and - * Cihan Topal's paper:"EDLines: A real-time line segment detector with a false detection control", 2011 - * PS: The linking step of edge detection has a little bit difference with the Edge drawing algorithm - * described in the paper. The edge chain doesn't stop when the pixel direction is changed. - */ -class EDLineDetector -{ - public: - EDLineDetector(); - EDLineDetector( EDLineParam param ); - ~EDLineDetector(); - - /*extract edges from image - *image: In, gray image; - *edges: Out, store the edges, each edge is a pixel chain - *return -1: error happen - */ - int EdgeDrawing( cv::Mat &image, EdgeChains &edgeChains ); - - /*extract lines from image - *image: In, gray image; - *lines: Out, store the extracted lines, - *return -1: error happen - */ - int EDline( cv::Mat &image, LineChains &lines ); - - /* extract line from image, and store them */ - int EDline( cv::Mat &image ); - - cv::Mat dxImg_; //store the dxImg; - - cv::Mat dyImg_; //store the dyImg; - - cv::Mat gImgWO_; //store the gradient image without threshold; - - LineChains lines_; //store the detected line chains; - - //store the line Equation coefficients, vec3=[w1,w2,w3] for line w1*x + w2*y + w3=0; - std::vector > lineEquations_; - - //store the line endpoints, [x1,y1,x2,y3] - std::vector > lineEndpoints_; - - //store the line direction - std::vector lineDirection_; - - //store the line salience, which is the summation of gradients of pixels on line - std::vector lineSalience_; - - // image sizes - unsigned int imageWidth; - unsigned int imageHeight; - - /*The threshold of line fit error; - *If lineFitErr is large than this threshold, then - *the pixel chain is not accepted as a single line segment.*/ - double lineFitErrThreshold_; - - /*the threshold of pixel gradient magnitude. - *Only those pixel whose gradient magnitude are larger than this threshold will be - *taken as possible edge points. Default value is 36*/ - short gradienThreshold_; - - /*If the pixel's gradient value is bigger than both of its neighbors by a - *certain threshold (ANCHOR_THRESHOLD), the pixel is marked to be an anchor. - *Default value is 8*/ - unsigned char anchorThreshold_; - - /*anchor testing can be performed at different scan intervals, i.e., - *every row/column, every second row/column etc. - *Default value is 2*/ - unsigned int scanIntervals_; - - int minLineLen_; //minimal acceptable line length - - private: - void InitEDLine_(); - - /*For an input edge chain, find the best fit line, the default chain length is minLineLen_ - *xCors: In, pointer to the X coordinates of pixel chain; - *yCors: In, pointer to the Y coordinates of pixel chain; - *offsetS:In, start index of this chain in vector; - *lineEquation: Out, [a,b] which are the coefficient of lines y=ax+b(horizontal) or x=ay+b(vertical); - *return: line fit error; -1:error happens; - */ - double LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, std::vector &lineEquation ); - - /*For an input pixel chain, find the best fit line. Only do the update based on new points. - *For A*x=v, Least square estimation of x = Inv(A^T * A) * (A^T * v); - *If some new observations are added, i.e, [A; A'] * x = [v; v'], - *then x' = Inv(A^T * A + (A')^T * A') * (A^T * v + (A')^T * v'); - *xCors: In, pointer to the X coordinates of pixel chain; - *yCors: In, pointer to the Y coordinates of pixel chain; - *offsetS:In, start index of this chain in vector; - *newOffsetS: In, start index of extended part; - *offsetE:In, end index of this chain in vector; - *lineEquation: Out, [a,b] which are the coefficient of lines y=ax+b(horizontal) or x=ay+b(vertical); - *return: line fit error; -1:error happens; - */ - double LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int newOffsetS, unsigned int offsetE, - std::vector &lineEquation ); - - /* Validate line based on the Helmholtz principle, which basically states that - * for a structure to be perceptually meaningful, the expectation of this structure - * by chance must be very low. - */ - bool LineValidation_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int offsetE, std::vector &lineEquation, - float &direction ); - - bool bValidate_; //flag to decide whether line will be validated - - int ksize_; //the size of Gaussian kernel: ksize X ksize, default value is 5. - - float sigma_; //the sigma of Gaussian kernal, default value is 1.0. - - /*For example, there two edges in the image: - *edge1 = [(7,4), (8,5), (9,6),| (10,7)|, (11, 8), (12,9)] and - *edge2 = [(14,9), (15,10), (16,11), (17,12),| (18, 13)|, (19,14)] ; then we store them as following: - *pFirstPartEdgeX_ = [10, 11, 12, 18, 19];//store the first part of each edge[from middle to end] - *pFirstPartEdgeY_ = [7, 8, 9, 13, 14]; - *pFirstPartEdgeS_ = [0,3,5];// the index of start point of first part of each edge - *pSecondPartEdgeX_ = [10, 9, 8, 7, 18, 17, 16, 15, 14];//store the second part of each edge[from middle to front] - *pSecondPartEdgeY_ = [7, 6, 5, 4, 13, 12, 11, 10, 9];//anchor points(10, 7) and (18, 13) are stored again - *pSecondPartEdgeS_ = [0, 4, 9];// the index of start point of second part of each edge - *This type of storage order is because of the order of edge detection process. - *For each edge, start from one anchor point, first go right, then go left or first go down, then go up*/ - - //store the X coordinates of the first part of the pixels for chains - unsigned int *pFirstPartEdgeX_; - - //store the Y coordinates of the first part of the pixels for chains - unsigned int *pFirstPartEdgeY_; - - //store the start index of every edge chain in the first part arrays - unsigned int *pFirstPartEdgeS_; - - //store the X coordinates of the second part of the pixels for chains - unsigned int *pSecondPartEdgeX_; - - //store the Y coordinates of the second part of the pixels for chains - unsigned int *pSecondPartEdgeY_; - - //store the start index of every edge chain in the second part arrays - unsigned int *pSecondPartEdgeS_; - - //store the X coordinates of anchors - unsigned int *pAnchorX_; - - //store the Y coordinates of anchors - unsigned int *pAnchorY_; - - //edges - cv::Mat edgeImage_; - - cv::Mat gImg_; //store the gradient image; - - cv::Mat dirImg_; //store the direction image - - double logNT_; - - cv::Mat_ ATA; //the previous matrix of A^T * A; - - cv::Mat_ ATV; //the previous vector of A^T * V; - - cv::Mat_ fitMatT; //the matrix used in line fit function; - - cv::Mat_ fitVec; //the vector used in line fit function; - - cv::Mat_ tempMatLineFit; //the matrix used in line fit function; - - cv::Mat_ tempVecLineFit; //the vector used in line fit function; - - /** 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 accumulated rounding error in the chain of - computation. 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 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 the natural logarithm of the absolute value of - the gamma function of x using the Lanczos approximation. - See http://www.rskey.org/gamma.htm - The formula used is - @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] - and - q0 = 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 of - the gamma function of x using Windschitl method. - See http://www.rskey.org/gamma.htm - The 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 -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 Tests - The computation is based in the gamma function by the following - relation: - @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 of - the gamma function. - To make the computation faster, not all the sum is computed, part - of the terms are neglected based on a bound to the error obtained - (an error of 10% in the result is accepted). - */ - static double nfa( int n, int k, double p, double logNT ) - { - 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 ) - { - std::cout << "nfa: wrong n, k or p values." << std::endl; - exit( 0 ); - } - /* 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. - But - bincoef(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 / MLN10 - logNT; /* end: use just the first term */ - else - return -logNT; /* begin: the tail is roughly 1 */ - } - - /* compute more terms if needed */ - bin_tail = term; - for ( i = k + 1; i <= n; i++ ) - { - /* As - term_i = bincoef(n,i) * p^i * (1-p)^(n-i) - and - bincoef(n,i)/bincoef(n,i-1) = n-i+1 / i, - then, - term_i / term_i-1 = (n-i+1)/i * p/(1-p) - and - term_i = term_i-1 * (n-i+1)/i * p/(1-p). - p/(1-p) is computed only once and stored in 'p_term'. - */ - bin_term = (double) ( n - i + 1 ) / (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_ji. - Then, the error on the binomial tail when truncated at - the i term can be bounded by a geometric series of form - term_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 is - given by tolerance*final_result divided by the derivative - of -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; - } -}; - -#endif /* EDLINEDETECTOR_HH_ */ diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/line_structure.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/line_structure.hpp deleted file mode 100644 index 78f4db3db..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/line_structure.hpp +++ /dev/null @@ -1,123 +0,0 @@ -/*IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - - By downloading, copying, installing or using the software you agree to this license. - If you do not agree to this license, do not download, install, - copy or use the software. - - - License Agreement - For Open Source Computer Vision Library - - Copyright (C) 2011-2012, Lilian Zhang, all rights reserved. - Copyright (C) 2013, Manuele Tamburrano, Stefano Fabri, all rights reserved. - Copyright (C) 2014, Biagio Montesano, all rights reserved. - Third party copyrights are property of their respective owners. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * The name of the copyright holders may not be used to endorse or promote products - derived from this software without specific prior written permission. - - This software is provided by the copyright holders and contributors "as is" and - any express or implied warranties, including, but not limited to, the implied - warranties of merchantability and fitness for a particular purpose are disclaimed. - In no event shall the Intel Corporation or contributors be liable for any direct, - indirect, incidental, special, exemplary, or consequential damages - (including, but not limited to, procurement of substitute goods or services; - loss of use, data, or profits; or business interruption) however caused - and on any theory of liability, whether in contract, strict liability, - or tort (including negligence or otherwise) arising in any way out of - the use of this software, even if advised of the possibility of such damage. - */ - -#ifndef __OPENCV_LINE_STRUCTURE_HH_ -#define __OPENCV_LINE_STRUCTURE_HH_ - -#include - -/* struct to represent lines extracted from an octave */ -struct OctaveLine -{ - unsigned int octaveCount; //the octave which this line is detected - unsigned int lineIDInOctave; //the line ID in that octave image - unsigned int lineIDInScaleLineVec; //the line ID in Scale line vector - float lineLength; //the length of line in original image scale -}; - -// A 2D line (normal equation parameters). -struct SingleLine -{ - //note: rho and theta are based on coordinate origin, i.e. the top-left corner of image - double rho; //unit: pixel length - double theta; //unit: rad - double linePointX; // = rho * cos(theta); - double linePointY; // = rho * sin(theta); - //for EndPoints, the coordinate origin is the top-left corner of image. - double startPointX; - double startPointY; - double endPointX; - double endPointY; - //direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis. - double direction; - //mean gradient magnitude - double gradientMagnitude; - //mean gray value of pixels in dark side of line - double darkSideGrayValue; - //mean gray value of pixels in light side of line - double lightSideGrayValue; - //the length of line - double lineLength; - //the width of line; - double width; - //number of pixels - int numOfPixels; - //the decriptor of line - std::vector descriptor; -}; - -// Specifies a vector of lines. -typedef std::vector Lines_list; - -struct OctaveSingleLine -{ - /*endPoints, the coordinate origin is the top-left corner of the original image. - *startPointX = sPointInOctaveX * (factor)^octaveCount; */ - float startPointX; - float startPointY; - float endPointX; - float endPointY; - //endPoints, the coordinate origin is the top-left corner of the octave image. - float sPointInOctaveX; - float sPointInOctaveY; - float ePointInOctaveX; - float ePointInOctaveY; - //direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis. - float direction; - //the summation of gradient magnitudes of pixels on lines - float salience; - //the length of line - float lineLength; - //number of pixels - unsigned int numOfPixels; - //the octave which this line is detected - unsigned int octaveCount; - //the decriptor of line - std::vector descriptor; -}; - -// Specifies a vector of lines. -typedef std::vector LinesVec; - -// each element in ScaleLines is a vector of lines -// which corresponds the same line detected in different octave images. -typedef std::vector ScaleLines; - -#endif diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/mihasher.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/mihasher.hpp deleted file mode 100644 index aa64565f8..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/mihasher.hpp +++ /dev/null @@ -1,130 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#ifndef __OPENCV_MIHASHER_HPP -#define __OPENCV_MIHASHER_HPP - -#ifdef _WIN32 -#pragma warning( disable : 4267 ) -#endif - -#include "types.hpp" -#include "bitops.hpp" -#include "sparse_hashtable.hpp" -#include "bitarray.hpp" - -#include -#include - -class Mihasher -{ - private: - - /* Bits per code */ - int B; - - /* B/8 */ - int B_over_8; - - /* Bits per chunk (must be less than 64) */ - int b; - - /* Number of chunks */ - int m; - - /* Number of chunks with b bits (have 1 bit more than others) */ - int mplus; - - /* Maximum hamming search radius (we use B/2 by default) */ - int D; - - /* Maximum hamming search radius per substring */ - int d; - - /* Maximum results to return */ - int K; - - /* Number of codes */ - UINT64 N; - - /* Table of original full-length codes */ - cv::Mat codes; - - /* Counter for eliminating duplicate results (it is not thread safe) */ - bitarray *counter; - - /* Array of m hashtables */ - SparseHashtable *H; - - /* Volume of a b-bit Hamming ball with radius s (for s = 0 to d) */ - UINT32 *xornum; - - /* Used within generation of binary codes at a certain Hamming distance */ - int power[100]; - - public: - - /* constructor */ - Mihasher(); - - /* desctructor */ - ~Mihasher(); - - /* constructor 2 */ - Mihasher( int B, int m ); - - /* K setter */ - void setK( int K ); - - /* populate tables */ - void populate( cv::Mat & codes, UINT32 N, int dim1codes ); - - /* execute a batch query */ - void batchquery( UINT32 * results, UINT32 *numres/*, qstat *stats*/, const cv::Mat & q, UINT32 numq, int dim1queries ); - - private: - - /* execute a single query */ - void query( UINT32 * results, UINT32* numres/*, qstat *stats*/, UINT8 *q, UINT64 * chunks, UINT32 * res ); -}; - -#endif diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/sparse_hashtable.hpp b/modules/line_descriptor/include/opencv2/line_descriptor/sparse_hashtable.hpp deleted file mode 100644 index 98586a3a9..000000000 --- a/modules/line_descriptor/include/opencv2/line_descriptor/sparse_hashtable.hpp +++ /dev/null @@ -1,89 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#ifndef __OPENCV_SPARSE_HASHTABLE_HPP -#define __OPENCV_SPARSE_HASHTABLE_HPP - -#ifdef _WIN32 -#pragma warning( disable : 4267 ) -#endif - -#include "types.hpp" -#include "bucket_group.hpp" - -class SparseHashtable -{ - - private: - - /* Maximum bits per key before folding the table */ - static const int MAX_B; - - /* Bins (each bin is an Array object for duplicates of the same key) */ - BucketGroup *table; - - public: - - /* constructor */ - SparseHashtable(); - - /* destructor */ - ~SparseHashtable(); - - /* initializer */ - int init( int _b ); - - /* insert data */ - void insert( UINT64 index, UINT32 data ); - - /* query data */ - UINT32* query( UINT64 index, int* size ); - - /* Bits per index */ - int b; - - /* Number of bins */ - UINT64 size; - -}; - -#endif diff --git a/modules/line_descriptor/samples/radius_matching.cpp b/modules/line_descriptor/samples/radius_matching.cpp index 91da9f5ec..fe71457e9 100644 --- a/modules/line_descriptor/samples/radius_matching.cpp +++ b/modules/line_descriptor/samples/radius_matching.cpp @@ -62,22 +62,22 @@ static const char* keys = static void help() { std::cout << "\nThis example shows the functionalities of radius matching " << "Please, run this sample using a command in the form\n" - << "./example_line_descriptor_radius_matching /" << std::endl; + << "./example_line_descriptor_radius_matching /" << std::endl; } int main( int argc, char** argv ) { /* get parameters from comand line */ CommandLineParser parser( argc, argv, keys ); - String pathToImages = parser.get( 0 ); + String pathToImages = parser.get < String > ( 0 ); /* create structures for hosting KeyLines and descriptors */ int num_elements = sizeof ( images ) / sizeof ( images[0] ); - std::vector descriptorsMat; - std::vector > linesMat; + std::vector < Mat > descriptorsMat; + std::vector < std::vector > linesMat; /*create a pointer to a BinaryDescriptor object */ - Ptr bd = BinaryDescriptor::createBinaryDescriptor(); + Ptr < BinaryDescriptor > bd = BinaryDescriptor::createBinaryDescriptor(); /* compute lines and descriptors */ for ( int i = 0; i < num_elements; i++ ) @@ -97,7 +97,7 @@ int main( int argc, char** argv ) } /* compute lines and descriptors */ - std::vector lines; + std::vector < KeyLine > lines; Mat computedDescr; bd->detect( loadedImage, lines ); bd->compute( loadedImage, lines, computedDescr ); @@ -121,19 +121,19 @@ int main( int argc, char** argv ) std::cout << "It has been generated a matrix of " << queries.rows << " descriptors" << std::endl; /* create a BinaryDescriptorMatcher object */ - Ptr bdm = BinaryDescriptorMatcher::createBinaryDescriptorMatcher(); + Ptr < BinaryDescriptorMatcher > bdm = BinaryDescriptorMatcher::createBinaryDescriptorMatcher(); /* populate matcher */ bdm->add( descriptorsMat ); /* compute matches */ - std::vector > matches; + std::vector < std::vector > matches; bdm->radiusMatch( queries, matches, 30 ); std::cout << "size matches sample " << matches.size() << std::endl; - for ( int i = 0; i < matches.size(); i++ ) + for ( int i = 0; i < (int) matches.size(); i++ ) { - for ( int j = 0; j < matches[i].size(); j++ ) + for ( int j = 0; j < (int) matches[i].size(); j++ ) { std::cout << "match: " << matches[i][j].queryIdx << " " << matches[i][j].trainIdx << " " << matches[i][j].distance << std::endl; } diff --git a/modules/line_descriptor/src/array32.cpp b/modules/line_descriptor/src/array32.cpp deleted file mode 100644 index cc42dd601..000000000 --- a/modules/line_descriptor/src/array32.cpp +++ /dev/null @@ -1,189 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#include "precomp.hpp" - -/* no need for the static keyword in the definition */ -double Array32::ARRAY_RESIZE_FACTOR = 1.1; // minimum is 1.0 -double Array32::ARRAY_RESIZE_ADD_FACTOR = 4; // minimum is 1 - -/* set ARRAY_RESIZE_FACTOR */ -void Array32::setArrayResizeFactor( double arf ) -{ - ARRAY_RESIZE_FACTOR = arf; -} - -/* constructor */ -Array32::Array32() -{ - arr = NULL; -} - -/* definition of operator = - Array32& Array32::operator = (const Array32 &rhs) */ -void Array32::operator =( const Array32 &rhs ) -{ - if( &rhs != this ) - this->arr = rhs.arr; -} - -/* destructor */ -Array32::~Array32() -{ - cleanup(); -} - -/* cleaning function used in destructor */ -void Array32::cleanup() -{ - free( arr ); -} - -/* push data */ -void Array32::push( UINT32 Data ) -{ - if( arr ) - { - if( arr[0] == arr[1] ) - { - arr[1] = (UINT32) std::max( ceil( arr[1] * ARRAY_RESIZE_FACTOR ), arr[1] + ARRAY_RESIZE_ADD_FACTOR ); - UINT32* new_Data = static_cast( realloc( arr, sizeof(UINT32) * ( 2 + arr[1] ) ) ); - if( new_Data == NULL ) - { - /* could not realloc, but orig still valid */ - std::cout << "ALERT!!!! Not enough memory, operation aborted!" << std::endl; - exit( 0 ); - } - else - { - arr = new_Data; - } - - } - - arr[2 + arr[0]] = Data; - arr[0]++; - - } - - else - { - arr = (UINT32*) malloc( (size_t) ( 2 + ARRAY_RESIZE_ADD_FACTOR ) * sizeof(UINT32) ); - arr[0] = 1; - arr[1] = 1; - arr[2] = Data; - } -} - -/* insert data at given index */ -void Array32::insert( UINT32 index, UINT32 Data ) -{ - if( arr ) - { - if( arr[0] == arr[1] ) - { - arr[1] = (UINT32) ceil( arr[0] * 1.1 ); - UINT32* new_data = static_cast( realloc( arr, sizeof(UINT32) * ( 2 + arr[1] ) ) ); - if( new_data == NULL ) - { - // could not realloc, but orig still valid - std::cout << "ALERT!!!! Not enough memory, operation aborted!" << std::endl; - exit( 0 ); - } - else - { - arr = new_data; - } - } - - memmove( arr + ( 2 + index ) + 1, arr + ( 2 + index ), ( arr[0] - index ) * sizeof ( *arr ) ); - - arr[2 + index] = Data; - arr[0]++; - } - - else - { - arr = (UINT32*) malloc( 3 * sizeof(UINT32) ); - arr[0] = 1; - arr[1] = 1; - arr[2] = Data; - } -} - -/* return data */ -UINT32* Array32::data() -{ - return arr ? arr + 2 : NULL; -} - -/* return data size */ -UINT32 Array32::size() -{ - return arr ? arr[0] : 0; -} - -/* return capacity */ -UINT32 Array32::capacity() -{ - return arr ? arr[1] : 0; -} - -/* print data */ -void Array32::print() -{ - for ( int i = 0; i < (int) size(); i++ ) - printf( "%d, ", arr[i + 2] ); - - printf( "\n" ); -} - -/* initializer */ -void Array32::init( int Size ) -{ - if( arr == NULL ) - { - arr = (UINT32*) malloc( ( 2 + Size ) * sizeof(UINT32) ); - arr[0] = 0; - arr[1] = Size; - } -} diff --git a/modules/line_descriptor/src/binary_descriptor.cpp b/modules/line_descriptor/src/binary_descriptor.cpp index a133a6eee..5ea1313f2 100644 --- a/modules/line_descriptor/src/binary_descriptor.cpp +++ b/modules/line_descriptor/src/binary_descriptor.cpp @@ -42,6 +42,14 @@ #include "precomp.hpp" #define NUM_OF_BANDS 9 +#define Horizontal 255//if |dx|<|dy|; +#define Vertical 0//if |dy|<=|dx|; +#define UpDir 1 +#define RightDir 2 +#define DownDir 3 +#define LeftDir 4 +#define TryTime 6 +#define SkipEdgePoint 2 //using namespace cv; namespace cv @@ -119,7 +127,7 @@ void BinaryDescriptor::setWidthOfBand( int width ) images_sizes.resize( params.numOfOctave_ ); for ( int i = 0; i < params.numOfOctave_; i++ ) - edLineVec_[i] = Ptr( new EDLineDetector() ); + edLineVec_[i] = Ptr < EDLineDetector > ( new EDLineDetector() ); /* prepare a vector to host local weights F_l*/ gaussCoefL_.resize( params.widthOfBand_ * 3 ); @@ -184,12 +192,12 @@ void BinaryDescriptor::Params::write( cv::FileStorage& fs ) const Ptr BinaryDescriptor::createBinaryDescriptor() { - return Ptr( new BinaryDescriptor() ); + return Ptr < BinaryDescriptor > ( new BinaryDescriptor() ); } Ptr BinaryDescriptor::createBinaryDescriptor( Params parameters ) { - return Ptr( new BinaryDescriptor( parameters ) ); + return Ptr < BinaryDescriptor > ( new BinaryDescriptor( parameters ) ); } /* construct a BinaryDescrptor object and compute external private parameters */ @@ -201,7 +209,7 @@ BinaryDescriptor::BinaryDescriptor( const BinaryDescriptor::Params ¶meters ) images_sizes.resize( params.numOfOctave_ ); for ( int i = 0; i < params.numOfOctave_; i++ ) - edLineVec_[i] = Ptr( new EDLineDetector() ); + edLineVec_[i] = Ptr < EDLineDetector > ( new EDLineDetector() ); /* prepare a vector to host local weights F_l*/ gaussCoefL_.resize( params.widthOfBand_ * 3 ); @@ -239,7 +247,7 @@ BinaryDescriptor::BinaryDescriptor( const BinaryDescriptor::Params ¶meters ) /* definition of operator () */ void BinaryDescriptor::operator()( InputArray image, InputArray mask, CV_OUT std::vector& keylines, OutputArray descriptors, - bool useProvidedKeyLines, bool returnFloatDescr ) const + bool useProvidedKeyLines, bool returnFloatDescr ) const { /* create some matrix objects */ @@ -258,7 +266,7 @@ void BinaryDescriptor::operator()( InputArray image, InputArray mask, CV_OUT std bn->edLineVec_.resize( params.numOfOctave_ ); for ( int i = 0; i < params.numOfOctave_; i++ ) - bn->edLineVec_[i] = Ptr( new EDLineDetector() ); + bn->edLineVec_[i] = Ptr( new EDLineDetector() ); detectImpl( imageMat, keylines, maskMat ); @@ -270,10 +278,10 @@ void BinaryDescriptor::operator()( InputArray image, InputArray mask, CV_OUT std //descrMat = descriptors.getMat(); /* compute descriptors */ if( !useProvidedKeyLines ) - computeImpl( imageMat, keylines, descrMat, returnFloatDescr, true ); + computeImpl( imageMat, keylines, descrMat, returnFloatDescr, true ); else - computeImpl( imageMat, keylines, descrMat, returnFloatDescr, false ); + computeImpl( imageMat, keylines, descrMat, returnFloatDescr, false ); descrMat.copyTo( descriptors ); } @@ -400,10 +408,10 @@ void BinaryDescriptor::detect( const Mat& image, CV_OUT std::vector& ke } if( mask.data != NULL && ( mask.size() != image.size() || mask.type() != CV_8UC1 ) ) - throw std::runtime_error( "Mask error while detecting lines: please check its dimensions and that data type is CV_8UC1" ); + throw std::runtime_error( "Mask error while detecting lines: please check its dimensions and that data type is CV_8UC1" ); else - detectImpl( image, keylines, mask ); + detectImpl( image, keylines, mask ); } /* requires line detection (more than one image) */ @@ -491,7 +499,7 @@ void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector& ke for ( size_t keyCounter = 0; keyCounter < keylines.size(); keyCounter++ ) { KeyLine kl = keylines[keyCounter]; - if( mask.at( (int) kl.startPointY, (int) kl.startPointX ) == 0 && mask.at( (int) kl.endPointY, (int) kl.endPointX ) == 0 ) + if( mask.at < uchar > ( (int) kl.startPointY, (int) kl.startPointX ) == 0 && mask.at < uchar > ( (int) kl.endPointY, (int) kl.endPointX ) == 0 ) keylines.erase( keylines.begin() + keyCounter ); } } @@ -500,7 +508,7 @@ void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector& ke /* requires descriptors computation (only one image) */ void BinaryDescriptor::compute( const Mat& image, CV_OUT CV_IN_OUT std::vector& keylines, CV_OUT Mat& descriptors, - bool returnFloatDescr ) const + bool returnFloatDescr ) const { computeImpl( image, keylines, descriptors, returnFloatDescr, false ); } @@ -703,7 +711,7 @@ int BinaryDescriptor::OctaveKeyLines( cv::Mat& image, ScaleLines &keyLines ) } /* end of loop over number of octaves */ /* prepare a vector to store octave information associated to extracted lines */ - std::vector octaveLines( numOfFinalLine ); + std::vector < OctaveLine > octaveLines( numOfFinalLine ); /* set lines' counter to 0 for reuse */ numOfFinalLine = 0; @@ -1346,5 +1354,1385 @@ int BinaryDescriptor::computeLBD( ScaleLines &keyLines, bool useDetectionData ) return 1; } + +BinaryDescriptor::EDLineDetector::EDLineDetector() +{ + //set parameters for line segment detection + ksize_ = 15; //15 + sigma_ = 30.0; //30 + gradienThreshold_ = 80; // ***** ORIGINAL WAS 25 + anchorThreshold_ = 8; //8 + scanIntervals_ = 2; //2 + minLineLen_ = 15; //15 + lineFitErrThreshold_ = 1.6; //1.4 + InitEDLine_(); +} +BinaryDescriptor::EDLineDetector::EDLineDetector( EDLineParam param ) +{ + //set parameters for line segment detection + ksize_ = param.ksize; + sigma_ = param.sigma; + gradienThreshold_ = (short) param.gradientThreshold; + anchorThreshold_ = (unsigned char) param.anchorThreshold; + scanIntervals_ = param.scanIntervals; + minLineLen_ = param.minLineLen; + lineFitErrThreshold_ = param.lineFitErrThreshold; + InitEDLine_(); +} +void BinaryDescriptor::EDLineDetector::InitEDLine_() +{ + bValidate_ = true; + ATA = cv::Mat_( 2, 2 ); + ATV = cv::Mat_( 1, 2 ); + tempMatLineFit = cv::Mat_( 2, 2 ); + tempVecLineFit = cv::Mat_( 1, 2 ); + fitMatT = cv::Mat_( 2, minLineLen_ ); + fitVec = cv::Mat_( 1, minLineLen_ ); + for ( int i = 0; i < minLineLen_; i++ ) + { + fitMatT[1][i] = 1; + } + dxImg_.create( 1, 1, CV_16SC1 ); + dyImg_.create( 1, 1, CV_16SC1 ); + gImgWO_.create( 1, 1, CV_8SC1 ); + pFirstPartEdgeX_ = NULL; + pFirstPartEdgeY_ = NULL; + pFirstPartEdgeS_ = NULL; + pSecondPartEdgeX_ = NULL; + pSecondPartEdgeY_ = NULL; + pSecondPartEdgeS_ = NULL; + pAnchorX_ = NULL; + pAnchorY_ = NULL; +} + +BinaryDescriptor::EDLineDetector::~EDLineDetector() +{ + if( pFirstPartEdgeX_ != NULL ) + { + delete[] pFirstPartEdgeX_; + delete[] pFirstPartEdgeY_; + delete[] pSecondPartEdgeX_; + delete[] pSecondPartEdgeY_; + delete[] pAnchorX_; + delete[] pAnchorY_; + } + if( pFirstPartEdgeS_ != NULL ) + { + delete[] pFirstPartEdgeS_; + delete[] pSecondPartEdgeS_; + } +} + +int BinaryDescriptor::EDLineDetector::EdgeDrawing( cv::Mat &image, EdgeChains &edgeChains ) +{ + imageWidth = image.cols; + imageHeight = image.rows; + unsigned int pixelNum = imageWidth * imageHeight; + + unsigned int edgePixelArraySize = pixelNum / 5; + unsigned int maxNumOfEdge = edgePixelArraySize / 20; + //compute dx, dy images + if( gImg_.cols != (int) imageWidth || gImg_.rows != (int) imageHeight ) + { + if( pFirstPartEdgeX_ != NULL ) + { + delete[] pFirstPartEdgeX_; + delete[] pFirstPartEdgeY_; + delete[] pSecondPartEdgeX_; + delete[] pSecondPartEdgeY_; + delete[] pFirstPartEdgeS_; + delete[] pSecondPartEdgeS_; + delete[] pAnchorX_; + delete[] pAnchorY_; + } + + dxImg_.create( imageHeight, imageWidth, CV_16SC1 ); + dyImg_.create( imageHeight, imageWidth, CV_16SC1 ); + gImgWO_.create( imageHeight, imageWidth, CV_16SC1 ); + gImg_.create( imageHeight, imageWidth, CV_16SC1 ); + dirImg_.create( imageHeight, imageWidth, CV_8UC1 ); + edgeImage_.create( imageHeight, imageWidth, CV_8UC1 ); + pFirstPartEdgeX_ = new unsigned int[edgePixelArraySize]; + pFirstPartEdgeY_ = new unsigned int[edgePixelArraySize]; + pSecondPartEdgeX_ = new unsigned int[edgePixelArraySize]; + pSecondPartEdgeY_ = new unsigned int[edgePixelArraySize]; + pFirstPartEdgeS_ = new unsigned int[maxNumOfEdge]; + pSecondPartEdgeS_ = new unsigned int[maxNumOfEdge]; + pAnchorX_ = new unsigned int[edgePixelArraySize]; + pAnchorY_ = new unsigned int[edgePixelArraySize]; + } + cv::Sobel( image, dxImg_, CV_16SC1, 1, 0, 3 ); + cv::Sobel( image, dyImg_, CV_16SC1, 0, 1, 3 ); + + //compute gradient and direction images + cv::Mat dxABS_m = cv::abs( dxImg_ ); + cv::Mat dyABS_m = cv::abs( dyImg_ ); + cv::Mat sumDxDy; + cv::add( dyABS_m, dxABS_m, sumDxDy ); + + cv::threshold( sumDxDy, gImg_, gradienThreshold_ + 1, 255, cv::THRESH_TOZERO ); + gImg_ = gImg_ / 4; + gImgWO_ = sumDxDy / 4; + cv::compare( dxABS_m, dyABS_m, dirImg_, cv::CMP_LT ); + + short *pgImg = gImg_.ptr(); + unsigned char *pdirImg = dirImg_.ptr(); + + //extract the anchors in the gradient image, store into a vector + memset( pAnchorX_, 0, edgePixelArraySize * sizeof(unsigned int) ); //initialization + memset( pAnchorY_, 0, edgePixelArraySize * sizeof(unsigned int) ); + unsigned int anchorsSize = 0; + int indexInArray; + unsigned char gValue1, gValue2, gValue3; + for ( unsigned int w = 1; w < imageWidth - 1; w = w + scanIntervals_ ) + { + for ( unsigned int h = 1; h < imageHeight - 1; h = h + scanIntervals_ ) + { + indexInArray = h * imageWidth + w; + //gValue1 = pdirImg[indexInArray]; + if( pdirImg[indexInArray] == Horizontal ) + { //if the direction of pixel is horizontal, then compare with up and down + //gValue2 = pgImg[indexInArray]; + if( pgImg[indexInArray] >= pgImg[indexInArray - imageWidth] + anchorThreshold_ + && pgImg[indexInArray] >= pgImg[indexInArray + imageWidth] + anchorThreshold_ ) + { // (w,h) is accepted as an anchor + pAnchorX_[anchorsSize] = w; + pAnchorY_[anchorsSize++] = h; + } + } + else + { // if(pdirImg[indexInArray]==Vertical){//it is vertical edge, should be compared with left and right + //gValue2 = pgImg[indexInArray]; + if( pgImg[indexInArray] >= pgImg[indexInArray - 1] + anchorThreshold_ && pgImg[indexInArray] >= pgImg[indexInArray + 1] + anchorThreshold_ ) + { // (w,h) is accepted as an anchor + pAnchorX_[anchorsSize] = w; + pAnchorY_[anchorsSize++] = h; + } + } + } + } + if( anchorsSize > edgePixelArraySize ) + { + std::cout << "anchor size is larger than its maximal size. anchorsSize=" << anchorsSize << ", maximal size = " << edgePixelArraySize << std::endl; + return -1; + } + + //link the anchors by smart routing + edgeImage_.setTo( 0 ); + unsigned char *pEdgeImg = edgeImage_.data; + memset( pFirstPartEdgeX_, 0, edgePixelArraySize * sizeof(unsigned int) ); //initialization + memset( pFirstPartEdgeY_, 0, edgePixelArraySize * sizeof(unsigned int) ); + memset( pSecondPartEdgeX_, 0, edgePixelArraySize * sizeof(unsigned int) ); + memset( pSecondPartEdgeY_, 0, edgePixelArraySize * sizeof(unsigned int) ); + memset( pFirstPartEdgeS_, 0, maxNumOfEdge * sizeof(unsigned int) ); + memset( pSecondPartEdgeS_, 0, maxNumOfEdge * sizeof(unsigned int) ); + unsigned int offsetPFirst = 0, offsetPSecond = 0; + unsigned int offsetPS = 0; + + unsigned int x, y; + unsigned int lastX = 0; + unsigned int lastY = 0; + unsigned char lastDirection; //up = 1, right = 2, down = 3, left = 4; + unsigned char shouldGoDirection; //up = 1, right = 2, down = 3, left = 4; + int edgeLenFirst, edgeLenSecond; + for ( unsigned int i = 0; i < anchorsSize; i++ ) + { + x = pAnchorX_[i]; + y = pAnchorY_[i]; + indexInArray = y * imageWidth + x; + if( pEdgeImg[indexInArray] ) + { //if anchor i is already been an edge pixel. + continue; + } + /*The walk stops under 3 conditions: + * 1. We move out of the edge areas, i.e., the thresholded gradient value + * of the current pixel is 0. + * 2. The current direction of the edge changes, i.e., from horizontal + * to vertical or vice versa.?? (This is turned out not correct. From the online edge draw demo + * we can figure out that authors don't implement this rule either because their extracted edge + * chain could be a circle which means pixel directions would definitely be different + * in somewhere on the chain.) + * 3. We encounter a previously detected edge pixel. */ + pFirstPartEdgeS_[offsetPS] = offsetPFirst; + if( pdirImg[indexInArray] == Horizontal ) + { //if the direction of this pixel is horizontal, then go left and right. + //fist go right, pixel direction may be different during linking. + lastDirection = RightDir; + while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) + { + pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel + pFirstPartEdgeX_[offsetPFirst] = x; + pFirstPartEdgeY_[offsetPFirst++] = y; + shouldGoDirection = 0; //unknown + if( pdirImg[indexInArray] == Horizontal ) + { //should go left or right + if( lastDirection == UpDir || lastDirection == DownDir ) + { //change the pixel direction now + if( x > lastX ) + { //should go right + shouldGoDirection = RightDir; + } + else + { //should go left + shouldGoDirection = LeftDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == RightDir || shouldGoDirection == RightDir ) + { //go right + if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the right and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-right + x = x + 1; + y = y + 1; + } + else + { //straight-right + x = x + 1; + } + lastDirection = RightDir; + } + else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) + { //go left + if( x == 0 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the left and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + gValue2 = (unsigned char) pgImg[indexInArray - 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-left + x = x - 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-left + x = x - 1; + } + lastDirection = LeftDir; + } + } + else + { //should go up or down. + if( lastDirection == RightDir || lastDirection == LeftDir ) + { //change the pixel direction now + if( y > lastY ) + { //should go down + shouldGoDirection = DownDir; + } + else + { //should go up + shouldGoDirection = UpDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == DownDir || shouldGoDirection == DownDir ) + { //go down + if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the down and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //down-right + x = x + 1; + y = y + 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-down + y = y + 1; + } + lastDirection = DownDir; + } + else if( lastDirection == UpDir || shouldGoDirection == UpDir ) + { //go up + if( x == 0 || x == imageWidth - 1 || y == 0 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the up and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //up-left + x = x - 1; + y = y - 1; + } + else + { //straight-up + y = y - 1; + } + lastDirection = UpDir; + } + } + indexInArray = y * imageWidth + x; + } //end while go right + //then go left, pixel direction may be different during linking. + x = pAnchorX_[i]; + y = pAnchorY_[i]; + indexInArray = y * imageWidth + x; + pEdgeImg[indexInArray] = 0; //mark the anchor point be a non-edge pixel and + lastDirection = LeftDir; + pSecondPartEdgeS_[offsetPS] = offsetPSecond; + while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) + { + pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel + pSecondPartEdgeX_[offsetPSecond] = x; + pSecondPartEdgeY_[offsetPSecond++] = y; + shouldGoDirection = 0; //unknown + if( pdirImg[indexInArray] == Horizontal ) + { //should go left or right + if( lastDirection == UpDir || lastDirection == DownDir ) + { //change the pixel direction now + if( x > lastX ) + { //should go right + shouldGoDirection = RightDir; + } + else + { //should go left + shouldGoDirection = LeftDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == RightDir || shouldGoDirection == RightDir ) + { //go right + if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the right and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-right + x = x + 1; + y = y + 1; + } + else + { //straight-right + x = x + 1; + } + lastDirection = RightDir; + } + else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) + { //go left + if( x == 0 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the left and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + gValue2 = (unsigned char) pgImg[indexInArray - 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-left + x = x - 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-left + x = x - 1; + } + lastDirection = LeftDir; + } + } + else + { //should go up or down. + if( lastDirection == RightDir || lastDirection == LeftDir ) + { //change the pixel direction now + if( y > lastY ) + { //should go down + shouldGoDirection = DownDir; + } + else + { //should go up + shouldGoDirection = UpDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == DownDir || shouldGoDirection == DownDir ) + { //go down + if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the down and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //down-right + x = x + 1; + y = y + 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-down + y = y + 1; + } + lastDirection = DownDir; + } + else if( lastDirection == UpDir || shouldGoDirection == UpDir ) + { //go up + if( x == 0 || x == imageWidth - 1 || y == 0 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the up and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //up-left + x = x - 1; + y = y - 1; + } + else + { //straight-up + y = y - 1; + } + lastDirection = UpDir; + } + } + indexInArray = y * imageWidth + x; + } //end while go left + //end anchor is Horizontal + } + else + { //the direction of this pixel is vertical, go up and down + //fist go down, pixel direction may be different during linking. + lastDirection = DownDir; + while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) + { + pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel + pFirstPartEdgeX_[offsetPFirst] = x; + pFirstPartEdgeY_[offsetPFirst++] = y; + shouldGoDirection = 0; //unknown + if( pdirImg[indexInArray] == Horizontal ) + { //should go left or right + if( lastDirection == UpDir || lastDirection == DownDir ) + { //change the pixel direction now + if( x > lastX ) + { //should go right + shouldGoDirection = RightDir; + } + else + { //should go left + shouldGoDirection = LeftDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == RightDir || shouldGoDirection == RightDir ) + { //go right + if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the right and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-right + x = x + 1; + y = y + 1; + } + else + { //straight-right + x = x + 1; + } + lastDirection = RightDir; + } + else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) + { //go left + if( x == 0 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the left and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + gValue2 = (unsigned char) pgImg[indexInArray - 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-left + x = x - 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-left + x = x - 1; + } + lastDirection = LeftDir; + } + } + else + { //should go up or down. + if( lastDirection == RightDir || lastDirection == LeftDir ) + { //change the pixel direction now + if( y > lastY ) + { //should go down + shouldGoDirection = DownDir; + } + else + { //should go up + shouldGoDirection = UpDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == DownDir || shouldGoDirection == DownDir ) + { //go down + if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the down and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //down-right + x = x + 1; + y = y + 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-down + y = y + 1; + } + lastDirection = DownDir; + } + else if( lastDirection == UpDir || shouldGoDirection == UpDir ) + { //go up + if( x == 0 || x == imageWidth - 1 || y == 0 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the up and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //up-left + x = x - 1; + y = y - 1; + } + else + { //straight-up + y = y - 1; + } + lastDirection = UpDir; + } + } + indexInArray = y * imageWidth + x; + } //end while go down + //then go up, pixel direction may be different during linking. + lastDirection = UpDir; + x = pAnchorX_[i]; + y = pAnchorY_[i]; + indexInArray = y * imageWidth + x; + pEdgeImg[indexInArray] = 0; //mark the anchor point be a non-edge pixel and + pSecondPartEdgeS_[offsetPS] = offsetPSecond; + while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) + { + pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel + pSecondPartEdgeX_[offsetPSecond] = x; + pSecondPartEdgeY_[offsetPSecond++] = y; + shouldGoDirection = 0; //unknown + if( pdirImg[indexInArray] == Horizontal ) + { //should go left or right + if( lastDirection == UpDir || lastDirection == DownDir ) + { //change the pixel direction now + if( x > lastX ) + { //should go right + shouldGoDirection = RightDir; + } + else + { //should go left + shouldGoDirection = LeftDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == RightDir || shouldGoDirection == RightDir ) + { //go right + if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the right and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-right + x = x + 1; + y = y + 1; + } + else + { //straight-right + x = x + 1; + } + lastDirection = RightDir; + } + else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) + { //go left + if( x == 0 || y == 0 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the left and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + gValue2 = (unsigned char) pgImg[indexInArray - 1]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-left + x = x - 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-left + x = x - 1; + } + lastDirection = LeftDir; + } + } + else + { //should go up or down. + if( lastDirection == RightDir || lastDirection == LeftDir ) + { //change the pixel direction now + if( y > lastY ) + { //should go down + shouldGoDirection = DownDir; + } + else + { //should go up + shouldGoDirection = UpDir; + } + } + lastX = x; + lastY = y; + if( lastDirection == DownDir || shouldGoDirection == DownDir ) + { //go down + if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the down and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //down-right + x = x + 1; + y = y + 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //down-left + x = x - 1; + y = y + 1; + } + else + { //straight-down + y = y + 1; + } + lastDirection = DownDir; + } + else if( lastDirection == UpDir || shouldGoDirection == UpDir ) + { //go up + if( x == 0 || x == imageWidth - 1 || y == 0 ) + { //reach the image border + break; + } + // Look at 3 neighbors to the up and pick the one with the max. gradient value + gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; + gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; + gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; + if( gValue1 >= gValue2 && gValue1 >= gValue3 ) + { //up-right + x = x + 1; + y = y - 1; + } + else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) + { //up-left + x = x - 1; + y = y - 1; + } + else + { //straight-up + y = y - 1; + } + lastDirection = UpDir; + } + } + indexInArray = y * imageWidth + x; + } //end while go up + } //end anchor is Vertical + //only keep the edge chains whose length is larger than the minLineLen_; + edgeLenFirst = offsetPFirst - pFirstPartEdgeS_[offsetPS]; + edgeLenSecond = offsetPSecond - pSecondPartEdgeS_[offsetPS]; + if( edgeLenFirst + edgeLenSecond < minLineLen_ + 1 ) + { //short edge, drop it + offsetPFirst = pFirstPartEdgeS_[offsetPS]; + offsetPSecond = pSecondPartEdgeS_[offsetPS]; + } + else + { + offsetPS++; + } + } + //store the last index + pFirstPartEdgeS_[offsetPS] = offsetPFirst; + pSecondPartEdgeS_[offsetPS] = offsetPSecond; + if( offsetPS > maxNumOfEdge ) + { + std::cout << "Edge drawing Error: The total number of edges is larger than MaxNumOfEdge, " + "numofedge = " << offsetPS << ", MaxNumOfEdge=" << maxNumOfEdge << std::endl; + return -1; + } + if( offsetPFirst > edgePixelArraySize || offsetPSecond > edgePixelArraySize ) + { + std::cout << "Edge drawing Error: The total number of edge pixels is larger than MaxNumOfEdgePixels, " + "numofedgePixel1 = " << offsetPFirst << ", numofedgePixel2 = " << offsetPSecond << ", MaxNumOfEdgePixel=" << edgePixelArraySize << std::endl; + return -1; + } + + /*now all the edge information are stored in pFirstPartEdgeX_, pFirstPartEdgeY_, + *pFirstPartEdgeS_, pSecondPartEdgeX_, pSecondPartEdgeY_, pSecondPartEdgeS_; + *we should reorganize them into edgeChains for easily using. */ + int tempID; + edgeChains.xCors.resize( offsetPFirst + offsetPSecond ); + edgeChains.yCors.resize( offsetPFirst + offsetPSecond ); + edgeChains.sId.resize( offsetPS + 1 ); + unsigned int *pxCors = edgeChains.xCors.data(); + unsigned int *pyCors = edgeChains.yCors.data(); + unsigned int *psId = edgeChains.sId.data(); + offsetPFirst = 0; + offsetPSecond = 0; + unsigned int indexInCors = 0; + unsigned int numOfEdges = 0; + for ( unsigned int edgeId = 0; edgeId < offsetPS; edgeId++ ) + { + //step1, put the first and second parts edge coordinates together from edge start to edge end + psId[numOfEdges++] = indexInCors; + indexInArray = pFirstPartEdgeS_[edgeId]; + offsetPFirst = pFirstPartEdgeS_[edgeId + 1]; + for ( tempID = offsetPFirst - 1; tempID >= indexInArray; tempID-- ) + { //add first part edge + pxCors[indexInCors] = pFirstPartEdgeX_[tempID]; + pyCors[indexInCors++] = pFirstPartEdgeY_[tempID]; + } + indexInArray = pSecondPartEdgeS_[edgeId]; + offsetPSecond = pSecondPartEdgeS_[edgeId + 1]; + for ( tempID = indexInArray + 1; tempID < (int) offsetPSecond; tempID++ ) + { //add second part edge + pxCors[indexInCors] = pSecondPartEdgeX_[tempID]; + pyCors[indexInCors++] = pSecondPartEdgeY_[tempID]; + } + } + psId[numOfEdges] = indexInCors; //the end index of the last edge + edgeChains.numOfEdges = numOfEdges; + + return 1; +} + +int BinaryDescriptor::EDLineDetector::EDline( cv::Mat &image, LineChains &lines ) +{ + + //first, call EdgeDrawing function to extract edges + EdgeChains edges; + if( ( EdgeDrawing( image, edges ) ) != 1 ) + { + std::cout << "Line Detection not finished" << std::endl; + return -1; + } + + //detect lines + unsigned int linePixelID = edges.sId[edges.numOfEdges]; + lines.xCors.resize( linePixelID ); + lines.yCors.resize( linePixelID ); + lines.sId.resize( 5 * edges.numOfEdges ); + unsigned int *pEdgeXCors = edges.xCors.data(); + unsigned int *pEdgeYCors = edges.yCors.data(); + unsigned int *pEdgeSID = edges.sId.data(); + unsigned int *pLineXCors = lines.xCors.data(); + unsigned int *pLineYCors = lines.yCors.data(); + unsigned int *pLineSID = lines.sId.data(); + logNT_ = 2.0 * ( log10( (double) imageWidth ) + log10( (double) imageHeight ) ); + double lineFitErr = 0; //the line fit error; + std::vector lineEquation( 2, 0 ); + lineEquations_.clear(); + lineEndpoints_.clear(); + lineDirection_.clear(); + unsigned char *pdirImg = dirImg_.data; + unsigned int numOfLines = 0; + unsigned int newOffsetS = 0; + unsigned int offsetInEdgeArrayS, offsetInEdgeArrayE; //start index and end index + unsigned int offsetInLineArray = 0; + float direction; //line direction + + for ( unsigned int edgeID = 0; edgeID < edges.numOfEdges; edgeID++ ) + { + offsetInEdgeArrayS = pEdgeSID[edgeID]; + offsetInEdgeArrayE = pEdgeSID[edgeID + 1]; + while ( offsetInEdgeArrayE > offsetInEdgeArrayS + minLineLen_ ) + { //extract line segments from an edge, may find more than one segments + //find an initial line segment + while ( offsetInEdgeArrayE > offsetInEdgeArrayS + minLineLen_ ) + { + lineFitErr = LeastSquaresLineFit_( pEdgeXCors, pEdgeYCors, offsetInEdgeArrayS, lineEquation ); + if( lineFitErr <= lineFitErrThreshold_ ) + break; //ok, an initial line segment detected + offsetInEdgeArrayS += SkipEdgePoint; //skip the first two pixel in the chain and try with the remaining pixels + } + if( lineFitErr > lineFitErrThreshold_ ) + break; //no line is detected + //An initial line segment is detected. Try to extend this line segment + pLineSID[numOfLines] = offsetInLineArray; + double coef1 = 0; //for a line ax+by+c=0, coef1 = 1/sqrt(a^2+b^2); + double pointToLineDis; //for a line ax+by+c=0 and a point(xi, yi), pointToLineDis = coef1*|a*xi+b*yi+c| + bool bExtended = true; + bool bFirstTry = true; + int numOfOutlier; //to against noise, we accept a few outlier of a line. + int tryTimes = 0; + if( pdirImg[pEdgeYCors[offsetInEdgeArrayS] * imageWidth + pEdgeXCors[offsetInEdgeArrayS]] == Horizontal ) + { //y=ax+b, i.e. ax-y+b=0 + while ( bExtended ) + { + tryTimes++; + if( bFirstTry ) + { + bFirstTry = false; + for ( int i = 0; i < minLineLen_; i++ ) + { //First add the initial line segment to the line array + pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; + pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; + } + } + else + { //after each try, line is extended, line equation should be re-estimated + //adjust the line equation + lineFitErr = LeastSquaresLineFit_( pLineXCors, pLineYCors, pLineSID[numOfLines], newOffsetS, offsetInLineArray, lineEquation ); + } + coef1 = 1 / sqrt( lineEquation[0] * lineEquation[0] + 1 ); + numOfOutlier = 0; + newOffsetS = offsetInLineArray; + while ( offsetInEdgeArrayE > offsetInEdgeArrayS ) + { + pointToLineDis = fabs( lineEquation[0] * pEdgeXCors[offsetInEdgeArrayS] - pEdgeYCors[offsetInEdgeArrayS] + lineEquation[1] ) * coef1; + pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; + pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; + if( pointToLineDis > lineFitErrThreshold_ ) + { + numOfOutlier++; + if( numOfOutlier > 3 ) + break; + } + else + { //we count number of connective outliers. + numOfOutlier = 0; + } + } + //pop back the last few outliers from lines and return them to edge chain + offsetInLineArray -= numOfOutlier; + offsetInEdgeArrayS -= numOfOutlier; + if( offsetInLineArray - newOffsetS > 0 && tryTimes < TryTime ) + { //some new pixels are added to the line + } + else + { + bExtended = false; //no new pixels are added. + } + } + //the line equation coefficients,for line w1x+w2y+w3 =0, we normalize it to make w1^2+w2^2 = 1. + std::vector lineEqu( 3, 0 ); + lineEqu[0] = lineEquation[0] * coef1; + lineEqu[1] = -1 * coef1; + lineEqu[2] = lineEquation[1] * coef1; + if( LineValidation_( pLineXCors, pLineYCors, pLineSID[numOfLines], offsetInLineArray, lineEqu, direction ) ) + { //check the line + //store the line equation coefficients + lineEquations_.push_back( lineEqu ); + /*At last, compute the line endpoints and store them. + *we project the first and last pixels in the pixelChain onto the best fit line + *to get the line endpoints. + *xp= (w2^2*x0-w1*w2*y0-w3*w1)/(w1^2+w2^2) + *yp= (w1^2*y0-w1*w2*x0-w3*w2)/(w1^2+w2^2) */ + std::vector lineEndP( 4, 0 ); //line endpoints + double a1 = lineEqu[1] * lineEqu[1]; + double a2 = lineEqu[0] * lineEqu[0]; + double a3 = lineEqu[0] * lineEqu[1]; + double a4 = lineEqu[2] * lineEqu[0]; + double a5 = lineEqu[2] * lineEqu[1]; + unsigned int Px = pLineXCors[pLineSID[numOfLines]]; //first pixel + unsigned int Py = pLineYCors[pLineSID[numOfLines]]; + lineEndP[0] = (float) ( a1 * Px - a3 * Py - a4 ); //x + lineEndP[1] = (float) ( a2 * Py - a3 * Px - a5 ); //y + Px = pLineXCors[offsetInLineArray - 1]; //last pixel + Py = pLineYCors[offsetInLineArray - 1]; + lineEndP[2] = (float) ( a1 * Px - a3 * Py - a4 ); //x + lineEndP[3] = (float) ( a2 * Py - a3 * Px - a5 ); //y + lineEndpoints_.push_back( lineEndP ); + lineDirection_.push_back( direction ); + numOfLines++; + } + else + { + offsetInLineArray = pLineSID[numOfLines]; // line was not accepted, the offset is set back + } + } + else + { //x=ay+b, i.e. x-ay-b=0 + while ( bExtended ) + { + tryTimes++; + if( bFirstTry ) + { + bFirstTry = false; + for ( int i = 0; i < minLineLen_; i++ ) + { //First add the initial line segment to the line array + pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; + pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; + } + } + else + { //after each try, line is extended, line equation should be re-estimated + //adjust the line equation + lineFitErr = LeastSquaresLineFit_( pLineXCors, pLineYCors, pLineSID[numOfLines], newOffsetS, offsetInLineArray, lineEquation ); + } + coef1 = 1 / sqrt( 1 + lineEquation[0] * lineEquation[0] ); + numOfOutlier = 0; + newOffsetS = offsetInLineArray; + while ( offsetInEdgeArrayE > offsetInEdgeArrayS ) + { + pointToLineDis = fabs( pEdgeXCors[offsetInEdgeArrayS] - lineEquation[0] * pEdgeYCors[offsetInEdgeArrayS] - lineEquation[1] ) * coef1; + pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; + pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; + if( pointToLineDis > lineFitErrThreshold_ ) + { + numOfOutlier++; + if( numOfOutlier > 3 ) + break; + } + else + { //we count number of connective outliers. + numOfOutlier = 0; + } + } + //pop back the last few outliers from lines and return them to edge chain + offsetInLineArray -= numOfOutlier; + offsetInEdgeArrayS -= numOfOutlier; + if( offsetInLineArray - newOffsetS > 0 && tryTimes < TryTime ) + { //some new pixels are added to the line + } + else + { + bExtended = false; //no new pixels are added. + } + } + //the line equation coefficients,for line w1x+w2y+w3 =0, we normalize it to make w1^2+w2^2 = 1. + std::vector lineEqu( 3, 0 ); + lineEqu[0] = 1 * coef1; + lineEqu[1] = -lineEquation[0] * coef1; + lineEqu[2] = -lineEquation[1] * coef1; + + if( LineValidation_( pLineXCors, pLineYCors, pLineSID[numOfLines], offsetInLineArray, lineEqu, direction ) ) + { //check the line + //store the line equation coefficients + lineEquations_.push_back( lineEqu ); + /*At last, compute the line endpoints and store them. + *we project the first and last pixels in the pixelChain onto the best fit line + *to get the line endpoints. + *xp= (w2^2*x0-w1*w2*y0-w3*w1)/(w1^2+w2^2) + *yp= (w1^2*y0-w1*w2*x0-w3*w2)/(w1^2+w2^2) */ + std::vector lineEndP( 4, 0 ); //line endpoints + double a1 = lineEqu[1] * lineEqu[1]; + double a2 = lineEqu[0] * lineEqu[0]; + double a3 = lineEqu[0] * lineEqu[1]; + double a4 = lineEqu[2] * lineEqu[0]; + double a5 = lineEqu[2] * lineEqu[1]; + unsigned int Px = pLineXCors[pLineSID[numOfLines]]; //first pixel + unsigned int Py = pLineYCors[pLineSID[numOfLines]]; + lineEndP[0] = (float) ( a1 * Px - a3 * Py - a4 ); //x + lineEndP[1] = (float) ( a2 * Py - a3 * Px - a5 ); //y + Px = pLineXCors[offsetInLineArray - 1]; //last pixel + Py = pLineYCors[offsetInLineArray - 1]; + lineEndP[2] = (float) ( a1 * Px - a3 * Py - a4 ); //x + lineEndP[3] = (float) ( a2 * Py - a3 * Px - a5 ); //y + lineEndpoints_.push_back( lineEndP ); + lineDirection_.push_back( direction ); + numOfLines++; + } + else + { + offsetInLineArray = pLineSID[numOfLines]; // line was not accepted, the offset is set back + } + } + //Extract line segments from the remaining pixel; Current chain has been shortened already. + } + } //end for(unsigned int edgeID=0; edgeID &lineEquation ) +{ + + float * pMatT; + float * pATA; + double fitError = 0; + double coef; + unsigned char *pdirImg = dirImg_.data; + unsigned int offset = offsetS; + /*If the first pixel in this chain is horizontal, + *then we try to find a horizontal line, y=ax+b;*/ + if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Horizontal ) + { + /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec + * [x0,1] [y0] + * [x1,1] [a] [y1] + * . [b] = . + * [xn,1] [yn]*/ + pMatT = fitMatT.ptr(); //fitMatT = [x0, x1, ... xn; 1,1,...,1]; + for ( int i = 0; i < minLineLen_; i++ ) + { + //*(pMatT+minLineLen_) = 1; //the value are not changed; + * ( pMatT++ ) = (float) xCors[offsetS]; + fitVec[0][i] = (float) yCors[offsetS++]; + } + ATA = fitMatT * fitMatT.t(); + ATV = fitMatT * fitVec.t(); + /* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */ + pATA = ATA.ptr(); + coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); + // lineEquation = svd.Invert(ATA) * matT * vec; + lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); + lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); + /*compute line fit error */ + for ( int i = 0; i < minLineLen_; i++ ) + { + //coef = double(yCors[offset]) - double(xCors[offset++]) * lineEquation[0] - lineEquation[1]; + coef = double( yCors[offset] ) - double( xCors[offset] ) * lineEquation[0] - lineEquation[1]; + offset++; + fitError += coef * coef; + } + return sqrt( fitError ); + } + /*If the first pixel in this chain is vertical, + *then we try to find a vertical line, x=ay+b;*/ + if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Vertical ) + { + /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec + * [y0,1] [x0] + * [y1,1] [a] [x1] + * . [b] = . + * [yn,1] [xn]*/ + pMatT = fitMatT.ptr(); //fitMatT = [y0, y1, ... yn; 1,1,...,1]; + for ( int i = 0; i < minLineLen_; i++ ) + { + //*(pMatT+minLineLen_) = 1;//the value are not changed; + * ( pMatT++ ) = (float) yCors[offsetS]; + fitVec[0][i] = (float) xCors[offsetS++]; + } + ATA = fitMatT * ( fitMatT.t() ); + ATV = fitMatT * fitVec.t(); + /* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */ + pATA = ATA.ptr(); + coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); + // lineEquation = svd.Invert(ATA) * matT * vec; + lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); + lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); + /*compute line fit error */ + for ( int i = 0; i < minLineLen_; i++ ) + { + //coef = double(xCors[offset]) - double(yCors[offset++]) * lineEquation[0] - lineEquation[1]; + coef = double( xCors[offset] ) - double( yCors[offset] ) * lineEquation[0] - lineEquation[1]; + offset++; + fitError += coef * coef; + } + return sqrt( fitError ); + } + return 0; +} +double BinaryDescriptor::EDLineDetector::LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, + unsigned int newOffsetS, unsigned int offsetE, std::vector &lineEquation ) +{ + int length = offsetE - offsetS; + int newLength = offsetE - newOffsetS; + if( length <= 0 || newLength <= 0 ) + { + std::cout << "EDLineDetector::LeastSquaresLineFit_ Error:" + " the expected line index is wrong...offsetE = " << offsetE << ", offsetS=" << offsetS << ", newOffsetS=" << newOffsetS << std::endl; + return -1; + } + if( lineEquation.size() != 2 ) + { + std::cout << "SHOULD NOT BE != 2" << std::endl; + } + cv::Mat_ matT( 2, newLength ); + cv::Mat_ vec( newLength, 1 ); + float * pMatT; + float * pATA; + double coef; + unsigned char *pdirImg = dirImg_.data; + /*If the first pixel in this chain is horizontal, + *then we try to find a horizontal line, y=ax+b;*/ + if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Horizontal ) + { + /*Build the new system,and solve it using least square regression: mat * [a,b]^T = vec + * [x0',1] [y0'] + * [x1',1] [a] [y1'] + * . [b] = . + * [xn',1] [yn']*/ + pMatT = matT.ptr(); //matT = [x0', x1', ... xn'; 1,1,...,1] + for ( int i = 0; i < newLength; i++ ) + { + * ( pMatT + newLength ) = 1; + * ( pMatT++ ) = (float) xCors[newOffsetS]; + vec[0][i] = (float) yCors[newOffsetS++]; + } + /* [a,b]^T = Inv(ATA + mat^T * mat) * (ATV + mat^T * vec) */ + tempMatLineFit = matT * matT.t(); + tempVecLineFit = matT * vec; + ATA = ATA + tempMatLineFit; + ATV = ATV + tempVecLineFit; + pATA = ATA.ptr(); + coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); + lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); + lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); + + return 0; + } + /*If the first pixel in this chain is vertical, + *then we try to find a vertical line, x=ay+b;*/ + if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Vertical ) + { + /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec + * [y0',1] [x0'] + * [y1',1] [a] [x1'] + * . [b] = . + * [yn',1] [xn']*/ + pMatT = matT.ptr(); //matT = [y0', y1', ... yn'; 1,1,...,1] + for ( int i = 0; i < newLength; i++ ) + { + * ( pMatT + newLength ) = 1; + * ( pMatT++ ) = (float) yCors[newOffsetS]; + vec[0][i] = (float) xCors[newOffsetS++]; + } + /* [a,b]^T = Inv(ATA + mat^T * mat) * (ATV + mat^T * vec) */ +// matT.MultiplyWithTransposeOf(matT, tempMatLineFit); + tempMatLineFit = matT * matT.t(); + tempVecLineFit = matT * vec; + ATA = ATA + tempMatLineFit; + ATV = ATV + tempVecLineFit; +// pATA = ATA.GetData(); + pATA = ATA.ptr(); + coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); + lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); + lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); + + } + return 0; +} + +bool BinaryDescriptor::EDLineDetector::LineValidation_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int offsetE, + std::vector &lineEquation, float &direction ) +{ + if( bValidate_ ) + { + int n = offsetE - offsetS; + /*first compute the direction of line, make sure that the dark side always be the + *left side of a line.*/ + int meanGradientX = 0, meanGradientY = 0; + short *pdxImg = dxImg_.ptr(); + short *pdyImg = dyImg_.ptr(); + double dx, dy; + std::vector pointDirection; + int index; + for ( int i = 0; i < n; i++ ) + { + index = yCors[offsetS] * imageWidth + xCors[offsetS]; + offsetS++; + meanGradientX += pdxImg[index]; + meanGradientY += pdyImg[index]; + dx = (double) pdxImg[index]; + dy = (double) pdyImg[index]; + pointDirection.push_back( atan2( -dx, dy ) ); + } + dx = fabs( lineEquation[1] ); + dy = fabs( lineEquation[0] ); + if( meanGradientX == 0 && meanGradientY == 0 ) + { //not possible, if happens, it must be a wrong line, + return false; + } + if( meanGradientX > 0 && meanGradientY >= 0 ) + { //first quadrant, and positive direction of X axis. + direction = (float) atan2( -dy, dx ); //line direction is in fourth quadrant + } + if( meanGradientX <= 0 && meanGradientY > 0 ) + { //second quadrant, and positive direction of Y axis. + direction = (float) atan2( dy, dx ); //line direction is in first quadrant + } + if( meanGradientX < 0 && meanGradientY <= 0 ) + { //third quadrant, and negative direction of X axis. + direction = (float) atan2( dy, -dx ); //line direction is in second quadrant + } + if( meanGradientX >= 0 && meanGradientY < 0 ) + { //fourth quadrant, and negative direction of Y axis. + direction = (float) atan2( -dy, -dx ); //line direction is in third quadrant + } + /*then check whether the line is on the border of the image. We don't keep the border line.*/ + if( fabs( direction ) < 0.15 || M_PI - fabs( direction ) < 0.15 ) + { //Horizontal line + if( fabs( lineEquation[2] ) < 10 || fabs( imageHeight - fabs( lineEquation[2] ) ) < 10 ) + { //upper border or lower border + return false; + } + } + if( fabs( fabs( direction ) - M_PI * 0.5 ) < 0.15 ) + { //Vertical line + if( fabs( lineEquation[2] ) < 10 || fabs( imageWidth - fabs( lineEquation[2] ) ) < 10 ) + { //left border or right border + return false; + } + } + //count the aligned points on the line which have the same direction as the line. + double disDirection; + int k = 0; + for ( int i = 0; i < n; i++ ) + { + disDirection = fabs( direction - pointDirection[i] ); + if( fabs( 2 * M_PI - disDirection ) < 0.392699 || disDirection < 0.392699 ) + { //same direction, pi/8 = 0.392699081698724 + k++; + } + } + //now compute NFA(Number of False Alarms) + double ret = nfa( n, k, 0.125, logNT_ ); + + return ( ret > 0 ); //0 corresponds to 1 mean false alarm + } + else + { + return true; + } +} + +int BinaryDescriptor::EDLineDetector::EDline( cv::Mat &image ) +{ + if( ( EDline( image, lines_/*, smoothed*/) ) != 1 ) + { + return -1; + } + lineSalience_.clear(); + lineSalience_.resize( lines_.numOfLines ); + unsigned char *pgImg = gImgWO_.ptr(); + unsigned int indexInLineArray; + unsigned int *pXCor = lines_.xCors.data(); + unsigned int *pYCor = lines_.yCors.data(); + unsigned int *pSID = lines_.sId.data(); + for ( unsigned int i = 0; i < lineSalience_.size(); i++ ) + { + int salience = 0; + for ( indexInLineArray = pSID[i]; indexInLineArray < pSID[i + 1]; indexInLineArray++ ) + { + salience += pgImg[pYCor[indexInLineArray] * imageWidth + pXCor[indexInLineArray]]; + } + lineSalience_[i] = (float) salience; + } + return 1; +} + } } diff --git a/modules/line_descriptor/src/binary_descriptor_matcher.cpp b/modules/line_descriptor/src/binary_descriptor_matcher.cpp index 7f044a4f4..268a4dc4b 100644 --- a/modules/line_descriptor/src/binary_descriptor_matcher.cpp +++ b/modules/line_descriptor/src/binary_descriptor_matcher.cpp @@ -41,6 +41,10 @@ #include "precomp.hpp" +#define MAX_B 37 +double ARRAY_RESIZE_FACTOR = 1.1; // minimum is 1.0 +double ARRAY_RESIZE_ADD_FACTOR = 4; // minimum is 1 + //using namespace cv; namespace cv { @@ -59,7 +63,7 @@ BinaryDescriptorMatcher::BinaryDescriptorMatcher() /* constructor with smart pointer */ Ptr BinaryDescriptorMatcher::createBinaryDescriptorMatcher() { - return Ptr( new BinaryDescriptorMatcher() ); + return Ptr < BinaryDescriptorMatcher > ( new BinaryDescriptorMatcher() ); } /* store new descriptors to be inserted in dataset */ @@ -132,7 +136,7 @@ void BinaryDescriptorMatcher::match( const Mat& queryDescriptors, std::vectorsecond << " in knnMatch function " << "should have " << queryDescriptors.rows << " and " - << "1 column. Program will be terminated"; + << "1 column. Program will be terminated"; //throw std::runtime_error( ss.str() ); } /* create a DMatch object if required by mask or if there is no mask at all */ - else if( masks.empty() || masks[itup->second].at( counter ) != 0 ) + else if( masks.empty() || masks[itup->second].at < uchar > ( counter ) != 0 ) { std::vector k_distances; checkKDistances( numres, 1, k_distances, counter, 256 ); @@ -227,7 +231,7 @@ void BinaryDescriptorMatcher::match( const Mat& queryDescriptors, const Mat& tra { /* create a DMatch object if required by mask or if there is no mask at all */ - if( mask.empty() || ( !mask.empty() && mask.at( counter ) != 0 ) ) + if( mask.empty() || ( !mask.empty() && mask.at < uchar > ( counter ) != 0 ) ) { std::vector k_distances; checkKDistances( numres, 1, k_distances, counter, 256 ); @@ -291,10 +295,10 @@ void BinaryDescriptorMatcher::knnMatch( const Mat& queryDescriptors, const Mat& for ( int counter = 0; counter < queryDescriptors.rows; counter++ ) { /* initialize a vector of matches */ - std::vector tempVec; + std::vector < DMatch > tempVec; /* chech whether query should be ignored */ - if( !mask.empty() && mask.at( counter ) == 0 ) + if( !mask.empty() && mask.at < uchar > ( counter ) == 0 ) { /* if compact result is not requested, add an empty vector */ if( !compactResult ) @@ -346,7 +350,7 @@ void BinaryDescriptorMatcher::knnMatch( const Mat& queryDescriptors, std::vector if( masks.size() != 0 && (int) masks.size() != numImages ) { std::cout << "Error: the number of images in dataset is " << numImages << " but knnMatch function received " << masks.size() - << " masks. Program will be terminated" << std::endl; + << " masks. Program will be terminated" << std::endl; return; } @@ -369,7 +373,7 @@ void BinaryDescriptorMatcher::knnMatch( const Mat& queryDescriptors, std::vector for ( int counter = 0; counter < queryDescriptors.rows; counter++ ) { /* create a void vector of matches */ - std::vector tempVector; + std::vector < DMatch > tempVector; /* loop over k results returned for every query */ for ( int j = index; j < index + k; j++ ) @@ -384,14 +388,14 @@ void BinaryDescriptorMatcher::knnMatch( const Mat& queryDescriptors, std::vector if( !masks.empty() && ( masks[itup->second].rows != queryDescriptors.rows || masks[itup->second].cols != 1 ) ) { std::cout << "Error: mask " << itup->second << " in knnMatch function " << "should have " << queryDescriptors.rows << " and " - << "1 column. Program will be terminated" << std::endl; + << "1 column. Program will be terminated" << std::endl; return; } /* decide if, according to relative mask, returned match should be considered */ - else if( masks.size() == 0 || masks[itup->second].at( counter ) != 0 ) + else if( masks.size() == 0 || masks[itup->second].at < uchar > ( counter ) != 0 ) { std::vector k_distances; checkKDistances( numres, k, k_distances, counter, 256 ); @@ -465,13 +469,13 @@ void BinaryDescriptorMatcher::radiusMatch( const Mat& queryDescriptors, const Ma std::vector k_distances; checkKDistances( numres, trainDescriptors.rows, k_distances, i, 256 ); - std::vector tempVector; + std::vector < DMatch > tempVector; for ( int j = index; j < index + trainDescriptors.rows; j++ ) { // if( numres[j] <= maxDistance ) if( k_distances[j - index] <= maxDistance ) { - if( mask.empty() || mask.at( i ) != 0 ) + if( mask.empty() || mask.at < uchar > ( i ) != 0 ) { DMatch dm; dm.queryIdx = i; @@ -515,7 +519,7 @@ void BinaryDescriptorMatcher::radiusMatch( const Mat& queryDescriptors, std::vec if( masks.size() != 0 && (int) masks.size() != numImages ) { std::cout << "Error: the number of images in dataset is " << numImages << " but radiusMatch function received " << masks.size() - << " masks. Program will be terminated" << std::endl; + << " masks. Program will be terminated" << std::endl; return; } @@ -537,7 +541,7 @@ void BinaryDescriptorMatcher::radiusMatch( const Mat& queryDescriptors, std::vec int index = 0; for ( int counter = 0; counter < queryDescriptors.rows; counter++ ) { - std::vector tempVector; + std::vector < DMatch > tempVector; for ( int j = index; j < index + descrInDS; j++ ) { std::vector k_distances; @@ -554,13 +558,13 @@ void BinaryDescriptorMatcher::radiusMatch( const Mat& queryDescriptors, std::vec if( !masks.empty() && ( masks[itup->second].rows != queryDescriptors.rows || masks[itup->second].cols != 1 ) ) { std::cout << "Error: mask " << itup->second << " in radiusMatch function " << "should have " << queryDescriptors.rows << " and " - << "1 column. Program will be terminated" << std::endl; + << "1 column. Program will be terminated" << std::endl; return; } /* add match if necessary */ - else if( masks.empty() || masks[itup->second].at( counter ) != 0 ) + else if( masks.empty() || masks[itup->second].at < uchar > ( counter ) != 0 ) { DMatch dm; @@ -587,5 +591,477 @@ void BinaryDescriptorMatcher::radiusMatch( const Mat& queryDescriptors, std::vec delete numres; } + +/* execute a batch query */ +void BinaryDescriptorMatcher::Mihasher::batchquery( UINT32 * results, UINT32 *numres, const cv::Mat & queries, UINT32 numq, int dim1queries ) +{ + /* create and initialize a bitarray */ + counter = new bitarray; + counter->init( N ); + + UINT32 *res = new UINT32[K * ( D + 1 )]; + UINT64 *chunks = new UINT64[m]; + UINT32 * presults = results; + UINT32 *pnumres = numres; + + /* make a copy of input queries */ + cv::Mat queries_clone = queries.clone(); + + /* set a pointer to first query (row) */ + UINT8 *pq = queries_clone.ptr(); + + /* loop over number of descriptors */ + for ( size_t i = 0; i < numq; i++ ) + { + /* for every descriptor, query database */ + query( presults, pnumres, pq, chunks, res ); + + /* move pointer to write next K indeces */ + presults += K; + pnumres += B + 1; + + /* move forward pointer to current row in descriptors matrix */ + pq += dim1queries; + + } + + delete[] res; + delete[] chunks; + + delete counter; +} + +/* execute a single query */ +void BinaryDescriptorMatcher::Mihasher::query( UINT32* results, UINT32* numres, UINT8 * Query, UINT64 *chunks, UINT32 *res ) +{ + /* if K == 0 that means we want everything to be processed. + So maxres = N in that case. Otherwise K limits the results processed */ + UINT32 maxres = K ? K : (UINT32) N; + + /* number of results so far obtained (up to a distance of s per chunk) */ + UINT32 n = 0; + + /* number of candidates tested with full codes (not counting duplicates) */ + UINT32 nc = 0; + + /* counting everything retrieved (duplicates are counted multiple times) + number of lookups (and xors) */ + UINT32 nl = 0; + + UINT32 nd = 0; + UINT32 *arr; + int size = 0; + UINT32 index; + int hammd; + + counter->erase(); + memset( numres, 0, ( B + 1 ) * sizeof ( *numres ) ); + + split( chunks, Query, m, mplus, b ); + + /* the growing search radius per substring */ + int s; + + /* current b: for the first mplus substrings it is b, for the rest it is (b-1) */ + int curb = b; + + for ( s = 0; s <= d && n < maxres; s++ ) + { + for ( int k = 0; k < m; k++ ) + { + if( k < mplus ) + curb = b; + else + curb = b - 1; + UINT64 chunksk = chunks[k]; + /* number of bit-strings with s number of 1s */ + nl += xornum[s + 1] - xornum[s]; + + /* the bit-string with s number of 1s */ + UINT64 bitstr = 0; + for ( int i = 0; i < s; i++ ) + /* power[i] stores the location of the i'th 1 */ + power[i] = i; + /* used for stopping criterion (location of (s+1)th 1) */ + power[s] = curb + 1; + + /* bit determines the 1 that should be moving to the left */ + int bit = s - 1; + + /* start from the left-most 1, and move it to the left until + it touches another one */ + + /* the loop for changing bitstr */ + bool infiniteWhile = true; + while ( infiniteWhile ) + { + if( bit != -1 ) + { + bitstr ^= ( power[bit] == bit ) ? (UINT64) 1 << power[bit] : (UINT64) 3 << ( power[bit] - 1 ); + power[bit]++; + bit--; + } + + else + { /* bit == -1 */ + /* the binary code bitstr is available for processing */ + arr = H[k].query( chunksk ^ bitstr, &size ); // lookup + if( size ) + { /* the corresponding bucket is not empty */ + nd += size; + for ( int c = 0; c < size; c++ ) + { + index = arr[c]; + if( !counter->get( index ) ) + { /* if it is not a duplicate */ + counter->set( index ); + hammd = cv::line_descriptor::match( codes.ptr() + (UINT64) index * ( B_over_8 ), Query, B_over_8 ); + + nc++; + if( hammd <= D && numres[hammd] < maxres ) + res[hammd * K + numres[hammd]] = index + 1; + + numres[hammd]++; + } + } + } + + /* end of processing */ + while ( ++bit < s && power[bit] == power[bit + 1] - 1 ) + { + bitstr ^= (UINT64) 1 << ( power[bit] - 1 ); + power[bit] = bit; + } + if( bit == s ) + break; + } + } + + n = n + numres[s * m + k]; + if( n >= maxres ) + break; + } + } + + n = 0; + for ( s = 0; s <= D && (int) n < K; s++ ) + { + for ( int c = 0; c < (int) numres[s] && (int) n < K; c++ ) + results[n++] = res[s * K + c]; + } + +} + +/* constructor 2 */ +BinaryDescriptorMatcher::Mihasher::Mihasher( int _B, int _m ) +{ + B = _B; + B_over_8 = B / 8; + m = _m; + b = (int) ceil( (double) B / m ); + + /* assuming that B/2 is large enough radius to include + all of the k nearest neighbors */ + D = (int) ceil( B / 2.0 ); + d = (int) ceil( (double) D / m ); + + /* mplus is the number of chunks with b bits + (m-mplus) is the number of chunks with (b-1) bits */ + mplus = B - m * ( b - 1 ); + + xornum = new UINT32[d + 2]; + xornum[0] = 0; + for ( int i = 0; i <= d; i++ ) + xornum[i + 1] = xornum[i] + (UINT32) choose( b, i ); + + H = new SparseHashtable[m]; + + /* H[i].init might fail */ + for ( int i = 0; i < mplus; i++ ) + H[i].init( b ); + for ( int i = mplus; i < m; i++ ) + H[i].init( b - 1 ); +} + +/* K setter */ +void BinaryDescriptorMatcher::Mihasher::setK( int _K ) +{ + K = _K; +} + +/* desctructor */ +BinaryDescriptorMatcher::Mihasher::~Mihasher() +{ + delete[] xornum; + delete[] H; +} + +/* populate tables */ +void BinaryDescriptorMatcher::Mihasher::populate( cv::Mat & _codes, UINT32 _N, int dim1codes ) +{ + N = _N; + codes = _codes; + UINT64 * chunks = new UINT64[m]; + + UINT8 * pcodes = codes.ptr(); + for ( UINT64 i = 0; i < N; i++, pcodes += dim1codes ) + { + split( chunks, pcodes, m, mplus, b ); + + for ( int k = 0; k < m; k++ ) + H[k].insert( chunks[k], (UINT32) i ); + + if( i % (int) ceil( N / 1000.0 ) == 0 ) + fflush (stdout); + } + + delete[] chunks; +} + +/* constructor */ +BinaryDescriptorMatcher::SparseHashtable::SparseHashtable() +{ + table = NULL; + size = 0; + b = 0; +} + +/* initializer */ +int BinaryDescriptorMatcher::SparseHashtable::init( int _b ) +{ + b = _b; + + if( b < 5 || b > MAX_B || b > (int) ( sizeof(UINT64) * 8 ) ) + return 1; + + size = UINT64_1 << ( b - 5 ); // size = 2 ^ b + table = (BucketGroup*) calloc( size, sizeof(BucketGroup) ); + + return 0; + +} + +/* destructor */ +BinaryDescriptorMatcher::SparseHashtable::~SparseHashtable() +{ + free( table ); +} + +/* insert data */ +void BinaryDescriptorMatcher::SparseHashtable::insert( UINT64 index, UINT32 data ) +{ + table[index >> 5].insert( (int) ( index % 32 ), data ); +} + +/* query data */ +UINT32* BinaryDescriptorMatcher::SparseHashtable::query( UINT64 index, int *Size ) +{ + return table[index >> 5].query( (int) ( index % 32 ), Size ); +} + +/* constructor */ +BinaryDescriptorMatcher::BucketGroup::BucketGroup() +{ + empty = 0; + group = NULL; +} + +/* destructor */ +BinaryDescriptorMatcher::BucketGroup::~BucketGroup() +{ + if( group != NULL ) + delete group; +} + +/* insert data into the bucket */ +void BinaryDescriptorMatcher::BucketGroup::insert( int subindex, UINT32 data ) +{ + if( group == NULL ) + { + group = new Array32(); + group->push( 0 ); + } + + UINT32 lowerbits = ( (UINT32) 1 << subindex ) - 1; + int end = popcnt( empty & lowerbits ); + + if( ! ( empty & ( (UINT32) 1 << subindex ) ) ) + { + group->insert( end, group->arr[end + 2] ); + empty |= (UINT32) 1 << subindex; + } + + int totones = popcnt( empty ); + group->insert( totones + 1 + group->arr[2 + end + 1], data ); + for ( int i = end + 1; i < totones + 1; i++ ) + group->arr[2 + i]++; +} + +/* perform a query to the bucket */ +UINT32* BinaryDescriptorMatcher::BucketGroup::query( int subindex, int *size ) +{ + if( empty & ( (UINT32) 1 << subindex ) ) + { + UINT32 lowerbits = ( (UINT32) 1 << subindex ) - 1; + int end = popcnt( empty & lowerbits ); + int totones = popcnt( empty ); + *size = group->arr[2 + end + 1] - group->arr[2 + end]; + return group->arr + 2 + totones + 1 + group->arr[2 + end]; + } + + else + { + *size = 0; + return NULL; + } +} + +/* set ARRAY_RESIZE_FACTOR */ +void BinaryDescriptorMatcher::Array32::setArrayResizeFactor( double arf ) +{ + ARRAY_RESIZE_FACTOR = arf; +} + +/* constructor */ +BinaryDescriptorMatcher::Array32::Array32() +{ + arr = NULL; + ARRAY_RESIZE_FACTOR = 1.1; // minimum is 1.0 + ARRAY_RESIZE_ADD_FACTOR = 4; // minimum is 1 + +} + +/* definition of operator = + Array32& Array32::operator = (const Array32 &rhs) */ +void BinaryDescriptorMatcher::Array32::operator =( const Array32 &rhs ) +{ + if( &rhs != this ) + this->arr = rhs.arr; +} + +/* destructor */ +BinaryDescriptorMatcher::Array32::~Array32() +{ + cleanup(); +} + +/* cleaning function used in destructor */ +void BinaryDescriptorMatcher::Array32::cleanup() +{ + free( arr ); +} + +/* push data */ +void BinaryDescriptorMatcher::Array32::push( UINT32 Data ) +{ + if( arr ) + { + if( arr[0] == arr[1] ) + { + arr[1] = (UINT32) std::max( ceil( arr[1] * ARRAY_RESIZE_FACTOR ), arr[1] + ARRAY_RESIZE_ADD_FACTOR ); + UINT32* new_Data = static_cast( realloc( arr, sizeof(UINT32) * ( 2 + arr[1] ) ) ); + if( new_Data == NULL ) + { + /* could not realloc, but orig still valid */ + std::cout << "ALERT!!!! Not enough memory, operation aborted!" << std::endl; + exit( 0 ); + } + else + { + arr = new_Data; + } + + } + + arr[2 + arr[0]] = Data; + arr[0]++; + + } + + else + { + arr = (UINT32*) malloc( ( size_t )( 2 + ARRAY_RESIZE_ADD_FACTOR ) * sizeof(UINT32) ); + arr[0] = 1; + arr[1] = 1; + arr[2] = Data; + } +} + +/* insert data at given index */ +void BinaryDescriptorMatcher::Array32::insert( UINT32 index, UINT32 Data ) +{ + if( arr ) + { + if( arr[0] == arr[1] ) + { + arr[1] = (UINT32) ceil( arr[0] * 1.1 ); + UINT32* new_data = static_cast( realloc( arr, sizeof(UINT32) * ( 2 + arr[1] ) ) ); + if( new_data == NULL ) + { + // could not realloc, but orig still valid + std::cout << "ALERT!!!! Not enough memory, operation aborted!" << std::endl; + exit( 0 ); + } + else + { + arr = new_data; + } + } + + memmove( arr + ( 2 + index ) + 1, arr + ( 2 + index ), ( arr[0] - index ) * sizeof ( *arr ) ); + + arr[2 + index] = Data; + arr[0]++; + } + + else + { + arr = (UINT32*) malloc( 3 * sizeof(UINT32) ); + arr[0] = 1; + arr[1] = 1; + arr[2] = Data; + } +} + +/* return data */ +UINT32* BinaryDescriptorMatcher::Array32::data() +{ + return arr ? arr + 2 : NULL; +} + +/* return data size */ +UINT32 BinaryDescriptorMatcher::Array32::size() +{ + return arr ? arr[0] : 0; +} + +/* return capacity */ +UINT32 BinaryDescriptorMatcher::Array32::capacity() +{ + return arr ? arr[1] : 0; +} + +/* print data */ +void BinaryDescriptorMatcher::Array32::print() +{ + for ( int i = 0; i < (int) size(); i++ ) + printf( "%d, ", arr[i + 2] ); + + printf( "\n" ); +} + +/* initializer */ +void BinaryDescriptorMatcher::Array32::init( int Size ) +{ + if( arr == NULL ) + { + arr = (UINT32*) malloc( ( 2 + Size ) * sizeof(UINT32) ); + arr[0] = 0; + arr[1] = Size; + } } + + } +} + diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/bitarray.hpp b/modules/line_descriptor/src/bitarray.hpp similarity index 100% rename from modules/line_descriptor/include/opencv2/line_descriptor/bitarray.hpp rename to modules/line_descriptor/src/bitarray.hpp diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/bitops.hpp b/modules/line_descriptor/src/bitops.hpp similarity index 98% rename from modules/line_descriptor/include/opencv2/line_descriptor/bitops.hpp rename to modules/line_descriptor/src/bitops.hpp index 80ed7dab9..4187b7528 100644 --- a/modules/line_descriptor/include/opencv2/line_descriptor/bitops.hpp +++ b/modules/line_descriptor/src/bitops.hpp @@ -43,6 +43,8 @@ #ifndef __OPENCV_BITOPTS_HPP #define __OPENCV_BITOPTS_HPP +#include "precomp.hpp" + #ifdef _WIN32 # include # define popcnt __popcnt @@ -63,6 +65,10 @@ const int lookup[] = 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8 }; +namespace cv +{ +namespace line_descriptor +{ /*matching function */ inline int match( UINT8*P, UINT8*Q, int codelb ) { @@ -163,5 +169,7 @@ inline UINT64 choose( int n, int r ) return nchooser; } +} +} #endif diff --git a/modules/line_descriptor/src/bucket_group.cpp b/modules/line_descriptor/src/bucket_group.cpp deleted file mode 100644 index ed44b0ea4..000000000 --- a/modules/line_descriptor/src/bucket_group.cpp +++ /dev/null @@ -1,100 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#include "precomp.hpp" - -/* constructor */ -BucketGroup::BucketGroup() -{ - empty = 0; - group = NULL; -} - -/* destructor */ -BucketGroup::~BucketGroup() -{ - if( group != NULL ) - delete group; -} - -/* insert data into the bucket */ -void BucketGroup::insert( int subindex, UINT32 data ) -{ - if( group == NULL ) - { - group = new Array32(); - group->push( 0 ); - } - - UINT32 lowerbits = ( (UINT32) 1 << subindex ) - 1; - int end = popcnt( empty & lowerbits ); - - if( ! ( empty & ( (UINT32) 1 << subindex ) ) ) - { - group->insert( end, group->arr[end + 2] ); - empty |= (UINT32) 1 << subindex; - } - - int totones = popcnt( empty ); - group->insert( totones + 1 + group->arr[2 + end + 1], data ); - for ( int i = end + 1; i < totones + 1; i++ ) - group->arr[2 + i]++; -} - -/* perform a query to the bucket */ -UINT32* BucketGroup::query( int subindex, int *size ) -{ - if( empty & ( (UINT32) 1 << subindex ) ) - { - UINT32 lowerbits = ( (UINT32) 1 << subindex ) - 1; - int end = popcnt( empty & lowerbits ); - int totones = popcnt( empty ); - *size = group->arr[2 + end + 1] - group->arr[2 + end]; - return group->arr + 2 + totones + 1 + group->arr[2 + end]; - } - - else - { - *size = 0; - return NULL; - } -} diff --git a/modules/line_descriptor/src/ed_line_detector.cpp b/modules/line_descriptor/src/ed_line_detector.cpp deleted file mode 100644 index f1fbe81c6..000000000 --- a/modules/line_descriptor/src/ed_line_detector.cpp +++ /dev/null @@ -1,1432 +0,0 @@ -/*IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - - By downloading, copying, installing or using the software you agree to this license. - If you do not agree to this license, do not download, install, - copy or use the software. - - - License Agreement - For Open Source Computer Vision Library - - Copyright (C) 2011-2012, Lilian Zhang, all rights reserved. - Copyright (C) 2013, Manuele Tamburrano, Stefano Fabri, all rights reserved. - Third party copyrights are property of their respective owners. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * The name of the copyright holders may not be used to endorse or promote products - derived from this software without specific prior written permission. - - This software is provided by the copyright holders and contributors "as is" and - any express or implied warranties, including, but not limited to, the implied - warranties of merchantability and fitness for a particular purpose are disclaimed. - In no event shall the Intel Corporation or contributors be liable for any direct, - indirect, incidental, special, exemplary, or consequential damages - (including, but not limited to, procurement of substitute goods or services; - loss of use, data, or profits; or business interruption) however caused - and on any theory of liability, whether in contract, strict liability, - or tort (including negligence or otherwise) arising in any way out of - the use of this software, even if advised of the possibility of such damage. - */ - -#include "precomp.hpp" - -#define Horizontal 255//if |dx|<|dy|; -#define Vertical 0//if |dy|<=|dx|; -#define UpDir 1 -#define RightDir 2 -#define DownDir 3 -#define LeftDir 4 -#define TryTime 6 -#define SkipEdgePoint 2 - -using namespace std; -EDLineDetector::EDLineDetector() -{ - //set parameters for line segment detection - ksize_ = 15; //15 - sigma_ = 30.0; //30 - gradienThreshold_ = 80; // ***** ORIGINAL WAS 25 - anchorThreshold_ = 8; //8 - scanIntervals_ = 2; //2 - minLineLen_ = 15; //15 - lineFitErrThreshold_ = 1.6; //1.4 - InitEDLine_(); -} -EDLineDetector::EDLineDetector( EDLineParam param ) -{ - //set parameters for line segment detection - ksize_ = param.ksize; - sigma_ = param.sigma; - gradienThreshold_ = (short) param.gradientThreshold; - anchorThreshold_ = (unsigned char) param.anchorThreshold; - scanIntervals_ = param.scanIntervals; - minLineLen_ = param.minLineLen; - lineFitErrThreshold_ = param.lineFitErrThreshold; - InitEDLine_(); -} -void EDLineDetector::InitEDLine_() -{ - bValidate_ = true; - ATA = cv::Mat_( 2, 2 ); - ATV = cv::Mat_( 1, 2 ); - tempMatLineFit = cv::Mat_( 2, 2 ); - tempVecLineFit = cv::Mat_( 1, 2 ); - fitMatT = cv::Mat_( 2, minLineLen_ ); - fitVec = cv::Mat_( 1, minLineLen_ ); - for ( int i = 0; i < minLineLen_; i++ ) - { - fitMatT[1][i] = 1; - } - dxImg_.create( 1, 1, CV_16SC1 ); - dyImg_.create( 1, 1, CV_16SC1 ); - gImgWO_.create( 1, 1, CV_8SC1 ); - pFirstPartEdgeX_ = NULL; - pFirstPartEdgeY_ = NULL; - pFirstPartEdgeS_ = NULL; - pSecondPartEdgeX_ = NULL; - pSecondPartEdgeY_ = NULL; - pSecondPartEdgeS_ = NULL; - pAnchorX_ = NULL; - pAnchorY_ = NULL; -} - -EDLineDetector::~EDLineDetector() -{ - if( pFirstPartEdgeX_ != NULL ) - { - delete[] pFirstPartEdgeX_; - delete[] pFirstPartEdgeY_; - delete[] pSecondPartEdgeX_; - delete[] pSecondPartEdgeY_; - delete[] pAnchorX_; - delete[] pAnchorY_; - } - if( pFirstPartEdgeS_ != NULL ) - { - delete[] pFirstPartEdgeS_; - delete[] pSecondPartEdgeS_; - } -} - -int EDLineDetector::EdgeDrawing( cv::Mat &image, EdgeChains &edgeChains ) -{ - imageWidth = image.cols; - imageHeight = image.rows; - unsigned int pixelNum = imageWidth * imageHeight; - - unsigned int edgePixelArraySize = pixelNum / 5; - unsigned int maxNumOfEdge = edgePixelArraySize / 20; - //compute dx, dy images - if( gImg_.cols != (int) imageWidth || gImg_.rows != (int) imageHeight ) - { - if( pFirstPartEdgeX_ != NULL ) - { - delete[] pFirstPartEdgeX_; - delete[] pFirstPartEdgeY_; - delete[] pSecondPartEdgeX_; - delete[] pSecondPartEdgeY_; - delete[] pFirstPartEdgeS_; - delete[] pSecondPartEdgeS_; - delete[] pAnchorX_; - delete[] pAnchorY_; - } - - dxImg_.create( imageHeight, imageWidth, CV_16SC1 ); - dyImg_.create( imageHeight, imageWidth, CV_16SC1 ); - gImgWO_.create( imageHeight, imageWidth, CV_16SC1 ); - gImg_.create( imageHeight, imageWidth, CV_16SC1 ); - dirImg_.create( imageHeight, imageWidth, CV_8UC1 ); - edgeImage_.create( imageHeight, imageWidth, CV_8UC1 ); - pFirstPartEdgeX_ = new unsigned int[edgePixelArraySize]; - pFirstPartEdgeY_ = new unsigned int[edgePixelArraySize]; - pSecondPartEdgeX_ = new unsigned int[edgePixelArraySize]; - pSecondPartEdgeY_ = new unsigned int[edgePixelArraySize]; - pFirstPartEdgeS_ = new unsigned int[maxNumOfEdge]; - pSecondPartEdgeS_ = new unsigned int[maxNumOfEdge]; - pAnchorX_ = new unsigned int[edgePixelArraySize]; - pAnchorY_ = new unsigned int[edgePixelArraySize]; - } - cv::Sobel( image, dxImg_, CV_16SC1, 1, 0, 3 ); - cv::Sobel( image, dyImg_, CV_16SC1, 0, 1, 3 ); - - //compute gradient and direction images - cv::Mat dxABS_m = cv::abs( dxImg_ ); - cv::Mat dyABS_m = cv::abs( dyImg_ ); - cv::Mat sumDxDy; - cv::add( dyABS_m, dxABS_m, sumDxDy ); - - cv::threshold( sumDxDy, gImg_, gradienThreshold_ + 1, 255, cv::THRESH_TOZERO ); - gImg_ = gImg_ / 4; - gImgWO_ = sumDxDy / 4; - cv::compare( dxABS_m, dyABS_m, dirImg_, cv::CMP_LT ); - - short *pgImg = gImg_.ptr(); - unsigned char *pdirImg = dirImg_.ptr(); - - //extract the anchors in the gradient image, store into a vector - memset( pAnchorX_, 0, edgePixelArraySize * sizeof(unsigned int) ); //initialization - memset( pAnchorY_, 0, edgePixelArraySize * sizeof(unsigned int) ); - unsigned int anchorsSize = 0; - int indexInArray; - unsigned char gValue1, gValue2, gValue3; - for ( unsigned int w = 1; w < imageWidth - 1; w = w + scanIntervals_ ) - { - for ( unsigned int h = 1; h < imageHeight - 1; h = h + scanIntervals_ ) - { - indexInArray = h * imageWidth + w; - //gValue1 = pdirImg[indexInArray]; - if( pdirImg[indexInArray] == Horizontal ) - { //if the direction of pixel is horizontal, then compare with up and down - //gValue2 = pgImg[indexInArray]; - if( pgImg[indexInArray] >= pgImg[indexInArray - imageWidth] + anchorThreshold_ - && pgImg[indexInArray] >= pgImg[indexInArray + imageWidth] + anchorThreshold_ ) - { // (w,h) is accepted as an anchor - pAnchorX_[anchorsSize] = w; - pAnchorY_[anchorsSize++] = h; - } - } - else - { // if(pdirImg[indexInArray]==Vertical){//it is vertical edge, should be compared with left and right - //gValue2 = pgImg[indexInArray]; - if( pgImg[indexInArray] >= pgImg[indexInArray - 1] + anchorThreshold_ && pgImg[indexInArray] >= pgImg[indexInArray + 1] + anchorThreshold_ ) - { // (w,h) is accepted as an anchor - pAnchorX_[anchorsSize] = w; - pAnchorY_[anchorsSize++] = h; - } - } - } - } - if( anchorsSize > edgePixelArraySize ) - { - cout << "anchor size is larger than its maximal size. anchorsSize=" << anchorsSize << ", maximal size = " << edgePixelArraySize << endl; - return -1; - } - - //link the anchors by smart routing - edgeImage_.setTo( 0 ); - unsigned char *pEdgeImg = edgeImage_.data; - memset( pFirstPartEdgeX_, 0, edgePixelArraySize * sizeof(unsigned int) ); //initialization - memset( pFirstPartEdgeY_, 0, edgePixelArraySize * sizeof(unsigned int) ); - memset( pSecondPartEdgeX_, 0, edgePixelArraySize * sizeof(unsigned int) ); - memset( pSecondPartEdgeY_, 0, edgePixelArraySize * sizeof(unsigned int) ); - memset( pFirstPartEdgeS_, 0, maxNumOfEdge * sizeof(unsigned int) ); - memset( pSecondPartEdgeS_, 0, maxNumOfEdge * sizeof(unsigned int) ); - unsigned int offsetPFirst = 0, offsetPSecond = 0; - unsigned int offsetPS = 0; - - unsigned int x, y; - unsigned int lastX = 0; - unsigned int lastY = 0; - unsigned char lastDirection; //up = 1, right = 2, down = 3, left = 4; - unsigned char shouldGoDirection; //up = 1, right = 2, down = 3, left = 4; - int edgeLenFirst, edgeLenSecond; - for ( unsigned int i = 0; i < anchorsSize; i++ ) - { - x = pAnchorX_[i]; - y = pAnchorY_[i]; - indexInArray = y * imageWidth + x; - if( pEdgeImg[indexInArray] ) - { //if anchor i is already been an edge pixel. - continue; - } - /*The walk stops under 3 conditions: - * 1. We move out of the edge areas, i.e., the thresholded gradient value - * of the current pixel is 0. - * 2. The current direction of the edge changes, i.e., from horizontal - * to vertical or vice versa.?? (This is turned out not correct. From the online edge draw demo - * we can figure out that authors don't implement this rule either because their extracted edge - * chain could be a circle which means pixel directions would definitely be different - * in somewhere on the chain.) - * 3. We encounter a previously detected edge pixel. */ - pFirstPartEdgeS_[offsetPS] = offsetPFirst; - if( pdirImg[indexInArray] == Horizontal ) - { //if the direction of this pixel is horizontal, then go left and right. - //fist go right, pixel direction may be different during linking. - lastDirection = RightDir; - while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) - { - pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel - pFirstPartEdgeX_[offsetPFirst] = x; - pFirstPartEdgeY_[offsetPFirst++] = y; - shouldGoDirection = 0; //unknown - if( pdirImg[indexInArray] == Horizontal ) - { //should go left or right - if( lastDirection == UpDir || lastDirection == DownDir ) - { //change the pixel direction now - if( x > lastX ) - { //should go right - shouldGoDirection = RightDir; - } - else - { //should go left - shouldGoDirection = LeftDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == RightDir || shouldGoDirection == RightDir ) - { //go right - if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the right and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-right - x = x + 1; - y = y + 1; - } - else - { //straight-right - x = x + 1; - } - lastDirection = RightDir; - } - else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) - { //go left - if( x == 0 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the left and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - gValue2 = (unsigned char) pgImg[indexInArray - 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-left - x = x - 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-left - x = x - 1; - } - lastDirection = LeftDir; - } - } - else - { //should go up or down. - if( lastDirection == RightDir || lastDirection == LeftDir ) - { //change the pixel direction now - if( y > lastY ) - { //should go down - shouldGoDirection = DownDir; - } - else - { //should go up - shouldGoDirection = UpDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == DownDir || shouldGoDirection == DownDir ) - { //go down - if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the down and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //down-right - x = x + 1; - y = y + 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-down - y = y + 1; - } - lastDirection = DownDir; - } - else if( lastDirection == UpDir || shouldGoDirection == UpDir ) - { //go up - if( x == 0 || x == imageWidth - 1 || y == 0 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the up and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //up-left - x = x - 1; - y = y - 1; - } - else - { //straight-up - y = y - 1; - } - lastDirection = UpDir; - } - } - indexInArray = y * imageWidth + x; - } //end while go right - //then go left, pixel direction may be different during linking. - x = pAnchorX_[i]; - y = pAnchorY_[i]; - indexInArray = y * imageWidth + x; - pEdgeImg[indexInArray] = 0; //mark the anchor point be a non-edge pixel and - lastDirection = LeftDir; - pSecondPartEdgeS_[offsetPS] = offsetPSecond; - while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) - { - pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel - pSecondPartEdgeX_[offsetPSecond] = x; - pSecondPartEdgeY_[offsetPSecond++] = y; - shouldGoDirection = 0; //unknown - if( pdirImg[indexInArray] == Horizontal ) - { //should go left or right - if( lastDirection == UpDir || lastDirection == DownDir ) - { //change the pixel direction now - if( x > lastX ) - { //should go right - shouldGoDirection = RightDir; - } - else - { //should go left - shouldGoDirection = LeftDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == RightDir || shouldGoDirection == RightDir ) - { //go right - if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the right and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-right - x = x + 1; - y = y + 1; - } - else - { //straight-right - x = x + 1; - } - lastDirection = RightDir; - } - else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) - { //go left - if( x == 0 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the left and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - gValue2 = (unsigned char) pgImg[indexInArray - 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-left - x = x - 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-left - x = x - 1; - } - lastDirection = LeftDir; - } - } - else - { //should go up or down. - if( lastDirection == RightDir || lastDirection == LeftDir ) - { //change the pixel direction now - if( y > lastY ) - { //should go down - shouldGoDirection = DownDir; - } - else - { //should go up - shouldGoDirection = UpDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == DownDir || shouldGoDirection == DownDir ) - { //go down - if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the down and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //down-right - x = x + 1; - y = y + 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-down - y = y + 1; - } - lastDirection = DownDir; - } - else if( lastDirection == UpDir || shouldGoDirection == UpDir ) - { //go up - if( x == 0 || x == imageWidth - 1 || y == 0 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the up and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //up-left - x = x - 1; - y = y - 1; - } - else - { //straight-up - y = y - 1; - } - lastDirection = UpDir; - } - } - indexInArray = y * imageWidth + x; - } //end while go left - //end anchor is Horizontal - } - else - { //the direction of this pixel is vertical, go up and down - //fist go down, pixel direction may be different during linking. - lastDirection = DownDir; - while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) - { - pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel - pFirstPartEdgeX_[offsetPFirst] = x; - pFirstPartEdgeY_[offsetPFirst++] = y; - shouldGoDirection = 0; //unknown - if( pdirImg[indexInArray] == Horizontal ) - { //should go left or right - if( lastDirection == UpDir || lastDirection == DownDir ) - { //change the pixel direction now - if( x > lastX ) - { //should go right - shouldGoDirection = RightDir; - } - else - { //should go left - shouldGoDirection = LeftDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == RightDir || shouldGoDirection == RightDir ) - { //go right - if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the right and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-right - x = x + 1; - y = y + 1; - } - else - { //straight-right - x = x + 1; - } - lastDirection = RightDir; - } - else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) - { //go left - if( x == 0 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the left and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - gValue2 = (unsigned char) pgImg[indexInArray - 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-left - x = x - 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-left - x = x - 1; - } - lastDirection = LeftDir; - } - } - else - { //should go up or down. - if( lastDirection == RightDir || lastDirection == LeftDir ) - { //change the pixel direction now - if( y > lastY ) - { //should go down - shouldGoDirection = DownDir; - } - else - { //should go up - shouldGoDirection = UpDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == DownDir || shouldGoDirection == DownDir ) - { //go down - if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the down and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //down-right - x = x + 1; - y = y + 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-down - y = y + 1; - } - lastDirection = DownDir; - } - else if( lastDirection == UpDir || shouldGoDirection == UpDir ) - { //go up - if( x == 0 || x == imageWidth - 1 || y == 0 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the up and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //up-left - x = x - 1; - y = y - 1; - } - else - { //straight-up - y = y - 1; - } - lastDirection = UpDir; - } - } - indexInArray = y * imageWidth + x; - } //end while go down - //then go up, pixel direction may be different during linking. - lastDirection = UpDir; - x = pAnchorX_[i]; - y = pAnchorY_[i]; - indexInArray = y * imageWidth + x; - pEdgeImg[indexInArray] = 0; //mark the anchor point be a non-edge pixel and - pSecondPartEdgeS_[offsetPS] = offsetPSecond; - while ( pgImg[indexInArray] > 0 && !pEdgeImg[indexInArray] ) - { - pEdgeImg[indexInArray] = 1; // Mark this pixel as an edge pixel - pSecondPartEdgeX_[offsetPSecond] = x; - pSecondPartEdgeY_[offsetPSecond++] = y; - shouldGoDirection = 0; //unknown - if( pdirImg[indexInArray] == Horizontal ) - { //should go left or right - if( lastDirection == UpDir || lastDirection == DownDir ) - { //change the pixel direction now - if( x > lastX ) - { //should go right - shouldGoDirection = RightDir; - } - else - { //should go left - shouldGoDirection = LeftDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == RightDir || shouldGoDirection == RightDir ) - { //go right - if( x == imageWidth - 1 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the right and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-right - x = x + 1; - y = y + 1; - } - else - { //straight-right - x = x + 1; - } - lastDirection = RightDir; - } - else if( lastDirection == LeftDir || shouldGoDirection == LeftDir ) - { //go left - if( x == 0 || y == 0 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the left and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - gValue2 = (unsigned char) pgImg[indexInArray - 1]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-left - x = x - 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-left - x = x - 1; - } - lastDirection = LeftDir; - } - } - else - { //should go up or down. - if( lastDirection == RightDir || lastDirection == LeftDir ) - { //change the pixel direction now - if( y > lastY ) - { //should go down - shouldGoDirection = DownDir; - } - else - { //should go up - shouldGoDirection = UpDir; - } - } - lastX = x; - lastY = y; - if( lastDirection == DownDir || shouldGoDirection == DownDir ) - { //go down - if( x == 0 || x == imageWidth - 1 || y == imageHeight - 1 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the down and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray + imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray + imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray + imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //down-right - x = x + 1; - y = y + 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //down-left - x = x - 1; - y = y + 1; - } - else - { //straight-down - y = y + 1; - } - lastDirection = DownDir; - } - else if( lastDirection == UpDir || shouldGoDirection == UpDir ) - { //go up - if( x == 0 || x == imageWidth - 1 || y == 0 ) - { //reach the image border - break; - } - // Look at 3 neighbors to the up and pick the one with the max. gradient value - gValue1 = (unsigned char) pgImg[indexInArray - imageWidth + 1]; - gValue2 = (unsigned char) pgImg[indexInArray - imageWidth]; - gValue3 = (unsigned char) pgImg[indexInArray - imageWidth - 1]; - if( gValue1 >= gValue2 && gValue1 >= gValue3 ) - { //up-right - x = x + 1; - y = y - 1; - } - else if( gValue3 >= gValue2 && gValue3 >= gValue1 ) - { //up-left - x = x - 1; - y = y - 1; - } - else - { //straight-up - y = y - 1; - } - lastDirection = UpDir; - } - } - indexInArray = y * imageWidth + x; - } //end while go up - } //end anchor is Vertical - //only keep the edge chains whose length is larger than the minLineLen_; - edgeLenFirst = offsetPFirst - pFirstPartEdgeS_[offsetPS]; - edgeLenSecond = offsetPSecond - pSecondPartEdgeS_[offsetPS]; - if( edgeLenFirst + edgeLenSecond < minLineLen_ + 1 ) - { //short edge, drop it - offsetPFirst = pFirstPartEdgeS_[offsetPS]; - offsetPSecond = pSecondPartEdgeS_[offsetPS]; - } - else - { - offsetPS++; - } - } - //store the last index - pFirstPartEdgeS_[offsetPS] = offsetPFirst; - pSecondPartEdgeS_[offsetPS] = offsetPSecond; - if( offsetPS > maxNumOfEdge ) - { - cout << "Edge drawing Error: The total number of edges is larger than MaxNumOfEdge, " - "numofedge = " - << offsetPS << ", MaxNumOfEdge=" << maxNumOfEdge << endl; - return -1; - } - if( offsetPFirst > edgePixelArraySize || offsetPSecond > edgePixelArraySize ) - { - cout << "Edge drawing Error: The total number of edge pixels is larger than MaxNumOfEdgePixels, " - "numofedgePixel1 = " - << offsetPFirst << ", numofedgePixel2 = " << offsetPSecond << ", MaxNumOfEdgePixel=" << edgePixelArraySize << endl; - return -1; - } - - /*now all the edge information are stored in pFirstPartEdgeX_, pFirstPartEdgeY_, - *pFirstPartEdgeS_, pSecondPartEdgeX_, pSecondPartEdgeY_, pSecondPartEdgeS_; - *we should reorganize them into edgeChains for easily using. */ - int tempID; - edgeChains.xCors.resize( offsetPFirst + offsetPSecond ); - edgeChains.yCors.resize( offsetPFirst + offsetPSecond ); - edgeChains.sId.resize( offsetPS + 1 ); - unsigned int *pxCors = edgeChains.xCors.data(); - unsigned int *pyCors = edgeChains.yCors.data(); - unsigned int *psId = edgeChains.sId.data(); - offsetPFirst = 0; - offsetPSecond = 0; - unsigned int indexInCors = 0; - unsigned int numOfEdges = 0; - for ( unsigned int edgeId = 0; edgeId < offsetPS; edgeId++ ) - { - //step1, put the first and second parts edge coordinates together from edge start to edge end - psId[numOfEdges++] = indexInCors; - indexInArray = pFirstPartEdgeS_[edgeId]; - offsetPFirst = pFirstPartEdgeS_[edgeId + 1]; - for ( tempID = offsetPFirst - 1; tempID >= indexInArray; tempID-- ) - { //add first part edge - pxCors[indexInCors] = pFirstPartEdgeX_[tempID]; - pyCors[indexInCors++] = pFirstPartEdgeY_[tempID]; - } - indexInArray = pSecondPartEdgeS_[edgeId]; - offsetPSecond = pSecondPartEdgeS_[edgeId + 1]; - for ( tempID = indexInArray + 1; tempID < (int) offsetPSecond; tempID++ ) - { //add second part edge - pxCors[indexInCors] = pSecondPartEdgeX_[tempID]; - pyCors[indexInCors++] = pSecondPartEdgeY_[tempID]; - } - } - psId[numOfEdges] = indexInCors; //the end index of the last edge - edgeChains.numOfEdges = numOfEdges; - - return 1; -} - -int EDLineDetector::EDline( cv::Mat &image, LineChains &lines ) -{ - - //first, call EdgeDrawing function to extract edges - EdgeChains edges; - if( ( EdgeDrawing( image, edges ) ) != 1 ) - { - cout << "Line Detection not finished" << endl; - return -1; - } - - //detect lines - unsigned int linePixelID = edges.sId[edges.numOfEdges]; - lines.xCors.resize( linePixelID ); - lines.yCors.resize( linePixelID ); - lines.sId.resize( 5 * edges.numOfEdges ); - unsigned int *pEdgeXCors = edges.xCors.data(); - unsigned int *pEdgeYCors = edges.yCors.data(); - unsigned int *pEdgeSID = edges.sId.data(); - unsigned int *pLineXCors = lines.xCors.data(); - unsigned int *pLineYCors = lines.yCors.data(); - unsigned int *pLineSID = lines.sId.data(); - logNT_ = 2.0 * ( log10( (double) imageWidth ) + log10( (double) imageHeight ) ); - double lineFitErr = 0; //the line fit error; - std::vector lineEquation( 2, 0 ); - lineEquations_.clear(); - lineEndpoints_.clear(); - lineDirection_.clear(); - unsigned char *pdirImg = dirImg_.data; - unsigned int numOfLines = 0; - unsigned int newOffsetS = 0; - unsigned int offsetInEdgeArrayS, offsetInEdgeArrayE; //start index and end index - unsigned int offsetInLineArray = 0; - float direction; //line direction - - for ( unsigned int edgeID = 0; edgeID < edges.numOfEdges; edgeID++ ) - { - offsetInEdgeArrayS = pEdgeSID[edgeID]; - offsetInEdgeArrayE = pEdgeSID[edgeID + 1]; - while ( offsetInEdgeArrayE > offsetInEdgeArrayS + minLineLen_ ) - { //extract line segments from an edge, may find more than one segments - //find an initial line segment - while ( offsetInEdgeArrayE > offsetInEdgeArrayS + minLineLen_ ) - { - lineFitErr = LeastSquaresLineFit_( pEdgeXCors, pEdgeYCors, offsetInEdgeArrayS, lineEquation ); - if( lineFitErr <= lineFitErrThreshold_ ) - break; //ok, an initial line segment detected - offsetInEdgeArrayS += SkipEdgePoint; //skip the first two pixel in the chain and try with the remaining pixels - } - if( lineFitErr > lineFitErrThreshold_ ) - break; //no line is detected - //An initial line segment is detected. Try to extend this line segment - pLineSID[numOfLines] = offsetInLineArray; - double coef1 = 0; //for a line ax+by+c=0, coef1 = 1/sqrt(a^2+b^2); - double pointToLineDis; //for a line ax+by+c=0 and a point(xi, yi), pointToLineDis = coef1*|a*xi+b*yi+c| - bool bExtended = true; - bool bFirstTry = true; - int numOfOutlier; //to against noise, we accept a few outlier of a line. - int tryTimes = 0; - if( pdirImg[pEdgeYCors[offsetInEdgeArrayS] * imageWidth + pEdgeXCors[offsetInEdgeArrayS]] == Horizontal ) - { //y=ax+b, i.e. ax-y+b=0 - while ( bExtended ) - { - tryTimes++; - if( bFirstTry ) - { - bFirstTry = false; - for ( int i = 0; i < minLineLen_; i++ ) - { //First add the initial line segment to the line array - pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; - pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; - } - } - else - { //after each try, line is extended, line equation should be re-estimated - //adjust the line equation - lineFitErr = LeastSquaresLineFit_( pLineXCors, pLineYCors, pLineSID[numOfLines], newOffsetS, offsetInLineArray, lineEquation ); - } - coef1 = 1 / sqrt( lineEquation[0] * lineEquation[0] + 1 ); - numOfOutlier = 0; - newOffsetS = offsetInLineArray; - while ( offsetInEdgeArrayE > offsetInEdgeArrayS ) - { - pointToLineDis = fabs( lineEquation[0] * pEdgeXCors[offsetInEdgeArrayS] - pEdgeYCors[offsetInEdgeArrayS] + lineEquation[1] ) * coef1; - pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; - pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; - if( pointToLineDis > lineFitErrThreshold_ ) - { - numOfOutlier++; - if( numOfOutlier > 3 ) - break; - } - else - { //we count number of connective outliers. - numOfOutlier = 0; - } - } - //pop back the last few outliers from lines and return them to edge chain - offsetInLineArray -= numOfOutlier; - offsetInEdgeArrayS -= numOfOutlier; - if( offsetInLineArray - newOffsetS > 0 && tryTimes < TryTime ) - { //some new pixels are added to the line - } - else - { - bExtended = false; //no new pixels are added. - } - } - //the line equation coefficients,for line w1x+w2y+w3 =0, we normalize it to make w1^2+w2^2 = 1. - std::vector lineEqu( 3, 0 ); - lineEqu[0] = lineEquation[0] * coef1; - lineEqu[1] = -1 * coef1; - lineEqu[2] = lineEquation[1] * coef1; - if( LineValidation_( pLineXCors, pLineYCors, pLineSID[numOfLines], offsetInLineArray, lineEqu, direction ) ) - { //check the line - //store the line equation coefficients - lineEquations_.push_back( lineEqu ); - /*At last, compute the line endpoints and store them. - *we project the first and last pixels in the pixelChain onto the best fit line - *to get the line endpoints. - *xp= (w2^2*x0-w1*w2*y0-w3*w1)/(w1^2+w2^2) - *yp= (w1^2*y0-w1*w2*x0-w3*w2)/(w1^2+w2^2) */ - std::vector lineEndP( 4, 0 ); //line endpoints - double a1 = lineEqu[1] * lineEqu[1]; - double a2 = lineEqu[0] * lineEqu[0]; - double a3 = lineEqu[0] * lineEqu[1]; - double a4 = lineEqu[2] * lineEqu[0]; - double a5 = lineEqu[2] * lineEqu[1]; - unsigned int Px = pLineXCors[pLineSID[numOfLines]]; //first pixel - unsigned int Py = pLineYCors[pLineSID[numOfLines]]; - lineEndP[0] = (float) ( a1 * Px - a3 * Py - a4 ); //x - lineEndP[1] = (float) ( a2 * Py - a3 * Px - a5 ); //y - Px = pLineXCors[offsetInLineArray - 1]; //last pixel - Py = pLineYCors[offsetInLineArray - 1]; - lineEndP[2] = (float) ( a1 * Px - a3 * Py - a4 ); //x - lineEndP[3] = (float) ( a2 * Py - a3 * Px - a5 ); //y - lineEndpoints_.push_back( lineEndP ); - lineDirection_.push_back( direction ); - numOfLines++; - } - else - { - offsetInLineArray = pLineSID[numOfLines]; // line was not accepted, the offset is set back - } - } - else - { //x=ay+b, i.e. x-ay-b=0 - while ( bExtended ) - { - tryTimes++; - if( bFirstTry ) - { - bFirstTry = false; - for ( int i = 0; i < minLineLen_; i++ ) - { //First add the initial line segment to the line array - pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; - pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; - } - } - else - { //after each try, line is extended, line equation should be re-estimated - //adjust the line equation - lineFitErr = LeastSquaresLineFit_( pLineXCors, pLineYCors, pLineSID[numOfLines], newOffsetS, offsetInLineArray, lineEquation ); - } - coef1 = 1 / sqrt( 1 + lineEquation[0] * lineEquation[0] ); - numOfOutlier = 0; - newOffsetS = offsetInLineArray; - while ( offsetInEdgeArrayE > offsetInEdgeArrayS ) - { - pointToLineDis = fabs( pEdgeXCors[offsetInEdgeArrayS] - lineEquation[0] * pEdgeYCors[offsetInEdgeArrayS] - lineEquation[1] ) * coef1; - pLineXCors[offsetInLineArray] = pEdgeXCors[offsetInEdgeArrayS]; - pLineYCors[offsetInLineArray++] = pEdgeYCors[offsetInEdgeArrayS++]; - if( pointToLineDis > lineFitErrThreshold_ ) - { - numOfOutlier++; - if( numOfOutlier > 3 ) - break; - } - else - { //we count number of connective outliers. - numOfOutlier = 0; - } - } - //pop back the last few outliers from lines and return them to edge chain - offsetInLineArray -= numOfOutlier; - offsetInEdgeArrayS -= numOfOutlier; - if( offsetInLineArray - newOffsetS > 0 && tryTimes < TryTime ) - { //some new pixels are added to the line - } - else - { - bExtended = false; //no new pixels are added. - } - } - //the line equation coefficients,for line w1x+w2y+w3 =0, we normalize it to make w1^2+w2^2 = 1. - std::vector lineEqu( 3, 0 ); - lineEqu[0] = 1 * coef1; - lineEqu[1] = -lineEquation[0] * coef1; - lineEqu[2] = -lineEquation[1] * coef1; - - if( LineValidation_( pLineXCors, pLineYCors, pLineSID[numOfLines], offsetInLineArray, lineEqu, direction ) ) - { //check the line - //store the line equation coefficients - lineEquations_.push_back( lineEqu ); - /*At last, compute the line endpoints and store them. - *we project the first and last pixels in the pixelChain onto the best fit line - *to get the line endpoints. - *xp= (w2^2*x0-w1*w2*y0-w3*w1)/(w1^2+w2^2) - *yp= (w1^2*y0-w1*w2*x0-w3*w2)/(w1^2+w2^2) */ - std::vector lineEndP( 4, 0 ); //line endpoints - double a1 = lineEqu[1] * lineEqu[1]; - double a2 = lineEqu[0] * lineEqu[0]; - double a3 = lineEqu[0] * lineEqu[1]; - double a4 = lineEqu[2] * lineEqu[0]; - double a5 = lineEqu[2] * lineEqu[1]; - unsigned int Px = pLineXCors[pLineSID[numOfLines]]; //first pixel - unsigned int Py = pLineYCors[pLineSID[numOfLines]]; - lineEndP[0] = (float) ( a1 * Px - a3 * Py - a4 ); //x - lineEndP[1] = (float) ( a2 * Py - a3 * Px - a5 ); //y - Px = pLineXCors[offsetInLineArray - 1]; //last pixel - Py = pLineYCors[offsetInLineArray - 1]; - lineEndP[2] = (float) ( a1 * Px - a3 * Py - a4 ); //x - lineEndP[3] = (float) ( a2 * Py - a3 * Px - a5 ); //y - lineEndpoints_.push_back( lineEndP ); - lineDirection_.push_back( direction ); - numOfLines++; - } - else - { - offsetInLineArray = pLineSID[numOfLines]; // line was not accepted, the offset is set back - } - } - //Extract line segments from the remaining pixel; Current chain has been shortened already. - } - } //end for(unsigned int edgeID=0; edgeID &lineEquation ) -{ - - float * pMatT; - float * pATA; - double fitError = 0; - double coef; - unsigned char *pdirImg = dirImg_.data; - unsigned int offset = offsetS; - /*If the first pixel in this chain is horizontal, - *then we try to find a horizontal line, y=ax+b;*/ - if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Horizontal ) - { - /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec - * [x0,1] [y0] - * [x1,1] [a] [y1] - * . [b] = . - * [xn,1] [yn]*/ - pMatT = fitMatT.ptr(); //fitMatT = [x0, x1, ... xn; 1,1,...,1]; - for ( int i = 0; i < minLineLen_; i++ ) - { - //*(pMatT+minLineLen_) = 1; //the value are not changed; - * ( pMatT++ ) = (float) xCors[offsetS]; - fitVec[0][i] = (float) yCors[offsetS++]; - } - ATA = fitMatT * fitMatT.t(); - ATV = fitMatT * fitVec.t(); - /* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */ - pATA = ATA.ptr(); - coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); - // lineEquation = svd.Invert(ATA) * matT * vec; - lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); - lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); - /*compute line fit error */ - for ( int i = 0; i < minLineLen_; i++ ) - { - //coef = double(yCors[offset]) - double(xCors[offset++]) * lineEquation[0] - lineEquation[1]; - coef = double( yCors[offset] ) - double( xCors[offset] ) * lineEquation[0] - lineEquation[1]; - offset++; - fitError += coef * coef; - } - return sqrt( fitError ); - } - /*If the first pixel in this chain is vertical, - *then we try to find a vertical line, x=ay+b;*/ - if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Vertical ) - { - /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec - * [y0,1] [x0] - * [y1,1] [a] [x1] - * . [b] = . - * [yn,1] [xn]*/ - pMatT = fitMatT.ptr(); //fitMatT = [y0, y1, ... yn; 1,1,...,1]; - for ( int i = 0; i < minLineLen_; i++ ) - { - //*(pMatT+minLineLen_) = 1;//the value are not changed; - * ( pMatT++ ) = (float) yCors[offsetS]; - fitVec[0][i] = (float) xCors[offsetS++]; - } - ATA = fitMatT * ( fitMatT.t() ); - ATV = fitMatT * fitVec.t(); - /* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */ - pATA = ATA.ptr(); - coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); - // lineEquation = svd.Invert(ATA) * matT * vec; - lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); - lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); - /*compute line fit error */ - for ( int i = 0; i < minLineLen_; i++ ) - { - //coef = double(xCors[offset]) - double(yCors[offset++]) * lineEquation[0] - lineEquation[1]; - coef = double( xCors[offset] ) - double( yCors[offset] ) * lineEquation[0] - lineEquation[1]; - offset++; - fitError += coef * coef; - } - return sqrt( fitError ); - } - return 0; -} -double EDLineDetector::LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int newOffsetS, - unsigned int offsetE, std::vector &lineEquation ) -{ - int length = offsetE - offsetS; - int newLength = offsetE - newOffsetS; - if( length <= 0 || newLength <= 0 ) - { - cout << "EDLineDetector::LeastSquaresLineFit_ Error:" - " the expected line index is wrong...offsetE = " - << offsetE << ", offsetS=" << offsetS << ", newOffsetS=" << newOffsetS << endl; - return -1; - } - if( lineEquation.size() != 2 ) - { - std::cout << "SHOULD NOT BE != 2" << std::endl; - } - cv::Mat_ matT( 2, newLength ); - cv::Mat_ vec( newLength, 1 ); - float * pMatT; - float * pATA; - double coef; - unsigned char *pdirImg = dirImg_.data; - /*If the first pixel in this chain is horizontal, - *then we try to find a horizontal line, y=ax+b;*/ - if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Horizontal ) - { - /*Build the new system,and solve it using least square regression: mat * [a,b]^T = vec - * [x0',1] [y0'] - * [x1',1] [a] [y1'] - * . [b] = . - * [xn',1] [yn']*/ - pMatT = matT.ptr(); //matT = [x0', x1', ... xn'; 1,1,...,1] - for ( int i = 0; i < newLength; i++ ) - { - * ( pMatT + newLength ) = 1; - * ( pMatT++ ) = (float) xCors[newOffsetS]; - vec[0][i] = (float) yCors[newOffsetS++]; - } - /* [a,b]^T = Inv(ATA + mat^T * mat) * (ATV + mat^T * vec) */ - tempMatLineFit = matT * matT.t(); - tempVecLineFit = matT * vec; - ATA = ATA + tempMatLineFit; - ATV = ATV + tempVecLineFit; - pATA = ATA.ptr(); - coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); - lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); - lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); - - return 0; - } - /*If the first pixel in this chain is vertical, - *then we try to find a vertical line, x=ay+b;*/ - if( pdirImg[yCors[offsetS] * imageWidth + xCors[offsetS]] == Vertical ) - { - /*Build the system,and solve it using least square regression: mat * [a,b]^T = vec - * [y0',1] [x0'] - * [y1',1] [a] [x1'] - * . [b] = . - * [yn',1] [xn']*/ - pMatT = matT.ptr(); //matT = [y0', y1', ... yn'; 1,1,...,1] - for ( int i = 0; i < newLength; i++ ) - { - * ( pMatT + newLength ) = 1; - * ( pMatT++ ) = (float) yCors[newOffsetS]; - vec[0][i] = (float) xCors[newOffsetS++]; - } - /* [a,b]^T = Inv(ATA + mat^T * mat) * (ATV + mat^T * vec) */ -// matT.MultiplyWithTransposeOf(matT, tempMatLineFit); - tempMatLineFit = matT * matT.t(); - tempVecLineFit = matT * vec; - ATA = ATA + tempMatLineFit; - ATV = ATV + tempVecLineFit; -// pATA = ATA.GetData(); - pATA = ATA.ptr(); - coef = 1.0 / ( double( pATA[0] ) * double( pATA[3] ) - double( pATA[1] ) * double( pATA[2] ) ); - lineEquation[0] = coef * ( double( pATA[3] ) * double( ATV[0][0] ) - double( pATA[1] ) * double( ATV[0][1] ) ); - lineEquation[1] = coef * ( double( pATA[0] ) * double( ATV[0][1] ) - double( pATA[2] ) * double( ATV[0][0] ) ); - - } - return 0; -} - -bool EDLineDetector::LineValidation_( unsigned int *xCors, unsigned int *yCors, unsigned int offsetS, unsigned int offsetE, - std::vector &lineEquation, float &direction ) -{ - if( bValidate_ ) - { - int n = offsetE - offsetS; - /*first compute the direction of line, make sure that the dark side always be the - *left side of a line.*/ - int meanGradientX = 0, meanGradientY = 0; - short *pdxImg = dxImg_.ptr(); - short *pdyImg = dyImg_.ptr(); - double dx, dy; - std::vector pointDirection; - int index; - for ( int i = 0; i < n; i++ ) - { - index = yCors[offsetS] * imageWidth + xCors[offsetS]; - offsetS++; - meanGradientX += pdxImg[index]; - meanGradientY += pdyImg[index]; - dx = (double) pdxImg[index]; - dy = (double) pdyImg[index]; - pointDirection.push_back( atan2( -dx, dy ) ); - } - dx = fabs( lineEquation[1] ); - dy = fabs( lineEquation[0] ); - if( meanGradientX == 0 && meanGradientY == 0 ) - { //not possible, if happens, it must be a wrong line, - return false; - } - if( meanGradientX > 0 && meanGradientY >= 0 ) - { //first quadrant, and positive direction of X axis. - direction = (float) atan2( -dy, dx ); //line direction is in fourth quadrant - } - if( meanGradientX <= 0 && meanGradientY > 0 ) - { //second quadrant, and positive direction of Y axis. - direction = (float) atan2( dy, dx ); //line direction is in first quadrant - } - if( meanGradientX < 0 && meanGradientY <= 0 ) - { //third quadrant, and negative direction of X axis. - direction = (float) atan2( dy, -dx ); //line direction is in second quadrant - } - if( meanGradientX >= 0 && meanGradientY < 0 ) - { //fourth quadrant, and negative direction of Y axis. - direction = (float) atan2( -dy, -dx ); //line direction is in third quadrant - } - /*then check whether the line is on the border of the image. We don't keep the border line.*/ - if( fabs( direction ) < 0.15 || M_PI - fabs( direction ) < 0.15 ) - { //Horizontal line - if( fabs( lineEquation[2] ) < 10 || fabs( imageHeight - fabs( lineEquation[2] ) ) < 10 ) - { //upper border or lower border - return false; - } - } - if( fabs( fabs( direction ) - M_PI * 0.5 ) < 0.15 ) - { //Vertical line - if( fabs( lineEquation[2] ) < 10 || fabs( imageWidth - fabs( lineEquation[2] ) ) < 10 ) - { //left border or right border - return false; - } - } - //count the aligned points on the line which have the same direction as the line. - double disDirection; - int k = 0; - for ( int i = 0; i < n; i++ ) - { - disDirection = fabs( direction - pointDirection[i] ); - if( fabs( 2 * M_PI - disDirection ) < 0.392699 || disDirection < 0.392699 ) - { //same direction, pi/8 = 0.392699081698724 - k++; - } - } - //now compute NFA(Number of False Alarms) - double ret = nfa( n, k, 0.125, logNT_ ); - - return ( ret > 0 ); //0 corresponds to 1 mean false alarm - } - else - { - return true; - } -} - -int EDLineDetector::EDline( cv::Mat &image ) -{ - if( ( EDline( image, lines_/*, smoothed*/ ) ) != 1 ) - { - return -1; - } - lineSalience_.clear(); - lineSalience_.resize( lines_.numOfLines ); - unsigned char *pgImg = gImgWO_.ptr(); - unsigned int indexInLineArray; - unsigned int *pXCor = lines_.xCors.data(); - unsigned int *pYCor = lines_.yCors.data(); - unsigned int *pSID = lines_.sId.data(); - for ( unsigned int i = 0; i < lineSalience_.size(); i++ ) - { - int salience = 0; - for ( indexInLineArray = pSID[i]; indexInLineArray < pSID[i + 1]; indexInLineArray++ ) - { - salience += pgImg[pYCor[indexInLineArray] * imageWidth + pXCor[indexInLineArray]]; - } - lineSalience_[i] = (float) salience; - } - return 1; -} - diff --git a/modules/line_descriptor/src/mihasher.cpp b/modules/line_descriptor/src/mihasher.cpp deleted file mode 100644 index 90f1cc71a..000000000 --- a/modules/line_descriptor/src/mihasher.cpp +++ /dev/null @@ -1,270 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Mohammad Norouzi, Ali Punjani, David J. Fleet, - // all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#include "precomp.hpp" - -/* execute a batch query */ -void Mihasher::batchquery( UINT32 * results, UINT32 *numres, const cv::Mat & queries, UINT32 numq, int dim1queries ) -{ - /* create and initialize a bitarray */ - counter = new bitarray; - counter->init( N ); - - UINT32 *res = new UINT32[K * ( D + 1 )]; - UINT64 *chunks = new UINT64[m]; - UINT32 * presults = results; - UINT32 *pnumres = numres; - - /* make a copy of input queries */ - cv::Mat queries_clone = queries.clone(); - - /* set a pointer to first query (row) */ - UINT8 *pq = queries_clone.ptr(); - - /* loop over number of descriptors */ - for ( size_t i = 0; i < numq; i++ ) - { - /* for every descriptor, query database */ - query( presults, pnumres, pq, chunks, res ); - - /* move pointer to write next K indeces */ - presults += K; - pnumres += B + 1; - - /* move forward pointer to current row in descriptors matrix */ - pq += dim1queries; - - } - - delete[] res; - delete[] chunks; - - delete counter; -} - -/* execute a single query */ -void Mihasher::query( UINT32* results, UINT32* numres, UINT8 * Query, UINT64 *chunks, UINT32 *res ) -{ - /* if K == 0 that means we want everything to be processed. - So maxres = N in that case. Otherwise K limits the results processed */ - UINT32 maxres = K ? K : (UINT32) N; - - /* number of results so far obtained (up to a distance of s per chunk) */ - UINT32 n = 0; - - /* number of candidates tested with full codes (not counting duplicates) */ - UINT32 nc = 0; - - /* counting everything retrieved (duplicates are counted multiple times) - number of lookups (and xors) */ - UINT32 nl = 0; - - UINT32 nd = 0; - UINT32 *arr; - int size = 0; - UINT32 index; - int hammd; - - counter->erase(); - memset( numres, 0, ( B + 1 ) * sizeof ( *numres ) ); - - split( chunks, Query, m, mplus, b ); - - /* the growing search radius per substring */ - int s; - - /* current b: for the first mplus substrings it is b, for the rest it is (b-1) */ - int curb = b; - - for ( s = 0; s <= d && n < maxres; s++ ) - { - for ( int k = 0; k < m; k++ ) - { - if( k < mplus ) - curb = b; - else - curb = b - 1; - UINT64 chunksk = chunks[k]; - /* number of bit-strings with s number of 1s */ - nl += xornum[s + 1] - xornum[s]; - - /* the bit-string with s number of 1s */ - UINT64 bitstr = 0; - for ( int i = 0; i < s; i++ ) - /* power[i] stores the location of the i'th 1 */ - power[i] = i; - /* used for stopping criterion (location of (s+1)th 1) */ - power[s] = curb + 1; - - /* bit determines the 1 that should be moving to the left */ - int bit = s - 1; - - /* start from the left-most 1, and move it to the left until - it touches another one */ - - /* the loop for changing bitstr */ - bool infiniteWhile = true; - while ( infiniteWhile ) - { - if( bit != -1 ) - { - bitstr ^= ( power[bit] == bit ) ? (UINT64) 1 << power[bit] : (UINT64) 3 << ( power[bit] - 1 ); - power[bit]++; - bit--; - } - - else - { /* bit == -1 */ - /* the binary code bitstr is available for processing */ - arr = H[k].query( chunksk ^ bitstr, &size ); // lookup - if( size ) - { /* the corresponding bucket is not empty */ - nd += size; - for ( int c = 0; c < size; c++ ) - { - index = arr[c]; - if( !counter->get( index ) ) - { /* if it is not a duplicate */ - counter->set( index ); - hammd = match( codes.ptr() + (UINT64) index * ( B_over_8 ), Query, B_over_8 ); - - nc++; - if( hammd <= D && numres[hammd] < maxres ) - res[hammd * K + numres[hammd]] = index + 1; - - numres[hammd]++; - } - } - } - - /* end of processing */ - while ( ++bit < s && power[bit] == power[bit + 1] - 1 ) - { - bitstr ^= (UINT64) 1 << ( power[bit] - 1 ); - power[bit] = bit; - } - if( bit == s ) - break; - } - } - - n = n + numres[s * m + k]; - if( n >= maxres ) - break; - } - } - - n = 0; - for ( s = 0; s <= D && (int) n < K; s++ ) - { - for ( int c = 0; c < (int) numres[s] && (int) n < K; c++ ) - results[n++] = res[s * K + c]; - } - -} - -/* constructor 2 */ -Mihasher::Mihasher( int _B, int _m ) -{ - B = _B; - B_over_8 = B / 8; - m = _m; - b = (int) ceil( (double) B / m ); - - /* assuming that B/2 is large enough radius to include - all of the k nearest neighbors */ - D = (int) ceil( B / 2.0 ); - d = (int) ceil( (double) D / m ); - - /* mplus is the number of chunks with b bits - (m-mplus) is the number of chunks with (b-1) bits */ - mplus = B - m * ( b - 1 ); - - xornum = new UINT32[d + 2]; - xornum[0] = 0; - for ( int i = 0; i <= d; i++ ) - xornum[i + 1] = xornum[i] + (UINT32) choose( b, i ); - - H = new SparseHashtable[m]; - - /* H[i].init might fail */ - for ( int i = 0; i < mplus; i++ ) - H[i].init( b ); - for ( int i = mplus; i < m; i++ ) - H[i].init( b - 1 ); -} - -/* K setter */ -void Mihasher::setK( int _K ) -{ - K = _K; -} - -/* desctructor */ -Mihasher::~Mihasher() -{ - delete[] xornum; - delete[] H; -} - -/* populate tables */ -void Mihasher::populate( cv::Mat & _codes, UINT32 _N, int dim1codes ) -{ - N = _N; - codes = _codes; - UINT64 * chunks = new UINT64[m]; - - UINT8 * pcodes = codes.ptr(); - for ( UINT64 i = 0; i < N; i++, pcodes += dim1codes ) - { - split( chunks, pcodes, m, mplus, b ); - - for ( int k = 0; k < m; k++ ) - H[k].insert( chunks[k], (UINT32) i ); - - if( i % (int) ceil( N / 1000.0 ) == 0 ) - fflush( stdout ); - } - - delete[] chunks; -} - diff --git a/modules/line_descriptor/src/precomp.hpp b/modules/line_descriptor/src/precomp.hpp index 024e89060..c9406b7e1 100644 --- a/modules/line_descriptor/src/precomp.hpp +++ b/modules/line_descriptor/src/precomp.hpp @@ -68,6 +68,10 @@ #include #include +#include "bitarray.hpp" +#include "bitops.hpp" +#include "types.hpp" + #include "opencv2/line_descriptor.hpp" #endif diff --git a/modules/line_descriptor/src/sparse_hashtable.cpp b/modules/line_descriptor/src/sparse_hashtable.cpp deleted file mode 100644 index 625c1e377..000000000 --- a/modules/line_descriptor/src/sparse_hashtable.cpp +++ /dev/null @@ -1,85 +0,0 @@ -/*M/////////////////////////////////////////////////////////////////////////////////////// - // - // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. - // - // By downloading, copying, installing or using the software you agree to this license. - // If you do not agree to this license, do not download, install, - // copy or use the software. - // - // - // License Agreement - // For Open Source Computer Vision Library - // - // Copyright (C) 2014, Biagio Montesano, all rights reserved. - // Third party copyrights are property of their respective owners. - // - // Redistribution and use in source and binary forms, with or without modification, - // are permitted provided that the following conditions are met: - // - // * Redistribution's of source code must retain the above copyright notice, - // this list of conditions and the following disclaimer. - // - // * Redistribution's in binary form must reproduce the above copyright notice, - // this list of conditions and the following disclaimer in the documentation - // and/or other materials provided with the distribution. - // - // * The name of the copyright holders may not be used to endorse or promote products - // derived from this software without specific prior written permission. - // - // This software is provided by the copyright holders and contributors "as is" and - // any express or implied warranties, including, but not limited to, the implied - // warranties of merchantability and fitness for a particular purpose are disclaimed. - // In no event shall the Intel Corporation or contributors be liable for any direct, - // indirect, incidental, special, exemplary, or consequential damages - // (including, but not limited to, procurement of substitute goods or services; - // loss of use, data, or profits; or business interruption) however caused - // and on any theory of liability, whether in contract, strict liability, - // or tort (including negligence or otherwise) arising in any way out of - // the use of this software, even if advised of the possibility of such damage. - // - //M*/ - -#include "precomp.hpp" - -const int SparseHashtable::MAX_B = 37; - -/* constructor */ -SparseHashtable::SparseHashtable() -{ - table = NULL; - size = 0; - b = 0; -} - -/* initializer */ -int SparseHashtable::init( int _b ) -{ - b = _b; - - if( b < 5 || b > MAX_B || b > (int) ( sizeof(UINT64) * 8 ) ) - return 1; - - size = UINT64_1 << ( b - 5 ); // size = 2 ^ b - table = (BucketGroup*) calloc( size, sizeof(BucketGroup) ); - - return 0; - -} - -/* destructor */ -SparseHashtable::~SparseHashtable() -{ - free( table ); -} - -/* insert data */ -void SparseHashtable::insert( UINT64 index, UINT32 data ) -{ - table[index >> 5].insert( (int) ( index % 32 ), data ); -} - -/* query data */ -UINT32* SparseHashtable::query( UINT64 index, int *Size ) -{ - return table[index >> 5].query( (int) ( index % 32 ), Size ); -} diff --git a/modules/line_descriptor/include/opencv2/line_descriptor/types.hpp b/modules/line_descriptor/src/types.hpp similarity index 100% rename from modules/line_descriptor/include/opencv2/line_descriptor/types.hpp rename to modules/line_descriptor/src/types.hpp