Optional standard gradient calculation in Sobel filter

[GoogleCodeExporter]

Instead

        Math.Sqrt( gx * gx + gy * gy )

, framework evaluates

        Math.Min( 255, Math.Abs(gx) + Math.Abs(gy) )

This produces differences in two ways

1.- Simpliffied expression is always a bit bigger. Differences are noticiable 
but not important.

gx + gy    <=>    Sqrt( (gx + gy)^2 )    <=>    Sqrt( gx^2 + 2*gx*gy + gy^2 )


2.- Evaluating "Min" before normalization, produces important differences in 
result.

We can see differences in this image that combines 4 results. Results has been 
inverted and thresholded (235) to clearly show differences  

* top-left: Min + Sum (current AForge implementation)
* top-right: Min + Sqrt
* bottom-left: Sum
* bottom-right: Sqrt

https://dl.dropboxusercontent.com/u/49778953/test_sobel_results.JPG

Differences between results clearly show that smoothing effect is lost when 
using "Min"

This is my test image
https://dl.dropboxusercontent.com/u/49778953/test_sobel.jpg


I suggest to introduce a new Property "StandardGradientCalculation" with 
default value = False. This way, behaviour will not change in production 
systems.


This is the function i've used. Note the use of 

float[,] gradients = new float[source.Height, source.Width];

, to avoid byte overflow


public unsafe void ProcessFilter(AForge.Imaging.UnmanagedImage source, 
AForge.Imaging.UnmanagedImage destination, Rectangle rect,
                                 bool takeMin, bool useSum)
{
    int[,] kernelX = {{-1, +0, +1},
                      {-2, +0, +2},
                      {-1, +0, +1}};

    int[,] kernelY = {{-1, -2, -1},
                      {+0, +0, +0},
                      {+1, +2, +1}};

    // 3x3 kernel => kernelX.GetUpperBound(0) = 2
    int kernelSize = kernelX.GetUpperBound(0) + 1;
    if (kernelSize % 2 == 0)
    {
        throw new ArgumentException("Size of kernel matrix must have odd");
    }
    int kernelRadius = (kernelSize - 1) / 2;
    byte srcBytesPerPixel = (byte)(System.Drawing.Bitmap.GetPixelFormatSize(source.PixelFormat) / 8);
    byte dstBytesPerPixel = (byte)(System.Drawing.Bitmap.GetPixelFormatSize(destination.PixelFormat) / 8);

    // processing start and stop X,Y positions                
    int startX = rect.Left + kernelRadius;
    int startY = rect.Top + kernelRadius;
    int stopX = startX + rect.Width - 2 * kernelRadius;
    int stopY = startY + rect.Height - 2 * kernelRadius;

    int dstStride = destination.Stride;
    int srcStride = source.Stride;

    byte* srcCurrent = (byte*)source.ImageData.ToPointer() + srcStride * startY + srcBytesPerPixel * startX;
    byte* dstCurrent = (byte*)destination.ImageData.ToPointer() + dstStride * startY + dstBytesPerPixel * startX;

    byte* src = null;
    byte* dst = null;

    // variables for gradient calculation
    float g = 0;
    float max = 0;
    float[,] gradients = new float[source.Height, source.Width];

    // for each line
    for (int y = startY; y < stopY; y++)
    {
        // for each pixel
        src = srcCurrent;
        dst = dstCurrent;
        for (int x = startX; x < stopX; x++)
        {
            //All cells for clarity: kernelX[0, 1], kernelX[1, 1] and kernelX[2, 1] are 0 and can be removed in expression
            int gx = kernelX[0, 0] * src[-srcStride - 1] + kernelX[0, 1] * src[-srcStride]  + kernelX[0, 2] * src[-srcStride + 1] +
                     kernelX[1, 0] * src[-1]             + kernelX[1, 1] * src[0]           + kernelX[1, 2] * src[1] +
                     kernelX[2, 0] * src[srcStride - 1]  + kernelX[2, 1] * src[srcStride]   + kernelX[2, 2] * src[srcStride + 1];

            //All cells for clarity: kernelY[1, 0], kernelY[1, 1] and kernelY[1, 2] are 0 and can be removed in expression
            int gy = kernelY[0, 0] * src[-srcStride - 1] + kernelY[0, 1] * src[-srcStride] + kernelY[0, 2] * src[-srcStride + 1] +
                     kernelY[1, 0] * src[-1]             + kernelY[1, 1] * src[0]          + kernelY[1, 2] * src[+1] +
                     kernelY[2, 0] * src[srcStride - 1]  + kernelY[2, 1] * src[srcStride]  + kernelY[2, 2] * src[srcStride + 1];

            if (useSum)
                g = Math.Abs(gy) + Math.Abs(gx);
            else
                g = (float)Math.Sqrt(gx * gx + gy * gy);

            if (takeMin)
                g = Math.Min(255, g);

            if (g > max)
                max = g;

            gradients[y, x] = g;
            src += srcBytesPerPixel;
            dst += dstBytesPerPixel;
        }
        srcCurrent += srcStride;
        dstCurrent += dstStride;
    }

    // make the second pass for intensity scaling and byte casting
    float factor = (float)255.0 / max;
    dstCurrent = (byte*)destination.ImageData.ToPointer() + dstStride * startY + dstBytesPerPixel * startX;
    // for each line
    for (int y = startY; y < stopY; y++)
    {
        // for each pixel
        dst = dstCurrent;
        for (int x = startX; x < stopX; x++)
        {
            *dst = (byte)(factor * gradients[y, x]);
            dst += dstBytesPerPixel;
        }
        dstCurrent += dstStride;
    }

    // draw black rectangle to remove those pixels, which were not processed
    // (this needs to be done for those cases, when filter is applied "in place" -
    // source image is modified instead of creating new copy)
    AForge.Imaging.Drawing.Rectangle(destination, rect, Color.Black);
}

Original issue reported on code.google.com by [email protected] on 15 Jun 2015 at 2:08