イメージの歪み補正アルゴリズムの開発

この例では、Symbolic Math Toolbox™ を使用して、イメージの歪みを補正する数学モデルを開発し、ライブスクリプトでローカル関数の特徴を示します。

問題の背景

現実世界の任意の点 $P (X, Y, Z)$ は、ある 3 次元世界の原点に対して定義できます。

カメラのレンズとの相関において、この 3 次元の点は $P$ の回転と平行移動によって得られる $p_{0}$ として定義できます。

$p_{0} = (x_{0}, y_{0}, z_{0})$ = $R P + t$

次に、3 次元の点 $p_{0}$ をカメラのイメージ平面に 2 次元の点 ( $x_{1}$ , $y_{1}$ ) として投影します。

$x_{1} = \frac{x_{0}}{z_{0}}$ , $y_{1} = \frac{y_{0}}{z_{0}}$

カメラがイメージをキャプチャする場合、実際の点を正確にキャプチャするのではなく、むしろ実際の点が少し歪んだものをキャプチャします。これは ( $x_{2}$ , $y_{2}$ ) で表すことができます。歪んだ点は、以下の関数で記述できます。

$x_{2} = x_{1} (1 + k_{1} r^{2} + k_{2} r^{4}) + 2 p_{1} x_{1} y_{1} + p_{2} (r^{2} + 2 {x_{1}}^{2})$

$y_{2} = y_{1} (1 + k_{1} r^{2} + k_{2} r^{4}) + 2 p_{2} x_{1} y_{1} + p_{1} (r^{2} + 2 {y_{1}}^{2})$

ここで

$k_{1}$ 、 $k_{2}$ = レンズの半径方向歪み係数

$p_{1}$ 、 $p_{2}$ = レンズの円周方向歪み係数

$r = \sqrt{{x_{1}}^{2} + {y_{1}}^{2}}$

歪み

レンズ歪みの例 (以下を参照): 元の歪んだイメージ (左) と歪み補正後のイメージ (右)

最初のイメージのエッジに向かうラインの湾曲に注目してください。イメージの再構成やトラッキングなどのアプリケーションでは、現実世界の点の位置を理解することが重要になります。歪んだイメージがある場合、歪んだピクセルの位置 ( $x_{2}$ , $y_{2}$ ) はわかっています。( $x_{2}$ , $y_{2}$ ) と特定のレンズの歪み係数を与えて、歪み補正後のピクセル位置 ( $x_{1}$ , $y_{1}$ ) を決定することが目標です。

特に複雑な部分はないにしても、レンズ歪みの非線形特性によって難易度の高い問題となっています。

歪みモデルの定義

歪みモデルの定義から開始

% 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} (3 x^{2} + y^{2}) + x (k_{2} {(x^{2} + y^{2})}^{2} + k_{1} (x^{2} + y^{2}) + 1) + 2 p_{1} x 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} (x^{2} + 3 y^{2}) + y (k_{2} {(x^{2} + y^{2})}^{2} + k_{1} (x^{2} + y^{2}) + 1) + 2 p_{2} x y$

半径方向歪み $k_{1} = 0$

レンズの半径方向歪み係数を $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$

半径方向歪み $k_{1} = 0.15$

$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 (\frac{3 x^{2}}{20} + \frac{3 y^{2}}{20} + 1)$

distortionY = 
 $y (\frac{3 x^{2}}{20} + \frac{3 y^{2}}{20} + 1)$

逆歪みモデルの計算

カメラのレンズ歪み係数と一連の歪んだピクセル位置 ( $x_{2}$ , $y_{2}$ ) を与えて、歪み補正後のピクセル位置 ( $x_{1}$ , $y_{1}$ ) を計算できるようにします。0.2 に等しい $k_{1}$ を除き、すべての歪み係数が 0 の特殊なケースを調べます。

まず、歪み係数を定義します。

syms X Y positive
eq1 = X == distortionX

eq1 =  $X = p_{2} (3 x^{2} + y^{2}) + x (k_{2} {(x^{2} + y^{2})}^{2} + k_{1} (x^{2} + y^{2}) + 1) + 2 p_{1} x y$

eq2 = Y == distortionY

eq2 =  $Y = p_{1} (x^{2} + 3 y^{2}) + y (k_{2} {(x^{2} + y^{2})}^{2} + k_{1} (x^{2} + y^{2}) + 1) + 2 p_{2} x y$

与えられた歪み係数に対する歪み方程式を定義して、歪み補正後のピクセル位置 ( $x_{1}$ , $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 y^{2}}{5} + x$

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

eq2 = 
 $Y = \frac{x^{2} 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^6)/(27*(X^2 + Y^2)^3))^(1/2))^(1/3) - (5*Y^2)/(3*(X^2 + Y^2)*((5*Y^3)/(2*X^2 + 2*Y^2) + ((25*Y^6)/(4*(X^2 + Y^2)^2) + (125*Y^6)/(27*(X^2 + Y^2)^3))^(1/2))^(1/3))))/Y
    y: ((5*Y^3)/(2*(X^2 + Y^2)) + ((25*Y^6)/(4*(X^2 + Y^2)^2) + (125*Y^6)/(27*(X^2 + Y^2)^3))^(1/2))^(1/3) - (5*Y^2)/(3*(X^2 + Y^2)*((5*Y^3)/(2*X^2 + 2*Y^2) + ((25*Y^6)/(4*(X^2 + Y^2)^2) + (125*Y^6)/(27*(X^2 + Y^2)^3))^(1/2))^(1/3))

要素 1 は実数解のみであるため、その式を自身の変数へ抽出します。

[Result.x Result.y]

ans = 
 $\begin{array}{l} (\begin{array}{c} \frac{X σ_{1}}{Y} & σ_{1} \end{array}) \\ where \\ σ_{1} = σ_{2} - \frac{5 Y^{2}}{3 (X^{2} + Y^{2}) σ_{2}} \\ σ_{2} = {(\frac{5 Y^{3}}{2 (X^{2} + Y^{2})} + \sqrt{\frac{25 Y^{6}}{4 {(X^{2} + Y^{2})}^{2}} + \frac{125 Y^{6}}{27 {(X^{2} + Y^{2})}^{3}}})}^{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