概要

Cannyでエッジ検出するプログラムをJavaで書いてみました。OpenCVを使えば一発なのでしょうが、趣味のコーディングなので、そこは頑張って自分でひと通り実装します。大まかな流れは以下の通り。

1. 入力画像をグレースケールに変換
2. ガウシアンフィルタでノイズ除去
3. ソーベルフィルタで1階微分
4. 細線化
5. ヒステリシス処理

Filter.java

public static BufferedImage canny(BufferedImage input, boolean debug) throws IOException {
    BufferedImage result = Color.grayScaleNTSC(input, true);
    result = Filter.gaussianBlur(result, (float) 1.3, 5, true);

    float[][] tan = new float[input.getWidth()][input.getHeight()];
    result = Filter.sobel(result, tan, true);
    result = Filter.nonMaximumSupression(result, tan, true);

    return hysteresis(result, (float) 20, (float) 10);
}

個別の手法の詳細は他のサイトや書籍を参考にしてください。掲載しているコードは、何となく動くことは確認していますが、バグがあったらごめんなさい。最後に追記あり。

それでは、処理の順番に沿って見ていきます。

グレースケール変換

RGBそれぞれの値にどんな数字を掛け合わせるかは諸論あるみたいですが、今回はあまり拘りがないのでネットで見かけた数字を使ってみます。ツール的なものは自作のColorクラスにまとめてあるので、グレースケール変換もそちらに実装してみます。気分の問題以外に別のクラスに実装する理由はなさそう。

Color.java

public static BufferedImage grayScaleNTSC(BufferedImage input, boolean debug) throws IOException {
    BufferedImage result = new BufferedImage(input.getWidth(), input.getHeight(), BufferedImage.TYPE_INT_RGB);
    for(int y=0; y<input.getHeight(); y++){
        for(int x=0; x<input.getWidth(); x++){
            int pixel = input.getRGB(x, y);
            float gray = (float)0.298839*r(pixel) + (float)0.586811*g(pixel) + (float)0.114350*b(pixel);
            result.setRGB(x, y, rgb(round(gray), round(gray), round(gray)));
        }
    }

    if(debug){
        ImageIO.write(result, "jpg", new File("grayscale.jpg"));
    }
    return result;
}

自作のColorクラスでは、グレースケール変換の他に、BufferedImageのgetRGBメソッドが返す整数値からR/G/Bそれぞれの値を取り出すメソッドや、R/G/Bの値からBufferedImageのsetRGBメソッドに渡すことができる整数値を作るメソッド、R/G/Bの値を整数に丸める時に0〜255に収めるメソッドなども実装しており、これらは以下のあちこちで使っています。

Color.java

public static int r(int pixel){
    return pixel >> 16 & 0xff;
}

public static int g(int pixel){
    return pixel >> 8 & 0xff;
}

public static int b(int pixel){
    return pixel & 0xff;
}

public static int rgb(int r, int g, int b){
    return 0xff000000 | r << 16 | g << 8 | b;
}

public static int round(float x) {
    if (x > 255) {
        return 255;
    } else if (x < 0) {
        return 0;
    } else {
        return Math.round(x);
    }
}

レナ画像で実行してみると、こんな感じに。それっぽい。

ガウシアンフィルタ

シグマの値とカーネルサイズに合わせてフィルタを作成、畳み込みします。畳み込み計算はあちこちで使いそうなので、別のメソッドとして独立させます。

Filter.java

public static BufferedImage gaussianBlur(BufferedImage input, float sigma, int size, boolean debug) throws IOException {
    BufferedImage result = new BufferedImage(input.getWidth(), input.getHeight(), BufferedImage.TYPE_INT_RGB);

    float[][] kernel = new float[size][size];
    float weight = gaussianKernel(sigma, kernel);

    for (int y = 0; y < input.getHeight(); y++) {
        for (int x = 0; x < input.getWidth(); x++) {
            result.setRGB(x, y, convolution(input, kernel, weight, x, y));
        }
    }

    if (debug) {
        ImageIO.write(result, "jpg", new File("gaussian.jpg"));
    }
    return result;
}

Filter.java

private static float gaussianKernel(float sigma, float[][] kernel) {
    int size = kernel.length;
    int mid = (size - 1) / 2;
    float tmp = 0;
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            kernel[i][j] = gaussian(i - mid, j - mid, sigma);
            tmp += kernel[i][j];
        }
    }
    return tmp;
}

Filter.java

private static int convolution(BufferedImage input, float[][] kernel, float weight, int x, int y) {
    int size = kernel.length;
    int mid = (size - 1) / 2;
    float tmp_r = 0;
    float tmp_g = 0;
    float tmp_b = 0;

    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            if (x - mid + i >= 0 && x - mid + i < input.getWidth() && y - mid + j >= 0 && y - mid + j < input.getHeight()) {
                int org = input.getRGB(x - mid + i, y - mid + j);
                tmp_r += kernel[i][j] * Color.r(org);
                tmp_g += kernel[i][j] * Color.g(org);
                tmp_b += kernel[i][j] * Color.b(org);
            }
        }
    }

    return Color.rgb(Color.round(tmp_r / weight), Color.round(tmp_g / weight), Color.round(tmp_b / weight));
}

Filter.java

private static float gaussian(int x, int y, float sigma) {
    return (float) (Math.exp(-(x * x + y * y) / (2 * sigma * sigma)) / (2 * Math.PI * sigma * sigma));
}

ここまで動かすとこんな感じ。画像サイズを小さくして表示しているので、ガウシアンフィルタの効果がわかりづらいですね。クリックして大きな画像で見ると、ガウシアンフィルタ適用後は、ぼかしが効いているのがわかります。

ソーベルフィルタ

これはもうフィルタが決め打ちなので、決まり切った数値でカーネルを作って、先ほどの畳み込みメソッドに叩き込むだけです。

Filter.java

public static BufferedImage sobel(BufferedImage input, float[][] tan, boolean debug) throws IOException {
    BufferedImage result = new BufferedImage(input.getWidth(), input.getHeight(), BufferedImage.TYPE_INT_RGB);
    float diffKernelVertical[][] = { { -1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } };
    float diffKernelHorizontal[][] = { { -1, -2, -1 }, { 0, 0, 0 }, { 1, 2, 1 } };

    for (int y = 0; y < input.getHeight(); y++) {
        for (int x = 0; x < input.getWidth(); x++) {
            int gy = convolution(input, diffKernelVertical, (float) 1.0, x, y);
            int gx = convolution(input, diffKernelHorizontal, (float) 1.0, x, y);

            float tmp_x = Color.r(gx);
            float tmp_y = Color.r(gy);
            float strength_r = (float)Math.sqrt(tmp_x*tmp_x + tmp_y*tmp_y);

            tmp_x = Color.g(gx);
            tmp_y = Color.g(gy);
            float strength_g = (float) Math.sqrt(tmp_x * tmp_x + tmp_y * tmp_y);

            tmp_x = Color.b(gx);
            tmp_y = Color.b(gy);
            float strength_b = (float) Math.sqrt(tmp_x * tmp_x + tmp_y * tmp_y);

            result.setRGB(x, y, Color.rgb(Color.round(strength_r), Color.round(strength_g), Color.round(strength_b)));
            if(Double.isNaN(gy/gx) || Double.isInfinite(gy/gx)){
                tan[x][y] = 0;
            }
            else{
                tan[x][y] = (float)Math.tan(gy/gx);
            }
        }
    }

    if (debug) {
        ImageIO.write(result, "jpg", new File("sobel.jpg"));
    }
    return result;
}

ここまでの経過はこんな感じ。なにやらエッジっぽいものが出始めました。

細線化

場合分けが面倒ですが、素直に書くだけです。

Filter.java

public static BufferedImage nonMaximumSupression(BufferedImage input, float[][] tan, boolean debug) throws IOException {
    BufferedImage result = new BufferedImage(input.getWidth(), input.getHeight(), BufferedImage.TYPE_INT_RGB);
    for (int y = 0; y < input.getHeight(); y++) {
        for (int x = 0; x < input.getWidth(); x++) {
            int s1 = 0;
            int s2 = 0;
            if (tan[x][y] > Math.tan(-Math.PI / 8) && Math.tan(Math.PI / 8) >= tan[x][y]) {
                // theta: -PI/8 to PI/8 -> 0
                if (x - 1 >= 0) {
                    s1 = input.getRGB(x-1, y);
                }

                if (x + 1 < input.getWidth()) {
                    s2 = input.getRGB(x+1, y);
                }
            } else if (tan[x][y] > Math.tan(Math.PI / 8) && Math.tan(3 * Math.PI / 8) >= tan[x][y]) {
                // theta: PI/8 to 3PI/8 -> PI/4
                if (x - 1 >= 0 && y + 1 < input.getHeight()) {
                    s1 = input.getRGB(x-1, y+1);
                }

                if (x + 1 < input.getWidth() && y - 1 >= 0) {
                    s2 = input.getRGB(x+1, y-1);
                }
            } else if (Math.abs(tan[x][y]) > Math.tan(3 * Math.PI / 8)) {
                // theta: 3PI/8 to 5PI/8 -> PI/2
                if (y - 1 >= 0) {
                    s1 = input.getRGB(x, y-1);
                }

                if (y + 1 < input.getHeight()) {
                    s2 = input.getRGB(x, y+1);
                }
            } else {
                // theta: 5PI/8 to 7PI/8-> 3PI/4
                if (x - 1 >= 0 && y - 1 >= 0) {
                    s1 = input.getRGB(x-1, y-1);
                }

                if (x + 1 < input.getWidth() && y + 1 < input.getHeight()) {
                    s2 = input.getRGB(x+1, y+1);
                }
            }

            int target = input.getRGB(x, y);
            int strength_t = Color.r(target) + Color.g(target) + Color.b(target);
            int strength_s1 = Color.r(s1) + Color.g(s1) + Color.b(s1);
            int strength_s2 = Color.r(s2) + Color.g(s2) + Color.b(s2);
            if (strength_s1 >= strength_t || strength_s2 >= strength_t) {
                result.setRGB(x, y, Color.rgb(0, 0, 0));
            } else {
                result.setRGB(x, y, input.getRGB(x, y));
            }
        }
    }

    if (debug) {
        ImageIO.write(result, "jpg", new File("nonMaximumSupression.jpg"));
    }
    return result;
}

実行してみると、こんな感じ。確かに線が細くなっています。

ヒステリシス処理

これも素直にそのまま書いています。ふたつの閾値の間にあるピクセルについて、今回は近傍にエッジなピクセルがひとつでもあるとエッジであると判定してしまっています（checkNeighborsメソッド）。面倒だったのでそうしてしまったのですが、ここはもう少し真面目にクラスタリングをした方が良かったな。まあ今後の課題ってことで。

Filter.java

public static BufferedImage hysteresis(BufferedImage input, float higherThreshold, float lowerThreshold) {
    BufferedImage buf = new BufferedImage(input.getWidth(), input.getHeight(), BufferedImage.TYPE_INT_RGB);
    for (int y = 0; y < input.getHeight(); y++) {
        for (int x = 0; x < input.getWidth(); x++) {
            int pixel = input.getRGB(x, y);
            float target = (float) ((Color.r(pixel) + Color.g(pixel) + Color.b(pixel)) / 3.0);
            if (target >= higherThreshold) {
                buf.setRGB(x, y, Color.rgb(255, 255, 255));
            } else if (target < lowerThreshold) {
                buf.setRGB(x, y, Color.rgb(0, 0, 0));
            } else {
                buf.setRGB(x, y, Color.rgb(255, 0, 0));
            }
        }
    }

    BufferedImage ret = new BufferedImage(buf.getWidth(), buf.getHeight(), BufferedImage.TYPE_INT_RGB);
    for (int y = 0; y < buf.getHeight(); y++) {
        for (int x = 0; x < buf.getWidth(); x++) {
            int pixel = buf.getRGB(x, y);
            if (pixel != Color.rgb(255, 0, 0)) {
                ret.setRGB(x, y, pixel);
                continue;
            }

            int white = Color.rgb(255, 255, 255);
            if (checkNeighbors(buf, x, y, white)) {
                ret.setRGB(x, y, white);
            } else {
                ret.setRGB(x, y, Color.rgb(0, 0, 0));
            }
        }
    }
    return ret;
}

Filter.java

private static boolean checkNeighbors(BufferedImage input, int x, int y, int color) {

    if (x - 1 >= 0 && input.getRGB(x - 1, y) == color) {
        return true;
    }

    if (x - 1 >= 0 && y - 1 >= 0 && input.getRGB(x - 1, y - 1) == color) {
        return true;
    }

    if (x - 1 >= 0 && y + 1 < input.getHeight() && input.getRGB(x - 1, y + 1) == color) {
        return true;
    }

    if (y - 1 >= 0 && input.getRGB(x, y - 1) == color) {
        return true;
    }

    if (y + 1 < input.getHeight() && input.getRGB(x, y + 1) == color) {
        return true;
    }

    if (x + 1 < input.getWidth() && input.getRGB(x + 1, y) == color) {
        return true;
    }

    if (x + 1 < input.getWidth() && y - 1 >= 0 && input.getRGB(x + 1, y - 1) == color) {
        return true;
    }

    if (x + 1 < input.getWidth() && y + 1 < input.getHeight() && input.getRGB(x + 1, y + 1) == color) {
        return true;
    }

    return false;
}

これで完成。

パラメータを調整したり、テキトーに済ませたところをちゃんと実装すれば、もう少し良い結果が出ると思いますが、手法をひと通り実装して動いたので、とりあえず満足。でもあちこちでグルングルン同じようなループが回っているのは非常に気にくわないので、何かしたい気になってしまう。あまり大幅にループを減らすのは難しいだろうけど。

(2019/08/22 追記）
グレースケール変換のバグを直したら、だいぶそれっぽくなってきました。ヒステリシスの閾値を20と60にしてみるとこんな感じ。

あまり褒められたものではないオマケ

レナ様は、背景ボケてるし、リアルな写真だし、アニメとかの方が取りやすいんじゃないのかなあってことで、拾い物のリヴァイ様で試した結果がこちら。グレースケールにした時点で、リヴァイ様の美しいお顔が、のっぺらぼうに…やっぱり難しい。

(2019/08/22 追記）
こちらもグレースケール変換のバグを直したら、リヴァイ様のお顔が消えなくなり、それっぽい結果が出るようになりました。

追記

2019/08/22
Color.grayScaleNTSC()にバグがあるとのご指摘をいただきました。ありがとうございます。typoしていてGとBの値を掛け合わせるという暴挙に出ていたので、修正しました。結果もだいぶ良くなった！

tan(x)のxが大きくなってしまうとどこに行くかわからないという懸念を持たれている方がおられるようですが、tan(theta) = tan(theta + nPI) (n:整数)ですから、実質、昨日の追記にある図の黒じゃないところを行ったり来たりするのと変わらないですね。細線化処理の中で、tan(x)の値は注目画素を含む3点のうちの強度が最大の点を求める際に使われているだけ。角度がthetaの時もtheta+PIの時も、比較に用いる3点の組み合わせは変わらないので、実はPIの範囲だけ考えればOKなのでは。tan(theta)=Math.Infiniteの時も、細線化処理の中でPI/2のケースにキチンと落ちるので問題なさそうに思います。

何よりも、個人的には、tan(theta)を使うと細線化処理の離散化における場合分けを綺麗に連続して取ることができる点が良いと感じます。クリティカルな問題がなければこちらにしたい。Math.atan2()の返り値は-PIからPIまでですから、1周全部考えなくてはいけないですし、-PIから-7PI/8、7PI/8からPIみたいな分断された場合ができてしまうのは可能であれば避けたいところです。見通しの良さは大事だと思うので。
https://docs.oracle.com/javase/jp/8/docs/api/java/lang/Math.html#atan2-double-double-

この話とは別に、下記をみると、Math.tan(double a)は「引数がNaNまたは無限大の場合、結果はNaN」とあります。引数チェックを忘れていることに気がついたので、修正しています。
https://docs.oracle.com/javase/jp/8/docs/api/java/lang/Math.html#tan-double-

この引数チェックに引っかかる事って、そうそうない気がするのですが、どういう時かを考えてみます。gxとgyは上位8ビットが同じと考えれば、大雑把なオーダーは同じ。しかも両方int。そうするとgy/gxがNaNになるのはgx=0のときでしょうか。それはつまり、水平方向に取るべきものがないとき。縦方向のエッジを強調する方向に倒しました。

ヒステリシス処理について、checkNeighbors()を再帰で呼ぶと良いのでは？というコメントを頂きました。ちょっと考えてみたのですが、注目画素の8近傍全部について再帰的にcheckNeighbors()を呼ぶという意味なら、スキャンの向きを変えているだけなのでは。さらに、ループの向きを変えて書く場合と違って、入力画像の画素数-1回分の再帰呼び出しがスタックに積まれる状態なのでは。例えば、縦に3ピクセル、横にNピクセル（N:任意の正の整数）な画像を、左から2ピクセル、上から2ピクセル目を注目画素として再帰的にcheckNeighbors()をかけていくと、ヒステリシス処理でみている3x3の範囲が、右端まで結果を保留しながら行って、右端から値を確定させながら帰ってきて、最後に注目画素をどうするか決まるイメージ。これだと、右端にエッジがあって、そこから左に向かって中間値が並ぶ場合は確かにエッジが引けるけど、逆に左端にエッジがあって、そこから右に向かって中間値が並んで、さらにヒステリシス処理の小さい方の閾値を下回る値が並ぶ場合にはエッジが一切引かれないのでは。スタックに再帰呼び出しが山ほど積まれるリスクの割には嬉しさが小さすぎる気がします。

2019/08/21
ソーベルフィルタの中でtanではなくatan2を使うべきというコメントを頂きました。ありがとうございます。定義に添えばその通り。ただタンジェントって図のような感じで（手抜きしたので軸系が情報ないですが、真ん中のグレーの線がtheta=0）、青が0、緑がPI/4、赤がPI/2、水色が3PI/4という扱い（のつもり）で書いていて、基本単位は満たしているし、角度とタンジェントの値は基本単位の中では1対1対応しているんじゃないのかな？という疑問が消えなかったので、そのままです。

一方で、この件を見直している時に、ソーベルフィルタに関する理解が微妙だったことに気が付いたので修正しました。具体的には、強度を別配列として持っていたのですが、それ画像だろってことに気がついて。道理で強度の値が馬鹿でかいはずだ…。だったらちゃんとRGBそれぞれで計算した方が良くない…？とか。ちょっとググってみたりしたけれども微妙にスッキリしなかったのも理由のひとつ。芋づる式に細線化の強度比較も書き直しました。勘違いバージョンよりは良さげな結果が出るようになりました。

ヒステリシスで、ラスタスキャンしてるから、右にエッジがあって左向きに中間値が並んでいる時には線が引かれなくなるよ、という指摘もいただきました。ありがとうございます。その通りかと。で、その問題を解決しようと思って考えてみると、スキャンの方向を変えたところで、同様の問題が向きを変えて起こるだけなんですよね。そこを解決するにはランダムにアクセスして誤魔化すとか、逆方向からもスキャンして合成するとか？ちょっとググった程度だとあまりその点を工夫されている方は見つからず。研究レベルではおられるのかもしれませんが。結論として「エッジって強くて長い線」だと思うので、中間値が並んでいるような状況では本来は線は引かれなくても良いのでは？という結論になって、保留したままです。