计算几何学 | 点的内包 | Polygon-Point Containment | C/C++实现

问题描述

对于多边形g与点p，p位于g内时输出2，p位于g的边上时输出1，其余情况输出0。

多边形g由其顶点的序列$p_1,p_2,...,p_n$p1​,p2​,...,pn​表示，相邻两点$p_i$pi​、$p_{i+1}$pi+1​（1 ≤ i ≤ n-1）相连构成g的边。另外，点$p_n$pn​和$p_1$p1​相连的线也是g的边。

请注意，g不一定是凸多边形。

输入：
输入按照以下格式给出：
$g$g（构成多边形的点的序列）
$q$q（问题数）
1st query
2nd query
…
qth query
$g$g为点的序列$p_1,...,p_n$p1​,...,pn​按照以下格式给出：
$n$n
$x_1$x1​ $y_1$y1​
$x_2$x2​ $y_2$y2​
…
$x_n$xn​ $n_y$ny​
第1行的n表示点的数量。点$p_i$pi​的坐标以2个整数$x_i$xi​、$y_i$yi​的形式给出。各点按照逆时针访问多边形相邻顶点的顺序排列。
输出：
根据各问题输出2、1或0，每个问题占1行。
限制：
3 ≤ n ≤ 100
1 ≤ q ≤ 1000
-10000 ≤ $x_i,y_i,x,y$xi​,yi​,x,y ≤ 10000
多边形各顶点坐标均不相同
多边形各边只在公共端点处相交。

输入示例

4
0 0
3 1
2 3
0 3
3
2 1
0 2
3 2

输出示例

2
1
0

讲解

只要检查以p为端点且平行于x轴的射线与多边形g的边的相交次数，我们就能判断出给定点p是否内包于多边形g。而且这条射线不必实际生成，只需要通过下面的算法即可完成判断。

对于构成多边形各边的线段$g_ig_{i+1}$gi​gi+1​，设$g_i-p$gi​−p与$g_{i+1}-p$gi+1​−p分别为向量a和向量b。点p是否位于$g_ig_{i+1}$gi​gi+1​上可通过类似ccw的方法（关于ccw的详细解释参见：计算几何学 | 逆时针方向 | Counter-Clockwise | C/C++实现）检查，即检查a和b是否在同一直线且方向相反。如果a和b外积大小为0且内积小于等于0，则点p位于$g_ig_{i+1}$gi​gi+1​上。

至于射线与线段$g_ig{i+1}$gi​gi+1是否相交（是否能更新内包状态），可以通过a、b构成的平行四边形的面积正负，即a、b外积的大小来判断。首先我们调整向量a、b重新将y值较小的向量定为a。在这个状态下，如果a和b的外积大小为正（a到b为逆时针）且a、b（的终点）位于射线两侧，则可确定射线与边交叉。这里要注意设置边界条件，避免射线与$g_ig_{i+1}$gi​gi+1​端点相交时对结果造成影响。

相交次数为奇数时表示“内包”，为偶数时表示“不内包”。此外，在判断射线与各线段是否相交时，一旦发现点p位于线段上，需要立刻返回“在线段上”。
判断多边形与点的内包关系的程序如下所示

点的内包：

/*
	IN 2
	ON 1
	OUT 0
*/
int contains(Polygon g, Point p) {
	int n = g.size();
	bool x = false;
	for(int i = 0; i < n; i++) {
		Point a = g[i] - p, b = g[(i+1) % n] - p;
		if( abs( cross(a,b) ) < EPS && dot(a, b) < EPS ) return 1;
		if( a.y > b.y ) swap(a, b);
		if( a.y < EPS && EPS < b.y && cross(a, b) > EPS ) x = !x;
	}
	return (x ? 2 : 0);
}

AC代码如下

#include<iostream>
#include<cmath>
#include<vector>
using namespace std;

#define EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS)

class Point {//Point类，点 
	public:
		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 a) { return Point(a * x, a * y); }
		Point operator / (double a) { return Point(x / a, y / a); }

		double abs() { return sqrt(norm()); }
		
		double norm() { return x * x + y * y; }
		
		bool operator < (const Point &p) const {
			return x != p.x ? x < p.x : y < p.y;
		}

		bool operator == (const Point &p) const {
			return fabs(x - p.x) < EPS && fabs(y - p.y) < EPS;
		}
};

typedef Point Vector;//Vector类，向量 

typedef vector<Point> Polygon;

double dot(Vector a, Vector b) {//内积 
	return a.x * b.x + a.y * b.y;
}

double cross(Vector a, Vector b) {//外积 
	return a.x*b.y - a.y*b.x;
}

/*
	IN 2
	ON 1
	OUT 0
*/
int contains(Polygon g, Point p) {
	int n = g.size();
	bool x = false;
	for(int i = 0; i < n; i++) {
		Point a = g[i] - p, b = g[(i+1) % n] - p;
		if( abs(cross(a, b) ) < EPS && dot(a, b) < EPS ) return 1;
		if( a.y > b.y ) swap(a, b);
		if( a.y < EPS && EPS < b.y && cross(a, b) > EPS ) x = !x;
	}
	return (x ? 2 : 0);
}

int main(){
	int n;
	cin>>n;
	
	Polygon g;
	Point temp;
	
	while(n--){
		cin>>temp.x>>temp.y;
		g.push_back(temp);
	}
	
	cin>>n;
	
	while(n--){
		cin>>temp.x>>temp.y;
		cout<<contains(g, temp)<<endl;
	}
} 

注：以上本文未涉及代码的详细解释参见：计算几何学