こんにちは、つまみ (@TrpFrog) です。生まれてから一度も歩いたことがありません。

大学の課題でタイトルの内容をどうしても計算したくなる機会があったのですが僕のググり力不足かヒットせず、 泣く泣く自力で計算する羽目になったので覚書を残します。 数学的な議論が怪しいのはご容赦ください。(は？)(ツッコミがあればお願いします)

あとProcessing.jsを使ってみよう！ ということで、記事の最後に触って確かめられるProcessingのコードを載せました。

前提知識

この記事における前提知識は

• 三角関数 ($\sin, \cos, \tan, \mathrm{arc\hspace{0px}tan})$
• 高校数学レベルのベクトル
• 行列の積と回転行列

です。結果を知りたいだけならばこの知識は不要です。

解きたい問題

次のような問題が解きたいです。

問題

$xy$ 平面上に点 $A$ から点 $B$ 方向に伸びる半直線 $L$ と、点 $P$ があります。 このとき、点 $P$ から最も近い $L$ 上の点 $Q$ を求めてください。

直線で考える

さて、いきなり半直線で考えるのは大変そうなので直線で考えてみましょう。すると次のような問題になります。

この問題は簡単です。直線 $AB$ に向けて点 $P$ から下ろした垂線とその交点が $R$ となります。さてそのような $R$ を実際に求めてみましょう。

ベクトル方程式で表す

原点を $O$ としたとき, $\vec{a} = \overrightarrow{OA}, \vec{b} = \overrightarrow{OB}, \vec{v} = \overrightarrow{AB}$ とします。また、直線 $AB$ 上の点 $R$ の位置ベクトルを $\vec r$ とします。このとき直線 $AB$ のベクトル方程式はパラメータ $t$ を用いて次のように表すことができます。

また、$\vec v$ の同じ大きさの法線ベクトルを $\vec u$ とします。すなわち

となるようなベクトル $\vec u$ です。このとき、平面上の任意の点 $P$ の位置ベクトル $\vec p$ はパラメータ $s, t$ を用いて次の式で表すことができます。

これらの情報から点 $R$ を求めてみます。$R$ は直線 $AB$ に向けた垂線と直線の交点でした。$\vec u$ と $\vec v$ は直交しますから、求める点 $R$ の位置ベクトル $\vec r$ は式 (2) におけるパラメータをそのまま使って次のように表すことができます。

これで問題(2) の答えがわかりました！ めでたしめでたし……ではなく $s, t$ を求める必要があります。

Pの座標からパラメータを求める

式(2) を変形して次のような式を作ります。

$\vec u$ と $\vec v$ は直行していますから、$\vec p - \vec a$ を $\vec v$ が $x$ 軸と平行になるように回転してあげて、そのときの $x$ 座標を $|\vec v|$ で割った値を読むと $u$ がわかります。これは $\vec v$ の $x$ 軸となす角 (反時計回りを正とする) を $\theta$ として次のように表すことができます。

従って $R$ を求めるのに必要だったパラメータ $s$ は

であるとわかりました。

$\vec v$ の $x$ 軸となす角 $\theta$ を求める

次に $s$ を求めるのに必要な $\theta$ を計算しましょう。$\vec v = \overrightarrow{AB} = \vec b - \vec a$ であったことを思い出してください。つまり直線 $AB$ の傾きが求まれば良いです。ここで $A(x_A, y_A), B(x_B, y_B)$ とすると

です。ただし $x_B - x_A \neq 0$ 。$x_B - x_A = 0$ の時は今回は議論しません。ここについては各プログラミング言語の標準ライブラリとして存在するはずの便利関数 $\mathrm{atan\hspace{0px}2}(y, x)$ を用いて

を使ってください。(は？)

以上の議論より、問題(2)の答えの点 $R$ の位置ベクトル $\vec r$ は

であると分かりました。

半直線の問題を解く

さて、ここまでで元の問題の9割は解き終わりました。ここで元の問題についてもう一度考えてみましょう。

問題

$xy$ 平面上に点 $A$ から点 $B$ 方向に伸びる半直線 $L$ と、点 $P$ があります。 このとき、点 $P$ から最も近い $L$ 上の点 $Q$ を求めてください。

結論から述べると、点 $Q$ の位置ベクトルを $\vec q$ としてこの問題の答えは次の通りです。

$\vec v$ の係数を $\max(0, s)$ としただけです。なぜこうなるかは簡単です。パラメータ $s$ は直線のときと同様の議論をすることで

を満たすような $s$ であることがわかります。このとき半直線 $AB$ に乗らないように、$R$ を点 $A$ から点 $B$ と反対方向に動かしていくと $s < 0$ となってしまいます。よって負の数だけをカットすれば良く、$\vec v$ の係数は $\max(0, s)$ であるとわかります。

本当にそうなるの？

数学が下手くそな人間が数字をこねくり回してるだけでは信憑性に欠けるので(？)確認用のProcessingのコードを置いておきます。また、下の画像をマウスオーバーすることでそのような点 $Q$ が描画されます。実際に上で示した式が正しいことを動かして確認してみてください。

PVector A, B;

void halfStraightLine(float x1, float y1, float x2, float y2) {
    final int INF = 10000;
    line(x1, y1, x1 + (x2 - x1) * INF, y1 + (y2 - y1) * INF);
}

PVector calcQ(float x, float y) {
    float theta = atan2(B.y - A.y, B.x - A.x);
    float d = dist(A.x, A.y, B.x, B.y);
    float s = (cos(theta) * (x - A.x) + sin(theta) * (y - A.y)) / d;
    x = A.x + Math.max(0, s) * (B.x - A.x);
    y = A.y + Math.max(0, s) * (B.y - A.y);
    return new PVector(x, y);
}

void setup() {
    size(300, 400);
    A = new PVector(100, 100);
    B = new PVector(200, 300);
}

void draw() {
    background(245);

    stroke(0);
    halfStraightLine(A.x, A.y, B.x, B.y);
    noStroke();

    final PVector p = new PVector(mouseX, mouseY);
    final PVector q = calcQ(p.x, p.y);
    final int r = 10;

    textSize(16);
    fill(200);
    text("A", A.x + r, A.y + r);
    text("B", B.x + r, B.y + r);
    ellipse(A.x, A.y, r, r);
    ellipse(B.x, B.y, r, r);

    fill(0);
    text("P", p.x + r, p.y + r);
    text("Q", q.x + r, q.y + r);

    fill(#dbbb92);
    ellipse(p.x, p.y, r, r);

    fill(#90e200);
    ellipse(q.x, q.y, r, r);
}