最近点对点问题.cpp

递归和分治法解一系列经典算法

//最接近点对问题

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define TOTAL 10000

struct PointX
{
	bool operator <= (PointX a) const
	{
		return (x <= a.x);
	}
	float x, y;
};

struct PointY
{
	bool operator <= (PointY a) const
	{
		return (y <= a.y);
	}
	int p;  // 同一点在数组X中的坐标
	float x, y;
};

template <class Type>
inline double distance(const Type &u, const Type &v)
{
	float dx = u.x - v.x;
	float dy = u.y - v.y;
	return sqrt(dx * dx + dy * dy);
}

template <class Type>
void Merge(Type Y[], Type Z[], int l, int m, int r)
{
	Type *a = &Z[l];
	Type *b = &Z[m + 1];
	Type *c = &Y[l];
	int anum = m - l + 1, ap = 0;
	int bnum = r - m, bp = 0;
	int cnum = r - l + 1, cp = 0;
	while (ap < anum && bp < bnum)
	{
		if (a[ap] <= b[bp]) c[cp++] = a[ap++];
		else c[cp++] = b[bp++];
	}
	while (ap < anum) c[cp++] = a[ap++];
	while (bp < bnum) c[cp++] = b[bp++];
}

template <class Type>
void MergeSort(Type a[], int n)
{
	if (n <= 1) return;
	int anum = n / 2;
	Type *b = a + anum;
	int bnum = n - anum;
	MergeSort(a, anum);
	MergeSort(b, bnum);
	Type *c = new Type[n];
	Merge(c, a, 0, anum - 1, n - 1);
	for (int i=0; i<n; i++) a[i] = c[i];
	delete []c;
}

void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX &a, PointX &b, float &d)
{
	if (r - l == 1)  // two node
	{
		a = X[l];
		b = X[r];
		d = (float) distance(X[l], X[r]);
		return;
	}
	if (r - l == 2)
	{
		float d1 =(float) distance(X[l], X[l + 1]);
		float d2 =(float) distance(X[l + 1], X[r]);
		float d3 =(float) distance(X[l], X[r]);
		if (d1 <= d2 && d1 <= d3)
		{
			a = X[l];
			b = X[l + 1];
			d = d1;
		}
		else if (d2 <= d3)
		{
			a = X[l + 1];
			b = X[r];
			d = d2;
		}
		else
		{
			a = X[l];
			b = X[r];
			d = d3;
		}
		return;
	}
	int mid = (l + r) / 2;
	int f = l;
	int g = mid + 1;
	for (int i=l; i<=r; i++)
	{
		if (Y[i].p > mid) Z[g++] = Y[i];
		else Z[f++] = Y[i];
	}
	closest(X, Z, Y, l, mid, a, b, d);
	float dr;
	PointX ar, br;
	closest(X, Z, Y, mid + 1, r, ar, br, dr);
	if (dr < d)
	{
		a = ar;
		b = br;
		d = dr;
	}
	Merge(Y, Z, l, mid, r);  // 重构数组Y
	int k = l;
	for (i=l; i<=r; i++)
		if (fabs(Y[mid].x - Y[i].x) < d) Z[k++] = Y[i];
	int j;
	for (i=l; i<k; i++)
	{
		for (j=i+1; j<k && Z[j].y-Z[i].y<d; j++)
		{
			float dp =(float) distance(Z[i], Z[j]);
			if (dp < d)
			{
				d = dp;
				a = X[Z[i].p];
				b = X[Z[j].p];
			}
		}
	}
} 

bool Cpair2(PointX X[], int n, PointX &a, PointX &b, float &d)
{
	if (n < 2) return false;
	MergeSort(X, n);
	PointY *Y = new PointY[n];
	for (int i=0; i<n; i++)
	{
		Y[i].p = i;  // 同一点在数组X中的坐标
		Y[i].x = X[i].x;
		Y[i].y = X[i].y;
	}
	MergeSort(Y, n);
	PointY *Z = new PointY[n];
	closest(X, Y, Z, 0, n - 1, a, b, d);
	delete []Y;
	delete []Z;
	return true;
}

void main()
{
	srand(time(0));
	PointX X[TOTAL];
	int i, j, n = TOTAL;
	PointX a, b;
	float d;
	for (i=0; i<n; i++)
	{
		X[i].x = (float)rand()/1000;
		X[i].y = (float)rand()/1000;
	}
	a = X[0];
	b = X[1];
	d =(float) distance(X[0], X[1]);
	for (i=0; i<n; i++)
	{
		for (j=i+1; j<n; j++)
		{
			if (distance(X[i], X[j]) < d)
			{
				a = X[i];
				b = X[j];
				float d =(float) distance(X[i], X[j]);
			}
		}
	}
	printf("(%.1f, %.1f) (%.1f, %.1f)  %.3f\n",a.x,a.y,b.x,b.y,d);
	Cpair2(X, n, a, b, d);
	printf("(%.1f, %.1f) (%.1f, %.1f)  %.3f\n",a.x,a.y,b.x,b.y,d);
}


/*
begin
	if |S|=2 then δ:=S中这2点的距离
	else
		if |S|=0 then δ:=∞
		else
		begin
		1. m:=S中各点x坐标值的中位数;
		构造S1和S2, 使S1={p∈S|px≤m}和S2={p∈S|px>m};
		2.δ1:=CPAIR2(S1);δ2:=CPAIR2(S2);
		3.δm:=min(δ1,δ2);
		4.设P1是S1中距垂直分割线l的距离在δm之内的所有点组成的集合,
		P2是S2中距分割线l的距离在δm之内所有点组成的集合。将P1和P2中
		的点依其y坐标值从小到大排序,并设P1*和P2*是相应的已排好序的点列;
		5.通过扫描P1*以及对于P1*中每个点检查P2*中与其距离在δm之内的所
		有点(最多6个)可以完成合并。当P1*中的扫描指针逐次向上移动时,P2*
 		中的扫描指针可在宽为2δm的一个区间内移动。设δl是按这种扫描方
		式找到的点对间的最小距离;
		6.δ=min(δm,δl);
		end;
	return(δ);
end;
*/