完全に実力不足な理系大学生からの成長ブログ

プログラミング能力皆無、でも頑張ります。

OpenCV + NumPy で画像の周波数フィルタ

プログラムが正しく動いていないことを確認しました。
修正まで少々お待ち下さい。2016/07/29

こんにちは、よりです。
今回は画像とフィルターをフーリエ変換し、周波数空間でフィルター処理を行います。

流れ

  1. 元画像をフーリエ変換
  2. フィルターをフーリエ変換
  3. 画素ごとに乗算
  4. 逆変換

大まかにはこんな感じです。
細かい点を補足するとフィルターをフーリエ変換するときの画像サイズですね。
元画像と同じサイズにしてからフーリエ変換します。


今回の入力画像はこちら。
f:id:yori1029:20160705105925j:plain

いつも通りです。

ではコードを見てみましょう。(めっちゃ汚いけど。)

import cv2
import numpy as np

image = cv2.imread("lena.png", 0)
cv2.imshow("in",image)

#画像のフーリエ変換
fimage = np.fft.fft2(image)
"""
outfimage = np.log(np.abs(fimage) + 1)
outfimage = outfimage / np.amax(outfimage) * 255
outfimage = outfimage.astype(np.uint8)
cv2.imshow("fftimage",outfimage)
cv2.imwrite("fftimage.png",outfimage)
"""
#単純平滑化フィルタ
a = np.array([[1/25.0, 1/25.0, 1/25.0, 1/25.0, 1/25.0],
              [1/25.0, 1/25.0, 1/25.0, 1/25.0, 1/25.0],
              [1/25.0, 1/25.0, 1/25.0, 1/25.0, 1/25.0],
              [1/25.0, 1/25.0, 1/25.0, 1/25.0, 1/25.0],
              [1/25.0, 1/25.0, 1/25.0, 1/25.0, 1/25.0]])

#画像サイズを0で初期化した配列を用意
A = np.zeros((512,512)) 

#用意した配列にフィルターをコピー
A[0:5, 0:5] = a[0:5, 0:5]

#フィルターのフーリエ変換
fa = np.fft.fft2(A)
"""
outfa = np.log(np.abs(fa) + 1)
outfa = outfa / np.amax(outfa) * 255
outfa = outfa.astype(np.uint8)
cv2.imshow("fftfa",outfa)
cv2.imwrite("fftfa.png",outfa)
"""
#要素ごとの乗算
outimg = fa * image
"""
foutimg = np.log(np.abs(outimg) + 1)
foutimg = foutimg / np.amax(foutimg) * 255
foutimg = foutimg.astype(np.uint8)
cv2.imshow("aftimg",foutimg)
cv2.imwrite("aftimg.png",foutimg)
"""

#逆変換
fout = np.fft.ifft2(outimg)
out = fout.real + 128l
out = out.astype(np.int8) 
cv2.imshow("out",out)

cv2.waitKey(0)

コメントにしている部分は出力ようなのでスルーしてください。
この処理で周波数空間でのフィルター処理ができます。

5×5の単純平滑化フィルタなのでボケ画像が出力されますね。
f:id:yori1029:20160713011136p:plain

ボケてるのがわかりますよね?


一応周波数空間に変換した時の画像も載せておきます。

まず画像のフーリエ変換
f:id:yori1029:20160713011233p:plain


次にフィルターのフーリエ変換
f:id:yori1029:20160713011259p:plain


最後にこれら二つを要素ごとに乗算したもの(逆変換前)
f:id:yori1029:20160713011336p:plain


他のフィルターも試してみたいですね!
(試さずに終わり)

OpenCVでフーリエ変換(dft)2 表示編

こんにちは、よりです。
今日は下記の記事で紹介したフーリエ変換をした後の画像の表示方法を紹介したいと思います。
yori1029.hatenablog.com

とはいえ手探りで頑張ってるので私のメモ程度のレベルです。
今回はさくさく行きましょう、表示するだけやし。

前回の記事を参照してもらうとわかるんですけど dft 後のデータって2チャンネルになってるし虚部があるから扱いにくいよなーって感じです。
それを何とかして、よく見るフーリエ変換後の画像にしていきたいと思います。


全体の流れとしては

  1. フーリエ変換したデータを実部と虚部に分ける
  2. すべて実数にする(二乗して平方根
  3. 表示用に対数に変換する
  4. 表示位置を変更する(よく見るフーリエ変換後の画像にする)

って感じです。

今回フーリエ変換する画像はこちら。
f:id:yori1029:20160705105925j:plain



正直全然大したことないので参考にしたサイトのURLとコードを載せます。
下記のサイトの関数 create_fourier_magnitude_image_from_complex を参考に書きました。

opencv_sample_list_jp/main.cpp at master · YusukeSuzuki/opencv_sample_list_jp · GitHub


それでは私のわかりにくいコードも載せます。

#include <iostream> 
#include <cmath>
#include <string>
#include <sstream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <ctype.h>

using namespace cv;
using namespace std;

int main(int argc, char **argv)
{
	//入力
	Mat image = cv::imread("lena.jpg", cv::IMREAD_GRAYSCALE);
	
	//フーリエ変換用Mat
	Mat furimg;

	//実部のみのimageと虚部を0で初期化したMatをRealImaginary配列に入れる
	Mat RealIamginary[] = { Mat_<float>(image), Mat::zeros(image.size(), CV_32F) };
	
	//配列を合成
	merge(RealIamginary, 2, furimg);

	//フーリエ変換
	dft(furimg, furimg);

	//表示用
	Mat divdisplay[2];

	//フーリエ後を実部と虚部に分ける
	split(furimg, divdisplay);

	//表示用にすべて実数に
	Mat display;
	magnitude(divdisplay[0], divdisplay[1], display);

	//対数に変換する(そのため各ピクセルに1を加算)
	display += Scalar::all(1);
	log(display, display);

	//表示用に正規化
	Mat outdisplay;
	normalize(display, outdisplay, 0, 1, CV_MINMAX);

	namedWindow("aftdft");
	imshow("aftdft", outdisplay);

	waitKey(-1);

	return 0;

}


実行すると下図が出てきます。


f:id:yori1029:20160705093146p:plain


これは画像をフーリエ変換し、それを可視化したものです。
しかし普段みるフーリエ変換後の画像とは異なってますよね?
普段見るのは白い部分が中心にきてると思いますが、それは見やすいように変換している画像です。
では変換してみましょう。

追加した部分がわかるように書いたのでそこだけ見てもらえれば大丈夫です。

#include <iostream> 
#include <cmath>
#include <string>
#include <sstream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <ctype.h>

using namespace cv;
using namespace std;

int main(int argc, char **argv)
{
	//入力
	Mat image = cv::imread("lena.jpg", cv::IMREAD_GRAYSCALE);
	
	//フーリエ変換用Mat
	Mat furimg;

	//実部のみのimageと虚部を0で初期化したMatをRealImaginary配列に入れる
	Mat RealIamginary[] = { Mat_<float>(image), Mat::zeros(image.size(), CV_32F) };
	
	//配列を合成
	merge(RealIamginary, 2, furimg);

	//フーリエ変換
	dft(furimg, furimg);

	//表示用
	Mat divdisplay[2];

	//フーリエ後を実部と虚部に分ける
	split(furimg, divdisplay);

	//表示用にすべて実数に
	Mat display;
	magnitude(divdisplay[0], divdisplay[1], display);

	//対数に変換する(そのため各ピクセルに1を加算)
	display += Scalar::all(1);
	log(display, display);

	//ここから下を追加
	//___________________________________________

	const int halfW = display.cols / 2;
	const int halfH = display.rows / 2;

	Mat tmp;

	Mat q0(display,
		Rect(0, 0, halfW, halfH));
	Mat q1(display,
		Rect(halfW, 0, halfW, halfH));
	Mat q2(display,
		Rect(0, halfH, halfW, halfH));
	Mat q3(display,
		Rect(halfW, halfH, halfW, halfH));

	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);

	//____________________________________________
	//ここから上を追加

	//表示用に正規化
	Mat outdisplay;
	normalize(display, outdisplay, 0, 1, CV_MINMAX);

	namedWindow("aftdft");
	imshow("aftdft", outdisplay);

	waitKey(-1);

	return 0;

}

上記のものを実行すると下図がでてきます。
f:id:yori1029:20160705094156p:plain


これで画像のフーリエ変換を可視化することができました!!

OpenCVでのフーリエ変換(dft)


はじめまして、よりです。
今回ブログを書く決意をしたのはやったことを忘れないようにメモしようと思ったからです。


OpenCVでは離散フーリエ変換用の関数 dft が用意されています。
これをつかってフーリエ変換したいと思います。
下記のサイトに dft の解説がありますが基本的には dft(Mat(入力),Mat(出力)) って感じです。
配列操作 — opencv 2.2 documentation


さっそくフーリエ変換するぞー!!

int main(int argc, char **argv)
{
        //入力
	Mat image = cv::imread("lena.jpg", cv::IMREAD_GRAYSCALE); 

    //出力配列の用意
	Mat imgout;

    //フーリエ変換
	dft(image, imgout);

	return 0;

}

はい、エラー出ました。
出力結果
OpenCV Error: Assertion failed (type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2) in cv::dft, file ..\..\..\modules\core\src\dxt.cpp, line 2506

英語が読めないのでよくわからないんですけど dft に入れるMatのtypeがよくないんですかね。
というわけで現在入力している画像の詳細を見てみます。

void checkMat(Mat m1) {
	// 行数
	std::cout << "rows:" << m1.rows << std::endl;
	// 列数
	std::cout << "cols:" << m1.cols << std::endl;
	// 次元数
	std::cout << "dims:" << m1.dims << std::endl;
	// サイズ(2次元の場合)
	std::cout << "size[]:" << m1.size().width << "," << m1.size().height << std::endl;
	// ビット深度ID
	std::cout << "depth (ID):" << m1.depth() << "(=" << CV_64F << ")" << std::endl;
	// チャンネル数
	std::cout << "channels:" << m1.channels() << std::endl;
	// (複数チャンネルから成る)1要素のサイズ [バイト単位]
	std::cout << "elemSize:" << m1.elemSize() << "[byte]" << std::endl;
	// 1要素内の1チャンネル分のサイズ [バイト単位]
	std::cout << "elemSize1 (elemSize/channels):" << m1.elemSize1() << "[byte]" << std::endl;
	// 要素の総数
	std::cout << "total:" << m1.total() << std::endl;
	// ステップ数 [バイト単位]
	std::cout << "step:" << m1.step << "[byte]" << std::endl;
	// 1ステップ内のチャンネル総数
	std::cout << "step1 (step/elemSize1):" << m1.step1() << std::endl;
	// データは連続か?
	std::cout << "isContinuous:" << (m1.isContinuous() ? "true" : "false") << std::endl;
	// 部分行列か?
	std::cout << "isSubmatrix:" << (m1.isSubmatrix() ? "true" : "false") << std::endl;
	// データは空か?
	std::cout << "empty:" << (m1.empty() ? "true" : "false") << std::endl;

	cout << endl;
}

この関数を使えば中身がだいたい分かりそうですね。
下記のサイトのcv::Matの様々なプロパティに載ってます。
cv::Matの基本処理 — OpenCV-CookBook

出力結果
rows:512
cols:512
dims:2
size[]:512,512
depth (ID):0(=6)
channels:1
elemSize:1[byte]
elemSize1 (elemSize/channels):1[byte]
total:262144
step:512[byte]
step1 (step/elemSize1):512
isContinuous:true
isSubmatrix:false
empty:false

いろいろ調べた結果、フーリエ変換をした時に出てくる複素数を扱うためにはMatのチャンネルを2にする必要があるそうです。
今回imageに画像を入れましたがこのままだと1チャンネルしかないようなので2チャンネルにします。
mergeという関数が適しているようなのでこれを用いましょう。
mergeさせるのはimageと同じ大きさの複素数のMatです。

int main(int argc, char **argv)
{
	//入力
	Mat image = cv::imread("lena.jpg", cv::IMREAD_GRAYSCALE);
	
	//フーリエ変換用Mat
	Mat furimg;

	//実部のみのimageと虚部を0で初期化したMatをRealImaginary配列に入れる
	Mat RealIamginary[] = { Mat_<float>(image), Mat::zeros(image.size(), CV_32F) };
	
	//配列を合成
	merge(RealIamginary, 2, furimg);

	//フーリエ変換
	dft(furimg, furimg);

	return 0;

}

今回はエラーが出ませんでした!
おそらくfurimgにフーリエ変換後の値が入っていると思います。

今日はこれからバイトなのでここまで。