計算幾何學 | 線段的交點 | Cross Point | C/C++實作

問題描述

輸出線段s1、s2交點的坐标。

設s1的端點為p0、p1，s2的端點為p2、p3。

輸入：

第1行輸入問題數q。接下來q行給出q個問題。各問題線段s1、s2的坐标按照以下格式給出：

x p 0 x_{p0} xp0​ y p 0 y_{p0} yp0​ x p 1 x_{p1} xp1​ y p 1 y_{p1} yp1​ x p 2 x_{p2} xp2​ y p 2 y_{p2} yp2​ x p 3 x_{p3} xp3​ y p 3 y_{p3} yp3​

輸出：

根據各問題輸出交點坐标(x, y)，每個問題占1行。輸出允許誤差不超過0.00000001

限制：

1 ≤ q ≤ 100

-10000 ≤ x p i , y p i 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。

d 1 = ∣ b a s e × h y p o ∣ / ∣ b a s e ∣ d1=|base×hypo|/|base| d1=∣base×hypo∣/∣base∣

d 2 d2 d2同理

d 1 = ∣ b a s e × ( s 1. p 1 − s 2. p 1 ) ∣ / ∣ b a s e ∣ d1=|base×(s1.p1-s2.p1)|/|base| d1=∣base×(s1.p1−s2.p1)∣/∣base∣

d 2 = ∣ b a s e × ( s 1. p 2 − s 2. p 1 ) ∣ / ∣ b a s e ∣ d2=|base×(s1.p2-s2.p1)|/|base| d2=∣base×(s1.p2−s2.p1)∣/∣base∣

然後，設線段s1的長度與點s1.p1到交點x的距離之比為t，則有

d 1 : d 2 = t : ( 1 − t ) d1:d2=t:(1-t) d1:d2=t:(1−t)

由此可得

t = d 1 / ( d 1 + d 2 ) t=d1/(d1+d2) t=d1/(d1+d2)

是以交點x為

x = s 1. p 1 + ( s 1. p 2 − s 1. p 1 ) × t 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);
	}
} 
           

注：以上本文未涉及代碼的詳細解釋參見：計算幾何學