9
$\begingroup$

I'm trying to calculate higher order (e.g., m=0, n=46) Zernike moments for some image. However, I'm running into a problem regarding the radial polynomial (see wikipedia). This is a polynomial defined on the interval [0 1]. See the MATLAB code below

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

However, this obviously runs into numerical issues near RHO > 0.9. A polynomial with lots of noise

I tried refactoring it to polyval thinking it might have some better behind-the-scenes algorithms but that didn't didn't solve anything. Converting it to a symbolic calculation did create the desired graph but was mindbogglingly slow even for a simple graph such as shown.

Is there a numerically stable way of evaluating such high-order polynomials?

Improve this question
$\endgroup$
  • 3
    $\begingroup$ Often it is better to use orthogonal polynomials, here Jacobi polynomials. Have you tried mathworks.com/help/symbolic/jacobip.html and the relation $$ R_n^{\,m}(r) = (-1)^{(n-m)/2}\,r^m\,P_{(n-m)/2}^{\,(m,0)}(1-2r^2)? $$ $\endgroup$ – gammatester Dec 13 '18 at 11:27
  • $\begingroup$ @gammatester That works! Could you perhaps elaborate in an answer on why this would be the case? $\endgroup$ – Sanchises Dec 13 '18 at 12:37
  • $\begingroup$ Nice hear that it works. Unfortunately I cannot give a decicated answer for two reasons. First: although it is commonly known that orthogonal polynomials have better stability properties than the standard form, I do not know a formal proof (especially in this case). Second I do not use Matlab and cannot give data for the implemented Jacobi polynomials. $\endgroup$ – gammatester Dec 13 '18 at 14:09
  • 1
    $\begingroup$ @Sanchises There's no free lunch here: just because something is a polynomial doesn't mean the direct formula in terms of powers is the right way to compute it, and computing Jacobi polynomials accurately is not itself a trivial matter—you don't do it through the coefficients, so it's not as cheap. $\endgroup$ – Kirill Dec 13 '18 at 23:52
  • 2
    $\begingroup$ The reason it works to use Jacobi polynomials is that you get rid of the catastrophic cancellation in your formula (look at all those oscillating factors with very large coefficients!), and the default Jacobi polynomial evaluation procedure is implemented carefully in a library so is guaranteed to be accurate. Most of the work here is done in making sure the Jacobi polynomials are evaluated accurately. $\endgroup$ – Kirill Dec 13 '18 at 23:56
7
$\begingroup$

In this paper, Honarvar and Paramesran derive an interesting method to compute the radial Zernike polynomials in a very nice recursive way. The recursion formula is surprisingly straightforward, without division or multiplication by large integers: $$ R^m_n(\rho) = \rho \left(R^{|m-1|}_{n-1}(\rho)+R^{m+1}_{n-1}(\rho)\right) - R^{m}_{n-2}(\rho)$$ I'd recommend to have a look at figure 1 in the Honarvar and Paramesran paper, which clearly illustrates the dependencies between the different Zernike polynomials.

This is implemented in the following Octave script:

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)

For example, the figure produced by this code shows that with $m = 22$, and $n = 112$, catastrophic cancellation occurs near $\rho = 0.7$, if Zernike radial polynomials are computed via Jacobi polynomials. Therefore, one also has to worry about the accuracy of the lower-degree Zernike polynomials.

enter image description here

The recursive method seems to be much more suitable for computing these higher-order Zernike polynomials in a stable way. Nevertheless, for $m = 0$ and $n = 46$, the maximum difference between the Jacobi and the recursive method is (only?) 1.4e-10, which might be accurate enough for your application.

Improve this answer
$\endgroup$
  • $\begingroup$ Your plot looks like a bug in Matlab's jacobiPD, not like any generic catastrophic cancellation. $\endgroup$ – Kirill Dec 17 '18 at 18:33
  • $\begingroup$ @Kiril: I used Sanchises' JacobiPD from his answer. This works well for low-order polynomials. For example, with $n=30$, arbitrary $m$, and arbitrary $\rho$, the difference between the two methods is less than 6.9e-13. Although the individual terms in the summation of JacobiPD are small, they may become large after multiplying by factorial(n+a) * factorial(n+b). Moreover they have alternating sign, which is a perfect recipe for catastrophic cancellation. $\endgroup$ – wim Dec 17 '18 at 20:53
  • $\begingroup$ (continued) E.g. with $m=22$ and $n=112$ the expression 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), may become as large as 1.4e18, while the sum is only -2.1 eventually. You can call this a bug, but with infinite precision, the answer would have been correct. Can you explain what you mean by "no generic catastrophic cancellation"? $\endgroup$ – wim Dec 17 '18 at 20:55
  • 1
    $\begingroup$ @wim I didn't notice it's not Matlab's. If somebody's Jacobi polynomial implementation is good enough for their purpose, that's no problem. I only said it's a bug because I misunderstood and thought it's a built-in function (I expect library functions to be very solid). By "generic" I meant that if you don't know how the function is implemented, you can't call incorrect outputs "catastrophic cancellation" like a catch-all term for all kinds of errors, but that was just my misunderstanding of what the code was doing. $\endgroup$ – Kirill Dec 17 '18 at 21:01
  • 1
    $\begingroup$ To be clear: my code is not recursive. It is the iterative standard three term recurrence relation (similiar to Chebyshev polynomials) which is supposed to be normally more stable than e.g. the Horner form for polynomials. $\endgroup$ – gammatester Dec 17 '18 at 21:49
8
$\begingroup$

A possible solution (suggested by @gammatester) is to use Jacobi polynomials. This circumvents the problem of catastrophic cancellation in adding the large polynomial coefficients by 'naive' polynomial evaluation.

The radial Zernike polynomial can be expressed by Jacobi polynomials as follows (see equation (6))

$$ R^m_n(\rho) = (-1)^{(n-m)/2}\rho^m \cdot P^{(m,0)}_{(n-m)/2} \Big(1-2\rho^2 \Big)$$

In MATLAB however, the use of jacobiP(n,a,b,x) is unacceptably slow for large vectors/matrices of x=rho. The jacobiP function is actually part of the Symbolic Toolbox, and evaluation of the polynomial is deferred to the symbolic engine, which trades speed for arbitrary precision. A manual implementation of the Jacobi polynomials is thus necessary.

Since the parameters to the Jacobi function are all nonnegative ($\alpha=m$, $\beta=0$, $n^*=(n-m/2)$), we can use the following expression (see Wikipedia, note that I filled in the values for $s$) \begin{multline} P_n^{(\alpha,\beta)}(\rho) = (n+\alpha)!(n+\beta)! \quad \cdot\\ \sum_{s={0}}^n \left[\frac{1}{s!(n+\alpha-s)!(\beta+s)!(n-s)!} \left(\frac{x-1}{2}\right)^{n-s}\left(\frac{x+1}{2}\right)^{s} \right]\end{multline}

In MATLAB, this translates to (Jacobi police department Polynomial, 'Double' implementation)

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

The actual radial Zernike polynomial is thus (for m=abs(m))

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

Note: this self-answer is only a practical solution; feel free to tag on another answer that explains why this works.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.