(For more resources related to this topic, see here.)

>>> data = scipy.stats.randint.rvs(0.4,10,size=(10,2))
>>> triangulation = scipy.spatial.Delaunay(data)

Any Delaunay class has the basic search attributes such as points (to obtain the set of points in the triangulation), vertices (that offers the indices of vertices forming simplices in the triangulation), neighbors (for the indices of neighbor simplices for each simplex—with the convention that “-1” indicates no neighbor for simplices at the boundary).

More advanced attributes, for example convex_hull, indicate the indices of the vertices that form the convex hull of the given points. If we desire to search for the simplices that share a given vertex, we may do so with the vertex_to_simplex method. If, instead, we desire to locate the simplices that contain any given point in the space, we do so with the find_simplex method.

At this stage we would like to point out the intimate relationship between triangulations and Voronoi diagrams, and offer a simple coding exercise. Let us start by choosing first a random set of points, and obtaining the corresponding triangulation.

>>> locations=scipy.stats.randint.rvs(0,511,size=(2,8))
>>> triangulation=scipy.spatial.Delaunay(locations.T)

We may use the matplotlib.pyplot routine triplot to obtain a graphical representation of this triangulation. We first need to obtain the set of computed simplices. Delaunay offers us this set, but by means of the indices of the vertices instead of their coordinates. We thus need to map these indices to actual points before feeding the set of simplices to the triplot routine:

>>>assign_vertex = lambda index: triangulation.points[index]
>>>triangle_set = map(assign_vertex, triangulation.vertices)
>>>matplotlib.pyplot.triplot(locations[1], locations[0], 
... triangles=triangle_set, color='r')

We will now obtain the edge map of the Voronoi diagram in a similar fashion as we did before, and plot it below the triangulation (since the former needs to be with either a pcolormesh or imshow command).

Note how the triangulation and the corresponding Voronoi diagrams are dual of each other; each edge in the triangulation (red) is perpendicular with an edge in the Voronoi diagram (white). How should we use this observation to code an actual Voronoi diagram for a cloud of points? The actual Voronoi diagram is the set of vertices and edges that composes it, rather than a binary image containing an approximation to the edges as we have computed.

Let us finish this Article with two applications to scientific computing that use these techniques extensively, in combination with routines from other SciPy modules.

Structural model of oxides

In this example we will cover the extraction of the structural model of a molecule of a bronze-type Niobium oxide, from HAADF-STEM micrographs.

The following diagram shows HAADF-STEM micrograph of a bronze-type Niobium oxide (taken from http://www.microscopy.ethz.ch/BFDF-STEM.htm, courtesy of ETH Zurich):

For pedagogical purposes, we took the following approach to solving this problem:

1. Segmentation of the atoms by thresholding and morphological operations.

2. Connected component labeling to extract each single atom for posterior examination.

3. Computation of the centers of mass of each label identified as an atom. This presents us with a lattice of points in the plane that shows a first insight in the structural model of the oxide.

4. Computation of the Voronoi diagram of the previous lattice of points. The combination of information with the output of the previous step will lead us to a decent (approximation of the actual) structural model of our sample.

Let us proceed in this direction.

Once retrieved, our HAADF-STEM images will be stored as big matrices with float32 precision. For this project, it is enough to retrieve some tools from the scipy.ndimage module, and some procedures from the matplotlib library. The preamble then looks like the following code:

import numpy
import scipy
from scipy.ndimage import *
from scipy.misc import imfilter
import matplotlib.pyplot as plt

The image is loaded with the imread(filename) command. This stores the image as a numpy.array with dtype = float32. Notice that the maxima and minima are 1.0 and 0.0, respectively. Other interesting information about the image can be retrieved:

img=imread('/Users/blanco/Desktop/NbW-STEM.png')
print "Image dtype: %s"%(img.dtype)
print "Image size: %6d"%(img.size)
print "Image shape: %3dx%3d"%(img.shape[0],img.shape[1])
print "Max value %1.2f at pixel %6d"%(img.max(),img.argmax())
print "Min value %1.2f at pixel %6d"%(img.min(),img.argmin())
print "Variance: %1.5fnStandard deviation:
%1.5f"%(img.var(),img.std())

This outputs the following information:

Image dtype: float32
Image size: 87025
Image shape: 295x295
Max value 1.00 at pixel 75440
Min value 0.00 at pixel 5703
Variance: 0.02580
Standard deviation: 0.16062

We perform thresholding by imposing an inequality in the array holding the data. The output is a Boolean array where True (white) indicates that the inequality is fulfilled, and False (black) otherwise. We may perform at this point several thresholding operations and visualize them to obtain the best threshold for segmentation purposes. The following images show several examples (different thresholdings applied to the oxide image):

By visual inspection of several different thresholds, we choose 0.62 as one that gives us a good map showing what we need for segmentation. We need to get rid of “outliers”, though; small particles that might fulfill the given threshold but are small enough not to be considered as an actual atom. Therefore, in the next step we perform a morphological operation of opening to get rid of those small particles. We decided that anything smaller than a square of size 2 x 2 is to be eliminated from the output of thresholding:

BWatoms = (img> 0.62)
BWatoms = binary_opening(BWatoms,structure=numpy.ones((2,2)))

We are ready for segmentation, which will be performed with the label routine from the scipy.ndimage module. It collects one slice per segmented atom, and offers the number of slices computed. We need to indicate the connectivity type. For example, in the following toy example, do we want to consider that situation as two atoms or one atom?

It depends; we would rather have it now as two different connected components, but for some other applications we might consider that they are one. The way we indicate the connectivity to the label routine is by means of a structuring element that defines feature connections. For example, if our criterion for connectivity between two pixels is that they are in adjacent edges, and then the structuring element looks like the image shown on the left-hand side from the images shown next. If our criterion for connectivity between two pixels is that they are also allowed to share a corner, then the structuring element looks like the image on the right-hand side. For each pixel we impose the chosen structuring element and count the intersections; if there are no intersections, then the two pixels are not connected. Otherwise, they belong to the same connected component.

We need to make sure that atoms that are too close in a diagonal direction are counted as two, rather than one, so we chose the structuring element on the left. The script then reads as follows:

structuring_element = [[0,1,0],[1,1,1],[0,1,0]]
segmentation,segments = label(BWatoms,structuring_element)

The segmentation object contains a list of slices, each of them with a Boolean matrix containing each of the found atoms of the oxide. We may obtain for each slice a great deal of useful information. For example, the coordinates of the centers of mass of each atom can be retrieved with the following commands:

coords = center_of_mass(img, segmentation, range(1,segments+1))
xcoords = numpy.array([x[1] for x in coords])
ycoords = numpy.array([x[0] for x in coords])

Note that, because of the way matrices are stored in memory, there is a transposition of the x and y coordinates of the locations of the pixels. We need to take it into account.

Notice the overlap of the computed lattice of points over the original image (the left-hand side image from the two images shown next). We may obtain it with the following commands:

>>>plt.imshow(img); plt.gray(); plt.axis('off')
>>>plt.plot(xcoords,ycoords,'b.')

We have successfully found the centers of mass for most atoms, although there are still about a dozen regions where we are not too satisfied with the result. It is time to fine-tune by the simple method of changing the values of some variables; play with the threshold, with the structuring element, with different morphological operations, and so on. We can even add all the obtained information for a wide range of those variables, and filter out outliers. An example with optimized segmentation is shown, as follows (look at the right-hand side image):

For the purposes of this exposition, we are happy to keep it simple and continue working with the set of coordinates that we have already computed. We will be now offering an approximation to the lattice of the oxide, computed as the edge map of the Voronoi diagram of the lattice.

L1,L2 = distance_transform_edt(segmentation==0,
return_distances=False,
return_indices=True)
Voronoi = segmentation[L1,L2]
Voronoi_edges= imfilter(Voronoi,'find_edges')
Voronoi_edges=(Voronoi_edges>0)

Let us overlay the result of Voronoi_edges with the locations of the found atoms:

>>>plt.imshow(Voronoi_edges); plt.axis('off'); plt.gray()
>>>plt.plot(xcoords,ycoords,'r.',markersize=2.0)

This gives the following output, which represents the structural model we were searching for: