Skip to content

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

Original code + comments

parent 079ff5c0
Show whitespace changes
Inline Side-by-side
Showing with 3368 additions and 0 deletions
modules/line_descriptor/CMakeLists.txt 0 → 100644
View file @ f2937652
#INCLUDE_DIRECTORIES(/data1/lz/LilianTests/ARPACKLAB/arpack++/include/)
#INCLUDE_DIRECTORIES(/data1/lz/LilianTests/ARPACKLAB/arpack++)
INCLUDE_DIRECTORIES(/usr/include/arpack++/include/)
INCLUDE_DIRECTORIES(/usr/include/arpack++)
#project name
PROJECT(LILIANTESTS)
cmake_minimum_required(VERSION 2.8)
if(COMMAND cmake_policy)
cmake_policy(SET CMP0003 NEW)
endif(COMMAND cmake_policy)
SET(BUILD_SHARED_LIBS ON)
## where are user-specific cmake modules
SET(CMAKE_MODULE_PATH $ENV{CMAKE_MODULE_PATH})
OPTION(USE_BIAS OFF)
#IF(USE_BIAS)
#FIND_PACKAGE(BIAS)
#IF(BIAS_FOUND)
# INCLUDE(${BIAS_USE_FILE})
#ELSE(BIAS_FOUND)
# MESSAGE(SEND_ERROR "BIAS lib not found.")
#ENDIF(BIAS_FOUND)
#ENDIF(USE_BIAS)
IF(USE_BIAS)
INCLUDE_DIRECTORIES(${BIAS_INCLUDE_DIR} /usr/local/include/BIAS /usr/include/ImageMagick /usr/local/include/opencv2 ${WXWIDGETS_INCLUDE_DIR})
LINK_DIRECTORIES( ${BIAS_LINK_DIRECTORIES})
ENDIF(USE_BIAS)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIAS_CXX_FLAGS} -std=c++0x")
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${BIAS_C_FLAGS}")
# source files of library "LineMatchingLib" to be created
SET(LineMatchingLib_SRCS
#src/PairwiseLineMatching.cpp
src/LineDescriptor.cpp
src/EDLineDetector.cpp
)
# header files to be installed
SET(LineMatchingLib_HEADER
#include/PairwiseLineMatching.hh
include/LineDescriptor.hh
include/EDLineDetector.hh
include/LineStructure.hh
)
ADD_LIBRARY(LineMatchingLib
${LineMatchingLib_SRCS}
${LineMatchingLib_HEADER})
TARGET_LINK_LIBRARIES(LineMatchingLib /usr/lib/libsuperlu.so opencv_core opencv_highgui opencv_imgproc)
#ADD_EXECUTABLE(TestLineMatchingAlgorithm src/TestLineMatchingAlgorithm.cpp)
#TARGET_LINK_LIBRARIES(TestLineMatchingAlgorithm LineMatchingLib )
modules/line_descriptor/include/opencv2/line_descriptor/EDLineDetector.hh 0 → 100644
View file @ f2937652
/*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.
*/
#ifndef EDLINEDETECTOR_HH_
#define EDLINEDETECTOR_HH_
#include <vector>
#include <iostream>
#include <list>
#include <opencv2/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/utility.hpp>
#include <array>
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 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
struct EDLineParam{
int ksize;
float sigma;
float gradientThreshold;
float anchorThreshold;
int scanIntervals;
int minLineLen;
double lineFitErrThreshold;
};
#define RELATIVE_ERROR_FACTOR 100.0
#define M_LN10 2.30258509299404568402
#define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
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();
}
/* This class is used to detect lines from input image.
* First, edges are extracted from input image following the method presented in Cihan Topal and
* Cuneyt Akinlar's paper:"Edge Drawing: A Heuristic Approach to Robust Real-Time Edge Detection", 2010.
* Then, lines are extracted from the edge image following the method presented in Cuneyt Akinlar and
* Cihan Topal's paper:"EDLines: A real-time line segment detector with a false detection control", 2011
* PS: The linking step of edge detection has a little bit difference with the Edge drawing algorithm
* described in the paper. The edge chain doesn't stop when the pixel direction is changed.
*/
class EDLineDetector
{
public:
EDLineDetector();
EDLineDetector(EDLineParam param);
~EDLineDetector();
/*extract edges from image
*image: In, gray image;
*edges: Out, store the edges, each edge is a pixel chain
*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::array<double, 3> > lineEquations_;
//store the line endpoints, [x1,y1,x2,y3]
std::vector<std::array<float, 4> > 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_;
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 array;
*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::array<double, 2> &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 array;
*newOffsetS: In, start index of extended part;
*offsetE:In, end index of this chain in array;
*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::array<double, 2> &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::array<double, 3> &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*/
unsigned int *pFirstPartEdgeX_;//store the X coordinates of the first part of the pixels for chains
unsigned int *pFirstPartEdgeY_;//store the Y coordinates of the first part of the pixels for chains
unsigned int *pFirstPartEdgeS_;//store the start index of every edge chain in the first part arrays
unsigned int *pSecondPartEdgeX_;//store the X coordinates of the second part of the pixels for chains
unsigned int *pSecondPartEdgeY_;//store the Y coordinates of the second part of the pixels for chains
unsigned int *pSecondPartEdgeS_;//store the start index of every edge chain in the second part arrays
unsigned int *pAnchorX_;//store the X coordinates of anchors
unsigned int *pAnchorY_;//store the Y coordinates of anchors
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/include/opencv2/line_descriptor/LineDescriptor.hh 0 → 100644
View file @ f2937652
/*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.
*/
#ifndef LINEDESCRIPTOR_HH_
#define LINEDESCRIPTOR_HH_
#include "EDLineDetector.hh"
#include "LineStructure.hh"
#include <map>
struct OctaveLine{
unsigned int octaveCount;//the octave which this line is detected
unsigned int lineIDInOctave;//the line ID in that octave image
unsigned int lineIDInScaleLineVec;//the line ID in Scale line vector
float lineLength; //the length of line in original image scale
};
/* This class is used to generate the line descriptors from multi-scale images */
class LineDescriptor
{
public:
LineDescriptor(int n_octave);
LineDescriptor(unsigned int numOfBand, unsigned int widthOfBand, int n_octave);
~LineDescriptor();
enum{
NearestNeighbor=0, //the nearest neighbor is taken as matching
NNDR=1//nearest/next ratio
};
/*This function is used to detect lines from multi-scale images.*/
int OctaveKeyLines(cv::Mat & image, ScaleLines &keyLines);
int GetLineDescriptor(cv::Mat & image,
ScaleLines &keyLines);
int GetLineDescriptor(cv::Mat & image,
ScaleLines &keyLines, bool lsd);
void findNearestParallelLines(ScaleLines & keyLines);
void GetLineBinaryDescriptor(std::vector<cv::Mat> & oct_binaryDescMat, ScaleLines & keyLines);
int MatchLineByDescriptor(ScaleLines &keyLinesLeft, ScaleLines &keyLinesRight,
std::vector<short> &matchLeft, std::vector<short> &matchRight,
int criteria=NNDR);
float LowestThreshold;//global threshold for line descriptor distance, default is 0.35
float NNDRThreshold;//the NNDR threshold for line descriptor distance, default is 0.6
int ksize_; //the size of Gaussian kernel: ksize X ksize, default value is 5.
unsigned int numOfOctave_;//the number of image octave
unsigned int numOfBand_;//the number of band used to compute line descriptor
unsigned int widthOfBand_;//the width of band;
std::vector<float> gaussCoefL_;//the local gaussian coefficient apply to the orthogonal line direction within each band;
std::vector<float> gaussCoefG_;//the global gaussian coefficient apply to each Row within line support region
private:
void sample(float *igray,float *ogray, float factor, int width, int height)
{
int swidth = (int)((float) width / factor);
int sheight = (int)((float) height / factor);
for(int j=0; j < sheight; j++)
for(int i=0; i < swidth; i++)
ogray[j*swidth + i] = igray[(int)((float) j * factor) * width + (int) ((float) i*factor)];
}
void sampleUchar(uchar *igray,uchar *ogray, float factor, int width, int height)
{
int swidth = (int)((float) width / factor);
int sheight = (int)((float) height / factor);
for(int j=0; j < sheight; j++)
for(int i=0; i < swidth; i++)
ogray[j*swidth + i] = igray[(int)((float) j * factor) * width + (int) ((float) i*factor)];
}
/*Compute the line descriptor of input line set. This function should be called
*after OctaveKeyLines() function; */
int ComputeLBD_(ScaleLines &keyLines);
/*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_;
};
#endif /* LINEDESCRIPTOR_HH_ */
modules/line_descriptor/include/opencv2/line_descriptor/LineStructure.hh 0 → 100644
View file @ f2937652
/*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.
*/
#ifndef LINESTRUCTURE_HH_
#define LINESTRUCTURE_HH_
#include <vector>
// A 2D line (normal equation parameters).
struct SingleLine
{
//note: rho and theta are based on coordinate origin, i.e. the top-left corner of image
double rho;//unit: pixel length
double theta;//unit: rad
double linePointX;// = rho * cos(theta);
double linePointY;// = rho * sin(theta);
//for EndPoints, the coordinate origin is the top-left corner of image.
double startPointX;
double startPointY;
double endPointX;
double endPointY;
//direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis.
double direction;
//mean gradient magnitude
double gradientMagnitude;
//mean gray value of pixels in dark side of line
double darkSideGrayValue;
//mean gray value of pixels in light side of line
double lightSideGrayValue;
//the length of line
double lineLength;
//the width of line;
double width;
//number of pixels
int numOfPixels;
//the decriptor of line
std::vector<double> descriptor;
};
// Specifies a vector of lines.
typedef std::vector<SingleLine> Lines_list;
struct OctaveSingleLine
{
/*endPoints, the coordinate origin is the top-left corner of the original image.
*startPointX = sPointInOctaveX * (factor)^octaveCount; */
float startPointX;
float startPointY;
float endPointX;
float endPointY;
//endPoints, the coordinate origin is the top-left corner of the octave image.
float sPointInOctaveX;
float sPointInOctaveY;
float ePointInOctaveX;
float ePointInOctaveY;
//direction of a line, the angle between positive line direction (dark side is in the left) and positive X axis.
float direction;
//the summation of gradient magnitudes of pixels on lines
float salience;
//the length of line
float lineLength;
//number of pixels
unsigned int numOfPixels;
//the octave which this line is detected
unsigned int octaveCount;
//the decriptor of line
std::vector<float> descriptor;
};
// Specifies a vector of lines.
typedef std::vector<OctaveSingleLine> LinesVec;
typedef std::vector<LinesVec> ScaleLines;//each element in ScaleLines is a vector of lines which corresponds the same line detected in different octave images.
#endif /* LINESTRUCTURE_HH_ */
modules/line_descriptor/src/EDLineDetector.cpp 0 → 100644
View file @ f2937652
/*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 "../include/EDLineDetector.hh"
#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_ = 60; // ***** 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);
#endif
if(!smoothed){//input image hasn't been smoothed.
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;
//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::array<double, 2> lineEquation;
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::array<double, 3> lineEqu = {lineEquation[0]*coef1, -1*coef1, 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::array<float, 4> lineEndP;//line endpoints
double a1 = lineEqu[1]*lineEqu[1];
double a2 = lineEqu[0]*lineEqu[0];
double a3 = lineEqu[0]*lineEqu[1];
double a4 = lineEqu[2]*lineEqu[0];
double a5 = lineEqu[2]*lineEqu[1];
unsigned int Px = pLineXCors[pLineSID[numOfLines] ];//first pixel
unsigned int Py = pLineYCors[pLineSID[numOfLines] ];
lineEndP[0] = 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::array<double, 3> lineEqu= {1*coef1,-lineEquation[0]*coef1, -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::array<float, 4> lineEndP;//line endpoints
double a1 = lineEqu[1]*lineEqu[1];
double a2 = lineEqu[0]*lineEqu[0];
double a3 = lineEqu[0]*lineEqu[1];
double a4 = lineEqu[2]*lineEqu[0];
double a5 = lineEqu[2]*lineEqu[1];
unsigned int Px = pLineXCors[pLineSID[numOfLines] ];//first pixel
unsigned int Py = pLineYCors[pLineSID[numOfLines] ];
lineEndP[0] = 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++)
std::cout<<"numOfLines: "<<numOfLines<<std::endl;
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::array<double, 2> &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];
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];
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::array<double, 2> &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::array<double, 3> &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++];
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;
}
modules/line_descriptor/src/LineDescriptor.cpp 0 → 100644
View file @ f2937652
/*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 "../include/LineDescriptor.hh"
#define SalienceScale 0.9//0.9
//test pairs (0,7),(0,8),(1,7),(1,8) are excluded to get a descriptor size of 256 bit
static const int combinations [32][2] = {{0,1},{0,2},{0,3},{0,4},{0,5},{0,6},/*{0,7},{0,8},*/{1,2},{1,3},{1,4},{1,5},{1,6},/*{1,7},{1,8},*/{2,3},{2,4},{2,5},{2,6},{2,7},{2,8},{3,4},{3,5},{3,6},{3,7},{3,8},{4,5},{4,6},{4,7},{4,8},{5,6},{5,7},{5,8},{6,7},{6,8},{7,8}};
static inline int get2Pow(int i) {
switch (i) {
case 0:
return 1;
case 1:
return 2;
case 2:
return 4;
case 3:
return 8;
case 4:
return 16;
case 5:
return 32;
case 6:
return 64;
case 7:
return 128;
default:
std::cout<<"Invalid pow"<<std::endl;
CV_Assert(false);
return -1; //avoid warning
}
}
//#define DEBUGLinesInOctaveImages
using namespace std;
LineDescriptor::LineDescriptor(int n_octave)
{
/* generate a random number by using time as a seed */
srand ( time(NULL) );
/* set gaussian kernel's size */
ksize_ = 5;
/* set number of octaves (default = 5) */
numOfOctave_ = n_octave;//5
/* create as many LineVecs as the number of octaves */
edLineVec_.resize(numOfOctave_);
for(unsigned int i=0; i<numOfOctave_; i++){
edLineVec_[i] = new EDLineDetector;
}
/* set number of bands */
numOfBand_ = 9; // suggested by authors
/* set bands' width */
widthOfBand_ = 7; //widthOfBand_%3 must equal to 0; 7 is a good value.
/* prepare a vector to host local weights F_l*/
gaussCoefL_.resize(widthOfBand_*3);
/* compute center of central band (every computation involves 2-3 bands) */
double u = (widthOfBand_*3-1)/2;
/* compute exponential part of F_l */
double sigma = (widthOfBand_*2+1)/2;// (widthOfBand_*2+1)/2;
double invsigma2 = -1/(2*sigma*sigma);
/* compute all local weights */
double dis;
for(int i=0; i<widthOfBand_*3; i++){
dis = i-u;
gaussCoefL_[i] = exp(dis*dis*invsigma2);
}
/* prepare a vector for global weights F_g*/
gaussCoefG_.resize(numOfBand_*widthOfBand_);
/* compute center of LSR */
u = (numOfBand_*widthOfBand_-1)/2;
/* compute exponential part of F_g */
sigma = u;
invsigma2 = -1/(2*sigma*sigma);
for(int i=0; i<numOfBand_*widthOfBand_; i++){
dis = i-u;
gaussCoefG_[i] = exp(dis*dis*invsigma2);
}
/* THRESHOLDS
2 is used to show recall ratio;
0.2 is used to show scale space results;
0.35 is used when verify geometric constraints */
LowestThreshold = 0.3;
NNDRThreshold = 0.6;
}
LineDescriptor::LineDescriptor(unsigned int numOfBand, unsigned int widthOfBand, int n_octave){
/* generate a random number by using time as a seed */
srand ( time(NULL) );
/* set gaussian kernel's size */
ksize_ = 5;
/* set number of octaves (default = 5) */
numOfOctave_ = n_octave;
/* create as many LineVecs as the number of octaves */
edLineVec_.resize(numOfOctave_);
for(unsigned int i=0; i<numOfOctave_; i++){
edLineVec_[i] = new EDLineDetector;
}
/* set number of bands */
numOfBand_ = numOfBand;
/* set bands' width */
widthOfBand_ = widthOfBand;
/* prepare a vector to host local weights F_l*/
gaussCoefL_.resize(widthOfBand_*3);
/* compute center of central band (every computation involves 2-3 bands) */
double u = (widthOfBand_*3-1)/2;
/* compute exponential part of F_l */
double sigma = (widthOfBand_*2+1)/2;// (widthOfBand_*2+1)/2;
double invsigma2 = -1/(2*sigma*sigma);
/* compute all local weights */
double dis;
for(int i=0; i<widthOfBand_*3; i++){
dis = i-u;
gaussCoefL_[i] = exp(dis*dis*invsigma2);
}
/* prepare a vector for global weights F_g */
gaussCoefG_.resize(numOfBand_*widthOfBand_);
/* compute center of LSR */
u = (numOfBand_*widthOfBand_-1)/2;
/* compute exponential part of F_g */
sigma = u;
invsigma2 = -1/(2*sigma*sigma);
for(int i=0; i<numOfBand_*widthOfBand_; i++){
dis = i-u;
gaussCoefG_[i] = exp(dis*dis*invsigma2);
}
/* THRESHOLDS
2 is used to show recall ratio;
0.2 is used to show scale space results;
0.35 is used when verify geometric constraints */
LowestThreshold = 0.35;//0.35;
NNDRThreshold = 0.2;//0.6
}
/* destructor */
LineDescriptor::~LineDescriptor(){
for(unsigned int i=0; i<numOfOctave_; i++){
if(edLineVec_[i] !=NULL){
delete edLineVec_[i];
}
}
}
/*Line detection method: element in keyLines[i] includes a set of lines which is the same line
* detected in different octave images.
*/
int LineDescriptor::OctaveKeyLines(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(unsigned int octaveCount = 0; octaveCount<numOfOctave_; octaveCount++){
/* matrix storing results from blurring processes */
cv::Mat blur;
/* Compute sigma value for Gaussian filter for each level by adding incremental blur from previous level.
* curSigma = [sqrt(2)]^octaveCount;
* increaseSigma^2 = curSigma^2 - preSigma^2 */
float increaseSigma = sqrt(curSigma2-preSigma2);
/* apply Gaussian blur */
cv::GaussianBlur(image, blur, cv::Size(ksize_,ksize_), increaseSigma);
/* 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[numOfOctave_];
scale[0] = 1;
for(unsigned int octaveCount = 1; octaveCount<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(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(unsigned int octaveCount = 1; octaveCount<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(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;
}
/*The definitions of line descriptor,mean values of {g_dL>0},{g_dL<0},{g_dO>0},{g_dO<0} of each row in band
*and std values of sum{g_dL>0},sum{g_dL<0},sum{g_dO>0},sum{g_dO<0} of each row in band.
* With overlap region. */
int LineDescriptor::ComputeLBD_(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 = widthOfBand_*numOfBand_;//the height of line support region;
short descriptorSize = numOfBand_ * 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[numOfBand_];//the summation of {g_dL |g_dL>0 } for each band of the region;
float *ngdLBandSum = new float[numOfBand_];//the summation of {g_dL |g_dL<0 } for each band of the region;
float *pgdL2BandSum = new float[numOfBand_];//the summation of {g_dL^2 |g_dL>0 } for each band of the region;
float *ngdL2BandSum = new float[numOfBand_];//the summation of {g_dL^2 |g_dL<0 } for each band of the region;
float *pgdOBandSum = new float[numOfBand_];//the summation of {g_dO |g_dO>0 } for each band of the region;
float *ngdOBandSum = new float[numOfBand_];//the summation of {g_dO |g_dO<0 } for each band of the region;
float *pgdO2BandSum = new float[numOfBand_];//the summation of {g_dO^2 |g_dO>0 } for each band of the region;
float *ngdO2BandSum = new float[numOfBand_];//the summation of {g_dO^2 |g_dO<0 } for each band of the region;
short numOfBitsBand = numOfBand_*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/widthOfBand_;
coefInGaussion = gaussCoefL_[hID%widthOfBand_+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%widthOfBand_ + 2*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<numOfBand_){/*the band below the current band */
coefInGaussion = gaussCoefL_[hID%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/(widthOfBand_ * 2.0);
float invN3 = 1.0/(widthOfBand_ * 3.0);
float invN, temp;
for(bandID = 0; bandID<numOfBand_; bandID++){
if(bandID==0||bandID==numOfBand_-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<numOfBand_*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<numOfBand_*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++) */
}/* 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;
}
unsigned char binaryTest(float* f1, float*f2)
{
uchar result=0;
for(int i = 0; i<8; i++)
{
// std::cout<<"f1[: "<<i<<"]: "<<f1[i];
// std::cout<<" -- f2[: "<<i<<"]: "<<f2[i]<<std::endl;
if(f1[i]>f2[i])
{
// std::cout<< " ------ 1 ------- "<<std::endl;
result+=get2Pow(i);
}
else
{
// std::cout<< " ------ 0 ------- "<<std::endl;
}
}
return result;
}
//unsigned char binaryTest2(float* f1, float*f2)
//{
// uchar result=0;
// for(int i = 0; i<8; i++)
// {
//// std::cout<<"f1[: "<<i<<"]: "<<f1[i];
//// std::cout<<" -- f2[: "<<i<<"]: "<<f2[i]<<std::endl;
// if(f1[i]>f2[i])
// {
//// std::cout<< " ------ 1 ------- "<<std::endl;
// result+=get2Pow(i);
// }
// else
// {
//// std::cout<< " ------ 0 ------- "<<std::endl;
// }
// }
//
// return result;
//
//}
unsigned char binaryIndexTest(int f1, int f2, std::vector<float>& desc)
{
uchar result=0;
for(int i = 0; i<8; i++)
{
// std::cout<<"f1[: "<<i<<"]: "<<f1[i];
// std::cout<<" -- f2[: "<<i<<"]: "<<f2[i]<<std::endl;
if(desc[i+f1]>desc[i+f2])
{
// std::cout<< " ------ 1 ------- "<<std::endl;
result+=get2Pow(i);
}
else
{
// std::cout<< " ------ 0 ------- "<<std::endl;
}
}
return result;
}
/* create a vector of matrices: every matrix represents i-th octave and it has a (binary) descriptor
of one of lines extracted from i-th octave on each row */
void LineDescriptor::GetLineBinaryDescriptor(std::vector<cv::Mat> & oct_binaryDescMat, ScaleLines & keyLines)
{
/* std::cout<<"numOfOctave: "<<numOfOctave_<<std::endl; */
/* create an int vector, whose i-th position stores how many lines where extracted from
i-th octave.
At the beginning, initialize vector as a sequence of counters set to 0 */
std::vector<int> rows_size;
for(int i = 0; i<numOfOctave_; i++)
rows_size.push_back(0);
/* modify counters in rows_size in order to reflect
the number of lines extracted from each octave */
for(int offsetKeyLines = 0; offsetKeyLines<keyLines.size(); offsetKeyLines++)
{
for(int offsetOctave = 0; offsetOctave<keyLines[offsetKeyLines].size(); offsetOctave++)
{
rows_size[keyLines[offsetKeyLines][offsetOctave].octaveCount]++;
}
}
/* prepare a vector of pointers to the first rows of matrices
that are going to be created in next loop */
std::vector<uchar *> vec_binaryMat_p;
/* loop on the number of the octaves */
for(int i = 0; i<numOfOctave_; i++)
{
/* create a matrix having as many rows as the number of lines
extracted from i-th octave and 32 columns (each row stores
a 32 bit string) */
cv::Mat mat_binary(rows_size[i], 32, CV_8UC1);
/* add created matrix to list passed as an argument */
oct_binaryDescMat.push_back(mat_binary);
/* store a pointer to the first row of the matrix that has been
just created */
vec_binaryMat_p.push_back(oct_binaryDescMat.at(i).ptr());
}
/* loop over the number of LineVecs */
for(int offsetKeyLines = 0; offsetKeyLines<keyLines.size(); offsetKeyLines++)
{
/* loop over the number of lines inside each LineVec */
for(int offsetOctave = 0; offsetOctave<keyLines[offsetKeyLines].size(); offsetOctave++)
{
/* binaryMat_p is a pointer to a pointer,
because we must increment the content of a vector,
in such a way that pointers are updated with the right row of matrix */
uchar ** binaryMat_p = &vec_binaryMat_p.at(keyLines[offsetKeyLines][offsetOctave].octaveCount);
/* get descriptor associated to i-th line (as a sequence of floats) */
float * desVec = keyLines[offsetKeyLines][offsetOctave].descriptor.data();
for(int comb = 0; comb < 32; comb++)
{
*(*binaryMat_p) = binaryTest(&desVec[8*combinations[comb][0]], &desVec[8*combinations[comb][1]]);
(*binaryMat_p)++;
}
}
}
//writeMat(oct_binaryDescMat[0], "binaryMat_deb",0);
}
/* compute LBD descriptors */
int LineDescriptor::GetLineDescriptor(cv::Mat & image, ScaleLines & keyLines)
{
/*check whether image depth is different from 0 */
if(image.depth() != 0)
{
std::cout << "Warning, depth image!= 0" << std::endl;
CV_Assert(false);
}
/* get clock's TIC */
double t = (double)cv::getTickCount();
/* compute LineVecs extraction and in case of failure, return an error code */
if((OctaveKeyLines(image,keyLines))!=true){
cout << "OctaveKeyLines failed" << endl;
return -1;
}
/* get time lapse */
t = ((double)cv::getTickCount() - t)/cv::getTickFrequency();
std::cout << "Time of line extraction: "<< t << "s" <<std::endl;
// for(int j = 0; j<keyLines.size(); j++)
// {
// for(int k = 0; k<keyLines[j].size(); k++)
// {
// OctaveSingleLine singleLine = keyLines[j][k];
// std::cout<<"-----------["<<j<<"]["<<k<<"]--------------"<<std::endl;
// std::cout<<"singleLine.octaveCount :"<<singleLine.octaveCount<<std::endl;
// std::cout<<"singleLine.direction :"<<singleLine.direction<<std::endl;
// std::cout<<"singleLine.lineLength :"<<singleLine.lineLength<<std::endl;
// std::cout<<"singleLine.salience :"<<singleLine.salience<<std::endl;
// std::cout<<"singleLine.numOfPixels :"<<singleLine.numOfPixels<<std::endl;
// std::cout<<"singleLine.sPointInOctaveX :"<<singleLine.sPointInOctaveX<<std::endl;
// std::cout<<"singleLine.sPointInOctaveY :"<<singleLine.sPointInOctaveY<<std::endl;
// std::cout<<"singleLine.ePointInOctaveX :"<<singleLine.ePointInOctaveX<<std::endl;
// std::cout<<"singleLine.ePointInOctaveY :"<<singleLine.ePointInOctaveY<<std::endl;
// std::cout<<"singleLine.startPointX :"<<singleLine.startPointX<<std::endl;
// std::cout<<"singleLine.startPointY :"<<singleLine.startPointY<<std::endl;
// std::cout<<"singleLine.endPointX :"<<singleLine.endPointX<<std::endl;
// std::cout<<"singleLine.endPointY :"<<singleLine.endPointY<<std::endl;
// std::cout<<"--------------------------------"<<std::endl;
// }
// }
// t = (double)cv::getTickCount();
/* compute LBD descriptors */
ComputeLBD_(keyLines);
// t = ((double)cv::getTickCount() - t)/cv::getTickFrequency();
// std::cout<<"time descriptor extraction: "<<t<<"s"<<std::endl;
//
// for(int j = 0; j<keyLines.size(); j++)
// {
// for(int k = 0; k<keyLines[j].size(); k++)
// {
// for(int i = 0; i<keyLines[j][k].descriptor.size(); i++)
// std::cout<<"keylines["<<j<<"]["<<k<<"].descriptor["<<i<<"]: "<<keyLines[j][k].descriptor[i]<<std::endl;
// }
// }
// srand((unsigned)time(0));
// int lowest=100, highest=25555;
// int range=(highest-lowest)+1;
// int r;
// r = lowest+int(rand()%range);
// string folder_debug = "/home/manuele/src/imgretrieval/mat_debug/";
// std::cout<<"LINDESCRIPTOR INIZIO"<<std::endl;
// cv::Mat desc_float(keyLines.size(), keyLines[0][0].descriptor.size(), CV_32FC1);
// for(int i = 0; i<keyLines.size(); i++)
// {
// for(int j = 0; j<keyLines[i][0].descriptor.size(); j++)
// {
// //std::cout<<keyLines[i][0].descriptor[j];
// desc_float.at<float>(i, j) = (float)keyLines[i][0].descriptor[j];
// }
// std::cout<<std::endl;
// }
//// std::cout<<desc_float<<std::endl;
// std::cout<<"LINDESCRIPTOR FINE"<<std::endl<<std::endl;
// writeMat(desc_float, folder_debug+"linedesc", r);
return 1;
}
/* check whether two lines are parallel, using their directions */
bool areParallels(float direction1, float direction2)
{
if(abs(abs(direction1) - abs(direction2)) <= 0.02)
return true;
if(abs(direction1) + abs(direction2) >= 3.12)
return true;
return false;
}
void LineDescriptor::findNearestParallelLines(ScaleLines & keyLines)
{
std::cout<<"PARALLELLINES: size: "<<keyLines.size()<<std::endl;
std::map<float, OctaveSingleLine> parallels;
/* loop over LineVecs */
for(int j = 0; j<keyLines.size(); j++)
{
/* loop over lines in current LineVec */
for(int k = 0; k<keyLines[j].size(); k++)
{
/* get current line */
OctaveSingleLine singleLine = keyLines[j][k];
/* get an iterator to map of lines */
std::map<float, OctaveSingleLine>::iterator it;
/* scan map to searck for a line parallel to current one */
bool foundParallel = false;
for(it = parallels.begin(); it != parallels.end(); it++) {
if(!areParallels(it->first, singleLine.direction))
{
foundParallel = true;
break;
}
}
/* if a parallel line has not been found, add current line
to map, using its direction as a key */
if(!foundParallel)
parallels[singleLine.direction] = singleLine;
}
}
/* create a vector of LineVecs, each one containing a line that is
not parallel to any other one inside a different LineVec */
ScaleLines newKeyLines;
std::map<float, OctaveSingleLine>::iterator it;
for(it = parallels.begin(); it != parallels.end(); it++) {
LinesVec lineScaleLs;
lineScaleLs.push_back(it->second);
newKeyLines.push_back(lineScaleLs);
}
keyLines = newKeyLines;
}
int LineDescriptor::GetLineDescriptor(cv::Mat & image, ScaleLines & keyLines, bool lsd)
{
/*check whether image depth is different from 0 */
if(image.depth() != 0)
{
std::cout << "Warning, depth image!= 0" << std::endl;
CV_Assert(false);
}
/* make a clone of input image */
cv::Mat original = image.clone();
/* get clock's TIC */
double t = (double)cv::getTickCount();
/* compute LineVecs extraction and in case of failure, return an error code */
if((OctaveKeyLines(image,keyLines))!=true){
cout << "OctaveKeyLines failed" << endl;
return -1;
}
/* get time lapse */
t = ((double)cv::getTickCount() - t)/cv::getTickFrequency();
std::cout << "time line extraction: " << t<< "s" <<std::endl;
// for(int j = 0; j<keyLines.size(); j++)
// {
// for(int k = 0; k<keyLines[j].size(); k++)
// {
// OctaveSingleLine singleLine = keyLines[j][k];
// std::cout<<"-----------["<<j<<"]["<<k<<"]--------------"<<std::endl;
//// std::cout<<"singleLine.octaveCount :"<<singleLine.octaveCount<<std::endl;
// std::cout<<"singleLine.direction :"<<singleLine.direction<<" -- "<<atan2((singleLine.startPointY - singleLine.endPointY),(singleLine.startPointX - singleLine.endPointX))<<" --- "<<atan2((singleLine.endPointY - singleLine.startPointY),(singleLine.endPointX - singleLine.startPointX))<<std::endl;
//// std::cout<<"singleLine.lineLength :"<<singleLine.lineLength<<std::endl;
//// std::cout<<"singleLine.salience :"<<singleLine.salience<<std::endl;
//// std::cout<<"singleLine.numOfPixels :"<<singleLine.numOfPixels<<std::endl;
//// std::cout<<"singleLine.sPointInOctaveX :"<<singleLine.sPointInOctaveX<<std::endl;
//// std::cout<<"singleLine.sPointInOctaveY :"<<singleLine.sPointInOctaveY<<std::endl;
//// std::cout<<"singleLine.ePointInOctaveX :"<<singleLine.ePointInOctaveX<<std::endl;
//// std::cout<<"singleLine.ePointInOctaveY :"<<singleLine.ePointInOctaveY<<std::endl;
//// std::cout<<"singleLine.startPointX :"<<singleLine.startPointX<<std::endl;
//// std::cout<<"singleLine.startPointY :"<<singleLine.startPointY<<std::endl;
//// std::cout<<"singleLine.endPointX :"<<singleLine.endPointX<<std::endl;
//// std::cout<<"singleLine.endPointY :"<<singleLine.endPointY<<std::endl;
//// std::cout<<"--------------------------------"<<std::endl;
// }
// }
if(lsd == true)
{
std::cout<<"GetLineDescriptor LSD"<<std::endl;
/* create an LSD detector and store a pointer to it */
cv::Ptr<cv::LineSegmentDetector> ls = cv::createLineSegmentDetector(cv::LSD_REFINE_STD);
/* prepare a vectore to host extracted segments */
std::vector<cv::Vec4i> lines_std;
ScaleLines scaleLs;
/* use detector to extract segments */
ls->detect(original, lines_std);
/* make a copy of input image */
cv::Mat drawnLines(original);
/* draw extracted segments on the copy of image */
ls->drawSegments(drawnLines, lines_std);
/* show extracted segments on image */
cv::imshow("Standard refinement", drawnLines);
/* get a matrix representation of segments */
cv::Mat _lines;
_lines = ((cv::InputArray)lines_std).getMat();
// Draw segments
for(int i = 0; i < _lines.size().width; ++i)
{
/* get current segments and store its extremes */
const cv::Vec4i& v = _lines.at<cv::Vec4i>(i);
cv::Point b(v[0], v[1]);
cv::Point e(v[2], v[3]);
/* create an object to store line information */
OctaveSingleLine osl;
osl.startPointX = b.x;
osl.startPointY = b.y;
osl.endPointX = e.x;
osl.endPointY = e.y;
osl.sPointInOctaveX = b.x;
osl.sPointInOctaveY = b.y;
osl.ePointInOctaveX = e.x;
osl.ePointInOctaveY = e.y;
osl.direction = 0;
osl.salience = 0;
osl.lineLength = 0;
osl.numOfPixels = std::sqrt((b.x-e.x)*(b.x-e.x) + (b.y-e.y)*(b.y-e.y));
osl.octaveCount = 0;
/* store information about line */
LinesVec lineScaleLs;
lineScaleLs.push_back(osl);
scaleLs.push_back(lineScaleLs);
}
keyLines = scaleLs;
}
else
std::cout<<"GetLineDescriptor EDLINE"<<std::endl;
//findNearestParallelLines(keyLines);
ComputeLBD_(keyLines);
return 1;
}
/*Match lines by their descriptors.
*The function will use opencv FlannBasedMatcher to mathc lines. */
int LineDescriptor::MatchLineByDescriptor(ScaleLines &keyLinesLeft, ScaleLines &keyLinesRight,
std::vector<short> &matchLeft, std::vector<short> &matchRight,
int criteria)
{
/* check whether any input is void */
int leftSize = keyLinesLeft.size();
int rightSize = keyLinesRight.size();
if(leftSize<1||rightSize<1){
return -1;
}
/* void vectors */
matchLeft.clear();
matchRight.clear();
int desDim = keyLinesLeft[0][0].descriptor.size();
float *desL, *desR, *desMax, *desOld;
if(criteria==NearestNeighbor)
{
float minDis,dis,temp;
int corresId;
/* loop over left list of LineVecs */
for(int idL=0; idL<leftSize; idL++)
{
/* get size of current left LineVec */
short sameLineSize = keyLinesLeft[idL].size();
/* initialize a distance threshold */
minDis = 100;
/* loop over lines inside current left LineVec */
for(short lineIDInSameLines = 0; lineIDInSameLines<sameLineSize; lineIDInSameLines++)
{
/* get current left line */
desOld = keyLinesLeft[idL][lineIDInSameLines].descriptor.data();
/* loop over right list of LineVecs */
for(int idR=0; idR<rightSize; idR++)
{
/* get size of current right LineVec */
short sameLineSizeR = keyLinesRight[idR].size();
/* loop over lines inside current right LineVec */
for(short lineIDInSameLinesR = 0; lineIDInSameLinesR<sameLineSizeR; lineIDInSameLinesR++)
{
/* store temporay descriptor of left line */
desL = desOld;
/* get descriptor of right line */
desR = keyLinesRight[idR][lineIDInSameLinesR].descriptor.data();
desMax = desR+desDim;
dis = 0;
/* compute euclidean distance between the two descriptors */
while(desR<desMax)
{
temp = *(desL++) - *(desR++);
dis += temp*temp;
}
dis = sqrt(dis);
/* if euclidean distance is below current threshold,
keep track of current match */
if(dis<minDis)
{
minDis = dis;
corresId = idR;
}
}
}//end for(int idR=0; idR<rightSize; idR++)
}//end for(short lineIDInSameLines = 0; lineIDInSameLines<sameLineSize; lineIDInSameLines++)
/* if minimum found distance is below fixed threshold, cnfirm match
by storing corresponding matcjed lines */
if(minDis<LowestThreshold)
{
matchLeft.push_back(idL);
matchRight.push_back(corresId);
}
}// end for(int idL=0; idL<leftSize; idL++)
}
}
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