关于半平面交问题

原帖地址

真是抱歉,拖了这么久才回复,确实有点忙啊。

先说点题外话。因为我是学医的(我专业是临床医学而不是计算机科学),在大三(有可能是大四,具体由于时间比较久,记不清楚了)的时候,老师开始教授《病理生理学》,这门课的目的是了解人体在疾病状态下组织和细胞所呈现的形态。奇怪的是,老师在最初的几周并不给我们看疾病状态的组织和细胞的切片,而是反复给我们看正常情形下人体组织的细胞切片(比如肝细胞、上皮细胞等)。我们很纳闷,不禁向老师提出疑问。老师解释说,给您们看正常情形下的切片是为了让您们熟悉正常情况下人体的细胞是什么样子,当您们对正常情况的细胞是什么样子了然如胸时,看疾病状态的切片就易如反掌了。

我想,如老师说的那样,要对半平面交算法的实现有比较深刻的理解,必须先看看标准情况下,它是如何实现的。

朱泽园关于求半平面交的增量排序算法是2006年写的。说实话,他的论文在叙述上很糟糕,有些地方晦涩难懂,给出的图示和文章内容也有不相符合的地方,最后也没有给出参考代码,使得很多初次看论文的人都会产生疑惑。在论文中,朱使用了四个步骤来介绍排序增量法,第一个步骤,将所有半平面分成两个部分,其中一部分的极角为(-1/2π,1/2π],另外一部分为(-π,-1/2π]∪(1/2π,π],这样划分,如果最终的交集不为空,那么从最终效果来看,极角为(-1/2π,1/2π]的半平面构成了结果凸包的右侧部分,极角为(-π,-1/2π]∪(1/2π,π]的半平面构成了结果凸包的左半部分。在第二个步骤中,朱使用图示说明了求(-1/2π,1/2π]的半平面的交,但是图示实际上展示的半平面已经包括了一部分属于(1/2π,π]的半平面。在第四个步骤中,朱写道:“相邻半平面的交点组成半个凸多边形,我们有两个点集,(-1/2π,1/2π]给出上半个,(-π,-1/2π]∪(1/2π,π]给出下半个”。这样描述并不准确,按照论文中的设定,(-1/2π,1/2π]给出的是右半凸包,(-π,-1/2π]∪(1/2π,π]给出的是左半凸包。论文中也是自相矛盾的,因为紧接着,朱马上提到:“初始时候,四个指针p1,p2,p3和p4指向上/下凸壳的最左最右边”。在论文中,朱还遗漏了一个小细节,按照他的描述,得到的是两个半凸包,在最后进行处理时,从论文的表述来看,似乎是分别对两部分凸包分开进行处理,但是这样可能会出现一些特殊情形,使得无法得到正确的凸包(就像您在洛谷P4196中所使用的方法),需要将左右凸包合并起来后进行检查,这样才能够避免两部分凸包最终是有效的。

理解排序增量法的核心在于定义凸多边形边的极角,并对极角进行排序,从凸多边形的交的几何图像不难看出,交一定是由更靠近“中心”的那些边构成的。想象一下,假设有多个凸多边形相交,它们的交集不为空,根据几何性质可以知道相交区域也是一个凸多边形。把相交后构成的凸多边形称之为“核心”,那么直观上不难理解(尽管严格的证明或许并不容易):“核心”是由构成这n个凸多边形的线段的延长而成的直线相交而成的,而且最终的“核心”是由这n个凸多边形的边线中的某一些边线构成。如果有两条边线,它们的极角相同,不难理解,距离“核心”越近的边线应该被选择。这可通过想象用刀切蛋糕的情形来感性的理解。在代码上,如果两条有向直线的极角相同,那么沿着极角方向,靠近左侧的直线应该被选择。

以下是洛谷P4196的参考解题代码,我使用标准的双端队列写的,加了比较详细的注释。

#include <bits/stdc++.h>

using namespace std;

const int MAXV = 1100;
const double EPSILON = 1E-7;

// 使用向量来表示点。
struct point {
    double x, y;
    point (double x = 0, double y = 0): x(x), y(y) {}
    point operator + (point p) { return point(x + p.x, y + p.y); };
    point operator - (point p) { return point(x - p.x, y - p.y); };
    point operator * (double k) { return point(x * k, y * k); };
    point operator / (double k) { return point(x / k, y / k); };
};

// 定义凸包为点数组。
typedef vector<point> polygon;

// 向量的范数。
double norm(point a)
{
    return a.x * a.x + a.y * a.y;
}

// 向量的模(长度)。
double abs(point a)
{
    return sqrt(norm(a));
}

// 向量的点积。
double dot(point a, point b)
{
    return a.x * b.x + a.y * b.y;
}

// 向量的叉积。
double cross(point a, point b)
{
    return a.x * b.y - a.y * b.x;
}

// 直线的定义,使用两个点来表示一条直线,angle表示直线的极角。
struct line
{
    point a, b;
    double angle;
};

// 叉积。
double cp(point a, point b, point c)
{
    return cross(b - a, c - a);
}

// 顺时针旋转判断谓词。如果函数返回真,表明从点a向点b望去,
// 点c在a到b构成的有向直线的右侧,亦即从a走到b,再走到c构成一个右转。
bool cw(point a, point b, point c)
{
    return cp(a, b, c) < -EPSILON;
}

// 直线排序比较函数。当两条直线的极角相同时,沿着极角方向位于左侧的
// 那条直线排列在前,否则按极角升序排列。
bool cmpLine(line p, line q)
{
    if (fabs(p.angle - q.angle) <= EPSILON) return cw(p.a, p.b, q.a);
    return p.angle < q.angle;
}

// 用于对有向直线进行去重的比较函数。
bool cmpAngle(line p, line q)
{
    return fabs(p.angle - q.angle) <= EPSILON;
}

// 判断两条直线是否平行。
bool parallel(line p, line q)
{
    return fabs((p.a.x - p.b.x) * (q.a.y - q.b.y) -
        (q.a.x - q.b.x) * (p.a.y - p.b.y)) <= EPSILON;
}

// 计算两条直线的交点。
point getIntersection(line p, line q)
{
    point p1 = p.a;
    double scale =
        ((p.a.x - q.a.x) * (q.a.y - q.b.y) - (p.a.y - q.a.y) * (q.a.x - q.b.x)) /
        ((p.a.x - p.b.x) * (q.a.y - q.b.y) - (p.a.y - p.b.y) * (q.a.x - q.b.x));
    p1.x += (p.b.x - p.a.x) * scale;
    p1.y += (p.b.y - p.a.y) * scale;
    return p1;
}

// 将给定的两个点转换为一条直线。
line pointToLine(point a, point b)
{
    line lr;
    lr.a = a, lr.b = b, lr.angle = atan2(b.y - a.y, b.x - a.x);
    return lr;
}

// 计算凸包的面积。
double getArea(polygon pg)
{
    if (pg.size() < 3) return 0.0;
    double A = 0.0;
    int n = pg.size();
    for (int i = 0, j = (i + 1) % n; i < n; i++, j = (i + 1) % n)
        A += (pg[i].x * pg[j].y - pg[j].x * pg[i].y);
    return fabs(A / 2.0);
}

// 排序增量法确定半平面交。sides为有向直线数组,nLine为有向直线的数量。
polygon halfPlaneIntersection(line *sides, int nLine)
{
    // pg为最终凸包。
    polygon pg;
    // 使用双端队列来实现。
    line deq[MAXV];

    // 先将给定的有向直线按极角升序排列。
    sort(sides, sides + nLine, cmpLine);
    // 对于极角相同的有向直线只保留一条,即沿着极角方向位于左侧的那条有向直线。
    nLine = unique(sides, sides + nLine, cmpAngle) - sides;
    
    // 初始时,将排序后的前两条直线置入双端队列中。
    int btm = 0, top = 1;
    deq[0] = sides[0], deq[1] = sides[1];

    // 从第三条有向直线开始,按指定规则对有向直线进行筛选。
    for (int i = 2; i < nLine; i++)
    {
        // 如果发现双端队列的前端或者后端的两条有向直线是平行的,就可以直接退出。
        // 在此种情况下,给定的有向直线不存在合法的交。为什么?因为两条有向直线平行,
        // 说明它们是方向相反的(之前通过处理,已经使得不存在两条极角相同的有向直线),
        // 两条方向相反且相邻的有向直线无法构成有效的凸多边形交。
        if (parallel(deq[top], deq[top - 1]) || parallel(deq[btm], deq[btm + 1]))
            return pg;

        // 对于给定的一条直线sides[i],如果双端队列首端的两条直线的交点位于
        // sides[i]的右侧(即构成一个右转),那么说明sides[i]更靠近“核心”,
        // 因此队首的那条有向直线需要剔除。
        while (btm < top &&
            cw(sides[i].a, sides[i].b, getIntersection(deq[top], deq[top - 1])))
            top--;

        // 同理,对于给定的一条直线sides[i],如果双端队列末端的两条直线的
        // 交点位于sides[i]的右侧(即构成一个右转),那么说明sides[i]更靠近“核心”,
        // 因此队末的那条有向直线需要剔除。
        while (btm < top &&
            cw(sides[i].a, sides[i].b, getIntersection(deq[btm], deq[btm + 1])))
            btm++;

	// 将有向直线sides[i]置入双端队列的首端。
        deq[++top] = sides[i];
    }

初步筛选之后,剩下的有向直线如下图所示(以你给的hack数据作为输入数据):
在这里插入图片描述由几何图像可知,并不存在相交区域,这可以通过后续的检查予以确定。

    // 此时双端队列存储的是结果凸包的下半部分和上半部分。但是这两部分可能并不存在交集,
    // 因此需要通过进一步的检查来予以确定。具体方法和前述的思路是类似的,
    // 即对于队列末端的一条直线dep[btm]来说,如果队列首端的两条直线deq[top]和deq[top - 1]的
    // 交点在dep[btm]的右侧,那么显然,队列首端的有向直线dep[top]需要剔除。
    while (btm < top &&
        cw(deq[btm].a, deq[btm].b, getIntersection(deq[top], deq[top - 1])))
        top--;

    // 反过来也是类似的。
    while (btm < top &&
        cw(deq[top].a, deq[top].b, getIntersection(deq[btm], deq[btm + 1])))
        btm++;

    // 退出条件,只剩下两条有向直线时,无法构成凸包。
    if (top <= (btm + 1)) return pg;

    // 此时经过检查和剔除,双端队列中所存储的有向直线就是构成凸多边形交的那些有向直线,
    // 确定其交点即可。
    for (int i = btm; i < top; i++)
        pg.push_back(getIntersection(deq[i], deq[i + 1]));
    // 注意不能忘了首尾两端直线的交点。
    if (btm < (top + 1))
        pg.push_back(getIntersection(deq[btm], deq[top]));

    return pg;
}

int main(int argc, char *argv[])
{
    int n, m;
    int fx, fy, sx, sy, nx, ny;
    int cnt = 0;
    line sides[MAXV];

    cin >> n;
    for (int i = 0; i < n; i++)
    {
        cin >> m;
        cin >> fx >> fy;
        sx = fx, sy = fy;
        for (int j = 1; j < m; j++)
        {
            cin >> nx >> ny;
            sides[cnt++] = pointToLine(point(sx, sy), point(nx, ny));
            sx = nx, sy = ny;   
        }
        sides[cnt++] = pointToLine(point(nx, ny), point(fx, fy));
    }

    cout << fixed << setprecision(3) << getArea(halfPlaneIntersection(sides, cnt)) << '\n';
}

现在来看您的代码。为了便于说明您代码实现的问题,以您给的那组hack数据作为代码的输入:

3
4
0 0
1 0
2 1
0 2
4
0 7
1 4
4 3
5 4
5
-8 -6
-10 -10
-5 -15
5 -10
6 -4

其几何图像为:
Hack数据所对应的几何图像

总的来说,您的实现在整体上是正确的,至不过在细节上出了小问题。代码也可以进行进一步简化和优化。您先将给定的凸多边形的边转换为有向直线,然后对有向直线进行排序,之后按照朱的论文中所提到的方法构建凸包的交,即两个部分的凸壳。我的注释使用“//”开始,与您原来的注释容易区分。

/*这种做法结尾合并要求两个壳必须有交才能保证正确...*/
/*然而我不知道怎么 On 判 | upd:其实可以用双端队列做法暴力判...但那样就没意义了...*/
/**/
/*直线左侧平面交*/
#include <cstdio>
#include <cmath>
#include <algorithm>
using std::sort;
using std::isnan;

/*------------------------------Computational geometry------------------------------*/

const double pi =acos(-1), eps =1e-6;

struct vect{
    /*用极坐标存的话另外一些存坐标方便运算的地方也会变麻烦*/
    double x, y;
    //double m()/*magnitude*/ { return x/cos(atan(y/x)); };
    double m()/*magnitude*/{ return sqrt(x*x+y*y); };
    vect(){}
    vect(double xx, double yy):x(xx), y(yy){}
    vect operator + (vect b){ return vect(x+b.x, y+b.y); }
    vect operator - (vect b){ return vect(x-b.x, y-b.y); }
    vect operator * (double u){ return vect(u*x, u*y); }
    double operator * (vect b){ return x*b.x+y*b.y; }/*点积*/
    double operator / (vect b){ return x*b.y-y*b.x; }/*叉积*/
};

struct line{
    vect a, b;
    double angle;
    line(){}
    line(vect aa, vect bb):a(aa), b(bb){ angle =atan2(bb.y-aa.y, bb.x-aa.x); }
};

/**/
inline short gtr/*less*/(double a, double b){ return (a-b > eps); }

inline bool eq(double a, double b){ return (a-b < eps && a-b > -eps); }

inline bool inright/*点是否在有向直线右侧*/(line f, vect x){ return (gtr((x-f.a)/(f.b-f.a), 0)); }

int cmp(line A, line B){
    if(eq(A.angle, B.angle)) return inright(B, A.a);/*有向直线最右的在最后面,会被保留*/
    else return (gtr(A.angle, B.angle));
}

/*double 0 不会 RE,会 nan;平行时*/
vect getIntersection(line f, line g){
    double x =((g.a-f.a)/(f.a-f.b))/((f.a-f.b)/(g.a-g.b));
    return g.a+(g.a-g.b)*x;
}

/*------------------------------Main------------------------------*/

line li[550];/*排序后:前:右半圆正(向内/左),后:左半圆负;按角降序排*/
int tot;
line stkl[550], stkr[550], stks[550];
int totl, totr, tots;

inline int read(){
    int x =0; bool f =0; char c =getchar();
    while(c < '0' || c > '9') (c == '-') ? f =1, c =getchar() : c =getchar();
    while(c >= '0' && c <= '9') x = (x<<3) + (x<<1) + (48^c), c =getchar();
    return (f) ? -x : x;
}

inline void out(line *l, line *r){
    printf("------------------------------\n");
    for(; l !=r; ++l)
        printf("(%.2lf, %.2lf) -> (%.2lf, %.2lf)\n", l->a.x, l->a.y, l->b.x, l->b.y);
    printf("------------------------------\n");
}

/*上半凸壳下半凸壳太难实现,我就干脆左右搞了*/
/*右半圆从下往上转,左半圆从上往下转(小到大排序)*/

// 实际上,您的代码实现就是上下凸壳而不是左右凸壳,原因见后续代码我的注释。

int main(){
    //freopen("P4196.in", "r", stdin);
    int n =read();
    for(int i =0; i < n; ++i){
        int J =read(), fstx =read(), fsty =read(), lstx =fstx, lsty =fsty;
        for(int j =0; j < J-1; ++j){
            int x =read(), y =read();
            li[tot++] =line(vect(lstx, lsty), vect(x, y));
            lstx =x, lsty =y;
        }
        li[tot++] =line(vect(lstx, lsty), vect(fstx, fsty));
    }
    li[tot++] =line(vect(1500, -1500), vect(1500, 1500));
    li[tot++] =line(vect(1500, 1500), vect(-1500, 1500));
    li[tot++] =line(vect(-1500, 1500), vect(-1500, -1500));
    li[tot++] =line(vect(-1500, -1500), vect(1500, -1500));

    //out(&li[0], &li[tot]);/**/
    sort(li, li+tot, cmp);
    /*V 求左半右半凸壳 V*/
    //out(&li[0], &li[tot]);/**/

    // 在您的实现中,如果两条有向直线的极角不同,极角较大的有向直线排列在前,极角较小的有向直线排列在后。
    // 如果两条有向之间的极角相同,则远离“核心”的有向直线排列在前,靠近“核心”的有向直线排列在后。
    // 与标准的实现不太一样。不过,这并没有什么问题。只不过在后续处理时,注意选择靠近“核心”的有向直线即可。
    // 在后续代码中,您用以下语句来实现:
    // 
    // while(i < tot-1 && eq(li[i].angle, li[i+1].angle)) ++i;
    // 
    // 更好的办法就是在初始时直接将这些“无用的”有向直线直接剔除掉。因为按照排序增量法,当两天有向直线
    // 的极角相同时,应该选择极角方向更靠近左侧的那条直线。

    // 以下是排序后的有向直线,最后的浮点数表示有向直线的极角,单位为弧度。

    // (1500.00, 1500.00) -> (-1500.00, 1500.00): 3.14
    // (2.00, 1.00) -> (0.00, 2.00): 2.68
    // (5.00, 4.00) -> (0.00, 7.00): 2.60
    // (1500.00, -1500.00) -> (1500.00, 1500.00): 1.57
    // (5.00, -10.00) -> (6.00, -4.00): 1.41
    // (1.00, 0.00) -> (2.00, 1.00): 0.79
    // (4.00, 3.00) -> (5.00, 4.00): 0.79
    // (-5.00, -15.00) -> (5.00, -10.00): 0.46
    // (-1500.00, -1500.00) -> (1500.00, -1500.00): 0.00
    // (0.00, 0.00) -> (1.00, 0.00): 0.00
    // (1.00, 4.00) -> (4.00, 3.00): -0.32
    // (-10.00, -10.00) -> (-5.00, -15.00): -0.79
    // (0.00, 7.00) -> (1.00, 4.00): -1.25
    // (-1500.00, 1500.00) -> (-1500.00, -1500.00): -1.57
    // (0.00, 2.00) -> (0.00, 0.00): -1.57
    // (-8.00, -6.00) -> (-10.00, -10.00): -2.03
    // (6.00, -4.00) -> (-8.00, -6.00): -3.00


    // 以下语句实际构建的是上凸包,因为您加了语句: gtr(li[i].angle, 0),也就是说,只有当极角大于 0 的有向直线才进入这个 for 循环。
    // 极角大于 0 的有向直线构成的就是上凸包了,而不是右侧凸包。

    int i =0;
    for(; i < tot && gtr(li[i].angle, 0); ++i){
        while(i < tot-1 && eq(li[i].angle, li[i+1].angle)) ++i;
        while(totr >= 2 && !inright(stkr[totr-1], getIntersection(li[i], stkr[totr-2]))) --totr;
        stkr[totr++] =li[i];
    }

    // 以下语句实际构建的是下凸包,因为您之前排序是按极角递减排序的,前一个for循环只是取用极角大于0的有向直线,一旦遇到极角小于等于
    // 0的有向之间就退出for循环,转而使用极角小于0的有向之间构建凸包,而极角小于0的有向直线构建得到的是下凸包,不是您之前注释所写的
    // 构建的是左侧凸包。

    for(; i < tot; ++i){
        while(i < tot-1 && eq(li[i].angle, li[i+1].angle)) ++i;
        while(totl >= 2 && !inright(stkl[totl-1], getIntersection(li[i], stkl[totl-2]))) --totl;
        stkl[totl++] =li[i];
    }

    // 输出左侧凸壳,实际上是下凸壳,因为这些有向直线的极角均小于等于0。其结果为:
    // (0.00, 0.00) -> (1.00, 0.00): 0.00
    // (6.00, -4.00) -> (-8.00, -6.00): -3.00

    out(&stkl[0], &stkl[totl]);/**/

其几何图像为:
下凸壳

    // 输出右侧凸壳,实际上是上凸壳,因为这些有向直线的极角均大于0。其结果为:
    // (1500.00, 1500.00) -> (-1500.00, 1500.00): 3.14
    // (2.00, 1.00) -> (0.00, 2.00): 2.68
    // (4.00, 3.00) -> (5.00, 4.00): 0.79
    // (-5.00, -15.00) -> (5.00, -10.00): 0.46

    out(&stkr[0], &stkr[totr]);/**/

其几何图像为:
上凸壳

    // 合并上下凸包,这里是关键,如果操作不慎,就很容易出现Bug。
    // 由前述过程可知,得到的上下凸壳可能并不存在相交区域,那么如何判断不存在相交区域?标准的做法还是使用类似的思路
    // 检查栈顶的两条直线的交点是否在栈底直线的右侧,即构成栈底直线的两个点a,b,与栈顶两条直线的交点c,是否构成
    // 一个右转,如果构成右转,说明栈底的直线更靠近“核心”,栈顶的直线更加远离“核心”,因此栈顶的那条直线应该被剔除。
    // 类似的,需要使用栈顶的直线来测试位于栈底的两条直线交点,如果存在同样的情形,那么栈底的直线需要被剔除。

    /*V 合并凸壳 V*/
    int lp1 =totl-1, rp1 =0, lp2 =0, rp2 =totr-1;/*注意这里顺序,因为要求同组指针斜率相邻*/
    for(bool upd =1/*是否有新的更新*/; upd;){
        upd =0;
        if(lp1 != lp2 && inright(stkr[rp1], getIntersection(stkl[lp1], stkl[lp1-1]))) --lp1, upd =1;
        if(rp1 != rp2 && inright(stkl[lp1], getIntersection(stkr[rp1], stkr[rp1+1]))) ++rp1, upd =1;
    }
    for(bool upd =1; upd;){
        upd =0;
        if(lp1 != lp2 && inright(stkr[rp2], getIntersection(stkl[lp2], stkl[lp2+1]))) ++lp2, upd =1;
        if(rp1 != rp2 && inright(stkl[lp2], getIntersection(stkr[rp2], stkr[rp2-1]))) --rp2, upd =1;
    }
    /*如果没交(点)这里一定可以查出来,但如果有交点其实还会到下面,不过答案是对的( 0.000)| upd:没交可能会炸...这里有交点指的是交为一个点*/
    if(lp1 == lp2 && rp1 == rp2) return puts("0.000") && 0;
    
    // 这里的判断条件并不能保证没有交就一定可以检查出来,为什么?因为您是分别对上半凸包和下半凸包进行检查,很有可能出现
    // 这样一种情形:经过检查,上半凸壳只剩下一条直线,下半凸壳剩下两条直线(或者相反,上半凸壳剩下两条直线,下半凸壳剩
    // 下一条直线),按照您使用的“更新”规则,当前已经不存在“更新”,因此退出了两个for循环,但是这三条直线是有可能并不存在交的。
    // 您自己列出的Hack数据已经证明了这一点,对于这组hack数据来说,只满足lp1==lp2的条件,但不满足rp1==rp2的条件。
    // 要想避免这种情况,应该先将上下凸壳合并到一个数组中,然后再使用您所写的“更新”检查方法。从本质上来说,您所写的“更新”
    // 检查方法和标准的检查方法是相同的,只不过把问题复杂化了,因为只需要两端相互检查即可。

    /*如果题目不保证面积不无穷这里还要全遍历一遍看看有没有用到边界线*/
    /*V 求面积 V*/

    // 输出左侧凸壳,实际上是下凸壳,包含一条直线:
    // (0.00, 0.00) -> (1.00, 0.00): 0.00

    out(&stkl[lp2], &stkl[lp1+1]);/**/
    
    // 输出右侧凸壳,实际上是上凸壳,包含两条直线:
    // (2.00, 1.00) -> (0.00, 2.00): 2.68
    // (4.00, 3.00) -> (5.00, 4.00): 0.79
    // 从几何图像可以容易看出,有向直线(0.00, 0.00) -> (1.00, 0.00)和(2.00, 1.00) -> (0.00, 2.00)的交点
    // 在有向直线(4.00, 3.00) -> (5.00, 4.00)的右侧,即构成一个右转,按照规则,有向直线(0.00, 0.00) -> (1.00, 0.00)
    // 应该被剔除。但您的代码无法进行这样的检查,因为您代码的判断逻辑并未包含这样的情形。改进的方法一个是先将凸壳
    // 合并再判断,或者将判断代码扩展包含这种情形,明显第一种方法更为简便。

    out(&stkr[rp1], &stkr[rp2+1]);/**/

    double S =0;
    for(i =lp1; i >= lp2; --i) stks[tots++] =stkl[i];/*栈内都是降序*/
    for(i =rp2; i >= rp1; --i) stks[tots++] =stkr[i];

    // 将构成凸多边形交的直线合并为一个数组,其结果为:
    // (0.00, 0.00) -> (1.00, 0.00): 0.00
    // (4.00, 3.00) -> (5.00, 4.00): 0.79
    // (2.00, 1.00) -> (0.00, 2.00): 2.68
    // 实际上应该是先合并再检查,否则就可能出现Bug。

    int p =0;
    vect v0 =getIntersection(stks[p], stks[p+1]), v1, v2 =getIntersection(stks[p+1], stks[p+2]);
    for(; p+3 < tots; ++p){
        v1 =v2, v2 =getIntersection(stks[p+2], stks[p+3]);
        S +=(v1-v0)/(v2-v1);
    }
    v1 =v2, v2 =getIntersection(stks[p+2], stks[0]);
    S +=(v1-v0)/(v2-v1);
    printf("%.3lf", S/2);
}

您所给出的CF上的题目,大意是给定一个凸多边形,要求在凸多边形内放置两个同样大小的圆形,且两个圆覆盖的面积要尽可能的大。首先需要得到圆心可能所在的区域,这个可以通过将给定的凸多边形的每条边沿着凸多边形内部法线方向移动距离r,然后求这些有向直线的交所得到的凸多边形核心,这个凸多边形就是圆心可能的位置。题目要求两个圆覆盖的面积尽可能的大(可以重叠),等价于两个圆的圆心距离尽可能的远,也就是说,可以将问题转化为确定凸多边形核心所围区域(包括边界)内的两个点,这两个点具有最大的距离,显然,这是求凸多边形核心的直径,不过由于题目的数据规模较小,可以直接枚举凸多边形的两个顶点来确定其直径,而不需要使用旋转卡壳算法。

下述实现之所以是正确的,是因为在构建凸包的过程中,所有有向直线一直都是在同一个栈中,这样您在使用“更新”迭代剔除无效直线的时候,能够保证将所有不符合要求的直线全部剔除掉,而不像之前那样的实现,只能将两部分凸壳的无效直线去除,而当这两部分凸壳合并后,仍然存在无效的直线,使得算法失效。

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1005;
const double eps = 1e-8;
typedef struct point {
    double x, y;
}vect;
point p[MAXN];
struct line {
    vect p;
    point u, v;
}Sline[MAXN];
vect operator - (const point a, const point b) {
    return (vect){a.x - b.x, a.y - b.y};
}
int pnz(double x) {
    return x < -eps ? -1 : x > eps ? 1 : 0;
}
double dist(const point a, const point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double Cross(const vect a, const vect b) {
    return a.x * b.y - a.y * b.x;
}
int Onleft(const point p, const line l) {
    return pnz(Cross(l.p, p - l.u));
}
point Inter(const line a, const line b) {
    double t = Cross(b.p, a.u - b.u) / Cross(a.p, b.p);
    return (point){a.u.x + t * a.p.x, a.u.y + t * a.p.y};
}
void Mobile(line l[], int n, int r) {
    for (int i = 0; i < n; i++) {
        double dis = dist(l[i].u, l[i].v);
        double dx = (l[i].v.y - l[i].u.y) / dis * r;
        double dy = (l[i].v.x - l[i].u.x) / dis * r;
        l[i].u.x -= dx, l[i].u.y += dy;
        l[i].v.x -= dx, l[i].v.y += dy;
    }
}
void Farthest(point p[], int n) {
    int x = 0, y = 0;
    double max_ = 0;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            double dis = dist(p[i], p[j]);
            if (pnz(dis - max_) > 0) {
                max_ = dis;
                x = i, y = j;
            }
        }
    }
    printf("%.5f %.5f %.5f %.5f\n", p[x].x, p[x].y, p[y].x, p[y].y);
}

int Q[10000];

void Halfplane(int n) {
    point Spt[n];
    memset(Q, 0, sizeof(Q));
    int l = 0, r = 0;
    for (int i = 1; i < n; i++) {
        while (l < r && Onleft(Spt[r - 1], Sline[i]) < 0)
            r--;
        //while (l < r && Onleft(Spt[l], Sline[i]) < 0)
            //l++;
        Q[++r] = i;
        if (!pnz(Cross(Sline[i].p, Sline[Q[r - 1]].p))) {
            r--;
            if (Onleft(Sline[i].u, Sline[Q[r]]) > 0)
                Q[r] = i;
        }
        if (l < r) Spt[r - 1] = Inter(Sline[Q[r]], Sline[Q[r - 1]]);
    }
    
    // 为什么这里是对的?因为在队列Q中,已经包含了经过筛选的有向直线,这些直线构成结果凸包的整体,但是
    // 需要进行进一步的检查,也就是类似于标准的做法,检查队列首端(假定l为尾端,r为首端)的两条直线的
    // 交点是否位于队列尾端直线Sline[Q[l]]的右侧,如果Spt[r - 1]位于Sline[Q[l]]的右侧,
    // 说明Sline[Q[l]]更靠近“核心”,交点Spt[r - 1]需要丢弃。

    for(bool upd =1; upd && l < r;){
        upd =0;
        if(l < r && Onleft(Spt[r - 1], Sline[Q[l]]) < 0) r--, upd =1;
        if(l < r && Onleft(Spt[l], Sline[Q[r]]) < 0) l++, upd =1;
    }
    Spt[r] = Inter(Sline[Q[l]], Sline[Q[r]]);
    Farthest(Spt + l, r - l + 1);
}
int main() {
    int n, r;
    while (~scanf("%d%d", &n, &r)) {
        for (int i = n - 1; i >= 0; i--)
            scanf("%lf%lf", &p[i].x, &p[i].y);
        p[n] = p[0];
        for (int i = 0; i < n; i++)
            Sline[i] = (line){p[i + 1] - p[i], p[i], p[i + 1]};
        Mobile(Sline, n, r);
        Halfplane(n);
    }
    return 0;
}

通过阅读您的代码,发现您在写计算几何时并未形成自己固定的代码模块,似乎是随用随写,这并不是一个好的习惯。因为计算几何对代码的鲁棒性(robustness,健壮性)要求比较高,对于基本的操作,如果将其编写为单个经过充分测试的小模块,在解题时能够避免错误。建议您可以按照我写的书中的那样,将基本的向量操作、点积、叉积、判断转向都写成一个小函数(可以将之视为基本的“元件”,使用类似于搭积木的方法将这些“元件”组合成功能更复杂的模块),然后再整合起来,这样在参加网络比赛时能够随时使用,不需要太多更改,而且也不容易出错。

综上,我想我已经清楚地回答了您的问题。如果还有进一步的疑问,您可以@我。

  • 2
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值