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:

800px-ColloidCrystal_10xBrightField_GlassInWater.jpg 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. fig1.png

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

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

fig4.png

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. ダウンロード.png

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)

fig5.png

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: fig9.png

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.

img_m1303_02.gif (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.

fig6.png

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

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! fig8.png 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 ~.

Recommended Posts

You can do it with Python! Structural analysis of two-dimensional colloidal crystals
Python | What you can do with Python
If you write TinderBot in Python, she can do it
[Python] What do you do with visualization of 4 or more variables?
You can try it with copy! Let's draw a cool network diagram with networkx of Python
What you can do with API vol.1
Two-dimensional unsteady heat conduction analysis with Python
What you can do with programming skills
Until you can use opencv with python
You will be an engineer in 100 days --Day 35 --Python --What you can do with Python
You can easily create a GUI with Python
You can see it! Comparison of optimization methods (2020)
If you guys in the scope kitchen can do it with a margin ~ ♪
Static analysis of Python code with GitLab CI
Two-dimensional elastic skeleton geometric nonlinear analysis with Python
[For beginners] You can do it from scratch! Creating APIs with AWS SAM and outputting OpenAPI documentation in Python
Consideration when you can do a good job in 10 years with Python3 and Scala3.
Until you can install blender and run it with python for the time being
[OpenCV / Python] I tried image analysis of cells with OpenCV
Calculate the regression coefficient of simple regression analysis with python
Challenge principal component analysis of text data with Python
Planar skeleton analysis with Python (4) Handling of forced displacement
What you can and cannot do with Tensorflow 2.x
Until you can do simple image recognition with Jupyter
How much do you know the basics of Python?
Voice analysis with python
Do Houdini with Python3! !! !!
Voice analysis with python
Data analysis with Python
3. Natural language processing with Python 5-1. Concept of sentiment analysis [AFINN-111]
Plane skeleton analysis with Python (3) Creation of section force diagram
Until you can install your own Python library with pip
Perform isocurrent analysis of open channels with Python and matplotlib
Do Django with CodeStar (Python3.6.8, Django2.2.9)
[Python] Morphological analysis with MeCab
[Co-occurrence analysis] Easy co-occurrence analysis with Python! [Python]
Do Django with CodeStar (Python3.8, Django2.1.15)
Sentiment analysis with Python (word2vec)
Static analysis of Python programs
Planar skeleton analysis with Python
Japanese morphological analysis with Python
Muscle jerk analysis with Python
Made it possible to convert PNG to JPG with Pillow of Python
It seems that you can now write gate books with blueqat
If you know Python, you can make a web application with Django
Don't write Python if you want to speed it up with Python
What to do if you can't install pyaudio with pip #Python
Practical exercise of data analysis with Python ~ 2016 New Coder Survey Edition ~
Introduction to Structural Equation Modeling (SEM), Covariance Structure Analysis with Python
IOS application development with Python (kivy) that even monkeys can do
From the introduction of JUMAN ++ to morphological analysis of Japanese with Python
Cosine similarity matrix? You can get it right away with NumPy