HRSTEM images characterization

In this project, we have used multiple Python libraries to do a full analysis of High Resolution Scanning Transmission Electron Microscope (HRSTEM) images of LSMO/STO Perovskites.

In this particular project, we will use atomap module for the detection of atoms in a HRSTEM image and then connect the atoms with colored lines that highlights the distance between them. This script is usefull to highlight any local distortions in the image.

Atomic columns detection

To use the following script, we assume that the user has a high resolution dark-field or bright-field image that had its scan drift distortions already corrected.

Let’s start with a HRSTEM bright field image of a LSMO/STO interface:

First, we import all the modules:

import atomap.api as am
import numpy as np
import hyperspy.api as hs
from atomap.tools import remove_atoms_from_image_using_2d_gaussian

Then Atomap finds initial positions for the atomic columns by using a peak finding algorithm. This algorithm needs to know the smallest peak separation (minimum pixel separation of the features):

s = hs.load(r'C:/Users/adchmielews/Documents/HRSTEM analysis/LSMO-STO bulk.dm3')
s.data=np.divide(1, s.data)
s_peaks = am.get_feature_separation(s, separation_range=(20,30))
s_peaks.plot(norm = 'linear')

Here 22 pixels seems to be a good peak separation, so we can store it in atom_positions:

atom_positions = am.get_atom_positions(s, separation = 22)

On the top left corner of the right image, we notice that few atoms are remaining, we can remove them by using an implement GUI function, in which we can add/remove atoms by mouse clicking directly on the image:

atom_positions = am.add_atoms_with_gui(s, atom_positions,distance_threshold=8,norm='linear')

Now that we have a rough estimate of the positions of our first type of atomic columns (usually composed by the heaviest type of atoms), we will refine it with a center of mass algorithm followed by a 2D Gaussian fitting:

sublattice_A = am.Sublattice(atom_positions, image=s.data)
sublattice_A.find_nearest_neighbors()
sublattice_A.refine_atom_positions_using_center_of_mass()
sublattice_A.refine_atom_positions_using_2d_gaussian(mask_radius=None)
sublattice_A.construct_zone_axes()

Now that we have the atomic columns of atom A, we will plot the planes and apply « sublattice_A.find_missing_atoms_from_zone_vector() » to have a rough estimate of the positions of atoms B.

direction_001 = sublattice_A.zones_axis_average_distances[1]
B_positions = sublattice_A.find_missing_atoms_from_zone_vector(direction_001,vector_fraction=0.5)

Finally, we will refine the position with the center of mass + 2D Gaussians fitting. Important: for a better fitting with the Gaussians, we can remove atoms A from the image, so it is not overlapping with the signal of atoms B.

image_without_A = remove_atoms_from_image_using_2d_gaussian(sublattice_A.image, sublattice_A)
sublattice_B = am.Sublattice(B_positions, image_without_A, color='blue', name='B')
sublattice_B.construct_zone_axes()
sublattice_B.refine_atom_positions_using_center_of_mass()
sublattice_B.refine_atom_positions_using_2d_gaussian(mask_radius=0.5)

We use the exact same workflow for the detection of atoms C, first a rough estimate by finding the missing atoms out of the planes:

image_without_A = remove_atoms_from_image_using_2d_gaussian(sublattice_A.image, sublattice_A)
sublattice_B = am.Sublattice(B_positions, image_without_A, color='blue', name='B')
sublattice_B.construct_zone_axes()
sublattice_B.refine_atom_positions_using_center_of_mass()
sublattice_B.refine_atom_positions_using_2d_gaussian(mask_radius=0.5)

Then removing the positions B of the image using ‘remove_atoms_from_image_using_2d_gaussian(s_image_without_A.data, sublattice_B)’ and refinement with center of mass + 2D Gaussians:

direction_001 = sublattice_B.zones_axis_average_distances[0]
C_positions = sublattice_B.find_missing_atoms_from_zone_vector(direction_001)
s_image_without_A = hs.signals.Signal2D(image_without_A)
image_without_B = remove_atoms_from_image_using_2d_gaussian(s_image_without_A.data, sublattice_B)
s_image_without_B = hs.signals.Signal2D(image_without_B)
sublattice_C = am.Sublattice(C_positions, s_image_without_B, color='green', name='C')
sublattice_C.find_nearest_neighbors()
sublattice_C.refine_atom_positions_using_center_of_mass()
sublattice_C.refine_atom_positions_using_2d_gaussian(mask_radius=0.5)
sublattice_C.construct_zone_axes()

We now have all the 3 types of atomic columns detected !

Distance between atoms calculation

Here’s the high-level functionality of the script:

1. It loads an HRSTEM image and positions of atoms (presumably determined using some other method) from files.
2. It calculates pairwise Euclidean distances between the positions of atoms. For each atom, it finds the k closest atoms (neighbors). The first neighbor is the atom itself. It then creates a plot where atoms are connected by lines. The color of the line depends on the distance between atoms and the colormap chosen. Only lines corresponding to distances within the provided thresh_min and thresh_max are shown on the plot. The plot also includes a scale bar and colorbar for reference, and the image is displayed in grayscale in the background. Finally,the function returns a list of line coordinates that met the distance thresholds.
3. This distance_plot function is then called twice. First, to simply display the plot, and second, to store the line coordinates that meet the distance thresholds.
4. Finally, the script removes duplicates from the line coordinates and sorts the coordinates in each line. It then prints the unique lines (pairs of atom positions) that meet the distance thresholds.

import numpy as np
from matplotlib_scalebar.scalebar import ScaleBar
import math
from scipy.spatial.distance import pdist, squareform
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
import hyperspy.api as hs
from matplotlib.colors import Normalize
from matplotlib.cm import ScalarMappable

# Load data
s = hs.load(r'E:/POSTDOC/2019-05-15 Al20Ga2O3-001-zone axis O imaging/HAADF-ABF Batch 1/HAADF-ABF Batch 1/ZOOM of 88.tif')
X = np.loadtxt(r'E:/POSTDOC/2019-05-15 Al20Ga2O3-001-zone axis O imaging/HAADF-ABF Batch 1/HAADF-ABF Batch 1/ZOOM of 88 O positions.txt')

def distance_plot(k, Calibration, color, thresh_max, thresh_min):
    """
    This function will display all the atoms in a HRSTEM image and connect them with colored lines.
    
    Parameters
    ----------
    k : int, number of neighbors, the first neigbor is the atom itself, so put 2 to have at least one connexion to a nearest
        neighbor.
    Calibration : float, in nanometer. For your image calibration.
    color : chose a color out of matplotlib.colors.ListedColormap
    thresh_max: float, maximum distance threshold that will be displayed on the image (in pm)
    thresh_min: float, minimum distance threshold that will be displayed on the image (in pm)
    
    Examples
    ----------
    distance_plot(k=5, Calibration=0.010863, color='magma',thresh_max=417, thresh_min=350)
    
    """
    
    # Colormap
    cmap=plt.get_cmap(color)

    # Number of atoms in the data
    N=len(X)

    # Calculate pairwise Euclidean distances
    distmat = squareform(pdist(X, 'euclidean'))

    # Function for calculating distance between two points
    def dist(x1, y1, x2, y2):
        return math.sqrt((x2-x1)**2 + (y2-y1)**2)

    # Select the kNN for each datapoint
    neighbors = np.sort(np.argsort(distmat, axis=1)[:, 0:k])

    # Get edge coordinates and distance
    coordinates = np.zeros((N, k, 2, 2))
    distance = []

    for i in np.arange(N):
        for j in np.arange(k):
            coordinates[i, j, :, 0] = np.array([X[i,:][0], X[neighbors[i, j], :][0]])
            coordinates[i, j, :, 1] = np.array([X[i,:][1], X[neighbors[i, j], :][1]])
            distance.append(dist(coordinates[i][j][0][0], coordinates[i][j][0][1], coordinates[i][j][1][0], coordinates[i][j][1][1]))

    # Threshold to keep only coordinates below a certain value
    coordinates = coordinates.tolist()
    new_coordinates = []

    fig, ax = plt.subplots(1, 1)

    for i in np.arange(len(coordinates)):
        for j in np.arange(k):
            if dist(coordinates[i][j][0][0], coordinates[i][j][0][1], coordinates[i][j][1][0], coordinates[i][j][1][1]) > thresh_max/(Calibration*1000):
                new_coordinates.append(coordinates[i][j])
                distance.remove(dist(coordinates[i][j][0][0], coordinates[i][j][0][1], coordinates[i][j][1][0], coordinates[i][j][1][1]))
            if dist(coordinates[i][j][0][0], coordinates[i][j][0][1], coordinates[i][j][1][0], coordinates[i][j][1][1]) < thresh_min/(Calibration*1000):
                new_coordinates.append(coordinates[i][j])
                distance.remove(dist(coordinates[i][j][0][0], coordinates[i][j][0][1], coordinates[i][j][1][0], coordinates[i][j][1][1]))

    distance = [i*Calibration*1000 for i in distance]
    new_coordinates = np.round_(new_coordinates, decimals=8)
    new_coordinates = new_coordinates.tolist()
    coordinates = np.asarray(coordinates)
    
    #Plotting
    
    norm = Normalize(vmin=thresh_min, vmax=thresh_max)
    line_colors = cmap(norm(distance))
    coco = coordinates.reshape((N*k, 2, 2)).tolist()
    coco = np.round_(coco, decimals=8)
    coco = coco.tolist()
    result = np.array([position for position in coco if position not in new_coordinates])
    
    lines = LineCollection(result, linewidth=3, colors=line_colors)
    
    ax.scatter(X[:, 0], X[:, 1], c='black', s=5)
    
    ax.add_collection(lines)
    
    mappable = ScalarMappable(cmap=cmap, norm=norm)
    cb = fig.colorbar(mappable)
    cb.set_label('Distance (pm)')
    lines.set_clim(vmin=thresh_min, vmax=thresh_max)
    scalebar = ScaleBar(Calibration*1000, 'pm', location=4)
    plt.gca().add_artist(scalebar)
    
    plt.imshow(s, cmap="gray")
    plt.sci(lines)
    plt.show()
    
    return result

Highlighting defects

In the following part, we will look at the Aluminum-Gallium-Oxide/Gallium-Oxide (AlGaO/GaO) interface. In this system, point defects induced during the growth process are often created. It is then necessary to detect and highlight them in order to see their distribution in the sample. One defect type that is seen in AlGaO is the interstitial atom that will usually sit in between two atoms (see the yellow circle on the figure below). To detect its position, we have to play around the thresh_min and thresh_max distance thresholds as the distance between an interstitial atom and its nearest neighbord will be very small (around 100 pm in the case shown below).

################# Defects detection ###########################################

# Sort the points in each line
sorted_coordinates = np.sort(coordinates, axis=1)

# Reshape the coordinates to a 2D array, where each row represents a line
reshaped_coordinates = sorted_coordinates.reshape(-1, sorted_coordinates.shape[-1]*2)

# Remove duplicate rows (lines)
unique_lines = np.unique(reshaped_coordinates, axis=0)

# Reshape the unique lines back to the original shape
unique_lines = unique_lines.reshape(-1, 2, 2)

print("Unique lines:")
print(unique_lines)

def find_duplicates(arr):
    seen = set()
    duplicates = set()

    for line in arr:
        for point in line:
            # Convert coordinate to a tuple so it can be inserted into a set
            point_tuple = tuple(point)
            if point_tuple in seen:
                duplicates.add(point_tuple)
            seen.add(point_tuple)
    
    return duplicates

duplicated_coordinates = find_duplicates(unique_lines)

print("Duplicated coordinates:")
for coord in duplicated_coordinates:
    print(coord)

from matplotlib.patches import Circle

# Your code before this is unchanged ...

print("Duplicated coordinates:")
for coord in duplicated_coordinates:
    print(coord)

fig, ax = plt.subplots(1, 1)
# Draw circles around duplicated coordinates
for coord in duplicated_coordinates:
    circle = Circle(coord, radius=10, fill=False, edgecolor='yellow', linewidth=2)
    ax.add_patch(circle)
    plt.imshow(s, cmap="gray")

# Adjust the plot limits if necessary. This might not be needed depending on your data and plot settings.
# ax.set_xlim(center[0] - 10, center[0] + 10)
# ax.set_ylim(center[1] - 10, center[1] + 10)

    plt.show()

Polarization

To be continued..