CodeIQ：円と線分の位置関係
私自身が表題の問題を解いた時のプログラムについて解説します。
問題の詳細は「円と線分の位置関係」（CodeIQ）を参照してください。
問題の概要
円と線分があります。
円と線分の関係は、下表のいずれかになります：
記号 説明 状態の例 A 円が線分を完全に含んでいます B 線分の端の一方が円周上に、他方は円内にあります。 C 線分の端の両方が円周上にあります。 D 線分の端点のひとつは円周上にあり、そこから円内を通って円の外まで突き抜けます。 E 線分の端点は両方とも円外にありますが、円の中を通ります。 F 線分の端点の一方は円内に、他方は円外にあります。 G 線分の端点のひとつが円周上にあり、それ以外は円外です。 H 線分は、端点以外の１点で、円に接しています。 I 線分と円は共有点を持ちません。
どの関係になるのかを求めるプログラムを書いてください。
【入出力】
入力は
(1,1)1/(20,20)(20,21) (10,20)10/(100,100)(100,200) (11,22)100/(22,33)(25,35)
のような感じです。
「(中心のX座標,中心のY座標)半径/(端点のX座標,端点のY座標)(端点のX座標,端点のY座標)」という形式で書かれた円と線分の組が、空白区切りで並んでいます。
出力は、各ペアが前述の A〜I のどれに該当するのかを、区切り文字なしで並べたものです。
先程の例だと
最初のペア(1,1)1/(20,20)(20,21)は、共有点を持たないので I
二番目のペア(10,20)10/(100,100)(100,200)も、共有点を持たないので I
最後のペア(11,22)100/(22,33)(25,35)は、線分が円内に含まれているのでA
ということで、
IIA
と出力すれば正解となります。
末尾の改行はあってもなくても構いません。
私のプログラム
Rubyで解答しています。
#!/usr/bin/ruby
# epsilon
Epsilon = Float::EPSILON * 1000
# 円クラス
class Circle
def initialize(x, y, r)
@x = x
@y = y
@r = r
end
attr_accessor :x, :y, :r
end
# 線クラス
class Line
def initialize(x1,y1,x2,y2)
@x1 = x1
@y1 = y1
@x2 = x2
@y2 = y2
end
attr_accessor :x1, :y1, :x2, :y2
end
# 入力値をパースする
# [[a1,a2],[b1,b2], ...]でa,bはCircleのインスタンス
def parseInput(line)
arr = line.split(" ")
ret = []
for pair in arr
vals=[]
splited = pair.split("/")
# 円
x,y,r = splited[0].gsub("(", "").gsub(")", ",").split(",").map{|v| v.to_i}
vals.push(Circle.new(x,y,r))
# 線
splited[1].match(/\((\d+),(\d+)\)\((\d+),(\d+)\)/){|md|
vals.push(Line.new(md[1].to_i, md[2].to_i, md[3].to_i, md[4].to_i))
}
ret.push(vals)
end
return ret
end
# cの中心からcとlの接点までの線とlの成す角度を求める
# 円と線が接していなければならない
# 返り値はcosθの値なので-1から1の間になる。
# 負数が帰った場合はcの中心とlの両端が鈍角、正の値が帰った場合は鋭角を成す
def angle(c, l)
ol1 = (c.x - l.x1)**2 + (c.y - l.y1)**2 # 円の中心と線の左端の距離の2乗
ol2 = (c.x - l.x2)**2 + (c.y - l.y2)**2 # 円の中心と線の右端の距離の2乗
if ol1 == c.r**2 then
co = [(c.x-l.x1), (c.y-l.y1)] # 接点から円の中心までのベクトル
ce = [(l.x2-l.x1), (l.y2-l.y1)] # 接点から線のもう一端までのベクトル
elsif ol2 == c.r**2 then
co = [(c.x-l.x2), (c.y-l.y2)] # 接点から円の中心までのベクトル
ce = [(l.x1-l.x2), (l.y1-l.y2)] # 接点から線のもう一端までのベクトル
else
puts "Angle Error"
return nil
end
d = (co[0]*ce[0] + co[1]*ce[1])/Math.sqrt((co[0]**2 + co[1]**2) * (ce[0]**2 + ce[1]**2))
return d
end
# 円の中心から線に伸ばした垂線の長さを求める
def vline(c,l)
oa = Math.sqrt((l.x1 - c.x)**2 + (l.y1 - c.y)**2)
ob = Math.sqrt((l.x2 - c.x)**2 + (l.y2 - c.y)**2)
ab = Math.sqrt((l.x2 - l.x1)**2 + (l.y2-l.y1)**2)
area = Math.sqrt((oa+ob+ab)*(ob+ab-oa)*(oa-ob+ab)*(oa+ob-ab))/4
return 2*area/ab
end
# 線の延長常に円の中心がある場合、線が円の中心を通っているかを判定する
# ただし、線の両端は円の外になければならない
# 線が円の中心を通っている場合true、通っていない場合falseを返す
def isLineAcrossCenter(c, l)
oa = Math.sqrt((l.x1 - c.x)**2 + (l.y1 - c.y)**2)
ob = Math.sqrt((l.x2 - c.x)**2 + (l.y2 - c.y)**2)
ab = Math.sqrt((l.x2 - l.x1)**2 + (l.y2-l.y1)**2)
# 円の中心から線の両端それぞれまでの線分の長さの和が線の長さよりも長い場合、線は円の外側にある
if ((oa + ob) - ab).abs >= Epsilon then true
else return false
end
end
# 円と直線の距離をチェックする
# cはCircleのインスタンス
# lは線のインスタンス
def checkDistance(c, l)
ol1 = (c.x - l.x1)**2 + (c.y - l.y1)**2 # 円の中心と線の左端の距離の2乗
ol2 = (c.x - l.x2)**2 + (c.y - l.y2)**2 # 円の中心と線の右端の距離の2乗
# 円の中心から線の両端までの距離が円の半径よりも小さい（A）
if (ol1 < c.r**2) && (ol2 < c.r**2) then
return "A"
# 円の中心から線の一方の端までの距離が円の半径と同じて他方が半径よりも小さい（B）
elsif ((ol1 == c.r**2) && (ol2 < c.r**2)) || (ol1 < c.r**2) && (ol2 == c.r**2)
return "B"
# 円の中心から線の一方の端までの距離が円の半径と同じて他方が半径よりも大きい（D,G）
elsif ((ol1 == c.r**2) && (ol2 > c.r**2)) || (ol1 > c.r**2) && (ol2 == c.r**2)
if angle(c,l) > 0 then return "D" # 鋭角を成す場合は円を横切る
else return "G" # 直角よりも鈍角の場合は円を横切らない
end
# 円の中心から線の一方の端までの距離が円の半径より小さく他方が半径よりも大きい（F）
elsif ((ol1 < c.r**2) && (ol2 > c.r**2)) || (ol1 > c.r**2) && (ol2 < c.r**2)
return "F"
# 円の中心から線の両方の端までの距離が円の半径と同じ（C）
elsif (ol1 == c.r**2) && (ol2 == c.r**2) then
return "C"
# 円の中心から線の両方の端までの距離が円の半径より長い(E,H,I)
else
d = vline(c,l)
if d <= Epsilon then # 線と円の中心の間の距離が0（中心を通る）
if isLineAcrossCenter(c,l) then #円の中心と線の両端までの線分の長さの和と線の長さが等しい
return "I"
else
return "E"
end
else
if (d - c.r).abs <= Epsilon then # 線と円の中心が円の半径に等しい
return "H"
elsif (d - c.r) > Epsilon then # 線と円の中心が円の半径より遠い
return "I"
else # 線と円の中心が円の半径より近い
return "E"
end
end
end
end
#main
while line = gets
params = parseInput(line.strip)
end
ret = ""
for cl in params
ret += checkDistance(cl[0], cl[1])
end
puts ret
解説
一つ前の問題「２つの円の位置関係」に比べるとかなり難しいと思います。ですが高校レベルの数学の知識でできるので頑張れば回答できます。
CircleクラスとLineクラス
円の中心座標と半径を保持するクラスと線の両端の座標を保持するクラスです。メンバ関数を1つも持っていないのでC言語で言うところの構造体的な使い方をします。Pythonのようにタプルがあればタプルを使いましたがRubyの場合、配列で表現しなければならず[カッコ]がいっぱいになるのでクラスを作りました。
checkDistance()
円と線分の関係をチェックする関数でプログラムの本体です。
まず、後で計算に使うため円の中心と線分の両端までの距離の2乗の値を計算します。2乗の値を使うのは整数のまま計算できるようにです。この問題では三角関数を使ったりしているので結局浮動小数点を使っているのですが、浮動小数点誤差の問題は結構面倒臭いのです。なので私は整数で計算できるところはなるべく整数で計算するようにしています。
あとは場合分けして解いて行きます。
円の中心から線の両端までの距離が円の半径よりも小さい場合、線分は円の中にあるのでAを返します。
円の中心から線の一方の端までの距離が円の半径と同じて他方が半径よりも小さい場合、半径と同じ距離の端で線分は円と接しており、小さい方の端は円の内側にあるのでBを返します。
円の中心から線の一方の端までの距離が円の半径と同じて他方が半径よりも大きい場合は2パターン考えられます。DとGのいずれの場合も円と線分が接している点と円の中心の距離は円の半径と同じですが、反対側は円の半径より大きくなります。
線分が円の内側に入っているかどうかを判断する必要があるのですが、これは円の中心をO、円と線分が接する点をA、線分の反対側をBとすると∠OACが鋭角なら円の内側を通り、90度以上開いていれば円の中を通りません。
これをangle()で計算します（後述）。
angle()はcosθの値を返すので正の値なら鋭角なのでD、0以下なら90度以上の鈍角なのでGを返します。
円の中心から線の一方の端までの距離が円の半径より小さく他方が半径よりも大きい場合、距離が小さい方は縁の内側にあり、大きい方は外側にあるのでFを返します。
円の中心から線の両方の端までの距離が円の半径と同じ場合は線分の両端が円周上にあるのでCを返します。
ここまででパターンわけできなかったものがE、H、Iになるのですが、これが面倒です。基本的な考え方は縁の中心から線分（の延長上）に垂直な線を引いた場合、その線の円の中心から線分と交わる点までの長さが円の半径より長いか（円の外側を通っている）、短いか（円の内側を通っている）、同じか（円と接している）で判断します（141〜148行）。ただし、線分（の延長線上）に円の中心がある場合、円の中心と線分までの距離が0になるので別に考えなければなりません（135〜140行）。
まず（141〜148行）の方を考えます。
vline()関数（）後述は線分と円の中心の距離を求める関数です。先に述べた通り、この距離が長さが円の半径より長い場合は円の外側を通っているのでI、短いか場合は円の内側を通っているのでE、同じか場合は円と接しているのでIを返します。
vline()が返す値はどうしても浮動小数点になります。なので比較は浮動小数点誤差を考慮しなければならないので半径と円の中心から線分までの距離の差をとってEpsilonと比較します。Epsilonは組み込み定数Float::EPSILONから適当に決めた値で、これはローカルで試験して決定しました。浮動小数点誤差は掛け算や割り算をすると蓄積するのでFloat::EPSILONより十分大きな値にすべきなのですが、妥当な値の決め方はない（実際にテストして決めるしかない）ようなのです。これが面倒なので私はできる限り浮動小数点で計算したくないのです。
（135〜140行）の方は円の中心と線の両端までの線分の長さの和と線の長さが等しいかどうかで判断できます。この条件に入ってきた場合、線分は円の中心を通っているので円の中心から線分の両端までの距離が等しい場合、線分は円の内側を通っているためE、そうでなければIです。
これを判断するのがisAcrossCenter()です（後述）が、コメントも関数名も間違っていて中心を通っていない場合true、中心を通っている場合falseを返します。
以上で全パターンを網羅したので入力値のペアに対して本関数を呼び出してやれば結果がわかります。
angle()
これは高校でやった内積で求めることができます。ベクトルAOをa、ベクトルABをbとするとcosθ = a・b / |a||b|と言うやつです。この関数はまさしくこの計算を行います。
まず、線分のどちらの点で円と接しているかを判断し（59〜60行目と62行目、65行目）、接している点から縁の中心方向のベクトル（変数co）と線の反対側の点まで（変数ce）を求めます。そして内積を計算し（73行目）その値を返します。
内積の値はcosθで角度ではありませんが、-1〜0の値を取り、負数なら鈍角、0なら直角、正数なら鋭角を意味するのでわざわざ角度にする必要がありません。cosθのまま使えば±符号と0かどうかだけを判断すれば良いため浮動小数点誤差を気にしなくて良いのが優れています。
vline()
円の中心から線分（の延長線上）まで垂線を引いた場合、その垂線の長さを求める関数です。
高校の教科書なら線分の方程式と円の中心座標から垂線の方程式を求め、その後線分の方程式と垂線の方程式から交点の座標を求め、その座標と円の中心の座標から長さを求めることになりますが複雑すぎます。
人間だとこの方法の方が計算しやすいですが、コンピュータの場合単に計算でできるならその計算が人間にとっては面倒くさくてもロジック的にはそちらの方が楽なことは多いです。
なので、この関数では三角形の面積から垂線の長さを計算しています。
まず、円の中心Oと線分の両端a,bの座標はわかっているので線分oa、ob、abの長さは計算できます（80〜82行目）。3辺の長さがわかればヘロンの公式で面積がわかります（84行目）。
そして三角形の底辺をabとすると高さは2×面積÷abで、高さ＝円の中心から線分までの垂線の長さになります（85行目）。
isAcrossCenter()
先に書いた通り、関数名もコメントも実装の結果と逆になっています。
この関数は線分が円の中心を通っているか、通っていない（延長線上に円の中心がある）のかを判断します。
ロジックは簡単です。線分が実際に円の中心を通っている場合、円の中心から線分のそれぞれの端までの長さ|oa|と|ob|の和が|ab|と同じです。線分が実際に円の中心を通っていなくて線分の延長線上に円の中心がある場合、|oa|+|ob|>|ab|です。これを計算して結果を返すだけです。
雑感
まさしく高校数学でやったような問題です。ベクトルの内積とかヘロンの公式とかの細部は忘れてしまっていたので調べなおしましたが、やり方自体は全て自力で考えつきました。
この問題、うまくやれば浮動小数点を使わずにできるような気がするのですがどうでしょうか？