2D Deconvolution for Image Reconstruction

"真理本身之所以是真理，就在于它穿透了语言的有限性而将人带入到对真实世界的直观把握中。"

 ——http://my1510.cn/article.php?id=69054

最近在做sparse coding, 用Bruno Olshausen最原始的方法, 因此却发现了一些背后直感上更接近真理的东西。中间有一步需要通过得到的sparse响应重建输入图像，之前一直是用Matlab for循环直接解决的，但是速度奇慢，于是今天尝试了一下把这一部分加速。

首先考虑的是fft2->ifft2的方法，代码如下(参考了http://www.mathworks.com/matlabcentral/newsreader/author/121671)：

[plain] view plaincopy

1. % Load image  
2. image = im2double(imread('./data/lena.png'));  
3. % image=image(201:240,201:240);  
4. [M N] = size(image);  
5. figure(101),subplot(1,2,1),imshow(image,[]);  
6.   
7. % Here filter should be a 7x7 patch  
8. filter = reshape(A(:,2),[7,7]);  
9.   
10. % Get the filtered response in fft2  
11. F=fft2(image);  
12. H=fft2(filter,M,N);  
13. G=H.*F;  
14.   
15. Gnew= G./H;  
16. gnew=real(ifft2(double(Gnew)));  
17.   
18. figure(101),subplot(1,2,2),imshow(gnew,[]);  

出来的结果是正常。

但由于sparse coding的方法并不是直接利用滤波后响应而是将这个响应离散化并且在时间上积累(可以看成是神经元

的membrane potential随着spike的输入而增高)

我实际implement的code是

[plain] view plaincopy

1. % Load image  
2. image = im2double(imread('./data/lena.png'));  
3. % image=image(201:240,201:240);  
4. [M N] = size(image);  
5. figure(101),subplot(1,2,1),imshow(image,[]);  
6.   
7. % Here filter should be a 7x7 patch  
8. filter = reshape(A(:,2),[7,7]);  
9.   
10. % Get the filtered response in fft2  
11. F=fft2(image);  
12. H=fft2(filter,M,N);  
13. G=H.*F;  
14.   
15. % for ISTA or LCA algorithm  
16. g=real(ifft2(double(G)));  
17. g=g/10;  
18. g=max(abs(g)-0.01,0).*sign(g);  
19. G=fft2(g);  
20.   
21. Gnew= G./H;  
22. gnew=real(ifft2(double(Gnew)));  
23.   
24. figure(101),subplot(1,2,2),imshow(gnew,[]);  

结果出来的结果就如下图中间一幅所示，与用filter乘以sparse响应重建的有一定差距

图1，原图(左), fft2->ifft2结果(中), sparse coding直接重建结果(右)

所以就不得不用c+mex的方法写了一个2d image的deconvolution 程序，mex文件代码如下：

[cpp] view plaincopy

1. #include <stdio.h>  
2. #include <math.h>  
3. #include "mex.h"  
4.   
5. #define a_IN        prhs[0]     /* sparse coefficient matrix */  
6. #define f_IN        prhs[1]     /* filter matrix*/  
7. #define fb_OUT      plhs[0]     /* reconstructed image */  
8.   
9. #define a_(i,n)     a[(i) + (n)*ra]     /* A is L x M */  
10. #define f_(i,n)     f[(i) + (n)*rf]     /* X is L x npats */  
11. #define fb_(i,n)    fb[(i) + (n)*rfb]   /* S is M x npats */  
12.   
13. static int ra;          /* row number of a */  
14. static int ca;          /* column number of a */  
15. static int rf;          /* row number of f */  
16. static int cf;          /* column number of f */  
17.   
18. static int rfb;         /* row number of fb */  
19. static int cfb;         /* column number of fb */  
20. // const mwSize N_DIM=2;  
21.   
22. void mexFunction(int nlhs, mxArray *plhs[],  
23.          int nrhs, const mxArray *prhs[])  
24. {  
25.     if (mxGetNumberOfDimensions(a_IN)!=2 || mxGetNumberOfDimensions(f_IN)!=2){  
26.         mexErrMsgTxt("Both inputs should be 2D matrix.");  
27.     }  
28.       
29.     double *a, *f, *fb;  
30.     a = mxGetPr(a_IN);  
31.     f = mxGetPr(f_IN);  
32.     ra = (int)mxGetM(a_IN);  
33.     ca = (int)mxGetN(a_IN);  
34.     rf = (int)mxGetM(f_IN);  
35.     cf = (int)mxGetN(f_IN);    
36.     rfb = ra+rf-1;  
37.     cfb = ca+cf-1;  
38.       
39.     fb_OUT = mxCreateDoubleMatrix(rfb, cfb, mxREAL);  
40.     fb = mxGetPr(fb_OUT);  
41.     //printf("a: %d, %d\n f: %d, %d\n fb: %d, %d\n", ra,ca,rf,cf,rfb,cfb);  
42.     //initialize fb to zero;  
43.       
44.     for (int i=0;i<rfb;i++){  
45.         for (int j=0;j<cfb;j++){  
46.             fb_(i,j)=0;  
47.         }  
48.     }  
49.       
50.     for (int i=0;i<ra;i++){  
51.         for (int j=0;j<ca;j++){  
52.             if (a_(i,j)!=0){  
53.                 for (int i1=0;i1<rf;i1++){  
54.                     for (int j1=0;j1<cf;j1++){  
55.                         fb_(i+i1,j+j1)+=a_(i,j)*f_(i1,j1);  
56.   
57.                     }  
58.                 }  
59.             }  
60.   
61.         }  
62.     }  
63.   
64. }  

把这个文件命名成deconv2.cpp 再在matlab里面mex deconv2.cpp编译一下，就是2d deconvolution函数了。运行速度比原来matlab code提高了1000倍左右

后记：后来发现原来这样implement没有必要，matlab里面conv2就能实现这个功能，只不过要得到正确结果的话要先把filter行列反转一下