Deblurring with OpenCV: Weiner filter reloaded


The problem of deblurring has really caught my fancy though I have only had partial success with it. Deblurring is basically an ill-posed problem where there are 2 unknowns namely the original image and a blurring function. There are many solutions to this problem involving a fair amount of mathematics.

Every now and then I will sneak into some white paper on this topic only to beat a hasty retreat gulping for air as I get drowned in the abstract math. Anyway my search led me to the following presentation “Deblurring in CT (Computer Tomography)” by Kriti Sen Sharma which seemed to make a lot of sense. This presentation shows how a PSF (Point Spread Function or the blur kernel) can blur an image as shown below.

Further it can be shown that the Wiener filter can be represented as
W(u) = P(u*)
——-
|P(u)|^2 + K

Where P(u) is PSD (Point Spread Distribution) of the blur kernel. K is S/N or signal to noise ratio.

For the PSF I took the Gaussian distribution given in Wikipedia – Gaussian blur given by

I tried with various values of SIGMA. The best value was SIGMA = 0.014089642. I get just one pixel with a value of > 0.  This is shown below.

I also tried various values of K. The larger the K the clearer was the DFT INVERSE. This was the best I got.

I am still nowhere near removing the blur but I will get there sometime in the future.
Any and all comments welcome.

Note: You can clone the code from GitHub: Weiner filter reloaded

I am including the code executed with Visual Studio 2010 express

//======================================================================================================================
// Wiener filter implemention using Gaussian blur kernel
// Developed by: Tinniam V Ganesh
// Date: 11 May 2012
//======================================================================================================================
#include “stdafx.h”
#include “math.h”
#include <cxcore.h>
#include <cv.h>
#include <highgui.h>

#define kappa 10000
int main(int argc, char ** argv)
{
int height,width,step,channels,depth;
uchar* data1;
CvMat *dft_A;
CvMat *dft_B;
CvMat *dft_C;
IplImage* im;
IplImage* im1;
IplImage* image_ReB;
IplImage* image_ImB;

IplImage* image_ReC;
IplImage* image_ImC;
IplImage* complex_ImC;
CvScalar val;
IplImage* k_image_hdr;
int i,j,k;

FILE *fp;
fp = fopen(“test.txt”,”w+”);
int dft_M,dft_N;
int dft_M1,dft_N1;

CvMat* cvShowDFT1(IplImage*, int, int,char*);
void 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);
IplImage* k_image;
int rowLength= 11;
long double kernels[11*11];
CvMat kernel;
int x,y;
long double PI_F=3.14159265358979;

//long double SIGMA = 0.84089642;
long double SIGMA = 0.014089642;
//long double SIGMA = 0.00184089642;
long double EPS = 2.718;
long double numerator,denominator;
long double value,value1;
long double a,b,c,d;

numerator = (pow((float)-3,2) + pow((float) 0,2))/(2*pow((float)SIGMA,2));
printf(“Numerator=%f\n”,numerator);
denominator = sqrt((float) (2 * PI_F * pow(SIGMA,2)));
printf(“denominator=%1.8f\n”,denominator);

value = (pow((float)EPS, (float)-numerator))/denominator;
printf(“Value=%1.8f\n”,value);
for(x = -5; x < 6; x++){
for (y = -5; y < 6; y++)
{
//numerator = (pow((float)x,2) + pow((float) y,2))/(2*pow((float)SIGMA,2));
numerator = (pow((float)x,2) + pow((float)y,2))/(2.0*pow(SIGMA,2));
denominator = sqrt((2.0 * 3.14159265358979 * pow(SIGMA,2)));
value = (pow(EPS,-numerator))/denominator;
printf(” %1.8f “,value);
kernels[x*rowLength +y+55] = (float)value;

}
printf(“\n”);
}
printf(“———————————\n”);
for (i=-5; i < 6; i++){
for(j=-5;j < 6;j++){
printf(” %1.8f “,kernels[i*rowLength +j+55]);
}
printf(“\n”);
}
kernel= cvMat(rowLength, // number of rows
rowLength, // number of columns
CV_32FC1, // matrix data type
&kernels);
k_image_hdr = cvCreateImageHeader( cvSize(rowLength,rowLength), IPL_DEPTH_32F,1);
k_image = cvGetImage(&kernel,k_image_hdr);

height = k_image->height;
width = k_image->width;
step = k_image->widthStep/sizeof(float);
depth = k_image->depth;
channels = k_image->nChannels;
//data1 = (float *)(k_image->imageData);
data1 = (uchar *)(k_image->imageData);
cvNamedWindow(“blur kernel”, 0);
cvShowImage(“blur kernel”, k_image);

dft_M = cvGetOptimalDFTSize( im->height – 1 );
dft_N = cvGetOptimalDFTSize( im->width – 1 );
//dft_M1 = cvGetOptimalDFTSize( im->height+99 – 1 );
//dft_N1 = cvGetOptimalDFTSize( im->width+99 – 1 );
dft_M1 = cvGetOptimalDFTSize( im->height+3 – 1 );
dft_N1 = cvGetOptimalDFTSize( im->width+3 – 1 );
printf(“dft_N1=%d,dft_M1=%d\n”,dft_N1,dft_M1);

// Perform DFT of original image
dft_A = cvShowDFT1(im, dft_M1, dft_N1,”original”);
//Perform inverse (check)
//cvShowInvDFT1(im,dft_A,dft_M1,dft_N1, “original”); – Commented as it overwrites the DFT
// Perform DFT of kernel
dft_B = cvShowDFT1(k_image,dft_M1,dft_N1,”kernel”);
//Perform inverse of kernel (check)
//cvShowInvDFT1(k_image,dft_B,dft_M1,dft_N1, “kernel”);- Commented as it overwrites the DFT
// Multiply numerator with complex conjugate
dft_C = cvCreateMat( dft_M1, dft_N1, CV_64FC2 );
printf(“%d %d %d %d\n”,dft_M,dft_N,dft_M1,dft_N1);

// Multiply DFT(blurred image) * complex conjugate of blur kernel
cvMulSpectrums(dft_A,dft_B,dft_C,CV_DXT_MUL_CONJ);
//cvShowInvDFT1(im,dft_C,dft_M1,dft_N1,”blur1?);

// Split Fourier in real and imaginary parts
image_ReC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
image_ImC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
complex_ImC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 2);
printf(“%d %d %d %d\n”, dft_M,dft_N,dft_M1,dft_N1);
//cvSplit( dft_C, image_ReC, image_ImC, 0, 0 );
cvSplit( dft_C, image_ReC, image_ImC, 0, 0 );

// Compute A^2 + B^2 of denominator or blur kernel
image_ReB = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
image_ImB = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);

// Split Real and imaginary parts
cvSplit( dft_B, image_ReB, image_ImB, 0, 0 );
cvPow( image_ReB, image_ReB, 2.0);
cvPow( image_ImB, image_ImB, 2.0);
cvAdd(image_ReB, image_ImB, image_ReB,0);
val = cvScalarAll(kappa);
cvAddS(image_ReB,val,image_ReB,0);
//Divide Numerator/A^2 + B^2
cvDiv(image_ReC, image_ReB, image_ReC, 1.0);
cvDiv(image_ImC, image_ReB, image_ImC, 1.0);

// Merge Real and complex parts
cvMerge(image_ReC, image_ImC, NULL, NULL, complex_ImC);
// Perform Inverse
cvShowInvDFT1(im, (CvMat *)complex_ImC,dft_M1,dft_N1,”Weiner o/p k=10000 SIGMA=0.014089642″);
cvWaitKey(-1);
return 0;
}

CvMat* cvShowDFT1(IplImage* im, int dft_M, int dft_N,char* src)
{
IplImage* realInput;
IplImage* imaginaryInput;
IplImage* complexInput;
CvMat* dft_A, tmp;
IplImage* image_Re;
IplImage* image_Im;
char str[80];
double m, M;
realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);
cvScale(im, realInput, 1.0, 0.0);
cvZero(imaginaryInput);
cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);

dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );
image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);

// copy A to dft_A and pad dft_A with zeros
cvGetSubRect( dft_A, &tmp, cvRect(0,0, im->width, im->height));
cvCopy( complexInput, &tmp, NULL );
if( dft_A->cols > im->width )
{
cvGetSubRect( dft_A, &tmp, cvRect(im->width,0, dft_A->cols – im->width, im->height));
cvZero( &tmp );
}
// no need to pad bottom part of dft_A with zeros because of
// use nonzero_rows parameter in cvDFT() call below

cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );
strcpy(str,”DFT -“);
strcat(str,src);
cvNamedWindow(str, 0);

// Split Fourier in real and imaginary parts
cvSplit( dft_A, image_Re, image_Im, 0, 0 );
// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re, NULL);
cvPow( image_Re, image_Re, 0.5 );

// Compute log(1 + Mag)
cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)
cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
cvShowImage(str, image_Re);
return(dft_A);
}

void cvShowInvDFT1(IplImage* im, CvMat* dft_A, int dft_M, int dft_N,char* src)
{
IplImage* realInput;
IplImage* imaginaryInput;
IplImage* complexInput;
IplImage * image_Re;
IplImage * image_Im;
double m, M;
char str[80];
realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);
image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);

//cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, complexInput->height );
cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, dft_M);
strcpy(str,”DFT INVERSE – “);
strcat(str,src);
cvNamedWindow(str, 0);
// Split Fourier in real and imaginary parts
cvSplit( dft_A, image_Re, image_Im, 0, 0 );
// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re, NULL);
cvPow( image_Re, image_Re, 0.5 );

// Compute log(1 + Mag)
cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)
cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
//cvCvtColor(image_Re, image_Re, CV_GRAY2RGBA);
cvShowImage(str, image_Re);
}

See also
– Experiments with deblurring using OpenCV
– De-blurring revisited with Wiener filter using OpenCV
–  Dabbling with Wiener filter using OpenCV
– Re-working the Lucy-Richardson Algorithm in 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

Find me on Google+

Experiments with deblurring using OpenCV


Deblurring refers to the removal of the blur from blurred images.  The removal of blur is extremely important in the fields of medical imaging, astronomy etc. In medical imaging this is also known as denoising and finds extensive applications in ultra sonic and CT images.  Similarly in astronomy there is a need to denoise and deblur images that are taken by space telescopes for e.g. the Hubble telescope.

Deblurring is really a tough problem to solve though it has been around for ages. While going through some of the white papers on deblurring I have been particularly fascinated by the results. The blurred images are restored back to its pristine, original sharp image. It is almost magical and is amazing to know that a computer program is able to detect and remove patches, almost bordering on “computer perception”.

However the solution to the problem is very involved. Almost every white paper I read in deblurring techniques involves abstruse mathematical concepts, enough to make one break into ‘cold sweat’ as I did several times. There are several well known methods of removing blur from images. The chief among them are the Richardson-Lucy method, the Weiner filter (do read my post “Dabbling with Wiener filter using OpenCV“), Radon transform or some form Bayesian filtering etc. All of these methods perform the converse of convolution known as de-convolution.  Almost all approaches are based on the following equation

b(x) = k(x) * i(x) + n(x)                           – (1)

Where b(x) – is the blurred image, i(x) is the latent (good) image, k(x) is the blur kernel also known as the Point Spread Function (PSF) and n(x) is random noise.

The key fact to note in the above equation is that the blurred image (b) is a convolution of a good latent image with a blur kernel plus some additive noise. So the good latent image is convolved by a blurring function or the points spread function and results in the blurred image.

Now there are 2 unknowns in the above equation, both the blur kernel and the latent image. To obtain the latent image a process known as de-convolution has to be performed.

If equation (1) is converted to the frequency domain using the Fourier transform on equation (1) we have

B(w) =  K(w) * I(w) + N(w)
Ignoring noise we have
I(w)  = B(w)/K(w)
or
I(w)  = DFT [B(w)]/DFT[K(w)]
Or i(t) = Inverse DFT {[DFT B(w)]/DFT[K(w)]}
If DFT[K(w)] can be represented as A+iB then the above can be represented as
i(t) = Inverse DFT { DFT [B(w)] * (A – iB)}
—————————
A**2 + B**2
where A-iB is the complex conjugate of the DFT of the blur kernel

This is known as de-convolution. There are two types of de-convolution blind and non-blind. In the non-blind de-convolution method an initial blur kernel is assumed and is used to de-blur the image. In the blind convolution an initial estimate of the blur kernel is made and the latent image is successively iterated to obtain the best image. This is based on a method known as maximum-a-posteriori (MAP) to iterate through successive estimations of better and better blur kernels. The Richardson-Lucy algorithm also tries to estimate the blur kernel from the blurred image.

In this post I  perform de-convolution using an estimated blur kernel. As can be seen there is a slight reduction of the blur. By choosing a large kernel probably a 100 x 100 matrix would give a better result.

You can clone this code from GitHub at “Deblurring by deconvolution in OpenCV

Kindly do share any thoughts, ideas that you have. I have included the full code for this. Feel free to use the code and tweak it as you see fit and please share any findings you come across.

Steps in the de-convolution process

  1. Perform DFT on blurred image. Also perform Inverse DFT to verify whether the process of DFT is correct or not. Make sure that the line for performing the inverse is commented out as it overwrites the DFT array.
           // Perform DFT of original image
		dft_A = cvShowDFT(im, dft_M1, dft_N1,"original");
		//Perform inverse (check)
		cvShowInvDFT(im,dft_A,dft_M1,dft_N1,fp, "original");

2. Perform DFT on blur kernel. Also perform inverse DFT to get back original contents. Make sure that the line for performing the inverse is commented out as it overwrites the DFT array.

		// Perform DFT of kernel
		dft_B = cvShowDFT(k_image,dft_M1,dft_N1,"kernel");
		//Perform inverse of kernel (check)
		cvShowInvDFT(k_image,dft_B,dft_M1,dft_N1,fp, "kernel");


    3. Multiply the DFT of image with the complex conjugate of the DFT of the blur kernel
		// Multiply numerator with complex conjugate
		dft_C = cvCreateMat( dft_M1, dft_N1, CV_64FC2 );
     4. Compute A**2 + B**2
		// Split Real and imaginary parts
		cvSplit( dft_B, image_ReB, image_ImB, 0, 0 );
    	        cvPow( image_ReB, image_ReB, 2.0);
		cvPow( image_ImB, image_ImB, 2.0);
		cvAdd(image_ReB, image_ImB, image_ReB,0);
     5. Divide numerator with A**2 + B**2
		//Divide Numerator/A^2 + B^2
		cvDiv(image_ReC, image_ReB, image_ReC, 1.0);
		cvDiv(image_ImC, image_ReB, image_ImC, 1.0);
                  6.Merge real and imaginary parts
		// Merge Real and complex parts
		cvMerge(image_ReC, image_ImC, NULL, NULL, complex_ImC);
         7.Finally perform Inverse DFT
		cvShowInvDFT(im, complex_ImC,dft_M1,dft_N1,fp,"deblur");

I have included the complete re-worked code below for the process of de-convolution. There was a small bug in the initial version. This has been fixed and the code below should work as is.

Note: You can clone this code from GitHub at “Deblurring by deconvolution in OpenCV
Do also take a look at my post “Reworking the Lucy Richardson algorithm in OpenCV
/*
============================================================================
Name : deconv.c
Author : Tinniam V Ganesh
Version :
Copyright :
Description : Deconvolution in OpenCV
============================================================================
*/

#include
#include
#include “cxcore.h”
#include “cv.h”
#include “highgui.h”

int main(int argc, char ** argv)
{
int height,width,step,channels;
uchar* data;
uchar* data1;
int i,j,k;
float s;

CvMat *dft_A;
CvMat *dft_B;
CvMat *dft_C;
IplImage* im;
IplImage* im1;
IplImage* image_ReB;
IplImage* image_ImB;
IplImage* image_ReC;
IplImage* image_ImC;
IplImage* complex_ImC;

IplImage* image_ReDen;
IplImage* image_ImDen;

FILE *fp;
fp = fopen(“test.txt”,”w+”);

int dft_M,dft_N;
int dft_M1,dft_N1;
CvMat* cvShowDFT();
void cvShowInvDFT();

im1 = cvLoadImage( “kutty-1.jpg”,1 );
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);
// Create blur kernel (non-blind)
//float vals[]={.000625,.000625,.000625,.003125,.003125,.003125,.000625,.000625,.000625};
//float vals[]={-0.167,0.333,0.167,-0.167,.333,.167,-0.167,.333,.167};

float vals[]={.055,.055,.055,.222,.222,.222,.055,.055,.055};
CvMat kernel = cvMat(3, // number of rows
3, // number of columns
CV_32FC1, // matrix data type
vals);
IplImage* k_image_hdr;
IplImage* k_image;

k_image_hdr = cvCreateImageHeader(cvSize(3,3),IPL_DEPTH_64F,2);
k_image = cvCreateImage(cvSize(3,3),IPL_DEPTH_64F,1);
k_image = cvGetImage(&kernel,k_image_hdr);

/*IplImage* k_image;
k_image = cvLoadImage( “kernel4.bmp”,0 );*/
cvNamedWindow(“blur kernel”, 0);

height = k_image->height;
width = k_image->width;
step = k_image->widthStep;

channels = k_image->nChannels;
//data1 = (float *)(k_image->imageData);
data1 = (uchar *)(k_image->imageData);

cvShowImage(“blur kernel”, k_image);

dft_M = cvGetOptimalDFTSize( im->height – 1 );
dft_N = cvGetOptimalDFTSize( im->width – 1 );

//dft_M1 = cvGetOptimalDFTSize( im->height+99 – 1 );
//dft_N1 = cvGetOptimalDFTSize( im->width+99 – 1 );

dft_M1 = cvGetOptimalDFTSize( im->height+3 – 1 );
dft_N1 = cvGetOptimalDFTSize( im->width+3 – 1 );

// Perform DFT of original image
dft_A = cvShowDFT(im, dft_M1, dft_N1,”original”);
//Perform inverse (check & comment out) – Commented as it overwrites dft_A
//cvShowInvDFT(im,dft_A,dft_M1,dft_N1,fp, “original”);

// Perform DFT of kernel
dft_B = cvShowDFT(k_image,dft_M1,dft_N1,”kernel”);
//Perform inverse of kernel (check & comment out) – commented as it overwrites dft_B
//cvShowInvDFT(k_image,dft_B,dft_M1,dft_N1,fp, “kernel”);

// Multiply numerator with complex conjugate
dft_C = cvCreateMat( dft_M1, dft_N1, CV_64FC2 );

printf(“%d %d %d %d\n”,dft_M,dft_N,dft_M1,dft_N1);

// Multiply DFT(blurred image) * complex conjugate of blur kernel
cvMulSpectrums(dft_A,dft_B,dft_C,CV_DXT_MUL_CONJ);

// Split Fourier in real and imaginary parts
image_ReC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
image_ImC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
complex_ImC = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 2);

printf(“%d %d %d %d\n”,dft_M,dft_N,dft_M1,dft_N1);

//cvSplit( dft_C, image_ReC, image_ImC, 0, 0 );
cvSplit( dft_C, image_ReC, image_ImC, 0, 0 );

// Compute A^2 + B^2 of denominator or blur kernel
image_ReB = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);
image_ImB = cvCreateImage( cvSize(dft_N1, dft_M1), IPL_DEPTH_64F, 1);

// Split Real and imaginary parts
cvSplit( dft_B, image_ReB, image_ImB, 0, 0 );
cvPow( image_ReB, image_ReB, 2.0);
cvPow( image_ImB, image_ImB, 2.0);
cvAdd(image_ReB, image_ImB, image_ReB,0);

//Divide Numerator/A^2 + B^2
cvDiv(image_ReC, image_ReB, image_ReC, 1.0);
cvDiv(image_ImC, image_ReB, image_ImC, 1.0);

// Merge Real and complex parts
cvMerge(image_ReC, image_ImC, NULL, NULL, complex_ImC);

// Perform Inverse
cvShowInvDFT(im, complex_ImC,dft_M1,dft_N1,fp,”deblur”);
cvWaitKey(-1);
return 0;
}

CvMat* cvShowDFT(im, dft_M, dft_N,src)
IplImage* im;
int dft_M, dft_N;
char* src;
{
IplImage* realInput;
IplImage* imaginaryInput;
IplImage* complexInput;

CvMat* dft_A, tmp;

IplImage* image_Re;
IplImage* image_Im;

char str[80];

double m, M;

realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);

cvScale(im, realInput, 1.0, 0.0);
cvZero(imaginaryInput);
cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);

dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );
image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);

// copy A to dft_A and pad dft_A with zeros
cvGetSubRect( dft_A, &tmp, cvRect(0,0, im->width, im->height));
cvCopy( complexInput, &tmp, NULL );
if( dft_A->cols > im->width )
{
cvGetSubRect( dft_A, &tmp, cvRect(im->width,0, dft_A->cols – im->width, im->height));
cvZero( &tmp );
}

// no need to pad bottom part of dft_A with zeros because of
// use nonzero_rows parameter in cvDFT() call below

cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );

strcpy(str,”DFT -“);
strcat(str,src);
cvNamedWindow(str, 0);

// Split Fourier in real and imaginary parts
cvSplit( dft_A, image_Re, image_Im, 0, 0 );

// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re, NULL);
cvPow( image_Re, image_Re, 0.5 );

// Compute log(1 + Mag)
cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)

cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
cvShowImage(str, image_Re);
return(dft_A);
}

void cvShowInvDFT(im, dft_A, dft_M, dft_N,fp, src)
IplImage* im;
CvMat* dft_A;
int dft_M,dft_N;
FILE *fp;
char* src;
{

IplImage* realInput;
IplImage* imaginaryInput;
IplImage* complexInput;

IplImage * image_Re;
IplImage * image_Im;

double m, M;
char str[80];

realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);

image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);

//cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, complexInput->height );
cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, dft_M);
strcpy(str,”DFT INVERSE – “);
strcat(str,src);
cvNamedWindow(str, 0);

// Split Fourier in real and imaginary parts
cvSplit( dft_A, image_Re, image_Im, 0, 0 );

// Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re, NULL);
cvPow( image_Re, image_Re, 0.5 );

// Compute log(1 + Mag)
cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
cvLog( image_Re, image_Re ); // log(1 + Mag)

cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
//cvCvtColor(image_Re, image_Re, CV_GRAY2RGBA);

cvShowImage(str, image_Re);
}


See also

My book ‘Practical Machine Learning with R and Python’ on Amazon
– De-blurring revisited with Wiener filter using OpenCV
–  Dabbling with Wiener filter using OpenCV
– Deblurring with OpenCV: Wiener filter reloaded
– Re-working the Lucy-Richardson Algorithm in 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

Find me on Google+

OpenCV: Fun with filters and convolution


My initial encounter with filters, convolution and correlation in OpenCV made me play around with the filters for Gaussian smooth, erosion and dilation operations on random image files. However I found the experience rather unsatisfactory and I wanted to get a real feel for the working of these operations. Suddenly a thought struck me. Could I restore an old family photograph of my parents? The photograph has areas of white patches that had to be removed.

So I started to dig a little more into the filters, convolution and correlation matrices to get a better understanding.

This is the original photograph.

My Mom & late Dad
As can be seen there are large patches in several places in the photograph. So I decided to use cvFloodFill to fill these areas.

a)      cvFloodFill: Since I had to identify the spots where these patches occurred I took a dump of the cvMat of the image which I resized to about 29 x 42. By inspecting the data I could see that the patches typically corresponded to intensity values that were greater than 170. So I decided that the cvFloodFill should happen with the seed around these parts.   So the code checks the intensity values > 170 and calls cvFloodFill. After much tweaking I could see that the white patches were now filled with gray (newval intensity= 150.0) So I was able to get rid of the white patches.

 

b)      cvSmooth: The next step that I took was to perform a Gaussian smooth of the picture. This smoothed out the filled parts

 

c)       cvErode & cvDilate: I followed this with cvErode to smooth out the dark areas  and cvDilate to smooth out the bright areas.

 

d)      cvFilter2D:  I wanted to now sharpen the image. I did a lot of experiments with different kernel values but I found this to be extremely difficult to work with. After much trial and error I came with a kernel values of
double a[9]={-1,20,1,-1,20,1,-1,20,1};
The sharpening was reasonable but there are areas where there are white streaks. I still need to figure out a kernel that can sharpen images. For this also to understand what was happening I tried to dump the values of the image to get a feel of where the values lay.

 

e)      cvSmooth: Finally I performed a cvSmooth of the filtered output.

While I have had fair success there is still a lot more left to be desired from the final version.
The complete process flow

The code is included below

#include “cv.h”
#include “highgui.h”
#include “stdio.h”
int main(int argc, char** argv)
{
IplImage* img = cvLoadImage(“dad_mom.jpg”,0);
IplImage* dst;
IplImage* dst1;
IplImage* dst2;
IplImage* dst3;
IplImage* dst4;
int i,j,k;
int height,width,step,channels;
uchar* data;
uchar* data1;
uchar* outdata;
CvScalar s;
CvScalar lodiff,highdiff,newval;
// get the image data
height    = img->height;
width     = img->width;
step      = img->widthStep;
channels  = img->nChannels;
data      = (uchar *)img->imageData;
double a[9]={-1,20,1,-1,20,1,-1,20,1};
//double a[9]={1/16,1/8,1/16,1/8,1/4,1/8,1/16,1/8,1/16};
double values[9]={1/16,0,-1/16,2/16,0,-2/16,1/16,0,-1/16};
CvPoint seed;
CvMat kernel= cvMat(3,3,CV_32FC1,a);
printf(“Processing a %d x %d image with %d channels\n”,height,width,channels);
// Create windows
cvNamedWindow(“Original”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Flood Fill”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Smooth”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Erode”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Dilate”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Filter”,CV_WINDOW_AUTOSIZE);
cvNamedWindow(“Smooth1”,CV_WINDOW_AUTOSIZE);
// Original image
cvShowImage(“Original”, img);
/* Flood fill in white patches intensity > 170.0 */
highdiff=cvRealScalar(5.0);
lodiff=cvRealScalar(5.0);
newval=cvRealScalar(150.0);
for(i=0;i<height;i++) {
for(j=0;j<width;j++) {
for(k=0;k<channels;k++) {
//printf(“Data = %d \n”,data[i*step+j*channels+k]);
if((data[i*step+j*channels+k]) > 170.0)
{
seed=cvPoint(j,i);
//printf(“data=%dFlood
seed=%d,%d\n”,data[i*step+j*channels+k],i,j);
cvFloodFill(img,seed,newval,lodiff,highdiff, NULL,CV_FLOODFILL_FIXED_RANGE,NULL);                                                                              }
else
{
;
}
}
}
//printf(“\n”);
}
cvShowImage(“Flood Fill”,img);
// Gaussian smooth
dst = cvCloneImage(img);
cvSmooth( img, dst, CV_GAUSSIAN, 3, 3, 0, 0 );
cvShowImage(“Smooth”,dst);
// Erode the image
dst1 = cvCloneImage(img);
IplConvKernel* kern = cvCreateStructuringElementEx(3,3,1,1,CV_SHAPE_RECT,values);
cvErode(dst,dst1,kern,1);
cvShowImage(“Erode”,dst1);
// Perform dilation operation
dst2 = cvCloneImage(img);
cvDilate(dst1,dst2,kern,1);
cvShowImage(“Dilate”,dst2);
// Filter the image with convolution kernel. Sharpen the image
dst3 = cvCloneImage(img);
printf(“reached here\n”);
cvFilter2D(dst2,dst3,&kernel,cvPoint(-1,-1));
cvShowImage(“Filter”,dst3);
// Smoothen the image
dst4 = cvCloneImage(img);
cvSmooth( dst3, dst4, CV_MEDIAN, 3, 0, 0, 0 );
cvShowImage(“Smooth1”,dst4);
// Cleanup
cvWaitKey(0);
cvReleaseImage(&img);
cvReleaseImage(&dst);
cvDestroyWindow(“Original”);
vDestroyWindow(“Restore”);

}


Find me on Google+