FFTW图像处理之频域滤波和频域分析

一、 频域滤波实例

1、低通滤波（模糊图像/去噪）

#include <opencv2/opencv.hpp>
#include <fftw3.h>
#include <cmath>

void lowpass_filter(const char* input_path, const char* output_path, int cutoff_radius) {
    // 使用 OpenCV 加载图像并转换为灰度格式
    cv::Mat image = cv::imread(input_path, cv::IMREAD_GRAYSCALE);
    if (image.empty()) {
        std::cerr << "Error loading image" << std::endl;
        return;
    }

    int width = image.cols;
    int height = image.rows;

    // 分配 FFTW 数组
    fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);
    fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);
    fftw_complex *filtered = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);

    // 初始化输入数据
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            in[y * width + x][0] = image.at<uchar>(y, x) / 255.0; // 归一化到 [0,1]
            in[y * width + x][1] = 0;
        }
    }

    // 正向 FFT
    fftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan_forward);

    // 频域中心坐标
    int center_x = width / 2;
    int center_y = height / 2;

    // 创建低通滤波器
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            // 计算到中心的距离
            int dx = x - center_x;
            int dy = y - center_y;
            double distance = sqrt(dx*dx + dy*dy);
            
            if (distance <= cutoff_radius) {
                // 保留低频成分
                filtered[y * width + x][0] = out[y * width + x][0];
                filtered[y * width + x][1] = out[y * width + x][1];
            } else {
                // 滤除高频成分
                filtered[y * width + x][0] = 0;
                filtered[y * width + x][1] = 0;
            }
        }
    }

    // 逆 FFT
    fftw_plan plan_backward = fftw_plan_dft_2d(height, width, filtered, in, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(plan_backward);

    // 保存结果并归一化
    cv::Mat result(height, width, CV_8UC1);
    double max_val = 0;
    double min_val = 0;
    
    // 首先找到最小和最大值
    for (int i = 0; i < width * height; i++) {
        double val = in[i][0] / (width * height); // FFTW需要归一化
        if (val > max_val) max_val = val;
        if (val < min_val) min_val = val;
    }
    
    // 线性映射到[0,255]
    double range = max_val - min_val;
    if (range == 0) range = 1; // 避免除以零
    
    for (int i = 0; i < width * height; i++) {
        double val = in[i][0] / (width * height); // 归一化
        val = (val - min_val) / range * 255;
        result.data[i] = static_cast<uchar>(fmin(fmax(val, 0), 255));
    }

    cv::imwrite(output_path, result);

    // 清理资源
    fftw_destroy_plan(plan_forward);
    fftw_destroy_plan(plan_backward);
    fftw_free(in);
    fftw_free(out);
    fftw_free(filtered);
}

注意：

1. 尝试不同的cutoff_radius值，通常从图像尺寸的1/10开始尝试。

2. 对于512x512的图像，可以尝试50-100作为初始cutoff_radius

 

2、高通滤波（边缘增强）

#include <opencv2/opencv.hpp>
#include <fftw3.h>
#include <cmath>

void highpass_filter(const char* input_path, const char* output_path, int cutoff_radius) {
    // 使用 OpenCV 加载图像并转换为灰度格式
    cv::Mat image = cv::imread(input_path, cv::IMREAD_GRAYSCALE);
    if (image.empty()) {
        std::cerr << "Error loading image" << std::endl;
        return;
    }

    int width = image.cols;
    int height = image.rows;

    // 分配 FFTW 数组
    fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);
    fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);
    fftw_complex *filtered = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);

    // 初始化输入数据
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            in[y * width + x][0] = image.at<uchar>(y, x) / 255.0; // 归一化到 [0,1]
            in[y * width + x][1] = 0;
        }
    }

    // 正向 FFT
    fftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan_forward);

    // 频域中心坐标
    int center_x = width / 2;
    int center_y = height / 2;


    // 高通滤波：去除中心区域
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            int dx = x - center_x;
            int dy = y - center_y;
            double distance = sqrt(dx*dx + dy*dy);
            if (distance <= cutoff_radius) {
                //滤除高频成分
                filtered[y * width + x][0] = 0;
                filtered[y * width + x][1] = 0;
            }
            else{
                //保留低频成分
                filtered[y * width + x][0] = out[y * width + x][0];
                filtered[y * width + x][1] = out[y * width + x][1];
            }
        }
    }

    // 逆 FFT
    fftw_plan plan_backward = fftw_plan_dft_2d(height, width, filtered, in, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(plan_backward);

    // 保存结果并归一化
    cv::Mat result(height, width, CV_8UC1);
    double max_val = 0;
    double min_val = 0;

    // 首先找到最小和最大值
    for (int i = 0; i < width * height; i++) {
        double val = in[i][0] / (wi