OpenCV日記（2）カスタム線形フィルタを画像に適用する

　線形フィルタを画像に適用して、適用結果を表示するプログラムを書いてみました。
以下のサンプルプログラムを参考にしました。

filter2D_demo.cpp

ソースコード

3 * 3の行列を標準入力から受け取って、その行列をカーネルとして画像に適用するようにしました。

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <cstdlib>
#include <iostream>

using namespace cv;

static const char* window_name = "filter2D Demo";

/**
 * 画像ファイルに3*3のカーネルを適用し、適用結果をウィンドウに表示する。
 */
int main (int argc, char** argv )
{
    /// Check command args
    if (argc != 2) {
        std::cout << "specify an image file path at the first argument." << std::endl;
        return -1;
    }

    /// Load an image
    Mat src = imread(argv[1]);
    if (!src.data) { 
        return -1; 
    }

    /// Initialize arguments for the filter
    Point anchor = Point(-1, -1);
    double delta = 0;
    int ddepth = -1;
    
    /// Define a kernel matrix
    double kernel[3][3];
    std::cout << "input a 3*3 filter matrix." << std::endl;
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            std::cin >> kernel[i][j];
    
    /// Apply filter
    Mat dst;
    filter2D(src, dst, ddepth, Mat(3, 3, CV_64F, &kernel), anchor, delta, BORDER_DEFAULT);
    
    /// Show the result
    imshow(window_name, dst);
    waitKey(0);

    return 0;
}

結果

lena.jpgにいくつかのカーネルフィルタを適用してみました。

まず、Original filter。Identity functionです。

0 0 0

0 1 0

0 0 0

次に、Edge-Detect filter。８近傍に対して微分しています。

-1 -1 -1

-1  8  -1

-1 -1 -1 

Blur filter。ぼかしです。ぼかし効果をあげたい場合は、行列のサイズを大きくするとよいです。

1/9 1/9 1/9

1/9 1/9 1/9

1/9 1/9 1/9

Contrast filter（このフィルタ名は勝手につけました）。セルの輝度を定数倍してみます。

0 0 0

0 2 0

0 0 0

OpenCV日記（1）インストール

　とりあえずインストールして動かしてみた。

インストール

Ubuntu Community Helpにインストールの説明があったので、それを参考にインストールを実施（といってもシェルスクリプトを実行するだけ）。
最新版(ver 2.4.8)のインストールスクリプトはここのを使用。

動作確認

同じくUbuntu Community Helpにサンプルソース、コンパイルコマンドが載っていたので、そのとおりに実行。

動いた！

コンパイルコマンドを毎回打つのは面倒なので、シェルスクリプトにしてopencvというaliasを設定した。

サンプルソースをチラ見

opencv-2.4.8/samples/cpp にエッジ検出とか、SVMとか、K-meanとか、その他いろいろのサンプルソースがある。なかなか面白そうだ。

クラスローダーからリソースを取得する

　JavaでOpenCVを動かして遊んでいたら、

getClass().getResource("ファイル名").getPath()

のようにしてファイルのパスを取得しているサンプルがあった。これ何だろうなーと思って調べてみた。

Mat image = Highgui.imread("lena.png");

でいいんじゃないの？と思った。しかし、上のように直接ファイルパスを書いてしまうと、このファイルパスはカレントディレクトリに依存してしまう。
よって、想定外のディレクトリからプログラムを実行する(※)と、ファイルを見つけることが出来なくなってしまう。

※一個上のディレクトリ階層からプログラムを実行する場合など。より実用的な例としては、jarにまとめて実行したり、サーブレット上でプログラムを実行したりする場合。

これに対して、

Mat image = Highgui.imread(getClass().getResource("/lena.png").getPath());

のようにしておけば、クラスパスからリソースが検索されるため、カレントディレクトリには影響しない。

std::swapで配列をスワップする

 C++のstd::swapを使うと

std::swap(a, b)

のようにして、変数aとbの値を入れ替えることが出来る。

これは配列に対しても適用出来る。例えば、以下のようにして配列a[]と配列b[]をスワップ出来る。

#include <iostream>

using namespace std;

void dump(int arr[], int sz) {
    for (int i = 0; i < sz; i++)
        cout << arr[i] << " ";
    cout << endl;
}

int main(int argc, char **argv) {
    int a[] = {1,2,3};
    int b[] = {4,5,6};

    swap(a, b);

    dump(a, 3);   // 4 5 6
    dump(b, 3);   // 1 2 3
    
    return 0;
}

また、行列（二次元配列）の行の入れ替え操作を行うことも出来る。例えば、以下のようにして行列aの1行目と3行目を入れ替えることが出来る。

#include <iostream>

using namespace std;

void dump(int arr[][3], int sz) {
    for (int i = 0; i < sz; i++) {
        for (int j = 0; j < 3; j++) {
            cout << arr[i][j] << " ";
        }
        cout << endl;
    }
}

int main(int argc, char **argv) {
    
    int a[][3] = {
        {1,2,3},
        {4,5,6},
        {7,8,9}
    };

    swap(a[0], a[2]);

    dump(a, 3);
    // 7 8 9 
    // 4 5 6
    // 1 2 3
    
    return 0;
}

なぜこのようなことが出来るのか納得できなかったので、swap関数の中身を見てみた。どうやら、配列に対してオーバーロードされているようだ。gccのstd::swapの定義は以下のような感じだった。

template<typename _Tp, size_t _Nm>
inline void swap(_Tp (&__a)[_Nm], _Tp (&__b)[_Nm])
{
    for (size_t __n = 0; __n < _Nm; ++__n)
        swap(__a[__n], __b[__n]);
}

ほうほう、テンプレート型のところに配列のサイズを指定できるのか。
何かテキトーな関数を書いて実験してみた。

#include <iostream>

using namespace std;

template<typename T>
inline void dump(const T &a) {
    cout << a << " ";
}

template <typename T, size_t N>
inline void dump(T const (&a)[N]) {
    for (size_t i = 0; i < N; i++)
        dump(a[i]);
    cout << endl;
}

int main(int argc, char **argv) {

    int x = 10;
    
    int a[] = {1,2,3,4,5,6,7,8,9,10};
    
    int b[][3] = {
        {1,2,3},
        {4,5,6},
        {7,8,9}
    };


    // x = 10 
    cout << "x = ";
    dump(x);
    cout << endl;

    // a = 1 2 3 4 5 6 7 8 9 10 
    cout << "a = ";
    dump(a);
    cout << endl;
    
    //b = 1 2 3 
    //4 5 6 
    //7 8 9 
    cout << "b = ";
    dump(b);
    cout << endl;
    
    return 0;
}

予想どおりに動いた。テンプレートはこんなことが出来るのかー。面白い。。

コルセルトの判定法を使ったカーマイケル数の検出

　コルセルトの判定法を使って[0, 1000000)内のカーマイケル数を列挙してみました。

コルセルトの判定法とは、以下のような条件です。

正の合成数nがカーマイケル数であるための必要十分条件は以下のとおりである。
・nはsquare-free（平方数で割り切れない）である。
・nのすべての素因数pについて、p-1 | n-1が成り立つ。

与えられた数の素因数分解が出来ればいいので、O(√n)で判定可能です。

　[0, 1000000)までの各整数に対してコルセルトの条件が成り立つかどうかを調べると、全体でO(n√n)かかってしまいます。少し工夫してエラトステネスの篩を応用すると、O(n log log n)でカーマイケル数の列挙が出来ます。

　以下ソースコードです。

#include <iostream>
#include <cassert>
#include <map>

using namespace std;

bool carmichael[1000000];
bool prime[1000000];

#ifdef debug
map<long long, int> factorize(long long n) {
    map<long long, int> ret;
    
    for (long long i = 2; i * i <= n; i++) {
        while (n % i == 0) {
            ++ret[i];
            n /= i;
        }
    }
    
    if (n != 1)
        ret[n] = 1;
    
    return ret;
}

void korseltCriterion(int x) {
    map<long long, int> fact = factorize(x);

    for (const auto &pr : fact) {
        assert(pr.second == 1);
        assert((x-1) % (pr.first-1) == 0);
    }
}
#endif

int main(int argc, char **argv) {
    fill(carmichael, carmichael+1000000, true);
    fill(prime, prime+1000000, true);

    carmichael[0] = carmichael[1] = false;
    prime[0] = prime[1] = false;

    for (int i = 2; i < 1000000; i++) {
        if (!prime[i])
            continue;

        carmichael[i] = false;   // i is a prime number

        long long sq = (long long)i * i;
        for (int j = 2 * i; j < 1000000; j += i) {
            prime[j] = false;
            if (j % sq == 0)
                carmichael[j] = false;  // j is not square-free
            if ((j - 1) % (i - 1))
                carmichael[j] = false;  // (i-1) | (j-1) does not hold
        }
    }
    
    for (int i = 0; i < 1000000; i++) {
        if (carmichael[i]) {
            cout << i << endl;
#ifdef debug
            korseltCriterion(i);
#endif
        }
    }

    return 0;
}

以下、[0, 1000000)のカーマイケル数のリストです。
561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361
101101
115921
126217
162401
172081
188461
252601
278545
294409
314821
334153
340561
399001
410041
449065
488881
512461
530881
552721
656601
658801
670033
748657
825265
838201
852841
997633

Solovay–Strassenの素数判定法

Solovay–Strassenの素数判定法をC++で書いて見ました。

Solovay–Strassenの素数判定法とは？

与えられた数が素数かどうかを判定する乱択アルゴリズムです。RSA暗号の実現可能性を示すために重要な役割を担った歴史的価値の高いアルゴリズムです[1]。

Solovay–Strassenの素数判定法は、オイラーの基準

a^(p-1)/2 = (a / p)  (mod p)  (1)

を利用します。ただし、pは素数です。(a / p)はルジャンドルの記号で、その値は

0　 a = 0 (mod p)の場合
1     aが法pにおける平方剰余の場合
-1　aが法pにおける平方剰余ではない場合

のように定義されます。

与えられたnが式(1)を満たせば素数、満たさなければ合成数と判定します。
左辺は繰り返し二乗法を使えば高速（nを多倍長整数とするとnの桁数に比例する時間）に計算できます。右辺は平方剰余の相互法則[2]を使えば高速に計算できます。

フェルマーテストのときと同様に、この判定法では疑陽性が発生します。
ただしaをランダムにt回選んで判定を行うことで、疑陽性の発生確率は (1/2)^t 以下になることが知られています。またこの判定法では、フェルマーテストのときのカーマイケル数のように、疑陽性の発生確率が大きくなるような特別な数はありません。

C++で書いたスパゲッティなソース

以下ソースコードです。

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;

typedef boost::multiprecision::cpp_int bigint;

bigint mpow(bigint x, bigint p, bigint mod) {
    bigint ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

bigint jacobi(bigint a, bigint n) {
    if (a == 0) return 0;
    if (a == 1) return 1;

    if ((a & 1) == 0) {
        if (n % 8 == 1 || n % 8 == 7)
            return jacobi(a>>1, n);
        else
            return -jacobi(a>>1, n);
    }

    if (a >= n)
        return jacobi(a % n, n);
    
    if (a % 4 == 3 && n % 4 == 3)
        return -jacobi(n, a);
    else
        return jacobi(n, a);
}

bool SolovayStrassen(bigint n, int t = 30) {
    if (n <= 2) 
        return n == 2;

    for (int i = 0; i < t; i++) {
        bigint a = rand() % n;
        if (mpow(a, (n-1)/2, n) != (jacobi(a, n) + n) % n)
            return false;
    }

    return true;
}

int main(int argc, char **argv) {
    srand((unsigned)time(NULL));

    bigint n;
    cin >> n;
    cout << SolovayStrassen(n) << endl;

    return 0;
}

参考URL

[1] http://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test
[2] http://math.fau.edu/richman/jacobi.htm

フェルマーテストで素数判定

　ちょっと気晴らしにフェルマーテストを書いてみました。

この前boostの多倍長ライブラリを入れたので、C++で書いてみました。テストの式の微妙な違いで大きな差が出たのでちょっとメモしておきます。

a^p = a (mod p)    (1)

のようなテストをすると、大量のカーマイケル数が陽性になりました。ちなみに(1)は正しいフェルマーテストのチェック式ではなくて、何も見ずにプログラムを書いてみたときに使った式です。

で疑陽性多すぎるなと思って、よくよく調べてみると正しい式は(2)でした。

 a^p-1 = 1 (mod p)    (2)

上式でテストすると、疑陽性は少なくなりました。[1, 1000000)まででテストしましたが、疑陽性は出ませんでした（各数に対するテスト反復回数=100）。(1)と(2)って同じように見えるけど、(2)の方が厳しい条件ですね。

例えばカーマイケル数561は、
237⁵⁶¹ = 237 (mod 561) のようになります。
しかし、 237⁵⁶⁰ = 375 (mod 561)
なので、(1)のようなテストをすると疑陽性が発生し、(2)のようなテストをすると正しく陰性になるというわけです。

んー、ちゃんと定義は見た方がいいな。ということで、カーマイケル数の定義を確認してみました。

wikipediaを見てみると、
合成数 n がカーマイケル数であるとは、自身と互いに素である任意の自然数 a に対し、
a^n-1 = 1 (mod n)    (3)
を満たすことをいう。

とあります。

カーマイケル数は互いに素でないaに対しては(3)を満たすとは限らないということか。。

確かに、例で見た237と561はどちらも3で割れるので、互いに素ではないです。

んーということは、カーマイケル数nに対するフェルマーテストを行う際に、nと互いに素ではないaがうまく選ばれれば疑陽性の発生頻度が少なくなる可能性があると考えてもよさそうです。

Fermat witnessの濃度を考えると実験で設定した反復回数=100ってちょっと多すぎやろと思いましたが、カーマイケル数に対するテストでカーマイケル数と互いに素ではない数がたくさん選ばれるようにするという意味で100という数はなかなかいいチョイスだったのかもしれません。。

今回のテストで使った範囲[0, 1000000)内のカーマイケル数のtotient functionを計算してみると、
Φ(561) = 320
Φ(1105) = 768

Φ(1729) = 1296

Φ(2465) = 1792
......
Φ(488881) = 466560
......

のように。

488881の場合は、coprimeじゃない数は5%くらいの確率でしか出ない。んー、やはり100という設定値は絶妙なチョイスだったようだ。

と、意外と発見があったフェルマーテストの実装でした。
コードを載せておきます。

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;

typedef boost::multiprecision::cpp_int bigint;

bigint pow(bigint x, bigint p, bigint mod) {
    bigint ret = 1;

    while (p > 0) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }

    return ret;
}

bool fermatTest(bigint p, int k = 100) {
    if (p < 2)
        return false;
    
    srand((unsigned)time(NULL));
    for (int i = 0; i < k; i++) {
        bigint a = rand() % (p-1) + 1;
        if (pow(a, p-1, p) != 1)  // this condition is better than pow(a, p, p) != a
            return false;
    }

    return true;
}

#ifndef TEST
int main() {
    bigint n;
    cin >> n;

    cout << fermatTest(n) << endl;
    
    return 0;
}
#else

bigint primalityTest(bigint x) {
    for (bigint i = 2; i * i <= x; i++) {
        if (x % i == 0)
            return false;
    }

    return x > 1;
}

int main() {
    bool res = true;
    for (int i = 1; i < 1000000; i++) {
        if (fermatTest(i) != primalityTest(i)) {
            res = false;
            cout << i << endl;
        }
    }
    if (res)
        cout << "passed." << endl;

    return 0;
}
#endif