ps:이론:원
원
삼각형의 외접원 (Circumscribed circle)
- Circumscribed circle
- 삼각형의 세 꼭짓점을 지나는 원은 유일하게 정의된다.
- 이 원을 찾는 것은, 원의 중심에서 세 꼭짓점까지의 거리가 모두 동일하다는 것을 가지고 연립방정식을 만들어서 풀면 된다.
- 식 정리가 복잡하긴 한데, 정리된 결과는 위키 링크를 참고하자.
- 처음에 꼭짓점 하나를 원점으로 평행이동시킨 다음에 중심점을 계산하고, 다시 중심점을 원래 좌표로 평행이동시키는 것이 계산이 쉽다.
- 이 식을 그대로 구현한 코드는 다음과 같다
def circumcircle_of_triangle(p, q, r): px, py = p bx, by = q[0] - px, q[1] - py cx, cy = r[0] - px, r[1] - py B = bx * bx + by * by C = cx * cx + cy * cy D = 2 * (bx * cy - by * cx) x, y = (cy * B - by * C) / D, (bx * C - cx * B) / D return px + x, py + y, math.hypot(x, y)
최소 외접원 (Smallest Enclosing Circle)
- 참고
- Smallest-circle problem
- 보통 최소 외접원이라고 번역하는데, 이게 정확한 번역인지는 의문이다. 외접원이면 모든 점과 접하는 원이어야 할텐데, 여기서 구하려는 것은 모든 점을 포함하는 원이기 때문이다. '모든 점을 포함하는 가장 작은 원'이 정확한 번역이긴 한데, 너무 길다보니 그냥 최소 외접원이라고 부르는것 같다. 마음에 드는 것은 아니지만 대안이 없으니 여기서도 그냥 그렇게 부르겠다.
- 위에서 설명한대로, 2차원 평면상에 주어진 N개의 점을 모두 포함하는 가장 작은 원을 찾는 문제이다
- 가장 먼 점과의 거리가 최소화되는 위치를 찾는 minimax problem 도 동일한 문제이다. 최소 외접원의 중심이 그 위치에 해당된다.
알고리즘
- 간단한 방법은 Gradient descent방식으로, 적당한 후보에서 출발해서 수렴할때까지 이터레이션을 도는 방식이다. 오차함수가 컨벡스하기 때문에, 글로벌 최적해를 찾을 수 있다.
- 좀더 PS적인 방식들은 위키 링크에 언급이 되어있다.
- Megiddo's algorithm 은 deterministic 하게 O(n)의 시간 복잡도를 갖는 알고리즘이다
- 하지만 그냥 위키에서 요약한 설명만 봐도 매우 지저분해보인다. 이터레이션을 돌때마다 후보 점의 갯수을 15/16 배로 줄여나간다고 한다;
- 자세한 설명은 https://greimul.tistory.com/32 에 있다.
- Welzl's algorithm 이 흔히 쓰이는 알고리즘이다. 랜덤에 기반한 알 고리즘이긴 하지만, 평균적으로 O(n)의 성능에 돌아간다.
- 삼소멤 링크 에 매우 친절하게 설명이 되어있다, 다만 여기에서 설명하는 알고리즘은 원래의 recursive한 Welzl's algorithm을, 조금 변형한 버전이다.
- 우리도 이 알고리즘을 사용하겠다
- Welzl's algorithm 의 원본 형태의 구현이 필요하면 https://www.geeksforgeeks.org/minimum-enclosing-circle-using-welzls-algorithm/ 을 참고.
구현
- 설명되어있는 대로 구현하면 되긴 하는데. 좀 헷갈릴 여지가 있기는 하니 주의..
- 대강 이런 식이 된다
def _is_in_circle(circle, p): *center, radius = circle dist = math.dist(center, p) return dist <= radius * MULTIPLICATIVE_EPS def _circle_from_pq(points, p, q): circle = ((p[0] + q[0]) * 0.5, (p[1] + q[1]) * 0.5, math.dist(p, q) * 0.5) cw_circle, ccw_circle = None, None ccw_area, cw_area = -INF, INF for r in points: if _is_in_circle(circle, r): continue ori = orientation(p, q, r) if ori == Orientation.COLLINEAR: continue *center, radi = circle = circumcircle_of_triangle(p, q, r) area = signed_parallelogram_area(p, q, center) if ori == Orientation.CCW and area > ccw_area: ccw_area, ccw_circle = area, circle elif ori == Orientation.CW and area < cw_area: cw_area, cw_circle = area, circle if cw_circle is None and ccw_circle is None: return circle elif cw_circle is None: return ccw_circle elif ccw_circle is None: return cw_circle else: return cw_circle if cw_circle[2] < ccw_circle[2] else ccw_circle def smallest_enclosing_circle(orig_points: Iterable[PointF]) -> CircleF: points = list(orig_points) if len(points) == 1: return (*points[0], 0) if len(points) == 2: p, q = points return (p[0] + q[0]) * 0.5, (p[1] + q[1]) * 0.5, math.dist(p, q) * 0.5 random.shuffle(points) circle = (0, 0, 0) for i, p in enumerate(points): if not _is_in_circle(circle, p): circle = (0, 0, 0) for j, q in enumerate(points[:i]): if not _is_in_circle(circle, q): circle = _circle_from_pq(points[:j], p, q) return circle- 코드를 실행시키면, 실수값을 계속 사용하게 된다.
- 실제 teflib에 구현된 코드는 위의 코드를 살짝 더 최적화한 버전이다
관련 문제
최대 내접원
- 반평면 교집합을 이용해서 해결 가능하다. 설명은 링크 참고.
ps/이론/원.txt · 마지막으로 수정됨: 2023/04/28 16:02 저자 teferi
토론