MATLAB,Python,Scilab,Julia比較 第3章 その41【非極大値抑制⑥】

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

はじめに
登場人物
動作確認用処理【再掲】
MATLABコード
処理結果
考察
まとめ

はじめに

非極大値抑制をプログラムで実現してみる。
今回はMATLABで実施。

登場人物

博識フクロウのフクさん

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

動作確認用処理【再掲】

太郎くん

で、今回から非極大値抑制をプログラミングで実現だね。

フクさん

まずはMATLABだな。
処理手順を再掲しよう。

Sobelフィルタ等の微分フィルタで以下を推定
- x軸、y軸の濃淡変化量
- 変化強度(ノルム)
「x軸、y軸の濃淡変化量」から勾配方向角を推定
- arntan関数を利用
勾配方向を垂直(UD)、水平(LR)、斜め右上から右下(RULD)、斜め左上から右下の4パターンに丸め。
勾配方向角に応じて極大値評価をして非極大値だったら「変化強度(ノルム)」を0値埋め
画像出力

MATLABコード

フクさん

MATLABコードは以下になる。

convolution2d.m

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

non_maximum_suppression.m

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
end

non_maximum_suppression_test.m

function []=non_maximum_suppression_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 = 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 = atan2(Gy, Gx) * 180 / pi;
    imwrite(uint8(abs(G)),'dog_nms_G.jpg'); % 255オーバーあり
    imwrite(uint8(abs(theta)),'dog_nms_theta.jpg'); % -180～180
    
    % 極大値以外を除去(Non maximum Suppression)
    G_nms = non_maximum_suppression(G, theta);
    imwrite(uint8(G_nms),'dog_nms.jpg'); % 255オーバーあり 
end