6. 距離

距離の求めかた。

パラメータで表される $\vec{\bm{X}}(\vec{s})$ と $\vec{Y}(\vec{t})$ , $\vec{s}\in A\subset \mathbb{R}^{n}$ , $\vec{t}\in B\subset \mathbb{R}^{m}$ の二乗距離を求めたいなら、

$\mathcal{Q}(\vec{s},\vec{t})=\left|\vec{X}(\vec{s})-\vec{Y}(\vec{t})\right|^{2} , (\vec{s},\vec{t})\in A\times B\subset \mathbb{R}^{n}\times \mathbb{R}^{m}$

この関数の $A\times B$ での極小値を考えれば良い。これは連続で微分可能な関数なので

$\vec{\nabla}(\mathcal{Q})=0$ か、さもなくば $A\times B$ の境界での点ある。

6.1 点からリニアなコンポーネント

ここで言うリニアコンポーネントとは、直線,線分の両方を含めたものである。

以下の議論は3次元に限らずいかなる次元でも成り立つ。

点を $\vec{P}$ とする。線はパラメータを使って $\vec{L}(t)=\vec{B}+t\vec{M}$ と書く。 $\vec{B}$ は線上の点, $\vec{M}$ は線の方向。 $t\in \mathbb{R}$ 。

半直線(ray)も同じ形式で、 $t\geq 0$ という制限がつくだけだ。

線分(line segment)の場合、 $t\in[0,1]$ という制限になる。

$\vec{P}$ に最も近い線上の点は、 $\vec{P}$ のこの線への投影地点であるから、

$t_{0}=\displaystyle \frac{proj_{M}\vec{BP}}{M}=\frac{\vec{M}\text{・}(\vec{P}-\vec{B})}{\left|\vec{M}\right|^{2}}$

Pから線分までの距離は、

$d = \left|\vec{P}-\vec{L}(t_{0})\right|$

言うまでもなく半直線と点との距離を求めていて、 $t\geq 0$ という制限があるときに、 $t_{0}<0$ なる答えになったら、

$d = \left|\vec{P}-\vec{L}(0)\right|$

ですよ、と。以下、線分の時も同様。

$\left|\vec{M}\right|^{2}$ での割り算は結構な計算コストがかかるので、 $t_{0}\geq 0$ かそういうの先に確認してから、本当に必要になってから計算するのが常道である。

もし可能ならば $\vec{M}$ を単位ベクトルとなるように正規化しておいても良い。

6.2 リニアコンポーネントからリニアコンポーネント

$L_{0}: \vec{L}_{0}(s)=\vec{B}_{0}+t\vec{M}_{0}$ , $s\in I\subset \mathbb{R}$

$L_{1}: \vec{L}_{1}(t)=\vec{B}_{1}+t\vec{M}_{1}$ , $t\in J\subset \mathbb{R}$

の２つのリニアコンポーネント間の距離を考えてみる。

言うまでもなく $L_{0}$ は $I=\mathbb{R}$ ならば、直線 , $I=[0,\infty)$ ならばray , $I=[0,1]$ ならば線分。 $L_{1}$ に関しても同様。

二乗距離 $\mathcal{Q}(s,t)=\left|\vec{L}_{0}(s)-\vec{L}_{1}(t)\right|^{2} , (s,t)\in I\times J$

なので、これ展開して係数うざいんで置き換えたら

$\mathcal{Q}(s,t)=as^{2}+2bst+ct^{2}+2ds+2et+f$

こんな感じ。

s,tで偏微分 $\displaystyle \vec{\nabla}\mathcal{Q}=(\frac{\partial \mathcal{Q}}{\partial s}$ , $\displaystyle \frac{\partial \mathcal{Q}}{\partial t})=2(as+bt+d,bs+ct+e)$ が(0,0)となる点か、領域の境界かが最小。この前者の方程式をs,tについて解くには、行列表現

$\left(\begin{array}{l}a b\\b c\end{array}\right)\left(\begin{array}{l}s\\t\end{array}\right)=-\left(\begin{array}{l}d\\e\end{array}\right)$ を考え、両辺の左から $\left(\begin{array}{l}a b\\b c\end{array}\right)^{-1}$ を掛けるとてっとり早い。

$(s,t)=-\displaystyle \frac{1}{ac-b^{2}}\left(\begin{array}{l}c -b\\-b a\end{array}\right)\left(\begin{array}{l}d\\e\end{array}\right)=(\frac{be-cd}{ac-b^{2}},\frac{bd-ae}{ac-b^{2}})$

また $ac-b^{2} =0$ のときは逆行列は存在しないが、これは

$ac-b^{2}=\vec{M}_{0}^{2}\vec{M}_{1}^{2}-(\vec{M}_{0}\text{・}\vec{M}_{1})^{2}=\left|\vec{M}_{0}\times\vec{M}_{1}\right|^{2}$ で、これが0ということは、 $\vec{M}_{0}$ , $\vec{M}_{1}$ が平行であることを意味する。この場合は、 $L_{0}$ 上のどこか一点と $L_{1}$ との距離を求めると良い。これは $t=0,s=-d/a$ にしとけばシンプルでよい。

直線と直線との距離

上で求めたので割愛。

線からrayか線分

これも同様なので割愛。

rayからrayか線分,線分から線分

これ実は少し面倒。

例えば線分と線分との交点。線分は、 $(s,t)\in I\times J=[0,1]^{2}$ だが、

解 $(\displaystyle \overline{s},\overline{t})=(\frac{be-cd}{ac-b^{2}},\frac{bd-ae}{ac-b^{2}})$ が、 $(\overline{s},\overline{t})\in[0,1]^{2}$ なら何の問題もないが、解がこの範囲外だとしたらどうすれば良いのか？

たとえば、 $(\overline{s},\overline{t})=(2,0.5)$ だったとしよう。 $\overline{t}$ は範囲内だが、 $\overline{s}$ は範囲外である。

$\mathcal{Q}$ は、s,tによる2次方程式だから、(x,y,z)=(s,t, $\mathcal{Q}$ (s,t))の描く図形はパラボロイド(paraboloid:放物線を回転させたような図形)である。この図形をz= $\mathcal{Q}_{0}$ でスライスしたらこれは $(\overline{s},\overline{t})$ を中心とする楕円である。そのへんを考えて $\mathcal{Q}_{0}$ を大きくしていき、この楕円が大きくなっていくときに、 $I\times J$ と最初に接するのはどこかを考える。

これは、 $(1,t_{0}) , t\in[0,1]$ で、 $t_{0}\neq\overline{t}$ である。これが $\mathcal{Q}$ を最小化する境界の端点。楕円が傾いているかも知れないので $t_{0}\neq\overline{t}$ になる。 $\overline{s},\overline{t}$ のうち片方が領域からはみ出している場合については、これと同様。

$\overline{s},\overline{t}$ の両方が領域からはみ出している場合、たとえば、 $(\overline{s},\overline{t})=(2,3)$ のような場合は、 $(1,t_{0}) , t\in[0,1]$ と $(1,s_{0}),s\in[0,1]$ と、どちらかが $\mathcal{Q}$ を最小化する境界の端点。このどちらが解なるかと言うと．．

$\vec{\nabla}\mathcal{Q}=(\mathcal{Q}_{s},\mathcal{Q}_{t})$ を考えてみると、パラボロイドなので、(s,t)=(1,1)の地点では $\mathcal{Q}_{s}$ , $\mathcal{Q}_{t}$ の両方が負になることはない。

$\mathcal{Q}_{s}(1,1)>0$ なら $s<1$ , $\mathcal{Q}_{s}(s,1)<\mathcal{Q}_{s}(1,1)$ , なので $t=1$ 上に解。

$\mathcal{Q}_{t}(1,1)>0$ なら $t<1$ , $\mathcal{Q}_{t}(t,1)<\mathcal{Q}_{t}(1,1)$ , なので $s=1$ 上に解。

以下、 $\overline{s},\overline{t}$ の両方が領域からはみ出している場合は、これと同様。

そんなわけで全部で3×3=9の場合分け。

6.3 点から3角形

点: $\vec{P}$

3角形: $\vec{T}(s,t)=\vec{B}+s\vec{E}_{0}+t\vec{E}_{1}, (s,t)\in D=\{(s,t):s\in[0,1],t\in[0,1],s+t\leq 1\}$

との距離を求めてみる。

例によって二乗距離 $\mathcal{Q}(s,t)=\left|\vec{T}(s,t)-\vec{P}\right|^{2} , (s,t)\in D$

これ展開して係数うざいんで置き換えたら

$\mathcal{Q}(s,t)=as^{2}+2bst+ct^{2}+2ds+2et+f$

とまあ、さっきと同じ式が出てくるではありませんか。

じゃあ例によって、 $\vec{\nabla}\mathcal{Q}=(0,0)$ を解くと

解 $(\displaystyle \overline{s},\overline{t})=(\frac{be-cd}{ac-b^{2}},\frac{bd-ae}{ac-b^{2}})$

解 $(\overline{s},\overline{t})\in D$ なら何も問題はないんですけど、これが $D$ の外にある場合はどうなるのかと。

これは、例によって、パラボロイド考えて、そのスライスがどこの境界に接するか考えていけばok。

そんなわけで場合分けは7個。

6.4 リニアコンポーネントから三角形

直線： $\vec{L}(r)=\vec{B}_{0}+r\vec{M}_{0}$ , $r\in I\subset \mathbb{R}$

3角形: $\vec{T}(s,t)=\vec{B}+s\vec{E}_{0}+t\vec{E}_{1}, (s,t)\in D=\{(s,t):s\in[0,1],t\in[0,1],s+t\leq 1\}$

例によって二乗距離 $\mathcal{Q}(s,t,r)=\left|\vec{T}(s,t)-\vec{L}(r)\right|^{2} , (s,t,r)\in D\times I$

$=a_{00}s^{2}+a_{11}t^{2}+a_{22}r^{2}+2a_{01}st+2a_{02}sr+2a_{12}tr+2b_{0}s+2b_{1}t+2b_{2}r+c$

(係数は例によって適当に置き換えた)

$\vec{\nabla}\mathcal{Q}=(0,0,0)$ で最小。これは
$\left(\begin{array}{l}a_{00} a_{01} a_{02}\\a_{01} a_{11} a_{12}\\a_{02} a_{12} a_{22}\end{array}\right)\left(\begin{array}{l}s\\t\\r\end{array}\right)=-\left(\begin{array}{l}b_{0}\\b_{1}\\b_{2}\end{array}\right)$

となるので、逆行列使えば簡単に求まる。

この解 $(\overline{s},\overline{t},\overline{r})\in D\times I$ ならその点で最小。これがこの領域外だと、例によって、もよりの境界上で最小。

線から三角形

同様なので割愛。

半直線から三角形、線分から三角形

同様なので割愛。

6.5 点から四角形

点から3角形の時と同様。求まった解が領域外の時についての場合分けが増えるだけ。

6.6 リニアコンポーネントから四角形

リニアコンポーネントから三角形の時と同様。求まった解が領域外の時についての場合分けが増えるだけ。

半直線から四角形,線分から四角形

これも同様。

6.7 三角形から三角形

二乗距離 $\mathcal{Q}(s_{0},t_{0},s_{1},t_{1})=\left|\vec{T}_{0}(s_{0},t_{0})-\vec{T}_{1}(s_{1},t_{1})\right|^{2}$ , $(s_{i},t_{i})\in D for 0\leq i\leq 1$

$\mathcal{Q}$ の領域はデカルト積(Cartesian product) $D\times D\in \mathbb{R}^{4}$ になる。

場合分けは、3角形によって7つの場合分けが発生するので、7×7=49個の場合分けが必要。

6.8 三角形から四角形

これまた同様。 $\mathcal{Q}$ の領域はデカルト積 $D\times[0,1]^{2}\in \mathbb{R}^{4}$ 。7×9=63の場合分けが必要。

6.9 四角形から四角形

これまた同様。9×9=81の場合分けが必要。

6.10 点から方位性のあるボックス(point to oriented box)

ボックス：中心 $\vec{C}$ ,正規直交単位軸 $\vec{U}_{i}$ ,範囲 $e_{i}$

これらを用いて点を表す：

$\vec{P}=\vec{C}+s_{0}\vec{U}_{0}+s_{1}\vec{U}_{1}+s_{2}\vec{U}_{2}$

係数 $s_{i}$ は $s_{i}=\vec{U}_{i}(\vec{P}-\vec{C})$ で求まる。

$\vec{P}$ がボックス内部にあれば、距離は0。ボックス外部にあれば、この座標系 $(s_{0},s_{1},s_{2})$ において、方位性のないボックスと点との距離を求めるようにすれば求まる。

6.11 その他

点から楕円(point to ellipse)

考えるのは軸に沿った(axis−aligned)楕円だけで十分。

方位性をもった楕円(oriented ellipses)は回転させて考えれば良い。

楕円は円を軸方向に引き延ばしただけなので、それを考慮に入れると解きやすい。

点P： $(u,v)$

楕円: $(x/a)^{2}+(y/b)^{2}=1$

この楕円上でこの点へもっとも近い点を $(u,v)$ とする。

このとき、 $(x-u,y-v)$ は、楕円の法線である。 $\vec{\nabla}((x/a)^{2}+(y/b)^{2})=(x/a^{2},y/b^{2})$

$x-u=t*x/a^{2},v-y=t*y/b^{2} , t$ はパラメータ。

これを解いて $x=a^{2}u/(t+a^{2}),y=b^{2}v/(t+b^{2})$

んで、楕円の方程式にこれ代入
$\left(\frac{au}{ t+a^{2}}\right)^{2}+\left(\frac{bv}{ t+b^{2}}\right)^{2}=1$

$\text{[math]}$ 分母である $(t+a^{2})^{2}(t+b^{2})^{2}$ を両辺に掛けてみる。

$F(t)=(t+a^{2})^{2}(t+b^{2})^{2}-a^{2}u^{2}(t+b^{2})^{2}-b^{2}v^{2}(t+a^{2})^{2}=0$

この解が点Pに対して楕円上で一番近い点でのtの値を表している。ニュートン法で解くと良い。

$t_{i+1}=t_{i}-F(t_{i})/F^{\prime}(t_{i}) , i\geq 0$

初期値は(u,v)が楕円内部なら $t_{0}=0$ , 楕円外部なら $t_{0}=nax\{a,b\}\sqrt{u^{2}+v^{2}}$ が良い。

点から楕円体(point to ellipsoid)

楕円体: $(x/a)^{2}+(y/b)^{2}+(z/c)^{2}=1$

までの距離。以下、楕円の時と同様。

点から二次曲線or二次曲面(point to quadratic curve or quadric surface)

手順は点から楕円の時と同様。

一般的な2次方程式：

$\mathcal{Q}(\vec{X})=\vec{X}^{T}A\vec{X}+\vec{b}^{T}\vec{X}+c=0$

$A$ は $N\times N$ の対称行列 , $N=2 or 3$

$b$ は $N\times 1$ のベクトル。

$c$ はスカラー。

と表される。そこで。

平面: $\mathcal{Q}(\vec{X})=0$

点: $\vec{Y}$

との距離を調べてみる。この平面上の $\vec{Y}$ に最も近い点を $\vec{X}_{c}$ とする。

楕円の時同様、 $\vec{Y}-\vec{X}_{c}$ は $\vec{X}_{c}$ での平面の法線になっているはずである。平面の法線は $\nabla \mathcal{Q}(\vec{X})$ である。

$\vec{Y}-\vec{X}_{c}=t\nabla \mathcal{Q}(\vec{X}_{c})=t(2A\vec{X}+\vec{b}) , t$ はパラメータ

$\vec{X}_{c}=(I+2tA)^{-1}(\vec{Y}-t\vec{b})$

これを元の方程式に代入して、tについて解くのだけど、その前に $A=RDR^{T}$ と対角化してからのほうが良い。

$\vec{X}_{c}=(I+2tA)^{-1}(\vec{Y}-t\vec{b})=...=R(I+2tD)^{-1}R^{T}(\vec{Y}-t\vec{b})=R(I+2tD)^{-1}(\vec{\alpha}-t\vec{\beta})$

$\vec{\alpha}$ と $\vec{\beta}$ は簡便化のための置換。

2次元のとき: $(I+2tD)^{-1}=diag\{1/(1+2td_{0}) , 1/(1+2td_{1})\}$

3次元のとき: $(I+2tD)^{-1}=diag\{1/(1+2td_{0}) , 1/(1+2td_{1}) , 1/(1+2td_{2})\}$

この逆行列求めるときに2次元のときは、tについての4次方程式,3次元のときはtについての6次方程式になる。

3次元で点から円(point to circle in 3D)

点P: $\vec{P}$

いま問題とする円を含む平面を考える。

平面： $\vec{N}\text{・}\left(\vec{X}-\vec{C}\right)=0$ , 法線を表す単位ベクトル $\vec{N}$ , $\vec{C}$ は円の中心

この平面上の円から点までの距離を考える。

$\vec{N}$ と正規直交になる $\vec{U}$ , $\vec{V}$ をとり円を

$\vec{X}=\vec{C}+R(cos(\theta)\vec{U}+sin(\theta)\vec{V})=\text{：}\vec{C}+R\vec{W}(\theta)$ , $R$ は円の半径

と表す。

距離の二乗関数は $F(\theta)=\left|\vec{C}+R\vec{W}(\theta)-\vec{P}\right|^{2}=R^{2}+\left|\vec{C}-\vec{P}\right|^{2}+2R(\vec{C}-\vec{P})\text{・}\vec{W}$

この最小値を求める。

$0=F^{\prime}(\theta)=2R(\vec{C}-\vec{P})\text{・}\vec{W}^{\prime}(\theta)$

いま $\vec{W}\text{・}\vec{W}=1$ なのでこれを微分(積の微分法)すれば、 $\vec{W}$ ・ $\vec{W}^{\prime}=0$ がわかる。

また、 $\vec{W}^{\prime\prime}=-\vec{W}$ で $0=\vec{W}$ ・ $\vec{W}^{\prime}$ これを微分すれば $0=\vec{W}\text{・}\vec{W}^{\prime\prime}+\vec{W}^{\prime}\text{・}\vec{W}^{\prime}=-1+\vec{W}^{\prime}\vec{W}^{\prime}$