---

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

一、算法原理

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

---

二、数学推导

1. 基本模型

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

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

2. 多项式展开

假设 $(f(x)=\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}=\left[\begin{array}{c}\mathbf{p}(x_1{)}^T\\ \mathbf{p}(x_2{)}^T\\ ⋮\\ \mathbf{p}(x_n{)}^T\end{array}\right],\mathbf{W}=\text{diag}(w(\Vert x-x_1\Vert ),\dots ,w(\Vert x-x_n\Vert ))$$
目标函数转化为：
$$\underset{\mathbf{a}}{\min }(\mathbf{y}-\mathbf{Pa}{)}^T\mathbf{W}(\mathbf{y}-\mathbf{Pa})$$
解得系数：
$$\mathbf{a}=({\mathbf{P}}^T\mathbf{WP}{)}^{-1}{\mathbf{P}}^T\mathbf{Wy}$$

---

三、应用场景

1. 点云平滑与重建
   • 去除噪声，生成连续曲面（如3D扫描数据处理）。
2. 图像变形与插值
   • 实现非刚性形变（如图像配准、医学图像处理）。
3. 曲线/曲面拟合
   • 从离散数据点生成光滑曲线（如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, yint m = 3; // 线性基: 1, x, yEigen::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. 求解系数aEigen::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. 代码说明

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

---

五、参数选择与优化

• 邻域半径 $(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 = bEigen::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;
}

---

通过合理选择参数与高效实现，移动最小二乘法可广泛用于：

• 点云平滑（如 PCL 中的 MLS 滤波器）
• 曲面重建（如基于 MLS 的表面重建）
• 数值计算（如无网格法中的插值）