【入門】Hysteresis Threshold(Julia)【数値計算】

【入門】Hysteresis Threshold(Julia)【数値計算】 数値計算
【入門】Hysteresis Threshold(Julia)【数値計算】

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

はじめに

の、
MATLAB,Python,Scilab,Julia比較 第3章 その56【Hysteresis Threshold⑦】

非極大値抑制にHysteresis Thresholdを加えた、Canny法を実施して、2値化を行う。
今回はJuliaで実施する。

【再掲】Hysteresis Thresholdを実施するための手順

まずは手順を再掲。

  • Sobelフィルタ等の微分フィルタで以下を推定
    • x軸、y軸の濃淡変化量
    • 変化強度(ノルム)
  • 「x軸、y軸の濃淡変化量」から勾配方向角を推定
    • arntan関数を利用
  • 勾配方向を垂直(UD)、水平(LR)、斜め右上から右下(RULD)、斜め左上から右下の4パターンに丸め。
  • 勾配方向角に応じて極大値評価をして非極大値だったら「変化強度(ノルム)」 を0値埋め
  • Hysteresis Threshold
    • High以上は白
    • Low未満は黒
    • High-Lowの間の場合は周辺を探査し、エッジが居れば白、いなければ黒
  • 画像出力

今回は、これをJuliaで実現する。

Juliaコード

Juliaコードは以下になる。

using Images

# 畳み込み演算
function convolution2d(img, kernel)
    (n, m) = size(kernel); # カーネルサイズ取得
    
    # カーネル中心からみた幅
    dy = Int64((n-1)/2);  # カーネル上下幅
    dx = Int64((m-1)/2);  # カーネル左右幅
    
    (h, w) = size(img);   # イメージサイズ
    out = zeros(h, w);    # 出力用イメージ
    
    # 畳み込み
    for y = dy+1:(h - dy)
        for x = dx+1:(w-dx)
            out[y, x] = sum( img[y-dy:y+dy, x-dx:x+dx].*kernel );
        end
    end
    return out;
end

# Non maximum Suppression
function non_maximum_suppression(G, theta)
    (h, w) = size(G);
    out = copy(G);

    # 勾配方向を4方向(LR,UD,RULD,LURD)に近似
    theta[ -22.5 .<= theta .<   22.5] .= 0;   # LR    ─
    theta[  22.5 .<= theta .<   67.5] .= 45;  # RULD  /
    theta[  67.5 .<= theta .<  112.5] .= 90;  # UD    │
    theta[ 112.5 .<= theta .<  157.5] .= 135; # LURD  \
    theta[ 157.5 .<= theta .<  180.0] .= 0;   # LR    ─
    theta[-180.0 .<= theta .< -157.5] .= 0;   # LR    ─
    theta[-157.5 .<= theta .< -112.5] .= 45;  # RULD  /
    theta[-112.5 .<= theta .<  -67.5] .= 90;  # UD    │
    theta[ -67.5 .<= theta .<  -22.5] .= 135; # LURD  \

    # 現画素の勾配方向に接する2つの画素値を比較し、現画素が極大値でなければ0にする。
    for y = 2:(h - 1)
        for x = 2:(w - 1)
            if theta[y,x] == 0    # LR    ─
                if (G[y,x] < G[y,x+1]) || (G[y,x] < G[y,x-1])
                    out[y,x] = 0;
                end
            elseif theta[y,x] == 45 # RULD  /
                if (G[y,x] < G[y-1,x+1]) || (G[y,x] < G[y+1,x-1])
                    out[y,x] = 0;
                end
            elseif theta[y,x] == 90 # UD    │
                if (G[y,x] < G[y+1,x]) || (G[y,x] < G[y-1,x])
                    out[y,x] = 0;
                end
            else                  # LURD  \
                if (G[y,x] < G[y+1,x+1]) || (G[y,x] < G[y-1,x-1])
                    out[y,x] = 0;
                end
            end
        end
    end
    return out
end

function hysteresis_threshold(img, low, high, r)
    (h, w) = size(img);
    out = copy(img);
    
    for y = (1+r):(h-r-1)
        for x = (1+r):(w-r-1)
           # 最大閾値より大きければ「エッジ」
           if img[y,x] >= high
               out[y,x] = 1.0;
           # 最小閾値より小さければ「非エッジ」
           elseif img[y,x] < low
               out[y,x] = 0;
           # 最小閾値と最大閾値の間で、半径rの範囲内に「エッジ」が1つでもあればエッジと判定
           else
               if maximum(img[y-r:y+r+1, x-r:x+r+1]) >= high
                   out[y,x] = 1.0;
               else
                   out[y,x] = 0;
               end
           end
        end
    end
    return out;
end

function canny_test()
    # 入力画像の読み込み
    img = channelview(load("dog.jpg"));
    r = img[1,:,:];
    g = img[2,:,:];
    b = img[3,:,:];
    
    # ガウシアンフィルタ用のkernel
    kernel_gauss = [ 1/16  2/16  1/16; 
                     2/16  4/16  2/16; 
                     1/16  2/16  1/16];
    # Sobelフィルタ用のKernel
    kernel_sx = [-1  0  1;
                 -2  0  2;
                 -1  0  1];
    kernel_sy = kernel_sx';
    
    # SDTVグレースケール
    gray_sdtv = 0.2990 * r + 0.5870 * g + 0.1140 * b;

    # ガウシアンフィルタ
    img_g = convolution2d(gray_sdtv, kernel_gauss);

    # Sobelフィルタ
    Gx = convolution2d(img_g, kernel_sx);
    Gy = convolution2d(img_g, kernel_sy);

    # 勾配強度と角度
    G = sqrt.( Gx.^2 + Gy.^2 );
    theta = atan.(Gy, Gx) * 180 / pi;
    
    #  極大値以外を除去(Non maximum Suppression)
    G_nms = non_maximum_suppression(G, theta);
    
    # Hysteresis Threshold(最小閾値除去、最大閾値以上を残し、且つそこと繋がっている最小閾値以上を残す)
    G_canny = hysteresis_threshold(G_nms, 0.1176, 0.2549, 1); # 0 or 1に2値化
    save("dog_canny_j.jpg",colorview(Gray, min.(G_canny,1)));
    
    return;
end

canny_test();

処理結果

処理結果は以下。

犬と自転車(Canny法)Julia

考察

これも2値化できてるし大丈夫そう。

ちなみに、一度処理を走らせてJITコンパイルが通り切った後はJuliaがもっとも高速。
この手の処理はMATLABが一番早いイメージがあるが、Juliaの方が速くなることも多い。
当然、環境依存はあるだろうが、
仕組み的にJuliaはJITコンパイル後はネイティブコードで動作するから、
速くなる理屈ははっきりはしてる。

まとめ

  • 非極大値抑制にHysteresis Thresholdを加えた、Canny法による2値化をJuliaで実施。
  • 環境依存はあるかもしれないが、処理としてはMATLABよりも高速。
    • JITコンパイル後はネイティブコードで動作するため。

MATLAB、Python、Scilab、Julia比較ページはこちら

コメント

タイトルとURLをコピーしました