Фон

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

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

$$p_0=(x_0,y_0,z_0)$$ = $$RP+t$$

Точка 3-D ${\mathit{p}}_0$ затем проецируется в плоскость изображения камеры как точка 2D, (${\mathit{x}}_1$,${\mathit{y}}_1$).

$$x_1=\frac{x_0}{z_0}$$ , $$y_1=\frac{y_0}{z_0}$$

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

где:

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

${\mathit{p}}_1$, ${\mathit{p}}_2$ = тангенциальные коэффициенты искажения объектива

Искажение

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

Обратите внимание на кривизну линий к ребрам первого изображения. Для таких приложений, как реконструкция изображений и отслеживание, важно знать реальное расположение точек в мире. Когда у нас есть искаженное изображение, мы знаем искаженные пиксельные местоположения (${\mathit{x}}_2$,${\mathit{y}}_2$). Наша цель - определить неискаженные пиксельные местоположения (${\mathit{x}}_1$,${\mathit{y}}_1$) данный (${\mathit{x}}_2$,${\mathit{y}}_2$) и коэффициенты искажения конкретной линзы.

В противном случае нелинейная природа искажения объектива делает задачу сложной.

Задайте модель искажения

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

% 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{k}}_1=0$

Мы строим сетку пиксельных местоположений, принимая, что наша линза имеет радиальный коэффициент искажения ${\mathit{k}}_1=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{k}}_1=0.15$

Исследуйте чувствительность к изменениям в ${\mathit{k}}_1$.

% 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{x}}_2$,${\mathit{y}}_2$), мы хотим иметь возможность вычислить неискаженные пиксельные местоположения (${\mathit{x}}_1$,${\mathit{y}}_1$). Мы рассмотрим конкретный случай, когда все коэффициенты искажения равны нулю, кроме ${\mathit{k}}_1$ который равен 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{x}}_1$,${\mathit{y}}_1$).

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