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

with open('Penobscot_HorB.txt') as f:
    data = f.read()

print data[:100]

data = np.loadtxt('Penobscot_HorB.txt')
data[:5]

inlines = 1597 - 1003 + 1
xlines = 1471 - 1009 + 1

print inlines, 'inlines'
print xlines, 'crosslines'

data[:,0] = inlines - (data[:,0] - 1002) 
data[:,1] -= 1008 # same as data[:,1] = data[:,1] - 1008

data[:5]

horizon = np.zeros((inlines, xlines))
print horizon
print horizon.shape

for sample in data:
    inline = sample[0]
    xline = sample[1]
    z_value = sample[2]
    
    # We have to subtract 1 to allow for 0-based indexing
    horizon[inline - 1, xline - 1] = z_value

horizon

np.save('Penobscot_HorB.npy', horizon)

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)

ch = plt.get_cmap('cubehelix_r')

plt.imshow(horizon, aspect=0.5, cmap=ch, vmin=vmin, vmax=vmax)
plt.colorbar(shrink=0.75)  # shrink makes the colourbar a bit shorter
plt.show()

import scipy.signal

laplacian = np.array([[0, -1,  0],
                      [-1, 4, -1],
                      [0, -1,  0]])

robinson = np.array([[-1, 0, 1],
                    [-1, 0, 1],
                    [-1, 0, 1]])

sobel = np.array([[-1, 0, 1],
                  [-2, 0, 2],
                  [-1, 0, 1]])

plt.figure(figsize=(8,4))

plt.subplot(121)
plt.imshow(robinson, interpolation='nearest')
plt.set_cmap('jet')

plt.subplot(122)
plt.imshow(sobel, interpolation='nearest')
plt.set_cmap('jet')

plt.show()

sobel_out = scipy.signal.convolve2d(horizon, sobel)

plt.imshow(-np.power(sobel_out, 1), aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=-50, vmax=50)
plt.show()

import scipy.ndimage

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

plt.figure(figsize=(12,8))
plt.imshow(mag, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.show()

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)
mag2 = do_sobel(mag)
mag3 = do_sobel(mag2)
mag4 = do_sobel(mag3)
mag5 = do_sobel(mag4)
mag6 = do_sobel(mag5)

m = 3

plt.figure(figsize=(18,12))
plt.subplot(231)
plt.imshow(mag[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.subplot(232)
plt.imshow(mag2[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.subplot(233)
plt.imshow(mag3[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.subplot(234)
plt.imshow(mag4[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.subplot(235)
plt.imshow(mag5[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.subplot(236)
plt.imshow(mag6[200:400,200:300], aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=m)
plt.show()

plt.figure(figsize=(18,5))

plt.subplot(121)
plt.imshow(mag2, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.axvline(190)

plt.subplot(122)
plt.plot(np.clip(mag[300:400,190], 0, 3))
plt.plot(np.clip(mag2[300:400,190], 0, 3), color='b', alpha=0.6)
plt.plot(np.clip(mag3[300:400,190], 0, 3), color='b', alpha=0.4)
plt.plot(np.clip(mag4[300:400,190], 0, 3), color='b', alpha=0.3)
plt.plot(np.clip(mag5[300:400,190], 0, 3), color='b', alpha=0.2)
plt.show()

plt.figure(figsize=(18,5))

plt.subplot(121)
plt.imshow(mag2, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.axvline(190)

plt.subplot(122)
plt.plot(np.clip(mag[300:400,190], 0, 3))
plt.plot(np.clip(mag6[300:400,190], 0, 3), color='lightblue')
plt.show()

b = np.array([3,2,3,4,3,19,2,3,2]) 
print "mean =", np.mean(b)
print "conservative mean =", np.mean([i for i in b if i < 2 * np.mean(b)])
print np.std(b)

m = b.size/2
print "middle index", m
print "middle value", b[m]
print "rest of values", np.concatenate((b[:m], b[m+1:]))

def conservative(a):
    m = a.size/2
    c = a[m] # middle value
    b = np.concatenate((a[:m], a[m+1:]))

    return min(np.amax(b), max(np.amin(b), c))

horizon_sp = np.copy(horizon)
horizon_sp[200,200] = 5000

cons = scipy.ndimage.generic_filter(horizon_sp, conservative, size=(5,5))
print cons[200,200]

plt.figure(figsize=(12,6))
plt.subplot(121)
plt.imshow(horizon_sp[180:220,180:220], cmap=ch, interpolation="nearest", vmin=vmin, vmax=vmax)
plt.subplot(122)
plt.imshow(cons[180:220,180:220], cmap=ch, interpolation="nearest", vmin=vmin, vmax=vmax)
plt.show()

%timeit scipy.ndimage.generic_filter(horizon, conservative, size=(5,5))

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()

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=(12,8))
plt.imshow(edge, aspect=0.5, cmap=plt.get_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.subplot(222)
plt.imshow(edge_4, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.subplot(223)
plt.imshow(edge_8, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.subplot(224)
plt.imshow(edge_16, aspect=0.5, cmap=plt.get_cmap('Greys'), vmin=0, vmax=1)
plt.show()

print np.sum(edge_8), np.sum(edge_16)

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

plt.figure(figsize=(12,8))
plt.imshow(result, aspect=0.5, cmap=plt.get_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=ch, 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_threshold=5, high_threshold=10)

from ipywidgets import StaticInteract, RangeWidget

from IPython.html import widgets

image = np.copy(horizon)

def canny_demo(**kwargs):
    edges = filter.canny(image, **kwargs)
    fig = plt.figure(figsize=(9,6))
    plt.imshow(horizon, aspect=0.5, cmap=ch, vmin=vmin, vmax=vmax)
    plt.imshow(edges, aspect=0.5, cmap='Greys', alpha=0.6)
    return fig

StaticInteract(canny_demo, sigma=RangeWidget(1,8,0.5))

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()

import mpld3
mpld3.enable_notebook()

fig, ax = plt.subplots(subplot_kw=dict(axisbg='#EEEEEE'))
ax.grid(color='white', linestyle='solid')

N = 50
scatter = ax.scatter(np.random.normal(size=N),
                     np.random.normal(size=N),
                     c=np.random.random(size=N),
                     s = 1000 * np.random.random(size=N),
                     alpha=0.3,
                     cmap=plt.cm.jet)

ax.set_title("D3 Scatter Plot", size=18);

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