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.

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 DFT(blurred image) * complex conjugate of blur kernel

      cvMulSpectrums(dft_A,dft_B,dft_C,CV_DXT_MUL_CONJ);

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

6.  Finally perform Inverse DFT

  // Perform Inverse

  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.


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

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.

  

1) Download OpenCV2.2  from Source Forge

Download the OpenCV-2.2.0-win32-vs2010.exe as it already contains the libraries and the binaries required for running OpenCV.

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

2) Download and install Visual Studio 2010 express (free) version and install on your PC.

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

4) Include any relevant OpenCV code in your first.cpp ( e.g. snippet given below)

5) Include the following 3 directories under

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

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) Now you are good to go. Build by choosing Debug->Build Solution

It should through  fine when you 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

    IplImage* img = cvLoadImage("baboon.jpg",1);

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

}

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 1^st order or 2^nd 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.

A code snippet (Complete code 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 2^nd 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* img = cvLoadImage("gazelle.jpg",1);
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 and 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_8uc1 = cvLoadImage( argv[1], CV_LOAD_IMAGE_GRAYSCALE );

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; itotal; ++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;

}

http://www.eecs.berkeley.edu/~trevor/CS280Notes/