Фон

Любая точка реального мира $\mathit{P}\mathit{(}X\mathit{,}\mathit{Y},$Z) может быть определена относительно некоторого 3-D мирового происхождения.

Относительно объектива камеры эта точка 3-D может быть определена как ${\mathit{}}_{\mathrm{p0}}$, которая получается вращением и переводом $\mathit{P}$.

$${}_{\mathrm{p0}}={}_{\mathrm{(}}{\mathrm{x0}}_{\mathrm{,}}{}_{\mathrm{y0}},$$ z0) $$=$$RP + t

Точка 3-D ${\mathit{}}_{\mathrm{p0}}$ затем проецируется в плоскость изображения камеры как точка 2D (${\mathit{}}_{\mathrm{x1}}$, ${\mathit{}}_{\mathrm{}}$y1).

$${}_{\mathrm{}}\frac{{}_{\mathrm{}}}{{}_{\mathrm{x1=x0z0}}}$$, $${}_{\mathrm{}}\frac{{}_{\mathrm{}}}{{}_{\mathrm{y1=y0z0}}}$$

Когда камера захватывает изображение, она точно не захватывает реальные точки, а, скорее, слегка искаженную версию реальных точек, которые можно обозначить (${\mathit{}}_{\mathrm{x2}}$, ${\mathit{}}_{\mathrm{}}$y2). Искаженные точки могут быть описаны с помощью следующей функции:

где:

${\mathit{}}_{\mathrm{k1}}$, ${\mathit{}}_{\mathrm{k2}}$ = коэффициенты радиального искажения объектива

${\mathit{}}_{\mathrm{p1}}$, ${\mathit{}}_{\mathrm{p2}}$ = коэффициенты тангенциального искажения линзы

Искажение

Пример искажения линзы показан ниже; исходное искаженное изображение (слева) и неискаженное изображение (справа).

Обратите внимание на кривизну линий по направлению к краям первого изображения. Для таких приложений, как восстановление и отслеживание изображений, важно знать реальное местоположение точек. Когда мы имеем искаженное изображение, мы знаем искаженные местоположения пикселей (${\mathit{}}_{\mathrm{x2}}$, ${\mathit{}}_{\mathrm{}}$y2). Наша цель - определить неискаженные местоположения пикселей ${\mathit{(}}_{\mathrm{}}$x1${\mathit{,}}_{\mathrm{}}$y1), заданные ${\mathit{}}_{\mathrm{(}}$${\mathit{x2}}_{\mathrm{,}}$y2), и коэффициенты искажения конкретной линзы.

Хотя в остальном это просто, нелинейный характер искажения линзы делает проблему сложной.

Определение модели искажения

Мы начинаем с определения нашей модели искажения:

% Parameters
syms k_1 k_2 p_1 p_2 real
syms r x y
distortionX = subs(x * (1 + k_1 * r^2 + k_2 * r^4) + 2 * p_1 * x * y + p_2 * (r^2 + 2 * x^2), r, sqrt(x^2 + y^2))

distortionX = $$p_2\hspace{0.17em}\left(3\hspace{0.17em}{x}^2+y^2\right)+x\hspace{0.17em}\left(k_2\hspace{0.17em}{\left(x^2+y^2\right)}^2+k_1\hspace{0.17em}\left(x^2+y^2\right)+1\right)+2\hspace{0.17em}{p}_1\hspace{0.17em}x\hspace{0.17em}y$$

distortionY = subs(y * (1 + k_1 * r^2 + k_2 * r^4) + 2 * p_2 * x * y + p_1 * (r^2 + 2 * y^2), r, sqrt(x^2 + y^2))

distortionY = $$p_1\hspace{0.17em}\left(x^2+3\hspace{0.17em}{y}^2\right)+y\hspace{0.17em}\left(k_2\hspace{0.17em}{\left(x^2+y^2\right)}^2+k_1\hspace{0.17em}\left(x^2+y^2\right)+1\right)+2\hspace{0.17em}{p}_2\hspace{0.17em}x\hspace{0.17em}y$$

Радиальное искажение ${\mathit{}}_{\mathrm{k1}}\mathrm{=}$0

Мы строим сетку местоположений пикселей, предполагая, что наша линза имеет коэффициент радиального искажения ${\mathit{}}_{\mathrm{k1}}\mathrm{=}$0. Заметим, что искажение является наименьшим вблизи центра изображения и наибольшим вблизи краев.

% Set Parameters
parameters = [k_1 k_2 p_1 p_2];
parameterValues = [0 0 0 0];
plotLensDistortion(distortionX,distortionY,parameters,parameterValues)

spacing = 0.2000

distortionX = $$x$$

distortionY = $$y$$

Радиальное искажение ${\mathit{}}_{\mathrm{k1}}\mathrm{=}\mathrm{}$0,15

Изучите чувствительность к изменениям в ${\mathit{}}_{\mathrm{k1}}$.

% Set Parameters
parameters = [k_1 k_2 p_1 p_2];
parameterValues = [0.15 0 0 0];
plotLensDistortion(distortionX,distortionY,parameters,parameterValues)

spacing = 0.2000

distortionX = 
$$x\hspace{0.17em}\left(\frac{3\hspace{0.17em}{x}^2}{20}+\frac{3\hspace{0.17em}{y}^2}{20}+1\right)$$

distortionY = 
$$y\hspace{0.17em}\left(\frac{3\hspace{0.17em}{x}^2}{20}+\frac{3\hspace{0.17em}{y}^2}{20}+1\right)$$

Расчет модели обратного искажения

Учитывая коэффициенты искажения объектива камеры и набор искаженных местоположений пикселей (${\mathit{}}_{\mathrm{x2}}$, ${\mathit{}}_{\mathrm{}}$y2), мы хотим иметь возможность вычислить неискаженные местоположения пикселей ${\mathit{(}}_{\mathrm{}}$x1${\mathit{,}}_{\mathrm{}}$y1). Рассмотрим конкретный случай, когда все коэффициенты искажения равны нулю, за исключением ${\mathit{}}_{\mathrm{}}$k1, который равен 0,2.

Мы начинаем с определения коэффициентов искажения

syms X Y positive
eq1 = X == distortionX

eq1 = $$X=p_2\hspace{0.17em}\left(3\hspace{0.17em}{x}^2+y^2\right)+x\hspace{0.17em}\left(k_2\hspace{0.17em}{\left(x^2+y^2\right)}^2+k_1\hspace{0.17em}\left(x^2+y^2\right)+1\right)+2\hspace{0.17em}{p}_1\hspace{0.17em}x\hspace{0.17em}y$$

eq2 = Y == distortionY

eq2 = $$Y=p_1\hspace{0.17em}\left(x^2+3\hspace{0.17em}{y}^2\right)+y\hspace{0.17em}\left(k_2\hspace{0.17em}{\left(x^2+y^2\right)}^2+k_1\hspace{0.17em}\left(x^2+y^2\right)+1\right)+2\hspace{0.17em}{p}_2\hspace{0.17em}x\hspace{0.17em}y$$

Мы определяем уравнения искажений для заданных коэффициентов искажений и решаем для неискаженных местоположений пикселей (${\mathit{}}_{\mathrm{x1}}$, ${\mathit{}}_{\mathrm{}}$y1).

parameters = [k_1 k_2 p_1 p_2];
parameterValues = [0.2 0 0 0];
eq1 = expand(subs(eq1, parameters, parameterValues))

eq1 = 
$$X=\frac{x^3}{5}+\frac{x\hspace{0.17em}{y}^2}{5}+x$$

eq2 = expand(subs(eq2, parameters, parameterValues))

eq2 = 
$$Y=\frac{x^2\hspace{0.17em}y}{5}+\frac{y^3}{5}+y$$

Result = solve([eq1, eq2], [x,y], 'MaxDegree', 3,'Real',true)

Result = struct with fields:
    x: [1x1 sym]
    y: [1x1 sym]

Поскольку элемент 1 является единственным реальным решением, мы извлекаем это выражение в его собственную переменную.

[Result.x Result.y]

ans = 
$$\begin{array}{l}\left(\begin{array}{cc}\frac{X\hspace{0.17em}{\sigma }_1}{Y}& {\sigma }_1\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1={\sigma }_2-\frac{5\hspace{0.17em}{Y}^2}{3\hspace{0.17em}\left(X^2+Y^2\right)\hspace{0.17em}{\sigma }_2}\\ \\ \mathrm{ }{\sigma }_2={\left(\frac{5\hspace{0.17em}{Y}^3}{2\hspace{0.17em}\left(X^2+Y^2\right)}+\sqrt{\frac{25\hspace{0.17em}{Y}^6}{4\hspace{0.17em}{\left(X^2+Y^2\right)}^2}+\frac{125\hspace{0.17em}{Y}^6}{27\hspace{0.17em}{\left(X^2+Y^2\right)}^3}}\right)}^{1/3}\end{array}$$

Теперь у нас есть аналитические выражения для местоположений пикселей X и Y, которые мы можем использовать для отмены искажения наших изображений.

Функция рисования искажения объектива

function plotLensDistortion(distortionX,distortionY,parameters,parameterValues)
% distortionX is the expression describing the distorted x coordinate
% distortionY is the expression describing the distorted y coordinate
% k1 and k2 are the radial distortion coefficients 
% p1 and p2 are the tangential distortion coefficients 

syms x y

% This is the grid spacing over the image
spacing = 0.2

% Inspect and parametrically substitute in the values for k_1 k_2 p_1 p_2
distortionX = subs(distortionX,parameters,parameterValues)
distortionY = subs(distortionY,parameters,parameterValues)

% Loop over the grid
for x_i = -1:spacing:1
    for y_j = -1:spacing:1
        
        % Compute the distorted location 
        xout = subs(distortionX, {x,y}, {x_i,y_j});
        yout = subs(distortionY, {x,y}, {x_i,y_j});
        
        % Plot the original point
        plot(x_i,y_j, 'o', 'Color', [1.0, 0.0, 0.0])
        hold on
       
        % Draw the distortion direction with Quiver
        p1 = [x_i,y_j];                     % First Point
        p2 = [xout,yout];                   % Second Point
        dp = p2-p1;                         % Difference
        quiver(p1(1),p1(2),dp(1),dp(2),'AutoScale','off','MaxHeadSize',1,'Color',[0 0 1])
        
    end
end
hold off
grid on

end