Computer vision in geoscience: recover seismic data from images – part 2

Posted on September 18, 2017 by matteomycarta

In part 1 of this short series I demonstrated how to detect the portion occupied by the seismic section in an image (Figure 1).

Figure 1

The result was a single binary image with the white object representing the pixels occupied by the seismic section (Figure 2).

Figure 2

You can download from GitHub all the tools for the automated workflow (including both part 1 and part 2, and some of the optional features outlined in the introduction) in the module mycarta.py, as well as an example Jupyter Notebook showing how to run it.

Next I want to use this binary object to derive a transformation function to rectify to a rectangle the seismic section in the input image.

The first step is to detect the contour of the object. Notice that because we used morphological operations it is not a perfect quadrilateral: it has rounded corners and some of the sides are bent, therefore the second step will be to approximate the contour with a polygon with enough tolerance to ensure it has 4 sides only(this took some trial and error but 25 turned out to be a good value for the parameter for a whole lot of test images I tried).

In reality, the two steps are performed together using the functions find_contours (there is only one to find, really) and approximate_polygon from the skimage.measure module, as below:

contour = np.squeeze(find_contours(enhanced, 0))
coords = approximate_polygon(contour, tolerance=25)

The variable coords contains the coordinate for the corner points of the polygon (the first point is repeated last to close the polygon), which in Figure 3 I plotted superimposed to the input binary object.

Figure 3 – approximated polygon

A problem with the output of approximate_polygon is that the points are not ordered; to solve it I adapted a function from a Stack Overflow answer to sort them based on the angle from their centroid:

def ordered(points):
  x = points[:,0]
  y = points[:,1]
  cx = np.mean(x)
  cy = np.mean(y)
  a = np.arctan2(y - cy, x - cx)
  order = a.ravel().argsort()
  x = x[order]
  y = y[order]
  return np.vstack([x,y])

I call the function as below to get the corners in the contour without the last one (repetition of the first point).

sortedCoords = ordered(coords[:-1]).T

I can then plot them using colors in a predefined order to convince myself the indeed are sorted:

plt.scatter(sortedCoords[:, 1], sortedCoords[:, 0], s=60, 
 color=['magenta', 'cyan', 'orange', 'green'])

Figure 4 – corners sorted in counter-clockwise order

The next bit of code may seem a bit complicated but it is not. With coordinates of the corners known, and their order as well, I can calculate the largest width and height of the input seismic section, and I use them to define the size of the registered output section, which is to be of rectangular shape:

w1 = np.sqrt(((sortedCoords[0, 1]-sortedCoords[3, 1])**2)
  +((sortedCoords[0, 0]-sortedCoords[3, 0])**2))
w2 = np.sqrt(((sortedCoords[1, 1]-sortedCoords[2, 1])**2)
  +((sortedCoords[1, 0]-sortedCoords[2, 0])**2))

h1 = np.sqrt(((sortedCoords[0, 1]-sortedCoords[1, 1])**2)
  +((sortedCoords[0, 0]-sortedCoords[1, 0])**2))
h2 = np.sqrt(((sortedCoords[3, 1]-sortedCoords[2, 1])**2)
  +((sortedCoords[3, 0]-sortedCoords[2, 0])**2))

w = max(int(w1), int(w2))
h = max(int(h1), int(h2))

and with those I define the coordinates of the output corners used to derive the transformation function:

dst = np.array([
  [0, 0],
  [h-1, 0],
  [h-1, w-1],
  [0, w-1]], dtype = 'float32')

Now I have everything I need to rectify the seismic section in the input image: it is warped using homologous points (the to sets of four corners) and a transformation function.

dst[:,[0,1]] = dst[:,[1,0]]
sortedCoords[:,[0,1]] = sortedCoords[:,[1,0]]
tform = skimage.transform.ProjectiveTransform()
tform.estimate(dst,sortedCoords)
warped =skimage.transform.warp(img, tform, output_shape=(h-1, w-1))

Notice that I had to swap the x and y coordinates to calculate the transformation function. The result is shown in Figure 5: et voilà!

Figure 5 – rectified seismic section

You can download from GitHub the code to try this yourself (both part 1 and part 2, and some of the optional features outlined in the introduction, like removing the rectangle with label inside the section) as well as an example Jupyter Notebook showing how to run it.

Computer vision in geoscience: recover seismic data from images – part 1

Posted on September 10, 2017 by matteomycarta

As anticipated in the introductory post of this short series I am going to demonstrate how to automatically detect where a seismic section is located in an image (be it a picture taken from your wall, or a screen capture from a research paper), rectify any distortions that might be present, and remove all sorts of annotations and trivia around and inside the section.

In this part one I will be focusing on the image preparation and enhancement, and the automatic detection of the seismic section (all done using functions from numpy, scipy, and scikit-image). In order to do that, first I convert the input image (Figure 1) containing the seismic section to grayscale and then enhance it by increasing the image contrast (Figure 2).

Figure 1 – input image

Figure 2 – grayscale image

All it takes to do that is three lines of code as follows:

gry = skimage.color.rgb2gray(img);
p2, p95 = numpy.percentile(gry, (2, 95))
rescale = exposure.rescale_intensity(gry, in_range=(p2, p95))

For a good visual intuition of what actually is happening during the contrast stretching, check my post sketch2model – sketch image enhancements: in there I show intensity profiles taken across the same image before and after the process.

Finding the seismic section in this image involve four steps:

converting the grayscale image to binary with a threshold (in this example a global threshold with the Otsu method)
finding and retaining only the largest object in the binary image (heuristically assumed to be the seismic section)
filling its holes
applying morphological operations to remove minutiae (tick marks and labels)

Below I list the code, and show the results.

global_thresh = threshold_otsu(rescale)
binary_global = rescale < global_thresh

Figure 3 – binary image

# (i) label all white objects (the ones in the binary image).
# scipy.ndimage.label actually labels 0s (the background) as 0 and then
# every non-connected, nonzero object as 1, 2, ... n.
label_objects, nb_labels = scipy.ndimage.label(binary_global)

# (ii) calculate every labeled object's binary size (including that 
# of the background)
sizes = numpyp.bincount(label_objects.ravel())

# (3) set the size of the background to 0 so that if it happened to be 
# larger than the largest white object it would not matter
sizes[0] = 0

# (4) keep only the largest object
binary_objects = remove_small_objects(binary_global, max(sizes))

Figure 4 – isolated seismic section

# Remove holes (black regions inside white object)
binary_holes = scipy.ndimage.morphology.binary_fill_holes(binary_objects)

Figure 5 – holes removed

enhanced = opening(binary_holes, disk(7))

Figure 6 – removed residual tick marks and labels

That’s it!!!

In the next post, we will use this polygonal binary object both as a basis to capture the actual coloured seismic section from the input image and to derive a transformation to rectify it to a rectangle.

Computer vision in geoscience: recover seismic data from images, introduction

Posted on July 30, 2017 by matteomycarta

In a recent post titled Unweaving the rainbow, Matt Hall described our joint attempt (partly successful) to create a Python tool to enable recovery of digital data from any pseudo-colour scientific image (and a seismic section in particular, like the one in Figure 1), without any prior knowledge of the colormap.

Figure 1. Test image: a photo of a distorted seismic section on my wall.

Please check our GitHub repository for the code and slides and watch Matt’s talk (very insightful and very entertaining) from the 2017 Calgary Geoconvention below:

In the next two post, coming up shortly, I will describe in greater detail my contribution to the project, which focused on developing a computer vision pipeline to automatically detect where the seismic section is located in the image, rectify any distortions that might be present, and remove all sorts of annotations and trivia around and inside the section. The full workflow is included below (with sections I-VI developed to date):

I – Image preparation, enhancement:
1. Convert to gray scale
2. Optional: smooth or blur to remove high frequency noise
3. Enhance contrast

II – Find seismic section:
1. Convert to binary with adaptive or other threshold method
2. Find and retain only largest object in binary image
3. Fill its holes
4. Apply opening and dilation to remove minutiae (tick marks and labels)

III – Define rectification transformation
1. Detect contour of largest object find in (2). This should be the seismic section.
2. Approximate contour with polygon with enough tolerance to ensure it has 4 sides only
3. Sort polygon corners using angle from centroid
4. Define new rectangular image using length of largest long and largest short sides of initial contour
5. Estimate and output transformation to warp polygon to rectangle

IV – Warp using transformation

V – Blanking annotations inside seismic section (if rectangular):
1. Start with output of (4)
2. Pre-process and apply canny filter
3. Find contours in the canny filter smaller than input size
4. Sort contours (by shape and angular relationships or diagonal lengths)
5. Loop over contours:
  1. Approximate contour
  2. If approximation has 4 points AND the 4 semi-diagonals are of same length: fill contour and add to mask

VI – Use mask to remove text inside rectangle in the input and blank (NaN) the whole rectangle.

VII – Optional: tools to remove arrows and circles/ellipses:
1. For arrows – contours from (4) find ones with 7 sizes and low convexity (concave) or alternatively Harris corner and count 7 corners, or template matching
2. For ellipses – template matching or regionprops

VIII – Optional FFT filters to remove timing lines and vertical lines

You can download from GitHub all the tools for the automated workflow (parts I-VI) in the module mycarta.py, as well as an example Jupyter Notebook showing how to run it.

The first post focuses on the image pre-processing and enhancement, and the detection of the seismic line (sections I and II, in green); the second one deals with the rectification of the seismic (sections IV to V, in blue). They are not meant as full tutorials, rather as a pictorial road map to (partial) success, but key Python code snippets will be included and discussed.

An introduction to color palettes for seismic amplitude – teaser

Posted on April 20, 2013 by matteomycarta

Introduction

In a future posts I will take a look at some of the color palettes used for seismic amplitude display, and discuss ways we can design more perceptual and more efficient ones.

For now, I would like to ask readers to look at two sets of seismic images and answer the survey questions in each section. Far from being exhaustive sets, these are meant as a teaser to get a conversation started and exchange opinions and preferences.

Stratigraphic interpretation

The seismic line below is inline 424 from the F3 dataset, offshore Netherlands from the Open Seismic Repository (licensed CC-BY-SA).

I generated an animation, played at 0.5 frames/second, where 8 different color palette are alternated in sequence. Please click on the image to see a full resolution animation. I also generated a 0.25 frame/second version and a 1 frame/second version.

Fault interpretation

The images used to create the panel below are portions of seismic displays kindly provided by Steve Lynch of 3rd Science Solutions, generated using data released by PeruPetro. I am grateful to both.

Thanks to Matt Hall and Evan Bianco of Agile Geoscience for their suggestions.

MyCarta

A blog about Technology, Geoscience, Visualization, Data Science

Tag Archives: seismic colormap

Computer vision in geoscience: recover seismic data from images – part 2

Like this:

Computer vision in geoscience: recover seismic data from images – part 1

Like this:

Computer vision in geoscience: recover seismic data from images, introduction

Like this:

An introduction to color palettes for seismic amplitude – teaser

Introduction

Stratigraphic interpretation

Fault interpretation

Like this:

Share this:

Like this:

Share this:

Like this:

Share this:

Like this:

Introduction

Stratigraphic interpretation

Fault interpretation

Share this:

Like this: