根据一组已知旧坐标系和新坐标系坐标推求七参数并转换未知旧坐标
%不同空间直角坐标系转换
X=input("请输入旧坐标X值");
Y=input("请输入旧坐标Y值");
Z=input("请输入旧坐标Z值");
P0=[X;Y;Z];
%计算七参数
%已知旧坐标p1
p1=[ -1964734.9635 4484768.5466 4075386.7697;
-1967174.8023 4490401.5079 4067948.1663;
-1958198.6099 4481934.1931 4081954.0888;
-1958489.2289 4485256.3985 4077866.1343;
-1953456.7885 4481362.6794 4084841.9629;
-1958020.9950 4492625.1514 4069911.5369];
%已知新坐标p2
p2=[-1964642.8359 4484908.5860 4075486.8981;
-1967082.7160 4490541.6460 4068048.1509;
-1958106.3701 4482074.1789 4082054.3211;
-1958396.9949 4485396.4450 4077966.2971;
-1953364.4591 4481502.6550 4084942.2649;
-1957928.7551 4492765.3050 4070011.5630];
B=-[1 0 0 p1(1,1) 0 -p1(1,3) p1(1,2) ;%已知旧坐标矩阵
0 1 0 p1(1,2) p1(1,3) 0 -p1(1,1) ;
0 0 1 p1(1,3) -p1(1,2) p1(1,1) 0 ;
1 0 0 p1(2,1) 0 -p1(2,3) p1(2,2) ;
0 1 0 p1(2,2) p1(2,3) 0 -p1(2,1) ;
0 0 1 p1(2,3) -p1(2,2) p1(2,1) 0 ;
1 0 0 p1(3,1) 0 -p1(3,3) p1(3,2) ;
0 1 0 p1(3,2) p1(3,3) 0 -p1(3,1) ;
0 0 1 p1(3,3) -p1(3,2) p1(3,1) 0 ;
1 0 0 p1(4,1) 0 -p1(4,3) p1(4,2) ;
0 1 0 p1(4,2) p1(4,3) 0 -p1(4,1) ;
0 0 1 p1(4,3) -p1(4,2) p1(4,1) 0 ;
1 0 0 p1(5,1) 0 -p1(5,3) p1(5,2) ;
0 1 0 p1(5,2) p1(5,3) 0 -p1(5,1) ;
0 0 1 p1(5,3) -p1(5,2) p1(5,1) 0 ];
L=[p2(1,1) ; p2(1,2) ; p2(1,3);%已知新坐标L矩阵
p2(2,1) ; p2(2,2) ; p2(2,3);
p2(3,1) ; p2(3,2) ; p2(3,3);
p2(4,1) ; p2(4,2) ; p2(4,3);
p2(5,1) ; p2(5,2) ; p2(5,3)];
P=eye(15);
EX=-inv((B.'*P*B))*B.'*P*L;%参数矩阵
detaX=EX(1);detaY=EX(2);detaZ=EX(3);%平移参数
m=EX(4)-1;%尺度参数
eX=EX(5)/EX(4);eY=EX(6)/EX(4);eZ=EX(7)/EX(4);
P1=(1+m)*[1 eZ -eY;-eZ 1 eX;eY -eX 1]*P0+[detaX;detaY;detaZ];
disp("尺度变化参数为"),m;
disp("平移参数为")
[detaX detaY detaZ]
disp("旋转参数为")
[eX eY eZ]
disp("新坐标为");
P1.'