标准椭圆方程推导

　

椭圆方程

初衷

用opencv拟合椭圆后，想评估一下拟合的质量，即被拟合点与拟合结果的接近程度。我首先想到的办法是将被拟合点带入椭圆方程 f(x,y)=Ax2+Bxy+Cy2+Dx+Ey+F" role="presentation">f(x,y)=Ax2+Bxy+Cy2+Dx+Ey+F$f(x,y)=A{x}^2+Bxy+C{y}^2+Dx+Ey+F$，如果一个点正好在椭圆上，那么 f(x,y)=0" role="presentation">f(x,y)=0$f(x,y)=0$，而一个点偏离椭圆越多，则其计算值越大，因此可以用来评估被被拟合点的偏差。 opencv 中 fitEllipse() 函数计算得到的是椭圆的中心坐标、长短轴的长度和长轴与 +x" role="presentation">+x$+x$ 轴的夹角（需要注意的是，拟合结果的 height 才是长轴），需要根据这些信息推导椭圆的标准方程。

---

推导

中心在原点且长轴在x轴上的椭圆方程为：

(1)x2a2+y2b2=1" role="presentation">x2a2+y2b2=1(1)$$\begin{array}{}\text{(1)}& {\displaystyle \frac{x^2}{a^2}}+{\displaystyle \frac{y^2}{b^2}}=1\end{array}$$

其中，a" role="presentation">a$a$和b" role="presentation">b$b$分别是长半轴和短半轴。若椭圆长轴与正x轴夹角为θ" role="presentation">θ$\theta$(逆时针)，且中心坐标为(x0,y0)" role="presentation">(x0,y0)$(x_0,y_0)$。其任意一点平移(−x0,−y0)" role="presentation">(−x0,−y0)$(-x_0,-y_0)$再旋转−θ" role="presentation">−θ$-\theta$即满足标准椭圆方程。

首先考虑点(x,y)" role="presentation">(x,y)$(x,y)$绕原点顺时针旋转角度θ" role="presentation">θ$\theta$到(x′,y′)" role="presentation">(x′,y′)$(x^{\prime },y^{\prime })$。

用参数方程表示点(x,y)" role="presentation">(x,y)$(x,y)$和点(x′,y′)" role="presentation">(x′,y′)$(x^{\prime },y^{\prime })$：

(2){x=tcos⁡ϕy=tsin⁡ϕ" role="presentation">{x=tcosϕy=tsinϕ(2)$$\begin{array}{}\text{(2)}& \{\begin{array}{l}x=t\cos \varphi \\ y=t\sin \varphi \end{array}\end{array}$$

(3){x′=tcos⁡(ϕ−θ)=t(cos⁡ϕcos⁡θ+sin⁡ϕsin⁡θ)y′=tsin⁡(ϕ−θ)=t(sin⁡ϕcos⁡θ−cos⁡ϕsin⁡θ)" role="presentation">{x′=tcos(ϕ−θ)=t(cosϕcosθ+sinϕsinθ)y′=tsin(ϕ−θ)=t(sinϕcosθ−cosϕsinθ)(3)$$\begin{array}{}\text{(3)}& \{\begin{array}{l}{x}^{\prime }=t\cos (\varphi -\theta )=t(\cos \varphi \cos \theta +\sin \varphi \sin \theta )\\ y^{\prime }=t\sin (\varphi -\theta )=t(\sin \varphi \cos \theta -\cos \varphi \sin \theta )\end{array}\end{array}$$

因此：

(4){x′=xcos⁡θ+ysin⁡θy′=−xsin⁡θ+ycos⁡θ" role="presentation">{x′=xcosθ+ysinθy′=−xsinθ+ycosθ(4)$$\begin{array}{}\text{(4)}& \{\begin{array}{l}{x}^{\prime }=x\cos \theta +y\sin \theta \\ y^{\prime }=-x\sin \theta +y\cos \theta \end{array}\end{array}$$

继而可得一般椭圆方程为：

[(x−x0)cos⁡θ+(y−y0)sin⁡θ]2a2+[(x−x0)sin⁡θ−(y−y0)cos⁡θ]2b2=1" role="presentation">[(x−x0)cosθ+(y−y0)sinθ]2a2+[(x−x0)sinθ−(y−y0)cosθ]2b2=1$${\displaystyle \frac{{\left[(x-x_0)\cos \theta +(y-y_0)\sin \theta \right]}^2}{a^2}}+{\displaystyle \frac{{\left[(x-x_0)\sin \theta -(y-y_0)\cos \theta \right]}^2}{b^2}}=1$$

化简为 Ax2+Bxy+Cy2+Dx+Ey+F=0" role="presentation">Ax2+Bxy+Cy2+Dx+Ey+F=0$A{x}^2+Bxy+C{y}^2+Dx+Ey+F=0$ 的形式可得：

(5)A=cos2⁡θa2+sin2⁡θb2B=2sin⁡θcos⁡θ(1a2−1b2)C=sin2⁡θa2+cos2⁡θb2D=−2[cos⁡θ(x0cos⁡θ+y0sin⁡θ)a2+sin⁡θ(x0sin⁡θ−y0cos⁡θ)b2]E=−2[sin⁡θ(x0cos⁡θ+y0sin⁡θ)a2−cos⁡θ(x0sin⁡θ−y0cos⁡θ)b2]F=(x0cos⁡θ+y0sin⁡θ)2a2+(x0sin⁡θ−y0cos⁡θ)2b2−1" role="presentation">A=B=C=D=E=F=cos2θa2+sin2θb22sinθcosθ(1a2−1b2)sin2θa2+cos2θb2−2[cosθ(x0cosθ+y0sinθ)a2+sinθ(x0sinθ−y0cosθ)b2]−2[sinθ(x0cosθ+y0sinθ)a2−cosθ(x0sinθ−y0cosθ)b2](x0cosθ+y0sinθ)2a2+(x0sinθ−y0cosθ)2b2−1(5)$$\begin{array}{}\text{(5)}& \begin{array}{rl}A=& {\displaystyle \frac{{\cos }^2\theta }{a^2}}+{\displaystyle \frac{{\sin }^2\theta }{b^2}}\\ B=& 2\sin \theta \cos \theta \left(\frac{1}{a^2}-\frac{1}{b^2}\right)\\ C=& {\displaystyle \frac{{\sin }^2\theta }{a^2}}+{\displaystyle \frac{{\cos }^2\theta }{b^2}}\\ D=& -2\left[{\displaystyle \frac{\cos \theta (x_0\cos \theta +y_0\sin \theta )}{a^2}}+{\displaystyle \frac{\sin \theta (x_0\sin \theta -y_0\cos \theta )}{b^2}}\right]\\ E=& -2\left[{\displaystyle \frac{\sin \theta (x_0\cos \theta +y_0\sin \theta )}{a^2}}-{\displaystyle \frac{\cos \theta (x_0\sin \theta -y_0\cos \theta )}{b^2}}\right]\\ F=& {\displaystyle \frac{(x_0\cos \theta +y_0\sin \theta {)}^2}{a^2}}+{\displaystyle \frac{(x_0\sin \theta -y_0\cos \theta {)}^2}{b^2}}-1\end{array}\end{array}$$

---

代码

在 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 六个参数。

---

后记

测试后发现一个问题，计算椭圆方程的方法计算比较复杂，而且计算结果受椭圆大小的影响，难以直观地反应点的偏差值。后来我发现其实可以利用“椭圆上的点到其两个焦点的距离之和不变”这一性质来判断一个点偏离椭圆的程度。

1. 椭圆的两个焦点在其长轴上，且与中心的距离为 a2+b2" role="presentation">a2+b2−−−−−−√$\sqrt{a^2+b^2}$;
2. 椭圆上的点到两个焦点的距离之和为 2a" role="presentation">2a$2a$.

计算椭圆焦点的函数为：

std::vector<cv::Point2f> calcFocal(cv::RotatedRect& ellipse)
{
    if (ellipse.size.width < ellipse.size.height) {
        float tmp = ellipse.size.width;
        ellipse.size.width = ellipse.size.height;
        ellipse.size.height = tmp;
        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));
    std::vector<cv::Point2f> rst;
    rst.push_back(ellipse.center + delta);
    rst.push_back(ellipse.center - delta);
    return rst;
}

对于任意一个被拟合的点，只要计算其到两焦点的距离之和与 2a" role="presentation">2a$2a$ 的差值，即可知道它与椭圆的像素偏差大小。

相关阅读

以贴吧和头条为例，为什么产品都有极速版和标准版

当用户需要从应用市场上下载一款软件的时候，往往可以发现不止一个版本。没有任何标志的为标准版本，而除此之外往往还有一个极速版本

标准的液压传动系统怎么构成？工作原理是什么？

液压传动中由液压泵、液压控制阀、液压执行元件(液压缸和液压马达等)和液压辅件(管道和蓄能器等)组成的液压系统。液

如何利用matlab求解方程

如何利用matlab求解方程1.    前言作为三大数学软件之一，matlab在数值计算方法的能力首屈一指。求解方程是工科学习和工程计算

Normalization(标准化)的原理和实现详解

本系列文章为原创，转载请注明出处。 作者：Tom Bai 邮箱： [email protected] 若您觉得本博文对您有帮助，请您为我点赞并关注我，以鼓

微分方程模型

在数学建模中，当要描述某事件的数量变化对时间或者其它东西的变化时，可以考虑采用微分方程模型。根据函数及其变化率之间的关系确定