UVa12565 Lovely Magical Curves 题解

· · 题解

题意简述

给定两条 NURBS 曲线，求它们的交点。

题目思路

首先，需要注意的是种种奇奇怪怪东西的定义。

第一就是公式：

C(u)=\dfrac{\sum_{i=1}^nw_iN_{i,k}(u)P_i}{\sum_{i=1}^nw_iN_{i,k}(u)}

它的意思就是说，我们可以把如上的这些点和它的函数以及权重相乘，而后我们需要一些比较新奇的理解，即把所有的 P_i 的两个坐标加起来，而后除以下面的一大坨东西。

而后，我们观察一下这个多项式函数：

N_{i,k}(u)=\frac{u-t_i}{t_{i+k}-t_i}N_{i,k-1}(u)+\frac{t_{i+k+1}-u}{t_{i+k+1}-t_{i+1}}N_{i+1,k-1}(u) N_{i,0}(u)=[t_i\le u<t_{i+1}]

它的意思，已经被直白的表述在了公式之中，因此不解释。

因此，这是第一版代码中 NURBS 曲线的定义：

struct nurbs{
    int n; // numbers of control points
    int k; // degree
    int m; // number of knots
    point P[25]; // control points
    double w[25]; // weight of points
    double t[25]; // knot vector
    double N(int i,int k,double u){ // Function N
        if(k==0) return (t[i]<=u&&u<t[i+1]?1.0:0.0);
        double co0,co1;
        if(fabs(t[i+k]-t[i])<EPS||fabs(u-t[i])<EPS) co0=0;
        else co0=(u-t[i])/(t[i+k]-t[i]);
        if(fabs(t[i+k+1]-u)<EPS||fabs(t[i+k+1]-t[i+1])<EPS) co1=0;
        else co1=(t[i+k+1]-u)/(t[i+k+1]-t[i+1]);
        return co0*N(i,k-1,u)+co1*N(i+1,k-1,u);
    }
    point C(double u){ // Function C (for a single curve)
        point num=point(0,0);
        double dem=0;
        for(int i=1;i<=n;i++){
            num=num+w[i]*N(i,k,u)*P[i];
            dem+=w[i]*N(i,k,u);
        }
        return num/dem;
    }
    void clear(){
        n=k=m=0;
        for(int i=0;i<25;i++) P[i].x=P[i].y=w[i]=t[i]=0;
    }
}curv[3];

当然，多测不清空的悲剧还是要尽量避免，因此有一个 clear() 在这里。

另外，题目加粗强调了 0/0=0，但是我要告诉你的是，x/y 中只要 x=0 或者 y=0 那么这个式子就等于 0，不等于任何别的数！

于是，我们当然可以把第一版程序中点结构体拿出来：

const double EPS=8e-5;
struct point{
    double x,y;
    point(double cx=0,double cy=0):x(cx),y(cy){}
};
point operator*(double a,point p){
    return point(p.x*a,p.y*a);
}
point operator+(point a,point b){
    return point(a.x+b.x,a.y+b.y);
}
point operator/(point a,double b){
    double x=a.x/b,y=a.y/b;
    if(fabs(a.x)<EPS||fabs(b)<EPS) x=0;
    if(fabs(a.y)<EPS||fabs(b)<EPS) y=0;
    return point(x,y);
}

好的，我们就可以通过暴力枚举 [0,1) 里的非常多的点来找到交点，如果要找到交点，可以直接用 set 里的 intersection。

但是，如果这样的话，你大概率会那道一个 WA 或者 RE 或者 TLE。

那么我们应该怎么办？

众所周知，我们可以通过把曲线近似为很多条线段（亦即折线）来达到相同的结果，因此，我们可以设一个非常大的正整数 R，而后设 s=1/R，随后对于 k=0\cdots(R-1)，或者说只要 ks\neq 1，我们就能找到点 C(ks)，并且把相邻的两个点用线段连接起来。

这是不需要多说的，因为下面就要开始让人心肺骤停了。

首先的问题是：如何求得两条线段的交点，而它的简化方式，就是求出直线的交点。我们可以设两条直线分别为 P+t\overrightarrow{v} 和 Q+t\overrightarrow{w}，而后设 \overrightarrow{u}=P-Q，然后通过解方程（过程略），我们就能得出交点的位置在

P+\left(\frac{\overrightarrow{w}\times\overrightarrow{u}}{\overrightarrow{v}\times\overrightarrow{w}}\right)\overrightarrow{v}

这样我们就可以给出这样的一份代码求两条直线的交点：

point getintersection(point p,vecto v,point q,vecto w){
    vecto u=p-q;
    double t=cross(w,u)/cross(v,w);
    return p+v*t;
}

那么，我们现在研究线段的相交，它的充分必要条件就是，每条线段的两个端点都在另一条线段的两侧（即叉积的符号不同）。因此我们这样写出：

bool isintersected(point a1,point a2,point b1,point b2){
    double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}

至少到这里为止的计算几何是不需要图片作为补充的，然而下面我们就需要了。

当然，这里需要先给出通过计算求近似折线的算法（也不知道第几版代码了）：

vector<pair<point,point>> aprlin[3];
void approx(int num){
    double st=1/ROX;
    vector<point> aprdot;
    for(double i=0;i*st<1;i++) aprdot.push_back(curv[num].C(st*i));
    for(int i=1;i<aprdot.size();i++) aprlin[num].push_back(make_pair(aprdot[i-1],aprdot[i]));
}

而后，我们就要求折线的全部交点了——这就是我们的任务。

首先上场的是 \mathrm O(n^2) 的暴力判断，这哥们把所有的线段两两互相判断，代码如下：

vector<point> ans,shans;
for(auto it:aprlin[1]) for(auto jt:aprlin[2]){
    if(isintersected(it.first,it.second,jt.first,jt.second)){
        point pi=getintersection(it.first,it.second-it.first,jt.first,jt.second-jt.first);
        ans.push_back(pi);
    }
}

这里的 shans 是为了精度和去重用的，暂且不需要它。

然而不幸的是，如果你使用如上的代码，你会 TLE。怎么办？我们需要更快的算法求任意两条线段的交点（当然，如果两条线段不分属两个集合，我们可以特判）。

在多方查找后，我在这里找到了一个好的算法，在这里叙述一下。

这个算法是扫描线的另一种奇妙的实现。算法要求我们首先开一个链表，用以保存所有的线段端点以及找到的交点，按照 y 坐标从大到小，x 坐标从小到大的递增顺序建立链表；以及一棵二叉查找树，记录与扫描线相交的线段的编号，按照所在线段的上端点的 x 坐标递增顺序存储。

在这里我采用的实现方式是指针链表和 set 实现的二叉查找树（Treap,Splay 太烦了不想写），当然用自定义的 cmp 来排列元素，因此代码是这样写的：

struct cmp{
    bool operator()(int a,int b){
        point at,bt;
        if(aprlin[a].first.p1.y<aprlin[a].first.p2.y) at=aprlin[a].first.p2;
        else at=aprlin[a].first.p1;
        if(aprlin[b].first.p1.y<aprlin[b].first.p2.y) bt=aprlin[b].first.p2;
        else bt=aprlin[b].first.p1;
        return at.x<bt.x;
    }
};
set<int,cmp> bst;

我们还得探究一下，如何保留点所在线段的信息（交点需要两个！），点的类型（分为上点、下点和交点，字面意思，与扫描线平行的强行区分），以及它们之间的比较，因此我又定义了个结构体 endpoint 用来记录它（下面的 node 是链表节点，不论）：

const int BOTTOM=0,TOP=1,INTERSECT=-1;
struct endpoint{
    point p;
    int seg,ges,st;
    endpoint(point p=point(0,0),int seg=0,int ges=0,int st=0):p(p),seg(seg),ges(ges),st(st){}
};
struct node{
    endpoint ep;
    node *next;
    node(endpoint ep=endpoint(),node *next=nullptr):ep(ep),next(next){}
}*head;

那么初始化这条链表的代码如下（一开始我们存了所有线段的端点）：

vector<endpoint> edps;
int n=aprlin.size();
for(int i=0;i<n;i++){
    if(aprlin[i].first.p1.y>aprlin[i].first.p2.y){
        edps.push_back(endpoint(aprlin[i].first.p1,i,0,TOP));
        edps.push_back(endpoint(aprlin[i].first.p2,i,0,BOTTOM));
    }else{
        edps.push_back(endpoint(aprlin[i].first.p1,i,0,BOTTOM));
        edps.push_back(endpoint(aprlin[i].first.p2,i,0,TOP));
    }
}
sort(edps.begin(),edps.end(),[](endpoint a,endpoint b)->bool {return a.p.y!=b.p.y?a.p.y>b.p.y:a.p.x<b.p.x;});
head=new node;
node *cur=head,*co=head;
for(int i=0;i<2*n;i++){
    cur->ep=edps[i];
    if(cur->next==nullptr) cur->next=new node;
    co=cur,cur=cur->next;
}
delete cur,co->next=nullptr;
cur=head;

下面的过程非常之恐怖，因此略去代码，专门讲解几何部分。

扫描线按照链表行进，可能会遇到如下的三种情况：

  1. 扫描线碰上了某条线段的上点。这时我们设该线段编号为 \mathfrak{it}，左边线段为 \mathfrak{lt}，右边线段为 \mathfrak{rt}，左边线段 \mathfrak{lt} 和右边线段 \mathfrak{rt} 指在二叉查找树中与 \mathfrak{it} 相邻的两条线段，那么如图所示： 我们就可以把 \mathfrak{it} 塞进二叉查找树，而后判断 \mathfrak{it} 和 \mathfrak{lt} 以及 \mathfrak{rt} 是否相交。相交的节点一定在扫描线的下方（因为如果它在上方，说明我们已经搜索到了它并且已经用它插入过线段了），而后我们就把交点插入链表。
  2. 扫描线碰到了某条线段的下点，设的同上，我们就有 这时我们把 \mathfrak{lt} 和 \mathfrak{rt} 判断一下是否相交，如果相交就把交点插入链表，而后把 \mathfrak{it} 从扫描线中删除。当然，交点还是在扫描线的下方。
  3. 扫描线碰到了两条线段的交点，这时我们设这两条相交的线段编号分别为 \mathfrak{it} 和 \mathfrak{jt}（这里 \mathfrak{it} 在二叉查找树中相对于 \mathfrak{jt} 是在左边的！），并且设 \mathfrak{it} 的左相邻线段为 \mathfrak{lt}，\mathfrak{jt} 的右相邻线段为 \mathfrak{rt}，如图所示： 我们只需要判定 \mathfrak{it} 和 \mathfrak{rt} 是否相交以及 \mathfrak{lt} 和 \mathfrak{jt} 是否相交，而其他的交点我们已经找过了。

当然，这里需要注意的是，由于采用的是 set 来实现二叉查找树，所以它的 iterator 中寻找该元素要用 lower_bound，它的左邻居需要先判断是否是 begin() 而后左移，右邻居需要先右移而后判断是不是 end()，因为 end() 指向最后一个元素的后一个元素。最后释放整个链表。

本题完整版代码如下（为防止抄袭，R 和 \epsilon 的数值更改）：

#include <iostream>
#include <cmath>
#include <iomanip>
#include <set>
#include <algorithm>
#include <vector>
#include <map>
#define fs(X) fixed<<setprecision(X)
using namespace std;
const double EPS=0,ROX=0;
int dcmp(double x){
    if(fabs(x)<EPS) return 0;
    else return x<0?-1:1;
}
struct point{
    double x,y;
    point(double cx=0,double cy=0):x(cx),y(cy){}
    point operator+(point b){
        return point(x+b.x,y+b.y);
    }
    point operator-(point b){
        return point(x-b.x,y-b.y);
    }
    point operator*(double cof){
        return point(x*cof,y*cof);
    }
    point operator/(double b){
        double xx=x/b,yy=y/b;
        if(dcmp(x)==0||dcmp(b)==0) xx=0;
        if(dcmp(y)==0||dcmp(b)==0) yy=0;
        return point(xx,yy);
    }
};
typedef point vecto;
double cross(vecto a,vecto b){
    return a.x*b.y-a.y*b.x;
}
struct segment{
    point p1,p2;
    segment(point p1,point p2):p1(p1),p2(p2){}
};
bool isintersected(point a1,point a2,point b1,point b2){
    double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
point getintersection(point p,vecto v,point q,vecto w){
    vecto u=p-q;
    double t=cross(w,u)/cross(v,w);
    return p+v*t;
}
struct nurbs{
    int n,k,m;
    point P[25];
    double w[25],t[25];
    double N(int i,int k,double u){
        if(k==0) return (t[i]<=u&&u<t[i+1]?1.0:0.0);
        double co0,co1;
        if(dcmp(t[i+k]-t[i])==0||dcmp(u-t[i])==0) co0=0;
        else co0=(u-t[i])/(t[i+k]-t[i]);
        if(dcmp(t[i+k+1]-u)==0||dcmp(t[i+k+1]-t[i+1])==0) co1=0;
        else co1=(t[i+k+1]-u)/(t[i+k+1]-t[i+1]);
        return co0*N(i,k-1,u)+co1*N(i+1,k-1,u);
    }
    point C(double u){
        point num=point(0,0);
        double dem=0;
        for(int i=1;i<=n;i++){
            double coef=w[i]*N(i,k,u);
            num=num+P[i]*coef;
            dem+=coef;
        }
        return num/dem;
    }
    void clear(){
        n=k=m=0;
        for(int i=0;i<25;i++) P[i].x=P[i].y=w[i]=t[i]=0;
    }
}curv[3];
vector<pair<segment,int>> aprlin;
const int BOTTOM=0,TOP=1,INTERSECT=-1;
struct endpoint{
    point p;
    int seg,ges,st;
    endpoint(point p=point(0,0),int seg=0,int ges=0,int st=0):p(p),seg(seg),ges(ges),st(st){}
};
struct node{
    endpoint ep;
    node *next;
    node(endpoint ep=endpoint(),node *next=nullptr):ep(ep),next(next){}
}*head;
void approx(int num){
    double st=1/ROX;
    vector<point> aprdot;
    for(double i=0;i*st<1;i++) aprdot.push_back(curv[num].C(st*i));
    for(int i=1;i<aprdot.size();i++) aprlin.push_back(make_pair(segment(aprdot[i-1],aprdot[i]),num));
}
void intersections(vector<point> &ans){
    vector<endpoint> edps;
    int n=aprlin.size();
    for(int i=0;i<n;i++){
        if(aprlin[i].first.p1.y>aprlin[i].first.p2.y){
            edps.push_back(endpoint(aprlin[i].first.p1,i,0,TOP));
            edps.push_back(endpoint(aprlin[i].first.p2,i,0,BOTTOM));
        }else{
            edps.push_back(endpoint(aprlin[i].first.p1,i,0,BOTTOM));
            edps.push_back(endpoint(aprlin[i].first.p2,i,0,TOP));
        }
    }
    sort(edps.begin(),edps.end(),[](endpoint a,endpoint b)->bool {return a.p.y!=b.p.y?a.p.y>b.p.y:a.p.x<b.p.x;});
    head=new node;
    node *cur=head,*co=head;
    for(int i=0;i<2*n;i++){
        cur->ep=edps[i];
        if(cur->next==nullptr) cur->next=new node;
        co=cur,cur=cur->next;
    }
    delete cur,co->next=nullptr;
    cur=head;
    struct cmp{
        bool operator()(int a,int b){
            point at,bt;
            if(aprlin[a].first.p1.y<aprlin[a].first.p2.y) at=aprlin[a].first.p2;
            else at=aprlin[a].first.p1;
            if(aprlin[b].first.p1.y<aprlin[b].first.p2.y) bt=aprlin[b].first.p2;
            else bt=aprlin[b].first.p1;
            return at.x<bt.x;
        }
    };
    set<int,cmp> bst;
    while(cur!=nullptr){
        if(cur->ep.st==TOP){
            bst.insert(cur->ep.seg);
            auto it=bst.lower_bound(cur->ep.seg),lt=it,rt=it;
            rt++;
            if(rt!=bst.end()){
                if(isintersected(aprlin[*it].first.p1,aprlin[*it].first.p2,aprlin[*rt].first.p1,aprlin[*rt].first.p2)){
                    point pi=getintersection(aprlin[*it].first.p1,aprlin[*it].first.p2-aprlin[*it].first.p1,aprlin[*rt].first.p1,aprlin[*rt].first.p2-aprlin[*rt].first.p1);
                    node *nx=cur->next,*add=new node;
                    add->ep=endpoint(pi,*it,*rt,INTERSECT);
                    add->next=nx,cur->next=add;
                }
            }
            if(lt!=bst.begin()){
                lt--;
                if(isintersected(aprlin[*it].first.p1,aprlin[*it].first.p2,aprlin[*lt].first.p1,aprlin[*lt].first.p2)){
                    point pi=getintersection(aprlin[*it].first.p1,aprlin[*it].first.p2-aprlin[*it].first.p1,aprlin[*lt].first.p1,aprlin[*lt].first.p2-aprlin[*lt].first.p1);
                    node *nx=cur->next,*add=new node;
                    add->ep=endpoint(pi,*lt,*it,INTERSECT);
                    add->next=nx,cur->next=add;
                }
            }
        }else if(cur->ep.st==BOTTOM){
            auto it=bst.lower_bound(cur->ep.seg),lt=it,rt=it;
            ++rt;
            if(lt!=bst.begin()&&rt!=bst.end()){
                lt--;
                if(isintersected(aprlin[*rt].first.p1,aprlin[*rt].first.p2,aprlin[*lt].first.p1,aprlin[*lt].first.p2)){
                    point pi=getintersection(aprlin[*rt].first.p1,aprlin[*rt].first.p2-aprlin[*rt].first.p1,aprlin[*lt].first.p1,aprlin[*lt].first.p2-aprlin[*lt].first.p1);
                    node *nx=cur->next,*add=new node;
                    add->ep=endpoint(pi,*lt,*rt,INTERSECT);
                    add->next=nx,cur->next=add;
                }
            }
            bst.erase(cur->ep.seg);
        }else if(cur->ep.st==INTERSECT){
            if(aprlin[cur->ep.seg].second+aprlin[cur->ep.ges].second==3) ans.push_back(cur->ep.p);
            auto it=bst.lower_bound(cur->ep.seg),jt=bst.lower_bound(cur->ep.ges);
            auto lt=it,rt=jt;
            ++rt;
            if(rt!=bst.end()){
                if(isintersected(aprlin[*jt].first.p1,aprlin[*jt].first.p2,aprlin[*rt].first.p1,aprlin[*rt].first.p2)){
                    point pi=getintersection(aprlin[*jt].first.p1,aprlin[*jt].first.p2-aprlin[*jt].first.p1,aprlin[*rt].first.p1,aprlin[*rt].first.p2-aprlin[*rt].first.p1);
                    node *nx=cur->next,*add=new node;
                    add->ep=endpoint(pi,*jt,*rt,INTERSECT);
                    add->next=nx,cur->next=add;
                }
            }
            if(lt!=bst.begin()){
                lt--;
                if(isintersected(aprlin[*lt].first.p1,aprlin[*lt].first.p2,aprlin[*it].first.p1,aprlin[*it].first.p2)){
                    point pi=getintersection(aprlin[*lt].first.p1,aprlin[*lt].first.p2-aprlin[*lt].first.p1,aprlin[*it].first.p1,aprlin[*it].first.p2-aprlin[*it].first.p1);
                    node *nx=cur->next,*add=new node;
                    add->ep=endpoint(pi,*lt,*it,INTERSECT);
                    add->next=nx,cur->next=add;
                }
            }
        }
        cur=cur->next;
    }
    cur=head;
    do{
        node *co=cur->next;
        delete cur;
        cur=co;
    }while(cur!=nullptr);
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    int T;
    cin>>T;
    for(int tid=1;tid<=T;tid++){
        if(tid>1) cout<<endl;
        aprlin.clear();
        for(int i=1;i<=2;i++) curv[i].clear();
        for(int i=1;i<=2;i++){
            cin>>curv[i].n>>curv[i].m,curv[i].k=curv[i].m-curv[i].n-1;
            for(int j=1;j<=curv[i].n;j++) cin>>curv[i].P[j].x>>curv[i].P[j].y>>curv[i].w[j];
            for(int j=1;j<=curv[i].m;j++) cin>>curv[i].t[j];
            approx(i);
        }
        vector<point> ans;
        intersections(ans);
        cout<<"Case "<<tid<<": ";
        if(ans.empty()){
            cout<<0<<endl;
            continue;
        }
        sort(ans.begin(),ans.end(),[](point a,point b)->bool {return a.x!=b.x?a.x<b.x:a.y<b.y;});
        cout<<ans.size()<<endl;
        for(auto it:ans) cout<<"("<<fs(3)<<it.x<<", "<<fs(3)<<it.y<<")"<<endl;
    }
    return 0;
}

真的好长。当然，对于交点，我们可以在特判后（采用在双向广搜中使用的技巧，编号和为 3）把答案压入 vector 而后排序。

写这道题是为了计算几何入门的，结果一不小心就变成了这样：

可能前两个还是较为有趣的，但是第三个就有些恐怖了——自 2012 年出的题目，竟然到现在 A 掉它的不超过 5 个人！

好吧，最后贡献一组 Hack（从比赛的官方文件中薅下来的）。

EOF