盲去卷积原理及在图像复原的应用

前言

       之前写过一篇维纳滤波在图像复原中的作用，讲述了图像退化模型以及维纳滤波的作用。维纳滤波使用的前提是知道信号和噪声的功率谱，但在实际应用中较难得到，只能根据先验知识进行估计。

       本文介绍盲去卷积复原算法，并在MATLAB中进行实验，和维纳滤波的复原效果进行一个对比。盲去卷积的方法有多种，本文主要介绍由fish提出的基于露西-理查德森(Richardson-Lucy)的盲去卷积算法。

盲去卷积原理

     露西-理查德森算法属于图像复原中的非线性算法，与维纳滤波这种较为直接的算法不同，该算法使用非线性迭代技术，在计算量、性能方面都有了一定提升。

     露西-理查德森算法是由贝叶斯公式推导而来，因为使用了条件概率（即算法考虑了信号的固有波动，因此具有复原噪声图像的能力）。贝叶斯公式如下：

      结合图像退化/复原模型，可以得到迭代函数：

  

    其中 fi 就是第i轮迭代复原图像，对应贝叶斯公式中的p(x)，g是退化函数，对应贝叶斯公式的p(y|x)，c为退化图像（c(y)dy意为在退化图像上积分），如果满足等晕条件，即图像各区域的模糊函数相同，则迭代公式可化简如下：

     这就是路西-理查德森迭代公式，其中c是退化图像，g是退化函数，f是第k轮复原图像。如果系统的退化函数PSF(g(x))已知，只要有一个初始估计f就可以进行迭代求解了。在开始迭代后，由于算法的形式，估计值会与真实值的差距迅速减小，从而后续迭代过程f的更新速度会逐渐变慢，直至收敛。算法的另一优点就是初始值f>0，后续迭代值均会保持非负性，并且能量不会发散。

     盲去卷积需要两步进行复原，原因是我们既不知道原始图像f，也不知道退化函数g。求解过程示意图如下：

      即在第k轮迭代，我们假a设原始图像已知，即k-1轮得到的fk-1，再通过R-L公式求解gk，随后，再用gk求解fk，反复迭代，最后求得最终f和g。因此，在求解最初，我们需要同时假设一个复原图像f0和一个退化函数g0。迭代公式如下：

      此外，有人采用这种盲去卷积方法进行了相关实验，下图为实验原图：

下图(a)左、右分别为附加标准差1.5%、10%泊松噪声的退化图像，图(b)左、右分别为1.5%的复原图像和PSF,图(c)为10%对应结果

     

盲去卷积MATLAB实验

     将一幅原始图像，进行模糊处理（模拟大气湍流），分别使用维纳滤波（由于没加噪声，就是逆滤波）和盲去卷积进行复原，复原结果如下：

      下面是网上找到的有关盲去卷积的MATLAB程序，可以加深理解。

[plain]  view plain  copy

1. %% Deblurring Images Using the Blind Deconvolution Algorithm   
2. %%盲反卷积算法复原图像  
3. % The Blind Deconvolution Algorithm can be used effectively when no  
4. % information about the distortion (blurring and noise) is known. The  
5. % algorithm restores the image and the point-spread function (PSF)  
6. % simultaneously. The accelerated, damped Richardson-Lucy algorithm is used  
7. % in each iteration. Additional optical system (e.g. camera)  
8. % characteristics can be used as input parameters that could help to  
9. % improve the quality of the image restoration. PSF constraints can be  
10. % passed in through a user-specified function  
11. %在不知道图像失真信息(模糊和噪声)信息情况下，盲反卷积算法可以有效地加以利用。该算法  
12. %对图像和点扩展函数（PSF）的同时进行复原。每次迭代都使用加速收敛Richardson-Lucy   
13. %算法。额外的光学系统（如照相机）的特性可作为输入参数，帮助改善图像复原质量。可以通  
14. %过用户指定的函数对PSF进行限制  
15. % Copyright 2004-2005 The MathWorks, Inc.  
16.    
17. %% Step 1: Read Image  
18. %%第一步：读取图像  
19. % The example reads in an intensity image. The |deconvblind| function can  
20. % handle arrays of any dimension.  
21. %该示例读取一个灰度图像。| deconvblind |函数可以处理任何维数组。  
22. I = imread('view.tif');  
23. figure;imshow(I);title('Original Image');  
24. %text(size(I,2),size(I,1)+15, ...  
25. %    'Image courtesy of Massachusetts Institute of Technology', ...  
26. %'FontSize',7,'HorizontalAlignment','right');    
27.      
28.   
29.    
30. %% Step 2: Simulate a Blur  
31. %%第二步：模拟一个模糊  
32. % Simulate a real-life image that could be blurred (e.g., due to camera  
33. % motion or lack of focus). The example simulates the blur by convolving a  
34. % Gaussian filter with the true image (using |imfilter|). The Gaussian filter  
35. % then represents a point-spread function, |PSF|.  
36.  %模拟一个现实中存在的模糊图像（例如，由于相机抖动或对焦不足）。这个例子通过对真实  
37. %图像进行高斯滤波器模拟图像模糊（使用|imfilter|）。高斯滤波器是一个点扩展函数，  
38. %|PSF|。  
39. PSF=fspecial('gaussian',7,10);  
40. Blurred=imfilter(I,PSF,'symmetric','conv');  %对图像I进行滤波处理；  
41. figure;imshow(Blurred);title('Blurred Image');    
42.   
43.      
44.   
45.    
46. %% Step 3: Restore the Blurred Image Using PSFs of Various Sizes  
47. %%第三步：使用不同的点扩展函数复原模糊图像  
48. % To illustrate the importance of knowing the size of the true PSF, this  
49. % example performs three restorations. Each time the PSF reconstruction  
50. % starts from a uniform array--an array of ones.  
51. %为了说明知道真实PSF的大小的重要性，这个例子执行三个修复。PSF函数重建每次都是从统一  
52. %的全一数组开始。  
53. %%  
54. % The first restoration, |J1| and |P1|, uses an undersized array, |UNDERPSF|, for  
55. % an initial guess of the PSF. The size of the UNDERPSF array is 4 pixels  
56. % shorter in each dimension than the true PSF.   
57. %第一次复原，|J1|和|P1|，使用一个较小数组，| UNDERPSF |，来对PSF的初步猜测。该  
58. %UNDERPSF数组每维比真实PSF少4个元素。  
59. UNDERPSF = ones(size(PSF)-4);  
60. [J1 P1] = deconvblind(Blurred,UNDERPSF);  
61. figure;imshow(J1);title('Deblurring with Undersized PSF');   
62.   
63.      
64.   
65. %%  
66. % The second restoration, |J2| and |P2|, uses an array of ones, |OVERPSF|, for an  
67. % initial PSF that is 4 pixels longer in each dimension than the true PSF.  
68. %第二次复原，|J2|和|P2|，使用一个元素全为1的数组，| OVERPSF|，初始PSF每维比真  
69. %实PSF多4个元素。  
70. OVERPSF = padarray(UNDERPSF,[4 4],'replicate','both');  
71. [J2 P2] = deconvblind(Blurred,OVERPSF);  
72. figure;imshow(J2);title('Deblurring with Oversized PSF');    
73.      
74.   
75.    
76. %%  
77. % The third restoration, |J3| and |P3|, uses an array of ones, |INITPSF|, for an  
78. % initial PSF that is exactly of the same size as the true PSF.  
79. %第三次复原，|J3|和|P3|，使用一个全为一的数组| INITPSF |作为初次PSF，每维与真正  
80. %的PSF相同。  
81. INITPSF = padarray(UNDERPSF,[2 2],'replicate','both');  
82. [J3 P3] = deconvblind(Blurred,INITPSF);  
83. figure;imshow(J3);title('Deblurring with INITPSF');    
84.   
85.      
86.   
87.    
88. %% Step 4: Analyzing the Restored PSF  
89. %%第四步：分析复原函数PSF  
90. % All three restorations also produce a PSF. The following pictures show  
91. % how the analysis of the reconstructed PSF might help in guessing the  
92. % right size for the initial PSF. In the true PSF, a Gaussian filter, the  
93. % maximum values are at the center (white) and diminish at the borders (black).  
94. %所有这三个复原也产生PSF。以下图片显示对PSF重建分析的如何可能有助于猜测最初PSF的大  
95. %小。在真正的PSF中，高斯滤波器的最高值在中心（白），到边界消失（黑）。  
96. figure;  
97. subplot(221);imshow(PSF,[],'InitialMagnification','fit');  
98. title('True PSF');  
99. subplot(222);imshow(P1,[],'InitialMagnification','fit');  
100. title('Reconstructed Undersized PSF');  
101. subplot(223);imshow(P2,[],'InitialMagnification','fit');  
102. title('Reconstructed Oversized PSF');  
103. subplot(224);imshow(P3,[],'InitialMagnification','fit');  
104. title('Reconstructed true PSF');    
105.   
106.      
107.   
108.    
109. %%   
110. % The PSF reconstructed in the first restoration, |P1|, obviously does not  
111. % fit into the constrained size. It has a strong signal variation at the  
112. % borders. The corresponding image, |J1|, does not show any improved clarity  
113. % vs. the blurred image,.  
114.  %第一次复原的PSF，|P1|，显然不适合大小的限制。它在边界有一个强烈的变化信号。  
115. %相应的图片|J1|，与模糊图像|Blurred|比没有表现出清晰度提高。  
116. %%  
117. % The PSF reconstructed in the second restoration, |P2|, becomes very smooth  
118. % at the edges. This implies that the restoration can handle a PSF of a  
119. % smaller size. The corresponding image, |J2|, shows some deblurring but it  
120. % is strongly corrupted by the ringing.  
121.  %第二次复原的PSF，|P2|，边缘变得非常平滑。这意味着复原可以处理一个更细致的  
122. %PSF。相应的图片|J2|，显得清晰了，但被一些“振铃”强烈破坏。  
123. %%  
124. % Finally, the PSF reconstructed in the third restoration, |P3|, is somewhat  
125. % intermediate between |P1| and |P2|. The array, |P3|, resembles the true PSF  
126. % very well. The corresponding image, |J3|, shows significant improvement;  
127. % however it is still corrupted by the ringing.  
128.  %最后，第三次复原的PSF，|P3|，介于|P1|和|P2|之间。该阵列|P3|，非常接近真  
129. %正的PSF。相应的图片，|J3|，显示了显着改善，但它仍然被一些“振铃”破坏。  
130.   
131.   
132. %% Step 5: Improving the Restoration  
133. %%第五步：改善图像复原  
134. % The ringing in the restored image, |J3|, occurs along the areas of sharp  
135. % intensity contrast in the image and along the image borders. This example  
136. % shows how to reduce the ringing effect by specifying a weighting  
137. % function. The algorithm weights each pixel according to the |WEIGHT| array  
138. % while restoring the image and the PSF. In our example, we start by  
139. % finding the "sharp" pixels using the edge function. By trial and error,  
140. % we determine that a desirable threshold level is 0.3.  
141. %在复原图像|J3|内部灰度对比鲜明的地方和图像边界都出现了“振铃”。这个例子说明了如何  
142. %通过定义一个加权函数来减少图像中的“振铃”。该算法是在对图像和PSF进行复原时，对每个  
143. %像元根据|WEIGHT|数组进行加权计算。在我们的例子，我们从用边缘函数查找“鲜明”像元  
144. %开始。通过反复试验，我们确定理想的阈值为0.3。  
145.   
146. %WEIGHT = edge(I,'sobel',.3);    
147.  WEIGHT = edge(Blurred,'sobel',.3);   
148. %%  
149. % To widen the area, we use |imdilate| and pass in a structuring element, |se|.  
150. %为了拓宽领域，我们使用|imdilate|并传递一个结构元素|se|。  
151. se = strel('disk',2);  
152. WEIGHT = 1-double(imdilate(WEIGHT,se));    
153.    
154. %%  
155. % The pixels close to the borders are also assigned the value 0.  
156. %在边界附近像素的值也被分配为0。  
157. WEIGHT([1:3 end-[0:2]],:) = 0;  
158. WEIGHT(:,[1:3 end-[0:2]]) = 0;  
159. figure;imshow(WEIGHT);title('Weight array');    
160.   
161.      
162.   
163.  %%  
164. % The image is restored by calling deconvblind with the |WEIGHT| array and an  
165. % increased number of iterations (30). Almost all the ringing is suppressed.  
166. %该图像通过|WEIGHT|数组和增加重复次数（30）调用deconvblind函数来复原。几乎所  
167. %有的“振铃”被抑制。  
168. [J P] = deconvblind(Blurred,INITPSF,30,[],WEIGHT);  
169. figure;imshow(J);title('Deblurred Image');    
170.      
171.   
172.    
173. %% Step 6: Using Additional Constraints on the PSF Restoration  
174. %第六步：使用附加约束对PSF复原  
175. % The example shows how you can specify additional constraints on the PSF.  
176. %这个例子说明了如何在PSF上指定额外的限制。  
177. % The function, |FUN|, below returns a modified PSF array which deconvblind  
178. % uses for the next iteration.   
179. %函数|FUN|返还一个修改了的PSF数组，用作deconvblind函数的下一次重复。  
180. % In this example, |FUN| modifies the PSF by cropping it by |P1| and |P2| number  
181. % of pixels in each dimension, and then padding the array back to its  
182. % original size with zeros. This operation does not change the values in  
183. % the center of the PSF, but effectively reduces the PSF size by |2*P1| and  
184. % |2*P2| pixels.   
185. %在这个例子中，通过对PSF数组各维数剪切|P1|和|P2|个值实现对PSF的修改，对数组填充  
186. %回零。此操作不会改变在PSF中心的值，而且有效地在各维减少了|2*P1|和| 2*P2|元  
187. %素。  
188.    
189. P1 = 2;  
190. P2 = 2;  
191. FUN = @(PSF) padarray(PSF(P1+1:end-P1,P2+1:end-P2),[P1 P2]);    
192.    
193. %%  
194. % The anonymous function, |FUN|, is passed into |deconvblind| last.  
195. %该匿名函数|FUN|，最后传递给| deconvblind |。  
196. %%  
197. % In this example, the size of the initial PSF, |OVERPSF|, is 4 pixels larger  
198. % than the true PSF. Setting P1=2 and P2=2 as parameters in |FUN|  
199. % effectively makes the valuable space in |OVERPSF| the same size as the true  
200. % PSF. Therefore, the outcome, |JF| and |PF|, is similar to the result of  
201. % deconvolution with the right sized PSF and no |FUN| call, |J| and |P|, from  
202. % step 4.  
203. %在这个例子中，初始PSF，|OVERPSF|，每维比真正的PSF多4个像素，。设置P1=2和P2=2作  
204. %为|FUN|的参数，可有效地使|OVERPSF|与真正的PSF的大小相同。因此，得到的结果|JF|  
205. %和|PF|，与第四步不使用|FUN|而仅用正确尺寸PSF盲反卷积得到的结果|J|和|P|类似。  
206. [JF PF] = deconvblind(Blurred,OVERPSF,30,[],WEIGHT,FUN);  
207. figure;imshow(JF);title('Deblurred Image');    
208.   
209.      
210.   
211.    
212. %%  
213. % If we had used the oversized initial PSF, |OVERPSF|, without the  
214. % constraining function, |FUN|, the resulting image would be similar to the  
215. % unsatisfactory result, |J2|, achieved in Step 3.  
216. %  
217. % Note, that any unspecified parameters before |FUN| can be omitted, such as  
218. % |DAMPAR| and |READOUT| in this example, without requiring a place holder,  
219. % ([]).  
220.  %如果我们使用了没有约束的函数|FUN|的较大的初始PSF，| OVERPSF |，所得图像将类  
221. %似第3步得到的效果并不理想的|J2|。   
222. %注意，任何在|FUN|之前未指定参数都可以省略，如|DAMPAR|和|READOUT|在这个例子中，而不需要指示他们的位置，([])。  
223.    
224. displayEndOfDemoMessage(mfilename)