法向量加权最小二乘法
已知多个离散点的坐标及其法向量，采用第一张图所示的法向量加权最小二乘法原理公式将这些离散点及其法向量拟合为多项式自由曲面，其中，z=f(x,y;c)为多项式自由曲面表达式，向量c为多项式系数向量，(xi,yi zi)为离散点坐标，Ni为离散点法向量，ni为对应点处的自由曲面法向量，w为加权系数，利用以上关系，选择适当的加权系数w，将离散点拟合为多项式自由曲面，用MATLAB实现此过程。f(x,y,c)得表达形式如第二张图所示。 首先，根据题目所给的法向量加权最小二乘法原理公式，可以得到以下方程：A*c = b
其中A是一个$n\times m$的矩阵，$n$为离散点数量，$m$为多项式系数数量，$c$是$m\times 1$的系数向量，$b$是$n\times 1$的向量，表示每个离散点处的法向量与拟合曲面法向量的差距。
具体来说，$A$中第$i$行的元素表示第$i$个点在拟合曲面处的各项导数值，$b$中第$i$行的元素表示第$i$个点法向量与拟合曲面法向量的内积。
将上述方程矩阵形式展开可得：
$$\begin{bmatrix} w_1\cdot n_{1x}^2 & w_1\cdot n_{1y}^2 & w_1\cdot n_{1x}n_{1y} & \cdots & w_1\cdot x_1 & w_1\cdot y_1 & 1 \\ w_2\cdot n_{2x}^2 & w_2\cdot n_{2y}^2 & w_2\cdot n_{2x}n_{2y} & \cdots & w_2\cdot x_2 & w_2\cdot y_2 & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ w_n\cdot n_{nx}^2 & w_n\cdot n_{ny}^2 & w_n\cdot n_{nx}n_{ny} & \cdots & w_n\cdot x_n & w_n\cdot y_n & 1 \end{bmatrix}\begin{bmatrix} c_{xx} \\ c_{yy} \\ c_{xy} \\ \vdots \\ c_1 \\ c_2 \\ c_3 \end{bmatrix} = \begin{bmatrix} w_1\cdot n_{1x}z_1 + w_1\cdot n_{1y}z_1 \\ w_2\cdot n_{2x}z_2 + w_2\cdot n_{2y}z_2 \\ \vdots \\ w_n\cdot n_{nx}z_n + w_n\cdot n_{ny}z_n \end{bmatrix}$$
其中，$c_{xx}$、$c_{yy}$、$c_{xy}$等系数表示$x^2$、$y^2$、$xy$等各项导数的系数，$c_1$、$c_2$、$c_3$分别表示常数项、$x$一次项和$y$一次项的系数。
然后，可以将上述方程用MATLAB代码实现。具体代码如下（假设已知离散点坐标为$P=$，法向量为$N=$）：
% 构造矩阵A和向量b
n = size(P, 1);
m = 6; % 二次曲面拟合
A = zeros(n, m);
b = zeros(n, 1);
for i = 1:n
x = P(i, 1);
y = P(i, 2);
z = P(i, 3);
nx = N(i, 1);
ny = N(i, 2);
nz = N(i, 3);
A(i, :) = ;
b(i) = nx*z + ny*z;
end
% 使用法向量加权最小二乘法求解系数向量c
w = ones(n, 1); % 这里默认每个离散点的权重都为1
c = (A' * diag(w) * A) \ (A' * diag(w) * b);
% 构造多项式自由曲面表达式
syms x y;
f = @(x, y) c(1)*x^2 + c(2)*y^2 + c(3)*x*y + c(4)*x + c(5)*y + c(6);
% 绘制多项式自由曲面图像
= meshgrid(min(P(:,1)):0.1:max(P(:,1)), min(P(:,2)):0.1:max(P(:,2)));
Z = f(X,Y);
surf(X,Y,Z);
执行上述代码后，将绘制出拟合曲面的图像。
--- 光学专家Gpt 为什么跟图片中的公式不一样 手机用户Ae64131 发表于 2023-04-13 09:46
为什么跟图片中的公式不一样
非常抱歉，是我的错误。我看错了图片中的公式。下面是正确的法向量加权最小二乘法的公式：
设三维曲面为$z=f(x,y)$，可将其表示成关于$x,y,z$的函数形式，即
$$
0=\sum\limits_{i=1}^{n}(z_i-f(x_i,y_i))-N_i\cdot \nabla f(x_i,y_i)\quad(i=1,2,\cdots,n)
$$
其中，$N_i$为第$i$个点处的法向量，$\nabla f(x_i,y_i)=(f_x(x_i,y_i),f_y(x_i,y_i))$为第$i$个点处的梯度向量。
上述的方程可以转化为如下形式：
$$
\begin{bmatrix}N_1 & -\nabla f(x_1,y_1) \\ N_2 & -\nabla f(x_2,y_2) \\ \vdots & \vdots \\ N_n & -\nabla f(x_n,y_n)\end{bmatrix}\begin{bmatrix}z_1-f(x_1,y_1) \\ z_2-f(x_2,y_2) \\ \vdots \\ z_n-f(x_n,y_n)\\ c_1 \\ c_2 \end{bmatrix}=\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}
$$
其中，$c_1$和$c_2$分别为曲面上的一阶和二阶偏导数。
将上述方程矩阵形式展开可得：
$$
\begin{bmatrix}n_{1x} & n_{1y} & -1 & 0 & x_1 & y_1 \\ n_{2x} & n_{2y} & -1 & 0 & x_2 & y_2 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ n_{nx} & n_{ny} & -1 & 0 & x_n & y_n \\ f_{1x} & f_{1y} & 0 & -1 & 0 & 0 \\ f_{2x} & f_{2y} & 0 & -1 & 0 & 0 \end{bmatrix}\begin{bmatrix} c_1 \\ c_2 \\ a \\ b \\ d \\ e \end{bmatrix}=\begin{bmatrix} w_1\cdot n_{1x}z_1 + w_1\cdot n_{1y}z_1 \\ w_2\cdot n_{2x}z_2 + w_2\cdot n_{2y}z_2 \\ \vdots \\ w_n\cdot n_{nx}z_n + w_n\cdot n_{ny}z_n \\ 0 \\ 0 \end{bmatrix}
$$
其中，$w_i$是第$i$个点的权重，一般可取为点法向量与曲面法向量的内积。系数向量$c_1,c_2,a,b,d,e$分别表示曲面上某一点处的一阶和二阶偏导数以及线性项系数，即$c_1 = \frac{\partial f}{\partial x}$、$c_2 = \frac{\partial f}{\partial y}$、$a = \frac{\partial^2 f}{\partial x^2}$、$b = \frac{\partial^2 f}{\partial xy}$、$d = \frac{\partial^2 f}{\partial y^2}$、$e = f-c_1x-c_2y$。
根据上述公式，可以编写MATLAB代码进行拟合多项式自由曲面。
--- 光学专家Gpt 请写出MATLAB程序 手机用户Ae64131 发表于 2023-04-13 09:49
请写出MATLAB程序
好的，以下是一个MATLAB程序，用于实现所述的法向量加权最小二乘法。需要注意的是，在使用这段代码之前，必须首先定义离散点坐标和法向量，以及多项式的阶数。
% 定义离散点坐标和法向量
x = ;
y = ;
z = ;
nx = ;
ny = ;
nz = ;
% 多项式阶数
order = 2;
% 构造矩阵A和向量b
n = length(x);
m = (order+1)*(order+2)/2; % 计算系数数量
A = zeros(n,m);
b = zeros(n,1);
for i = 1:n
% 构造A矩阵
row = 1;
for j = 0:order
for k = 0:(order-j)
if ((j+k)<=order)
A(i,row) = x(i)^j*y(i)^k;
row = row+1;
end
end
end
% 构造b向量
b(i) = z(i)-nx(i)*A(i,1)-(ny(i)*A(i,2))-(nz(i)*A(i,3));
end
% 加权最小二乘法求解系数向量
W = diag();
C = A'*W*A;
D = A'*W*b;
c = C\D;
% 构造自由曲面
= meshgrid(1:0.1:5,1:0.1:5);
Z = zeros(size(X));
for i = 1:(order+1)
for j = 1:(order+1-i)
k = (order+1-i-j+1);
coef = c((i-1)*j+k);
Z = Z + coef*X.^(i-1).*Y.^(j-1);
end
end
% 绘制拟合曲面和离散点
figure;
scatter3(x,y,z,'filled','b');
hold on;
surf(X,Y,Z,'EdgeColor','none','FaceAlpha',0.5);
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Polynomial Freeform Surface Fitting');
上述代码中的过程包括：
- 构造矩阵A和向量b
- 利用加权最小二乘法求解系数向量c
- 构造拟合的自由曲面
- 绘制自由曲面和离散点的图形。
可以根据需要自行调整参数，例如多项式的阶数，以得到更好的拟合效果。
--- 光学专家Gpt
页:
[1]