기하 - 직선의 표현, 직선의 교점 구하기, 반평면 교집합

 

 

기하에 관하여 다루는 알고리즘은 문제의 답을 직접 구할 때와 원리는 비슷하나 구현에서 상당히 다른 측면이 있다. 부동소수점 오차를 최소화하고 예외(ex. 기울기가 0 혹은 무한대) 처리의 상황을 가급적 방지하고자 하기 때문이다. 때문에 대체로 선형대수학적 접근을 통해 구현한다. 이번에는 알고리즘에서 직선을 다루는 방식에 대해 정리하였다.

 

 

atan2를 이용한 직선의 표현

직선의 표현에 관한 직관적인 접근은 일차함수 $y=ax+b$, 즉 기울기와 절편으로 표현하는 것이다. 그러나 기울기가 무한대가 되는 예외 상황이 발생할 수 있어 좋지 않다. 따라서 두 점 p1, p2로 표현하는 것이 좀 더 적합하다.

 

정렬 등을 위해 기울기가 필요할 때는, 대신 두 점이 이루는 벡터의 부호 있는 각도를 사용할 수 있다. 벡터의 각도는 atan2 함수로 구한다. tan의 역함수로 기울기 l를 $(-\pi/2, \pi/2)$ 범위의 각도로 변환해 주는 atan(l)와 달리, atan2(y, x)는 y, x 각각을 입력받아 벡터 (x, y)가 양의 x축과 이루는 각도 를 범위 $(-\pi, \pi]$ 내에서 구한다. 예를 들어, $atan2(2, -2) = \frac{3}{4}\pi$, $atan2(-2, 2) = -\frac{1}{4}\pi$이다. 따라서 다음과 같이 직선에 관한 클래스를 구현할 수 있다. (Point 역시 클래스로 구현할 수도 있지만, 간단하게 (x,y)의 튜플 혹은 리스트를 가정했다)

 

class Line:
    def __init__(self, p1, p2):
        self.p1 = p1
        self.p2 = p2
        self.angle = math.atan2(p2[1]-p1[1], p2[0]-p1[0])

 

 

직선의 교점 구하기

직선의 교점을 구하는 것은 직선(선분)의 교차 판정만큼이나 생각보다 까다롭다. 먼저 p1,p2와 p3, p4를 지나는 각 직선의 일반형을 생각하자.

 

$$ (y_2 - y_1 )x - (x_2 - x_1 ) y = x_1 y_2 - x_2 y_1 $$

$$ (y_4 - y_3 )x - (x_4 - x_3 ) y = x_3 y_4 - x_4 y_3 $$

 

만약 두 직선이 평행하다면(det=0), 교점은 없거나 무수히 많다. 그렇지 않은 경우, Cramer's Rule을 활용해 x와 y를 구할 수 있다. 연립방정식의 상수항을 각각 C1, C2라 놓으면 다음과 같은 결과를 얻는다.

 

$$ x = \frac{det(A_1 ) }{det} = \frac{C_1 (x_3 - x_4 ) - C_2 (x_1 - x_2) }{det} $$

$$ y = \frac{det(A_2 ) }{det} = \frac{C_1 (y_3 - y_4 ) - C_2 (y_1 - y_2) }{det} $$

$$ det = (x_1 - x_2 ) (y_3 - y_4) - (y_1 - y_2 ) (x_3 - x_4 ) $$

 

점들의 좌표가 실수 형태로 주어진다면, 소수점 오차를 고려하여 코드를 아래와 같이 작성할 수 있다.

EPS = 1e-9
def get_intersection(l1, l2):
    p1, p2 = l1.p1, l1.p2
    p3, p4 = l2.p1, l2.p2
    det = (p1[0]-p2[0])*(p3[1]-p4[1]) - (p1[1]-p2[1])*(p3[0]-p4[0])
    if abs(det) < EPS: return None

    c1 = p1[0]*p2[1] - p1[1]*p2[0]
    c2 = p3[0]*p4[1] - p3[1]*p4[0]
    x = (c1*(p3[0]-p4[0]) - c2*(p1[0]-p2[0])) / det
    y = (c1*(p3[1]-p4[1]) - c2*(p1[1]-p2[1])) / det
    return [x, y]

 

 

반평면 교집합

반평면 교집합(Half Plane Intersection)은 여러 반평면들의 교집합 부분을 찾는 알고리즘이다. 편의상 반평면은 경계에 해당하는 직선(p1p2의 벡터로 표현)의 왼쪽 편으로 가정하자. 볼록 껍질 알고리즘을 점이 아닌 직선 단위로 수행하여, 큐에 필요한 직선만을 남겨 각각의 교점을 구하면 답을 얻을 수 있다. 시간복잡도는 볼록 껍질과 마찬가지로 정렬로 인해 $O(n \log n)$이다. 

 

1. 직선을 기울기(각도) 순으로 정렬해 차례로 큐에 넣는다.
2. 큐에 새로 넣을 직선을 기준으로, 직전 두 직선의 교점이 직선의 왼쪽 편, 즉 반평면에 속하는지 판정한다. 반평면에 속하지 않는다면 현재 큐에 들어 있는 마지막 직선은 필요 없으므로 pop한다.
3. 볼록 껍질과 차이점은, 첫 직선(들)이 반평면 교집합에 속한다는 보장은 없다는 것이다. 따라서 q[0], q[1]의 교점 역시 반평면에 속하는지 검사해야 하고, 속하지 않는다면 popleft한다. 남은 모든 직선들에 대해 2, 3을 반복한다.
4. append해준 마지막 직선(들) 역시 반평면 교집합에 속한다는 보장이 없다. 이번에는 q[0]를 기준으로, q[-1], q[-2]의 교점을 검사해서 반평면에 속하지 않는 경우 pop한다.

 

def not_in_hp(l, p):
    if p is None: return False
    return ccw(l.p1, l.p2, p) < 0

def hpi(lines):
    lines.sort(key=lambda l:l.angle)

    q = deque()
    for l in lines:
        while len(q) >= 2 and not_in_hp(l, get_intersection(q[-1], q[-2])):
            q.pop()
        while len(q) >= 2 and not_in_hp(l, get_intersection(q[0], q[1])):
            q.popleft()
        q.append(l)
    
    while len(q) >= 3 and not_in_hp(q[0], get_intersection(q[-1], q[-2])):
        q.pop()
    
    if len(q) < 3: return []
    res = []
    for i in range(len(q)):
        inter = get_intersection(q[i], q[(i+1)%len(q)])
        if inter is not None:
            res.append(inter)
    return res