26_geometry.cpp

Download
cpp 353 lines 11.1 KB
  1/*
  2 * 기하 알고리즘 (Computational Geometry)
  3 * CCW, Convex Hull, Line Intersection, Point in Polygon
  4 *
  5 * 점, 선, 다각형 등 기하학적 문제를 해결합니다.
  6 */
  7
  8#include <iostream>
  9#include <vector>
 10#include <cmath>
 11#include <algorithm>
 12#include <stack>
 13
 14using namespace std;
 15
 16const double EPS = 1e-9;
 17const double PI = acos(-1.0);
 18
 19// =============================================================================
 20// 1. 점과 벡터
 21// =============================================================================
 22
 23struct Point {
 24    double x, y;
 25
 26    Point(double x = 0, double y = 0) : x(x), y(y) {}
 27
 28    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
 29    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
 30    Point operator*(double t) const { return Point(x * t, y * t); }
 31    Point operator/(double t) const { return Point(x / t, y / t); }
 32
 33    bool operator<(const Point& p) const {
 34        if (abs(x - p.x) > EPS) return x < p.x;
 35        return y < p.y;
 36    }
 37
 38    bool operator==(const Point& p) const {
 39        return abs(x - p.x) < EPS && abs(y - p.y) < EPS;
 40    }
 41
 42    double norm() const { return sqrt(x * x + y * y); }
 43    double norm2() const { return x * x + y * y; }
 44};
 45
 46// 내적
 47double dot(const Point& a, const Point& b) {
 48    return a.x * b.x + a.y * b.y;
 49}
 50
 51// 외적
 52double cross(const Point& a, const Point& b) {
 53    return a.x * b.y - a.y * b.x;
 54}
 55
 56// 두 점 사이 거리
 57double dist(const Point& a, const Point& b) {
 58    return (a - b).norm();
 59}
 60
 61// =============================================================================
 62// 2. CCW (Counter-Clockwise)
 63// =============================================================================
 64
 65// 양수: 반시계, 음수: 시계, 0: 일직선
 66double ccw(const Point& a, const Point& b, const Point& c) {
 67    return cross(b - a, c - a);
 68}
 69
 70int ccwSign(const Point& a, const Point& b, const Point& c) {
 71    double v = ccw(a, b, c);
 72    if (v > EPS) return 1;   // 반시계
 73    if (v < -EPS) return -1; // 시계
 74    return 0;                 // 일직선
 75}
 76
 77// =============================================================================
 78// 3. 볼록 껍질 (Convex Hull)
 79// =============================================================================
 80
 81vector<Point> convexHull(vector<Point> points) {
 82    int n = points.size();
 83    if (n < 3) return points;
 84
 85    sort(points.begin(), points.end());
 86
 87    vector<Point> hull;
 88
 89    // 하단 껍질
 90    for (int i = 0; i < n; i++) {
 91        while (hull.size() >= 2 &&
 92               ccw(hull[hull.size()-2], hull[hull.size()-1], points[i]) <= 0) {
 93            hull.pop_back();
 94        }
 95        hull.push_back(points[i]);
 96    }
 97
 98    // 상단 껍질
 99    int lower = hull.size();
100    for (int i = n - 2; i >= 0; i--) {
101        while (hull.size() > lower &&
102               ccw(hull[hull.size()-2], hull[hull.size()-1], points[i]) <= 0) {
103            hull.pop_back();
104        }
105        hull.push_back(points[i]);
106    }
107
108    hull.pop_back();  // 마지막 점(시작점과 동일) 제거
109    return hull;
110}
111
112// =============================================================================
113// 4. 선분 교차 판정
114// =============================================================================
115
116bool onSegment(const Point& p, const Point& q, const Point& r) {
117    return min(p.x, r.x) <= q.x && q.x <= max(p.x, r.x) &&
118           min(p.y, r.y) <= q.y && q.y <= max(p.y, r.y);
119}
120
121bool segmentIntersect(const Point& a, const Point& b,
122                      const Point& c, const Point& d) {
123    int d1 = ccwSign(c, d, a);
124    int d2 = ccwSign(c, d, b);
125    int d3 = ccwSign(a, b, c);
126    int d4 = ccwSign(a, b, d);
127
128    if (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) &&
129        ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0))) {
130        return true;
131    }
132
133    if (d1 == 0 && onSegment(c, a, d)) return true;
134    if (d2 == 0 && onSegment(c, b, d)) return true;
135    if (d3 == 0 && onSegment(a, c, b)) return true;
136    if (d4 == 0 && onSegment(a, d, b)) return true;
137
138    return false;
139}
140
141// =============================================================================
142// 5. 직선 교차점
143// =============================================================================
144
145// ax + by + c = 0 형태의 직선
146struct Line {
147    double a, b, c;
148
149    Line(const Point& p1, const Point& p2) {
150        a = p2.y - p1.y;
151        b = p1.x - p2.x;
152        c = -a * p1.x - b * p1.y;
153    }
154};
155
156bool lineIntersection(const Line& l1, const Line& l2, Point& intersection) {
157    double det = l1.a * l2.b - l2.a * l1.b;
158    if (abs(det) < EPS) return false;  // 평행
159
160    intersection.x = (l1.b * l2.c - l2.b * l1.c) / det;
161    intersection.y = (l2.a * l1.c - l1.a * l2.c) / det;
162    return true;
163}
164
165// =============================================================================
166// 6. 점과 다각형
167// =============================================================================
168
169// 점이 다각형 내부에 있는지 (Ray Casting)
170bool pointInPolygon(const Point& p, const vector<Point>& polygon) {
171    int n = polygon.size();
172    int count = 0;
173
174    for (int i = 0; i < n; i++) {
175        Point a = polygon[i];
176        Point b = polygon[(i + 1) % n];
177
178        if ((a.y <= p.y && p.y < b.y) || (b.y <= p.y && p.y < a.y)) {
179            double x = a.x + (p.y - a.y) * (b.x - a.x) / (b.y - a.y);
180            if (p.x < x) count++;
181        }
182    }
183
184    return count % 2 == 1;
185}
186
187// =============================================================================
188// 7. 다각형 넓이
189// =============================================================================
190
191double polygonArea(const vector<Point>& polygon) {
192    int n = polygon.size();
193    double area = 0;
194
195    for (int i = 0; i < n; i++) {
196        area += cross(polygon[i], polygon[(i + 1) % n]);
197    }
198
199    return abs(area) / 2.0;
200}
201
202// =============================================================================
203// 8. 최근접 점 쌍
204// =============================================================================
205
206double closestPair(vector<Point>& points, int lo, int hi) {
207    if (hi - lo <= 3) {
208        double minDist = 1e18;
209        for (int i = lo; i < hi; i++) {
210            for (int j = i + 1; j < hi; j++) {
211                minDist = min(minDist, dist(points[i], points[j]));
212            }
213        }
214        sort(points.begin() + lo, points.begin() + hi,
215             [](const Point& a, const Point& b) { return a.y < b.y; });
216        return minDist;
217    }
218
219    int mid = (lo + hi) / 2;
220    double midX = points[mid].x;
221
222    double d = min(closestPair(points, lo, mid),
223                   closestPair(points, mid, hi));
224
225    // 병합 (y 좌표 기준)
226    vector<Point> merged(hi - lo);
227    merge(points.begin() + lo, points.begin() + mid,
228          points.begin() + mid, points.begin() + hi,
229          merged.begin(),
230          [](const Point& a, const Point& b) { return a.y < b.y; });
231    copy(merged.begin(), merged.end(), points.begin() + lo);
232
233    // 스트립 내 점들 확인
234    vector<Point> strip;
235    for (int i = lo; i < hi; i++) {
236        if (abs(points[i].x - midX) < d) {
237            strip.push_back(points[i]);
238        }
239    }
240
241    for (int i = 0; i < (int)strip.size(); i++) {
242        for (int j = i + 1; j < (int)strip.size() &&
243             strip[j].y - strip[i].y < d; j++) {
244            d = min(d, dist(strip[i], strip[j]));
245        }
246    }
247
248    return d;
249}
250
251double closestPairDistance(vector<Point> points) {
252    sort(points.begin(), points.end());
253    return closestPair(points, 0, points.size());
254}
255
256// =============================================================================
257// 9. 회전 캘리퍼스
258// =============================================================================
259
260double convexDiameter(vector<Point>& hull) {
261    int n = hull.size();
262    if (n == 1) return 0;
263    if (n == 2) return dist(hull[0], hull[1]);
264
265    int j = 1;
266    double maxDist = 0;
267
268    for (int i = 0; i < n; i++) {
269        while (true) {
270            int next = (j + 1) % n;
271            Point v1 = hull[(i + 1) % n] - hull[i];
272            Point v2 = hull[next] - hull[j];
273
274            if (cross(v1, v2) > 0) {
275                j = next;
276            } else {
277                break;
278            }
279        }
280        maxDist = max(maxDist, dist(hull[i], hull[j]));
281    }
282
283    return maxDist;
284}
285
286// =============================================================================
287// 테스트
288// =============================================================================
289
290int main() {
291    cout << "============================================================" << endl;
292    cout << "기하 알고리즘 예제" << endl;
293    cout << "============================================================" << endl;
294
295    // 1. CCW
296    cout << "\n[1] CCW (Counter-Clockwise)" << endl;
297    Point a(0, 0), b(1, 1), c(2, 0);
298    cout << "    A(0,0), B(1,1), C(2,0)" << endl;
299    int ccwResult = ccwSign(a, b, c);
300    cout << "    CCW: " << (ccwResult > 0 ? "반시계" : ccwResult < 0 ? "시계" : "일직선") << endl;
301
302    // 2. 볼록 껍질
303    cout << "\n[2] 볼록 껍질" << endl;
304    vector<Point> points = {{0, 0}, {1, 1}, {2, 2}, {0, 2}, {2, 0}, {1, 3}};
305    auto hull = convexHull(points);
306    cout << "    입력 점: (0,0), (1,1), (2,2), (0,2), (2,0), (1,3)" << endl;
307    cout << "    볼록 껍질: ";
308    for (auto& p : hull) {
309        cout << "(" << p.x << "," << p.y << ") ";
310    }
311    cout << endl;
312
313    // 3. 선분 교차
314    cout << "\n[3] 선분 교차" << endl;
315    Point p1(0, 0), p2(2, 2), p3(0, 2), p4(2, 0);
316    cout << "    선분1: (0,0)-(2,2), 선분2: (0,2)-(2,0)" << endl;
317    cout << "    교차: " << (segmentIntersect(p1, p2, p3, p4) ? "예" : "아니오") << endl;
318
319    // 4. 다각형 넓이
320    cout << "\n[4] 다각형 넓이" << endl;
321    vector<Point> polygon = {{0, 0}, {4, 0}, {4, 3}, {0, 3}};
322    cout << "    사각형 (0,0)-(4,0)-(4,3)-(0,3)" << endl;
323    cout << "    넓이: " << polygonArea(polygon) << endl;
324
325    // 5. 점의 포함 판정
326    cout << "\n[5] 점의 다각형 포함 판정" << endl;
327    Point inside(2, 1), outside(5, 5);
328    cout << "    사각형: (0,0)-(4,0)-(4,3)-(0,3)" << endl;
329    cout << "    (2,1) 포함: " << (pointInPolygon(inside, polygon) ? "예" : "아니오") << endl;
330    cout << "    (5,5) 포함: " << (pointInPolygon(outside, polygon) ? "예" : "아니오") << endl;
331
332    // 6. 최근접 점 쌍
333    cout << "\n[6] 최근접 점 쌍" << endl;
334    vector<Point> pts = {{0, 0}, {3, 4}, {1, 1}, {5, 5}, {2, 1}};
335    cout << "    점: (0,0), (3,4), (1,1), (5,5), (2,1)" << endl;
336    cout << "    최근접 거리: " << closestPairDistance(pts) << endl;
337
338    // 7. 복잡도 요약
339    cout << "\n[7] 복잡도 요약" << endl;
340    cout << "    | 알고리즘       | 시간복잡도    |" << endl;
341    cout << "    |----------------|---------------|" << endl;
342    cout << "    | CCW            | O(1)          |" << endl;
343    cout << "    | 볼록 껍질      | O(n log n)    |" << endl;
344    cout << "    | 선분 교차      | O(1)          |" << endl;
345    cout << "    | 다각형 넓이    | O(n)          |" << endl;
346    cout << "    | 최근접 점 쌍   | O(n log n)    |" << endl;
347    cout << "    | 회전 캘리퍼스  | O(n)          |" << endl;
348
349    cout << "\n============================================================" << endl;
350
351    return 0;
352}