최근접 점의 쌍 찾기

#include <stdio.h>
#include <math.h>

# define SWAP(x, y, temp) ( (temp)=(x), (x)=(y), (y)=(temp) )

void merging(double* X, unsigned* Xid, unsigned* TMP, unsigned left, 
unsigned mid, unsigned right) {
int i, j, k, l;
i = left;
j = mid + 1;
k = left;

/* 분할 정렬된 list의 합병 */
while (i <= mid && j <= right) {
if (X[TMP[i]] <= X[TMP[j]])
Xid[k++] = TMP[i++];
else
Xid[k++] = TMP[j++];
}

// 남아 있는 값들을 일괄 복사
if (i > mid) {
for (l = j; l <= right; l++)
Xid[k++] = TMP[l];
}
// 남아 있는 값들을 일괄 복사
else {
for (l = i; l <= mid; l++)
Xid[k++] = TMP[l];
}

// 배열 sorted[](임시 배열)의 리스트를 배열 list[]로 재복사
for (l = left; l <= right; l++) {
TMP[l] = Xid[l];
}
}

void merge_sort(double* X, unsigned* Xid, unsigned* TMP, unsigned left, unsigned right) {
int mid;

if (left < right) {
mid = (left + right) / 2; // 중간 위치를 계산하여 리스트를 균등 분할 -분할(Divide)
merge_sort(X, Xid, TMP, left, mid); // 앞쪽 부분 리스트 정렬 -정복(Conquer)
merge_sort(X, Xid, TMP, mid + 1, right); // 뒤쪽 부분 리스트 정렬 -정복(Conquer)
merging(X, Xid, TMP, left, mid, right); // 정렬된 2개의 부분 배열을 합병하는 과정 -결합(Combine)
}
}

//x좌표에 비내림차순으로 정렬해서 그 인덱스를 xid에 저장
void   sortXid(double* X, unsigned* Xid, unsigned* TMP, unsigned N)
// input  : X[]   (x position array of N points)
//          Xid[] (index array of X (initially [0,1,...,N-1])
//          TMP[] (temporary array of size N. may or may not be used)
//          N   number of points to sort
// output : Xid[] (sorted indices by x position in X[])
{
for (int i = 0; i < N; i++) {
TMP[i] = i;
}
merge_sort(X, Xid, TMP, 0, N - 1);

for (int i = 0; i < N; i++) {
Xid[i] = TMP[i];
}
for (int i = 0; i < N; i++) {
TMP[i] = 0;
}
}

//y위치의 오름차순 정렬. 
void selection_sort(double * Y, unsigned* Yid, unsigned L, unsigned R, int n) {
int i, j, least, temp;

// 마지막 숫자는 자동으로 정렬되기 때문에 (숫자 개수-1) 만큼 반복한다.
for (i = L; i < R; i++) {
least = i;

// 최솟값을 탐색한다.
for (j = i + 1; j <= R; j++) {
if (Y[Yid[j]] < Y[Yid[least]])
least = j;
}

// 최솟값이 자기 자신이면 자료 이동을 하지 않는다.
if (i != least) {
SWAP(Yid[i], Yid[least], temp);
}
}
}

void merge(unsigned * Yid, unsigned* TMP, int left, int mid, int right, double *Y) {
int i, j, k, l;
i = left;
j = mid + 1;
k = left;

/* 분할 정렬된 list의 합병 */
while (i <= mid && j <= right) {
if (Y[Yid[i]] <= Y[Yid[j]])
TMP[k++] = Yid[i++];
else
TMP[k++] = Yid[j++];
}

// 남아 있는 값들을 일괄 복사
if (i > mid) {
for (l = j; l <= right; l++)
TMP[k++] = Yid[l];
}
// 남아 있는 값들을 일괄 복사
else {
for (l = i; l <= mid; l++)
TMP[k++] = Yid[l];
}

// 배열 sorted[](임시 배열)의 리스트를 배열 list[]로 재복사
for (l = left; l <= right; l++) {
Yid[l] = TMP[l];
}
}

double  combine(unsigned* Yid, unsigned L, unsigned R, unsigned M,
double dlr, double* X, unsigned* Xid, unsigned* TMP, double* Y,
unsigned* point1, unsigned* point2) {
double xmid = (X[Xid[M]] + X[Xid[M + 1]]) / 2;
/*yid의 l부터 r까지 스캔, 만일 x좌표가 dlr 내에 속하면 tmp로 카피
* xmid기준으로, xid의 xmid-dlr, xmid+dlr까지 내에 속할 경우 tmp로 카피
*/
int k = 0;
double cpc = 1.797679e+308;
double dist;
for (int i = L; i <= R; i++) {
if(X[Yid[i]] >= xmid-dlr && X[Yid[i]] <= xmid + dlr) {
//y방향 정렬 순서가 유지됨
TMP[k] = Yid[i];

//tmp배열에 2개이상 삽입시, 비교시작. 만일 하나만 충족 시, 하나만 리턴할 듯
if (k > 0) {
dist= sqrt((X[TMP[k]] - X[TMP[k-1]]) * (X[TMP[k]] - X[TMP[k - 1]]) 
+ (Y[TMP[k]] - Y[TMP[k-1]]) * (Y[TMP[k]] - Y[TMP[k - 1]]));
if (dist < cpc) {
cpc = dist;
*point1 = (unsigned)TMP[k];
*point2 = (unsigned)TMP[k-1];
}
}
k++;
}
}

return cpc;
}

double closestPairDC(
unsigned L, unsigned R,   // current leftmost and rightmost indices
unsigned* pt1, unsigned* pt2, // closest pair points indices
double* X, double* Y,         // (x,y) positions array(input)
unsigned* Xid,  // point index array(sorted by x coordinates)
unsigned* Yid,  // point index array(gradually sorted by y coordinates)
unsigned* TMP,  // temporal index array
unsigned THR // threshold value
) {
double dl, dr, dmin;
unsigned M;
double dlr;
unsigned point1 = 0;
unsigned point2 = 0;
double d = 1.797679e+308;

for (int i = 0; i <= THR; i++) {
TMP[i] = 0;
}

if (R - L + 1 <= THR) {
for (int i = L; i <= R; i++) {
Yid[i] = Xid[i];
}

//y위치의 오름차순 정렬
selection_sort(Y, Yid, L, R, R - L + 1);

double dij;
for (int i = L; i <= R-1; i++) {
for (int j = i + 1; j <= R; j++) {
// if (Y[Yid[j]] - Y[Yid[i]] >= d) 
// break;

dij= sqrt((X[Yid[j]] - X[Yid[i]]) * (X[Yid[j]] - X[Yid[i]])
+ (Y[Yid[j]] - Y[Yid[i]]) * (Y[Yid[j]] - Y[Yid[i]]));

if (dij < d) {
d = dij;
*pt1 = (unsigned)Yid[i];
*pt2 = (unsigned)Yid[j];
}
}
}
return d;
//ok
}
else {

M = (L + R) / 2;

//지점 바꾸기용
unsigned ptr1 = 0; unsigned ptr2 = 0;
unsigned ptr3 = 0; unsigned ptr4 = 0;

dl = closestPairDC(L, M, pt1, pt2, X, Y, Xid, Yid, TMP, THR);
ptr1 = *pt1;
ptr2 = *pt2;
dr = closestPairDC(M+1, R, pt1, pt2, X, Y, Xid, Yid, TMP, THR);
ptr3 = *pt1;
ptr4 = *pt2;

//merge
merge(Yid, TMP, L, M, R, Y);

//dlr = min(dl, dr)
if (dl > dr) {
dlr = dr;
*pt1 = ptr3;
*pt2 = ptr4;
}
else {
dlr = dl;
*pt1 = ptr1;
*pt2 = ptr2;
}

//closest pair distnace dmin = combine(Yid, L, R, dlr,...)
dmin = combine(Yid, L, R, M, dlr, X, Xid, TMP, Y, &point1, &point2);

if (dmin > dlr) {
dmin = dlr;
}
else {
*pt1 = point1;
*pt2 = point2;
}
}
//dl, dr, cpc중 거리가 가장 짧은 쌍

return dmin;
}