2次元フーリエ変換,フーリエ逆変換
YouTube で見た以下のビデオが面白かったので,似たようなものを OpenCV のサンプルとして作ってみました.cv::dft,cv::idft を利用しています.処理の詳細は,ビデオをご覧ください.
C++
#include <iostream> #include <cv.h> #include <highgui.h> using namespace cv; using namespace std; // class declaration class DftIdftApp { public: DftIdftApp(const string filename); ~DftIdftApp(){}; void calcMagImage(); void showOrgImage(){imshow(org_win, src_img);} void showMagImage(){imshow(mag_win, mag_img);} void showIDFTImage(){imshow(idft_win, idft_img);} void shiftDFT(Mat &src, Mat &dst); void calcIDFT(bool all=false); void clear(){idft_img.setTo(0);} static void onMouse(int event, int x, int y, int flags, void* param); private: Mat src_img, mag_img; Mat Re_img, Im_img, Complex_img; Mat zero, dft_src, dft_dst, dft_dst_p; Mat idft_img; vector<Mat> mv; string org_win, mag_win, idft_win; int src_cols, src_rows, dft_cols, dft_rows; }; // constructor DftIdftApp::DftIdftApp(const string filename) :org_win("Original"), mag_win("Magnitude"), idft_win("IDTF") { src_img = imread(filename, 0); if(!src_img.data) return; Size s_size = src_img.size(); src_cols = s_size.width; src_rows = s_size.height; Re_img = Mat(s_size, CV_64F); Im_img = Mat::zeros(s_size, CV_64F); Complex_img = Mat(s_size, CV_64FC2); src_img.convertTo(Re_img, CV_64F); mv.push_back(Re_img); mv.push_back(Im_img); merge(mv, Complex_img); idft_img = zero = Mat::zeros(s_size, CV_64F); namedWindow(org_win, CV_WINDOW_AUTOSIZE); namedWindow(mag_win, CV_WINDOW_AUTOSIZE); namedWindow(idft_win, CV_WINDOW_AUTOSIZE); cvSetMouseCallback(mag_win.c_str(), DftIdftApp::onMouse, this); } // calc magnitude image void DftIdftApp::calcMagImage() { dft_rows = getOptimalDFTSize(src_rows); dft_cols = getOptimalDFTSize(src_cols); dft_src = Mat::zeros(dft_rows, dft_cols, CV_64FC2); Mat roi(dft_src, Rect(0, 0, src_cols, src_rows)); Complex_img.copyTo(roi); dft(dft_src, dft_dst); dft_dst_p = Mat::zeros(dft_dst.size(), dft_dst.type()); //split(dft_dst.mul(dft_dst), mv); //sqrt(mv[0]+mv[1], mv[0]); split(dft_dst, mv); magnitude(mv[0], mv[1], mv[0]); log(mv[0]+1, mv[0]); // for ver. 2.1 or later shiftDFT(mv[0], mv[0]); //double min, max; //minMaxLoc(mv[0], &min, &max); //mag_img = mv[0]*1.0/(max-min) - 1.0*min/(max-min); normalize(mv[0], mag_img, 0, 1, CV_MINMAX); } // swap 1,3 and 2,4 quadrants respectively void DftIdftApp::shiftDFT(Mat &src, Mat& dst) { Mat tmp; int cx = src.cols/2; int cy = src.rows/2; for(int i=0; i<=cx; i+=cx) { Mat qs(src, Rect(i^cx,0,cx,cy)); Mat qd(dst, Rect(i,cy,cx,cy)); qs.copyTo(tmp); qd.copyTo(qs); tmp.copyTo(qd); } } // mouse event callback void DftIdftApp::onMouse(int event, int x, int y, int flags, void* param) { DftIdftApp *app = static_cast<DftIdftApp*>(param); Mat &dft_dst = app->dft_dst; Mat &dft_dst_p = app->dft_dst_p; Mat &mag_img = app->mag_img; int cx = app->src_cols/2; int cy = app->src_rows/2; int mx=x, my=y; int w = app->dft_cols; int h = app->dft_rows; switch(event) { case CV_EVENT_MOUSEMOVE: if((flags&CV_EVENT_FLAG_LBUTTON)==0) return; case CV_EVENT_LBUTTONUP: if(flags&CV_EVENT_FLAG_CTRLKEY) { // LeftButton+CTRL int step = 5; for(int j=-step; j<=step; j++) { my = y+j + ((y+j)<cy ? cy:-cy); if(y+j<0 || y+j>=h) continue; double *from = dft_dst.ptr<double>(my); double *to = dft_dst_p.ptr<double>(my); double *mag = mag_img.ptr<double>(y+j); for(int i=-step; i<=step; i++) { mx = x+i + ((x+i)<cx ? cx:-cx); if(x+i<0 || x+i>=w) break; to[(mx)*2+0] = from[(mx)*2+0]; to[(mx)*2+1] = from[(mx)*2+1]; mag[x+i] = 0; } } } else { // LeftButton mx += x<cx ? cx:-cx; my += y<cy ? cy:-cy; double *from = dft_dst.ptr<double>(my); double *to = dft_dst_p.ptr<double>(my); double *mag = mag_img.ptr<double>(y); to[(mx)*2+0] = from[(mx)*2+0]; to[(mx)*2+1] = from[(mx)*2+1]; mag[x] = 0; } break; default: return; } app->calcIDFT(); } // Inverse Discrete Fourier Transforma void DftIdftApp::calcIDFT(bool all) { if(all) { dft_dst.copyTo(dft_dst_p); mag_img.setTo(0); } double min, max; idft(dft_dst_p, dft_src); split(dft_src, mv); minMaxLoc(mv[0], &min, &max); idft_img = Mat(mv[0]*1.0/max, Rect(0, 0, src_cols, src_rows)); } int main(int argc, char *argv[]) { const string filename = argc > 1 ? argv[1] : "../image/Fourier.png"; DftIdftApp app(filename); cout << "Usage: click or drag on Magnitude Image.\n" << "Hot keys: \n" "\tESC - quit the program\n" "\ta - select all pixels\n" "\tr - restore original images\n" "\tleft mouse button - select one pixel\n" "\tCTRL+left mouse button - select a pixel and neighbors\n"; app.showOrgImage(); app.calcMagImage(); while(1) { int c = waitKey(10); switch((char)c) { case '\x1b': // exit return 0; case 'a': // IDFT all app.calcIDFT(true); break; case 'r': // reset app.clear(); app.calcMagImage(); break; } app.showIDFTImage(); app.showMagImage(); } return 0; }
実行結果例
(1)動画(サンプル実行の例)
(2)静止画
(左)対数パワースペクトル, (右)IDFT結果 左画像の黒い領域は,選択された(IDFTの対象となる)周波数,方向の正弦波を表します.この領域は,対数パワースペクトル画像上をマウスでクリック,ドラッグなどすることで選択できます.多くの正弦波を重ね合わせることで,より元の画像に近づく様子が分かります.