2直線の交点座標を求めるプログラム

相変わらず，この記事が参照されるようで。

以下のプログラムを見て分かる人はそもそも，こんな記事を読みに来るとも思えないので，一応，式を書いておきましょうかねぇ。

二点 (x1, y1) と (x2, y2) を通る直線と，もう二組みの点 (x3,y3) と (x4,y4) を通る直線の交点の座標 （x，y） は，まず，

a1 = (y2-y1)/(x2-x1)
a3 = (y4-y3)/(x4-x3)

を計算して，その後，

x = (a1*x1-y1-a3*x3+y3)/(a1-a3)
y = (y2-y1)/(x2-x1)*(x-x1)+y1

ということです。下のプログラムそのまんまですね。。。下の方に追記

(x1,y1)と(x2,y2)を通る直線と，(x3,y3)と(x4,y4)を通る直線の交点の座標を求める
prog <- function(x1, y1, x2, y2, x3, y3, x4, y4)
{
a1 <- (y2-y1)/(x2-x1)
a3 <- (y4-y3)/(x4-x3)
x <- (a1*x1-y1-a3*x3+y3)/(a1-a3)
y <- (y2-y1)/(x2-x1)*(x-x1)+y1
return(c(x=x, y=y))
}
以下は，利用法と説明図を描くためのプログラム
line2 <- function(x1, y1, x2, y2)
{
abline(y1-(y2-y1)/(x2-x1)*x1, (y2-y1)/(x2-x1), col="green")
}
x1 <- 4
y1 <- 8
x2 <- 8
y2 <- 0
x3 <- 2
y3 <- 2
x4 <- 12
y4 <- 7
(ans <- prog(x1, y1, x2, y2, x3, y3, x4, y4))
par(xpd=TRUE, mar=c(5, 5, 2, 5))
plot(c(x1, x2, x3, x4), c(y1, y2, y3, y4), pch=19, xlab="x", ylab="y", asp=1)
line2(x1, y1, x2, y2)
line2(x3, y3, x4, y4)
abline(v=ans["x"], h=ans["y"], col="blue", lty=3)
text(c(x1, x2, x3, x4), c(y1, y2, y3, y4), paste("(x", 1:4, ", y", 1:4, ")", sep=""), pos=4)

実行結果として　（6, 4） を得る

追記：ずいぶん前の記事にコメントが付いたので，プログラムの補完

prog = function(x1, y1, x2, y2, x3, y3, x4, y4) {
  if (x1 == x2 && x3 == x4) { # 追加1
    x = y = NA
  } else if (x1 == x2) { # 追加2
    x = x1
    y = (y4 - y3) / (x4 - x3) * (x1 - x3) + y3
  } else if (x3 == x4) { # 追加3
    x = x3
    y = (y2 - y1) / (x2 - x1) * (x3 - x1) + y1
  } else {
    a1 = (y2 - y1) / (x2 - x1)
    a3 = (y4 - y3) / (x4 - x3)
    if (a1 == a3) { # 追加4
      x = y = NA
    } else {
      x = (a1 * x1 - y1 - a3 * x3 + y3) / (a1 - a3)
      y = (y2 - y1) / (x2 - x1) * (x - x1) + y1
    }
  }
  return(c(x = x, y = y))
}

line2 = function(x1, y1, x2, y2) {
  if (x1 == x2) { # 追加5
    abline(v = x1)
  } else {
    abline(y1 - (y2 - y1) / (x2 - x1) * x1, (y2 - y1) / (x2 - x1), col = "green")
  }
}

check.program = function(x1, y1, x2, y2, x3, y3, x4, y4) {
  ans = prog(x1, y1, x2, y2, x3, y3, x4, y4)
  print(ans)
  par(xpd = TRUE, mar = c(5, 5, 2, 5))
  plot(c(x1, x2, x3, x4), c(y1, y2, y3, y4),
    pch = 19, xlab = "x", ylab = "y", asp = 1)
  line2(x1, y1, x2, y2)
  line2(x3, y3, x4, y4)
  abline(v = ans["x"], h = ans["y"], col = "blue", lty = 3)
  text(c(x1, x2, x3, x4), c(y1, y2, y3, y4),
    paste("(x", 1:4, ", y", 1:4, ")", sep = ""), pos = 4)
}

check.program(4, 8, 8, 0, 2, 2, 12, 7)
check.program(4, 1, 4, 5, -1, -1, 6, 3)
check.program(-1, -1, 6, 3, 4, 1, 4, 5)
check.program(-1, -1, 6, 3, 4, 1, 4+1e-3, 5)
check.program(1, 2, 5, 2, 0, 3, 8, 3)
check.program(1, 1, 5, 5, -1, 0, 6, 7)
check.program(-5, 2, -5, 8, 2, 3, 2, 7)