import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

horizon = np.load('Penobscot_HorB.npy')

plt.figure()
plt.imshow(-1 * horizon, aspect=0.5)
plt.show()

vmin = np.percentile(horizon[horizon > 0], 1)
vmax = np.percentile(horizon[horizon > 0], 99)
print 'min {0}, max {1}'.format(vmin, vmax)

import scipy.ndimage

def do_sobel(image):
    dx = scipy.ndimage.sobel(image, 0)  # horizontal derivative
    dy = scipy.ndimage.sobel(image, 1)  # vertical derivative
    mag = np.hypot(dx, dy)      # magnitude
    mag *= 255.0 / np.max(mag)  # normalize
    return mag

mag = do_sobel(horizon)

plt.imshow(mag, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=3)
plt.show()

plt.figure(figsize=(15,10))
plt.imshow(-1*horizon, aspect=0.5, cmap="cubehelix", vmin=-vmax, vmax=-vmin)
plt.colorbar(shrink=0.75)  # shrink makes the colourbar a bit shorter
plt.show()

# You might need to install this first: pip install scikit-image
from skimage import filter

# Generate noisy image of a square
im = np.zeros((128, 128))
im[32:-32, 32:-32] = 1

im = scipy.ndimage.rotate(im, 15, mode='constant')
im = scipy.ndimage.gaussian_filter(im, 4)
im += 0.2 * np.random.random(im.shape)

# Compute the Canny filter for two values of sigma
edges1 = filter.canny(im)
edges2 = filter.canny(im, sigma=3)

# display results
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(12, 5))

ax1.imshow(im, cmap=plt.cm.jet)
ax1.axis('off')
ax1.set_title('noisy image', fontsize=20)

ax2.imshow(edges1, cmap=plt.cm.gray)
ax2.axis('off')
ax2.set_title('Canny filter, $\sigma=1$', fontsize=20)

ax3.imshow(edges2, cmap=plt.cm.gray)
ax3.axis('off')
ax3.set_title('Canny filter, $\sigma=3$', fontsize=20)

fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
                    bottom=0.02, left=0.02, right=0.98)

plt.show()

edge = filter.canny(horizon)

plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.imshow(horizon, aspect=0.5, cmap='cubehelix_r', vmin=vmin, vmax=vmax)
plt.subplot(122)
plt.imshow(edge, aspect=0.5, cmap='Greys', vmin=0, vmax=1)
plt.show()

edge_2  = filter.canny(horizon, sigma=2)
edge_4  = filter.canny(horizon, sigma=4)
edge_8  = filter.canny(horizon, sigma=8)
edge_16 = filter.canny(horizon, sigma=16)

plt.figure(figsize=(20,12))
plt.subplot(221)
plt.imshow(edge_2, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.text(10, 50, "sigma = 2", size=24, color="blue")
plt.subplot(222)
plt.imshow(edge_4, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.text(10, 50, "sigma = 4", size=24, color="blue")
plt.subplot(223)
plt.imshow(edge_8, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.text(10, 50, "sigma = 8", size=24, color="blue")
plt.subplot(224)
plt.imshow(edge_16, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.text(10, 50, "sigma = 16", size=24, color="blue")
plt.show()

pixels = float(horizon.size)
print np.sum(edge_8)/pixels, np.sum(edge_16)/pixels

images = []
for i in range(16):
    images.append(filter.canny(horizon, sigma=i))
    if np.sum(images[-1]) / pixels < 0.028:
        result = images[-1]
        break

plt.figure(figsize=(8,5))
plt.imshow(result, aspect=0.5, cmap='Greys', vmin=0, vmax=1)
plt.title('Canny filter, sigma ={0}'.format(i), size=20)
plt.show()

from IPython.html import widgets

image = np.copy(horizon)

def canny_demo(**kwargs):
    edges = filter.canny(image, **kwargs)
    plt.figure(figsize=(9,6))
    plt.imshow(horizon, aspect=0.5, cmap='cubehelix_r', vmin=vmin, vmax=vmax)
    plt.imshow(edges, aspect=0.5, cmap='Greys', alpha=0.6)
    plt.show()
    
# Add widgets with keyword arguments for `canny`
widgets.interactive(canny_demo, sigma=8)

widgets.interactive(canny_demo, sigma=4, low=5, high=10)

images = []
for i in range(16):
    images.append(filter.canny(horizon, sigma=i))
out = np.dstack(images)
result = np.mean(out, axis=2)

plt.figure(figsize=(18,6))
plt.subplot(121)
plt.imshow(result, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=0.4)
plt.subplot(122)
plt.imshow(mag, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.show()