约束最小平方滤波

Halcom

约束最小平方滤波
% 约束最小平方滤波
Xd = im2double(im);
HSIZE = [3 3];   % 模板窗口大小
SIGMA = 0.5;     % 标准差
H = fspecial('gaussian',HSIZE,SIGMA);
noise_power = noise_var * prod(size(Xd));   % prod(size(Xd))=65536;噪声的功率
[Zd, LAGRA] = deconvreg_filter(im,H,noise_power);  % 应用约束最小平方滤波

1. function [J, LAGRA] = deconvreg_filter(I,PSF,NP)
   
2. % 约束最小平方滤波器
   
3. %函数输入:
   
4. %        I:   输入的二维图像矩阵
   
5. %        PSF: 退化函数的空域模板
   
6. %        NP:  表示噪声的功率
   
7. %函数输出:
   
8. %        J: 约束最小平方滤波图像
   
9. %        LAGRA: 可以为一个数值，表示指定约束最小平方的最佳复原参数y，
   
10. %               也可以为[min,max]形式，表示参数y的搜索范围
    
11. %               若此参数省略，则表示搜索范围为[1e-9,1e9] 。
    
12. 
    
13. % 约束最小平方滤波
    
14. if ~isa(I,'double')
    
15.     I = double(I)/255;
    
16. end
    
17. LR = [1e-9 1e9];    % 复原参数搜索范围
    
18. % 求解输入图像维数
    
19. sizeI = size(I);
    
20. % PSF 矩阵
    
21. sizePSF = size(PSF);
    
22. % 输入图像的维数求解
    
23. numNSdim = find(sizePSF~=1);
    
24. NSD = length(numNSdim);
    
25. % 转化PSF函数到期望的维数 光传递函数OTF
    
26. otf = psf2otf(PSF,sizeI);
    
27. % regop：通过计算拉普拉斯算子计算图像的平滑性
    
28. % 具体见表达式（10.23）
    
29. if NSD == 1,
    
30.     regop = [1 -2 1];
    
31. else % 二维矩阵
    
32.     % 3x3 Laplacian matrixes
    
33.     regop = repmat(zeros(3),[1 1 3*ones(1,NSD-2)]);
    
34.     % center matrix of Laplacian
    
35.     idx = [{':'}, {':'}, repmat({2},[1 NSD-2])];
    
36.     regop(idx{:}) = [0 1 0;1 -NSD*2 1;0 1 0];  %　模板算子
    
37. end
    
38. %   regop =
    
39. %      0     1     0
    
40. %      1    -4     1
    
41. %      0     1     0
    
42. % 改变REGOP折返回原始维数
    
43. idx1 = repmat({1},[1 length(sizePSF)]);
    
44. idx1(numNSdim) = repmat({':'},[1 NSD]);
    
45. REGOP(idx1{:}) = regop;
    
46. % 转化PSF函数到期望的维数 光传递函数OTF
    
47. REGOP = psf2otf(REGOP,sizeI);
    
48. 
    
49. fftnI = fftn(I);
    
50. R2 = abs(REGOP).^2;
    
51. H2 = abs(otf).^2;
    
52. 
    
53. % 计算LAGRA值
    
54. if (numel(LR)==1) || isequal(diff(LR),0),% LAGRA is given
    
55.     LAGRA = LR(1);
    
56. else % 采用fminbnd在[1e-9,1e9]优化，加速计算
    
57.     R4G2 = (R2.*abs(fftnI)).^2;
    
58.     H4 = H2.^2;
    
59.     R4 = R2.^2;
    
60.     H2R22 = 2*H2.*R2;
    
61.     ScaledNP = NP*prod(sizeI);
    
62.     LAGRA = fminbnd(@ResOffset,LR(1),LR(2),[],R4G2,H4,R4,H2R22,ScaledNP);
    
63. end;
    
64. 
    
65. % 重构图像
    
66. Denom = H2 + LAGRA*R2;
    
67. Nomin = conj(otf).*fftnI;
    
68. 
    
69. % 保证Denom中的最小值取为sqrt(eps)
    
70. Denom = max(Denom, sqrt(eps));
    
71. J = real(ifftn(Nomin./Denom));   % 取实部
    
72. end
    
73. 
    
74. function f = ResOffset(LAGRA,R4G2,H4,R4,H2R22,ScaledNP)
    
75. % 计算反向卷积残差-留数计算
    
76. % Parseval's theorem
    
77. Residuals = R4G2./(H4 + LAGRA*H2R22 + LAGRA^2*R4 + sqrt(eps));
    
78. f = abs(LAGRA^2*sum(Residuals(:)) - ScaledNP);
    
79. end

复制代码