기하 관련 프로그래밍을 하다 보면 평면 상의 한 점이 다각형 내의 점인가 아니면 밖의 점인가를 판단해야 할 필요성이 많이 생긴다. 이 문제에 대한 답을 생각해보기로 하자. $N$개의 꼭짓점 $\{ x_i, y_i| i=0,...., N-1\}$으로 구성이 된 다각형을 놓고 생각을 하도록 하자. 평면상의 주어진 점 $(x_p, y_p)$에서 오른쪽으로 출발하는 반직선을 쭉 그으면 다각형과 교차하는 점들이 생기게 된다. 교차하는 상황을 잘 헤아려 보면 해당점이 다각형의 내부에 들어 있는 경우에는 항상 홀수번의 교차점이 생기고, 외부에 있는 경우에는 짝수번(0번 포함)의 교차점이 생김을 알 수 있다. 따라서 주어진 점에서 시작하는 반직선과 다각형의 교차점 개수의 홀짝 여부를 판별하면 내부/외부점인지에 대한 문제를 풀 수 있다. 여기서 한 가지 주의할 사항은, 주어짐 점에서 그은 반직선이 다각형의 변(edge)을 포함하는 경우는 제외해야 하고, 꼭짓점이 반직선에 걸리는 경우는 꼭짓점을 공유하는 두 변과 만나지만 한 번으로 카운트해야 한다.
** 주의: 경계점(+꼭지점)에 대한 처리는 일관성이 있는 결과를 주지 않는다.
BOOL is_inside_polygon(POINT p, POINT polygon[], int N) {
int counter = 0;
POINT* q1 = &polygon[0];
for (int i = 1; i <= N; i++) {
POINT *q2 = &polygon[i % N];
if (p.y > min(q1->y, q2->y)) { //y-intersect.;
if (p.y <= max(q1->y, q2->y)) { //y-intersect(=single);
if (p.x <= max(q1->x, q2->x)) //x-intersect;
if (q1->y != q2->y) { //not parallel;
double xq = double(p.y - q1->y) * (q2->x - q1->x) / (q2->y - q1->y) + q1->x;
// x-intersect;
if (q1->x == q2->x || p.x <= xq)
counter++;
}
}
}
}
q1 = q2; // move to next edge;
}
return (counter % 2);
}
원점을 중심으로 분산이 1 인 정규분포를 갖는 평면상의 점을 발생시킬 필요가 있는 경우에 아래의 함수를 이용한다: Box-Muller method.
내장 함수 rand()를 이용하면 $[0,1]$ 구간에서 균일한 분포를 만들 수 있다. $[0,1]\times[0,1]$ 영역에서 균일한 분포를 갖는 두 개의 1차원 랜덤 변수 $x_1$, $x_2$을 이용해서 확률 변수 $y_1$, $y_2$의 2차원 Gaussian 분포를 주는 변환을 구하려면, 우선 확률 보존에 의해서
두 개의 점이 주어지는 경우 이들을 잇는 가장 단순한 곡선은 직선이다. 또 직선은 두 점을 잇는 가장 짧은 거리를 갖는 곡선이기도 하다. 그러면 $N$개의 점이 주어지는 경우는 이들 모두 지나는 곡선을 어떻게 찾을까? 구간별로 직선으로 연결을 시켜서 하나의 곡선을 구성할 수 있으나 부드럽게 이어지는 형태의 곡선은 아니다. 곡선이 smooth 함은 그 곡선의 곡률이 잘 정의된다는 의미이다(곡률은 곡선의 이차 미분에 의존한다). 주어진 점들을 연결하는 충분히 짧고 smooth 한 곡선을 찾는 문제는 곡률에 의한 에너지를 최소화시키는 해를 구하는 문제로 환원할 수 있다. 이 문제의 해는 잘 알려진 cubic spline이다. 그러나 cubic spline은 지나는 점(컨트롤 점)의 변동에 대해서 안정적이 아니다. 컨트롤 점 중에서 한 점만 변해도 곡선이 전역적(global)으로 변하는 특성을 나타낸다. 컨트롤 점의 변동이 있을 때 그 점 주위에서만 변화하는 국소 특성을 갖는 곡선을 찾기 위해서는 smoothness 조건을 완화해야 한다. 그렇지만 충분히 부드러운 곡선이어야 하므로 적어도 일차의 미분 값의 연속성 정도는 보장되어야 한다. 이런 특성을 갖는 삼차 곡선이 Catmull-Rom spline이다. 이름은 Edwin Catmull와 Raphie Rom의 이름에서 유래한다.
열린 곡선으로 보간을 하는 경우 처음과 끝 구간에는 미분 값을 구할 수 없으므로 정의가 안되지만, 나머지 구간에서는 주어진 점들을 연결하는 충분히 부드러운 곡선이 된다. 양 끝 구간의 보간은 동일점을 반복함으로써 전체적인 구간에서 잘 정의되게 만들 수 있다(끝점에서 기울기를 그 점과 인접점과의 기울기/2로 잡는 것과 같은 효과임). 닫힌 곡선인 경우에는 모든 구간에서 잘 정의가 된다. Catmull-Rom spline은 Bezier나 B-Spline 곡선의 경우와 다르게 컨트롤 점의 convex hull 내부에서만 정의되는 특성을 따르지 않는다.
CPoint CRSpline(double t, CPoint p1, CPoint p2, CPoint p3, CPoint p4) {
double tt = t * t ;
double ttt = tt * t ;
double x = 0.5 * ((-p1.x + 3 * p2.x - 3 * p3.x + p4.x) * ttt
+ (2 * p1.x - 5 * p2.x + 4 * p3.x - p4.x) * tt
+ (-p1.x + p3.x) * t
+ 2 * p2.x);
double y = 0.5 * ((-p1.y + 3 * p2.y - 3 * p3.y + p4.y) * ttt
+ (2 * p1.y - 5 * p2.y + 4 * p3.y - p4.y) * tt
+ (-p1.y + p3.y) * t
+ 2 * p2.y);
return CPoint(int(x + .5), int(y + .5)) ;
}
//open spline의 drawing;
void DrawCatmullRom(std::vector<CPoint>& Q, CDC* pDC) {
#define STEPS (20)
if (Q.size() < 4) return ;
CPen red(PS_SOLID, 1, RGB(0xFF, 0, 0));
CPen *pOld = pDC->SelectObject(&red);
const int n = Q.size();
for (int i = 0; i < n - 1; i++) {
pDC->MoveTo(Q[i]);
// i = 0 인 경우에는 처음 점을 반복, i = n - 2인 경우에는 마지막 점을 반복..
int ip = max(i - 1, 0);
int inn = min(i + 2, n - 1);
for (int it = 1; it <= STEPS; it++)
pDC->LineTo(CRSpline(double(it)/STEPS, Q[ip], Q[i], Q[i + 1], Q[inn]));
};
pDC->SelectObject(pOld);
}
#define REMOVE_FLAG (-1)
struct TRIANGLE {
int a, b, c;
TRIANGLE() {}
TRIANGLE(int va, int va, int vc) : a(va), b(vb), c(vc) {}
};
struct EDGE {
int a, b;
EDGE() {}
EDGE(int va, int vb) : a(va), b(vb) {}
void Unflag() {a = REMOVE_FLAG; b = REMOVE_FLAG;}
bool Unflagged() {a == REMOVE_FLAG || b == REMOVE_FLAG;}
bool IsSame(EDGE e) {
return ((a == e.a && b == e.b) || (a == e.b && b == e.a)) ;
}
};
// 수퍼 삼각형은 입력 점집합을 모두 포함하도록 충분히 큰 삼각형을 선택한다;
// 삼각화한 후 최외각이 convex가 아니면 수퍼 삼각형을 더 크게 다시 잡아야 함.
// 마지막 수정: 2021.4.7, 코드 간결화; 2024.4.5, 단순화;
// max num of triangles = 2 * (N - 2) - 1;
// max num of edegs = 3 * (N - 2);
int inc_delaunay_triangulation(std::vector<CPoint>& pts, std::vector<TRIANGLE>& V) {
const int N = pts.size()
if (N < 3) return 0;
int tmax = 2 * ((N + 3) - 2);
int emax = 3 * ((N + 3) - 2);
std::vector<EDGE> E(emax);
V.resize(tmax);
std::vector<CPoint> P(pts.size());
for (int i = points.size(); i-->0; ) P[i] = pts[i];
set_supertriangle(P); // P.size()은 +3 만큼 증가함;
// First element of triangle index is supertriangle vertex;
int nt = 0;
V[nt++] = TRIANGLE(N, N + 1, N + 2);
for (int i = 0; i < N; i++) {
int ne = 0;
/* 추가점(P[i])를 외접원에 포함시키는 삼각형을 모두 찾아 제거한다.
** 이들 삼각형의 공유되는 내부 edge를 제외한 나머지 외부 edge와
** 추가점으로 새로운 삼각형을 만들어 리스트에 추가한다. */
for (int j = 0; j < nt; j++) {
TRIANGLE &t = V[j];
if (inCircle(P[t.a], P[t.b], P[t.c], P[i])) {
// 제거되는 삼각형의 edge는 백업한다; 내부 edge는 중복된다.
E[ne++] = EDGE(t.a, t.b);
E[ne++] = EDGE(t.b, t.c);
E[ne++] = EDGE(t.c, t.a);
// P[i]를 포함하는 외접원 삼각형은 제거;
// swap(j, nt-1) to remove triangle at j
V[j--] = V[--nt];
}
}
// 중복되어 있는 내부 edge를 제거;
remove_double_edges(E, ne);
// 외부 edge와 P[i]가 만드는 새 삼각형을 리스트에 추가;
for (int j = 0; j < ne; j++) {
if (E[j].Unflagged()) continue;
V[nt++] = TRIANGLE(E[j].a, E[j].b, i);
}
}
// done with super_triangle;
V.resize(nt);
return remove_supertriangle_relics(N, V) ;
}
// double edge 제거;
void remove_double_edges(std::vector<EDGE>& E, int ne) {
for (int j = 0; j < ne - 1; j++)
for (int k = j + 1; k < ne; k++)
if (E[j].IsSame(E[k])) {
E[j].Unflag(); E[k].Unflag();
}
}
// 삼각형 중에서 수퍼 삼각형의 꼭지점(N,N+1,N+2)을 가지고 있는 것을 제거;
int remove_supertriangle_relics(int N, std::vector<TRIANGLE>& V) {
int nt = V.size();
for (int i = 0; i < nt; i++)
if (V[i].a >= N || V[i].b >= N || V[i].c >= N)
// swap(i, nt - 1) to remove triangle at j
V[i--] = V[--nt];
V.resize(nt);
return nt;
}
// A,B,C가 만드는 외접원에 D가 들어가는가?;
int inCircle(CPoint A, CPoint B, CPoint C, CPoint D) {
CPoint AD = A - D, BD = B - D, CD = C - D;
double ccw = CCW(A, B, C);
double det = norm2(AD) * cross(BD, CD)
+ norm2(BD) * cross(CD, AD)
+ norm2(CD) * cross(AD, BD)
if (ccw > 0) return det > 0; //CCW인 경우에..
else return det < 0; //CW인 경우에..
}
