Numerische Stabilität von Zernike-Polynomen höherer Ordnung

Ich versuche m=0, n=46Zernike-Momente höherer Ordnung (z. B. ) für ein Bild zu berechnen . Ich habe jedoch ein Problem mit dem radialen Polynom (siehe Wikipedia ). Dies ist ein Polynom, das im Intervall [0 1] definiert ist. Siehe den MATLAB-Code unten

function R = radial_polynomial(m,n,RHO)
    R = 0;
    for k = 0:((n-m)/2)        
        R = R + (-1).^k.*factorial(n-k) ...
            ./ ( factorial(k).*factorial((n+m)./2-k) .* factorial((n-m)./2-k) ) ...
            .*RHO.^(n-2.*k);
    end
end

Dies stößt jedoch offensichtlich auf numerische Probleme in der Nähe RHO > 0.9.

Ich habe versucht, es so umzugestalten, dass ich polyvaldachte, es könnte einige bessere Algorithmen hinter den Kulissen geben, aber das hat nichts gelöst. Die Konvertierung in eine symbolische Berechnung hat zwar das gewünschte Diagramm erstellt, war jedoch selbst für ein einfaches Diagramm wie das gezeigte erstaunlich langsam.

Gibt es eine numerisch stabile Methode zur Bewertung solcher Polynome höherer Ordnung?

Dies ist im folgenden Octave-Skript implementiert:

clear                                     % Tested with Octave instead of Matlab
N = 120;
n_r = 1000;
R = cell(N+1,N+1);
rho = [0:n_r]/n_r;
rho_x_2 = 2*[0:n_r]/n_r;

R{0+1,0+1} = ones(1,n_r+1);               % R^0_0  Unfortunately zero based cell indexing is not possible
R{1+1,1+1} = R{0+1,0+1}.*rho;             % R^1_1  ==>  R{...+1,...+1} etc.
for n = 2:N,
    if bitget(n,1) == 0,                  % n is even
        R{0+1,n+1} = -R{0+1,n-2+1}+rho_x_2.*R{1+1,n-1+1};                % R^0_n
        m_lo = 2;
        m_hi = n-2;
    else
        m_lo = 1;
        m_hi = n-1;
    end
    for m = m_lo:2:m_hi,
        R{m+1,n+1} = rho.*(R{m-1+1,n-1+1}+R{m+1+1,n-1+1})-R{m+1,n-2+1};  % R^m_n
    end
    R{n+1,n+1} = rho.*R{n-1+1,n-1+1};                                    % R^n_n
end;


Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);
m = 22;
n = 112;
figure
plot(rho,Z(m,n,rho))
hold on
plot(rho,R{m+1,n+1},'r');
xlabel("rho")
ylabel("R^{22}_{112}(rho)")
legend("via Jacobi","recursive");
%print -djpg plt.jpg

m = 0;
n = 46;
max_diff_m_0_n_46 = norm(Z(m,n,rho)-R{m+1,n+1},inf)

$m = 22$$n = 112$$\rho = 0.7$

$m = 0$$n = 46$1.4e-10

Eine mögliche Lösung (von @gammatester vorgeschlagen) ist die Verwendung von Jacobi-Polynomen. Dies umgeht das Problem der katastrophalen Aufhebung beim Addieren der großen Polynomkoeffizienten durch "naive" Polynombewertung.

Das radiale Zernike-Polynom kann durch Jacobi-Polynome wie folgt ausgedrückt werden (siehe Gleichung (6) )

In MATLAB ist die Verwendung von jacobiP(n,a,b,x)jedoch für große Vektoren / Matrizen von unannehmbar langsam x=rho. Die jacobiPFunktion ist tatsächlich Teil der Symbolic Toolbox, und die Auswertung des Polynoms wird auf die symbolische Engine verschoben, die Geschwindigkeit gegen willkürliche Präzision eintauscht. Eine manuelle Implementierung der Jacobi-Polynome ist daher erforderlich.

$\alpha=m$$\beta=0$$n^*=(n-m/2)$$s$

In MATLAB bedeutet dies (Jacobi p olice d epartment P olynomial, ' D ouble' Implementierung)

function P = jacobiPD(n,a,b,x)
    P = 0;
    for  s  0:n
        P = P + ...
            1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ...
            ((x-1)/2).^(n-s).*((x+1)/2).^s;
    end
    P = P*factorial(n+a) * factorial(n+b);
end

Das tatsächliche radiale Zernike-Polynom ist also (für m=abs(m))

Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);

Hinweis: Diese Selbstantwort ist nur eine praktische Lösung. Fühlen Sie sich frei, eine andere Antwort zu markieren, die erklärt, warum dies funktioniert.