Re-working the Lucy Richardson algorithm in OpenCV

Here is my latest attempt at deblurring using the Lucy-Richardson algorithm. For this I looked up the chapter on Iterative deconvolution and the Lucy Richardson algorithm in scribd.

As mentioned in my previous posts the blurred image can be represented as
We can represent the ill-posed blurring problem as
b(x,y) = i(x,y) ** k(x,y) + n(x,y)
where b(x,y) is the blurred image, i(x,y) the original image, k(x,y) the blur kernel and n(x,y) the noise function. If our estimate of the original image is good then n(x,y) = 0

Hence b(x,y) – i(x,y) ** k(x,y) = 0
If we add i(x,y) to both sides of the equation we have
i(x,y) = i(x,y) + b(x,y) – i(x,y) ** k(x,y)
This can be represented iteratively as
ik+1(x,y) = ik(x,y) + b(x,y) – ik(x,y) ** k(x,y) (1)

The underlined terms is the error correction.
We have to add the previous estimate with the error correction to get the new estimate.
Now we can seed this by setting ik(x,y) with the blurred image.
Hence our iteration 1 we would substitute
ik(x,y) = b(x,y) in Eqn (1)
So I have done this as follows
I have chosen a blur kernel
double a[9] = {0,40,0,0,40,0,0,40,0};

In the 1^st iteration I convolve the blurred image with the kernel
cvFilter2D(im,im_conv_kernel,&kernel1,cvPoint(-1,-1)); – A
To get the error correction I subtract with the convolved term
cvSub(im,im_conv_kernel,im_correction, 0); – B
Now I add the previous estimate with the error correction to get the new estimate
cvAdd(im,im_correction,im_new_est,NULL); – C

Finally I repeat the process
im = im_new_est;
im = cvCloneImage(im_new_est); – D
The convolved image, the error correction and the estimates of the nth iteration is shown below

The 7th,8th and 9th iteration are shown below

Note: You can clone the code from GitHub – An implementation of Lucy-Richardson algorithm in OpenCV

The complete code is given below
// deconvlucy.cpp : Defines the entry point for the console application.
//
// ===================================================================================================================================
// ========================================================Lucy-Richardson algorithm ===================================
//
// Author: Tinniam V Ganesh
// Developed 14 May 2012
// File: deconvlucy.cpp
//=====================================================================================================================================
#include “stdafx.h”
#include “math.h”
#include <cxcore.h>
#include <cv.h>
#include <highgui.h>

#define kappa 10000
int main(int argc, char ** argv)
{
IplImage* im;
IplImage* im_conv_kernel;
IplImage* im_correction;
IplImage* im_new;
IplImage* im_new_est;
IplImage* im1;

char str[80];
int i;
CvMat* cvShowDFT1(IplImage*, int, int,char*);
IplImage* cvShowInvDFT1(IplImage*, CvMat*, int, int,char*);

im1 = cvLoadImage(“kutty-1.jpg”);
cvNamedWindow(“Original-Color”, 0);
cvShowImage(“Original-Color”, im1);
im = cvLoadImage(“kutty-1.jpg”, CV_LOAD_IMAGE_GRAYSCALE );
if( !im )
return -1;

cvNamedWindow(“Original-Gray”, 0);
cvShowImage(“Original-Gray”, im);

// fk+1(x,y) = fk(x,y)

for(i=0;i < 10;i++) {

// Convolve f0(x,y)= g(x,y) with blur kernel
// f0(x,y) ** kernel

// Create a blur kernel
//double a[9]={-1,200,1,-1,200,1,-1,200,1};
//double a[9]={0,-1,0,-1,4,-1,0,-1,0};
//double a[9]={-4,40,4,-4,40,4,-4,40,4};
//double a[9]={-1,2,-1,-1,2,-1,-1,2,-1};
double a[9] = {0,40,0,0,40,0,0,40,0};
CvMat kernel1 = cvMat(3,3,CV_32FC1,a);

// Convolve the kernel with the blurred image as the seed i0(x,y) ** k(x,y)
im_conv_kernel= cvCloneImage(im);
cvFilter2D(im,im_conv_kernel,&kernel1,cvPoint(-1,-1));

cvNamedWindow(“conv”, 0);
cvShowImage(“conv”, im_conv_kernel);

// Subtract from blurred image. Error correction = b(x,y) – ik(x,y) ** k(x.y)
im_correction = cvCreateImage(cvSize(383,357),8,1);;
cvSub(im,im_conv_kernel,im_correction, 0);
cvNamedWindow(“Sub”, 0);
cvShowImage(“Sub”, im_correction);

// Add ik(x,y) with imCorrection – ik(x,y) + b(x,y) – ik(x,y) ** k(x,y)
im_new_est = cvCreateImage(cvSize(383,357),8,1);;
cvAdd(im,im_correction,im_new_est,NULL);

cvNamedWindow(“Add”, 0);
cvShowImage(“Add”, im_new_est);
sprintf(str,”Iteration – %d”,i);
cvNamedWindow(str, 0);
cvShowImage(str, im_new_est);

//Set the estimate as the previous estimate and repeat
im = im_new_est;
im = cvCloneImage(im_new_est);
}
cvWaitKey(-1);
return 0;
}