MATLAB,Python,Scilab,Julia比較第3章その116【射影変換、アフィン変換合成⑧】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia3-backnumber/

はじめに
登場人物
【再掲】数式とパラメータ
Juliaコード
処理結果
考察
まとめ

はじめに

射影変換とアフィン変換の合成をプログラムで実現する。
今回はJulia。

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

【再掲】数式とパラメータ

太郎くん

まずは、数式とパラメータを再掲。

射影(逆)変換

\(
\begin{eqnarray}
\begin{bmatrix}
x\prime\\
y\prime\\
1
\end{bmatrix}&=&
\begin{bmatrix}
a&b&c\\
d&e&f\\
g&h&1
\end{bmatrix}^{-1}
\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(\theta)&-\sin(\theta)&0\\
\sin(\theta)&\cos(\theta)&0\\
0&0&1
\end{bmatrix}^{-1}
\begin{bmatrix}
x\\
y\\
1
\end{bmatrix}
\end{eqnarray}
\)

パラメータ

射影変換

\(
\begin{eqnarray}
(-1,-1)→&(-0.5,-0.8)\\
(-1,1)→&(-0.8,0.8)\\
(1,1)→&(1,1)\\
(1,-1)→&(0.4,-1)\\
\end{eqnarray}
\)

伸縮アフィン
　縦横を1/2へ縮小

移動アフィン
　x軸方向+0.5

回転アフィン
　+60°

フクさん

これを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 homography(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_homography = matrix * points;
    
    # 画像と同一形状の2次元配列に変換元座標配列を生成
    dx = reshape(points_homography[1,:],hight,width);
    dy = reshape(points_homography[2,:],hight,width);
    ds = reshape(points_homography[3,:],hight,width);
    dx = dx./ds;
    dy = dy./ds;
    
    # 変換元座標をピクセル位置に変換
    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 homography_affine()
    # 入力画像の読み込み
    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, 100, 100);
    
    sx = 0.5;
    sy = 0.5;
    tx = 0.5;
    ty = 0;
    theta = 60/180*pi;
    
    x0=-1; y0=-1;   # 左上
    x1=-1; y1= 1;   # 左下
    x2= 1; y2= 1;   # 右下
    x3= 1; y3=-1;   # 右上

    x0t=-0.5; y0t=-0.8;	# 左上変換先
    x1t=-0.8; y1t= 0.8;	# 左下変換先
    x2t= 1; y2t= 1;	# 右下変換先
    x3t= 0.4; y3t=-1;	# 右上変換先
    
    mat = [x0 y0 1  0  0 0 -x0*x0t -y0*x0t;
            0  0 0 x0 y0 1 -x0*y0t -y0*y0t;
           x1 y1 1  0  0 0 -x1*x1t -y1*x1t;
            0  0 0 x1 y1 1 -x1*y1t -y1*y1t;
           x2 y2 1  0  0 0 -x2*x2t -y2*x2t;
            0  0 0 x2 y2 1 -x2*y2t -y2*y2t;
           x3 y3 1  0  0 0 -x3*x3t -y3*x3t;
            0  0 0 x3 y3 1 -x3*y3t -y3*y3t];
    dst = [x0t y0t x1t y1t x2t y2t x3t y3t]';
    res = inv(mat)*dst;

    homo_matrix = inv([ res[1] res[2] res[3];
                        res[4] res[5] res[6];
                        res[7] res[8] 1]);
    
    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];
    
    matrix = (homo_matrix*scaling_matrix
        *translation_matrix*rotation_matrix);
    
    homography_img = homography(img, matrix);
    
    save("dog_homography_affine_j.jpg",
        colorview(Gray, min.(abs.(homography_img),1)));
end

homography_affine();