FFTW图像处理入门

FFTW (Fastest Fourier Transform in the West) 是一个高效的离散傅里叶变换(DFT)库，特别适合图像处理应用。本指南将带你从零开始学习如何使用FFTW进行基本的图像处理操作。

1. 安装FFTW

Linux (Ubuntu/Debian)

bash

sudo apt-get install libfftw3-dev libfftw3-doc

MacOS (使用Homebrew)

bash

brew install fftw

Windows

1. 从官网下载预编译库: FFTW Installation on Windows

2. 或将源代码编译为DLL使用

2. 基本概念

• 傅里叶变换：将图像从空间域转换到频率域

• 频域表示：图像中的低频对应整体形状，高频对应细节和噪声

• 复数表示：FFTW使用复数数组表示频域数据

3. 第一个FFTW图像处理程序

以下是一个简单的程序，演示如何对图像进行傅里叶变换和逆变换：

c

#include <fftw3.h>
#include <stdlib.h>
#include <math.h>#define WIDTH 512
#define HEIGHT 512int main() {// 1. 分配输入输出数组fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * WIDTH * HEIGHT);fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * WIDTH * HEIGHT);// 2. 初始化输入数据(这里用简单的图案作为示例)for (int y = 0; y < HEIGHT; y++) {for (int x = 0; x < WIDTH; x++) {in[y*WIDTH + x][0] = sin(0.1*x) + cos(0.05*y); // 实部in[y*WIDTH + x][1] = 0;                        // 虚部(图像数据通常为实数)}}// 3. 创建正向变换计划fftw_plan plan_forward = fftw_plan_dft_2d(HEIGHT, WIDTH, in, out, FFTW_FORWARD, FFTW_ESTIMATE);// 4. 执行正向变换fftw_execute(plan_forward);// 5. 频域处理可以在这里进行// 6. 创建逆向变换计划fftw_plan plan_backward = fftw_plan_dft_2d(HEIGHT, WIDTH, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);// 7. 执行逆向变换fftw_execute(plan_backward);// 8. 归一化(因为FFTW的逆变换不自动缩放)for (int i = 0; i < WIDTH*HEIGHT; i++) {in[i][0] /= (WIDTH * HEIGHT);in[i][1] /= (WIDTH * HEIGHT);}// 9. 清理资源fftw_destroy_plan(plan_forward);fftw_destroy_plan(plan_backward);fftw_free(in);fftw_free(out);return 0;
}

4. 实际图像处理示例

加载真实图像(使用stb_image.h)

c

#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"void process_image(const char* input_path, const char* output_path) {int width, height, channels;unsigned char *img = stbi_load(input_path, &width, &height, &channels, 1); // 加载为灰度if (!img) {printf("Error loading image\n");return;}// 分配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);// 初始化输入数据for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {in[y*width + x][0] = img[y*width + x] / 255.0; // 归一化到[0,1]in[y*width + x][1] = 0;}}// 创建和执行正向变换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;int radius = 30; // 截止频率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  > radius) { // 超出截止频率out[y*width + x][0] = 0;out[y*width + x][1] = 0;}}}// 逆变换fftw_plan plan_backward = fftw_plan_dft_2d(height, width, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);fftw_execute(plan_backward);// 保存结果unsigned char *result = (unsigned char*)malloc(width * height);for (int i = 0; i < width*height; i++) {double val = in[i][0] / (width * height); // 归一化result[i] = (unsigned char)(fmin(fmax(val * 255, 0), 255));}stbi_write_png(output_path, width, height, 1, result, width);// 清理fftw_destroy_plan(plan_forward);fftw_destroy_plan(plan_backward);fftw_free(in);fftw_free(out);free(result);stbi_image_free(img);
}

加载真实图像(使用QImage)

#include <QImage>
#include <fftw3.h>
#include <cmath>
#include <algorithm>void process_image(const QImage& inputImg, QImage& resultImg) {// 加载图像if (inputImg.isNull()) {qDebug() << "Error loading image";return;}// 转换为灰度图像QImage grayImg = inputImg.convertToFormat(QImage::Format_Grayscale8);int width = grayImg.width();int height = grayImg.height();// 分配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);// 初始化输入数据for (int y = 0; y < height; y++) {const uchar* scanLine = grayImg.scanLine(y);for (int x = 0; x < width; x++) {in[y * width + x][0] = scanLine[x] / 255.0; // 归一化到[0,1]in[y * width + x][1] = 0;}}// 创建和执行正向变换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;int radius = std::min(width, height) / 4; // 动态截止频率double filter_factor = 0.1; // 外部频率衰减因子//低通滤波for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {// 计算相对于频谱中心的距离（考虑周期性）int dx = (x + center_x) % width - center_x;int dy = (y + center_y) % height - center_y;double distance = sqrt(dx * dx + dy * dy);// 应用平滑过渡的滤波if (distance > radius) {// 高斯衰减而不是硬截断double attenuation = exp(-(distance - radius) * (distance - radius) / 100.0);out[y * width + x][0] *= (filter_factor + (1 - filter_factor) * attenuation);out[y * width + x][1] *= (filter_factor + (1 - filter_factor) * attenuation);}}}// 逆变换fftw_plan plan_backward = fftw_plan_dft_2d(height, width, out, in,FFTW_BACKWARD, FFTW_ESTIMATE);fftw_execute(plan_backward);// 归一化因子double scale = 1.0 / (width * height);// 创建结果图像resultImg = QImage(width, height, QImage::Format_Grayscale8);// 首先找到最小和最大值以便对比度拉伸double min_val = 1.0, max_val = 0.0;for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {double val = in[y * width + x][0] * scale;min_val = std::min(min_val, val);max_val = std::max(max_val, val);}}// 防止除以零if (max_val <= min_val) {max_val = min_val + 1.0;}// 应用对比度拉伸并转换为8位图像for (int y = 0; y < height; y++) {uchar* scanLine = resultImg.scanLine(y);for (int x = 0; x < width; x++) {double val = (in[y * width + x][0] * scale - min_val) / (max_val - min_val);val = (val < 0.0) ? 0.0 : ((val > 1.0) ? 1.0 : val);scanLine[x] = static_cast<uchar>(val * 255.0);}}// 清理fftw_destroy_plan(plan_forward);fftw_destroy_plan(plan_backward);fftw_free(in);fftw_free(out);
}

 

5. 常见图像处理操作

频域滤波

• 低通滤波：去除高频噪声，模糊图像

• 高通滤波：增强边缘，保留细节

• 带通滤波：保留特定频率范围

频域分析

• 幅度谱分析：sqrt(re^2 + im^2)

• 相位谱分析：atan2(im, re)

• 功率谱分析：re^2 + im^2

6. 编译命令

使用gcc编译FFTW程序：

bash

gcc -o fftw_image fftw_image.c -lfftw3 -lm -lstb

7. 学习资源

1. FFTW官方文档

2. FFTW教程