Skip to content

O
Opencv

Greenplum / Opencv
11 个月 前同步成功

通知 7
Star 0
Fork 0
O
Opencv

体验新版 GitCode，发现更多精彩内容 >>

提交 d25ee8a2 编写于 作者: V Vadim Pisarevsky
浏览文件
操作

Merge pull request #9761 from Jazmann:ellipseFitAMS&Direct

上级 6c6900a6 800da724
展开全部 隐藏空白更改
内联 并排
Showing with 1596 addition and 37 deletion
doc/opencv.bib
浏览文件 @ d25ee8a2
......@@ -295,6 +295,66 @@
pages = {513--522},
organization = {BMVA Press}
}
@ARTICLE{fitzgibbon1999,
abstract = {This work presents a new efficient method for fitting
ellipses to scattered data. Previous algorithms either
fitted general conics or were computationally expensive. By
minimizing the algebraic distance subject to the constraint
4ac-b<sup>2</sup>=1, the new method incorporates the
ellipticity constraint into the normalization factor. The
proposed method combines several advantages: It is
ellipse-specific, so that even bad data will always return
an ellipse. It can be solved naturally by a generalized
eigensystem. It is extremely robust, efficient, and easy to
implement},
author = {Fitzgibbon, Andrew and Pilu, Maurizio and Fisher, Robert B.},
doi= {10.1109/34.765658},
isbn= {0162-8828},
issn= {01628828},
journal = {IEEE Transactions on Pattern Analysis and Machine
Intelligence},
number = {5},
pages= {476--480},
pmid= {708},
title= {{Direct least square fitting of ellipses}},
volume = {21},
year= {1999}
}
@Article{taubin1991,
abstract = {The author addresses the problem of parametric
representation and estimation of complex planar curves in
2-D surfaces in 3-D, and nonplanar space curves in 3-D.
Curves and surfaces can be defined either parametrically or
implicitly, with the latter representation used here. A
planar curve is the set of zeros of a smooth function of
two variables <e1>x</e1>-<e1>y</e1>, a surface is the set
of zeros of a smooth function of three variables
<e1>x</e1>-<e1>y</e1>-<e1>z</e1>, and a space curve is the
intersection of two surfaces, which are the set of zeros of
two linearly independent smooth functions of three
variables <e1>x</e1>-<e1>y</e1>-<e1>z</e1> For example, the
surface of a complex object in 3-D can be represented as a
subset of a single implicit surface, with similar results
for planar and space curves. It is shown how this unified
representation can be used for object recognition, object
position estimation, and segmentation of objects into
meaningful subobjects, that is, the detection of `interest
regions' that are more complex than high curvature regions
and, hence, more useful as features for object
recognition},
author = {Taubin, Gabriel},
doi= {10.1109/34.103273},
isbn= {0162-8828},
issn= {01628828},
journal = {IEEE Transactions on Pattern Analysis and Machine Intelligence},
number = {11},
pages= {1115--1138},
title= {{Estimation of planar curves, surfaces, and nonplanar
space curves defined by implicit equations with
applications to edge and range image segmentation}},
volume = {13},
year= {1991}
}
@INPROCEEDINGS{G11,
author = {Grundmann, Matthias and Kwatra, Vivek and Essa, Irfan},
title = {Auto-directed video stabilization with robust l1 optimal camera paths},
......
modules/core/include/opencv2/core/utils/logger.hpp
浏览文件 @ d25ee8a2
......@@ -16,6 +16,7 @@
//
//! @{
namespace cv {
namespace utils {
namespace logging {
......@@ -77,7 +78,7 @@ enum LogLevel {
#endif
}} // namespace
}}} // namespace
//! @}
......
modules/imgproc/include/opencv2/imgproc.hpp
浏览文件 @ d25ee8a2
......@@ -4066,6 +4066,88 @@ border of the containing Mat element.
*/
CV_EXPORTS_W RotatedRect fitEllipse( InputArray points );
/** @brief Fits an ellipse around a set of 2D points.
The function calculates the ellipse that fits a set of 2D points.
It returns the rotated rectangle in which the ellipse is inscribed.
The Approximate Mean Square (AMS) proposed by @cite Taubin1991 is used.
For an ellipse, this basis set is \f$ \chi= \left(x^2, x y, y^2, x, y, 1\right) \f$,
which is a set of six free coefficients \f$ A^T=\left\{A_{\text{xx}},A_{\text{xy}},A_{\text{yy}},A_x,A_y,A_0\right\} \f$.
However, to specify an ellipse, all that is needed is five numbers; the major and minor axes lengths \f$ (a,b) \f$,
the position \f$ (x_0,y_0) \f$, and the orientation \f$ \theta \f$. This is because the basis set includes lines,
quadratics, parabolic and hyperbolic functions as well as elliptical functions as possible fits.
If the fit is found to be a parabolic or hyperbolic function then the standard fitEllipse method is used.
The AMS method restricts the fit to parabolic, hyperbolic and elliptical curves
by imposing the condition that \f$ A^T ( D_x^T D_x + D_y^T D_y) A = 1 \f$ where
the matrices \f$ Dx \f$ and \f$ Dy \f$ are the partial derivatives of the design matrix \f$ D \f$ with
respect to x and y. The matrices are formed row by row applying the following to
each of the points in the set:
\f{align*}{
D(i,:)&=\left\{x_i^2, x_i y_i, y_i^2, x_i, y_i, 1\right\} &
D_x(i,:)&=\left\{2 x_i,y_i,0,1,0,0\right\} &
D_y(i,:)&=\left\{0,x_i,2 y_i,0,1,0\right\}
\f}
The AMS method minimizes the cost function
\f{equation*}{
\epsilon ^2=\frac{ A^T D^T D A }{ A^T (D_x^T D_x + D_y^T D_y) A^T }
\f}
The minimum cost is found by solving the generalized eigenvalue problem.
\f{equation*}{
D^T D A = \lambda \left( D_x^T D_x + D_y^T D_y\right) A
\f}
@param points Input 2D point set, stored in std::vector\<\> or Mat
*/
CV_EXPORTS_W RotatedRect fitEllipseAMS( InputArray points );
/** @brief Fits an ellipse around a set of 2D points.
The function calculates the ellipse that fits a set of 2D points.
It returns the rotated rectangle in which the ellipse is inscribed.
The Direct least square (Direct) method by @cite Fitzgibbon1999 is used.
For an ellipse, this basis set is \f$ \chi= \left(x^2, x y, y^2, x, y, 1\right) \f$,
which is a set of six free coefficients \f$ A^T=\left\{A_{\text{xx}},A_{\text{xy}},A_{\text{yy}},A_x,A_y,A_0\right\} \f$.
However, to specify an ellipse, all that is needed is five numbers; the major and minor axes lengths \f$ (a,b) \f$,
the position \f$ (x_0,y_0) \f$, and the orientation \f$ \theta \f$. This is because the basis set includes lines,
quadratics, parabolic and hyperbolic functions as well as elliptical functions as possible fits.
The Direct method confines the fit to ellipses by ensuring that \f$ 4 A_{xx} A_{yy}- A_{xy}^2 > 0 \f$.
The condition imposed is that \f$ 4 A_{xx} A_{yy}- A_{xy}^2=1 \f$ which satisfies the inequality
and as the coefficients can be arbitrarily scaled is not overly restrictive.
\f{equation*}{
\epsilon ^2= A^T D^T D A \quad \text{with} \quad A^T C A =1 \quad \text{and} \quad C=\left(\begin{matrix}
0 & 0 & 2 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 & 0 & 0 \\
2 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0
\end{matrix} \right)
\f}
The minimum cost is found by solving the generalized eigenvalue problem.
\f{equation*}{
D^T D A = \lambda \left( C\right) A
\f}
The system produces only one positive eigenvalue \f$ \lambda\f$ which is chosen as the solution
with its eigenvector \f$\mathbf{u}\f$. These are used to find the coefficients
\f{equation*}{
A = \sqrt{\frac{1}{\mathbf{u}^T C \mathbf{u}}} \mathbf{u}
\f}
The scaling factor guarantees that \f$A^T C A =1\f$.
@param points Input 2D point set, stored in std::vector\<\> or Mat
*/
CV_EXPORTS_W RotatedRect fitEllipseDirect( InputArray points );
/** @brief Fits a line to a 2D or 3D point set.
The function fitLine fits a line to a 2D or 3D point set by minimizing \f$\sum_i \rho(r_i)\f$ where
......
modules/imgproc/src/shapedescr.cpp
浏览文件 @ d25ee8a2
......@@ -39,7 +39,6 @@
//
//M*/
#include "precomp.hpp"
namespace cv
{
......@@ -454,6 +453,329 @@ cv::RotatedRect cv::fitEllipse( InputArray _points )
return box;
}
cv::RotatedRect cv::fitEllipseAMS( InputArray _points )
{
Mat points = _points.getMat();
int i, n = points.checkVector(2);
int depth = points.depth();
CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
RotatedRect box;
if( n < 5 )
CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
Point2f c(0,0);
bool is_float = depth == CV_32F;
const Point* ptsi = points.ptr<Point>();
const Point2f* ptsf = points.ptr<Point2f>();
Mat A( n, 6, CV_64F);
Matx<double, 6, 6> DM;
Matx<double, 5, 5> M;
Matx<double, 5, 1> pVec;
Matx<double, 6, 1> coeffs;
double x0, y0, a, b, theta;
for( i = 0; i < n; i++ )
{
Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
c += p;
}
c.x /= (float)n;
c.y /= (float)n;
for( i = 0; i < n; i++ )
{
Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
p -= c;
A.at<double>(i,0) = (double)(p.x)*(p.x);
A.at<double>(i,1) = (double)(p.x)*(p.y);
A.at<double>(i,2) = (double)(p.y)*(p.y);
A.at<double>(i,3) = (double)p.x;
A.at<double>(i,4) = (double)p.y;
A.at<double>(i,5) = 1.0;
}
cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
DM *= (1.0/n);
double dnm = ( DM(2,5)*(DM(0,5) + DM(2,5)) - (DM(1,5)*DM(1,5)) );
double ddm = (4.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
double ddmm = (2.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
M(0,0)=((-DM(0,0) + DM(0,2) + DM(0,5)*DM(0,5))*(DM(1,5)*DM(1,5)) + (-2*DM(0,1)*DM(1,5) + DM(0,5)*(DM(0,0) \
- (DM(0,5)*DM(0,5)) + (DM(1,5)*DM(1,5))))*DM(2,5) + (DM(0,0) - (DM(0,5)*DM(0,5)))*(DM(2,5)*DM(2,5))) / ddm;
M(0,1)=((DM(1,5)*DM(1,5))*(-DM(0,1) + DM(1,2) + DM(0,5)*DM(1,5)) + (DM(0,1)*DM(0,5) - ((DM(0,5)*DM(0,5)) + 2*DM(1,1))*DM(1,5) + \
(DM(1,5)*DM(1,5)*DM(1,5)))*DM(2,5) + (DM(0,1) - DM(0,5)*DM(1,5))*(DM(2,5)*DM(2,5))) / ddm;
M(0,2)=(-2*DM(1,2)*DM(1,5)*DM(2,5) - DM(0,5)*(DM(2,5)*DM(2,5))*(DM(0,5) + DM(2,5)) + DM(0,2)*dnm + \
(DM(1,5)*DM(1,5))*(DM(2,2) + DM(2,5)*(DM(0,5) + DM(2,5))))/ddm;
M(0,3)=(DM(1,5)*(DM(1,5)*DM(2,3) - 2*DM(1,3)*DM(2,5)) + DM(0,3)*dnm) / ddm;
M(0,4)=(DM(1,5)*(DM(1,5)*DM(2,4) - 2*DM(1,4)*DM(2,5)) + DM(0,4)*dnm) / ddm;
M(1,0)=(-(DM(0,2)*DM(0,5)*DM(1,5)) + (2*DM(0,1)*DM(0,5) - DM(0,0)*DM(1,5))*DM(2,5))/ddmm;
M(1,1)=(-(DM(0,1)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,2)*DM(1,5)) + 2*DM(1,1)*DM(2,5)))/ddmm;
M(1,2)=(-(DM(0,2)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,2)) + 2*DM(1,2)*DM(2,5)))/ddmm;
M(1,3)=(-(DM(0,3)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,3)) + 2*DM(1,3)*DM(2,5)))/ddmm;
M(1,4)=(-(DM(0,4)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,4)) + 2*DM(1,4)*DM(2,5)))/ddmm;
M(2,0)=(-2*DM(0,1)*DM(0,5)*DM(1,5) + (DM(0,0) + (DM(0,5)*DM(0,5)))*(DM(1,5)*DM(1,5)) + DM(0,5)*(-(DM(0,5)*DM(0,5)) \
+ (DM(1,5)*DM(1,5)))*DM(2,5) - (DM(0,5)*DM(0,5))*(DM(2,5)*DM(2,5)) + DM(0,2)*(-(DM(1,5)*DM(1,5)) + DM(0,5)*(DM(0,5) + DM(2,5)))) / ddm;
M(2,1)=((DM(0,5)*DM(0,5))*(DM(1,2) - DM(1,5)*DM(2,5)) + (DM(1,5)*DM(1,5))*(DM(0,1) - DM(1,2) + DM(1,5)*DM(2,5)) \
+ DM(0,5)*(DM(1,2)*DM(2,5) + DM(1,5)*(-2*DM(1,1) + (DM(1,5)*DM(1,5)) - (DM(2,5)*DM(2,5))))) / ddm;
M(2,2)=((DM(0,5)*DM(0,5))*(DM(2,2) - (DM(2,5)*DM(2,5))) + (DM(1,5)*DM(1,5))*(DM(0,2) - DM(2,2) + (DM(2,5)*DM(2,5))) + \
DM(0,5)*(-2*DM(1,2)*DM(1,5) + DM(2,5)*((DM(1,5)*DM(1,5)) + DM(2,2) - (DM(2,5)*DM(2,5))))) / ddm;
M(2,3)=((DM(1,5)*DM(1,5))*(DM(0,3) - DM(2,3)) + (DM(0,5)*DM(0,5))*DM(2,3) + DM(0,5)*(-2*DM(1,3)*DM(1,5) + DM(2,3)*DM(2,5))) / ddm;
M(2,4)=((DM(1,5)*DM(1,5))*(DM(0,4) - DM(2,4)) + (DM(0,5)*DM(0,5))*DM(2,4) + DM(0,5)*(-2*DM(1,4)*DM(1,5) + DM(2,4)*DM(2,5))) / ddm;
M(3,0)=DM(0,3);
M(3,1)=DM(1,3);
M(3,2)=DM(2,3);
M(3,3)=DM(3,3);
M(3,4)=DM(3,4);
M(4,0)=DM(0,4);
M(4,1)=DM(1,4);
M(4,2)=DM(2,4);
M(4,3)=DM(3,4);
M(4,4)=DM(4,4);
if (fabs(cv::determinant(M)) > 1.0e-10) {
Mat eVal, eVec;
eigenNonSymmetric(M, eVal, eVec);
// Select the eigen vector {a,b,c,d,e} which has the lowest eigenvalue
int minpos = 0;
double normi, normEVali, normMinpos, normEValMinpos;
normMinpos = sqrt(eVec.at<double>(minpos,0)*eVec.at<double>(minpos,0) + eVec.at<double>(minpos,1)*eVec.at<double>(minpos,1) + \
eVec.at<double>(minpos,2)*eVec.at<double>(minpos,2) + eVec.at<double>(minpos,3)*eVec.at<double>(minpos,3) + \
eVec.at<double>(minpos,4)*eVec.at<double>(minpos,4) );
normEValMinpos = eVal.at<double>(minpos,0) * normMinpos;
for (i=1; i<5; i++) {
normi = sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + \
eVec.at<double>(i,2)*eVec.at<double>(i,2) + eVec.at<double>(i,3)*eVec.at<double>(i,3) + \
eVec.at<double>(i,4)*eVec.at<double>(i,4) );
normEVali = eVal.at<double>(i,0) * normi;
if (normEVali < normEValMinpos) {
minpos = i;
normMinpos=normi;
normEValMinpos=normEVali;
}
};
pVec(0) =eVec.at<double>(minpos,0) / normMinpos;
pVec(1) =eVec.at<double>(minpos,1) / normMinpos;
pVec(2) =eVec.at<double>(minpos,2) / normMinpos;
pVec(3) =eVec.at<double>(minpos,3) / normMinpos;
pVec(4) =eVec.at<double>(minpos,4) / normMinpos;
coeffs(0) =pVec(0) ;
coeffs(1) =pVec(1) ;
coeffs(2) =pVec(2) ;
coeffs(3) =pVec(3) ;
coeffs(4) =pVec(4) ;
coeffs(5) =-pVec(0) *DM(0,5)-pVec(1) *DM(1,5)-coeffs(2) *DM(2,5);
// Check that an elliptical solution has been found. AMS sometimes produces Parabolic solutions.
bool is_ellipse = (coeffs(0) < 0 && \
coeffs(2) < (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
coeffs(5) > (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4) - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
((coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) )) || \
(coeffs(0) > 0 && \
coeffs(2) > (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
coeffs(5) < (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4) - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
( (coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) ));
if (is_ellipse) {
double u1 = pVec(2) *pVec(3) *pVec(3) - pVec(1) *pVec(3) *pVec(4) + pVec(0) *pVec(4) *pVec(4) + pVec(1) *pVec(1) *coeffs(5) ;
double u2 = pVec(0) *pVec(2) *coeffs(5) ;
double l1 = sqrt(pVec(1) *pVec(1) + (pVec(0) - pVec(2) )*(pVec(0) - pVec(2) ));
double l2 = pVec(0) + pVec(2) ;
double l3 = pVec(1) *pVec(1) - 4.0*pVec(0) *pVec(2) ;
double p1 = 2.0*pVec(2) *pVec(3) - pVec(1) *pVec(4) ;
double p2 = 2.0*pVec(0) *pVec(4) -(pVec(1) *pVec(3) );
x0 = p1/l3 + c.x;
y0 = p2/l3 + c.y;
a = sqrt(2)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3));
b = sqrt(2)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)));
if (pVec(1) == 0) {
if (pVec(0) < pVec(2) ) {
theta = 0;
} else {
theta = CV_PI/2.;
}
} else {
theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0) - pVec(2) ));
}
box.center.x = (float)x0; // +c.x;
box.center.y = (float)y0; // +c.y;
box.size.width = (float)(2.0*a);
box.size.height = (float)(2.0*b);
if( box.size.width > box.size.height )
{
float tmp;
CV_SWAP( box.size.width, box.size.height, tmp );
box.angle = (float)(90 + theta*180/CV_PI);
} else {
box.angle = (float)(fmod(theta*180/CV_PI,180.0));
};
} else {
box = cv::fitEllipseDirect( points );
}
} else {
box = cv::fitEllipse( points );
}
return box;
}
cv::RotatedRect cv::fitEllipseDirect( InputArray _points )
{
Mat points = _points.getMat();
int i, n = points.checkVector(2);
int depth = points.depth();
CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
RotatedRect box;
if( n < 5 )
CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
Point2f c(0,0);
bool is_float = (depth == CV_32F);
const Point* ptsi = points.ptr<Point>();
const Point2f* ptsf = points.ptr<Point2f>();
Mat A( n, 6, CV_64F);
Matx<double, 6, 6> DM;
Matx33d M, TM, Q;
Matx<double, 3, 1> pVec;
double x0, y0, a, b, theta, Ts;
for( i = 0; i < n; i++ )
{
Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
c += p;
}
c.x /= (float)n;
c.y /= (float)n;
for( i = 0; i < n; i++ )
{
Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
p -= c;
A.at<double>(i,0) = (double)(p.x)*(p.x);
A.at<double>(i,1) = (double)(p.x)*(p.y);
A.at<double>(i,2) = (double)(p.y)*(p.y);
A.at<double>(i,3) = (double)p.x;
A.at<double>(i,4) = (double)p.y;
A.at<double>(i,5) = 1.0;
}
cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
DM *= (1.0/n);
TM(0,0) = DM(0,5)*DM(3,5)*DM(4,4) - DM(0,5)*DM(3,4)*DM(4,5) - DM(0,4)*DM(3,5)*DM(5,4) + \
DM(0,3)*DM(4,5)*DM(5,4) + DM(0,4)*DM(3,4)*DM(5,5) - DM(0,3)*DM(4,4)*DM(5,5);
TM(0,1) = DM(1,5)*DM(3,5)*DM(4,4) - DM(1,5)*DM(3,4)*DM(4,5) - DM(1,4)*DM(3,5)*DM(5,4) + \
DM(1,3)*DM(4,5)*DM(5,4) + DM(1,4)*DM(3,4)*DM(5,5) - DM(1,3)*DM(4,4)*DM(5,5);
TM(0,2) = DM(2,5)*DM(3,5)*DM(4,4) - DM(2,5)*DM(3,4)*DM(4,5) - DM(2,4)*DM(3,5)*DM(5,4) + \
DM(2,3)*DM(4,5)*DM(5,4) + DM(2,4)*DM(3,4)*DM(5,5) - DM(2,3)*DM(4,4)*DM(5,5);
TM(1,0) = DM(0,5)*DM(3,3)*DM(4,5) - DM(0,5)*DM(3,5)*DM(4,3) + DM(0,4)*DM(3,5)*DM(5,3) - \
DM(0,3)*DM(4,5)*DM(5,3) - DM(0,4)*DM(3,3)*DM(5,5) + DM(0,3)*DM(4,3)*DM(5,5);
TM(1,1) = DM(1,5)*DM(3,3)*DM(4,5) - DM(1,5)*DM(3,5)*DM(4,3) + DM(1,4)*DM(3,5)*DM(5,3) - \
DM(1,3)*DM(4,5)*DM(5,3) - DM(1,4)*DM(3,3)*DM(5,5) + DM(1,3)*DM(4,3)*DM(5,5);
TM(1,2) = DM(2,5)*DM(3,3)*DM(4,5) - DM(2,5)*DM(3,5)*DM(4,3) + DM(2,4)*DM(3,5)*DM(5,3) - \
DM(2,3)*DM(4,5)*DM(5,3) - DM(2,4)*DM(3,3)*DM(5,5) + DM(2,3)*DM(4,3)*DM(5,5);
TM(2,0) = DM(0,5)*DM(3,4)*DM(4,3) - DM(0,5)*DM(3,3)*DM(4,4) - DM(0,4)*DM(3,4)*DM(5,3) + \
DM(0,3)*DM(4,4)*DM(5,3) + DM(0,4)*DM(3,3)*DM(5,4) - DM(0,3)*DM(4,3)*DM(5,4);
TM(2,1) = DM(1,5)*DM(3,4)*DM(4,3) - DM(1,5)*DM(3,3)*DM(4,4) - DM(1,4)*DM(3,4)*DM(5,3) + \
DM(1,3)*DM(4,4)*DM(5,3) + DM(1,4)*DM(3,3)*DM(5,4) - DM(1,3)*DM(4,3)*DM(5,4);
TM(2,2) = DM(2,5)*DM(3,4)*DM(4,3) - DM(2,5)*DM(3,3)*DM(4,4) - DM(2,4)*DM(3,4)*DM(5,3) + \
DM(2,3)*DM(4,4)*DM(5,3) + DM(2,4)*DM(3,3)*DM(5,4) - DM(2,3)*DM(4,3)*DM(5,4);
Ts=(-(DM(3,5)*DM(4,4)*DM(5,3)) + DM(3,4)*DM(4,5)*DM(5,3) + DM(3,5)*DM(4,3)*DM(5,4) - \
DM(3,3)*DM(4,5)*DM(5,4) - DM(3,4)*DM(4,3)*DM(5,5) + DM(3,3)*DM(4,4)*DM(5,5));
M(0,0) = (DM(2,0) + (DM(2,3)*TM(0,0) + DM(2,4)*TM(1,0) + DM(2,5)*TM(2,0))/Ts)/2.;
M(0,1) = (DM(2,1) + (DM(2,3)*TM(0,1) + DM(2,4)*TM(1,1) + DM(2,5)*TM(2,1))/Ts)/2.;
M(0,2) = (DM(2,2) + (DM(2,3)*TM(0,2) + DM(2,4)*TM(1,2) + DM(2,5)*TM(2,2))/Ts)/2.;
M(1,0) = -DM(1,0) - (DM(1,3)*TM(0,0) + DM(1,4)*TM(1,0) + DM(1,5)*TM(2,0))/Ts;
M(1,1) = -DM(1,1) - (DM(1,3)*TM(0,1) + DM(1,4)*TM(1,1) + DM(1,5)*TM(2,1))/Ts;
M(1,2) = -DM(1,2) - (DM(1,3)*TM(0,2) + DM(1,4)*TM(1,2) + DM(1,5)*TM(2,2))/Ts;
M(2,0) = (DM(0,0) + (DM(0,3)*TM(0,0) + DM(0,4)*TM(1,0) + DM(0,5)*TM(2,0))/Ts)/2.;
M(2,1) = (DM(0,1) + (DM(0,3)*TM(0,1) + DM(0,4)*TM(1,1) + DM(0,5)*TM(2,1))/Ts)/2.;
M(2,2) = (DM(0,2) + (DM(0,3)*TM(0,2) + DM(0,4)*TM(1,2) + DM(0,5)*TM(2,2))/Ts)/2.;
if (fabs(cv::determinant(M)) > 1.0e-10) {
Mat eVal, eVec;
eigenNonSymmetric(M, eVal, eVec);
// Select the eigen vector {a,b,c} which satisfies 4ac-b^2 > 0
double cond[3];
cond[0]=(4.0 * eVec.at<double>(0,0) * eVec.at<double>(0,2) - eVec.at<double>(0,1) * eVec.at<double>(0,1));
cond[1]=(4.0 * eVec.at<double>(1,0) * eVec.at<double>(1,2) - eVec.at<double>(1,1) * eVec.at<double>(1,1));
cond[2]=(4.0 * eVec.at<double>(2,0) * eVec.at<double>(2,2) - eVec.at<double>(2,1) * eVec.at<double>(2,1));
if (cond[0]<cond[1]) {
i = (cond[1]<cond[2]) ? 2 : 1;
} else {
i = (cond[0]<cond[2]) ? 2 : 0;
}
double norm = std::sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + eVec.at<double>(i,2)*eVec.at<double>(i,2));
if (((eVec.at<double>(i,0)<0.0 ? -1 : 1) * (eVec.at<double>(i,1)<0.0 ? -1 : 1) * (eVec.at<double>(i,2)<0.0 ? -1 : 1)) <= 0.0) {
norm=-1.0*norm;
}
pVec(0) =eVec.at<double>(i,0)/norm; pVec(1) =eVec.at<double>(i,1)/norm;pVec(2) =eVec.at<double>(i,2)/norm;
// Q = (TM . pVec)/Ts;
Q(0,0) = (TM(0,0)*pVec(0) +TM(0,1)*pVec(1) +TM(0,2)*pVec(2) )/Ts;
Q(0,1) = (TM(1,0)*pVec(0) +TM(1,1)*pVec(1) +TM(1,2)*pVec(2) )/Ts;
Q(0,2) = (TM(2,0)*pVec(0) +TM(2,1)*pVec(1) +TM(2,2)*pVec(2) )/Ts;
// We compute the ellipse properties in the shifted coordinates as doing so improves the numerical accuracy.
double u1 = pVec(2)*Q(0,0)*Q(0,0) - pVec(1)*Q(0,0)*Q(0,1) + pVec(0)*Q(0,1)*Q(0,1) + pVec(1)*pVec(1)*Q(0,2);
double u2 = pVec(0)*pVec(2)*Q(0,2);
double l1 = sqrt(pVec(1)*pVec(1) + (pVec(0) - pVec(2))*(pVec(0) - pVec(2)));
double l2 = pVec(0) + pVec(2) ;
double l3 = pVec(1)*pVec(1) - 4*pVec(0)*pVec(2) ;
double p1 = 2*pVec(2)*Q(0,0) - pVec(1)*Q(0,1);
double p2 = 2*pVec(0)*Q(0,1) - pVec(1)*Q(0,0);
x0 = p1/l3 + c.x;
y0 = p2/l3 + c.y;
a = sqrt(2)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3));
b = sqrt(2)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)));
if (pVec(1) == 0) {
if (pVec(0) < pVec(2) ) {
theta = 0;
} else {
theta = CV_PI/2.;
}
} else {
theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0) - pVec(2) ));
}
box.center.x = (float)x0;
box.center.y = (float)y0;
box.size.width = (float)(2.0*a);
box.size.height = (float)(2.0*b);
if( box.size.width > box.size.height )
{
float tmp;
CV_SWAP( box.size.width, box.size.height, tmp );
box.angle = (float)(fmod((90 + theta*180/CV_PI),180.0)) ;
} else {
box.angle = (float)(fmod(theta*180/CV_PI,180.0));
};
} else {
box = cv::fitEllipse( points );
}
return box;
}
namespace cv
{
......@@ -1080,5 +1402,4 @@ cvBoundingRect( CvArr* array, int update )
return rect;
}
/* End of file. */
modules/imgproc/test/test_fitellipseAMS.cpp 0 → 100644
浏览文件 @ d25ee8a2
此差异已折叠。
modules/imgproc/test/test_fitellipseDirect.cpp 0 → 100644
浏览文件 @ d25ee8a2
此差异已折叠。
samples/cpp/fitellipse.cpp
浏览文件 @ d25ee8a2
/********************************************************************************
*
*
* This program is demonstration for ellipse fitting. Program finds
* contours and approximate it by ellipses.
*
* Trackbar specify threshold parametr.
*
* White lines is contours. Red lines is fitting ellipses.
*
*
* Autor: Denis Burenkov.
*
*
*
********************************************************************************/
*
*
* This program is demonstration for ellipse fitting. Program finds
* contours and approximate it by ellipses using three methods.
* 1: OpenCV's original method fitEllipse which implements Fitzgibbon 1995 method.
* 2: The Approximate Mean Square (AMS) method fitEllipseAMS proposed by Taubin 1991
* 3: The Direct least square (Direct) method fitEllipseDirect proposed by Fitzgibbon1999.
*
* Trackbar specify threshold parameter.
*
* White lines is contours/input points and the true ellipse used to generate the data.
* 1: Blue lines is fitting ellipses using openCV's original method.
* 2: Green lines is fitting ellipses using the AMS method.
* 3: Red lines is fitting ellipses using the Direct method.
*
*
* Original Author: Denis Burenkov
* AMS and Direct Methods Autor: Jasper Shemilt
*
*
********************************************************************************/
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
......@@ -22,26 +28,171 @@
using namespace cv;
using namespace std;
class canvas{
public:
bool setupQ;
cv::Point origin;
cv::Point corner;
int minDims,maxDims;
double scale;
int rows, cols;
cv::Mat img;
void init(int minD, int maxD){
// Initialise the canvas with minimum and maximum rows and column sizes.
minDims = minD; maxDims = maxD;
origin = cv::Point(0,0);
corner = cv::Point(0,0);
scale = 1.0;
rows = 0;
cols = 0;
setupQ = false;
}
void stretch(cv::Point2f min, cv::Point2f max){
// Stretch the canvas to include the points min and max.
if(setupQ){
if(corner.x < max.x){corner.x = (int)(max.x + 1.0);};
if(corner.y < max.y){corner.y = (int)(max.y + 1.0);};
if(origin.x > min.x){origin.x = (int) min.x;};
if(origin.y > min.y){origin.y = (int) min.y;};
} else {
origin = cv::Point((int)min.x, (int)min.y);
corner = cv::Point((int)(max.x + 1.0), (int)(max.y + 1.0));
}
int c = (int)(scale*((corner.x + 1.0) - origin.x));
if(c<minDims){
scale = scale * (double)minDims/(double)c;
} else {
if(c>maxDims){
scale = scale * (double)maxDims/(double)c;
}
}
int r = (int)(scale*((corner.y + 1.0) - origin.y));
if(r<minDims){
scale = scale * (double)minDims/(double)r;
} else {
if(r>maxDims){
scale = scale * (double)maxDims/(double)r;
}
}
cols = (int)(scale*((corner.x + 1.0) - origin.x));
rows = (int)(scale*((corner.y + 1.0) - origin.y));
setupQ = true;
}
void stretch(vector<Point2f> pts)
{ // Stretch the canvas so all the points pts are on the canvas.
cv::Point2f min = pts[0];
cv::Point2f max = pts[0];
for(size_t i=1; i < pts.size(); i++){
Point2f pnt = pts[i];
if(max.x < pnt.x){max.x = pnt.x;};
if(max.y < pnt.y){max.y = pnt.y;};
if(min.x > pnt.x){min.x = pnt.x;};
if(min.y > pnt.y){min.y = pnt.y;};
};
stretch(min, max);
}
void stretch(cv::RotatedRect box)
{ // Stretch the canvas so that the rectangle box is on the canvas.
cv::Point2f min = box.center;
cv::Point2f max = box.center;
cv::Point2f vtx[4];
box.points(vtx);
for( int i = 0; i < 4; i++ ){
cv::Point2f pnt = vtx[i];
if(max.x < pnt.x){max.x = pnt.x;};
if(max.y < pnt.y){max.y = pnt.y;};
if(min.x > pnt.x){min.x = pnt.x;};
if(min.y > pnt.y){min.y = pnt.y;};
}
stretch(min, max);
}
void drawEllipseWithBox(cv::RotatedRect box, cv::Scalar color, int lineThickness)
{
if(img.empty()){
stretch(box);
img = cv::Mat::zeros(rows,cols,CV_8UC3);
}
box.center = scale * cv::Point2f(box.center.x - origin.x, box.center.y - origin.y);
box.size.width = (float)(scale * box.size.width);
box.size.height = (float)(scale * box.size.height);
ellipse(img, box, color, lineThickness, LINE_AA);
Point2f vtx[4];
box.points(vtx);
for( int j = 0; j < 4; j++ ){
line(img, vtx[j], vtx[(j+1)%4], color, lineThickness, LINE_AA);
}
}
void drawPoints(vector<Point2f> pts, cv::Scalar color)
{
if(img.empty()){
stretch(pts);
img = cv::Mat::zeros(rows,cols,CV_8UC3);
}
for(size_t i=0; i < pts.size(); i++){
Point2f pnt = scale * cv::Point2f(pts[i].x - origin.x, pts[i].y - origin.y);
img.at<cv::Vec3b>(int(pnt.y), int(pnt.x))[0] = (uchar)color[0];
img.at<cv::Vec3b>(int(pnt.y), int(pnt.x))[1] = (uchar)color[1];
img.at<cv::Vec3b>(int(pnt.y), int(pnt.x))[2] = (uchar)color[2];
};
}
void drawLabels( std::vector<std::string> text, std::vector<cv::Scalar> colors)
{
if(img.empty()){
img = cv::Mat::zeros(rows,cols,CV_8UC3);
}
int vPos = 0;
for (size_t i=0; i < text.size(); i++) {
cv::Scalar color = colors[i];
std::string txt = text[i];
Size textsize = getTextSize(txt, FONT_HERSHEY_COMPLEX, 1, 1, 0);
vPos += (int)(1.3 * textsize.height);
Point org((img.cols - textsize.width), vPos);
cv::putText(img, txt, org, FONT_HERSHEY_COMPLEX, 1, color, 1, LINE_8);
}
}
};
static void help()
{
cout <<
"\nThis program is demonstration for ellipse fitting. The program finds\n"
"contours and approximate it by ellipses.\n"
"Call:\n"
"./fitellipse [image_name -- Default ../data/stuff.jpg]\n" << endl;
"\nThis program is demonstration for ellipse fitting. The program finds\n"
"contours and approximate it by ellipses. Three methods are used to find the \n"
"elliptical fits: fitEllipse, fitEllipseAMS and fitEllipseDirect.\n"
"Call:\n"
"./fitellipse [image_name -- Default ../data/stuff.jpg]\n" << endl;
}
int sliderPos = 70;
Mat image;
bool fitEllipseQ, fitEllipseAMSQ, fitEllipseDirectQ;
cv::Scalar fitEllipseColor = Scalar(255, 0, 0);
cv::Scalar fitEllipseAMSColor = Scalar( 0,255, 0);
cv::Scalar fitEllipseDirectColor = Scalar( 0, 0,255);
cv::Scalar fitEllipseTrueColor = Scalar(255,255,255);
void processImage(int, void*);
int main( int argc, char** argv )
{
cv::CommandLineParser parser(argc, argv,
"{help h||}{@image|../data/stuff.jpg|}"
);
fitEllipseQ = true;
fitEllipseAMSQ = true;
fitEllipseDirectQ = true;
cv::CommandLineParser parser(argc, argv,"{help h||}{@image|../data/ellipses.jpg|}");
if (parser.has("help"))
{
help();
......@@ -56,10 +207,11 @@ int main( int argc, char** argv )
}
imshow("source", image);
namedWindow("result", 1);
namedWindow("result", CV_WINDOW_NORMAL );
// Create toolbars. HighGUI use.
createTrackbar( "threshold", "result", &sliderPos, 255, processImage );
processImage(0, 0);
// Wait for a key stroke; the same function arranges events processing
......@@ -71,13 +223,35 @@ int main( int argc, char** argv )
// draw it and approximate it by ellipses.
void processImage(int /*h*/, void*)
{
RotatedRect box, boxAMS, boxDirect;
vector<vector<Point> > contours;
Mat bimage = image >= sliderPos;
findContours(bimage, contours, RETR_LIST, CHAIN_APPROX_NONE);
Mat cimage = Mat::zeros(bimage.size(), CV_8UC3);
canvas paper;
paper.init(int(0.8*MIN(bimage.rows, bimage.cols)), int(1.2*MAX(bimage.rows, bimage.cols)));
paper.stretch(cv::Point2f(0.0f, 0.0f), cv::Point2f((float)(bimage.cols+2.0), (float)(bimage.rows+2.0)));
std::vector<std::string> text;
std::vector<cv::Scalar> color;
if (fitEllipseQ) {
text.push_back("OpenCV");
color.push_back(fitEllipseColor);
}
if (fitEllipseAMSQ) {
text.push_back("AMS");
color.push_back(fitEllipseAMSColor);
}
if (fitEllipseDirectQ) {
text.push_back("Direct");
color.push_back(fitEllipseDirectColor);
}
paper.drawLabels(text, color);
int margin = 2;
vector< vector<Point2f> > points;
for(size_t i = 0; i < contours.size(); i++)
{
size_t count = contours[i].size();
......@@ -86,19 +260,57 @@ void processImage(int /*h*/, void*)
Mat pointsf;
Mat(contours[i]).convertTo(pointsf, CV_32F);
RotatedRect box = fitEllipse(pointsf);
if( MAX(box.size.width, box.size.height) > MIN(box.size.width, box.size.height)*30 )
vector<Point2f>pts;
for (int j = 0; j < pointsf.rows; j++) {
Point2f pnt = Point2f(pointsf.at<float>(j,0), pointsf.at<float>(j,1));
if ((pnt.x > margin && pnt.y > margin && pnt.x < bimage.cols-margin && pnt.y < bimage.rows-margin)) {
if(j%20==0){
pts.push_back(pnt);
}
}
}
points.push_back(pts);
}
for(size_t i = 0; i < points.size(); i++)
{
vector<Point2f> pts = points[i];
if (pts.size()<=5) {
continue;
drawContours(cimage, contours, (int)i, Scalar::all(255), 1, 8);
}
if (fitEllipseQ) {
box = fitEllipse(pts);
if( MAX(box.size.width, box.size.height) > MIN(box.size.width, box.size.height)*30 ||
MAX(box.size.width, box.size.height) <= 0 ||
MIN(box.size.width, box.size.height) <= 0){continue;};
}
if (fitEllipseAMSQ) {
boxAMS = fitEllipseAMS(pts);
if( MAX(boxAMS.size.width, boxAMS.size.height) > MIN(boxAMS.size.width, boxAMS.size.height)*30 ||
MAX(box.size.width, box.size.height) <= 0 ||
MIN(box.size.width, box.size.height) <= 0){continue;};
}
if (fitEllipseDirectQ) {
boxDirect = fitEllipseDirect(pts);
if( MAX(boxDirect.size.width, boxDirect.size.height) > MIN(boxDirect.size.width, boxDirect.size.height)*30 ||
MAX(box.size.width, box.size.height) <= 0 ||
MIN(box.size.width, box.size.height) <= 0 ){continue;};
}
ellipse(cimage, box, Scalar(0,0,255), 1, LINE_AA);
ellipse(cimage, box.center, box.size*0.5f, box.angle, 0, 360, Scalar(0,255,255), 1, LINE_AA);
Point2f vtx[4];
box.points(vtx);
for( int j = 0; j < 4; j++ )
line(cimage, vtx[j], vtx[(j+1)%4], Scalar(0,255,0), 1, LINE_AA);
if (fitEllipseQ) {
paper.drawEllipseWithBox(box, fitEllipseColor, 3);
}
if (fitEllipseAMSQ) {
paper.drawEllipseWithBox(boxAMS, fitEllipseAMSColor, 2);
}
if (fitEllipseDirectQ) {
paper.drawEllipseWithBox(boxDirect, fitEllipseDirectColor, 1);
}
paper.drawPoints(pts, cv::Scalar(255,255,255));
}
imshow("result", cimage);
imshow("result", paper.img);
}
samples/data/ellipses.jpg 0 → 100644
浏览文件 @ d25ee8a2
此差异由.gitattributes 抑制。
Markdown is supported
0% 或 .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑！
想要评论请 注册 或