基于MATLAB的两平面间8参数变换求解函数设计（附MATLAB函数）

最新推荐文章于 2023-08-21 17:02:05 发布

原创最新推荐文章于 2023-08-21 17:02:05 发布 · 1.9k 阅读

15 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #平面 #算法

摄影测量同时被 2 个专栏收录

3 篇文章

订阅专栏

MATLAB

3 篇文章

订阅专栏

本文详细介绍了如何使用MATLAB构建函数来求解摄影测量中两平面间的8参数变换关系。通过控制点的像素坐标和指定坐标，运用间接平差原理计算转换参数，并计算拟合数据在X轴和Y轴上的中误差（RMSE），以评估变换的精度。

在摄影测量中，为了方便数据后续解算，会存在将单幅影像中像素点坐标转换为规定平面坐标的处理。利用两平面的8参数变化，可以快速的计算出由像素点坐标解算到指定坐标的具体转换关系。本文将从构建MATLAB函数的角度，对两平面间8参数变换参数的求解进行分析。

具体算法原理可以参照我之前的文章：两平面间8参数变换参数求解简单原理解析

参数矩阵的计算

已知某一影像上控制点像素坐标数据与指定坐标系坐标数据ControlPointsCoordinate

为了方便后续计算，将其分为4个n行1列的数组x,y,X,Y，n为数组行数

x = ControlPointsCoordinate(:,1);
y = ControlPointsCoordinate(:,2);
X = ControlPointsCoordinate(:,3);
Y = ControlPointsCoordinate(:,4);
[n,~] = size(ControlPointsCoordinate);

用每一个控制点的坐标求解其对应在参数阵 $A$ 和 $A_{0}$ 的部分，并创建参数阵 $B$ 和 $P$

B = eye(8); P = eye(8);
A = zeros(2*n,8); A0 = zeros(2*n,1); 
for i = 1:n
Ai = [x(i,:) , 0 , -x(i,:).*X(i,:) , y(i,:) , 0 , -y(i,:).*X(i,:) , 1 , 0;...
      0 , x(i,:) , -x(i,:).*Y(i,:) , 0 , y(i,:) , -y(i,:).*Y(i,:) , 0 , 1];
A((2*i-1):(2*i),:) = Ai;
A0((2*i-1):(2*i),:) = [X(i);Y(i)];
end

待求参数数平差值计算

利用间接平差原理方法，设此时计算转换参数近似值 $X^{0}_{efg}$ 还未求解并视为0，计算转换参数平差值 $\widehat{X}_{efg}$

l = pinv(A)*A0;
Nbb = B' * P * B;
U = B' * P * l;
X_efg = inv(Nbb) * U;

在X轴与Y轴上的中误差计算

依据中误差计算原理，计算依照转换参数的模型拟合出的数据在指定坐标系X轴与Y轴上的中误差RMSE_X和RMSE_Y

X_Fitted = zeros(n,1);
Y_Fitted = zeros(n,1);
residual_X = zeros(n,1);
residual_Y = zeros(n,1);
RMSE_X = 0;RMSE_Y = 0;
for i =1:n
    X_Fitted(i,1) = (X_efg(1,1)*x(i,1)+X_efg(4,1)*y(i,1)+X_efg(7,1))/(X_efg(3,1)*x(i,1)+X_efg(6,1)*y(i,1)+1);
    Y_Fitted(i,1) = (X_efg(2,1)*x(i,1)+X_efg(5,1)*y(i,1)+X_efg(8,1))/(X_efg(3,1)*x(i,1)+X_efg(6,1)*y(i,1)+1);
    residual_X(i,1) = X(i,1)-X_Fitted(i,1);
    residual_Y(i,1) = Y(i,1)-Y_Fitted(i,1);
    RMSE_X = RMSE_X+(residual_X(i,1))^2/(n-1);
    RMSE_Y = RMSE_Y+(residual_Y(i,1))^2/(n-1);
end
RMSE_X = sqrt(RMSE_X);
RMSE_Y = sqrt(RMSE_Y);

整体代码

function [tarns_parameters,RMSE] = get8parameters(ControlPointsCoordinate)
% Deal with input data
x = ControlPointsCoordinate(:,1);
y = ControlPointsCoordinate(:,2);
X = ControlPointsCoordinate(:,3);
Y = ControlPointsCoordinate(:,4);
[n,~] = size(ControlPointsCoordinate);

% Build parameter matrix
B = eye(8); P = eye(8);
A = zeros(2*n,8); A0 = zeros(2*n,1); 
for i = 1:n
Ai = [x(i,:) , 0 , -x(i,:).*X(i,:) , y(i,:) , 0 , -y(i,:).*X(i,:) , 1 , 0;...
      0 , x(i,:) , -x(i,:).*Y(i,:) , 0 , y(i,:) , -y(i,:).*Y(i,:) , 0 , 1];
A((2*i-1):(2*i),:) = Ai;
A0((2*i-1):(2*i),:) = [X(i);Y(i)];
end

% Compute adjustment value
l = pinv(A)*A0;
Nbb = B' * P * B;
U = B' * P * l;
X_efg = inv(Nbb) * U;

% Compute RMSE on X-axis and Y-axis
X_Fitted = zeros(n,1);
Y_Fitted = zeros(n,1);
residual_X = zeros(n,1);
residual_Y = zeros(n,1);
RMSE_X = 0;RMSE_Y = 0;
for i =1:n
    X_Fitted(i,1) = (X_efg(1,1)*x(i,1)+X_efg(4,1)*y(i,1)+X_efg(7,1))/(X_efg(3,1)*x(i,1)+X_efg(6,1)*y(i,1)+1);
    Y_Fitted(i,1) = (X_efg(2,1)*x(i,1)+X_efg(5,1)*y(i,1)+X_efg(8,1))/(X_efg(3,1)*x(i,1)+X_efg(6,1)*y(i,1)+1);
    residual_X(i,1) = X(i,1)-X_Fitted(i,1);
    residual_Y(i,1) = Y(i,1)-Y_Fitted(i,1);
    RMSE_X = RMSE_X+(residual_X(i,1))^2/(n-1);
    RMSE_Y = RMSE_Y+(residual_Y(i,1))^2/(n-1);
end
RMSE_X = sqrt(RMSE_X);
RMSE_Y = sqrt(RMSE_Y);

% Output result
tarns_parameters = table(X_efg(1,1),X_efg(2,1),X_efg(3,1),X_efg(4,1),X_efg(5,1),X_efg(6,1),X_efg(7,1),X_efg(8,1),...
    'VariableNames',{'e1','e2','e3','f1','f2','f3','g1','g2'});
RMSE = table(RMSE_X,RMSE_Y,'VariableNames',{'RMSE X','RMSE Y'});

end

运行示例