You can do it with Python! Structural analysis of two-dimensional colloidal crystals

Introduction

In this article, we will perform the following analysis using a two-dimensional colloidal crystal image.

1. Detection of colloidal position
2. Derivation of radial distribution function (like)
3. Counting and plotting the first neighboring atom
4. Visualization of polycrystalline grains using the coordination angle of nearby atoms as an index

I used the following colloidal images: yum:

Author Zephyris [CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0) or GFDL (http://www.gnu.org/copyleft/fdl.html)], [Wikimedia From Commons (link)

List of libraries to use

It is as follows.

import cv2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy import ndimage as ndi
from sklearn.neighbors import NearestNeighbors

Detection of colloidal position

A maximum filter is used to detect the colloidal position. This method does much the same as ImageJ's Find Maxima (https://imagej.nih.gov/ij/docs/menus/process.html#find-maxima).

First, load the image as shown below and convert it to grayscale. Also get the image size.

#Image loading
img_raw = cv2.imread('ColloidCrystal_10xBrightField_GlassInWater.jpg', 1)
#Grayscale
img = cv2.cvtColor(img_raw, cv2.COLOR_BGR2GRAY)
#Image height,Get width
h, w = img.shape[:2]

Particle detection using maximum filter

As the name implies, the maximum value filter is a filter that refers to the maximum brightness within a certain range around the pixel of interest. By using this, it is possible to detect the local maximum value of the image, that is, the position of colloidal particles.

Let's apply the filter with the size of the attention range (kernel) set to 9x9. (Noise is removed by smoothing before application)

#Blur(Smoothing)Noise removal by
img = cv2.blur(img,(3,3))
#Apply maximum filter
img_max = ndi.maximum_filter(img, size=9, mode='constant')

Click here for a comparison of images before and after applying the filter.

In these two images, ** the position where the brightness value is the same ** is the local maximum value, that is, the position of the particle. (Reference) This can be easily found using numpy.where.

#Extraction of local maximum value
pts_list = np.where(img == img_max)
pts_list = np.array(pts_list)

#Plot on the original image
fig = plt.figure(dpi=150,facecolor='white')
plt.gray()
plt.imshow(img)
plt.scatter(pts_list[1],pts_list[0],s=3,c='red',alpha=0.7)
plt.xlim(0,300)
plt.ylim(0,300)

The points obtained in this way are displayed above the original image and are shown below.

You can roughly detect the position of the particles. However, some particles have been detected twice. In order to remove this, we will summarize the points that are close to each other with their center of gravity as a representative.

Remove duplicate points

Use sklearn.neighbors.NearestNeighbors to calculate the distance between points.

nbrs = NearestNeighbors(n_neighbors=10, algorithm='ball_tree').fit(pts_list.T)
distances, indices = nbrs.kneighbors(pts_list.T)

Nearest Neighbor returns the distance and index to the nearest neighbor by specifying the nearest neighbor.

This time, we will calculate that the points existing within 3 pixels are derived from the same particle.

#Proximity threshold
th = 3

#Calculate the center of gravity with the proximity point below the threshold.
center_list = []
for d,i in zip(distances,indices):
    i = i[np.where(d < 3)]
    pts = pts_list.T[i]
    center = np.mean(pts,axis=0)
    center_list.append(center)
center_list = np.array(center_list).astype(np.int32)

#Remove duplicates
center_list = np.unique(center_list,axis=0)
center_list = center_list.T

The points obtained in this way are shown below in the same way. The points that overlapped in one particle are summarized. This completes particle detection:: thumbs up :: thumbs up:

Derivation of radial distribution function (like)

Next, we will derive the Radial Distribution Function (RDF). RDF indicates the probability that another particle will exist with respect to the distance from one particle, and is often used to describe the periodic structure of a crystal.

See below for the derivation of RDF and its calculation method. Calculation of radial distribution function --Koishi's Page

Use sklearn.neighbors.NearestNeighbors used earlier to calculate the proximity atoms.

#Calculate the position and distance of nearby atoms up to 40
nbrs = NearestNeighbors(n_neighbors=40, algorithm='ball_tree').fit(center_list.T)
distances, indices = nbrs.kneighbors(center_list.T)

#Calculate the distance to other particles existing within 50 pixels and create a histogram
dist_list = []
for d,i in zip(distances,indices):
    d = d[np.where(d < 50)][1:]
    dist_list.extend(d)
dist_list = np.array(dist_list)
dist_hist = np.histogram(dist_list,range=(0,50),bins=50)

#RDF calculation
rdf = dist_hist[0]/(4*np.pi*np.power(dist_hist[1][1:],2))

This is all right. This time I calculated RDF in the range up to 50 pixels. (Strictly different because it is not divided by the density at the end, but this time it is OK if the peak value is known)

The result is shown below.

You can see the peaks up to the 1st, 2nd, and 3rd proximity. The colloidal crystals this time are densely packed, but there are some defects, so the peaks are split or broad.

For the time being, it turns out that the distance to the first proximity particle is about 10 to 18 pixels.

Next, let's count the number of these first proximity particles.

Counting and plotting of the first neighboring atom

When the spherical colloid is packed in two dimensions, the coordination number with the first proximity particle is ** 6 **. The area below this is where the colloids are not neatly arranged. In other words, plotting the number of first proximity particles can illustrate defects in colloidal crystals.

Let's do it: v :: v: I also use sklearn.neighbors.NearestNeighbors. The threshold value to be regarded as the first neighboring atom is 16 px.

#Counting of first neighboring atoms
co_list = []
for d,i in zip(distances,indices):
    d = d[np.where(d < 16)][1:]
    co_list.append(d.shape[0])

In addition, the figure below shows the number of first proximity particles in a histogram.

The number of closest particles is generally distributed between 3 and 7. For now, let's plot in the range of (3,6).

fig = plt.figure(dpi=150,facecolor='white')
plt.imshow(img)
plt.scatter(center_list[1],center_list[0],s=3,c=co_list,cmap='rainbow',alpha=0.7)
plt.colorbar()
plt.clim(3,6)
plt.xlim(0,1000)
plt.ylim(0,1000)

The original image is shown on the left, and the plot of each particle color-coded by coordination number is shown on the right. Disordered part = The location of the defect can be illustrated: heart_eyes :: heart_eyes:

Finally, here is a plot over the entire original image: point_down:

Visualization of polycrystalline grains using the coordination angle of neighboring atoms as an index

Finally, let's visualize the polycrystalline grains of colloidal crystals. Specifically, two-dimensional mapping of crystal grains (see the figure below) as obtained by EBSD (Electron Backscatter Diffraction) is also performed for colloidal crystals.

(Source: https://www.ube-ind.co.jp/usal/documents/m1303_550.htm)

Such an analysis is also performed in a paper on colloidal crystals. For example, in Fig. 1 of the following paper, the vector formed by the particle of interest ⇒ the first proximity particle is classified according to the angle formed by the X-axis, and the crystal grain We are visualizing.

Reference: [Ramananarivo, S., Ducrot, E. & Palacci, J. Activity-controlled annealing of colloidal monolayers. Nat Commun 10, 3380 (2019).](Https://doi.org/10.1038/s41467-019- 11362-y)

Let's visualize with the same method!

theta_list = []
for d,i in zip(distances,indices):
    #Extract proximity particles excluding itself
    idx = i[np.where(d < 16)][1:]
    cnts = center_list.T[idx]
    #Get the position of the particle with the largest X value
    top = cnts[np.argmin(cnts[:,1])]
    #Get the position of the particle of interest
    center = center_list.T[i[0]]
    #Calculate the angle with the X axis
    axis = np.array([0,center[1]])
    u = top - center
    v = axis - center
    c = np.inner(u, v) / (np.linalg.norm(u) * np.linalg.norm(v) + 1e-99)
    a = np.rad2deg(np.arccos(np.clip(c, -1.0, 1.0))) - 90
    theta_list.append(a)

The histogram of the angles obtained in this way is shown below.

The angle seems to be within the range of ± 30 deg. Let's plot!

The original image is shown on the left, and the angle-coded plot of the first neighboring atom is shown on the right. Finally, let's do this color coding for the entire image! You can clearly see the grain of the colloidal crystal: yum :: yum: It seems to be noisy compared to the mapping of the papers I referred to, but it should be enough for a defeat job.

Appeal to rivals with colorful figures!

Finally

Recently, I received a request for image analysis from a person who is doing colloid-related research, and I was doing various analyzes while staying at home. For the time being, I've just completed a paragraph, so it feels like I've put together my know-how based on the images of colloidal crystals that were rolling on the net.

In fact, the images obtained in experiments are often very noisy and may require proper filtering and denoising. Regarding this, it is basically a game that comes out, or a part that requires trial and error. I would like to summarize useful methods soon. It's going to be a messy article.

Well then ~.