Hi,

Deine Umsetzung ähnelt sehr der von Lehmann (2014): Transformation model selection by multiple hypotheses testing., welcher auf Andrei (2006): 3D affine coordinate transformations. verweist. Das Paper von Lehmann habe ich bzgl. der Wahl der Restriktionen kurz überflogen und bei meiner Implementierung herangezogen. Dennoch konnte ich auf die schnell keine zufriedenstellende Lösung finden, die auch bei "großen" Koordinaten numerisch stabil bleibt. Im Anhang findest Du meinen Versuch aus Octave/Matlab. Vielleicht hilft es aber weiter in dieser Diskussion.

Dennoch zwei Hinweise, die Dir ggf. bei Deinem Problem weiterhelfen können, wenn Du auch für andere Lösungswege offen bist. Zum Bestimmen der 7 Parameter der Helmert-Transformation existiert eine geschlossene Lösung, die keine Iterationen benötigt. Siehe hierzu auch den Artikel: Bestimmung der 7 Parameter einer Helmerttransformation mit Quaternionen.

Solltest Du vorhaben, Dein Modell auf mehr als 7 Parameter zu erweitern bspw. Affintransformation, beachte bitte, dass in der Geodäsie üblicherweise folgendes Modell verwendet wird:

X = T + RSMx

worin T die Translation, R die Drehung, S die Scherung und M die Matrix der Maßstäbe ist. Um dieses Modell zu reduzieren, könnte bspw. S = I gesetzt werden, sodass keine Scherungen auftreten. Die Maßstäbe in diesem Modell sind dem Startsystem zugeordnet und stehen somit nicht vor der Drehung, vgl. auch Koch & Fröhlich (2000): Transformation räumlicher variabler Koordinaten. Für dieses allg. Transformationsmodell und die nötigen Bedingungen zur Reduktion der Parameter findest Du entsprechende Gleichungen in Lösler & Eschelbach (2014): Zur Bestimmung der Parameter einer räumlichen Affintransformation.

Ich hoffe, dass hilft Dir weiter.

Viele Grüße
Micha

clear all
close all
format long g
clc
 
 
%% Punkte im Startsystem
PunkteS = [
    304.507 2262.732 8226.374
    3267.980 5992.254 4529.311
    568.348 4081.212 6567.109
];
 
%% Punkte im Zielsystem
PunkteZ = [
    275.525 1668.545 7810.831
    3239.034 5398.100 4113.819
    539.378 3487.020 6151.622
];
 
schwerpunktS = mean(PunkteS);
schwerpunktZ = mean(PunkteZ);
PunkteS = PunkteS - (schwerpunktS'*ones(1,size(PunkteS,1)))';
PunkteZ = PunkteZ - (schwerpunktZ'*ones(1,size(PunkteZ,1)))';
 
nop = size(PunkteS,1);
iter = 0;
%x=(f11 f12 f13 f21 f22 f23 f31 f32 f33)
X0=[1 0 0   0 1 0   0 0 1 ]';
dx = inf;
while(iter < 200 && max(abs(dx)) > eps^(0.75))
    iter = iter + 1;
    A = zeros(size(PunkteS,1)*size(PunkteS,2), length(X0));
    w = zeros(size(A,1),1);
    R = zeros(5, length(X0));
    r = zeros(size(R,1), 1);
 
    f11 = X0(1);
    f12 = X0(2);
    f13 = X0(3);
 
    f21 = X0(4);
    f22 = X0(5);
    f23 = X0(6);
 
    f31 = X0(7);
    f32 = X0(8);
    f33 = X0(9);
 
    for i=1:nop
        xP = PunkteS(i,1);
        yP = PunkteS(i,2);
        zP = PunkteS(i,3);
 
        XP = PunkteZ(i,1);
        YP = PunkteZ(i,2);
        ZP = PunkteZ(i,3);
 
        A(i,:)       = [xP yP zP   0  0  0    0  0  0 ];
        A(i+nop, :)  = [0  0  0   xP yP zP    0  0  0 ];
        A(i+2*nop,:) = [0  0  0    0  0  0    xP yP zP];
 
        w(i, 1)      = (f11*xP+f12*yP+f13*zP) - XP;
        w(i+nop, 1)  = (f21*xP+f22*yP+f23*zP) - YP;
        w(i+2*nop,1) = (f31*xP+f32*yP+f33*zP) - ZP;
    end
    R(1,:) = [f21 f22 f23   f11 f12 f13   0 0 0 ];
    R(2,:) = [0 0 0   f31 f32 f33   f21 f22 f23 ];
    R(3,:) = [f31 f32 f33   0 0 0   f11 f12 f13 ];
 
    r(1,1) = f11*f21 + f12*f22 + f13*f23;
    r(2,1) = f21*f31 + f22*f32 + f23*f33;
    r(3,1) = f11*f31 + f12*f32 + f13*f33;
 
    R(4,:) = [2*f11 2*f12 2*f13   -2*f21 -2*f22 -2*f23   0 0 0 ];
    R(5,:) = [-f11 -f12 -f13   -f21 -f22 -f23   2*f31 2*f32 2*f33 ];
 
    r(4,1) = f11*f11 + f12*f12 + f13*f13 - (f21*f21 + f22*f22 + f23*f23);
    r(5,1) = f31*f31 + f32*f32 + f33*f33 - 0.5*(f11*f11 + f12*f12 + f13*f13  + f21*f21 + f22*f22 + f23*f23);
 
    N = [A'*A  R'
        R     zeros(length(r))];
 
    n = -[A'*w; r];
 
    dxk = N\n;
    dx = dxk(1:length(X0));
    X0 = X0 + dx;
end
 
F     = [f11 f12 f13; f21 f22 f23; f31 f32 f33];
scale = mean([norm(F(1,:)) norm(F(2,:)) norm(F(3,:))]);
R     = 1/scale.*F;
 
% Pruefung auf orth. Matrix R'*R == I
disp(R'*R);
 
disp('Lösung')
disp(R);
disp(scale);
 

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

Tags:
Transformation, Helmert, Affin