星图识别5-小结验证识别、求指向、求滚转角

总结前面几个博客，首先main4.m计算出是哪些星，再带入main10.m求出滚转角。

%  main10.m

clc;clear;close all;
index = [518
472
537
428
491
469
482
503
478
499
547
460
507
556
447
479
];
load("附件2  简易星表.mat");

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

newMatrix= [A1,repmat(f,size(A1,1),1)];
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:size(index,1),1).^2 + newMatrix(1:size(index,1),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 
    alpha=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    alpha = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end

beta=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))  ];


syms a b fai;
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)
];

x_ba = f .* ( R(1,1).*corWorld(:,1) +  R(1,2).*corWorld(:,2)  +  R(1,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));
y_ba = f .* ( R(2,1).*corWorld(:,1) +  R(2,2).*corWorld(:,2)  +  R(2,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));


A = [diff(x_ba,a) , diff(x_ba,b) , diff(x_ba,fai)];
B = [diff(y_ba,a) , diff(y_ba,b) , diff(y_ba,fai)];

iter=[alpha;beta;0];
times = 10;
iterRecord = zeros(3, times);
for i = 1 : times
    
    p = [double(subs(x_ba , {a,b,fai},{iter(1),iter(2),iter(3)}))-xingtu1(:,1); double(subs(y_ba , {a,b,fai},{iter(1),iter(2),iter(3)})) - xingtu1(:,2)];
    M = [double(subs(A , {a,b,fai},{iter(1),iter(2),iter(3)})); double(subs(B , {a,b,fai},{iter(1),iter(2),iter(3)}))];
    wow = inv((M.') * M) * (M.')*p;
    iter = iter - wow
    iterRecord(:,i) = iter;

end
x=1:1:times;
figure(1);
plot(x,iterRecord(3,:),'-or' , x,iterRecord(1,:),'-*b' , x,iterRecord(2,:),'-*b' ); %线性，颜色，标记
disp('maxp=');
disp(max(p));
RMS
%  最后看各点的差值p即可

x_real = double(subs(x_ba , {a,b,fai},{iter(1),iter(2),iter(3)}));
y_real = double(subs(y_ba , {a,b,fai},{iter(1),iter(2),iter(3)}));
figure(2);
scatter(x_real,y_real,140,'dk');
hold on;
scatter(xingtu1(:,1), xingtu1(:,2), 'red','filled');

发现星图4和星图6有所偏差。

针对星图1的点

1670
1477
1502
1631
1603
1453
1432
1492
1488
1648
1646
1566
1688
1655
1505

iter =

109.6000
34.5000
-89.9997

maxp=
0.0048 表征最大的x、y值的变差绝对值。

RMS =

1.1755e-07

星图2

518
472
537
428
491
469
482
503
478
499
547
460
507
556
447
479

iter =

43.0000
18.0000
-29.9994

maxp=
0.0054

修改main6.m得到上图

星图3

1864
1942
1825
1941
1943
1900
1722

iter =

124.9771
43.9942
-90.0708

maxp=
2.3836

星图4（指向117，-28 (后查看应该为243，-28)

3249
3346
3364
3275
3421
3319
3283
3370
3265
3261
3321
3309

iter =

61.7471
28.6652
28.7760

maxp=
441.2296

由main6看上去是旋转270得到的结果。

main8得到

很奇怪，明明一致。

经过 360-赤经修改后

iter =

243.0000
-28.0000
-90.0002

RMS =

5.8798e-08

星图5

1230
1150
1017
1033
1014
1223
1008
1208
1201

iter =

86.5836
44.4766
270.0121

maxp=
0.0058

星图6

1670
1675
1502
1631
1655
1692
1492
1576
1646
1566
1655
1505

iter =

113.7861
35.4645
280.3664

maxp=
258.0423

查看问题

逆时针旋转90度，并且参考053.pdf的F02/F05/F08无解，则圈红，猜测大致是如此。

main6中用到了错误的点，去除错误点得到赤经赤纬为112.6，34.5

%  main10.m

clc;clear;close all;
index = [1670
1675
1502
1631
1655
1692
1492
1576
1646
1566
1655
1505
];
load("附件2  简易星表.mat");

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

newMatrix= [A1,repmat(f,size(A1,1),1)];
% 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:size(index,1),1).^2 + newMatrix(1:size(index,1),2).^2))  )];
A  =  [  cosd(star_data(index([1,3,4,6,7,9,10,11,12]),3)).*cosd(star_data(index([1,3,4,6,7,9,10,11,12]),2))   ,   cosd(star_data(index([1,3,4,6,7,9,10,11,12]),3)).*sind(star_data(index([1,3,4,6,7,9,10,11,12]),2))  ,  sind(star_data(index([1,3,4,6,7,9,10,11,12]),3))  ];
Y = [   (f  ./  (sqrt(f^2 + newMatrix([1,3,4,6,7,9,10,11,12],1).^2 + newMatrix([1,3,4,6,7,9,10,11,12],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 
    alpha=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    alpha = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end
beta=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))  ];
corWorld =  [  cosd(star_data(index([1,3,4,6,7,9,10,11,12]),3)).*cosd(star_data(index([1,3,4,6,7,9,10,11,12]),2))   ,   cosd(star_data(index([1,3,4,6,7,9,10,11,12]),3)).*sind(star_data(index([1,3,4,6,7,9,10,11,12]),2))  ,  sind(star_data(index([1,3,4,6,7,9,10,11,12]),3))  ];


syms a b fai;
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)
];

x_ba = f .* ( R(1,1).*corWorld(:,1) +  R(1,2).*corWorld(:,2)  +  R(1,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));
y_ba = f .* ( R(2,1).*corWorld(:,1) +  R(2,2).*corWorld(:,2)  +  R(2,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));


AA = [diff(x_ba,a) , diff(x_ba,b) , diff(x_ba,fai)];
B = [diff(y_ba,a) , diff(y_ba,b) , diff(y_ba,fai)];

iter=[alpha;beta;0];
times = 10;
iterRecord = zeros(3, times);
for i = 1 : times
    
    p = [double(subs(x_ba , {a,b,fai},{iter(1),iter(2),iter(3)}))-xingtu1([1,3,4,6,7,9,10,11,12],1); double(subs(y_ba , {a,b,fai},{iter(1),iter(2),iter(3)})) - xingtu1([1,3,4,6,7,9,10,11,12],2)];
    disp(max(p));
    M = [double(subs(AA , {a,b,fai},{iter(1),iter(2),iter(3)})); double(subs(B , {a,b,fai},{iter(1),iter(2),iter(3)}))];
    wow = inv((M.') * M) * (M.')*p;
    iter = iter - wow
    iterRecord(:,i) = iter;

end
x=1:1:times;
figure(1);
plot(x,iterRecord(3,:),'-or' , x,iterRecord(1,:),'-*b' , x,iterRecord(2,:),'-*b' ); %线性，颜色，标记
disp('maxp=');
disp(max(p));
RMS
%  最后看各点的差值p即可

x_real = double(subs(x_ba , {a,b,fai},{iter(1),iter(2),iter(3)}));
y_real = double(subs(y_ba , {a,b,fai},{iter(1),iter(2),iter(3)}));
figure(2);
scatter(x_real,y_real,140,'dk');
hold on;
scatter(xingtu1([1,3,4,6,7,9,10,11,12],1), xingtu1([1,3,4,6,7,9,10,11,12],2), 'red','filled');

iter =

112.6000
34.5000
-89.9997

maxp=
0.0048

RMS =

4.6872e-08

星图7

1525
1572
1443
1748
1780
1675
1720
1503
1634
1577
1757
1586
1536
1610
1681
1606
1670
1477
1790
1502
1631
1603
1415
1692
1453
1432
1492
1488
1648
1646
1566
1688
1655
1505
1373
1576
1545
1424
1375
1825
1401

iter =

112.6000
34.5000
269.9999

maxp=
0.0060

星图8

1572
1387
1675
1634
1586
1610
1354
1681
1606
1670
1390
1502
1631
1603
1692
1432
1488
1648
1646
1566
1688
1359
1505
1373
1576
1424
1375
1401
1385

iter =

106.0000
35.0000
270.0000

maxp=
0.0081

针对星图4，需要考察一下问题。

正向验证

%  main11.m
%  正向求解指向区域内的星图和题图的对比。

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;
xingtu1=[-A1(:,1), A1(:,2) ];

newMatrix= [A1,repmat(f,size(A1,1),1)];
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:size(index,1),1).^2 + newMatrix(1:size(index,1),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 
    alpha=acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
else 
    alpha = 360-acosd( zhixiang(1)  /  sqrt(1-((zhixiang(3))^2) )   );
end
beta=asind( zhixiang(3) );
corWorld =  [  cosd(star_data(:,3)).*cosd(star_data(:,2))   ,   cosd(star_data(:,3)).*sind(star_data(:,2))  ,  sind(star_data(:,3))  ];


% syms a b fai;
a = alpha;
b = beta;
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)
];

x_ba = f .* ( R(1,1).*corWorld(:,1) +  R(1,2).*corWorld(:,2)  +  R(1,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));
y_ba = f .* ( R(2,1).*corWorld(:,1) +  R(2,2).*corWorld(:,2)  +  R(2,3).*corWorld(:,3)) ./ (R(3,1).*corWorld(:,1) +  R(3,2).*corWorld(:,2)  +  R(3,3).*corWorld(:,3));

indexIn = find(x_ba(:,1)<256 & x_ba(:,1)>-256 & y_ba(:,1)<256 & y_ba(:,1)>-256);
idx = kmeans(indexIn,2);
idx1 = find(idx == 1);
idx2 = find(idx == 2);
all(ismember(index,indexIn))

CorIn = [x_ba(indexIn(idx2)) , y_ba(indexIn(idx2))];
figure(1);  
plot(CorIn(:,1),CorIn(:,2),'ok'  ); %线性，颜色，标记
% plot(x,iterRecord(3,:),'-or' , x,iterRecord(1,:),'-*b' , x,iterRecord(2,:),'-*b' ); %线性，颜色，标记
text(CorIn(:,1)-20,CorIn(:,2)-17,string(indexIn(idx2)));
hold on
scatter(xingtu1(:,1),xingtu1(:,2),140,'or'  ); %线性，颜色，标记

figure(2);
plot(star_data(:,2),star_data(:,3),'ok'  ); %线性，颜色，标记

xlim([109.6-6  109.6+6]);
ylim([34.5-6  34.5+6]);
text(star_data(:,2),star_data(:,3),string(star_data(:,1)) );


figure(3);
scatter3(corWorld(indexIn,1),corWorld(indexIn,2),corWorld(indexIn,3));