⼆维坐标系空间变换（详细解读，附MATLAB代码）
⼆维坐标系空间变换
假如存在任意两个⼆维坐标系，如下图所⽰：
⽬的：将xoy坐标系经过处理变换到XOY坐标系。
经过分析可知：对于⼆维平⾯上的变换需要 x y ⽅向上两个平移参数dx，dy，以及⼀个绕着垂直于xy平⾯的旋转参数θ，此外，还存在⼀个缩放尺度因⼦m（可选）。
求解出以上四个参数就能够实现任意两个⼆维坐标系之间的相互转换，转换公式如下。
下⾯分两种⽅法求解以上四个参数：
1 直接法
福尔摩斯 演绎法前提：已知两对公共点分别在两个坐标系下的坐标，如上图中待转换坐标系中的o,p 和在参考坐标系中的O,P。
求解过程：
1）对于平移参数，简单认为 第⼀对公共点（如o与O）的差值就是所有点的平移量，即 O 点在两个坐标系下的坐标值直接相减得到dx，dy。
2）对于旋转参数，认为op连线在xoy坐标系与XOY坐标系下的⽅位⾓差值为旋转⾓θ。
spend⾸先计算线段op在两个坐标系下的向量表⽰：
基于上述向量，计算⽅位⾓：
注意：⽅位⾓A和α计算时要区分不同象限，判断△x,△y的正负。这⾥简单附上代码，就不过多解释了。
azimuth.m
global warming
function[outputArg]= azimuth(dx,dy)神探夏洛克第五季
%AZIMUTH xiaochen 2021/07/17
% ⽅位⾓计算（没有考虑dx,dy为0的情况）
% 通过dx 与 dy，计算⽅位⾓判断属于哪个象限
if(dy >0&& dx >0) % 第⼀象限锐⾓0-90°，直接求解
a_orb = atan(dy / dx);
elif(dy >0&& dx <0) % 第⼆象限钝⾓90-180°，180°- XX
a_orb = pi - atan(dy / abs(dx));
ttlementelif(dy <0&& dx <0) % 第三象限180-270° 180°+ XX
a_orb = pi + atan(dy / dx);
elif(dy <0&& dx >0) % 第四象限270-360° 360°- XX
a_orb =2*pi - atan(abs(dy) / dx);
end
outputArg = a_orb;
end
3）对于尺度参数，认为op连线在xoy坐标系与XOY坐标系下的长度的⽐值即为对应的缩放参数。长度计算：
参数m计算：
⾄此，我们就通过直接法计算出了⼆维坐标系转换所有的4个转换参数。
代码如下：（azimuth函数在上⽅）
%% xiaochen 2021/07/17
% ⼆维坐标转换3参数计算直接参数法，⽆最⼩⼆乘优化
clc
clear
clo all
%% 任选两对公共点
xoy = load(''); % 读取坐标，xoy_Cor 待转换坐标系登船
XOY_C = load('XOY_'); % XOY_C数据，参考坐标系
2010年12月四级真题%% 计算平移量
xoy_p1 =[xoy(1,1), xoy(1,2)]; % xoy 点1
xoy_p2 =[xoy(150,1), xoy(150,2)]; % xoy 点2
XOY_C_p1 =[XOY_C(1,1), XOY_C(1,2)]; % XOY_C 点1
ugg是什么意思netantXOY_C_p2 =[XOY_C(150,1), XOY_C(150,2)]; % XOY_C 点2
%% 计算平移量
dx = XOY_C_p1(1) - xoy_p1(1);
dy = XOY_C_p1(2) - xoy_p1(2);
%% 计算旋转⾓度
xoy_dy = xoy_p2(2)-xoy_p1(2); xoy_dx = xoy_p2(1)-xoy_p1(1);
XOY_C_dy = XOY_C_p2(2)-XOY_C_p1(2); XOY_C_dx = XOY_C_p2(1)-XOY_C_p1(1);
% xoy⽅位⾓计算判断
a_xoy = azimuth(xoy_dx, xoy_dy);
% XOY_C⽅位⾓计算判断
a_XOY_C = azimuth(XOY_C_dx, XOY_C_dy);
% a_xoy = atan(xoy_dy / xoy_dx);
% a_XOY_C = atan(XOY_C_dy / XOY_C_dx);
xita = a_XOY_C - a_xoy; % 旋转⾓度
rotation =[cos(xita), sin(xita); -sin(xita), cos(xita)]; % 旋转矩阵
xoy_aft =[xoy(:,1),xoy(:,2)]* rotation + [dx, dy]; % 使⽤计算的三个参数进⾏变换
xoy_aft =[xoy_aft, xoy(:,3)]; % 与最后⼀列⾼程进⾏组合
dlmwrite('',xoy_aft,'delimiter','\t','precision',10); % 写出数据⽂件
%% 查看效果
figure;
plot(XOY_C(:,1),XOY_C(:,2),'.-');
hold on;
plot(xoy(:,1),xoy(:,2),'.-');
hold on;
plot(xoy_aft(:,1),xoy_aft(:,2),'.-');
% legend('XOY_C', 'xoy', 'xoyaft');
fprintf('translation: [%f, %f], rotation: %f° \n', dx, dy, rad2deg(xita));
disp('the rotation matrix is:');
disp(rotation);
disp('finish!');
%% 两点之间距离，计算尺度时使⽤，这⾥没有进⾏使⽤
s_xoy = sqrt(power((xoy_p1(1)-xoy_p2(1)),2) + power((xoy_p1(2)-xoy_p2(2)),2));
s_XOY_C = sqrt(power((XOY_C_p1(1)-XOY_C_p2(1)),2) + power((XOY_C_p1(2)-XOY_C_p2(2)),2)); scale =(s_XOY_C - s_xoy)/s_xoy;
2 最⼩⼆乘优化法
对于多公共点（已知多于两对公共点分别在两个坐标系下的坐标）的坐标变换，需利⽤最⼩⼆乘原理，将模型参数⼀起作为未知参数列⽴误差⽅程进⾏解算。
转换公式仍与上述⽅法相同：
写成线性公式的形式：
现金流量表的编制方法
每⼀个（x，y）可由上式计算得（X’，Y’），将该点在XOY坐标系下的（X，Y）当做观测值，则有：
以dx、dy、m、θ为未知数，将以上两式线性化展开，保留⼀次项，可得：
上式写成矩阵形式：
求解参数 X：
编程实现过程：
1）未知数近似值的求解，可选取任意两点，采⽤上述直接法求得；
2）B矩阵元素的计算，采⽤1）求得的近似值，代⼊线性化后的误差⽅程式计算系数项，每对公共点可列两个误差⽅程式，可计算得B矩阵中的两⾏，如下图；
3）L矩阵元素的计算，采⽤1）求得的近似值，将（x，y）计算出（X’，Y），⽤公共点坐标减去计算坐标即得L矩阵，每对公共点可计算得L矩阵中的两个元素；
