上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
以前、Matlab(Octave)で実装した。
位相限定相関法(POC)




位相画像をOpenCVSharpで実装した
位相画像をC#とOpenCVSharpで実装

の続きです。

ソース

ソースはクラスもきちんと作れていないし、汚いソースなので紙面の都合上、肝の部分だけ抜粋します。

あとバグがありそうですので・・・見つけたらオシエテクダサイ。

// DFT用の最適サイズを計算し,そのサイズを返却する。
public static int GetOptimalDFTSize(int size)
{
    return Cv.GetOptimalDFTSize(size);
}

public static CvMat dft(IplImage img)
{
    IplImage reimg = Cv.CreateImage(img.Size, BitDepth.F64, 1);
    IplImage imimg = Cv.CreateImage(img.Size, BitDepth.F64, 1);
    IplImage cmpimg = Cv.CreateImage(img.Size, BitDepth.F64, 2);

    // (1)入力画像を実数配列にコピーし,虚数配列とマージして複素数平面を構成
    Cv.Scale(img, reimg, 1.0, 0.0);
    Cv.Zero(imimg);
    Cv.Merge(reimg, imimg, null, null, cmpimg);

    // (2)DFT用の最適サイズを計算し,そのサイズで行列を確保する
    int dft_M = GetOptimalDFTSize(img.Height);
    int dft_N = GetOptimalDFTSize(img.Width);

    CvMat dft_A = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C2);

    // (3)複素数平面をdft_Aにコピーし,残りの行列右側部分を0で埋める
    CvMat tmp;
    Cv.GetSubRect(dft_A, out tmp, new CvRect(0, 0, img.Width, img.Height));
    Cv.Copy(cmpimg, tmp, null);
    if (dft_A.Cols > img.Width)
    {
        Cv.GetSubRect(dft_A, out tmp, new CvRect(img.Width, 0, dft_A.Cols - img.Width, img.Height));
        Cv.Zero(tmp);
    }

    // (4)離散フーリエ変換を行う
    Cv.DFT(dft_A, dft_A, DFTFlag.Forward, cmpimg.Height);
    //Cv.Split(dft_A, image_Re, image_Im, null, null);

    return dft_A;
}

public static IplImage idft(CvMat dft)
{
    // DFT用の最適サイズを計算し,そのサイズで行列を確保する
    int dft_M = GetOptimalDFTSize(dft.Height);
    int dft_N = GetOptimalDFTSize(dft.Width);
    IplImage img_Re = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1);
    IplImage img_Im = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1);

    // 逆変換
    Cv.DFT(dft, dft, DFTFlag.Inverse, dft_N);

    // 逆変換した画像を表示用の変数に格納
    Cv.Split(dft, img_Re, img_Im, null, null);

    return img_Re;
}

public static CvMat Conj(CvMat dft)
{
    int dft_M = Image.GetOptimalDFTSize(dft.Height);
    int dft_N = Image.GetOptimalDFTSize(dft.Width);

    using (CvMat srcRe = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C1))
    using (CvMat srcIm = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C1))
    {
        // 実部と虚部に分解
        Cv.Split(dft, srcRe, srcIm, null, null);
        for (int i = 0; i < dft_M; i++)
        {
            for (int j = 0; j < dft_N; j++)
            {
                srcIm[i, j] = -1 * srcIm[i, j];
            }
        }
        // 実部・虚部で計算した値をマージする
        Cv.Merge(srcRe, srcIm, null, null, dft);
    }
    return dft;
}

public static IplImage PhaseOnlyCorrelation(IplImage img1, IplImage img2)
{
    int dft_M = Image.GetOptimalDFTSize(img1.Height);
    int dft_N = Image.GetOptimalDFTSize(img1.Width);

    IplImage img = Cv.CreateImage(new CvSize(dft_M, dft_N), BitDepth.F64, 1);
    CvMat dft1 = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C2);
    CvMat dft2 = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C2);
    CvMat result = Cv.CreateMat(dft_M, dft_N, MatrixType.F64C2);

    // フーリエ変換
    dft1 = Image.dft(img1);
    // フーリエ変換後に複素共役にする
    dft2 = Image.Conj(Image.dft(img2));

    using (IplImage srcRe1 = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    using (IplImage srcIm1 = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    using (IplImage srcRe2 = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    using (IplImage srcIm2 = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    using (IplImage dstRe = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    using (IplImage dstIm = Cv.CreateImage(new CvSize(dft_N, dft_M), BitDepth.F64, 1))
    {
        // フーリエ変換した画像を実部と虚部に分解
        Cv.Split(dft1, srcRe1, srcIm1, null, null);
        Cv.Split(dft2, srcRe2, srcIm2, null, null);

        for (int i = 0; i < dft_M; i++)
        {
            for (int j = 0; j < dft_N; j++)
            {
                // 画像をかけあわせる F×G*
                dstRe[i, j] = srcRe1[i, j] * srcRe2[i, j] - srcIm1[i, j] * srcIm2[i, j];
                dstIm[i, j] = srcRe1[i, j] * srcIm2[i, j] + srcRe2[i, j] * srcIm1[i, j];

                // 絶対値を計算しその値で割る |F×G*|
                double spectrum = Math.Sqrt(dstRe[i, j] * dstRe[i, j] + dstIm[i, j] * dstIm[i, j]);
                dstRe[i, j] = dstRe[i, j] / spectrum;
                dstIm[i, j] = dstIm[i, j] / spectrum;
            }
        }
        // 計算結果を統合する
        Cv.Merge(dstRe, dstIm, null, null, result);
    }

    // 逆フーリエ変換
    img = idft(result);

    // データを解放
    dft1.Dispose();
    dft2.Dispose();
    result.Dispose();

    return img;
}



しかし、OpenCVって複素数をサポートしていないのかなぁ・・・
チャンネル毎に行うのはちょっと面倒で計算ミスが起こりそう・・・


結果

今回も画像はnashruddin.com - Phase Correlation in OpenCVさんから拝借したもので、行いました。


結果は以下のようになるようです。

相関値:0.6344
移動距離:(6,15)



1.jpg 2.jpg

POCを行ったグラフ
左上に白い点がありますが、そこがピーク値です。
POC1.png


ちなみに、私の計算結果は・・・
相関値:0.5904
移動距離:(6,15)



・・・あれ?ちょっと相関値が違う・・・
移動させた距離はあっているのだけどなぁ・・・


関連記事
スポンサーサイト
コメント
ありがとうございます!
とても参考になりました!ただ、以下の点が気になりました
dstRe[i, j] = srcRe1[i, j] * srcRe2[i, j] - srcIm1[i, j] * srcIm2[i, j];
dstIm[i, j] = srcRe1[i, j] * srcIm2[i, j] + srcRe2[i, j] * srcIm1[i, j];

dstRe[i, j] = srcRe1[i, j] * srcRe2[i, j] + srcIm1[i, j] * srcIm2[i, j];
dstIm[i, j] = srcRe1[i, j] * srcIm2[i, j] - srcRe2[i, j] * srcIm1[i, j];
ではないでしょうか?勘違いだったらもうしわけありません。
Re: ありがとうございます!
>かず様
コメントいただきありがとうございます。

画像2に対して複素共役をとったものを画像1と
掛け合わせるため上記のようになるのではないかということでしょうか。

コード内の94行目にて以下のように複素共役にしております。
dft2 = Image.Conj(Image.dft(img2));

そのため、該当箇所は通常の複素数の積として扱ってよいのではないかと考えました。
いかがでしょうか。
もし、考え違いなどございましたら、ご指摘お願いします。
勘違いでした。
申し訳ありません、十分に確認していませんでした。
丁寧な解説ありがとうございます、複素共役できちんととっておられたのですね。




Re: 勘違いでした。
>かず様
こちらこそ、ご指摘ありがとうございます。
コード内のコメントや記載方法がわかりにくいため誤解を与えてしまいました。

若干バグが潜んでいそうな気がしているため、
また何かありましたらよろしくお願いいたしますm(_ _)m
OpenCVで・・・
こんにちは.
POCについて詳しく解説してくださり,大変参考になりました.

1つ質問があるのですが,OpenCVで実装されているPOCは返り値がcv::Point2dなのですが,管理人様が実装したように,OpenCVで相関値を求めることはできるのでしょうか?

お時間ある時にご回答頂けると幸いです.
Re: OpenCVで・・・
コメントいただきありがとうございます。
そして返信が遅くなり申し訳ありません。

modules\imgproc\src\phasecorr.cpp内の関数を見ましたが、
どうやら相関値はなさそうです。。。

これは予想なのですが、
そもそもPOCの場合の相関値は通常はデルタ相関に近い形となるようです。
(以前、数式で証明を見た気がいたします)
そのため、そもそも複数のピークを想定しなくてもよいという理由のためではないでしょうか。

あとは、C言語の複数のreturn値を持てないため切り捨てられてしまったなどです。。。
もし、相関値を必要とされるのでしたら、ライブラリを修正してビルドするか自前の関数内に同様の実装を実現する方がよいかもしれません。。。

お答えになっていれば幸いです。
コメントの投稿
Contributor
Comment
トラックバック URL
URL
このエントリーの固定リンク
トラックバック
FC2カウンター
プロフィール

Author:詠み人知らず
プログラム好きな名もなき凡人がお送りしています。(得意とは言っていない
最近の興味はPython、C#、Matlab(Octave)、画像処理、AR(拡張現実)、統計などなど・・・

気分で思いついたことを書くため話題に一貫性がないかもしれません。

カレンダー
08 | 2017/03 | 09
- - - 1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31 -
最新記事
タグクラウドとサーチ

カテゴリ
最新コメント
最新トラックバック
月別アーカイブ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。