星图识别3-旋转矩阵与坐标转换

这篇文章将介绍世界坐标转为星敏坐标的过程。

根据专利 https://patentimages.storage.googleapis.com/da/9a/a0/6568cfa38683e7/CN1923621A.pdf 的矩阵转换描述，再结合003.pdf文中的图来理解矩阵变换的算法过程。

（请注意，这里的R已是可以直接乘的姿态变换矩阵，为旋转矩阵的逆。）

搜索旋转矩阵的基本原理，

对本文的旋转应是左乘还是右乘？上文是绕固定轴旋转，为左乘。

再看一个定理 绕固定坐标轴旋转与绕自身坐标轴旋转一致性证明 https://blog.csdn.net/jiongjiongxia123/article/details/90236737

绕固定轴左乘的公式为

而绕自身轴右乘的公式为

和专利的对应，但差了个转置。

想了很久这个转置的问题，其实是旋转矩阵与坐标变换的关系问题。从原坐标系到新坐标系的旋转矩阵为M。原坐标系下的坐标a，新坐标系下的坐标b，若a为行向量，则aM=b，若a为列向量，则m’a=b，公式应为b=m’a，举例可知。

在这需要重述一下坐标系的定义

（1）天球坐标系。以天赤道为基圈，过春分点的时圈为主圈，春分点为主点。天球坐标系采用赤经、赤纬作为坐标量。参见附件1相关叙述。

图1 星敏感器坐标系、图像坐标系及前视投影成像示意图

（2）星敏感器坐标系xyz。以投影中心O（光轴上与感光面距离为f的点，即光心，参见图1）为坐标原点，以光轴为z轴（后面的讨论中，光轴OO’与天球面的交点记为D点），过O点平行于感光面两边的直线作为x轴和y轴。图1为星敏感器坐标系、图像坐标系及前视投影成像示意图。

（3）图像坐标系XY。以感光面的中心O’（O点在该平面上的投影点）为坐标原点，平行于感光面两边的直线为X轴和Y轴的平面坐标系，参见图1。

clc;clear;close all;
index = [1670
1477
1502
1631
1603
1453
1432
1492
1488
1648
1646
1566
1688
1655
1505
];
load("附件2  简易星表.mat");
f = 256/(tand(6));
A1 = xlsread('D:\文献\starimage\npmcm2019b\npmcm2019-B\附件3  8幅星图相关数据\xingtu01.xls') - 256;
newMatrix= [A1,repmat(f,size(A1,1),1)];
for j = 1:size(newMatrix,1)
    ll(j,1) = norm(newMatrix(j,:));
end
newNewMatrix = newMatrix./ll;

A  =  [  cosd(star_data(index,3)).*cosd(star_data(index,2))   ,   cosd(star_data(index,3)).*sind(star_data(index,2))  ,  sind(star_data(index,3))  ];
Y = [   (f  ./  (sqrt(f^2 + newMatrix(1:15,1).^2 + newMatrix(1:15,2).^2))  )];
zhixiang = inv(A'*A)*A'*Y;
E = (A*zhixiang - Y);
RMS = sqrt(mse(E));
if   zhixiang(2)/sqrt(1-zhixiang(3)^2) > 0 
    a=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    a = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end
b=asind( zhixiang(3) );


% syms fai x;
% assume(fai>=0 & fai<=2*pi);

fai = 0;
R = [
-sind(a)*cosd(fai)-cosd(a)*sind(b)*sind(fai)      cosd(a)*cosd(fai)-sind(a)*sind(b)*sind(fai)      cosd(b)*sind(fai)
 sind(a)*sind(fai)-cosd(a)*sind(b)*cosd(fai)     -cosd(a)*sind(fai)-sind(a)*sind(b)*cosd(fai)      cosd(b)*cosd(fai)
 cosd(a)*cosd(b)                                  sind(a)*cosd(b)                                  sind(b)
];

inn = 1670;
star1 =  [  cosd(star_data(inn,3)).*cosd(star_data(inn,2))   ;   cosd(star_data(inn,3)).*sind(star_data(inn,2))  ;  sind(star_data(inn,3))  ];
R * star1


subplot(2,2,1); 
scatter(star_data(index,2),star_data(index,3))
title('天球坐标系，从外向内看')

subplot(2,2,2); 
scatter(A1(:,1),A1(:,2))
title('星图1')

subplot(2,2,3); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
title('星敏坐标系，从内往外看')

subplot(2,2,4); 
scatter(A1(:,1),A1(:,2))
view(-90,90);  % view(az,el) 为当前坐标区设置相机视线的方位角(顺时针为正)和仰角。
title('星图1逆时针旋转90')

针对星图1的第一个点，即第1670个点，它的赤经赤纬为 115.828 , 28.8835，用旋转矩阵乘得到星敏坐标（0.0950，-0.0949， 0.9909）与左下图对应，第二个点，即第1477个点，赤经赤纬为105.88,29.33 ，星敏坐标为（-0.0566，-0.0889，0.9944）

这个感光面，究竟对于图像坐标系来说，是哪一面呢？不管哪一面，都不影响已经确定了的XY轴和点坐标。但我画星图1的意义就不大，不好直接比对。

但在这里，经过很长时间的考虑，如题所述的（图像坐标系与星敏坐标系的XY对应）经过验证后是错误的（倘若如上上上图的XY轴示意图和上图的文字描述，第二列和第三列没有负坐标）。

%  星图1各个坐标

clc;clear;close all;
index = [1670
1477
1502
1631
1603
1453
1432
1492
1488
1648
1646
1566
1688
1655
1505
];
load("附件2  简易星表.mat");

Attribute_Set = {'LineWidth',1.5}; 

f = 256/(tand(6));
A1 = xlsread('D:\文献\starimage\npmcm2019b\npmcm2019-B\附件3  8幅星图相关数据\xingtu01.xls') - 256;
newMatrix= [A1,repmat(f,size(A1,1),1)];

for j = 1:size(newMatrix,1)
    ll(j,1) = norm(newMatrix(j,:));
end
newNewMatrix = newMatrix./ll;

corWorld =  [  cosd(star_data(index,3)).*cosd(star_data(index,2))   ,   cosd(star_data(index,3)).*sind(star_data(index,2))  ,  sind(star_data(index,3))  ];

subplot(3,4,1); 
scatter(star_data(index,2),star_data(index,3))
title('天球坐标系，从外向内看')
axis equal;hold on;
quiver(109.6,34.5,6,0,'LineWidth',2,'Color','r');
text(109.6+6,34.5,'x');
text(109.6,34.5+6,'y');
text(115.828,28.8835,'第一点（正负）');
text(105.88,29.33,'第二点（负负）');
quiver(109.6,34.5,0,6,'LineWidth',2,'Color','g');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]);

subplot(3,4,5); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
axis equal;hold on;
title('星敏坐标系，从内往外看')
quiver(109.6,34.5,6,0,'LineWidth',2,'Color','r');
quiver(109.6,34.5,0,6,'LineWidth',2,'Color','g');
text(109.6+6,34.5,'x');
text(109.6,34.5+6,'y');
text(115.828,28.8835,'第一点（正负）');
text(105.88,29.33,'第二点（负负）');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]);

subplot(3,4,9); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
axis equal;hold on;
title('上图旋转')
quiver(109.6,34.5,6,0,'LineWidth',2,'Color','r');
quiver(109.6,34.5,0,6,'LineWidth',2,'Color','g');
text(109.6+6,34.5,'x');
text(109.6,34.5+6,'y');
text(115.828,28.8835,'第一点（正负）');
text(105.88,29.33,'第二点（负负）');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]);
view(180,90); 




subplot(3,4,6); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
axis equal;hold on;
title('星敏坐标系，旋转90')
quiver(109.6,34.5,0,6,'LineWidth',2,'Color','r');
quiver(109.6,34.5,-6,0,'LineWidth',2,'Color','g');
text(109.6-6,34.5,'y');
text(109.6,34.5+6,'x');
text(115.828,28.8835,'第一点（负负）');
text(105.88,29.33,'第二点（负正）');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]); 
view(-90,90); 

subplot(3,4,7); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
axis equal;hold on;
title('星敏坐标系，旋转180')
quiver(109.6,34.5,-6,0,'LineWidth',2,'Color','r');
quiver(109.6,34.5,0,-6,'LineWidth',2,'Color','g');
text(109.6-6,34.5,'y');
text(109.6,34.5-6,'x');
text(115.828,28.8835,'第一点（负正）');
text(105.88,29.33,'第二点（正正）');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]);
view(180,90); 



subplot(3,4,8); 
scatter(star_data(index,2),star_data(index,3))
set(gca,'XDir','reverse')%对X方向反转
axis equal;hold on;
title('星敏坐标系，旋转270')
quiver(109.6,34.5,0,-6,'LineWidth',2,'Color','r');
quiver(109.6,34.5,6,0,'LineWidth',2,'Color','g');
text(109.6+6,34.5,'y');
text(109.6,34.5-6,'x');
text(115.828,28.8835,'第一点（正正）');
text(105.88,29.33,'第二点（正负）');
xlim([109.6-8 109.6+8]);
ylim([34.5-8 34.5+8]);
view(90,90);


A  =  [  cosd(star_data(index,3)).*cosd(star_data(index,2))   ,   cosd(star_data(index,3)).*sind(star_data(index,2))  ,  sind(star_data(index,3))  ];
Y = [   (f  ./  (sqrt(f^2 + newMatrix(1:15,1).^2 + newMatrix(1:15,2).^2))  )];
zhixiang = inv(A'*A)*A'*Y;
E = (A*zhixiang - Y);
RMS = sqrt(mse(E));
if   zhixiang(2)/sqrt(1-zhixiang(3)^2) > 0 
    a=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    a = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end
b=asind( zhixiang(3) );

% syms fai x;
% assume(fai>=0 & fai<=2*pi);

fai = 90;
R = [
-sind(a)*cosd(fai)-cosd(a)*sind(b)*sind(fai)      cosd(a)*cosd(fai)-sind(a)*sind(b)*sind(fai)      cosd(b)*sind(fai)
 sind(a)*sind(fai)-cosd(a)*sind(b)*cosd(fai)     -cosd(a)*sind(fai)-sind(a)*sind(b)*cosd(fai)      cosd(b)*cosd(fai)
 cosd(a)*cosd(b)                                  sind(a)*cosd(b)                                  sind(b)
];
corSensor =  corWorld * R';

subplot(3,4,2); 
scatter(corSensor(:,1),corSensor(:,2))
% set(gca,'XDir','reverse')%对X方向反转
title('右旋90度')
axis equal;hold on;
xlim([-0.1 0.1]);
ylim([-0.1 0.1]);
quiver(0,0,0.1,0,'LineWidth',2,'Color','r');
quiver(0,0,0,0.1,'LineWidth',2,'Color','g');
text(0,0.1,'y');
text(0.1,0,'x');
text(corSensor(1,1),corSensor(1,2),'第一点（负负）');
text(corSensor(2,1),corSensor(2,2),'第二点（负正）');
% view(0,90);

fai = 180;
R = [
-sind(a)*cosd(fai)-cosd(a)*sind(b)*sind(fai)      cosd(a)*cosd(fai)-sind(a)*sind(b)*sind(fai)      cosd(b)*sind(fai)
 sind(a)*sind(fai)-cosd(a)*sind(b)*cosd(fai)     -cosd(a)*sind(fai)-sind(a)*sind(b)*cosd(fai)      cosd(b)*cosd(fai)
 cosd(a)*cosd(b)                                  sind(a)*cosd(b)                                  sind(b)
];
corSensor =  corWorld * R';

subplot(3,4,3); 
scatter(corSensor(:,1),corSensor(:,2))
% set(gca,'XDir','reverse')%对X方向反转
title('右旋180度')
axis equal;hold on;
xlim([-0.1 0.1]);
ylim([-0.1 0.1]);
quiver(0,0,0.1,0,'LineWidth',2,'Color','r');
quiver(0,0,0,0.1,'LineWidth',2,'Color','g');
text(0,0.1,'y');
text(0.1,0,'x');
text(corSensor(1,1),corSensor(1,2),'第一点（负正）');
text(corSensor(2,1),corSensor(2,2),'第二点（正正）');


fai = 270;
R = [
-sind(a)*cosd(fai)-cosd(a)*sind(b)*sind(fai)      cosd(a)*cosd(fai)-sind(a)*sind(b)*sind(fai)      cosd(b)*sind(fai)
 sind(a)*sind(fai)-cosd(a)*sind(b)*cosd(fai)     -cosd(a)*sind(fai)-sind(a)*sind(b)*cosd(fai)      cosd(b)*cosd(fai)
 cosd(a)*cosd(b)                                  sind(a)*cosd(b)                                  sind(b)
];
corSensor =  corWorld * R';

subplot(3,4,4); 
scatter(corSensor(:,1),corSensor(:,2))
% set(gca,'XDir','reverse')%对X方向反转
title('右旋270度')
axis equal;hold on;
xlim([-0.1 0.1]);
ylim([-0.1 0.1]);
quiver(0,0,0.1,0,'LineWidth',2,'Color','r');
quiver(0,0,0,0.1,'LineWidth',2,'Color','g');
text(0,0.1,'y');
text(0.1,0,'x');
text(corSensor(1,1),corSensor(1,2),'第一点（正正）');
text(corSensor(2,1),corSensor(2,2),'第二点（正负）');


subplot(3,4,[9,10]); 
scatter(A1(:,1),A1(:,2))
% set(gca,'YDir','reverse')%对y方向反转
title('星图1')
axis equal;hold on;
quiver(0,0,0,200,'LineWidth',2,'Color','g');
quiver(-256,-256,0,512,'LineWidth',2,'Color','g');
quiver(0,0,200,0,'LineWidth',2,'Color','r');
quiver(-256,-256,512,0,'LineWidth',2,'Color','r');
text(200,0,'x');
text(0,200,'y');
text(-233,233,'第一点（负正）');
text(-217,-138,'第二点（负负）');
xlim([-260 260]);
ylim([-260 260]);

subplot(3,4,[11,12]); 
scatter(A1(:,1),A1(:,2))
set(gca,'YDir','reverse')%对y方向反转
title('星图1')
axis equal;hold on;
quiver(0,0,0,200,'LineWidth',2,'Color','g');
quiver(-256,-256,0,512,'LineWidth',2,'Color','g');
quiver(0,0,200,0,'LineWidth',2,'Color','r');
quiver(-256,-256,512,0,'LineWidth',2,'Color','r');
text(200,0,'x');
text(0,200,'y');
text(-233,233,'第一点（负正）');
text(-217,-138,'第二点（负负）');
xlim([-260 260]);
ylim([-260 260]);

我又细细画了一幅图，图2是 将 图1的xyz轴滚转角fai = 90（图为从外往内看，xy轴逆时针旋转90°，相当于图顺时针旋转了90度），绕着朝向自己的拇指方向，右旋90°，相当于将图1顺时针旋转90°。依次画其他图。

根据第一点、第二点的正负可见，没有一幅图上的星点坐标值是符合星图1的描述，而图8和星图1是一致的，可以猜测出结论。

因此正确的理解应该是将天球坐标转为星敏坐标，然后将x轴反向 x’ = -x 。

星图1就是由星敏坐标经过滚转角270度旋转后得到图4，再将x轴反向 x’ = -x得到图像坐标。

如何得到图像坐标呢？

（注意，本题与上文所描述不同，连相机坐标系都不是右手坐标系）

（https://qiy.net/2020/06/09/IAP-Cam-calibration/）

f = 256/(tand(6)); % f = 2435.677300280981像素

x = f * Xc / Zc ;

该公式一步到像素坐标。

下面进行验证。

%  正向验证星图1从世界坐标到星敏坐标到像素坐标并画图

clc;clear;close all;
index = [1670
1477
1502
1631
1603
1453
1432
1492
1488
1648
1646
1566
1688
1655
1505
];
load("附件2  简易星表.mat");

f = 256/(tand(6));
A1 = xlsread('D:\文献\starimage\npmcm2019b\npmcm2019-B\附件3  8幅星图相关数据\xingtu01.xls') - 256;
newMatrix= [A1,repmat(f,size(A1,1),1)];

for j = 1:size(newMatrix,1)
    ll(j,1) = norm(newMatrix(j,:));
end
newNewMatrix = newMatrix./ll;


A  =  [  cosd(star_data(index,3)).*cosd(star_data(index,2))   ,   cosd(star_data(index,3)).*sind(star_data(index,2))  ,  sind(star_data(index,3))  ];
Y = [   (f  ./  (sqrt(f^2 + newMatrix(1:15,1).^2 + newMatrix(1:15,2).^2))  )];
zhixiang = inv(A'*A)*A'*Y;
E = (A*zhixiang - Y);
RMS = sqrt(mse(E));

if   zhixiang(2)/sqrt(1-zhixiang(3)^2) > 0 
    a=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    a = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end
b=asind( zhixiang(3) );

corWorld =  [  cosd(star_data(index,3)).*cosd(star_data(index,2))   ,   cosd(star_data(index,3)).*sind(star_data(index,2))  ,  sind(star_data(index,3))  ];


fai = 270;
R = [
-sind(a)*cosd(fai)-cosd(a)*sind(b)*sind(fai)      cosd(a)*cosd(fai)-sind(a)*sind(b)*sind(fai)      cosd(b)*sind(fai)
 sind(a)*sind(fai)-cosd(a)*sind(b)*cosd(fai)     -cosd(a)*sind(fai)-sind(a)*sind(b)*cosd(fai)      cosd(b)*cosd(fai)
 cosd(a)*cosd(b)                                  sind(a)*cosd(b)                                  sind(b)
];
corSensor =  corWorld * R';

figure(1);
scatter(corSensor(:,1),corSensor(:,2))
% set(gca,'XDir','reverse')%对X方向反转
title('右旋270度')
axis equal;hold on;
xlim([-0.1 0.1]);
ylim([-0.1 0.1]);
quiver(0,0,0.1,0,'LineWidth',2,'Color','r');
quiver(0,0,0,0.1,'LineWidth',2,'Color','g');
text(0,0.1,'y');
text(0.1,0,'x');
text(corSensor(1,1),corSensor(1,2),'第一点（正正）');
text(corSensor(2,1),corSensor(2,2),'第二点（正负）');

corSensor1=[-corSensor(:,1), corSensor(:,2:3) ];
corPic = [f .* corSensor1(:,1) ./corSensor1(:,3)  ,  f .* corSensor1(:,2) ./corSensor1(:,3) ];
figure(2);
scatter(corPic(:,1),corPic(:,2))

figure(3);
scatter(A1(:,1),A1(:,2));