The problem

Let’s consider a polymer sheet, 50 µm thick, where very small holes (~ 50 nm) are randomly distributed. When these holes are filled, this system becomes an array of nanowires embedded in a polymer template.

If the holes are filled with a magnetic materiel (e.g. Cobalt, Iron) the array can be analyzed by ferromagnetic resonance (FMR). The experiment consists in injecting a high frequency (1 GHz – 65 GHz) microwave into the nanowires, while sweeping an external magnetic field and measuring the absorption of the microwave.

For the sake of simplicity, just recall that the FMR analysis yields an absoption curve that follows a gaussian distribution with a maximum occuring at the resonance field.

According to theory for a single, isolated nanowire the absorption curve should be perfectly gaussian with a linewidth about 150-300 Oe (depending on the type of metal, the shape of the wire, etc.). For a system perfectly distributed and wires isolated, the total absorption curve should be roughly equivalent to that of a single nanowire.

However, in the real world, holes are not always perfectly distributed and wires have neighbors. This results in interactions, leading to an enlargement of the absorption curves mainly due to varying dipolar interactions between wires (basically the magnetic force responsible for the attractive/repulsive force).

The objective

The goal here is to analyse pictures taken with an electron microscope, perform an analysis of distribution and input the physical parameters to try and reproduce the experimental measurements.

Image analysis

To setup the script, we’ll start with a small picture of a template with few holes rather regularly distributed. Note from the scale on the picture below that the average diameter is about 50 nm. Some holes are not nice regular circles, this is due to complex imaging conditions (We’re using an electron microscope on an insulating surface: electrons accumulate on the surface, leading to distortions).

The process flow is rather straightforward :

Open file and crop to the region of interest
Adjust contrast/exposition to get a picture as smooth as possible
Identify shapes and segment picture
Get the geometrical parameters : mean diameter, distance between holes
Perform Delaunay triangulation

A top-view of a nanoporous template showing holes of ~ 50 nm diameter, distributed as a honeycomb.

First, we’ll crop the image to remove the text information at its bottom.

from skimage import io
figName = "template2.jpg"
img_init = io.imread(figName, as_grey=True)
img_crop = img_init[0:302,0:432]
io.imshow(img_crop)

Image cropped

Now, we need to improve the quality of the image. The constrast should be smoothed to help identify all holes, in particular the blurry ones.

#Import modules
from scipy.ndimage import gaussian_filter
from skimage import img_as_float
from skimage.morphology import reconstruction

#Apply a gaussian filter to smooth things
image = img_as_float(img_crop)
image = gaussian_filter(image,0.9)

#Smooth contrast by subtracting a mask
seed = np.copy(image)
seed[1:-1, 1:-1] = image.min()
mask = image
dilated = reconstruction(seed-0.15, mask, method='dilation')

The result is satisfying : the area between the holes is a smooth grey with no variation in constrast.

We have fixed the variations of contrast : black holes on a grey, uniform background.

Now, we need to identify the holes. First, I’ll use the Sobel transform to find the hole’s edges then use a watershed transform to make sure all holes are properly separated.

We’ll use markers to separate the two regions : holes and inter-hole : these regions are separated by intensity value manually chosen from the numpy array representing the image.

from skimage.filters import sobel
from skimage.morphology import watershed

elevation = sobel(image)
markers = np.zeros_like(img_crop)
markers[img_crop < 70] = 1
markers[img_crop > 120] = 2
segmentation = watershed(elevation,markers)
io.imshow(segmentation)

The result is pretty nice : only one hole is missing at the bottom right, and some fine tuning could help but let’s keep it as is.

Before performing the pattern recognition and Delaunay triangulation, let’s remove (or rather, ignore) all holes at the borders. These act as outliers diameter or distances calculations.

Unwanted objects at the borders are removed. We can proceed to further analysis.

OK, that’s perfect now. We have a nice netword of isolated circles on a uniform background.

Let’s proceed to pattern recognition using scikit’s image measure.label and measure.regionprops. Also, to check if this worked, I’ll superimpose red patches around each hole.

from skimage import measure
import matplotlib.patches as mpatches

img_label = measure.label(segmentation) 
img_regions = measure.regionprops(img_label)#Identified regions
image_label_overlay = color.label2rgb(img_label, image=image)
fig, ax = plt.subplots(figsize=(10, 6))
ax.imshow(image_label_overlay)
for region in img_regions:
    # take regions with large enough areas
    if region.area >= 10:
        # draw rectangle around segmented coins
        minr, minc, maxr, maxc = region.bbox
        rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                  fill=False, edgecolor='red', linewidth=2)
        ax.add_patch(rect)

ax.set_axis_off()
plt.tight_layout()
plt.show()

All holes have been identified

Let’s perform the triangulation. We will need :

The centroids, i.e. the coordinates of each hole.
The diameter of each hole
The porosity, i.e. in each cell (or triangle), the ratio of the area occupied by a hole.

from scipy.spatial ConvexHull
maxdim = len(img_regions)-1 #Number of objects identified
centroids = np.zeros((maxdim,2))
diameters = np.zeros((maxdim))
porosite = np.zeros((nb_triangles))

Now, we get the coordinates of all centroids and perform the triangulation on these elements. The diameter is measured with the “equivalent_diameter” property in skimage.measure.regionprops that determines the diameter of a circle fitting into the square surrounding the item.

#We can use built-in functions to get coordinates and diameter of each holes
for i in range(0,len(img_regions)-1):
    centroids[i] = np.array(img_regions[i+1].centroid) 
    diameters[i] = np.array(img_regions[i+1].equivalent_diameter)

#Triangulation
from scipy.spatial import Delaunay
triangles = Delaunay(centroids) #This is a qhull.Delaunay object
#Plot
plt.triplot(centroids[:,1], centroids[:,0], triangles.simplices.copy())

Here is a nice triangulation

The triangulation stores a lot of information: each hole and triangle are indexed their coordinates are stored. Then, we have attributes and methods that allow for useful operations, for example:

triangle.vertices returns an array containing the indices of each triangle’s points
triangle.points returns an array with the coordinates of each hole (i.e. each identified circle).
triangle.neighbors returns an array with the indices of neighboring triangles

triangles.vertices
>>>array([[ 66,  45,  82],
>>>       [ 66,  48,  45],
>>>       [ 48,  66,  67],
>>>       ..., 
>>>       [207, 224, 204],
>>>       [204, 224, 261],
>>>       [224, 243, 261]], dtype=int32)

triangles.points
>>>array([[   8.63492063,   77.06349206],
>>>       [   8.60305344,  243.76335878],
>>>       [   7.84375   ,  265.6015625 ],
>>>       ..., 
>>>       [ 292.30952381,   89.92857143],
>>>       [ 293.78333333,  112.44166667],
>>>       [ 294.15833333,  284.2       ]])

triangles.neighbors
>>>array([[ -1,  35,   1],
>>>       [ 38,   0,   2],
>>>       [ 36,  42,   1],
>>>       ..., 
>>>       [528, 518, 526],
>>>       [529,  -1, 527],
>>>       [277, 528, 525]], dtype=int32)

To calculate areas, we could use the above functions to take the coordinates of each triangle, then compute distances, use Heron’s formula etc. Hopefully, we can use scipy’s ConvexHull module to get the same information.

Be careful here, there is a trick. ConvexHull is meant initially for 3D objects and the methods (area, volume) are named accordingly.

When used in 2D : volume will return the area and are will return the perimeter ! Be aware of this 😉

nb_triangles = triangles.nsimplex #nsimplex = number of triangles
for i in range(nb_triangles): #Area calculation
    # WARNING ConvexHull.volume yields a surface in our case
    porosite[i] = (mean_area_px/2)/(ConvexHull(triangles.points[triangles.vertices[i]]).volume)  
#From the scale we know that 1 px is 4.16 nm
diameters_nm = (diameters)*4.16
mean_diameter_px = np.mean(diameters) #Mean diameter
mean_diameter_nm = np.mean(diameters_nm)
mean_area_px = np.pi*(mean_diameter_px/2)**2
mean_area_nm = np.pi*(mean_diameter_nm/2)**2 #Conversion into nanometers
#Switch to a list just for the fun of using a list comprehension
porosite = list(porosite)
porosite = [_ for _ in porosite if _<1] #A porosity > 1 has no physical meaning
print("Mean porosity", np.mean(porosite))

Let’s proceed to some geometrical analysis. The mean diameter is 12.8 px, which corresponds to 53.2 nm.

The average porosity is about 30%, which is great since it is the expected value !

And now for some fun physics

The magnetic susceptibility is written as a non-symetric tensor :

$\bar{\bar\chi} = \begin{vmatrix} \chi & i \chi_a & 0 \\ -i \chi_a & \chi & 0 \\ 0 & 0 & 0 \end{vmatrix}$

The diagonal and off-diagonal terms are :

$\chi = \frac{\gamma M_0 \omega}{{\omega_H^2-\omega^2}}$ and $\chi_a = \frac{\gamma M_0 \omega}{{\omega_H^2-\omega^2}}$

With : $\chi = \chi^' - i\chi''$ and $\chi_a = \chi_a^' - i\chi_a''$

Where :

$\chi^' = \frac{1}{D}\gamma M_0 \omega_H \left[ \omega_H^2 - \left( 1-\alpha^2\right)\omega^2 \right]$

$\chi^'' = \frac{1}{D}\alpha\gamma M_0 \omega \left[ \omega_H^2 + \left( 1+\alpha^2\right)\omega^2 \right]$

$\chi^'_a = \frac{1}{D}\gamma M_0 \omega \left[ \omega_H^2 - \left( 1+\alpha^2\right)\omega^2 \right]$

$\chi^''_a = \frac{1}{D}2\alpha\gamma M_0 \omega^2\omega_H$

$D = \left[ \omega_H^2 - \left( 1+\alpha^2\right)\omega^2\right] + 4 \alpha^2\omega^2\omega_H^2$

Where : $\alpha$ , $\gamma$ , $M_0$ are constants and $\omega$ is the excitation frequency and $\omega_H$ that is linear with the applied magnetic field.

The porosity, which we calculate from the images, comes into play via the following equation :

$\omega = 2\pi\gamma M \times \left( 1 - 3 P\right)$

We will input these equations into the image analysis code. What I do is simply loop over all the holes detected into the image, calculate the local porosity and and calculate the susceptibility $\chi$ . Note I could get better performances by vectorizing the code.

def calcHeff(poro): 
    return 2*np.pi*Ms*(1-3*poro)

def calc_Chi(Heff,Omega): #,porosity,field
    Omega_H = ((f/gamma)-Heff)/(gamma)
    Chi = (alpha*gamma*Ms*Omega*(Omega_H**2- ((1+alpha**2)*Omega**2)))/((((Omega_H**2)-((1+alpha**2)*Omega**2))**2)+4*alpha**2*Omega**2*Omega_H**2)
    Chi_a = -(2*Omega**2*Ms*gamma*Omega_H*alpha)/((((Omega_H**2)-((1+alpha**2)*Omega**2))**2)+4*alpha**2*Omega**2*Omega_H**2)
    Chi_2 =Chi+Chi_a
    return Chi_2

Omega = H/gamma
for h_ in calcHeff(porosite):
    for o_ in Omega:
        Chi_2 = np.append(Chi_2,calc_Chi(h_,o_))

Let’s plot the corresponding graphs. First, the individual absorption curve for each individual wire. As we can see below, since the network is nicely ordered the majority of wire have an absorption field close to 5.8 kOe.

Plot of the absorption curves for each individual wire.

If we plot the sum, we have a slightly skewed curve as expected.

A typical absorption curve for an array of nanowires into an ordered template.

The peaks between 2-3 kOe are due to border effects in the triangulation: the code could be optimized by ignoring distances larger than a given threshold. A wider area should also be taken to average things, which I’ll do in the next example.

A more complex sample

Now we have a working algorithm, let’s use it on a larger image of a different template. As shown below the distribution of the holes is random: there are packed wires as well as (almost) isolated one.

Image constructed from multiple acquisitions. Here, the distribution is random.

The image correction and analysis steps are exactly as before, so I’ll jump to the triangulation :

Result of the triangulation, there are about 1600 holes.

The algorithm counts around 1600 holes and the estimated porosity is 7.8%, in good agreement with the expected value.

The calculation of the absorption curves yields the result below:

Plot of the absorption curves for each individual wire.

We can see that the distribution is much wider than for the previous case.

Indeed, ploting the sum shows a wide curve, skewed towards high field values (i.e. packed wires).

Sum of absorption curves for a random distribution.

I unfortunately do not have the original experimental data to plot the comparison properly. However, if you are interested, I suggest you read that paper or this one.

Final thoughts

That’s been a long post, so let’s summarize.

I have built an image analysis algorithm to analyse the distribution of holes in a porous template, that could be extend to similar problems : counting particles, objects on a map, etc.

We have seen how an image should be equalized and cleaned for a correct analysis using scikit-image and scipy’s image analysis modules. Note that matplotlib also has similar tools.

Finally, a physical model has been implemented into the image analysis algorithm, that reproduces experimental data pretty well.For the physicists reading this post, I’ll mention that this is a simplified model that does not take into account effects due to the microstructure of the wires, their orientation into the matrix and that is limited to first neighbors. So it could definitely be improved.

As usual, do not hesitate to leave a comment below 🙂

Image analysis : Delaunay triangulation to predict FMR linewidth