Thursday, October 25, 2007

Canny Edge Detection

The process of finding edges in the image is accomplished by using the Canny Edge Detection method. This method is well documented on the internet and an excellent tutorial has been written by Bill Green.

The three steps to a Canny Edge Detection are:
1. Blur the image
2. Calculate image gradients
3. Create a hysteresis

Blurring
The image is blurred using a Gaussian Mask in order to filter out noise such as dust particles and faint lines (which we are not interested in). A normal probability function in 2 dimensions is used to calculate the Mask. The mask is applied to every pixel of the image for which the edges of the mask still fall on pixels (i.e. there will be some space on the edges that the mask will not be applied to). The following equation calculates the normal probability function used in creating the mask:

p = exp(-(x^2 + y^2)/(2 * sigma^2)) / (2 * Pi * sigma^2)

The mask is created in the following code:

win = 1 + 2 * Ceiling(3 * sigma)
center = Floor(win / 2)
ReDim gk(win - 1, win - 1)
xsum = 0
For i = 0 To win - 1
*****For j = 0 To win - 1
**********b = i - center
**********d = j - center
**********fx = Exp(-0.5 * (b ^ 2 + d ^ 2) / (sigma * sigma)) / _
************(sigma * sigma * 2 * pi)
**********gk(i, j) = fx
**********xsum = xsum + fx
*****Next j
Next i
For i = 0 To win - 1
*****For j = 0 To win - 1
**********gk(i, j) = gk(i, j) / xsum
*****Next j
Next i

Sigma is generally a user input. Good values of sigma range from 0.5 to 1.4. The larger the sigma, the greater the blurring, but the time required to apply the mask increases exponentially. I recommend a value of 1.0 as it seems to give adequate blurring and can be done in little time.

The mask is applied in the following code (note: idat() is the array containing the grayscale pixels values where each pixel is represented by a single byte):

For i = 0 To numberofpixels - 1
*****xdot = 0.0
*****xsum = 0.0
*****For cc = -center To center
**********For j = -center To center
***************If (((i + row) Mod row) + cc) >= 0 And _
*****************(((i + row) Mod row) + cc) _
*****************And (i + j * row) > 0 And (i + j * row)
********************xdot = xdot + idat(i + cc + j * row) * _
**********************gk(center + cc, center + j)
********************xsum = xsum + gk(center + cc, center + j)
***************End If
**********Next j
*****
Next cc
*****gdat(i) = (xdot / xsum) * boostblurfactor + 0.5
Next i
For i = 1 To numberofpixels - 1
*****idat(i) = gdat(i)
Next i

The boostblurfactor in the code above is just a constant for changing the amount of blur effect. I use a value of 0.9 because that is the value that was used in the source code I modeled my code from. The constant "row" is equal to the number of pixels in each row of the original image.

Gradients
To find the edges, it is not enough to know where the image is intensely dark or light, but rather to know where there is a very sudden change from dark to light (or light to dark). The gradient is a measure of the amount of change around each pixel. The following code calculates the gradients at each pixel (note that derx() and dery() are arrays that hold the values of the gradients in those directions and the magnitude() array holds the combined gradient for the pixel):

For i = 0 To numberofpixels - 1
*****'Compute x derivative
*****If i = 0 Or (i Mod row) = 0 Then
**********c = idat(i + 1)
**********d = idat(i)
**********delx(i) = c - d
*****ElseIf ((i + 1) Mod row) = 0 Then
**********c = idat(i)
**********d = idat(i - 1)
**********delx(i) = c - d
*****Else
**********c = idat(i + 1)
**********d = idat(i - 1)
**********delx(i) = c - d
*****End If

*****'Compute y derivative
*****If i <>**********c = idat(i + row)
**********d = idat(i)
**********dely(i) = c - d
*****ElseIf i > numberofpixels - row - 1 Then
**********c = idat(i)
**********d = idat(i - row)
**********dely(i) = c - d
*****Else
**********c = idat(i + row)
**********d = idat(i - row)
**********dely(i) = c - d
*****End If
*****If i Mod (467 * row) = 0 Then Application.DoEvents()
Next i

'Calculate the magnitude of the gradient
For i = 0 To numberofpixels - 1
*****sq1 = delx(i) * delx(i)
*****sq2 = dely(i) * dely(i)
*****magnitude(i) = 0.5 + Sqrt(sq1 + sq2)
Next i

After calculating the gradient, we need to suppress any pixel that is not a maximum. We will use three values to indicate pixel status as:
1. No Edge (0)
2. Possible Edge (128)
3. Definite Edge (255)

The pixels around the edges of the image will be set to No Edge. For the remaining pixels, a simple test shows whether the pixel in question is a possible edge or not, but the eight different directions from which the gradient may be coming create a complicated set of conditions for the test. The code is:

Dim i As Integer
Dim xperp, yperp, z1, z2, mag1, mag2 As Double

'Initialize parameters
ReDim idat(bmsize - 1) 'idat() is now the result array
noedge = 0
possedge = 128
edge = 255

'Zero the Result() Image Edges
For i = 0 To bmsize - 1
***** If i <> (bmsize - row - 1) Or (i Mod row) = 0 _
**********Or ((i + 1) Mod row) = 0 Or magnitude(i) = 0 Then
********** idat(i) = noedge
***** Else
********** xperp = -delx(i) / magnitude(i)
********** yperp = dely(i) / magnitude(i)
***** End If

***** If delx(i) >= 0 And dely(i) >= 0 And delx(i) >= dely(i) _
**********And i > row And i < (bmsize - row - 1) And (i Mod row) <> 0 _
**********And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i - 1)
********** z2 = magnitude(i - row - 1)
********** mag1 = (magnitude(i) - z1) * xperp + (z2 - z1) * yperp

********** 'Right Point
********** z1 = magnitude(i + 1)
********** z2 = magnitude(i + row + 1)
********** mag2 = (magnitude(i) - z1) * xperp + (z2 - z1) * yperp
***** ElseIf delx(i) >= 0 And dely(i) >= 0 And delx(i) <> row _
**********And i < (bmsize - row - 1) And (i Mod row) <> 0 _
**********And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i - row)
********** z2 = magnitude(i - row - 1)
********** mag1 = (z1 - z2) * xperp + (z1 - magnitude(i)) * yperp

**********'Right Point
********** z1 = magnitude(i + row)
********** z2 = magnitude(i + row + 1)
********** mag2 = (z1 - z2) * xperp + (z1 - magnitude(i)) * yperp
***** ElseIf delx(i) >= 0 And dely(i) <>= -dely(i) And i > row _
**********And i < (bmsize - row - 1) And (i Mod row) <> 0 _
**********And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i - 1)
********** z2 = magnitude(i + row - 1)
********** mag1 = (magnitude(i) - z1) * xperp + (z1 - z2) * yperp

**********'Right Point
********** z1 = magnitude(i + 1)
********** z2 = magnitude(i - row + 1)
********** mag2 = (magnitude(i) - z1) * xperp + (z1 - z2) * yperp
***** ElseIf delx(i) >= 0 And dely(i) <> row And i < (bmsize - row - 1) _ **********And (i Mod row) <> 0 And ((i + 1) Mod row) <> 0 Then
**********'Left Point
**********z1 = magnitude(i + row)
**********z2 = magnitude(i + row - 1)
**********mag1 = (z1 - z2) * xperp + (magnitude(i) - z1) * yperp

********** 'Right Point
********** z1 = magnitude(i - row)
********** z2 = magnitude(i - row + 1)
********** mag2 = (z1 - z2) * xperp + (magnitude(i) - z1) * yperp
***** ElseIf delx(i) <>= 0 And -delx(i) >= dely(i) And i > row _
**********And i < (bmsize - row - 1) And (i Mod row) <> 0 _
**********And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i + 1)
********** z2 = magnitude(i - row + 1)
********** mag1 = (z1 - magnitude(i)) * xperp + (z2 - z1) * yperp

**********'Right Point
********** z1 = magnitude(i - 1)
********** z2 = magnitude(i + row - 1)
********** mag2 = (z1 - magnitude(i)) * xperp + (z2 - z1) * yperp
***** ElseIf delx(i) <>= 0 And -delx(i) <> row And i < (bmsize - row - 1) **********And (i Mod row) <> 0 And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i - row)
********** z2 = magnitude(i - row - 1)
********** mag1 = (z2 - z1) * xperp + (z1 - magnitude(i)) * yperp

********** 'Right Point
********** z1 = magnitude(i + row)
********** z2 = magnitude(i + row - 1)
********** mag2 = (z2 - z1) * xperp + (z1 - magnitude(i)) * yperp
***** ElseIf delx(i) <>= -dely(i) And i > row And i < (bmsize - row - 1) _ **********And (i Mod row) <> 0 And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i + 1)
********** z2 = magnitude(i + row + 1)
********** mag1 = (z1 - magnitude(i)) * xperp + (z1 - z2) * yperp

********** 'Right Point
********** z1 = magnitude(i - 1)
********** z2 = magnitude(i - row - 1)
********** mag2 = (z1 - magnitude(i)) * xperp + (z1 - z2) * yperp
***** ElseIf i > row And i < (bmsize - row - 1) And (i Mod row) <> 0 _
**********And ((i + 1) Mod row) <> 0 Then
********** 'Left Point
********** z1 = magnitude(i + row)
********** z2 = magnitude(i + row + 1)
********** mag1 = (z2 - z1) * xperp + (magnitude(i) - z1) * yperp

**********'Right Point
********** z1 = magnitude(i - row)
********** z2 = magnitude(i - row - 1)
********** mag2 = (z2 - z1) * xperp + (magnitude(i) - z1) * yperp
***** End If

***** 'Now determine if the point is a maximum
*****If mag1 > 0.0 Or mag2 > 0.0 Then
********** idat(i) = noedge
***** ElseIf mag2 = 0.0 Then
********** idat(i) = noedge
***** Else
********** idat(i) = possedge
***** End If
***** If i Mod (467 * row) = 0 Then Application.DoEvents()
Next i

After filtering out all of the points that probably aren't edges, a hysteresis is created on the magnitude of the gradients for the remaining points. We choose a low and high threshold. If the magnitude of the gradient is lower than the low threshold, it cannot be an edge, but if it is above the high threshold, we assume it is an edge. Anything in between is considered a possible edge. As points are eliminated from the list of possible edge points, the image changes. We have to continue our comparison of gradient magnitudes and process of eliminating points until the image stops changing.

The code for generating the hysteresis and creating the edge data is given here:
Dim i, j, highc, hist(), edges, hight, lowt, tpos As Integer
Dim thigh As Single
Dim tlow As Single
Dim mmag As Short
Dim x() As Integer = {1, 1, 0, -1, -1, -1, 0, 1}
Dim y() As Integer = {0, 1, 1, 1, 0, -1, -1, -1}
Dim unchanged As Boolean

'Initialize parameters
ReDim hist(32768)
thigh = 0.9
tlow = 0.5
hist.Initialize()

For i = 0 To bmsize - 1
***** 'Initialize edges of image to noedge
***** If idat(i) = possedge Then idat(i) = possedge
**********Else idat(i) = noedge
*****If i <> bmsize - 1 - row Or ((i + row) Mod row) = 0
**********Or ((i + row + 1) Mod row) = 0 Then idat(i) = noedge

***** 'Creates the histogram of the magnitudes in the image
***** If idat(i) = possedge Then _
**********hist(magnitude(i)) = hist(magnitude(i)) + 1
Next i

'Calculate the number of pixels that pass the suppression
For i = 1 To 32767
***** If hist(i) <> 0 Then
********** mmag = i
********** edges = edges + hist(i)
***** End If
Next i

highc = edges * thigh + 0.5

I found this comment from Jim Carol's source code useful: Compute the high threshold value as the (100 * thigh) percentage point in the magnitude of the gradient histogram of all the pixels that passes non-maximal suppression. Then calculate the low threshold as a fraction of the computed high threshold value. John Canny said in his paper "A Computational Approach to Edge Detection" that "The ratio of the high to low threshold in the implementation is in the range two or three to one." That means that in terms of this implementation, we should choose tlow ~= 0.5 or 0.33333.

i = 1
edges = hist(1)
Do While i < (mmag - 1) And edges <>
*****i = i + 1
*****edges = edges + hist(i)
Loop
hight = i
lowt = hight * tlow + 0.5

'Looks for pixels above the hight to locate edges we will need to run over the entire image until it stops changing. In this first step we seed the image with everything that is clearly an edge by being above the high threshold.

For i = 0 To bmsize - 1
*****If idat(i) = possedge And magnitude(i) >= hight Then idat(i) = edge
Next i

'run over image until it stops changing
unchanged = False
Do While unchanged = False
***** unchanged = True
*****For j = 0 To bmsize - 1
**********If idat(j) = edge Then
***************For i = 0 To 7
******************** tpos = j - y(i) * row + x(i)
********************If idat(tpos) = possedge _
*************************And magnitude(tpos) > lowt Then
************************* unchanged = False
******************** ***** idat(tpos) = edge
********************End If
***************Next i
**********End If
*****Next j
Loop

'Set any remaining possible edges to non-edge
For i = 0 To bmsize - 1
*****If idat(i) <> edge Then idat(i) = noedge
Next i

'idat() now has the edge data stored to it

An edge image can be exported to bmp format by assigning the value of each point in idat() to all three of the color values for each pixel of a bitmap image with the original dimensions as explained in the bitmap tutorial. Edge images have exactly 2 colors: black and white. I will occasionally use an edge image for troubleshooting purposes.

With the edge image, it is now possible to look for shapes in our bitmap. We will use this ability to find the sprockets holes in the film. Once the sprockets are known, the frames are easily picked off.

Go on to Sprockets

No comments: