バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia3-backnumber/
はじめに
非極大値抑制にHysteresis Thresholdを加えた、Canny法を実施して、2値化を行う。
今回は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
【再掲】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();
処理結果
フクさん
処理結果は以下。
考察
太郎くん
これも2値化できてるし大丈夫そうだね。
フクさん
ちなみに、一度処理を走らせてJITコンパイルが通り切った後はJuliaがもっとも高速だったな。
太郎くん
ほう。
この手の処理はMATLABが一番早いかと思ってたけど、Juliaの方が速くなるのか。
フクさん
当然、環境依存はあるだろうが、
仕組み的にJuliaはJITコンパイル後はネイティブコードで動作するからね。
速くなる理屈ははっきりはしてるね。
まとめ
フクさん
まとめだよ。
- 非極大値抑制にHysteresis Thresholdを加えた、Canny法による2値化をJuliaで実施。
- 環境依存はあるかもしれないが、処理としてはMATLABよりも高速。
- JITコンパイル後はネイティブコードで動作するため。
バックナンバーはこちら。
コメント