【图像处理】MATLAB：退化函数及多种复原方法

退化函数建模

  通过点扩散函数 PSF 进行图像复原的实验，添加适当的噪声。

  fspecial函数中，’motion’指的是运动模糊算子。

代码示例

% checkerboard产生测试板图像，第一个参数是每个正方形一边的像素数，第二个参数行数，第三为列数（缺省则等于行数）
f = checkerboard(8);                                     % 产生一个一面为8个正方形的测试板
PSF = fspecial('motion',7,45);                           % 运动模糊，PSF刚好为空间滤波器
gb = imfilter(f,PSF,'circular');                         % 减少边界效应
noise = imnoise(zeros(size(f)),'gaussian',0,0.001);      % 高斯噪声
g = gb + noise;                                          % 添加高斯噪声构造退化的图像模型

subplot(2,2,1);imshow(pixeldup(f,8),[ ]);title('原图像'); % 大图像运算过慢，故选用小图像来节省时间，   
subplot(2,2,2);imshow(gb);title('运动模糊图像');          % 以显示为目的，可通过像素赋值来放大图像。      
subplot(2,2,3);imshow(noise,[ ]);title('高斯噪声图像');
subplot(2,2,4);imshow(g);title('运动模糊+高斯噪声');

运行结果

---

频域滤波降低周期噪声

  周期噪声本身表现为类似冲击的串，这种串在傅里叶频谱中通常是可见的。滤除这些成分的主要途径是通过陷波滤波。n阶巴特沃兹陷波滤波器的传递函数如下：

---

直接逆滤波

---

维纳滤波

代码示例

f = checkerboard(8);                                     % 产生一个一面为8个正方形的测试板
PSF = fspecial('motion',7,45);                           % 运动模糊，PSF刚好为空间滤波器
gb = imfilter(f,PSF,'circular');                         % 减少边界效应
noise = imnoise(zeros(size(f)),'gaussian',0,0.001);      % 高斯噪声
g = gb + noise;                                          % 添加高斯噪声构造退化的图像模型

fr1 = deconvwnr(g,PSF);                                  % 直接逆滤波

Sn = abs(fft(noise)).^2;                                 % 噪声功率谱
nA = sum(Sn(:))/prod(size(noise));                       % 平均噪声功率，prod计算数组元素的连乘积。
Sf = abs(fft2(f)).^2;                                    % 图像功率谱
fA = sum(Sf(:))/prod(size(f));                           % 平均图像功率
R = nA/fA;                                               % 噪信功率比
fr2 = deconvwnr(g,PSF,R);                                % 参数维纳滤波器

NCORR = fftshift(real(ifft2(Sn)));                       % 自相关函数
ICORR = fftshift(real(ifft2(Sf)));  
fr3 = deconvwnr(g,PSF,NCORR,ICORR);

subplot(2,2,1);imshow(g,[]);title('噪声图像'); 
subplot(2,2,2);imshow(fr1);title('逆滤波结果');          
subplot(2,2,3);imshow(fr2);title('使用常数比率的维纳滤波的结果');
subplot(2,2,4);imshow(fr3);title('使用自相关函数的维纳滤波的结果');

运行结果

---

正则滤波(约束的最小二乘方滤波)

算法另寻…

代码示例

f = checkerboard(8);                                     % 产生一个一面为8个正方形的测试板
PSF = fspecial('motion',7,45);                           % 运动模糊，PSF刚好为空间滤波器
gb = imfilter(f,PSF,'circular');                         % 减少边界效应
noise = imnoise(zeros(size(f)),'gaussian',0,0.001);      % 高斯噪声
g = gb + noise;                                          % 添加高斯噪声构造退化的图像模型

fr1 = deconvreg(g,PSF,4);                                % 正则滤波器
fr2 = deconvreg(g,PSF,0.4,[1e-7 1e7]);

subplot(1,3,1);imshow(g);title('运动模糊+高斯噪声(g)');
subplot(1,3,2);imshow(fr1);title('正则滤波器');
subplot(1,3,3);imshow(fr2);title('正则滤波器(缩小范围)');

运行结果

---

使用L-R算法的迭代非线性复原

  前面所讨论的图像复原方法都是线性的，而LR算法是非线性方法，又称Lucy-Richardson算法。

代码示例

f = checkerboard(8);                                % 产生64×64像素的方形图像
PSF = fspecial('gaussian',7,10);                    % 产生一个大小为7×7且标准偏差为10的高斯PSF
SD = 0.01;                                          % 标准偏差
g = imnoise(imfilter(f,PSF),'gaussian',0,SD^2);     % 添加均值为0、标准偏差为0.01的高斯噪声

DAMPAR = 10*SD;                                     % 结果图像与原图像的偏离阈值
LIM = ceil(size(PSF,1)/2);
WEIGHT = zeros(size(g));                            % WEIGHT数组大小64×64
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1;              % WEIGHT数组有值为0的4像素宽的边界，其余像素为1

f1 = deconvlucy(g,PSF,5,DAMPAR,WEIGHT);             % 迭代5次
f2 = deconvlucy(g,PSF,10,DAMPAR,WEIGHT);            % 迭代10次
f3 = deconvlucy(g,PSF,20,DAMPAR,WEIGHT);            % 迭代20次
f4 = deconvlucy(g,PSF,100,DAMPAR,WEIGHT);           % 迭代100次

subplot(2,3,1);imshow(pixeldup(f,8));title('原图像');
subplot(2,3,2);imshow(g);title('高斯噪声+运动模糊');    
subplot(2,3,3);imshow(f1);title('L-R算法迭代5次');   
subplot(2,3,4);imshow(f2);title('L-R算法迭代10次');  % 图像虽然改进，但依然模糊
subplot(2,3,5);imshow(f3);title('L-R算法迭代20次');  % 迭代20次为合理复原
subplot(2,3,6);imshow(f4);title('L-R算法迭代100次'); % 除了稍微清晰明亮，并无显著改进

运行结果

---

盲去卷积

  不以PSF知识为基础的图像复原方法统称为盲去卷积算法。

代码示例

PSF = fspecial('gaussian',7,10);
SD = 0.01;
g = imnoise(imfilter(f,PSF),'gaussian',0,SD^2);

INITPSF = ones(size(PSF));                                              % 点扩散函数的初始估计
DAMPAR = 10*SD;                                                         % 结果图像与原图像的偏离阈值
LIM = ceil(size(PSF,1)/2);
WEIGHT = zeros(size(g));                                                % WEIGHT数组大小64×64
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1;                                  % WEIGHT数组有值为0的4像素宽的边界，其余像素为1
[fr1,PSFe1] = deconvblind(g,INITPSF,5,DAMPAR,WEIGHT);                   % 盲去卷积迭代5次
[fr2,PSFe2] = deconvblind(g,INITPSF,10,DAMPAR,WEIGHT);                  % 盲去卷积迭代10次
[fr3,PSFe3] = deconvblind(g,INITPSF,20,DAMPAR,WEIGHT);                  % 盲去卷积迭代20次

subplot(2,2,1);imshow(pixeldup(PSF,73),[ ]);title('退化函数的PSF');        
subplot(2,2,2);imshow(pixeldup(PSFe1,73),[ ]);title('盲去卷积迭代5次');
subplot(2,2,3);imshow(pixeldup(PSFe2,73),[ ]);title('盲去卷积迭代10次');
subplot(2,2,4);imshow(pixeldup(PSFe3,73),[ ]);title('盲去卷积迭代20次'); % 最接近真正PSF

运行结果

---