Last night I stayed up late wrestling with an artificial vision project. The goal was simple on paper: isolate the silhouette of some mechanical parts passing on a conveyor belt. The usual. The problem is that the lighting in the industrial warehouse where this is set up is an absolute disaster. Using my rudimentary threshold edge detector was no good. It either gave me noise everywhere (it looked like the image had measles), or it lost the actual outlines of the part if I raised the filter's strictness too much.

After reviewing old notes and diving through forums for a while, I decided to implement the famous Canny edge detector. We already saw in the previous article how to fight with BMP headers in C, so this time we're going straight to the fun part: the pure mathematical algorithm in C.

Why Canny and not a simple Sobel?

Up until now I had been playing with the Sobel operator, which is not bad to start with, but it leaves you with edge lines thicker than a permanent marker. Canny, developed by John F. Canny back in 1986 (yes, it's old but still the king), solves this brilliantly and, besides, handles noise like a champ.

The Canny algorithm is not a simple function to which you pass an image and it returns another. It is a process in several very distinct stages:

  1. Gaussian Filter: Smooths the image through spatial convolution to remove high-frequency noise (the hateful salt and pepper effect).
  2. Gradient Calculation: Finds the intensity and direction of the edge. This is where we usually put in the Sobel masks to calculate the derivatives on the X and Y axes.
  3. Non-maximum suppression: Thins the detected edges. It basically ensures that the edge is only one pixel thick by solely selecting the local maximum in the direction of the gradient. Pure magic.
  4. Hysteresis Thresholding: This is the part I like the most. It uses two thresholds (high and low) to decide which pixels are a real edge and which are just connected noise.

Implementing Canny from scratch in C

I'm not going to paste the entire source code of the project because it would blow up your browser, but I want to show you the core of the most complex part: calculating the gradient and the directions.

Assuming we already have our image in memory, this would be the snippet where we calculate the magnitude and direction of the edges:

// We assume that imagen_suavizada is a 1D array for our 2D grayscale image
// ancho and alto are the dimensions obtained from the header
int Gx, Gy;
double magnitud;
double direccion;

// We traverse the image avoiding the external borders so we don't go out of bounds
for (int y = 1; y < alto - 1; y++) {
    for (int x = 1; x < ancho - 1; x++) {
        // Sobel operator on the X axis
        Gx = (imagen_suavizada[(y-1)*ancho + (x+1)] + 2*imagen_suavizada[y*ancho + (x+1)] + imagen_suavizada[(y+1)*ancho + (x+1)])
           - (imagen_suavizada[(y-1)*ancho + (x-1)] + 2*imagen_suavizada[y*ancho + (x-1)] + imagen_suavizada[(y+1)*ancho + (x-1)]);

        // Sobel operator on the Y axis
        Gy = (imagen_suavizada[(y-1)*ancho + (x-1)] + 2*imagen_suavizada[(y-1)*ancho + x] + imagen_suavizada[(y-1)*ancho + (x+1)])
           - (imagen_suavizada[(y+1)*ancho + (x-1)] + 2*imagen_suavizada[(y+1)*ancho + x] + imagen_suavizada[(y+1)*ancho + (x+1)]);

        // We calculate the magnitude of the vector (using a fast approximation to save CPU cycles)
        magnitud = abs(Gx) + abs(Gy); 

        // Direction (angle in radians)
        direccion = atan2((double)Gy, (double)Gx);

        // Here we would save the results in auxiliary matrices
        // matriz_magnitud[y*ancho + x] = magnitud;
        // matriz_angulo[y*ancho + x] = direccion;
    }
}

That code is a nested for loop that doesn't hold much mystery, but it's the heart that makes the system beat. You have to be very careful with the indices.

After this step would come the non-maximum suppression. That is, look in the direction that atan2 has given us and check if the current pixel has a greater value than its two immediate neighbors on that axis. If it does not, we suppress it and set it to zero (black). This is what achieves that fine, sharp, professional edge. Finally, we apply hysteresis to it and cross our fingers that we don't pick up noise.

Compiling the creature

I pull the GNU compiler (GCC) directly from the terminal, remembering to link the math library:

gcc -Wall -O2 canny.c -o canny -lm
./canny entrada.bmp salida.bmp 50 150

Note for the clueless: Those 50 and 150 are my low and high threshold values for the hysteresis algorithm. You're going to have to play with them depending on the contrast, the lighting, and the quality of the camera you take the photos with.

Reflections on artificial vision

Programming these things from scratch based on one-dimensional arrays in C is painful, I won't lie to you, but it turns out to be tremendously educational. I'm aware that there are very powerful libraries in development, like OpenCV (which just released a version a couple of years ago), that solve all this hassle for you with a single call to the cvCanny function. But, honestly, understanding what happens under the hood of memory and dealing with derivatives gives me much more peace of mind when the control system later fails in production.

It's fascinating what we can achieve by brute-force processing enormous matrices of pixels. If processors keep increasing in speed at this hellish rate (we're already bordering on 3GHz in the brand new dual-core Intel Core 2 Duo processors), I wonder if in a few years we'll be able to apply these mathematical algorithms to real-time, high-resolution video streams, right out of the box and without burning the CPU. Only time will tell.