diff --git a/py4DSTEM/process/__init__.py b/py4DSTEM/process/__init__.py index 0df11ef01..6f0019019 100644 --- a/py4DSTEM/process/__init__.py +++ b/py4DSTEM/process/__init__.py @@ -1,7 +1,6 @@ from py4DSTEM.process.polar import PolarDatacube from py4DSTEM.process.strain.strain import StrainMap -from py4DSTEM.process import latticevectors from py4DSTEM.process import phase from py4DSTEM.process import calibration from py4DSTEM.process import utils diff --git a/py4DSTEM/process/latticevectors/__init__.py b/py4DSTEM/process/latticevectors/__init__.py deleted file mode 100644 index cda4f91e5..000000000 --- a/py4DSTEM/process/latticevectors/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from py4DSTEM.process.latticevectors.initialguess import * -from py4DSTEM.process.latticevectors.index import * -from py4DSTEM.process.latticevectors.fit import * diff --git a/py4DSTEM/process/latticevectors/fit.py b/py4DSTEM/process/latticevectors/fit.py deleted file mode 100644 index d36b10bca..000000000 --- a/py4DSTEM/process/latticevectors/fit.py +++ /dev/null @@ -1,71 +0,0 @@ -# Functions for fitting lattice vectors to measured Bragg peak positions - -import numpy as np -from numpy.linalg import lstsq - -from emdfile import tqdmnd, PointList, PointListArray -from py4DSTEM.data import RealSlice - - -def fit_lattice_vectors_masked(braggpeaks, mask, x0=0, y0=0, minNumPeaks=5): - """ - Fits lattice vectors g1,g2 to each diffraction pattern in braggpeaks corresponding - to a scan position for which mask==True. - - Args: - braggpeaks (PointList): A 6 coordinate PointList containing the data to fit. - Coords are 'qx','qy' (the bragg peak positions), 'intensity' (used as a - weighting factor when fitting), 'h','k' (indexing). May optionally also - contain 'index_mask' (bool), indicating which peaks have been successfully - indixed and should be used. - mask (boolean array): real space shaped (R_Nx,R_Ny); fit lattice vectors where - mask is True - x0 (float): x-coord of the origin - y0 (float): y-coord of the origin - minNumPeaks (int): if there are fewer than minNumPeaks peaks found in braggpeaks - which can be indexed, return None for all return parameters - - Returns: - (RealSlice): A RealSlice ``g1g2map`` containing the following 8 arrays: - - * ``g1g2_map.get_slice('x0')`` x-coord of the origin of the best fit lattice - * ``g1g2_map.get_slice('y0')`` y-coord of the origin - * ``g1g2_map.get_slice('g1x')`` x-coord of the first lattice vector - * ``g1g2_map.get_slice('g1y')`` y-coord of the first lattice vector - * ``g1g2_map.get_slice('g2x')`` x-coord of the second lattice vector - * ``g1g2_map.get_slice('g2y')`` y-coord of the second lattice vector - * ``g1g2_map.get_slice('error')`` the fit error - * ``g1g2_map.get_slice('mask')`` 1 for successful fits, 0 for unsuccessful - fits - """ - assert isinstance(braggpeaks, PointListArray) - assert np.all( - [name in braggpeaks.dtype.names for name in ("qx", "qy", "intensity")] - ) - - # Make RealSlice to contain outputs - slicelabels = ("x0", "y0", "g1x", "g1y", "g2x", "g2y", "error", "mask") - g1g2_map = RealSlice( - data=np.zeros((braggpeaks.shape[0], braggpeaks.shape[1], 8)), - slicelabels=slicelabels, - name="g1g2_map", - ) - - # Fit lattice vectors - for Rx, Ry in tqdmnd(braggpeaks.shape[0], braggpeaks.shape[1]): - if mask[Rx, Ry]: - braggpeaks_curr = braggpeaks.get_pointlist(Rx, Ry) - qx0, qy0, g1x, g1y, g2x, g2y, error = fit_lattice_vectors( - braggpeaks_curr, x0, y0, minNumPeaks - ) - # Store data - if g1x is not None: - g1g2_map.get_slice("x0").data[Rx, Ry] = qx0 - g1g2_map.get_slice("y0").data[Rx, Ry] = qx0 - g1g2_map.get_slice("g1x").data[Rx, Ry] = g1x - g1g2_map.get_slice("g1y").data[Rx, Ry] = g1y - g1g2_map.get_slice("g2x").data[Rx, Ry] = g2x - g1g2_map.get_slice("g2y").data[Rx, Ry] = g2y - g1g2_map.get_slice("error").data[Rx, Ry] = error - g1g2_map.get_slice("mask").data[Rx, Ry] = 1 - return g1g2_map diff --git a/py4DSTEM/process/latticevectors/index.py b/py4DSTEM/process/latticevectors/index.py deleted file mode 100644 index 03cdf07ce..000000000 --- a/py4DSTEM/process/latticevectors/index.py +++ /dev/null @@ -1,147 +0,0 @@ -# Functions for indexing the Bragg directions - -import numpy as np -from numpy.linalg import lstsq - -from emdfile import tqdmnd, PointList, PointListArray - - -def get_selected_lattice_vectors(gx, gy, i0, i1, i2): - """ - From a set of reciprocal lattice points (gx,gy), and indices in those arrays which - specify the center beam, the first basis lattice vector, and the second basis lattice - vector, computes and returns the lattice vectors g1 and g2. - - Args: - gx (1d array): the reciprocal lattice points x-coords - gy (1d array): the reciprocal lattice points y-coords - i0 (int): index in the (gx,gy) arrays specifying the center beam - i1 (int): index in the (gx,gy) arrays specifying the first basis lattice vector - i2 (int): index in the (gx,gy) arrays specifying the second basis lattice vector - - Returns: - (2-tuple of 2-tuples) A 2-tuple containing - - * **g1**: *(2-tuple)* the first lattice vector, (g1x,g1y) - * **g2**: *(2-tuple)* the second lattice vector, (g2x,g2y) - """ - for i in (i0, i1, i2): - assert isinstance(i, (int, np.integer)) - g1x = gx[i1] - gx[i0] - g1y = gy[i1] - gy[i0] - g2x = gx[i2] - gx[i0] - g2y = gy[i2] - gy[i0] - return (g1x, g1y), (g2x, g2y) - - -def generate_lattice(ux, uy, vx, vy, x0, y0, Q_Nx, Q_Ny, h_max=None, k_max=None): - """ - Returns a full reciprocal lattice stretching to the limits of the diffraction pattern - by making linear combinations of the lattice vectors up to (±h_max,±k_max). - - This can be useful when there are false peaks or missing peaks in the braggvectormap, - which can cause errors in the strain finding routines that rely on those peaks for - indexing. This allows us to create a reference lattice that has all combinations of - the lattice vectors all the way out to the edges of the frame, and excluding any - erroneous intermediate peaks. - - Args: - ux (float): x-coord of the u lattice vector - uy (float): y-coord of the u lattice vector - vx (float): x-coord of the v lattice vector - vy (float): y-coord of the v lattice vector - x0 (float): x-coord of the lattice origin - y0 (float): y-coord of the lattice origin - Q_Nx (int): diffraction pattern size in the x-direction - Q_Ny (int): diffraction pattern size in the y-direction - h_max, k_max (int): maximal indices for generating the lattice (the lattive is - always trimmed to fit inside the pattern so you can overestimate these, or - leave unspecified and they will be automatically found) - - Returns: - (PointList): A 4-coordinate PointList, ('qx','qy','h','k'), containing points - corresponding to linear combinations of the u and v vectors, with associated - indices - """ - - # Matrix of lattice vectors - beta = np.array([[ux, uy], [vx, vy]]) - - # If no max index is specified, (over)estimate based on image size - if (h_max is None) or (k_max is None): - (y, x) = np.mgrid[0:Q_Ny, 0:Q_Nx] - x = x - x0 - y = y - y0 - h_max = np.max(np.ceil(np.abs((x / ux, y / uy)))) - k_max = np.max(np.ceil(np.abs((x / vx, y / vy)))) - - (hlist, klist) = np.meshgrid( - np.arange(-h_max, h_max + 1), np.arange(-k_max, k_max + 1) - ) - - M_ideal = np.vstack((hlist.ravel(), klist.ravel())).T - ideal_peaks = np.matmul(M_ideal, beta) - - coords = [("qx", float), ("qy", float), ("h", int), ("k", int)] - - ideal_data = np.zeros(len(ideal_peaks[:, 0]), dtype=coords) - ideal_data["qx"] = ideal_peaks[:, 0] - ideal_data["qy"] = ideal_peaks[:, 1] - ideal_data["h"] = M_ideal[:, 0] - ideal_data["k"] = M_ideal[:, 1] - - ideal_lattice = PointList(data=ideal_data) - - # shift to the DP center - ideal_lattice.data["qx"] += x0 - ideal_lattice.data["qy"] += y0 - - # trim peaks outside the image - deletePeaks = ( - (ideal_lattice.data["qx"] > Q_Nx) - | (ideal_lattice.data["qx"] < 0) - | (ideal_lattice.data["qy"] > Q_Ny) - | (ideal_lattice.data["qy"] < 0) - ) - ideal_lattice.remove(deletePeaks) - - return ideal_lattice - - -def bragg_vector_intensity_map_by_index(braggpeaks, h, k, symmetric=False): - """ - Returns a correlation intensity map for an indexed (h,k) Bragg vector - Used to obtain a darkfield image corresponding to the (h,k) reflection - or a bightfield image when h=k=0 - - Args: - braggpeaks (PointListArray): must contain the coordinates 'h','k', and - 'intensity' - h, k (int): indices for the reflection to generate an intensity map from - symmetric (bool): if set to true, returns sum of intensity of (h,k), (-h,k), - (h,-k), (-h,-k) - - Returns: - (numpy array): a map of the intensity of the (h,k) Bragg vector correlation. - Same shape as the pointlistarray. - """ - assert isinstance(braggpeaks, PointListArray), "braggpeaks must be a PointListArray" - assert np.all([name in braggpeaks.dtype.names for name in ("h", "k", "intensity")]) - intensity_map = np.zeros(braggpeaks.shape, dtype=float) - - for Rx in range(braggpeaks.shape[0]): - for Ry in range(braggpeaks.shape[1]): - pl = braggpeaks.get_pointlist(Rx, Ry) - if pl.length > 0: - if symmetric: - matches = np.logical_and( - np.abs(pl.data["h"]) == np.abs(h), - np.abs(pl.data["k"]) == np.abs(k), - ) - else: - matches = np.logical_and(pl.data["h"] == h, pl.data["k"] == k) - - if len(matches) > 0: - intensity_map[Rx, Ry] = np.sum(pl.data["intensity"][matches]) - - return intensity_map diff --git a/py4DSTEM/process/latticevectors/initialguess.py b/py4DSTEM/process/latticevectors/initialguess.py deleted file mode 100644 index d8054143f..000000000 --- a/py4DSTEM/process/latticevectors/initialguess.py +++ /dev/null @@ -1,229 +0,0 @@ -# Obtain an initial guess at the lattice vectors - -import numpy as np -from scipy.ndimage import gaussian_filter -from skimage.transform import radon - -from py4DSTEM.process.utils import get_maxima_1D - - -def get_radon_scores( - braggvectormap, - mask=None, - N_angles=200, - sigma=2, - minSpacing=2, - minRelativeIntensity=0.05, -): - """ - Calculates a score function, score(angle), representing the likelihood that angle is - a principle lattice direction of the lattice in braggvectormap. - - The procedure is as follows: - If mask is not None, ignore any data in braggvectormap where mask is False. Useful - for removing the unscattered beam, which can dominate the results. - Take the Radon transform of the (masked) Bragg vector map. - For each angle, get the corresponding slice of the sinogram, and calculate its score. - If we let R_theta(r) be the sinogram slice at angle theta, and where r is the - sinogram position coordinate, then the score of the slice is given by - score(theta) = sum_i(R_theta(r_i)) / N_i - Here, r_i are the positions r of all local maxima in R_theta(r), and N_i is the - number of such maxima. Thus the score is large when there are few maxima which are - high intensity. - - Args: - braggvectormap (ndarray): the Bragg vector map - mask (ndarray of bools): ignore data in braggvectormap wherever mask==False - N_angles (int): the number of angles at which to calculate the score - sigma (float): smoothing parameter for local maximum identification - minSpacing (float): if two maxima are found in a radon slice closer than - minSpacing, the dimmer of the two is removed - minRelativeIntensity (float): maxima in each radon slice dimmer than - minRelativeIntensity compared to the most intense maximum are removed - - Returns: - (3-tuple) A 3-tuple containing: - - * **scores**: *(ndarray, len N_angles, floats)* the scores for each angle - * **thetas**: *(ndarray, len N_angles, floats)* the angles, in radians - * **sinogram**: *(ndarray)* the radon transform of braggvectormap*mask - """ - # Get sinogram - thetas = np.linspace(0, 180, N_angles) - if mask is not None: - sinogram = radon(braggvectormap * mask, theta=thetas, circle=False) - else: - sinogram = radon(braggvectormap, theta=thetas, circle=False) - - # Get scores - N_maxima = np.empty_like(thetas) - total_intensity = np.empty_like(thetas) - for i in range(len(thetas)): - theta = thetas[i] - - # Get radon transform slice - ind = np.argmin(np.abs(thetas - theta)) - sinogram_theta = sinogram[:, ind] - sinogram_theta = gaussian_filter(sinogram_theta, 2) - - # Get maxima - maxima = get_maxima_1D(sinogram_theta, sigma, minSpacing, minRelativeIntensity) - - # Calculate metrics - N_maxima[i] = len(maxima) - total_intensity[i] = np.sum(sinogram_theta[maxima]) - scores = total_intensity / N_maxima - - return scores, np.radians(thetas), sinogram - - -def get_lattice_directions_from_scores( - thetas, scores, sigma=2, minSpacing=2, minRelativeIntensity=0.05, index1=0, index2=0 -): - """ - Get the lattice directions from the scores of the radon transform slices. - - Args: - thetas (ndarray): the angles, in radians - scores (ndarray): the scores - sigma (float): gaussian blur for local maxima identification - minSpacing (float): minimum spacing for local maxima identification - minRelativeIntensity (float): minumum intensity, relative to the brightest - maximum, for local maxima identification - index1 (int): specifies which local maximum to use for the first lattice - direction, in order of maximum intensity - index2 (int): specifies the local maximum for the second lattice direction - - Returns: - (2-tuple) A 2-tuple containing: - - * **theta1**: *(float)* the first lattice direction, in radians - * **theta2**: *(float)* the second lattice direction, in radians - """ - assert len(thetas) == len(scores), "Size of thetas and scores must match" - - # Get first lattice direction - maxima1 = get_maxima_1D( - scores, sigma, minSpacing, minRelativeIntensity - ) # Get maxima - thetas_max1 = thetas[maxima1] - scores_max1 = scores[maxima1] - dtype = np.dtype( - [("thetas", thetas.dtype), ("scores", scores.dtype)] - ) # Sort by intensity - ar_structured = np.empty(len(thetas_max1), dtype=dtype) - ar_structured["thetas"] = thetas_max1 - ar_structured["scores"] = scores_max1 - ar_structured = np.sort(ar_structured, order="scores")[::-1] - theta1 = ar_structured["thetas"][index1] # Get direction 1 - - # Apply sin**2 damping - scores_damped = scores * np.sin(thetas - theta1) ** 2 - - # Get second lattice direction - maxima2 = get_maxima_1D( - scores_damped, sigma, minSpacing, minRelativeIntensity - ) # Get maxima - thetas_max2 = thetas[maxima2] - scores_max2 = scores[maxima2] - dtype = np.dtype( - [("thetas", thetas.dtype), ("scores", scores.dtype)] - ) # Sort by intensity - ar_structured = np.empty(len(thetas_max2), dtype=dtype) - ar_structured["thetas"] = thetas_max2 - ar_structured["scores"] = scores_max2 - ar_structured = np.sort(ar_structured, order="scores")[::-1] - theta2 = ar_structured["thetas"][index2] # Get direction 2 - - return theta1, theta2 - - -def get_lattice_vector_lengths( - u_theta, - v_theta, - thetas, - sinogram, - spacing_thresh=1.5, - sigma=1, - minSpacing=2, - minRelativeIntensity=0.1, -): - """ - Gets the lengths of the two lattice vectors from their angles and the sinogram. - - First, finds the spacing between peaks in the sinogram slices projected down the u- - and v- directions, u_proj and v_proj. Then, finds the lengths by taking:: - - |u| = v_proj/sin(u_theta-v_theta) - |v| = u_proj/sin(u_theta-v_theta) - - The most important thresholds for this function are spacing_thresh, which discards - any detected spacing between adjacent radon projection peaks which deviate from the - median spacing by more than this fraction, and minRelativeIntensity, which discards - detected maxima (from which spacings are then calculated) below this threshold - relative to the brightest maximum. - - Args: - u_theta (float): the angle of u, in radians - v_theta (float): the angle of v, in radians - thetas (ndarray): the angles corresponding to the sinogram - sinogram (ndarray): the sinogram - spacing_thresh (float): ignores spacings which are greater than spacing_thresh - times the median spacing - sigma (float): gaussian blur for local maxima identification - minSpacing (float): minimum spacing for local maxima identification - minRelativeIntensity (float): minumum intensity, relative to the brightest - maximum, for local maxima identification - - Returns: - (2-tuple) A 2-tuple containing: - - * **u_length**: *(float)* the length of u, in pixels - * **v_length**: *(float)* the length of v, in pixels - """ - assert ( - len(thetas) == sinogram.shape[1] - ), "thetas must corresponding to the number of sinogram projection directions." - - # Get u projected spacing - ind = np.argmin(np.abs(thetas - u_theta)) - sinogram_slice = sinogram[:, ind] - maxima = get_maxima_1D(sinogram_slice, sigma, minSpacing, minRelativeIntensity) - spacings = np.sort(np.arange(sinogram_slice.shape[0])[maxima]) - spacings = spacings[1:] - spacings[:-1] - mask = ( - np.array( - [ - max(i, np.median(spacings)) / min(i, np.median(spacings)) - for i in spacings - ] - ) - < spacing_thresh - ) - spacings = spacings[mask] - u_projected_spacing = np.mean(spacings) - - # Get v projected spacing - ind = np.argmin(np.abs(thetas - v_theta)) - sinogram_slice = sinogram[:, ind] - maxima = get_maxima_1D(sinogram_slice, sigma, minSpacing, minRelativeIntensity) - spacings = np.sort(np.arange(sinogram_slice.shape[0])[maxima]) - spacings = spacings[1:] - spacings[:-1] - mask = ( - np.array( - [ - max(i, np.median(spacings)) / min(i, np.median(spacings)) - for i in spacings - ] - ) - < spacing_thresh - ) - spacings = spacings[mask] - v_projected_spacing = np.mean(spacings) - - # Get u and v lengths - sin_uv = np.sin(np.abs(u_theta - v_theta)) - u_length = v_projected_spacing / sin_uv - v_length = u_projected_spacing / sin_uv - - return u_length, v_length