OpenCV DFT
(OpenCV) DFT Discite Fourier Transform
離散フーリエ変換
cv::dft, cvDFTを使った2次元フーリエ。簡単に出来るかと思ったら死ぬほど面倒だった。
※OpenCV2.xにはsamples/c/dft.cがあり、そっちの方が分かり易い。
OpenCV 窓関数も参照のこと。
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <stdio.h>
using namespace cv;
using namespace std;
static void help()
{
printf("\nThis program demonstrated the use of the discrete Fourier transform (dft)\n"
"The dft of an image is taken and it's power spectrum is displayed.\n"
"Usage:\n"
"./dft [image_name -- default lena.jpg)]\n");
}
const char* keys =
{
"{1| |lena.jpg)|input image file}"
};
int main(int argc, const char ** argv)
{
help();
CommandLineParser parser(argc, argv, keys);
string filename = parser.get<string>("1");
//画像をグレースケールで読み込む+エラー処理
Mat img = imread(filename.c_str(), CV_LOAD_IMAGE_GRAYSCALE);
if( img.empty() )
{
help();
printf("Cannot read image file: %s\n", filename.c_str());
return -1;
}
//DFTに最適なサイズを取得
int M = getOptimalDFTSize( img.rows );
int N = getOptimalDFTSize( img.cols );
//最適なサイズはimgより大きいので周辺を0で埋めたMatを作成
Mat padded;
copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));
//複素数画像(実部チャネルと虚部チャネルの2チャネル)の準備
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
//2枚のMatをmergeしてcomplexImgへ収める
Mat complexImg;
merge(planes, 2, complexImg);
//DFT処理
dft(complexImg, complexImg);
// DFTの結果は値の幅が大きいのでlogにするが，値が0だとエラーになるので全体に1を加えてから
// DFTの結果は√((実部)^2+(虚部)^2) とする
// compute log(1 + sqrt(Re(DFT(img))**2 + Im(DFT(img))**2))
split(complexImg, planes);
//要素を2乗したものをplanes[0]に上書きで収める
magnitude(planes[0], planes[1], planes[0]);
Mat mag = planes[0];
// 1足してlog
mag += Scalar::all(1);
log(mag, mag);
// crop the spectrum, if it has an odd number of rows or columns
// 大きさが奇数の場合削る
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));
int cx = mag.cols/2;
int cy = mag.rows/2;
// rearrange the quadrants of Fourier image
// so that the origin is at the image center
//低周波が四隅になっているので入れ替えて中心を直流に
// MATLABならfftshiftで済むが
// OpenCVの場合はMatにRectを指定してROI(Region of Interest)
// をq0-3の4つ作成し入れ替える
Mat tmp;
Mat q0(mag, Rect(0, 0, cx, cy));
Mat q1(mag, Rect(cx, 0, cx, cy));
Mat q2(mag, Rect(0, cy, cx, cy));
Mat q3(mag, Rect(cx, cy, cx, cy));
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//最小最大を0-1に割当てるよう正規化
normalize(mag, mag, 0, 1, CV_MINMAX);
imshow("spectrum magnitude", mag);
waitKey();
return 0;
}C++のサンプルのdft.cppを基に作成
- Matは無く，取得される画像はNumPyのndarrayそのものなので，ncolやnrowは無い
- merge，clone，ROI等の行列操作はNumPyのやり方に従う
- 筆者はNumPy素人
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# dft.py
# 2011/8/29 ScyPyを使った方が速そうだが samples/cpp/dft.cpp に似せて作成
# Eiichiro Momma
import cv2 as cv
import numpy as np
#今のところcv２ではMatはndarrayそのもの
if __name__ == '__main__':
img = cv.imread("cos2img.jpg)",0)
M = cv.getOptimalDFTSize(img.shape[0])
N = cv.getOptimalDFTSize(img.shape[1])
padded = np.zeros([M, N], dtype=float)
complexImg = np.zeros([M, N, 2], dtype=float)
#copyMakeBorderは無駄にややこしいので使わない
complexImg[ (M-img.shape[0])/2: (M-img.shape[0])/2+img.shape[0], \
(N-img.shape[1])/2: (N-img.shape[1])/2+img.shape[1], 0 ] = img/255.0
cv.dft(complexImg, complexImg)
#planesに分ける作業は不要
padded = cv.magnitude(complexImg[:,:,0], complexImg[:,:,1])
padded = cv.log(padded+1)
padded = padded[0:padded.shape[0] & -2, 0:padded.shape[1] & -2]
cv.normalize(padded, padded, 0, 1, cv.NORM_MINMAX)
#象限の入れ替え。恐らくもっとエレガントなやり方がある
cr = padded.shape[0]/2
cc = padded.shape[1]/2
tmp = padded[0:cr, 0:cc].copy()
padded[0:cr, 0:cc] = padded[cr: cr*2, cc: cc*2].copy()
padded[cr: cr*2, cc: cc*2] = tmp.copy()
tmp = padded[0:cr, cc: cc*2].copy()
padded[0:cr, cc: cc*2] = padded[cr: cr*2, 0: cc].copy()
padded[cr: cr*2, 0:cc] = tmp.copy()
cv.imshow("org", img)
cv.imshow("DFT",padded)
cv.waitKey(0)OpenCV 窓関数から流用
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 );
- src
- -元の画像、実数か複素数。縦横は偶数であること。
- dst
- -出力。srcと同じ形式にする。
- flags
- -動作のフラグ
- --CV_DXT_FORWARD 1Dか2Dの順変換。出力は生の値(正規化されない)
- --CV_DXT_INVERSE 逆変換
- --CV_DXT_SCALE 順変換。出力を行列の総数で割る。
- --CV_DXT_INV_SCALE 上の逆変換版
- --CV_DXT_ROWS 行ごとに独立して変換。3D以上をやりたい場合に使える(らしい)
- nonzero_rows
- -行の指定。指定しないか0なら全行変換。
実数のみの場合は1channelとなるのだが、出力が下記のように変態的な並びになる。
Re Y0,0 Re Y0,1 Im Y0,1 Re Y0,2 Im Y0,2 ... Re Y0,N/2-1 Im Y0,N/2-1 Re Y0,N/2
Re Y1,0 Re Y1,1 Im Y1,1 Re Y1,2 Im Y1,2 ... Re Y1,N/2-1 Im Y1,N/2-1 Re Y1,N/2
Im Y1,0 Re Y2,1 Im Y2,1 Re Y2,2 Im Y2,2 ... Re Y2,N/2-1 Im Y2,N/2-1 Im Y2,N/2
............................................................................................
Re YM/2-1,0 Re YM-3,1 Im YM-3,1 Re YM-3,2 Im YM-3,2 ... Re YM-3,N/2-1 Im YM-3,N/2-1 Re YM-3,N/2
Im YM/2-1,0 Re YM-2,1 Im YM-2,1 Re YM-2,2 Im YM-2,2 ... Re YM-2,N/2-1 Im YM-2,N/2-1 Im YM-2,N/2
Re YM/2,0 Re YM-1,1 Im YM-1,1 Re YM-1,2 Im YM-1,2 ... Re YM-1,N/2-1 Im YM-1,N/2-1 Im YM-1,N/2
後処理が洒落にならないので2channel画像にして、1channel目をReal、2channelをImaginaryにする。
よくある中心が原点の画像にはならず、(0,0)が原点でx,yが入れ替わっている。
1channel画像
↓
Imaginaryが0の2channel画像
↓
cvDFT
↓
原点を中心に変換
↓
sqrt(real^2+image^2)
↓
最大最小値を使ってスケーリング
という手順が必要になる
Yahoo! Groups(Yahoo.com)のOpenCV MLで投稿されていたのを流用してよくある出力になる関数を作成した。
IplImage *myDFT(IplImage *src)
{
IplImage* fftImg=0;
IplImage* realImg=0;
IplImage* imageImg=0;
IplImage* prealImg=0;
IplImage* pimageImg=0;
IplImage* modulo=0;
//for only 1 channel image
if (src->nChannels > 1)
return NULL;
//for FFT
realImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
imageImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
fftImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,2);
prealImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
pimageImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
modulo=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
cvConvertScale(src,realImg,1,0); // uchar image to double image
cvZero(imageImg);
cvMerge(realImg,imageImg,NULL,NULL,fftImg);
cvDFT(fftImg,fftImg,CV_DXT_FORWARD);
cvSplit(fftImg,realImg,imageImg,NULL,NULL);
cvPow(realImg,prealImg,2); //prealImg=realImg^2
cvPow(imageImg,pimageImg,2);//pimageImg=imageImg^2
cvAdd(prealImg,pimageImg,modulo); //modulo = prealImg+pimageImg
cvPow(modulo,modulo,0.5); //modulo = sqrt(modulo)
//rearrange quadrant
double tmp13,tmp24;
int nRow = srcImg->height;
int nCol = srcImg->width;
int cy = nRow/2; // image center
int cx = nCol/2;
int i, j;
for( j = 0; j < cy; j++ ){
for( i = 0; i < cx; i++ ){
tmp13 = CV_IMAGE_ELEM( modulo, double, j, i);
CV_IMAGE_ELEM( modulo, double, j, i) = CV_IMAGE_ELEM(modulo, double, j+cy, i+cx);
CV_IMAGE_ELEM( modulo, double, j+cy, i+cx) = tmp13;
tmp24 = CV_IMAGE_ELEM( modulo, double, j, i+cx);
CV_IMAGE_ELEM( modulo, double, j, i+cx) = CV_IMAGE_ELEM(modulo, double, j+cy, i);
CV_IMAGE_ELEM( modulo, double, j+cy, i) = tmp24;
}
}
//normalize
CvPoint maxPoint,minPoint;
double max,min;
cvMinMaxLoc(modulo,&min,&max,&minPoint,&maxPoint);
double scale=255.0/(max-min);
double shift=min*scale;
cvConvertScale(modulo,modulo,scale,shift);
//clean up
cvReleaseImage(&fftImg);
cvReleaseImage(&realImg);
cvReleaseImage(&imageImg);
cvReleaseImage(&prealImg=0);
cvReleaseImage(&pimageImg=0);
return modulo;
}使い方は IplImage *dst; //受取用ポインタ dst = myDFT(1channel画像へのポインタ); で良い。
変換前
Green成分の変換後
OpenCV Library Wikiを参考に pixvalue = ((double*)(dst->imageData+y*dst->widthStep))[x]; のように取得する。
cvMinMaxLocの使い方が異なる。
from opencv.cv import *
from opencv.highgui import *
import sys
def main():
src = cvLoadImage(sys.argv[1],0)
realImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1)
imageImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1)
fftImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2)
prealImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1)
pimageImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1)
modulo = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1)
cvConvertScale(src,realImg,1,0)
cvZero(imageImg)
cvMerge(realImg,imageImg,None,None,fftImg)
cvDFT(fftImg,fftImg,CV_DXT_FORWARD)
cvSplit(fftImg,realImg,imageImg,None,None)
cvPow(realImg,prealImg,2)
cvPow(imageImg,pimageImg,2)
cvAdd(prealImg,pimageImg,modulo)
cvPow(modulo,modulo,0.5)
#rearrange quadrant
nRow = src.height
nCol = src.width
cy = nRow/2
cx = nCol/2
for j in range(cy):
for i in range(cx):
tmp13 = cvGet2D(modulo,j,i)
tmp31 = cvGet2D(modulo,j+cy, i+cx)
cvSet2D(modulo,j,i,tmp31)
cvSet2D(modulo,j+cy,i+cx,tmp13)
tmp24 = cvGet2D(modulo,j,i+cx)
tmp42 = cvGet2D(modulo,j+cy,i)
cvSet2D(modulo,j,i+cx, tmp42)
cvSet2D(modulo,j+cy,i, tmp24)
minval=1000.0
maxval=0.0
minPoint = cvPoint(0,0)
maxPoint = cvPoint(0,0)
minval,maxval= cvMinMaxLoc(modulo,minPoint,maxPoint,None )
scale = 255.0/(maxval-minval)
shift = minval*scale
cvConvertScale(modulo,modulo,scale,shift)
cvNamedWindow("src",1)
cvNamedWindow("modulo",1)
cvShowImage("src",src)
cvShowImage("modulo",modulo)
cvWaitKey(0)
cvDestroyWindow("src")
cvDestroyWindow("modulo")
if __name__ == '__main__':
main()1ch画像は扱いが面倒なので上と同様に2chとする。 その際は1chスペクトラム画像用のcvMulSpectrumsではなく、多チャンネルの要素ごとの乗算であるcvMulを使う(多分)。 ※ざっと書いたのでいい加減
from opencv.cv import *
from opencv.highgui import *
FLT_RADIUS = 20
#原画像
Aorg = cvLoadImage("lena.jpg)",0)
#複素画像用
Areal = cvCreateImage(cvGetSize(Aorg),IPL_DEPTH_64F,1)
Aimg = cvCloneImage(Areal)
Breal = cvCloneImage(Areal)
Bimg = cvCloneImage(Areal)
#複素画像
FFTimg = cvCreateImage(cvGetSize(Aorg),IPL_DEPTH_64F,2)
FLTimg = cvCloneImage(FFTimg)
cvConvertScale(Aorg,Areal,1.0/255.0,0)
cvSetZero(Aimg)
cvSetZero(Breal)
cvSetZero(Bimg)
#各パートをFFTimg, FLTimgそれぞれにマージ
cvMerge(Areal,Aimg,None,None,FFTimg)
cvMerge(Breal,Bimg,None,None,FLTimg)
#FLTimgに白の円を四隅に書く(簡易的ながらLPFとなる)。正確には1pixelずらして書く必要あり
cvCircle(FLTimg,cvPoint(0,0),FLT_RADIUS,cvScalar(1.0,1.0),-1)
cvCircle(FLTimg,cvPoint(0,FLTimg.height - 1),FLT_RADIUS,cvScalar(1.0,1.0),-1)
cvCircle(FLTimg,cvPoint(FLTimg.width - 1,0),FLT_RADIUS,cvScalar(1.0,1.0),-1)
cvCircle(FLTimg,cvPoint(FLTimg.width - 1,FLTimg.height - 1),FLT_RADIUS,cvScalar(1.0,1.0),-1)
cvDFT(FFTimg,FFTimg,CV_DXT_FORWARD)
#cvMulSpectrumではなくcvMul
cvMul(FFTimg,FLTimg,FFTimg)
cvDFT(FFTimg,FFTimg,CV_DXT_INVERSE)
fmin = 0.0
fmax = 0.0
cvSplit(FFTimg,Breal,Bimg,None,None)
#8bitへ正規化のため最小、最大値の取得
fmin, fmax = cvMinMaxLoc(Breal,None,None,None)
cvConvertScaleAbs(Breal,Aorg,255/fmax,0)
cvNamedWindow("test",1)
cvShowImage("test",Aorg)
cvWaitKey(0)原画像
FLT_RADIUS = 10の結果
- double型配列のf64scanedLine[width]のデータを解析する。(データ数width)
- 結果を簡易折れ線グラフで表示する。但し、IGNORE_FREQ[lp/pixel]より低周波は飽和させる。
- グラフを保存する。
※生データの書き出しはOpenCVに関係無いので省略
#define IGNORE_FREQ 3
void myDisplayDFT(int width, double *f64scanedLine)
{
//2チャンネルの1次元マトリクス
CvMat *dftData = cvCreateMat(width, 1, CV_64FC2);
//RealとIm用の1チャンネルの1次元マトリクス
CvMat *dftReData = cvCreateMat(width, 1, CV_64FC1);
CvMat *dftImData = cvCreateMat(width, 1, CV_64FC1);
//表示用ダミーデータ
//(この関数を繰り返し呼ぶ場合は下のmalloc,freeを繰り返すのは得策ではない。
//mainなどの呼び出し元で確保してポインタを渡すのがスマート)
double *f64DispDFTLine;
//グラフ用画像
//グラフは扱い易いように横軸をwidth[pixel]の幅にし、縦は255。上下左右のマージンは20とする。
IplImage *graphImage=cvCreateImage(cvSize(width+40,296),IPL_DEPTH_8U, 3);
int x;
double max=0.0;
f64DispDFTLine = (double*)malloc(sizeof(double)*width);
//貰ったデータをdftDataの1チャネルに格納。アクセスの仕方はOpenCV wikiのFAQ参照
//2チャネルは0で初期化(追加)
for (x=0; x<width; x++){
((double*)(dftData->data.ptr+dftData->step*x))[0] = f64scanedLine[x];
((double*)(dftData->data.ptr+dftData->step*x))[1] = 0.0;
}
cvDFT(dftData,dftData,CV_DXT_FORWARD);
//1,2チャンネルにそれぞれReal,Imパートのデータが入っているので分ける
cvSplit(dftData,dftReData,dftImData,NULL,NULL);
//二乗和のルート。マトリクスは使い回し。
cvPow(dftReData,dftReData,2);
cvPow(dftImData,dftImData,2);
cvAdd(dftReData,dftImData,dftReData);
cvPow(dftReData,dftReData,0.5);
//マトリクスから扱い易いただの配列にデータを移す。同時に描画用に最大値を調べる
for (x=0; x<width; x++){
f64DFTLine[x]= *((double*)(dftReData->data.ptr+dftReData->step*x));
if(x>=IGNORE_FREQ && x<width/2 && f64DFTLine[x]>max){
max=f64DFTLine[x];
}
}
//表示用配列に正規化した値を代入
for (x=0; x<width; x++){
f64DispDFTLine[x]=(255.0*f64DFTLine[x]/max)>255.0?255.0:(255.0*f64DFTLine[x]/max);
}
//グラフ作成
cvSet(graphImage,CV_RGB(255,255,255)); //白塗り
cvLine(graphImage,cvPoint(20,20),cvPoint(20,276),CV_RGB(0,0,0)); //軸
cvLine(graphImage,cvPoint(20,276),cvPoint(width+20,276),CV_RGB(0,0,0)); //軸
//スペクトルの範囲はwidth/2になるのでグラフを描画する際は横軸を倍にしている
for (x=0; x<width/2; x++){
cvLine(graphImage,cvPoint(20+x*2,276-(int)f64DispDFTLine[x]),
cvPoint(20+(x+1)*2,276-(int)f64DispDFTLine[x+1]),CV_RGB(255,0,0));
}
cvNamedWindow("Graph",1);
cvShowImage("Graph",graphImage);
cvWaitKey(0);
cvSaveImage("graph.png)",graphImage);
cvDestroyWindow("Graph");
cvReleaseImage(&graphImage);
free(f64DispDFTLine);
cvReleaseMat(&dftData);
cvReleaseMat(&dftReData);
cvReleaseMat(&dftImData);
} //2チャンネルの1次元マトリクス
CvMat *dftData = cvCreateMat(width, 1, CV_64FC2);
//RealとIm用の1チャンネルの1次元マトリクス
CvMat *dftReData = cvCreateMat(width, 1, CV_64FC1);
CvMat *dftImData = cvCreateMat(width, 1, CV_64FC1);
cvCreateMatで作成。基本的な考え方はcvCreateImageと同じ。width*1の1次元データを作る。
//貰ったデータをdftDataの1チャネルに格納。アクセスの仕方はOpenCV wikiのFAQ参照
for (x=0; x<width; x++){
((double*)(dftData->data.ptr+dftData->step*x))[0] = f64scanedLine[x];
((double*)(dftData->data.ptr+dftData->step*x))[1] = 0.0;
}
画像同様ただの配列ではないのでテクニックが必要になる。CvMat構造体のdata構造体のptr(ややこしい)を参照して、データ幅がstepとなる。ケツの[0]は1チャンネルを示す。 2チャンネルを初期化しておかないと処理系によってはデータが爆発する。
cvDFT(dftData,dftData,CV_DXT_FORWARD);
//1,2チャンネルにそれぞれReal,Imパートのデータが入っているので分ける
cvSplit(dftData,dftReData,dftImData,NULL,NULL);
//二乗和のルート。マトリクスは使い回し。
cvPow(dftReData,dftReData,2);
cvPow(dftImData,dftImData,2);
cvAdd(dftReData,dftImData,dftReData);
cvPow(dftReData,dftReData,0.5);
2次元フーリエ変換と同様。1チャネルだと面倒なので元データは2チャンネルマトリクスとして扱う。
//マトリクスから扱い易いただの配列にデータを移す。同時に描画用に最大値を調べる
for (x=0; x<width; x++){
f64DFTLine[x]= *((double*)(dftReData->data.ptr+dftReData->step*x));
if(x>=IGNORE_FREQ && x<width/2 && f64DFTLine[x]>max){
max=f64DFTLine[x];
}
}
これも格納と同様にポインタを取得してやる必要がある。
//表示用配列に正規化した値を代入
for (x=0; x<width; x++){
f64DispDFTLine[x]=(255.0*f64DFTLine[x]/max)>255.0?255.0:(255.0*f64DFTLine[x]/max);
}
//グラフ作成
cvSet(graphImage,CV_RGB(255,255,255)); //白塗り
cvLine(graphImage,cvPoint(20,20),cvPoint(20,276),CV_RGB(0,0,0)); //軸
cvLine(graphImage,cvPoint(20,276),cvPoint(width+20,276),CV_RGB(0,0,0)); //軸
//スペクトルの範囲はwidth/2になるのでグラフを描画する際は横軸を倍にしている
for (x=0; x<width/2; x++){
cvLine(graphImage,cvPoint(20+x*2,276-(int)f64DispDFTLine[x]),
cvPoint(20+(x+1)*2,276-(int)f64DispDFTLine[x+1]),CV_RGB(255,0,0));
}
ここら辺はMS-DOSの時代からプログラムをやっている人間には常套手段の描画法。
free(f64DispDFTLine);
cvReleaseMat(&dftData);
cvReleaseMat(&dftReData);
cvReleaseMat(&dftImData);
IplImageと同様にReleaseする。