计算几何 - 凸包 - 判断某点是否可以被一点集中的某三个点围成的三角形包围

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N = 300005;

// 存储二维平面点
struct ip
{
    int x, y;
    ip(int x = 0, int y = 0) : x(x), y(y) { }
    void ri() { scanf("%d %d", &x, &y); }
};

// iv表示一个向量类型，其存储方式和ip一样
typedef ip iv;

// 先比较x轴再比较y轴
bool operator < (const ip &a, const ip &b)
{
    return a.x == b.x ? a.y < b.y : a.x < b.x;
} 

// 两点相减得到的向量
iv operator - (const ip &a, const ip &b)
{
    return iv(a.x - b.x, a.y - b.y);
}

// 计算a和b的叉积（外积）
ll operator ^ (const iv &a, const iv &b)
{
    return (ll)a.x * b.y - (ll)a.y * b.x;
}

// 计算a和b的点积（内积） 
ll operator * (const iv &a, const iv &b)
{
    return (ll)a.x * b.x + (ll) a.y * b.y;
}

class myComp
{
private:
    long long ltl_x, ltl_y;    // 存储左下角点的坐标 
public:
    myComp() : ltl_x(0), ltl_y(0) { }
    myComp(int x, int y){ ltl_x = (long long)x; ltl_y = (long long)y; }
    bool operator()(const ip & p1, const ip & p2)
    {
        long long x1 = (long long) p1.x - ltl_x, y1 = (long long) p1.y - ltl_y;
        long long x2 = (long long) p2.x - ltl_x, y2 = (long long) p2.y - ltl_y;
        long long outter_product = x1 * y2 - x2 * y1;
        if ( outter_product > 0 )
            return true;
        else if ( outter_product == 0 )    // 若LTL、p1和p2共线 
        {
            long long x3 = (long long) p2.x - (long long) p1.x;
            long long y3 = (long long) p2.y - (long long) p1.y;
            if ( x3 * x1 + y3 * y1 > 0 )    // 若 LTL->p1 与 p1->p2 方向相同（用内积判断），则p1到LTL距离更短 
                return true;
        }
        return false;
    }
};

/* 判断点 c 是否在 ab 向量的左侧。若 c 在 ab 上，或者 a→b 正向延长线上，或者 c 与 a 或 b 重合，返回皆为 true。 */
bool toLeft(ip & c, ip & a, ip & b)
{
    long long x_ab = (long long) b.x - (long long) a.x; 
    long long y_ab = (long long) b.y - (long long) a.y;
    long long x_ac = (long long) c.x - (long long) a.x; 
    long long y_ac = (long long) c.y - (long long) a.y;
    long long outter_product = x_ab * y_ac - x_ac * y_ab;
    long long inner_product = x_ab * x_ac + y_ab * y_ac;
    if ( outter_product > 0 || (outter_product == 0 && inner_product >= 0) )
        return true;
    return false;
} 

/* 计算二维点数组a的凸包，将凸包放入b数组中，下标均从0开始
   a, b：如上 
   n：表示a中元素个数
   返回：凸包元素个数 
*/
int convex(ip * a, ip * b, int n)
{
    // 升序排序
    sort(a, a+n);
    // n = unique(a, a+n) - a;    // 若题目中有重复点，必须去重
    
    int m = 0;
    
    // 求下凸壳
    for ( int i = 0; i < n; ++i )
    {
        for ( ; m > 1 && ((b[m-1] - b[m-2]) ^ (a[i] - b[m-2])) < 0; --m );
        b[m++] = a[i];
    }
    
    // 求上凸壳
    for ( int i = n - 2, t = m; i >= 0; --i )
    {
        for ( ; m > t && ((b[m-1] - b[m-2]) ^ (a[i] - b[m-2])) < 0; --m );
        b[m++] = a[i];
    }
     
    return m - 1;
}

/* 对第n_b个点构成的点集B(坐标从b开始存储)中每个点，判断是否能被n_a个点构成的点集A(坐标从a开始存储)中的三个点包围。 */
vector<int> isSurrounded(ip * a, ip * b, int n_a, int n_b);

ip A[N], B[N];

int main()
{   
    int n = 0, m = 0;    // 点集A中点的数量n和点集B中点的数量m 
    
    // 读入各点坐标 
    scanf("%d", &n);    
    for ( int i = 0; i < n; ++i )
        A[i].ri();
    scanf("%d", &m);    
    for ( int i = 0; i < m; ++i )
        B[i].ri();
    
    // 求解 
    vector<int> ans = isSurrounded(A, B, n, m);
    
    // 打印结果 
    for ( int i = 0; i < m; ++i )
        printf("%d\n", ans[i]);
    
    return 0;
}

vector<int> isSurrounded(ip * a, ip * b, int n_a, int n_b)
{
    vector<int> ans; ans.resize(n_b);
    
    // 求a的凸包 (如果不找到凸包，直接在整个点集A上执行下面的判断是不正确的)
    ip * conv = new ip[n_a];
    int n_conv = convex(a, conv, n_a);
    
    // 找到A中的LTL，并将它换到conv指向的位置 (即conv数组的第0个元素) 
    ip ltl(conv->x, conv->y);
    ip * p = conv;
    for ( int i = 1; i < n_conv; ++i )
    {
        if ( (conv+i)->y < ltl.y || ((conv+i)->y == ltl.y && (conv+i)->x < ltl.x ) )
        {
            ltl.x = (conv+i)->x;
            ltl.y = (conv+i)->y;
            p = conv+i;
        }
    }
    p->x = conv->x; p->y = conv->y;
    conv->x = ltl.x; conv->y = ltl.y;
    
    // 对conv中各点关于LTL进行极角排序
    sort(conv+1, conv+n_conv, myComp(ltl.x, ltl.y)); 
     
    // 对B中各点判断其是否能被A中的三个点包围
    for ( int i = 0; i < n_b; ++i )
    {
        if ( b[i].y < ltl.y || (b[i].y == ltl.y && b[i].x < ltl.x ) )
        {
            ans[i] = 0;
            continue;
        }
        
        // 找到A凸包中不小于b[i]的最小点a1
        ip * a1 = lower_bound(conv+1, conv+n_conv, b[i], myComp(ltl.x, ltl.y));
        ip * a0 = a1 - 1;
        if ( (a1 - a) == n_conv )    // 若找不到这样的a1，则不能被包围 
            ans[i] = 0;
        else
        {
            if ( toLeft(b[i], ltl, *a0) && toLeft(b[i], *a0, *a1) && toLeft(b[i], *a1, ltl) )
                ans[i] = 1;
            else
                ans[i] = 0;
        }
    } 
    
    return ans;
}