标准椭圆方程推导

初衷

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

---

推导

中心在原点且长轴在x轴上的椭圆方程为：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\text{\hspace{1em}(1)}$$其中，$a$和$b$分别是长半轴和短半轴。若椭圆长轴与正x轴夹角为$\theta$(逆时针)，且中心坐标为$(x_0,y_0)$。其任意一点平移$(-x_0,-y_0)$再旋转$-\theta$即满足标准椭圆方程。

首先考虑点$(x,y)$绕原点顺时针旋转角度$\theta$到$(x^{\prime },y^{\prime })$。
用参数方程表示点$(x,y)$和点$(x^{\prime },y^{\prime })$：
$$\{\begin{array}{l}x=t\cos \varphi \\ y=t\sin \varphi \end{array}\text{\hspace{1em}(2)}$$$$\{\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}\text{\hspace{1em}(3)}$$因此：
$$\{\begin{array}{l}{x}^{\prime }=x\cos \theta +y\sin \theta \\ y^{\prime }=-x\sin \theta +y\cos \theta \end{array}\text{\hspace{1em}(4)}$$继而可得一般椭圆方程为：
$$\frac{{\left[(x-x_0)\cos \theta +(y-y_0)\sin \theta \right]}^2}{a^2}+\frac{{\left[(x-x_0)\sin \theta -(y-y_0)\cos \theta \right]}^2}{b^2}=1$$化简为 $A{x}^2+Bxy+C{y}^2+Dx+Ey+F=0$ 的形式可得：
$$\begin{array}{rl}A=& \frac{{\cos }^2\theta }{a^2}+\frac{{\sin }^2\theta }{b^2}\\ B=& 2\sin \theta \cos \theta \left(\frac{1}{a^2}-\frac{1}{b^2}\right)\\ C=& \frac{{\sin }^2\theta }{a^2}+\frac{{\cos }^2\theta }{b^2}\\ D=& -2\left[\frac{\cos \theta (x_0\cos \theta +y_0\sin \theta )}{a^2}+\frac{\sin \theta (x_0\sin \theta -y_0\cos \theta )}{b^2}\right]\\ E=& -2\left[\frac{\sin \theta (x_0\cos \theta +y_0\sin \theta )}{a^2}-\frac{\cos \theta (x_0\sin \theta -y_0\cos \theta )}{b^2}\right]\\ F=& \frac{(x_0\cos \theta +y_0\sin \theta {)}^2}{a^2}+\frac{(x_0\sin \theta -y_0\cos \theta {)}^2}{b^2}-1\end{array}\text{\hspace{1em}(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 六个参数。

---

后记

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

1. 椭圆的两个焦点在其长轴上，且与中心的距离为 $\sqrt{a^2-b^2}$;
2. 椭圆上的点到两个焦点的距离之和为 $2a$.

计算椭圆焦点的函数为：

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};
}

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