I want to investigate the assertion I read recently that doing Sobel filtering on a coherency slice is a good idea. I highly doubt it, since we are interested in edges, not in the edges of edges.
We will use the open-access Penobscot dataset from offshore Nova Scotia, quite near my house. This is a great dataset for lots of reasons, but especially because it is freely available and can be shared in almost any way you like. You can download it from the Open Seismic Repository. The data is owned by the Nova Scotia Deptartment of Energy and licensed CC-BY. I used OpendTect software to load the data, then export one of the horizons that comes with the data, a shallow horizon called Horizon B. I exported the whole thing as an inline,crossline (as opposed to x,y) ASCII file.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
with open('Penobscot_HorB.txt') as f:
data = f.read()
print data[:100]
1003 1010 926.391959190369 1003 1011 925.13644695282 1003 1012 924.199283123016 1003 1013 923.783898
The square bracket notation is called a 'slice'. The data object is a string. The slice gives us the first 100 characters of it. Slicing, and the closely related indexing action, are very useful — and you'll see them a lot in Python code.
As expected, we have inline, xline, and two-way time. There is no header. We only really need the time values, so let's use slicing and indexing some more
We can actually make our live a bit easier by re-reading the file into a list, one line at a time. This is easy:
data = np.loadtxt('Penobscot_HorB.txt')
data[:5]
array([[ 1003. , 1010. , 926.39195919], [ 1003. , 1011. , 925.13644695], [ 1003. , 1012. , 924.19928312], [ 1003. , 1013. , 923.78389835], [ 1003. , 1014. , 923.74902964]])
inlines = 1597 - 1003 + 1
xlines = 1471 - 1009 + 1
print inlines, 'inlines'
print xlines, 'crosslines'
595 inlines 463 crosslines
data[:,0] = inlines - (data[:,0] - 1002)
data[:,1] -= 1008 # same as data[:,1] = data[:,1] - 1008
data[:5]
array([[ 594. , 2. , 926.39195919], [ 594. , 3. , 925.13644695], [ 594. , 4. , 924.19928312], [ 594. , 5. , 923.78389835], [ 594. , 6. , 923.74902964]])
horizon = np.zeros((inlines, xlines))
print horizon
print horizon.shape
[[ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] ..., [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.]] (595, 463)
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
array([[ 0. , 0. , 0. , ..., 1088.39404583, 1088.08231354, 1087.8816843 ], [ 0. , 0. , 0. , ..., 1088.29820156, 1087.96453476, 1087.78274059], [ 0. , 0. , 0. , ..., 1088.17863464, 1087.8431797 , 1087.68224716], ..., [ 0. , 926.46729946, 925.56166649, ..., 972.23240137, 972.9629159 , 974.19649363], [ 0. , 926.39195919, 925.13644695, ..., 973.05631638, 973.32763672, 974.58219528], [ 0. , 0. , 0. , ..., 1088.48810196, 1088.18411827, 1087.98325062]])
Aside: save this out as a file to use in other examples.
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)
min 891.804814339, max 1065.95427275
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()
Let's try some classic edge detecting filters.
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()
Let's try it on some real data.
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()
Let's use the ndimage.sobel
implementation to compute the full (non-oriented) result:
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()
First we'll turn the Sobel thing into a function, to make it easier to iterate.
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()
scipy.ndimage.generic_filter()
seems like it can be used as a non-linear filter, because you just provide a function to work on the 1D array that ndimage
passes.
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:]))
mean = 4.55555555556 conservative mean = 2.75 5.14481640124 middle index 4 middle value 3 rest of values [ 3 2 3 4 19 2 3 2]
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]
919.426202774
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))
1 loops, best of 3: 4 s per loop
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()
Example from SciKits website.
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)
8864 5060
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)
This is more elegant, and works in nbviewer.org because everything is precomputed...
from ipywidgets import StaticInteract, RangeWidget
Need a different function, slightly.
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))
Now that we have everything in code, we can also get creative... how about summing all sigmas up to 20? Let's try that, and compare to the Sobel edge map we created earlier:
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()
Warning This will turn all subsequent plots into mplD3 plots.
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()