diff --git a/py4DSTEM/process/diffraction/flowlines.py b/py4DSTEM/process/diffraction/flowlines.py index 5af64f73e..07f8d93d7 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_analysis.py b/py4DSTEM/process/polar/polar_analysis.py index dbf16ff97..0104e7ba8 100644 --- a/py4DSTEM/process/polar/polar_analysis.py +++ b/py4DSTEM/process/polar/polar_analysis.py @@ -1,12 +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 -from scipy.ndimage import gaussian_filter +import matplotlib.pyplot as plt +import numpy as np from emdfile import tqdmnd +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 def calculate_radial_statistics( @@ -1017,3 +1021,723 @@ 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, + distance_tolerance_px=1, + plot_grain_clusters=True, + returncalc=False, + **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 + 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 + 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: + 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 + vec_search = np.arange( + -distance_tolerance_px, distance_tolerance_px + 1, dtype="int" + ) + 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 + 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) + + 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) + + if returncalc: + return self.cluster_sizes, self.cluster_sig, self.cluster_inds + + +def plot_clusters( + self, + area_min=2, + outline_grains=True, + outline_thickness=1, + fill_grains=0.25, + weight_by_area=False, + smooth_grains=1.0, + cmap="viridis", + figsize=(8, 8), + progress_bar=True, + returncalc=True, + returnfig=False, +): + """ + 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 domains with outlines + outline_thickness: int (optional) + Thickenss of the domain outline + fill_grains: float (optional) + 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) + Domain boundaries are smoothed by this value in pixels. + figsize: tuple + Size of the figure panel + returncalc: bool + Return the domain 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", + ) + + 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 domain 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 * weights[a0] + # 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" + ) + * weights[a0] + ) + + # 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: + if returnfig: + return im_plot, fig, ax + else: + return im_plot + else: + if returnfig: + return fig, ax + + +def plot_clusters_area( + self, + area_min=None, + area_max=None, + area_step=1, + weight_intensity=False, + weight_area=True, + 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. + weight_area: 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 + + """ + + 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) + 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 + ")") + 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_clusters_diameter( + self, + cluster_sizes=None, + diameter_min=None, + diameter_max=None, + diameter_step=1, + weight_intensity=False, + weight_diameter=False, + weight_area=False, + pixel_size=1.0, + pixel_area_units="px", + figsize=(8, 6), + returnfig=False, +): + """ + Plot the cluster sizes + + 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) + 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 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 + 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_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) + if diameter_min is None: + sub = cluster_diam < diameter_max + else: + sub = np.logical_and( + 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"), + minlength=diameter.shape[0], + ) + + # plotting + fig, ax = plt.subplots(figsize=figsize) + if weight_diameter: + ax.bar( + diameter, + hist * diameter, + width=0.8 * diameter_step, + ) + else: + 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") + elif weight_diameter: + ax.set_ylabel("Total Domain Diameter") + else: + 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 + + +# 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..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 @@ -105,7 +105,13 @@ def __init__( plot_reduced_pdf, plot_pdf, background_pca, + cluster_grains, + plot_clusters, + plot_clusters_area, + plot_clusters_diameter, + plot_clusters_max_length, ) + from py4DSTEM.process.polar.polar_peaks import ( find_peaks_single_pattern, find_peaks, @@ -114,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 16d0531db..64f730581 100644 --- a/py4DSTEM/process/polar/polar_peaks.py +++ b/py4DSTEM/process/polar/polar_peaks.py @@ -34,11 +34,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 ): @@ -81,16 +84,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. @@ -161,6 +170,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, @@ -247,6 +258,25 @@ 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)) @@ -303,19 +333,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( @@ -741,6 +786,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( @@ -758,6 +807,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 v_lines is not None: y_min, y_max = ax.get_ylim() @@ -778,8 +832,10 @@ def model_radial_background( ring_sigma=None, ring_int=None, refine_model=True, + maxfev = None, plot_result=True, figsize=(8, 4), + returnfig=True, ): """ User provided radial background model, of the form: @@ -806,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: @@ -877,18 +934,232 @@ 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 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 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_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)) + + # 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_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_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), + returnfig=False, +): + """ + 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: + cbar = fig.colorbar( + h_im, + ax=ax, + # orientation='horizontal', + # fraction=0.1, + ) + cbar.ax.set_ylabel("More Amorphous <----> More Crystalline") + + if returnfig: + return fig,ax def refine_peaks( self, @@ -1235,25 +1506,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 @@ -1261,10 +1550,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) @@ -1327,7 +1617,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