Hej David,

In deinem Dokument in Abschnitt 2 ist mir bis zur symetrischen Matrix N auch alles klar. Aber bei "..zum grössten Eigenwert gehörenden Eigenvektor der symetrischen Matrix N" komm ich nicht so ganz weiter. Wie komme ich über N zu meinen Unbekannten der Quaternion q0, q1, q2 und q3?

Wenn Du die Singulärwertzerlegung auf die Matrix N anwendest, erhältst Du die Singulärwerte und die Singulärvektoren. Alternativ kannst Du auch eine Eigenwertzerlegung in diesem Fall durchführen und würdest die Eigenwerte und -vektoren bekommen. Zu jedem Eigenwert existier also ein Eigenvektor. Nun suchst Du den größten Eigenwert von N. Der korrespondierende Eigenvektor ist dann die Quaternion.

Hier mal ein Matlab-Code, der Dir sicher weiterhilft und auch in Octave lauffähig sein sollte.

%%
% Direkte Bestimmung der Transformationsparameter einer raeumlichen
% Helmert-Transformation mittels Quaternion zur Beschreibung
% der Rotationssequenz.
%
% Michael Loesler, 16.06.2018
 
clear all;
close all;
format long g;
clc;
 
%% Punkte im Startsystem
PunkteS = [
4157222.543 664789.307 4774952.099
4149043.336 688836.443 4778632.188
4172803.511 690340.078 4758129.701
4177148.376 642997.635 4760764.800
4137012.190 671808.029 4791128.215
4146292.729 666952.887 4783859.856
4138759.902 702670.738 4785552.196
];
 
%% Punkte im Zielsystem
PunkteZ = [
4157870.237 664818.678 4775416.524
4149691.049 688865.785 4779096.588
4173451.354 690369.375 4758594.075
4177796.064 643026.700 4761228.899
4137659.549 671837.337 4791592.531
4146940.228 666982.151 4784324.099
4139407.506 702700.227 4786016.645
];
 
%% Ermittlung des Schwerpunktes beider System
cX = mean(PunkteS);
cY = mean(PunkteZ);
 
%% Schwerpunktreduktion
X = PunkteS - (cX'*ones(1,size(PunkteS,1)))';
Y = PunkteZ - (cY'*ones(1,size(PunkteZ,1)))';
 
S = X'*Y;
 
Sxx = S(1,1);
Sxy = S(1,2);
Sxz = S(1,3);
 
Syx = S(2,1);
Syy = S(2,2);
Syz = S(2,3);
 
Szx = S(3,1);
Szy = S(3,2);
Szz = S(3,3);
 
N = [ Sxx+Syy+Szz      Syz-Szy         Szx-Sxz       Sxy-Syx
        Syz-Szy      Sxx-Syy-Szz       Sxy+Syx       Szx+Sxz
        Szx-Sxz        Sxy+Syx      -Sxx+Syy-Szz     Syz+Szy
        Sxy-Syx        Szx+Sxz         Syz+Szy    -Sxx-Syy+Szz];
 
% groesster positiver Eigenwert suchen und Eigenvektor nehmen == Rotationsquaternion
[V,D] = eig(N);
[maxD, idxD] = max(diag(D));
%% Quaternion der Drehung
q = V(:,idxD);
 
q0 = q(1);
qq = q(2:end);
 
Cq = [  0    -qq(3)  qq(2)
       qq(3)   0    -qq(1)
      -qq(2)  qq(1)   0   ];
 
%% Ableitung der Rotationsmatrix R
R = (q0^2 - qq'*qq)*eye(3) + 2*(qq*qq' + q0*Cq );
 
 
%% Maßstab
xTx = trace(X' * X);
m = maxD/xTx;
T = cY' - m*R*cX';
 
disp('Massstab');
disp(m);
 
disp('Translation');
disp(T);
 
disp('Rotationsmatrix');
disp(R);
 
disp('Quaternion');
disp(q);
 
disp('Omgea');
V = (m * R * PunkteS' + T * ones(1,size(PunkteS,1)))' - PunkteZ;
v = V(:);
disp(v'*v);

Viele Grüße
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Matlab, Helmert-Transformation, Singulärwertzerlegung, Quaternion, Eigenwertzerlegung