# Mini-Spline

**Repository Path**: null_552_0547/mini-spline

## Basic Information

- **Project Name**: Mini-Spline
- **Description**: Mini-Spline 通过 Qt 下的 OpenGL 编程实现样条曲线曲面生成、数据拟合、点云参数化、正交投影等功能。其中整理了一个小型的数值计算库，一个半边结构的网格数据类型，并实现网格模型的基本处理方法。
- **Primary Language**: C++
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 1
- **Created**: 2025-06-01
- **Last Updated**: 2025-06-01

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# Mini-Spline

#### 介绍
Mini-Spline 通过 Qt 下的 OpenGL 编程实现样条曲线曲面生成、数据拟合、点云参数化、正交投影等功能。其中整理了一个小型的数值计算库，一个半边结构的网格数据类型，并实现网格模型的基本处理方法。



#### 开发进度

3.9 

* 整理线性代数库，修改了断言的方式为标准 assert；
* 使用 catch2 开源库进行数值测试
    * 完成向量操作测试
    * 完成矩阵基本操作测试
    * 对基本矩阵求解器进行测试



3.10

* 完成隐式 QR 方法求解器的调试；
* 搭建窗口程序；
* 建立抽象类 MiniObject 以及 MiniCloud, MiniCurve 子类；
* 完成一般着色器配置，可以绘制 MiniCloud；



3.11

* 优化片段着色器，可以修改颜色 UserColor；
* 顶点着色器添加视图、模型、投影矩阵；
* 之后合并了模型视图矩阵，去除了单独的模型矩阵；
* 完成相机类的设计 MiniCamera；
* 完成鼠标滚动调整距离，鼠标点击拖拽调整视角；
* 完成基本绘制框架：坐标轴 drawAxis，物体绘制流程 drawObject；
* 将半边网格整合为 MiniMesh 类型；
* 重新调整了 point3f, vector3f, rgb_space 结构，整合为 Vector3f 统一表示三维点/向量/RGB 空间类型；
* 整合 Quaternion 四元数类型；
* 去除了包围盒和 QEM 简化算法；
* 完成兔子模型的读取显示；
* **注意到 SpLinearSolver 的 Jacobi 迭代没有用到权 omega；暂时搁置稀疏矩阵类型和求解器；**
* 预计使用 Catch2 检测 Vector3f, Quaternion 类型的算法，还未考虑如何测试 MiniMesh 半边网格类型；



3.12

* 移除增量凸包和平面凸包算法；
* 移除 MiniMesh 中的颜色列表；
* 修改了网格材质结构；
* 完成 Phong 光照着色器：包括漫反射、镜面反射、环境光、距离衰减效果；



3.13

* 将 Test 测试部分整合进来，调整 makefile 文件，使得测试程序可以直接调用 Qt 程序相关的头文件和源文件；
* 测试 Vector3f 类型及相关算法；



3.14

* 测试 Quaternion 类型及相关算法；
* 调整了光照系数，将其与渲染颜色合并；
* 将光照和材质结构单独封装；
* 调整继承关系：
    * MiniObject -> MiniCurve -> MiniPolyline
    * MiniObject -> MiniCurve -> MiniBSpline

* 规划文件数据格式



3.15

* 暂时确定项目类图

    <img src="README.assets/MinSpline.jpg" alt="MinSpline" style="zoom:80%;" />

* 完善这些类的成员信息

* 修改成员继承关系

    <img src="README.assets/MiniSpline.jpg" alt="MiniSpline" style="zoom: 67%;" />

* 完成 MiniPolyline MiniLine 类；

* 整合 B 样条曲线类型及其算法，并用 Catch2 进行测试：

    * 寻找 u 参数区间
    * 计算 u 参数下的非零基函数
    * 计算 u 参数对应的曲线值



3.16

* 整合 B 样条曲线类型及其算法，并用 Catch2 进行测试：
    * 计算 u 参数对应的导数值
    * 计算导数曲线控制点
* 整合 B 样条曲面及其算法：
    * 计算 u,v 参数对应的点




3.17

* 修改继承关系：简化图元的继承；
* 重新整理图形绘制，将相关方法封装在 MiniHandler 中，通过 handler 获得图元的绘制信息；
* 整合 B 样条曲面及其算法，并用 Catch2 进行测试：
    * 计算 u,v 参数对应的偏导曲线的控制点
    * 计算 u,v 参数对应的偏导




3.18

* 整合 B 样条曲线类型及其算法，并用 Catch2 进行测试：
    * 计算 u,v 参数对应的非零基函数及其导数
* 整合 B 样条曲面及其算法，并用 Catch2 进行测试：
    * 验证计算 u,v 参数对应的偏导及其控制点
* 实现 Handler 类的继承关系，分别保管不同类型的图元；
* 实现对点云的拾取操作；



3.19

* 实现对样条曲线、样条曲面的基本拾取操作（基于点的拾取，因此对于物体中的“空隙”不能拾取）；
* 实现样条曲面的光照绘制；
* 实现样条曲线控制多边形的绘制；
* 实现样条曲面控制点的绘制；
* 增加对点云的顶点修改交互方法；
* 取消了对 B 样条曲线曲面的直接拾取，改为可以拾取控制顶点；
* 增加 B 样条曲线和曲面的交互方法（修改控制顶点），按照三个给定方向计算鼠标移动的投影，沿着投影最大的方向移动；
* 增加线框模式和填充模式的切换；



3.21

* 将 Polyline 类型重新分离出来，并创建对应的处理类；
* 增加 B 样条曲面显示等参线；



3.22

* 优化修改控制顶点的方法，对 x,y,z 坐标轴和鼠标移动都进行投影到投影平面上，然后计算鼠标移动在投影后的 x,y,z 轴上的移动；

* 增加显示/隐藏控制点的功能；

* 测试样条曲线最小二乘拟合点云的方法（给定参数曲线表达式，手动输入数据点进行拟合）；

* 将曲线拟合方法合并到 MiniCloud 中；

* 测试样条曲面最小二乘拟合点云的方法；

* 将曲面拟合方法合并到 MiniCloud 中；

* 现在的类继承关系

    <img src="README.assets/MiniSpline.png" alt="MiniSpline" style="zoom:80%;" />

* **增加对 MiniObject 类的全部 assert 断言（MiniBSplineSurface 没有设置断言）**；



3.24

* 实现 PCA 方法计算离散点集的 3D 包围盒算法；
* 实现 B 样条曲线的节点插入算法；
* 实现 B 样条曲面的节点插入算法（基于 B 样条曲线的节点插入，效率较低）；
* **预期使用 Catch2 测试这些算法**；



3.25

* 实现最近点的牛顿迭代计算点到 B 样条曲线投影的算法；
* 实现最近点的牛顿迭代计算点到 B 样条曲面投影的算法；
* 实现点到 B 样条曲线的一阶正交投影算法；
* 实现点到 B 样条曲面的二阶正交投影算法（从下向上投影还有问题，可能和曲率圆的方向有关系，之后测试）；



3.26

* 实现一个曲面拟合曲线的组合算法；
* 整理曲线、曲面拟合以及最近点投影、正交投影算法的测试；
* 发现一个 bug，点云拖动会消失（修复，由于测试代码中 cloud 指针被删除）；
* 修改曲面的最近点和正交投影算法的输入：增加了选择初始迭代点的参数；
* 测试空间曲线的曲面拟合方法；



4.13

* 修改图形界面：
    * 增加显示图元信息的表格控件
    * 增加输入命令名的输入框
    * 增加执行命令的按钮
    * 增加填写参数、显示算法结果（误差）的表格控件
* 增加对 B 样条曲线曲面的节点插入算法的测试程序
* 增加 B 样条曲线曲线节点优化算法以及测试程序



4.14

* 实现最小二乘渐进迭代逼近算法 LSPIA（曲线、曲面）
* 实现增量迭代计算（不断插入新的节点）
* 产生的问题
    * 不使用增量迭代，计算效果不佳，因为曲线围成的空洞里的控制点没有办法移动
    * 新的节点插入在曲线附近，导致插入节点在这里非常密集，产生曲面的质量差
    * 好处是只要限制可移动的控制点，就可以避免奇异情况



4.20

* 并入多项式类型 Polynomial 和数值积分方法 Integrator

* 提取 FindSpan 等基函数的函数到 MiniObject 中定义
  
* MiniBSpline 重命名为 MiniBSplineCurve
  
* 修改了 MiniBSplineCurve，MiniBSplineSurface 中关于 FindSpan 等基函数的部分
  
* 对应调整了 Test 程序中的测试代码



4.21

* 增加能量项约束
    $$
    E_v=E_{LS}+E_{bending}+E_{streching}
    $$

* 使用均匀样条拟合，节点插入时也均匀插入，迭代两三次

* 增加能量项，系数和为 1

* 增加参数值变化量作为终止准则

* 针对参数化算法优化 Cholesky 分解过程

* 测试 CPU 时间



4.22

* 修改并测试稀疏矩阵类型，加速最小二乘系数的计算速度；
* 修改基本数据类型
    * MiniBSpline -> MiniBSplineCurve
    * MiniPolyline -> MiniDiscreteCurve
* 删除 MiniLine 类型
* 增加 MiniParaCurve, MiniParaSurface 类型



4.23

* 增加测试稀疏矩阵的 Cholesky 分解算法，并验证其复杂度过高，最后删除；
* 实现 MiniParaCurve, MiniParaSurface 类型的基本成员函数；
* 实现加权中传播的 LSPIA 方法：使用高斯加权，即用一个矩阵作为核心平均加权，虽然能算，但是效果不佳；
* 考虑将控制点向外移动，避免折叠的方法



4.24

* 增加 RungeKutta 法和 ODE 求解器并进行测试
    * 经典 RK 方法
    * ESDI RK 方法
    * Gauss-Legendre RK 方法
* 增加了产生首一平移 Gauss-Legendre 多项式的函数
* 为只投影曲线的参数化方法增加曲线正交投影算法；



4.25

* 修改图元显示颜色的交互模式：
    * 取消了曲面的选取操作
    * 增加了 basisColor 变量，控制不同图元的显示颜色

* 增加了参量曲线、曲面的求导和求偏导操作
    * 曲面求偏导设置精度为 1e-6，因为 1e-8 精度下函数差为零（不稳定现象）

* 增加了 MiniParaCurve, MiniParaSurface 的显示控制句柄
* 测试曲线投影下的参数化方法：效果不佳；
* 增加步进修正的正交投影算法；



4.26

* 完成步进修正的正交投影方法；
* 修正有向包围盒投影方向的问题；
* 增加基于 QR 分解求解拟合线性系统的方法；
* 修改生成均匀样条的高效构造函数；
* 增加使用环面正交投影的点投影算法；
* 修正了法曲率圆的正交点投影算法；



4.27

* 增加迭代生成初始 DBS 的算法；
    * 网格单点距离最小化投影算法
    * 网格等参线点弦长最小化 + 距离最小化加权投影算法
    * 两个算法效果相近
    * 后者的系数矩阵是三对角阵
        * **考虑增加三对角阵的上下二对角分解**
        * **考虑增加追赶法求解方程组**
* 统一能量项拟合与 LSPIA 函数的参数输入；
* 发现一个问题：如果控制网格在 uv 方向相同，产生的包围盒的坐标会沿着对角线，导致一般 LSPIA 与最小二乘拟合都不能计算；
    * 修改了冒泡排序的 bug：同时调整了 Vector, MiniCloud 中的 sort 方法；（之后又进行了调整）
    * **暂时搁置包围盒的问题**



5.2

* 移除 QR 分解求解拟合线性系统的方法；

    ```cpp
    MiniBSplineSurface *MiniCloud::EnergyFittingSurfaceQR(MiniArray<Real> &u, MiniArray<Real> &v, int Nu, int Nv, int p, int q) const
    {
        // N 个控制点
        int nu = Nu - 1;
        int nv = Nv - 1;
        int m = m_points.size() - 1;
    
        // 构造样条
        MiniBSplineSurface *spline = new MiniBSplineSurface(nu + 1, nv + 1, p, q);
        MiniArray2D<Vector3f> &cpts = spline->GetPoints();
    
        // 标准化到 [0,1] 上，由于这时 u,v 参数不是分别递增的，所以需要寻找最大最小值
        Real umin = u[0];
        Real umax = u[0];
        for (int i = 0; i < m + 1; i++)
        {
            if (umin > u[i])
                umin = u[i];
            if (umax < u[i])
                umax = u[i];
        }
    
        Real du = umax - umin;
        for (int i = 0; i < m + 1; i++)
            u[i] = (u[i] - umin) / du;
    
        Real vmin = v[0];
        Real vmax = v[0];
        for (int i = 0; i < m + 1; i++)
        {
            if (vmin > v[i])
                vmin = v[i];
            if (vmax < v[i])
                vmax = v[i];
        }
    
        Real dv = vmax - vmin;
        for (int i = 0; i < m + 1; i++)
            v[i] = (v[i] - vmin) / dv;
    
        // 计算系数矩阵的 CPU 时间
        clock_t start = clock();
    
        // 填充最小二乘系数矩阵
        Matrix A(m + 1 + (nu + 1) * (nv + 1), (nu + 1) * (nv + 1));
        for (int k = 0; k < m + 1; k++)
        {
            int uSpan = FindSpan(u[k], spline->GetUDegree(), spline->GetUKnots());
            int vSpan = FindSpan(v[k], spline->GetVDegree(), spline->GetVKnots());
    
            MiniArray<Real> Bu = BasisFuncs(uSpan, u[k], spline->GetUDegree(), spline->GetUKnots());
            MiniArray<Real> Bv = BasisFuncs(vSpan, v[k], spline->GetVDegree(), spline->GetVKnots());
    
            // 每次填充 (p + 1) * (q + 1) 规模的元素
            int uIndex = (uSpan - p) * (nv + 1);
            for (int i = 0; i < p + 1; i++)
            {
                int vIndex = vSpan - q;
                for (int j = 0; j < q + 1; j++)
                   A(k, uIndex + vIndex + j) = Bu[i] * Bv[j];
                uIndex += nv + 1;
            }
        }
    
        // 构造含有能量项约束的系数矩阵，先计算出 Gauss 积分权和节点
        GaussIntegrator I(3);
        Vector nodes = I.nodes();
        Vector weights = I.weights();
    
        // 分别沿着 u,v 方向计算积分，然后保存
        Matrix Bu0(nu + 1, nu + 1);
        Matrix Bu1(nu + 1, nu + 1);
        Matrix Bu2(nu + 1, nu + 1);
        Matrix Bv0(nv + 1, nv + 1);
        Matrix Bv1(nv + 1, nv + 1);
        Matrix Bv2(nv + 1, nv + 1);
    
        MiniArray<Real> &uknots = spline->GetUKnots();
        for (int i = p; i < nu + 1; i++)
        {
            Real mid = (uknots[i + 1] + uknots[i]) / 2;
            Real ratio = (uknots[i + 1] - uknots[i]) / 2;
    
            // 三个积分节点分别计算值相加
            for (int j = 0; j < 3; j++)
            {
                MiniArray2D<Real> Bu = DerivBasisFuncs(2, i, mid + nodes[j] * ratio, p, uknots);
                for (int ui = 0; ui < p + 1; ui++)
                {
                    for (int uj = 0; uj < p + 1; uj++)
                    {
                        Bu0(i - p + ui, i - p + uj) += Bu[0][ui] * Bu[0][uj] * weights[j];
                        Bu1(i - p + ui, i - p + uj) += Bu[1][ui] * Bu[1][uj] * weights[j];
                        Bu2(i - p + ui, i - p + uj) += Bu[2][ui] * Bu[2][uj] * weights[j];
                    }
                }
            }
        }
    
        MiniArray<Real> &vknots = spline->GetVKnots();
        for (int i = q; i < nv + 1; i++)
        {
            Real mid = (vknots[i + 1] + vknots[i]) / 2;
            Real ratio = (vknots[i + 1] - vknots[i]) / 2;
    
            // 三个积分节点分别计算值相加
            for (int j = 0; j < 3; j++)
            {
                MiniArray2D<Real> Bv = DerivBasisFuncs(2, i, mid + nodes[j] * ratio, q, vknots);
                for (int vi = 0; vi < q + 1; vi++)
                {
                    for (int vj = 0; vj < q + 1; vj++)
                    {
                        Bv0(i - q + vi, i - q + vj) += Bv[0][vi] * Bv[0][vj] * weights[j];
                        Bv1(i - q + vi, i - q + vj) += Bv[1][vi] * Bv[1][vj] * weights[j];
                        Bv2(i - q + vi, i - q + vj) += Bv[2][vi] * Bv[2][vj] * weights[j];
                    }
                }
            }
        }
    
        double time = double(clock() - start) / CLOCKS_PER_SEC;
    //    qDebug() << "Prepare Matrix: " << time;
    
        // 计算获得最终矩阵的 CPU 时间
        start = clock();
    
        // 配置权重
        int Nuv = (nu + 1) * (nv + 1);
        Real alpha = 1e-6;
        Real beta = 1e-6;
        Real gamma = 1 - alpha - beta;
        for (int i = 0; i < Nuv; i++)
        {
            for (int j = 0; j < Nuv; j++)
            {
                // 获得 B_{ui}'(u) * B_{vi}(v) * B_{uj}'(u) * B_{vj}(v) 对应的索引
                int ui = i / (nv + 1);
                int uj = j / (nv + 1);
                int vi = i % (nv + 1);
                int vj = j % (nv + 1);
    
                // 累计在 B 矩阵上
                A(m + 1 + i, j) += Bu1(ui, uj) * Bv0(vi, vj) * alpha / gamma;    // BuBu^T
                A(m + 1 + i, j) += Bu0(ui, uj) * Bv1(vi, vj) * alpha / gamma;    // BvBv^T
                A(m + 1 + i, j) += Bu2(ui, uj) * Bv0(vi, vj) * beta / gamma;     // BuuBuu^T
                A(m + 1 + i, j) += Bu0(ui, uj) * Bv2(vi, vj) * beta / gamma;     // BvvBvv^T
                A(m + 1 + i, j) += 2 * Bu1(ui, uj) * Bv1(vi, vj) * beta / gamma; // 2 * BuvBuv^T
            }
        }
    
        time = double(clock() - start) / CLOCKS_PER_SEC;
    //    qDebug() << "Final Matrix: " << time;
    
        // 计算求解整个系统的 CPU 时间
        start = clock();
    
        // 求解最小二乘+能量约束系统
        Vector bx(m + 1 + Nuv), by(m + 1 + Nuv), bz(m + 1 + Nuv);
        for (int j = 0; j < m + 1; j++)
        {
            bx[j] = m_points[j].x();
            by[j] = m_points[j].y();
            bz[j] = m_points[j].z();
        }
    
        // 作 QR 分解
        int row = A.row();
        int col = A.col();
        Vector d(col);
        Matrix B = A.QR(d);
    
        // Hb = b - beta * v * v^T * b = b - (beta * v^T * b) * v
        for (int i = 0; i < col; i++)
        {
            // Householder 变换作用于 bx
            Real sum = bx[i];
            for (int j = i + 1; j < row; j++)
                sum += B(j, i) * bx[j];
    
            sum *= d[i];
            bx[i] -= sum;
            for (int j = i + 1; j < row; j++)
                bx[j] -= B(j, i) * sum;
    
            // Householder 变换作用于 by
            sum = by[i];
            for (int j = i + 1; j < row; j++)
                sum += B(j, i) * by[j];
    
            sum *= d[i];
            by[i] -= sum;
            for (int j = i + 1; j < row; j++)
                by[j] -= B(j, i) * sum;
    
            // Householder 变换作用于 bz
            sum = bz[i];
            for (int j = i + 1; j < row; j++)
                sum += B(j, i) * bz[j];
    
            sum *= d[i];
            bz[i] -= sum;
            for (int j = i + 1; j < row; j++)
                bz[j] -= B(j, i) * sum;
        }
    
        // 获得一个更小规模的矩阵
        Vector b1(col), b2(col), b3(col);
        for (int i = 0; i < col; i++)
        {
            b1[i] = bx[i];
            b2[i] = by[i];
            b3[i] = bz[i];
        }
    
        // 求解上半部分的线性方程组
        LinearSolver solver;
        Vector X = solver.solveUpper(B, b1);
        Vector Y = solver.solveUpper(B, b2);
        Vector Z = solver.solveUpper(B, b3);
    
        for (int i = 0; i < nu + 1; i++)
        {
            for (int j = 0; j < nv + 1; j++)
            {
                int index = i * (nv + 1) + j;
                cpts[i][j] = {X[index], Y[index], Z[index]};
            }
        }
    
        time = double(clock() - start) / CLOCKS_PER_SEC;
    //    qDebug() << "Solving Time: " << time;
    
        return spline;
    }
    ```

    

* 将有向包围盒和归一化算法提取出来单独作为外部函数；



5.3

* 增加累积弦长最小化求解拟合线性系统的方法；

* 增加平滑化函数约束求解拟合线性系统的算法；

* 移除 WLSPIA 方法

    ```cpp
    MiniBSplineSurface *MiniCloud::WLSPIASurface(MiniArray<Real> &u, MiniArray<Real> &v, int Nu, int Nv, int p, int q, Real tol) const
    {
        // N 个控制点，m + 1 个数据点和对应参数
        int nu = Nu - 1;
        int nv = Nv - 1;
        int m = m_points.size() - 1;
    
        // 构造样条
        MiniBSplineSurface *spline = new MiniBSplineSurface(nu + 1, nv + 1, p, q);
        MiniArray2D<Vector3f> &P = spline->GetPoints();
        const MiniArray<Vector3f> &Q = GetPoints();
    
        // 几何迭代法
        Real e = tol * 10;
        int times = 0;
        while (e > tol && times < ITERATION)
        {
            MiniArray2D<Real> basisSum(nu + 1, MiniArray<Real>(nv + 1));         // 记录每个控制点对应的那组基函数和
            MiniArray2D<Vector3f> deltaSum(nu + 1, MiniArray<Vector3f>(nv + 1)); // 记录每个控制点对应的那组偏差向量加权和
    
            // 对每个数据点计算增量
            for (int i = 0; i < m + 1; i++)
            {
                Vector3f delta = Q[i] - spline->SurfacePoint(u[i], v[i]);
    
                // 计算此参数对应的非零基函数值
                int ju = FindSpan(u[i], spline->GetUDegree(), spline->GetUKnots());
                int jv = FindSpan(v[i], spline->GetVDegree(), spline->GetVKnots());
                MiniArray<Real> basisFuncU = BasisFuncs(ju, u[i], spline->GetUDegree(), spline->GetUKnots());
                MiniArray<Real> basisFuncV = BasisFuncs(jv, v[i], spline->GetVDegree(), spline->GetVKnots());
    
                // 累计参数
                for (int ku = 0; ku < int(basisFuncU.size()); ku++)
                {
                    for (int kv = 0; kv < int(basisFuncV.size()); kv++)
                    {
                        basisSum[ju - p + ku][jv - q + kv] += basisFuncU[ku] * basisFuncV[kv];
                        deltaSum[ju - p + ku][jv - q + kv] = deltaSum[ju - p + ku][jv - q + kv] + delta * basisFuncU[ku] * basisFuncV[kv];
                    }
                }
            }
    
            // 得到控制点增量的平均值并加上增量，注意由于存在空洞，不是所有控制点都可以移动
            for (int i = 0; i < nu + 1; i++)
            {
                for (int j = 0; j < nv + 1; j++)
                {
                    if (fabs(basisSum[i][j]) > FLT_MIN)
                        deltaSum[i][j] = deltaSum[i][j] / basisSum[i][j];
                }
            }
    
            // 加权矩阵规模是控制点网格的一半以上
            int wu = (nu + 1) / 4 + 1;
            int wv = (nv + 1) / 4 + 1;
            int cu = wu / 2;
            int cv = wv / 2;
            Matrix weight(wu, wv);
            for (int i = 0; i < wu; i++)
            {
                int di = cu - std::abs(i - cu);
                for (int j = 0; j < wv; j++)
                {
                    int dj = cv - std::abs(j - cv);
    //                weight(i, j) = pow(2, std::min(di, dj));
                    weight(i, j) = pow(1.5, di + dj);
                }
            }
    //        cout << weight << endl;
    
            for (int i = 0; i < nu + 1; i++)
            {
                for (int j = 0; j < nv + 1; j++)
                {
                    Real s = 0;
                    Vector3f delta = {0, 0, 0};
                    for (int ki = 0; ki < wu; ki++)
                    {
                        for (int kj = 0; kj < wv; kj++)
                        {
                            int wi = i + ki - wu / 2;
                            int wj = j + kj - wv / 2;
                            if (wi >= 0 && wi < nu + 1 && wj >= 0 && wj < nv + 1)
                            {
                                delta = delta + deltaSum[wi][wj] * weight(ki, kj);
                                s += weight(ki, kj);
                            }
                        }
                    }
                    // 取平均值
                    P[i][j] = P[i][j] + delta / s;
                }
            }
    
            if (times < 10)
            {
                // 加权光滑化
                for (int i = 0; i < wu; i++)
                {
                    int di = cu - std::abs(i - cu);
                    for (int j = 0; j < wv; j++)
                    {
                        int dj = cv - std::abs(j - cv);
                        weight(i, j) = pow(2, std::min(di, dj));
    //                    weight(i, j) = pow(1.2, di + dj);
                    }
                }
    
                MiniArray2D<Vector3f> Q(nu + 1, MiniArray<Vector3f>(nv + 1));
                for (int i = 0; i < nu + 1; i++)
                {
                    for (int j = 0; j < nv + 1; j++)
                    {
                        Real s = 0;
                        Vector3f delta = {0, 0, 0};
    
                        for (int ki = 0; ki < wu; ki++)
                        {
                            for (int kj = 0; kj < wv; kj++)
                            {
                                int wi = i + ki - wu / 2;
                                int wj = j + kj - wv / 2;
                                if (wi >= 0 && wi < nu + 1 && wj >= 0 && wj < nv + 1)
                                {
                                    delta = delta + P[wi][wj] * weight(ki, kj);
                                    s += weight(ki, kj);
                                }
                            }
                        }
    
                        // 取平均值
                        Q[i][j] = delta / s;
                    }
                }
    
                for (int i = 0; i < nu + 1; i++)
                    for (int j = 0; j < nv + 1; j++)
                        P[i][j] = Q[i][j];
            }
    
            // 计算区间误差
            e = 0;
            for (int i = 0; i < m + 1; i++)
            {
                Real err = length(spline->SurfacePoint(u[i], v[i]) - m_points[i]);
                e += err * err / (m + 1);
            }
            times++;
        }
        return spline;
    }
    ```

    

* **考虑优化 SSDE 初始化 DBS 的算法：最小化横向与纵向累计弦长函数之和，修改论文中的矩阵；**移除 SSDE 生成 DBS 的方法

    ```cpp
    MiniBSplineSurface *MiniCloud::SSDEInitialDBS(int Nu, int Nv, int p, int q, Real tol, Real ru, Real rv) const
    {
        // 计算一个包围盒
        MiniArray<Vector3f> box = OrientBoundingBox(m_points);
    
        // N 个控制点，m + 1 个数据点和对应参数
        int nu = Nu - 1;
        int nv = Nv - 1;
        int m = m_points.size() - 1;
    
        // 构造样条
        MiniBSplineSurface *spline = new MiniBSplineSurface(nu + 1, nv + 1, p, q);
        MiniArray2D<Vector3f> &P = spline->GetPoints();
        const MiniArray<Vector3f> &Q = GetPoints();
    
        // 调整为右手坐标系，两个方向的先后顺序会决定光照时的法向
        if (dot(cross(box[0], box[2]), box[4]) < 0)
        {
            Vector3f tx = box[0];
            box[0] = box[2];
            box[2] = tx;
    
            // 注意正反方向的都要交换
            tx = box[1];
            box[1] = box[3];
            box[3] = tx;
        }
    
        // 稍微放大包围盒，不要放大原点位置
        box[0] = box[0] * ru;
        box[1] = box[1] * ru;
        box[2] = box[2] * rv;
        box[3] = box[3] * rv;
    
        // 根据包围盒设置初始控制点
        MiniArray<Vector3f> oxy(3);
        oxy[0] = box.back() + box[1] + box[3] + box[5];
        oxy[1] = normalize(box[0]);
        oxy[2] = normalize(box[2]);
        Real su = (length(box[0]) + length(box[1])) / nu;
        Real sv = (length(box[2]) + length(box[3])) / nv;
    
        for (int i = 0; i < nu + 1; i++)
            for (int j = 0; j < nv + 1; j++)
                P[i][j] = oxy[0] + oxy[1] * su * i + oxy[2] * sv * j;
    
        // 设置平均权重
        MiniArray<Real> w(m + 1);
        for (int i = 0; i < m + 1; i++)
            w[i] = 0.5;
    
        // 加细分割曲面网格
        nu *= 2, nv *= 2;
    
        LinearSolver solver;
        int times = 0;
        Real e = 0, e1 = tol * 10;
        Real kappa = 1e-3;
        while (fabs(e - e1) > tol && times < ITERATION)
        {
            // 记录新的误差
            e = e1;
    
            // 构造相交网格及其参数
            MiniArray<Real> u((nu + 1) * (nv + 1));
            MiniArray<Real> v((nu + 1) * (nv + 1));
            MiniArray2D<Vector3f> R(nu + 1, MiniArray<Vector3f>(nv + 1));
            MiniArray2D<Vector3f> N(nu + 1, MiniArray<Vector3f>(nv + 1));
            for (int i = 1; i < nu; i++)
            {
                for (int j = 1; j < nv; j++)
                {
                    u[i * (nv + 1) + j] = 1.0 / nu * i;
                    v[i * (nv + 1) + j] = 1.0 / nv * j;
    
                    // 内部点计算法向
                    MiniArray2D<Vector3f> S = spline->SurfaceDerivPoint(1, 1.0 / nu * i, 1.0 / nv * j);
                    R[i][j] = S[0][0];
                    N[i][j] = normalize(cross(S[1][0], S[0][1]));
                }
            }
    
            // 外插计算边界点处的法向
            // 左右边界外插
            for (int i = 1; i < nu; i++)
            {
                Real t = 1.0 / nu * i;
                u[i * (nv + 1)] = t;
                v[i * (nv + 1)] = 0;
                R[i][0] = spline->SurfacePoint(t, 0);
                N[i][0] = normalize(N[i][1] * 2 - N[i][2]);
    
                u[i * (nv + 1) + nv] = t;
                v[i * (nv + 1) + nv] = 1;
                R[i][nv] = spline->SurfacePoint(t, 1);
                N[i][nv] = normalize(N[i][nv-1] * 2 - N[i][nv-2]);
            }
    
            // 上下边界外插
            for (int i = 1; i < nv; i++)
            {
                Real t = 1.0 / nv * i;
                u[i] = 0;
                v[i] = t;
                R[0][i] = spline->SurfacePoint(0, t);
                N[0][i] = normalize(N[1][i] * 2 - N[2][i]);
    
                u[nu * (nv + 1) + i] = 1;
                v[nu * (nv + 1) + i] = t;
                R[nu][i] = spline->SurfacePoint(1, t);
                N[nu][i] = normalize(N[nu-1][i] * 2 - N[nu-2][i]);
            }
    
            // 四个角点外插
            u[0] = 0, v[0] = 0;
            u[nv] = 0, v[nv] = 1;
            u[nu * (nv + 1)] = 1, v[nu * (nv + 1)] = 0;
            u[nu * (nv + 1) + nv] = 1, v[nu * (nv + 1) + nv] = 1;
    
            R[0][0] = spline->SurfacePoint(0, 0);
            R[0][nv] = spline->SurfacePoint(0, 1);
            R[nu][0] = spline->SurfacePoint(1, 0);
            R[nu][nv] = spline->SurfacePoint(1, 1);
    
            N[0][0] = normalize(N[0][1] * 2 - N[0][2]);
            N[0][nv] = normalize(N[0][nv-1] * 2 - N[0][nv-2]);
            N[nu][0] = normalize(N[nu-1][0] * 2 - N[nu-2][0]);
            N[nu][nv] = normalize(N[nu-1][nv] * 2 - N[nu-2][nv]);
    
            // 遍历每行控制点构造线性系统求解
            e1 = 0;
            for (int i = 0; i < nu + 1; i++)
            {
                // 先计算出两个目标向量
                Vector b(nv + 1);
                for (int j = 0; j < nv + 1; j++)
                {
                    Vector c(5);
    
                    // 对每个点云中的点计算误差
                    for (int mu = 0; mu < m + 1; mu++)
                    {
                        Real a = 1 / pow(length(R[i][j] - Q[mu]), 4);
                        c[0] += a;
                        c[1] += a * Q[mu].x();
                        c[2] += a * Q[mu].y();
                        c[3] += a * Q[mu].z();
                        c[4] += a * (Q[mu].x() * Q[mu].x() + Q[mu].y() * Q[mu].y() + Q[mu].z() * Q[mu].z());
                    }
    
                    // 计算出极小值参数
                    Real lambda = (c[1] * N[i][j].x() + c[2] * N[i][j].y() + c[3] * N[i][j].z()) / c[0];
                    b[j] = lambda - dot(R[i][j], N[i][j]);
                }
    
                // 开始加权
                b[0] *= 1 - kappa;
                b[nv] *= 1 - kappa;
                for (int j = 1; j < nv; j++)
                    b[j] = b[j] * (1 - kappa) + dot(R[i][j - 1] - R[i][j] * 2 + R[i][j + 1], N[i][j]) * kappa;
    
                // 然后计算系数矩阵
                Matrix A(nv + 1, nv + 1);
                A(0, 0) = kappa;
                A(nv, nv) = kappa;
                for (int j = 1; j < nv; j++)
                {
                    A(j, j - 1) = -dot(N[j][j - 1], N[j][j]) * kappa;
                    A(j, j) = 2 * kappa;
                    A(j, j + 1) = -dot(N[j][j], N[j][j + 1]) * kappa;
                }
                A(1, 0) = A(nv - 1, nv) = 0;
    
                for (int j = 0; j < nv + 1; j++)
                    A(j, j) += 1 - kappa;
    
                // 求解方程得到参数 t，然后修改这一行
                Vector t = solver.solve(A, b, true);
                for (int j = 0; j < nv + 1; j++)
                {
                    R[i][j] = R[i][j] + N[i][j] * t[j];
                    e1 += length(spline->SurfacePoint(1.0 / nu * i, 1.0 / nv * j) - R[i][j]);
                }
            }
            e1 /= (nu + 1) * (nv + 1);
    
            MiniCloud cloud(MiniArray<Vector3f>(0));
            MiniArray<Vector3f> &final = cloud.GetPoints();
            for (int i = 0; i < nu + 1; i++)
                for (int j = 0; j < nv + 1; j++)
                    final.push_back(R[i][j]);
    
            // 计算完成后即可移除原先曲面
            delete spline;
    
            // 重新拟合曲面
            spline = cloud.ChordFittingSurface(u, v, Nu, Nv, p, q, w, 0, 0);
            times++;
    
            // 逐渐减小光滑化操作
            kappa *= 1;
        }
        return spline;
    }
    ```

    

* 增加使用稀疏矩阵求解的 Jacobi, GS, SOR 迭代方法；

* 引入 Eigen3 库到 MiniObject.h 头文件中；

* 用 Eigen3 中的算法替换掉 FittingSurface 的两个方法中稀疏矩阵和 Cholesky 分解求解方程的算法；



5.4

* 修改了 B 样条曲面和参量曲面中的计算光照法向时 Su, Sv 的叉乘顺序；
* 向有向包围盒算法中加入了调整坐标系的方法，使得包围盒的最短正向轴总是沿着 y 轴正向，其余坐标轴与其构成右手坐标系；
* 封装了毕业论文中的边界曲线参数化算法；





* 引入 OCC 的数据结构：熟悉裁剪曲面、基本的样条数据结构等；
* 通过 OCC 架构实现场景中图元属性的保存和加载：为 MiniObject 构成的场景设置文件数据，能够同时将多条曲线、网格保存在一个文件中，提供读取文件的方法；



* **由于加入了 Eigen 库，Test 程序无法编译，考虑使用 makefile 协调**；
* **考虑将 Basic 部分的测试和 MiniObject 部分的测试分离**；
* **实现对所有算法的测试程序**；



* 整理当前数据结构
    * 统一 MiniDiscreteCurve, MiniParaCurve, MiniBSplineCurve 的继承关系：都继承自 MiniCurve 
    * 统一 MiniParaSurface, MiniBSplineSurface 的继承关系：都继承自 MiniSurface
* 移除网格中的材质信息；
* 添加菜单功能：打开和保存文件；
* 增加创建 B 样条曲线曲面的菜单功能（设置中心位置、控制点数和节点数（自动均匀节点）以及是否是周期样条，曲线设置长度，曲面设置长和宽，控制点自动均匀排列，可双击曲线、曲面设置节点向量）；
* 封装各种不同的算法
    * 完善菜单功能，并合并之前的各种算法；
    * 封装空间封闭曲线参数化函数；
* 增加图元信息显示
* 增加坐标轴显示（虚线坐标，可以显示或隐藏）
* 增加观察者坐标显示
* 增加纹理映射（斑马纹）
* 优化图像数据显示形式（修改当前的数据点光照格式）
* 增加框选功能



* 增加嵌入 RungeKutta 方法（允许自适应步长）*
* 增加对稀疏矩阵迭代求解算法的测试；