Фон

Любая точка реального мира $\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: (X*(((5*Y^3)/(2*(X^2 + Y^2)) + ((25*Y^6)/(4*(X^2 + Y^2)^2) + (125*...
    y: ((5*Y^3)/(2*(X^2 + Y^2)) + ((25*Y^6)/(4*(X^2 + Y^2)^2) + (125*Y^6)...

Поскольку элементом 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