2차원 점집합의 median
'Computational Geometry' 카테고리의 다른 글
Bresenham's Line Algorithm (0) | 2021.04.07 |
---|---|
Rotating Calipers (3) | 2021.03.31 |
Graham Scan (0) | 2021.03.28 |
Jarvis March (0) | 2021.03.26 |
Approximate Minimum Enclosing Circle (1) | 2021.03.18 |
2차원 점집합의 median
Bresenham's Line Algorithm (0) | 2021.04.07 |
---|---|
Rotating Calipers (3) | 2021.03.31 |
Graham Scan (0) | 2021.03.28 |
Jarvis March (0) | 2021.03.26 |
Approximate Minimum Enclosing Circle (1) | 2021.03.18 |
Convex hull을 찾을 때, 중복점과 일직선상 위에 있는 점들에 대해서 적절한 처리를 하지 않으면 효율적으로 작성된 알고리즘이라도 정상적으로 동작하지 않을 수 있다. 그림은 Graham Scan 알고리즘을 이용해서 50,000개 점들의 convex hull을 구한 결과다. qsort()을 호출할 때 compare()에서 global 변수를 쓰지 않기 위해서 pivot point를 뺀 후 사용한다. 정렬 후에는 다시 더해준다.
inline int LeftSide(CPoint A, CPoint B, CPoint C) {
return (B.x - A.x) * (C.y - A.y) - (C.x - A.x) * (B.y - A.y); // cross(AB, AC);
}
static int cross(CPoint *p1, CPoint *p2) {
return p1->x * p2->y - p1->y * p2->x;
}
static int norm2(CPoint *p) {
return p->x * p->x + p->y * p->y;
}
static int compare(const void *a, const void *b) {
CPoint *p1 = (CPoint *)a, *p2 = (CPoint *)b;
int ccw = cross(p2, p1);
// 거리가 먼 곳 --> 가까운 곳;
if (ccw == 0) return norm2(p2) > norm2(p1) ? 1 : -1;
return ccw; // 반시계방향 정렬;
}
static void angularSort(std::vector<CPoint> & pts) {
if ((pts.size()) < 2) return;
int id = 0;
for (int i = pts.size() - 1; i >= 1; --i)
if ((pts[id].y > pts[i].y) || (pts[id].y == pts[i].y && pts[i].x > pts[id].x))
id = i; //lowest--y-value; //largest--x-value;
std::swap(pts[0], pts[id]);
for (int i = pts.size() - 1; i >= 1; --i) pts[i] -= pts[0];
qsort(&pts[1], pts.size() - 1, sizeof(CPoint), compare);
for (int i = pts.size() - 1; i >= 1; --i) pts[i] += pts[0];
};
// AB와 BC가 나란할 때, AB 방향과 BC방향이 같은가(0) 반대면 (1): dot(AB, BC);
static int folded(CPoint A, CPoint B, CPoint C) {
return (B.x - A.x) * (C.x - B.x) + (B.y - A.y) * (C.y - B.y) <= 0;
}
int GrahamScan(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
if (pts.size() < 3) return 0;
angularSort(pts);
hull.resize(pts.size());
hull[0] = pts[0];
int i = 1;
while (i < pts.size()) {
// peek out points degenerate with pts[0];
if (pts[i] != hull[0]) {
hull[1] = pts[i++]; break; // (i++)
}
i++;
}
if (i == pts.size()) return 0;
int k = 2; // 새로 추가되는 hull index;
while (i < pts.size()) {
// (i)가 hull의 left에 있으면(1) hull의 꼭지점에 추가;
int ccw = LeftSide(hull[k - 2], hull[k - 1], pts[i]);
if (ccw > 0) hull[k++] = pts[i++];
// (i)가 hull의 직전 변과 collinear한 데 반대방향을 향하면 skip;
else if (ccw == 0 && folded(hull[k - 2], hull[k - 1], pts[i])) ++i;
// (i)가 hull의 right 또는 위에 있으면(직전 변과 같은 방향)
// hull의 직전 꼭지점을 순차적으로 제거한 후 다시 검사;
else --k;
ASSERT(k >= 2);
}
hull.resize(k);
return hull.size();
};
Rotating Calipers (3) | 2021.03.31 |
---|---|
Convex Hull Peeling (0) | 2021.03.29 |
Jarvis March (0) | 2021.03.26 |
Approximate Minimum Enclosing Circle (1) | 2021.03.18 |
Minimum Bounding Rectangle (3) | 2021.03.15 |
기하 알고리즘에서 입력점의 개수가 많아지면 중복점이 생기고, 이 때문에 최적화에 주안점을 둔 알고리즘은 중복점에 대한 예외처리를 하지 않으면 정상적으로 동작하지 않는 경우가 종종 생긴다. 안정성과 구현의 간결함을 유지하고 싶으면 좀 더 brute force method에 가까운 알고리즘을 사용하는 것도 고려해 볼 수 있다. 다음은 Jarvis March 알고리즘을 써서 2차원 convex hull을 구하는 코드이다. (그림은 50,000개의 입력점을 사용함.)
BOOL ccw(CPoint a, CPoint b, CPoint c){
// cross(ac, ab) < 0: ab가 ac보다 오른쪽에 있음 = b가 ac변의 오른쪽에 있음.
return ((c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x)) < 0;
}
int JarvisMarch(std::vector<CPoint> &pts, std::vector<CPoint> &hull) {
int n = pts.size();
if (n < 3) return 0;
int left = 0;
for (int i = 1; i < n; i++) //convex hull의 첫번째는 맨 왼쪽 점;
if (pts[i].x < pts[left].x) left = i;
hull.clear();
int turns = 0; //debug 용도
int prev = left;
while (1) {
hull.push_back(pts[prev]);
int next = (prev + 1) % n;
// prev-next 변이 모든 점의 오른쪽에 있게 하는 next를 찾으면 됨;
// 없으면 next가 convex hull 위의 점임;
for (int i = 0; i < n; i++) {
if (i == next || i == prev) continue;
if (ccw(pts[prev], pts[i], pts[next]))
next = i;
}
prev = next;
turns++; //debug 용도;
if (prev == left) break;
}
TRACE("number of turns = %d\n", turns);
remove_collinear(hull); // removes collinear vertices(optional);
return hull.size();
}
int remove_collinear(std::vector<CPoint>& V) {
int N = V.size();
if (N < 3) return 0;
// 먼저, 일직선상에 있지 않는 vertex을 찾음;
int start;
for (start = 0; start < N; start++)
if (ccw(V[(start - 1 + N) % N], V[start], V[(start + 1) % N]))
break;
std::vector<CPoint> H;
H.push_back(V[start]);
int prev = start;
int curr = (prev + 1) % N;
int next = (curr + 1) % N;
while (curr != start) {
if (ccw(V[prev], V[curr], V[next])) {
H.push_back(V[curr]);
prev = curr;
}
curr = (curr + 1) % N;
next = (curr + 1) % N;
}
std::swap(V, H);
return V.size();
};
Convex Hull Peeling (0) | 2021.03.29 |
---|---|
Graham Scan (0) | 2021.03.28 |
Approximate Minimum Enclosing Circle (1) | 2021.03.18 |
Minimum Bounding Rectangle (3) | 2021.03.15 |
Minimum Enclosing Circle (0) | 2021.03.01 |
주어진 점집합을 포함하는 최소 면적의 rectangle을 찾는 문제는 점집합의 convex hull을 포함하는 rectangle을 찾는 문제와 같다. convex hull의 한 변을 공유하면서 전체를 포함하는 rectangle 중 최소 면적인 것을 찾으면 된다. Rotating Caliper를 쓰면 O(n) step이 요구된다.
// 방향이 dir이고 A점을 통과하는 직선과 방향이 dir에 수직이고 B점을 통과하는 직선의 교점을 반환;
CPoint LineIntersect(CPoint A, CPoint dir, CPoint B);
//
#define DOT(A,B) ((A).x*(B).x + (A).y*(B).y))
void BoundingRect(std::vector<CPoint>& hull, CPoint vtx[4]) {//non-optimal.
if (hull.size() < 3) return;
double minArea = DBL_MAX;
int bIndex = 0, tIndex = 0, lIndex = 0, rIndex = 0;
for (int j = 0, i = hull.size() - 1; j < hull.size(); i = j++) {
CPoint &A = hull[i], &B = hull[j];
CPoint BA = B - A; // convex hull의 한 기준변;
CPoint BAnormal = CPoint(-BA.y, BA.x); //(B-A)에 반시계방향으로 수직인 벡터;
double BAsq = DOT(BA, BA);
//기준변과 대척점(antipodal point)을 찾는다: BAnormal방향으로 가장 멀리 떨어진 점
int id = FindAntipodal(BAnormal, hull) ;
CPoint QA = hull[id] - hull[i];
//기준변과 대척점을 통과하는 직선과의 사이거리;
double D1 = fabs(double(DOT(BAnormal, QA))) / sqrt(BAsq);
//left_side_end_point;//기준변 반대방향으로 가장 멀리 떨어진 꼭지점;
int id1 = FindAntipodal(-BA, hull);
//right_side_end_point;//기준변 방향으로 가장 멀리 떨어진 꼭지점;
int id2 = FindAntipodal(BA, hull);
///가장 왼쪽과 가장 오른쪽 꼭지점을 연결하는 선분
CPoint H = hull[id1] - hull[id2];
//기준변 방향으로 정사영하면 두 평행선 사이거리를 줌;
double D2 = fabs(double(DOT(BA,H))) / sqrt(BAsq);
double area = D1 * D2;
if (area < minArea) {
minArea = area;
bIndex = i ; tIndex = id;
lIndex = id1; rIndex = id2;
}
}
//directional vector for base_edge;
CPoint dir = hull[(bIndex + 1) % hull.size()] - hull[bIndex] ;
vtx[0] = LineIntersect(hull[bIndex], dir, hull[lIndex]);
vtx[1] = LineIntersect(hull[tIndex], dir, hull[lIndex]);
vtx[2] = LineIntersect(hull[tIndex], dir, hull[rIndex]);
vtx[3] = LineIntersect(hull[bIndex], dir, hull[rIndex]);
}
Jarvis March (0) | 2021.03.26 |
---|---|
Approximate Minimum Enclosing Circle (1) | 2021.03.18 |
Minimum Enclosing Circle (0) | 2021.03.01 |
Creating Simple Polygons (0) | 2021.01.25 |
단순 다각형의 Convex hull (0) | 2021.01.24 |
평면 상에 주어진 점들의 convex hull을 구하는 알고리즘 중의 하나인 Graham scan에서는 먼저 주어진 점들을 한 점을 기준으로 각도로 정렬하는 과정이 필요했다. 그러면 점들이 순차적으로 연결된 단순 다각형에서는 sorting과정이 없이 Graham scan 알고리즘을 바로 적용이 가능한 것일까? 이에 대한 counter example은 쉽게 찾을 수 있다. 단순 폴리곤에 대해서도 항상 각도에 대한 정렬 과정이 필요한 것인가? 답은 아니다. 정렬 과정이 없이도 단순 폴리곤에 대해서는 쉽게 convex hull을 구할 수 있는 알고리즘이 존재한다. 정렬 과정이 없이 단순 폴리곤의 convex hull을 찾는 알고리즘에 대해서 알아보자.
Melkman Algorithm;
우선 폴리곤의 이웃하는 세 꼭짓점을 잡아서, 반시계 방향의 삼각형을 구성한다. 이 삼각형을 deque에 넣는다 (bottom = top). 폴리곤을 순환하면서 새 점이 들어올 때 이미 만들어진 convex hull의 내부점이 아니면 이 점이 포함되도록 convex hull을 업데이트한다: Graham scan 알고리즘의 scanning 과정을 bottom을 기준으로 반시계 방향으로 convexity를 만족시킬 때까지 bottom을 제거하고, top을 기준으로는 시계방향으로 convexity를 만족시킬 때까지 top을 제거한다. 이 과정이 끝나면 새 점을 deque에 추가한다. 이 과정을 나머지 모든 점들에 대해서도 진행한다.
int LEFTSIDE(CPoint A, CPoint B, CPoint C){ // cross(AB, AC) > 0: C가 AB 대해서 왼쪽
return ((B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x)) > 0;
}
int convexhull_spolygon(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
int N = pts.size();
if (N < 3) return 0;
std::vector<CPoint> Q(2 * N + 1); // deque;
int bot = N - 2, top = bot + 3;
// |0|1|......|N-2|N-1|N+0|N+1|N+2|..............|2N|;
// bot top
Q[bot] = Q[top] = pts[2];
if (LEFTSIDE(pts[0], pts[1], pts[2])) { //2->0->1->2; Q에 ccw로 입력;
Q[bot + 1] = pts[0]; Q[bot + 2] = pts[1];
} else {
Q[bot + 1] = pts[1]; Q[bot + 2] = pts[0];
}
int i = 2; // last touched index;
while (++i < N) {
// 기존의 convex_hull에 들어있으면 제외.
if (LEFTSIDE(Q[bot + 0], Q[bot + 1], pts[i]) &&
LEFTSIDE(Q[top - 1], Q[top + 0], pts[i]))
continue;
// bottom에 대해서 ccw 방향으로 체크(bot 증가 방향)
// pts[i]가 (bot)->(bot+1)라인의 오른쪽에 있으면, 기존의 bottom을 제외;
while (!LEFTSIDE(Q[bot], Q[bot + 1], pts[i]) && (bot < top)) bot++;
Q[--bot] = pts[i];
// 추가점에서 top->top-1의 ccw 방향으로 convexity 체크하여 만족하지 않은 경우
// top을 감소
while (!LEFTSIDE(Q[top - 1], Q[top], pts[i]) && (top > bot)) top-- ;
Q[++top] = pts[i];
}
hull.resize(top - bot);
for (int i = hull.size() - 1; i >= 0; --i) hull[i] = Q[i + bot];
return hull.size();
};
**네이버 블로그 이전;
Minimum Enclosing Circle (0) | 2021.03.01 |
---|---|
Creating Simple Polygons (0) | 2021.01.25 |
단순 다각형의 무게중심 (0) | 2021.01.24 |
단순 다각형의 면적(2D) (0) | 2021.01.23 |
삼각형 외접원의 Inclusion Test (0) | 2020.12.30 |