Commit 75ff3686 authored by Mukherjee, Debangshu's avatar Mukherjee, Debangshu
Browse files

Accurate DPC potential calculation, units may still be off

parent 289f6e74
This diff is collapsed.
......@@ -88,6 +88,10 @@ class atomic_dpc(object):
Shiteng Zhao, Philipp M. Pelz, Edward S. Barnard et al. "py4DSTEM: a software package for
multimodal analysis of four-dimensional scanning transmission electron microscopy datasets."
arXiv preprint arXiv:2003.09523 (2020).
.. [3] Ishizuka, Akimitsu, Masaaki Oka, Takehito Seki, Naoya Shibata,
and Kazuo Ishizuka. "Boundary-artifact-free determination of
potential distribution from differential phase contrast signals."
Microscopy 66, no. 6 (2017): 397-405.
"""
def __init__(self, Data_4D, Data_ADF, calib_pm, voltage, aperture):
......@@ -113,6 +117,12 @@ class atomic_dpc(object):
self.planck = 6.62607004 * (10 ** (-34))
self.epsilon0 = 8.85418782 * (10 ** (-12))
self.e_charge = (-1) * 1.60217662 * (10 ** (-19))
e_mass = 9.109383 * (10 ** (-31))
c = 299792458
self.sigma = (
(2 * np.pi / (self.wavelength * self.voltage))
* ((e_mass * (c ** 2)) + (self.e_charge * self.voltage))
) / ((2 * e_mass * (c ** 2)) + (self.e_charge * self.voltage))
def show_ADF_BF(self, imsize=(20, 10)):
"""
......@@ -181,7 +191,7 @@ class atomic_dpc(object):
plt.gca().add_artist(scalebar)
plt.axis("off")
def initial_dpc(self, imsize=(30, 17)):
def initial_dpc(self, imsize=(30, 17), normalize=True):
"""
"""
......@@ -199,6 +209,9 @@ class atomic_dpc(object):
self.XCom[yy[ii], xx[ii]] = self.inverse * (
(np.sum(np.multiply(pp, pattern)) / np.sum(pattern)) - self.beam_x
)
if normalize:
self.YCom = self.YCom - np.mean(self.YCom)
self.XCom = self.XCom - np.mean(self.XCom)
vm = (np.amax(np.abs(np.concatenate((self.XCom, self.YCom), axis=1)))) / (
10 ** 9
......@@ -384,7 +397,7 @@ class atomic_dpc(object):
ax1 = plt.subplot(gs[0:15, 0:15])
ax2 = plt.subplot(gs[15:17, :])
ax1.imshow(self.charge, vmin=-cm, vmax=cm, cmap="seismic")
ax1.imshow(self.charge, vmin=-cm, vmax=cm, cmap="RdBu_r")
scalebar = mpss.ScaleBar(self.calib / 1000, "nm")
scalebar.location = "lower right"
scalebar.box_alpha = 1
......@@ -420,7 +433,10 @@ class atomic_dpc(object):
def show_potential(self, imsize=(10, 10)):
fontsize = int(np.amax(np.asarray(imsize)))
self.potential = st.dpc.integrate_dpc(self.e_fieldX, self.e_fieldY)
self.phase = st.dpc.integrate_dpc(
self.XComC * self.wavelength, self.YComC * self.wavelength
)
self.potential = self.phase / self.sigma
cm = np.amax(np.abs(self.potential))
plt.figure(figsize=imsize)
plt.imshow(self.potential, vmin=-cm, vmax=cm, cmap="RdBu_r")
......
......@@ -4,6 +4,7 @@ import numpy as np
import numba
import warnings
import stemtool as st
import pyfftw.interfaces.numpy_fft as pfft
def cart2pol(xx, yy):
......@@ -55,28 +56,40 @@ def data_rotator(cbed_pattern, rotangle, xcenter, ycenter, data_radius):
return rotated_cbed
def integrate_dpc(xshift, yshift, fourier_calibration=1):
def integrate_dpc(xshift, yshift, padf=4, lP=0.5, hP=100, stepsize=0.2, iter_count=100):
"""
Integrate DPC shifts using Fourier transforms and
preventing edge effects
Parameters
----------
xshift: ndarray
Beam shift in the X dimension
yshift: ndarray
Beam shift in the X dimensions
fourier_calibration: float
Pixel size of the Fourier space
xshift: ndarray
Beam shift in the X dimension
yshift: ndarray
Beam shift in the X dimensions
padf: int, optional
padding factor for accurate FFT,
default is 4
lP: float, optional
low pass filter, default is 0.5
hP: float, optional
high pass filter, default is 100
stepsize: float, optional
fraction of phase differnce to update every
iteration. Default is 0.5. This is a dynamic
factor, and is reduced if the error starts
increasing
iter_count: int, optional
Number of iterations to run. Default is 100
Returns
-------
integrand: ndarray
Integrated DPC
phase_final: ndarray
Phase of the matrix that leads to the displacement
Notes
-----
This is based on two ideas - noniterative complex plane
This is based on two ideas - iterative complex plane
integration and antisymmetric mirror integration. First
two antisymmetric matrices are generated for each of the
x shift and y shifts. Then they are integrated in Fourier
......@@ -86,74 +99,65 @@ def integrate_dpc(xshift, yshift, fourier_calibration=1):
References
----------
Bon, Pierre, Serge Monneret, and Benoit Wattellier.
"Noniterative boundary-artifact-free wavefront
reconstruction from its derivatives." Applied optics
51, no. 23 (2012): 5698-5704.
.. [1] Ishizuka, Akimitsu, Masaaki Oka, Takehito Seki, Naoya Shibata,
and Kazuo Ishizuka. "Boundary-artifact-free determination of
potential distribution from differential phase contrast signals."
Microscopy 66, no. 6 (2017): 397-405.
:Authors:
Debangshu Mukherjee <mukherjeed@ornl.gov>
"""
# Initialize matrices
size_array = np.asarray(np.shape(xshift))
x_mirrored = np.zeros(2 * size_array, dtype=np.float64)
y_mirrored = np.zeros(2 * size_array, dtype=np.float64)
# Generate antisymmetric X arrays
x_mirrored[0 : size_array[0], 0 : size_array[1]] = np.fliplr(np.flipud(0 - xshift))
x_mirrored[0 : size_array[0], size_array[1] : (2 * size_array[1])] = np.fliplr(
0 - xshift
)
x_mirrored[size_array[0] : (2 * size_array[0]), 0 : size_array[1]] = np.flipud(
xshift
)
x_mirrored[
size_array[0] : (2 * size_array[0]), size_array[1] : (2 * size_array[1])
] = xshift
# Generate antisymmetric Y arrays
y_mirrored[0 : size_array[0], 0 : size_array[1]] = np.fliplr(np.flipud(0 - yshift))
y_mirrored[0 : size_array[0], size_array[1] : (2 * size_array[1])] = np.fliplr(
yshift
)
y_mirrored[size_array[0] : (2 * size_array[0]), 0 : size_array[1]] = np.flipud(
0 - yshift
)
y_mirrored[
size_array[0] : (2 * size_array[0]), size_array[1] : (2 * size_array[1])
] = yshift
# Calculated Fourier transform of antisymmetric matrices
x_mirr_ft = np.fft.fft2(x_mirrored)
y_mirr_ft = np.fft.fft2(y_mirrored)
# Calculated inverse Fourier space calibration
qx = np.mean(
np.diff(
(np.arange(-size_array[1], size_array[1], 1))
/ (2 * fourier_calibration * size_array[1])
imshape = np.asarray(xshift.shape)
padshape = (imshape * padf).astype(int)
qx = np.fft.fftfreq(padshape[0])
qy = np.fft.rfftfreq(padshape[1])
qr2 = qx[:, None] ** 2 + qy[None, :] ** 2
denominator = qr2 + hP + ((qr2 ** 2) * lP)
_ = np.seterr(divide="ignore")
denominator = 1.0 / denominator
denominator[0, 0] = 0
_ = np.seterr(divide="warn")
f = 1j * 0.25 * stepsize
qxOperator = f * qx[:, None] * denominator
qyOperator = f * qy[None, :] * denominator
padded_phase = np.zeros(padshape)
update = np.zeros_like(padded_phase)
dx = np.zeros_like(padded_phase)
dy = np.zeros_like(padded_phase)
error = np.zeros(iter_count)
mask = np.zeros_like(padded_phase, dtype=bool)
mask[
int(0.5 * (padshape[0] - imshape[0])) : int(0.5 * (padshape[0] + imshape[0])),
int(0.5 * (padshape[1] - imshape[1])) : int(0.5 * (padshape[1] + imshape[1])),
] = True
maskInv = mask == False
for ii in range(iter_count):
dx[mask] -= xshift.ravel()
dy[mask] -= yshift.ravel()
dx[maskInv] = 0
dy[maskInv] = 0
update = pfft.irfft2(pfft.rfft2(dx) * qxOperator + pfft.rfft2(dy) * qyOperator)
padded_phase += scnd.gaussian_filter((stepsize * update), 1)
dx = (
np.roll(padded_phase, (-1, 0), axis=(0, 1))
- np.roll(padded_phase, (1, 0), axis=(0, 1))
) / 2.0
dy = (
np.roll(padded_phase, (0, -1), axis=(0, 1))
- np.roll(padded_phase, (0, 1), axis=(0, 1))
) / 2.0
xDiff = dx[mask] - xshift.ravel()
yDiff = dy[mask] - yshift.ravel()
error[ii] = np.sqrt(
np.mean((xDiff - np.mean(xDiff)) ** 2 + (yDiff - np.mean(yDiff)) ** 2)
)
)
qy = np.mean(
np.diff(
(np.arange(-size_array[0], size_array[0], 1))
/ (2 * fourier_calibration * size_array[0])
)
)
# Calculate mirrored CPM integrand
mirr_ft = (x_mirr_ft + ((1j) * y_mirr_ft)) / (qx + ((1j) * qy))
mirr_int = np.fft.ifft2(mirr_ft)
# Select integrand from antisymmetric matrix
integrand = np.abs(
mirr_int[
size_array[0] : (2 * size_array[0]), size_array[1] : (2 * size_array[1])
]
)
return integrand
if ii > 0:
if error[ii] > error[ii - 1]:
stepsize /= 2
phase_final = np.reshape(padded_phase[mask], imshape)
return phase_final
def potential_dpc(x_dpc, y_dpc, angle=0):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment