#include <opencv2/opencv.hpp>
#include <cmath>
#include <complex>
#include <vector>
#include <iostream>
using namespace cv;
using namespace std;
// 位反转函数（用于蝶形运算准备）
int bitReverse(int num, int log2n) {
int result = 0;
for (int i = 0; i < log2n; ++i) {
if (num & (1 << i))
result |= 1 << (log2n - 1 - i);
}
return result;
}
// 1D FFT实现（基2-FFT蝶形运算）
void fft(vector<complex<float>>& data, bool inverse = false) {
const int n = data.size();
const int log2n = static_cast<int>(log2(n));
// 位反转重排
vector<complex<float>> temp(n);
for (int i = 0; i < n; ++i) {
int rev = bitReverse(i, log2n);
temp[i] = data[rev];
}
// 蝶形运算
for (int s = 1; s <= log2n; ++s) {
int m = 1 << s; // 当前处理块大小
float ang = 2 * CV_PI / m * (inverse ? 1 : -1);
complex<float> wm(cos(ang), sin(ang));
for (int k = 0; k < n; k += m) {
complex<float> w(1.0f);
for (int j = 0; j < m / 2; ++j) {
complex<float> t = w * temp[k + j + m / 2];
complex<float> u = temp[k + j];
temp[k + j] = u + t;
temp[k + j + m / 2] = u - t;
w *= wm;
}
}
}
// 逆变换缩放
if (inverse) {
for (auto& x : temp) x /= n;
}
data = temp;
}
// 2D FFT实现
void fft2d(Mat& complexImage, bool inverse = false) {
int rows = complexImage.rows;
int cols = complexImage.cols;
// 行变换
for (int y = 0; y < rows; ++y) {
vector<complex<float>> rowData;
for (int x = 0; x < cols; ++x) {
rowData.push_back(complexImage.at<Vec2f>(y, x).val[0] +
complexImage.at<Vec2f>(y, x).val[1] * 1if);
}
fft(rowData, inverse);
for (int x = 0; x < cols; ++x) {
complexImage.at<Vec2f>(y, x) = Vec2f(rowData[x].real(), rowData[x].imag());
}
}
// 列变换
for (int x = 0; x < cols; ++x) {
vector<complex<float>> colData;
for (int y = 0; y < rows; ++y) {
colData.push_back(complexImage.at<Vec2f>(y, x).val[0] +
complexImage.at<Vec2f>(y, x).val[1] * 1if);
}
fft(colData, inverse);
for (int y = 0; y < rows; ++y) {
complexImage.at<Vec2f>(y, x) = Vec2f(colData[y].real(), colData[y].imag());
}
}
}
// 频域中心化 (fftshift)
void fftShift(Mat& mag) {
int cx = mag.cols / 2;
int cy = mag.rows / 2;
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));
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
// 计算MSE
double calculateMSE(const Mat& img1, const Mat& img2) {
if (img1.size() != img2.size() || img1.type() != img2.type()) {
cerr << "Images must have same size and type" << endl;
return -1;
}
Mat diff;
absdiff(img1, img2, diff);
diff = diff.mul(diff);
Scalar s = sum(diff);
return s[0] / (img1.total() * img1.channels());
}
// 计算PSNR
double calculatePSNR(const Mat& img1, const Mat& img2) {
double mse = calculateMSE(img1, img2);
if (mse <= 1e-10) return 100; // 避免除以0
return 10.0 * log10((255.0 * 255.0) / mse);
}
// 创建低通滤波器掩膜
Mat createLowPassFilter(Size size, double D0) {
Mat mask = Mat::zeros(size, CV_32F);
Point center(size.width / 2, size.height / 2);
for (int y = 0; y < size.height; y++) {
for (int x = 0; x < size.width; x++) {
double d = sqrt(pow(x - center.x, 2) + pow(y - center.y, 2));
if (d <= D0) {
mask.at<float>(y, x) = 1.0f;
}
}
}
return mask;
}
int main() {
// 读取灰度图像
Mat img = imread("input.jpg", IMREAD_GRAYSCALE);
if (img.empty()) {
cerr << "Could not read image" << endl;
return -1;
}
// 调整大小为2的幂次（FFT要求）
int newRows = pow(2, ceil(log2(img.rows)));
int newCols = pow(2, ceil(log2(img.cols)));
resize(img, img, Size(newCols, newRows));
// 转换为浮点型并归一化
Mat floatImg;
img.convertTo(floatImg, CV_32F, 1.0 / 255.0);
// 准备复数矩阵（实部=图像，虚部=0）
Mat planes[] = { Mat_<float>(floatImg), Mat::zeros(floatImg.size(), CV_32F) };
Mat complexImg;
merge(planes, 2, complexImg);
// ===== 执行FFT =====
fft2d(complexImg);
// 分离实部和虚部
split(complexImg, planes);
// 计算幅度谱 (log scale)
Mat mag;
magnitude(planes[0], planes[1], mag);
mag += Scalar::all(1); // 避免log(0)
log(mag, mag);
normalize(mag, mag, 0, 1, NORM_MINMAX);
// 频域中心化
fftShift(mag);
// 显示幅度谱
imshow("Frequency Spectrum", mag);
imwrite("spectrum.jpg", mag * 255);
waitKey(0);
// ===== 重建原始图像 (IFFT) =====
fft2d(complexImg, true); // 逆变换
// 获取实部作为重建图像
split(complexImg, planes);
Mat reconstructed;
normalize(planes[0], reconstructed, 0, 1, NORM_MINMAX);
// 显示重建图像
imshow("Reconstructed Image", reconstructed);
imwrite("reconstructed.jpg", reconstructed * 255);
// 计算重建MSE/PSNR
double mseOriginal = calculateMSE(floatImg, reconstructed);
double psnrOriginal = calculatePSNR(floatImg, reconstructed * 255);
cout << "MSE (Original): " << mseOriginal << endl;
cout << "PSNR (Original): " << psnrOriginal << " dB" << endl;
// ===== 低通滤波处理 =====
// 创建低通滤波器（截止频率=图像尺寸的15%）
double D0 = min(img.rows, img.cols) * 0.45;
Mat lpfMask = createLowPassFilter(complexImg.size(), D0);
// 应用滤波器（在频域中心化状态）
fftShift(mag); // 确保频域在中心
Mat filteredPlanes[] = { planes[0].mul(lpfMask), planes[1].mul(lpfMask) };
Mat filteredComplex;
merge(filteredPlanes, 2, filteredComplex);
// 逆变换得到滤波后图像
fft2d(filteredComplex, true);
split(filteredComplex, filteredPlanes);
Mat filteredImg;
normalize(filteredPlanes[0], filteredImg, 0, 1, NORM_MINMAX);
// 显示滤波后图像
imshow("Low-Pass Filtered", filteredImg);
imwrite("filtered.jpg", filteredImg * 255);
// 计算滤波后PSNR
double mseFiltered = calculateMSE(floatImg, filteredImg);
double psnrFiltered = calculatePSNR(floatImg, filteredImg * 255);
cout << "MSE (Filtered): " << mseFiltered << endl;
cout << "PSNR (Filtered): " << psnrFiltered << " dB" << endl;
waitKey(0);
return 0;
}在这个代码上修改并给出完整代码
这篇博客是学习OpenCV的笔记，重点介绍了如何使用Python实现图像处理的几个关键函数，包括：copyMakeBorder用于添加图像边框，merge用于通道合并，dft进行离散傅里叶变换，log计算数组元素的自然对数，以及normalize进行图像数据的归一化操作。
扫一扫
803
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
