なめこ備忘録

プログラミングに関する備忘録や経験したこと,考えたことなど好き勝手書きます.

[tex: ]

画像処理.js -エッジ検出 3

こんにちは,なめこです.

jsでの画像処理3回目です.
前回の記事はこちらになります.

nanameko.hatenablog.com

今回はガウス関数を使ったエッジ検出を2種類実装していきます.

2つ目は普通のエッジ検出とはちょっと毛色が違う感じです.

DoG(Difference of Gaussian)

アルゴリズム

その名の通りGaussianの差分を元にエッジを検出します.

概要としては簡単で,分散の違う2つのガウシアンフィルタを適用した画像の差分を取り,その差が大きいものをエッジとします.

参考までに以下にグラフを示します.

f:id:NAMEKO:20190224181241j:plain

赤と青が分散の違う2つのガウシアン関数,緑がその差をとったグラフです.

ガウシアンフィルタはフィルタの総和が1になるように調整されているので,差分を取ると総和は0となります.
周辺画素の値と対象画素の値の差が小さい(エッジがない)と0に近づき,大きい(エッジがある)と適用後の値も大きくなります.

実装

let sig1 = 1.3;
let sig2 = 3.2;

const filterSize = parseInt(sig2*4+1);
const filterSizeHalf = (filterSize-1)/2;

function gaussian(sig, x, y){
    return math.exp(-(x*x+y*y)/2/sig/sig)/2/math.PI/sig/sig;
}

// DoGフィルタの作成
const filter1 = Array.from(new Array(filterSize)).map((v, i) =>
    Array.from(new Array(filterSize)).map((v, j) => 
        gaussian(sig1,i-filterSizeHalf,j-filterSizeHalf)
        ));
const filter2 = Array.from(new Array(filterSize)).map((v, i) =>
    Array.from(new Array(filterSize)).map((v, j) => 
        gaussian(sig2,i-filterSizeHalf,j-filterSizeHalf)
        ));

// 各フィルタの総和計算
const sumFilter1 = filter1.reduce((pre ,cur) => pre+cur.reduce((pre, cur) => pre+cur, 0),0);
const sumFilter2 = filter2.reduce((pre ,cur) => pre+cur.reduce((pre, cur) => pre+cur, 0),0);

// フィルタを正規化しつつ合成
const filter = map(filter1, (v, i, j) => v/sumFilter1-filter2[j][i]/sumFilter2);

// DoGフィルタ適用
const imgDog = conv.conv(img, image.sizes, filter, [filterSize,filterSize], {mode:conv.EXPAND});

// abs
const maxVal = max(imgDog);
const imgDogAbs = map(imgDog, v => v<0?0:(v/maxVal*255));

// 2値化
const imgDogAbsBin = bin(imgDogAbs, {threshold:20, maxVal:255});

2つのフィルタの正規化処理があるので少し長いです.

計算量を減らすために,別々にフィルタをかけてから差分を取るような方法ではなく,合成したフィルタをかけるようにしてあります.

フィルタ適用後の負の値を使用すると一つのエッジに対して複数回反応する(フィルタの両側に負値があるため)ので,負の値は切り捨てています.
また最大値が255となるような正規化も同時に行なっています.

結果を見やすくするために2値化処理も加えました.

実行結果

DoGフィルタの適用結果です.

計算量や実装の難しさはSobelフィルタ程度ですが,陰影への反応も少なく,主要なエッジは綺麗に検出できています.
ただしCanny法と同様にジグザグやグラデーションがかかってるようなエッジに関しては検出できていません.

f:id:NAMEKO:20190224160226p:plain
結果

f:id:NAMEKO:20190224160303p:plain
2値化したもの

FDoG(Flow-Based Difference of Gaussian)

FDoGはこちらの論文にある手法の一部です.

http://www.cs.umsl.edu/~kang/Papers/kang_tvcg09.pdf

アルゴリズム

DoGは等方性によって,エッジが途切れやすいなどの問題があります.
FDoGはエッジのベクトル方向を加味することでより綺麗なエッジ検出を可能とします(Sobelフィルタに対するCanny法のようなものです.たぶん).

ちょっと長くなるので今回は詳細は省きますが一連の流れとしては以下のような感じになります(詳細説明は全体の実装をした際にまとめてします).

  1. ノイズ除去
  2. エッジの接ベクトルを元にした流れ場(ETF)を求める
  3. 直交ベクトル方向に1次元DoGフィルタを適用
  4. ETFの曲線に沿って1次元Gaussianフィルタを適用

ノイズ除去にはGaussianフィルタを使用します.

2ではSobelフィルタでETFの元となるベクトル場を求め,ETFを求める処理を3回ほど繰り返します.

接ベクトルの直交方向にDoGフィルタをかけることでより効率的なエッジの検出ができます.
また,接ベクトル方向の曲線に沿ってGaussianフィルタをかけることで途切れているエッジをつなげることができます.

3,4の処理は複数回繰り返すことで結果が良くなるのでこちらも3回ほど繰り返します.

イラスト風の画像となるという特徴を持っており,処理後の画像も通常のエッジ検出とは違い,白地に黒で描かれたものになります.

実装

// 関数定義
const G = (r, sig) => math.exp(-(r*r)/2/sig/sig)/math.sqrt(2*math.PI)/sig;
const Wm = (g, g_) => (1 + math.tanh(g_ - g))/2;
const Wd = (tx, ty, tx_, ty_) => tx*tx_ + ty*ty_;

// 定数定義
let sigM = 2.5;
let sigC = 1.1;
let tau = 0.5;
let range = 3;

const sigS = sigC*1.6;
const rho = 0.997;

const T = parseInt(sigS*2+0.5);
const S = parseInt(sigM*2+0.5);

// 画像データ取得
const image = cv.imread(argv[0], cv.CV_8UC1);
const img = image.getDataAsArray();

const [height, width] = image.sizes;

// ノイズ除去
const imgGaus = gaussianBlur(img, {kSize:3});

// sobelフィルタ
let {sobelImg:g, sobelX:gx, sobelY:gy} = sobel(imgGaus);
sobelImg = normalize(g);
map(g, (v, i, j) => {
    if(v == 0)return;
    gx[j][i] /= v;
    gy[j][i] /= v;
    return;
});

// gx,gyを90度回転したものをtx,tyの初期値とする
let tx = map(gy, v => -v);
let ty = gx.concat();

// 正規化
g = normalize(g);

// ETF計算を複数回回すループ
for(let n=0;n<3;n++){
    console.log(`Edge Tangent Flow ${n} iteration`);

    const tx_cur = tx.concat();
    const ty_cur = ty.concat();

    // ETF計算
    map(g, (v, i, j) => {
        let tx_sum = 0;
        let ty_sum = 0;

        const tx_num = tx_cur[j][i];
        const ty_num = ty_cur[j][i];

        for(let l=-range; l<=range; l++){
            const j_ = j+l;
            if(j_ < 0 || j_ >= height)continue;
            for(let k=-range; k<=range; k++){
                const i_ = i+k;
                if(i_ < 0 || i_ >= width)continue;
                
                // 対象外領域は弾く
                const r = math.sqrt(k*k+l*l);
                if(r > range)continue;

                // 重み計算
                const w = Wm(v, g[j_][i_])*Wd(tx_num, ty_num, tx_cur[j_][i_], ty_cur[j_][i_]);

                tx_sum += tx_cur[j_][i_]*w;
                ty_sum += ty_cur[j_][i_]*w;
            }
        }

        // 正規化
        const r = math.sqrt(tx_sum*tx_sum + ty_sum*ty_sum);
        if(r != 0){
            tx_sum /= r;
            ty_sum /= r;
        }

        // tx,tyを更新
        tx[j][i] = tx_sum;
        ty[j][i] = ty_sum;
    });
}

// gx,gyを更新
gx = ty.concat();
gy = map(tx, v => -v);

// 定数定義
const THRESHOLD = math.cos(math.PI/8);

let He = img.concat();

// FDoGを複数回回すループ
for(let n=0;n<3;n++){
    console.log(`Flow Based DoG ${n} iteration`);

    // 直交方向に1次元DoGフィルタをかける
    const Hg = map(He, (_, i, j, img) => {
        const x = gx[j][i];
        const y = gy[j][i];
        if(x==0 && y==0)return 0;

        let Csum = 0;
        let Ssum = 0;
        let GCsum = 0;
        let GSsum = 0;

        // 角度に応じて探索方向を定める
        let dk = 0;
        let dl = 0;
        if(math.abs(x) < THRESHOLD){
            dl = 1;
        }
        if(math.abs(y) < THRESHOLD){
            dk = 1;
        }
        if(x*y < 0){
            dl = -1;
        }

        // DoG計算
        for(let l=-T; l<=T; l++){
            const di = l*dk;
            const dj = l*dl;
            const i_ = i + di;
            const j_ = j + dj;
            if(j_ < 0 || j_ >= height || i_ < 0 || i_ >= width)continue;

            const r = math.sqrt(di*di + dj*dj);
            gc = G(r, sigC);
            gs = G(r, sigS);
            GCsum += gc;
            GSsum += gs;
            Csum += img[j_][i_] * gc;
            Ssum += img[j_][i_] * gs;
        }
        return Csum/GCsum - rho*Ssum/GSsum;
    });

    // 流れ場に沿ってGaussianフィルタをかける
    He = map(Hg, (v, i, j, img) => {
        let x = tx[j][i];
        let y = ty[j][i];
        // 流れがなければ白地
        if(x==0 && y==0)return 255;

        sum = v;
        Gsum = G(0, sigM);

        // 接ベクトルの順方向と逆方向に探索
        for(let m=-1; m<2; m+=2){
            let r = 0;
            let i_ = i;
            let j_ = j;

            while(r < S){
                let di = 0;
                let dj = 0;
                if(math.abs(x) < THRESHOLD){
                    dj = m;
                }
                if(math.abs(y) < THRESHOLD){
                    di = m;
                }
                if(x*y < 0){
                    dj *= -1;
                }
                j_ += dj;
                i_ += di;
                if(j_ < 0 || j_ >= height || i_ < 0 || i_ >= width)break;

                r += math.sqrt(di*di + dj*dj);
                const g = G(r, sigM);
                Gsum += g;
                sum += img[j_][i_]*g;

                x = tx[j_][i_];
                y = ty[j_][i_];
                if(x==0 && y==0)break;
            }
        }

        sum /= Gsum;

        // 閾値処理
        return (sum<0 && (1 + math.tanh(sum)) < tau) ? 0 : 255;
    });
}

実行結果

実装時の定数をそのまま使用した場合の結果が以下になります.
DoGフィルタでは検出できなかったジグザグとしたエッジやグラデーションのかかった部分もこの手法では綺麗に検出できています.
またイラスト風というだけあって,陰影も表現されているのが面白いです.

f:id:NAMEKO:20190224160047p:plain
tau=0.5

比較のためtauを0.1, 0.9としたものを示します.tauを小さくすると陰影などの線が少なくなります.

f:id:NAMEKO:20190224160121p:plain
tau=0.1

f:id:NAMEKO:20190224160139p:plain
tau=0.9

まとめ

手軽に綺麗なエッジを取れるDoGとその改良版のFDoGを実装しました.

FDoGは処理は若干複雑ですが,ETF計算以外はガウス関数しか使っていないのが面白いです. またFDoGは自然画像に対して実行すると他の手法とは少し違った特徴が出てくるので試してみるといいかもしれません.

今回の実装含め,これまでのエッジ検出のプログラムは以下のGitHubのedgeディレクトリに置いてあります.

github.com

今回で一旦エッジ検出はクローズとなります.また面白いものが見つかったら書くかもしれません.

次回からは劣化画像修復に移る予定です. 今回までは畳み込みや対象画素付近の探索がメインでしたが,次からは行列計算がメインになります.

参考文献

【画像処理】DoGフィルタの原理・特徴・計算式 | アルゴリズム雑記