POJ 4048 Chinese Repeating Crossbow（点积加叉积和极角排序的两种做法）

题意：平面上给你若干条线段和一个点o，从点o为端点，确定一条射线，最多能穿过多少条线段（刚好穿过线段的端点也算）。

易知，我们可以把点o向每条线段的每个端点划一条射线，然后取某一条穿过最多线段的射线就可以了。可以用两种做法。

1.点积加叉积

思路是：通过点o和某线段的一个端点得到一条射线（就是将连接两点得到的线段按向量的方向延长，得到近似无限长的射线），然后判断是否与线段相交即可。这里的相交是非规范相交。

这部分的大部份内容参考了黑书的内容。

先从线段的规范相交说起。规范相交<==>两条线段恰有一个不是端点的公共点<==>每条线段两个端点都在另一条线段所在直线的异端（充要条件）。

  

如上图，要判断线段p11p12和线段p21p22是否相交，可以转化成，判断点p21和p22是否在向量(p11p12)的异侧。如果我们能够判断一个点是在向量的左侧和右侧，那么就可以解决这个问题，注意，这里说的向量的左边还是右边不是相对于绝对的坐标系的，而是相对于它本身的方向的。如下图所示：

要判断点p22是在向量(p11p12)的哪边也很简单，只要添加一条辅助向量(p11p22)即可，现在我们只要判断对于公共起点p11，向量(p11p22)相对于向量(p11p12)是顺时针还是逆时针方向，或者说，从向量(p11p12)到向量(p11p22)是左手螺旋还是右手螺旋方向。

方法是叉积，向量a叉乘向量b的结果是：xa*yb-xb*ya。当向量a逆时针旋转到向量b的角度小于180度（即向量b在向量a的逆时针方向，这时向量a到向量b是右手螺旋，成右手系），叉积的结果是正的，当超过180度（即向量b在向量a的顺时针方向，这时候向量a到向量b是左手螺旋，成左手系），结果是负的，而当向量a和向量b共线时候，不管是正向还是反向，其结果总为零。

容易证明叉积的结果的绝对值等于以a和b的边的平行四边形面积。如果保留叉积的结果的正负号，可以得到有向面积。即当向量b在向量a的逆时针方向时，叉积的结果是正的为正面积，反之叉积的结果是负的为负面积。如果定义角度Ɵ为从向量a逆时针旋转到向量b的角度，可以得到叉积的几何形式：|a|*|b|*sinƟ。

那么到现在就大致可以判断两条线段是否规范相交了，主要的代码是：

struct Point {double x,y;};
double det(double x1,double y1,double x2,double y2)
{
	return x1*y2-x2*y1;
}
double cross(Point a,Point b,Point c)
{
	return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}

但是计算机在处理实数的时候是用浮点数表示的，精度受到了限制，特别是结过相当的运算（如乘除）后，误差会累积起来，这时候，判断“<0"就容易出错。所以要改写"<0"的判断，允许一定的误差，即可以认识在0附近的某个小范围内（如10^-6）都可以看作0。由此可以写一个三出口的判断函数。

#define eps (1e-6)
int dcmp(double d)
{
	if(fabs(d)<eps) return 0;
	return d>0?1:-1;
}

判断两线段是否相交的关键代码段如下：

int segCross(Point a,Point b,Point c,Point d)
{
	return (dcmp(cross(a,c,d))^dcmp(cross(b,c,d))==-2) &&
		(dcmp(cross(c,a,b))^dcmp(cross(d,a,b))==-2);
}

非规范相交在规范相交的基础上允许线段的端点在另一条线段上（内部或者端点重合）。那么我们只要判断是否至少有一端点A在另一条线段BC上就可以了。而这无非是以下两点：1.A在BC所在的直线上（即向量(AB)叉乘(AC)的结果为0）。2.A在BC范围内。

对于2可以借助点积，如果说叉积是判断的是左右，那么点积判断的就是前后，或者说，点积是考察两个向量在方向上的一致性。如果已经点A是在线段BC上，那么通过向量(CA)点乘(CB)的结果就可以判断，如果结果等于0，说明C与线段AB的端点A或B重合，如果结果大于0说明在直线外，如果小于0说明在直线内。这样就可以判断出非规范相交了，下面是黑书上的判断非规范相交的代码：

double dot(double x1,double y1,double x2,double y2)//点积
{
    return x1*x2+y1*y2;
}
double dotdet(Point a,Point b,Point c)
{
    return dot(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
int betweenCmp(Point a,Point b,Point c)//通过点积判断点a是否在线段bc上
{
    return dblcmp(dotdet(a,b,c));
}
double det(double x1,double y1,double x2,double y2)//叉积
{
    return x1*y2-x2*y1;
}
double cross(Point a,Point b,Point c)
{
    return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
int segCross(Point a,Point b,Point c,Point d,Point &p)
{
    double s1,s2,s3,s4;
    int d1,d2,d3,d4;
    d1=dblcmp(s1=cross(a,b,c));
    d2=dblcmp(s2=cross(a,b,d));
    d3=dblcmp(s3=cross(c,d,a));
    d4=dblcmp(s4=cross(c,d,b));
    if((d1^d2)==-2&&(d3^d4)==-2)//如果是规范相交则通过定比分点求其交点
    {
        p.x=(c.x*s2-d.x*s1)/(s2-s1);
        p.y=(c.y*s2-d.y*s1)/(s2-s1);
        return 1;
    }
    if(d1==0&&betweenCmp(c,a,b)<=0||
        d2==0&&betweenCmp(d,a,b)<=0||
        d3==0&&betweenCmp(a,c,d)<=0||
        d4==0&&betweenCmp(b,c,d)<=0)
        return 2;
    return 0;
}

到这里，就可以解决这题了，下面是我的代码。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define eps (1e-6)
int dblcmp(double a)
{
    if(fabs(a)<eps) return 0;
    return a>0?1:-1;
}
struct Point
{
    double x,y;
    Point(){}
    Point(double a,double b){x=a;y=b;}
};
struct Line
{
    Point x,y;
    Line(){}
    Line(Point a,Point b){x=a;y=b;}
}data[2000];
double dot(double x1,double y1,double x2,double y2)
{
    return x1*x2+y1*y2;
}
double dotdet(Point a,Point b,Point c)
{
    return dot(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
int betweenCmp(Point a,Point b,Point c)
{
    return dblcmp(dotdet(a,b,c));
}
double det(double x1,double y1,double x2,double y2)
{
    return x1*y2-x2*y1;
}
double cross(Point a,Point b,Point c)
{
    return det(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
int segCross(Point a,Point b,Point c,Point d,Point &p)
{
    double s1,s2,s3,s4;
    int d1,d2,d3,d4;
    d1=dblcmp(s1=cross(a,b,c));
    d2=dblcmp(s2=cross(a,b,d));
    d3=dblcmp(s3=cross(c,d,a));
    d4=dblcmp(s4=cross(c,d,b));
    if((d1^d2)==-2&&(d3^d4)==-2)
    {
        p.x=(c.x*s2-d.x*s1)/(s2-s1);
        p.y=(c.y*s2-d.y*s1)/(s2-s1);
        return 1;
    }
    if(d1==0&&betweenCmp(c,a,b)<=0||
        d2==0&&betweenCmp(d,a,b)<=0||
        d3==0&&betweenCmp(a,c,d)<=0||
        d4==0&&betweenCmp(b,c,d)<=0)
        return 2;
    return 0;
}
bool judge(Point a,Point b,Line line,Point &p)
{
    Point tmp;
    tmp.x=b.x-a.x;    tmp.y=b.y-a.y;
    tmp.x=tmp.x*10000+a.x;
    tmp.y=tmp.y*10000+a.y;
    return segCross(a,tmp,line.x,line.y,p);
}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        int n;
        Point o;
        scanf("%d",&n);
        for(int i=0;i<n;i++)
        {
            double a,b,c,d;
            scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
            data[i]=Line(Point(a,b),Point(c,d));
        }
        scanf("%lf%lf",&o.x,&o.y);
        int res=0;
        for(int i=0;i<n;i++)
        {
            Point p;
            int cnt1=0,cnt2=0;
            for(int j=0;j<n;j++)
            {
                if(judge(o,data[i].x,data[j],p)) cnt1++;
                if(judge(o,data[i].y,data[j],p)) cnt2++;
            }
            res=max(res,max(cnt1,cnt2));
        }
        printf("%d\n",res);
    }
    return 0;
}

2.极角排序。

引用这里的一段话http://www.cnblogs.com/karlvin/archive/2012/05/23/poj4048.html。

我们来看另外一个情景：某个电影院门口有个光门，如果有人经过那个门口时就会立即触发光门记录一个时间值，若是离开电影院标记时间属性为out，若是进入电影就标记时间属性为in。求电影院里最大人流量。
     其实把光门记录的时间按从小到大排序，时间值相同时out属性的时间值排在in后面，我们将所有时间值属性按照时间值从小到大遍历一遍，in加1，out减1，中间出现的最大数字即为所求。
      所以本题最重要的是要将线段的端点进行排序。如果排序（即极角排序）解决了，按照上面的方法求解就行。

为处理方便可以将o点移到原点，当然直接调用atan2可以直接得到极角方便地处理，但是会有些精度的误差，对于这题影响不大。这里要做的也是完全用整数来实现的极角排序。我们规定如果向量a在向量b的顺时针方向，那么向量a的极角就比较小——即借助叉积，同时，当y等于0时，x大于0的向量最小。有一种特殊的情况是，如果有线段跨越了x的正半轴，这时候将这个线段拆分成两段，我的做法是将与x正半轴相交的那一点拆成两点，一点还是在x轴上，另一点向下移一些，将其y坐标赋值成-1。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
#include <algorithm>
using namespace std;
#define INF (1<<30)
const int N=2005;
struct Point
{
    int x,y,valu;
    Point(){}
    Point(int a,int b,int c){x=a;y=b;valu=c;}
    Point operator-(const Point &b){return Point(x-b.x,y-b.y,0);}
    int operator*(const Point &b){return x*b.y-y*b.x;}
    friend bool operator<( Point p1,Point p2)
    {
        if(p1.y==0&&p1.x>0)
        {
            if(p2.y==0&&p2.x>0) return p1.valu>p2.valu;
            return true;
        }
        else if(p2.y==0&&p2.x>0) return false;
        if(p1.y*p2.y<0) return p1.y>p2.y;
        int tmp=p1*p2;
        return tmp>0||(tmp==0&&p1.valu>p2.valu);
    }
};
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        vector<Point> p1,p2,polar;
        int n;
        scanf("%d",&n);
        for(int i=0;i<n;i++)
        {
            int a,b,c,d;
            scanf("%d%d%d%d",&a,&b,&c,&d);
            p1.push_back(Point(a,b,0));
            p2.push_back(Point(c,d,0));
        }
        Point o;
        scanf("%d%d",&o.x,&o.y);
        int iadd=0;
        for(int i=0;i<n;i++)
        {
            p1[i]=p1[i]-o;
            p2[i]=p2[i]-o;
            if(p2[i]*p1[i]<0) swap(p1[i],p2[i]);
            bool flag1=(p1[i].x==0&&p1[i].y==0);
            bool flag2=(p2[i].x==0&&p2[i].y==0);
            bool flag3=(p1[i].x*p2[i].x<0)||(p1[i].y*p2[i].y<0);
            if(flag1||flag2||(p1[i]*p2[i]==0&&flag3)) iadd++;
            else 
            {
                if(p1[i]<p2[i])
                {
                    polar.push_back(Point(p1[i].x,p1[i].y,-1));
                    polar.push_back(Point(1,0,1));
                    polar.push_back(Point(p2[i].x,p2[i].y,1));
                    polar.push_back(Point(INF,-1,-1));
                }
                else 
                {
                    polar.push_back(Point(p1[i].x,p1[i].y,-1));
                    polar.push_back(Point(p2[i].x,p2[i].y,1));
                }
            }
        }
        sort(polar.begin(),polar.end());
        int tmp=0,imax=0;
        for(int i=0;i<(int)polar.size();i++)
        {
            tmp+=polar[i].valu;
            imax=max(imax,tmp);
        }
        printf("%d\n",imax+iadd);
    }
    return 0;
}