一、前人種樹
二、知識梳理
曼哈頓距離：給定二維平面上的N個點，在兩點之間連邊的代價。（即distance(P1,P2) = |x1－x2|+|y1－y2|）
曼哈頓距離最小生成樹問題求什么？求使所有點連通的最小代價。
最小生成樹的“環切”性質：在圖G = (V, E)中，如果存在一個環，那么把環上的最大邊e刪除后得到的圖G’ = (V, E- {e})的最小生成樹的邊權和與G相同。
三、難點剖析
【廢話定理神馬的，很難懂只要記住就是了】
朴素的算法可以用O(N2)的Prim，或者處理出所有邊做Kruskal，但在這里總邊數有O(N2)條，所以Kruskal的復雜度變成了O(N2logN)。
但是事實上，真正有用的邊遠沒有O(N2)條。我們考慮每個點會和其他一些什么樣的點連邊。
可以得出這樣一個結論：以一個點為原點建立直角坐標系，在每45度內只會向距離該點最近的一個點連邊。
證明結論：假設我們以點A為原點建系，考慮在y軸向右45度區域內的任意兩點B(x1,y1)和C(x2,y2)，不妨設|AB|≤|AC|（這里的距離為曼哈頓距離），如下圖：
|AB|=x1+y1，|AC|=x2+y2，|BC|=|x1-x2|+|y1-y2|。而由於B和C都在y軸向右45度的區域內，有y-x>0且x>0。下面我們分情況討論：
- x1>x2且y1>y2。這與|AB|≤|AC|矛盾；
- x1≤x2且y1>y2。此時|BC|=x2-x1+y1-y2，|AC|-|BC|=x2+y2-x2+x1-y1+y2=x1-y1+2*y2。由前面各種關系可得y1>y2>x2>x1。假設|AC|<|BC|即y1>2*y2+x1，那么|AB|=x1+y1>2*x1+2*y2，|AC|=x2+y2<2*y2<|AB|與前提矛盾，故|AC|≥|BC|；
- x1>x2且y1≤y2。與2同理；
- x1≤x2且y1≤y2。此時顯然有|AB|+|BC|=|AC|，即有|AC|>|BC|。
綜上有|AC|≥|BC|，也即在這個區域內只需選擇距離A最近的點向A連邊。
這種連邊方式可以保證邊數是O(N)的，那么如果能高效處理出這些邊，就可以用Kruskal在O(NlogN)的時間內解決問題。下面我們就考慮怎樣高效處理邊。
我們只需考慮在一塊區域內的點，其他區域內的點可以通過坐標變換“移動”到這個區域內。為了方便處理，我們考慮在y軸向右45度的區域。在某個點A(x0,y0)的這個區域內的點B(x1,y1)滿足x1≥x0且y1-x1>y0-x0。這里對於邊界我們只取一邊，但是操作中兩邊都取也無所謂。那么|AB|=y1-y0+x1-x0=(x1+y1)-(x0+y0)。在A的區域內距離A最近的點也即滿足條件的點中x+y最小的點。因此我們可以將所有點按x坐標排序，再按y-x離散，用線段樹或者樹狀數組維護大於當前點的y-x的最小的x+y對應的點。時間復雜度O(NlogN)。
至於坐標變換，一個比較好處理的方法是第一次直接做；第二次沿直線y=x翻轉，即交換x和y坐標；第三次沿直線x=0翻轉，即將x坐標取相反數；第四次再沿直線y=x翻轉。注意只需要做4次，因為邊是雙向的。
至此，整個問題就可以在O(NlogN)的復雜度內解決了。
【回到正題】
一個點把平面分成了8個部分：
由上面的廢話可知，我們只需要讓這個點與每個部分里距它最近的點連邊。
拿R1來說吧：
如圖，i的R1區域里距i最近的點是j。也就是說，其他點k都有：
xj + yj <= xk + yk
那么k將落在如下陰影部分：
顯然，邊(i,j), (j,k), (i,k)構成一個環<i,j,k>，而(i,k)一定是最長邊，可以被刪去。所以我們只連邊(i,j)。
為了避免重復加邊，我們只考慮R1~R4這4個區域。（總共加了4N條邊）
這4個區域的點(x,y)要滿足什么條件？
- 如果點(x,y)在R1，它要滿足：x ≥ xi ,y – x ≥ yi – xi（最近點的x + y最小）
- 如果點(x,y)在R2，它要滿足：y ≥ yi ,y – x ≤ yi – xi（最近點的x + y最小）
- 如果點(x,y)在R3，它要滿足：y ≤ yi ,y + x ≥ yi + xi（最近點的y – x最小）
- 如果點(x,y)在R4，它要滿足：x ≥ xi ,y + x ≤ yi – xi（最近點的y – x最小）
其中一個條件用排序，另一個條件用數據結構（這種方法很常用），在數據結構上詢問，找最近點。因為詢問總是前綴或后綴，所以可以用樹狀數組。
四、代碼模板
//離散化：
scanf("%d", &N);
for (int i=1; i<=N; ++i)
{
scanf("%d%d", &P[i].x, &P[i].y);
P[i].id = i;
P[i].d = P[i].y - P[i].x;
P[i].s = P[i].y + P[i].x;
}
//對x,y離散化
int totxy = 0;
for (int i=1; i<=N; ++i)
{
xy[totxy++] = P[i].x;
xy[totxy++] = P[i].y;
}
sort(xy, xy+totxy);
for (int i=1; i<=N; ++i)
{
P[i].idx = lower_bound(xy, xy+totxy, P[i].x) - xy + 1;
P[i].idy = lower_bound(xy, xy+totxy, P[i].y) - xy + 1;
}
//樹狀數組：
struct BIT
{
pii a[maxN * 2];
int N;
void Init(int _N)
{
N = _N;
for (int i=0; i<=N; ++i) a[i] = pii(oo, 0);
}
pii ask(int x)
{
return x == 0 ? pii(oo, 0) : min(a[x], ask(x - (x & (-x))));
}
void update(int x, const pii &v)
{
if (x > N) return ;
a[x] = min(a[x], v);
update(x + (x & (-x)), v);
}
pii ask_front(int x) {return ask(x);}
pii ask_back(int x) {return ask(N - x + 1);}
void update_front(int x, const pii &v) {update(x, v);}
void update_back(int x, const pii &v) {update(N - x + 1, v);}
} tree;
//構圖：
bool cmp1(const Tpoint &A, const Tpoint &B)
{
//return A.x < B.x || (A.x == B.x && A.y < B.y);
return (A.y - A.x > B.y - B.x || A.y - A.x == B.y - B.x && A.x > B.x);
}
bool cmp2(const Tpoint &A, const Tpoint &B)
{
//return A.x < B.x || (A.x == B.x && A.y > B.y);
return (A.y + A.x < B.y + B.x || A.y + A.x == B.y + B.x && A.x > B.x);
}
bool cmp3(const Tpoint &A, const Tpoint &B)
{
//return A.y < B.y || (A.y == B.y && A.x < B.x);
return A.y - A.x < B.y - B.x || A.y - A.x == B.y - B.x && A.y > B.y;
}
bool cmp4(const Tpoint &A, const Tpoint &B)
{
//return A.y < B.y || (A.y == B.y && A.x > B.x);
return A.s > B.s || A.s == B.s && A.y < B.y;;
}
bool cmpE(const E_arr &A, const E_arr &B) {return A.v < B.v;}
void Make_Graph()
{
#define Connect(i,j) E[++tot_E].Init(P[i].id,P[j].id,getdis(i,j))
int LL, RR;
tree.Init(2 * N);
sort(P+1, P+N+1, cmp1);
for (int i=1; i<=N; ++i)
{
pii tmp = tree.ask_back(P[i].idx);
if (tmp.first < oo) Connect(i, tmp.second);
tree.update_back(P[i].idx, pii(P[i].x + P[i].y, i));
}
sort(P+1, P+N+1, cmp2);
tree.Init(2 * N);
for (int i=1; i<=N; ++i)
{
pii tmp = tree.ask_back(P[i].idx);
if (tmp.first < oo) Connect(i, tmp.second);
tree.update_back(P[i].idx, pii(P[i].x - P[i].y, i));
}
sort(P+1, P+N+1, cmp3);
tree.Init(2 * N);
for (int i=1; i<=N; ++i)
{
pii tmp = tree.ask_back(P[i].idy);
if (tmp.first < oo) Connect(i, tmp.second);
tree.update_back(P[i].idy, pii(P[i].x + P[i].y, i));
}
sort(P+1, P+N+1, cmp4);
tree.Init(2 * N);
for (int i=1; i<=N; ++i)
{
pii tmp = tree.ask_front(P[i].idy);
if (tmp.first < oo) Connect(i, tmp.second);
tree.update_front(P[i].idy, pii(P[i].x - P[i].y, i));
}
}
五、沙場練兵
POJ 3241 Object Clustering 求曼哈頓距離最小生成樹上第k大的邊
//POJ3241; Object Clustering; Manhattan Distance MST
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#define N 100000
#define INFI 123456789
struct point
{
int x, y, n;
bool operator < (const point &p) const
{ return x == p.x ? y < p.y : x < p.x; }
}p[N + 1];
struct inedge
{
int a, b, w;
bool operator < (const inedge &x) const
{ return w < x.w; }
}e[N << 3 | 1];
struct BITnode
{
int w, p;
}arr[N + 1];
int n, k, tot = 0, f[N + 1], a[N + 1], *l[N + 1], ans;
template <typename T>
inline T abs(T x)
{ return x < (T)0 ? -x : x; }
int find(int x)
{ return x == f[x] ? x : f[x] = find(f[x]); }
inline bool cmp(int *a, int *b)
{ return *a < *b; }
inline int query(int x)
{
int r = INFI, p = -1;
for (; x <= n; x += x & -x)
if (arr[x].w < r) r = arr[x].w, p = arr[x].p;
return p;
}
inline void modify(int x, int w, int p)
{
for (; x > 0; x -= x & -x)
if (arr[x].w > w) arr[x].w = w, arr[x].p = p;
}
inline void addedge(int a, int b, int w)
{
++tot;
e[tot].a = a, e[tot].b = b, e[tot].w = w;
// printf("%d %d %d\n", a, b, w);
}
inline int dist(point &a, point &b)
{ return abs(a.x - b.x) + abs(a.y - b.y); }
int main()
{
//Initialize
scanf("%d%d", &n, &k);
for (int i = 1; i <= n; ++i)
{
scanf("%d%d", &p[i].x, &p[i].y);
p[i].n = i;
}
//Solve
for (int dir = 1; dir <= 4; ++dir)
{
//Coordinate transform - reflect by y=x and reflect by x=0
if (dir == 2 || dir == 4)
for (int i = 1; i <= n; ++i) p[i].x ^= p[i].y ^= p[i].x ^= p[i].y;
else if (dir == 3)
for (int i = 1; i <= n; ++i) p[i].x = -p[i].x;
//Sort points according to x-coordinate
std::sort(p + 1, p + n + 1);
//Discretize
for (int i = 1; i <= n; ++i) a[i] = p[i].y - p[i].x, l[i] = &a[i];
std::sort(l + 1, l + n + 1, cmp);
/*
int cnt = 1;
for (int i = 2; i <= n; ++i)
if (*l[i] != *l[i - 1]) *l[i - 1] = cnt++;
else *l[i - 1] = cnt;
*l[n] = cnt;
*/
for (int i = 1; i <= n; ++i) *l[i] = i;
//Initialize BIT
for (int i = 1; i <= n; ++i) arr[i].w = INFI, arr[i].p = -1;
//Find points and add edges
for (int i = n; i > 0; --i)
{
int pos = query(a[i]);
if (pos != -1)
addedge(p[i].n, p[pos].n, dist(p[i], p[pos]));
modify(a[i], p[i].x + p[i].y, i);
}
}
//Kruskal
std::sort(e + 1, e + tot + 1);
for (int i = 1; i <= n; ++i) f[i] = i;
for (int i = 1, ec = n; ec > k && i <= tot; ++i)
if (find(e[i].a) != find(e[i].b))
{
f[find(e[i].a)] = find(e[i].b);
if (--ec == k) ans = e[i].w;
}
printf("%d\n", ans);
return 0;
}