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 1st 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

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*);

cvNamedWindow(“Original-Color”, 0);
cvShowImage(“Original-Color”, im1);
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);;

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;
}

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.

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*);

cvNamedWindow(“Original-Color”, 0);
cvShowImage(“Original-Color”, im1);
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 = 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);
val = cvScalarAll(kappa);
//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);
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);
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);
}

De-blurring revisited with Wiener filter using OpenCV

In this post I continue to experiment with the de-blurring of images using the Wiener filter. For details on the Wiener filter, please look at my earlier post “Dabbling with Wiener filter using OpenCV”.  Thanks to Egli Simon, Switzerland for pointing out a bug in the earlier post which I have now fixed.  The Wiener filter attempts to de-blur by assuming that the source signal is convolved with a blur kernel in the presence of noise.   I have also included the blur kernel as estimated by E. Simon in the code. I am including the de-blurring with 3 different blur kernel radii and different values for the Wiener constant K.  While the de-blurring is still a long way off there is some success.

One of the reasons I have assumed a non-blind blur kernel and try to de-convolve with that. The Wiener filter tries to minimize the Mean Square Error (MSE)  which can be expressed as
e(f) = E[X(f) – X1(f)]^2                                – (1)

where e(f) is the Mean Square Error(MSE) in the frequency domain, X(f) is the original image and X1(f) is the estimated signal which we get by de-convolving the Wiener filter with the observed blurred image i.e. and E[] is the expectation

X1(f) = G(f) * Y(f)                                        -(2)
where G(f) is the Wiener de-convolution filter and Y(f) is the observed blurred image
Substituting (2) in (1) we get
e(f) = E[X(f) – G(f) *Y(f)]^2

If the above equation is solved we can effectively remove the blur.
Feel free to post comments, opinions or ideas.

Watch this space …
I will be back! Hasta la Vista!

Note: You can clone the code below from Git Hib – Another implementation of Weiner filter in OpenCV

The complete code is included below and should work as is

/*
============================================================================
Name : deblur_wiener.c
Author : Tinniam V Ganesh & Egli Simon
Version :
Description : Implementation of Wiener filter in OpenCV
============================================================================
*/

#include <stdio.h>
#include <stdlib.h>
#include “cxcore.h”
#include “cv.h”
#include “highgui.h”

#define kappa 10.0

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;
char str[80];
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*);

cvNamedWindow(“Original-color”, 0);
cvShowImage(“Original-color”, im1);
if( !im )
return -1;

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

// Create a random noise matrix
fp = fopen(“test.txt”,”w+”);
int val_noise[357*383];
for(i=0; i <im->height;i++){
for(j=0;j<im->width;j++){
fprintf(fp, “%d “,(383*i+j));
val_noise[383*i+j] = rand() % 128;
}
fprintf(fp, “\n”);
}

CvMat noise = cvMat(im->height,im->width, CV_8UC1,val_noise);

// Add the random noise matric to the image

cvNamedWindow(“Original + Noise”, 0);
cvShowImage(“Original + Noise”, im);

cvSmooth( im, im, CV_GAUSSIAN, 7, 7, 0.5, 0.5 );
cvNamedWindow(“Gaussian Smooth”, 0);
cvShowImage(“Gaussian Smooth”, im);

// Create a blur kernel
IplImage* k_image;

printf(“rowlength %d\n”,rowLength);
float kernels[rowLength*rowLength];
printf(“rowl: %i”,rowLength);
int norm=0; //Normalization factor
int x,y;
CvMat kernel;
for(x = 0; x < rowLength; x++)
for (y = 0; y < rowLength; y++)
norm++;
// Populate matrix
for (y = 0; y < rowLength; y++) //populate array with values
{
for (x = 0; x < rowLength; x++) {
//kernels[y * rowLength + x] = 255;
kernels[y * rowLength + x] =1.0/norm;
printf(“%f “,1.0/norm);
}
else{
kernels[y * rowLength + x] =0;
}
}
}

/*for (i=0; i < rowLength; i++){
for(j=0;j < rowLength;j++){
printf(“%f “, kernels[i*rowLength +j]);
}
}*/

kernel= cvMat(rowLength, // number of rows
rowLength, // number of columns
CV_32FC1, // matrix data type
&kernels);
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);
val = cvScalarAll(kappa);

//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,str);

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);
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);
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, CV_GRAY2RGBA);
cvShowImage(str, image_Re);
}
– 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

Dabbling with Wiener filter using OpenCV

The technique of reduction of blur and restoration of images is an extremely important field of study and finds numerous applications in medical imaging and astronomy.  One such technique is Wiener filter named after the originator of the technique Prof. Norbert Wiener (1894-1964).  It is usually important to consider these techniques in the frequency domain.

This can be represented diagrammatically as follows

Given that we have
b(x,y) = i(x,y) * k(x,y) + n(x,y)                                       – (1)
where b(x,y) is the blurred image, i(x,y) is the latent image, k(x,y) is the blur kernel and n(x,y) is random noise.
The above in frequency domain can be written as
B(u,v) = I(u,v) * K(u,v) + N(u,v)                                     –  (2)
Inv K(u,v) represents the inverse of the blur kernel and can be determined  from equation (2) by disregarding the noise.  However the inverse filter has high gain at low frequencies and hence tends to amplify noise at these low frequencies.
Wiener filter and Butterworth filter try to minimize these
Wiener filter attempts to minimize the Mean Squared Error (MSE) of the following equation
e^2 = E[( i(x,y)   – mean (i (x,y))^2]

The frequency domain expression for this
T(u,v)           =
K*(u,v)
——–                                                   -(3)
| K(u,v)|^2 + Sn(u,v)/Sf(u,v)
Where T(u,v) is the Wiener filter
K*(u,v) is the complex conjugate of the blur kernel
|K(u,v)|^2  = k^2 + m^2 where k and m are the coefficients of the real and imaginary parts of the blur kernel
Sn(u,v) is the power spectra of noise
Sf(u,v) is the power spectra of the signal
Equation (3) can be represented as
T(u,v)            =
K*(u,v)
——–                                                   –   (4)
| K(u,v)|^2 +  L
Where L is a small positive value
Eqn(4) represents Wiener filter

The lecture  “Image deblurring by Frequency Domain Operations” by Prof. Harvey Rhody is worth a read.
One thing that puzzles me is the Wiener filter is the inverse of the blur kernel with the additional factor of ‘L’ in the denominator. The contribution of the factor ‘L’ will be small, so the Wiener filter appears to be the same as a regular inverse filter. I am probably missing something here(??).

You can clone the code below from Git Hub – Weiner filter with OpenCV

The steps to create the Wiener filter were as follow

1. Create the gray image of the blurred original

// Add the random noise matric to the image
3.Gaussian smooth the noisy blurred image
cvSmooth( im, im, CV_GAUSSIAN, 7, 7, 0.5, 0.5 );
4.Perform DFT and inverse DFT (to check) of the original image (noisy, blurred image)
// 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");
5.Perform DFT and inverse DFT of blur kernel
// 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");

6.Multiply DFT of original with complex conjugate of blur kernel
// Multiply numerator with complex conjugate
dft_C = cvCreateMat( dft_M1, dft_N1, CV_64FC2 );
7.Calculate modulus of blur kernel (denominator)
// 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);
8.Add Wiener constant to modulus of blur kernel (denominator)
val = cvScalarAll(0.00001);
9.Divide the numerator with the denominator
//Divide Numerator/A^2 + B^2
cvDiv(image_ReC, image_ReB, image_ReC, 1.0);
cvDiv(image_ImC, image_ReB, image_ImC, 1.0);
10.Perform inverse DFT of the filtered signal

// Perform Inverse

cvShowInvDFT1(im, (CvMat *)complex_ImC,dft_M1,dft_N1,"Output - Wiener filter");

There is still a residual amount of blur (do also read De-blurring revisited with Wiener filter using OpenCV). This can be tried for different values of the Wiener introduced noise.
Note: I would like to give special thanks to Egli Simon who identified the bug in my previous version of the code. The line to display the inverse DFT of the original image & the blur kernel overwrote the DFT contents and I was not getting the Wiener filter to work
Note: You can clone the code below from Git Hub – Weiner filter with OpenCV
The complete code is given below (re-worked and should work as is)
/*
============================================================================
Name : wiener1.c
Author : Tinniam V Ganesh & Egli Simon
Version :
Description : Wiener filter implementation in OpenCV (22 Nov 2011)
============================================================================
*/
#include <stdio.h>
#include <stdlib.h>
#include “cxcore.h”
#include “cv.h”
#include “highgui.h”

#define kappa 1
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*);
cvNamedWindow(“Original-color”, 0);
cvShowImage(“Original-color”, im1);
if( !im )
return -1;
cvNamedWindow(“Original-gray”, 0);
cvShowImage(“Original-gray”, im);

// Create a random noise matrix
fp = fopen(“test.txt”,”w+”);
int val_noise[357*383];
for(i=0; i <im->height;i++){
for(j=0;j<im->width;j++){
fprintf(fp, “%d “,(383*i+j));
val_noise[383*i+j] = rand() % 128;
}
fprintf(fp, “\n”);
}

CvMat noise = cvMat(im->height,im->width, CV_8UC1,val_noise);

// Add the random noise matric to the image
cvNamedWindow(“Original + Noise”, 0);
cvShowImage(“Original + Noise”, im);
cvSmooth( im, im, CV_GAUSSIAN, 7, 7, 0.5, 0.5 );
cvNamedWindow(“Gaussian Smooth”, 0);
cvShowImage(“Gaussian Smooth”, im);

// Create a blur kernel
IplImage* k_image;
printf(“rowlength %d\n”,rowLength);
float kernels[rowLength*rowLength];
printf(“rowl: %i”,rowLength);
int norm=0; //Normalization factor
int x,y;
CvMat kernel;
for(x = 0; x < rowLength; x++)
for (y = 0; y < rowLength; y++)
norm++;
// Populate matrix
for (y = 0; y < rowLength; y++) //populate array with values
{
for (x = 0; x < rowLength; x++) {
//kernels[y * rowLength + x] = 255;
kernels[y * rowLength + x] =1.0/norm;
printf(“%f “,1.0/norm);
}
else{
kernels[y * rowLength + x] =0;
}
}
}

/*for (i=0; i < rowLength; i++){
for(j=0;j < rowLength;j++){
printf(“%f “, kernels[i*rowLength +j]);
}
}*/

kernel= cvMat(rowLength, // number of rows
rowLength, // number of columns
CV_32FC1, // matrix data type
&kernels);
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);

val = cvScalarAll(kappa);

//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,”O/p Wiener k=1 rad=2″);

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);
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);
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);

}

Checkout my book ‘Deep Learning from first principles Second Edition- In vectorized Python, R and Octave’.  My book is available on Amazon  as paperback (\$18.99) and in kindle version(\$9.99/Rs449).

You may also like my companion book “Practical Machine Learning with R and Python:Second Edition- Machine Learning in stereo” available in Amazon in paperback(\$12.99) and Kindle(\$9.99/Rs449) versions.

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);
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 :
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();

cvNamedWindow(“original-color”, 0);
cvShowImage(“original-color”, im1);
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 = cvCreateImage(cvSize(3,3),IPL_DEPTH_64F,1);
k_image = cvGetImage(&kernel,k_image_hdr);

/*IplImage* k_image;
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);

//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);
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);
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);
}

Checkout my book ‘Deep Learning from first principles Second Edition- In vectorized Python, R and Octave’.  My book is available on Amazon  as paperback (\$18.99) and in kindle version(\$9.99/Rs449).

You may also like my companion book “Practical Machine Learning with R and Python:Second Edition- Machine Learning in stereo” available in Amazon in paperback(\$12.99) and Kindle(\$9.99/Rs449) versions.

Installing and using OpenCV with Visual Studio 2010 express

Here’s the low-down on installing and getting started with OpenCV on Visual Studio Express 2010 on Windows.

You can untar the file into a directory named /opencv2.2/

3) Once Visual Studio 2010  is installed open VC++ and select  “New Project” and choose Win32 Console Application for e.g. first

4) Include the any relevant OpenCV code into first.cpp (snippet given below)

5) Include the following 3 directories under
Browse to the folder containing the files below and select them
Project->first properties->VC++ Directories->Include Directories. Click Edit
.\opencv2.2\include
.\opencv2.2\include\opencv
.\opencv2.2\include\opencv2
Click Apply and OK

6) Now include the library path
Browse to the folder below

Project->first properties->Library Directories
.\opencv2.2\lib
Click Apply and OK

7) Now the last step is to include all the necessary OpenCV libraries during the linking phase
For this go to Project->first properties->Linker->Input->Additional Dependencies and cut and paste all the following libraries by clicking the “Edit”

opencv_calib3d220.lib
opencv_calib3d220d.lib
opencv_contrib220.lib
opencv_contrib220d.lib
opencv_core220.lib
opencv_core220d.lib
opencv_features2d220.lib
opencv_features2d220d.lib
opencv_ffmpeg220.lib
opencv_ffmpeg220d.lib
opencv_flann220.lib
opencv_flann220d.lib
opencv_gpu220.lib
opencv_gpu220d.lib
opencv_highgui220.lib
opencv_highgui220d.lib
opencv_imgproc220.lib
opencv_imgproc220d.lib
opencv_legacy220.lib
opencv_legacy220d.lib
opencv_ml220.lib
opencv_ml220d.lib
opencv_objdetect220.lib
opencv_objdetect220d.lib
opencv_ts220.lib
opencv_video220.lib
opencv_video220d.lib

Click Apply and Ok.

8) Copy the baboon.jpg from the .\opencv2.2\samples\cpp to your projects\first\first folder

9) Now you are good to go. Build by choosing Debug->Build Solution
It should through fine  when  you now run the code
You should see

Here is a sample snippet

#include "stdafx.h"
#include "cv.h"
#include "highgui.h"
int main( int argc, char** argv ) {
// cvLoadImage determines an image type and creates datastructure with appropriate size
// create a window. Window name is determined by a supplied argument
cvNamedWindow( "test", CV_WINDOW_AUTOSIZE );
// Apply Gaussian smooth
//cvSmooth( img, img, CV_GAUSSIAN, 9, 9, 0, 0 );
cvErode (img, img,NULL,2);
// Display an image inside and window. Window name is determined by a supplied argument
cvShowImage( "test", img );
//Save image
cvSaveImage( "c:\\shanthi\\baboon2.png", img, 0);
// wait indefinitely for keystroke
cvWaitKey(0);
// release pointer to an object
cvReleaseImage( &img );
// Destroy a window
cvDestroyWindow( argv[1] );
}

Checkout my book ‘Deep Learning from first principles Second Edition- In vectorized Python, R and Octave’.  My book is available on Amazon  as paperback (\$18.99) and in kindle version(\$9.99/Rs449).

You may also like my companion book “Practical Machine Learning with R and Python:Third Edition- Machine Learning in stereo” available in Amazon in paperback(\$12.99) and Kindle(\$9.99/Rs449) versions.

Computer Vision: Ramblings on derivatives, histograms and contours

Images can be visualized to be functions of the form f(x,y) where f(x,y) represents the intensity at the pixel position ,y. However images can be grayscale, color or four channels and each channel may consist of integers or floating point numbers. However the changes in the values can be viewed as a continuous function. Here is a nice representation of an image as a continuous function (courtesy: Prof Darrell’s lecture at Berkeley on Filters)

Given that the image can be viewed as a continuous function in 2 or 3 axes we have derivatives that can be taken of the images. The derivates determine the maximum and minimum of this changing function. The key derivatives in image processing are the Sobel, the Scharr and the Laplacian filters. These provide the 1st order or 2nd order derivative and hence can be used for determining edges of an image.

I was keen on playing around with derivatives and also understanding how the histograms look like.

Here is the original image and its histogram. Clearly there is a nice spread of the values.

Sobel filter and its histogram

The output of the Sobel filter on the original image is shown. The edges with Sobel’s derivative somehow are not too pronounced. The Sobel derivate can be used for obtaining the gradient of the image. The corresponding histogram of the Sobel’s gradient is shown.

The code snippet ( The complete code is given below)

IplImage* out_sobel = cvCreateImage( cvSize(img->width, img->height), IPL_DEPTH_16S, 1);
cvSobel(in_gray, out_sobel, 1,1,7);
cvShowImage("Sobel", out_sobel);

//create an image to hold the histogram
IplImage* histImage_Sobel = cvCreateImage(cvSize(300,400), 8, 1);

//create a histogram to store the information from the image
CvHistogram* histSobel = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);

//calculate the histogram and apply to hist
cvCalcHist( &histImage_Sobel, histSobel, 0, NULL );

//grab the min and max values and their indeces
cvGetMinMaxHistValue( histSobel, &min_value, &max_value, &min_idx, &max_idx);

//scale the bin values so that they will fit in the image representation
cvScale( histSobel->bins, histSobel->bins, ((double)histImage_Sobel->height)/max_value, 0 );

//set all histogram values to 255
cvSet( histImage_Sobel, cvScalarAll(255), 0 );

//create a factor for scaling along the width
bin_w = cvRound((double)histImage_Sobel->width/hist_size);

for( i = 0; i < hist_size; i++ ) {
//draw the histogram data onto the histogram image
cvRectangle( histImage_Sobel, cvPoint(i*bin_w, histImage_Sobel->height),
cvPoint((i+1)*bin_w,
histImage_Sobel->height - cvRound(cvGetReal1D(histSobel->bins,i))),
cvScalarAll(0), -1, 8, 0 );
//get the value at the current histogram bucket
float* bins = cvGetHistValue_1D(histSobel,i);
//increment the mean value
mean += bins[0];
}

cvShowImage("Hist Sobel",histImage_Sobel);
...
...
(Please see Gavin S. Page's tutorial(vast.uccs.edu/~tboult/CS330/NOTES/OpenCVTutorial_II.ppt) on histograms)

Laplacian and its histogram

The Laplacian provides the 2nd order derivative and hence can be used to determine local maxima and local minima. The Laplacian provides for much more pronounced edges and can be used to extract features of an object of interest. Its corresponding histogram is also included.

Canny filter and Contours

The third filter is cvCanny which is most suitable for obtaining clear edges in an image. The canny is usually used along with cvFindContours to determine the general shape of an object. I used the canny filter which I passed to a contour detecting function. However the contour detecting function identified more than 228 contours most of which were useless except for 1 which had included the complete contour of the hand as shown.

However when I increased the max_depth to 1 I found that it was immediately able to get the complete contour of the hand besides a lot of extraneous contours.

I guess the challenge with the contour function is being able to programmatically reject all those contours which of lesser importance (possibly a future post).

Code for Sobel, Laplacian and Histograms

#include "cv.h"
#include "highgui.h"
#include "stdio.h"

int main(int argc, char** argv)
{
IplImage* dst;
IplImage* in_gray;
int hist_size=30;
float gray_ranges[] = { 0, 255 };
float* ranges[]     = { gray_ranges};
int min_idx,max_idx;
float min_value,max_value;
int bin_w;
int i;
float mean,variance;

cvNamedWindow("Original",CV_WINDOW_AUTOSIZE);
cvNamedWindow("histogram",CV_WINDOW_AUTOSIZE);
cvNamedWindow("Sobel",CV_WINDOW_AUTOSIZE);
cvNamedWindow("Hist Sobel",CV_WINDOW_AUTOSIZE);

cvNamedWindow("Laplacian",CV_WINDOW_AUTOSIZE);
cvNamedWindow("Hist Laplace",CV_WINDOW_AUTOSIZE);

in_gray = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_8U, 1);
cvCvtColor(img, in_gray, CV_BGR2GRAY);
cvShowImage("Original", in_gray);

//create a rectangular area to evaluate
CvRect rect = cvRect(0, 0, 300, 400 );
//apply the rectangle to the image and establish a region of interest
cvSetImageROI(in_gray, rect);

//create an image to hold the histogram
IplImage* histImage = cvCreateImage(cvSize(300,400), 8, 1);

//create a histogram to store the information from the image
CvHistogram* hist = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);

//calculate the histogram and apply to hist
cvCalcHist( &in_gray, hist, 0, NULL );

//grab the min and max values and their indeces
cvGetMinMaxHistValue( hist, &min_value, &max_value, &min_idx, &max_idx);

//scale the bin values so that they will fit in the image representation
cvScale( hist->bins, hist->bins, ((double)histImage->height)/max_value, 0 );

//set all histogram values to 255
cvSet( histImage, cvScalarAll(255), 0 );

//create a factor for scaling along the width
bin_w = cvRound((double)histImage->width/hist_size);

for( i = 0; i < hist_size; i++ ) {
//draw the histogram data onto the histogram image
cvRectangle( histImage, cvPoint(i*bin_w, histImage->height),
cvPoint((i+1)*bin_w,
histImage->height - cvRound(cvGetReal1D(hist->bins,i))),
cvScalarAll(0), -1, 8, 0 );
//get the value at the current histogram bucket
float* bins = cvGetHistValue_1D(hist,i);
//increment the mean value
mean += bins[0];
}

//finish mean calculation
mean /= hist_size;

//go back through now that mean has been calculated in order to calculate variance
for( i = 0; i < hist_size; i++ ) {
float* bins = cvGetHistValue_1D(hist,i);
variance += pow((bins[0] - mean),2);
}
//finish variance calculation
variance /= hist_size;

cvShowImage("histogram",histImage);

IplImage* out_sobel = cvCreateImage( cvSize(img->width, img->height), IPL_DEPTH_16S, 1);
cvSobel(in_gray, out_sobel, 1,1,7);
cvShowImage("Sobel", out_sobel);

//create an image to hold the histogram
IplImage* histImage_Sobel = cvCreateImage(cvSize(300,400), 8, 1);

//create a histogram to store the information from the image
CvHistogram* histSobel = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);

//calculate the histogram and apply to hist
cvCalcHist( &histImage_Sobel, histSobel, 0, NULL );

//grab the min and max values and their indeces
cvGetMinMaxHistValue( histSobel, &min_value, &max_value, &min_idx, &max_idx);

//scale the bin values so that they will fit in the image representation
cvScale( histSobel->bins, histSobel->bins, ((double)histImage_Sobel->height)/max_value, 0 );

//set all histogram values to 255
cvSet( histImage_Sobel, cvScalarAll(255), 0 );

//create a factor for scaling along the width
bin_w = cvRound((double)histImage_Sobel->width/hist_size);

for( i = 0; i < hist_size; i++ ) {
//draw the histogram data onto the histogram image
cvRectangle( histImage_Sobel, cvPoint(i*bin_w, histImage_Sobel->height),
cvPoint((i+1)*bin_w,
histImage_Sobel->height - cvRound(cvGetReal1D(histSobel->bins,i))),
cvScalarAll(0), -1, 8, 0 );
//get the value at the current histogram bucket
float* bins = cvGetHistValue_1D(histSobel,i);
//increment the mean value
mean += bins[0];
}

cvShowImage("Hist Sobel",histImage_Sobel);

// Create Laplacian and the histogram for it
IplImage *output=cvCreateImage( cvSize(img->width, img->height), IPL_DEPTH_16S, 1);
cvLaplace(in_gray, output, 7);
cvShowImage("Laplacian", output);

//create an image to hold the histogram
IplImage* histImage_Laplace = cvCreateImage(cvSize(300,400), 8, 1);

//create a histogram to store the information from the image
CvHistogram* histLaplace = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);

//calculate the histogram and apply to hist
cvCalcHist( &histImage_Laplace, histLaplace, 0, NULL );

//grab the min and max values and their indeces
cvGetMinMaxHistValue( histLaplace, &min_value, &max_value, &min_idx, &max_idx);

//scale the bin values so that they will fit in the image representation
cvScale( histLaplace->bins, histLaplace->bins, ((double)histImage_Laplace->height)/max_value, 0 );

//set all histogram values to 255
cvSet( histImage_Laplace, cvScalarAll(255), 0 );

//create a factor for scaling along the width
bin_w = cvRound((double)histImage_Laplace->width/hist_size);

for( i = 0; i < hist_size; i++ ) {
//draw the histogram data onto the histogram image
cvRectangle( histImage_Laplace, cvPoint(i*bin_w, histImage_Laplace->height),
cvPoint((i+1)*bin_w,
histImage_Laplace->height - cvRound(cvGetReal1D(histLaplace->bins,i))),
cvScalarAll(0), -1, 8, 0 );
//get the value at the current histogram bucket
float* bins = cvGetHistValue_1D(histLaplace,i);
//increment the mean value
mean += bins[0];
}

cvShowImage("Hist Laplace",histImage_Laplace);

cvWaitKey(0);

printf("Mean= %f\n",mean);
printf("variance=%f\n",variance);

//clean up images
cvReleaseImage(&histImage_Laplace);
cvReleaseImage(&histImage_Sobel);
cvReleaseImage(&histImage);
cvReleaseImage(&in_gray);
cvReleaseImage(&img);

//remove windows
cvDestroyWindow("Original");

cvDestroyWindow("histogram");
}

Code for Canny & Contours
#include "cv.h"
#include "highgui.h"

#define CVX_RED		CV_RGB(0xff,0x00,0x00)
#define CVX_GREEN	CV_RGB(0x00,0xff,0x00)
#define CVX_BLUE	CV_RGB(0x00,0x00,0xff)

int main(int argc, char* argv[])
{
CvSeq* c;
int i;
cvNamedWindow("Original", 1 );
cvNamedWindow("Canny_Edge", 1 );
cvNamedWindow("Contours", 1 );
IplImage* img_edge = cvCreateImage( cvGetSize(img_8uc1), 8, 1 );
IplImage* img_8uc3 = cvCreateImage( cvGetSize(img_8uc1), 8, 3 );
cvThreshold( img_8uc1, img_edge, 128, 255, CV_THRESH_BINARY );
CvMemStorage* storage =cvCreateMemStorage(0);

CvSeq* first_contour = NULL;

int Nc;
int n=0;

cvShowImage("Original", img_8uc1);

IplImage *out_canny=cvCreateImage( cvSize(img_8uc1->width, img_8uc1->height), IPL_DEPTH_8U, 1);
cvCanny(img_8uc1, out_canny, 50.0 ,100.0, 3);
cvShowImage("Canny_Edge", out_canny);

/*	Nc = cvFindContours(
img_edge,
storage,
&first_contour,
sizeof(CvContour),
CV_RETR_LIST,
CV_CHAIN_APPROX_SIMPLE,
cvPoint(0,0)// Try all four values and see what happens
);*/

Nc = cvFindContours(
out_canny,
storage,
&first_contour,
sizeof(CvContour),
CV_RETR_TREE,
CV_CHAIN_APPROX_SIMPLE,
cvPoint(0,0)// Try all four values and see what happens
);

printf("Total contours detected: %d\n",Nc);

for(c=first_contour; c!=NULL; c=c->h_next )
{
cvCvtColor( img_8uc1, img_8uc3, CV_GRAY2BGR );
cvDrawContours(
img_8uc3,
c,
CVX_RED,
CVX_BLUE,
1,        // Try different values of max_level, and see what happens
2,
8,
cvPoint(0,0));
printf("Contour #%d\n",n);

cvShowImage("Contours", img_8uc3 );
printf(" %d elements: \n",c->total);

for(i=0; i<c->total; ++i ) {
CvPoint* p = CV_GET_SEQ_ELEM( CvPoint, c, i );
printf("(%d,%d)\n",p->x,p->y);

}
cvWaitKey(0);
n++;
}
printf("Finished all contours\n");
cvCvtColor( img_8uc1, img_8uc3, CV_GRAY2BGR );
cvShowImage( argv[0], img_8uc3 );
cvWaitKey(0);
cvDestroyWindow( argv[0] );
cvReleaseImage( &img_8uc1 );
cvReleaseImage( &img_8uc3 );
cvReleaseImage( &img_edge );
return 0;
}

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.

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* 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”);

}

Checkout my book ‘Deep Learning from first principles Second Edition- In vectorized Python, R and Octave’.  My book is available on Amazon  as paperback (\$18.99) and in kindle version(\$9.99/Rs449).

You may also like my companion book “Practical Machine Learning with R and Python:Second Edition- Machine Learning in stereo” available in Amazon in paperback(\$12.99) and Kindle(\$9.99/Rs449) versions.

Hand detection through Haartraining: A hands-on approach

Detection of objects from images or from video is no trivial task. We need to use some type of machine learning algorithm and train it to detect features and also identify misses and false positives.  The haartraining algorithm does just this. It creates a series of haarclassifiers which ensure that non-features are quickly rejected as the object is identified.

This post will highlight the necessary steps required to build a haarclassifier for detection a hand or any object of interest. This post is sequel to my earlier post (OpenCV: Haartraining and all that jazz!) and has a lot more detail. In order to train the haarclassifier, it is suggested, that at least 1000 positive samples (images with the object of interest- hand in this case) and 2000 negative samples (any other image) is required.

As before for performing haartraining the following 3 steps have to be performed
1)      Create samples (createsamples.cpp)
2)      Haar Training (haartraining.cpp)
3)      Performance testing of the classifier (performance.cpp)

In order to build the above 3 utilities please refer to my earlier post OpenCV: Haartraining and all that jazz!

The steps required for training a haarcascade to recognize on normal open palm is given below

Create Samples: This step can be further broken down into the following 3 steps
a)      Creation of positive samples
b)      Superimposing the positive sample on the negative sample
c)      Merging of vector files of samples.
A)      Creation of positive samples:

a)Get a series of images with objects of interest (positive samples):For this step take photos using the webcam (or camera) of the objects that you are interested. I had taken several snapshots of my hand ( I later simply downloaded images from Google images for because of the excessive clutter in the snaps taken).

b) Crop Images:In this step you need to crop the images such that it only contains the object of interest. You can use any photo editing tool of your choice and save all the positive images in a directory for e.g. ./hands
c) Mark the object:Now the image with the positive sample has to be marked. This will be used for creating the samples.
The tool that is to be used for marking is objectmarker.cpp. You can downloaded the source code for this from Achu Wilson’s blogNow build objectmarker.cpp with the include directories and libraries of OpenCV. Once you have successfully built objectmarker you are ready to mark the positive samples. Samples have to be marked  because the description file for the positive samples file must be in the following format

[filename] [# of objects] [[x y width height] [… 2nd object] …]

This will be used for generating the positive training samples.  This file is will be used with the createsamples utility to create positive training samples.
The command to use to run objectmarker is as follows
\$ objectmarker <output file> <dir>
for e.g.
\$objectmarker pos_desc.txt ./images
where dir is the directory containing the positive images
and output file will contain the positions of the objects marked as follows
pos_desc.txt
/images/hand1.bmp 1 0 0 246 50
/images/hand2.bmp 1 187 26 333 400
….
The use of the utility objectmarker is an art 😉 and you need to be trained to get the proper data. Make sure you get sensible widths and heights. There are times when you get -ve widths and -ve heights which are clearly wrong.

An alternative easy way is to open the jpg/bmp in a photo editor and check the width and height in pixels (188 x 200 say) and create the description file as
/images/hand1.bmp 1 0 0 188 200
Ideally the objectmarker should give you values close to this if the sample image file has only the object of interest (hand)
Create positive training samples
The createsamples utility can now be used for creating positive training sample files. The command is
createsamples -info pos_desc.txt -vec pos.vec -w 20 -h 20

This will create positive training samples with the object of interest, in our case the “hand”.  Now, you should verify that the tool has done something sensible by checking the training samples generated.  The command to use to display the positive training samples captured in the pos.vec is to run
createsamples -vec pos.vec -w 20 -h 20
If the objects have been marked accurately you should see a series of hands (positive samples). If the samples have not been marked correctly the positive training file will give incorrect results.
Create negative background training samples
This step is used to create training samples with one positive image superimposed against a set of negative background samples. The command to use is
createsamples -img ./image/hand_1.BMP -num 9 -bg bg.txt -vec neg1.vec -maxxangle 0.6 -maxyangle 0 -maxzangle 0.3 -maxidev 100 -bgcolor 0 -bgthresh 0 -w 20 -h 20
where bg.txt will contain the list of negative samples in the following format
./negative/airplane.jpg
./negative/baboon.jpg
./negative/cat.jpg

The positive hand images will be superimposed against the negative background in various angles of distortion.

All the negative training samples will be collected in the training file neg1.vec as above. As before you can verify that the createsamples utility has done something reasonable by executing

createsamples -vec neg1.vec -w 20 -h 20.
This should show a series of images of hands in various angles of distortion against the negative image background.
Creating several negative training samples with all positive samples
The createsamples utility takes one positive sample and superimposes it against the negative samples in bg.txt file. However we need to repeat this process for each of the positive sample (hand) that we have. Since this can be laborious process you can use the following shell script, I wrote,  to repeat the negative training sample with every positive sample with create_neg_training.sh

create_neg_training.sh
#!/bin/bash
let j=1
for  i in `cat hands.txt`
do
createsamples -img ./image/\$i -num 9 -bg bg.txt -vec “neg_training\$j.vec” -maxxangle 0.6 -maxyangle 0 -maxzangle 0.3 -maxidev 100 -bgcolor 0 -bgthresh 0 -w 20 -h 20
let j=j+1
done
where the positive samples are under ./image directory. For each positive hand sample a negative training sample file is create namely neg_training1.vec, neg_training2.vec etc.

Merging samples:
As seen above the createsamples utility superimposes only positive sample (hand) against a series of negative samples. Since we would like to repeat the utility for multiple positive samples there needs to be a way for merging all the training samples. A useful utility (mergevec.cpp) has been created  and is available for download in Natoshi Seo’s blog
Download mergevec.cpp and build it along with (cvboost.cpp, cvhaarclassifier.cpp, cvhaartarining.cpp, cvsamples.cpp, cvcommon.cpp) with usual include files and the link libraries.
Once built it can be executed by executing
mergevec pos_neg.txt pos_neg.vec – w 20 -h 20
where pos_neg.txt will contain both the positive and negative training sample files as follows
pos_neg.txt
./vec/pos.vec
./vec/neg_training1.vec
./vec/neg_training2.vec
….

As before you can verify that the entire training file pos_neg.vec is sensible by executing
createsamples -vec pos_neg.vec -w 20 -h 20
Now the pos_neg.vec will contain all the training samples that are required for the haartraining process.

HaarTraining
The haartraining can be run with the training samples generated from the mergevec utility described above.
The command is
./haartrainer -data haarcascade -vec pos_neg.vec -bg bg.txt -nstages 20 -nsplits 2 -minhitrate 0.999 -maxfalsealarm 0.5 -npos 7 -nneg 9 -w 20 -h 20 -nonsym -mem 512 -mode ALL
Several posts have suggested that nstages should be ideally 20 and splits should be 2.
npos indicates the number of positive samples and -nneg the number of negative samples. In my case I had just used 7 positive and 9 negative samples. This step is extremely CPU intensive and can take several hours/days to complete. I had reduced the number of stages to 14. The haartraining utility will create a haarcascade directory and a haarcascade.xml.

Performance testing : The first way to test the integrity of the haarcascades is to run the performance utility described in my earlier post OpenCV:Haartraining and all that jazz! If you want to use the performance utility you should also create test samples which can be used for testing with the command

createsamples -img ./image/hand-1.bmp -num 10 -bg bg.txt -info test.dat -maxxangle 0.6 -maxyangle 0 -maxzangle 0.3 -maxidev 100 -bgcolor 0 -bgthresh 0
The test.dat can be used with the performance utility as follows
./performance -data haarcascade.xml -info test.dat -w 20 -h 20

I wanted something that would be more visually satisfying that seeing the output of the  performance utility. This utility just spews out textual information of hits, misses.
So I decided to use the facedetect.c with the haarcascade trained by my positive and negative samples to check whether it was working.  The test below describes using the facedetect.c for detecting the hand.

Compile and link the facedetect.c as handetect with the usual suspects. Once this builds successfully you can use

where hands.jpg was my test image. If the handdetect does indeed detect the hand in the image it will l enclose the object detected with a rectangle. See the output below

While my haarcascade does detect 5 hands it does appear shifted. This could be partly because the number of positive and negative training samples used were 7 & 9 respectively which is very low. Also I used only 14 stages in the haartraining. As mentioned in the beginning there is a need for at least 1000 positive and 2000 negative samples which has to be used. Also haartraining should have 20 stages and 2 nsplits.  Hopefully if you follow this you would have developed a fairly true haarcascade.
If you are adventurous enough you could run the above with the webcam as
where 0 indicates the webcam. Assuming that your haarcascade is perfect you should be able to track your hand in real time.
Haarpy training!
– OpenCV: Haartraing and all that jazz!
– 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

OpenCV: Haartraining and all that jazz!

Object detection in OpenCV can be done through Haartraining. OpenCV haartraining applications provide the ability to detect objects of interest to us like faces, eyes, moving cars etc. What is required is that we need to identify positive and negative samples and use the utilities that OpenCV provides us to train the application to be able to recognize the objects we want. The positive and negative samples are used to create classifiers which are then utilized to detect the objects that we intend to.

Fortunately the tar file (OpenCV-2.3.1a.tar.bz2) which can be downloaded from the OpenCV site already comes bundled with all the necessary utilities to create a fully trained haar classifier. We just to have to build and run the commands to create our our classifiers.

This post will look at the necessary steps for creating a haarclassifier. To get started with OpenCV please look at my earlier post Computer Vision: Getting started with OpenCV.

Once you have installed OpenCV look under modules/haartraining. All the necessary files are included.

There are 3 steps in this process

1) Create samples   2) Train and create a classified haar 3) Performance testing

Create Samples (createsamples.cpp) : This is the first utility that has to be executed. This utility takes as input positive samples (images of the object that we are interested in) and negative samples (images that do not include the object that we are interested in). The createsamples utility superimposes the objects that we want to recognize in various degrees of rotation against the background of the negative samples. The composite images file is then used to train the haar application.

The first step is to build createsamples.cpp. Make sure you include ../OpenCV-2.3./modules/haartraining. Build the following files

(createsamples.cpp, cvboost.cpp, cvhaarclassifier.cpp, cvhaartarining.cpp, cvsamples.cpp, cvcommon.cpp)

Once createsamples.cpp successfully builds we can create samples required  for the training.

The command is

\$./myhaartraining -img logo.png -vec samples1.jpg -bg bg.txt  -w 20 -h 20

Where I have chosen the OpenCV logo as my object of interest.

The samples are created in samples1.jpg

The bg.txt contains a list of the negative samples included as below

./img/airplane.jpg

./img/baboon.jpg

./img/kid.jpg

./img/lena.jpg

-w stands for the width and -h stands for the height of the samples/

This will create superimposed positive samples in the negative sample background. A few of these a are shown below

Make sure that the width and height are small i.e < 50 otherwise the haartraining application core dumps because of lack of memory.

To check if the samples are created properly run a test round as follows run the command with the following options

\$./myhaartraining -img logo.png -vec samples1.jpg -bg bg.txt -n 10 – show  -w 20 -h 20

where -n is the number of samples to be generated (default is 1000) -show will show a series of images with the positive samples superimposed on the negative samples. This can be used to check if createsamples utility is working properly. For a more thorough and detailed explanation see my post Hand detection through Haartraining: A hands-on approach

Haartrainer (haartrainer.cpp): This utility takes as input the samples from the createsamples utility and creates a trained haar classifier. To build the haartrained files, build the haartraining.cpp. As before make sure you include the appropriate files along with the all the opencv libraries. Buid the following files

(haartraining.cpp, cvboost.cpp, cvhaarclassifier.cpp, cvhaartarining.cpp, cvsamples.cpp, cvcommon.cpp)

The command to use is

./haartrainer -data test2 -vec samples1.jpg -bg bg.txt -npos 1 -nneg 4 -nstages 20 -mem 500 -w 20 -h 20

In the above command test2 is the directory name in which the trained classifier is stored. The -vec option denotes the samples that were captured by the createsamples utility above. The bg.txt contains the negative samples file. The width and height have to be the same as used in the create samples utility. As mentioned before if the width and height are too large you the haartrainer will bail out with “Insufficient memory”

If the haartrainer executes successfully the test2 directory will have all the trained files and the directory will also have the haarclassifier as a test2.xml.

Once these two steps go through successfully we have to run the performance step

Performance (performance.cpp): This step is run to ensure that we have trained our application to properly recognize the object we intended to. As before the necessary file to build is performance.cpp. The files to build are

(performance.cpp, cvboost.cpp, cvhaarclassifier.cpp, cvhaartarining.cpp, cvsamples.cpp, cvcommon.cpp)

Once this built successfuly you can now run using the command

./performance -data test2.xml -info bg.txt -w 20 -h 20

Make sure that the bg.txt contains the images from which the object has to be detected in the following

format

[positive filename] [# of objects] [[x y width height]

bg.txt

/img/logo5.jpg 1 145 100 20 20

./img/logo6.jpg 1 145 100 20 20

./img/logo3.jpg 1 145 100 45 45

./img/logo4.jpg 1 145 100 45 45

./img/airplane.jpg 1 145 100 45 45

./img/baboon.jpg 1 145 100 45 45

./img/kid.jpg 1 145 100 45 45

./img/lena.jpg 1 145 100 45 45

./img/opencv-logo2.png 1 145 100 35 35

./img/logo.png 1 145 100 45 45

When this is run we get

ganesh@localhost Debug]\$ ./performance -data test3.xml -info bg.txt -w 20 -h 20

+================================+======+======+======+

|            File Name           | Hits |Missed| False|

+================================+======+======+======+

|                 ./img/logo5.jpg|     0|     1|    43|

+——————————–+——+——+——+

|                 ./img/logo6.jpg|     0|     1|    51|

+——————————–+——+——+——+

|                 ./img/logo3.jpg|     1|     0|    37|

+——————————–+——+——+——+

|                 ./img/logo4.jpg|     0|     1|     7|

+——————————–+——+——+——+

|              ./img/airplane.jpg|     1|     0|   226|

+——————————–+——+——+——+

|                ./img/baboon.jpg|     0|     1|   236|

+——————————–+——+——+——+

|                   ./img/kid.jpg|     0|     1|  1291|

+——————————–+——+——+——+

|                  ./img/lena.jpg|     0|     1|   188|

+——————————–+——+——+——+

|          ./img/opencv-logo2.png|     0|     1|     3|

+——————————–+——+——+——+

|                  ./img/logo.png|     0|     1|    33|

+——————————–+——+——+——+

|                           Total|     2|     8|  2115|

+================================+======+======+======+

As can be seen the training has been too accurate. There are hits and misses along with a false positive. The number of positive and negative samples, the co-ordinates and the number of stages all have to fine tuned to get the correct result.

Watch this space!

I will be back! Hasta la vista!

Please do take a look at my sequel to this post Hand detection through Haartraining: A Hands-on approach for a more robust haartraining method