USACO 6.5.3 Closed Fences 计算几何

http://train.usaco.org/usacoprob2?a=AFsbkYO7d5R&S=fence4

题目大意：按点的顺时针顺序给出一个由N条线段（N个点，3<N<200）组成的闭合图形，给出一个观察点（不在图形边界上），求从观察点望去可以被看到的所有边（边没有厚度，如果侧对观察点视为不可见），并按顺序输出。用两个点的四个坐标数字描述一条边，同时在输入文件里先出现的点在前，在边之间排序时，要求按照第二个点在输入文件中的顺序排序，如果相同再按照第一个点的顺序排序；

思路很直观，遍历每一条边看能不能找到该边上的一点和观察点的连线与其他任意一条边不相交。

然后因为这是一个封闭图形，所以对于某一条边的遮挡必然是从两边开始连续的向中间延伸的，所以我本能地想到了二分搜索，因为一旦确定一个点被遮挡（该点与观察点的连线与别的边相交）那么从这个点开始到端点的一整段都必然被遮挡。看上去稳得不行，但在移动边界时就发生了问题，因为不知道该移动哪一个边界——所以说使用这样的二分搜索的条件是单调性的原因不在于能够批量地去除错误解而在于能够有效地确定合法解的位置，后者是算法合理性的基础，前者是算法高效性的保证。显然在这个问题下是不能确定合法解位置的，所以我考虑使用三分的方法，因为中心合法嘛，但是没法比较两个位置的优劣，所以也不行。

所以我看了别人的题解才恍然大悟可以使用分治法搜索合法解，当尝试区间足够小（一定要足够小，我就因为区间设置的太大在其中一个样例上卡了好久）时就返回失败，一旦有一个区间内找到了一个合法解则成功。我以前一直搞不清楚分治法和二分搜索有什么区别，尤其是很多场合下分治法确实就是二分的，现在看来二分搜索应该是特殊的分治法因为它的解的位置是能够（大致）确定的，可以说是剪了枝的分治法，而像这种几乎没有特殊性的题目由于没法剪枝所以只能采用分治法而不能二分搜索。

参考资料：http://blog.youkuaiyun.com/lyd_7_29/article/details/48998035

在计算几何学方面，我操作还是很不熟练，有时候想偷懒就出了错，比如在叉积为零时判断点是否在线段上，我就想当然的仅仅判断了一个坐标，太naive了。

顺便：我存图形时是存的点，所以开了250够用，但是在存答案时我是存的线段，所以一开始只开了200就越界了（笑

/*
ID: frontie1
TASK: fence4
LANG: C++
*/
#include <iostream>
#include <cstdio>

using namespace std;

struct point
{
    double x;
    double y;
};

double crossproduct(point p1, point p2, point p3, point p4) //p1p2 X p3p4
{
    double a_x = p2.x-p1.x, a_y = p2.y-p1.y;
    double b_x = p4.x-p3.x, b_y = p4.y-p3.y;

    return (a_x*b_y-a_y*b_x);
}

int direction(point p1, point p2, point p3) // p1 p2p3
{
    double tem = crossproduct(p2, p3, p2, p1);
    if(tem > 0) return 1;       //left
    else if(tem < 0) return -1;  //right
    else return 0;              //on
}

bool onsegment(point p1, point p2, point p3)
{
    return (p1.x >= min(p2.x, p3.x) && p1.x <= max(p2.x, p3.x) && p1.y >= min(p2.y, p3.y) && p1.y <= max(p2.y, p3.y));
}

bool intersect(point p1, point p2, point p3, point p4) //p1p2 p3p4
{
    int dir1 = direction(p1, p3, p4);
    int dir2 = direction(p2, p3, p4);
    int dir3 = direction(p3, p1, p2);
    int dir4 = direction(p4, p1, p2);

    if(dir1*dir2*dir3*dir4 != 0 && dir1 != dir2 && dir3 != dir4) return true;

    if(dir1 == 0 && onsegment(p1, p3, p4)) return true;
    if(dir2 == 0 && onsegment(p2, p3, p4)) return true;
    if(dir3 == 0 && onsegment(p3, p1, p2)) return true;
    if(dir4 == 0 && onsegment(p4, p1, p2)) return true;

    return false;
}

double dist(point p1, point p2)
{
    double x = p1.x - p2.x;
    double y = p1.y - p2.y;
    return (x*x + y*y);
}

point arr[250];
int N;
point st;
int cnt = 0;
point output[500];

bool visible(point lo, point hi, int cur)
{
    if(dist(lo, hi) < 0.0001) return false;
    if(direction(st, lo, hi) == 0) return false;
    point mid;
    bool flag = true;
    mid.x = (lo.x + hi.x) / 2;
    mid.y = (lo.y + hi.y) / 2;
    for(int i = 0; i < N; ++i){
        if(i == cur) continue;
        if(intersect(st, mid, arr[i], arr[i+1])){
            flag = false;
            break;
        }
    }
    if(!flag) return (visible(lo, mid, cur) || visible(mid, hi, cur));
    else return true;
}

int main()
{
    freopen("fence4.in", "r", stdin);
    freopen("fence4.out", "w", stdout);

    cin >> N;
    cin >> st.x >> st.y;
    for(int i = 0; i < N; ++i){
        cin >> arr[i].x >> arr[i].y;
    }
    arr[N] = arr[0];

    for(int i = 0; i < N-2; ++i){
        if(visible(arr[i], arr[i+1], i)){
            output[2*cnt] = arr[i];
            output[2*cnt+1] = arr[i+1];
            ++cnt;
        }
    }
    if(visible(arr[0], arr[N-1], N-1)){
        output[2*cnt] = arr[0];
        output[2*cnt+1] = arr[N-1];
        ++cnt;
    }
    if(visible(arr[N-2], arr[N-1], N-2)){
        output[2*cnt] = arr[N-2];
        output[2*cnt+1] = arr[N-1];
        ++cnt;
    }

    cout << cnt << endl;
    for(int i = 0; i < cnt; ++i){
        cout << output[2*i].x << ' ' << output[2*i].y << ' ' << output[2*i+1].x << ' ' << output[2*i+1].y << endl;
    }

    return 0;
}