#就很固执#对上一篇的退化图像的最小二乘方滤波

最小二乘方算法写得头疼，头绪很足但是效果不尽人意，希望记录下来日后常进行思考。

---

还是只谈干货

原理简单一提

---

退化模型的得出是从基本原理推导出一个数学模型的结果，假设图像在平面进行运动，通过对时间间隔内瞬时曝光量进行积分得到曝光总数，再通过傅里叶变换求出传递函数H。

约束最小二乘方滤波的核心是H对噪声的敏感问题，通过建立约束条件形成方程并求解得出最佳滤波方案，并可以迭代优化求解最优参数。

---

想做到的事情是：

• 利用所给H（u，v）公式和原图手动做出a=0.1，b=0.1，T=1时的运动退化图像（不使用自带fspecial函数）
• 对运动退化图像加一定的高斯噪声，并利用约束最小二乘方法进行滤波处理，目的得到“原图”。
• 在进行滤波过程中，使用手动设定gama和迭代求出最优gama两种方法。

  ---

  clc,clear;
  [I,map]=imread('im.bmp'); % 原始图像X数据类型为uint8
  f=double(rgb2gray(I));
  [M,N]=size(f);
  % imshow(X,[]);title('ImageOrg.bmp');
  H = zeros(M,N);
  T=1;
  a=0.1;
  b=0.1;
  for i = 1:M
     for k = 1:N
         ii=i-M/2+0.000001;
         kk=k-N/2+0.000001;
         H(i,k) = (T/(pi*(ii*a+kk*b))) * sin(pi*(ii*a+kk*b))*exp(-j*pi*(ii*a+kk*b));
     end
  end
  F=fft2(f);
  F=fftshift(F);
  G=F.*H;
  g1=ifftshift(G);
  gg=ifft2(g1);
  ggg=abs(gg)/200;
  %gg=abs(g)/3;
  %imshow(ggg);title('result.bmp');
  %imwrite(ggg,'ims.bmp')
  %gsnoise
  noise_mean = 0;% 噪声均值
  noise_var = 0.00001;% 噪声方差
  GSnoisy = imnoise(ggg,'gaussian',...
      noise_mean, noise_var);% 加噪的运动模糊图像
  %%
  %约束最小二乘滤波
  %手动给gama赋值
  g=GSnoisy;
  %[g,map]=imread('ims.bmp');
  %g=double(rgb2gray(g));
  [M,N]=size(g);
  gama=0.00006;                  %手动赋值gama
  p=[0,-1,0;-1,4,-1;0,-1,0];   %p时域
  P=fft2(p,M,N);               %P频域
  G=fft2(g);
  Gs=fftshift(G);
  Fs=(conj(H).*Gs)./(abs(H).^2+gama*(abs(P).^2));
  f=ifft2(ifftshift(Fs));      %滤波结果(手动)
  figure(1)
  imshow((f));
  figure(2)
  imshow(ggg);
  
  
  %%
  %约束最小二乘滤波
  %自动迭代选取最优gama
  g=GSnoisy;
  G=fft2(g);
  Gs=fftshift(G)       %再写一遍防止出现混乱
  gap=1e-6;
  gama=1e-6;
  a=0.5;
  yita=M*N*(noise_mean^2+noise_var);
  iteration=0;
  while (true)
      iteration=iteration+1;
      if iteration>1000
          break;
      end
      Fs=(conj(H).*Gs)./(abs(H).^2+gama*(abs(P).^2));
      R=Gs-H.*Fs;
      r=ifftshift(R);
      r_=ifft2(r);
      rr=sum(sum(r_.^2));
      if rr<yita+a
          gama=gama+gap;
      elseif rr>yita+a
          gama=gama-gap;
      else
          gama=gama;
          break
      end
  end
  f=ifft2(ifftshift(Fs)); 
  figure(3)
  imshow(f);

  如果使用自带fspecial的H函数得到效果甚好，我猜是约束最小二乘方方法非常贴合fs函数的原理而不一定贴合冈萨雷斯里按照时域里积分对应到频域进行变换的方法，但是其运动退化去除效果还是很明显的，总之还是蛮玄学的，值得推敲的一个算法。

• 下图是本文的效果图：

• 退化：

• 

• 自己的H滤波（手动gama）：

• 

• 自动gama：

• 

• 使用fspecial，滤波方法相同：

• 

•