计算几何摸黑(3)-圆及有关计算(上)

计算几何摸黑(3)-圆及有关计算(上)

圆的表示方法

​ 和直线一样，圆也可以使用参数式表示。显而易见，一个唯一的圆可以用一个点表示圆心，和一个实数表示半径。所以，我们可以这样定义一个圆。

class circle_base
{
public:
    double r;
	Point c;
	Circle(Point C = Point(0.0, 0.0), R = 0.0)
		: c(C), r(R) { }
}C;

​ 约定：这里用C.c表示圆C的圆心，C.r表示圆C的半径。

通过圆心角求点的坐标

​ 我们作一条过圆心，平行于\(x\)轴的直线\(lx\)，再作一条过圆心，平行于\(y\)轴的直线\(ly\)，那么这又构成了一个坐标系。这样，在圆上且和\(lx\)正半轴有一个夹角的点唯一。我们把那个角称为圆心角。所以，我们可以用三角函数和圆的成分来实现这个功能。注意：在这里我用到了类继承，也可以直接将Circle类声明中添加point函数。

class Circle : public circle_base
{
public:
	Circle(circle_base &t)
		: circle_base(t) { }
	inline Point point(const double theta)
	{
		return Point(x + std::cos(theta) * r, y + std::sin(theta) * r);
	}
}

直线和圆的交点

​ 我们可以通过解方程法求出直线和圆的交点。 若直线为\(AB\)，我们可以设交点为\(P=A+t(B-A)\)。这里要用到高中课本上的圆方程公式：

\[(x - x_0) ^ 2 + (y - y_0) ^ 2 = r ^ 2 \]

​ 由上面的设出的式子，我们可以这样推导：

\(P = (ax, ay) + t(bx - ax, by - ay)\)

$ = (ax - tax + tbx, ay - tay + tby)$

\(\because P \in \odot C\)

\(\therefore (ax - tax + tbx - C.c.x)^2 +(ay - tay + tby - C.c.y)^2=r^2\)

\(\because bx - ax = \vec {AB}.x, by - ay = \vec {AB}.y\)

\(整理得，(t(\vec {AB}.x) + (ax - C.c.x))^2 + (t(\vec {AB}.y) + (ay - C.c.y))^2 - r^2 = 0\)

\(令a = (\vec {AB}.x), b = (ax - C.c.x), c = (\vec{AB}.y), d = (ay - C.c.y)\)

\(\therefore (at+b)^2+(ct+d)^2 - r^2 = 0\)

\((a^2+c^2)t^2+(ab+cd)t+b^2+d^2-r^2=0\)

\(令e = (a^2 + c^2), f = (ab + cd), g = (b^2 + d^2 - r^2)\)

\(\therefore et^2 + ft + g = 0\)

\(令\Delta= f^2 - 4eg\)

​ 易知， 如果\(\Delta > 0\)，那么直线将与圆有两个交点，\(\Delta < 0\)，直线将与圆有相离，\(\Delta = 0\)，直线将与圆有一个交点。而把这个一元二次方程的解带入，并用三角函数计算就可以求出角了。实现如下。

int getLineCircInter(const Line &L, const Circle &C,
	double &t1, double &t2,
	double &s1, double &s2)
	// 意义: L交C于圆心角为t1和t2的点，分别为s1, s2
{
	double a = L.v.x, b = L.p.x - C.c.x,
		   c = L.v.y, d = L.p.y - C.c.y,
		   e = a * a + c * c, f = 2 * (a * b + c * d),
		   g = (b * b + d * d - C.r * C.r); // 详细推导过程见上文
	double delta = f * f - 4 * a * c; // 判别式
	if (dcmp(delta) < 0)
		return 0;
	if (dcmp(delta) == 0)
	{
		t1 = t2 = -(f / (2 * e));
		s1 = s2 = C.point(t1);
		return 1;
	}
	t1 = (-f + std::sqrt(delta)) / (2 * e);
	t2 = (-f - std::sqrf(delta)) / (2 * e);
	s1 = C.point(t1);
	s2 = C.point(t2);
	return 2;
} // 注意，这里的返回值的含义是交点个数。