This Notebook accompanies an Agile Geoscience blog post on the same subject.
A quick look at a progression of ways to parameterize a seismic attribute. We'll use the Canny filter, since I was playing with it anyway (see Sobel_filtering_horizons.ipynb).
Some ways to parameterize, in roughly increasing order of sophistication:
We could use all of the above approaches in a Notebook too — but why would we? There are better ways.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Next we need some data.
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 saved this horizon as a NumPy object in the Notebook Sobel_filtering_horizons.ipynb.
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)
min 891.804814339, max 1065.95427275
We will compare our new attribute with this one later.
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()
The nicest way to figure out how it works is to look at the example from SciKits website. Copy and paste and it runs perfectly...
# 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()
Now we can get to work. Let's apply some different ways to investigate the parameters of skimage.filter.canny
.
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()
This is essentially the same as stepping over values and comparing side-by-side in PowerPoint, but about 100 times faster.
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()
I might decide I think something about halfway between those last two is about right. So as a very crude proxy for 'contrast' of the result, let's compute the average pixel, and aim for something in between.
pixels = float(horizon.size)
print np.sum(edge_8)/pixels, np.sum(edge_16)/pixels
0.0321759805434 0.0183676062218
I'm going to arbitrarily say that I want an average pixel of 0.028, so we'll compute increasing sigmas and stop when it falls below that value.
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()
One of the great things about computers is the interactivity. So why not replace the clumsy parameter setting in little bits of code with an interactive slider widget? We can also do cool things like overlay it on a map of the horizon, as a sort of quality check.
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)
I love that the slider 'incorrectly' lets us choose negative numbers, thereby letting us explore parts of the parameter space we might not even have bothered thinking of. (In this case, this is quite proper — negative sigmas make no sense and do not compute properly). But it's the thought that counts.
Even cooler is that this is the first method we've tried that makes it reasonable to explore multiple parameters...
widgets.interactive(canny_demo, sigma=4, low=5, high=10)
OK, I have explored my data, but now I'm confused, where can I go from here?
I could try to extend the 'goal seek' option, and come up with a more sophisticated cost function. If I could define something well enough — for edge detection, like coherence, I might be interested in contrast — then I could potentially just find the best answers, in much the same way that a digital camera autofocuses (indeed, many of them look for the highest contrast image). But goal seeking, if the cost function is too literal, in a way begs the question. I mean, you sort of have to know the answer — or something about the answer — before you find it.
Perhaps I can turn to other humans, in my social and professional networks. I could...
What if I could combine the best of all these approaches? Interactive exploration, with guided optimization, constrained by some cost function or other expectation. That could be interesting.
Unfortunately, I have absolutely no idea how that would work. But I think the optimization of the future will contain all of these elements.
Oh well, we can think of lots of easy things to do in code... things that might be hard to do in desktop software.
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()