在摄影测量中,为了方便数据后续解算,会存在将单幅影像中像素点坐标转换为规定平面坐标的处理。利用两平面的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);
用每一个控制点的坐标求解其对应在参数阵和
的部分,并创建参数阵
和
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
待求参数数平差值计算
利用间接平差原理方法,设此时计算转换参数近似值还未求解并视为0,计算转换参数平差值
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
运行示例