Translated 출처: e-maxx.ru

두 선분의 교점 찾기

두 선분 AB와 CD가 각각 두 끝점의 쌍으로 주어집니다. 두 끝점이 같으면 선분이 하나의 점일 수도 있습니다. 이 두 선분의 교집합을 구해야 합니다. 교집합은 비어 있거나(선분이 교차하지 않는 경우), 하나의 점이거나, 선분(주어진 선분이 겹치는 경우)일 수 있습니다.

풀이

선분의 교점은 직선의 교점과 같은 방법으로 구할 수 있습니다. 선분의 끝점으로부터 직선의 방정식을 다시 세우고 두 직선이 평행한지 확인합니다.

두 직선이 평행하지 않다면 교점을 구한 뒤 그 교점이 두 선분 모두에 속하는지 확인해야 합니다(이를 위해서는 교점이 각 선분을 X축과 Y축에 투영한 구간에 속하는지 확인하는 것으로 충분합니다). 이 경우 답은 "교점 없음"이거나 두 직선의 교점 하나입니다.

두 직선이 평행한 경우는 조금 더 복잡합니다(하나 이상의 선분이 하나의 점인 경우도 여기에 해당합니다). 이 경우 두 선분이 같은 직선 위에 있는지 확인해야 합니다. 같은 직선 위에 있지 않다면 답은 "교점 없음"입니다. 같은 직선 위에 있다면 답은 같은 직선 위에 있는 두 선분의 교집합입니다. 이는 두 선분의 끝점을 특정 좌표의 오름차순으로 정렬한 뒤, 왼쪽 끝점 중 가장 오른쪽에 있는 점과 오른쪽 끝점 중 가장 왼쪽에 있는 점을 취하여 구합니다.

두 선분이 모두 하나의 점이라면 두 점은 같아야 하며, 이 확인은 따로 수행하는 것이 좋습니다.

알고리즘의 시작 부분에 바운딩 박스 검사를 추가합시다. 이 검사는 두 선분이 같은 직선 위에 있는 경우에 반드시 필요하며, (가벼운 검사이므로) 무작위 테스트에서 알고리즘이 평균적으로 더 빠르게 동작하도록 해 줍니다.

구현

다음은 직선과 선분을 처리하는 모든 보조 함수를 포함한 구현입니다.

메인 함수 intersect는 두 선분의 교집합이 비어 있지 않으면 true를 반환하고, 교집합 선분의 끝점을 인자 left와 right에 저장합니다. 답이 하나의 점이라면 left와 right에 기록되는 값은 같습니다.

const double EPS = 1E-9;

struct pt {
    double x, y;

    bool operator<(const pt& p) const
    {
        return x < p.x - EPS || (abs(x - p.x) < EPS && y < p.y - EPS);
    }
};

struct line {
    double a, b, c;

    line() {}
    line(pt p, pt q)
    {
        a = p.y - q.y;
        b = q.x - p.x;
        c = -a * p.x - b * p.y;
        norm();
    }

    void norm()
    {
        double z = sqrt(a * a + b * b);
        if (abs(z) > EPS)
            a /= z, b /= z, c /= z;
    }

    double dist(pt p) const { return a * p.x + b * p.y + c; }
};

double det(double a, double b, double c, double d)
{
    return a * d - b * c;
}

inline bool betw(double l, double r, double x)
{
    return min(l, r) <= x + EPS && x <= max(l, r) + EPS;
}

inline bool intersect_1d(double a, double b, double c, double d)
{
    if (a > b)
        swap(a, b);
    if (c > d)
        swap(c, d);
    return max(a, c) <= min(b, d) + EPS;
}

bool intersect(pt a, pt b, pt c, pt d, pt& left, pt& right)
{
    if (!intersect_1d(a.x, b.x, c.x, d.x) || !intersect_1d(a.y, b.y, c.y, d.y))
        return false;
    line m(a, b);
    line n(c, d);
    double zn = det(m.a, m.b, n.a, n.b);
    if (abs(zn) < EPS) {
        if (abs(m.dist(c)) > EPS || abs(n.dist(a)) > EPS)
            return false;
        if (b < a)
            swap(a, b);
        if (d < c)
            swap(c, d);
        left = max(a, c);
        right = min(b, d);
        return true;
    } else {
        left.x = right.x = -det(m.c, m.b, n.c, n.b) / zn;
        left.y = right.y = -det(m.a, m.c, n.a, n.c) / zn;
        return betw(a.x, b.x, left.x) && betw(a.y, b.y, left.y) &&
               betw(c.x, d.x, left.x) && betw(c.y, d.y, left.y);
    }
}