三角剖分就是能在O(n)或者O(nlogn)的时间內将平面上的点集用许多互不相交的三角形连起来，这些三角形有一个共同的性质：过该三角形的外接圆内不含其他任何点，称为空圆。

分治法求解三角剖分：

1. 按x-y排序并去重

2. 递归做 work(l,mid),work(mid+1,r);

3. 合并(l,mid)~(mid+1,r)

递归边界条件为当前处理的点数<=3，此时直接两两连边

合并第一步为找到两块点集的下公切线，这一步直接用凸包求就可以了。

然后就是把相邻两个点集连在一起。

假设当前刚刚连的一条边左端点为第ld号点，右端点为第rd号点。

我们遍历所有与当前ld,rd相邻的点，这些点中必定存在一个点id使得以三点(ld,rd,id)构成的圆不包含当前其它任何点。于是直接扫一遍，如果新加的点id'在圆(ld,rd,id)内则id=id'

然后就是连接这个点id，如果id原来在左边就连边(id,rd),删掉之前端点为ld，与(id,rd)有交的边，然后更新ld=id。id在右边同理。

删边只要用双向链表删就可以了。

由于每一步操作，边数都是O(n)的，所以整个复杂度是O(nlogn)的。

学习链接：

http://www.geom.uiuc.edu/~samuelp/del_project.html

http://blog.csdn.net/donkeylong/article/details/4909774

以BZOJ-3911为例，三角剖分出最小生成树之后，就是跟noip2013货车运输一样的了，lca倍增询问。

/**************************************************************
    Problem: 3911
    User: wjy1998
    Language: C++
    Result: Accepted
    Time:11336 ms
    Memory:93412 kb
****************************************************************/
 
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 100010
#define sqr(x) ((x)*(x))
#define max(a,b) (a>b?a:b)
int aa;char ch;int F(){
    while(ch=getchar(),ch<'0'||ch>'9');aa=ch-'0';
    while(ch=getchar(),ch>='0'&&ch<='9')aa=aa*10+ch-'0';return aa;
}
typedef double ld;
struct P{
    ld x,y;
#define PP const P&
    bool operator<(PP a)const {return x<a.x||x==a.x&&y<a.y;}
    P operator-(PP a)const {return (P){x-a.x,y-a.y};}
    ld operator&(PP a)const {return x*a.y-y*a.x;}
    ld operator|(PP a)const {return x*a.x+y*a.y;}
}p[N];
#define check(a,b,c) ((b-a)&(c-a))
ld dis2(PP a){return a.x*a.x+a.y*a.y;}
#define cross(a,b,c,d) (check(p[a],p[c],p[d])*check(p[b],p[c],p[d])<0&&check(p[c],p[a],p[b])*check(p[d],p[a],p[b])<0)
struct P3{
    ld x,y,z;
    bool operator<(const P3&a)const {return x<a.x||x==a.x&&y<a.y;}
    P3 operator-(const P3&a)const {return (P3){x-a.x,y-a.y,z-a.z};}
    ld operator|(const P3&a)const {return x*a.x+y*a.y+z*a.z;}
    P3 operator&(const P3&a)const {return (P3){y*a.z-z*a.y,z*a.x-x*a.z,x*a.y-y*a.x};}
}ori[N];
#define gp3(a) (P3){a.x,a.y,a.x*a.x+a.y*a.y}
bool incir(int a,int b,int c,int d){
    P3 aa=gp3(p[a]),bb=gp3(p[b]),cc=gp3(p[c]),dd=gp3(p[d]);
    if(check(p[a],p[b],p[c])<0)std::swap(bb,cc);
    return (check(aa,bb,cc)|(dd-aa))<0;
}
int n,m,i,j,et=1,la[N],ts,xx,yy,fa[N][20],tot,cnt,dep[N],l,r,q[N<<2],ufs[N];ld mx[N][20];
struct E{int to,l,r;}e[N<<5];
void add(int x,int y){
    e[++et]=(E){y,la[x]},e[la[x]].r=et,la[x]=et;
    e[++et]=(E){x,la[y]},e[la[y]].r=et,la[y]=et;
}
void del(int x){
    e[e[x].r].l=e[x].l,e[e[x].l].r=e[x].r,la[e[x^1].to]==x?la[e[x^1].to]=e[x].l:1;
}
void delaunay(int l,int r){
    if(r-l<=2){
        for(int i=l;i<r;i++)
        for(int j=i+1;j<=r;j++)add(i,j);
        return;
    }
    int i,j,mid=l+r>>1,ld=0,rd=0,id,op;
    delaunay(l,mid),delaunay(mid+1,r);
    for(tot=0,i=l;i<=r;q[++tot]=i++)
    while(tot>1&&check(p[q[tot-1]],p[q[tot]],p[i])<0)tot--;
    for(i=1;i<tot&&!ld;i++)if(q[i]<=mid&&mid<q[i+1])ld=q[i],rd=q[i+1];
    for(;add(ld,rd),1;){
        id=op=0;
        for(i=la[ld];i;i=e[i].l)if(check(p[ld],p[rd],p[e[i].to])>0)
        if(!id||incir(ld,rd,id,e[i].to))op=-1,id=e[i].to;
        for(i=la[rd];i;i=e[i].l)if(check(p[rd],p[ld],p[e[i].to])<0)
        if(!id||incir(ld,rd,id,e[i].to))op=1,id=e[i].to;
        if(op==0)break;
        if(op==-1){
            for(i=la[ld];i;i=e[i].l)
            if(cross(rd,id,ld,e[i].to))del(i),del(i^1),i=e[i].r;
            ld=id;
        }else{
            for(i=la[rd];i;i=e[i].l)
            if(cross(ld,id,rd,e[i].to))del(i),del(i^1),i=e[i].r;
            rd=id;
        }
    }
}
struct D{int x,y;ld v;}d[N<<3];
bool operator<(const D&i,const D&j){return i.v<j.v;}
int gf(int x){return ufs[x]==x?x:ufs[x]=gf(ufs[x]);}
struct G{int to;double v;int nxt;}g[N<<3];
#define addg(x,y,v) (g[++et]=(G){y,v,la[x]},la[x]=et)
ld query(int x,int y){
    int k,i;ld ans=0;
    if(dep[x]<dep[y])k=x,x=y,y=k;
    for(k=dep[x]-dep[y],i=0;k;k>>=1,i++)if(k&1)
    ans=max(ans,mx[x][i]),x=fa[x][i];
    if(x==y)return ans;
     
    for(i=0;fa[x][i]!=fa[y][i];i++);
    for(i--;~i;i--)if(fa[x][i]!=fa[y][i])
    ans=max(ans,max(mx[x][i],mx[y][i])),x=fa[x][i],y=fa[y][i];
    ans=max(ans,max(mx[x][0],mx[y][0]));return ans;
}
int main(){
    for(n=F(),i=1;i<=n;i++)xx=F(),yy=F(),p[i]=(P){xx,yy},ori[i]=(P3){xx,yy,i},ufs[i]=i;
    std::sort(p+1,p+1+n);std::sort(ori+1,ori+1+n);delaunay(1,n);
     
    for(i=1;i<=n;i++)
    for(j=la[i];j;j=e[j].l)xx=ori[i].z,yy=ori[e[j].to].z,
    d[++tot]=(D){xx,yy,dis2(p[i]-p[e[j].to])};
    std::sort(d+1,d+1+tot);
     
    memset(la,0,sizeof(la)),et=0;
    for(i=1;i<=tot&&cnt<n-1;i++)if(gf(d[i].x)!=gf(d[i].y))
    cnt++,ufs[ufs[d[i].x]]=ufs[d[i].y],
    addg(d[i].x,d[i].y,d[i].v),addg(d[i].y,d[i].x,d[i].v);
     
    for(q[l=r=1]=dep[1]=1;l<=r;l++)
    for(i=la[q[l]];i;i=g[i].nxt)if(!dep[g[i].to])
    for(dep[q[++r]=g[i].to]=dep[q[l]]+1,fa[g[i].to][j=0]=q[l],mx[g[i].to][0]=g[i].v;fa[g[i].to][j];j++)
    fa[g[i].to][j+1]=fa[fa[g[i].to][j]][j],mx[g[i].to][j+1]=max(mx[g[i].to][j],mx[fa[g[i].to][j]][j]);
     
    for(m=F();m--;printf("%.6lf\n",sqrt(query(F(),F()))));
}

一个Linux下可视化小程序