c++ 点群 相切连接 源代码

以下是C++实现的点群相切连接的源代码：

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

const double eps = 1e-8;
const double inf = 1e20;

// 二维向量
struct Vector {
    double x, y;
    Vector() {}
    Vector(double x, double y) : x(x), y(y) {}
    Vector operator + (const Vector &v) const { return Vector(x + v.x, y + v.y); }
    Vector operator - (const Vector &v) const { return Vector(x - v.x, y - v.y); }
    Vector operator * (double d) const { return Vector(x * d, y * d); }
    double dot(const Vector &v) const { return x * v.x + y * v.y; } // 向量点积
    double cross(const Vector &v) const { return x * v.y - y * v.x; } // 向量叉积
    double length() const { return sqrt(x * x + y * y); } // 向量长度
    Vector unit() const { return *this * (1 / length()); } // 单位向量
};

// 圆
struct Circle {
    Vector c; // 圆心
    double r; // 半径
    Circle() {}
    Circle(Vector c, double r) : c(c), r(r) {}
    Vector point(double a) const { return Vector(c.x + cos(a) * r, c.y + sin(a) * r); } // 圆上某点
};

// 判定两个圆是否相交
int getCircleCircleIntersection(const Circle& c1, const Circle& c2, vector<Vector>& sol) {
    double d = (c1.c - c2.c).length();
    if (d > c1.r + c2.r + eps || d < fabs(c1.r - c2.r) - eps) return 0;
    double a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d));
    Vector p1 = c1.point(a), p2 = c1.point(-a);
    sol.push_back(p1);
    if (p1.x == p2.x && p1.y == p2.y) return 1;
    sol.push_back(p2);
    return 2;
}

// 判定三个圆是否相交
int getCircleCircleIntersection(const Circle& c1, const Circle& c2, const Circle& c3, vector<Vector>& sol) {
    vector<Vector> p12;
    if (getCircleCircleIntersection(c1, c2, p12) == 0) return 0;
    vector<Vector> p13;
    if (getCircleCircleIntersection(c1, c3, p13) == 0) return 0;
    vector<Vector> p23;
    if (getCircleCircleIntersection(c2, c3, p23) == 0) return 0;
    for (int i = 0; i < p12.size(); i++) {
        for (int j = 0; j < p13.size(); j++) {
            for (int k = 0; k < p23.size(); k++) {
                if ((p12[i] - p13[j]).length() < eps && (p12[i] - p23[k]).length() < eps) {
                    sol.push_back(p12[i]);
                    return 1;
                }
            }
        }
    }
    return 0;
}

// 计算两圆相切时的切点
int getTangents(const Circle& a, const Circle& b, Vector* v) {
    int cnt = 0;
    if (a.r < b.r) { swap(a, b); }
    double d2 = (a.c - b.c).lengthSquared();
    double rdif = a.r - b.r, rsum = a.r + b.r;
    if (d2 < rdif * rdif - eps) return 0; // 内含
    double base = atan2(b.c.y - a.c.y, b.c.x - a.c.x);
    if (d2 < rsum * rsum + eps) { // 相交
        double ang = acos((a.r - b.r) / sqrt(d2));
        v[cnt++] = a.point(base + ang);
        if (fabs(ang) > eps) v[cnt++] = a.point(base - ang);
        v[cnt++] = b.point(base + ang);
        if (fabs(ang) > eps) v[cnt++] = b.point(base - ang);
    } else { // 相离
        double ang = acos((a.r + b.r) / sqrt(d2));
        v[cnt++] = a.point(base + ang);
        v[cnt++] = a.point(base - ang);
        v[cnt++] = b.point(base + ang + M_PI);
        v[cnt++] = b.point(base - ang + M_PI);
    }
    return cnt;
}

// 计算三圆相切时的切线
int getTangents(Circle a, Circle b, Circle c, Vector* v) {
    if (b.r > c.r) { swap(b, c); }
    if (a.r > c.r) { swap(a, c); }
    if (a.r > b.r) { swap(a, b); }
    if (fabs((a.c - b.c).length() - a.r - b.r) < eps) { // 两圆相切
        Vector u[2];
        int cnt = getTangents(a, b, u);
        for (int i = 0; i < cnt; i++) {
            if (fabs((u[i] - c.c).length() - c.r) < eps) {
                v[0] = u[i];
                return 1;
            }
        }
        return 0;
    } else if (fabs((a.c - b.c).length() + a.r + b.r) < eps) { // 两圆相离
        Vector u[4], t[2];
        int cnt = getTangents(a, b, u);
        int k = 0;
        for (int i = 0; i < cnt; i++) {
            if (getCircleCircleIntersection(c, Circle(u[i], (u[i] - a.c).length())) > 0) {
                t[k++] = u[i];
            }
            if (getCircleCircleIntersection(c, Circle(u[i], (u[i] - b.c).length())) > 0) {
                t[k++] = u[i];
            }
        }
        if (k < 2) return 0;
        double ang = acos((c.r / (t[1] - c.c).length()));
        v[0] = t[1];
        v[1] = c.point(atan2(t[1].y - c.c.y, t[1].x - c.c.x) + ang);
        v[2] = c.point(atan2(t[1].y - c.c.y, t[1].x - c.c.x) - ang);
        if (k > 2) {
            ang = acos((c.r / (t[0] - c.c).length()));
            v[3] = t[0];
            v[4] = c.point(atan2(t[0].y - c.c.y, t[0].x - c.c.x) + ang);
            v[5] = c.point(atan2(t[0].y - c.c.y, t[0].x - c.c.x) - ang);
            return 6;
        }
        return 3;
    } else { // 两圆相交
        Vector u[4], t[2];
        int cnt = getTangents(a, b, u);
        int k = 0;
        for (int i = 0; i < cnt; i++) {
            if (getCircleCircleIntersection(c, Circle(u[i], (u[i] - a.c).length())) > 0) {
                t[k++] = u[i];
            }
            if (getCircleCircleIntersection(c, Circle(u[i], (u[i] - b.c).length())) > 0) {
                t[k++] = u[i];
            }
        }
        if (k < 2) return 0;
        if (k == 2) {
            v[0] = t[0];
            v[1] = t[1];
            return 2;
        }
        double ang = acos((c.r / (t[1] - c.c).length()));
        v[0] = t[1];
        v[1] = c.point(atan2(t[1].y - c.c.y, t[1].x - c.c.x) + ang);
        v[2] = c.point(atan2(t[1].y - c.c.y, t[1].x - c.c.x) - ang);
        ang = acos((c.r / (t[0] - c.c).length()));
        v[3] = t[0];
        v[4] = c.point(atan2(t[0].y - c.c.y, t[0].x - c.c.x) + ang);
        v[5] = c.point(atan2(t[0].y - c.c.y, t[0].x - c.c.x) - ang);
        return 6;
    }
}

int main() {
    int T;
    cin >> T;
    while (T--) {
        Circle a, b, c;
        cin >> a.c.x >> a.c.y >> a.r;
        cin >> b.c.x >> b.c.y >> b.r;
        cin >> c.c.x >> c.c.y >> c.r;
        Vector v[6];
        int cnt = getTangents(a, b, c, v);
        if (cnt == 0) {
            cout << "NO\n";
        } else if (cnt == 2) {
            cout << "YES\n";
            printf("%.6lf %.6lf\n", v[0].x, v[0].y);
            printf("%.6lf %.6lf\n", v[1].x, v[1].y);
        } else {
            cout << "YES\n";
            printf("%.6lf %.6lf\n", v[0].x, v[0].y);
            printf("%.6lf %.6lf\n", v[2].x, v[2].y);
            printf("%.6lf %.6lf\n", v[4].x, v[4].y);
        }
    }
    return 0;
}

该程序实现了计算三个圆的相切连接，并输出连接的切点坐标。