From d1933b7e4ad84e3f6af0ec4a9e1706af97f3bf26 Mon Sep 17 00:00:00 2001 From: Stephanie Ribet Date: Mon, 24 Jun 2024 19:07:06 -0700 Subject: [PATCH 01/15] first pass at grain clustering --- py4DSTEM/process/polar/polar_analysis.py | 270 ++++++++++++++++++++++- py4DSTEM/process/polar/polar_datacube.py | 3 + 2 files changed, 269 insertions(+), 4 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index 78a95c24a..80ec4bd93 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1,14 +1,16 @@ # Analysis scripts for amorphous 4D-STEM data using polar transformations. +import sys +import time from typing import Tuple -import numpy as np + import matplotlib.pyplot as plt -from scipy.optimize import curve_fit +import numpy as np +from emdfile import tqdmnd from scipy.ndimage import gaussian_filter +from scipy.optimize import curve_fit from sklearn.decomposition import PCA -from emdfile import tqdmnd - def calculate_radial_statistics( self, @@ -1016,3 +1018,263 @@ def background_pca( ax.plot(coef_pca[:, 0], coef_pca[:, 1]) return im_pca, coef_pca + + +def cluster_grains( + self, + orientation_histogram=None, + radial_range=None, + threshold_add=0.05, + threshold_grow=0.2, + angle_tolerance_deg=5, + plot_grain_clusters=True, + **kwargs +): + """ + Cluster grains from a specific radial bin + + Parameters + -------- + orientation_histogram: np.ndarray + Must be a 3D array of size Rx, Ry, and polar_shape[0] + Can be used to bypass step to calculate orientation histogram for + grain clustering + radial_range: int + Size (2,) array for bins + threshold_add: float + Minimum signal required for a probe position to initialize a cluster. + threshold_grow: float + Minimum signal required for a probe position to be added to a cluster. + angle_tolerance_deg: float + Rotation rolerance for clustering grains. + plot_grain_clusters: bool + If True, plots clusters. **kwargs passed to `plot_grain_clusters` + + """ + + if orientation_histogram is None: + if hasattr(self, "peaks_refine"): + use_refined_peaks = True + else: + use_refined_peaks = False + sig = np.squeeze( + self.make_orientation_histogram( + [ + radial_range, + ], + use_refined_peaks=use_refined_peaks, + upsample_factor=1, + progress_bar=False, + ) + ) + else: + sig = np.squeeze(orientation_histogram.copy()) + assert sig.shape == ( + self.peaks.shape[0], + self.peaks.shape[1], + self.polar_shape[0], + ), "orientation_histogram must have an `upsample_factor` of 1" + + sig_init = sig.copy() + mark = sig >= threshold_grow + sig[np.logical_not(mark)] = 0 + + phi = np.arange(sig.shape[-1]) / sig.shape[-1] * np.pi + + # init + self.cluster_sizes = np.array((), dtype="int") + self.cluster_sig = np.array(()) + self.cluster_inds = [] + inds_all = np.zeros_like(sig, dtype="int") + inds_all.ravel()[:] = np.arange(inds_all.size) + + # Tolerance + tol = np.deg2rad(angle_tolerance_deg) + + # Main loop + search = True + while search is True: + inds_grain = np.argmax(sig) + val = sig.ravel()[inds_grain] + + if val < threshold_add: + search = False + + else: + # progressbar + comp = 1 - np.mean(np.max(mark, axis=2)) + update_progress(comp) + + # Start cluster + x, y, z = np.unravel_index(inds_grain, sig.shape) + mark[x, y, z] = False + sig[x, y, z] = 0 + phi_cluster = phi[z] + + # Neighbors to search + xr = np.clip(x + np.arange(-1, 2, dtype="int"), 0, sig.shape[0] - 1) + yr = np.clip(y + np.arange(-1, 2, dtype="int"), 0, sig.shape[1] - 1) + inds_cand = inds_all[xr[:, None], yr[None], :].ravel() + inds_cand = np.delete(inds_cand, mark.ravel()[inds_cand] == False) + + if inds_cand.size == 0: + grow = False + else: + grow = True + + # grow the cluster + while grow is True: + inds_new = np.array((), dtype="int") + + keep = np.zeros(inds_cand.size, dtype="bool") + for a0 in range(inds_cand.size): + xc, yc, zc = np.unravel_index(inds_cand[a0], sig.shape) + + phi_test = phi[zc] + dphi = ( + np.mod(phi_cluster - phi_test + np.pi / 2.0, np.pi) + - np.pi / 2.0 + ) + + if np.abs(dphi) < tol: + keep[a0] = True + + sig[xc, yc, zc] = 0 + mark[xc, yc, zc] = False + + xr = np.clip( + xc + np.arange(-1, 2, dtype="int"), 0, sig.shape[0] - 1 + ) + yr = np.clip( + yc + np.arange(-1, 2, dtype="int"), 0, sig.shape[1] - 1 + ) + inds_add = inds_all[xr[:, None], yr[None], :].ravel() + inds_new = np.append(inds_new, inds_add) + + inds_grain = np.append(inds_grain, inds_cand[keep]) + inds_cand = np.unique( + np.delete(inds_new, mark.ravel()[inds_new] == False) + ) + + if inds_cand.size == 0: + grow = False + + # convert grain to x,y coordinates, add = list + xg, yg, zg = np.unravel_index(inds_grain, sig.shape) + xyg = np.unique(np.vstack((xg, yg)), axis=1) + sig_mean = np.mean(sig_init.ravel()[inds_grain]) + self.cluster_sizes = np.append(self.cluster_sizes, xyg.shape[1]) + self.cluster_sig = np.append(self.cluster_sig, sig_mean) + self.cluster_inds.append(xyg) + + # finish progressbar + update_progress(1) + + if plot_grain_clusters: + self.plot_grain_clusters(**kwargs) + + return + + +def plot_grain_clusters( + self, + area_min=None, + area_max=None, + area_step=1, + weight_intensity=False, + pixel_area=1.0, + pixel_area_units="px^2", + figsize=(8, 6), + returnfig=False, +): + """ + Plot the cluster sizes + + Parameters + -------- + area_min: int (optional) + Min area bin in pixels + area_max: int (optional) + Max area bin in pixels + area_step: int (optional) + Step size of the histogram bin + weight_intensity: bool + Weight histogram by the peak intensity. + pixel_area: float + Size of pixel area unit square + pixel_area_units: string + Units of the pixel area + figsize: tuple + Size of the figure panel + returnfig: bool + Setting this to true returns the figure and axis handles + + Returns + -------- + fig, ax (optional) + Figure and axes handles + + """ + + if area_max is None: + area_max = np.max(self.cluster_sizes) + area = np.arange(0, area_max, area_step) + if area_min is None: + sub = self.cluster_sizes.astype("int") < area_max + else: + sub = np.logical_and( + self.cluster_sizes.astype("int") >= area_min, + self.cluster_sizes.astype("int") < area_max, + ) + if weight_intensity: + hist = np.bincount( + self.cluster_sizes[sub] // area_step, + weights=self.cluster_sig[sub], + minlength=area.shape[0], + ) + else: + hist = np.bincount( + self.cluster_sizes[sub] // area_step, + minlength=area.shape[0], + ) + + # plotting + fig, ax = plt.subplots(figsize=figsize) + ax.bar( + area * pixel_area, + hist, + width=0.8 * pixel_area * area_step, + ) + ax.set_xlim((0, area_max * pixel_area)) + ax.set_xlabel("Grain Area [" + pixel_area_units + "]") + if weight_intensity: + ax.set_ylabel("Total Signal [arb. units]") + else: + ax.set_ylabel("Number of Grains") + + if returnfig: + return fig, ax + + +# Progressbar taken from stackexchange: +# https://stackoverflow.com/questions/3160699/python-progress-bar +def update_progress(progress): + barLength = 60 # Modify this to change the length of the progress bar + status = "" + if isinstance(progress, int): + progress = float(progress) + if not isinstance(progress, float): + progress = 0 + status = "error: progress var must be float\r\n" + if progress < 0: + progress = 0 + status = "Halt...\r\n" + if progress >= 1: + progress = 1 + status = "Done...\r\n" + block = int(round(barLength * progress)) + text = "\rPercent: [{0}] {1}% {2}".format( + "#" * block + "-" * (barLength - block), np.round(progress * 100, 4), status + ) + sys.stdout.write(text) + sys.stdout.flush() diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index de3c463b8..869e62072 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -105,7 +105,10 @@ def __init__( plot_reduced_pdf, plot_pdf, background_pca, + cluster_grains, + plot_grain_clusters, ) + from py4DSTEM.process.polar.polar_peaks import ( find_peaks_single_pattern, find_peaks, From 263e0f43055103bc312dc7f191bfb4b38a3ecebf Mon Sep 17 00:00:00 2001 From: Stephanie Ribet Date: Wed, 26 Jun 2024 07:30:23 +0100 Subject: [PATCH 02/15] return calc --- py4DSTEM/process/polar/polar_analysis.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index 80ec4bd93..ec78c505e 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1028,6 +1028,7 @@ def cluster_grains( threshold_grow=0.2, angle_tolerance_deg=5, plot_grain_clusters=True, + returncalc=False, **kwargs ): """ @@ -1049,7 +1050,8 @@ def cluster_grains( Rotation rolerance for clustering grains. plot_grain_clusters: bool If True, plots clusters. **kwargs passed to `plot_grain_clusters` - + return_calc: bool + If True, return cluster_sizes, cluster_sig, and cluster_inds """ if orientation_histogram is None: @@ -1173,7 +1175,8 @@ def cluster_grains( if plot_grain_clusters: self.plot_grain_clusters(**kwargs) - return + if returncalc: + return self.cluster_sizes, self.cluster_sig, self.cluster_inds def plot_grain_clusters( From f66001623cc5728a0b41117569cbec9cdfde1e08 Mon Sep 17 00:00:00 2001 From: Colin Date: Sun, 30 Jun 2024 15:08:44 -0700 Subject: [PATCH 03/15] updates to polar peak detection / plotting --- py4DSTEM/process/polar/polar_peaks.py | 59 ++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 5 deletions(-) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 535ae7143..58194cbdf 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -35,11 +35,14 @@ def find_peaks_single_pattern( remove_masked_peaks=False, scale_sigma_annular=0.5, scale_sigma_radial=0.25, + refine_subpixel = True, return_background=False, plot_result=True, + plot_smoothed_image=False, plot_power_scale=1.0, plot_scale_size=10.0, figsize=(12, 6), + figax = None, returnfig=False, **kwargs ): @@ -82,16 +85,22 @@ def find_peaks_single_pattern( Scaling of the estimated annular standard deviation. scale_sigma_radial: float Scaling of the estimated radial standard deviation. + refine_subpixel: bool + Use parabolic fit to find subpixel positions. return_background: bool Return the background signal. - plot_result: - Plot the detector peaks + plot_result: bool + Plot the detected peaks. + plot_smoothed_image: bool + Plot the image after smoothing is applied. plot_power_scale: float Image intensity power law scaling. plot_scale_size: float Marker scaling in the plot. figsize: 2-tuple Size of the result plotting figure. + fig,ax: 2-tuple + Matplotlib figure and axes handles to plot result in. returnfig: bool Return the figure and axes handles. @@ -160,6 +169,8 @@ def find_peaks_single_pattern( im_polar_sm[sub] /= im_mask[sub] # Find local maxima + # with warnings.catch_warnings(): + # warnings.simplefilter("ignore") peaks = peak_local_max( im_polar_sm, num_peaks=num_peaks_max, @@ -246,6 +257,20 @@ def find_peaks_single_pattern( axis=0, ) + # If needed, refine peaks using parabolic fit + peaks = peaks.astype('float') + if refine_subpixel: + for a0 in range(peaks.shape[0]): + if peaks[a0,1] > 0 and peaks[a0,1] < im_polar.shape[1]-1: + im_crop = im_polar_sm[ + (peaks[a0,0] + np.arange(-1,2)).astype('int')[:,None] % im_polar.shape[0], + (peaks[a0,1] + np.arange(-1,2)).astype('int'), + ] + dx = (im_crop[2,1] - im_crop[0,1]) / (4*im_crop[1,1] - 2*im_crop[2,1] - 2*im_crop[0,1]) + dy = (im_crop[1,2] - im_crop[1,0]) / (4*im_crop[1,1] - 2*im_crop[1,2] - 2*im_crop[1,0]) + peaks[a0,0] += dx + peaks[a0,1] += dy + # combine peaks into one array peaks_all = np.column_stack((peaks, peaks_int, peaks_prom)) @@ -302,19 +327,34 @@ def find_peaks_single_pattern( if plot_result: # init - im_plot = im_polar.copy() + if plot_smoothed_image: + im_plot = im_polar_sm.copy() + else: + im_plot = im_polar.copy() im_plot = np.maximum(im_plot, 0) ** plot_power_scale t = np.linspace(0, 2 * np.pi, 180 + 1) ct = np.cos(t) st = np.sin(t) - fig, ax = plt.subplots(figsize=figsize) + if figax is None: + fig, ax = plt.subplots(figsize=figsize) + else: + fig = figax[0] + ax = figax[1] cmap = kwargs.pop("cmap", "gray") vmax = kwargs.pop("vmax", 1) vmin = kwargs.pop("vmin", 0) - show(im_plot, figax=(fig, ax), cmap=cmap, vmax=vmax, vmin=vmin, **kwargs) + show( + im_plot, + figax=(fig, ax), + cmap=cmap, + vmax=vmax, + vmin=vmin, + **kwargs, + ) + ax.set_aspect('auto') # peaks ax.scatter( @@ -723,6 +763,10 @@ def plot_radial_peaks( color="r", linewidth=2, ) + ax.yaxis.get_ticklocs(minor=True) + ax.minorticks_on() + ax.yaxis.set_tick_params(which='minor',left=False) + ax.xaxis.grid(True) ax.set_xlim((q_bins[0], q_bins[-1])) if q_pixel_units: ax.set_xlabel( @@ -740,6 +784,11 @@ def plot_radial_peaks( ) if not label_y_axis: ax.tick_params(left=False, labelleft=False) + # ax.tick_params( + # axis='y', + # which='minor', + # bottom=True, + # ) if returnfig: return fig, ax From 25cc5ac57fb8635c841a213268239305f8687213 Mon Sep 17 00:00:00 2001 From: Colin Date: Sun, 30 Jun 2024 19:23:50 -0700 Subject: [PATCH 04/15] adding grain diameter plots --- py4DSTEM/process/polar/polar_analysis.py | 121 +++++++++++++++++++++-- py4DSTEM/process/polar/polar_datacube.py | 1 + 2 files changed, 115 insertions(+), 7 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index ec78c505e..3b01dc11f 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1185,6 +1185,7 @@ def plot_grain_clusters( area_max=None, area_step=1, weight_intensity=False, + weight_area=True, pixel_area=1.0, pixel_area_units="px^2", figsize=(8, 6), @@ -1203,6 +1204,8 @@ def plot_grain_clusters( Step size of the histogram bin weight_intensity: bool Weight histogram by the peak intensity. + weight_area: bool + Weight histogram by the area of each bin. pixel_area: float Size of pixel area unit square pixel_area_units: string @@ -1243,15 +1246,119 @@ def plot_grain_clusters( # plotting fig, ax = plt.subplots(figsize=figsize) - ax.bar( - area * pixel_area, - hist, - width=0.8 * pixel_area * area_step, - ) + if weight_area: + ax.bar( + area * pixel_area, + hist * area, + width=0.8 * pixel_area * area_step, + ) + else: + ax.bar( + area * pixel_area, + hist, + width=0.8 * pixel_area * area_step, + ) ax.set_xlim((0, area_max * pixel_area)) - ax.set_xlabel("Grain Area [" + pixel_area_units + "]") + ax.set_xlabel("Grain Area (" + pixel_area_units + ")") + if weight_intensity: + ax.set_ylabel("Total Signal (arb. units)") + elif weight_area: + ax.set_ylabel("Area-Weighted Signal (arb. units)") + else: + ax.set_ylabel("Number of Grains") + + if returnfig: + return fig, ax + + +def plot_grain_clusters_diameter( + self, + diameter_min=None, + diameter_max=None, + diameter_step=1, + weight_intensity=False, + weight_diameter=True, + pixel_area=1.0, + pixel_area_units="px", + figsize=(8, 6), + returnfig=False, +): + """ + Plot the cluster sizes + + Parameters + -------- + diameter_min: int (optional) + Min area bin in pixels + diameter_max: int (optional) + Max area bin in pixels + diameter_step: int (optional) + Step size of the histogram bin + weight_intensity: bool + Weight histogram by the peak intensity. + weight_diameter: bool + Weight histogram by the area of each bin. + pixel_area: float + Size of pixel area unit square + pixel_area_units: string + Units of the pixel area + figsize: tuple + Size of the figure panel + returnfig: bool + Setting this to true returns the figure and axis handles + + Returns + -------- + fig, ax (optional) + Figure and axes handles + + """ + + cluster_diam = 0.5 * np.sqrt(self.cluster_sizes.astype('float') / np.pi) + + + if diameter_max is None: + diameter_max = np.max(cluster_diam) + diameter = np.arange(0, diameter_max, diameter_step) + if diameter_min is None: + sub = cluster_diam < diameter_max + else: + sub = np.logical_and( + cluster_diam >= diameter_min, + cluster_diam < diameter_max, + ) + if weight_intensity: + hist = np.bincount( + np.round(cluster_diam[sub] / diameter_step).astype('int'), + weights=self.cluster_sig[sub], + minlength=diameter.shape[0], + ) + else: + hist = np.bincount( + np.round(cluster_diam[sub] / diameter_step).astype('int'), + minlength=diameter.shape[0], + ) + + # plotting + fig, ax = plt.subplots(figsize=figsize) + if weight_diameter: + ax.bar( + diameter * pixel_area, + hist * diameter, + width=0.8 * pixel_area * diameter_step, + ) + else: + ax.bar( + diameter * pixel_area, + hist, + width=0.8 * pixel_area * diameter_step, + ) + ax.set_xlim((0, diameter_max * pixel_area)) + ax.set_xlabel("Domain Size (" + pixel_area_units + ")") if weight_intensity: - ax.set_ylabel("Total Signal [arb. units]") + ax.set_ylabel("Total Signal (arb. units)") + # elif weight_area: + # ax.set_ylabel("Area-Weighted Signal (arb. units)") else: ax.set_ylabel("Number of Grains") diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index 869e62072..e249d2877 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -107,6 +107,7 @@ def __init__( background_pca, cluster_grains, plot_grain_clusters, + plot_grain_clusters_diameter, ) from py4DSTEM.process.polar.polar_peaks import ( From d7fe1ea7e4b9e8daefbfa8480f7817a2829aacf3 Mon Sep 17 00:00:00 2001 From: Colin Date: Mon, 1 Jul 2024 09:40:55 -0700 Subject: [PATCH 05/15] Adding grain plotting --- py4DSTEM/process/polar/polar_analysis.py | 152 ++++++++++++++++++++++- py4DSTEM/process/polar/polar_datacube.py | 3 +- 2 files changed, 151 insertions(+), 4 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index 3b01dc11f..bb349ebfa 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -7,7 +7,8 @@ import matplotlib.pyplot as plt import numpy as np from emdfile import tqdmnd -from scipy.ndimage import gaussian_filter +from scipy.ndimage import gaussian_filter, distance_transform_edt +from skimage.morphology import dilation, erosion from scipy.optimize import curve_fit from sklearn.decomposition import PCA @@ -1035,7 +1036,7 @@ def cluster_grains( Cluster grains from a specific radial bin Parameters - -------- + ---------- orientation_histogram: np.ndarray Must be a 3D array of size Rx, Ry, and polar_shape[0] Can be used to bypass step to calculate orientation histogram for @@ -1052,6 +1053,18 @@ def cluster_grains( If True, plots clusters. **kwargs passed to `plot_grain_clusters` return_calc: bool If True, return cluster_sizes, cluster_sig, and cluster_inds + + + Returns + ---------- + cluster_sizes: np.array + (N,) shape vector giving the size of N grains as the # of probe positions + cluster_sig: np.array + (N,) shape vector giving the mean signal of each of N grains. + cluster_inds: list of np.array + Length N list of probe positions for each grain, stored as arrays with + shape (2,M) for M probe positions where rows correspond to the (x,y) indices. + """ if orientation_histogram is None: @@ -1179,7 +1192,140 @@ def cluster_grains( return self.cluster_sizes, self.cluster_sig, self.cluster_inds -def plot_grain_clusters( +def plot_clusters( + self, + area_min=2, + outline_grains=True, + outline_thickness=1, + fill_grains=0.25, + smooth_grains=1.0, + cmap="viridis", + figsize=(8, 8), + progress_bar=True, + returncalc=True, + returnfig=False, +): + """ + Plot the clusters as an image. + + Parameters + -------- + area_min: int (optional) + Min cluster size to include, in units of probe positions. + outline_grains: bool (optional) + Set to True to draw grains with outlines + outline_thickness: int (optional) + Thickenss of the grain outline + fill_grains: float (optional) + Outlined grains are filled with this value in pixels. + smooth_grains: float (optional) + Grain boundaries are smoothed by this value in pixels. + figsize: tuple + Size of the figure panel + returncalc: bool + Return the grain image. + returnfig: bool, optional + Setting this to true returns the figure and axis handles + + Returns + -------- + im_plot: np.array + The plotting image + fig, ax (optional) + Figure and axes handles + + """ + + + # init + im_plot = np.zeros( + ( + self.data_raw.shape[0], + self.data_raw.shape[1], + ) + ) + im_grain = np.zeros( + ( + self.data_raw.shape[0], + self.data_raw.shape[1], + ), + dtype="bool", + ) + + # make plotting image + for a0, ry in tqdmnd( + self.cluster_sizes.shape[0]. + desc="Generating grain image", + unit=" grains", + disable=not progress_bar, + ): + if self.cluster_sizes[a0] >= area_min: + if outline_grains: + im_grain[:] = False + im_grain[ + self.cluster_inds[a0][0, :], + self.cluster_inds[a0][1, :], + ] = True + + im_dist = distance_transform_edt( + erosion( + np.invert(im_grain), footprint=np.ones((3, 3), dtype="bool") + ) + ) - distance_transform_edt(im_grain) + im_dist = gaussian_filter(im_dist, sigma=smooth_grains, mode="nearest") + im_add = np.exp(im_dist**2 / (-0.5 * outline_thickness**2)) + + if fill_grains > 0: + im_dist = distance_transform_edt( + erosion( + np.invert(im_grain), footprint=np.ones((3, 3), dtype="bool") + ) + ) + im_dist = gaussian_filter( + im_dist, sigma=smooth_grains, mode="nearest" + ) + im_add += fill_grains * np.exp( + im_dist**2 / (-0.5 * outline_thickness**2) + ) + + # im_add = 1 - np.exp( + # distance_transform_edt(im_grain)**2 \ + # / (-2*outline_thickness**2)) + im_plot += im_add + # im_plot = np.minimum(im_plot, im_add) + else: + # xg,yg = np.unravel_index(self.cluster_inds[a0], im_plot.shape) + im_grain[:] = False + im_grain[ + self.cluster_inds[a0][0, :], + self.cluster_inds[a0][1, :], + ] = True + im_plot += gaussian_filter( + im_grain.astype("float"), sigma=smooth_grains, mode="nearest" + ) + + # im_plot[ + # self.cluster_inds[a0][0,:], + # self.cluster_inds[a0][1,:], + # ] += 1 + + if outline_grains: + im_plot = np.clip(im_plot, 0, 2) + + # plotting + fig, ax = plt.subplots(figsize=figsize) + ax.imshow( + im_plot, + # vmin = -3, + # vmax = 3, + cmap=cmap, + ) + + if returncalc: + return im_plot + + +def plot_grain_clusters_area( self, area_min=None, area_max=None, diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index e249d2877..60cd2ca94 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -106,7 +106,8 @@ def __init__( plot_pdf, background_pca, cluster_grains, - plot_grain_clusters, + plot_clusters, + plot_grain_clusters_area, plot_grain_clusters_diameter, ) From 99078f0073fe24aff772f07da8d5e66fd173461c Mon Sep 17 00:00:00 2001 From: Colin Date: Mon, 1 Jul 2024 09:45:08 -0700 Subject: [PATCH 06/15] typo --- py4DSTEM/process/polar/polar_analysis.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index bb349ebfa..84bae3903 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1253,8 +1253,8 @@ def plot_clusters( ) # make plotting image - for a0, ry in tqdmnd( - self.cluster_sizes.shape[0]. + for a0 in tqdmnd( + self.cluster_sizes.shape[0], desc="Generating grain image", unit=" grains", disable=not progress_bar, @@ -1322,7 +1322,13 @@ def plot_clusters( ) if returncalc: - return im_plot + if returnfig: + return im_plot, fig, ax + else: + return im_plot + else: + if returnfig: + return fig, ax def plot_grain_clusters_area( From 084dd1bcd2a3fd5438b950b699ce4ad4a8a0e5ca Mon Sep 17 00:00:00 2001 From: Colin Date: Tue, 2 Jul 2024 17:21:30 -0700 Subject: [PATCH 07/15] More updates to the polar branch --- py4DSTEM/process/polar/polar_analysis.py | 250 ++++++++++++++++++++--- py4DSTEM/process/polar/polar_datacube.py | 5 +- py4DSTEM/process/polar/polar_peaks.py | 8 +- 3 files changed, 232 insertions(+), 31 deletions(-) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index 84bae3903..19e630d8d 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1028,6 +1028,7 @@ def cluster_grains( threshold_add=0.05, threshold_grow=0.2, angle_tolerance_deg=5, + distance_tolerance_px = 1, plot_grain_clusters=True, returncalc=False, **kwargs @@ -1048,7 +1049,9 @@ def cluster_grains( threshold_grow: float Minimum signal required for a probe position to be added to a cluster. angle_tolerance_deg: float - Rotation rolerance for clustering grains. + Angular rotation tolerance for clustering grains. + distance_tolerance_px: int + Distance tolerance for clustering grains. plot_grain_clusters: bool If True, plots clusters. **kwargs passed to `plot_grain_clusters` return_calc: bool @@ -1107,6 +1110,7 @@ def cluster_grains( tol = np.deg2rad(angle_tolerance_deg) # Main loop + vec_search = np.arange(-distance_tolerance_px, distance_tolerance_px+1, dtype="int") search = True while search is True: inds_grain = np.argmax(sig) @@ -1127,8 +1131,8 @@ def cluster_grains( phi_cluster = phi[z] # Neighbors to search - xr = np.clip(x + np.arange(-1, 2, dtype="int"), 0, sig.shape[0] - 1) - yr = np.clip(y + np.arange(-1, 2, dtype="int"), 0, sig.shape[1] - 1) + xr = np.clip(x + vec_search, 0, sig.shape[0] - 1) + yr = np.clip(y + vec_search, 0, sig.shape[1] - 1) inds_cand = inds_all[xr[:, None], yr[None], :].ravel() inds_cand = np.delete(inds_cand, mark.ravel()[inds_cand] == False) @@ -1198,6 +1202,7 @@ def plot_clusters( outline_grains=True, outline_thickness=1, fill_grains=0.25, + weight_by_area=False, smooth_grains=1.0, cmap="viridis", figsize=(8, 8), @@ -1206,24 +1211,26 @@ def plot_clusters( returnfig=False, ): """ - Plot the clusters as an image. + Plot the clusters / domains as an image. Parameters -------- area_min: int (optional) Min cluster size to include, in units of probe positions. outline_grains: bool (optional) - Set to True to draw grains with outlines + Set to True to draw domains with outlines outline_thickness: int (optional) - Thickenss of the grain outline + Thickenss of the domain outline fill_grains: float (optional) - Outlined grains are filled with this value in pixels. + Outlined domains are filled with this value in pixels. + weight_by_area: bool (optional) + Weight the domain fill and edges by each domain's area. smooth_grains: float (optional) - Grain boundaries are smoothed by this value in pixels. + Domain boundaries are smoothed by this value in pixels. figsize: tuple Size of the figure panel returncalc: bool - Return the grain image. + Return the domain image. returnfig: bool, optional Setting this to true returns the figure and axis handles @@ -1252,10 +1259,15 @@ def plot_clusters( dtype="bool", ) + if weight_by_area: + weights = self.cluster_sizes / np.max(self.cluster_sizes) * (1-fill_grains) + fill_grains + else: + weights = np.ones_like(self.cluster_sizes) + # make plotting image for a0 in tqdmnd( self.cluster_sizes.shape[0], - desc="Generating grain image", + desc="Generating domain image", unit=" grains", disable=not progress_bar, ): @@ -1291,7 +1303,7 @@ def plot_clusters( # im_add = 1 - np.exp( # distance_transform_edt(im_grain)**2 \ # / (-2*outline_thickness**2)) - im_plot += im_add + im_plot += im_add * weights[a0] # im_plot = np.minimum(im_plot, im_add) else: # xg,yg = np.unravel_index(self.cluster_inds[a0], im_plot.shape) @@ -1302,7 +1314,7 @@ def plot_clusters( ] = True im_plot += gaussian_filter( im_grain.astype("float"), sigma=smooth_grains, mode="nearest" - ) + ) * weights[a0] # im_plot[ # self.cluster_inds[a0][0,:], @@ -1331,7 +1343,7 @@ def plot_clusters( return fig, ax -def plot_grain_clusters_area( +def plot_clusters_area( self, area_min=None, area_max=None, @@ -1423,14 +1435,16 @@ def plot_grain_clusters_area( return fig, ax -def plot_grain_clusters_diameter( +def plot_clusters_diameter( self, + cluster_sizes=None, diameter_min=None, diameter_max=None, diameter_step=1, weight_intensity=False, - weight_diameter=True, - pixel_area=1.0, + weight_diameter=False, + weight_area=False, + pixel_size=1.0, pixel_area_units="px", figsize=(8, 6), returnfig=False, @@ -1440,6 +1454,8 @@ def plot_grain_clusters_diameter( Parameters -------- + cluster_sizes: np.array + Size in pixels^2 of all clusters. diameter_min: int (optional) Min area bin in pixels diameter_max: int (optional) @@ -1449,8 +1465,10 @@ def plot_grain_clusters_diameter( weight_intensity: bool Weight histogram by the peak intensity. weight_diameter: bool - Weight histogram by the area of each bin. - pixel_area: float + Weight histogram by the diameter of each cluster. + weight_diameter: bool + Weight histogram by the area of each cluster. + pixel_size: float Size of pixel area unit square pixel_area_units: string Units of the pixel area @@ -1466,9 +1484,14 @@ def plot_grain_clusters_diameter( """ - cluster_diam = 0.5 * np.sqrt(self.cluster_sizes.astype('float') / np.pi) + if cluster_sizes is None: + cluster_sizes = self.cluster_sizes.copy().astype('float') + cluster_diam = 0.5 * np.sqrt(cluster_sizes / np.pi) + # units + cluster_diam *= pixel_size + # subset if diameter_max is None: diameter_max = np.max(cluster_diam) diameter = np.arange(0, diameter_max, diameter_step) @@ -1479,12 +1502,26 @@ def plot_grain_clusters_diameter( cluster_diam >= diameter_min, cluster_diam < diameter_max, ) + + # histogram if weight_intensity: hist = np.bincount( np.round(cluster_diam[sub] / diameter_step).astype('int'), weights=self.cluster_sig[sub], minlength=diameter.shape[0], ) + elif weight_area: + hist = np.bincount( + np.round(cluster_diam[sub] / diameter_step).astype('int'), + weights=cluster_sizes[sub], + minlength=diameter.shape[0], + ) + elif weight_diameter: + hist = np.bincount( + np.round(cluster_diam[sub] / diameter_step).astype('int'), + weights=cluster_diam[sub], + minlength=diameter.shape[0], + ) else: hist = np.bincount( np.round(cluster_diam[sub] / diameter_step).astype('int'), @@ -1495,24 +1532,181 @@ def plot_grain_clusters_diameter( fig, ax = plt.subplots(figsize=figsize) if weight_diameter: ax.bar( - diameter * pixel_area, + diameter, hist * diameter, - width=0.8 * pixel_area * diameter_step, + width=0.8 * diameter_step, ) else: ax.bar( - diameter * pixel_area, + diameter, hist, - width=0.8 * pixel_area * diameter_step, + width=0.8 * diameter_step, ) - ax.set_xlim((0, diameter_max * pixel_area)) + ax.set_xlim((0, diameter_max)) ax.set_xlabel("Domain Size (" + pixel_area_units + ")") + if weight_intensity: - ax.set_ylabel("Total Signal (arb. units)") - # elif weight_area: - # ax.set_ylabel("Area-Weighted Signal (arb. units)") + ax.set_ylabel("Intensity-Weighted Domain Size") + elif weight_area: + ax.set_ylabel("Total Domain Area") + elif weight_diameter: + ax.set_ylabel("Total Domain Diameter") else: - ax.set_ylabel("Number of Grains") + ax.set_ylabel("Number of Domains") + + if weight_intensity or weight_area or weight_diameter: + ax.set_yticks([]) + + if returnfig: + return fig, ax + + + +def plot_clusters_max_length( + self, + cluster_inds=None, + length_min=None, + length_max=None, + length_step=1, + weight_intensity=False, + weight_length=False, + weight_area=False, + pixel_size=1.0, + pixel_area_units="px", + figsize=(8, 6), + returnfig=False, +): + """ + Plot the histogram of the max length of each cluster (over all angles). + + Parameters + -------- + cluster_inds: list of np.array + List of clusters, with the indices given for each cluster. + length_min: int (optional) + Min area bin in pixels + length_max: int (optional) + Max area bin in pixels + length_step: int (optional) + Step size of the histogram bin + weight_intensity: bool + Weight histogram by the peak intensity. + weight_length: bool + Weight histogram by the length of each cluster. + weight_diameter: bool + Weight histogram by the area of each cluster. + pixel_size: float + Size of pixel area unit square + pixel_area_units: string + Units of the pixel area + figsize: tuple + Size of the figure panel + returnfig: bool + Setting this to true returns the figure and axis handles + + Returns + -------- + fig, ax (optional) + Figure and axes handles + + """ + + if cluster_inds is None: + cluster_inds = self.cluster_inds.copy().astype('float') + + # init + num_clusters = len(cluster_inds) + cluster_sizes = np.zeros(num_clusters) + cluster_length = np.zeros(num_clusters) + t_all = np.linspace(0,np.pi,45,endpoint=False) + ct = np.cos(t_all) + st = np.sin(t_all) + + # calculate size and max lengths of each cluster + for a0 in range(num_clusters): + cluster_sizes[a0] = cluster_inds[a0].shape[1] + + x0 = cluster_inds[a0][0] + y0 = cluster_inds[a0][1] + + length_current_max = 0 + for a1 in range(t_all.size): + p = x0*ct[a1] + y0*st[a1] + length_current_max = np.maximum( + length_current_max, + np.max(p) - np.min(p) + 1, + ) + cluster_length[a0] = length_current_max + + + # units + cluster_length *= pixel_size + + # subset + if length_max is None: + length_max = np.max(cluster_length) + length = np.arange(0, length_max, length_step) + if length_min is None: + sub = cluster_length < length_max + else: + sub = np.logical_and( + cluster_length >= length_min, + cluster_length < length_max, + ) + + # histogram + if weight_intensity: + hist = np.bincount( + np.round(cluster_length[sub] / length_step).astype('int'), + weights=self.cluster_sig[sub], + minlength=length.shape[0], + ) + elif weight_area: + hist = np.bincount( + np.round(cluster_length[sub] / length_step).astype('int'), + weights=cluster_sizes[sub], + minlength=length.shape[0], + ) + elif weight_diameter: + hist = np.bincount( + np.round(cluster_length[sub] / length_step).astype('int'), + weights=cluster_diam[sub], + minlength=length.shape[0], + ) + else: + hist = np.bincount( + np.round(cluster_length[sub] / length_step).astype('int'), + minlength=length.shape[0], + ) + + # plotting + fig, ax = plt.subplots(figsize=figsize) + if weight_length: + ax.bar( + length, + hist * length, + width=0.8 * length_step, + ) + else: + ax.bar( + length, + hist, + width=0.8 * length_step, + ) + ax.set_xlim((0, length_max)) + ax.set_xlabel("Maximum Domain Length (" + pixel_area_units + ")") + + if weight_intensity: + ax.set_ylabel("Intensity-Weighted Domain Lengths") + elif weight_area: + ax.set_ylabel("Total Domain Area") + elif weight_diameter: + ax.set_ylabel("Total Domain Diameter") + else: + ax.set_ylabel("Number of Domains") + + if weight_intensity or weight_area or weight_length: + ax.set_yticks([]) if returnfig: return fig, ax diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index 60cd2ca94..87baba756 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -107,8 +107,9 @@ def __init__( background_pca, cluster_grains, plot_clusters, - plot_grain_clusters_area, - plot_grain_clusters_diameter, + plot_clusters_area, + plot_clusters_diameter, + plot_clusters_max_length, ) from py4DSTEM.process.polar.polar_peaks import ( diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 58194cbdf..cbf29ba7d 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -802,6 +802,7 @@ def model_radial_background( refine_model=True, plot_result=True, figsize=(8, 4), + returnfig = True, ): """ User provided radial background model, of the form: @@ -905,12 +906,17 @@ def background_model(q, *coefs): # plotting if plot_result: - self.plot_radial_background( + fig, ax = self.plot_radial_background( q_pixel_units=False, plot_background_model=True, figsize=figsize, + returnfig=returnfig, ) + if returnfig: + return fig,ax + + def refine_peaks( self, From 78b67d3ab6126d2efd4167cf4a747ca3f749fa1c Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 6 Jul 2024 08:51:47 -0700 Subject: [PATCH 08/15] updating docstring --- py4DSTEM/process/polar/polar_peaks.py | 56 ++++++++++++++++++--------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index cbf29ba7d..7bc9637d7 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -1263,25 +1263,43 @@ def make_orientation_histogram( NOTE - currently assumes two fold rotation symmetry TODO - add support for non two fold symmetry polardatacube - Args: - radial_ranges (np array): Size (N x 2) array for N radial bins, or (2,) for a single bin. - orientation_flip_sign (bool): Flip the direction of theta - orientation_offset_degrees (float): Offset for orientation angles - orientation_separate_bins (bool): whether to place multiple angles into multiple radial bins. - upsample_factor (float): Upsample factor - use_refined_peaks (float): Use refined peak positions - use_peak_sigma (float): Spread signal along annular direction using measured std. - theta_step_deg (float): Step size along annular direction in degrees - sigma_x (float): Smoothing in x direction before upsample - sigma_y (float): Smoothing in x direction before upsample - sigma_theta (float): Smoothing in annular direction (units of bins, periodic) - normalize_intensity_image (bool): Normalize to max peak intensity = 1, per image - normalize_intensity_stack (bool): Normalize to max peak intensity = 1, all images - progress_bar (bool): Enable progress bar - - Returns: - orient_hist (array): 4D array containing Bragg peak intensity histogram - [radial_bin x_probe y_probe theta] + Parameters + ---------- + radial_ranges (np array) + Size (N x 2) array for N radial bins, or (2,) for a single bin. + orientation_flip_sign (bool) + Flip the direction of theta + orientation_offset_degrees (float) + Offset for orientation angles. + This value is a rotation of Q space with respect to R space. + orientation_separate_bins (bool) + whether to place multiple angles into multiple radial bins. + upsample_factor (float) + Upsample factor + use_refined_peaks (float) + Use refined peak positions + use_peak_sigma (float) + Spread signal along annular direction using measured std. + theta_step_deg (float) + Step size along annular direction in degrees + sigma_x (float) + Smoothing in x direction before upsample + sigma_y (float) + Smoothing in x direction before upsample + sigma_theta (float) + Smoothing in annular direction (units of bins, periodic) + normalize_intensity_image (bool) + Normalize to max peak intensity = 1, per image + normalize_intensity_stack (bool) + Normalize to max peak intensity = 1, all images + progress_bar (bool) + Enable progress bar + + Returns + ---------- + orient_hist (array) + 4D array containing Bragg peak intensity histogram + [radial_bin x_probe y_probe theta] """ # coordinates From d8bd92e0d3cce845cbdc2042bfb2fa5933fc2d6c Mon Sep 17 00:00:00 2001 From: Colin Date: Sun, 7 Jul 2024 13:41:27 -0700 Subject: [PATCH 09/15] Fixing the orientation offset, angular legend plotting --- py4DSTEM/process/diffraction/flowlines.py | 47 +++++++++++++++-------- py4DSTEM/process/polar/polar_peaks.py | 9 +++-- 2 files changed, 35 insertions(+), 21 deletions(-) diff --git a/py4DSTEM/process/diffraction/flowlines.py b/py4DSTEM/process/diffraction/flowlines.py index 5af64f73e..177baf8b2 100644 --- a/py4DSTEM/process/diffraction/flowlines.py +++ b/py4DSTEM/process/diffraction/flowlines.py @@ -800,7 +800,7 @@ def make_flowline_rainbow_image( def make_flowline_rainbow_legend( im_size=np.array([256, 256]), sym_rotation_order=2, - theta_offset=0.0, + theta_offset_degrees=0.0, white_background=False, return_image=False, radial_range=np.array([0.45, 0.9]), @@ -810,26 +810,39 @@ def make_flowline_rainbow_legend( """ This function generates a legend for a the rainbow colored flowline maps, and returns it as an RGB image. - Args: - im_size (np.array): Size of legend image in pixels. - sym_rotation_order (int): rotational symmety for colouring - theta_offset (float): Offset the anglular coloring by this value in radians. - white_background (bool): For either color or greyscale output, switch to white background (from black). - return_image (bool): Return the image array. - radial_range (np.array): Inner and outer radius for the legend ring. - plot_legend (bool): Plot the generated legend. - figsize (tuple or list): Size of the plotted legend. + Parameters + ---------- + im_size (np.array): + Size of legend image in pixels. + sym_rotation_order (int): + rotational symmety for colouring + theta_offset_degrees (float): + Offset the anglular coloring by this value in degrees. + Rotation is Q with respect to R, in the positive (counter clockwise) direction. + white_background (bool): + For either color or greyscale output, switch to white background (from black). + return_image (bool): + Return the image array. + radial_range (np.array): + Inner and outer radius for the legend ring. + plot_legend (bool): + Plot the generated legend. + figsize (tuple or list): + Size of the plotted legend. - Returns: - im_legend (array): Image array for the legend. + Returns + ---------- + + im_legend (array): + Image array for the legend. """ # Coordinates x = np.linspace(-1, 1, im_size[0]) y = np.linspace(-1, 1, im_size[1]) - ya, xa = np.meshgrid(-y, x) + ya, xa = np.meshgrid(y, x) ra = np.sqrt(xa**2 + ya**2) - ta = np.arctan2(ya, xa) + theta_offset + ta = np.arctan2(ya, xa) + np.deg2rad(theta_offset_degrees) ta_sym = ta * sym_rotation_order # mask @@ -837,12 +850,12 @@ def make_flowline_rainbow_legend( # rgb image z = mask * np.exp(1j * ta_sym) - hue_start = -90 + # hue_offset = 0 amp = np.abs(z) vmin = np.min(amp) vmax = np.max(amp) - ph = np.angle(z, deg=1) + hue_start - h = (ph % 360) / 360 + ph = np.angle(z)# + hue_offset + h = np.mod(ph/(2*np.pi), 1) s = 0.85 * np.ones_like(h) v = (amp - vmin) / (vmax - vmin) im_legend = hsv_to_rgb(np.dstack((h, s, v))) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 7bc9637d7..8442c4950 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -1307,10 +1307,11 @@ def make_orientation_histogram( # Get angles from polardatacube theta = self.tt else: - theta = np.arange(0, 180, theta_step_deg) * np.pi / 180.0 - dtheta = theta[1] - theta[0] - dtheta_deg = dtheta * 180 / np.pi + theta = np.deg2rad(np.arange(0, 180, theta_step_deg)) num_theta_bins = np.size(theta) + dtheta = theta[1] - theta[0] + dtheta_deg = np.rad2deg(dtheta) + orientation_offset = np.deg2rad(orientation_offset_degrees) # Input bins radial_ranges = np.array(radial_ranges) @@ -1373,7 +1374,7 @@ def make_orientation_histogram( theta = self.peaks[rx, ry]["qt"][sub] * self._annular_step if orientation_flip_sign: theta *= -1 - theta += orientation_offset_degrees + theta += orientation_offset t = theta / dtheta From 5f124d1e2444b79ac09cdb7223fa8468d5a4dede Mon Sep 17 00:00:00 2001 From: Colin Date: Tue, 16 Jul 2024 18:40:18 -0700 Subject: [PATCH 10/15] Adding integrated crystalline and amorphous signals --- py4DSTEM/process/polar/polar_datacube.py | 2 + py4DSTEM/process/polar/polar_peaks.py | 190 +++++++++++++++++++++++ 2 files changed, 192 insertions(+) diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index 87baba756..05c80d2c7 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -120,6 +120,8 @@ def __init__( plot_radial_peaks, plot_radial_background, model_radial_background, + fit_crystal_amorphous_fraction, + plot_crystal_amorphous_fraction, make_orientation_histogram, ) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 8442c4950..31238580b 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -916,6 +916,196 @@ def background_model(q, *coefs): if returnfig: return fig,ax +def fit_crystal_amorphous_fraction( + self, + fitting_range_sigma = 2.0, + plot_result = True, + figsize = (4,4), + progress_bar = True, + ): + """ + Fit an amorphous halo and the crystalline peak signals from each + polar transformed image, in order to estimate the fraction of + crystalline and amorphous signal. + + Parameters + ---------- + + plot_result: bool + + progress_bar: bool + + + Returns + ---------- + + + """ + + # Basis for amorphous fitting function + num_rings = np.round((self.background_coefs.shape[0] - 3) / 3).astype("int") + num_bins = self.polar_shape[1] + basis = np.ones(( + num_bins, + num_rings + 2, + )) + + # init + self.signal_amorphous_peaks = np.zeros(( + num_rings, + self._datacube.shape[0], + self._datacube.shape[1], + )) + self.signal_crystal = np.zeros(( + num_rings, + self._datacube.shape[0], + self._datacube.shape[1], + )) + crystal_fitting_mask = np.zeros(( + num_rings, + num_bins, + ), dtype = 'bool') + + # background model + sigma_0 = self.background_coefs[2] + basis[:,1] = np.exp(self.qq**2/(-2.0*sigma_0**2)) + + # amorphous halos / rings + for a0 in range(num_rings): + ring_sigma = self.background_coefs[3 * a0 + 4] + ring_position = self.background_coefs[3 * a0 + 5] + basis[:,2+a0] = np.exp( + (self.qq - ring_position)**2 \ + /(-2.0*ring_sigma**2) + ) + + sub = np.logical_and( + self.qq > ring_position - ring_sigma*fitting_range_sigma, + self.qq <= ring_position + ring_sigma*fitting_range_sigma, + ) + crystal_fitting_mask[a0,sub] = True + + # main loop over probe positions + for rx, ry in tqdmnd( + # np.arange(60,100), + # np.arange(290,351), + self._datacube.shape[0], + self._datacube.shape[1], + desc="Refining peaks ", + unit=" probe positions", + disable=not progress_bar, + ): + # polar signal + im_polar = self.data[rx,ry] + im_polar_med = np.ma.median(im_polar, axis = 0) + + # fit amorphous background coefficients + coefs = np.linalg.lstsq( + basis, + im_polar_med, + rcond=None, + )[0] + + # background estimate + # im_polar_bg = basis @ coefs + im_polar_sub = im_polar - (basis @ coefs) + + # Output signals + self.signal_amorphous_peaks[:,rx,ry] = coefs[2:] + for a0 in range(num_rings): + self.signal_crystal[a0,rx,ry] = np.sum( + im_polar_sub[:,crystal_fitting_mask[a0]], + ) + + self.signal_amorphous_peaks = np.maximum(self.signal_amorphous_peaks,0.0) + self.signal_crystal = np.maximum(self.signal_crystal,0.0) + + # convert amorphous signal from peaks to integrated intensity + self.signal_amorphous = self.signal_amorphous_peaks.copy() + for a0 in range(num_rings): + ring_sigma = self.background_coefs[3 * a0 + 4] + self.signal_amorphous[a0] *= ring_sigma * 2.5069 + + # fractional crystal signal + sig_sum = self.signal_crystal + self.signal_amorphous + sub = sig_sum > 0.0 + self.signal_fractional_crystal = self.signal_crystal.copy() + self.signal_fractional_crystal[sub] /= sig_sum[sub] + + # plotting + if plot_result: + fig,ax = plt.subplots(figsize=figsize) + ax.imshow( + self.signal_fractional_crystal[0], + vmin = 0, + vmax = 1, + cmap = 'turbo', + ) + +def plot_crystal_amorphous_fraction( + self, + index_ring_plot = 0, + plot_range = (0.0,1.0), + sigma = 0.0, + cmap = 'PiYG', + legend = True, + ticks = False, + figsize = (5,4), + ): + """ + Plotting function for the crystal / amorphous fraction image. + + """ + + + sig_crys = self.signal_crystal[index_ring_plot].copy() + sig_amor = self.signal_amorphous[index_ring_plot].copy() + + if sigma > 0.0: + sig_crys = gaussian_filter( + sig_crys, + sigma, + mode = 'nearest', + ) + sig_amor = gaussian_filter( + sig_amor, + sigma, + mode = 'nearest', + ) + sig_sum = sig_crys + sig_amor + sub = sig_sum > 0.0 + signal_fractional_crystal = sig_crys.copy() + signal_fractional_crystal[sub] /= sig_sum[sub] + + + # im_plot = self.signal_fractional_crystal[index_ring_plot].copy() + # if sigma > 0.0: + # im_plot = gaussian_filter( + # im_plot, + # sigma, + # mode = 'nearest', + # ) + + fig,ax = plt.subplots(figsize=figsize) + h_im = ax.imshow( + signal_fractional_crystal, + vmin = plot_range[0], + vmax = plot_range[1], + cmap = cmap, + ) + if ticks is False: + ax.set_xticks([]) + ax.set_yticks([]) + + if legend is True: + fig.colorbar( + h_im, + ax=ax, + # orientation='horizontal', + # fraction=0.1, + ) + + def refine_peaks( From 77b21dfd80f80cbd60cc229dc31f7e6a7f609444 Mon Sep 17 00:00:00 2001 From: Colin Date: Tue, 16 Jul 2024 20:11:25 -0700 Subject: [PATCH 11/15] updating plotting --- py4DSTEM/process/polar/polar_peaks.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 31238580b..4d1cc67f5 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -951,6 +951,11 @@ def fit_crystal_amorphous_fraction( )) # init + self.signal_background = np.zeros(( + 2, + self._datacube.shape[0], + self._datacube.shape[1], + )) self.signal_amorphous_peaks = np.zeros(( num_rings, self._datacube.shape[0], @@ -1011,6 +1016,7 @@ def fit_crystal_amorphous_fraction( im_polar_sub = im_polar - (basis @ coefs) # Output signals + self.signal_background[:,rx,ry] = coefs[:2] self.signal_amorphous_peaks[:,rx,ry] = coefs[2:] for a0 in range(num_rings): self.signal_crystal[a0,rx,ry] = np.sum( @@ -1098,12 +1104,13 @@ def plot_crystal_amorphous_fraction( ax.set_yticks([]) if legend is True: - fig.colorbar( + cbar = fig.colorbar( h_im, ax=ax, # orientation='horizontal', # fraction=0.1, ) + cbar.ax.set_ylabel('More Amorphous <----> More Crystalline') From e434dfe725d8440bf55c3f0ab7c4cba773d89c41 Mon Sep 17 00:00:00 2001 From: Colin Date: Fri, 23 Aug 2024 16:19:55 -0700 Subject: [PATCH 12/15] black formatting --- py4DSTEM/process/diffraction/flowlines.py | 22 +-- py4DSTEM/process/polar/polar_analysis.py | 67 ++++--- py4DSTEM/process/polar/polar_peaks.py | 221 ++++++++++++---------- 3 files changed, 164 insertions(+), 146 deletions(-) diff --git a/py4DSTEM/process/diffraction/flowlines.py b/py4DSTEM/process/diffraction/flowlines.py index 177baf8b2..07f8d93d7 100644 --- a/py4DSTEM/process/diffraction/flowlines.py +++ b/py4DSTEM/process/diffraction/flowlines.py @@ -812,28 +812,28 @@ def make_flowline_rainbow_legend( Parameters ---------- - im_size (np.array): + im_size (np.array): Size of legend image in pixels. - sym_rotation_order (int): + sym_rotation_order (int): rotational symmety for colouring - theta_offset_degrees (float): + theta_offset_degrees (float): Offset the anglular coloring by this value in degrees. Rotation is Q with respect to R, in the positive (counter clockwise) direction. - white_background (bool): + white_background (bool): For either color or greyscale output, switch to white background (from black). - return_image (bool): + return_image (bool): Return the image array. - radial_range (np.array): + radial_range (np.array): Inner and outer radius for the legend ring. - plot_legend (bool): + plot_legend (bool): Plot the generated legend. - figsize (tuple or list): + figsize (tuple or list): Size of the plotted legend. Returns ---------- - im_legend (array): + im_legend (array): Image array for the legend. """ @@ -854,8 +854,8 @@ def make_flowline_rainbow_legend( amp = np.abs(z) vmin = np.min(amp) vmax = np.max(amp) - ph = np.angle(z)# + hue_offset - h = np.mod(ph/(2*np.pi), 1) + ph = np.angle(z) # + hue_offset + h = np.mod(ph / (2 * np.pi), 1) s = 0.85 * np.ones_like(h) v = (amp - vmin) / (vmax - vmin) im_legend = hsv_to_rgb(np.dstack((h, s, v))) diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index 19e630d8d..1366b5ecf 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1028,7 +1028,7 @@ def cluster_grains( threshold_add=0.05, threshold_grow=0.2, angle_tolerance_deg=5, - distance_tolerance_px = 1, + distance_tolerance_px=1, plot_grain_clusters=True, returncalc=False, **kwargs @@ -1061,11 +1061,11 @@ def cluster_grains( Returns ---------- cluster_sizes: np.array - (N,) shape vector giving the size of N grains as the # of probe positions + (N,) shape vector giving the size of N grains as the # of probe positions cluster_sig: np.array (N,) shape vector giving the mean signal of each of N grains. cluster_inds: list of np.array - Length N list of probe positions for each grain, stored as arrays with + Length N list of probe positions for each grain, stored as arrays with shape (2,M) for M probe positions where rows correspond to the (x,y) indices. """ @@ -1110,7 +1110,9 @@ def cluster_grains( tol = np.deg2rad(angle_tolerance_deg) # Main loop - vec_search = np.arange(-distance_tolerance_px, distance_tolerance_px+1, dtype="int") + vec_search = np.arange( + -distance_tolerance_px, distance_tolerance_px + 1, dtype="int" + ) search = True while search is True: inds_grain = np.argmax(sig) @@ -1243,7 +1245,6 @@ def plot_clusters( """ - # init im_plot = np.zeros( ( @@ -1260,7 +1261,10 @@ def plot_clusters( ) if weight_by_area: - weights = self.cluster_sizes / np.max(self.cluster_sizes) * (1-fill_grains) + fill_grains + weights = ( + self.cluster_sizes / np.max(self.cluster_sizes) * (1 - fill_grains) + + fill_grains + ) else: weights = np.ones_like(self.cluster_sizes) @@ -1312,9 +1316,12 @@ def plot_clusters( self.cluster_inds[a0][0, :], self.cluster_inds[a0][1, :], ] = True - im_plot += gaussian_filter( - im_grain.astype("float"), sigma=smooth_grains, mode="nearest" - ) * weights[a0] + im_plot += ( + gaussian_filter( + im_grain.astype("float"), sigma=smooth_grains, mode="nearest" + ) + * weights[a0] + ) # im_plot[ # self.cluster_inds[a0][0,:], @@ -1417,11 +1424,11 @@ def plot_clusters_area( width=0.8 * pixel_area * area_step, ) else: - ax.bar( + ax.bar( area * pixel_area, hist, width=0.8 * pixel_area * area_step, - ) + ) ax.set_xlim((0, area_max * pixel_area)) ax.set_xlabel("Grain Area (" + pixel_area_units + ")") if weight_intensity: @@ -1485,7 +1492,7 @@ def plot_clusters_diameter( """ if cluster_sizes is None: - cluster_sizes = self.cluster_sizes.copy().astype('float') + cluster_sizes = self.cluster_sizes.copy().astype("float") cluster_diam = 0.5 * np.sqrt(cluster_sizes / np.pi) # units @@ -1506,25 +1513,25 @@ def plot_clusters_diameter( # histogram if weight_intensity: hist = np.bincount( - np.round(cluster_diam[sub] / diameter_step).astype('int'), + np.round(cluster_diam[sub] / diameter_step).astype("int"), weights=self.cluster_sig[sub], minlength=diameter.shape[0], ) elif weight_area: hist = np.bincount( - np.round(cluster_diam[sub] / diameter_step).astype('int'), + np.round(cluster_diam[sub] / diameter_step).astype("int"), weights=cluster_sizes[sub], minlength=diameter.shape[0], ) elif weight_diameter: hist = np.bincount( - np.round(cluster_diam[sub] / diameter_step).astype('int'), + np.round(cluster_diam[sub] / diameter_step).astype("int"), weights=cluster_diam[sub], minlength=diameter.shape[0], ) else: hist = np.bincount( - np.round(cluster_diam[sub] / diameter_step).astype('int'), + np.round(cluster_diam[sub] / diameter_step).astype("int"), minlength=diameter.shape[0], ) @@ -1537,18 +1544,18 @@ def plot_clusters_diameter( width=0.8 * diameter_step, ) else: - ax.bar( + ax.bar( diameter, hist, width=0.8 * diameter_step, - ) + ) ax.set_xlim((0, diameter_max)) ax.set_xlabel("Domain Size (" + pixel_area_units + ")") if weight_intensity: ax.set_ylabel("Intensity-Weighted Domain Size") elif weight_area: - ax.set_ylabel("Total Domain Area") + ax.set_ylabel("Total Domain Area") elif weight_diameter: ax.set_ylabel("Total Domain Diameter") else: @@ -1561,7 +1568,6 @@ def plot_clusters_diameter( return fig, ax - def plot_clusters_max_length( self, cluster_inds=None, @@ -1612,13 +1618,13 @@ def plot_clusters_max_length( """ if cluster_inds is None: - cluster_inds = self.cluster_inds.copy().astype('float') + cluster_inds = self.cluster_inds.copy().astype("float") # init num_clusters = len(cluster_inds) cluster_sizes = np.zeros(num_clusters) cluster_length = np.zeros(num_clusters) - t_all = np.linspace(0,np.pi,45,endpoint=False) + t_all = np.linspace(0, np.pi, 45, endpoint=False) ct = np.cos(t_all) st = np.sin(t_all) @@ -1631,14 +1637,13 @@ def plot_clusters_max_length( length_current_max = 0 for a1 in range(t_all.size): - p = x0*ct[a1] + y0*st[a1] + p = x0 * ct[a1] + y0 * st[a1] length_current_max = np.maximum( length_current_max, np.max(p) - np.min(p) + 1, ) cluster_length[a0] = length_current_max - # units cluster_length *= pixel_size @@ -1657,25 +1662,25 @@ def plot_clusters_max_length( # histogram if weight_intensity: hist = np.bincount( - np.round(cluster_length[sub] / length_step).astype('int'), + np.round(cluster_length[sub] / length_step).astype("int"), weights=self.cluster_sig[sub], minlength=length.shape[0], ) elif weight_area: hist = np.bincount( - np.round(cluster_length[sub] / length_step).astype('int'), + np.round(cluster_length[sub] / length_step).astype("int"), weights=cluster_sizes[sub], minlength=length.shape[0], ) elif weight_diameter: hist = np.bincount( - np.round(cluster_length[sub] / length_step).astype('int'), + np.round(cluster_length[sub] / length_step).astype("int"), weights=cluster_diam[sub], minlength=length.shape[0], ) else: hist = np.bincount( - np.round(cluster_length[sub] / length_step).astype('int'), + np.round(cluster_length[sub] / length_step).astype("int"), minlength=length.shape[0], ) @@ -1688,18 +1693,18 @@ def plot_clusters_max_length( width=0.8 * length_step, ) else: - ax.bar( + ax.bar( length, hist, width=0.8 * length_step, - ) + ) ax.set_xlim((0, length_max)) ax.set_xlabel("Maximum Domain Length (" + pixel_area_units + ")") if weight_intensity: ax.set_ylabel("Intensity-Weighted Domain Lengths") elif weight_area: - ax.set_ylabel("Total Domain Area") + ax.set_ylabel("Total Domain Area") elif weight_diameter: ax.set_ylabel("Total Domain Diameter") else: diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index 4d1cc67f5..cb89c58ad 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -35,14 +35,14 @@ def find_peaks_single_pattern( remove_masked_peaks=False, scale_sigma_annular=0.5, scale_sigma_radial=0.25, - refine_subpixel = True, + refine_subpixel=True, return_background=False, plot_result=True, plot_smoothed_image=False, plot_power_scale=1.0, plot_scale_size=10.0, figsize=(12, 6), - figax = None, + figax=None, returnfig=False, **kwargs ): @@ -258,18 +258,23 @@ def find_peaks_single_pattern( ) # If needed, refine peaks using parabolic fit - peaks = peaks.astype('float') + peaks = peaks.astype("float") if refine_subpixel: for a0 in range(peaks.shape[0]): - if peaks[a0,1] > 0 and peaks[a0,1] < im_polar.shape[1]-1: + if peaks[a0, 1] > 0 and peaks[a0, 1] < im_polar.shape[1] - 1: im_crop = im_polar_sm[ - (peaks[a0,0] + np.arange(-1,2)).astype('int')[:,None] % im_polar.shape[0], - (peaks[a0,1] + np.arange(-1,2)).astype('int'), + (peaks[a0, 0] + np.arange(-1, 2)).astype("int")[:, None] + % im_polar.shape[0], + (peaks[a0, 1] + np.arange(-1, 2)).astype("int"), ] - dx = (im_crop[2,1] - im_crop[0,1]) / (4*im_crop[1,1] - 2*im_crop[2,1] - 2*im_crop[0,1]) - dy = (im_crop[1,2] - im_crop[1,0]) / (4*im_crop[1,1] - 2*im_crop[1,2] - 2*im_crop[1,0]) - peaks[a0,0] += dx - peaks[a0,1] += dy + dx = (im_crop[2, 1] - im_crop[0, 1]) / ( + 4 * im_crop[1, 1] - 2 * im_crop[2, 1] - 2 * im_crop[0, 1] + ) + dy = (im_crop[1, 2] - im_crop[1, 0]) / ( + 4 * im_crop[1, 1] - 2 * im_crop[1, 2] - 2 * im_crop[1, 0] + ) + peaks[a0, 0] += dx + peaks[a0, 1] += dy # combine peaks into one array peaks_all = np.column_stack((peaks, peaks_int, peaks_prom)) @@ -347,14 +352,14 @@ def find_peaks_single_pattern( vmax = kwargs.pop("vmax", 1) vmin = kwargs.pop("vmin", 0) show( - im_plot, - figax=(fig, ax), - cmap=cmap, - vmax=vmax, - vmin=vmin, + im_plot, + figax=(fig, ax), + cmap=cmap, + vmax=vmax, + vmin=vmin, **kwargs, ) - ax.set_aspect('auto') + ax.set_aspect("auto") # peaks ax.scatter( @@ -765,7 +770,7 @@ def plot_radial_peaks( ) ax.yaxis.get_ticklocs(minor=True) ax.minorticks_on() - ax.yaxis.set_tick_params(which='minor',left=False) + ax.yaxis.set_tick_params(which="minor", left=False) ax.xaxis.grid(True) ax.set_xlim((q_bins[0], q_bins[-1])) if q_pixel_units: @@ -802,7 +807,7 @@ def model_radial_background( refine_model=True, plot_result=True, figsize=(8, 4), - returnfig = True, + returnfig=True, ): """ User provided radial background model, of the form: @@ -914,17 +919,18 @@ def background_model(q, *coefs): ) if returnfig: - return fig,ax + return fig, ax + def fit_crystal_amorphous_fraction( self, - fitting_range_sigma = 2.0, - plot_result = True, - figsize = (4,4), - progress_bar = True, - ): + fitting_range_sigma=2.0, + plot_result=True, + figsize=(4, 4), + progress_bar=True, +): """ - Fit an amorphous halo and the crystalline peak signals from each + Fit an amorphous halo and the crystalline peak signals from each polar transformed image, in order to estimate the fraction of crystalline and amorphous signal. @@ -938,57 +944,67 @@ def fit_crystal_amorphous_fraction( Returns ---------- - - """ + + """ # Basis for amorphous fitting function num_rings = np.round((self.background_coefs.shape[0] - 3) / 3).astype("int") num_bins = self.polar_shape[1] - basis = np.ones(( - num_bins, - num_rings + 2, - )) + basis = np.ones( + ( + num_bins, + num_rings + 2, + ) + ) # init - self.signal_background = np.zeros(( - 2, - self._datacube.shape[0], - self._datacube.shape[1], - )) - self.signal_amorphous_peaks = np.zeros(( - num_rings, - self._datacube.shape[0], - self._datacube.shape[1], - )) - self.signal_crystal = np.zeros(( - num_rings, - self._datacube.shape[0], - self._datacube.shape[1], - )) - crystal_fitting_mask = np.zeros(( - num_rings, - num_bins, - ), dtype = 'bool') + self.signal_background = np.zeros( + ( + 2, + self._datacube.shape[0], + self._datacube.shape[1], + ) + ) + self.signal_amorphous_peaks = np.zeros( + ( + num_rings, + self._datacube.shape[0], + self._datacube.shape[1], + ) + ) + self.signal_crystal = np.zeros( + ( + num_rings, + self._datacube.shape[0], + self._datacube.shape[1], + ) + ) + crystal_fitting_mask = np.zeros( + ( + num_rings, + num_bins, + ), + dtype="bool", + ) # background model sigma_0 = self.background_coefs[2] - basis[:,1] = np.exp(self.qq**2/(-2.0*sigma_0**2)) + basis[:, 1] = np.exp(self.qq**2 / (-2.0 * sigma_0**2)) # amorphous halos / rings for a0 in range(num_rings): ring_sigma = self.background_coefs[3 * a0 + 4] ring_position = self.background_coefs[3 * a0 + 5] - basis[:,2+a0] = np.exp( - (self.qq - ring_position)**2 \ - /(-2.0*ring_sigma**2) + basis[:, 2 + a0] = np.exp( + (self.qq - ring_position) ** 2 / (-2.0 * ring_sigma**2) ) sub = np.logical_and( - self.qq > ring_position - ring_sigma*fitting_range_sigma, - self.qq <= ring_position + ring_sigma*fitting_range_sigma, + self.qq > ring_position - ring_sigma * fitting_range_sigma, + self.qq <= ring_position + ring_sigma * fitting_range_sigma, ) - crystal_fitting_mask[a0,sub] = True + crystal_fitting_mask[a0, sub] = True # main loop over probe positions for rx, ry in tqdmnd( @@ -999,10 +1015,10 @@ def fit_crystal_amorphous_fraction( desc="Refining peaks ", unit=" probe positions", disable=not progress_bar, - ): + ): # polar signal - im_polar = self.data[rx,ry] - im_polar_med = np.ma.median(im_polar, axis = 0) + im_polar = self.data[rx, ry] + im_polar_med = np.ma.median(im_polar, axis=0) # fit amorphous background coefficients coefs = np.linalg.lstsq( @@ -1016,15 +1032,15 @@ def fit_crystal_amorphous_fraction( im_polar_sub = im_polar - (basis @ coefs) # Output signals - self.signal_background[:,rx,ry] = coefs[:2] - self.signal_amorphous_peaks[:,rx,ry] = coefs[2:] + self.signal_background[:, rx, ry] = coefs[:2] + self.signal_amorphous_peaks[:, rx, ry] = coefs[2:] for a0 in range(num_rings): - self.signal_crystal[a0,rx,ry] = np.sum( - im_polar_sub[:,crystal_fitting_mask[a0]], + self.signal_crystal[a0, rx, ry] = np.sum( + im_polar_sub[:, crystal_fitting_mask[a0]], ) - self.signal_amorphous_peaks = np.maximum(self.signal_amorphous_peaks,0.0) - self.signal_crystal = np.maximum(self.signal_crystal,0.0) + self.signal_amorphous_peaks = np.maximum(self.signal_amorphous_peaks, 0.0) + self.signal_crystal = np.maximum(self.signal_crystal, 0.0) # convert amorphous signal from peaks to integrated intensity self.signal_amorphous = self.signal_amorphous_peaks.copy() @@ -1040,30 +1056,30 @@ def fit_crystal_amorphous_fraction( # plotting if plot_result: - fig,ax = plt.subplots(figsize=figsize) + fig, ax = plt.subplots(figsize=figsize) ax.imshow( self.signal_fractional_crystal[0], - vmin = 0, - vmax = 1, - cmap = 'turbo', + vmin=0, + vmax=1, + cmap="turbo", ) + def plot_crystal_amorphous_fraction( self, - index_ring_plot = 0, - plot_range = (0.0,1.0), - sigma = 0.0, - cmap = 'PiYG', - legend = True, - ticks = False, - figsize = (5,4), - ): + index_ring_plot=0, + plot_range=(0.0, 1.0), + sigma=0.0, + cmap="PiYG", + legend=True, + ticks=False, + figsize=(5, 4), +): """ Plotting function for the crystal / amorphous fraction image. """ - sig_crys = self.signal_crystal[index_ring_plot].copy() sig_amor = self.signal_amorphous[index_ring_plot].copy() @@ -1071,19 +1087,18 @@ def plot_crystal_amorphous_fraction( sig_crys = gaussian_filter( sig_crys, sigma, - mode = 'nearest', + mode="nearest", ) sig_amor = gaussian_filter( sig_amor, sigma, - mode = 'nearest', + mode="nearest", ) sig_sum = sig_crys + sig_amor sub = sig_sum > 0.0 signal_fractional_crystal = sig_crys.copy() signal_fractional_crystal[sub] /= sig_sum[sub] - # im_plot = self.signal_fractional_crystal[index_ring_plot].copy() # if sigma > 0.0: # im_plot = gaussian_filter( @@ -1092,12 +1107,12 @@ def plot_crystal_amorphous_fraction( # mode = 'nearest', # ) - fig,ax = plt.subplots(figsize=figsize) + fig, ax = plt.subplots(figsize=figsize) h_im = ax.imshow( signal_fractional_crystal, - vmin = plot_range[0], - vmax = plot_range[1], - cmap = cmap, + vmin=plot_range[0], + vmax=plot_range[1], + cmap=cmap, ) if ticks is False: ax.set_xticks([]) @@ -1105,14 +1120,12 @@ def plot_crystal_amorphous_fraction( if legend is True: cbar = fig.colorbar( - h_im, - ax=ax, - # orientation='horizontal', + h_im, + ax=ax, + # orientation='horizontal', # fraction=0.1, ) - cbar.ax.set_ylabel('More Amorphous <----> More Crystalline') - - + cbar.ax.set_ylabel("More Amorphous <----> More Crystalline") def refine_peaks( @@ -1462,39 +1475,39 @@ def make_orientation_histogram( Parameters ---------- - radial_ranges (np array) + radial_ranges (np array) Size (N x 2) array for N radial bins, or (2,) for a single bin. - orientation_flip_sign (bool) + orientation_flip_sign (bool) Flip the direction of theta orientation_offset_degrees (float) Offset for orientation angles. This value is a rotation of Q space with respect to R space. orientation_separate_bins (bool) whether to place multiple angles into multiple radial bins. - upsample_factor (float) + upsample_factor (float) Upsample factor - use_refined_peaks (float) + use_refined_peaks (float) Use refined peak positions - use_peak_sigma (float) + use_peak_sigma (float) Spread signal along annular direction using measured std. - theta_step_deg (float) + theta_step_deg (float) Step size along annular direction in degrees - sigma_x (float) + sigma_x (float) Smoothing in x direction before upsample - sigma_y (float) + sigma_y (float) Smoothing in x direction before upsample - sigma_theta (float) + sigma_theta (float) Smoothing in annular direction (units of bins, periodic) - normalize_intensity_image (bool) + normalize_intensity_image (bool) Normalize to max peak intensity = 1, per image normalize_intensity_stack (bool) Normalize to max peak intensity = 1, all images - progress_bar (bool) + progress_bar (bool) Enable progress bar Returns ---------- - orient_hist (array) + orient_hist (array) 4D array containing Bragg peak intensity histogram [radial_bin x_probe y_probe theta] """ From 3e5844369e1086398caabb54523e0744441387af Mon Sep 17 00:00:00 2001 From: Colin Date: Fri, 23 Aug 2024 19:18:27 -0700 Subject: [PATCH 13/15] Allowing user to specify # of fitting iterations --- py4DSTEM/process/polar/polar_peaks.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index fdd0a6a2e..d2c45c530 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -832,6 +832,7 @@ def model_radial_background( ring_sigma=None, ring_int=None, refine_model=True, + maxfev = None, plot_result=True, figsize=(8, 4), returnfig=True, @@ -861,6 +862,7 @@ def model_radial_background( self.background_mask ] self.background_radial_mean[np.logical_not(self.background_mask)] = 0 + self.background_radial_mean = np.maximum(self.background_radial_mean, 0.0) # init if ring_position is not None: @@ -932,8 +934,9 @@ def background_model(q, *coefs): self.qq[self.background_mask], self.background_radial_mean[self.background_mask], p0=self.background_coefs, - xtol=1e-12, + xtol=1e-8, bounds=(lb, ub), + maxfev = maxfev, )[0] # plotting From 0ccf4dbf732a34dfc4f7f8d1baed332ca9da0a4d Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 24 Aug 2024 15:32:05 -0700 Subject: [PATCH 14/15] ability to return fig handles --- py4DSTEM/process/polar/polar_peaks.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/py4DSTEM/process/polar/polar_peaks.py b/py4DSTEM/process/polar/polar_peaks.py index d2c45c530..64f730581 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -1104,6 +1104,7 @@ def plot_crystal_amorphous_fraction( legend=True, ticks=False, figsize=(5, 4), + returnfig=False, ): """ Plotting function for the crystal / amorphous fraction image. @@ -1157,6 +1158,8 @@ def plot_crystal_amorphous_fraction( ) cbar.ax.set_ylabel("More Amorphous <----> More Crystalline") + if returnfig: + return fig,ax def refine_peaks( self, From 7060ebe64afdf67189713ab4eb0ead1616bb9fab Mon Sep 17 00:00:00 2001 From: cophus Date: Sun, 17 Nov 2024 14:23:19 -0800 Subject: [PATCH 15/15] small fixes --- py4DSTEM/process/polar/polar_datacube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py4DSTEM/process/polar/polar_datacube.py b/py4DSTEM/process/polar/polar_datacube.py index 05c80d2c7..f6d90df5e 100644 --- a/py4DSTEM/process/polar/polar_datacube.py +++ b/py4DSTEM/process/polar/polar_datacube.py @@ -55,7 +55,7 @@ def __init__( """ # attach datacube - assert isinstance(datacube, DataCube) + # assert isinstance(datacube, DataCube) self._datacube = datacube self._datacube.polar = self