Численная устойчивость полиномов Зернике высших порядков

9

Я пытаюсь вычислить моменты Zernike более высокого порядка (например m=0, n=46) для некоторого изображения. Тем не менее, я сталкиваюсь с проблемой, касающейся радиального полинома (см. Википедию ). Это полином, определенный на интервале [0 1]. Смотрите код MATLAB ниже

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

Тем не менее, это, очевидно, сталкивается с численными проблемами рядом RHO > 0.9. Многочлен с большим количеством шума

Я попытался реорганизовать его, polyvalдумая, что у него могут быть какие-то лучшие закулисные алгоритмы, но это ничего не решило. Преобразование его в символьное вычисление создало желаемый граф, но было ошеломительно медленным даже для простого графика, такого как показано.

Существует ли численно устойчивый способ оценки таких многочленов высокого порядка?

Sanchises
источник
3
Часто лучше использовать ортогональные многочлены, здесь многочлены Якоби . Вы пробовали mathworks.com/help/symbolic/jacobip.html и его связь
рNм(р)знак равно(-1)(N-м)/2рмп(N-м)/2(м,0)(1-2р2)?
Gammatester
@gammatester Это работает! Не могли бы вы уточнить, почему это так?
Санчиз
Приятно слышать, что это работает. К сожалению, я не могу дать точный ответ по двум причинам. Во-первых: хотя общеизвестно, что ортогональные многочлены обладают лучшими свойствами устойчивости, чем стандартная форма, я не знаю формального доказательства (особенно в этом случае). Во-вторых, я не использую Matlab и не могу дать данные для реализованных полиномов Якоби.
Gammatester
1
@Sanchises. Здесь нет бесплатного обеда: просто потому, что что-то является полиномом, не означает, что прямая формула в терминах степеней является правильным способом его вычисления, и точное вычисление полиномов Якоби само по себе не является тривиальным вопросом - вы не делаете это через коэффициенты, так что это не так дешево.
Кирилл
2
Причиной использования полиномов Якоби является то, что вы избавляетесь от катастрофической отмены в своей формуле (посмотрите на все эти осциллирующие факторы с очень большими коэффициентами!), А процедура полиномиальной оценки Якоби по умолчанию тщательно реализована в библиотеке, так что это гарантировано. быть точным Большая часть работы здесь сделана, чтобы гарантировать, что многочлены Якоби оценены точно.
Кирилл

Ответы:

7

В этой статье Хонарвар и Парамесран выводят интересный метод вычисления радиальных полиномов Цернике очень хорошим рекурсивным способом. Формула рекурсии удивительно проста, без деления или умножения на большие целые числа:

рNм(ρ)знак равноρ(рN-1|м-1|(ρ)+рN-1м+1(ρ))-рN-2м(ρ)
Я бы порекомендовал взглянуть на рисунок 1 в статье Хонарвара и Парамесрана, которая ясно иллюстрирует зависимости между различными полиномами Цернике.

Это реализовано в следующем скрипте Octave:

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)

Например, рисунок, созданный этим кодом, показывает, что с мзнак равно22, а также Nзнак равно112, катастрофическая отмена происходит вблизи ρзнак равно0.7, если радиальные полиномы Зернике вычисляются с помощью полиномов Якоби. Следовательно, нужно также беспокоиться о точности полиномов Цернике нижней степени.

введите описание изображения здесь

Рекурсивный метод, по-видимому, гораздо лучше подходит для стабильного вычисления этих полиномов Цернике высшего порядка. Тем не менее, длямзнак равно0 а также Nзнак равно46максимальная разница между Якоби и рекурсивным методом составляет (только?) 1.4e-10, которая может быть достаточно точной для вашего приложения.

Wim
источник
Ваш сюжет выглядит как ошибка в Matlab'е jacobiPD, а не как какая-либо общая катастрофическая отмена.
Кирилл
@Kiril: я использовал Sanchises ' JacobiPDиз его ответа . Это хорошо работает для полиномов низкого порядка. Например, сNзнак равно30произвольно ми произвольно ρразница между этими двумя методами меньше, чем 6.9e-13. Хотя отдельные суммы в сумме JacobiPDявляются маленькими, они могут стать большими после умножения на factorial(n+a) * factorial(n+b). Кроме того, у них есть знакопеременный знак, который является идеальным рецептом для катастрофической отмены.
Вим
(продолжение) Например, с мзнак равно22 а также Nзнак равно112выражение 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)может стать таким же большим 1.4e18, как и сумма, но только в -2.1конце концов. Вы можете назвать это ошибкой, но с бесконечной точностью ответ будет правильным. Можете ли вы объяснить, что вы имеете в виду под «отсутствием общей катастрофической отмены»?
Вим
1
@ Я не заметил, что это не Мэтлэб. Если чья-то полиномиальная реализация Якоби достаточно хороша для своих целей, это не проблема. Я только сказал, что это ошибка, потому что я неправильно понял и подумал, что это встроенная функция (я ожидаю, что библиотечные функции будут очень надежными). Под «универсальным» я имел в виду, что если вы не знаете, как реализована функция, вы не можете назвать неверные выходные данные «катастрофической отменой», как универсальный термин для всех видов ошибок, но это было только мое недопонимание того, что код делал.
Кирилл
1
Чтобы было ясно: мой код не является рекурсивным. Это итеративное стандартное трехчленное рекуррентное отношение (аналогичное полиномам Чебышева), которое обычно считается более устойчивым, чем, например, форма Хорнера для полиномов.
Gammatester
8

Возможное решение (предложенное @gammatester) заключается в использовании полиномов Якоби. Это позволяет обойти проблему катастрофической отмены при сложении больших полиномиальных коэффициентов путем «наивной» полиномиальной оценки.

Радиальный многочлен Зернике может быть выражен многочленами Якоби следующим образом (см. Уравнение (6) )

рNм(ρ)знак равно(-1)(N-м)/2ρмп(N-м)/2(м,0)(1-2ρ2)

В MATLAB, однако, использование jacobiP(n,a,b,x)неприемлемо медленно для больших векторов / матриц x=rho. Эта jacobiPфункция на самом деле является частью Symbolic Toolbox, и оценка полинома откладывается на символический механизм, который обменивает скорость на произвольную точность. Таким образом, необходима ручная реализация полиномов Якоби.

Поскольку все параметры функции Якоби неотрицательны (αзнак равном, βзнак равно0, N*знак равно(N-м/2)), мы можем использовать следующее выражение (см. Википедия , обратите внимание, что я заполнил значения дляs)

пN(α,β)(ρ)знак равно(N+α)!(N+β)!Σsзнак равно0N[1s!(N+α-s)!(β+s)!(N-s)!(Икс-12)N-s(Икс+12)s]

В MATLAB, это приводит к (Jacobi р olice д epartment Р olynomial, ' D ouble' реализации)

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

Таким образом, фактический радиальный полином Зернике является (для m=abs(m))

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

Примечание: этот ответ на себя - только практическое решение; не стесняйтесь отмечать другой ответ, который объясняет, почему это работает.

Sanchises
источник