# 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

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)

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.

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.

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

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

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 :

The diagonal and off-diagonal terms are :

and

With : and

Where :

Where : , , are constants and is the excitation frequency and that is linear with the applied magnetic field.

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

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 . 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.

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

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.

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

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:

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

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 🙂