Chapter 3: Imaging

Image Analysis¶

part of

MSE672: Introduction to Transmission Electron Microscopy

Spring 2025
by Gerd Duscher

Microscopy Facilities
Institute of Advanced Materials & Manufacturing
Materials Science & Engineering
The University of Tennessee, Knoxville

Background and methods to analysis and quantification of data acquired with transmission electron microscopes.

Load important packages¶

Check Installed Packages¶

import sys
import importlib.metadata
def test_package(package_name):
    """Test if package exists and returns version or -1"""
    try:
        version = importlib.metadata.version(package_name)
    except importlib.metadata.PackageNotFoundError:
        version = '-1'
    return version

if test_package('pyTEMlib') < '0.2026.1.0':
    print('installing pyTEMlib')
    !{sys.executable} -m pip install  --upgrade pyTEMlib -q

print('done')

done

Import all relevant libraries¶

Besides the ususal libraries we also load the blob detectors of the scipy-image package

%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
import sys
import os
if 'google.colab' in sys.modules:
    from google.colab import output
    from google.colab import drive
    output.enable_custom_widget_manager()

# our blob detectors from the scipy image package
import skimage    #.feature import blob_dog, blob_log, blob_doh

# Multidimensional Image library
import scipy

import time
# Import libraries from the book

import pyTEMlib

# it is a good idea to show the version numbers at this point for archiving reasons.
print('pyTEM version: ',pyTEMlib.__version__)

pyTEM version:  0.2026.1.0

Load an atomic resolution image¶

As an example we will use again p1-3-hr3.dm3 in the TEMdata directory. This image is of a nO particle on graphene.

# ---- Input ------
load_example = True
# -----------------
if not load_example:
    if 'google.colab' in sys.modules:
        drive.mount("/content/drive")

    fileWidget = pyTEMlib.file_tools.FileWidget()

# ---- Input ------
file_name = 'p1-3-hr3.dm3'
file_name = 'p1-hr4-ZnOonGraphite.dm3'
# -----------------
      
if load_example:
    if 'google.colab' in sys.modules:
      if not os.path.exists('./'+file_name):
        !wget  https://github.com/gduscher/MSE672-Introduction-to-TEM/raw/main/example_data/p1-3-hr3.dm3
        !wget  https://github.com/gduscher/MSE672-Introduction-to-TEM/raw/main/example_data/p1-hr4-ZnOonGraphite.dm3
    else:
        datasets = pyTEMlib.file_tools.open_file('../example_data/'+file_name)
        dataset = datasets['Channel_000']
      
else:
    datasets = fileWidget.datasets
    dataset = fileWidget.selected_dataset
# Load file

view = dataset.plot()

Fourier Transform of Image¶


fft_mag = pyTEMlib.image_tools.power_spectrum(dataset)
view = fft_mag.plot()

Spot Detection¶

# ---- Input ---------
spot_threshold = .2
# --------------------
## pixel_size in recipical space
rec_scale_x = 1/dataset.x[-1]  
rec_scale_y = 1/dataset.y[-1]

## Field of View (FOV) in recipical space please note: rec_FOV_x = 1/(scaleX*2)
rec_FOV_x = rec_scale_x * dataset.shape[0] /2.
rec_FOV_y = rec_scale_y * dataset.shape[1] /2.


## Field ofView (FOV) in recipical space
rec_extend = (-rec_FOV_x,rec_FOV_x,rec_FOV_y,-rec_FOV_y)

## Needed for conversion from pixel to Reciprocal space
rec_scale = np.array([rec_scale_x, rec_scale_y,1])
center = np.array([int(dataset.shape[0]/2), int(dataset.shape[1]/2),1] )


# spot detection ( for future referece there is no symmetry assumed here)

spots_random =  (skimage.feature.blob_log(fft_mag,  max_sigma= 5 , threshold=spot_threshold)-center)*rec_scale

print(f'Found {spots_random.shape[0]} reflections')
spots_random[:,2] = np.linalg.norm(spots_random[:,0:2], axis=1)
spots_index = np.argsort(spots_random[:,2])
spots = spots_random[spots_index]

## plot Fourier transform and found spots
view = fft_mag.plot()
plt.gca().scatter(spots[:,0], spots[:,1], c='Orange',  alpha = 0.2, label='spots');

Found 19 reflections

Adaptive Fourier Filtering¶

We mask the fourier transformed image so that the information can pass through is selected.

The information is in the spots and in the center of the Fourier transformed image,the rest is noise.

Please modify the radius of the mask of the reflections and the low-path area in the code below and notice the effects on the Fourier filtered image.

# Input 
reflection_radius = 0.3 # in 1/nm
low_pass = 1/.34 # in 1/nm diameter of mask for low pass filter
FOV_x = dataset.x[-1]
FOV_y = dataset.y[-1]


#prepare mask
pixels = (np.linspace(0,dataset.shape[0]-1,dataset.shape[0])-dataset.shape[0]/2)* rec_scale_x

x,y = np.meshgrid(pixels,pixels);
mask = np.zeros(dataset.shape)

# mask reflections
for spot in spots:
    mask_spot = (x-spot[0])**2+(y-spot[1])**2 < reflection_radius**2 # make a spot 
    mask = mask + mask_spot# add spot to mask
    
# mask zero region larger (low-pass filter = intensity variations)

mask_spot = x**2+y**2 < low_pass**2 
mask = mask + mask_spot

mask[np.where(mask>1)]=1 # just in case of overlapping disks

plt.figure()
ax1 = plt.subplot(1,2,1)
#ax1.imshow(mask)
fft_filtered = np.fft.fftshift(np.fft.fft2(np.array(dataset)))*mask.T
ax1.imshow(np.log(1+np.abs(fft_filtered)).real,extent=rec_extend, origin = 'upper')
plt.xlabel('reciprocal distance [1/nm]')
ax2 = plt.subplot(1,2,2)
filtered = np.fft.ifft2(np.fft.fftshift(fft_filtered).T)

real_extent = (0,FOV_x,FOV_y,0)
    
ax2.imshow(filtered.real.T,extent=real_extent, origin = 'upper')
plt.xlabel('distance [nm]');

Check on filtered images¶

We don’t want to filter anything out that caries information, or at least we want to be aware of that. An easy check is to subtract the filtered imag fromt he image and to determine that only noise is left.

Please note that any processing can be easily determined in the Fourier transformed, so be meticulous on reporting what you did to an image.

plt.figure()
ax1 = plt.subplot(1,2,1)
ax1.imshow(dataset-filtered.real.T,extent=real_extent, origin = 'upper')
plt.xlabel(' distance [nm]')

ax2 = plt.subplot(1,2,2)
fft_difference = np.fft.fftshift(np.fft.fft2(np.array(dataset-filtered.real.T)))
ax2.imshow(np.log(1+np.abs(fft_difference)).real,extent=rec_extend, origin = 'upper')
plt.xlabel('reciprocal distance [1/nm]');

Rotational Symmetry¶

In our context of symmetry, we just need to deal with the discrete values of Θ = 2π/n for the angle of rotation.

In two dimensios we have then the rotation matrix:

Undefined control sequence: \array at position 20: …-fold} = \left[\̲a̲r̲r̲a̲y̲{ \cos{( 2\pi/n…

C_{n-fold} = \left[\array{ \cos{( 2\pi/n)} & \sin{( 2\pi/n)}\\ -\sin{( 2\pi/n)} & \cos{ (2\pi/n)}\\}\right]

If we subtract all spots from all rotated spots we have a set of distances where for each spot there is minimal distance to the next spot. If we have a very small distance for each original spot, we have a found rotational symmetry operation.

from itertools import product

for n in [2,3,4,5,6]:
    C = np.array([[np.cos(2*np.pi/n), np.sin(2*np.pi/n),0],[-np.sin(2*np.pi/n), np.cos(2*np.pi/n),0], [0,0,1]])
    sym_spots = np.dot(spots,C)
    dif = []
    for p0, p1 in product(sym_spots[:,0:2], spots[:,0:2]):
        dif.append(np.linalg.norm(p0-p1))
    dif = np.array(sorted(dif))
    #print(dif[0:spots.shape[0]])
    if dif[int(spots.shape[0]*.7)] < 0.2:
        
        print(f'Found {n}-fold symmetry in diffraction pattern')

Found 2-fold symmetry in diffraction pattern
Found 3-fold symmetry in diffraction pattern
Found 6-fold symmetry in diffraction pattern

Mirror symmetry¶

Any mirror axis has to go through the origin.

Let’s consider the points as vectors and let’s condier the miror axis as a vector field that lays on that axis.

The mirror axis of any point has to go through that point plus the half of the vector to another point.

Reflection across a line through the origin in two dimensions can be described by the following formula

\operatorname{Ref}_l(v) = 2\frac{v \cdot l}{l \cdot l}l - v,

(2)

where $v$ is the vector of a point, while $l$ is a vector on the mirror axis on which $v$ is reflected on. $v \cdot l$ denotes the dot product of $v$ and $l$ .

mirror_axes = []
plt.figure()
plt.imshow(fft_mag.T,extent=rec_extend, origin = 'upper')
plt.xlabel('reciprocal distance [1/nm]')
plt.scatter(spots[:,0], spots[:,1], c='Red',  alpha = 0.5, label='spots');

for spot in spots:
    if spot[2] > .1:
        mirror_spots = []
        for spot2 in spots:
            if spot2[2]>.1:
                l = spot[0:2]+spot2[0:2]/2
                v = spot[0:2]
                ref = 2*np.dot(v,l) / np.dot(l,l)*l-v
                mirror_spots.append(ref)
        mirror_spots = np.array(mirror_spots)
        dif = []
        for p0, p1 in product(mirror_spots, spots[:,0:2]):
            dif.append(np.linalg.norm(p0-p1))
        dif = np.array(sorted(dif))
        #print(dif[0:25])
        #print(int(spots.shape[0]/2))
        #print(dif[int(spots.shape[0]/2)])
        if dif[int(spots.shape[0]/2)] < 0.5:
            #print(l,dif[0:25])
            mirror_axes.append(l)
            axis=l/np.linalg.norm(l)
            plt.plot([-l[0],l[0]],[-l[1],l[1]],c='yellow')
            print(f'Found mirror {axis} axis in diffraction pattern')
print(f'Found {len(mirror_axes)} mirror axes')

Found 0 mirror axes

Reference Crystals¶

graphite  = pyTEMlib.crystal_tools.structure_by_name('Graphite')
zno = pyTEMlib.crystal_tools.structure_by_name('ZnO')
### Define exxperimental parameters:
tags_experiment= {}
tags_experiment['acceleration_voltage'] = 200.0 *1000.0 #V
tags_experiment['new_figure'] = False
tags_experiment['plot FOV'] = 30
tags_experiment['convergence_angle'] = 0
tags_experiment['zone_hkl'] = np.array([0,0,1])  # incident neares zone axis: defines Laue Zones!!!!
tags_experiment['Sg_max'] = .02 # 1/Ang  maximum allowed excitation error ; This parameter is related to the thickness
tags_experiment['hkl_max'] = 2  # Highest evaluated Miller indices


######################################
# Diffraction Simulation of Graphite #
######################################

graphite_diff = {}
graphite_diff = pyTEMlib.diffraction_tools.get_bragg_reflections(graphite, tags_experiment, verbose=True) 

######################################
# Diffraction Simulation of ZnO      #
######################################
zno_diff = {}
zno_diff = pyTEMlib.diffraction_tools.get_bragg_reflections(zno, tags_experiment, verbose=True)

Of the 22 possible reflection 22 are allowed.
Of those, there are 22 in ZOLZ  and 0 in HOLZ
Of the 0 forbidden reflection in ZOLZ  0 can be dynamically activated.
Of the 24 possible reflection 24 are allowed.
Of those, there are 24 in ZOLZ  and 0 in HOLZ
Of the 0 forbidden reflection in ZOLZ  0 can be dynamically activated.

The diffraciton calculation is in Angstrom (or better 1/Ang) but the image is in nm. So we convert to 1/nm


spots_Graphite = pyTEMlib.diffraction_tools.plotting_coordinates(graphite_diff['allowed']['g'])  
spots_ZnO = pyTEMlib.diffraction_tools.plotting_coordinates(zno_diff['allowed']['g'])

resolution = 0.1#nm

plt.figure()
plt.imshow(fft_mag.T,extent=np.array(rec_extend), origin = 'upper')
plt.xlabel('reciprocal distance [1/nm]')
plt.scatter(spots_Graphite[:,0], spots_Graphite[:,1], c='Red',  alpha = 0.5, label='Graphite');
plt.scatter(spots_ZnO[:,0], spots_ZnO[:,1], c='orange',  alpha = 0.5, label='ZnO');
plt.scatter(spots[:,0], spots[:,1], c='green',  alpha = 0.5, label='exp.');

plt.legend(loc=1);

Reflections in Polar Coordinates¶

A more interesting way of comparing a simulation and experiment is to compare the spots in polar coordinates:

conversion to Euclidean space:

\begin{align} x &= r \cos\varphi \\ y &= r \sin\varphi \end{align}

(3)

conversion to polar coordinates:

\begin{align} r &= \sqrt{x^2 + y^2} \\ \varphi &= \operatorname{atan2}(y, x) \end{align}

(4)

Note:
The Bragg reflection from the simulation are already in polar coordinates.

##  Polar Coordinates of experimental and reference lattice

def cart2pol(points):
    rho = np.linalg.norm(points[:,0:2], axis=1)
    phi = np.arctan2(points[:,1], points[:,0])
    return(rho, phi)

def pol2cart(rho, phi):
    x = rho * np.cos(phi)
    y = rho * np.sin(phi)
    return(x, y)


def xy2polar(points, rounding = 1e-3):
    """
    Conversion from carthesian to polar coordinates
    
    the angles and distances are sorted by r and then phi
    The indices of this sort is also returned
    
    input points: numpy array with number of points in axis 0 first two elements in axis 1 are x and y
    
    optional rounding in significant digits 
    
    returns r,phi, sorted_indices
    """
    
    r,phi = cart2pol(points)
    
    phi = phi-phi.min() # only positive angles
    r = (np.floor(r/rounding) )*rounding # Remove rounding error differences

    sorted_indices = np.lexsort((phi,r) ) # sort first by r and then by phi
    r = r[sorted_indices]
    phi = phi[sorted_indices]
    
    return r, phi, sorted_indices

## Transfer everything to polar coordinates 
graphite_spherical = graphite_diff['allowed']['g']
zno_spherical = zno_diff['allowed']['g']

exp_r, exp_phi = cart2pol(spots) # just in polar coordinates

angleI = np.argmin(np.abs(exp_r-graphite_spherical[0, 0]) )
angle = exp_phi[angleI] - graphite_spherical[0, 1] ## Determine rotation angle

graphite_rotation = angle # rotation
print(f" Rotated Graphite SAD by {np.degrees(angle):.1f}°")


angleI2 = np.argmin(np.abs(exp_r-zno_spherical[0,0]) )
angle2 = exp_phi[angleI2] - zno_spherical[0, 1] ## Determine rotation angle


zno_rotation = angle2 # rotation
print(f" Rotated ZnO SAD by {np.degrees(angle2):.1f}°")


fft_mag.plot()

x, y= pol2cart(exp_r, exp_phi)
plt.scatter(x,y, c='green',  alpha = 0.2,label='exp')

x, y= pol2cart(graphite_spherical[:, 0]*10, graphite_spherical[:, 1] + graphite_rotation)
plt.scatter(x,y, c='red',  alpha = 0.4,label='Graphite')

x,y = pol2cart(zno_spherical[:, 0]*10, zno_spherical[:, 1] + zno_rotation)
plt.scatter(x,y, c='orange',  alpha = 0.4,label='ZnO')



plt.xlim(-8,8);plt.ylim(-8,8)
plt.legend(loc=1);

 Rotated Graphite SAD by 130.9°
 Rotated ZnO SAD by 120.0°

Calibrate Distortions with Reference Crystal¶

g = gx = gy = rec_scale_x

spots_reference = np.array( pol2cart(graphite_spherical[:, 0]*10, graphite_spherical[:, 1] + graphite_rotation)).T
#spots_reference = np.array(pol2cart(ZnO_r, ZnO_phi)).T

spots_experiment = spots[:,0:2]
dist_graphite = np.linalg.norm(spots_reference, axis=1)

distance_experiment = np.linalg.norm(spots_experiment, axis=1)
first_reflections = abs(distance_experiment - dist_graphite.min()) < 2
print('Evaluate ', first_reflections.sum(), 'reflections')
reference_reflections = spots_experiment[first_reflections]

import scipy.optimize as optimization
def func(params, xdata, ydata):
    dgx , dgy = params
    return (np.sqrt((xdata*dgx)**2 + (ydata*dgy)**2 ) - dist_graphite.min())

x0 = [1.,0.999]
[dgx, dgy], sig = optimization.leastsq(func, x0, args=(reference_reflections[:,0], reference_reflections[:,1]))

print('distortion x:', dgx, 'y: ', dgy) 
gx /=dgx
gy /=dgy

spots_experiment = spots_experiment*(dgx,dgy)

Evaluate  18 reflections
distortion x: 0.9281637216731751 y:  0.9157357742878829

Plot Corrected Image and Reference Lattice¶

spots_Graphite = pyTEMlib.diffraction_tools.plotting_coordinates(graphite_diff['allowed']['g'], rotation=graphite_rotation)  
spots_ZnO = pyTEMlib.diffraction_tools.plotting_coordinates(zno_diff['allowed']['g'], rotation= zno_rotation)

FOV_cor = (0,FOV_x/dgx,FOV_y/dgy,0)
rec_FOV_cor = (-rec_FOV_x*dgx,rec_FOV_x*dgx,rec_FOV_y*dgy,-rec_FOV_y*dgy)

plt.figure()
plt.title(' Graphite and ZnO in ' + str(tags_experiment['zone_hkl']), fontsize=20)     

plt.imshow(fft_mag.T,extent=rec_FOV_cor, origin = 'upper')
plt.xlabel('reciprocal distance [1/nm]')
plt.scatter(spots_Graphite[:,0], spots_Graphite[:,1], c='Red',  alpha = 0.5, label='Graphite');
plt.scatter(spots_ZnO[:,0], spots_ZnO[:,1], c='orange',  alpha = 0.5, label='ZnO');
plt.scatter(spots[:,0]*dgx, spots[:,1]*dgy, c='green',  alpha = 0.5, label='exp.');

plt.xlim(-7,7)
plt.ylim(-7,7)
plt.legend(loc=1)

Plot Histograms of Distances¶

Next we want to plot the radially summed profile of the diffractogram and compare it with the contrast transfer function.

from scipy.interpolate import interp1d
from scipy.ndimage import map_coordinates


def cartesian2polar(x, y, grid, r, t, order=3):

    R,T = np.meshgrid(r, t)

    new_x = R*np.cos(T)
    new_y = R*np.sin(T)

    ix = interp1d(x, np.arange(len(x)))
    iy = interp1d(y, np.arange(len(y)))

    new_ix = ix(new_x.ravel())
    new_iy = iy(new_y.ravel())

    
    return map_coordinates(grid, np.array([new_ix, new_iy]),
                            order=order).reshape(new_x.shape)

Plot Distances and CTF¶

please note that we use the skimage.transform.warp_polar function to obtain a radial average of the fourier transform.

# --------- Input -----------
Cs = 2.2
defocus = -180 # underfocus is negative
acceleration_voltage = 200*1e3
# ---------------------------------
wavelength = pyTEMlib.utilities.get_wavelength(acceleration_voltage, unit='nm') # in nm


polar_projection = skimage.transform.warp_polar(fft_mag)
below_zero = polar_projection<0.
polar_projection[below_zero]=0.

# Sum over all angles (axis 1)
profile = polar_projection.sum(axis=0)
profile = (profile-profile[1:400].min())/profile.max()*10
print(profile[-1])

k =np.linspace(1,len(profile),len(profile))*rec_scale_x
print (profile)

def calculateScherzer(wavelength, Cs):
    # Calculate the Scherzer defocus. Cs is in mm, lambda is in nm
    # Everything is trasnfered to m
    #The returned value is in nm
    Cs=Cs*1e-3 # now in m
    scherzer=-1.155*(Cs*wavelength*1e-9)**0.5 # in m
    scherzer=scherzer*1e9 # convert to nm
    return scherzer

scherzer =calculateScherzer(wavelength, Cs)

print(f'Scherzer defocus is {scherzer:.1f} nm')

def calculateCTF(wavelength, Cs, defocus,k3):
    # everything in nm
    Cs=Cs*10**6
    ctf=np.sin(np.pi*defocus*wavelength*k**2+0.5*np.pi*Cs*wavelength**3*k**4)
    return ctf
ctf = calculateCTF(wavelength/10, Cs, defocus,k)


fig = plt.figure()
plt.hist(exp_r,45,facecolor='green', alpha=0.75, label='experiment')
plt.hist(graphite_spherical[:,0]*10, 200, facecolor='red', alpha=0.75, label='graphite')
plt.xlim(0,exp_r.max()*1.05)

plt.hist(zno_spherical[:, 0]*10, 200, facecolor='blue', alpha=0.75, label='ZnO')
plt.plot(k,profile/profile.max()*20,c='orange',label='diffractogram' );
plt.plot(k,np.abs(ctf)*3, label = f'ctf df = {defocus:.1f}nm')
#plt.hist(np.linalg.norm(p_ZnO, axis=1),200,facecolor='orange', alpha=0.75, label='hexagonal ZnO')
plt.legend()
plt.xlabel('reciprocal lattice vector [1/nm]')
plt.ylabel('number of reflections')

plt.ylim(0,9)

plt.show()

-6.588371940729947
[ 3.34695425  3.32376075  3.24472325 ... -6.58839474 -6.58837951
 -6.58837194]
Scherzer defocus is -85.8 nm

Another Image¶

Go back to Load an atomic resolution image

Try the same procedure with image p1-hr4-ZnOonGraphite.dm3 in the TEMdata directory, which is taken with a different defocus (which one?)

Conclusion:¶

We see that we can calibrate an image with sub pm precission if the lattice of a reference crystal is resolved. The original image scale and distortion was accurate within less than 2%. A quite respectable result, if one takes into account that the objective stigmators and the objecitve focus vary from image to image. These electron optic values will change the magnification of course.

Up Chapter 2: Diffraction
Back: Image Processing
Next: Z-Contrast Imaging
List of Content: Front

Image Analysis