请用C++写出根据点集求NURBS控制点的代码

这是一个比较复杂的问题，需要考虑到很多细节。以下是一个简单的伪代码，用于描述计算NURBS控制点的基本流程：

// 假设已经有了点集和权重数组
vector<Point> points;
vector<double> weights;

// 初始化控制点数组
vector<Point> control_points;

// 设置NURBS的阶数和节点向量
int degree = 3; // 阶数为3
vector<double> knots; // 节点向量

// 计算节点向量
int n = points.size() + degree + 1;
knots.resize(n);
for (int i = 0; i < n; i++) {
    if (i < degree) {
        knots[i] = 0;
    } else if (i >= n - degree) {
        knots[i] = 1;
    } else {
        knots[i] = (double)(i - degree + 1) / (double)(n - 2 * degree);
    }
}

// 构建矩阵A和向量b
int m = points.size();
int k = degree + 1;
Matrix A(m, k + m);
Vector b(m);

for (int i = 0; i < m; i++) {
    double w = weights[i];
    Point p = points[i];

    for (int j = 0; j < k; j++) {
        A(i, j) = N(j, degree, knots, u) * w;
    }

    for (int j = 0; j < m; j++) {
        A(i, k + j) = 0;
    }

    b[i] = p * w;
}

// 解线性方程组
Vector x = A.solve(b);

// 将解向量转换为控制点
for (int i = 0; i < m; i++) {
    control_points.push_back(Point(x[i + k], x[i] / x[i + k]));
}

其中，`N(j, degree, knots, u)`是计算NURBS基函数的函数。它的实现可以参考以下代码：

double N(int i, int p, const vector<double>& U, double u) {
    if (p == 0) {
        if (u >= U[i] && u < U[i + 1]) {
            return 1;
        } else {
            return 0;
        }
    }

    double a = 0;
    double b = 0;

    if (U[i + p] - U[i] != 0) {
        a = (u - U[i]) / (U[i + p] - U[i]);
    }

    if (U[i + p + 1] - U[i + 1] != 0) {
        b = (U[i + p + 1] - u) / (U[i + p + 1] - U[i + 1]);
    }

    return a * N(i, p - 1, U, u) + b * N(i + 1, p - 1, U, u);
}

请注意，这只是一个简单的伪代码，实际实现中还需要考虑很多细节，如矩阵求解的算法、NURBS基函数的优化等等。