Automatic Threshold

Threshold plays an important role in image segmentation. The threshold value in an image can be manually or automatically set. Manual threshold is not always suitable for all situation because different image may have different intensity distribution.

There are some developed automatic threshold algorithm that have been developed. In this article we are going to use the triangle algorithm, one of the simplest method that I know. It's quite old but still effective to be used under certain image conditions. The one we are going to write here is a slightly modified version from the one written in this reference:

Zack GW, Rogers WE, Latt SA (1977), "Automatic measurement of sister chromatid exchange frequency", J. Histochem. Cytochem. 25 (7): 741–53, PMID 70454

In order to follow these steps, we need to set up our compiler, in this example I use Visual Studio 2008 and OpenCV 2.1.

Ok, let's get our hands dirty, just make sure that OpenCV has been setup properly in your system:
1. Declare some IplImage objects:

 IplImage* image(0);     
IplImage* imageGrey(0);     
IplImage* imageBin(0);

2. Load the image from a file:

 image     = cvLoadImage((char*)(void*)Marshal::StringToHGlobalAnsi(openFileDialog1->FileName), 1);    
imageGrey = cvCreateImage(cvSize(image->width, image->height), 8, 1);     imageBin  = cvCreateImage(cvSize(image->width, image->height), 8, 1);

3. Convert the color image into grayscale image, we're going to perform the threshold in grayscale value. We can also perform threshold for all the channels afterwards.

cvCvtColor(image, imageGrey, CV_BGR2GRAY );

4. Create the histogram of the image, there's a trick in histogram calculation here, hope you understand the trick:

int histogramGrey[256] = {0};    
for(int row=0;rowheight;row++)    
{       
 for(int col=0;colwidth;col++)       
 {          
   CvScalar m, n;          
   m = cvGet2D(imageGrey, row, col);          
   histogramGrey[(int)m.val[0]]++;       
 }              
}

5. There's an assumption in this algorithm, the background and the foreground must be separated, so we're going to find out the peak of the background and the foreground. The simplest way is to calculate the center of mass of the histogram and find the peak in each region (left a nd right).

 float Iz  = 0;     
 float I   = 0;     
 float peak = 0;       

 for(int i=0;i<256;i++)     
 {        
   Iz  += (i * histogramGrey[i]);        
   I   += histogramGrey[i];     
 }

 peak = Iz / I;

6. Then the next step is to find the left peak (background) and the right peak (foreground). A simple peak search will do.

// finding the left peak               
for(int i=0;i<(int)peak;i++)     
{        
 if(peakL < histogramGrey[i])        
 {           
   peakL = i;         
 }              
}         
//finding the right peak      
for(int j=(int)peak;j<256;j++)      
{         
 if(peakR < histogramGrey[j])         
 {            
   peakR = j;         
 }                        
}

7. Ok, now we have a line con necting the left and the right peak. Then the next thing to do is to find a point in the histogram with the maximum distance from t he line. The code is a little bit dirty with typecasting but it works.

for(int i=peakL+dist;i<peakR-dist;i++)
{
 float distance = (((peakR - peakL)*(histogramGrey[peakL] - histogramGrey[i])) - ((peakL - i)*(histogramGrey[peakR] - histogramGrey[peakL])) )/
    sqrt( ((float)peakR - (float)peakL) * ((float)peakR - (float)peakL) + (((float)histogramGrey[peakR] - (float)histogramGrey[peakL]) * ((float)histogramGrey[peakR] - (float)histogramGrey[peakL])));

 if(minDist > distance)
 {
   autoThresholdVal = i;
 }         
}

8. Finally we got the threshold valu e and then we can call cvThreshold() function.

cvThreshold( imageGrey, imageBin, autoThresholdVal, 255, CV_THRESH_BINARY );

With a little modification, the algorithm can be implemented real time by adding it into the webcam callback.

Et Voila! Here is the screen shot of the result. The left image is the live image from a webcam and the right image is the image after performing auto threshold.

Hope it's useful!

Sine Fitting

Sine is not a linear function, but fitting series of data to sine function is actually not a difficult task. I need to apply this fitting function so I was trying to find a way to do this.

Sine function can be represented in general form as:

- y(x) = A + C * sin(x + b)

this function can be rewritten as:

- y(x) = A + C * sin(b) * cos(x) + C * cos(b) * sin(x)

or can also be written as:

- y = b0 + b1 * x1 + b2 * x2

which is a general linear regression form and can be solved by using LM function to obtain b0, b1, and b2, then based on these values we can calculate A, b and C

1. A = b0

2. b1^2 + b2^2 = c^2 * (sin^2(b) + cos^2(b)) = c^2

C = sqrt(b1^2 + b2^2)

3. tan(b) = b1/b2

b = arctan(b1/b2)

C++ implementation:

1. Prepare the buffer for datax, datay for for cos(b), sin(b), and dataz for the input data ) and variables for the result, phase and mag:

float * datax = new float[numOfData]; float * datay = new float[numOfData]; float * dataz = new float[numOfData];  float phase(0); float mag(0);

2. Then fill the data with a sine or cosine function and the simulated input data:

for(int i=0;i

3. Prepare the variables for performing least square fitting.

float a(0); float b(0); float c(0); float p[3] = {0}; // product of fitting equation float XiYi(0); float XiZi(0); float YiZi(0); float XiXi(0); float YiYi(0); float Xi(0); float Yi(0); float Zi(0);  for(int i=0;i


4. Then calculate the inverse of the matrix as follows:

float A[3][3];
float B[3][3];
float C[3][3];
float X[3][3];
int i;
int j;
float x = 0;
float n = 0; //n is the determinant of A

A[0][0] = XiXi;
A[0][1] = XiYi;
A[0][2] = Xi;
A[1][0] = XiYi;
A[1][1] = YiYi;
A[1][2] = Yi;
A[2][0] = Xi;
A[2][1] = Yi;
A[2][2] = numOfData;

for( i=0;i<3;i++)
{
 for(int j=0;j<3;j++)
 {
   B[i][j] = 0;
   C[i][j] = 0;
 }
}

for( i=0, j=0;j<3;j++)
{
 if(j == 2)
 {
   n += A[i][j] * A[i+1][0] * A[i+2][1];
 }
 else if(j == 1)
 {
   n += A[i][j] * A[i+1][j+1] * A[i+2][0];
 }
 else
 {
   n += A[i][j] * A[i+1][j+1] * A[i+2][j+2];
 }
}

for( i=2, j=0;j<3;j++)
{
 if(j == 2)
 {
   n -= A[i][j] * A[i-1][0] * A[i-2][1];
 }
 else if(j == 1)
 {
   n -= A[i][j] * A[i-1][j+1] * A[i-2][0];
 }
 else
 {
   n -= A[i][j]*A[i-1][j+1]*A[i-2][j+2];
 }
}

// Check determinant n of matrix A
if (n)
{
 x = 1.0/n;

 for(i=0;i<3;i++)
 {
   for(j=0;j<3;j++)
   {
     B[i][j] = A[j][i];
   }
 }

 C[0][0] =  (B[1][1] * B[2][2]) - (B[2][1] * B[1][2]);
 C[0][1] = ((B[1][0] * B[2][2]) - (B[2][0] * B[1][2])) * (-1);
 C[0][2] =  (B[1][0] * B[2][1]) - (B[2][0] * B[1][1]);

 C[1][0] = ((B[0][1] * B[2][2]) - (B[2][1] * B[0][2])) * (-1);
 C[1][1] =  (B[0][0] * B[2][2]) - (B[2][0] * B[0][2]);
 C[1][2] = ((B[0][0] * B[2][1]) - (B[2][0] * B[0][1])) * (-1);

 C[2][0] =  (B[0][1] * B[1][2]) - (B[1][1] * B[0][2]);
 C[2][1] = ((B[0][0] * B[1][2]) - (B[1][0] * B[0][2])) * (-1);
 C[2][2] =  (B[0][0] * B[1][1]) - (B[1][0] * B[0][1]);

 for( i=0;i<3;i++)
 {
   for( j=0;j<3;j++)
   {
     X[i][j] = C[i][j] * x;
   }
 }

 p[0] = XiZi;
 p[1] = YiZi;
 p[2] = Zi;

 a = X[0][0] * p[0] + X[0][1] * p[1] + X[0][2] * p[2];
 b = X[1][0] * p[0] + X[1][1] * p[1] + X[1][2] * p[2];
 c = X[2][0] * p[0] + X[2][1] * p[1] + X[2][2] * p[2];
}
else  // determinant=0
{
 a = 1;
 b = 1;
 c = 0; 
}

5. Then we could obtain the magnitude:

mag = sqrt(a*a + b*b);

6. Finally we could calculate the phase of the signal:

float phi;

if((b == 0) && (a >= 0))
{
 phase = M_PI/2;
}
else if((b == 0) && (a >= 0))
{
 phase = 3*M_PI/2;
}
else if ((a >= 0) && (b > 0))
{
 phi = atan(fabs(a/b));
 phase = phi;
}
else if ((a >= 0) && (b < 0))
{
 phi = atan(fabs(a/b));
 phase = M_PI - phi;
}
else if ((a <  0) && (b > 0))
{
 phi = atan(fabs(a/b));
 phase = 2 * M_PI - phi;
}
else if ((a <  0) && (b < 0))
{
 phi = atan(fabs(a/b)); 
 phase = phi + M_PI;
}

Then it's time to visualize the result...

Good luck!

Multiple Webcam using OpenCV 2.1 and Visual Studio 2008

I just found out that it is not difficult to connect to two webcam using OpenCV and do image processing on the captured images. I used OpenCV 2.1 and Visual C++ 2008 for this program, so those who are familiar with Visual C++ 2008 should be familiar with this code. The installation and configuration of OpenCV is explained clearly in their wiki and you can follow the steps written there.

We need to create two capture objects, one for each camera:

CvCapture* capture;
CvCapture* capture2;

Also, create two IplImage objects for the two cameras:

IplImage* image;
IplImage* image2;

The CvCapture and IplImage objects have to be initialized before use, add a button to your form and add these initialization lines:

int w = 640;
int h = 480;

// creating the capture fwom webcam #1
capture = cvCreateCameraCapture(0);
capture2 = cvCreateCameraCapture(1);

// parameter setting
cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_WIDTH, w);
cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_HEIGHT, h);

cvSetCaptureProperty(capture2, CV_CAP_PROP_FRAME_WIDTH, w);
cvSetCaptureProperty(capture2, CV_CAP_PROP_FRAME_HEIGHT, h);

// initialization of iplimage object
image  = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 3);
image2 = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 3);

Then the next thing is to declare a function for capturing the images:

private: System::Void ProcessFrame(System::Object^ sender, System::EventArgs^ e)
{              // query the frame and capture
image = cvQueryFrame(capture);
image2 = cvQueryFrame(capture2);

cvShowImage("camera 1", image);
cvShowImage("camera 2", image2);
  
if(blnGrabImage)
{
cvSaveImage("image_left.bmp", image);
cvSaveImage("image_right.bmp", image2);
blnGrabImage = false;
}
}

We have to register this function to be called by the system, add these lines after the initialization of the capture and the images:

Application::Idle += gcnew EventHandler(this, &OpenCVImage::Form1::ProcessFrame);

There are two important parts in this line, the first one is

Application::Idle += gcnew EventHandler

which tell the system to add an event handler to the system when the system is idle and the other thing is this part:

&OpenCVImage::Form1::ProcessFrame

which tells the system to call this function every time there's a new event in the system, in this case OpenCVImage is the name of the project, Form1 is the name of the Form, and ProcessFrame is the function that is going to be called.

So, that's it, tidy up the code a little bit and ready to run. If you're successful you'll get something like this:
1. For single camera:
2. For multiple cameras:
- left camera:
- right camera:
Left and right are seen from the objects' viewpoint facing the camera...

Good luck!

Playing with webcam in OpenCV 2.1

I have posted an article about capturing and processing image using webcam in OpenCV 1.0. Right now I'm trying to upgrade my program to use OpenCV 2.1. Some functions are different but basically most of them are similar.

I'm using Microsoft Visual C++ 2008 as the compiler and Windows Form Application. It's quite straight away. This is an example of the user interface:

Basically, there are many ways of doing this, the one that I choose is using event trigger approach. So the image will be shown every time there's a new event from the camera received by the system.

Here's the details, I assume that you are already familiar with C++, Visual C++ 2008, and installing OpenCV 2.1 in C++ 2008:

1. include all the needed include files:
   

   #include < cv.h >
   #include < cxcore.h >
   #include < highgui.h >
   

2. Declare the needed objects, to make it simple just declare them as global variables:

   // object for capturing image   
   CvCapture*  capture;

   // iplimage objects   
   IplImage*  image;
   IplImage*  procimage;
   

   // image size
   double w;
   double h;
   

   // boolean variables for flag
   bool blnCapInProgress;
   bool blnGrabImage;
   

3. Initialize the capture when the start button is clicked.
   

4. // size of the image   
   w = 640;
   h = 480;
   

   // creating the capture object, connecting it to the camera, and set the size   
   capture = cvCreateCameraCapture(0);
   cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_WIDTH, w);
   cvSetCaptureProperty(capture, CV_CAP_PROP_FRAME_HEIGHT, h);
   

   // initialize the images and windows   
   image      = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 3);
   procimage = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
   cvNamedWindow("continuous", CV_WINDOW_AUTOSIZE);
   cvNamedWindow("gray", CV_WINDOW_AUTOSIZE);
   

   // declaration of event handle  
   if(capture != NULL)
   {
   this->label1->Text = "ok";
   
   if (blnCapInProgress)
   {
   // stop the capture     // the capture is stopped by removing the event handler..
   Application::Idle -= gcnew EventHandler(this, &OpenCVImage::Form1::ProcessFrame);
   buttonStart->Text = "Start Capture";
   }
   else
   {
   //start the capture     // the capture is stopped by adding the event handler..
   Application::Idle += gcnew EventHandler(this, &OpenCVImage::Form1::ProcessFrame);
   buttonStart->Text = "Stop";
   }
   
   blnCapInProgress = !blnCapInProgress;        
   }
   

   
   
5. The event handler is inside the ProcessFrame function, create a function like this:

   private: System::Void ProcessFrame(System::Object^ sender, System::EventArgs^ e) {
   // query the image from the camera   
   image = cvQueryFrame(capture);
   // show the image
   if(image != NULL)
   {                   // image processing algorithm can be performed here
   cvCvtColor(image, procimage, CV_RGB2GRAY);
   cvShowImage("continuous", image);
   cvShowImage("gray", procimage);
   }
   // grab the image
   if(blnGrabImage)
   {
   cvSaveImage("image.bmp", image);
   cvShowImage("grab image", image);
   blnGrabImage = false;
   }
   }
   

6. You can grab the image simply by setting the blnGrabImage variable to true.
7. And what you can get is something like this... :D

Plane Fitting using OpenCV

As what I've promised, I'm going to write about plane fitting using OpenCV. It's quite straight away and quite similar to my previous post, the difference is only about the physical representation of the variables in the equations. This time it's only a plane fitting, so it's a linear least square fitting. Later I would like to explore non-linear fitting as well and share it here when I have made it.

Let's say that we have a set of data that represents a plane in 3D coordinate X-Y-Z and modeled as Axi + Byi + C = zi. where i is the enumerator of the data (0~n).

Using the same approach, we can use matrix to represent the simultaneous equation and solve it.I will directly show you how.

1. prepare the matrices for the data.

CvMat *res  = cvCreateMat(3, 1, CV_32FC1);
CvMat *matX = cvCreateMat(10000, 3, CV_32FC1);
CvMat *matZ = cvCreateMat(10000, 1, CV_32FC1);

2. input the data into the matrices.

for(int row=0;row<100;row++)
{
 for(int col=0;col<100;col++)
 {
   int idx     = row * col;
   double val  = cvmGet(filteredMat, row, col);

   cvmSet(matX, idx, 0, col*pixScale);
   cvmSet(matX, idx, 1, row*pixScale);
   cvmSet(matX, idx, 2, 1);
   cvmSet(matZ, idx, 0, val);
 }
}

3. solve the equation

cvSolve(matX, matZ, res, CV_SVD);
A = cvmGet(res, 0, 0);
B = cvmGet(res, 1, 0);
C = cvmGet(res, 2, 0);

4. calculate the distance between the plane and the real data.

double sqrtc = sqrt(pow(A, 2) + pow(B, 2) + pow(C, 2));

float min      = 100;
float max      = -100;

// plane generation, Z = Ax + By + C
// distance calculation d = Ax + By - z + C / sqrt(A^2 + B^2 + C^2)
// pointlist generation
for(int row=0;row<filteredMat->rows;row++)
{
 for(int col=0;col<filteredMat->cols;col++)
 {
   double Zval     = cvmGet(filteredMat, row, col);
   double val      = col*pixScale*A + row*pixScale*B - Zval + C;
   double distance = val / sqrtc;

   if(min > distance)
     min = distance;

   if(max < distance)
     max = distance;

   cvmSet(planeMat, row, col, distance);
 }
}

Here is an example of this plane fitting, this is the surface before fitting, this is measurement result of a flat surface using white light interferometry:

and after fitting the result is:

The color represents the height, so we can see that after fitting the color of the surface is almost the same only a small variation of the surface. The sample should be flat and smooth, however it's not smooth anymore due to wear. However, the real scale is nanometer, the variation that we see in the result is around 100 nm.

Linear Least Square calculation using OpenCV

Least Square method is an important tool for fitting either a plane or surface.This article will explain briefly on how we can use OpenCV for performing Linear Least Square for solving simultaneous equation. This technique can be applied for line or surface fitting.

As an example, an equation: Z = alpha * X + beta + h * delta, which has three unknown variables alpha, beta and h. The known variables are Z, X, and delta.

We have to create a matrix for the known variables according to the number of available data.

CvMat * matZ = cvCreateMat((dataLength), 1, CV_32FC1);
CvMat * matX = cvCreateMat((dataLength), 3, CV_32FC1);
CvMat * matRes = cvCreateMat(3, 1, CV_32FC1);

Then the data can be put into the matrix, assume that the data for Z are obtained from an array data:

//fill data
for(int i=0;i<dataLength;i++)
{
 cvmSet(matZ, i, 0, data[i]);
 cvmSet(matX, i, 0, (i*OBJSCALE));
 cvmSet(matX, i, 1, 1);
 cvmSet(matX, i, 2, 1);
}

after that, call the least square solver using CvSolve and if nothing is wrong with the data, the result can be directly accessed by accessing the matRes matrix

cvSolve(matX, matZ, matRes, CV_SVD);
double h     = cvmGet(matRes, 2, 0);
double beta  = cvmGet(matRes, 1, 0);
double alpha = cvmGet(matRes, 0, 0);

Finally, clean the memory of the matrices:

cvReleaseMat(&matZ);
cvReleaseMat(&matX);
cvReleaseMat(&matRes);

That's it. Hope it helps :)

(the source code is formatted using http://formatmysourcecode.blogspot.com/, nice stuff!)

haven't written any code yet...

It has been almost three months in 2011, I haven't written any serious C++ code yet. Recently I've been fixing memory leak errors in a software, playing with some openGL stuffs and lots of matlab stuffs.

I'm planning to write something about simple image processing algorithms and examples with a webcam, the easiest and cheapest way of giving real time examples (which is quite impressive and efficient).

Another thing is about plane fitting, I'm planning to use OpenCV matrix library for doing plane fitting. This will be very useful for me and hopefully for all of you too. I use matlab to perform this I can imagine how to convert this into C++, what I need is some time and a clear mind. :)

So, stay tune...