学习内容：求一个点到椭球面的距离（上）

写在前面

最近的研究中大量的用到了椭圆、椭球等东西,而且我在学习过程中发现中文网络上居然查不到“点到椭球面的距离”这种看起来非常理所当然大家都会的东西,甚至连个推导过程都没有.在谷歌上查到了一篇相关的文章,在学习过程中,想把心得体会等记录下来,写在这里.那篇文章的链接先放在这里了.
下载原pdf

点到椭圆的距离

椭圆的一些基础知识

一个二维的普通椭圆可以由一个中心点$\mathbf{C}$,一组椭圆轴方向的标准正交的向量 { U 0 , U 1 } \{\mathbf{U}_0,\mathbf{U}_1\} {U0​,U1​},以及$e_i$表示(其中$e_0\ge e_1>0$).椭圆上的所有点可以表示为

$$\mathbf{P}=\mathbf{C}+x_0{\mathbf{U}}_0+x_1{\mathbf{U}}_1\text{\hspace{1em}(1)}$$

其中

$$(\frac{x_0}{e_0}{)}^2+(\frac{x_1}{e_1}{)}^2=1\text{\hspace{1em}(2)}$$

由于轴方向是正交的,则有$x_i={\mathbf{U}}_i\cdot (\mathbf{P}-\mathbf{C})$,带入(2)式中,可以得到
$$(\mathbf{P}-\mathbf{C}{)}^{T}M(\mathbf{P}-\mathbf{C})=1\text{\hspace{1em}(3)}$$

其中$M=R^{T}DR$,$R$是正交矩阵,每一列分别为 { U 0 , U 1 } \{\mathbf{U}_0,\mathbf{U}_1\} {U0​,U1​};$D$是一个对角阵,对角元素为$\frac{1}{e_0^2},\frac{1}{e_1^2}$.

想求出点到椭圆的距离,在椭圆对应的坐标系下进行求解就很重要.也就是说将任一点$Q$表示为$\mathbf{Q}=\mathbf{C}+y_0{\mathbf{U}}_0+y_1{\mathbf{U}}_1$的形式.由方程(3)定义的点$Q$到点$P$ 的距离,与由椭圆标准方程定义的$Y=(y_0,y_1)$到$X=(x_0,x_1)$的距离是相等的.

再利用椭圆的对称性简化我们的问题（这里就不再赘述）.我们只考虑在第一象限的点.

我们的问题是找到点到椭圆的距离,那么也就是要在椭圆上找到到目标点最近的点.

一个重要发现:椭圆上最近点处的法线指向目标点

标准椭圆的参数化表达式为$X(\theta )=(e_0\cos \theta ,e_1\sin \theta ),\theta \in \left[0,2\pi \right)$,点$\mathbf{Y}$ 到椭圆上任意一点的距离可表示为:

$$F(\theta )=|\mathbf{X}(\theta )-\mathbf{Y}{|}^2\text{\hspace{1em}(4)}$$

这是一个非负、周期性、可微的函数,因此它一定会有一个全局最小值,此时的一阶导数值为0,即

$$F^{\prime }(\theta )=2(\mathbf{X}(\theta )-\mathbf{Y})\ast {\mathbf{X}}^{\prime }(\theta )=0\text{\hspace{1em}(5)}$$

这是两个向量的内积,那么也就是说这两个向量是垂直的.${\mathbf{X}}^{\prime }(\theta )$是$\mathbf{X}(\theta )$处的切线,则说明点$\mathbf{Y}$ 到最近的点$\mathbf{X}$所连接成的向量就正是这里的曲线的法线.

有了这一发现,我们再利用椭圆的隐式表达式:$G(x_0,x_1)=(x_0/e_0{)}^2+(x_1/e_1{)}^2-1=0$,$G(x_0,x_1)$处的梯度就是点$(x_0,x_1)$ 处的法向量.则有$(y_0,y_1)-(x_0,x_1)=t\nabla G(x_0,x_1)/2=t(x_0/e_0^2,x_1/e_1^2)$, 也即

$$y_0=x_0(1+\frac{t}{e_0^2}),y_1=x_1(1+\frac{t}{e_1^2})\text{\hspace{1em}(6)}$$

如果$(y_0,y_1)$在椭圆外,则$t>0$;如果$(y_0,y_1)$在椭圆内,则$t<0$;如果$(y_0,y_1)$在椭圆上,则$t=0$,距离为0.

圆的情况

如果$e_0=e_1$,则椭圆为一个圆.原点$(0,0)$在圆上有无限个距离最近的点,显然最近距离为$e_0$. 非原点的任意一点$(y_0,y_1)$的距离最近的点为$(x_0,x_1)=e_0(y_0,y_1)/|(y_0,y_1)|$.
式(6)中的$t=-e_0^2+e_0\sqrt{y_0^2+y_1^2}$.

目标点为原点

$y_0=0,y_1=0$.式(6)化为$0=x_0\left(1+t/e_0^2\right),0=x_1\left(1+t/e_1^2\right)$, 考虑4个情况

• $x_0=0,x_1=0$,显然不在椭圆上,排除掉;
• $x_0=0,x_1\ne 0$,有$t=-e_1^2$,这个点为$(0,e_1)$
• $x_0\ne 0,x_1=0$,有$t=-e_0^2$,这个点为$(e_0,0)$
• $x_0\ne 0,x_1\ne 0$,有$t=-e_1^2$且$t=-e_1^2$,显然不成立,排除掉;

那么只有$(0,e_1)$和$(e_0,0)$这两个点满足条件,而显然$(0,e_1)$更近,则其距离为$d=e_1$

目标点在$y$轴上

$y_0=0,y_1>0$.式(6)化为$0=x_0\left(1+t/e_0^2\right),y_1=x_1\left(1+t/e_1^2\right)$, 考虑2 个情况

• $x_0=0$,为使其在椭圆上,这个点为$(0,e_1)$;
• $x_0\ne 0$,有$t=-e_0^2,y_1=x_1(1-e_0^2/e_1^2$.而由于$x_1\ge 0,e_0>e_1$,有$y_1<0$ 矛盾
  因此最近的点为$(0,e_1)$,距离为$d=|e_1-y_1|$

目标点在$x$轴上

$y_0>0,y_1=0$.式(6)化为$y_0=x_0\left(1+t/e_0^2\right),0=x_1\left(1+t/e_1^2\right)$, 考虑2 个情况

• $x_1=0$,为使其在椭圆上,这个点为$(e_0,0)$;距离的平方为

$$d_0^2=(y_0-e_0{)}^2\text{\hspace{1em}(7)}$$

• $x_1\ne 0$,有$t=-e_1^2$,把$t$代入可得$y_0=x_0(1-e_1^2/e_0^2)$.由于点$(x_0,x_1)$在椭圆上,$x_0=e_0^2{y}_0/(e_0^2-e_1^2)\le e_0$,可知道$y_0\le (e_0^2-e_1^2)/e_0<e_0.$因此这个情况只包含了$(y_0,0)$在椭圆中的情况.而当$y_0\ge (e_0^2-e_1^2)/e_0$时,最近的点一定为$(e_0,0)$.若$y_0\le (e_0^2-e_1^2)/e_0$,解得点$x_1=e_1\sqrt{1-(x_0/e_0{)}^2}$,距离的平方为

$$d_1^2={\left(x_0-y_0\right)}^2+x_1^2=e_1^2\left(1-\frac{y_0^2}{e_0^2-e_1^2}\right)\text{\hspace{1em}(8)}$$

两个候选点分别为$\left(e_0,0\right)\text{和}\left(x_0,x_1\right)$,我们需要判别哪个更接近$(y_0,0)$.

$$d_0^2-d_1^2=\frac{{\left(e_0{y}_0+e_1^2-e_0^2\right)}^2}{e_0^2-e_1^2}\ge 0\text{\hspace{1em}(9)}$$

因此$\left(x_0,x_1\right)$是更近的点.

总结:在第一象限中,距离$(y_0,0)$最近的点需要分类讨论:

1. $0<y_0<(e_0^2-e_1^2)/e_0$,最近的点为$(x_0,x_1)$其中$x_0=e_0^2{y}_0/(e_0^2-e_1^2),x_1=e_1\sqrt{1-(x_0/e_0{)}^2}$,距离为$d=\sqrt{{\left(x_0-y_0\right)}^2+x_1^2}=\sqrt{e_1^2\left(1-\frac{y_0^2}{e_0^2-e_1^2}\right)}$
2. $y_0\ge (e_0^2-e_1^2)/e_0$,最近的点为$(e_0,0)$,距离为$d=|y_0-e_0|$

目标点严格在第一象限

取$y_0>0,y_1>0$,则方程(6)的解为

$$x_0=\frac{e_0^2{y}_0}{t+e_0^2},x_1=\frac{e_1^2{y}_1}{t+e_1^2}\text{\hspace{1em}(10)}$$

由于只考虑第一象限,则有$x_0\ge 0,x_1\ge 0$,也即$t>-e_0^2,t>-e_1^2$,又由于$e_0>e_1$,则有$t>-e_1^2$.
将等式(10)代回到(2)中有

$$F(t)={\left(\frac{e_0{y}_0}{t+e_0^2}\right)}^2+{\left(\frac{e_1{y}_1}{t+e_1^2}\right)}^2-1=0\text{\hspace{1em}(11)}$$

此函数的定义域为$t\in (-e_1^2,\infty )$,距离目标点最近的点由此等式的根定义.

一阶导、二阶导为

$$F^{\prime }(t)=-\frac{2e_0^2{y}_0^2}{{\left(t+e_0^2\right)}^3}-\frac{2e_1^2{y}_1^2}{{\left(t+e_1^2\right)}^3},F^{\prime \prime }(t)=\frac{6e_0^2{y}_0^2}{{\left(t+e_0^2\right)}^4}+\frac{6e_1^2{y}_1^2}{{\left(t+e_1^2\right)}^4}\text{\hspace{1em}(12)}$$

在定义域上有$F^{\prime }(t)<0,F^{\prime \prime }(t)>0$,且

$$\underset{t\to -e_1^2+}{\lim }F(t)=+\infty ,\underset{t\to \infty }{\lim }F(t)=-1\text{\hspace{1em}(13)}$$
第一个表达式表示单侧极限($t\to -e_1^2,\text{且}t>-e_1^2$)
$F(t)$是一个关于$t$的严格递减函数,函数值首先为正、然后变成负数,因此必定有唯一零点.

因此只需要解出来这样的一个$\bar{t}$,代入到(6)中即可得出$(x_0,x_1)$,然后就可以求出来最小距离.

写在后面

要求点到椭球距离，需要先学习点到椭圆的距离，说实话，光打这篇文章就已经非常累了，还有些地方不是很能看得懂。
明天继续打后面的第二部分，也就是点到椭球面的距离。
而且在这里写博客的时候发现不会给方程标号，也不会引用。随后再研究研究，修改一下，尽量完善一下。