kdtree

注意：这是一篇个人学习笔记，如果有人因为某些原因点了进来并且要看一下，请一定谨慎地阅读，因为可能存在各种奇怪的错误，如果有人发现错误请指出谢谢！

资料待看：

https://www.cnblogs.com/yangsongyi/p/10000756.html

http://nfile.oss-cn-hangzhou.aliyuncs.com/KDTree.pdf

kdtree

将一些点组织成二叉树，对于每一个节点，有一个特定的划分维度，其左子树中点都在这一维上小于等于它，其右子树中点都在这一维上大于等于它。（等于的归于哪一边没有办法规定，解决问题的时候可能需要注意）这个划分维度一般是按层轮流划分（比如二维平面，树的第1,3,5,7,..层按x轴，第2,4,6,8,..层按y轴）

最近点对

例题：洛谷P1429平面最近点对

查询离某一个点p最近的点：在kdtree上dfs，进行一些剪枝。

到某个节点时，优先向它的两个子节点中，"估计"更优的一个点走。

到任何一个节点都可以用来更新当前最优解。

然后就不要走那些一定不会产生比当前最优解更优的解的区域（也就是那些离p的最小可能距离也要大于等于当前最近距离的区域）。

怎么样”估计“更优呢？网上有方法1，每次优先走包含p点的那个半空间（https://blog.csdn.net/u013534123/article/details/80952174，https://codeforces.com/blog/entry/54080）,然后如果p到用于划分的“面”的距离<=当前最短距离，才去走另外一个半空间（实际实现见代码）

 1 #include<cstdio>
 2 #include<ctime>
 3 #include<algorithm>
 4 #include<cstring>
 5 #include<vector>
 6 #include<cmath>
 7 using namespace std;
 8 #define fi first
 9 #define se second
10 #define mp make_pair
11 #define pb push_back
12 typedef long long ll;
13 typedef unsigned long long ull;
14 #define cmin(a,b) (((b)<(a))?void((a)=(b)):void())
15 template<class T>
16 inline T sqr(T x){return x*x;}
17 namespace K
18 {
19     struct N
20     {
21         double d[2];
22         int lc,rc;
23     }d[200011],p[200011];
24     int W;//W:当前维度
25     bool c1(const N &a,const N &b)
26     {
27         return a.d[W]<b.d[W];
28     }
29     double dis2(const N &a,const N &b)
30     {
31         return sqr(a.d[0]-b.d[0])+sqr(a.d[1]-b.d[1]);
32     }
33 #define LC (d[u].lc)
34 #define RC (d[u].rc)
35     int build(int l,int r,int w)//w:当前维度
36     {
37         W=w;
38         int u=(l+r)>>1;
39         nth_element(p+l,p+u,p+r+1,c1);
40         d[u]=p[u];
41         if(l!=u)    LC=build(l,u-1,w^1);
42         if(u!=r)    RC=build(u+1,r,w^1);
43         return u;
44     }
45     int x;
46     const double inf=1e18;
47     double ans;
48     void querymn(int u,int w)
49     {
50         if(!u)    return;
51         double tmp=dis2(d[u],d[x]);
52         if(u!=x)    cmin(ans,tmp);
53         int lc=LC,rc=RC;
54         if(d[x].d[w]>d[u].d[w])    swap(lc,rc);
55         querymn(lc,w^1);
56         if(sqr(d[x].d[w]-d[u].d[w])<=ans)    querymn(rc,w^1);
57     }
58 }
59 int n,rt;
60 double anss;
61 int main()
62 {
63     srand(time(0));
64     int i;
65     scanf("%d",&n);
66     {
67         using namespace K;
68         for(i=1;i<=n;++i)
69             scanf("%lf%lf",&p[i].d[0],&p[i].d[1]);
70         anss=inf;
71         rt=build(1,n,0);
72         for(i=1;i<=n;++i)
73         {
74             x=i;ans=inf;
75             querymn(rt,0);
76             anss=min(anss,ans);
77         }
78         printf("%.4f\n",sqrt(anss));
79     }
80     return 0;
81 }

View Code

还有方法2（http://hzwer.com/6954.html），维护每个区域的各维最小坐标和最大坐标，然后可以分别求出p到左右子节点的区域的某种特征（这个版本是将各个维度的最小差值直接加了起来，因为那个题是曼哈顿距离；如果是欧几里得距离，那么应该是各个维度的最小差值平方后加起来），先走较优的（这个版本就是指较小的）

可能需要特判同一个点/重合点

  1 #include<cstdio>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<vector>
  5 #include<cmath>
  6 using namespace std;
  7 #define fi first
  8 #define se second
  9 #define mp make_pair
 10 #define pb push_back
 11 typedef long long ll;
 12 typedef unsigned long long ull;
 13 template<class T1,class T2>
 14 inline T1 min1(T1 a,T2 b){return a<b?a:b;}
 15 template<class T1,class T2>
 16 inline T1 max1(T1 a,T2 b){return a<b?b:a;}
 17 #define min min1
 18 #define max max1
 19 #define cmin(a,b) (((b)<(a))?void((a)=(b)):void())
 20 #define cmax(a,b) (((a)<(b))?void((a)=(b)):void())
 21 template<class T>
 22 inline T sqr(T x){return x*x;}
 23 namespace K
 24 {
 25     struct N
 26     {
 27         double d[2],mx[2],mn[2];//自身，子树min,max
 28         int lc,rc;
 29     }d[200011],p[200011];
 30     int W;//W:当前维度
 31     bool c1(const N &a,const N &b)
 32     {
 33         return a.d[W]<b.d[W];
 34     }
 35     double dis2(const N &a,const N &b)
 36     {
 37         return sqr(a.d[0]-b.d[0])+sqr(a.d[1]-b.d[1]);
 38     }
 39 #define LC (d[u].lc)
 40 #define RC (d[u].rc)
 41     void upd(int u)
 42     {
 43         for(int i=0;i<2;++i)
 44         {
 45             d[u].mx[i]=d[u].mn[i]=d[u].d[i];
 46             if(LC)
 47             {
 48                 cmax(d[u].mx[i],d[LC].mx[i]);
 49                 cmin(d[u].mn[i],d[LC].mn[i]);
 50             }
 51             if(RC)
 52             {
 53                 cmax(d[u].mx[i],d[RC].mx[i]);
 54                 cmin(d[u].mn[i],d[RC].mn[i]);
 55             }
 56         }
 57     }
 58     int build(int l,int r,int w)//w:当前维度
 59     {
 60         W=w;
 61         int u=(l+r)>>1;
 62         nth_element(p+l,p+u,p+r+1,c1);
 63         d[u]=p[u];
 64         //for(int i=0;i<2;++i)
 65         //    d[u].mx[i]=d[u].mn[i]=d[u].d[i];
 66         if(l!=u)    LC=build(l,u-1,w^1);
 67         if(u!=r)    RC=build(u+1,r,w^1);
 68         upd(u);
 69         return u;
 70     }
 71     int x;
 72     const double inf=1e18;
 73     double ans;
 74     double getmn(int u)//x到u的子树区域的最小可能距离
 75     {
 76         double ans=0;
 77         for(int i=0;i<2;++i)
 78             ans+=sqr(max(d[x].d[i]-d[u].mx[i],0)
 79                 +max(d[u].mn[i]-d[x].d[i],0));
 80         return ans;
 81     }
 82     void querymn(int u)
 83     {
 84         double tmp=dis2(d[u],d[x]);
 85         if(u!=x)    cmin(ans,tmp);
 86         double dl=inf,dr=inf;
 87         if(LC)    dl=getmn(LC);
 88         if(RC)    dr=getmn(RC);
 89         if(dl<dr)
 90         {
 91             if(dl<ans)    querymn(LC);
 92             if(dr<ans)    querymn(RC);
 93         }
 94         else
 95         {
 96             if(dr<ans)    querymn(RC);
 97             if(dl<ans)    querymn(LC);
 98         }
 99     }
100 }
101 int n,rt;
102 double anss;
103 int main()
104 {
105     int i;
106     scanf("%d",&n);
107     {
108         using namespace K;
109         for(i=1;i<=n;++i)
110             scanf("%lf%lf",&p[i].d[0],&p[i].d[1]);
111         anss=inf;
112         rt=build(1,n,0);
113         for(i=1;i<=n;++i)
114         {
115             x=i;ans=inf;
116             querymn(rt);
117             anss=min(anss,ans);
118         }
119         printf("%.4f\n",sqrt(anss));
120     }
121     return 0;
122 }

View Code

（实际测试结果：方法1比方法2快很多（接近1/2），因为维护的东西少，且估价函数计算更简单；不清楚估价函数本身的优劣，应该是差不多的）

最远点对

这个一定要维护信息了...方法几乎和最近点对的方法2是一样的，只需要改一下估价函数（见下面K远点对）

K远点对

例题：[CQOI2016]K远点对洛谷P4357 bzoj4520

几乎和最近/远点对是一样的，就是用一个堆维护一下当前的前K远点对

（对于此题由于点对是无序的，实际要求出第2*K大的点对）

  1 #include<cstdio>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<vector>
  5 #include<cmath>
  6 #include<queue>
  7 using namespace std;
  8 #define fi first
  9 #define se second
 10 #define mp make_pair
 11 #define pb push_back
 12 typedef long long ll;
 13 typedef unsigned long long ull;
 14 template<class T1,class T2>
 15 inline T1 min1(T1 a,T2 b){return a<b?a:b;}
 16 template<class T1,class T2>
 17 inline T1 max1(T1 a,T2 b){return a<b?b:a;}
 18 template<class T>
 19 inline T abs(T a){return a>0?a:-a;}
 20 #define min min1
 21 #define max max1
 22 #define cmin(a,b) (((b)<(a))?void((a)=(b)):void())
 23 #define cmax(a,b) (((a)<(b))?void((a)=(b)):void())
 24 template<class T>
 25 inline T sqr(T x){return x*x;}
 26 int n,rt;unsigned K;
 27 namespace KK
 28 {
 29     struct N
 30     {
 31         ll d[2],mx[2],mn[2];//自身，子树min,max
 32         int lc,rc;
 33     }d[100011],p[100011];
 34     int W;//W:当前维度
 35     bool c1(const N &a,const N &b)
 36     {
 37         return a.d[W]<b.d[W];
 38     }
 39     ll dis2(const N &a,const N &b)
 40     {
 41         return sqr(a.d[0]-b.d[0])+sqr(a.d[1]-b.d[1]);
 42     }
 43 #define LC (d[u].lc)
 44 #define RC (d[u].rc)
 45     void upd(int u)
 46     {
 47         for(int i=0;i<2;++i)
 48         {
 49             d[u].mx[i]=d[u].mn[i]=d[u].d[i];
 50             if(LC)
 51             {
 52                 cmax(d[u].mx[i],d[LC].mx[i]);
 53                 cmin(d[u].mn[i],d[LC].mn[i]);
 54             }
 55             if(RC)
 56             {
 57                 cmax(d[u].mx[i],d[RC].mx[i]);
 58                 cmin(d[u].mn[i],d[RC].mn[i]);
 59             }
 60         }
 61     }
 62     int build(int l,int r,int w)//w:当前维度
 63     {
 64         W=w;
 65         int u=(l+r)>>1;
 66         nth_element(p+l,p+u,p+r+1,c1);
 67         d[u]=p[u];
 68         if(l!=u)    LC=build(l,u-1,w^1);
 69         if(u!=r)    RC=build(u+1,r,w^1);
 70         upd(u);
 71         return u;
 72     }
 73     int x;
 74     const ll inf=0x3f3f3f3f3f3f3f3f;
 75     ll getmx(int u)//x到u的子树区域的最大可能距离的平方
 76     {
 77         ll ans=0;
 78         for(int i=0;i<2;++i)
 79             ans+=sqr(max(abs(d[x].d[i]-d[u].mn[i]),
 80                 abs(d[u].mx[i]-d[x].d[i])));
 81         return ans;
 82     }
 83     priority_queue<ll,vector<ll>,greater<ll> > ans;
 84 #define OK(x) ( (ans.size()<K) || (ans.top()<(x)) )
 85     void querymx(int u)
 86     {
 87         ll tmp=dis2(d[u],d[x]);
 88         if(u!=x)
 89         {
 90             if(ans.size()<K)
 91                 ans.push(tmp);
 92             else if(ans.top()<tmp)
 93                 ans.pop(),ans.push(tmp);
 94             //printf("2t%d %lld %lld\n",u,tmp,ans.top());
 95         }
 96         ll dl=-inf,dr=-inf;
 97         if(LC)    dl=getmx(LC);
 98         if(RC)    dr=getmx(RC);
 99         if(dl>dr)
100         {
101             if(LC && OK(dl))    querymx(LC);
102             if(RC && OK(dr))    querymx(RC);
103         }
104         else
105         {
106             if(RC && OK(dr))    querymx(RC);
107             if(LC && OK(dl))    querymx(LC);
108         }
109     }
110     /*
111     void out()
112     {
113         printf("size:%d\n",int(ans.size()));
114             auto qqq=ans;
115             while(!qqq.empty())
116             {
117                 auto t=qqq.top();qqq.pop();
118                 printf("%lld\n",t);
119             }
120             {
121                 char s[233];
122                 scanf("%s",s);
123             }
124     }
125     */
126 }
127 int main()
128 {
129     int i;
130     scanf("%d%u",&n,&K);
131     K*=2;
132     {
133         using namespace KK;
134         for(i=1;i<=n;++i)
135             scanf("%lld%lld",&p[i].d[0],&p[i].d[1]);
136         rt=build(1,n,0);
137         for(i=1;i<=n;++i)
138         {
139             x=i;
140             querymx(rt);
141             //out();
142         }
143         printf("%lld\n",ans.top());
144     }
145     return 0;
146 }

View Code

改一下可以过：JZPFAR 清橙A1302 bzoj2626 洛谷P2093

（要求输出编号，注意细节）

  1 #include<cstdio>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<vector>
  5 #include<cmath>
  6 #include<queue>
  7 using namespace std;
  8 #define fi first
  9 #define se second
 10 #define mp make_pair
 11 #define pb push_back
 12 typedef long long ll;
 13 typedef unsigned long long ull;
 14 struct pli
 15 {
 16     ll fi;int se;
 17 };
 18 bool operator<(const pli &a,const pli &b)
 19 {
 20     return (a.fi<b.fi) || (a.fi==b.fi && a.se>b.se);
 21 }
 22 template<class T1,class T2>
 23 inline T1 min1(T1 a,T2 b){return a<b?a:b;}
 24 template<class T1,class T2>
 25 inline T1 max1(T1 a,T2 b){return a<b?b:a;}
 26 template<class T>
 27 inline T abs(T a){return a>0?a:-a;}
 28 #define min min1
 29 #define max max1
 30 #define cmin(a,b) (((b)<(a))?void((a)=(b)):void())
 31 #define cmax(a,b) (((a)<(b))?void((a)=(b)):void())
 32 template<class T>
 33 inline T sqr(T x){return x*x;}
 34 int n,m,rt;unsigned K;
 35 namespace KK
 36 {
 37     struct N
 38     {
 39         ll d[2],mx[2],mn[2];//自身，子树min,max
 40         int lc,rc,n;
 41     }d[100011],p[100011];
 42     int W;//W:当前维度
 43     bool c1(const N &a,const N &b)
 44     {
 45         return a.d[W]<b.d[W];
 46     }
 47     ll dis2(const N &a,const N &b)
 48     {
 49         return sqr(a.d[0]-b.d[0])+sqr(a.d[1]-b.d[1]);
 50     }
 51 #define LC (d[u].lc)
 52 #define RC (d[u].rc)
 53     void upd(int u)
 54     {
 55         for(int i=0;i<2;++i)
 56         {
 57             d[u].mx[i]=d[u].mn[i]=d[u].d[i];
 58             if(LC)
 59             {
 60                 cmax(d[u].mx[i],d[LC].mx[i]);
 61                 cmin(d[u].mn[i],d[LC].mn[i]);
 62             }
 63             if(RC)
 64             {
 65                 cmax(d[u].mx[i],d[RC].mx[i]);
 66                 cmin(d[u].mn[i],d[RC].mn[i]);
 67             }
 68         }
 69     }
 70     int build(int l,int r,int w)//w:当前维度
 71     {
 72         W=w;
 73         int u=(l+r)>>1;
 74         nth_element(p+l,p+u,p+r+1,c1);
 75         d[u]=p[u];
 76         if(l!=u)    LC=build(l,u-1,w^1);
 77         if(u!=r)    RC=build(u+1,r,w^1);
 78         upd(u);
 79         return u;
 80     }
 81     N x;
 82     const ll inf=0x3f3f3f3f3f3f3f3f;
 83     ll getmx(int u)//x到u的子树区域的最大可能距离的平方
 84     {
 85         ll ans=0;
 86         for(int i=0;i<2;++i)
 87             ans+=sqr(max(abs(x.d[i]-d[u].mn[i]),
 88                 abs(d[u].mx[i]-x.d[i])));
 89         return ans;
 90     }
 91     bool c2(const pli &a,const pli &b)
 92     {
 93         return (a.fi>b.fi) || (a.fi==b.fi && a.se<b.se);
 94     }
 95     pli ans[101];unsigned alen;
 96     inline void push(const pli &x)
 97     {
 98         ans[++alen]=x;
 99         push_heap(ans+1,ans+alen+1,c2);
100     }
101     inline void pop()
102     {
103         pop_heap(ans+1,ans+alen+1,c2);
104         --alen;
105     }
106 #define OK(x) ( (alen<K) || (ans[1].fi<=(x)) )
107     void querymx(int u)
108     {
109         pli tmp=(pli){dis2(d[u],x),d[u].n};
110         if(alen<K)
111             push(tmp);
112         else if(ans[1]<tmp)
113             pop(),push(tmp);
114         ll dl=-inf,dr=-inf;
115         if(LC)    dl=getmx(LC);
116         if(RC)    dr=getmx(RC);
117         if(dl>dr)
118         {
119             if(LC && OK(dl))    querymx(LC);
120             if(RC && OK(dr))    querymx(RC);
121         }
122         else
123         {
124             if(RC && OK(dr))    querymx(RC);
125             if(LC && OK(dl))    querymx(LC);
126         }
127     }
128     /*
129     void out()
130     {
131         printf("size:%d\n",int(ans.size()));
132             auto qqq=ans;
133             while(!qqq.empty())
134             {
135                 auto t=qqq.top();qqq.pop();
136                 printf("%lld\n",t);
137             }
138             {
139                 char s[233];
140                 scanf("%s",s);
141             }
142     }
143     */
144 }
145 int main()
146 {
147     int i;
148     scanf("%d",&n);
149     {
150         using namespace KK;
151         for(i=1;i<=n;++i)
152         {
153             scanf("%lld%lld",&p[i].d[0],&p[i].d[1]);
154             p[i].n=i;
155         }
156         rt=build(1,n,0);
157         scanf("%d",&m);
158         while(m--)
159         {
160             scanf("%lld%lld%u",&x.d[0],&x.d[1],&K);
161             alen=0;
162             querymx(rt);
163             //out();
164             printf("%d\n",ans[1].se);
165         }
166     }
167     return 0;
168 }

View Code

posted @ 2019-01-02 17:18 hehe_54321 阅读(256) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

哈

233

kdtree

公告