EE4208 Canny edge detector

canny.py

import numpy as np
import matplotlib.pyplot as plot
import matplotlib.cm as cm
import sys
import pylab
import math
import scipy

# Mostly implemented by Matt Cumming

range_inc = lambda start, end: range(start, end+1) #Because this is easier to write

def convolve(image, mask):
        
	# Convolute the mask with the image # Yes, technically this should be a vector operation, but I'm totally not hardcore enough for that yet
	width = image.shape[1]
	height = image.shape[0]
	w_range = int(math.floor(mask.shape[0]/2))

	res_image = np.zeros((height, width))


	# Iterate over every pixel that can be covered by the mask
	for i in range(w_range,width-w_range):
		for j in range(w_range,height-w_range):
			# Then convolute with the mask 
			for k in range_inc(-w_range,w_range):
				for h in range_inc(-w_range,w_range):
					res_image[j, i] += mask[w_range+h,w_range+k]*image[j+h,i+k]
	return res_image

def Gaussian(A, sigma,x,y):

        return(A*math.exp(-((x**2)+(y**2))/(2*sigma**2)))

def CreateMask(A,sigma):
        TotalA=0
        #width must be at least 5*sigma and must be an odd integer
        width=int(math.floor((sigma*5)+0.5))
        if width%2==0:
                width+=1

        w_range = int(math.floor(width/2))

        GaussianMask = np.zeros((width, width))
        min_array=99999 #value higher than possible minimum

        for x in range_inc(-w_range,w_range):
                for y in range_inc(-w_range,w_range):
                        GaussianMask[x,y]=(Gaussian(A,sigma,x,y))
                        if GaussianMask[x,y]<min_array:
                                min_array=GaussianMask[x,y]

        #to increase computational efficiecny, use approximate integer values by normalising lowest value to 1
        
        for x in range_inc(-w_range,w_range):
                for y in range_inc(-w_range,w_range):
                        GaussianMask[x,y]=int(math.floor(((GaussianMask[x,y])/min_array)+0.5))                  

        return(GaussianMask)

def Fat(Image):

        sobel_x = np.array([[1,0,-1], [2,0,-2], [1,0,-1]])
        sobel_y = sobel_x.T

        gradx=convolve(smoothed_image,sobel_x)
        grady=convolve(smoothed_image,sobel_y)

        sobeloutmag = scipy.hypot(gradx, grady)
        sobeloutdir = scipy.arctan2(grady, gradx)

        width = sobeloutdir.shape[0]
        height = sobeloutdir.shape[1]

        #notes taken from tutorial at http://pythongeek.blogspot.sg/2012/06/canny-edge-detection.html
        for x in range(width-1):
            for y in range(height-1):
                if (sobeloutdir[x][y]<22.5 and sobeloutdir[x][y]>=0) or (sobeloutdir[x][y]>=157.5 and sobeloutdir[x][y]<202.5) or (sobeloutdir[x][y]>=337.5 and sobeloutdir[x][y]<=360):
                    sobeloutdir[x][y]=0
                elif (sobeloutdir[x][y]>=22.5 and sobeloutdir[x][y]<67.5) or (sobeloutdir[x][y]>=202.5 and sobeloutdir[x][y]<247.5):
                    sobeloutdir[x][y]=135
                elif (sobeloutdir[x][y]>=67.5 and sobeloutdir[x][y]<112.5)or (sobeloutdir[x][y]>=247.5 and sobeloutdir[x][y]<292.5):
                    sobeloutdir[x][y]=90
                else:
                    sobeloutdir[x][y]=45

        return([sobeloutmag, sobeloutdir])

def Thin(ImageMag,ImageDir):

        thin_image=np.copy(ImageMag)

        width = ImageMag.shape[0]
        height = ImageMag.shape[1]
     
        for x in range(1, width-1):
            for y in range(1, height-1):
                if ImageDir[x][y]==0:
                    if (ImageMag[x][y]<=ImageMag[x][y+1]) or (ImageMag[x][y]<=ImageMag[x][y-1]):
                        thin_image[x][y]=0
                elif ImageDir[x][y]==45:
                    if (ImageMag[x][y]<=ImageMag[x-1][y+1]) or (ImageMag[x][y]<=ImageMag[x+1][y-1]):
                        thin_image[x][y]=0
                elif ImageDir[x][y]==90:
                    if (ImageMag[x][y]<=ImageMag[x+1][y]) or (ImageMag[x][y]<=ImageMag[x-1][y]):
                        thin_image[x][y]=0
                else:
                    if (ImageMag[x][y]<=ImageMag[x+1][y+1]) or (ImageMag[x][y]<=ImageMag[x-1][y-1]):
                        thin_image[x][y]=0

        return(thin_image)


print "start"
#Read in the user specified standard deviation and file name
filename = "cana.raw"
sigma=1.4
#filename=raw_input("Enter the file to open: ")
#sigma= float(raw_input("Enter the required standard deviation: "))

#Load the raw file
f = open(filename,'rb') # switch to command line args later
#Because the byte order is weird
a = f.read(1)
b = f.read(1)
#First line is rows
rows = int((b+a).encode('hex'), 16)
a = f.read(1)
b = f.read(1)
#Second line is columns
cols = int((b+a).encode('hex'), 16)
#Last byte is encoding, but we're just going to ignore it
f.read(1)
#And everything else is 8 bit encoded, so let's load it into numpy and display it with matplotlib
bin_image = np.fromstring(f.read(), dtype=np.uint8)
#Change the shape of the array to the actual shape of the picture
bin_image.shape = (cols, rows)

fig = pylab.figure()
#Display the original image
fig.add_subplot(1,4,1)
pylab.imshow(bin_image, cmap=cm.gray)

#smooth the image
GaussianMask=CreateMask(1,sigma) #A=1 as it is not needed


smoothed_image = convolve(bin_image, GaussianMask)

#Display the smoothed image
fig.add_subplot(1,4,2)
pylab.imshow(smoothed_image, cmap=cm.gray)

#get the fat edged image
FatImage=Fat(smoothed_image)
fat_image_mag=FatImage[0]
fat_image_dir=FatImage[1]

#Display Fat Image
fig.add_subplot(1,4,3)
pylab.imshow(fat_image_mag, cmap=cm.gray)

#get the thin edged image
thin_image = Thin(fat_image_mag,fat_image_dir)

#Display Thin Image
fig.add_subplot(1,4,4)
pylab.imshow(thin_image, cmap=cm.gray)

pylab.show()

print "done"