import numpy as np
import matplotlib.pyplot as plt

#Read image
I = plt.imread('umn.jpg')

#Convert to grayscale
G = 0.2989*I[:,:,0] + 0.5870*I[:,:,1] + 0.1140*I[:,:,2] #Gray = 0.2989*Red + 0.5870*Green + 0.1140*Blue

#Compute gradient
GE = G[2:,1:-1]     #Neighboring pixel east
GW = G[:-2,1:-1]    #Neighboring pixel west
GN = G[1:-1,2:]     #Neighboring pixel north
GS = G[1:-1,:-2]    #Neighboring pixel south
grad = (GE - GW)**2 + (GS - GN)**2

#Normalize gradient and boost the response (for visualization)
grad = grad/np.max(grad)
boost_factor = 15
edges =(255*np.minimum(grad,1/boost_factor)*boost_factor).astype(int)

#Splice together with the original image
J = I[1:-1,1:-1,:].copy()
h = int(J.shape[1]/2)
J[:,h:,:] = edges[:,h-1::-1,np.newaxis]

#Subsample (the image is rather high resolution)
J = J[::2,::2,:]

#Write the image to a new file 
plt.imsave('umn_edges.jpg',J)