2d deconvolution for Image Reconstruction

"Truth itself is truth, because it penetrate the limits of language and brings people into an intuitive grasp of the real world. "

-- Http://my1510.cn/article.php? Id = 69054

Recently, Sparse Coding uses the original method of Bruno olshausen to discover something that is closer to the truth. One step in the middle is to re-create the input image through the sparse response. Previously, it was directly solved using the MATLAB for loop, but the speed was very slow. So today I tried to accelerate this part.

First consider the fft2-> ifft2 method, the Code is as follows (refer to the http://www.mathworks.com/matlabcentral/newsreader/author/121671 ):

% Load imageimage = im2double(imread('./data/lena.png'));% image=image(201:240,201:240);[M N] = size(image);figure(101),subplot(1,2,1),imshow(image,[]);% Here filter should be a 7x7 patchfilter = reshape(A(:,2),[7,7]);% Get the filtered response in fft2F=fft2(image);H=fft2(filter,M,N);G=H.*F;Gnew= G./H;gnew=real(ifft2(double(Gnew)));figure(101),subplot(1,2,2),imshow(gnew,[]);

The result is normal.

However, since the Sparse Coding method does not directly use the filtered response, instead, it discretization the response and accumulates in time (it can be seen as a neuron ).

Membrane potential of increases with Spike input)

The actual implement code is

% Load imageimage = im2double(imread('./data/lena.png'));% image=image(201:240,201:240);[M N] = size(image);figure(101),subplot(1,2,1),imshow(image,[]);% Here filter should be a 7x7 patchfilter = reshape(A(:,2),[7,7]);% Get the filtered response in fft2F=fft2(image);H=fft2(filter,M,N);G=H.*F;% for ISTA or LCA algorithmg=real(ifft2(double(G)));g=g/10;g=max(abs(g)-0.01,0).*sign(g);G=fft2(g);Gnew= G./H;gnew=real(ifft2(double(Gnew)));figure(101),subplot(1,2,2),imshow(gnew,[]);

The result is shown in the middle one, which is different from the reconstruction of the response by multiplying filter by sparse.

Figure 1: source image (left), fft2-> ifft2 result (middle), sparse coding directly rebuilding result (right)

Therefore, you have to use the C + Mex method to write a 2D image deconvolution program. The code of the Mex file is as follows:

#include <stdio.h>#include <math.h>#include "mex.h"#definea_INprhs[0]/* sparse coefficient matrix */#definef_INprhs[1]/* filter matrix*/#definefb_OUT      plhs[0]/* reconstructed image */#define a_(i,n)a[(i) + (n)*ra]/* A is L x M */#define f_(i,n)f[(i) + (n)*rf]/* X is L x npats */#define fb_(i,n)fb[(i) + (n)*rfb]/* S is M x npats */static int ra;/* row number of a */static int ca;/* column number of a */static int rf;/* row number of f */static int cf;/* column number of f */static int rfb;/* row number of fb */static int cfb;/* column number of fb */// const mwSize N_DIM=2;void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){    if (mxGetNumberOfDimensions(a_IN)!=2 || mxGetNumberOfDimensions(f_IN)!=2){        mexErrMsgTxt("Both inputs should be 2D matrix.");    }        double *a, *f, *fb;    a = mxGetPr(a_IN);    f = mxGetPr(f_IN);    ra = (int)mxGetM(a_IN);    ca = (int)mxGetN(a_IN);    rf = (int)mxGetM(f_IN);    cf = (int)mxGetN(f_IN);      rfb = ra+rf-1;    cfb = ca+cf-1;        fb_OUT = mxCreateDoubleMatrix(rfb, cfb, mxREAL);    fb = mxGetPr(fb_OUT);    //printf("a: %d, %d\n f: %d, %d\n fb: %d, %d\n", ra,ca,rf,cf,rfb,cfb);    //initialize fb to zero;        for (int i=0;i<rfb;i++){        for (int j=0;j<cfb;j++){            fb_(i,j)=0;        }    }        for (int i=0;i<ra;i++){        for (int j=0;j<ca;j++){            if (a_(i,j)!=0){                for (int i1=0;i1<rf;i1++){                    for (int j1=0;j1<cf;j1++){                        fb_(i+i1,j+j1)+=a_(i,j)*f_(i1,j1);                    }                }            }        }    }}

Name this file deconv2.cpp and compile it in MATLAB. It is the 2D deconvolution function. The running speed is about 1000 times faster than the original Matlab code.

Note: It turns out that implement is unnecessary. in MATLAB, conv2 can implement this function, but to get the correct result, you must first reverse the filter row and column.