Próbuję obliczyć wyższego rzędu (np m=0
, n=46
) chwile Zernike jakiegoś obrazu. Mam jednak problem dotyczący wielomianu promieniowego (patrz wikipedia ). Jest to wielomian zdefiniowany w przedziale [0 1]. Zobacz kod MATLAB poniżej
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
To oczywiście napotyka na problemy numeryczne w pobliżu RHO > 0.9
.
Próbowałem przeredagować go, tak aby polyval
myślał, że może mieć lepsze algorytmy za kulisami, ale to niczego nie rozwiązało. Przekształcenie go w obliczenie symboliczne stworzyło pożądany wykres, ale było zadziwiająco powolne, nawet w przypadku prostego wykresu, takiego jak pokazano.
Czy istnieje stabilny numerycznie sposób oceny wielomianów wysokiego rzędu?
matlab
polynomials
numerical-limitations
Sanchises
źródło
źródło
Odpowiedzi:
W tym artykule Honarvar i Paramesran opracowali interesującą metodę obliczenia wielomianów promieniowych Zernike w bardzo przyjemny sposób rekurencyjny. Formuła rekurencyjna jest zaskakująco prosta, bez dzielenia lub mnożenia przez duże liczby całkowite:Rmn(ρ)=ρ(R|m−1|n−1(ρ)+Rm+1n−1(ρ))−Rmn−2(ρ)
Poleciłbym rzucić okiem na rysunek 1 w pracy Honarvara i Paramesrana, który wyraźnie ilustruje zależności między różnymi wielomianami Zernike.
Jest to zaimplementowane w następującym skrypcie Octave:
Na przykład liczba wygenerowana przez ten kod pokazuje, że za pomocąm=22 , i n=112 , katastrofalne anulowanie następuje w pobliżu ρ=0.7 , jeśli wielomiany promieniowe Zernike są obliczane za pomocą wielomianów Jacobiego. Dlatego należy również martwić się o dokładność wielomianów Zernike niższego stopnia.
Metoda rekurencyjna wydaje się znacznie bardziej odpowiednia do stabilnego obliczania wielomianów zerowych wyższego rzędu. Niemniej jednak dlam=0 i n=46 , maksymalna różnica między metodą Jacobi a metodą rekurencyjną wynosi (tylko?)
1.4e-10
, co może być wystarczająco dokładne dla Twojej aplikacji.źródło
jacobiPD
, a nie jak jakakolwiek katastroficzna anulacja.JacobiPD
z jego odpowiedzi . Działa to dobrze w przypadku wielomianów niskiego rzędu. Na przykład za pomocą6.9e-13
. Chociaż poszczególne warunki w podsumowaniuJacobiPD
są małe, mogą się zwiększyć po pomnożeniu przezfactorial(n+a) * factorial(n+b)
. Ponadto mają naprzemienny znak, który jest idealnym przepisem na katastrofalne anulowanie.1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ((x-1)/2).^(n-s).*((x+1)/2).^s * factorial(n+a) * factorial(n+b)
może stać się tak duże, jak1.4e18
, a suma jest dopiero w-2.1
końcu. Możesz nazwać to błędem, ale z nieskończoną precyzją odpowiedź byłaby poprawna. Czy możesz wyjaśnić, co rozumiesz przez „brak ogólnej katastrofalnej rezygnacji”?Możliwym rozwiązaniem (sugerowanym przez @gammatester) jest użycie wielomianów Jacobi. To omija problem katastroficznego anulowania przy dodawaniu dużych współczynników wielomianowych poprzez „naiwną” ocenę wielomianu.
Promieniowy wielomian Zernike może być wyrażony przez wielomiany Jacobiego w następujący sposób (patrz równanie (6) )
Jednak w MATLAB użycie
jacobiP(n,a,b,x)
jest niedopuszczalnie powolne w przypadku dużych wektorów / macierzyx=rho
. TajacobiP
funkcja jest właściwie częścią Symbolicznego zestawu narzędzi, a ocena wielomianu jest odraczana do silnika symbolicznego, który zamienia prędkość na dowolną precyzję. Konieczna jest zatem ręczna implementacja wielomianów Jacobi.Ponieważ wszystkie parametry funkcji Jacobi są nieujemne (α=m , β=0 , n∗=(n−m/2) ), możemy użyć następującego wyrażenia (patrz Wikipedia , zauważ, że wpisałem wartości dlas )
P(α,β)n(ρ)=(n+α)!(n+β)!⋅∑s=0n[1s!(n+α−s)!(β+s)!(n−s)!(x−12)n−s(x+12)s]
MATLAB przekłada się (Jacobim
P olice d epartmentP olynomial ' D ouble' Wykonanie)Rzeczywisty wielomian promieniowy Zernike jest zatem (dla
m=abs(m)
)Uwaga: ta odpowiedź jest tylko praktycznym rozwiązaniem; dodaj inną odpowiedź, która wyjaśnia, dlaczego to działa.
źródło