# How can I improve the Canny edge detection code?

I have converted this Matlab code for Canny edge detection into C# given that the Matlab version gives out a good output.

The output from Matlab code is as follows:

And my C# version's outputs are as follows:

Gaussian blur output.

Left: Horizontal Sobel output; Right: Vertical Sobel output.

Left: using Spatial domain Convolution; Right: using FFT Convolution (slower).

As we can see, my C# version's outputs are inferior. I couldn't fix them by varying the high and low threshold values.

On which step should I consider applying an improvement?

.

Source Code

public class Canny
{
public static double[,] Apply(double[,] image)
{
double dTresLo = 0;
double dTresHi = 0;

dTresLo = 0.075;
dTresHi = 0.175;

double[,] dGaussKernel = new double[,]
{   {2.0/159.0,  4.0/159.0,  5.0/159.0, 4.0/159.0, 2.0/159.0},
{4.0/159.0,  9.0/159.0, 12.0/159.0, 9.0/159.0, 4.0/159.0},
{5.0/159.0, 12.0/159.0, 15.0/159.0,12.0/159.0, 5.0/159.0},
{4.0/159.0,  9.0/159.0, 12.0/159.0, 9.0/159.0, 4.0/159.0},
{2.0/159.0,  4.0/159.0,  5.0/159.0, 4.0/159.0, 2.0/159.0} };

double[,] dSobelKernelX = new double[,]{{1, 0, -1},
{2, 0, -2},
{1, 0, -1} };

double[,] dSobelKernelY = new double[,]{{1,   2,   1},
{0,   0,   0},
{-1, -2,  -1}};

double[,] dGausOutput = CannyHelper.LinearConvolutionSpatial(image, dGaussKernel);
double[,] matSobelOutX = CannyHelper.LinearConvolutionSpatial(dGausOutput, dSobelKernelX);
double[,] matSobelOutY = CannyHelper.LinearConvolutionSpatial(dGausOutput, dSobelKernelY);

// Calculate directions/orientations
// Adjustment for negative directions, making all directions positive
double[,] matDirRad = new double[image.GetLength(0), image.GetLength(1)];
for (int x = 0; x < matDirRad.GetLength(0); x++)
{
for (int y = 0; y < matDirRad.GetLength(1); y++)
{
double atan = (matSobelOutY[x, y] / matSobelOutX[x, y]) * 180 / Math.PI;

if (atan < 0)
{
matDirRad[x, y] = atan + 360;
}
else
{
}
}
}

// Adjusting directions to nearest 0, 45, 90, or 135 degree
double[,] arah2 = new double[image.GetLength(0), image.GetLength(1)];
for (int i = 0; i < arah2.GetLength(0); i++)
{
for (int j = 1; j < arah2.GetLength(1); j++)
{
{
arah2[i, j] = 0;
}
{
arah2[i, j] = 45;
}
{
arah2[i, j] = 90;
}
{
arah2[i, j] = 135;
}
}
}

// Calculate magnitude
double[,] matSobelImg = new double[image.GetLength(0), image.GetLength(1)];
for (int x = 0; x < matSobelImg.GetLength(0); x++)
{
for (int y = 0; y < matSobelImg.GetLength(1); y++)
{
matSobelImg[x, y] = Math.Sqrt((matSobelOutX[x, y] * matSobelOutX[x, y]) + (matSobelOutY[x, y] * matSobelOutY[x, y]));
}
}

// Non-Maximum Supression
double[,] matSuppImg = new double[image.GetLength(0), image.GetLength(1)];
for (int i = 1; i < matSuppImg.GetLength(0)-1; i++)
{
for (int j = 1; j < matSuppImg.GetLength(1)-1; j++)
{
if (arah2[i, j]==0)
{
if (matSobelImg[i, j] == Math.Max(Math.Max(matSobelImg[i, j], matSobelImg[i, j + 1]), matSobelImg[i, j - 1]))
{
matSuppImg[i, j] = 1;
}
else
{
matSuppImg[i, j] = 0;
}
}
else if (arah2[i, j]==45)
{
if(matSobelImg[i, j] == Math.Max(Math.Max(matSobelImg[i, j], matSobelImg[i+1,j-1]), matSobelImg[i-1,j+1]))
{
matSuppImg[i, j] = 1;
}
else
{
matSuppImg[i, j] = 0;
}
}
else if (arah2[i, j]==90)
{
if(matSobelImg[i, j] == Math.Max(Math.Max(matSobelImg[i, j], matSobelImg[i+1,j]), matSobelImg[i-1,j]))
{
matSuppImg[i, j] = 1;
}
else
{
matSuppImg[i, j] = 0;
}
}
else if (arah2[i, j]==135)
{
if(matSobelImg[i, j] == Math.Max(Math.Max(matSobelImg[i, j], matSobelImg[i+1,j+1]), matSobelImg[i-1,j-1]))
{
matSuppImg[i, j] = 1;
}
else
{
matSuppImg[i, j] = 0;
}
}
}
}

for (int i = 0; i < matSuppImg.GetLength(0); i++)
{
for (int j = 0; j < matSuppImg.GetLength(1); j++)
{
matSuppImg[i, j] = matSuppImg[i,j] * matSobelImg[i,j];
}
}

// Hysteresis Thresholding
dTresLo = dTresLo * CannyHelper.Max(matSuppImg);
dTresHi = dTresHi * CannyHelper.Max(matSuppImg);

double[,] T_res = new double[image.GetLength(0), image.GetLength(1)];
for (int i = 1; i < T_res.GetLength(0) - 1; i++)
{
for (int j = 1; j < T_res.GetLength(1) - 1; j++)
{
if (matSuppImg[i, j] < dTresLo)
T_res[i, j] = 0;
else if (matSuppImg[i, j] > dTresHi)
T_res[i, j] = 1;
//%Using 8-connected components
else if (matSuppImg[i + 1, j] > dTresHi
|| matSuppImg[i - 1, j] > dTresHi
|| matSuppImg[i, j + 1] > dTresHi
|| matSuppImg[i, j - 1] > dTresHi
|| matSuppImg[i - 1, j - 1] > dTresHi
|| matSuppImg[i - 1, j + 1] > dTresHi
|| matSuppImg[i + 1, j + 1] > dTresHi
|| matSuppImg[i + 1, j - 1] > dTresHi)
T_res[i, j] = 1;
}
}

return T_res;
}
}


By the way, the following is my convolution function look like:

    public static double[,] ConvolutionSpatial(double[,] paddedImage, double[,] mask, double offset)
{
double min = 0.0;
double max = 1.0;

double[,] convolve = new double[imageWidth, imageHeight];

for (int x = 0; x < imageWidth; x++)
{
for (int y = 0; y < imageHeight; y++)
{

convolve[x, y] = Math.Min(Math.Max((sum / factor) + offset, min), max);

string str = string.Empty;
}
}

return convolve;
}

{
double sum = 0;

for (int x = startX; x < (startX + maskWidth); x++)
{
for (int y = startY; y < (startY + maskHeight); y++)
{
double msk = mask1[maskWidth - x + startX - 1, maskHeight - y + startY - 1];
sum = sum + (img * msk);
}
}

return sum;
}

public static double GetFactor(double[,] kernel)
{
double sum = 0.0;

int width = kernel.GetLength(0);
int height = kernel.GetLength(1);

for (int x = 0; x < width; x++)
{
for (int y = 0; y < height; y++)
{

sum += kernel[x, y];
}
}

return (sum == 0) ? 1 : sum;
}


Does anyone sense anything wrong in this?