SciPy for Computational Geometry

For solving complex problems in mathematics, science, or engineering, SciPy is the solution. Building on your basic knowledge of Python, and using a wealth of examples from many scientific fields, this book is your expert tutor.


>>> 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, courtesy of ETH Zurich):


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



  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:

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.5f\nStandard deviation:

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')

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,
Voronoi = segmentation[L1,L2]
Voronoi_edges= imfilter(Voronoi,'find_edges')

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

>>>plt.imshow(Voronoi_edges); plt.axis('off'); plt.gray()

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

A finite element solver for Poisson's equation

We use finite elements when the size of the data is so large that it results prohibitive to deal with finite differences. To illustrate this case, we would like to explore the potential flow over a wing, as a solution to the Laplace equation subjects to certain boundary conditions.

We wish to create a simple profile of a wing, and produce a mesh surrounding it. This will be our starting point to solve this problem using finite elements, as we will be placing on the domain a piecewise continuous function, whose pieces are linear and supported on each of the triangles.

import numpy
from numpy import pi, cos, sin, hstack, vstack, linspace, where
from numpy import ones, multiply, cross, array, mat, zeros, mgrid
import scipy
import matplotlib.pyplot as plt
from scipy.special import exp10
from scipy.linalg import norm
from scipy.sparse import dok_matrix
from scipy.sparse.linalg import spsolve
from scipy.interpolate import LinearNDInterpolator
from scipy.spatial import Delaunay

We will be using two functions to generate vertices of our triangulation:

paramtr=lambda s:linspace(0,1,s)
ellipse=lambda a,b,s:[a*cos(2*pi*paramtr(s)), b*sin(2*pi*paramtr(s))]

We will start with a grid of a sufficiently large domain where the wing profile is to be included. We will complement this basic grid with enough points on the wing profile, which is designed as an ellipse:

for k in range(16):

We will be restricting the domain to a small rectangular region. We wish to introduce enough points in that border:


Let us now perform the restriction of vertices, as follows:

inside_vertices=where( multiply(abs(vertices[0])<=200,
abs(vertices[1])<=100 ))

We may create now the triangulation, and erase from it all triangles that are inside of the wing profile, and outside the rectangle [-200,200]x[-100,100]. We do so by computing the center of mass for each triangle, and discarding those triangles whose centers are inside of the ellipse, or outside the rectangle:

triangulation = Delaunay(vertices.T)
index2point = lambda index: triangulation.points[index]
all_centers = index2point(triangulation.vertices).mean(axis=1)
not_in_wing = lambda pt: (pt[0]/128)**2+(pt[1]/16)**2>=1

We then have the following triangulation:


This produces the following graph:

In this case, the flow potential is the solution of the Laplace equation, with boundary conditions as follows:

Here, Cin is the set of vertical edges on the leftmost side of the rectangle. Cout is the set of vertical edges on the rightmost side of the rectangle. We code the solution in the usual fashion. We compute the stiff matrix A (which for obvious reasons need to be sparse), the matrix R and the vector r holding the Robin conditions. With them, the solution to the system comes from the solution X of the system (A + R) X = r. This should be no trouble for SciPy. Let us start with the stiff matrix:

for triangle in triangulation.vertices:

Note the cumbersome way to update the matrix stiff_matrix. This is due to the fact that the matrix is sparse, and the current choice of representation does not behave well with indexing.

To compute the Robin matrix and vector we need to collect all edges on the boundary first. We also need to define the kappa and gN functions to help us design the boundary conditions:

kappa=lambda pt: exp10(6)*(pt[0]>99.99)
gN=lambda pt:float(pt[0]<=99.99)
for edge in triangulation.convex_hull:
helper_matrix[edge,edge]= length*kappa(center)*array([2,1,1,2])

We are ready to solve the equation, precisely by computing the linear interpolant on the vertices of the triangulation, with the values obtained in our previous step:

>>> solution=LinearNDInterpolator(triangulation.points,sltn_v)


This article explored the construction of triangulation of points, convex hulls, Voronoi diagrams, and some applications.

Resources for Article :

Further resources on this subject:

Books to Consider

comments powered by Disqus