We'll start by importing the usual packages
# Import numpy and fft functions directly
import numpy as np
from numpy.fft import fft2, ifft2, fftshift
# Set some default parameters for the plots
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [16, 8]
plt.rcParams['figure.frameon'] = False
plt.rcParams['image.cmap'] = 'gray'
save_figures = False # Use this to save figures to png files
The data from our microscope is in the proprietary .dm3
format from Gatan, the maker of the microscope. Usually this would be opened in something like Gatan Digital Micrograph, but luckily the OpenNCEM project helps us open it in Python (look here if using Matlab)
from ncempy.io.dm import dmReader
# Load dm3 files for a particle and reference
dmdata1 = dmReader('data/h05 +30sat.dm3')
dmdata1r = dmReader('data/h06r +30sat.dm3')
# Extract image data
im1 = dmdata1['data']
im1r = dmdata1r['data']
plt.imshow(im1)
if save_figures:
plt.savefig('01-particle1.png', dpi=200, bbox_inches = 'tight')
There is a wealth of metadata available about the microscope conditions embedded in the .dm3
files, but the most relevant to us is simply the scale, i.e. how many nanometers does each pixel correspond to in the image?
print('Scale from metadata: 1 px = %s %s' %(dmdata1['pixelSize'][0], dmdata1['pixelUnit'][0]))
# Redefine to nm and correct for Lorentz lens factor 10 error
pxsize = dmdata1['pixelSize'][0] * 1000/10
pixel_unit = 'nm'
print('Corrected scale: 1 px = %.4f %s' %(pxsize, pixel_unit))
print('Field of view: %.2f nm x %.2f nm' %(im1.shape[1]*pxsize, im1.shape[0]*pxsize))
# Uncomment below to see all tags and their values in the .dm3 file
#from ncempy.io.dm import fileDM
#dmdata = fileDM(p1_path)
#[('%s: %s' % (key, dmdata.allTags[key])) for key in dmdata.allTags]
Scale from metadata: 1 px = 0.0037131761 µm Corrected scale: 1 px = 0.3713 nm Field of view: 760.46 nm x 760.46 nm
In the fourier transformed image we see the two distinct sidebands. fftshift
is used to center the main peak in the image, as it would otherwise be split in the corners.
# Transform image, ready for visualization
im1_fft = fftshift(fft2(im1))
im1_fft_abslog = np.log(np.abs(im1_fft))
# Adjust colormap to increase contrast
vmin_fft = np.amin(im1_fft_abslog)*3
vmax_fft = np.amax(im1_fft_abslog)*0.8
plt.imshow(im1_fft_abslog, vmax=vmax_fft, vmin=vmin_fft)
plt.xlim([600, 1448])
plt.ylim([1448, 600])
if save_figures:
plt.savefig('02-particle1-fft.png', dpi=200, bbox_inches = 'tight')
Note that since the fourier image consists of complex numbers, we visualize the absolute value of the image, i.e. $\sqrt{\mathrm{Re}^2+\mathrm{Im}^2}$ using numpy.abs()
. The central peak is much more intense than the sidebands, so we also calculate the logarithm of the image. In this way the lower values get more pronounced and easier to distinguish.
To select a sideband, we find coordinates of the peaks in the image using peak_local_max
from the scikit-image
package as in this example. The central peak will always be the most intense followed by the two sidebands. We just find the second strongest peak and check if it's a positive or negative phase change, otherwise we take the third strongest.
We make this into a function as we'll do it a few times
from scipy.ndimage import maximum_filter
from skimage.feature import peak_local_max
def peak_coordinate(image_fft, peak_number=0):
# Find the peaks
im = np.abs(image_fft)
image_max = maximum_filter(im, size=200, mode='constant')
coordinates = peak_local_max(im, min_distance=200)
# Get values of peaks and sort, get coordinates of second largest peak
peaks = [image_fft[coordinate[0], coordinate[1]] for coordinate in coordinates]
peaks_idx = np.argsort(peaks)
coordinates_sorted = [coordinates[peak_idx] for peak_idx in peaks_idx]
sideband_idx = coordinates_sorted[peak_number]
# Choose sideband on the right side of the image
if sideband_idx[1] < image_fft.shape[1]//2:
sideband_idx[0] -= 2*(sideband_idx[0]-image_fft.shape[0]//2)
sideband_idx[1] -= 2*(sideband_idx[1]-image_fft.shape[1]//2)
return sideband_idx
We'll also make a function to extract sidebands from an FFT image
def extract_sideband(image_fft, width=128):
sideband_idx = peak_coordinate(np.abs(image_fft))
width = width//2
x1 = sideband_idx[1]-width
x2 = sideband_idx[1]+width
y1 = sideband_idx[0]-width
y2 = sideband_idx[0]+width
sideband_fft = image_fft[y1:y2, x1:x2]
return sideband_fft
Let's extract a sideband and see it up close
sideband_size = 200
sideband1_fft = extract_sideband(im1_fft, width=sideband_size)
plt.imshow(np.log(np.abs(sideband1_fft)))
if save_figures:
plt.savefig('03-particle1-sideand.png', dpi=200, bbox_inches = 'tight')
We can now calculate the amplitude and phase of the image with numpy's abs
and angle
methods
# Calculate amplitude and phase of image
sideband1 = ifft2(fftshift(sideband1_fft))
sideband1_abs = np.abs(sideband1)
sideband1_phase = np.angle(sideband1)
plt.title('Absolute value of sideband')
plt.imshow(sideband1_abs)
if save_figures:
plt.savefig('04a-p1-amp.png', dpi=200, bbox_inches = 'tight')
plt.title('Phase of sideband')
plt.imshow(sideband1_phase)
if save_figures:
plt.savefig('04b-p1-phase.png', dpi=200, bbox_inches = 'tight')
There are two important things to notice here. First of all the numpy
implementation of np.angle
uses the arctan2
function, a numerical implementation of $\tan^{-1}$. There is a nice description of the reasons in the Wikipedia article on atan2. For us, this means that the angle, or phase, we calculate has values as $-\pi<\phi\leq\pi$.
Second, and much more obvious is: the phase image contains sudden jumps in intensity. This happens when the phase "wraps around" with a period of $2 \pi$ (or 360 degrees). In other words the phase angle $\phi=0$ is equivalent to $\phi=2\pi$. Thus we need to "unwrap" our phase image. We will use the scikit-image phase unwrapping implementation
from skimage.restoration import unwrap_phase
sideband1_phase = unwrap_phase(sideband1_phase)
plt.imshow(sideband1_phase)
plt.title('Unwrapped phase of sideband')
if save_figures:
plt.savefig('05-p1-unwrapped-phase.png', dpi=200, bbox_inches = 'tight')
There is some significant background variations in the phase shift across the image due to imperfections in the biprism, e.g. the bottom left corner is quite dark while the bottom right is lighter. This can be negated with a reference image taken somewhere with a constant background, like an area with an even thickness of carbon substrate:
im1r_fft = fftshift(fft2(im1r))
sideband1r_fft = extract_sideband(im1r_fft, width=sideband_size)
fig, (ax1, ax2, ax3) = plt.subplots(1,3)
ax1.imshow(im1r)
ax2.imshow(np.abs(np.log(im1r_fft)))
ax3.imshow(np.abs(np.log(sideband1r_fft)))
if save_figures:
fig.savefig('06-p1r-extraction.png', dpi=200, bbox_inches = 'tight')
This is the reference phase shift we correct with
sideband1r = ifft2(fftshift(sideband1r_fft))
sideband1r_phase = unwrap_phase(np.angle(sideband1r))
plt.title('Background reference phase image')
plt.imshow(sideband1r_phase)
if save_figures:
plt.savefig('06-p1r-phase.png', dpi=200, bbox_inches = 'tight')
Now we subtract the reference from our image
phi1 = sideband1_phase-sideband1r_phase
plt.title('Corrected phase image of particle')
plt.imshow(phi1)
if save_figures:
plt.savefig('08a-phi1-corrected.png', dpi=200, bbox_inches = 'tight')
So far so good! We have now extracted a phase image of a magnetic nanoparticle. In order to extract just the magnetic contribution $\phi_m$ to the phase shift we need an image of the same particle where the particle's magnetization, and thus $\phi_m$, is reversed as described in the blog post:
$\phi = \phi_e+\phi_m - (\phi_e - \phi_m) = 2\phi_m$
The magnetic reversal was done by applying a strong magnetic field inside the microscope in the opposite direction of the previous used to magnetize the particle in the first place, using one of the magnetic lenses. We load in the image and subtract a reference phase image just like before:
# Load data
dmdata2 = dmReader('data/h19 -30sat.dm3')
dmdata2r = dmReader('data/h17r.dm3')
im2 = dmdata2['data']
im2r = dmdata2r['data']
# Calculate FFT and extract sideband for image
im2_fft = fftshift(fft2(im2))
sideband2_fft = extract_sideband(im2_fft, width=sideband_size)
sideband2 = ifft2(fftshift(sideband2_fft))
sideband2_phase = unwrap_phase(np.angle(sideband2))
# Calculate FFT and extract sideband for reference image
im2r_fft = fftshift(fft2(im2r))
sideband2r_fft = extract_sideband(im2r_fft, width=sideband_size)
sideband2r = ifft2(fftshift(sideband2r_fft))
sideband2r_phase = unwrap_phase(np.angle(sideband2r))
# Subtract reference phase shift
phi2 = sideband2_phase-sideband2r_phase
plt.title('Corrected phase image of particle with reversed magnetization')
plt.imshow(phi2)
if save_figures:
plt.savefig('08b-phi2-corrected.png', dpi=200, bbox_inches = 'tight')
As expected it looks very similar to the first phase image as $\phi_e$ is dominating the signal. Note that the particle is slightly shifted relative to the first image and subtrating the two directly would not help much. We need to figure out the relative displacement of the two by matching the images.
We will use a numeric method phase cross correlation to find the translation difference between the two images. To start with, we apply a bandpass filter using gaussian blur to remove high and low frequencies, as this can mess up the matching
from skimage.filters import window, gaussian, difference_of_gaussians
from skimage.registration import phase_cross_correlation
from skimage import transform
# Bandpass with Gaussian blur
sigma1 = 1
sigma2 = 3
phi1_filt = difference_of_gaussians(phi1, sigma1, sigma2)
phi2_filt = difference_of_gaussians(phi2, sigma1, sigma2)
fig, ax = plt.subplots(1,2)
ax[0].imshow(phi1_filt)
ax[1].imshow(phi2_filt)
if save_figures:
fig.savefig('09a-phi12-blurred.png', dpi=200, bbox_inches = 'tight')
In the earlier FFT images you might have noticed horizontal and vertical streaks. These come from the fact that our images are finite with an abrupt cutoff at the edges, which can lead to artifacts (sometimes termed spectral leakage). The borders of the above phase images also show large positive and negative values in some places around the edges, which is related to this effect.
These artifacts are detrimental to our image matching, so to avoid this we apply a 2D window function with radial symmetry, in our case a Hann window, before trying to match the images.
# Apply Hann window to avoid border issues
window_type = 'hann'
window_image = window(window_type, phi1_filt.shape)
phi1_windowed = phi1_filt*window_image
phi2_windowed = phi2_filt*window_image
fig, ax = plt.subplots(1,2)
ax[0].imshow(window_image)
ax[1].imshow(phi1_windowed)
if save_figures:
fig.savefig('09b-phi-windowing.png', dpi=200, bbox_inches = 'tight')
Since the images are very similar with negligible rotation we can calculate the shift directly from the phase correlation. Knowing this shift we can ensure the images have a common center and easily subtract one from the other
# Calculate the shift using phase cross correlation on filtered images
shifts, error, phasediff = phase_cross_correlation(phi1_filt, phi2_filt, upsample_factor=100)
# Note that matrices are indexed as rows ("y") at [0] and columns ("x") at [1]
dx = shifts[1]
dy = shifts[0]
print('Translation: (dx, dy) = (%.2f, %.2f)' % (dx, dy))
# Translate the second unfiltered phase image for matching on top of the first
tform = transform.EuclideanTransform(translation=(dx, dy))
phi2_tf = transform.warp(phi2, tform.inverse)
# Subtract unfiltered phase images, remember division by 2!
phi_m_raw = (phi1-phi2_tf)/2
plt.title('Magnetic phase shift')
plt.imshow(phi_m_raw)
if save_figures:
plt.savefig('10-phim-raw.png', dpi=200, bbox_inches = 'tight')
Translation: (dx, dy) = (15.04, 13.14)
The magnetic phase shift has revealed itself! Note how there is a lighter region on one side of the particle and dark on the other. We'll enhance it by cutting out the relevant part and filter high frequency noise with a simple gaussian blur
def crop_from_shifts(image, dx, dy, border=0):
dx, dy = int(dx), int(dy)
image = image[:, dx:-1] if dx >= 0 else image[:, 0:dx]
image = image[dy:-1, :] if dy >= 0 else image[0:dy, :]
border = int(border)
if border > 0:
image = image[border:-border, border:-border]
return image
Now we can plot the cleaned up image for a nice view of the phase image. We also add a scalebar with matplotlib-scalebar
from matplotlib_scalebar.scalebar import ScaleBar
# Scale original pixel size to match resolution of phase image and add to scalebar
pxsize_phase = im1.shape[0]/sideband_size*pxsize
scalebar = ScaleBar(pxsize_phase, pixel_unit, length_fraction=0.25)
# Blur before cropping to minimize edge effects
phi_m_filt = gaussian(phi_m_raw, 3)
# Crop non-overlapping part and a bit of border
crop_border = 5
phi_m = crop_from_shifts(phi_m_filt, dx, dy, border=crop_border)
fig, ax = plt.subplots()
ax.axis('off')
ax.imshow(phi_m)
ax.add_artist(scalebar)
plt.title('Magnetic phase shift')
if save_figures:
fig.savefig('10-phim-filtered.png', dpi=200, bbox_inches = 'tight')
A profile plot across the particle to shows the value of the phase change. Now is also a good time to ensure the axes on the image match the physical size of the picture
from skimage.measure import profile_line
# Extract profile plot across particle
x1, y1 = 20, 110
x2, y2 = 160, 60
dist = np.sqrt((x2-x1)**2 + (y2-y1)**2)*pxsize_phase
profile = profile_line(phi_m, (y1, x1), (y2, x2), linewidth=3, mode='reflect')
profile_dist = np.linspace(0, dist, len(profile))
phi_min, phi_max = np.min(profile), np.max(profile)
delta_phi = phi_max-phi_min
print('Phase change across particle: %.2f radians' %(delta_phi))
fig, ax = plt.subplots(1,2)
ax[0].axis('off')
ax[0].imshow(phi_m)
ax[0].plot((x1, x2), (y1, y2), color='r', linewidth=4)
ax[0].scatter(x1, y1, color='r', marker='o', s=100)
scalebar = ScaleBar(pxsize_phase, pixel_unit, length_fraction=0.25)
ax[0].add_artist(scalebar)
ax[1].plot(profile_dist, profile, color='r')
ax[1].set_xlabel('Distance [nm]')
ax[1].set_ylabel('Phase $\phi$ [radian]')
ax[1].set_box_aspect(1)
if save_figures:
fig.savefig('11-phim-profile.png', dpi=200, bbox_inches = 'tight')
Phase change across particle: 1.05 radians
We can now visualize the magnetic field lines by displaying contour lines of the phase image (i.e. lines where the phase has the same value). By using a divergent colormap 'RdBu' which goes from red to blue with white for central values we can also easily distinguish values above and below the central values.
It can also be useful to overlay contourlines onto the bright field image to illustrate the magnetic field lines coming from the particle
# Apply a bit stronger blur to enhance contour lines
phi_m_extrafilt = crop_from_shifts(gaussian(phi_m_raw, 4), dx, dy, border=crop_border)
fig, ax = plt.subplots(1,2)
# Plot phase with divergent colormap, overlay contour lines
ax[0].imshow(phi_m, cmap='RdBu')
ax[0].contour(phi_m_extrafilt, 15, colors='k', alpha=0.5)
ax[0].axis('off')
scalebar = ScaleBar(pxsize_phase, pixel_unit, length_fraction=0.25)
ax[0].add_artist(scalebar)
# Crop BF image with same parameters as the phase image was cropped to get same center
bf_image = crop_from_shifts(sideband1_abs, dx, dy, border=crop_border)
# Overlay contours onto bright field image
ax[1].imshow(gaussian(bf_image, 1))
ax[1].contour(phi_m_extrafilt, 12, colors='white', linewidths=2, alpha=0.8)
ax[1].axis('off')
scalebar = ScaleBar(pxsize_phase, pixel_unit, length_fraction=0.25)
ax[1].add_artist(scalebar)
if save_figures:
fig.savefig('11-phim-contour-overlay.png', dpi=200, bbox_inches = 'tight')
Now that is looking like proper magnetic field lines!
We can also overlay the contours on our high resolution original image, we just need to careful to rescale and crop everything the same way. We will also overlay a scalebar for a nice, final image.
from skimage.transform import rescale
# Increase resolution of phase image
scaling = im1.shape[0]/sideband_size
phi_m_upscaled = rescale(phi_m_extrafilt, scaling)
# Crop original image the same way as the sideband images were, correcting for scale
im1_crop = crop_from_shifts(im1, dx*scaling, dy*scaling, border=crop_border*scaling)
im1_filt = gaussian(im1_crop, 3)
fig, ax = plt.subplots()
ax.axis("off")
ax.imshow(im1_filt)
ax.contour(phi_m_upscaled, 12, colors='white', linewidths=2, alpha=0.6)
scalebar = ScaleBar(pxsize, pixel_unit, length_fraction=0.25)
ax.add_artist(scalebar)
if save_figures:
plt.savefig('12-hires-contour-overlay.png', dpi=200, bbox_inches = 'tight')