最近对问题
最近点对问题要求在一个包含n个点的集合中,找出距离最近的两个点
蛮力法
很显然,蛮力法应该是这样,分别计算每一对点之间的距离,然后找出距离最小的那一对。时间复杂度为O(n2)
double BruteForceCloestPoints(Node a[])
{
int mindist = inf;
for(int i = 0;i<n-1;++i){
for(int j = i+1;j<n;++j){
double curdist = sqrt(pow(a[i].x-a[j].x,2)+pow(a[i].y-a[j].y,2));
mindist = min(mindist,curdist);
}
}
return mindist;
}
分治法
我们将点按照x轴排序后,就会发现,最左边的那些结点和最右边的那些结点显然不是最短距离,所以最短距离肯定存在于某个区间上面,且子问题没有交集,所以我们可以使用分治法,将区间一分为二,然后求解子问题。
当然有可能存在一个点在中线左边,一个点在中线右边,这时我们求完子问题后得到最小距离还有再求位于两端的距离小于d的点的最小距离与子问题得出的最小距离比较。当n<=3的时候,返回由蛮力算法求出的最小距离
typedef struct
{
int x;
int y;
}point;
double distance(const point& a, const point& b)
{
return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
}
double merge(point a[], int low, int high,double d)
{
double mindist = d;
vector<point> left, right;
int m = a[(low + high) / 2].x;
for (int i = low; i <= high; ++i){
if (a[i].x - m < 0 && a[i].x - m > -d)
left.push_back(a[i]);
if (a[i].x - m > 0 && a[i].x - m < d)
right.push_back(a[i]);
}
for (auto& p : left) {
for (auto& q : right) {
mindist = min(mindist, distance(p,q));
}
}
return mindist;
}
double EfficientClosestPair(point a[], int low, int high)
{
if (high - low < 2)
return distance(a[low],a[low + 1]);
if (high - low < 3)
return min(distance(a[low],a[low+1]), min(distance(a[low], a[low + 2]), distance(a[low+1], a[low + 2])));
int mid = (low + high) / 2;
double d1 = EfficientClosestPair(a, low, mid);
double d2 = EfficientClosestPair(a, mid + 1, high);
double d = min(d1, d2);
return merge(a, low, high,d);
}
凸包问题
凸包问题被视为计算几何中最重要的问题之一,定义如下:
在平面或者高维空间的一个给定点集合中寻找凸包
该问题之所以如此突出,是由于许多各种各样的应用要么本身就是凸包问题,要么其中一部分需要按照凸包问题来解决。大多数此类应用都是基于这样一个原理:凸包能方便地提供目标形状或给定数据集的一个近似。
凸集合的概念:
对于平面上的一个点集合(有限的或无限的),如果以集合中任意两点p和q为端点的线段上所有点都属于该集合,我们说这个集合是凸的。
下面这个凸包的正式定义可以应用于任意集合,包括那些正好位于一条直线上的点的集合。
一个点集合S的凸包是包含S的最小凸集合(“最小”意指,S的凸包一定是所有包含S的凸集合的子集)
如果S是凸的,它的凸包很明显是它本身。如果S是两个点组成的集合,它的凸包是连接着两个点的线段。如果S是由三个不同线的点组成的集合,它的凸包是以这三个点为顶点的三角形;如果三点同线,凸包是以距离最远的两个点为端点的线段。对于更大集合的凸包的例子,请见下图
凸包问题是为一个有n个点的集合构造凸包问题。为了解决该问题,需要找出某个点,它们将作为这个集合的凸多边形的顶点。
贪心法
对于一个n个点集合中的两个点pi和pj,当且仅当该集合中的其他点都位于穿过这两点的直线的同一边时,它们的连线是该集合凸包边界的一部分。对每一对点都做一遍检验之后,满足条件的线段构成了该凸包的边界
为了实现这个算法,需要用到一些解析集合的基本知识。首先,在坐标平面上穿过两个点(x1,y1),(x2,y2)的直线是由下列方程定义的:
ax+by = c
其中 a = y2 - y1,b = x1-x2 ,c = x1y2 - y1x2。
其次,这样一根直线把平面分为两个半平面:其中一个半平面中的点都满足ax+by>c,而另一个半平面中的点都满足ax+by<c。因此,为了检验某些点是否位于这条直线的同一边,只需要把每个点代入ax+by-c,检验这个表达式的符号是否相同。
typedef struct
{
int x;
int y;
}point;
typedef struct
{
point one;
point two;
}Edge;
void BruteForceClosure(const vector<point>& invec,vector<point>& outvec)
{
int a, b, c;
int last;
bool first = true;
bool can = true;
for (int i = 0; i < invec.size() - 1; ++i) {
for (int j = i + 1; j < invec.size(); ++j) {
a = invec[j].y - invec[i].y;
b = invec[i].x - invec[j].x;
c = invec[i].x * invec[j].y - invec[i].y * invec[j].x;
for (int k = 0; k < invec.size(); ++k) {
if (k != i && k != j) {
int ck = a * invec[k].x + b * invec[k].y - c;
if (first) {
first = false;
last = ck;
}
else {
if (last * ck < 0) {
can = false;
break;
}
}
}
}
if (can) {
Edge e;
e.one = invec[i];
e.two = invec[j];
outvec.emplace_back(e);
}
}
}
}
分治法(快包算法)
假设集合S是平面上上n>1个点构成的。我们还假设这些点按照他们的x轴坐标升序排列的,如果x轴坐标相同,则按照y轴坐标升序排列。不难证明这样一个几何学上的明显事实,最左面的点p1和最右面的点pn一定是该集合的凸包顶点。这条直线把点分为两个集合:S1是位于直线左侧或在直线上的点构成的集合,S2是位于直线右侧或在直线上的点构成的集合。
S的凸包的边界是由下面两个多角形链条构成的:一条“上”边界和一条“下”边界。“上”边界称为上包,是一系列线段的序列,“下”边界称为下包,是一系列线段的序列
快包算法:
到现在我们还需要具备一个几何知识, 即:
如果我们有三个点A(x1,y1),B(x2,y2), C(x3,y3)
我们可以用以下这个行列式求面积 (也就是我们常说的"叉积")
这个行列式对于平面上任意三角形, 求解面积都很方便, 所得结果是三角形ABC面积的两倍 (可以回忆高中知识~)
而且, 当A->B->C为顺时针顺序时, 该值为正, 反之则为负
struct point
{
int x;
int y;
}p[nmax];
int visited[nmax];
int Djudge(point a1, point a2, point a3)
{
int calculate = a1.x * a2.y + a3.x * a1.y + a2.x * a3.y - a3.x * a2.y - a2.x * a1.y - a1.x * a3.y;
return calculate;
}
bool cmp(const point& a, const point& b)
{
if (a.x != b.x)
return a.x < b.x;
else
return a.y < b.y;
}
void DealConvex(int first, int last)
{
int max = 0, index = -1;
//求上包
if (first < last) {
//注意两端没有必要计算
for (int i = first + 1; i < last; ++i) {
int calcu = Djudge(p[first], p[i], p[last]);
if (calcu == 0) visited[i] = 1;
if (calcu > nmax) {
max = calcu;
index = i;
}
}
}
//求下包
else {
for (int i = first - 1; i > last; --i) {
int calcu = Djudge(p[first], p[i], p[last]);
if (calcu == 0) visited[i] = 1;
if (calcu > nmax) {
max = calcu;
index = i;
}
}
}
if (index = -1) {
visited[index] = 1;
DealConvex(first, index);
DealConvex(index, last);
}
}
void PrintPoint(int n) {
int mark[nmax];
int t = 0;
point ans[nmax];
for (int i = 0; i < n; i++)
{
if (visited[i] == 1)
{
ans[t].x = p[i].x;
ans[t].y = p[i].y;
t++;
}
}
//顺时针输出
mark[0] = mark[t - 1] = 1; //数组mark避免重复检查降低效率
for (int i = 1; i < t - 1; i++)
{
mark[i] = 0;
}
cout << ans[0].x << " " << ans[0].y << endl;
for (int i = 1; i < t - 1; i++)
{
int d = Djudge(ans[0], ans[t - 1], ans[i]);
if (d >= 0)
{
cout << ans[i].x << " " << ans[i].y << endl;
mark[i] = 1;
}
}
cout << ans[t - 1].x << " " << ans[t - 1].y << endl;
for (int i = 1; i < t; i++)
{
if (mark[i] != 1)
{
int d = Djudge(ans[0], ans[t - 1], ans[i]);
if (d < 0)
{
cout << ans[i].x << " " << ans[i].y << endl;
}
}
}
}