よく位相限定相関(POC)をキーワードでブログにいらっしゃる方がいるのと、
自分のためにまとめておきます。


理論
位相限定相関法についての求め方を書いてあります。

また、4ではOpenCVで実装されているものの紹介。
5では従来の手法では検出できない回転に対して検出するための前処理であるLog-Polar変換を書いております。

  1. 位相画像
  2. 位相限定相関法(POC)
  3. 位相限定相関(POC)でサブピクセル精度を求める方法
  4. 位相限定相関(POC)がOpenCVで実装されている
  5. Log-Polar変換を画像に適用



OpenCvSharp
OpenCvSharpで実装してみたものです。。。があまり自信はありません。
  1. 位相画像をC#とOpenCVSharpで実装
  2. 位相限定相関法(POC)をOpenCVSharpで



画像補間について
回転がされている場合でも検出できるようにLog-Polar変換の時に必要であったため記載したものです。
ですが、有名な補間手法なので画像の回転・拡大・縮小に使えます。

  1. 画像を回転させる際の画素を補間するコード
  2. 回転した画像をLanczos補間する
  3. 回転した画像をバイキュービック補間する
  4. 回転した画像をバイリニア補間する
  5. 画像を回転させる



何か進展があればまた追記していきます。


以前、書いた投稿でSURFで二枚の画像での対応付けを行いました。

PythonからOpenCVのSURFを使う | 詠み人知らずの備忘録



SURFで抽出した対応する座標を用いれば、画像同士を合成することができます。
以下に参考にしたサイトです。


2枚の画像のモザイキング | Miyabiarts.net
前に紹介した「OpenCVを用いた特徴点対応付け」を応用することで複数の画像をモザイキングすることができます。 今回は、OpenCVを用いて下に示す2枚の画像をモザイキングして、1枚の大きな画像を作ります。 今回の例では、2枚の画像が大体半分程度重なり合ってないと上手くいかないです。 あと、画像合成の関係上、1枚目の画像に対して2枚目の画像が左側になるように撮影しないと正しい結果が得られませんので気をつけてください。



Panorama – Image Stitching in OpenCV | ramsrigoutham
he code snippet shown below is for simple image stitching of two images in OpenCV . It can easily be modified to stitch multiple images together and create a Panorama. OpenCV also has a stitching module which helps in achieving this task and which is more robust than this. The code presented here will help in understanding the major steps involved in image stitching algorithm. I am using OpenCV 2.4.3 and Visual studio 2010. This code is based on the openCV tutorial available here.



SURFを実行しても、surf結果画像のように完璧に対応点を抽出できません。 いくつか、誤検出が発生してしまいます。


surf結果画像
surf.png


そのため、誤検出を取り除く処理が必要になります。
RANSACが有名な方法のようです。

画像を合成する際に移動および回転などをさせるために射影変換を行いますが、OpenCVのfindhomographyでは、CV_RANSACを第3引数にすることにより、対応点をRANSACでふるいにかけてくれるようです。

ソース

途中まではSURFを求めるまでは前回までと同じソースを使用しています。
その座標を用いてHomography行列を求め合成しています。

# -*- coding: shift_jis -*-
# matching features of two images
import cv2
import sys
import scipy as sp
import numpy

if len(sys.argv) < 3:
    print 'usage: %s img1 img2' % sys.argv[0]
    sys.exit(1)

img1_path = sys.argv[1]
img2_path = sys.argv[2]

img1 = cv2.imread(img1_path, cv2.CV_LOAD_IMAGE_GRAYSCALE)
img2 = cv2.imread(img2_path, cv2.CV_LOAD_IMAGE_GRAYSCALE)

detector = cv2.FeatureDetector_create("SURF")
descriptor = cv2.DescriptorExtractor_create("BRIEF")
matcher = cv2.DescriptorMatcher_create("BruteForce-Hamming")

# detect keypoints
kp1 = detector.detect(img1)
kp2 = detector.detect(img2)

print '#keypoints in image1: %d, image2: %d' % (len(kp1), len(kp2))

# descriptors
k1, d1 = descriptor.compute(img1, kp1)
k2, d2 = descriptor.compute(img2, kp2)

print '#keypoints in image1: %d, image2: %d' % (len(d1), len(d2))

# match the keypoints
matches = matcher.match(d1, d2)

# visualize the matches
print '#matches:', len(matches)
dist = [m.distance for m in matches]

print 'distance: min: %.3f' % min(dist)
print 'distance: mean: %.3f' % (sum(dist) / len(dist))
print 'distance: max: %.3f' % max(dist)

# threshold: half the mean
thres_dist = (sum(dist) / len(dist)) * 0.9

# keep only the reasonable matches
sel_matches = [m for m in matches if m.distance < thres_dist]

print '#selected matches:', len(sel_matches)

point1 = [[k1[m.queryIdx].pt[0], k1[m.queryIdx].pt[1]] for m in sel_matches]
point2 = [[k2[m.trainIdx].pt[0], k2[m.trainIdx].pt[1]] for m in sel_matches]

point1 = numpy.array(point1)
point2 = numpy.array(point2)

H, Hstatus = cv2.findHomography(point2,point1,cv2.RANSAC)

# 移動量を算出
x=0
y=0
cnt=0
for i,v in enumerate(Hstatus):
    if v==1:
        x += point1[i][0]-point2[i][0]
        y += point1[i][1]-point2[i][1]
        cnt += 1

# カラー画像として改めて読み込む
img1 = cv2.imread(img1_path)
img2 = cv2.imread(img2_path)

x = abs(int(round(x/cnt)))
y = abs(int(round(y/cnt)))

# sizeを取得
h1, w1 = img1.shape[:2]
h2, w2 = img2.shape[:2]

dst = cv2.warpPerspective(img2,H,(w2+x,h2+y))

for i in xrange(w1):
    for j in xrange(h1):
        dst[j,i] = img1[j,i]

cv2.imshow("result", dst)
cv2.waitKey()



結果

入力画像
wall1.jpg


wall2.jpg


2枚目画像のHomography行列と演算結果
homography1.png


結果画像
homography2.png


2枚目画像のHomography行列と演算結果を見て分かるかもしれませんが、気づいたことを何点か・・・

このプログラムでは、1枚目の画像に合わせるように2枚目の画像を射影変換して合わせているため、 画像の左右の対応関係が逆だったりすると合成がうまくいかなかったり。
1枚目の画像より上部にある2枚目の画像が反映されていなかったりとあります。
#それは、理論というよりもプログラム上の問題かもしれませんが

ただ、複数枚つなげていくと誤差が積み重なっていきそうではあります。



お久しぶりです。
今回は普通(?)のOpencvとPythonです。
#OpenCVSharpとの兼ね合いのため、私の環境のOpenCVのバージョンが2.4.0となっておりますので、ご注意願いますm(_ _)m



導入

OpenCVはインストールされていると仮定してですが、
Pythonで動かすにはもう一手間必要です。

C:\opencv\build\python\2.7 にある
cv2.pyd というファイルを下記のフォルダにコピーします。

C:\Python27\Lib\site-packages
#フォルダ名などの違いは各自のフォルダに読み替えて実行願います。



ソース

SURFの実行方法を調べていたところ、下記サイトで相談&回答されています。


How to visualize descriptor matching using opencv module in python - Stack Overflow

ただ、このまま実行したところ、エラーが生じたため。その部分(65行目)を修正しました。
# matching features of two images
import cv2
import sys
import scipy as sp

if len(sys.argv) < 3:
    print 'usage: %s img1 img2' % sys.argv[0]
    sys.exit(1)

img1_path = sys.argv[1]
img2_path = sys.argv[2]

img1 = cv2.imread(img1_path, cv2.CV_LOAD_IMAGE_GRAYSCALE)
img2 = cv2.imread(img2_path, cv2.CV_LOAD_IMAGE_GRAYSCALE)

detector = cv2.FeatureDetector_create("SURF")
descriptor = cv2.DescriptorExtractor_create("BRIEF")
matcher = cv2.DescriptorMatcher_create("BruteForce-Hamming")

# detect keypoints
kp1 = detector.detect(img1)
kp2 = detector.detect(img2)

print '#keypoints in image1: %d, image2: %d' % (len(kp1), len(kp2))

# descriptors
k1, d1 = descriptor.compute(img1, kp1)
k2, d2 = descriptor.compute(img2, kp2)

print '#keypoints in image1: %d, image2: %d' % (len(d1), len(d2))

# match the keypoints
matches = matcher.match(d1, d2)

# visualize the matches
print '#matches:', len(matches)
dist = [m.distance for m in matches]

print 'distance: min: %.3f' % min(dist)
print 'distance: mean: %.3f' % (sum(dist) / len(dist))
print 'distance: max: %.3f' % max(dist)

# threshold: half the mean
thres_dist = (sum(dist) / len(dist)) * 0.5

# keep only the reasonable matches
sel_matches = [m for m in matches if m.distance < thres_dist]

print '#selected matches:', len(sel_matches)

# #####################################
# visualization
h1, w1 = img1.shape[:2]
h2, w2 = img2.shape[:2]
view = sp.zeros((max(h1, h2), w1 + w2, 3), sp.uint8)
view[:h1, :w1, 0] = img1
view[:h2, w1:, 0] = img2
view[:, :, 1] = view[:, :, 0]
view[:, :, 2] = view[:, :, 0]

for m in sel_matches:
    # draw the keypoints
    # print m.queryIdx, m.trainIdx, m.distance
    color = tuple([sp.random.randint(0, 255) for _ in xrange(3)])
    cv2.line(view, (int(k1[m.queryIdx].pt[0]),int(k1[m.queryIdx].pt[1])), (int(k2[m.trainIdx].pt[0] + w1), int(k2[m.trainIdx].pt[1])), color)

cv2.imshow("view", view)
cv2.waitKey()



実行結果

結果は以下のようになります。

result.png


無事結果が表示されました ε-(´∀`*)ホッ

Log-Polar変換なのですが、霊長類の網膜のモデルらしいです。
中心部は解像度が高いのですが、周辺部は低いという特性らしいです。

以下のような式で変換をします。(極座標変換でのrにlogをとったものです。)
logpolar.png


ですが、今まで行った回転のときのように、どの変換先の座標は元画像のどの位置を参照するかと考えます。
すなわち以下の逆変換を使います。
logpolar2.png



ソース


OpenCVのやりかたを真似したため、一部のみ抜粋です。
bmpが入力画像、lpcが出力画像です。
また、Mはformから入力された値を使用しました。256x256だと40くらいがよさそうです。
BiLinearInterpolationは以前の補間のさいに使用したものです。

画像を回転させる際の画素を補間するコード | 詠み人知らずの備忘録


Bitmap lpc = new Bitmap(bmp.Width, bmp.Height);

int cx = lpc.Width / 2;
int cy = lpc.Height / 2;

double M = Convert.ToDouble(textBox1.Text);

for (int i = 0; i < lpc.Width; i++)
{
    for (int j = 0; j < lpc.Height; j++)
    {
        double r = Math.Exp((double)i / (double)M);
        double x = r * Math.Cos(2 * Math.PI * j / lpc.Height) + cx;
        double y = r * Math.Sin(2 * Math.PI * j / lpc.Height) + cy;
        
        if (0 < x && x < bmp.Width - 1 && 0 < y && y < bmp.Height - 1)
        {
            Color bmpCol = BiLinearInterpolation(x, y, bmp);
            lpc.SetPixel(i, j, bmpCol);
        }

    }
}



結果

変換した結果、(x,y)と(rho,phi)の関係がどのような対応になるのかは、下記のサイトがイメージがつきやすいと思います。
Log polar transform - Rhea



※クリックで拡大します。

【元画像】
Lenna

【変換画像】
Lenna_log.png

【元画像】
Lenna_90.png

【変換画像】
Lenna_log90.png



画像の回転が、変換後の画像の下への移動としてあらわされています。
この考え方を、位相限定相関法(POC)と供に使用すると、角度のズレを判別できるためより便利な手法となります。(RIPOCという手法らしいです。)

位相限定相関法(POC) - スズメレンダラー・クマ将棋の開発日記

以前、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)



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


FC2カウンター
プロフィール

詠み人知らず

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

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

カレンダー
03 | 2017/04 | 05
- - - - - - 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 - - - - - -
最新記事
タグクラウドとサーチ

カテゴリ
最新コメント
最新トラックバック
月別アーカイブ