离散傅里叶变换DFT应用

#include <opencv2/opencv.hpp>

const int PIXEL_EXCLUDED = 64;
const cv::Size PATCH_SIZE = cv::Size(256, 256);
const float STEP = 0.5;

bool SortBySumFFT(const std::pair<cv::Point, double> &a,
	const std::pair<cv::Point, double> &b)
{
	return (a.second > b.second);
}

cv::Mat CvtImgs16To8(cv::Mat img)
{
	cv::Mat imgs_8bit;
	normalize(img, imgs_8bit, 0, 255, cv::NORM_MINMAX, CV_8U);
	return imgs_8bit;
}

void Recomb(cv::Mat &src, cv::Mat &dst)
{
	int cx = src.cols >> 1;
	int cy = src.rows >> 1;
	cv::Mat tmp;
	tmp.create(src.size(), src.type());
	src(cv::Rect(0, 0, cx, cy)).copyTo(tmp(cv::Rect(cx, cy, cx, cy)));
	src(cv::Rect(cx, cy, cx, cy)).copyTo(tmp(cv::Rect(0, 0, cx, cy)));
	src(cv::Rect(cx, 0, cx, cy)).copyTo(tmp(cv::Rect(0, cy, cx, cy)));
	src(cv::Rect(0, cy, cx, cy)).copyTo(tmp(cv::Rect(cx, 0, cx, cy)));
	dst = tmp;
}

void ForwardFFT(cv::Mat &src, cv::Mat *FImg, bool do_recomb = true)
{
	int M = cv::getOptimalDFTSize(src.rows);
	int N = cv::getOptimalDFTSize(src.cols);
	cv::Mat padded;
	copyMakeBorder(src, padded, 0, M - src.rows, 0, N - src.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
	cv::Mat planes[] = { cv::Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
	cv::Mat complexImg;
	merge(planes, 2, complexImg);
	cv::dft(complexImg, complexImg);
	split(complexImg, planes);
	planes[0] = planes[0](cv::Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2));
	planes[1] = planes[1](cv::Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2));
	if (do_recomb)
	{
		Recomb(planes[0], planes[0]);
		Recomb(planes[1], planes[1]);
	}
	planes[0] /= float(M*N);
	planes[1] /= float(M*N);
	FImg[0] = planes[0].clone();
	FImg[1] = planes[1].clone();
}

void AutoRegistration(cv::Mat img, std::vector<cv::Rect> &patches)
{
	if (img.type() == CV_16UC1)
		img = CvtImgs16To8(img);

	cv::Rect roi(PIXEL_EXCLUDED, PIXEL_EXCLUDED, img.cols - 2 * PIXEL_EXCLUDED, img.rows - 2 * PIXEL_EXCLUDED);
	cv::Mat image_roi = img(roi).clone();

	std::vector<std::pair<cv::Point, float>> Fourier_result;

	for (int j = 0; j < ceil(((float)image_roi.rows / PATCH_SIZE.height - 1) / STEP) + 1; j++)
	{
		for (int i = 0; i < ceil(((float)image_roi.cols / PATCH_SIZE.height - 1) / STEP) + 1; i++)
		{
			cv::Rect PATCH(i * PATCH_SIZE.width * STEP, j * PATCH_SIZE.height * STEP, PATCH_SIZE.width, PATCH_SIZE.height);
			cv::Mat img_rect = image_roi(PATCH).clone();

			cv::Mat F0[2];
			cv::Mat f0;
			ForwardFFT(img_rect, F0);
			cv::magnitude(F0[0], F0[1], f0);

			cv::Mat f0_LPF;
			cv::Mat GaussianX = cv::getGaussianKernel(256, 15, CV_32F);
			cv::Mat GaussianY = cv::getGaussianKernel(256, 15, CV_32F);
			cv::Mat Gaussian_LPF = GaussianX * GaussianY.t();
			float center = Gaussian_LPF.ptr<float>(127)[127];
			Gaussian_LPF /= center;
			f0_LPF = f0.mul(Gaussian_LPF);

			cv::Mat kernel = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(3, 3));
			cv::morphologyEx(f0_LPF, f0_LPF, cv::MORPH_OPEN, kernel);

			float sum_mag = cv::sum(f0_LPF)[0];
			cv::Point sum_mag_idx = cv::Point(i, j);
			Fourier_result.push_back(std::make_pair(cv::Point(i, j), sum_mag));
		}
	}

	std::sort(Fourier_result.begin(), Fourier_result.end(), SortBySumFFT);

	for (int i = 0; i < 3; i++)
	{
		patches.push_back(cv::Rect(Fourier_result[i].first.x * PATCH_SIZE.width * STEP + PIXEL_EXCLUDED,
			Fourier_result[i].first.y * PATCH_SIZE.height * STEP + PIXEL_EXCLUDED, PATCH_SIZE.width, PATCH_SIZE.height));
	}
}

int main()
{
	const char* filename = "C:\Zmg12161632-1-3.tif";
	cv::Mat img = imread(filename, cv::IMREAD_UNCHANGED);

	std::vector<cv::Rect> patches;
	AutoRegistration(img, patches);

	return 0;
}

推荐阅读：

OpenCV核心模块（7）DFT离散傅里叶变换