Computer Laboratory

Raspberry Pi

Edge Detection

0 Edge Detection

The algorithm we will look at in this tutorial is an edge detection algorithm, specifically an edge detection algorithm based on the Sobel operator. This algorithm works by calculating the gradient of the intensity of the image at each point, finding the direction of the change from light to dark and the magnitude of the change. This magnitude corresponds to how sharp the edge is.

1 The Sobel Operator

To calculate the gradient of each point in the image, the image is convolved with the Sobel Kernel. Convolution is done by moving the kernel across the image, one pixel at a time. At each pixel, the pixel and its neighbours are weighted by the corresponding value in the kernel, and summed to produce a new value. This operation is shown in the following diagram.

~~~~~~~~~~~~~~~~~~~            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
| k00 | k10 | k20 |            | i00 | i10 | i20 | i30 | ..
|~~~~~+~~~~~+~~~~~+            |~~~~~+~~~~~+~~~~~+~~~~~+~~~
| k01 | k11 | k21 |            | i01 | i11 | i21 | i31 | ..
|~~~~~+~~~~~+~~~~~+            |~~~~~+~~~~~+~~~~~+~~~~~+~~~
| k02 | k12 | k22 |            | i02 | i12 | i22 | i32 | ..
|~~~~~+~~~~~+~~~~~+            |~~~~~+~~~~~+~~~~~+~~~~~+~~~
                               | i03 | i13 | i23 | i33 | ..
Applying the kernel to the     |~~~~~+~~~~~+~~~~~+~~~~~+~~~
to the first pixel gives us    |  .. |  .. |  .. |  .. | ..	
an output value for that
pixel

     o00 = k11 * i00 + k12 * i01 + k21 * i10 + k22 * i11

we treat it as though the kernel has been overlaid onto the image, with the centre pixel of the kernel aligning we the first pixel in the image. Then we multiple each entry in the kernel by the value beneath it, and sum them to produce a single single output value from that pixel.

For the pixels on the boundary, we just ignore any entries in the kernel that fall outside.

Edge detection using the Sobel Operator applies two separate kernels to calculate the x and y gradients in the image. The length of this gradient is then calculated and normalised to produce a single intensity approximately equal to the sharpness of the edge at that position.

The kernels used for Sobel Edge Detection are shown below.

~~~~~~~~~~~~~
|-1 | 0 | 1 |
|~~~+~~~+~~~|
|-2 | 0 | 2 |
|~~~+~~~+~~~|
|-1 | 0 | 1 |
~~~~~~~~~~~~~
x gradient kernel

~~~~~~~~~~~~~
|-1 |-2 |-1 |
|~~~+~~~+~~~|
| 0 | 0 | 0 |
|~~~+~~~+~~~|
| 1 | 2 | 1 |
~~~~~~~~~~~~~
y gradient kernel

2 The Algorithm

Now the algorithm can be broken down into its constituent steps

1 Iterate over every pixel in the image
2 Apply the x gradient kernel
3 Apply the y gradient kernel
4 Find the length of the gradient using pythagoras' theorem
5 Normalise the gradient length to the range 0-255
6 Set the pixels to the new values

This illustrates the main steps, though we miss out some specifics such as what we do when we meet the boundary of the image. Now we will try implementing this in Python, and we will fill in any gaps in the algorithm as we go along

Create a new file for the python script, you can call it what you like. We will use the code from the first tutorial as a template for this one.

import imgproc
from imgproc import *

# import the maths module
import math


# open a webcam to take pictures, or load a sample image
camera = Camera(160, 120)

# Open a viewer window to display images
viewer = Viewer(160, 120, "Edge Detection")

# take a picture from the camera
img = camera.grabImage()

# display the image in the viewer
viewer.displayImage(img)


# iterate over each pixel in the image
for x in range(0, img.width):
	for y in range(0, img.height):

		# our code will go here


# display the image again
viewer.displayImage(img)

# delay for 5000 milliseconds to give us time to see the changes, then exit
waitTime(5000)

We need to do one tweak so that this will work for our edge detection. When we are iterating over the image, we will be changing value however when we move on to the next pixel, the just changed value will still be in the area of influence of our kernel. For this reason we need to make sure the image we read pixels from and the image we write pixels to are different, and we can do this by just created an empty image, which we will populate with the new values.

img_edge = Image(img.width, img.height)

Now we will read pixels from the original image, but write the edge values to img_edge

Before we start fiddling with the pixel values, we need to decide how we will handle the image boundaries. There are a few different ways to do this, one way is to check each pixel to see if it is on the boundary of the image, and ignore the non-existant pixels on that side. An alternative way would be to wrap the pixels, so that a non-existant pixel on the left edge actually maps to the pixel on the same row on the right edge. For simplicities sake I am simply going to ignore the boundary pixels, instead choosing to start one pixel in from the starting point, and 1 pixel away from the finishing edge. To do this, change the iterators to start at 1, and finish at img.width - 1 and img.height - 1 respectively. Now we can safely index adjacent pixels without falling off the edge of the image.

Right, lets move on to applying the gradients. For the x gradient, we need the column of pixels to the left and right, for each of the 6 pixels we index we will add up the red, green and blue intensities, multiply it by the kernel, and accumulated it into Gx. We will also do the same for the y gradient. As Gx and Gy both use the same pixels on occasions, it is best to apply the two kernels simultaneously.

		# initialise Gx to 0 and Gy to 0 for every pixel
		Gx = 0
		Gy = 0


		# top left pixel
		r, g, b = img[x - 1, y - 1]
		# intensity ranges from 0 to 765 (255 * 3)
		intensity = r + g + b   
		# accumulate the value into Gx, and Gy
		Gx += -intensity
		Gy += -intensity


		# now we do the same for the remaining pixels, left to right, top to bottom


		# remaining left column
		r, g, b = img[x - 1, y]
		Gx += -2 * (r + g + b)

		r, g, b = img[x - 1, y + 1]
		Gx += -(r + g + b)
		Gy += (r + g + b)


		# middle pixels
		r, g, b = img[x, y - 1]
		Gy += -2 * (r + g + b)

		r, g, b = img[x, y + 1]
		Gy += 2 * (r + g + b)


		# right column
		r, g, b = img[x + 1, y - 1]
		Gx += (r + g + b)
		Gy += -(r + g + b)

		r, g, b = img[x + 1, y]
		Gx += 2 * (r + g + b)

		r, g, b= img[x + 1, y + 1]
		Gx += (r + g + b)
		Gy += (r + g + b)

Now that the kernels are applied, Gx and Gy contain un-normalised values, which equal the relative length of the gradient in their respective axis. We wish to calculate the length of the gradient, and normalise it to a range suitable for displaying.

Lets calculate the relative length of each gradient. We can do this with pythagoras' theorem:

Pythagoras' theorem

sqr(a) + sqr(b) = sqr(c)
c = sqrt( sqr(a) + sqr(b) )

And in Python:

		length = math.sqrt((Gx * Gx) + (Gy * Gy))

This gives us an un-normalised length, to get the length normalised, we need to know the range of length

We know that each pixel's intensity value has a range of 0 to 765, due to it being the addition of 3 variables with range 0 to 255. Also, from looking at the kernels, we can see the maximum amount you can accumulate is +/- 4 intensities. Therefore the total range of Gx and Gy is between 0 and 3060. These values are then squared, added and square rooted, giving us a final range of 0 to 4328 ( sqrt(2 * sqr(3060)) ). We can renormalise by dividing by the upper bounds and multiplying by 255.

		# normalise the length
		length = length / 4328 * 255

		# convert the length to an integer
		length = int(length)

		# set edge image to a grayscale of the normalised length
		edge_img[x, y] = length, length, length

Right at the end of the code we need to change displayImage to display our new image rather than the original. Other than that, the script is finished, and we can test it to see the results. The script should produce a new grayscale image, highlighting edges from the original image.

This concludes this tutorial. Below are the code listings for this tutorial. The C code will actually continually capture images from the webcam, and you will find that on the Raspberry Pi it is capable of displaying the edge detected image in realtime.

3 Code Listing

3.1 Python Code

import imgproc
from imgproc import *

# also import the maths module for the square root function
import math


# open a webcam to take pictures
camera = Camera(160, 120)

# Open a viewer window to display images
viewer = Viewer(160, 120, "The first step")

# grab a couple of frames, just to flush the camera
camera.grabImage()
camera.grabImage()

# take a picture from the camera
img = camera.grabImage()

# create an empty image to show the edges
img_edge = Image(160, 120)


# display the image in the viewer
viewer.displayImage(img)


# iterate over each pixel in the image
for x in range(1, img.width - 1):
	for y in range(1, img.height - 1):
		# initialise Gx to 0 and Gy to 0 for every pixel
		Gx = 0
		Gy = 0


		# top left pixel
		r, g, b = img[x - 1, y - 1]
		# intensity ranges from 0 to 765 (255 * 3)
		intensity = (r + g + b)
		# accumulate the value into Gx, and Gy
		Gx += -intensity
		Gy += -intensity


		# now we do the same for the remaining pixels, left to right, top to bottom


		# remaining left column
		r, g, b = img[x - 1, y]
		Gx += -2 * (r + g + b)

		r, g, b = img[x - 1, y + 1]
		Gx += -(r + g + b)
		Gy += (r + g + b)


		# middle pixels
		r, g, b = img[x, y - 1]
		Gy += -2 * (r + g + b)

		r, g, b = img[x, y + 1]
		Gy += 2 * (r + g + b)


		# right column
		r, g, b = img[x + 1, y - 1]
		Gx += (r + g + b)
		Gy += -(r + g + b)

		r, g, b = img[x + 1, y]
		Gx += 2 * (r + g + b)

		r, g, b= img[x + 1, y + 1]
		Gx += (r + g + b)
		Gy += (r + g + b)


		# calculate the length of the gradient
		length = math.sqrt((Gx * Gx) + (Gy * Gy))

		# normalise the length of gradient to the range 0 to 255
		length = length / 4328 * 255

		# convert the length to an integer
		length = int(length)

		# draw the length in the edge image
		img_edge[x, y] = length, length, length


# display the edge image
viewer.displayImage(img_edge)

# delay for 5000 milliseconds to give us time to see the changes, then exit
waitTime(5000)


# end of the script

3.2 C Code

To compile the C code type this command in the terminal.

sudo gcc -O2 -o edge edge.c -limgproc

Now here is the C source


#include <stdlib.h>
#include <stdio.h>
#include "imgproc.h"
#include <math.h>

int main(int argc, char * argv[])
{
	// initialise the library
	init_imgproc();


	// open the webcam, with a capture resolution of width 320 and height 240
	Camera * cam = camOpen(320, 240);

	// create a new viewer of the same resolution with a caption
	Viewer * view = viewOpen(320, 240, "First");


	while(1){

		// capture an image from the webcam	
		Image * img = camGrabImage(cam);

		// Create an empty image
		Image * img_edge = imgNew(320, 240);


		// iterators
		unsigned int x, y;
		for(x = 1; x < img->width - 1; x++){
			for(y = 1; y < img->height - 1; y++){
				// initialise Gx and Gy to 0
				int Gx = 0;
				int Gy = 0;
				unsigned int intensity = 0;
				char * pixel;

				// Left column
				pixel = imgGetPixel(img, x - 1, y - 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += -intensity;
				Gy += -intensity;

				pixel = imgGetPixel(img, x - 1, y);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += -2 * intensity;

				pixel = imgGetPixel(img, x - 1, y + 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += -intensity;
				Gy += +intensity;

				// middle column
				pixel = imgGetPixel(img, x, y - 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gy += -2 * intensity;
			
				pixel = imgGetPixel(img, x, y + 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gy += +2 * intensity;

				// right column
				pixel = imgGetPixel(img, x + 1, y - 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += +intensity;
				Gy += -intensity;

				pixel = imgGetPixel(img, x + 1, y);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += +2 * intensity;
			
				pixel = imgGetPixel(img, x + 1, y + 1);
				intensity = pixel[0] + pixel[1] + pixel[2];
				Gx += +intensity;
				Gy += +intensity;

				// calculate the gradient length
				unsigned int length = (unsigned int)sqrt( (float)(Gx * Gx) + (float)(Gy * Gy) );

				// normalise the length to 0 to 255
				length = length / 17;

				// draw the pixel on the edge image
				imgSetPixel(img_edge, x, y, length, length, length);
			}
		}

		// display the edge image
		viewDisplayImage(view, img_edge);

		// now we will free the memory for the various objects
		imgDestroy(img_edge);
		imgDestroy(img);
	}

	// finally quit
	viewClose(view);
	camClose(cam);

	// finally we uninitialise the library
	quit_imgproc();
	return 0;
}