MATLAB,Python,Scilab,Julia比較 第3章 その55【Hysteresis Threshold⑥】

MATLAB,Python,Scilab,Julia比較 第3章 その55【Hysteresis Threshold⑥】 数値計算
MATLAB,Python,Scilab,Julia比較 第3章 その55【Hysteresis Threshold⑥】

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

はじめに

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

登場人物

博識フクロウのフクさん

指差しフクロウ

イラスト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

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

太郎くん
太郎くん

まずは手順を再掲。

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

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

Scilabコード

フクさん
フクさん

Scilabコードは以下になる。

convolution2d.sci

function out = convolution2d(img, Kernel)
    [m, n] = size(Kernel);  // カーネルサイズ取得
    
    // カーネル中心からみた幅
    dy = int32((m-1)/2);  // カーネル上下幅
    dx = int32((n-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( double(img(y-dy:y+dy, x-dx:x+dx)).*Kernel );
        end
    end
endfunction

non_maximum_suppression.sci

function out = non_maximum_suppression(G, theta)
    [h, w] = size(G);
    out = G;
    
    // 勾配方向を4方向(LR,UD,RULD,LURD)に近似
    theta( -22.5 <= theta & theta <   22.5) = 0;   // LR    ─
    theta(  22.5 <= theta & theta <   67.5) = 45;  // RULD  /
    theta(  67.5 <= theta & theta <  112.5) = 90;  // UD    │
    theta( 112.5 <= theta & theta <  157.5) = 135; // LURD  \
    theta( 157.5 <= theta & theta <  180.0) = 0;   // LR    ─
    theta(-180.0 <= theta & theta < -157.5) = 0;   // LR    ─
    theta(-157.5 <= theta & theta < -112.5) = 45;  // RULD  /
    theta(-112.5 <= theta & theta <  -67.5) = 90;  // UD    │
    theta( -67.5 <= theta & 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
endfunction

hysteresis_threshold.sci

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

canny_test.sci

function [] = canny_test()
    // 入力画像の読み込み
    img = imread('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 = uint8([0.2990 * double(r) ...
                + 0.5870 * double(g) + 0.1140 * double(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, 30, 65, 1); // 0 or 255に2値化
    imwrite(uint8(G_canny),'dog_canny.jpg'); 
endfunction

処理結果

フクさん
フクさん

処理結果は以下。

犬と自転車(Canny法)Scilab

メモリ不足問題

フクさん
フクさん

ちなみに環境依存ではあるようだが、Scilabでこの処理をさせる場合、
メモリ不足になって処理が完了できないことがある。
この場合、Scilabのメニューの「編集」から「設定」選んで、Javaピープメモリを増加方向に調整してみると解決することが多い。

Scilab Javaヒープメモリ調整、Javaで利用可能なメモリ(単位:MB)を指定
フクさん
フクさん

画面上にも書いてあるが、再起動で反映される設定である点に注意だ。

太郎くん
太郎くん

なんかマニアックな問題を抱えてるな・・・。

考察

太郎くん
太郎くん

なんかメモリ不足問題とかあるようだけど、
処理としてはちゃんと2値化出来てるね。
ソースコードもMATLABとほぼ一緒だし。

フクさん
フクさん

まぁ、今回も、コード自体はMATLABからコピペしたようなもんだしね。

まとめ

フクさん
フクさん

まとめだよ。

  • 非極大値抑制にHysteresis Thresholdを加えた、Canny法による2値化をScilabで実施。
  • 基本的にはいままでのコードを再利用。
  • 環境によってはメモリ不足に陥るのでJavaヒープメモリを調整する必要がある。

バックナンバーはこちら。

コメント

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