kd树【@Abandon】

刚开始学习，介绍先搁着~等理解透彻了再来写~~~

我是学习的mzry1992（UESTC_Izayoi ~）------http://www.mzry1992.com/blog/miao/kd%E6%A0%91.html

先去看mzry1992大牛博客里的讲解吧。。。

再附两篇论文：（看英文看得好爽。。。~@.@）

《An intoductory tutorial on kd-trees》　　★（里面就介绍了kd-tree和nearest neighbour algorithm（最*邻算法）、Q nearest neighbour（Q*邻））

《Range Searching Using Kd-Tree.》

模板：

kd-tree模板

/*
    HDOJ 2966
    KD-Tree模板
*/
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <stack>
#include <queue>
#include <map>
#include <algorithm>
#include <string>
#include <cstring>
#define MID(x, y) ( (x + y)>>1 )

using namespace std;
typedef long long LL;

//KD-Tree模板
const int N=100005;
LL res;
struct Point
{
    int x, y;        //点是二维的，此时是2D-Tree
};

LL dist2(const Point &a, const Point &b)            //距离的*方
{
    return LL(a.x - b.x) * LL(a.x - b.x) + LL(a.y - b.y) * LL(a.y - b.y);
}

bool cmpX(const Point &a, const Point &b)
{
    return a.x < b.x;
}
bool cmpY(const Point &a, const Point &b)
{
    return a.y < b.y;
}

struct KDTree        //很崇拜这种销魂的建树方法啊~0.0~很抽象很强大------p数组已经代表了KD-Tree了，神马左右子树全省了，OOOOOrrz!!!……
{
    Point p[N];        //空间内的点
    int Div[N];        //记录区间是按什么方式划分（分割线*行于x轴还是y轴, ==1*行y轴切;==0*行x轴切）

    void build(int l, int r)            //记得先把p备份一下。
    {
        if (l > r)    return;
        int mid=MID(l, r);
        int minX, minY, maxX, maxY;
        minX = min_element(p + l, p + r + 1, cmpX)->x;
        minY = min_element(p + l, p + r + 1, cmpY)->y;
        maxX = max_element(p + l, p + r + 1, cmpX)->x;
        maxY = max_element(p + l, p + r + 1, cmpY)->y;
        Div[mid] = (maxX - minX >= maxY - minY);
        nth_element(p + l, p + mid, p + r + 1, Div[mid] ? cmpX : cmpY);
        build(l, mid - 1);
        build(mid+1, r);
    }

    void find(int l, int r, Point a)                //查找最*点的*方距离
    {
        if (l > r)    return;
        int mid = MID(l, r);
        LL dist = dist2(a, p[mid]);
        if (dist > 0)   //如果有重点不能这么判断
            res = min(res, dist);
        LL d = Div[mid] ? (a.x - p[mid].x) : (a.y - p[mid].y);
        int l1, l2, r1, r2;
        l1 = l , l2 = mid + 1;
        r1 = mid - 1, r2 = r;
        if (d > 0)
            swap(l1, l2), swap(r1, r2);
        find(l1, r1, a);
        if (d * d < res)
            find(l2, r2, a);
    }
};

kd-tree入门题：

　　HDOJ 2966 In case of failure (最*邻，直接套模板~)

　　查找*面点最*点的距离（此题中是距离的*方）

HDOJ 2966

  1 /*
  2     HDOJ 2966
  3     KD-Tree模板
  4 */
  5 #include <iostream>
  6 #include <cstdio>
  7 #include <cstdlib>
  8 #include <cmath>
  9 #include <vector>
 10 #include <stack>
 11 #include <queue>
 12 #include <map>
 13 #include <algorithm>
 14 #include <string>
 15 #include <cstring>
 16 #define MID(x, y) ( (x + y)>>1 )
 17 
 18 using namespace std;
 19 typedef long long LL;
 20 
 21 //KD-Tree模板
 22 const int N=100005;
 23 LL res;
 24 struct Point
 25 {
 26     int x, y;        //点是二维的，此时是2D-Tree
 27 };
 28 
 29 LL dist2(const Point &a, const Point &b)            //距离的*方
 30 {
 31     return LL(a.x - b.x) * LL(a.x - b.x) + LL(a.y - b.y) * LL(a.y - b.y);
 32 }
 33 
 34 bool cmpX(const Point &a, const Point &b)
 35 {
 36     return a.x < b.x;
 37 }
 38 bool cmpY(const Point &a, const Point &b)
 39 {
 40     return a.y < b.y;
 41 }
 42 
 43 struct KDTree        //很崇拜这种销魂的建树方法啊~0.0~很抽象很强大------p数组已经代表了KD-Tree了，神马左右子树全省了，OOOOOrrz!!!……
 44 {
 45     Point p[N];        //空间内的点
 46     int Div[N];        //记录区间是按什么方式划分（分割线*行于x轴还是y轴, ==1*行y轴切;==0*行x轴切）
 47 
 48     void build(int l, int r)            //记得先把p备份一下。
 49     {
 50         if (l > r)    return;
 51         int mid=MID(l, r);
 52         int minX, minY, maxX, maxY;
 53         minX = min_element(p + l, p + r + 1, cmpX)->x;
 54         minY = min_element(p + l, p + r + 1, cmpY)->y;
 55         maxX = max_element(p + l, p + r + 1, cmpX)->x;
 56         maxY = max_element(p + l, p + r + 1, cmpY)->y;
 57         Div[mid] = (maxX - minX >= maxY - minY);
 58         nth_element(p + l, p + mid, p + r + 1, Div[mid] ? cmpX : cmpY);
 59         build(l, mid - 1);
 60         build(mid+1, r);
 61     }
 62 
 63     void find(int l, int r, Point a)                //查找最*点的*方距离
 64     {
 65         if (l > r)    return;
 66         int mid = MID(l, r);
 67         LL dist = dist2(a, p[mid]);
 68         if (dist > 0)   //如果有重点不能这么判断
 69             res = min(res, dist);
 70         LL d = Div[mid] ? (a.x - p[mid].x) : (a.y - p[mid].y);
 71         int l1, l2, r1, r2;
 72         l1 = l , l2 = mid + 1;
 73         r1 = mid - 1, r2 = r;
 74         if (d > 0)
 75             swap(l1, l2), swap(r1, r2);
 76         find(l1, r1, a);
 77         if (d * d < res)
 78             find(l2, r2, a);
 79     }
 80 };
 81 
 82 Point pp[N];
 83 KDTree kd;
 84 
 85 int main()
 86 {
 87     //freopen("test.txt","r+",stdin);
 88     int t;
 89     scanf("%d", &t);
 90     while(t--)
 91     {
 92         int n;
 93         scanf("%d", &n);
 94         for (int i = 0; i < n; i++)
 95         {
 96             scanf("%d%d", &pp[i].x, &pp[i].y);
 97             kd.p[i] = pp[i];
 98         }
 99 
100         kd.build(0, n-1);
101         for (int i = 0; i < n; i++)
102         {
103             res = 9223372036854775807LL;
104             kd.find(0, n - 1, pp[i]);
105             printf("%I64d\n", res);
106         }
107     }
108     return 0;
109 }

　　HDOJ 4347 The Closest M Points (Q*邻)

　　与上题不同的是，一是k维（这个好处理~）,二是求最*的m个点而不单是最*点了。这个也好处理~递归查找时处理的时候采取如下策略：如果当前找到的点小于k个，那么两个区间都要处理。。否则根据当前找到的第k个点决定是否去另外一个区间，如果目标点到分界线的距离大于等于已经找到的第k远的点，那么就不用查找另一个分界了。。。更新答案可以用一个大小为k的堆去维护（一个最大堆，一旦超过k个点就把最大的扔掉）。。。

　　2s，速度还是不错的^_^~

6811512

2012-09-21 16:58:20

Accepted

4347

2421MS

3848K

4152 B

G++

AbandonZHANG

HDOJ 4347

  1 /*
  2     HDOJ 2966
  3     KD-Tree模板
  4 */
  5 #include <iostream>
  6 #include <cstdio>
  7 #include <cstdlib>
  8 #include <cmath>
  9 #include <vector>
 10 #include <stack>
 11 #include <queue>
 12 #include <set>
 13 #include <map>
 14 #include <algorithm>
 15 #include <string>
 16 #include <cstring>
 17 #define MID(x,y) ( (x + y)>>1 )
 18 
 19 using namespace std;
 20 typedef long long LL;
 21 
 22 //KD-Tree模板
 23 const int N=50005;
 24 
 25 struct Point
 26 {
 27     int x[5];
 28     LL dis;
 29     Point()
 30     {
 31         for (int i = 0; i < 5; i++)
 32             x[i] = 0;
 33         dis = 9223372036854775807LL;
 34     }
 35     friend bool operator < (const Point &a, const Point &b)
 36     {
 37         return a.dis < b.dis;
 38     }
 39 };
 40 priority_queue <Point, vector<Point> > res;
 41 
 42 LL dist2(const Point &a, const Point &b, int k)            //距离的*方,开根号很耗时，能不开就不开
 43 {
 44     LL ans = 0;
 45     for (int i = 0; i < k; i++)         //一开始这儿写的i < 5，WA了N次。。。原本以为只要初始值设0就无所谓,
 46         ans += (a.x[i] - b.x[i]) * (a.x[i] - b.x[i]);   //但发现Point有个全局变量，在多case下会出错。。。
 47     return ans;
 48 }
 49 
 50 int ddiv;
 51 bool cmpX(const Point &a, const Point &b)
 52 {
 53     return a.x[ddiv] < b.x[ddiv];
 54 }
 55 
 56 struct KDTree        //很崇拜这种销魂的建树方法啊~0.0~很抽象很强大------p数组已经代表了KD-Tree了，神马左右子树全省了，OOOOOrrz!!!……
 57 {
 58     Point p[N];        //空间内的点
 59     int Div[N];        //记录区间是按什么方式划分（分割线*行于x轴还是y轴, ==1*行y轴切;==0*行x轴切）
 60     int k;          //维数
 61     int m;          //*邻
 62 
 63     int getDiv(int l, int r)        //寻找区间跨度最大的划分方式
 64     {
 65         map <int, int> ms;
 66         int minx[5],maxx[5];
 67         for (int i = 0; i < k; i++)
 68         {
 69             ddiv = i;
 70             minx[i] = min_element(p + l, p + r + 1, cmpX)->x[i];
 71             maxx[i] = max_element(p + l, p + r + 1, cmpX)->x[i];
 72             ms[maxx[i] - minx[i]] = i;
 73         }
 74         map <int ,int>::iterator pm = ms.end();
 75         pm--;
 76         return pm->second;
 77     }
 78 
 79     void build(int l, int r)            //记得先把p备份一下。
 80     {
 81         if (l > r)    return;
 82         int mid = MID(l,r);
 83         Div[mid] = getDiv(l,r);
 84         ddiv = Div[mid];
 85         nth_element(p + l, p + mid, p + r + 1, cmpX);
 86         build(l, mid - 1);
 87         build(mid + 1, r);
 88     }
 89 
 90     void findk(int l, int r, Point a)                //k(m)*邻，查找k*点的*方距离
 91     {
 92         if (l > r)    return;
 93         int mid = MID(l,r);
 94         LL dist = dist2(a, p[mid], k);
 95         if (dist >= 0)
 96         {
 97             p[mid].dis = dist;
 98             res.push(p[mid]);
 99             while ((int)res.size() > m)
100                 res.pop();
101         }
102         LL d = a.x[Div[mid]] - p[mid].x[Div[mid]];
103         int l1, l2, r1, r2;
104         l1 = l , l2 = mid + 1;
105         r1 = mid - 1, r2 = r;
106         if (d > 0)
107             swap(l1, l2), swap(r1, r2);
108         findk(l1, r1, a);
109         if ((int)res.size() < m || d*d < res.top().dis )
110             findk(l2, r2, a);
111     }
112 };
113 
114 Point pp[N];
115 KDTree kd;
116 
117 int main()
118 {
119 //    freopen("test.txt","r+",stdin);
120 //    freopen("ans.txt","w+",stdout);
121 
122     int n;
123     while(scanf("%d%d", &n, &kd.k)!=EOF)
124     {
125         for (int i = 0; i < n; i++)
126             for (int j = 0; j < kd.k; j++)
127             {
128                 scanf("%d", &pp[i].x[j]);
129                 kd.p[i] = pp[i];
130             }
131         kd.build(0, n - 1);
132         int t;
133         scanf("%d", &t);
134         while(t--)
135         {
136             Point a;
137             for (int i = 0; i < kd.k; i++)
138                 scanf("%d", &a.x[i]);
139             scanf("%d", &kd.m);
140             kd.findk(0, n - 1, a);
141             printf("the closest %d points are:\n", kd.m);
142             Point ans[11];
143             for (int i = 0; !res.empty(); i++)
144             {
145                 ans[i] = res.top();
146                 res.pop();
147             }
148             for (int i = kd.m - 1; i >= 0; i--)
149             {
150                 for (int j = 0; j < kd.k - 1; j++)
151                     printf("%d ", ans[i].x[j]);
152                 printf("%d\n", ans[i].x[kd.k - 1]);
153             }
154         }
155 
156     }
157     return 0;
158 }

(未完待续。。。)

举杯独醉，饮罢飞雪，茫然又一年岁。 ------AbandonZHANG

kd树 【@Abandon】

kd树【@Abandon】