当前位置：首页 > news >正文

移动最小二乘法（Moving Least Squares, MLS）原理和c++实现

news 来源：原创 2025/4/27 17:12:02

移动最小二乘法（Moving Least Squares, MLS）原理和c++实现

一、算法原理

移动最小二乘法（MLS） 是一种局部加权回归方法，通过对每个查询点邻域内的数据点进行加权最小二乘拟合，生成光滑的曲面或曲线。与传统最小二乘法不同，MLS的拟合函数随查询点位置动态变化，适用于非均匀采样数据与复杂几何建模。

二、数学推导

1. 基本模型

对于查询点 $(x)$ ，在其邻域内找到最优拟合函数 $(f (x))$ ，最小化加权残差平方和：
$\min_{f \in \Pi_m} \sum_{i=1}^n w(\|x - x_i\|) \left( f(x_i) - y_i \right)^2$
其中：

$\Pi_m )$ 为 $(m)$ 次多项式空间（如线性、二次多项式）。
$(w (d))$ 为权重函数（如高斯核 $w(d) = e^{-d^2/h^2} )$ ）。

2. 多项式展开

假设 $\mathbf{p}(x)^T \mathbf{a} )$ ，其中 $\mathbf{p}(x) )$ 为基函数向量，例如：

线性基： $\mathbf{p}(x) = [1, x, y]^T )$ （2D）
二次基： $\mathbf{p}(x) = [1, x, y, x^2, xy, y^2]^T )$ （2D）

3. 求解系数

构建加权设计矩阵 $\mathbf{P} )$ 和权重矩阵 $\mathbf{W} )$ ：
$\mathbf{P} = \begin{bmatrix} \mathbf{p}(x_1)^T \\ \mathbf{p}(x_2)^T \\ \vdots \\ \mathbf{p}(x_n)^T \end{bmatrix}, \quad \mathbf{W} = \text{diag}(w(\|x - x_1\|), \dots, w(\|x - x_n\|))$
目标函数转化为：
$\min_{\mathbf{a}} (\mathbf{y} - \mathbf{P}\mathbf{a})^T \mathbf{W} (\mathbf{y} - \mathbf{P}\mathbf{a})$
解得系数：
$\mathbf{a} = (\mathbf{P}^T \mathbf{W} \mathbf{P})^{-1} \mathbf{P}^T \mathbf{W} \mathbf{y}$

三、应用场景

点云平滑与重建
- 去除噪声，生成连续曲面（如3D扫描数据处理）。
图像变形与插值
- 实现非刚性形变（如图像配准、医学图像处理）。
曲线/曲面拟合
- 从离散数据点生成光滑曲线（如CAD建模）。

四、C++代码实现

1. 依赖库

Eigen：线性代数运算。
FLANN：快速邻域搜索（需安装 libflann-dev）。

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <flann/flann.hpp>

// 2D点结构体
struct Point2D {
    double x, y;
    Point2D(double x_, double y_) : x(x_), y(y_) {}
};

// MLS拟合函数
double mlsFit(const std::vector<Point2D>& points, 
              const Point2D& query, 
              double radius, 
              double h) {
    // 1. 邻域搜索
    flann::Matrix<double> dataset(new double[points.size()*2], points.size(), 2);
    for (size_t i = 0; i < points.size(); ++i) {
        dataset[i][0] = points[i].x;
        dataset[i][1] = points[i].y;
    }
    flann::Index<flann::L2<double>> index(dataset, flann::KDTreeIndexParams(4));
    index.buildIndex();

    // 查询邻域点
    double query_point[2] = {query.x, query.y};
    std::vector<int> indices;
    std::vector<double> dists;
    index.radiusSearch(flann::Matrix<double>(query_point, 1, 2), 
                     indices, dists, radius, flann::SearchParams(128));

    // 2. 构建矩阵P, W, y
    int m = 3; // 线性基: 1, x, y
    Eigen::MatrixXd P(indices.size(), m);
    Eigen::VectorXd y(indices.size());
    Eigen::MatrixXd W = Eigen::MatrixXd::Zero(indices.size(), indices.size());

    for (size_t i = 0; i < indices.size(); ++i) {
        int idx = indices[i];
        double dx = query.x - points[idx].x;
        double dy = query.y - points[idx].y;
        double weight = exp(-(dx*dx + dy*dy)/(h*h));
        
        P(i, 0) = 1.0;
        P(i, 1) = points[idx].x;
        P(i, 2) = points[idx].y;
        y(i) = points[idx].y; // 假设拟合y值
        W(i, i) = weight;
    }

    // 3. 求解系数a
    Eigen::MatrixXd PTW = P.transpose() * W;
    Eigen::MatrixXd PTWP = PTW * P;
    Eigen::VectorXd a = PTWP.inverse() * PTW * y;

    // 4. 预测query点的y值
    Eigen::VectorXd p_query(m);
    p_query << 1.0, query.x, query.y;
    return p_query.dot(a);
}

int main() {
    // 生成带噪声的测试数据 (y = 2x + 1 + noise)
    std::vector<Point2D> points;
    for (int i = 0; i < 100; ++i) {
        double x = i * 0.1;
        double y = 2*x + 1 + 0.1 * ((double)rand()/RAND_MAX - 0.5);
        points.emplace_back(x, y);
    }

    // MLS参数
    double radius = 1.0; // 邻域半径
    double h = 0.5;      // 高斯核带宽

    // 在x=5.0处预测y值
    Point2D query(5.0, 0.0);
    double y_pred = mlsFit(points, query, radius, h);
    std::cout << "Predicted y at x=5.0: " << y_pred << std::endl;

    return 0;
}

2. 代码说明

邻域搜索：使用FLANN库快速查找查询点邻域内的数据点。
矩阵构建：构建设计矩阵 $\mathbf{P} )$ 、权重矩阵 $\mathbf{W} )$ 和观测向量 $\mathbf{y} )$ 。
系数求解：通过加权最小二乘公式计算多项式系数 $\mathbf{a} )$ 。
预测值计算：利用拟合多项式预测查询点的函数值。

五、参数选择与优化

邻域半径 $(r)$ ：过小导致欠拟合，过大增加计算量。建议根据数据密度调整。
带宽 $(h)$ ：控制权重衰减速度，影响平滑程度。可通过交叉验证优化。
多项式阶数：高阶多项式拟合能力强但易过拟合，通常选择线性或二次。

六、扩展应用示例

点云平滑2D点云：对每个点应用MLS，用邻域点的加权平均更新其位置。

void smoothPointCloud(std::vector<Point2D>& points, double radius, double h) {
    std::vector<Point2D> smoothed = points;
    for (auto& p : smoothed) {
        p.y = mlsFit(points, p, radius, h);
    }
    points = smoothed;
}

点云平滑3D点云：

using namespace pcl;

// 高斯权重函数
double gaussianWeight(double x, double y, double x0, double y0, double h) {
    return exp(-((x - x0) * (x - x0) + (y - y0) * (y - y0)) / (h * h));
}

// 移动最小二乘计算
Eigen::VectorXd movingLeastSquares(const vector<PointXYZ>& points, const vector<double>& values, double x0, double y0, double h, int polyOrder = 2) {
    int N = points.size();
    int numCoeffs = (polyOrder + 1) * (polyOrder + 2) / 2; // 二维多项式项数

    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(numCoeffs, numCoeffs);
    Eigen::VectorXd b = Eigen::VectorXd::Zero(numCoeffs);

    for (int i = 0; i < N; ++i) {
        double w = gaussianWeight(points[i].x, points[i].y, x0, y0, h);

        // 构建多项式基矩阵
        Eigen::VectorXd p(numCoeffs);
        p(0) = 1;
        if (polyOrder >= 1) {
            p(1) = points[i].x;
            p(2) = points[i].y;
        }
        if (polyOrder >= 2) {
            p(3) = points[i].x * points[i].x;
            p(4) = points[i].x * points[i].y;
            p(5) = points[i].y * points[i].y;
        }

        A += w * (p * p.transpose());
        b += w * (p * values[i]);
    }

    // 解线性方程 A * a = b
    Eigen::VectorXd coeffs = A.ldlt().solve(b);
    return coeffs;
}

void TestMLS3D()
{
  // 生成示例点云
    vector<PointXYZ> points = { {0, 0, 1}, {1, 0, 2}, {0, 1, 2}, {1, 1, 3}, {0.5, 0.5, 2.5} };
    vector<double> values = {1, 2, 2, 3, 2.5};

    double x0 = 0.5, y0 = 0.5, h = 1.0;
    Eigen::VectorXd coeffs = movingLeastSquares(points, values, x0, y0, h);

    // 计算平滑后点值
    double smoothedValue = coeffs(0) + coeffs(1) * x0 + coeffs(2) * y0 + coeffs(3) * x0 * x0 + coeffs(4) * x0 * y0 + coeffs(5) * y0 * y0;
    
    cout << "MLS 估计的值：" << smoothedValue << endl;
}