初衷
用opencv拟合椭圆后，想评估一下拟合的质量，即被拟合点与拟合结果的接近程度。我首先想到的办法是将被拟合点带入椭圆方程
f
(
x
,
y
)
=
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
f(x, y) =Ax^2+Bxy+Cy^2+Dx+Ey+F
f
(
x
,
y
)
=
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
，如果一个点正好在椭圆上，那么
f
(
x
,
y
)
=
0
f(x, y) =0
f
(
x
,
y
)
=
0
，而一个点偏离椭圆越多，则其计算值越大，因此可以用来评估被被拟合点的偏差。 opencv 中 fitEllipse() 函数计算得到的是椭圆的中心坐标、长短轴的长度和长轴与
+
x
+x
+
x
轴的夹角（需要注意的是，拟合结果的 height 才是长轴），需要根据这些信息推导椭圆的标准方程。
推导
中心在原点且长轴在x轴上的椭圆方程为：
x
2
a
2
+
y
2
b
2
=
1
(1)
\dfrac{x^2}{a^2} + \dfrac{y^2}{b^2} = 1 \tag 1
a
2
x
2
+
b
2
y
2
=
1
(
1
)
其中，
a
a
a
和
b
b
b
分别是长半轴和短半轴。若椭圆长轴与正x轴夹角为
θ
\theta
θ
(逆时针)，且中心坐标为
(
x
0
,
y
0
)
(x_0,y_0)
(
x
0
,
y
0
)
。其任意一点平移
(
−
x
0
,
−
y
0
)
(-x_0,-y_0)
(
−
x
0
,
−
y
0
)
再旋转
−
θ
-\theta
−
θ
即满足标准椭圆方程。
首先考虑点
(
x
,
y
)
(x,y)
(
x
,
y
)
绕原点顺时针旋转角度
θ
\theta
θ
到
(
x
′
,
y
′
)
(x’,y’)
(
x
′
,
y
′
)
。
用参数方程表示点
(
x
,
y
)
(x,y)
(
x
,
y
)
和点
(
x
′
,
y
′
)
(x’,y’)
(
x
′
,
y
′
)
：
{
x
=
t
cos
ϕ
y
=
t
sin
ϕ
(2)
\begin{cases} x = t\cos\phi\\ y = t\sin\phi \end{cases} \tag 2
{
x
=
t
cos
ϕ
y
=
t
sin
ϕ
(
2
)
{
x
′
=
t
cos
(
ϕ
−
θ
)
=
t
(
cos
ϕ
cos
θ
+
sin
ϕ
sin
θ
)
y
′
=
t
sin
(
ϕ
−
θ
)
=
t
(
sin
ϕ
cos
θ
−
cos
ϕ
sin
θ
)
(3)
\begin{cases} x’ = t\cos(\phi-\theta)=t(\cos\phi\cos\theta+\sin\phi\sin\theta)\\ y’ = t\sin(\phi-\theta)=t(\sin\phi\cos\theta-\cos\phi\sin\theta) \end{cases} \tag 3
{
x
′
=
t
cos
(
ϕ
−
θ
)
=
t
(
cos
ϕ
cos
θ
+
sin
ϕ
sin
θ
)
y
′
=
t
sin
(
ϕ
−
θ
)
=
t
(
sin
ϕ
cos
θ
−
cos
ϕ
sin
θ
)
(
3
)
因此：
{
x
′
=
x
cos
θ
+
y
sin
θ
y
′
=
−
x
sin
θ
+
y
cos
θ
(4)
\begin{cases} x’ = x\cos\theta+y\sin\theta\\ y’ = -x\sin\theta+y\cos\theta \end{cases} \tag 4
{
x
′
=
x
cos
θ
+
y
sin
θ
y
′
=
−
x
sin
θ
+
y
cos
θ
(
4
)
继而可得一般椭圆方程为：
[
(
x
−
x
0
)
cos
θ
+
(
y
−
y
0
)
sin
θ
]
2
a
2
+
[
(
x
−
x
0
)
sin
θ
−
(
y
−
y
0
)
cos
θ
]
2
b
2
=
1
\dfrac{\left[(x-x_0)\cos\theta+(y-y_0)\sin\theta\right]^2}{a^2}+ \dfrac{\left[(x-x_0)\sin\theta-(y-y_0)\cos\theta\right]^2}{b^2}=1
a
2
[
(
x
−
x
0
)
cos
θ
+
(
y
−
y
0
)
sin
θ
]
2
+
b
2
[
(
x
−
x
0
)
sin
θ
−
(
y
−
y
0
)
cos
θ
]
2
=
1
化简为
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
=
0
Ax^2+Bxy+Cy^2+Dx+Ey+F=0
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
=
0
的形式可得：
A
=
cos
2
θ
a
2
+
sin
2
θ
b
2
B
=
2
sin
θ
cos
θ
(
1
a
2
−
1
b
2
)
C
=
sin
2
θ
a
2
+
cos
2
θ
b
2
D
=
−
2
[
cos
θ
(
x
0
cos
θ
+
y
0
sin
θ
)
a
2
+
sin
θ
(
x
0
sin
θ
−
y
0
cos
θ
)
b
2
]
E
=
−
2
[
sin
θ
(
x
0
cos
θ
+
y
0
sin
θ
)
a
2
−
cos
θ
(
x
0
sin
θ
−
y
0
cos
θ
)
b
2
]
F
=
(
x
0
cos
θ
+
y
0
sin
θ
)
2
a
2
+
(
x
0
sin
θ
−
y
0
cos
θ
)
2
b
2
−
1
(5)
\begin{aligned} A=&\dfrac{\cos^2\theta}{a^2}+\dfrac{\sin^2\theta}{b^2} \\ B=&2\sin\theta\cos\theta\left(\frac{1}{a^2}-\frac{1}{b^2}\right) \\ C=&\dfrac{\sin^2\theta}{a^2}+\dfrac{\cos^2\theta}{b^2} \\ D=&-2\left[\dfrac{\cos\theta(x_0\cos\theta+y_0\sin\theta)}{a^2}+ \dfrac{\sin\theta(x_0\sin\theta-y_0\cos\theta)}{b^2}\right] \\ E=&-2\left[\dfrac{\sin\theta(x_0\cos\theta+y_0\sin\theta)}{a^2}- \dfrac{\cos\theta(x_0\sin\theta-y_0\cos\theta)}{b^2}\right] \\ F=&\dfrac{(x_0\cos\theta+y_0\sin\theta)^2}{a^2}+ \dfrac{(x_0\sin\theta-y_0\cos\theta)^2}{b^2}-1 \end{aligned} \tag 5
A
=
B
=
C
=
D
=
E
=
F
=
a
2
cos
2
θ
+
b
2
sin
2
θ
2
sin
θ
cos
θ
(
a
2
1
−
b
2
1
)
a
2
sin
2
θ
+
b
2
cos
2
θ
−
2
[
a
2
cos
θ
(
x
0
cos
θ
+
y
0
sin
θ
)
+
b
2
sin
θ
(
x
0
sin
θ
−
y
0
cos
θ
)
]
−
2
[
a
2
sin
θ
(
x
0
cos
θ
+
y
0
sin
θ
)
−
b
2
cos
θ
(
x
0
sin
θ
−
y
0
cos
θ
)
]
a
2
(
x
0
cos
θ
+
y
0
sin
θ
)
2
+
b
2
(
x
0
sin
θ
−
y
0
cos
θ
)
2
−
1
(
5
)
代码
在 C++ 中实现上述计算的函数为：
cv::Mat normEllipseParams(cv::RotatedRect box)
{
double params[6];
cv::Mat rst(6, 1, CV_64FC1, params);
double theta = box.angle / 180 * CV_PI;
double st = sin(theta);
double ct = cos(theta);
double a = box.size.width / 2;
double b = box.size.height / 2;
double a2 = a * a;
double b2 = b * b;
double x0 = box.center.x;
double y0 = box.center.y;
double xcys = x0 * ct + y0 * st;
double xsyc = x0 * st - y0 * ct;
params[0] = ct * ct / a2 + st * st / b2;
params[1] = 2 * st * ct * (1 / a2 - 1 / b2);
params[2] = st * st / a2 + ct * ct / b2;
params[3] = -2 * (ct * xcys / a2 + st * xsyc / b2);
params[4] = -2 * (st * xcys / a2 - ct * xsyc / b2);
params[5] = xcys * xcys / a2 + xsyc * xsyc / b2 - 1;
return rst.clone();
}
返回的 Mat 即 A~F 六个参数。
后记
测试后发现一个问题，计算椭圆方程的方法计算比较复杂，而且计算结果受椭圆大小的影响，难以直观地反应点的偏差值。后来我发现其实可以利用“椭圆上的点到其两个焦点的距离之和不变”这一性质来判断一个点偏离椭圆的程度。
-
椭圆的两个焦点在其长轴上，且与中心的距离为
a2
−
b
2
\sqrt{a^2-b^2}
a
2
−
b
2
; -
椭圆上的点到两个焦点的距离之和为
2a
2a
2
a
.
计算椭圆焦点的函数为：
std::array<cv::Point2f, 2> calcFocal(cv::RotatedRect& ellipse)
{
if (ellipse.size.width < ellipse.size.height) {
swap(ellipse.size.width, ellipse.size.height);
ellipse.angle -= 90;
}
float a = ellipse.size.width / 2;
float b = ellipse.size.height / 2;
float c = sqrt(a * a - b * b);
float theta = ellipse.angle / 180 * CV_PI;
cv::Point2f delta(c * cos(theta), c * sin(theta));
return {ellipse.center + delta, ellipse.center - delta};
}
对于任意一个被拟合的点，只要计算其到两焦点的距离之和与
2
a
2a
2
a
的差值，即可知道它与椭圆的像素偏差大小。