Красиво описать дисторсию объектива

[ilghiz]

Добрый день,

имею очень китайский объектив и хочу затабулировать для него дисторсию, причем дисторсия не центрально симметрична, и имеет реально недецкие искажения. Автофокуса у объектива нет.

Я хочу отписать дисторсию массивом размера с число пикселей так, что $(n,m)$ пиксель, который видит камера, в реалии на плоскости матрицы сенсоров попадает на позицию $(p_{n,m}^{(x)}, p_{n,m}^{(y)}), n=1,...,N, m=1,...,M.$

Итак, у меня есть градуировочный функционал (набор картинок, по которым я градуирую дисторсию), который можно записать так

$f_i(s_j) \in R^2,$

$f_i(s_j) \simeq (p_{\alpha_{ij},\beta_{ij}}^{(x)}, p_{\alpha_{ij},\beta_{ij}}^{(y)})^T$

здесь

индекс $i=1,...,K_i$ бежит по реперным точкам,

индекс $j=1,...,K_j$ бежит по номеру картинки,

$s_j$ - параметры поворота картинки относительно камеры (я не могу гарантировать, что я всегда и супер точно располагаю картинку относительно камеры), реально там унитарная матрица поворотов и вектор линейного переноса, то есть 6 параметров,

$\alpha_{ij},\beta_{ij}$ то, куда в пикселах попала $i$-ая реперная точка для $j$-ой картинки,

$f_i(s) \in R^2$ - функция отображения $i$-ой точки для заданной ($s$) трансформации,

Дополнительно я хочу потребовать, что массивы $p_{n,m}^{(x)}$ и $p_{n,m}^{(y)}$ имеют по $x$ и $y$ очень маленькие отклонения от единицы, то есть я формулирую задачу минимизации функции вида

$$\min_{s_1,...,s_{K_j}, p_{1,1},...,p_{n,m}} \sum_{i=1,...,K_i} \sum_{j=1,...,K_j} ||f_i(s_j) - (p_{\alpha_{ij},\beta_{ij}}^{(x)}, p_{\alpha_{ij},\beta_{ij}}^{(y)})^T ||_2^2$$

$$ + \gamma \left[ \sum_{n=2,...,N; m=1,...,M} | p_{n,m}^{(x)} - p_{n-1,m}^{(x)} - 1|_2^2 + \sum_{n=2,...,N; m=1,...,M} | p_{n,m}^{(y)} - p_{n-1,m}^{(y)} - 1|_2^2 \right.$$

$$ \left. + \sum_{n=1,...,N; m=2,...,M} | p_{n,m}^{(x)} - p_{n,m-1}^{(x)} - 1|_2^2 + \sum_{n=1,...,N; m=2,...,M} | p_{n,m}^{(y)} - p_{n,m-1}^{(y)} - 1|_2^2 \right] $$

Не найдя ничего лучшего я вижу возможность решать эту задачу покомпонентной минимизацией, то есть на четном шаге мы фиксируем значения $s_1,...,s_K$ и ищем значения $p_{1,1},...,p_{n,m}$, а не нечетном, наоборот.

Задача минимизации $s_1,...,s_{K_j}$, ИМХО, довольно легко решается, так как она разваливается на $K_j$ независимых задач с 6 неизвестными, для которых можно посчитать аналитически (хоть и не сильно просто) и градиент, и Гессиан, тогда Ньютон в окрестности решения должен быстро сходиться. А вот задача по $p_{1,1},...,p_{n,m}$ меня немного озадачила.

Тут все линейно, все сводится к двум системам линейных уравнений с разреженной матрицей вида

$$(\gamma H_{N \times N} \times I_{M \times M} + \gamma I_{N \times N} \times H_{M \times M} + D) p = r + u_N \times z_M + z_N \times u_M,$$

$$H = \left( \begin{tabular}{ccccccc} 1 & -1 & 0 & 0 & ... & ... & 0 \\ -1 & 2 & -1 & 0 & ... & ... & 0 \\ 0 & -1 & 2 & -1 & ... & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & ... & ... & ... & ... & 0 \\ 0 & ... & ... & ... & -1 & 2 & -1 \\ 0 & ... & ... & ... & 0 & -1 & 1 \end{tabular} \right) \in \R^{N \times N},$$

$I_{M \times M} \in \R^{M \times M}$ - единичная матрица,

$D \in \R^{NM \times NM}$ - диагональная матрица, где на диагонали стоят единички в позициях, которые определены индексами $\alpha_{ij},\beta_{ij}$,

$u_N = (1,...,1)^T \in \R^N$,

$z_N = (-1, 0, ..., 0, 1) \in \R^N$,

$r \in \R^{NM}$ вектор, в котором ненулевые значения определены $\alpha_{ij},\beta_{ij}$.

Я понимаю как решить систему или перестановкой и полным $LDL^T$, или сопряженным градиентом с или без предобуславливателя, но мне не нравится ее размерность.

$K_i$ и $K_j$ у меня маленькие, около сотни, и все сопутствующие параметры $s_1,...,s_{K_j}, p_{1,1},...,p_{n,m}$ поместятся в память микроконтроллера, объем которого около 500 КБайт, а вот размерность матрицы (даже не матрицы) а просто решения $p$ у меня будет около 2-10 миллионов, причем будет две такие векторные системы, и, даже если рассмотреть итерационное решение методом сопряженного градиента без предобуславливателя, где умножение на матрицу будет происходить совсем без генерации самой матрицы, мне потребуется не менее 5 векторов по 10 миллионов, и я совсем не впишусь в память. У процессора конечно есть внешняя память около ГБайта, но это флеш и в ней не поитерируешь, но посчитанную таблицу хранить можно.

Понятно, можно задать $p_{1,1},...,p_{n,m}$ на сетке, в которой, например, задана только каждая 10-ая точка, и аппроксимировать все, что между ними либо линейно, либо кубическими сплайнами, размерность задачи, конечно, будет существенно меньше, и есть вероятность, что взяв билинейную интерполяцию на каждой 24-той точке на 10МПикселей мне потребуется только 500КБайт, но есть сомнения, что точности хватит, да и по памяти совсем впритык. Получаемые формулы будут очень похожими, но довольно громоздкими, поэтому я их тут не привожу.

Вопрос.

Скажите, пожалуйста, как еще можно регуляризовать эту задачу, чтобы натянуть сову на кактус, то есть чтобы выполнить условие, что две соседние точки в $p_{1,1},...,p_{n,m}$ отстоят друг от друга почти на единицу, но не генерировать явно весь тот массив точек?

Спасибо!

PS: Чтобы подтвердить мои слова, почему я выбрал такую сложную модель для дисторсии, прилагаю фото этой камеры листа миллиметровки. Дисторсия видна совершенно не вооруженным глазом.