# Sobel filtering horizons¶

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.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
with open('Penobscot_HorB.txt') as f:

In [3]:
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:

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

Out[4]:
array([[ 1003.        ,  1010.        ,   926.39195919],
[ 1003.        ,  1011.        ,   925.13644695],
[ 1003.        ,  1012.        ,   924.19928312],
[ 1003.        ,  1013.        ,   923.78389835],
[ 1003.        ,  1014.        ,   923.74902964]])
In [5]:
inlines = 1597 - 1003 + 1
xlines = 1471 - 1009 + 1

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

595 inlines
463 crosslines

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

data[:5]

Out[6]:
array([[ 594.        ,    2.        ,  926.39195919],
[ 594.        ,    3.        ,  925.13644695],
[ 594.        ,    4.        ,  924.19928312],
[ 594.        ,    5.        ,  923.78389835],
[ 594.        ,    6.        ,  923.74902964]])
In [7]:
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)

In [8]:
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

Out[8]:
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.

In [9]:
np.save('Penobscot_HorB.npy', horizon)

In [10]:
plt.figure()
plt.imshow(-1 * horizon, aspect=0.5)
plt.show()

In [11]:
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

In [12]:
ch = plt.get_cmap('cubehelix_r')

In [13]:
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()


# Edge detection by conovlution¶

Let's try some classic edge detecting filters.

In [14]:
import scipy.signal

In [15]:
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]])

In [16]:
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.

In [17]:
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:

In [18]:
import scipy.ndimage

In [19]:
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

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


## Iterative Sobel filtering¶

First we'll turn the Sobel thing into a function, to make it easier to iterate.

In [21]:
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

In [22]:
mag = do_sobel(horizon)
mag2 = do_sobel(mag)
mag3 = do_sobel(mag2)
mag4 = do_sobel(mag3)
mag5 = do_sobel(mag4)
mag6 = do_sobel(mag5)

In [23]:
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()