基于C++ Eigen FFT实现二维图像反卷积的技术求助

嘿，很高兴看到你已经摸透了傅里叶域反卷积的核心逻辑——没错，卷积在频域确实对应逐元素乘积，所以用频域除法还原原始图像O的思路完全站得住脚！不过这里有几个容易踩的小坑，我结合Eigen FFT的使用给你梳理清楚，再附上可直接参考的代码示例：

• 卷积类型对齐：要先确认你做正向卷积时用的是循环卷积还是线性卷积。如果是线性卷积，直接用原尺寸n×m做FFT相除会出现边界伪影，得先把I和K补零扩展到(n+n-1)×(m+m-1)的尺寸再处理；如果本身就是循环卷积场景，原尺寸就没问题。
• 数值稳定性：K的傅里叶变换里大概率会有接近0的数值，直接相除会导致结果爆炸。一定要加个极小的正则化项（比如1e-8），避免除以零的情况。
• Eigen FFT细节：Eigen的FFT模块处理2D数据时，要注意用复数矩阵存储变换结果，而且默认不会自动归一化，逆变换后记得手动调整幅值。

假设你已经用Eigen::MatrixXd存储了灰度图像I和核K，下面是一步步的代码：

1. 转成复数矩阵（FFT需要复数输入）

// 初始化和原图像同尺寸的复数矩阵
Eigen::MatrixXcd F_I(I.rows(), I.cols());
Eigen::MatrixXcd F_K(K.rows(), K.cols());

// 填充实部，虚部初始化为0
F_I.real() = I;
F_I.imag().setZero();
F_K.real() = K;
F_K.imag().setZero();

2. 执行2D正向傅里叶变换

Eigen::FFT<double> fft;
fft.fwd2(F_I, F_I);  // 对I做FFT，结果存回F_I
fft.fwd2(F_K, F_K);  // 对K做FFT，结果存回F_K

3. 频域除法（带正则化）

const double eps = 1e-8; // 防止除以零的小常数
Eigen::MatrixXcd F_O = F_I.array() / (F_K.array() + eps);

4. 逆傅里叶变换还原图像

Eigen::MatrixXcd O_complex(I.rows(), I.cols());
fft.inv2(O_complex, F_O);  // 执行逆FFT

// 提取实部作为最终图像（虚部是数值误差，可忽略）
Eigen::MatrixXd O = O_complex.real();

// 手动归一化（Eigen FFT默认不做归一化）
O /= I.rows() * I.cols();

• 线性卷积场景补零：如果你的原始卷积是线性卷积，先给I和K补零：

int new_rows = I.rows() + K.rows() - 1;
int new_cols = I.cols() + K.cols() - 1;
Eigen::MatrixXd I_padded = Eigen::MatrixXd::Zero(new_rows, new_cols);
Eigen::MatrixXd K_padded = Eigen::MatrixXd::Zero(new_rows, new_cols);

// 将原图像和核复制到补零矩阵的左上角
I_padded.block(0, 0, I.rows(), I.cols()) = I;
K_padded.block(0, 0, K.rows(), K.cols()) = K;

之后用补零后的矩阵执行上面的FFT流程即可。

• 噪声鲁棒性提升：如果图像I带有噪声，直接除法会放大噪声。可以换成维纳滤波的形式：

Eigen::MatrixXcd F_O = (F_I.array() * F_K.array().conjugate()) / (F_K.array().abs2() + eps);

这种方式对噪声的容忍度会高很多。

内容的提问来源于stack exchange，提问作者ABS