贝塞尔（Bezier）曲线与曲面详解

一、贝塞尔曲线基础

1.1 数学定义

先照搬教科书中的定义：

贝塞尔曲线由一组控制点定义，n次贝塞尔曲线的参数方程为：

其中PiP_iPi​是控制点坐标， Bi,n(t)B_{i,n}(t)Bi,n​(t) 是伯恩斯坦基函数：

其中CniC_n^iCni​ 为[t+(1−t)]n[t+(1-t)]^n[t+(1−t)]n二项展开式的二项式系数；

1.2 德卡斯特里奥算法

德卡斯特里奥算法是计算贝塞尔曲线点的高效递归算法：

#include <iostream>
#include <vector>
#include <cmath>
#include <memory>

// 二维点或二维向量
struct Point2D 
{
    double x, y;
    Point2D(double x = 0, double y = 0) : x(x), y(y) {}
    
    Point2D operator+(const Point2D& other) const {
        return Point2D(x + other.x, y + other.y);
    }
    
    Point2D operator*(double scalar) const {
        return Point2D(x * scalar, y * scalar);
    }
};

// 三维点或三维向量
struct Point3D 
{
    double x, y, z;
    Point3D(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
    
    Point3D operator+(const Point3D& other) const {
        return Point3D(x + other.x, y + other.y, z + other.z);
    }
    
    Point3D operator*(double scalar) const {
        return Point3D(x * scalar, y * scalar, z * scalar);
    }
};

class BezierCurve {
private:
    std::vector<Point2D> control_points_;
    
    // 计算二项式系数 C(n, k)
    int binomialCoefficient(int n, int k) const {
        if (k == 0 || k == n) return 1;
        if (k > n - k) k = n - k;
        
        int result = 1;
        for (int i = 1; i <= k; ++i) {
            result = result * (n - i + 1) / i;
        }
        return result;
    }
    
    // 伯恩斯坦基函数
    double bernstein(int n, int i, double t) const {
        return binomialCoefficient(n, i) * 
               std::pow(t, i) * 
               std::pow(1 - t, n - i);
    }

public:
    BezierCurve(const std::vector<Point2D>& control_points) 
        : control_points_(control_points) {}
    
    // 计算曲线上参数 t 对应的点
    Point2D evaluate(double t) const {
        if (control_points_.empty()) {
            return Point2D();
        }
        
        int n = control_points_.size() - 1;
        Point2D result(0, 0);
        
        for (int i = 0; i <= n; ++i) {
            double basis = bernstein(n, i, t);
            result = result + control_points_[i] * basis;
        }
        
        return result;
    }
    
    // 离散化曲线 - 生成采样点
    std::vector<Point2D> discretize(int num_samples = 100) const {
        std::vector<Point2D> samples;
        samples.reserve(num_samples);
        
        for (int i = 0; i < num_samples; ++i) {
            double t = static_cast<double>(i) / (num_samples - 1);
            samples.push_back(evaluate(t));
        }
        
        return samples;
    }
    
    // 德卡斯特里奥算法 - 更稳定的计算方法
    Point2D evaluateDeCasteljau(double t) const {
        if (control_points_.empty()) {
            return Point2D();
        }
        
        std::vector<Point2D> points = control_points_;
        int n = points.size();
        
        for (int k = 1; k < n; ++k) {
            for (int i = 0; i < n - k; ++i) {
                points[i] = points[i] * (1 - t) + points[i + 1] * t;
            }
        }
        
        return points[0];
    }
    
    // 获取控制点
    const std::vector<Point2D>& getControlPoints() const {
        return control_points_;
    }
    
    // 设置控制点
    void setControlPoints(const std::vector<Point2D>& points) {
        control_points_ = points;
    }
};

二、贝塞尔曲面

2.1 数学定义

贝塞尔曲面是贝塞尔曲线在二维参数空间的扩展，由控制点网格定义：

class BezierSurface {
private:
    std::vector<std::vector<Point3D>> control_net_;
    int u_degree_, v_degree_;
    
    // 伯恩斯坦基函数
    double bernstein(int n, int i, double t) const {
        if (n == 0) return 1.0;
        
        int coeff = 1;
        for (int k = 1; k <= i; ++k) {
            coeff = coeff * (n - k + 1) / k;
        }
        
        return coeff * std::pow(t, i) * std::pow(1 - t, n - i);
    }

public:
    BezierSurface(const std::vector<std::vector<Point3D>>& control_net)
        : control_net_(control_net) {
        u_degree_ = control_net.size() - 1;
        v_degree_ = control_net[0].size() - 1;
    }
    
    // 计算曲面上参数 (u,v) 对应的点
    Point3D evaluate(double u, double v) const {
        Point3D result(0, 0, 0);
        
        for (int i = 0; i <= u_degree_; ++i) {
            double basis_u = bernstein(u_degree_, i, u);
            
            for (int j = 0; j <= v_degree_; ++j) {
                double basis_v = bernstein(v_degree_, j, v);
                result = result + control_net_[i][j] * (basis_u * basis_v);
            }
        }
        
        return result;
    }
    
    // 离散化曲面 - 生成网格点
    std::vector<std::vector<Point3D>> discretize(int u_samples = 20, int v_samples = 20) const {
        std::vector<std::vector<Point3D>> mesh;
        mesh.resize(u_samples);
        
        for (int i = 0; i < u_samples; ++i) {
            double u = static_cast<double>(i) / (u_samples - 1);
            mesh[i].resize(v_samples);
            
            for (int j = 0; j < v_samples; ++j) {
                double v = static_cast<double>(j) / (v_samples - 1);
                mesh[i][j] = evaluate(u, v);
            }
        }
        
        return mesh;
    }
    
    // 德卡斯特里奥算法计算曲面点
    Point3D evaluateDeCasteljau(double u, double v) const {
        // 首先在 u 方向应用德卡斯特里奥
        std::vector<Point3D> temp_points;
        temp_points.reserve(v_degree_ + 1);
        
        for (int j = 0; j <= v_degree_; ++j) {
            // 对每一列应用德卡斯特里奥算法
            std::vector<Point3D> column;
            for (int i = 0; i <= u_degree_; ++i) {
                column.push_back(control_net_[i][j]);
            }
            
            // 在 u 方向计算
            for (int k = 1; k <= u_degree_; ++k) {
                for (int i = 0; i < u_degree_ - k + 1; ++i) {
                    column[i] = column[i] * (1 - u) + column[i + 1] * u;
                }
            }
            temp_points.push_back(column[0]);
        }
        
        // 然后在 v 方向应用德卡斯特里奥
        for (int k = 1; k <= v_degree_; ++k) {
            for (int j = 0; j < v_degree_ - k + 1; ++j) {
                temp_points[j] = temp_points[j] * (1 - v) + temp_points[j + 1] * v;
            }
        }
        
        return temp_points[0];
    }
};

三、完整的建模管理

#include <fstream>
#include <sstream>

class BezierModeler {
private:
    std::vector<std::shared_ptr<BezierCurve>> curves_;
    std::vector<std::shared_ptr<BezierSurface>> surfaces_;
    
public:
    // 创建贝塞尔曲线
    std::shared_ptr<BezierCurve> createCurve(const std::vector<Point2D>& control_points) {
        auto curve = std::make_shared<BezierCurve>(control_points);
        curves_.push_back(curve);
        return curve;
    }
    
    // 创建贝塞尔曲面
    std::shared_ptr<BezierSurface> createSurface(const std::vector<std::vector<Point3D>>& control_net) {
        auto surface = std::make_shared<BezierSurface>(control_net);
        surfaces_.push_back(surface);
        return surface;
    }
    
    // 导出曲线为CSV格式（用于可视化）
    void exportCurveToCSV(const BezierCurve& curve, const std::string& filename, int samples = 100) {
        std::ofstream file(filename);
        file << "x,y" << std::endl;
        
        auto discrete_points = curve.discretize(samples);
        for (const auto& point : discrete_points) {
            file << point.x << "," << point.y << std::endl;
        }
        
        // 同时导出控制多边形
        std::ofstream ctrl_file("control_" + filename);
        ctrl_file << "x,y" << std::endl;
        for (const auto& point : curve.getControlPoints()) {
            ctrl_file << point.x << "," << point.y << std::endl;
        }
    }
    
    // 导出曲面为OBJ格式
    void exportSurfaceToOBJ(const BezierSurface& surface, const std::string& filename, 
                           int u_samples = 20, int v_samples = 20) {
        std::ofstream file(filename);
        
        auto mesh = surface.discretize(u_samples, v_samples);
        
        // 写入顶点
        for (const auto& row : mesh) {
            for (const auto& vertex : row) {
                file << "v " << vertex.x << " " << vertex.y << " " << vertex.z << std::endl;
            }
        }
        
        // 写入面
        for (int i = 0; i < u_samples - 1; ++i) {
            for (int j = 0; j < v_samples - 1; ++j) {
                int v1 = i * v_samples + j + 1;
                int v2 = i * v_samples + j + 2;
                int v3 = (i + 1) * v_samples + j + 2;
                int v4 = (i + 1) * v_samples + j + 1;
                
                file << "f " << v1 << " " << v2 << " " << v3 << " " << v4 << std::endl;
            }
        }
    }
    
    // 曲线求导 - 返回导数的控制点
    std::vector<Point2D> differentiateCurve(const BezierCurve& curve) {
        const auto& points = curve.getControlPoints();
        int n = points.size() - 1;
        std::vector<Point2D> derivative_points;
        
        for (int i = 0; i < n; ++i) {
            Point2D diff = (points[i + 1] + points[i] * (-1)) * n;
            derivative_points.push_back(diff);
        }
        
        return derivative_points;
    }
};

// 示例和使用演示
int main() {
    BezierModeler modeler;
    
    // 创建二次贝塞尔曲线
    std::vector<Point2D> control_points = {
        Point2D(0, 0),
        Point2D(1, 2), 
        Point2D(3, 1)
    };
    
    auto curve = modeler.createCurve(control_points);
    
    // 测试曲线求值
    Point2D p = curve->evaluate(0.5);
    std::cout << "曲线在 t=0.5 的点: (" << p.x << ", " << p.y << ")" << std::endl;
    
    // 离散化并导出
    modeler.exportCurveToCSV(*curve, "bezier_curve.csv");
    
    // 创建贝塞尔曲面示例
    std::vector<std::vector<Point3D>> control_net = {
        { Point3D(0,0,0), Point3D(0,1,1), Point3D(0,2,0) },
        { Point3D(1,0,1), Point3D(1,1,2), Point3D(1,2,1) },
        { Point3D(2,0,0), Point3D(2,1,1), Point3D(2,2,0) }
    };
    
    auto surface = modeler.createSurface(control_net);
    
    // 测试曲面求值
    Point3D sp = surface->evaluate(0.5, 0.5);
    std::cout << "曲面在 (0.5,0.5) 的点: (" << sp.x << ", " << sp.y << ", " << sp.z << ")" << std::endl;
    
    // 导出曲面
    modeler.exportSurfaceToOBJ(*surface, "bezier_surface.obj");
    
    std::cout << "贝塞尔曲线曲面建模完成！" << std::endl;
    
    return 0;
}

四、关键特性说明

4.1 贝塞尔曲线的性质

4.2 应用场景

五、进阶

在实际几何内核开发中，曲线曲面类的实现还需要考虑以下几点：