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

Этот пример разрабатывает математическую модель с помощью Symbolic Math Toolbox, чтобы не исказить изображение и показывает локальную функцию в live скрипте.

Фон

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

Относительно объектива фотокамеры эта 3-D точка может быть задана как p0, который получен путем вращения и перевода P.

p0=(x0,y0,z0) = RP+t

3-D точка p0 затем спроектирован в плоскость изображения камеры как 2D точка, (x1,y1).

x1=x0z0 , y1=y0z0

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

x2=x1(1+k1r2+k2r4)+2p1x1y1+p2(r2+2x12)

y2=y1(1+k1r2+k2r4)+2p2x1y1+p1(r2+2y12)

где:

k1, k2 = радиальные коэффициенты искажения линзы

p1, p2 = тангенциальные коэффициенты искажения линзы

r=x12+y12

Искажение

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

Отметьте искривление строк к ребрам первого изображения. Для приложений, таких как реконструкция изображений и отслеживание, важно знать местоположение реального мира точек. Когда у нас есть искаженное изображение, мы знаем искаженные пиксельные местоположения (x2,y2). Это - наша цель определить неискаженные пиксельные местоположения (x1,y1) данный (x2,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 = p23x2+y2+xk2x2+y22+k1x2+y2+1+2p1xy
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 = p1x2+3y2+yk2x2+y22+k1x2+y2+1+2p2xy

Радиальное искажение k1=0

Мы строим сетку пиксельных местоположений, принимающих, что наша линза имеет радиальный коэффициент искажения k1=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

Радиальное искажение k1=0.15

Исследуйте чувствительность к изменениям в 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 = 

x3x220+3y220+1

distortionY = 

y3x220+3y220+1

Вычислите обратную модель искажения

Учитывая коэффициенты искажения объектива камеры и набор искаженных пиксельных местоположений (x2,y2), мы хотим смочь вычислить неискаженные пиксельные местоположения (x1,y1). Мы посмотрим на конкретный случай, где все коэффициенты искажения являются нулем за исключением k1 который равняется 0.2.

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

syms X Y positive
eq1 = X == distortionX
eq1 = X=p23x2+y2+xk2x2+y22+k1x2+y2+1+2p1xy
eq2 = Y == distortionY
eq2 = Y=p1x2+3y2+yk2x2+y22+k1x2+y2+1+2p2xy

Мы определяем уравнения искажения для данных коэффициентов искажения и решаем для неискаженных пиксельных местоположений (x1,y1).

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

X=x35+xy25+x

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

Y=x2y5+y35+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 = 

(Xσ1Yσ1)где  σ1=σ2-5Y23X2+Y2σ2  σ2=5Y32X2+Y2+25Y64X2+Y22+125Y627X2+Y231/3

Теперь у нас есть аналитические выражения для пиксельных 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
Для просмотра документации необходимо авторизоваться на сайте