/// Process image looking for corners.
/// </summary>
///
/// <param name="image">Source image data to process.</param>
///
/// <returns>Returns list of found corners (X-Y coordinates).</returns>
///
/// <exception cref="UnsupportedImageFormatException">
/// The source image has incorrect pixel format.
/// </exception>
///
public unsafe List<IntPoint> ProcessImage(UnmanagedImage image)
{
// check image format
if (
(image.PixelFormat != PixelFormat.Format8bppIndexed) &&
(image.PixelFormat != PixelFormat.Format24bppRgb) &&
(image.PixelFormat != PixelFormat.Format32bppRgb) &&
(image.PixelFormat != PixelFormat.Format32bppArgb)
)
{
throw new UnsupportedImageFormatException("Unsupported pixel format of the source image.");
}
// make sure we have grayscale image
UnmanagedImage grayImage = null;
if (image.PixelFormat == PixelFormat.Format8bppIndexed)
{
grayImage = image;
}
else
{
// create temporary grayscale image
grayImage = Grayscale.CommonAlgorithms.BT709.Apply(image);
}
// get source image size
int width = grayImage.Width;
int height = grayImage.Height;
int srcStride = grayImage.Stride;
int srcOffset = srcStride - width;
// 1. Calculate partial differences
float[,] diffx = new float[height, width];
float[,] diffy = new float[height, width];
float[,] diffxy = new float[height, width];
fixed (float* pdx = diffx, pdy = diffy, pdxy = diffxy)
{
byte* src = (byte*)grayImage.ImageData.ToPointer() + srcStride + 1;
// Skip first row and first column
float* dx = pdx + width + 1;
float* dy = pdy + width + 1;
float* dxy = pdxy + width + 1;
// for each line
for (int y = 1; y < height - 1; y++)
{
// for each pixel
for (int x = 1; x < width - 1; x++, src++, dx++, dy++, dxy++)
{
// Convolution with horizontal differentiation kernel mask
float h = ((src[-srcStride + 1] + src[+1] + src[srcStride + 1]) -
(src[-srcStride - 1] + src[-1] + src[srcStride - 1])) * 0.166666667f;
// Convolution vertical differentiation kernel mask
float v = ((src[+srcStride - 1] + src[+srcStride] + src[+srcStride + 1]) -
(src[-srcStride - 1] + src[-srcStride] + src[-srcStride + 1])) * 0.166666667f;
// Store squared differences directly
*dx = h * h;
*dy = v * v;
*dxy = h * v;
}
// Skip last column
dx++; dy++; dxy++;
src += srcOffset + 1;
}
// Free some resources which wont be needed anymore
if (image.PixelFormat != PixelFormat.Format8bppIndexed)
grayImage.Dispose();
}
// 2. Smooth the diff images
if (sigma > 0.0)
{
float[,] temp = new float[height, width];
// Convolve with Gaussian kernel
convolve(diffx, temp, kernel);
convolve(diffy, temp, kernel);
convolve(diffxy, temp, kernel);
}
// 3. Compute Harris Corner Response Map
float[,] map = new float[height, width];
fixed (float* pdx = diffx, pdy = diffy, pdxy = diffxy, pmap = map)
{
float* dx = pdx;
float* dy = pdy;
float* dxy = pdxy;
float* H = pmap;
float M, A, B, C;
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++, dx++, dy++, dxy++, H++)
{
A = *dx;
B = *dy;
C = *dxy;
if (measure == HarrisCornerMeasure.Harris)
{
// Original Harris corner measure
M = (A * B - C * C) - (k * ((A + B) * (A + B)));
}
else
{
// Harris-Noble corner measure
M = (A * B - C * C) / (A + B + Accord.Math.Special.SingleEpsilon);
}
if (M > threshold)
{
*H = M; // insert value in the map
}
}
}
}
// 4. Suppress non-maximum points
List<IntPoint> cornersList = new List<IntPoint>();
// for each row
for (int y = r, maxY = height - r; y < maxY; y++)
{
// for each pixel
for (int x = r, maxX = width - r; x < maxX; x++)
{
float currentValue = map[y, x];
// for each windows' row
for (int i = -r; (currentValue != 0) && (i <= r); i++)
{
// for each windows' pixel
for (int j = -r; j <= r; j++)
{
if (map[y + i, x + j] > currentValue)
{
currentValue = 0;
break;
}
}
}
// check if this point is really interesting
if (currentValue != 0)
{
cornersList.Add(new IntPoint(x, y));
}
}
}
return cornersList;
}
