B样条网格矫正相机畸变

这是我使用B样条网格去矫正相机的畸变，都是为什么我生成的网格无法覆盖整个图像。我的目标是通过这些已知点，使用B样条曲面去拟合，从而获得整个平面的分布，实现畸变矫正。

然后下一步使用多层三次B样条网格去进一步细化，实现矫正的准确性。都是现在卡在第一步了。

% 初始化相机参数
FOV = 15; % 视场角度
img_height = 1024; % 图像高度
img_width = 1024; 
% 划分的网格数
piece_u = 20;
piece_v = 20;
% 标是通过调整控制点来校正畸变坐标，使其与理论坐标一致
match_RA_DEC_points = load("match_RA_DEC_points.mat").match_RA_DEC_points;
X = groups(match_RA_DEC_points(:,8),1);
Y = groups(match_RA_DEC_points(:,9),2);% groups是对X,Y进行维度构建
Z = ones(size(X));
[M, N] = size(X);
k = 3; l = 3;
Surf_PlotSubMesh(M, N, k, l, X, Y, Z,piece_u,piece_v)

function Surf_PlotSubMesh(M, N, k, l, X, Y, Z,piece_u,piece_v)
% 绘制网格曲面，X Y Z是所绘制的M*N网格的坐标，k, l是网格的u向和v向的次数
n = N - 1;
m = M - 1;
Nik_u = zeros(n+1, 1);  % 基函数初始化
Nik_v = zeros(m+1, 1);

NodeVector_u = linspace(0, 1, n+k+2);     % 均匀B样条的u向节点矢量
u = linspace(k/(n+k+1), (n+1)/(n+k+1), piece_u);
X_M_piece = zeros(M, piece_u);    % 沿着u方向的网格，M * piece
Y_M_piece = zeros(M, piece_u);
Z_M_piece = zeros(M, piece_u);
for i = 1 : M
    for j = 1 : piece_u
        for ii = 0 : 1: n
            Nik_u(ii+1, 1) = BaseFunction(ii, k , u(j), NodeVector_u);
        end
        X_M_piece(i, j) = X(i, :) * Nik_u;
        Y_M_piece(i, j) = Y(i, :) * Nik_u;
        Z_M_piece(i, j) = Z(i, :) * Nik_u;
    end
end

NodeVector_v = linspace(0, 1, m+l+2);     % 均匀B样条的u向节点矢量
v = linspace(l/(m+l+1), (m+1)/(m+l+1), piece_v);
X_MN_piece = zeros(piece_v, piece_u);
Y_MN_piece = zeros(piece_v, piece_u);
Z_MN_piece = zeros(piece_v, piece_u);
for i = 1 : piece_u
    for j = 1 : piece_v
        for ii = 0 : 1 : m
            Nik_v(ii+1, 1) = BaseFunction(ii, l, v(j), NodeVector_v);
        end
    X_MN_piece(j, i) = Nik_v' * X_M_piece(:, i);
    Y_MN_piece(j, i) = Nik_v' * Y_M_piece(:, i);
    Z_MN_piece(j, i) = Nik_v' * Z_M_piece(:, i);
    end
end
  
% 根据得到的曲面数据点X_MN_piece... 绘制网格曲面
for j = 1 : piece_u
    for i = 1 : piece_v -1
        hold on;
        plot3([X_MN_piece(i, j) X_MN_piece(i+1, j)],...
            [Y_MN_piece(i, j) Y_MN_piece(i+1, j)],...
            [Z_MN_piece(i, j) Z_MN_piece(i+1, j)], 'b-');
    end
end
for i = 1 : piece_v
    for j = 1 : piece_u -1
        
        plot3([X_MN_piece(i, j) X_MN_piece(i, j+1)],...
            [Y_MN_piece(i, j) Y_MN_piece(i, j+1)],...
            [Z_MN_piece(i, j) Z_MN_piece(i, j+1)], 'b-');
    end
end
function Nik_u = BaseFunction(i, k , u, NodeVector)
    % 计算基函数Ni,k(u),NodeVector为节点向量
    
    if k == 0       % 0次B样条
        if (u >= NodeVector(i+1)) && (u < NodeVector(i+2))
            Nik_u = 1.0;
        else
            Nik_u = 0.0;
        end
    else
        Length1 = NodeVector(i+k+1) - NodeVector(i+1);
        Length2 = NodeVector(i+k+2) - NodeVector(i+2);      % 支撑区间的长度
        if Length1 == 0.0       % 规定0/0 = 0
            Length1 = 1.0;
        end
        if Length2 == 0.0
            Length2 = 1.0;
        end
        Nik_u = (u - NodeVector(i+1)) / Length1 * BaseFunction(i, k-1, u, NodeVector) ...
            + (NodeVector(i+k+2) - u) / Length2 * BaseFunction(i+1, k-1, u, NodeVector);
    end
end