My References: Visions Image Checker Development

Thursday, August 28, 2014

Pattern Recognition, Find a Mark in an Image

To find a specific mark in an image we use statistical correlation to compare a reference image that defines the mark to a test image that may contain the mark. The result image is filled with the correlation coefficient data. The pixel in the result image with the highest value is the most probable location for the mark.

This procedure supports grayscale images. Use the Victor Library conversion functions to convert color or black and white images to 8-bit grayscale. In general, the steps are:

determine the dimensions and pixel depth of the test and reference images
allocate buffer space for the images
load the images
fill the result image with correlation coefficient data calculated between the test image and the reference image
sort the correlation coefficient data to determine the most probable location of the mark


Test image	Reference image	Result image, mark found at (84, 213)

The result image is filled with pixel data that represent the correlation coefficient calculated between the reference image and the test image at each pixel location. The pixel values can range from zero (negative correlation) through 128 (no correlation) up to 255 (perfect positive correlation). A negative correlation would mean a good match with the negative image of the reference mark.
The brightest pixel in the result image represents the most probable location of the reference mark. In this case that is at location (84, 213).

Find a Mark - the VB.NET Source Code

Requires Victor Image Processing Library for 32-bit Windows v 6.0 or higher.

// Add viclib as a reference in your VB.NET project
   Private Function findmark(ByVal srcimg As vicwin.imgdes, ByVal oprimg As vicwin.imgdes, ByVal resimg As vicwin.imgdes, ByVal ptarray As vicwin.COORD_VAL(), ByVal numelem As Long) As Long
        Dim rcode As Long

        '   Fill the result image with correlation coefficient data
        rcode = vicwin.correlateimages(srcimg, oprimg, resimg)
        If (rcode = vicwin.NO_ERROR) Then
            '  The pixel data in the result image will be sorted
            '  and placed in the array in descending order
            rcode = vicwin.sortpixelsbyval(resimg, ptarray(0), numelem)

        End If

        findmark = rcode

    End Function


    Private Sub mnu_findmark_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles mnu_findmark.Click
        Dim rcode As Long
        Dim testimage As vicwin.imgdes
        Dim markimage As vicwin.imgdes
        Dim testfile As vicwin.JpegData
        Dim markfile As vicwin.JpegData
        Dim scols, srows, ocols, orows As Long
        Dim numelem As Long

        ' For easier understanding of the sample code no error checking is performed here
        ' but recognize that you should always test the return code to verify that each function is successful 
        rcode = vicwin.jpeginfo("testc.jpg", testfile)
        rcode = vicwin.allocimage(testimage, testfile.width, testfile.length, testfile.vbitcount)
        rcode = vicwin.loadjpg("testc.jpg", testimage)

        rcode = vicwin.jpeginfo("register.jpg", markfile)
        rcode = vicwin.allocimage(markimage, markfile.width, markfile.length, markfile.vbitcount)
        rcode = vicwin.loadjpg("register.jpg", markimage)

        scols = testimage.endx - testimage.stx + 1
        srows = testimage.endy - testimage.sty + 1
        ocols = markimage.endx - markimage.stx + 1
        orows = markimage.endy - markimage.sty + 1


        ' Allocate an array to hold the pixel data for the sorting
        numelem = scols * srows
        Dim ptarray(numelem) As vicwin.COORD_VAL

        rcode = findmark(testimage, markimage, testimage, ptarray, numelem)
        If (rcode = vicwin.NO_ERROR) Then
            MsgBox("Most probable location for mark: " & ptarray(0).val & "  (" & ptarray(0).x & "," & ptarray(0).y & ")")
        End If

        rcode = vicwin.savejpg("testr.jpg", testimage, 50)
        vicwin.freeimage(markimage)
        vicwin.freeimage(testimage)
    End Sub

Find a Mark - the C# Source Code

Requires Victor Image Processing Library for 32-bit Windows v 6.0 or higher.

                case "findmark":   // Assumes three existing images:
                                   // the subject image is in simage, the register mark is in oimage, and the result will be placed in rimage.

                    {   // Fill the result image with correlation coefficient data
                        rcode = vicwin.correlateimages(ref simage, ref oimage, ref rimage);
                        if (rcode == vicwin.NO_ERROR)
                        {
                            int numelem;
                            string message;
                            double coef;
                            int rwidth = 0;
                            int rlength = 0;

                            rwidth = rimage.endx - rimage.stx + 1;
                            rlength = rimage.endy - rimage.sty + 1;
                            numelem = rwidth * rlength;
                            
                            // Allocate an array to hold the pixel data for the sorting
                            vicwin.COORD_VAL[] ptarray = new vicwin.COORD_VAL[numelem];
                            
                            // The pixel data in the result image will be sorted
                            // and placed in the array in descending order

                            rcode = vicwin.sortpixelsbyval(ref rimage, ref ptarray[0], numelem);
                            if (rcode == vicwin.NO_ERROR)
                            {
                                coef = (double)ptarray[0].val / 255.0;
                                if (coef > 1.0)
                                    coef = (double)ptarray[0].val / 65535.0;
                                message = "Highest correlation value: " + coef.ToString("0.0000");
                                message = message + "       Most probable location for mark: " + "  (" + ptarray[0].x + "," + ptarray[0].y + ")";
                                MsgBox(0, message, "Find mark", 0);
                            }
                        }
                        break;
                    }






Find a Mark - the Visual Basic Source Code
Requires Victor Image Processing Library for 32-bit Windows v 5.4 or higher and the helper dll vicstats.dll.
' Add the data type definition and function declarations to vicdef32.bas ...............................................

Type COORD_VAL
   val As Long
   x As Long
   y As Long
End Type

Public Declare Function correlationcoef Lib "vicstats.dll" (ByRef srcimg As imgdes, ByRef oprimg As imgdes, ByRef pCoef As Double) As Long
Public Declare Function correlateimages Lib "vicstats.dll" (ByRef srcimg As imgdes, ByRef oprimg As imgdes, ByRef resimg As imgdes) As Long
Public Declare Function calcavglevelfloat Lib "vicstats.dll" (ByRef srcimg As imgdes, ByRef redavg As Double, ByRef grnavg As Double, ByRef bluavg As Double) As Long
Public Declare Function sortpixelsbyval Lib "vicstats.dll" (ByRef srcimg As imgdes, ByRef first_elem As COORD_VAL, ByVal nelem As Long) As Long

' The Functions .................................................... 
Private Function findmark(srcimg As imgdes, oprimg As imgdes, resimg As imgdes, ptarray As COORD_VAL, numelem As Long) As Long
Dim rcode As Long

'   Fill the result image with correlation coefficient data
    rcode = correlateimages(srcimg, oprimg, resimg)
                                                      
   If (rcode = NO_ERROR) Then
        '  The pixel data in the result image will be sorted
        '  and placed in the array in descending order
       rcode = sortpixelsbyval(resimg, ptarray, numelem)

   End If
    
   findmark = rcode

End Function


Private Sub mnufindmark_Click()
Dim rcode As Long
Dim testimage As imgdes
Dim markimage As imgdes
Dim testfile As JpegData
Dim markfile As JpegData
Dim scols, srows, ocols, orows As Long

Dim ptarray() As COORD_VAL
Dim numelem As Long

   ' For easier understanding of the sample code no error checking is performed here
   ' but recognize that you should always test the return code to verify that each function is successful 
   rcode = jpeginfo("testc.jpg", testfile)
   rcode = allocimage(testimage, testfile.width, testfile.length, testfile.vbitcount)
   rcode = loadjpg("testc.jpg", testimage)

   rcode = jpeginfo("register.jpg", markfile)
   rcode = allocimage(markimage, markfile.width, markfile.length, markfile.vbitcount)
   rcode = loadjpg("register.jpg", markimage)

   scols = testimage.endx - testimage.stx + 1
   srows = testimage.endy - testimage.sty + 1
   ocols = markimage.endx - markimage.stx + 1
   orows = markimage.endy - markimage.sty + 1

   numelem = scols * srows

   ' Allocate an array to hold the pixel data for the sorting
   ReDim ptarray(0 To numelem - 1) As COORD_VAL
   
   rcode = findmark(testimage, markimage, testimage, ptarray(0), numelem)
   If (rcode = NO_ERROR) Then
      MsgBox ("Most probable location for mark: " & ptarray(0).val & "  (" & ptarray(0).x & "," & ptarray(0).y & ")")
   End If

   rcode = savejpg("testr.jpg", testimage, 50)
   freeimage markimage
   freeimage testimage

End Sub






Find a Mark - the C Source Code
Requires Victor Image Processing Library for 32-bit Windows v 5.4 or higher and the helper dll vicstats.dll.
// Add these declarations to vicdefs.h
typedef struct { unsigned val, x, y; } COORD_VAL;
int WINAPI correlationcoef(imgdes *srcimg, imgdes *oprimg, double *pCoef);
int WINAPI correlateimages(imgdes *srcimg, imgdes *oprimg, imgdes *resimg);
int WINAPI calcavglevelfloat(imgdes *srcimg, double *redavg, double *grnavg, double *bluavg);
int WINAPI sortpixelsbyval(imgdes *srcimg, COORD_VAL *darray, int nelem);

int findmark(imgdes *srcimg, imgdes *oprimg, imgdes *resimg, COORD_VAL *ptarray, int numelem)  
{
   int rcode;

   // Fills the result image with correlation coefficient data
   rcode = correlateimages(srcimg, oprimg, resimg); 
                                                      
   if(rcode == NO_ERROR) {
      // The pixel data in the result image will be sorted
      // and placed in the array in descending order
      rcode = sortpixelsbyval(resimg, ptarray, numelem);
      }
   return (rcode);
}

void DoMenuFindmark(HWND hWnd)
{
   int rcode;
   imgdes testimage;
   imgdes markimage;
   JpegData testfile;
   JpegData markfile;
   int scols, srows, ocols, orows;

   COORD_VAL *ptarray;
   int numelem;

   // For easier understanding of the sample code no error checking is performed here
   // but recognize that you should always test the return code to verify that each function is successful 
   rcode = jpeginfo("testc.jpg", &testfile);
   rcode = allocimage(&testimage, testfile.width, testfile.length, testfile.vbitcount);
   rcode = loadjpg("testc.jpg", &testimage);

   rcode = jpeginfo("register.jpg", &markfile);
   rcode = allocimage(&markimage, markfile.width, markfile.length, markfile.vbitcount);
   rcode = loadjpg("register.jpg", &markimage);

   scols = CALC_WIDTH(&testimage);
   srows = CALC_HEIGHT(&testimage);
   ocols = CALC_WIDTH(&markimage);
   orows = CALC_HEIGHT(&markimage);

   numelem =  scols * srows;

   // Allocate an array to hold the pixel data for the sorting
   ptarray = calloc(numelem, sizeof(COORD_VAL));
   rcode = findmark(&testimage, &markimage, &testimage, ptarray, numelem);
   if(rcode == NO_ERROR) {
      char szBuff[128];

      // Display the values for the four most likely locations
      wsprintf(szBuff, "Most probable locations for mark:"
         "\n%3d: (%3d, %3d)"
         "\n%3d: (%3d, %3d)"
         "\n%3d: (%3d, %3d)"
         "\n%3d: (%3d, %3d)",
          ptarray[0].val, ptarray[0].x, ptarray[0].y, 
          ptarray[1].val, ptarray[1].x, ptarray[1].y,
          ptarray[2].val, ptarray[2].x, ptarray[2].y,
          ptarray[3].val, ptarray[3].x, ptarray[3].y
          );
      MessageBox(hWnd, szBuff, szAppName, MB_OK);
      }

   if(ptarray)
      free(ptarray);

   savejpg("testr.jpg", &testimage, 50);
   freeimage(&markimage);
   freeimage(&testimage);
}

Image processing C# tutorial 4 – Gaussian blur

No, we are not off edge-detection topic yet. However, in order to introduce some fancier algorithms, we have to do some ground work first. And to keep every step fruitful, we’ll start with Gaussian blur.
The math (as in code)
I don’t want to go into much details so here’s the brief version of Gaussian blur introduction – it blurs an image by adding weighted grayscales of surrounding pixels on to the pixel to be blurred. The center pixel will have biggest weight, and the surrounding pixels will have lesser and lesser weight as their distances to the center increase. The table of weights are calculated using Gaussian function, which is

in two dimensions. The standard deviation

controls how the blur the image will get. x,y are coordinates relative to the center pixel. For example, you have a 3X3 picture will grayscale values as shown in the following table:

200	200	100
130	200	120
30	130	30

The weight table when

=1 (or a sampled Gaussian kernel in fancier terms) is calculated:

G(-1,-1)=0.0585498	G(0,-1)=0.0965324	G(1,-1)=0.0585498
G(-1,0)=0.0965324	G(0,0)=0.159155	G(1,0)=0.0965324
G(-1,1)=0.0585498	G(0,1)=0.0965324	G(1,1)=0.0585498

Notice that the table (kernel) is constant for given size, so you don’t need to re-calculate it for each pixel. You may also notice that lots of cells hold the same value in above table. This is because Gaussian function is circularly symmetric. Of course I can’t get away without showing the pretty picture of a 2D Gaussian function:

What we need to do next is to apply the G(x,y) values to grayscale values and add them up as the new grayscale value at (0,0):
newValue = 200*G(-1,-1) + 200*G(0,-1) + 100*G(1,-1) + 130*G(-1,0) + …… 30*G(1,1)
Basically the function takes the pixel itself (at center) and its 8 neighbors and mix them together as the new value of center. The original center pixel has the most effect on the result because it has biggest weight, however the neighbors also get a share in the final result. Repeat the operation of each pixel of the image then you’ll get a blurred image.
Another thing to mention is that G(x,y) can be split into two steps by doing the 1D G(x) first and then applying 1D G(y) on the result. 1D Gaussian function looks like this:

Applying G(x) and G(y) separately saves lots of calculation because 1D Gaussian kernel is a 1D array, not a nXn matrix. In the code below we’ll use the 2-pass solution.
Also, in the above example we chose a 3X3 matrix, however it’s not enough. In theory the value of Gaussian function never falls down to zero so to really consider all surrounding pixels we need a infinitely large matrix. However, as Gaussian function decreases rapidly – in other worlds the G(x,y) values will be close to zero as the distance goes further from the center. Usually 3

can be considered a safe cut-off point, so for

=1, what we need is a 7X7 matrix (3

on each side plus center).
At last, as the values in the kernel doesn’t add up to 1, the result image will be darkened or brightened. To fix this the values in the matrix need to be normalized. This can be done by dividing each value by the sum of all values in the matrix. The following sample code implements such normalization.
The code
The Gaussian class

 public class Gaussian
{
        public static double[,] Calculate1DSampleKernel(double deviation, int size)
        {
            double[,] ret = new double[size, 1];
            double sum = 0;
            int half = size / 2;
            for (int i = 0; i < size; i++)
            {
                    ret[i, 0] = 1 / (Math.Sqrt(2 * Math.PI) * deviation) * Math.Exp(-(i - half) * (i - half) / (2 * deviation * deviation));
                    sum += ret[i, 0];
            }
            return ret;
        }
        public static double[,] Calculate1DSampleKernel(double deviation)
        {
            int size = (int)Math.Ceiling(deviation * 3) * 2 + 1;
            return Calculate1DSampleKernel(deviation, size);
        }
        public static double[,] CalculateNormalized1DSampleKernel(double deviation)
        {
            return NormalizeMatrix(Calculate1DSampleKernel(deviation));
        }
        public static double[,] NormalizeMatrix(double[,] matrix)
        {
            double[,] ret = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double sum = 0;
            for (int i = 0; i < ret.GetLength(0); i++)
            {
                for (int j = 0; j < ret.GetLength(1); j++)
                    sum += matrix[i,j];
            }
            if (sum != 0)
            {
                for (int i = 0; i < ret.GetLength(0); i++)
                {
                    for (int j = 0; j < ret.GetLength(1); j++)
                        ret[i, j] = matrix[i,j] / sum;
                }
            }
            return ret;
        }
}

The class defines a series of overloaded method to calculate a sample kernel base on given standard deviation. For example, to get a kernel for deviation 1, call:

double[,] ret = Gaussian.CalculateNormalized1DSampleKernel(1);

The returned array is a 7X1 array containing the kernel.
Then the code to apply the kernel is like this (also part of Gaussian class):

 public static double[,] GaussianConvolution(double[,] matrix, double deviation)
        {
            double[,] kernel = CalculateNormalized1DSampleKernel(deviation);
            double[,] res1 = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] res2 = new double[matrix.GetLength(0), matrix.GetLength(1)];
            //x-direction
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                    res1[i, j] = processPoint(matrix, i, j, kernel, 0);
            }
            //y-direction
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                    res2[i, j] = processPoint(res1, i, j, kernel, 1);
            }
            return res2;
        }
        private static double processPoint(double[,] matrix, int x, int y, double[,] kernel, int direction)
        {
            double res = 0;
            int half = kernel.GetLength(0) / 2;
            for (int i = 0; i < kernel.GetLength(0); i++)
            {
                int cox = direction == 0 ? x + i - half : x;
                int coy = direction == 1 ? y + i - half : y;
                if (cox >= 0 && cox < matrix.GetLength(0) && coy >= 0 && coy < matrix.GetLength(1))
                {
                    res += matrix[cox, coy] * kernel[i, 0];
                }
            }
            return res;
        }

At last, the code to process the Bitmap type:

private Color grayscale(Color cr)
{
    return Color.FromArgb(cr.A, (int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11),
       (int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11),
      (int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11));
}
public override Bitmap FilterProcessImage(double d, Bitmap image)
{
            Bitmap ret = new Bitmap(image.Width, image.Height);
            double[,] matrix = new double[image.Width, image.Height];
            for (int i = 0; i < image.Width; i++)
            {
                for (int j = 0; j < image.Height; j++)
                    matrix[i, j] = grayscale(image.GetPixel(i, j)).R;
            }
            matrix = Gaussian.GaussianConvolution(matrix, d);
            for (int i = 0; i < image.Width; i++)
            {
                for (int j = 0; j < image.Height; j++)
                {
                    int val = (int) Math.Min(255, matrix[i,j]);
                    ret.SetPixel(i, j, Color.FromArgb(255,val,val,val));
                }
            }
            return ret;
}

The result
The following is the screen shot of my toolkit. To the left is the original image. In the middle is the blurred (smoothed) image – you can’t easily tell on a natural photograph like this, though. To the right is the edge detected by Krisch detector based on the blurred image. You can already see the result is a great improvement from the previous post – the facial hair don’t cause false-positives; the branch and the arm start to take shape; eyes and nose are cleared marked, etc. However this is not the full story yet. Stay tuned for further discussions on edge detection!

And the following is an example of Gaussian blur with a big deviation (8). Note in the toolkit I implemented filters on all channels so the original color is preserved – all you need to do is to apply the filter on each channel and then combine them back together.

We can also have a little fun by applying the filter to only one channel. The left image shows a Star War-like sword and a boring red blade. Applying a Gaussian blur on red channel lights up the blade. Note here i cheated a little by multiplying Green channel by 2 to make it brighter in the result. Nothing significant here – this is just for fun ;).

Image processing C# tutorial 3 – edge detection: template-based algorithms

The math
Like derivation method, the template-based methods also try to detect edges by detecting sudden grayscale changes. As the name suggests, these methods use templates to model edges. The Sobel edge detector uses two templates, modeling changes along x-axis and y-axis. The Sobel templates are defined as following:

For each pixel on the image, both templates are applied to the 3X3 area around the pixel. The grayscale values of surrounding pixels are multiplied by corresponding matrix values and then summed together:

,where T is either Sx (for x-axis) or Sy (for y-axis). Then the absolute values of Ix and Iy are summed:

At last, I is compare with threshold and the pixel is marked as edge pixel if I is larger than threshold.
Another template-based method, Kirsch method, uses 8 templates instead of 2. The templates model 8 different possible edge directions: left-right, up-down, and upperleft-lowerright, etc. For each pixel of the image, all 8 templates are applied and maximum I is chosen to be compared with the threshold.

The code
The version (v0.2) of my image processing tool that include these two template-based methods can be downloaded here. The implementation of the two methods are very similar. Here I’m showing the Krisch code that works on the R channel of the image:

public override Bitmap FilterProcessImage(Bitmap image)
 {
   Bitmap ret = new Bitmap(image.Width, image.Height);
   for (int i = 1; i < image.Width - 1; i++)
   {
       for (int j = 1; j < image.Height - 1; j++)
       {
          Color cr = image.GetPixel(i + 1, j);
          Color cl = image.GetPixel(i - 1, j);
          Color cu = image.GetPixel(i, j - 1);
          Color cd = image.GetPixel(i, j + 1);
          Color cld = image.GetPixel(i - 1, j + 1);
          Color clu = image.GetPixel(i - 1, j - 1);
          Color crd = image.GetPixel(i + 1, j + 1);
          Color cru = image.GetPixel(i + 1, j - 1);
          int power = getMaxD(cr.R, cl.R, cu.R, cd.R, cld.R, clu.R, cru.R, crd.R);
           if (power > 50)
             ret.SetPixel(i, j, Color.Yellow);
           else
              ret.SetPixel(i, j, Color.Black);
          }
      }
      return ret;
}

And getMaxD is defined as:

private int getD(int cr, int cl, int cu, int cd, int cld, int clu, int cru, int crd, int[,] matrix)
{
   return Math.Abs(  matrix[0, 0]*clu + matrix[0, 1]*cu + matrix[0, 2]*cru
      + matrix[1, 0]*cl + matrix[1, 2]*cr
         + matrix[2, 0]*cld + matrix[2, 1]*cd + matrix[2, 2]*crd);
}
private int getMaxD(int cr, int cl, int cu, int cd, int cld, int clu, int cru, int crd)
{
   int max = int.MinValue;
   for (int i = 0; i < templates.Count; i++)
   {
      int newVal = getD(cr, cl, cu, cd, cld, clu, cru, crd, templates[i]);
      if (newVal > max)
         max = newVal;
    }
    return max;
}
private List<int[,]> templates = new List<int[,]> 
{
   new int[,] {{ -3, -3, 5 }, { -3, 0, 5 }, { -3, -3, 5 } },
   new int[,] {{ -3, 5, 5 }, { -3, 0, 5 }, { -3, -3, -3 } },
   new int[,] {{ 5, 5, 5 }, { -3, 0, -3 }, { -3, -3, -3 } },
   new int[,] {{ 5, 5, -3 }, { 5, 0, -3 }, { -3, -3, -3 } },
   new int[,] {{ 5, -3, -3 }, { 5, 0, -3 }, { 5, -3, -3 } },
   new int[,] {{ -3, -3, -3 }, { 5, 0, -3 }, { 5, 5, -3 } },
   new int[,] {{ -3, -3, -3 }, { -3, 0, -3 }, { 5, 5, 5 } },
   new int[,] {{ -3, -3, -3 }, { -3, 0, 5 }, { -3, 5, 5 } }
};

Note the templates variable contains all 8 Krisch templates. In the sample code above the threshold is hard-coded to 50.
The result
The following result is my implementation of Krisch detector with threshold 185.

So far all the 3 edge detectors we’ve implemented only consider local grayscale changes. This works fine for simpler images – as shown in the following example the result is very accurate:

However the algorithms are greatly challenged by some of natural images. In the following example the hairs of the koala are detected as edges as they do appear as sudden grayscale changes at pixel level. However, the edges around the arm and the tree branches are not detected because they are blurred in the picture, messing up our 3-by-3 local detection. To better simulate biological vision (like we have) more sophisticated methods are needed. In the coming posts I’ll introduce some of such algorithms.

Image processing C# tutorial 2 – edge detection: derivation method

The math (simplified)
An edge on a grayscale image is defined as sudden change in pixel brightness. Imagine a picture of a black square on a white background – the pixel right to the edge of the square will be pure white and one of its neighbor pixels will be pure black. When this happens we mark the black pixel as one of the edge pixel.
Of course, a natural picture rarely has such distinct feature. The edge is usually blurred by gray pixels. However the idea of detect sudden change still applies in finding the edges.
To detect such sudden change we use “the rate of change”, which is to find out how much brightness changes over a small interval – say 1 pixel. Let’s assume A(x,y) is the brightness of pixel (x,y), then the bright change “rate” at pixel (x,y) can be calculated as:

When the rate of change G is greater than a threshold then we say the pixel is an edge pixel. The implementation doesn’t consider noise. The threshold should be adjusted based on the noise level of the image, however here we’ll just pick one. (The tool allows you to dynamically change the threshold manually). Also note that the pixel at (x,y) is actually not used during the calculation for (x,y). This is because we assume the brightness levels vary linearly between the pixels.
The code
The following is a simplified implementation that works on the R channel of the image only. Note the threshold is not adaptive (fixed to 50 in the code) – we’ll discuss that in later posts.

public override Bitmap FilterProcessImage(Bitmap image)
{
    Bitmap ret = new Bitmap(image.Width, image.Height);
   for (int i = 1; i < image.Width-1; i++)
   {
      for (int j = 1; j < image.Height-1; j++)
      {
         Color cr = image.GetPixel(i+1, j);
         Color cl = image.GetPixel(i-1, j);
         Color cu = image.GetPixel(i, j-1);
         Color cd = image.GetPixel(i, j+1);
         int dx = cr.R - cl.R;
         int dy = cd.R - cu.R;
         double power = Math.Sqrt(dx * dx / 4 + dy * dy / 4);
         if (power > 14)
            ret.SetPixel(i, j, Color.Yellow);
         else
            ret.SetPixel(i, j, Color.Black);
      }
   }
   return ret;
}

The result
I started writing a image processing tool that provide a single UI for testing different image filters. Version 0.1 of the image processing test tool can be downloaded here. Note that I just started with the tool so there were some simplification/shortcuts taken and the framework is rough. I’ll keep improving the tool overtime. Following is a sample result of using the tool:
Original image:

Detected edge on grayscale image using threshold 17:

Detected edge on Red channel only using threshold 17:

Interestingly, the R channel result is visually superior than the all-channel result for this particular picture.

Reference : Image processing C# tutorial 1– find binary image boundary

Image processing C# tutorial 1– find binary image boundary

Finding the boundary of a binary image can be done in two steps 1) apply a binary erosion. The binary erosion “erodes” the edge of the image. 2) Subtract the eroded image from original image. What is left is the border. A binary erosion is to move a small structure element (a 3 by 3 image) across the target image, and the pixel on the target image is kept only if the pixels on the target image fully “cover” the structure element for the area covered by the structure element. In other words R AND S = S for any test region R and structure element S. The C# code looks like this:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace EdgeOfBinary
{
    class Program
    {
        static void Main(string[] args)
        {
            BinaryImage testImage = new BinaryImage(new int[20, 20]
                {
                    {0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0},
                    {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
                    {0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0},
                    {0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0},
                    {0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0},
                    {0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0},
                    {0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0}
                });
            BinaryImage structureElement = new BinaryImage(new int[3, 3]
            {
                {0,1,0},
                {1,1,1},
                {0,1,0}
            });
            BinaryImage resultImage = new BinaryImage(testImage.Boundary(structureElement));
            testImage.Print();
            resultImage.Print();
        }
    }
    class BinaryImage
    {
        public int Width { get; private set; }
        public int Height { get; private set; }
        public int[,] Pixels { get; set; }
        public BinaryImage(int width, int height)
        {
            Pixels = new int[width, height];
            Width = width;
            Height = height;
        }
        public BinaryImage(int[,] pixels)
        {
            Pixels = pixels;
            Width = pixels.GetLength(0);
            Height = pixels.GetLength(1);
        }
        public void Flush(int value)
        {
            for (int i =0;i < Width; i++)
                for (int j =0; j <Height; j++)
                    Pixels[i, j] = value;
        }
        public int[,] Erode(BinaryImage structureElement)
        {
            int[,] ret = new int[Width, Height];
            for (int i = 0; i < Width; i++)
            {
                for (int j = 0; j < Height; j++)
                {
                    bool keep = true;
                    for (int m = 0; m < structureElement.Width; m++)
                    {
                        for (int n = 0; n < structureElement.Height; n++)
                        {
                            if (Pixels[Math.Min(Math.Max(0, i + m - structureElement.Width / 2), Width-1),
                                        Math.Min(Math.Max(0, j + n - structureElement.Height / 2), Height-1)] != structureElement.Pixels[m, n] && structureElement.Pixels[m,n]!=0)
                            {
                                keep = false;
                                break;
                            }
                        }
                    }
                    ret[i, j] = keep ? 1 : 0;
                }
            }
            return ret;
        }
        public int[,] Boundary(BinaryImage structureElement)
        {
            int[,] result = Erode(structureElement);
            return Minus(new BinaryImage(result));
        }
        public int[,] Minus(BinaryImage image)
        {
            int[,] minus = new int[Width, Height];
            for (int i = 0; i < Width; i++)
                for (int j = 0; j < Height; j++)
                    minus[i, j] = Pixels[i, j] - image.Pixels[i, j];
            return minus;
        }
       
        public void Print()
        {
            for (int i = 0; i < Width; i++)
            {
                for (int j = 0; j < Height; j++)
                {
                    if (Pixels[i, j] != 0)
                        Console.Write("*");
                    else
                        Console.Write(" ");
                }
                Console.WriteLine();
            }
        }
    }
}

The result of above program looks like the following picture. Left is the original image, right is detected border.

The quality of the border is affected by the structure element. The 3 by 3 structure element used by the code can get fairly good result in many cases. The following example is the border detected using a 4 by 4 structure element:

  BinaryImage structureElement = new BinaryImage(new int[4, 4]
            {
                {0,1,1,0},
                {1,1,1,1},
                {1,1,1,1},
                {0,1,1,0}
            });

Reference_Draw Hyperbolic Geometry

Draw hyperbolic geometry (Poincare disc model) using C# WPF / GDI+

This post is about how to draw a hyperbolic plane using WPF and GDI+ in C#. The program uses Poincare disc (conformal model) representation. I’ll present a little math background and then the code to draw these two patterns:

Background

2D Hyperbolic geometry use an Euclidean circle to represent a hyperbolic plane. One interesting feature of this presentation is that although the hyperbolic plane is infinite, it’s represented in a finite Euclidean circle. You may have seen some tree structure presentations utilizing such features so user can navigate through a massive tree structure within a circular viewport. Poincare disc model uses arcs of circles, which meet the bounding circle orthogonally, to represent straight lines – yes the curves you see in above patterns are actually straight lines in hyperbolic plane. To illustrate above pictures, the core task is to find such circles. Specifically, the arc decided by an angle q1 and an angle q2 is provided by a circle shown in the following diagram. To find the circle we need to know both the center of the circle (cx, cy) and radius of the circle (r).

Now the question turns into several simple calculations in our familiar Euclidean plane. I won’t repeat the details here because the calculation will be in the source code. Now let’s start coding.
Drawing the pattern to the left using WPF
Create a new WPF application and update MainWindow.xaml:

<Window x:Class="HyperbolicGeometry.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        Title="MainWindow" Height="500" Width="525" Loaded="Window_Loaded">
    <StackPanel HorizontalAlignment="Stretch" VerticalAlignment="Stretch" Margin="10,10,10,10">
        <Image Name="space" Height="350" Width="350" HorizontalAlignment="Center" />
    </StackPanel>
</Window>

Nothing exciting here, we just created a image holder that we’ll draw on on the Window_Load event handler. The code-behind looks like this:

private static bool flip = true;
private void Window_Loaded(object sender, RoutedEventArgs e)
{
    int R = (int)(space.Width / 2);
    DrawingVisual visual = new DrawingVisual();
    using (DrawingContext context = visual.RenderOpen())
     {
        visual.Clip = new EllipseGeometry(new Point(R, R), R, R);
        context.DrawEllipse(Brushes.Transparent, new Pen(Brushes.Black,1), new Point(R,R),R,R);
        for (double s = 1.5; s <= 48; s = s * 2)
        {
            for (double q1 = 0; q1 < Math.PI * 2; q1 += Math.PI / s)
            {
                double q2 = q1 + Math.PI / s;
                drawLine(R, context, q1, q2);
            }
            flip = !flip;
        }
    }
    RenderTargetBitmap bmp = new RenderTargetBitmap(R*2,R*2, 96, 96, PixelFormats.Pbgra32);
    bmp.Render(visual);
    space.Source = bmp;
}
private static void drawLine(int R, DrawingContext context, double q1, double q2)
{
    double f = (q2 - q1) / 2;
    double dq = Math.Abs(f);
    double r = R * Math.Tan(dq);
    double rp = Math.Sqrt(r * r + R * R);
    double cx = R + rp * Math.Cos(q1 + f);
    double cy = R + rp * Math.Sin(q1 + f);
    double beta = Math.PI - dq * 2;
    double k = Math.PI / 2 + q2;
    context.DrawEllipse(flip?Brushes.Blue: Brushes.White, new Pen(Brushes.Blue, 1), new Point(cx, cy), r, r);
}

The algorithm of finding the circles is all in the drawLine method. Window_Load generates a series of angle pairs and feed them to drawLine method to draw lines at different locations. Note the visual.Clip call – we are constraining our painting within the bounding circle. This is necessary because when we call DrawEllipse method in drawLine parts of the ellipses will fall outside of the bounding circle. At last, the flip flag allows to vary between blue color and white color to create the final pattern.
Drawing the pattern to the right using WPF
It’s indeed very easy to modify above implementation to draw the pattern to the right but we’ll go one step further here: Instead of relying on the clipping region, we’ll figure out proper arc starting angle and ending angles so that we only paint the necessary portions:

Comment out the viual.Clip call – this makes sure we don’t cheat.
Add a drawArc method. WPF doesn’t have a native draw arc method so we have to fashion one:

private static void drawArc(DrawingContext drawingContext, Brush brush,
            Pen pen, Point start, Point end, Size radius)
{
    PathGeometry geometry = new PathGeometry();
    PathFigure figure = new PathFigure();
    geometry.Figures.Add(figure);
    figure.StartPoint = start;
    figure.Segments.Add(new ArcSegment(end, radius, 0, false, SweepDirection.Clockwise, true));
    drawingContext.DrawGeometry(brush, pen, geometry);
}

Modify drawLine method. Instead of calling DrawEllipse, call drawArc:

drawArc(context, Brushes.Transparent, new Pen(Brushes.Blue, 2),
        new Point(cx + r * Math.Cos(k), cy + r * Math.Sin(k)),
        new Point(cx + r * Math.Cos(k + beta), cy + r * Math.Sin(k + beta)), new Size(r,r));

Compile and run the code again, you should see the pattern to the right.

Drawing the patterns using GDI+
Above code can be easily modified to work with GDI+. The only place that needs special attention is the way how GDI+ draws arcs is different from how WPF does it. Under GDI+, the drawArc call should look like:

gc.DrawArc(new Pen(new SolidBrush(Color.Blue),2),
    new Rectangle((int)(cx-r), (int)(cy-r), (int)(r*2), (int)(r*2)),
    (float)(k * 180 / Math.PI), 
    (float)(beta * 180) / Math.PI));

(where gc is the Graphics instance you use to paint on). Consult .Net documentations for differences how arcs are represented differently in the two frameworks.

Wednesday, August 27, 2014

Detecting Simple Shapes C#

Detecting some simple shapes in images by Andrew Kirillov The article describes how to detect/check some simple shapes like circles, quadrilaterals, triangles, rectangles, etc.	Posted: June 30, 2010 Programming languages: C# AForge.NET framework: 2.1.3

Sample application (sources) - 215K
Sample application (biaries) - 206K

Introduction

Doing image processing and especially blob analysis it is often required to check some objects' shape and depending on it perform further processing of a particular object or not. For example, some applications may require finding only circles from all the detected objects, or quadrilaterals, rectangles, etc. The article will describe some basic techniques, which may allow detecting such shapes like circles, triangles and quadrilaterals (plus their subtypes, like rectangle, rhombus, rectangled triangle, etc.). Note: most of the code snippets provided in the article are based on the AForge.NET framework, but can be easily implemented using any other library providing similar image processing routines.

Prerequisites

Before we go into shape analysis, we need to locate objects we are going to check and their edge pixels, since they will be used as the input for shape analysis algorithms we are going to discuss. First thing to note is that we'll suppose our input images contain objects of different colors on a black background. If the background is not black, then it should be turned into black using, for example, some color filters, before we start with described below steps. Here is a sample source image we are going to process and detect shapes of objects on it:

The first step to do is to locate each separate object in the input image, which is done using BlobCounter:

// create instance of blob counter
BlobCounter blobCounter = new BlobCounter( );
// process input image
blobCounter.ProcessImage( bitmap );
// get information about detected objects
Blob[] blobs = blobCounter.GetObjectsInformation( );

The array with information about blobs contains such values for each detected blob like its rectangle, centre of gravity, mean color, standard deviation of color, etc. Some of these values can be use for filtering of blobs. For example, user may not need to do further processing of blobs which have too small width or height.
For analyzing of blobs' shapes we'll need to get their edge pixels. For this purpose the blob counting class provides 3 methods, which allow to get left/right edge pixels, top/bottom edge pixels and all edge pixels. The methods are GetBlobsLeftAndRightEdges, GetBlobsTopAndBottomEdges and GetBlobsEdgePoints. The image below demonstrates edge pixel detected by each of these methods).

Now we can start with developing a shape detection algorithm, which should detect type of a shape for a given set of shapes' edge pixels. For these purpose we'll use the GetBlobsEdgePoints method, which provides all edge points, so the detection could be more accurate.

Circle detection

The idea of circle detection algorithm is based on the assumption, that all circle's edge points have the same distance to its center, which equals to circle's radius. Of course doing different image processing tasks it may happen that circles may have some distortion of their shape, so some edge pixels may be closer or further to circle's center. But this variation in distance to center should not be very big and should reside in certain limits. If it is too big, than most probably the object has other than circle shape.
For a given set of circle's edge pixels it is quite easy to estimate its center and radius by finding bounding rectangle for the specified set of points (GetBoundingRectangle() method can be used, for example):

// get bounding rectangle of the points list
IntPoint minXY, maxXY;
PointsCloud.GetBoundingRectangle( edgePoints, out minXY, out maxXY );
// get cloud's size
IntPoint cloudSize = maxXY - minXY;
// calculate center point
DoublePoint center = minXY + (DoublePoint) cloudSize / 2;
// calculate radius
float radius = ( (float) cloudSize.X + cloudSize.Y ) / 4;

When we have estimated circle's radius and center, all we need to do is to go through the list of the shape's edge pixels, calculate their distance to the estimated center and check the difference with estimated radius - distance between provided edge points and estimated circle. Instead of checking each individual distance value for each edge pixel, we may just calculate mean distance:

// calculate mean distance between provided edge points
// and estimated circle's edge
float meanDistance = 0;

for ( int i = 0, n = edgePoints.Count; i < n; i++ )
{
    meanDistance += Math.Abs(
        (float) center.DistanceTo( edgePoints[i] ) - radius );
}
meanDistance /= edgePoints.Count;

Now, when we have calculated mean distance between provided shape's edge points and estimated circle, we just need to check if the value falls into certain range or not. If it is too big, then it means that the specified shape is not a circle, since its edge points are quite away on the average from the estimated circle. Ideally the value should be close to 0 meaning that all specified edge points fit very well the estimated circle. But, as we already mentioned, we may allow some distortion of the shape.
Which level of distortion do we want to allow for circle shapes? Definitely we can not set any fixed value. The distortion limit should depend on shape's size, so it would allow higher level of distortion for bigger shapes and lower value of distortion for smaller shapes. For example, distortion level may be calculated this way:

float maxDitance = ( (float) cloudSize.X + cloudSize.Y ) / 2 *
    relativeDistortionLimit;

The relativeDistortionLimit value can be made configurable, so user could affect acceptable distortion level and precision of circles' recognition. Experimentally I found that 0.03 value (3%) works quite well and allows to filter circles from other shapes.
The approach with relative distortion limit works well, but there is one issue, which we may need to resolve. In the case if we deal with small circles, like 10x10 pixels in size, the calculated distortion limit will be equal to 0.3. If a circle has some little distortion, then it may not be recognized as circle. For example, doing some image processing we may get circles which are 9x10 or 11x10 in size, which potentially may lead to higher calculated distortion than the specified limit. To avoid this we may add additional parameter, which is minimum acceptable distortion. So the resulting formula looks like this:

// some parameters, which also are exposed as properties
private float minAcceptableDistortion = 0.5f;
private float relativeDistortionLimit = 0.03f;

...

float maxDitance = Math.Max( minAcceptableDistortion,
    ( (float) cloudSize.X + cloudSize.Y ) / 2 * relativeDistortionLimit );

All we need to do now is to compare previously calculated meanDistance with the maxDitance. If it is less or equal to maxDistance, then we have a circle, otherwise we have too big distortion, which means the shape is not a circle most probably.
All the above code is also available as a single call in the AForge.NET framework: IsCircle(). Finally, here is a small sample showing how to use/test this method and some results:

// locate objects using blob counter
BlobCounter blobCounter = new BlobCounter( );
blobCounter.ProcessImage( bitmap );
Blob[] blobs = blobCounter.GetObjectsInformation( );
// create Graphics object to draw on the image and a pen
Graphics g = Graphics.FromImage( bitmap );
Pen redPen = new Pen( Color.Red, 2 );
// check each object and draw circle around objects, which
// are recognized as circles
for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );

    Point center;
    float radius;

    if ( shapeChecker.IsCircle( edgePoints, out center, out radius ) )
    {
        g.DrawEllipse( redPen,
            (int) ( center.X - radius ),
            (int) ( center.Y - radius ),
            (int) ( radius * 2 ),
            (int) ( radius * 2 ) );
    }
}

redPen.Dispose( );
g.Dispose( );

The first two sample images above demonstrate detection of more or less ideal circles. We can see there that ellipses and other shapes are not misclassified as circles. On the third image we have two circles, which have some distortion in their shape, but still are recognized successfully as circles.

Detection of quadrilaterals

Detection of quadrilaterals and triangles has pretty much the same idea - we are checking mean distance between provided shape's edge pixels and the edge of estimated quadrilateral/triangle. The only difference here is the way how to estimate parameters of the shape we want to recognize and how to check distance to the estimated shape.
Let's start with quadrilaterals first. For a given shape we can make an assumption that it is a quadrilateral, find its four corners (using FindQuadrilateralCorners, for example) and then using similar method as we've used for circles we can check if our assumption is correct or not - check how good the shape fits into the quadrilateral with assumed parameters. The below code and picture demonstrate the idea of finding quadrilateral points for different types of shapes:

// locate objects using blob counter
BlobCounter blobCounter = new BlobCounter( );
blobCounter.ProcessImage( bitmap );
Blob[] blobs = blobCounter.GetObjectsInformation( );
// create Graphics object to draw on the image and a pen
Graphics g = Graphics.FromImage( bitmap );
Pen bluePen = new Pen( Color.Blue, 2 );
// check each object and draw circle around objects, which
// are recognized as circles
for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );
    List<IntPoint> corners = PointsCloud.FindQuadrilateralCorners( edgePoints );
    
    g.DrawPolygon( bluePen, ToPointsArray( corners ) );    
}

bluePen.Dispose( );
g.Dispose( );

As we can see on the image above, we may have different objects and quadrilateral finder provides different results for them. For shapes which really look like quadrilateral, the quadrilateral finder is able to find their corners more or less correctly (problem may occur in the case if object has rounded corners). But for other types of objects (circles, starts, etc.), quadrilateral finder does not provide any good result. And this is correct, since this routine supposes the given shape is really quadrilateral.
Now, when we made the assumption that a subjected object is quadrilateral and we got its corner, we just need to see how good is the fit of the object into the quadrilateral with those found corner - we need to make sure there are no edge points of the given shape which are too far away from the edge of the assumed quadrilateral. And here we'll apply the same algorithm as we used for circle checking. The only difference is the way of calculating distance between point and the closest quadrilateral' edge ...

int cornersCount = corners.Count; // should be 4 for quadrilateral

// 1 - prepare for calculations of distance between a point
//     and the quadrilateral

// lines coefficients (for representation as y(x)=k*x+b)
double[] k = new double[cornersCount];
double[] b = new double[cornersCount];
double[] div = new double[cornersCount]; // precalculated divisor
bool[] isVert = new bool[cornersCount];

for ( int i = 0; i < cornersCount; i++ )
{
    IntPoint currentPoint = corners[i];
    IntPoint nextPoint = ( i + 1 == cornersCount ) ?
        corners[0] : corners[i + 1];

    if ( !( isVert[i] = nextPoint.X == currentPoint.X ) )
    {
        k[i] = (double) ( nextPoint.Y - currentPoint.Y ) /
                        ( nextPoint.X - currentPoint.X );
        b[i] = currentPoint.Y - k[i] * currentPoint.X;
        div[i] = Math.Sqrt( k[i] * k[i] + 1 );
    }
}

// 2 - calculate distances between edge points and quadrilateral sides
float meanDistance = 0;

for ( int i = 0, n = edgePoints.Count; i < n; i++ )
{
    float minDistance = float.MaxValue;

    for ( int j = 0; j < cornersCount; j++ )
    {
        float distance = 0;

        if ( !isVert[j] )
        {
            distance = (float) Math.Abs(
                ( k[j] * edgePoints[i].X + b[j] - edgePoints[i].Y )
                    / div[j] );
        }
        else
        {
            distance = Math.Abs( edgePoints[i].X - corners[j].X );
        }

        // finding distance to the closest side
        if ( distance < minDistance )
            minDistance = distance;
    }

    meanDistance += minDistance;
}
meanDistance /= edgePoints.Count;

The above code does the same as we did with circles - it goes through the list of all provided shape's edge points and finds distance between them and the assumed shape. Since we deal here with quadrilateral, the distance between point and the assumed quadrilateral is calculated as a distance to the closest quadrilateral's side. In the end of the above calculations we receive the same as before - mean distance between given edge points and the edge of the assumed quadrilateral. If the mean distance is not greater than a certain limit, then the shape seems to be a quadrilateral; otherwise it is not. This check is done in the similar way like before:

// get bounding rectangle of the corners list
IntPoint minXY, maxXY;
PointsCloud.GetBoundingRectangle( corners, out minXY, out maxXY );
IntPoint rectSize = maxXY - minXY;

float maxDitance = Math.Max( minAcceptableDistortion,
    ( (float) rectSize.X + rectSize.Y ) / 2 * relativeDistortionLimit );

return ( meanDistance <= maxDitance ); // true for quadrilateral

Same as with circles detection, the code for detection quadrilaterals was put into the AForge.NET framework: IsQuadrilateral(). And here is a small sample to test it:

// check each object and draw polygon around objects, which
// are recognized as quadrilaterals
for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );
    List<IntPoint> corners;

    if ( shapeChecker.IsQuadrilateral( edgePoints, out corners ) )
    {
        g.DrawPolygon( redPen, ToPointsArray( corners ) );
    }
}

As we can see on the picture above, only quadrilaterals were highlighted with a border. Other shapes, like circles, triangles, stars, etc., were filtered out, which means that although quadrilateral finder gave us some points as 4 corners, the further check for fitting assumed quadrilateral did not pass. Which is nice.

Detection of triangles

Idea of detecting triangles is almost the same like the idea of detecting quadrilaterals - just need to find 3 corners of the supposed triangle and then check the fitting of specified edge points into the assumed triangle. More of it, it is even based on the same API - FindQuadrilateralCorners. Why is it so? Because if we take a closer look into documentation of FindQuadrilateralCorners, we may find that the function returns only 3 corners, if it cannot find the 4^th one, which is the case with triangles. For example, using the code below we can make sure it really works - we'll highlight quadrilaterals and triangles with different colors:

for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );
    List<IntPoint> corners = PointsCloud.FindQuadrilateralCorners( edgePoints );

    g.DrawPolygon( ( corners.Count == 4 ) ? redPen : bluePen,
        ToPointsArray( corners ) );
}

Usage of quadrilateral finder - checking number of detected corners

Since the previously written code for checking fitting of provided edge points into assumed quadrilateral does not depend on number of corners (that code actually works not just for quadrilaterals, but for any convex polygon), we can use it for triangles as well. So nothing more needs to be written - check of triangles is the same as check for quadrilaterals; we just distinguish these shapes by number of corners found by quadrilateral finder.
Here is a test code with new framework's function - IsTriangle():

// check each object and draw triangle around objects, which
// are recognized as triangles
for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );
    List<IntPoint> corners;

    if ( shapeChecker.IsTriangle( edgePoints, out corners ) )
    {
        g.DrawPolygon( redPen, ToPointsArray( corners ) );
    }
}

Rectangles, squares, equilateral triangles, etc.

Once we made a check for quadrilateral shape and got its 4 corners, we may do some further checks to detect sub-type of the quadrilateral: trapezoid, parallelogram, rectangle, rhombus or square. These checks are based on checking angles between opposite/adjacent sides and their length. First we check two opposite sides - if they are parallel, then we get at least trapezoid. Then we check two other sides - if they are parallel two, then we have parallelogram. If two adjacent sides are perpendicular, then we got rectangle. The last check is to compare length of two adjacent sides - if they are equal, then parallelogram becomes rhombus and rectangle becomes square. And of course we need to add small possible error, so if angle between lines equals to 2 degrees, they are still treated as parallel. The below code shows implementation of the idea:

// get angles between 2 pairs of opposite sides
float angleBetween1stPair = Tools.GetAngleBetweenLines(
    corners[0], corners[1], corners[2], corners[3] );
float angleBetween2ndPair = Tools.GetAngleBetweenLines(
    corners[1], corners[2], corners[3], corners[0] );

// check 1st pair for parallelism
if ( angleBetween1stPair <= angleError )
{
    subType = PolygonSubType.Trapezoid;

    // check 2nd pair for parallelism
    if ( angleBetween2ndPair <= angleError )
    {
        subType = PolygonSubType.Parallelogram;

        // check angle between adjacent sides
        if ( Math.Abs( Tools.GetAngleBetweenVectors(
            corners[1], corners[0], corners[2] ) - 90 ) <= angleError )
        {
            subType = PolygonSubType.Rectangle;
        }

        // get length of 2 adjacent sides
        float side1Length = (float) corners[0].DistanceTo( corners[1] );
        float side2Length = (float) corners[0].DistanceTo( corners[3] );

        if ( Math.Abs( side1Length - side2Length ) <= maxLengthDiff )
        {
            subType = ( subType == PolygonSubType.Parallelogram ) ?
                PolygonSubType.Rhombus : PolygonSubType.Square;
        }
    }
}
else
{
    // check 2nd pair for parallelism - last chence to detect trapezoid
    if ( angleBetween2ndPair <= angleError )
    {
        subType = PolygonSubType.Trapezoid;
    }
}

Checking triangle subtype is even simpler. If it's all angles are near to 60 degrees, then we have equilateral triangle. If not then we check if some of two angles are equal - isosceles triangle. Also we may need to check if one of the angles equals to 90 degrees - this will give us rectangled triangle.

// get angles of the triangle
float angle1 = Tools.GetAngleBetweenVectors( corners[0], corners[1], corners[2] );
float angle2 = Tools.GetAngleBetweenVectors( corners[1], corners[2], corners[0] );
float angle3 = Tools.GetAngleBetweenVectors( corners[2], corners[0], corners[1] );

// check for equilateral triangle
if ( ( Math.Abs( angle1 - 60 ) <= angleError ) &&
     ( Math.Abs( angle2 - 60 ) <= angleError ) &&
     ( Math.Abs( angle3 - 60 ) <= angleError ) )
{
    subType = PolygonSubType.EquilateralTriangle;
}
else
{
    // check for isosceles triangle
    if ( ( Math.Abs( angle1 - angle2 ) <= angleError ) ||
         ( Math.Abs( angle2 - angle3 ) <= angleError ) ||
         ( Math.Abs( angle3 - angle1 ) <= angleError ) )
    {
        subType = PolygonSubType.IsoscelesTriangle;
    }

    // check for rectangled triangle
    if ( ( Math.Abs( angle1 - 90 ) <= angleError ) ||
         ( Math.Abs( angle2 - 90 ) <= angleError ) ||
         ( Math.Abs( angle3 - 90 ) <= angleError ) )
    {
        subType = ( subType == PolygonSubType.IsoscelesTriangle ) ?
            PolygonSubType.RectangledIsoscelesTriangle :
            PolygonSubType.RectangledTriangle;
    }
}

All the above code is available in AForge.NET Framework as CheckPolygonSubType(). Here is a quick test of it and some pictures with result:

// dictionary of color to highlight different shapes
Dictionary<PolygonSubType, Color> colors =
    new Dictionary<PolygonSubType, Color>( );

colors.Add( PolygonSubType.Unknown, Color.White );
colors.Add( PolygonSubType.Trapezoid, Color.Orange );
colors.Add( PolygonSubType.Parallelogram, Color.Red );
colors.Add( PolygonSubType.Rectangle, Color.Green );
colors.Add( PolygonSubType.Square, Color.Blue );
colors.Add( PolygonSubType.Rhombus, Color.Gray );

colors.Add( PolygonSubType.EquilateralTriangle, Color.Pink );
colors.Add( PolygonSubType.IsoscelesTriangle, Color.Purple );
colors.Add( PolygonSubType.RectangledTriangle, Color.SkyBlue );
colors.Add( PolygonSubType.RectangledIsoscelesTriangle, Color.SeaGreen );

for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints = blobCounter.GetBlobsEdgePoints( blobs[i] );
    List<IntPoint> corners;

    if ( shapeChecker.IsConvexPolygon( edgePoints, out corners ) )
    {
        // check subtype
        PolygonSubType subType = shapeChecker.CheckPolygonSubType( corners );

        using ( Pen pen = new Pen( colors[subType], 2 ) )
        {
            g.DrawPolygon( pen, ToPointsArray( corners ) );
        }
    }
}

Taking a look at the above pictures and colors of highlighting borders, we may see that all triangles' and quadrilaterals' sub-types were recognized quite successfully. You may ask - Why there is one triangle with white border, which corresponds to unknown sub-type? Because it is just a triangle - not equilateral, not isosceles and not rectangled - sub-type is unknown.

Test on a real image

Now, let's try all the above ideas on some real world image. Consider the image below. Our task is to find all coins (circles) and see what else we can detect.

As a first step we need to remove background - make it completely black, so we could easily locate stand alone objects. For this purpose we may try using ColorFiltering image processing routine:

// lock image
BitmapData bitmapData = bitmap.LockBits(
    new Rectangle( 0, 0, bitmap.Width, bitmap.Height ),
    ImageLockMode.ReadWrite, bitmap.PixelFormat );

// step 1 - turn background to black
ColorFiltering colorFilter = new ColorFiltering( );

colorFilter.Red   = new IntRange( 0, 64 );
colorFilter.Green = new IntRange( 0, 64 );
colorFilter.Blue  = new IntRange( 0, 64 );
colorFilter.FillOutsideRange = false;

colorFilter.ApplyInPlace( bitmapData );

Now we'll do detection of blobs. Although we have some noisy pixels, which were left from bright areas of background, we are not going to use any noise suppression image processing routines. We'll just instruct blob counter to filter objects by size, so all tiny objects are ignored:

// step 2 - locating objects
BlobCounter blobCounter = new BlobCounter( );

blobCounter.FilterBlobs = true;
blobCounter.MinHeight = 5;
blobCounter.MinWidth = 5;

blobCounter.ProcessImage( bitmapData );
Blob[] blobs = blobCounter.GetObjectsInformation( );
bitmap.UnlockBits( bitmapData );

And finally we'll check type of each object and highlight them with different color:

// step 3 - check objects' type and highlight
SimpleShapeChecker shapeChecker = new SimpleShapeChecker( );

Graphics g = Graphics.FromImage( bitmap );
Pen redPen = new Pen( Color.Red, 2 );
Pen yellowPen = new Pen( Color.Yellow, 2 );
Pen greenPen = new Pen( Color.Green, 2 );
Pen bluePen = new Pen( Color.Blue, 2 );

for ( int i = 0, n = blobs.Length; i < n; i++ )
{
    List<IntPoint> edgePoints =
        blobCounter.GetBlobsEdgePoints( blobs[i] );

    Point center;
    double radius;

    if ( shapeChecker.IsCircle( edgePoints, out center, out radius ) )
    {
        g.DrawEllipse( yellowPen,
            (float) ( center.X - radius ), (float) ( center.Y - radius ),
            (float) ( radius * 2 ), (float) ( radius * 2 ) );
    }
    else
    {
        List<IntPoint> corners;

        if ( shapeChecker.IsQuadrilateral( edgePoints, out corners ) )
        {
            if ( shapeChecker.CheckPolygonSubType( corners ) ==
                PolygonSubType.Rectangle )
            {
                g.DrawPolygon( greenPen, ToPointsArray( corners ) );
            }
            else
            {
                g.DrawPolygon( bluePen, ToPointsArray( corners ) );
            }
        }
        else
        {
            corners = PointsCloud.FindQuadrilateralCorners( edgePoints );
            g.DrawPolygon( redPen, ToPointsArray( corners ) );
        }
    }
}

redPen.Dispose( );
greenPen.Dispose( );
bluePen.Dispose( );
yellowPen.Dispose( );
g.Dispose( );

So, what did we get? As we can see from the above picture, all coins (circles) were detected successfully and highlighted with yellow borders. Also we got one rectangle in green border, one quadrilateral (it is actually trapezoid, we just did not make appropriate check) in blue border and 2 unknown shapes in red borders. Everything seems to be fine except of these two last shapes in red. Why they were not recognized? With gray box it is more or less obvious - since it has rounded corners, the quadrilateral finding routine was a bit confused and picked up those corners, which seem to be best from its algorithm's point of view. After that shape fitting has failed, since some edge points of the box are quite away from the supposed quadrilateral. With Panasonic battery the failure reason is caused by the "+" connector, which was chosen as wrong corner.
Although we got two unrecognized shapes, it does not mean that all of the above ideas about shape checking are absolutely bad and cannot be used. These two failures are actually caused not by failure of the described shape fitting algorithm, but by quadrilateral finder, which has failed to find proper corners. And it has a reason for this. So if we would use another algorithm for finding corners of the supposed quadrilateral, then things could work better for the last two objects.

Conclusion

As it was already mentioned in the article, all of the code became available as part of the AForge.NET framework. Except of the few mentioned functions, there is a SimpleShapeChecker class, which encapsulates them all, plus has some additional methods. For example, we already have seen, that there is a method, which checks if the specified set of edge points form a quadrilateral and provides its corners on success, which are detected by quadrilateral finder implemented in the framework. However, there is also another method, which allows user to specify his own quadrilateral corners. These gives user a bit more flexibility, since he may use another algorithm for finding corners of quadrilateral and then just use framework's method for checking fitting of edge points into the assumed quadrilateral.
In addition the framework provides a sample application, which demonstrates usage of the implemented functions and allow to do some quick experiments. The sample application is available as part of framework's distribution and also can be download using links in the top of the article.