图像复原原理及实现

补充知识

相关函数

• 信号$X$是随时间变化的随机变量序列，不同时间起始点的序列$X(s)$和$X(t)$的相关函数可以表示为

$$R(s,t) = \frac{{E[({X_t} - {\mu _t})({X_s} - {\mu _s})]}}{{{\sigma _t}{\sigma _s}}}$$
* 对于二阶稳态过程，即信号的均值和方差不随时间的变化而变化，则此时相关函数只是时间差 $\tau = s - t$的函数，可以简化为

$$R(\tau ) = \frac{{E[({X_t} - \mu )({X_s} - \mu )]}}{{{\sigma ^2}}}$$
这即是统计学函数里面的 自相关函数。
* 将相关函数离散化，可以得到

$$R(\tau ) = \frac{1}{N}\sum\limits_{i = 0}^{N - 1} {\frac{{{X_i} - \mu }}{\sigma }\frac{{{X_{i + \tau }} - \mu }}{\sigma }}$$
即两个向量的内积，它可以反映同一序列在不同时刻取值之间的相似程度。

标准互相关函数

• Normalized Cross-Correlation, NCC

• 两幅图像$f(x,y)$, $t(x,y)$，均值和方差分别为$({u_f},\sigma _f^2)$,$({u_t},\sigma _t^2)$，则NCC为

$$NCC = \frac{1}{{N - 1}}\sum\limits_{x,y} {\frac{{f(x,y) - {\mu _f}}}{{{\sigma _f}}}\frac{{t(x,y) - {\mu _t}}}{{{\sigma _t}}}}$$
对图像减去均值并除以标准差，可以保证图片的亮度和对比度不变性，即NCC具有在亮度和对比度下的稳定性。

功率谱密度

• power spectrum density, PSD

• 功率谱密度是信号在某个频率下拥有的能量，一个信号的功率谱密度就是该信号自相关函数的傅立叶变换，将一个信号的psd对所有频率进行积分，得到的值即为信号的总能量。功率谱密度计算方法为

$${P_{XX}}(s) = F(s){F^*}(s) = {\left| {F(s)} \right|^2}$$
其中 $F(s)$和 $F^*(s)$是 $f(t)$的傅立叶变换和 $F(s)$的复共轭。

• 每个信号的功率谱密度唯一，是信号的一种属性，在图像处理中，可以将图像看作信号，对图像的ROI进行标准化之后，再进行傅立叶变化，就可以得到图像的功率谱，在进行模板匹配等操作时，就有了一个可以匹配的模板属性，可以作为ROI的一个特征。
• FFT可以为计算psd提供一个高效的计算手段

白噪声

• 白噪声是一种功率谱密度为常数的随机信号或者随机过程，信号在各个频率上的能量相同，因此白噪声不具有周期性

参考内容

• http://skpsun.blog.163.com/blog/static/2760055201122043815725/

主要图像复原方法简介

图像复原处理的过程

• 设原始图像$f(x,y)$，则退化后的图像可以表示为

$$g(x,y)=H[f(x,y)]+\eta(x,y)$$
其中 $H$为退化函数，$\eta(x,y)$为加性噪声，复原的目标是得到基于退化图像的一个估计 $\hat f(x,y)$，使得图像能够尽可能地贴近原始图像。
* 若H是线性的、空间不变的过程，则退化图像在空间频域中可以表示为

$$g(x,y) = h(x,y) * f(x,y) + \eta(x,y)$$
进行傅立叶变换，有
$$G(u,v) = H(u,v)F(u,v) + N(u,v)$$
* 噪声模型有很多种，包括高斯噪声、椒盐噪声、泊松噪声等

直接逆滤波处理

• 用于复原一幅退化图像的最简方法是构成以下形式的估计

$$\hat F(u,v) = \frac{G(u,v)}{H(u,v)}$$
然后进行逆傅立叶变换得到原图的估计。

可以将上面的公式写为

$$\hat F(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)}$$
则即使知道了准确的 $H(u,v)$，也无法恢复 $F(u,v)$，因为还需要知道噪声信号。而且很多 $H(u,v)$为0的情况也是一个很大的问题。可以将 $H(u,v)=0$的点特殊处理一下，此时这种方法被称为 伪逆滤波。
* 很少有使用价值

维纳滤波

• 维纳滤波是一种线性图像复原方法，寻找一个使统计误差函数

$${e^2} = E\{ {(f - \hat f)^2}\}$$ 最小的估计 $\hat f$，在频域可以表示为

$${\rm{\hat F}}(u,v) = [\frac{1}{{H(u,v)}}\frac{{{{\left| {H(u,v)} \right|}^2}}}{{{{\left| {H(u,v)} \right|}^2} + {S_\eta }(u,v)/{S_f}(u,v)}}]G(u,v)$$

其中
$H(u,v)$为退化函数，$S_\eta$和$S_f$是噪声与未退化图像的功率谱。信噪比为

$$SNR=\frac{S_f(u,v)}{S_\eta(u,v)}$$
信噪比为无穷大时，表明没有噪声，则维纳滤波退化为逆滤波。

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

• 约束最小二乘滤波的核心是$H$对噪声敏感的问题。处理这个问题的一种方法是基于平滑度测量的最优性。复原过程中，需要寻找准则函数$C$的最小值。函数$C$定义为

$$C = \sum\limits_{x = 0}^{M - 1} {\sum\limits_{y = 0}^{N - 1} {{{[{\nabla ^2}f(x,y)]}^2}} }$$

函数的约束条件为

$$\left\| {g - H\hat f} \right\| = {\left\| \eta \right\|^2}$$

这个最优化问题的频域解决办法由下式给出

$$\hat F(u,v) = [\frac{{{H^*}(u,v)}}{{{{\left| {H(u,v)} \right|}^2} + \gamma {{\left| {P(u,v)} \right|}^2}}}]G(u,v)$$

其中$P(u,v)$是拉普拉斯算子$p(x,y)$的傅立叶变换，有

$$P = \left[ \begin{array}{l} 0{\kern 15pt} 1{\kern 15pt} 0\\ 1{\kern 3pt} -4{\kern 15pt} 1\\ 0{\kern 15pt} 1{\kern 15pt} 0 \end{array} \right]$$

未知量有$\gamma$和${\left\| \eta \right\|^2}$两个。若已知和噪声功率（标量）成比例的${\left\| \eta \right\|^2}$，则通过迭代，$\gamma$可以得出。

LR滤波

• 是一种非线性滤波方法，也是找到模型的极大似然函数

代码及实现部分

• 主要参考的是冈萨雷斯的数字图像处理matlab版

代码

%%
clc,clear,close all
% 原始图像
f = checkerboard(8);
% 噪声滤波器
PSF = fspecial('motion', 7, 45);
% 退化图像
gb = imfilter( f, PSF, 'circular' );
% 高斯滤波
noise = imnoise( zeros(size(f)), 'gaussian', 0, 0.001 );
% 将噪声加到原图上
g = gb + noise;

figure
subplot(221), imshow( f ), title('原图')
subplot(222), imshow( noise, [] ), title('高斯噪声')
subplot(223), imshow( gb ), title('退化图像')
subplot(224), imshow( g ), title('退化图像加高斯噪声')

%% 图像复原
% 噪信比默认为0，即信噪比相当于无穷大，即直接逆滤波的过程
fr1 = deconvwnr( g, PSF );

Sn = abs( fft2(noise) ) .^2;
nA = sum( Sn(:) ) / numel( noise );
Sf = abs( fft2(f) ) .^2;
fA = sum( Sf(:) ) / numel( f );
R = nA / fA;  % 平均噪信比计算
% 使用常数比例的维纳滤波进行复原
fr2 = deconvwnr( g, PSF, R );
% 使用自相关函数的维纳滤波进行复原
Ncorr = fftshift( real(ifft2(Sn)) );
Fcorr = fftshift( real(ifft2(Sf)) );
fr3 = deconvwnr( g, PSF, Ncorr, Fcorr );

figure
subplot(221), imshow( g ), title('模糊的运动图像')
subplot(222), imshow( fr1, [] ), title('直接逆滤波复原图像')
subplot(223), imshow( fr2, [] ), title('常数比例维纳滤波复原图像')
subplot(224), imshow( fr3, [] ), title('使用自相关函数的维纳滤波复原图像')

%% 约束最小二乘滤波
% 4 约等于  64*64*0.001
fr1 = deconvreg( g, PSF, 4 );

fr2 = deconvreg( g, PSF, 0.4, [1e-7, 1e7] );

figure
subplot(221), imshow( g ), title('模糊的运动图像')
subplot(222), imshow( fr1, [] ), title('仅噪声功率参数的正则滤波')
subplot(223), imshow( fr2, [] ), title('包含噪声和gamma范围的正则滤波')

%% Lucy-Richardson算法的迭代非线性复原
g = checkerboard(8);
ori = g;
PSF = fspecial( 'gaussian', 7, 10 );
SD = 0.01;
g = imnoise( imfilter(g, PSF), 'gaussian', 0, SD^2 );

damper = 10*SD;
lim = ceil( size(PSF, 1) / 2 );
weight = zeros( size(g) );
weight( lim+1:end-lim, lim+1:end-lim ) = 1;
numit = 5;
f5 = deconvlucy( g, PSF, numit, damper, weight );
numit = 20;
f20 = deconvlucy( g, PSF, numit, damper, weight );
numit = 50;
f50 = deconvlucy( g, PSF, numit, damper, weight );
numit = 100;
f100 = deconvlucy( g, PSF, numit, damper, weight );


figure
subplot(231), imshow( ori ), title('原始图像')
subplot(232), imshow( g, [] ), title('加两次高斯噪声的图')
subplot(233), imshow( f5, [] ), title('LR 5次迭代')
subplot(234), imshow( f20 ), title('LR 20次迭代')
subplot(235), imshow( f50, [] ), title('LR 50次迭代')
subplot(236), imshow( f100,[] ), title('LR 100次迭代')

结果

• 原图+噪声

  

• 逆滤波和维纳滤波

  

• 正则滤波

  

• LR滤波

  

任意模糊运动图像的复原

运动参数估计

• 主要包括运动的角度和长度，运动的角度可以用Radon方法进行估计，运动的长度可以用倒谱法进行估计

• 参考链接:http://www.cnblogs.com/yomman/p/3424494.html

• 代码

  function [beta, len] = getMotionPara( img, showImg )
  % getMotionPara : 得到运动模糊图像的的参数估计
  %   img   : 输入的灰度图像
  %   beta  : 运动参数的角度
  %   len   : 运动参数的长度
  % 用这两个参数可以对图像进行复原
      if nargin <= 1
          showImg = false;
      end
      fg = fft2( img );
      spec = log(1+abs(fg));
      specCenter = fftshift( spec );
      thres = graythresh( specCenter );
      bw = edge( specCenter, 'canny', thres );
      th = 1:180;
      R = radon( bw, th );
      maxRadon = max(R(:));
      [~, n] = find( R == maxRadon );
      [M, N] = size( specCenter );
      beta = atan( tan(n*pi/180)*M/N)*180/pi;
  
      %利用倒谱的方法估计模糊尺度
      theta = beta;
      %――――――――――――――――――――――――――――
      %转化为倒谱域：Cep(g(x,y)) = invFT{log(FT(g(x,y)))}
      fin = fft2(img); %转化为频域
  
      lgfin = abs(log(1 + abs(fin))); %转化为对数，abs是求复数的模
      cin = ifft2(lgfin); %得到倒谱域的图像
      cinrot = imrotate(cin, -theta); %旋转图像，使模糊方向水平
  
      %计算每列的平均值
      for i=1:size(cinrot, 2)   %列数
          avg(i) = 0;
          for j=1:size(cinrot, 1)  %行数
              avg(i) = avg(i) + cinrot(j, i);
          end
          avg(i) = avg(i)/size(cinrot, 1); %第i列的平均值
      end
      avgr = real(avg);
  
      %用第一个负值计算模糊尺度
      index = 0;
      for i = 1:round(size(avg,2)),
          if real(avg(i))<0,
              index = i;
              break;
          end
      end
  
      %如果找到了模糊尺度，则：
      if index~=0
          len = index;
      else
  
      %如果没有找到模糊尺度（无负峰），则利用查找最低峰的方法再查找
          index = 1;
          startval = avg(index);
          for i = 1 : round(size(avg, 2)/2),
              if startval>avg(i),
                  startval = avg(i);
                  index = i;
              end
          end
          len = index;
      end
  
      if showImg
          subplot(221), imshow( img ), title('原图')
          subplot(222), imshow( specCenter, [] ), title('居中后的压缩频谱')
          subplot(223), imshow( bw ), title('频谱二值化图像')
          subplot(224), imshow( R ), title('Radon变换后图像')
      end
  
  end

滤波效果对比

• 一些信噪比等的参数主要是通过尝试得出

• 代码(main.m)

  %% 
  clc,clear,close all
  img = imread( '1.jpg' );
  gb = double(rgb2gray( img )) / 255;
  %% 运动参数估计
  [beta, len] = getMotionPara( gb );
  fprintf('PSF parameters, angle : %.2f, length : %.2f\n', beta, len);
  %% 图像噪声去除
  % 运动噪声滤波器估计
  PSF = fspecial('motion', len, beta);
  % 直接逆滤波求得结果
  % 逆滤波得到的结果很差，说明原图中有较多噪声，需要通过其他方法去除
  res1 = deconvwnr( gb, PSF );
  % 维纳滤波
  % 随便给的一个信噪比的值
  nsr = 0.01;
  res2 = deconvwnr( gb, PSF, nsr );
  % 约束最小二乘滤波
  % 尝试给出的一个噪声能量值
  noisePower = 10;
  res3 = deconvreg( gb, PSF, noisePower );
  
  figure
  subplot(231), imshow( img ), title('原图')
  subplot(232), imshow( res1, [] ), title('逆滤波复原后图像')
  subplot(233), imshow( res2, [] ), title('维纳后图像')
  subplot(234), imshow( res3, [] ), title('约束最小二乘复原后图像')

• 效果图

  

  

• 注：其实图像复原最重要的还是参数估计和噪声估计，不然复原效果会一直很差