计算几何学 | 距离 | Distance | C/C++实现

问题描述

输出线段s1、s2之间的距离。设s1的端点为p0、p1，s2的端点为p2、p3。

输入：
第1行输入问题数q。接下来q行给出q个问题。各问题线段s1、s2的坐标按照以下格式给出：
$x_{p0}$ $y_{p0}$ $x_{p1}$ $y_{p1}$ $x_{p2}$ $y_{p2}$ $x_{p3}$ $y_{p3}$
输出：
根据各个问题输出距离，每个问题占1行。输出允许误差不超过0.00000001
限制：
1 ≤ q ≤ 1000
-10000 ≤ $x_{p_i},y_{p_i}$ ≤ 10000
p0、p1不是同一个点
p2、p3不是同一个点。

输入示例

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

输出示例

1.0000000000
1.4142135624
0.0000000000

讲解

我先不针对问题内容，而是对点、线段的相关距离进行讲解。

1.两点间的距离

点a与点b之间的距离等于向量a - b或b - a的绝对值。我们可以像下面这样设计程序来求点a和点b的距离。

点a和点b的距离：

double getDistance(Point a, Point b) {
	return abs(a - b);
}

2.点与直线的距离
设直线p1p2上的向量为a = p2 - p1，p与p1构成的向量为b = p - p1，则点p与直线p1p2的距离d就等于a、b构成的平行四边形的高。用a与b外积的大小（平行四边形的面积）除以a的大小|a|即可求出高d，于是有：
$d=|a\times b|/|a|=|(p2-p1)\times (p-p1)|/|p2-p1|$
我们可以这样设计程序来求直线$l$和点$p$的距离。

直线$l$和点$p$的距离：

double getDistanceLP(Line l, Point p) {
	return abs(cross(l.p2 - l.p1, p - l.p1) / abs(l.p2 - l.p1));
}

3.点与线段的距离
1.向量p2 - p1与向量p - p1的夹角θ大于90度（或小于-90度）时，d为点p到点p1的距离
2.向量p1 - p2与向量p - p2的夹角θ大于90度（或小于-90度）时，d为点p到点p2的距离
3.除上述两种情况外，d为点p到直线p1p2的距离。

由于θ大于90度时cosθ<0，所以可以通过2个向量的内积是否为负来判断第1、2种情况。我们可以这样设计程序来求线段s与点p的距离。

线段s与点p的距离：

double getDistanceSP(Segment s, Point p) {
	if( dot(s.p2 - s.p1, p - s.p1) < 0.0 ) return abs(p - s.p1);
	if( dot(s.p1 - s.p2, p - s.p2) < 0.0 ) return abs(p - s.p2);
	return getDistanceLP(s, p);
}

4.线段与线段的距离
线段s1与线段s2的距离为以下四个距离种最小的一个：

1.线段s1与线段s2的端点s2.p1的距离
2.线段s1与线段s2的端点s2.p2的距离
3.线段s2与线段s1的端点s1.p1的距离
4.线段s2与线段s1的端点s1.p2的距离

此外，如果两条线段相交，则他们的距离为0
我们可以这样设计程序来求线段s1与线段s2的距离。

线段s1与线段s2的距离：

double getDistance(Segment s1, Segment s2) {
	if(intersect(s1, s2) ) return 0.0;
	return min(min(getDistanceSP(s1, s2.p1), getDistanceSP(s1, s2.p2)), 
				min(getDistanceSP(s2, s1.p1), getDistanceSP(s2, s1.p2)));
}

上面程序种的intersect函数用于判断两条线段是否相交。

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;
};

typedef Segment Line;//Line 直线 

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;
}

double getDistance(Point a, Point b) {//点a和点b的距离
	return (a - b).abs();
}

double getDistanceLP(Line l, Point p) {//直线l和点p的距离
	return abs(cross(l.p2 - l.p1, p - l.p1) / (l.p2 - l.p1).abs() );
}

double getDistanceSP(Segment s, Point p) {//线段s与点p的距离
	if( dot(s.p2 - s.p1, p - s.p1) < 0.0 ) return (p - s.p1).abs();
	if( dot(s.p1 - s.p2, p - s.p2) < 0.0 ) return (p - s.p2).abs();
	return getDistanceLP(s, p);
}

/*ccw模块*/ 
static const int COUNTER_CLOCKWISE = 1;
static const int CLOCKWISE = -1;
static const int ONLINE_BACK = 2;
static const int ONLINE_FRONT = -2;
static const int ON_SEGMENT = 0;

int ccw(Point p0, Point p1, Point p2) {
	Vector a = p1 - p0;
	Vector b = p2 - p0;
	if( cross(a, b) > EPS ) return COUNTER_CLOCKWISE;
	if( cross(a, b) < -EPS ) return CLOCKWISE;
	if( dot(a, b) < -EPS ) return ONLINE_BACK;
	if( a.norm() < b.norm() ) return ONLINE_FRONT;

	return ON_SEGMENT;
}

bool intersect(Point p1, Point p2, Point p3, Point p4) {//判断线段p1p2和线段p3p4是否相交
	return ( ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 &&
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0 );
}

bool intersect(Segment s1, Segment s2) {//判断线段s1和线段s2是否相交
	return intersect(s1.p1, s1.p2, s2.p1, s2.p2);
}

double getDistance(Segment s1, Segment s2) {//线段与线段的距离
	if(intersect(s1, s2) ) return 0.0;
	return min(min(getDistanceSP(s1, s2.p1), getDistanceSP(s1, s2.p2)), 
				min(getDistanceSP(s2, s1.p1), getDistanceSP(s2, s1.p2)));
}

int main(){
	int q;
	cin>>q;
	
	Segment s1, s2;
	
	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;
		printf("%.10f\n", getDistance(s1, s2) );
	}
}

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