算法第一次考试 计算几何题解
题目描述
给定平面内一个圆(圆心为C,半径为R)和严格在圆外的两个点A, B,现要用一根曲线段连接两点,且其不可贯穿该圆。求该曲线段的最短长度。
输入
多组,每行7个浮点数,分别表示C的横纵坐标、R、A的横纵坐标、B的横纵坐标。
输出
对于每组输入,输出一行一个浮点数,表示曲线段的最短长度。进行spj,当绝对误差达到精度要求时你的答案即可AC。
输入样例
3 0 2 -1 0 6 0
3 0 2 -1 10 6 10
输出样例
9.88810282661858153119
7.00000000000000000000
解题思路
本题考察了计算几何基础。
首先可根据线段AB和圆C是否相交分为两种情况:
- 若不相交,答案显然即为线段AB长度。此处考察了判断圆和线段的位置关系,即求点到线段的距离。
- 若相交,考虑寻找最短长度的过程。假设用一根松弛的软绳连接AB,然后渐渐用力拉扯,可以预见的是,当拉扯到最紧时,曲线段达到最短,且其由两条切线段和一段圆弧组成。此处考察了求切点和求视角。
因此答案就是圆弧长+两条切线段长。为编程方便,我们只需暴力地考虑切点组合的4种组合情况即可,即求出$\min{\overset{\frown}{P_aP_b}, \overset{\frown}{P_aQ_b}, \overset{\frown}{Q_aP_b}, \overset{\frown}{Q_aQ_b}}$即为所选弧长:
至于求切点和求视角的细节:
-
对于求切点,可以使用直观的几何做法:考虑 $A$ 的切线,记两切点为 $P_a, Q_a$,先用反三角函数求出 $\theta = \angle CAP_a ($亦等于 $\angle CAQ_a)$,然后将 $\overrightarrow{AC}$ 以 $A$ 为心分别顺逆时针旋转 $\theta$,再缩放该向量长度为 $|AP_a|($亦等于 $|AQ_a|)$,则得到 $\overrightarrow{AP_a}$ 和 $\overrightarrow{AQ_a}$,再分别加上点 $A$ 坐标即得到 $P_a$ 和 $Q_a$ 的坐标。
-
对于求视角:求视角是为了求四个圆心角,进而求得最短的圆弧。利用叉积和反正切即可轻松求得视角。
因此本题即可得解。时间复杂度当然是 $O(n)$ 的啦。
代码
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using std::min;
#define IL inline
#define OP operator
#define CON constexpr
typedef double flot;
CON flot EPS = 1e-5, S_EPS = EPS / 16, INF = 1e22, PI = acos(-1), RAD = 180 / PI;
IL CON flot ABS(flot x) {return x<-EPS ? -x : (x>EPS ? x:0);}
IL CON int sign(flot x) {return x<-EPS ? -1 : x>EPS;}
struct Po
{
flot x, y;
Po(void) { }
Po(flot xx, flot yy) : x(xx), y(yy) { }
Po(const Po &a, const Po &b) : x(b.x-a.x), y(b.y-a.y) { }
Po OP +(const Po &o) const {return Po(x+o.x, y+o.y);}
Po OP -(const Po &o) const {return Po(x-o.x, y-o.y);}
flot OP *(const Po &o) const {return x*o.x + y*o.y;}
flot OP ^(const Po &o) const {return x*o.y - y*o.x;}
flot OP |(const Po &o) const {return sqrt((x-o.x)*(x-o.x)+ (y-o.y)*(y-o.y));}
flot len() const {return sqrt(x*x + y*y);}
Po zoom(flot r) const
{
flot l = len();
if (!sign(l))
return Po(0, 0);
r /= l;
return Po(r*x, r*y);
}
Po rot(const Po &o, flot a) const // 绕o旋转a,a是弧度制!
{
Po v = *this - o;
flot c = cos(a), s = sin(a);
return Po(o.x+ c*v.x - s*v.y, o.y+ s*v.x + c*v.y);
}
flot th(const Po &a, const Po &b) const
{
Po p = *this; // 在this观察ab的视角
return ABS(atan2( ABS((a-p)^(b-p)),(a-p)*(b-p) ));
}
};
flot po_d_seg(Po c, Po s, Po t)
{
Po ST(t-s), SP(c-s), TP(c-t);
if (sign(ST * SP) < 0) // c在反向延长线一侧
return SP.len();
if (sign(ST * TP) > 0) // c在延长线一侧
return TP.len();
return ABS((ST ^ SP) / ST.len()); // c在线段正两侧
}
void calc_tan(Po c, flot r, Po a, Po &pa, Po &qa)
{
flot dca = c | a;
flot alp = acos(r / dca);
pa = a + c.rot(a, alp).zoom(sqrt(dca * dca - r * r));
qa = a + c.rot(a, -alp).zoom(sqrt(dca * dca - r * r));
}
int main()
{
flot r;
Po c, a, b;
while (~scanf("%lf %lf %lf %lf %lf %lf %lf", &c.x, &c.y, &r, &a.x, &a.y, &b.x, &b.y))
{
Po pa, qa, pb, qb;
if (sign(po_d_seg(c, a, b) - r) >= 0)
{
printf("%.20f\n", a | b);
continue;
}
calc_tan(c, r, a, pa, qa);
calc_tan(c, r, b, pb, qb);
flot min_th = min(c.th(pa, pb), c.th(pa, qb));
min_th = min(min_th, c.th(qa, pb));
min_th = min(min_th, c.th(qa, qb));
printf("%.20f\n", min_th*r + (a | pa) + (b | pb));
}
return 0;
}