计算几何学 | 圆与直线的交点 | Cross Points of a Circle and a Line | C/C++实现

问题描述

求圆c与直线$l$的交点。

输入：
输入按照下述格式给出：
$cx$ $cy$ $r$
$q$
$Lin{e}_1$
$Lin{e}_2$
…
$Lin{e}_q$
第1行输入圆心坐标cx，cy以及半径r。第2行输入问题数q。
接下来q行按照下述格式输入q个直线$Lin{e}_i$作为问题。
$x_1$ $y_1$ $x_2$ $y_2$
各直线由其通过的2个p1、p2表示，$x_1$、$y_1$是p1的坐标，$x_2$、$y_2$是p2的坐标。以上输入均为整数。
输出：
根据各个问题输出交点坐标
请按照下述规则输出2个交点的坐标，相邻数值用空格隔开：
只有1个交点时输出2个相同的坐标
先输出x坐标较小的点。x坐标相同时先输出y坐标较小的点。
允许误差不超过0.000001
限制：
p1、p2不是同一个点
圆与直线必然存在交点
1 ≤ q ≤ 1000
-10000 ≤ $cx,cy,x_1,y_1,x_2,y_2$ ≤ 10000
1 ≤ r ≤ 10000

输入示例

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

输出示例

1.000000 1.000000 3.000000 1.000000
3.000000 1.000000 3.000000 1.000000

讲解

首先我们来考虑题目限制中“圆与直线存在交点”这一条件。实际上，只有检查圆心与直线的距离是否大于r（关于距离的问题请参见：计算几何学 | 距离 | Distance | C/C++实现），我们就能判断圆与直线是否存在交点。

求圆与直线的方法有许多方法，这里我们用向量来求解

首先求圆心c在直线$l$上的投影点pr。接下来，求直线$l$上的单位向量（大小为1的向量）e。如下所示，e可用向量除以其大小求出：
$e=(p2-p1)/|p2-p1|$
再接下来，根据半径r和向量pr的长度计算出圆内线段部分的一半base。然后把刚才求出的单位向量e乘以大小base，即可得出直线$l$上与base同样大小（长度）的向量。最后，以投影点pr为起点，向正/负方向加上该向量，我们就得到圆与直线的交点（的坐标）了。
求圆c与直线$l$交点的程序如下所示。

圆c与直线$l$的交点：

pair<Point, Point> getCrossPoints(Circle c, Line l) {
	assert(intersect(c, l));
	Vector pr = project(l, c.c);
	Vector e = (l.p2 - l.p1) / abs(l.p2 - l.p1);
	double base = sqrt(c.r * c.r - norm(pr - c.c));
	return make_pair(pr + e * base, pr - e * base);
}

AC代码如下

#include<stdio.h>
#include<iostream>
#include<cmath>
#include<assert.h>
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;

class Circle {//Circle 圆 
	public:
		Point c;
		double r;
		Circle(Point c = Point(), double r = 0.0): c(c), r(r) {}
};

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 project(Segment s, Point p) {//投影 对于给定的三个点p1、p2、p，从点p向通过
//p1、p2的直线引一条垂线，求垂足x的坐标。（点p在直线p1p2上的投影） 
	Vector base = s.p2 - s.p1;
	double r = dot(p - s.p1, base) / base.norm();
	return s.p1 + base * r;
} 

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

bool intersect(Circle c, Line l) {
	if(getDistanceLP(l, c.c) > c.r) {
		return false;
	} else {
		return true;
	}
}
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;
}

pair<Point, Point> getCrossPoints(Circle c, Line l) {
	assert(intersect(c, l));
	Vector pr = project(l, c.c);
	Vector e = (l.p2 - l.p1) / (l.p2 - l.p1).abs();
	double base = sqrt(c.r * c.r - (pr - c.c).norm() );
	return make_pair(pr + e * base, pr - e * base);
}

int main(){
	Circle c;
	cin>>c.c.x>>c.c.y>>c.r;
	
	int q;
	cin>>q;
	
	Line l;
	pair<Point, Point> p;
	
	while(q--){
		cin>>l.p1.x>>l.p1.y>>l.p2.x>>l.p2.y;
		p = getCrossPoints(c, l);
		if(p.first.x < p.second.x || (p.first.x == p.second.x && p.first.y <= p.second.y) ) {
			printf("%.6f %.6f %.6f %.6f\n", p.first.x, p.first.y, p.second.x, p.second.y);
		} else {
			printf("%.6f %.6f %.6f %.6f\n", p.second.x, p.second.y, p.first.x, p.first.y);
		}
	}
} 

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