二维坐标系空间变换(详细解读,附MATLAB代码)

简介: 二维坐标系空间变换(详细解读,附MATLAB代码)

二维坐标系空间变换

参考链接;

代码资源;

假如存在任意两个二维坐标系,如下图所示:

目的:将xoy坐标系经过处理变换到XOY坐标系。

经过分析可知:对于二维平面上的变换需要 x y 方向上两个平移参数dx,dy,以及一个绕着垂直于xy平面的旋转参数θ,此外,还存在一个缩放尺度因子m(可选)。

求解出以上四个参数就能够实现任意两个二维坐标系之间的相互转换,转换公式如下。

下面分两种方法求解以上四个参数:

1 直接法

前提:已知两对公共点分别在两个坐标系下的坐标,如上图中待转换坐标系中的o,p 和在参考坐标系中的O,P。

求解过程

1)对于平移参数,简单认为 第一对公共点(如o与O)的差值就是所有点的平移量,即 O 点在两个坐标系下的坐标值直接相减得到dx,dy。

2)对于旋转参数,认为op连线在xoy坐标系与XOY坐标系下的方位角差值为旋转角θ。

首先计算线段op在两个坐标系下的向量表示:

基于上述向量,计算方位角:

注意:方位角A和α计算时要区分不同象限,判断△x,△y的正负。这里简单附上代码,就不过多解释了

azimuth.m

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);
elseif(dy >0 && dx < 0) % 第二象限 钝角 90-180°,180°- XX
    a_orb = pi - atan(dy / abs(dx));
elseif(dy <0 && dx < 0) % 第三象限 180-270° 180°+ XX
    a_orb = pi + atan(dy / dx);
elseif(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
close all
%% 任选两对公共点
xoy = load('xoy_Cor.txt'); % 读取坐标,xoy_Cor 待转换坐标系
XOY_C = load('XOY_C_Cor.txt'); % XOY_C数据,参考坐标系
%% 计算平移量
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
XOY_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_Coor.txt',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矩阵中的两个元素;

4)如果未知数的初值不够准确,则需要迭代计算,即将解出的参数X

加到对应的初值上,形成新的初值,重复上述迭代计算,直到参数X的改正数接近于0,或一个可以接受的阈值,结束迭代。

代码如下:(init_cal函数如下)(azimuth函数在上方)

function [trans_dx,trans_dy, xita] = init_cal(xoy_p1,xoy_p2,XOY_C_p1,XOY_C_p2)
%INIT_CAL xiaochen 2021/07/17
%   两点法计算初值
%   先通过两点法计算转换参数,为后续进行优化提供初值
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);
xita =  a_XOY_C - a_xoy; % 旋转角度
trans_dx = dx; % x平移
trans_dy = dy; % y平移
end

最小二乘法求解脚本

%% xiaochen 2021/07/17
% 二维坐标转换  3参数计算 最小二乘优化
clc
clear
close all
%% 读取数据
xoy = load('xoy_Cor.txt'); % 读取坐标,xoy_Cor.txt为翻转180度后的数据
XOY_C = load('XOY_C_Cor.txt'); % XOY_C数据,参考坐标系
xoy_p1 = [xoy(100,1), xoy(100,2)]; % xoy 点1
xoy_p2 = [xoy(150,1), xoy(150,2)]; % xoy 点2
XOY_C_p1 = [XOY_C(100,1), XOY_C(100,2)]; % XOY_C 点1
XOY_C_p2 = [XOY_C(150,1), XOY_C(150,2)]; % XOY_C 点2
[init_x, init_y, init_xita] = init_cal(xoy_p1,xoy_p2,XOY_C_p1,XOY_C_p2); % 直接法计算参数初值
%% 最小二乘法优化 
xita = init_xita;
dx = init_x;
dy = init_y;
m = 0; % 参数赋初值
for i = 1:length(xoy)
   B = [1,0,cos(xita)*xoy(i,1)+sin(xita)*xoy(i,2),(1+m)*(-sin(xita)*xoy(i,1)+cos(xita)*xoy(i,2));...
       0,1, -sin(xita)*xoy(i,1)+cos(xita)*xoy(i,2),(1+m)*(-cos(xita)*xoy(i,1)-sin(xita)*xoy(i,2))];
   rotation = [cos(xita), sin(xita); -sin(xita), cos(xita)];
   xoy_trans = (1+m)*rotation*[xoy(i,1);xoy(i,2)] + [dx; dy];
   L = [XOY_C(i,1); XOY_C(i,2)] - xoy_trans;
   BB = (B'*B)\B'; % (B'*B)\B' 等同于 inv(B'*B)* B'
   para = BB * L;
   if (abs(para(1)) < 0.01 && abs(para(2)) < 0.01 && abs(para(3)) < 0.01 && abs(para(4)) < 0.01)
       fprintf('第 %d 次循环,达到收敛。\n', i);
       break;
   else
       dx = dx + para(1);
       dy = dy + para(2);
       m = m + para(3);
       xita = xita + para(4); % 改进参数
   end
   fprintf('第 %d 次循环。\n', i);
end
%% 转换数据并输出
fprintf('translation: [%f, %f], rotation: %f° \n', dx, dy, rad2deg(xita));
disp('the rotation matrix is:');
disp(rotation);
xoy_aft = [xoy(:,1),xoy(:,2)]* rotation  + [dx, dy]; % 使用计算的三个参数进行变换,没有使用尺度参数m
xoy_aft = [xoy_aft, xoy(:,3)]; % 与最后一列高程进行组合
dlmwrite('xoy_Opti.txt',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),'.-');
目录
相关文章
|
1月前
|
算法 5G 数据安全/隐私保护
MIMO系统中差分空间调制解调matlab误码率仿真
本项目展示了一种基于Matlab 2022a的差分空间调制(Differential Space Modulation, DMS)算法。DMS是一种应用于MIMO通信系统的信号传输技术,通过空间域的不同天线传输符号序列,并利用差分编码进行解调。项目包括算法运行效果图预览、核心代码及详细中文注释、理论概述等内容。在发送端,每次仅激活一个天线发送符号;在接收端,通过差分解调估计符号和天线选择。DMS在快速衰落信道中表现出色,尤其适用于高速移动和卫星通信系统。
|
1月前
|
存储 数据可视化 数据挖掘
使用Matlab绘制简单的二维与三维图形
【10月更文挑战第3天】本文详细介绍了如何在 Matlab 中绘制简单的二维和三维图形,包括曲线图、柱状图、散点图、网格图、表面图、等高线图、多边形填充图、切片图及矢量场等。文章提供了丰富的代码示例,如使用 `plot`、`bar`、`scatter`、`plot3`、`mesh`、`surf`、`contour` 等函数绘制不同类型图形的方法,并介绍了 `rotate3d`、`comet3` 和 `movie` 等工具实现图形的交互和动画效果。通过这些示例,读者可以轻松掌握 Matlab 的绘图技巧,并应用于数据可视化和分析中。
44 6
|
2月前
|
算法 5G 数据安全/隐私保护
SCM信道模型和SCME信道模型的matlab特性仿真,对比空间相关性,时间相关性,频率相关性
该简介展示了使用MATLAB 2022a进行无线通信信道仿真的结果,仿真表明信道的时间、频率和空间相关性随间隔增加而减弱,并且宏小区与微小区间的相关性相似。文中介绍了SCM和SCME模型,分别用于WCDMA和LTE/5G系统仿真,重点在于其空间、时间和频率相关性的建模。SCME模型在SCM的基础上进行了扩展,提供了更精细的参数化,增强了模型的真实性和复杂度。最后附上了MATLAB核心程序,用于计算不同天线间距下的空间互相关性。
76 0
|
3月前
|
安全
【2023高教社杯】D题 圈养湖羊的空间利用率 问题分析、数学模型及MATLAB代码
本文介绍了2023年高教社杯数学建模竞赛D题的圈养湖羊空间利用率问题,包括问题分析、数学模型建立和MATLAB代码实现,旨在优化养殖场的生产计划和空间利用效率。
198 6
【2023高教社杯】D题 圈养湖羊的空间利用率 问题分析、数学模型及MATLAB代码
|
3月前
|
存储 算法 搜索推荐
【2022年华为杯数学建模】B题 方形件组批优化问题 方案及MATLAB代码实现
本文提供了2022年华为杯数学建模竞赛B题的详细方案和MATLAB代码实现,包括方形件组批优化问题和排样优化问题,以及相关数学模型的建立和求解方法。
128 3
【2022年华为杯数学建模】B题 方形件组批优化问题 方案及MATLAB代码实现
|
3月前
|
数据采集 存储 移动开发
【2023五一杯数学建模】 B题 快递需求分析问题 建模方案及MATLAB实现代码
本文介绍了2023年五一杯数学建模竞赛B题的解题方法,详细阐述了如何通过数学建模和MATLAB编程来分析快递需求、预测运输数量、优化运输成本,并估计固定和非固定需求,提供了完整的建模方案和代码实现。
90 0
【2023五一杯数学建模】 B题 快递需求分析问题 建模方案及MATLAB实现代码
|
6月前
|
算法 数据安全/隐私保护 C++
基于二维CS-SCHT变换和扩频方法的彩色图像水印嵌入和提取算法matlab仿真
该内容是关于一个图像水印算法的描述。在MATLAB2022a中运行,算法包括水印的嵌入和提取。首先,RGB图像转换为YUV格式,然后水印通过特定规则嵌入到Y分量中,并经过Arnold置乱增强安全性。水印提取时,经过逆过程恢复,使用了二维CS-SCHT变换和噪声对比度(NC)计算来评估水印的鲁棒性。代码中展示了从RGB到YUV的转换、水印嵌入、JPEG压缩攻击模拟以及水印提取的步骤。
|
6月前
|
数据安全/隐私保护
耐震时程曲线,matlab代码,自定义反应谱与地震波,优化源代码,地震波耐震时程曲线
地震波格式转换、时程转换、峰值调整、规范反应谱、计算反应谱、计算持时、生成人工波、时频域转换、数据滤波、基线校正、Arias截波、傅里叶变换、耐震时程曲线、脉冲波合成与提取、三联反应谱、地震动参数、延性反应谱、地震波缩尺、功率谱密度
基于混合整数规划的微网储能电池容量规划(matlab代码)
基于混合整数规划的微网储能电池容量规划(matlab代码)
|
6月前
|
算法 调度
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)