KD-Tree 算法总结

An_Account · 2018-08-10 21:46:30 · 个人记录

KD-Tree 算法总结

KD-Tree 是什么

简而言之，KD-Tree是一种能维护高维数据空间的结构,主要支持几个操作：

1.插入点 2.进行距离查询（例如：查询距离某个点第k近的点）

建树（build）

考虑两个问题：

1.如何选择划分的维度，使得KD-Tree的结构尽可能更优秀 2.如何选择当前的根节点，使得子树的深度尽量最小

显然，按照1,2,3,\dots k维的顺序来划分并不一定最好，我们考虑这样一种情况：假如在一个二维平面上，这k个点排成一条与x轴平行的直线，那么按y的大小来划分就会出现一些很尴尬的情况。

其实还有一种划分方法，我们不按1,2,\dots k这样顺序划分，而是按方差最大的那一维划分，将点最分散的那一维化成两部分，这也是我们希望看到的结果。

其实在实际应用中，顺序划分是一种最常见的方式，因为求方差的时间复杂度很高，而顺序划分对于随机数据来说表现也很出色

至于第二个问题，很明显的一件事就是我们可以选择中位数，左边一半，右边一半，这样很平均。

因此上面那个图的树应该建成这样：

关于怎么求中位数，algorithm头文件中很贴心的为我们准备了一个函数nth_element

template<class _RanIt> inline  
    void nth_element(_RanIt _First, _RanIt _Nth, _RanIt _Last)

对于[First,Last)区间内的数重新排序，使得位置为Nth的值是第Nth小值

程序实现如下

bool del;
bool cmp(Point p1,Point p2) {
    if (!del) return (p1.x<p2.x||(p1.x==p2.x&&p1.y<p2.y));
    return (p1.y<p2.y||(p1.y==p2.y&&p1.x<p2.x));
}

两点间距离

double dis(Point p1,Point p2) {
    return (double)(p1.x-p2.x)*(p1.x-p2.x)+(double)(p1.y-p2.y)*(p1.y-p2.y);
}

void pushup(int rt) {
    T[rt].r1.x=min(min(T[T[rt].ch[0]].r1.x,T[T[rt].ch[1]].r1.x),T[rt].r1.x);
    T[rt].r1.y=min(min(T[T[rt].ch[0]].r1.y,T[T[rt].ch[1]].r1.y),T[rt].r1.y);
    T[rt].r2.x=max(max(T[T[rt].ch[0]].r2.x,T[T[rt].ch[1]].r2.x),T[rt].r2.x);
    T[rt].r2.y=max(max(T[T[rt].ch[0]].r2.y,T[T[rt].ch[1]].r2.y),T[rt].r2.y);
}

初始化

void init() {
    T[0].r1=Point(0x3f3f3f3f,0x3f3f3f3f),T[0].r2=Point(-0x3f3f3f3f,-0x3f3f3f3f);
}

建树

int build(int l,int r,int d) {
    if (l>r) return 0;
    del=d;int mid=(l+r)>>1,at=++ncnt;
    nth_element(ps+l,ps+mid,ps+r+1,cmp);
    T[at]=Tree(ps[mid],mid);
    T[at].ch[0]=build(l,mid-1,d^1),T[at].ch[1]=build(mid+1,r,d^1);
    pushup(at);return at;
}

优先队列中的点（注意优先队列是大根堆）

struct node {
    double dis;int id;
    node(double dis=0,int id=0):dis(dis),id(id) {}
    bool operator < (node b) const {
        return dis>b.dis||(dis==b.dis&&id<b.id);
    }
};

查询

priority_queue<node> q;
void query(int now,Point p) {
    if (!now) return;
    node st=node(dis(T[now].p,p),T[now].p.id);
    if (st<q.top()) q.pop(),q.push(st);
    double dis[2]={T[T[now].ch[0]].dis(p),T[T[now].ch[1]].dis(p)};
    int next=dis[0]<dis[1];query(T[now].ch[next],p);
    if (node(dis[next^1],T[T[now].ch[next^1]].id)<q.top()) query(T[now].ch[next^1],p);
}

总代码实现如下

#include <iostream>
#include <cstdio>
#include <cstring>
#include <queue>
#include <algorithm>
using namespace std;
#define N 100010
const double inf=1e300;
struct Point {
    int x,y,id;
    Point (int x=0,int y=0):x(x),y(y) {}
}ps[N];
bool del;
bool cmp(Point p1,Point p2) {
    if (!del) return (p1.x<p2.x||(p1.x==p2.x&&p1.y<p2.y));
    return (p1.y<p2.y||(p1.y==p2.y&&p1.x<p2.x));
}
double dis(Point p1,Point p2) {
    return (double)(p1.x-p2.x)*(p1.x-p2.x)+(double)(p1.y-p2.y)*(p1.y-p2.y);
}
struct Tree {
    int ch[2],id;Point p,r1,r2;
    Tree(Point p=Point(),int id=0):p(p),r1(p),r2(p),id(id) {}
    double dis(Point p) {
        if (!id) return -inf;
        return max(max(::dis(p,r1),::dis(p,r2)),max(::dis(p,Point(r1.x,r2.y)),::dis(p,Point(r2.x,r1.y))));
    }
}T[N];
void pushup(int rt) {
    T[rt].r1.x=min(min(T[T[rt].ch[0]].r1.x,T[T[rt].ch[1]].r1.x),T[rt].r1.x);
    T[rt].r1.y=min(min(T[T[rt].ch[0]].r1.y,T[T[rt].ch[1]].r1.y),T[rt].r1.y);
    T[rt].r2.x=max(max(T[T[rt].ch[0]].r2.x,T[T[rt].ch[1]].r2.x),T[rt].r2.x);
    T[rt].r2.y=max(max(T[T[rt].ch[0]].r2.y,T[T[rt].ch[1]].r2.y),T[rt].r2.y);
}
int ncnt;
void init() {
    T[0].r1=Point(0x3f3f3f3f,0x3f3f3f3f),T[0].r2=Point(-0x3f3f3f3f,-0x3f3f3f3f);
}
int build(int l,int r,int d) {
    if (l>r) return 0;
    del=d;int mid=(l+r)>>1,at=++ncnt;
    nth_element(ps+l,ps+mid,ps+r+1,cmp);
    T[at]=Tree(ps[mid],mid);
    T[at].ch[0]=build(l,mid-1,d^1),T[at].ch[1]=build(mid+1,r,d^1);
    pushup(at);return at;
}
struct node {
    double dis;int id;
    node(double dis=0,int id=0):dis(dis),id(id) {}
    bool operator < (node b) const {
        return dis>b.dis||(dis==b.dis&&id<b.id);
    }
};
priority_queue<node> q;
void query(int now,Point p) {
    if (!now) return;
    node st=node(dis(T[now].p,p),T[now].p.id);
    if (st<q.top()) q.pop(),q.push(st);
    double dis[2]={T[T[now].ch[0]].dis(p),T[T[now].ch[1]].dis(p)};
    int next=dis[0]<dis[1];query(T[now].ch[next],p);
    if (node(dis[next^1],T[T[now].ch[next^1]].id)<q.top()) query(T[now].ch[next^1],p);
}
int main() {
    init();int n,m;scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%d%d",&ps[i].x,&ps[i].y),ps[i].id=i;
    build(1,n,0),scanf("%d",&m);
    while (m--) {
        int x,y,k;scanf("%d%d%d",&x,&y,&k);
        while (!q.empty()) q.pop();
        for (int i=1;i<=k;i++) q.push(node(-inf));
        query(1,Point(x,y)),printf("%d\n",q.top().id);
    }
    return 0;
}