Applying the Watershed Transform to a Cuneiform Tablet

This post was entirely written using the IPython notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at 20160830_SegmentingAnImageWithWatershed.ipynb.

This week, I've been watching the Scipy 2016 tutorial for Scikit-image on YouTube. This tutorial has prompted me to write a little application of the technique of segmentation using the watershed algorithms. This post is essentially based on the material found on Github here and the official documentation for scikit-image.

Goal of this post

The goal of this post is to clearly explain the different steps used to provide a useful segmentation of an image using the watershed algorithm.

To this end, we're going to use an image from wikipedia dealing with one of my current interests: cuneiform tablets.

tablet

As one sees in the image above, there might be several steps involved before getting a clean image with signs in it:

  • isolate the tablet from the background
  • pick out individual cuneiform characters from the tablet
  • threshold these to produce a final image

Let's get started!

Loading the image

To work on the image, we first need to retrieve it as an array of pixels. Luckily for us, scikit-image allows you to download an image if you hand an url to the imread function:

In [1]:
from skimage import io
In [2]:
img = io.imread("https://upload.wikimedia.org/wikipedia/commons/thumb/f/f2/Tablet_with_Cuneiform_Inscription_LACMA_M.79.106.2_%284_of_4%29.jpg/800px-Tablet_with_Cuneiform_Inscription_LACMA_M.79.106.2_%284_of_4%29.jpg")

Let's check the properties of this image:

In [3]:
img.shape
Out[3]:
(531, 800, 3)

It's a 531 by 800 pixels image with a rgb tuple defining the color of each pixel.

In [4]:
img.dtype
Out[4]:
dtype('uint8')

The colors are defined by integers (values between 0 and 255), which we will convert to floats while working on them.

In [5]:
from skimage import img_as_float
img = img_as_float(img)

Segmenting the background

To cut out the background, we're going to apply the watershed algorithm while initializing the seeds for the algorithm with some foreground and some background pixels.

This works in the following way:

  • create the seeds using a contour filter of the tablet image
  • perform the segmentation

Creating the seeds

We first create a gray image, from which we will compute edges.

In [6]:
from skimage import color
gray = color.rgb2gray(img)

Let's display the gray image:

In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] ='nearest'
plt.style.use('seaborn-muted')
In [8]:
plt.imshow(gray)
plt.colorbar()
Out[8]:
<matplotlib.colorbar.Colorbar at 0x8667d68>

Now, we filter it using an edge filter.

In [9]:
from skimage import filters
elevation_map = filters.canny(gray)
C:\Anaconda3\lib\site-packages\skimage\filters\__init__.py:28: skimage_deprecation: Call to deprecated function ``canny``. Use ``skimage.feature.canny`` instead.
  def canny(*args, **kwargs):

Let's display this:

In [10]:
plt.imshow(elevation_map)
plt.colorbar()
Out[10]:
<matplotlib.colorbar.Colorbar at 0xbc48828>

The canny filter found some edges inside and outside the tablet. To really make the interior stick out, we average it using a gaussian filter:

In [11]:
filled_elevation_map = filters.gaussian(elevation_map, sigma=5)
plt.imshow(filled_elevation_map)
plt.colorbar()
Out[11]:
<matplotlib.colorbar.Colorbar at 0xc051a90>

We now initialize the seeds by thresholding the previous image:

In [12]:
threshold = filters.threshold_otsu(filled_elevation_map)
threshold
Out[12]:
0.13830405859398809

This allows us to make seeds for the background (1) and the tablet (2):

In [13]:
import numpy as np
seeds = np.zeros_like(gray)
seeds[filled_elevation_map < threshold / 5] = 1
seeds[filled_elevation_map > threshold] = 2

This yields three types of pixels:

  • 0: to be assigned by the watershed algorihtm
  • 1: background
  • 2: foreground
In [14]:
plt.imshow(seeds)
plt.colorbar()
Out[14]:
<matplotlib.colorbar.Colorbar at 0xc2bcb38>

Performing the watershed

Let's now flood the image:

In [15]:
from skimage.morphology import watershed
segmentation = watershed(filled_elevation_map, seeds)

What does the result look like?

In [16]:
plt.imshow(segmentation)
plt.colorbar()
Out[16]:
<matplotlib.colorbar.Colorbar at 0xcefd128>

Great! We can recognize our tablet!

Finally, we can extract the tablet itself from the contour:

In [17]:
mask = (segmentation == 2).reshape(tuple([*segmentation.shape, 1]))
tablet = mask * img
plt.imshow(tablet)
Out[17]:
<matplotlib.image.AxesImage at 0xdaceba8>

Cleaning the characters

Now that we have our tablet, we can start working on the characters itself. What assyriologist often do to work with tablets is to first produce a transcription such as this one:

transcription

The transcription can be helpful for reading a tablet when it is not available anymore or quite simply to make a clearer picture of the signs, as tablets are difficult to read.

To obtain a sort of automated transcription, let's apply a local histogram equalization to the tablet, to remove the artifacts from the lighting.

In [18]:
from skimage.morphology import disk
equalized = color.rgb2gray(tablet)
selem = disk(30)
equalized = filters.rank.equalize(equalized, selem=selem)
plt.imshow(equalized)
plt.colorbar()
C:\Anaconda3\lib\site-packages\skimage\util\dtype.py:110: UserWarning: Possible precision loss when converting from float64 to uint8
  "%s to %s" % (dtypeobj_in, dtypeobj))
Out[18]:
<matplotlib.colorbar.Colorbar at 0xdb86a90>

Let's have a look at the resulting histogram:

In [19]:
plt.hist(equalized.ravel(), bins=256);
plt.ylim(0, 2000);

Ideally, we would like to have some points at the bottom of the tablet and some at the surface and do the watershed on them.

In [20]:
seeds = np.zeros_like(equalized)
seeds[equalized < 10] = 1
seeds[(equalized > 140) & (equalized < 200)] = 2
seeds[(~mask[:, :, 0]).nonzero()] = 3
In [21]:
plt.imshow(seeds)
plt.colorbar()
Out[21]:
<matplotlib.colorbar.Colorbar at 0xeba2278>
In [22]:
edges = filters.sobel(equalized)
In [23]:
ws = watershed(edges, seeds)
In [24]:
plt.imshow(ws)
Out[24]:
<matplotlib.image.AxesImage at 0xef5a240>

Finally, we can produce our cleaned image and compare it to the original:

In [25]:
import scipy.ndimage as ndi
cleaned = np.ones_like(ws)
# signs
cleaned[(ws == 1)] = 0
cleaned = ndi.binary_opening(cleaned)
# border of tablet
cleaned[filters.canny(mask[:, :, 0]) == 1] = 0
C:\Anaconda3\lib\site-packages\skimage\filters\__init__.py:28: skimage_deprecation: Call to deprecated function ``canny``. Use ``skimage.feature.canny`` instead.
  def canny(*args, **kwargs):

Let's compare both images in full size:

In [26]:
fig, ax = plt.subplots(1, 2, figsize=(16, 10))
ax[0].imshow(tablet)
ax[1].imshow(cleaned)
plt.tight_layout()

We can also apply some zoom on this:

In [27]:
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

zoom = lambda ax: ax.set_xlim(300, 500) and ax.set_ylim(300, 100)
ax[0].imshow(tablet)
zoom(ax[0])
ax[1].imshow(cleaned)
zoom(ax[1])
plt.tight_layout()