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;

}

See also

1. Deblurring with OpenCV: Wiener filter reloaded

2. Dabbling with Wiener filter using OpenCV

3.Experiments with deblurring using OpenCV

4. De-blurring revisited with Wiener filter using OpenCV

You may also like

1. What’s up Watson? Using IBM Watson’s QAAPI with Bluemix, NodeExpress – Part 1

2. Bend it like Bluemix, MongoDB with autoscaling – Part 1

3. Informed choices through Machine Learning : Analyzing Kohli, Tendulkar and Dravid

4. A crime map of India in R: Crimes against women