Skip to content

O
opencv_contrib
Commit f3b4dd1f authored by biagio montesano's avatar biagio montesano
Browse files
Options

EDLine porting completed

parent b2484cfd
Hide whitespace changes
Inline Side-by-side
Showing with 2522 additions and 470 deletions
modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp
View file @ f3b4dd1f
......@@ -123,6 +123,9 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
CV_PROP_RW
int reductionRatio;
CV_PROP_RW
int ksize_;
/* read parameters from a FileNode object and store them (struct function) */
void read( const FileNode& fn );
......@@ -131,6 +134,15 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
};
struct CV_EXPORTS LineDetectionMode
{
enum
{
LSD_DETECTOR = 0, // detect lines using LSD
EDL_DETECTOR = 1 // detect lines using EDLines
};
};
/* constructor */
CV_WRAP
BinaryDescriptor( const BinaryDescriptor::Params &parameters = BinaryDescriptor::Params() );
......@@ -158,18 +170,21 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
/* requires line detection (only one image) */
CV_WRAP
void detect( const Mat& image, CV_OUT std::vector<KeyLine>& keypoints, const Mat& mask = Mat() );
void detect( const Mat& image, CV_OUT std::vector<KeyLine>& keypoints, const Mat& mask = Mat(), int flags = LineDetectionMode::LSD_DETECTOR );
/* requires line detection (more than one image) */
void detect( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, const std::vector<Mat>& masks =
std::vector<Mat>() ) const;
void detect( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, const std::vector<Mat>& masks = std::vector<Mat>(),
int flags = LineDetectionMode::LSD_DETECTOR ) const;
/* requires descriptors computation (only one image) */
CV_WRAP
void compute( const Mat& image, CV_OUT CV_IN_OUT std::vector<KeyLine>& keylines, CV_OUT Mat& descriptors, bool returnFloatDescr=false ) const;
void compute( const Mat& image, CV_OUT CV_IN_OUT std::vector<KeyLine>& keylines, CV_OUT Mat& descriptors, bool returnFloatDescr = false, int flags =
LineDetectionMode::LSD_DETECTOR ) const;
/* requires descriptors computation (more than one image) */
void compute( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, std::vector<Mat>& descriptors, bool returnFloatDescr=false ) const;
void compute( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, std::vector<Mat>& descriptors, bool returnFloatDescr =
false,
int flags = LineDetectionMode::LSD_DETECTOR ) const;
/*return descriptor size */
int descriptorSize() const;
......@@ -186,14 +201,14 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
/* definition of operator () */
CV_WRAP_AS(detectAndCompute)
virtual void operator()( InputArray image, InputArray mask, CV_OUT std::vector<KeyLine>& keylines, OutputArray descriptors,
bool useProvidedKeyLines = false, bool returnFloatDescr=false ) const;
bool useProvidedKeyLines = false, bool returnFloatDescr = false, int flags = LineDetectionMode::LSD_DETECTOR ) const;
protected:
/* implementation of line detection */
virtual void detectImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, const Mat& mask = Mat() ) const;
virtual void detectImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, int flags, const Mat& mask = Mat() ) const;
/* implementation of descriptors' computation */
virtual void computeImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, Mat& descriptors, bool returnFloatDescr ) const;
virtual void computeImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, Mat& descriptors, bool returnFloatDescr, int flags ) const;
/* function inherited from Algorithm */
AlgorithmInfo* info() const;
......@@ -203,7 +218,10 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
unsigned char binaryConversion( float* f1, float* f2 );
/* compute LBD descriptors */
int computeLBD( ScaleLines &keyLines );
int computeLBD( ScaleLines &keyLines, int flags );
/* compute LBD descriptors using EDLine extractor */
int computeLBD_EDL( ScaleLines &keyLines );
/* compute Gaussian pyramid of input image */
void computeGaussianPyramid( const Mat& image );
......@@ -212,6 +230,10 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
Each group contains the same line, detected in different octaves */
int OctaveKeyLines( ScaleLines &keyLines );
/* gather lines in groups using EDLine extractor.
Each group contains the same line, detected in different octaves */
int OctaveKeyLines_EDL( cv::Mat& image, ScaleLines &keyLines );
/* get coefficients of line passing by two points (in line_extremes) */
void getLineParameters( cv::Vec4i &line_extremes, cv::Vec3i &lineParams );
......@@ -239,6 +261,12 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
/* vector to store the Gaussian pyramid od an input image */
std::vector<cv::Mat> octaveImages;
/*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<EDLineDetector*> edLineVec_;
};
class CV_EXPORTS_W BinaryDescriptorMatcher : public Algorithm
......
modules/line_descriptor/include/opencv2/line_descriptor/ed_line_detector.hpp
View file @ f3b4dd1f
......@@ -5,38 +5,38 @@
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.
*/
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_
......@@ -50,34 +50,38 @@ the use of this software, even if advised of the possibility of such damage.
#include <opencv2/highgui.hpp>
#include <opencv2/core/utility.hpp>
struct Pixel{
unsigned int x;//X coordinate
unsigned int y;//Y coordinate
struct Pixel
{
unsigned int x; //X coordinate
unsigned int y; //Y coordinate
};
struct EdgeChains{
std::vector<unsigned int> xCors;//all the x coordinates of edge points
std::vector<unsigned int> yCors;//all the y coordinates of edge points
std::vector<unsigned int> 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 EdgeChains
{
std::vector<unsigned int> xCors; //all the x coordinates of edge points
std::vector<unsigned int> yCors; //all the y coordinates of edge points
std::vector<unsigned int> 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<unsigned int> xCors;//all the x coordinates of line points
std::vector<unsigned int> yCors;//all the y coordinates of line points
std::vector<unsigned int> 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;
struct LineChains
{
std::vector<unsigned int> xCors; //all the x coordinates of line points
std::vector<unsigned int> yCors; //all the y coordinates of line points
std::vector<unsigned int> 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<Pixel> PixelChain;//each edge is a pixel chain
typedef std::list<Pixel> PixelChain; //each edge is a pixel chain
struct EDLineParam{
int ksize;
float sigma;
float gradientThreshold;
float anchorThreshold;
int scanIntervals;
int minLineLen;
double lineFitErrThreshold;
struct EDLineParam
{
int ksize;
float sigma;
float gradientThreshold;
float anchorThreshold;
int scanIntervals;
int minLineLen;
double lineFitErrThreshold;
};
#define RELATIVE_ERROR_FACTOR 100.0
......@@ -94,374 +98,377 @@ struct EDLineParam{
*/
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
*smoothed: In, flag to mark whether the image has already been smoothed by Gaussian filter.
*return -1: error happen
*/
int EdgeDrawing(cv::Mat &image, EdgeChains &edgeChains, bool smoothed=false);
/*extract lines from image
*image: In, gray image;
*lines: Out, store the extracted lines,
*smoothed: In, flag to mark whether the image has already been smoothed by Gaussian filter.
*return -1: error happen
*/
int EDline(cv::Mat &image, LineChains &lines, bool smoothed=false);
/* extract line from image, and store them */
int EDline(cv::Mat &image, bool smoothed=false);
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<std::vector<double> > lineEquations_ ;
//store the line endpoints, [x1,y1,x2,y3]
std::vector<std::vector<float> > lineEndpoints_ ;
//store the line direction
std::vector<float> lineDirection_;
//store the line salience, which is the summation of gradients of pixels on line
std::vector<float> 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<double> &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<double> &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<double> &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_<float> ATA; //the previous matrix of A^T * A;
cv::Mat_<float> ATV; //the previous vector of A^T * V;
cv::Mat_<float> fitMatT; //the matrix used in line fit function;
cv::Mat_<float> fitVec; //the vector used in line fit function;
cv::Mat_<float> tempMatLineFit; //the matrix used in line fit function;
cv::Mat_<float> 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 / M_LN10 - 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_j<mult_term_i for j>i.
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;
}
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
*smoothed: In, flag to mark whether the image has already been smoothed by Gaussian filter.
*return -1: error happen
*/
int EdgeDrawing( cv::Mat &image, EdgeChains &edgeChains, bool smoothed );
/*extract lines from image
*image: In, gray image;
*lines: Out, store the extracted lines,
*smoothed: In, flag to mark whether the image has already been smoothed by Gaussian filter.
*return -1: error happen
*/
int EDline( cv::Mat &image, LineChains &lines, bool smoothed );
/* extract line from image, and store them */
int EDline( cv::Mat &image, bool smoothed );
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<std::vector<double> > lineEquations_;
//store the line endpoints, [x1,y1,x2,y3]
std::vector<std::vector<float> > lineEndpoints_;
//store the line direction
std::vector<float> lineDirection_;
//store the line salience, which is the summation of gradients of pixels on line
std::vector<float> 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<double> &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<double> &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<double> &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_<float> ATA; //the previous matrix of A^T * A;
cv::Mat_<float> ATV; //the previous vector of A^T * V;
cv::Mat_<float> fitMatT; //the matrix used in line fit function;
cv::Mat_<float> fitVec; //the vector used in line fit function;
cv::Mat_<float> tempMatLineFit; //the matrix used in line fit function;
cv::Mat_<float> 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 / M_LN10 - 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_j<mult_term_i for j>i.
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_ */
modules/line_descriptor/samples/compute_descriptors.cpp
View file @ f3b4dd1f
......@@ -61,6 +61,20 @@ static void help()
<< std::endl;
}
inline void writeMat(cv::Mat m, std::string name, int n)
{
std::stringstream ss;
std::string s;
ss << n;
ss >> s;
std::string fileNameConf = name + s;
cv::FileStorage fsConf(fileNameConf, cv::FileStorage::WRITE);
fsConf << "m" << m;
fsConf.release();
}
int main( int argc, char** argv )
{
/* get parameters from command line */
......@@ -88,7 +102,7 @@ int main( int argc, char** argv )
/* compute lines */
std::vector<KeyLine> keylines;
bd->detect( imageMat, keylines, mask );
bd->detect( imageMat, keylines, mask, 1 );
/* select only lines from first octave */
std::vector<KeyLine> octave0;
......@@ -101,6 +115,7 @@ int main( int argc, char** argv )
/* compute descriptors */
cv::Mat descriptors;
bd->compute( imageMat, octave0, descriptors );
bd->compute( imageMat, octave0, descriptors, false, 1 );
writeMat(descriptors, "bd_descriptors", 0);
}
modules/line_descriptor/samples/lines_extraction.cpp
View file @ f3b4dd1f
......@@ -90,10 +90,11 @@ int main( int argc, char** argv )
vector<KeyLine> lines;
/* extract lines */
bd->detect( imageMat, lines, mask );
cv::Mat output = imageMat.clone();
bd->detect( imageMat, lines, mask, 1 );
/* draw lines extracted from octave 0 */
cv::Mat output = imageMat.clone();
if( output.channels() == 1 )
cvtColor( output, output, COLOR_GRAY2BGR );
for ( size_t i = 0; i < lines.size(); i++ )
......
modules/line_descriptor/src/binary_descriptor.cpp
View file @ f3b4dd1f
......@@ -88,6 +88,7 @@ BinaryDescriptor::Params::Params()
numOfOctave_ = 1;
widthOfBand_ = 7;
reductionRatio = 2;
ksize_ = 5;
}
/* setters and getters */
......@@ -153,6 +154,12 @@ BinaryDescriptor::BinaryDescriptor( const BinaryDescriptor::Params &parameters )
params( parameters )
{
/* reserve enough space for EDLine objects */
edLineVec_.resize( params.numOfOctave_ );
images_sizes.resize( params.numOfOctave_ );
for ( unsigned int i = 0; i < params.numOfOctave_; i++ )
edLineVec_[i] = new EDLineDetector;
/* prepare a vector to host local weights F_l*/
gaussCoefL_.resize( params.widthOfBand_ * 3 );
......@@ -189,7 +196,7 @@ BinaryDescriptor::BinaryDescriptor( const BinaryDescriptor::Params &parameters )
/* definition of operator () */
void BinaryDescriptor::operator()( InputArray image, InputArray mask, CV_OUT std::vector<KeyLine>& keylines, OutputArray descriptors,
bool useProvidedKeyLines, bool returnFloatDescr ) const
bool useProvidedKeyLines, bool returnFloatDescr, int flags ) const
{
/* create some matrix objects */
......@@ -207,10 +214,10 @@ void BinaryDescriptor::operator()( InputArray image, InputArray mask, CV_OUT std
/* require drawing KeyLines detection if demanded */
if( !useProvidedKeyLines )
detectImpl( imageMat, keylines, maskMat );
detectImpl( imageMat, keylines, flags, maskMat );
/* compute descriptors */
computeImpl( imageMat, keylines, descrMat, returnFloatDescr );
computeImpl( imageMat, keylines, descrMat, returnFloatDescr, flags );
}
BinaryDescriptor::~BinaryDescriptor()
......@@ -382,7 +389,7 @@ inline void checkLineExtremes( cv::Vec4i& extremes, cv::Size imageSize )
}
/* requires line detection (only one image) */
void BinaryDescriptor::detect( const Mat& image, CV_OUT std::vector<KeyLine>& keylines, const Mat& mask )
void BinaryDescriptor::detect( const Mat& image, CV_OUT std::vector<KeyLine>& keylines, const Mat& mask, int flags )
{
if( mask.data != NULL && ( mask.size() != image.size() || mask.type() != CV_8UC1 ) )
{
......@@ -393,11 +400,12 @@ void BinaryDescriptor::detect( const Mat& image, CV_OUT std::vector<KeyLine>& ke
}
else
detectImpl( image, keylines, mask );
detectImpl( image, keylines, flags, mask );
}
/* requires line detection (more than one image) */
void BinaryDescriptor::detect( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, const std::vector<Mat>& masks ) const
void BinaryDescriptor::detect( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, const std::vector<Mat>& masks,
int flags ) const
{
/* detect lines from each image */
for ( size_t counter = 0; counter < images.size(); counter++ )
......@@ -410,11 +418,11 @@ void BinaryDescriptor::detect( const std::vector<Mat>& images, std::vector<std::
CV_Assert( false );
}
detectImpl( images[counter], keylines[counter], masks[counter] );
detectImpl( images[counter], keylines[counter], flags, masks[counter] );
}
}
void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, const Mat& mask ) const
void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, int flags, const Mat& mask ) const
{
cv::Mat image;
......@@ -434,11 +442,32 @@ void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector<KeyLine>& ke
BinaryDescriptor *bn = const_cast<BinaryDescriptor*>( this );
/* compute Gaussian pyramid */
bn->computeGaussianPyramid( image );
if( flags == 0 )
bn->computeGaussianPyramid( image );
/* detect and arrange lines across octaves */
ScaleLines sl;
bn->OctaveKeyLines( sl );
Mat m = image.clone();
cvtColor( m, m, COLOR_GRAY2BGR );
if( flags == 0 )
bn->OctaveKeyLines( sl );
else
bn->OctaveKeyLines_EDL( image, sl );
Mat temp = image.clone();
cvtColor( temp, temp, COLOR_GRAY2BGR );
for ( size_t i = 0; i < sl.size(); i++ )
{
for ( size_t j = 0; j < sl[i].size(); j++ )
{
OctaveSingleLine tempOSL = sl[i][j];
line( m, Point( tempOSL.startPointX, tempOSL.startPointY ), Point( tempOSL.endPointX, tempOSL.endPointY ), Scalar( 255, 0, 0 ), 5 );
}
}
imshow( "Immagine", m );
waitKey();
/* fill KeyLines vector */
for ( int i = 0; i < (int) sl.size(); i++ )
......@@ -452,14 +481,13 @@ void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector<KeyLine>& ke
KeyLine kl;
/* check data validity */
cv::Vec4i extremes( osl.startPointX, osl.startPointY, osl.endPointX, osl.endPointY );
checkLineExtremes( extremes, image.size() );
// cv::Vec4i extremes( osl.startPointX, osl.startPointY, osl.endPointX, osl.endPointY );
// checkLineExtremes( extremes, imageSrc.size() );
/* fill KeyLine's fields */
kl.startPointX = extremes[0];
kl.startPointY = extremes[1];
kl.endPointX = extremes[2];
kl.endPointY = extremes[3];
kl.startPointX = osl.startPointX; //extremes[0];
kl.startPointY = osl.startPointY; //extremes[1];
kl.endPointX = osl.endPointX; //extremes[2];
kl.endPointY = osl.endPointY; //extremes[3];
kl.sPointInOctaveX = osl.sPointInOctaveX;
kl.sPointInOctaveY = osl.sPointInOctaveY;
kl.ePointInOctaveX = osl.ePointInOctaveX;
......@@ -494,22 +522,22 @@ void BinaryDescriptor::detectImpl( const Mat& imageSrc, std::vector<KeyLine>& ke
}
/* requires descriptors computation (only one image) */
void BinaryDescriptor::compute( const Mat& image, CV_OUT CV_IN_OUT std::vector<KeyLine>& keylines, CV_OUT Mat& descriptors,
bool returnFloatDescr ) const
void BinaryDescriptor::compute( const Mat& image, CV_OUT CV_IN_OUT std::vector<KeyLine>& keylines, CV_OUT Mat& descriptors, bool returnFloatDescr,
int flags ) const
{
computeImpl( image, keylines, descriptors, returnFloatDescr );
computeImpl( image, keylines, descriptors, returnFloatDescr, flags );
}
/* requires descriptors computation (more than one image) */
void BinaryDescriptor::compute( const std::vector<Mat>& images, std::vector<std::vector<KeyLine> >& keylines, std::vector<Mat>& descriptors,
bool returnFloatDescr ) const
bool returnFloatDescr, int flags ) const
{
for ( size_t i = 0; i < images.size(); i++ )
computeImpl( images[i], keylines[i], descriptors[i], returnFloatDescr );
computeImpl( images[i], keylines[i], descriptors[i], returnFloatDescr, flags );
}
/* implementation of descriptors computation */
void BinaryDescriptor::computeImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, Mat& descriptors, bool returnFloatDescr ) const
void BinaryDescriptor::computeImpl( const Mat& imageSrc, std::vector<KeyLine>& keylines, Mat& descriptors, bool returnFloatDescr, int flags ) const
{
/* convert input image to gray scale */
cv::Mat image;
......@@ -615,7 +643,18 @@ void BinaryDescriptor::computeImpl( const Mat& imageSrc, std::vector<KeyLine>& k
}
/* compute LBD descriptors */
bn->computeLBD( sl );
if(flags == 0)
bn->computeLBD( sl, flags );
else
bn->computeLBD_EDL(sl);
/* resize output matrix */
if( !returnFloatDescr )
descriptors = cv::Mat( keylines.size(), 32, CV_8UC1 );
else
descriptors = cv::Mat( keylines.size(), NUM_OF_BANDS * 8, CV_32FC1 );
/* fill output matrix with descriptors */
for ( size_t k = 0; k < sl.size(); k++ )
......@@ -628,8 +667,6 @@ void BinaryDescriptor::computeImpl( const Mat& imageSrc, std::vector<KeyLine>& k
if( !returnFloatDescr )
{
/* resize output matrix */
descriptors = cv::Mat( keylines.size(), 32, CV_8UC1 );
/* get a pointer to correspondent row in output matrix */
uchar* pointerToRow = descriptors.ptr( originalIndex );
......@@ -648,8 +685,6 @@ void BinaryDescriptor::computeImpl( const Mat& imageSrc, std::vector<KeyLine>& k
else
{
/* resize output matrix */
descriptors = cv::Mat( keylines.size(), NUM_OF_BANDS * 8, CV_32FC1 );
/* get a pointer to correspondent row in output matrix */
uchar* pointerToRow = descriptors.ptr( originalIndex );
......@@ -689,7 +724,7 @@ int BinaryDescriptor::OctaveKeyLines( ScaleLines &keyLines )
std::vector<cv::Vec4i> lines_std;
/* use detector to extract segments */
ls->detect( currentScaledImage, lines_std, w_idth, prec/*, nfa*/ );
ls->detect( currentScaledImage, lines_std, w_idth, prec/*, nfa*/);
/* store lines extracted from current image */
extractedLines.push_back( lines_std );
......@@ -999,8 +1034,350 @@ int BinaryDescriptor::OctaveKeyLines( ScaleLines &keyLines )
}
int BinaryDescriptor::OctaveKeyLines_EDL( cv::Mat& image, ScaleLines &keyLines )
{
/* final number of extracted lines */
unsigned int numOfFinalLine = 0;
/* sigma values and reduction factor used in Gaussian pyramids */
float preSigma2 = 0; //orignal image is not blurred, has zero sigma;
float curSigma2 = 1.0; //[sqrt(2)]^0=1;
float factor = sqrt( 2 ); //the down sample factor between connective two octave images
/* loop over number of octaves */
for ( int octaveCount = 0; octaveCount < params.numOfOctave_; octaveCount++ )
{
/* matrix storing results from blurring processes */
cv::Mat blur;
/* apply Gaussian blur */
float increaseSigma = sqrt( curSigma2 - preSigma2 );
cv::GaussianBlur( image, blur, cv::Size( params.ksize_, params.ksize_ ), increaseSigma );
images_sizes[octaveCount] = blur.size();
/* for current octave, extract lines */
if( ( edLineVec_[octaveCount]->EDline( blur, true ) ) != true )
{
return -1;
}
/* update number of total extracted lines */
numOfFinalLine += edLineVec_[octaveCount]->lines_.numOfLines;
/* resize image for next level of pyramid */
cv::resize( blur, image, cv::Size(), ( 1.f / factor ), ( 1.f / factor ) );
/* update sigma values */
preSigma2 = curSigma2;
curSigma2 = curSigma2 * 2;
} /* end of loop over number of octaves */
/*lines which correspond to the same line in the octave images will be stored
in the same element of ScaleLines.*/
/* prepare a vector to store octave information associated to extracted lines */
std::vector<OctaveLine> octaveLines( numOfFinalLine );
/* set lines' counter to 0 for reuse */
numOfFinalLine = 0;
/* counter to give a unique ID to lines in LineVecs */
unsigned int lineIDInScaleLineVec = 0;
/* floats to compute lines' lengths */
float dx, dy;
/* loop over lines extracted from scale 0 (original image) */
for ( unsigned int lineCurId = 0; lineCurId < edLineVec_[0]->lines_.numOfLines; lineCurId++ )
{
/* FOR CURRENT LINE: */
/* set octave from which it was extracted */
octaveLines[numOfFinalLine].octaveCount = 0;
/* set ID within its octave */
octaveLines[numOfFinalLine].lineIDInOctave = lineCurId;
/* set a unique ID among all lines extracted in all octaves */
octaveLines[numOfFinalLine].lineIDInScaleLineVec = lineIDInScaleLineVec;
/* compute absolute value of difference between X coordinates of line's extreme points */
dx = fabs( edLineVec_[0]->lineEndpoints_[lineCurId][0] - edLineVec_[0]->lineEndpoints_[lineCurId][2] );
/* compute absolute value of difference between Y coordinates of line's extreme points */
dy = fabs( edLineVec_[0]->lineEndpoints_[lineCurId][1] - edLineVec_[0]->lineEndpoints_[lineCurId][3] );
/* compute line's length */
octaveLines[numOfFinalLine].lineLength = sqrt( dx * dx + dy * dy );
/* update counters */
numOfFinalLine++;
lineIDInScaleLineVec++;
}
/* create and fill an array to store scale factors */
float *scale = new float[params.numOfOctave_];
scale[0] = 1;
for ( int octaveCount = 1; octaveCount < params.numOfOctave_; octaveCount++ )
{
scale[octaveCount] = factor * scale[octaveCount - 1];
}
/* some variables' declarations */
float rho1, rho2, tempValue;
float direction, near, length;
unsigned int octaveID, lineIDInOctave;
/*more than one octave image, organize lines in scale space.
*lines corresponding to the same line in octave images should have the same index in the ScaleLineVec */
if( params.numOfOctave_ > 1 )
{
/* some other variables' declarations */
float twoPI = 2 * M_PI;
unsigned int closeLineID;
float endPointDis, minEndPointDis, minLocalDis, maxLocalDis;
float lp0, lp1, lp2, lp3, np0, np1, np2, np3;
/* loop over list of octaves */
for ( int octaveCount = 1; octaveCount < params.numOfOctave_; octaveCount++ )
{
/*for each line in current octave image, find their corresponding lines in the octaveLines,
*give them the same value of lineIDInScaleLineVec*/
/* loop over list of lines extracted from current octave */
for ( unsigned int lineCurId = 0; lineCurId < edLineVec_[octaveCount]->lines_.numOfLines; lineCurId++ )
{
/* get (scaled) known term from equation of current line */
rho1 = scale[octaveCount] * fabs( edLineVec_[octaveCount]->lineEquations_[lineCurId][2] );
/*nearThreshold depends on the distance of the image coordinate origin to current line.
*so nearThreshold = rho1 * nearThresholdRatio, where nearThresholdRatio = 1-cos(10*pi/180) = 0.0152*/
tempValue = rho1 * 0.0152;
float nearThreshold = ( tempValue > 6 ) ? ( tempValue ) : 6;
nearThreshold = ( nearThreshold < 12 ) ? nearThreshold : 12;
/* compute scaled lenght of current line */
dx = fabs( edLineVec_[octaveCount]->lineEndpoints_[lineCurId][0] - edLineVec_[octaveCount]->lineEndpoints_[lineCurId][2] ); //x1-x2
dy = fabs( edLineVec_[octaveCount]->lineEndpoints_[lineCurId][1] - edLineVec_[octaveCount]->lineEndpoints_[lineCurId][3] ); //y1-y2
length = scale[octaveCount] * sqrt( dx * dx + dy * dy );
minEndPointDis = 12;
/* loop over the octave representations of all lines */
for ( unsigned int lineNextId = 0; lineNextId < numOfFinalLine; lineNextId++ )
{
/* if a line from same octave is encountered,
a comparison with it shouldn't be considered */
octaveID = octaveLines[lineNextId].octaveCount;
if( (int) octaveID == octaveCount )
{ //lines in the same layer of octave image should not be compared.
break;
}
/* take ID in octave of line to be compared */
lineIDInOctave = octaveLines[lineNextId].lineIDInOctave;
/*first check whether current line and next line are parallel.
*If line1:a1*x+b1*y+c1=0 and line2:a2*x+b2*y+c2=0 are parallel, then
*-a1/b1=-a2/b2, i.e., a1b2=b1a2.
*we define parallel=fabs(a1b2-b1a2)
*note that, in EDLine class, we have normalized the line equations
*to make a1^2+ b1^2 = a2^2+ b2^2 = 1*/
direction = fabs( edLineVec_[octaveCount]->lineDirection_[lineCurId] - edLineVec_[octaveID]->lineDirection_[lineIDInOctave] );
/* the angle between two lines are larger than 10degrees
(i.e. 10*pi/180=0.1745), they are not close to parallel */
if( direction > 0.1745 && ( twoPI - direction > 0.1745 ) )
{
continue;
}
/*now check whether current line and next line are near to each other.
*If line1:a1*x+b1*y+c1=0 and line2:a2*x+b2*y+c2=0 are near in image, then
*rho1 = |a1*0+b1*0+c1|/sqrt(a1^2+b1^2) and rho2 = |a2*0+b2*0+c2|/sqrt(a2^2+b2^2) should close.
*In our case, rho1 = |c1| and rho2 = |c2|, because sqrt(a1^2+b1^2) = sqrt(a2^2+b2^2) = 1;
*note that, lines are in different octave images, so we define near = fabs(scale*rho1 - rho2) or
*where scale is the scale factor between to octave images*/
/* get known term from equation to be compared */
rho2 = scale[octaveID] * fabs( edLineVec_[octaveID]->lineEquations_[lineIDInOctave][2] );
/* compute difference between known ters */
near = fabs( rho1 - rho2 );
/* two lines are not close in the image */
if( near > nearThreshold )
{
continue;
}
/*now check the end points distance between two lines, the scale of distance is in the original image size.
* find the minimal and maximal end points distance*/
/* get the extreme points of the two lines */
lp0 = scale[octaveCount] * edLineVec_[octaveCount]->lineEndpoints_[lineCurId][0];
lp1 = scale[octaveCount] * edLineVec_[octaveCount]->lineEndpoints_[lineCurId][1];
lp2 = scale[octaveCount] * edLineVec_[octaveCount]->lineEndpoints_[lineCurId][2];
lp3 = scale[octaveCount] * edLineVec_[octaveCount]->lineEndpoints_[lineCurId][3];
np0 = scale[octaveID] * edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][0];
np1 = scale[octaveID] * edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][1];
np2 = scale[octaveID] * edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][2];
np3 = scale[octaveID] * edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][3];
/* get the distance between the two leftmost extremes of lines
L1(0,1)<->L2(0,1) */
dx = lp0 - np0;
dy = lp1 - np1;
endPointDis = sqrt( dx * dx + dy * dy );
/* set momentaneously min and max distance between lines to
the one between left extremes */
minLocalDis = endPointDis;
maxLocalDis = endPointDis;
/* compute distance between right extremes
L1(2,3)<->L2(2,3) */
dx = lp2 - np2;
dy = lp3 - np3;
endPointDis = sqrt( dx * dx + dy * dy );
/* update (if necessary) min and max distance between lines */
minLocalDis = ( endPointDis < minLocalDis ) ? endPointDis : minLocalDis;
maxLocalDis = ( endPointDis > maxLocalDis ) ? endPointDis : maxLocalDis;
/* compute distance between left extreme of current line and
right extreme of line to be compared
L1(0,1)<->L2(2,3) */
dx = lp0 - np2;
dy = lp1 - np3;
endPointDis = sqrt( dx * dx + dy * dy );
/* update (if necessary) min and max distance between lines */
minLocalDis = ( endPointDis < minLocalDis ) ? endPointDis : minLocalDis;
maxLocalDis = ( endPointDis > maxLocalDis ) ? endPointDis : maxLocalDis;
/* compute distance between right extreme of current line and
left extreme of line to be compared
L1(2,3)<->L2(0,1) */
dx = lp2 - np0;
dy = lp3 - np1;
endPointDis = sqrt( dx * dx + dy * dy );
/* update (if necessary) min and max distance between lines */
minLocalDis = ( endPointDis < minLocalDis ) ? endPointDis : minLocalDis;
maxLocalDis = ( endPointDis > maxLocalDis ) ? endPointDis : maxLocalDis;
/* check whether conditions for considering line to be compared
worth to be inserted in the same LineVec are satisfied */
if( ( maxLocalDis < 0.8 * ( length + octaveLines[lineNextId].lineLength ) ) && ( minLocalDis < minEndPointDis ) )
{ //keep the closest line
minEndPointDis = minLocalDis;
closeLineID = lineNextId;
}
}
/* add current line into octaveLines */
if( minEndPointDis < 12 )
{
octaveLines[numOfFinalLine].lineIDInScaleLineVec = octaveLines[closeLineID].lineIDInScaleLineVec;
}
else
{
octaveLines[numOfFinalLine].lineIDInScaleLineVec = lineIDInScaleLineVec;
lineIDInScaleLineVec++;
}
octaveLines[numOfFinalLine].octaveCount = octaveCount;
octaveLines[numOfFinalLine].lineIDInOctave = lineCurId;
octaveLines[numOfFinalLine].lineLength = length;
numOfFinalLine++;
}
} //end for(unsigned int octaveCount = 1; octaveCount<numOfOctave_; octaveCount++)
} //end if(numOfOctave_>1)
////////////////////////////////////
//Reorganize the detected lines into keyLines
keyLines.clear();
keyLines.resize( lineIDInScaleLineVec );
unsigned int tempID;
float s1, e1, s2, e2;
bool shouldChange;
OctaveSingleLine singleLine;
for ( unsigned int lineID = 0; lineID < numOfFinalLine; lineID++ )
{
lineIDInOctave = octaveLines[lineID].lineIDInOctave;
octaveID = octaveLines[lineID].octaveCount;
direction = edLineVec_[octaveID]->lineDirection_[lineIDInOctave];
singleLine.octaveCount = octaveID;
singleLine.direction = direction;
singleLine.lineLength = octaveLines[lineID].lineLength;
singleLine.salience = edLineVec_[octaveID]->lineSalience_[lineIDInOctave];
singleLine.numOfPixels = edLineVec_[octaveID]->lines_.sId[lineIDInOctave + 1] - edLineVec_[octaveID]->lines_.sId[lineIDInOctave];
//decide the start point and end point
shouldChange = false;
s1 = edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][0]; //sx
s2 = edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][1]; //sy
e1 = edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][2]; //ex
e2 = edLineVec_[octaveID]->lineEndpoints_[lineIDInOctave][3]; //ey
dx = e1 - s1; //ex-sx
dy = e2 - s2; //ey-sy
if( direction >= -0.75 * M_PI && direction < -0.25 * M_PI )
{
if( dy > 0 )
{
shouldChange = true;
}
}
if( direction >= -0.25 * M_PI && direction < 0.25 * M_PI )
{
if( dx < 0 )
{
shouldChange = true;
}
}
if( direction >= 0.25 * M_PI && direction < 0.75 * M_PI )
{
if( dy < 0 )
{
shouldChange = true;
}
}
if( ( direction >= 0.75 * M_PI && direction < M_PI ) || ( direction >= -M_PI && direction < -0.75 * M_PI ) )
{
if( dx > 0 )
{
shouldChange = true;
}
}
tempValue = scale[octaveID];
if( shouldChange )
{
singleLine.sPointInOctaveX = e1;
singleLine.sPointInOctaveY = e2;
singleLine.ePointInOctaveX = s1;
singleLine.ePointInOctaveY = s2;
singleLine.startPointX = tempValue * e1;
singleLine.startPointY = tempValue * e2;
singleLine.endPointX = tempValue * s1;
singleLine.endPointY = tempValue * s2;
}
else
{
singleLine.sPointInOctaveX = s1;
singleLine.sPointInOctaveY = s2;
singleLine.ePointInOctaveX = e1;
singleLine.ePointInOctaveY = e2;
singleLine.startPointX = tempValue * s1;
singleLine.startPointY = tempValue * s2;
singleLine.endPointX = tempValue * e1;
singleLine.endPointY = tempValue * e2;
}
tempID = octaveLines[lineID].lineIDInScaleLineVec;
keyLines[tempID].push_back( singleLine );
}
////////////////////////////////////
delete[] scale;
return 1;
}
/* compute LBD descriptors */
int BinaryDescriptor::computeLBD( ScaleLines &keyLines )
int BinaryDescriptor::computeLBD( ScaleLines &keyLines, int flags )
{
//the default length of the band is the line length.
short numOfFinalLine = keyLines.size();
......@@ -1056,11 +1433,18 @@ int BinaryDescriptor::computeLBD( ScaleLines &keyLines )
pSingleLine = & ( keyLines[lineIDInScaleVec][lineIDInSameLine] );
octaveCount = pSingleLine->octaveCount;
/* retrieve associated dxImg and dyImg
pdxImg = edLineVec_[octaveCount]->dxImg_.ptr<short>();
pdyImg = edLineVec_[octaveCount]->dyImg_.ptr<short>(); */
pdxImg = dxImg_vector[octaveCount].ptr<short>();
pdyImg = dyImg_vector[octaveCount].ptr<short>();
/* retrieve associated dxImg and dyImg*/
if( flags == 1 )
{
pdxImg = edLineVec_[octaveCount]->dxImg_.ptr<short>();
pdyImg = edLineVec_[octaveCount]->dyImg_.ptr<short>();
}
else
{
pdxImg = dxImg_vector[octaveCount].ptr<short>();
pdyImg = dyImg_vector[octaveCount].ptr<short>();
}
/* get image size to work on from real one
realWidth = edLineVec_[octaveCount]->imageWidth;
......@@ -1331,3 +1715,313 @@ int BinaryDescriptor::computeLBD( ScaleLines &keyLines )
return 1;
}
int BinaryDescriptor::computeLBD_EDL( ScaleLines &keyLines )
{
//the default length of the band is the line length.
short numOfFinalLine = keyLines.size();
float *dL = new float[2];//line direction cos(dir), sin(dir)
float *dO = new float[2];//the clockwise orthogonal vector of line direction.
short heightOfLSP = params.widthOfBand_*NUM_OF_BANDS;//the height of line support region;
short descriptorSize = NUM_OF_BANDS * 8;//each band, we compute the m( pgdL, ngdL, pgdO, ngdO) and std( pgdL, ngdL, pgdO, ngdO);
float pgdLRowSum;//the summation of {g_dL |g_dL>0 } for each row of the region;
float ngdLRowSum;//the summation of {g_dL |g_dL<0 } for each row of the region;
float pgdL2RowSum;//the summation of {g_dL^2 |g_dL>0 } for each row of the region;
float ngdL2RowSum;//the summation of {g_dL^2 |g_dL<0 } for each row of the region;
float pgdORowSum;//the summation of {g_dO |g_dO>0 } for each row of the region;
float ngdORowSum;//the summation of {g_dO |g_dO<0 } for each row of the region;
float pgdO2RowSum;//the summation of {g_dO^2 |g_dO>0 } for each row of the region;
float ngdO2RowSum;//the summation of {g_dO^2 |g_dO<0 } for each row of the region;
float *pgdLBandSum = new float[NUM_OF_BANDS];//the summation of {g_dL |g_dL>0 } for each band of the region;
float *ngdLBandSum = new float[NUM_OF_BANDS];//the summation of {g_dL |g_dL<0 } for each band of the region;
float *pgdL2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dL^2 |g_dL>0 } for each band of the region;
float *ngdL2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dL^2 |g_dL<0 } for each band of the region;
float *pgdOBandSum = new float[NUM_OF_BANDS];//the summation of {g_dO |g_dO>0 } for each band of the region;
float *ngdOBandSum = new float[NUM_OF_BANDS];//the summation of {g_dO |g_dO<0 } for each band of the region;
float *pgdO2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dO^2 |g_dO>0 } for each band of the region;
float *ngdO2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dO^2 |g_dO<0 } for each band of the region;
short numOfBitsBand = NUM_OF_BANDS*sizeof(float);
short lengthOfLSP; //the length of line support region, varies with lines
short halfHeight = (heightOfLSP-1)/2;
short halfWidth;
short bandID;
float coefInGaussion;
float lineMiddlePointX, lineMiddlePointY;
float sCorX, sCorY,sCorX0, sCorY0;
short tempCor, xCor, yCor;//pixel coordinates in image plane
short dx, dy;
float gDL;//store the gradient projection of pixels in support region along dL vector
float gDO;//store the gradient projection of pixels in support region along dO vector
short imageWidth, imageHeight, realWidth;
short *pdxImg, *pdyImg;
float *desVec;
short sameLineSize;
short octaveCount;
OctaveSingleLine *pSingleLine;
/* loop over list of LineVec */
for(short lineIDInScaleVec = 0; lineIDInScaleVec<numOfFinalLine; lineIDInScaleVec++){
sameLineSize = keyLines[lineIDInScaleVec].size();
/* loop over current LineVec's lines */
for(short lineIDInSameLine = 0; lineIDInSameLine<sameLineSize; lineIDInSameLine++){
/* get a line in current LineVec and its original ID in its octave */
pSingleLine = &(keyLines[lineIDInScaleVec][lineIDInSameLine]);
octaveCount = pSingleLine->octaveCount;
/* retrieve associated dxImg and dyImg */
pdxImg = edLineVec_[octaveCount]->dxImg_.ptr<short>();
pdyImg = edLineVec_[octaveCount]->dyImg_.ptr<short>();
/* get image size to work on from real one */
realWidth = edLineVec_[octaveCount]->imageWidth;
imageWidth = realWidth -1;
imageHeight = edLineVec_[octaveCount]->imageHeight-1;
/* initialize memory areas */
memset(pgdLBandSum, 0, numOfBitsBand);
memset(ngdLBandSum, 0, numOfBitsBand);
memset(pgdL2BandSum, 0, numOfBitsBand);
memset(ngdL2BandSum, 0, numOfBitsBand);
memset(pgdOBandSum, 0, numOfBitsBand);
memset(ngdOBandSum, 0, numOfBitsBand);
memset(pgdO2BandSum, 0, numOfBitsBand);
memset(ngdO2BandSum, 0, numOfBitsBand);
/* get length of line and its half */
lengthOfLSP = keyLines[lineIDInScaleVec][lineIDInSameLine].numOfPixels;
halfWidth = (lengthOfLSP-1)/2;
/* get middlepoint of line */
lineMiddlePointX = 0.5 * (pSingleLine->sPointInOctaveX + pSingleLine->ePointInOctaveX);
lineMiddlePointY = 0.5 * (pSingleLine->sPointInOctaveY + pSingleLine->ePointInOctaveY);
/*1.rotate the local coordinate system to the line direction (direction is the angle
between positive line direction and positive X axis)
*2.compute the gradient projection of pixels in line support region*/
/* get the vector representing original image reference system after rotation to aligh with
line's direction */
dL[0] = cos(pSingleLine->direction);
dL[1] = sin(pSingleLine->direction);
/* set the clockwise orthogonal vector of line direction */
dO[0] = -dL[1];
dO[1] = dL[0];
/* get rotated reference frame */
sCorX0= -dL[0]*halfWidth + dL[1]*halfHeight + lineMiddlePointX;//hID =0; wID = 0;
sCorY0= -dL[1]*halfWidth - dL[0]*halfHeight + lineMiddlePointY;
/* BIAS::Matrix<float> gDLMat(heightOfLSP,lengthOfLSP) */
for(short hID = 0; hID <heightOfLSP; hID++){
/*initialization */
sCorX = sCorX0;
sCorY = sCorY0;
pgdLRowSum = 0;
ngdLRowSum = 0;
pgdORowSum = 0;
ngdORowSum = 0;
for(short wID = 0; wID <lengthOfLSP; wID++){
tempCor = round(sCorX);
xCor = (tempCor<0)?0:(tempCor>imageWidth)?imageWidth:tempCor;
tempCor = round(sCorY);
yCor = (tempCor<0)?0:(tempCor>imageHeight)?imageHeight:tempCor;
/* To achieve rotation invariance, each simple gradient is rotated aligned with
* the line direction and clockwise orthogonal direction.*/
dx = pdxImg[yCor*realWidth+xCor];
dy = pdyImg[yCor*realWidth+xCor];
gDL = dx * dL[0] + dy * dL[1];
gDO = dx * dO[0] + dy * dO[1];
if(gDL>0){
pgdLRowSum += gDL;
}else{
ngdLRowSum -= gDL;
}
if(gDO>0){
pgdORowSum += gDO;
}else{
ngdORowSum -= gDO;
}
sCorX +=dL[0];
sCorY +=dL[1];
/* gDLMat[hID][wID] = gDL; */
}
sCorX0 -=dL[1];
sCorY0 +=dL[0];
coefInGaussion = gaussCoefG_[hID];
pgdLRowSum = coefInGaussion * pgdLRowSum;
ngdLRowSum = coefInGaussion * ngdLRowSum;
pgdL2RowSum = pgdLRowSum * pgdLRowSum;
ngdL2RowSum = ngdLRowSum * ngdLRowSum;
pgdORowSum = coefInGaussion * pgdORowSum;
ngdORowSum = coefInGaussion * ngdORowSum;
pgdO2RowSum = pgdORowSum * pgdORowSum;
ngdO2RowSum = ngdORowSum * ngdORowSum;
/* compute {g_dL |g_dL>0 }, {g_dL |g_dL<0 },
{g_dO |g_dO>0 }, {g_dO |g_dO<0 } of each band in the line support region
first, current row belong to current band */
bandID = hID/params.widthOfBand_;
coefInGaussion = gaussCoefL_[hID%params.widthOfBand_+params.widthOfBand_];
pgdLBandSum[bandID] += coefInGaussion * pgdLRowSum;
ngdLBandSum[bandID] += coefInGaussion * ngdLRowSum;
pgdL2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdL2RowSum;
ngdL2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdL2RowSum;
pgdOBandSum[bandID] += coefInGaussion * pgdORowSum;
ngdOBandSum[bandID] += coefInGaussion * ngdORowSum;
pgdO2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdO2RowSum;
ngdO2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdO2RowSum;
/* In order to reduce boundary effect along the line gradient direction,
* a row's gradient will contribute not only to its current band, but also
* to its nearest upper and down band with gaussCoefL_.*/
bandID--;
if(bandID>=0){/* the band above the current band */
coefInGaussion = gaussCoefL_[hID%params.widthOfBand_ + 2*params.widthOfBand_];
pgdLBandSum[bandID] += coefInGaussion * pgdLRowSum;
ngdLBandSum[bandID] += coefInGaussion * ngdLRowSum;
pgdL2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdL2RowSum;
ngdL2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdL2RowSum;
pgdOBandSum[bandID] += coefInGaussion * pgdORowSum;
ngdOBandSum[bandID] += coefInGaussion * ngdORowSum;
pgdO2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdO2RowSum;
ngdO2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdO2RowSum;
}
bandID = bandID+2;
if(bandID<NUM_OF_BANDS){/*the band below the current band */
coefInGaussion = gaussCoefL_[hID%params.widthOfBand_];
pgdLBandSum[bandID] += coefInGaussion * pgdLRowSum;
ngdLBandSum[bandID] += coefInGaussion * ngdLRowSum;
pgdL2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdL2RowSum;
ngdL2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdL2RowSum;
pgdOBandSum[bandID] += coefInGaussion * pgdORowSum;
ngdOBandSum[bandID] += coefInGaussion * ngdORowSum;
pgdO2BandSum[bandID] += coefInGaussion * coefInGaussion * pgdO2RowSum;
ngdO2BandSum[bandID] += coefInGaussion * coefInGaussion * ngdO2RowSum;
}
}
/* gDLMat.Save("gDLMat.txt");
return 0; */
/* construct line descriptor */
pSingleLine->descriptor.resize(descriptorSize);
desVec = pSingleLine->descriptor.data();
short desID;
/*Note that the first and last bands only have (lengthOfLSP * widthOfBand_ * 2.0) pixels
* which are counted. */
float invN2 = 1.0/(params.widthOfBand_ * 2.0);
float invN3 = 1.0/(params.widthOfBand_ * 3.0);
float invN, temp;
for(bandID = 0; bandID<NUM_OF_BANDS; bandID++){
if(bandID==0||bandID==NUM_OF_BANDS-1){ invN = invN2;
}else{ invN = invN3;}
desID = bandID * 8;
temp = pgdLBandSum[bandID] * invN;
desVec[desID] = temp;/* mean value of pgdL; */
desVec[desID+4] = sqrt(pgdL2BandSum[bandID] * invN - temp*temp);//std value of pgdL;
temp = ngdLBandSum[bandID] * invN;
desVec[desID+1] = temp;//mean value of ngdL;
desVec[desID+5] = sqrt(ngdL2BandSum[bandID] * invN - temp*temp);//std value of ngdL;
temp = pgdOBandSum[bandID] * invN;
desVec[desID+2] = temp;//mean value of pgdO;
desVec[desID+6] = sqrt(pgdO2BandSum[bandID] * invN - temp*temp);//std value of pgdO;
temp = ngdOBandSum[bandID] * invN;
desVec[desID+3] = temp;//mean value of ngdO;
desVec[desID+7] = sqrt(ngdO2BandSum[bandID] * invN - temp*temp);//std value of ngdO;
}
// normalize;
float tempM, tempS;
tempM = 0;
tempS = 0;
desVec = pSingleLine->descriptor.data();
int base = 0;
for(short i=0; i<NUM_OF_BANDS*8; ++base, i=base*8){
tempM += *(desVec+i) * *(desVec+i);//desVec[8*i+0] * desVec[8*i+0];
tempM += *(desVec+i+1) * *(desVec+i+1);//desVec[8*i+1] * desVec[8*i+1];
tempM += *(desVec+i+2) * *(desVec+i+2);//desVec[8*i+2] * desVec[8*i+2];
tempM += *(desVec+i+3) * *(desVec+i+3);//desVec[8*i+3] * desVec[8*i+3];
tempS += *(desVec+i+4) * *(desVec+i+4);//desVec[8*i+4] * desVec[8*i+4];
tempS += *(desVec+i+5) * *(desVec+i+5);//desVec[8*i+5] * desVec[8*i+5];
tempS += *(desVec+i+6) * *(desVec+i+6);//desVec[8*i+6] * desVec[8*i+6];
tempS += *(desVec+i+7) * *(desVec+i+7);//desVec[8*i+7] * desVec[8*i+7];
}
tempM = 1/sqrt(tempM);
tempS = 1/sqrt(tempS);
desVec = pSingleLine->descriptor.data();
base = 0;
for(short i=0; i<NUM_OF_BANDS*8; ++base, i=base*8){
*(desVec+i) = *(desVec+i) * tempM;//desVec[8*i] = desVec[8*i] * tempM;
*(desVec+1+i) = *(desVec+1+i) * tempM;//desVec[8*i+1] = desVec[8*i+1] * tempM;
*(desVec+2+i) = *(desVec+2+i) * tempM;//desVec[8*i+2] = desVec[8*i+2] * tempM;
*(desVec+3+i) = *(desVec+3+i) * tempM;//desVec[8*i+3] = desVec[8*i+3] * tempM;
*(desVec+4+i) = *(desVec+4+i) * tempS;//desVec[8*i+4] = desVec[8*i+4] * tempS;
*(desVec+5+i) = *(desVec+5+i) * tempS;//desVec[8*i+5] = desVec[8*i+5] * tempS;
*(desVec+6+i) = *(desVec+6+i) * tempS;//desVec[8*i+6] = desVec[8*i+6] * tempS;
*(desVec+7+i) = *(desVec+7+i) * tempS;//desVec[8*i+7] = desVec[8*i+7] * tempS;
}
/* In order to reduce the influence of non-linear illumination,
* a threshold is used to limit the value of element in the unit feature
* vector no larger than this threshold. In Z.Wang's work, a value of 0.4 is found
* empirically to be a proper threshold.*/
desVec = pSingleLine->descriptor.data();
for(short i=0; i<descriptorSize; i++ ){
if(desVec[i]>0.4){
desVec[i]=0.4;
}
}
//re-normalize desVec;
temp = 0;
for(short i=0; i<descriptorSize; i++){
temp += desVec[i] * desVec[i];
}
temp = 1/sqrt(temp);
for(short i=0; i<descriptorSize; i++){
desVec[i] = desVec[i] * temp;
}
}/* end for(short lineIDInSameLine = 0; lineIDInSameLine<sameLineSize;
lineIDInSameLine++) */
cv::Mat appoggio = cv::Mat(1, 32, CV_32FC1);
float* pointerToRow = appoggio.ptr<float>(0);
for(int g = 0; g<32; g++)
{
/* get LBD data */
float* desVec = keyLines[lineIDInScaleVec][0].descriptor.data();
*pointerToRow = desVec[g];
pointerToRow++;
}
}/* end for(short lineIDInScaleVec = 0;
lineIDInScaleVec<numOfFinalLine; lineIDInScaleVec++) */
delete [] dL;
delete [] dO;
delete [] pgdLBandSum;
delete [] ngdLBandSum;
delete [] pgdL2BandSum;
delete [] ngdL2BandSum;
delete [] pgdOBandSum;
delete [] ngdOBandSum;
delete [] pgdO2BandSum;
delete [] ngdO2BandSum;
return 1;
}
modules/line_descriptor/src/ed_line_detector.cpp 0 → 100644
View file @ f3b4dd1f
/*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
//#define DEBUGEdgeDrawing
//#define DEBUGEDLine
using namespace std;
EDLineDetector::EDLineDetector()
{
// cout<<"Call EDLineDetector constructor function"<<endl;
//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_ = param.gradientThreshold;
anchorThreshold_ = param.anchorThreshold;
scanIntervals_ = param.scanIntervals;
minLineLen_ = param.minLineLen;
lineFitErrThreshold_ = param.lineFitErrThreshold;
InitEDLine_();
}
void EDLineDetector::InitEDLine_()
{
bValidate_ = true;
ATA = cv::Mat_<int>(2,2);
ATV = cv::Mat_<int>(1,2);
tempMatLineFit = cv::Mat_<int>(2,2);
tempVecLineFit = cv::Mat_<int>(1,2);
fitMatT = cv::Mat_<int>(2,minLineLen_);
fitVec = cv::Mat_<int>(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, bool smoothed )
{
imageWidth = image.cols;
imageHeight= image.rows;
unsigned int pixelNum = imageWidth*imageHeight;
#ifdef DEBUGEdgeDrawing
cv::imshow("prima blur", image);
cv::waitKey();
#endif
if(!smoothed){//input image hasn't been smoothed.
std::cout << "Dentro smoothed " << std::endl;
cv::Mat InImage = image.clone();
cv::GaussianBlur(InImage, image, cv::Size(ksize_,ksize_), sigma_);
}
#ifdef DEBUGEdgeDrawing
cv::imshow("dopo blur", image);
cv::waitKey();
#endif
//Is this useful? Doesn't seem to have references
//else{
// unsigned char *pImage = image.GetImageData();
// memcpy(cvImage->data.ptr, pImage, pixelNum*sizeof(unsigned char));
// }
unsigned int edgePixelArraySize = pixelNum/5;
unsigned int maxNumOfEdge = edgePixelArraySize/20;
std::cout<<"imageHeight: "<<imageHeight<<" - imageWidth:"<<imageWidth<<std::endl;
//compute dx, dy images
if(gImg_.cols!= imageWidth||gImg_.rows!= 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);
#ifdef DEBUGEdgeDrawing
cv::imshow("dxImg_", dxImg_);
cv::imshow("dyImg_", dyImg_);
cv::waitKey();
#endif
//compute gradient and direction images
// double t = (double)cv::getTickCount();
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);
// t = ((double)cv::getTickCount() - t)/cv::getTickFrequency();
// std::cout<<"FOR ABS: "<<t<<"s"<<std::endl;
short *pdxImg = dxImg_.ptr<short>();
short *pdyImg = dyImg_.ptr<short>();
short *pgImg = gImg_.ptr<short>();
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;
}
#ifdef DEBUGEdgeDrawing
cout<<"Anchor point detection, anchors.size="<<anchorsSize<<endl;
#endif
//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, lastY;
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 = pgImg[indexInArray-imageWidth+1];
gValue2 = pgImg[indexInArray+1];
gValue3 = 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 = pgImg[indexInArray-imageWidth-1];
gValue2 = pgImg[indexInArray-1];
gValue3 = 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 = pgImg[indexInArray+imageWidth+1];
gValue2 = pgImg[indexInArray+imageWidth];
gValue3 = 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 = pgImg[indexInArray-imageWidth+1];
gValue2 = pgImg[indexInArray-imageWidth];
gValue3 = 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 = pgImg[indexInArray-imageWidth+1];
gValue2 = pgImg[indexInArray+1];
gValue3 = 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 = pgImg[indexInArray-imageWidth-1];
gValue2 = pgImg[indexInArray-1];
gValue3 = 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 = pgImg[indexInArray+imageWidth+1];
gValue2 = pgImg[indexInArray+imageWidth];
gValue3 = 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 = pgImg[indexInArray- imageWidth+1];
gValue2 = pgImg[indexInArray- imageWidth];
gValue3 = 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 = pgImg[indexInArray- imageWidth+1];
gValue2 = pgImg[indexInArray+1];
gValue3 = 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 = pgImg[indexInArray- imageWidth-1];
gValue2 = pgImg[indexInArray-1];
gValue3 = 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 = pgImg[indexInArray+ imageWidth+1];
gValue2 = pgImg[indexInArray+ imageWidth];
gValue3 = 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 = pgImg[indexInArray- imageWidth+1];
gValue2 = pgImg[indexInArray- imageWidth];
gValue3 = 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 = pgImg[indexInArray- imageWidth+1];
gValue2 = pgImg[indexInArray+1];
gValue3 = 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 = pgImg[indexInArray- imageWidth -1];
gValue2 = pgImg[indexInArray-1];
gValue3 = 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 = pgImg[indexInArray+ imageWidth +1];
gValue2 = pgImg[indexInArray+ imageWidth];
gValue3 = 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 = pgImg[indexInArray- imageWidth +1];
gValue2 = pgImg[indexInArray- imageWidth];
gValue3 = 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<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;
//#ifdef DEBUGEdgeDrawing
// cout<<"Edge Drawing, numofedge = "<<numOfEdges <<endl;
// /*Show the extracted edge cvImage in color. Each chain is in different color.*/
// IplImage* cvColorImg = cvCreateImage(cvSize(imageWidth,imageHeight),IPL_DEPTH_8U, 3);
// cvSet(cvColorImg, cvScalar(0,0,0));
// CvScalar s;
// srand((unsigned)time(0));
//
// int lowest=100, highest=255;
// int range=(highest-lowest)+1;
// int r, g, b; //the color of lines
// for(unsigned int i=0; i<edgeChains.numOfEdges; i++){
// r = lowest+int(rand()%range);
// g = lowest+int(rand()%range);
// b = lowest+int(rand()%range);
// s.val[0] = r; s.val[1] = g; s.val[2] = b;
// for(indexInCors = psId[i]; indexInCors<psId[i+1]; indexInCors++){
// cvSet2D(cvColorImg,pyCors[indexInCors],pxCors[indexInCors],s);
// }
// }
//
// cvNamedWindow("EdgeColorImage", CV_WINDOW_AUTOSIZE);
// cvShowImage("EdgeColorImage", cvColorImg);
// cvWaitKey(0);
// cvReleaseImage(&cvColorImg);
//#endif
return 1;
}
int EDLineDetector::EDline(cv::Mat &image, LineChains &lines, bool smoothed)
{
//first, call EdgeDrawing function to extract edges
EdgeChains edges;
if((EdgeDrawing(image, edges, smoothed))!=true){
cout<<"Line Detection not finished"<<endl;
return -1;
}
// bValidate_ =false;
#ifdef DEBUGEDLine
cv::imshow("EdgeDrawing", image);
cv::waitKey();
#endif
//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;//the line fit error;
std::vector<double> lineEquation(2, 0);
//lineEquation.reserve(2);
lineEquations_.clear();
lineEndpoints_.clear();
lineDirection_.clear();
unsigned char *pdirImg = dirImg_.data;
unsigned int numOfLines = 0;
unsigned int offsetInEdgeArrayS, offsetInEdgeArrayE, newOffsetS;//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;//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<double> lineEqu(3,0);
//lineEqu.reserve(3);
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<float> lineEndP(4, 0);//line endpoints
//lineEndP.resize(4);
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] = a1*Px-a3*Py-a4;//x
lineEndP[1] = a2*Py-a3*Px-a5;//y
Px = pLineXCors[offsetInLineArray -1 ];//last pixel
Py = pLineYCors[offsetInLineArray -1 ];
lineEndP[2] = a1*Px-a3*Py-a4;//x
lineEndP[3] = 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<double> lineEqu(3,0);
//lineEqu.reserve(3);
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<float> lineEndP(4,0);//line endpoints
//lineEndP.reserve(4);
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] = a1*Px-a3*Py-a4;//x
lineEndP[1] = a2*Py-a3*Px-a5;//y
Px = pLineXCors[offsetInLineArray -1 ];//last pixel
Py = pLineYCors[offsetInLineArray -1 ];
lineEndP[2] = a1*Px-a3*Py-a4;//x
lineEndP[3] = 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<edges.numOfEdges; edgeID++)
// if(numOfLines < 5)
// {
// std::cout<<"LINEE MINORI DI 5"<<std::endl;
// //return -1;
// }
pLineSID[numOfLines] = offsetInLineArray;
lines.numOfLines = numOfLines;
#ifdef DEBUGEDLine
// /*Show the extracted lines in color. Each line is in different color.*/
// IplImage* cvColorImg = cvCreateImage(cvSize(imageWidth,imageHeight),IPL_DEPTH_8U, 3);
// cvSet(cvColorImg, cvScalar(0,0,0));
// CvScalar s;
// srand((unsigned)time(0));
// int lowest=100, highest=255;
// int range=(highest-lowest)+1;
// // CvPoint point;
// // CvFont font;
// // cvInitFont(&font,CV_FONT_HERSHEY_SIMPLEX ,1.0,1.0,0,1);
// int r, g, b; //the color of lines
// for(unsigned int i=0; i<lines.numOfLines; i++){
// r = lowest+int(rand()%range);
// g = lowest+int(rand()%range);
// b = lowest+int(rand()%range);
// s.val[0] = b; s.val[1] = g; s.val[2] = r;
// for(offsetInLineArray = pLineSID[i]; offsetInLineArray<pLineSID[i+1]; offsetInLineArray++){
// cvSet2D(cvColorImg,pLineYCors[offsetInLineArray],pLineXCors[offsetInLineArray],s);
// }
// // iter = lines[i].begin();
// // point = cvPoint(iter->x,iter->y);
// // char buf[10];
// // sprintf( buf, "%d ", i);
// // cvPutText(cvColorImg,buf,point,&font,CV_RGB(r,g,b));
// }
// cvNamedWindow("LineColorImage", CV_WINDOW_AUTOSIZE);
// cvShowImage("LineColorImage", cvColorImg);
// cvWaitKey(0);
// cvReleaseImage(&cvColorImg);
#endif
return 1;
}
double EDLineDetector::LeastSquaresLineFit_( unsigned int *xCors, unsigned int *yCors,
unsigned int offsetS, std::vector<double> &lineEquation)
{
// if(lineEquation.size()!=2){
// //std::cout<<"SHOULD NOT BE != 2"<<std::endl;
// //CV_Assert(false);
// }
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<float>();//fitMatT = [x0, x1, ... xn; 1,1,...,1];
for(int i=0; i<minLineLen_; i++){
//*(pMatT+minLineLen_) = 1; //the value are not changed;
*(pMatT++) = xCors[offsetS];
fitVec[0][i] = yCors[offsetS++];
}
ATA = fitMatT * fitMatT.t();
ATV = fitMatT * fitVec.t();
/* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */
pATA = ATA.ptr<float>();
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<float>();//fitMatT = [y0, y1, ... yn; 1,1,...,1];
for(int i=0; i<minLineLen_; i++){
//*(pMatT+minLineLen_) = 1;//the value are not changed;
*(pMatT++) = yCors[offsetS];
fitVec[0][i] = xCors[offsetS++];
}
ATA = fitMatT*(fitMatT.t());
ATV = fitMatT * fitVec.t();
/* [a,b]^T = Inv(mat^T * mat) * mat^T * vec */
pATA = ATA.ptr<float>();
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<double> &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_Assert(false);
}
cv::Mat_<float> matT(2,newLength);
cv::Mat_<float> vec(newLength,1);
float * pMatT;
float * pATA;
// double fitError = 0;
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<float>();//matT = [x0', x1', ... xn'; 1,1,...,1]
for(int i=0; i<newLength; i++){
*(pMatT+newLength) = 1;
*(pMatT++) = xCors[newOffsetS];
vec[0][i] = 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<float>();
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]));
/*compute line fit error */
// for(int i=0; i<length; i++){
// coef = double(yCors[offsetS]) - double(xCors[offsetS++]) * lineEquation[0] - lineEquation[1];
// fitError += coef*coef;
// }
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.GetData();//matT = [y0', y1', ... yn'; 1,1,...,1]
pMatT = matT.ptr<float>();//matT = [y0', y1', ... yn'; 1,1,...,1]
for(int i=0; i<newLength; i++){
*(pMatT+newLength) = 1;
*(pMatT++) = yCors[newOffsetS];
vec[0][i] = 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<float>();
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]));
/*compute line fit error */
// for(int i=0; i<length; i++){
// coef = double(xCors[offsetS]) - double(yCors[offsetS++]) * lineEquation[0] - lineEquation[1];
// fitError += coef*coef;
// }
}
return 0;
}
bool EDLineDetector::LineValidation_( unsigned int *xCors, unsigned int *yCors,
unsigned int offsetS, unsigned int offsetE,
std::vector<double> &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>();
short *pdyImg = dyImg_.ptr<short>();
double dx, dy;
std::vector<double> pointDirection;
int index;
for(int i=0; i<n; i++){
//index = yCors[offsetS]*imageWidth + xCors[offsetS++];
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 = atan2(-dy,dx);//line direction is in fourth quadrant
}
if(meanGradientX<=0&&meanGradientY>0){//second quadrant, and positive direction of Y axis.
direction = atan2(dy,dx);//line direction is in first quadrant
}
if(meanGradientX<0&&meanGradientY<=0){//third quadrant, and negative direction of X axis.
direction = atan2(dy,-dx);//line direction is in second quadrant
}
if(meanGradientX>=0&&meanGradientY<0){//fourth quadrant, and negative direction of Y axis.
direction = 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, bool smoothed)
{
if((EDline(image, lines_, smoothed)) != true){
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;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment