基于Matlab的含噪图像恢复：已知幅度相位谱的PSF与噪声处理困惑

问题概述 #

已知目标图像的中心化且经过压缩的幅度谱与相位谱，通过逆傅里叶变换初步恢复图像后，在PSF（点扩散函数）确定、噪声处理及后续操作流程上存在瓶颈。当前实现的Matlab代码如下：

mag = load('esp_mag_g5.txt');   % 加载幅度谱
fase = load('esp_fase_g5.txt'); % 加载相位谱 
figure(1)
imshow(mag,[]);
mag1=exp(mag)-1; % 还原压缩的幅度谱
F = mag1 .* exp(1j * fase); % 构造复数频谱
F_shifted = ifftshift(F); % 逆中心化
imagen_restaurada = ifft2(F_shifted); % 逆傅里叶变换
imagen_restaurada = real(imagen_restaurada); % 取实部得到图像
figure(2);
imshow(imagen_restaurada, []);
title('Imagen Restaurada (Inversa de Fourier)');
LEN = 15;      % 假设的模糊长度
THETA = 106;    % 假设的模糊角度
PSF = fspecial('motion', LEN, THETA); % 构造运动模糊PSF
% 参数K调整噪声/锐度平衡
%K = 0:0.001:0.1;
K = 0.01;

% 应用Wiener滤波器
imagen_restaurada_wiener = deconvwnr(imagen_restaurada, PSF, K);

% 显示Wiener滤波后的图像
figure;
imshow(imagen_restaurada_wiener, []);
title(['Imagen Restaurada (Filtro Wiener, K = ', num2str(K), ')']);

针对性解决方案 #

一、PSF的合理确定 #

不要直接假设运动模糊类型的PSF，需结合数据特征或盲估计方法：

• 盲反卷积估计：使用deconvblind函数，从初始PSF开始迭代估计真实PSF与清晰图像，无需先验模糊信息：

  init_psf = ones(3,3); % 初始化小尺寸PSF
  [restored_img, estimated_psf] = deconvblind(imagen_restaurada, init_psf);
  

• 频谱特征推断：观察幅度谱的特征：运动模糊的频谱会沿模糊方向出现周期性条纹；高斯模糊的频谱呈径向衰减。根据这些特征构造对应PSF后，用deconvwnr验证效果。
• 参考区域匹配：若图像存在局部清晰区域（或同场景清晰参考图），可通过模糊区域与清晰区域的傅里叶比值计算PSF。

二、噪声处理优化 #

Wiener滤波器的参数K需根据实际噪声水平调整，而非固定值：

• 基于图像区域计算K：从图像平坦区域（如背景）估计噪声方差，结合全局信号方差计算K：

  % 选取图像左上角平坦区域估计噪声
  noise_region = imagen_restaurada(1:20,1:20);
  noise_var = var(noise_region(:));
  signal_var = var(imagen_restaurada(:));
  K = noise_var / signal_var; % 噪声/信号功率比
  

• 迭代筛选最优K：遍历K∈[0.001, 0.1]，通过图像清晰度指标（如梯度幅值总和、信息熵）选择最优值。
• 预处理去噪：反卷积前先对初步恢复的图像做去噪处理，抑制噪声干扰：

  imagen_denoised = wiener2(imagen_restaurada, [5 5]); % 自适应Wiener去噪
  

三、后续操作流程优化 #

1. 频谱还原验证：确认mag1=exp(mag)-1是否匹配实际的幅度谱压缩方式。如果是线性归一化压缩，需改为mag1 = mag * (max_mag - min_mag) + min_mag（需获取原始幅度的极值）。
2. 逆变换正确性验证：可尝试另一种逆变换流程（先ifft2再fftshift），对比结果是否一致，确保中心化处理正确：

   imagen_restaurada_alt = fftshift(ifft2(F));
   imagen_restaurada_alt = real(imagen_restaurada_alt);
   

3. 振铃效应抑制：反卷积后易出现振铃效应，用edgetaper函数弱化边缘振铃：

   imagen_restaurada_wiener = edgetaper(imagen_restaurada_wiener, PSF);
   

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