计算几何学 | 线段的交点 | Cross Point | C/C++实现

问题描述

输出线段s1、s2交点的坐标。
设s1的端点为p0、p1，s2的端点为p2、p3。

输入：
第1行输入问题数q。接下来q行给出q个问题。各问题线段s1、s2的坐标按照以下格式给出：
$x_{p0}$xp0​ $y_{p0}$yp0​ $x_{p1}$xp1​ $y_{p1}$yp1​ $x_{p2}$xp2​ $y_{p2}$yp2​ $x_{p3}$xp3​ $y_{p3}$yp3​
输出：
根据各问题输出交点坐标(x, y)，每个问题占1行。输出允许误差不超过0.00000001
限制：
1 ≤ q ≤ 100
-10000 ≤ $x_{p_i},y_{p_i}$xpi​​,ypi​​ ≤ 10000
p0、p1不是同一个点
p2、p3不是同一个点
s1、s2有交点，且量线段不重叠。

输入示例

3
0 0 2 0 1 1 1 -1
0 0 1 1 0 1 1 0
0 0 1 1 1 0 0 1

输出示例

1.0000000000 0.0000000000
0.5000000000 0.5000000000
0.5000000000 0.5000000000

讲解

两条线段s1、s2的交点坐标可以通过外积的大小求得。

我们用向量s2.p2 - s2.p1 = base来表示线段s2。接下来，分别求通过s2.p1和s2.p2的直线与线段s1两端点的距离d1、d2。举个例子，设s1.p1 - s2.p1为向量hypo，那么base与hypo构成的平行四边形的面积就是外积base×hypo的大小。这样一来，只要用面积处以base的大小即可求出d1。

$d1=|base×hypo|/|base|$d1=∣base×hypo∣/∣base∣
$d2$d2同理
$d1=|base×(s1.p1-s2.p1)|/|base|$d1=∣base×(s1.p1−s2.p1)∣/∣base∣
$d2=|base×(s1.p2-s2.p1)|/|base|$d2=∣base×(s1.p2−s2.p1)∣/∣base∣
然后，设线段s1的长度与点s1.p1到交点x的距离之比为t，则有
$d1:d2=t:(1-t)$d1:d2=t:(1−t)
由此可得
$t=d1/(d1+d2)$t=d1/(d1+d2)
所以交点x为
$x=s1.p1+(s1.p2-s1.p1)×t$x=s1.p1+(s1.p2−s1.p1)×t
求线段s1与线段s2交点的程序可以像下面这样写

线段s1与线段s2的交点：

Point getCrossPoint(Segment s1, Segment s2) {
	Vector base = s2.p2 - s2.p1;
	double d1 = abs(cross(base, s1.p1 - s2.p1));
	double d2 = abs(cross(base, s1.p2 - s2.p1));
	double t = d1 / (d1 + d2);
	return s1.p1 + (s1.p2 - s1.p1) * t;
}

上述程序在计算d1、d2的过程种会用到|base|，但这个数会在计算t时被约分消去。

AC代码如下

#include<stdio.h>
#include<iostream>
#include<cmath>
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类，向量 

struct Segment{//Segment 线段 
	Point p1, p2;
};
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;
}

Point getCrossPoint(Segment s1, Segment s2) {
	Vector base = s2.p2 - s2.p1;
	double d1 = abs(cross(base, s1.p1 - s2.p1));
	double d2 = abs(cross(base, s1.p2 - s2.p1));
	double t = d1 / (d1 + d2);
	return s1.p1 + (s1.p2 - s1.p1) * t;
}

int main(){
	int q;
	cin>>q;
	
	Segment s1, s2;
	Point p;
	
	while(q--){
		cin>>s1.p1.x>>s1.p1.y>>s1.p2.x>>s1.p2.y>>s2.p1.x>>s2.p1.y>>s2.p2.x>>s2.p2.y;
		p = getCrossPoint(s1, s2);
		printf("%.10f %.10f\n", p.x, p.y);
	}
} 

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