![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
eddy_em«Как она ни плакала…», а все равно нужно реализовать разложение волнового фронта (а затем — и градиентного поля для восстановления ВФ) по ортонормированным на кольце полиномам, линейно зависимым от полиномов Цернике (чтобы можно было восстановить по найденным коэффициентам разложения коэффициенты Цернике).
Пока — размышления.
Классикой в этом деле является статья Mahajan V.N. «Zernike annular polynomials for imaging systems with annular pupils», вот ее данные:
Mahajan предлагает довольно-таки громоздкую схему ортогонализации (по Грэму-Шмидту), достоинство которой в том, что вычисляются лишь нужные N полиномов, а недостаток в том, что делается это рекурсивно.
У Mahajan выводятся полиномы общего вида (к тому же непрерывные, а не дискретные). Учитывая то, что физического смысла в этих полиномах нет, можно попробовать реализовать какой-нибудь другой вариант.
Итак, мы имеем значения волнового фронта в M точках (вектор-столбец), W, и значения N полиномов Цернике (в нотации Нолля) в данных точках (матрица M строк на N столбцов), Z. Нам надо найти N коэффициентов Цернике (вектор-столбец), z, для однозначной аппроксимации данного волнового фронта, т.е. в идеале должно получиться
Z · z = W.
Так как Z не ортогональны на кольце, необходимо ввести промежуточный класс полиномов, U, по которым можно будет разложить волновой фронт: U·u=W, а сами U должны быть линейно зависимыми от Z (чтобы можно было однозначно восстановить z, зная b).
Одним словом, матричные методы базируются на следующих уравнениях:
U = Z · k, UT ≡ U-1 → u = WTU = UTW,
а так как U k-1 = Z, получим, что коэффициенты u очень просто преобразуются в z: z = u k-1. Следовательно, задача — правильно выбрать матрицу k, чтобы получить ортогональную матрицу U.
A = Q · R.
Матрица A может быть прямоугольной, N столбцов (количество полиномов) и M строк (количество точек с данными), в этом случае ортогональным базисом будут первые N столбцов квадратной матрицы Q.
Каждый столбец матрицы R будет представлять собой коэффициенты для выражения ортогональных полиномов через неортогональные. Восстановить затем коэффициенты оригинальных полиномов можно по строкам матрицы R: сумма элементов каждой строки и есть искомые коэффициенты.
В общем, в данном случае U ≡ Q и k ≡ R-1.
Если мы работаем с набором в несколько тысяч точек, этот метод будет пожирать очень много ресурсов! Однако, с точки зрения формальной реализации он, пожалуй, наипростейший.
AT·A = L·LT.
Ортогональный базис получаем, умножая оригинал на (L-1)T:
U ≡ A (L-1)T,
отсюда получаем: k ≡ (L-1)T.
Судя по тому, что во всех найденных мной статьях авторы ссылаются на оригинальную статью Mahajan, другие варианты если и есть, то «сокрыты они за семью печатями».
Кстати, я даже нагуглил готовый код для вычисления ортогональных на кольце полиномов Цернике по Mahajan, но, к сожалению, они написаны под Wolfram mathematica, т.е. толку от него не больше, чем от формул у Махаджи (хотя, нет: у Махаджи понятней).
Пока — размышления.
Классикой в этом деле является статья Mahajan V.N. «Zernike annular polynomials for imaging systems with annular pupils», вот ее данные:
@ARTICLE{1981JOSA...71...75M,
author = {{Mahajan}, V.~N.},
title = "{Zernike annular polynomials for imaging systems with annular pupils.}",
journal = {Journal of the Optical Society of America (1917-1983)},
keywords = {Astronomical Optics},
year = 1981,
volume = 71,
pages = {75-85},
adsurl = {http://adsabs.harvard.edu/abs/1981JOSA...71...75M},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
Mahajan предлагает довольно-таки громоздкую схему ортогонализации (по Грэму-Шмидту), достоинство которой в том, что вычисляются лишь нужные N полиномов, а недостаток в том, что делается это рекурсивно.
У Mahajan выводятся полиномы общего вида (к тому же непрерывные, а не дискретные). Учитывая то, что физического смысла в этих полиномах нет, можно попробовать реализовать какой-нибудь другой вариант.
Матричные методы
Эти методы в основном являются разными видами ортогонализации. К счастью, ортогонализация — очень популярный в алгебре процесс, поэтому в GSL все необходимое есть, т.е. одним велосипедом получается меньше.Итак, мы имеем значения волнового фронта в M точках (вектор-столбец), W, и значения N полиномов Цернике (в нотации Нолля) в данных точках (матрица M строк на N столбцов), Z. Нам надо найти N коэффициентов Цернике (вектор-столбец), z, для однозначной аппроксимации данного волнового фронта, т.е. в идеале должно получиться
Z · z = W.
Так как Z не ортогональны на кольце, необходимо ввести промежуточный класс полиномов, U, по которым можно будет разложить волновой фронт: U·u=W, а сами U должны быть линейно зависимыми от Z (чтобы можно было однозначно восстановить z, зная b).
Одним словом, матричные методы базируются на следующих уравнениях:
U = Z · k, UT ≡ U-1 → u = WTU = UTW,
а так как U k-1 = Z, получим, что коэффициенты u очень просто преобразуются в z: z = u k-1. Следовательно, задача — правильно выбрать матрицу k, чтобы получить ортогональную матрицу U.
QR-декомпозиция
Неортогональная матрица A (в нашем случае — полиномы Цернике для кольца) может быть разложена на ортогональную Q и треугольную (верхнюю) R:A = Q · R.
Матрица A может быть прямоугольной, N столбцов (количество полиномов) и M строк (количество точек с данными), в этом случае ортогональным базисом будут первые N столбцов квадратной матрицы Q.
Каждый столбец матрицы R будет представлять собой коэффициенты для выражения ортогональных полиномов через неортогональные. Восстановить затем коэффициенты оригинальных полиномов можно по строкам матрицы R: сумма элементов каждой строки и есть искомые коэффициенты.
В общем, в данном случае U ≡ Q и k ≡ R-1.
Если мы работаем с набором в несколько тысяч точек, этот метод будет пожирать очень много ресурсов! Однако, с точки зрения формальной реализации он, пожалуй, наипростейший.
Разложение Холецкого
(если бы не глянул на ссылку на русскую статью в википедии, так и думал бы дальше, что фамилия — Чолески, а не Холецкий!)
В данном случае используется тот факт, что ATA можно представить как произведение левой диагональной матрицы на ее сопряжение:AT·A = L·LT.
Ортогональный базис получаем, умножая оригинал на (L-1)T:
U ≡ A (L-1)T,
отсюда получаем: k ≡ (L-1)T.
Метод Mahajan
Конечно, в случае неудачи матричных способов можно будет попробовать и оригинальный метод Mahajan. К сожалению, он имеет зависимость от ε, поэтому затабулировать его заранее не получится, а рекурсия не даст полноценно распараллелить вычисления.Судя по тому, что во всех найденных мной статьях авторы ссылаются на оригинальную статью Mahajan, другие варианты если и есть, то «сокрыты они за семью печатями».
Кстати, я даже нагуглил готовый код для вычисления ортогональных на кольце полиномов Цернике по Mahajan, но, к сожалению, они написаны под Wolfram mathematica, т.е. толку от него не больше, чем от формул у Махаджи (хотя, нет: у Махаджи понятней).
