【入門】アフィン逆変換時の行列合成(Julia)【数値計算】

MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

はじめに
【再掲】プログラムで実現したいアフィン行列の合成
Juliaコード
処理結果
考察
まとめ

はじめに

の、

MATLAB,Python,Scilab,Julia比較第3章その86【アフィン行列の合成⑩】

を書き直したもの。

アフィン変換のアフィン行列の合成の話。
今回はJuliaで実施する。

【再掲】プログラムで実現したいアフィン行列の合成

まずは、実現したいアフィン行列の合成。

\(
\begin{eqnarray}
\begin{bmatrix}
x\prime\\
y\prime\\
1
\end{bmatrix}=
\begin{bmatrix}
1 & m_x & 0 \\
m_y & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\cos(30^\circ) & -\sin(30^\circ) & 0 \\
\sin(30^\circ) & \cos(30^\circ) & 0 \\
0 & 0 & 1
\end{bmatrix} \\
\begin{bmatrix}
1 & 0 & T_x \\
0 & 1 & T_y \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
S_x & 0 & 0 \\
0 & S_y & 0 \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
x\\
y\\
1
\end{bmatrix}
\end{eqnarray}
\)

そして、アフィン逆変換のアルゴリズム都合で実際に使用する数式がこれ。

\(
\begin{eqnarray}
\begin{bmatrix}
x\\
y\\
1
\end{bmatrix}=
\begin{bmatrix}
S_x & 0 & 0 \\
0 & S_y & 0 \\
0 & 0 & 1
\end{bmatrix}^{-1}
\begin{bmatrix}
1 & 0 & T_x \\
0 & 1 & T_y \\
0 & 0 & 1
\end{bmatrix}^{-1}
\begin{bmatrix}
\cos(30^\circ) & -\sin(30^\circ) & 0 \\
\sin(30^\circ) & \cos(30^\circ) & 0 \\
0 & 0 & 1
\end{bmatrix}^{-1}\\
\begin{bmatrix}
1 & m_x & 0 \\
m_y & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}^{-1}
\begin{bmatrix}
x\prime\\
y\prime\\
1
\end{bmatrix}
\end{eqnarray}
\)

これをJuliaで実現する。

Juliaコード

Juliaコードは以下になる。

using Images

function meshgrid(xin,yin)
    nx=length(xin)
    ny=length(yin)
    xout=zeros(ny,nx)
    yout=zeros(ny,nx)
    for jx=1:nx
        for ix=1:ny
            xout[ix,jx]=xin[jx]
            yout[ix,jx]=yin[ix]
        end
    end
    return (x=xout, y=yout)
end

# アフィン変換関数
function affine_transformation(img, matrix)
    # 画像サイズ取得
    (hight, width) = size(img);
    
    # 中心を0とした座標系を生成
    x_axis = range(-1, 1, length=width);
    y_axis = range(-1, 1, length=hight);
    (xim,yim) = meshgrid(x_axis, y_axis);
    
    # 座標x,y,1の3次元ベクトルの配列
    # n(:)表記で列ベクトル化したあとに転置して行ベクトル化
    points = [xim[:]';yim[:]'; ones(1, width*hight)];
    
    # 変換元座標算出(アフィン逆変換)
    points_affine = matrix * points;
    
    # 画像と同一形状の2次元配列に変換元座標配列を生成
    dx = reshape(points_affine[1,:],hight,width);
    dy = reshape(points_affine[2,:],hight,width);
    
    # 変換元座標をピクセル位置に変換
    v = UInt32.(round.(min.(max.((dx.+1)*width/2, 1), width )));
    h = UInt32.(round.(min.(max.((dy.+1)*hight/2, 1), hight )));
    
    # 元画像と変換元座標を元に変換先へコピー
    affine_img = img[h+(v.-1)*hight];
    
    return affine_img
end

# キャンパス拡張
function canvas_expansion(img, x, y)
    (H, W) = size(img);
    WID = W+x;
    HID = H+y;
    e_img = zeros(HID, WID);
    e_img[Int32((HID-H)/2)+1:Int32((HID+H)/2), Int32((WID-W)/2)+1:Int32((WID+W)/2)] = img;    
    img = e_img;
    
    return img
end

function affine_transformation_test()
    # 入力画像の読み込み
    img = channelview(load("dog.jpg"));
    
    r = img[1,:,:];
    g = img[2,:,:];
    b = img[3,:,:];

    # SDTVグレースケール
    img = 0.2990 * r + 0.5870 * g + 0.1140 * b;
    
    # キャンパス拡張
    img = canvas_expansion(img, 800, 800);
    
    sx = 1;
    sy = -1;
    tx = 0.5;
    ty = 0;
    theta = 150/180*pi;
    mx = tan(0/180*pi);
    my = tan(0/180*pi);
    
    scaling_matrix     = inv([ sx   0  0;
                                0  sy  0;
                                0   0  1]);
    
    translation_matrix = inv([ 1 0  tx;
                               0 1 -ty;
                               0 0   1]);
    
    rotation_matrix    = [ cos(theta) -sin(theta)  0;
                           sin(theta)  cos(theta)  0;
                                   0           0   1];
    
    shear_matrix       = inv([  1 -mx  0;
                              -my   1  0;
                                0   0  1]);
    
    matrix = scaling_matrix*translation_matrix*rotation_matrix*shear_matrix;
    
    affine_img = affine_transformation(img, matrix);
    
    save("dog_affine_j.jpg",colorview(Gray, min.(abs.(affine_img),1)));
end

affine_transformation_test();