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

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

写在前面

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

点到椭球的距离

椭球的基本知识

一个普通椭球面可以由一个中心点$\mathbf{C}$,一组椭球轴方向的标准正交的向量 { U 0 , U 1 , U 2 } \{\mathbf{U}_0,\mathbf{U}_1,\mathbf{U}_2\} {U0​,U1​,U2​},以及$e_i$表示(其中$e_0\ge e_1\ge e_2>0$).椭圆上的所有点可以表示为
$$\mathbf{P}=\mathbf{C}+x_0{\mathbf{U}}_0+x_1{\mathbf{U}}_1+x_2{\mathbf{U}}_2\text{\hspace{1em}(25)}$$

其中
$$(\frac{x_0}{e_0}{)}^2+(\frac{x_1}{e_1}{)}^2+(\frac{x_2}{e_2}{)}^2=1\text{\hspace{1em}(26)}$$

由于轴方向是正交的,则有$x_i={\mathbf{U}}_i\cdot (\mathbf{P}-\mathbf{C})$,带入(\ref{26})式中,可以得到

$$(\mathbf{P}-\mathbf{C}{)}^{T}M(\mathbf{P}-\mathbf{C})=1\text{\hspace{1em}(27)}$$

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

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

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

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

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

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

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

$$\frac{\partial F}{\partial \theta }=2(\mathbf{X}(\theta ,\varphi )-\mathbf{Y})\cdot \frac{\partial \mathbf{X}}{\partial \theta }=0,\frac{\partial F}{\partial \varphi }=2(\mathbf{X}(\theta ,\varphi )-\mathbf{Y})\cdot \frac{\partial \mathbf{X}}{\partial \varphi }=0\text{\hspace{1em}(29)}$$

一阶导数为0,则$(\mathbf{X}(\theta ,\varphi )-\mathbf{Y})$一定与$\frac{\partial \mathbf{X}}{\partial \theta },\frac{\partial \mathbf{X}}{\partial \varphi }$垂直.这说明点$\mathbf{Y}$ 到最近的点$\mathbf{X}$所连接成的向量是这里的平面的法向量.

使用椭球的隐函数形式,$G(x_0,x_1,x_2)=(x_0/e_0{)}^2+(x_1/e_1{)}^2(x_2/e_2{)}^2-1=0$,$G(x_0,x_1,x_2)$处的梯度就是点$(x_0,x_1,x_2)$ 处的法向量.则有$(y_0,y_1,y_2)-(x_0,x_1,x_2)=t\nabla G(x_0,x_1,x_2)/2=t(x_0/e_0^2,x_1/e_1^2,x_2/e_2^2)$, 也即

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

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

球的情况

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

Oblate的情况

如果$e_0=e_1>e_2$,则椭球为一个Oblate球.取$r=\sqrt{x_0^2+x_1^2}$,标准椭圆方程在$(r,x_2)-$平面上退化成为$(x_0^2+x_1^2)/e_0^2+x_2^2/e_2^2=1$.用映射$(y_0,y_1,y_2)\to (\sqrt{y_0^2+y_1^2},y_2)$并在椭圆上讨论这一问题.

Prolate的情况

如果$e_0>e_1=e_2$,则椭球为一个Prolate球.取$r=\sqrt{x_1^2+x_2^2}$,标准椭圆方程在$(x_0,r)-$平面上退化成为$x_0^2/e_0^2+(x_1^2+x_2^2)/e_1^2=1$.用映射$(y_0,y_1,y_2)\to (y_0,\sqrt{y_1^2+y_2^2})$并在椭圆上讨论这一问题.

本部分主要讨论椭球的问题.因此假设$e_0>e_1>e_2$,并在第一卦限讨论.

目标点是原点

方程(30)变为$0=x_i\left(1+t/e_i^2\right),i=0,1,2.$.与椭圆的讨论类似,共有8种情况需要考虑.首先排除掉$x_0=x_1=x_2=0$和同时满足两个$t=-e_i^2,t=-e_j^2,i\ne j$的情况,只剩下三个情况$(e_0,0,0),(0,e_1,0),(0,0,e_2)$,显然最近距离为$d=e_2$.

目标点$y_2>0$

方程(30)变为$y_2=x_2\left(1+t/e_2^2\right).$ 由于$y_2>0$,寻找的最近点在第一卦限,必须有$x_2>0,1+t/e_2^2>0$.考虑$y_0,y_1$分别为0或者正时,共考虑4种情况

1. $y_0=0,y_1=0\text{则有}0=x_0\left(1+t/e_2^2\right),0=x_1\left(1+t/e_1^2\right),\text{如果}{x}_0>0,\text{则有}t=-e_0^2,y_2=x_2(1-e_0^2/e_1^2)<0,\text{矛盾,相似的,}{x}_1>0\text{矛盾,因此最近点为}(x_0,x_1,x_2)=(0,0,e_2)$

2. $y_0>0,y_1=0\text{则有}{y}_0=x_0\left(1+t/e_2^2\right),\text{则}{x}_1=0,\text{则问题退化为在椭圆上的问题,即求点}(y_0,y_2)\text{到}(x_0/e_0{)}^2+(x_2/e_2{)}^2=1\text{的距离}$

3. $y_0=0,y_1>0\text{则有}{y}_1=x_1\left(1+t/e_2^2\right),\text{则}{x}_0=0,\text{则问题退化为在椭圆上的问题,即求点}(y_1,y_2)\text{到}(x_1/e_1{)}^2+(x_2/e_2{)}^2=1\text{的距离}$

4. $y_0>0,y_1>0$,则方程(30})解为$x_i=e_i^2{y}_i/(t+e_i^2),i=0,1,2$.与前面椭圆的解法类似.将$x_i$带入(26)得到
   $$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+{\left(\frac{e_2{y}_2}{t+e_2^2}\right)}^2-1=0\text{\hspace{1em}(31)}$$
   此函数$F(t)$的定义域为$t\in (-e_2^2,\infty )$.椭球上最近一点$(x_0,x_1,x_2)$由$F(t)=0$的根给出.

   一阶、二阶导数分别为
   $$\begin{gathered} 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}-\frac{2e_2^2{y}_2^2}{{\left(t+e_2^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}+\frac{6e_2^2{y}_2^2}{{\left(t+e_2^2\right)}^4}\text{\hspace{1em}(32)} \end{gathered}$$
   由$y_0>0,y_1>0,y_2>0,t>-e_2^2,F^{\prime }(t)<0,F^{\prime \prime }(t)>0$,又由于${\lim }_{t\to -e_2^{2+}}F(t)=+\infty ,{\lim }_{t\to \infty }F(t)=-1$
   $F(t)$是一个关于$t$的严格递减函数,函数值首先为正、然后变成负数,因此必定有唯一零点.
   
   可以发现定义域中包含了0,而且$F(0)=(y_0/e_0{)}^2+(y_1/e_1{)}^2+(y_2/e_2{)}^2$,$(y_0,y_1,y_2)$在椭球内时$F(0)<0,\bar{t}<0$,在椭球外时$F(0)>0,\bar{t}>0$,在椭球上时$F(0)=0,\bar{t}=0$

目标点$y_2=0$

与椭圆中“点在x轴上”类似.方程(30)转化为$y_0=x_0(1+t/e_0^2),y_1=x_1(1+t/e_1^2),0=x_2(1+t/e_2^2)$.最后一个方程有两种情况

• $x_2=0$.为了使$(x_0,x_1,0)$在椭球上,则有$(x_0/e_0{)}^2+(x_1/e_1{)}^2=1$,降维至二维椭圆问题.
• $x_2>0$,则有$t=-e_2^2$,带入可得:$y_0=x_0(1-e_2^2/e_0^2),y_1=x_0(1-e_2^2/e_1^2)$.
  由于$(x_0,x_1,x_2)$在椭球上,$$1\ge 1-(\frac{x_2}{e_2}{)}^2=(\frac{x_0}{e_0}{)}^2+(\frac{x_1}{e_1}{)}^2=(\frac{e_0{y}_0}{e_0^2-e_2^2}{)}^2+(\frac{e_1{y}_1}{e_1^2-e_2^2}{)}^2\text{\hspace{1em}(34)}$$此不等式表明$(y_0,y_1)$包含在一个$\text{半椭球上的}{x}_0{x}_1\_\text{域子椭圆}$中
  $(y_0,y_1)$在子椭圆外,最近的点需要在椭圆$(x_0/e_0{)}^2+(x_1/e_1{)}^2=1$上
  $(y_0,y_1)$在子椭圆内,为了使$(x_0,x_1,x_2)$在椭球上,解得$x_2=e_2\sqrt{1-(x_0/e_0{)}^2-(x_1/e_1{)}^2}$

至此 讨论结束。