低通滤波器和高通滤波器

应用于图像低通滤波器和高通滤波器的实现

需要用到傅里叶变换

#include <opencv2/opencv.hpp>
#include <Eigen>
#include <iostream>
#include <vector>
#include <cmath>
#include <complex>#define M_PI       3.14159265358979323846   // pi
// 对数幅度缩放
Eigen::MatrixXd logAmplitudeSpectrum(const Eigen::MatrixXd& spectrum) {return (spectrum.array() + 1).log();
}//  乘幂尺度变换
Eigen::MatrixXd powerLawScaling(const Eigen::MatrixXd& spectrum, double gamma) {return spectrum.array().pow(gamma);
}// 归一化 to [0, 1]
Eigen::MatrixXd normalize(const Eigen::MatrixXd& spectrum) {double minVal = spectrum.minCoeff();double maxVal = spectrum.maxCoeff();return (spectrum.array() - minVal) / (maxVal - minVal);
}// 增强频谱显示
Eigen::MatrixXd enhanceSpectrumDisplay(const Eigen::MatrixXd& spectrum, double gamma = 1) {Eigen::MatrixXd logSpectrum = logAmplitudeSpectrum(spectrum);Eigen::MatrixXd powerScaledSpectrum = powerLawScaling(logSpectrum, gamma);return normalize(powerScaledSpectrum);
}// 从复数矩阵中获取幅度谱和相位谱
void getAmplitudeAndPhaseSpectra(const Eigen::MatrixXcd& data, Eigen::MatrixXd& amplitude, Eigen::MatrixXd& phase) {amplitude = data.array().abs().matrix();phase = data.array().arg().matrix();
}// 从幅度谱和相位谱重构复数矩阵
void reconstructFromAmplitudeAndPhase(const Eigen::MatrixXd& amplitude,const Eigen::MatrixXd& phase,Eigen::MatrixXcd& data)
{data = (amplitude.array() * (phase.array().cos() + std::complex<double>(0, 1) * phase.array().sin())).matrix();}// 1：**振幅谱低频移动到中心（频率平移）**:方便操作,利用象限对称互换 fft之后
Eigen::MatrixXd fftShift(const Eigen::MatrixXd &F)
{int M = F.rows();int N = F.cols();Eigen::MatrixXd F_shifted(M, N);int mid_M = M >> 1;int mid_N = N >> 1;// 交换第一象限和第三象限F_shifted.block(0, 0, mid_M, mid_N) = F.block(mid_M, mid_N, mid_M, mid_N);F_shifted.block(mid_M, mid_N, mid_M, mid_N) = F.block(0, 0, mid_M, mid_N);// 交换第二象限和第四象限F_shifted.block(0, mid_N, mid_M, mid_N) = F.block(mid_M, 0, mid_M, mid_N);F_shifted.block(mid_M, 0, mid_M, mid_N) = F.block(0, mid_N, mid_M, mid_N);return F_shifted;
}// 2：振幅谱低频移动到中心  **图像进行-1幂操作**:然后经过fft变换后，低频会在振幅谱中间 fft之前
Eigen::MatrixXd imageShift(const Eigen::MatrixXd &image)
{int M = image.rows();int N = image.cols();Eigen::MatrixXd F_shifted(M, N);// 通过乘以 (-1)^(u+v) 来平移频率for (int u = 0; u < M; ++u){for (int v = 0; v < N; ++v){//(u + v) & 1 通过位的判断末尾如果是1 为奇数，为0，为偶数。 -1的奇次幂还是-1，偶次幂为1.F_shifted(u, v) = image(u, v) * ((u + v) & 1 ? -1 : 1);}}return F_shifted;
}Eigen::MatrixXd tMatrixXd(const cv::Mat &img)
{Eigen::MatrixXd image(img.rows, img.cols);for (int i = 0; i < img.rows; ++i){for (int j = 0; j < img.cols; ++j){image(i, j) = img.at<float>(i, j);}}return image;
}
cv::Mat tMat(const Eigen::MatrixXd &img)
{cv::Mat image(img.rows(), img.cols(), CV_32FC1);for (int i = 0; i < img.rows(); ++i){for (int j = 0; j < img.cols(); ++j){image.at<float>(i, j) = img(i, j);}}return image;
}//inv =true 为逆变换，false为 正变换
//FFT 使用Eigen库中的向量表示，方便二维计算
Eigen::VectorXcd FFT(const Eigen::VectorXcd& y, bool inv = false)
{Eigen::VectorXcd x = y;//数据的大小int N = x.size();// 按位反转for (int i = 1, j = 0; i < N; i++){//bit = N/2int bit = N >> 1;/** 我们正在查看j的二进制表示中的每一位。从左到右，我们检查每一位是否为1。对于每一个为1的位，我们将其反转为0。当我们遇到第一个为0的位时，循环终止。*/for (; j & bit; bit >>= 1){j ^= bit;}//进行异或运算（XOR运算的工作原理是：当两个比较的位相同时，结果是0；当两个比较的位不同时，结果是1。）j ^= bit;//当 i >j 表示前面已经替换了位置， i==j,表示位置不用变if (i < j){std::swap(x[i], x[j]);}}// 预先计算旋转因子Eigen::VectorXcd w(N >> 1);//逆变换 = 1；傅里叶变换  = -1；double imag_i = inv ? 1.0 : -1.0;std::complex<double> tempW = std::exp(std::complex<double>(0, imag_i * 2 * M_PI / N));w[0] = 1;for (int i = 1; i < (N >> 1); ++i){// w[i] = std::polar(1.0, imag_i * 2 * M_PI * i / N);//w[i] = std::pow(tempW, i);w[i] = w[i - 1] * tempW;}// 迭代FFTfor (int len = 2; len <= N; len <<= 1){int halfLen = len >> 1;int step = N / len;for (int i = 0; i < N; i += len){for (int j = 0; j < halfLen; ++j){std::complex<double> u = x[i + j];std::complex<double> v = x[i + j + halfLen] * w[j * step];x[i + j] = u + v;x[i + j + halfLen] = u - v;}}}//使用逆变换时if (inv){for (std::complex<double>& a : x) {a /= N;}}return x;
}//检查二维数据大小是否为2的幂数，不是则填充为0到2的幂数大小
Eigen::MatrixXd padToPowerOfTwo(const Eigen::MatrixXd& matrix) {int rows = matrix.rows();int cols = matrix.cols();// 确保输入的大小是2的幂int newRows = 1, newCols = 1;while (newRows < rows){newRows <<= 1;}while (newCols < cols){newCols <<= 1;}Eigen::MatrixXd paddedMatrix = Eigen::MatrixXd::Zero(newRows, newCols);paddedMatrix.block(0, 0, rows, cols) = matrix;return paddedMatrix;
}
// 离散傅里叶变换 -  二维
Eigen::MatrixXcd FFT2D(const Eigen::MatrixXcd& image, bool inv = false) {int rows = image.rows();int cols = image.cols();Eigen::MatrixXcd result(rows, cols);for (int i = 0; i < rows; i++) {Eigen::VectorXcd temp = image.row(i).transpose();result.row(i) = FFT(temp, inv);}for (int j = 0; j < cols; j++) {Eigen::VectorXcd temp = result.col(j);result.col(j) = FFT(temp, inv);}return result;
}//创建高频  radius越小，越减少高频的部分，越大，越还原图像
Eigen::MatrixXd highFrequency (const Eigen::MatrixXd & data,int radius)
{int rows = data.rows();int cols = data.cols();// 创建高通滤波器（圆形掩码）Eigen::MatrixXd mask = Eigen::MatrixXd::Ones(rows, cols);//中心坐标int centerRow = rows >> 1;int centerCol = cols >> 1;//找到以 radius 大小的矩形范围内//左边位置int left = centerCol - radius;//超过边界 为0left = left > 0 ? left : 0;//右边位置int right = centerCol + radius;//超过边界 为0right = right < cols ? right : cols;//上边位置int top = centerRow - radius;//超过边界 为0top = top > 0 ? top : 0;//下边位置int down = centerRow + radius;//超过边界 为0down = down < rows ? down : rows;//在正矩形内画最大的圆for (int i = top; i < down; ++i){for (int j = left; j < right; ++j){//在图像中心画圆，半径不能超过double distance = std::sqrt(std::pow(i - centerRow, 2) + std::pow(j - centerCol, 2));if (distance <= radius){mask(i, j) = 0.0;}}}return mask;}
//创建高通滤波器 -
//简单的说，就是靠近频谱图中心的低频部分给舍弃掉，远离频谱图中心的高频部分保留。通常会保留物体的边界。
Eigen::MatrixXd  highPassFilter(const Eigen::MatrixXd & image, int radius = 0)
{int rows = image.rows();int cols = image.cols();//大小变为2的幂数Eigen::MatrixXd data = padToPowerOfTwo(image);// 计算傅里叶变换Eigen::MatrixXcd transformed = FFT2D(data);//获取幅度谱和相位谱Eigen::MatrixXd amplitude, phase;getAmplitudeAndPhaseSpectra(transformed, amplitude, phase);//振幅谱移动到中心（频率平移）amplitude = fftShift(amplitude);/////增强振幅，用于观测 --  实际运算注释掉Eigen::MatrixXd amplitude1 = enhanceSpectrumDisplay(amplitude,1);cv::Mat highP = tMat(amplitude1);cv::imshow("highPassFilter",highP);/////根据radius创建高频掩码Eigen::MatrixXd mask = highFrequency(amplitude, radius);amplitude = amplitude.array() * mask.array();//振幅谱移动到中心（频率平移）反转换amplitude = fftShift(amplitude);// 从幅度谱和相位谱重构复数矩阵reconstructFromAmplitudeAndPhase(amplitude, phase, transformed);// 计算逆变换Eigen::MatrixXcd reconstructed = FFT2D(transformed, true);return reconstructed.real().block(0, 0, rows, cols);
}
//创建低频  radius越小，越还原图像，越大，减少低频的部分，
Eigen::MatrixXd lowFrequency(const Eigen::MatrixXd & data,int radius)
{int rows = data.rows();int cols = data.cols();// 创建低通滤波器（圆形掩码）Eigen::MatrixXd mask = Eigen::MatrixXd::Zero(rows, cols);//中心坐标int centerRow = rows >> 1;int centerCol = cols >> 1;//找到以 radius 大小的矩形范围内//左边位置int left = centerCol - radius;//超过边界 为0left = left > 0 ? left : 0;//右边位置int right = centerCol + radius;//超过边界 为0right = right < cols ? right : cols;//上边位置int top = centerRow - radius;//超过边界 为0top = top > 0 ? top : 0;//下边位置int down = centerRow + radius;//超过边界 为0down = down < rows ? down : rows;//在正矩形内画最大的圆for (int i = top; i < down; ++i){for (int j = left; j < right; ++j){//在图像中心画圆，半径不能超过double distance = std::sqrt(std::pow(i - centerRow, 2) + std::pow(j - centerCol, 2));if (distance <= radius){mask(i, j) = 1.0;}}}return mask;}
//创建低通滤波器 -
//简单的说，就是靠近频谱图中心的低频部分给保留，远离频谱图中心的高频部分给去除掉。但是这会影响图像的清晰度。
Eigen::MatrixXd  lowPassFilter(const Eigen::MatrixXd & image, int radius = 0)
{int rows = image.rows();int cols = image.cols();//大小变为2的幂数Eigen::MatrixXd data = padToPowerOfTwo(image);// 计算傅里叶变换Eigen::MatrixXcd transformed = FFT2D(data);//获取幅度谱和相位谱Eigen::MatrixXd amplitude, phase;getAmplitudeAndPhaseSpectra(transformed, amplitude, phase);//振幅谱移动到中心（频率平移）amplitude = fftShift(amplitude);///
//    //增强振幅，用于观测 --  实际运算注释掉
//    Eigen::MatrixXd amplitude1 = enhanceSpectrumDisplay(amplitude,1);
//    cv::Mat highP = tMat(amplitude1);
//    cv::imshow("lowPassFilter",highP);/////根据radius创建高频掩码Eigen::MatrixXd mask = lowFrequency(amplitude, radius);amplitude = amplitude.array() * mask.array();//振幅谱移动到中心（频率平移）反转换amplitude = fftShift(amplitude);// 从幅度谱和相位谱重构复数矩阵reconstructFromAmplitudeAndPhase(amplitude, phase, transformed);// 计算逆变换Eigen::MatrixXcd reconstructed = FFT2D(transformed, true);return reconstructed.real().block(0, 0, rows, cols);
}int main()
{cv::Mat img = cv::imread("193560523230866.png");if (img.empty()){std::cout << "请确定是否输入正确的图像文件" << std::endl;}cv::Mat gray;cvtColor(img, gray, cv::COLOR_BGR2GRAY);//图像转换CV_32F储存gray.convertTo(gray, CV_32F, 1 / 255.0, 0);//图像太大，用直接计算 耗时太长，缩小比例//resize(gray, gray, cv::Size(80, 80));//Mat  转 MatrixXdEigen::MatrixXd image = tMatrixXd(gray);// 记录开始时间auto start = std::chrono::high_resolution_clock::now();//低频Eigen::MatrixXd low = lowPassFilter(image,50);//高频Eigen::MatrixXd high = highPassFilter(image, 50);// 记录结束时间auto stop = std::chrono::high_resolution_clock::now();// 计算持续时间auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);qDebug() << "代码运行时长: " << duration.count() << " 微秒" ;cv::Mat lowP = tMat(low);//+0.5 增加显示效果high = high.array() + 0.5;cv::Mat highP = tMat(high);cv::imshow("lowP",lowP);cv::imshow("highP",highP);return 0;
}

• 振幅增强