📐

線分同士の交点の求め方

2023/09/10に公開

Ruby

目的

次の線分 ab cd の交点を求めたい。

a = V[288.0, 104.96]
b = V[496.0, 29.44]
c = V[344.0, 26.88]
d = V[472.0, 108.8]

ab = b - a  # => (208.0, -75.52)
cd = d - c  # => (128.0, 81.92)
ca = a - c  # => (-56.0, 78.08)
ac = c - a  # => (56.0, -78.08)

直線 ab 上の t を求める

ベクトル ab を t 倍進めたら直線に接触するとする。

t は、a からすぐのところに交点があれば 0.1 ぐらいになるだろうし、もうすぐ b に到着するところに交点があれば 0.9 ぐらいになる。なので最初の図から t は 0.6 ぐらいになりそうだと見当がつく。

交点は e とする。

次に外積の公式を思い出す。二つのベクトルが平行だと外積は 0 になるので

cd × ce = 0 ^[1]

と定義する。次に ce は遠回りして ca と ae を足したものなので、

ce = ca + ae

と定義する。続いて ae は ab を 1.0 としたときの割合 t で求まるので、

ae = ab * t

と定義する。それを一つ前の式に代入すると

ce = ca + ab * t

となり、さらに前の式に代入すると

cd × (ca + ab * t) = 0

となる。ここから t = の式にする。内積と外積は分配法則 A × (B + C) = A × B + A × C が成り立つ^[2]ので、

(cd × ca) + (cd × ab * t) = 0cd × ab * t = -cd × cat = -cd × ca / cd × ab

となり実際に計算すると

t = -cd.cross(ca) / cd.cross(ab)  # => 0.5460122699386504

t が求まる。ベクトル ab の 0.55 あたりに交点 e があるのがわかった。正確な交点 e の位置は

a.lerp(b, t)  # => (401.57055214723925, 63.725153374233116)

で求まる。lerp は線形補間のことで

a + (b - a) * t  # => (401.57055214723925, 63.725153374233116)

と等しい。

t の値で相対的な位置がわかる

直線 ab 上の線分 ab を電車とし、交点 e を駅と考えたとき、t が 1.1 だと車両の先頭が駅にまだ来てないことがわかる。一方 -0.1 だと車両の最後がもう駅を出ちゃってる。そういうのが t だけでわかる。実際に図で見てみる。

来てない

来た

通過した

したがって直線ではなく線分上にあるかどうかを調べるには

(0..1).cover?(t)  # => true

とする。

一方の線分が飛び出ている場合

今度は線分 cd の方が飛び出ている場合を見てみる。

その場合、次のように直線 ab 上の t を調べても意味がない。

今度は直線 cd 上の t を求める。

直線 cd 上の t を求める

いったん綺麗な形で交差している図に戻して考える。

さっきまでは ab 上で考えていたが今度は cd 上で考える。c から e までは半分より少し短かいので t = 0.4 ぐらいだろうか。さきほどとほぼ同じ計算を行う。過程は同じなので結果だけ書くと、

t = -ab.cross(ac) / ab.cross(cd)  # => 0.4497699386503068
c.lerp(d, t)                      # => (401.57055214723925, 63.72515337423313)
(0..1).cover?(t)                  # => true

となる。

判定コードまとめ

ab 上の t

t = -cd.cross(ca) / cd.cross(ab)  # => 0.5460122699386504
(0..1).cover?(t)                  # => true
a.lerp(b, t)                      # => (401.57055214723925, 63.725153374233116)

cd 上の t

t = -ab.cross(ac) / ab.cross(cd)  # => 0.4497699386503068
(0..1).cover?(t)                  # => true
c.lerp(d, t)                      # => (401.57055214723925, 63.72515337423313)

両方の t が 0..1 の範囲のとき重なっている
交点はどちら側から計算しても同じになる

やる必要のない高速化

外積の場合、右辺と左辺を入れ替えると符号が反転するので、直線 ab 上の t を求める部分は、次のように変形できる。

t = -cd.cross(ca) /  cd.cross(ab)  # => 0.5460122699386504
t = -cd.cross(ca) / -ab.cross(cd)  # => 0.5460122699386504
t =  cd.cross(ca) /  ab.cross(cd)  # => 0.5460122699386504

すると直線 cd 上の t を求めるときの分母 ab.cross(cd) と一致するようになる。

t = -ab.cross(ac) / ab.cross(cd)  # => 0.4497699386503068

そこで ab.cross(cd) の結果を使いまわせば処理は1回で済むようになる。とはいえ、外積は四則演算によって一瞬なので共通化できるはずのコードをわざわざ展開してまでやるメリットはない。

途中から内積を使って判定する方法

上の方法では

ab 上の t が 0..1 の範囲か？
cd 上の t が 0..1 の範囲か？

を順に調べる方法だった。

それを少し変えて

ab 上の t が 0..1 の範囲か？
t から交点 e を求める
その交点 e は線分 cd に含まれるか？

とする方法^[3]もある。

具体的には

t = -cd.cross(ca) / cd.cross(ab)  # => 0.5460122699386504
(0..1).cover?(t)                  # => true

ここまでは同じで、ここからいきなり交点を求めて

e = a.lerp(b, t)  # => (401.57055214723925, 63.725153374233116)

その交点は直線 cd 上にあるので、線分 cd の両端から交点に対してベクトルを作り、

ce = e - c  # => (57.570552147239255, 36.84515337423312)
de = e - d  # => (-70.42944785276075, -45.07484662576688)

内積を調べる。

ce.dot(de)  # => -5715.451837555045

もし交点が線分 cd の中にあるなら、その2つのベクトルはお互いに向き合うため内積は負になる。(内積は同じ方向を向くほど正の大きなスカラー値を返し、同じ方向を向かないほど負の大きなスカラー値を返す)

ce.dot(de).negative?  # => true

なので負であれば線分 cd の中にも交点 e があると証明できる。

参照

コード

class App < Base
  def initialize
    super

    self.width, self.height = V[800, 128]

    a = window_wh * V[0.36, 0.82]  # => (288.0, 104.96)
    b = window_wh * V[0.62, 0.23]  # => (496.0, 29.44)
    c = window_wh * V[0.43, 0.21]  # => (344.0, 26.88)
    d = window_wh * V[0.59, 0.85]  # => (472.0, 108.8)

    @points.concat([a, b, c, d])

    @mode = 0
  end

  def button_down(id)
    super

    if id == Gosu::KB_Z
      @mode = @mode.next.modulo(3)
    end
  end

  def draw
    super

    a, b, c, d = @points
    ab = b - a
    cd = d - c

    if @mode == 0
      vputs "a: #{a.round}"
      vputs "b: #{b.round}"
      vputs "c: #{c.round}"
      vputs "d: #{d.round}"
    end

    if @mode == 1
      ca = a - c
      vector_draw(c, a, " ", " ")
      vputs "cd×ce = 0"
      vputs "cd×(ca + ab*t) = 0"
      # (cd × ca) + (cd × ab*t) = 0
      # cd × ab*t = -cd × ca
      t = -cd.cross(ca) / cd.cross(ab)
      vputs "t = -cd×ca / cd×ab = #{t.round(2)}"
      vputs "e = a.lerp(b, t) = #{a.lerp(b, t).round}"
      vector_draw(c, a.lerp(b, t), " ", "e")
    end

    if @mode == 2
      ac = c - a
      vector_draw(a, c, " ", " ")
      vputs "ab×ae = 0"
      vputs "ab×(ac + cd*t) = 0"
      # (ab × ac) + (ab × cd*t) = 0
      # ab × cd*t = -ab × ac
      t = -ab.cross(ac) / ab.cross(cd)
      vputs "t = -ab×ac / ab×cd = #{t.round(2)}"
      vputs "e = c.lerp(d, t) = #{c.lerp(d, t).round}"
      vector_draw(a, c.lerp(d, t), " ", "e")
    end

    vector_draw(a, b, "a", "b")
    vector_draw(c, d, "c", "d")
  end

  show
end

脚注

×は掛け算ではなく外積のこと ↩︎
クロス積→分配律 ↩︎
レイと線分とで当たり判定を取るで紹介されていた方法 ↩︎