From 3b18a5f29f924232511f57f648f74fb00b639fe9 Mon Sep 17 00:00:00 2001 From: Gerd Duscher <50049264+gduscher@users.noreply.github.com> Date: Sat, 24 Feb 2024 11:53:25 -0500 Subject: [PATCH 1/4] Update version.py --- pyTEMlib/version.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyTEMlib/version.py b/pyTEMlib/version.py index e0fe706f..44efe4c8 100644 --- a/pyTEMlib/version.py +++ b/pyTEMlib/version.py @@ -1,6 +1,6 @@ """ version """ -_version = '0.2024.02.0' +_version = '0.2024.02.1' __version__ = _version -_time = '2024-02-08 19:58:26' +_time = '2024-02-18 19:58:26' From 3400d60f7bf468273c416f7fe85d694d0068a688 Mon Sep 17 00:00:00 2001 From: Gerd Duscher <50049264+gduscher@users.noreply.github.com> Date: Mon, 26 Feb 2024 14:48:31 -0500 Subject: [PATCH 2/4] Update version.py --- pyTEMlib/version.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyTEMlib/version.py b/pyTEMlib/version.py index 44efe4c8..ea5db57a 100644 --- a/pyTEMlib/version.py +++ b/pyTEMlib/version.py @@ -1,6 +1,6 @@ """ version """ -_version = '0.2024.02.1' +_version = '0.2024.02.2' __version__ = _version -_time = '2024-02-18 19:58:26' +_time = '2024-02-26 19:58:26' From a8479cf2f961358f4753088f0338814dd504f122 Mon Sep 17 00:00:00 2001 From: ahoust17 <88668350+ahoust17@users.noreply.github.com> Date: Mon, 11 Mar 2024 17:30:39 -0400 Subject: [PATCH 3/4] center_diffractogram function --- pyTEMlib/image_tools.py | 74 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 66 insertions(+), 8 deletions(-) diff --git a/pyTEMlib/image_tools.py b/pyTEMlib/image_tools.py index 3be384e9..fc545979 100644 --- a/pyTEMlib/image_tools.py +++ b/pyTEMlib/image_tools.py @@ -50,6 +50,11 @@ from collections import Counter +# center diff function +from skimage.filters import threshold_otsu, sobel +from scipy.optimize import leastsq +from sklearn.cluster import DBSCAN + _SimpleITK_present = True try: @@ -275,6 +280,61 @@ def diffractogram_spots(dset, spot_threshold, return_center=True, eps=0.1): return spots, center +def center_diffractogram(dset, return_plot = True, histogram_factor = None): + diff = np.array(dset).T.astype(np.float16) + diff[diff < 0] = 0 + + if histogram_factor is not None: + hist, bins = np.histogram(np.ravel(diff), bins=256, range=(0, 1), density=True) + threshold = threshold_otsu(diff, hist = hist * histogram_factor) + else: + threshold = threshold_otsu(diff) + binary = (diff > threshold).astype(float) + + # Find the edges using the Sobel operator + edges = sobel(binary) + edge_points = np.argwhere(edges) + + # Use DBSCAN to cluster the edge points + db = DBSCAN(eps=10, min_samples=100).fit(edge_points) + labels = db.labels_ + if len(set(labels)) == 1: + raise ValueError("DBSCAN clustering resulted in only one group, check the parameters.") + + # Get the largest group of edge points + unique, counts = np.unique(labels, return_counts=True) + counts = dict(zip(unique, counts)) + largest_group = max(counts, key=counts.get) + edge_points = edge_points[labels == largest_group] + + # Fit a circle to the diffraction ring + def calc_distance(c, x, y): + Ri = np.sqrt((x - c[0])**2 + (y - c[1])**2) + return Ri - Ri.mean() + x_m = np.mean(edge_points[:, 1]) + y_m = np.mean(edge_points[:, 0]) + center_guess = x_m, y_m + center, ier = leastsq(calc_distance, center_guess, args=(edge_points[:, 1], edge_points[:, 0])) + mean_radius = np.mean(calc_distance(center, edge_points[:, 1], edge_points[:, 0])) + np.sqrt((edge_points[:, 1] - center[0])**2 + (edge_points[:, 0] - center[1])**2).mean() + + if return_plot: + fig, ax = plt.subplots(1, 3, figsize=(10, 4)) + circle = plt.Circle(center, mean_radius, color='red', fill=False) + ax[0].set_title('Diffractogram') + ax[0].imshow(dset.T, cmap='viridis') + ax[1].set_title('Otsu Binary Image') + ax[1].imshow(binary, cmap='gray') + ax[2].set_title('Edge Detection and Fitting') + ax[2].imshow(edges, cmap='gray') + ax[2].scatter(center[0], center[1], c='r', s=10) + ax[2].add_artist(circle) + for axis in ax: + axis.axis('off') + fig.tight_layout() + + return center + + def adaptive_fourier_filter(dset, spots, low_pass=3, reflection_radius=0.3): """ Use spots in diffractogram for a Fourier Filter @@ -1070,9 +1130,8 @@ def cartesian2polar(x, y, grid, r, t, order=3): return ndimage.map_coordinates(grid, np.array([new_ix, new_iy]), order=order).reshape(new_x.shape) -def warp(diff): - """Takes a centered diffraction pattern (as a sidpy dataset)and warps it to a polar grid""" - """Centered diff can be produced with it.diffractogram_spots(return_center = True)""" +def warp(diff, center): + """Takes a diffraction pattern (as a sidpy dataset)and warps it to a polar grid""" # Define original polar grid nx = np.shape(diff)[0] @@ -1080,15 +1139,14 @@ def warp(diff): # Define center pixel pix2nm = np.gradient(diff.u.values)[0] - center_pixel = [abs(min(diff.u.values)), abs(min(diff.v.values))]//pix2nm - x = np.linspace(1, nx, nx, endpoint=True)-center_pixel[0] - y = np.linspace(1, ny, ny, endpoint=True)-center_pixel[1] + x = np.linspace(1, nx, nx, endpoint=True)-center[0] + y = np.linspace(1, ny, ny, endpoint=True)-center[1] z = diff # Define new polar grid - nr = int(min([center_pixel[0], center_pixel[1], diff.shape[0]-center_pixel[0], diff.shape[1]-center_pixel[1]])-1) - nt = 360*3 + nr = int(min([center[0], center[1], diff.shape[0]-center[0], diff.shape[1]-center[1]])-1) + nt = 360 * 3 r = np.linspace(1, nr, nr) t = np.linspace(0., np.pi, nt, endpoint=False) From e3b01e2d8ea6744d01f22a86b3eba66874e8879e Mon Sep 17 00:00:00 2001 From: ahoust17 <88668350+ahoust17@users.noreply.github.com> Date: Wed, 3 Apr 2024 22:11:53 -0400 Subject: [PATCH 4/4] center_diffractogram plotting improvement --- pyTEMlib/image_tools.py | 106 +++++++++++++++++++++------------------- 1 file changed, 56 insertions(+), 50 deletions(-) diff --git a/pyTEMlib/image_tools.py b/pyTEMlib/image_tools.py index fc545979..8abf9b12 100644 --- a/pyTEMlib/image_tools.py +++ b/pyTEMlib/image_tools.py @@ -280,57 +280,63 @@ def diffractogram_spots(dset, spot_threshold, return_center=True, eps=0.1): return spots, center -def center_diffractogram(dset, return_plot = True, histogram_factor = None): - diff = np.array(dset).T.astype(np.float16) - diff[diff < 0] = 0 +def center_diffractogram(dset, return_plot = True, histogram_factor = None, smoothing = 1, min_samples = 100): + try: + diff = np.array(dset).T.astype(np.float16) + diff[diff < 0] = 0 + + if histogram_factor is not None: + hist, bins = np.histogram(np.ravel(diff), bins=256, range=(0, 1), density=True) + threshold = threshold_otsu(diff, hist = hist * histogram_factor) + else: + threshold = threshold_otsu(diff) + binary = (diff > threshold).astype(float) + smoothed_image = ndimage.gaussian_filter(binary, sigma=smoothing) # Smooth before edge detection + smooth_threshold = threshold_otsu(smoothed_image) + smooth_binary = (smoothed_image > smooth_threshold).astype(float) + # Find the edges using the Sobel operator + edges = sobel(smooth_binary) + edge_points = np.argwhere(edges) + + # Use DBSCAN to cluster the edge points + db = DBSCAN(eps=10, min_samples=min_samples).fit(edge_points) + labels = db.labels_ + if len(set(labels)) == 1: + raise ValueError("DBSCAN clustering resulted in only one group, check the parameters.") + + # Get the largest group of edge points + unique, counts = np.unique(labels, return_counts=True) + counts = dict(zip(unique, counts)) + largest_group = max(counts, key=counts.get) + edge_points = edge_points[labels == largest_group] + + # Fit a circle to the diffraction ring + def calc_distance(c, x, y): + Ri = np.sqrt((x - c[0])**2 + (y - c[1])**2) + return Ri - Ri.mean() + x_m = np.mean(edge_points[:, 1]) + y_m = np.mean(edge_points[:, 0]) + center_guess = x_m, y_m + center, ier = leastsq(calc_distance, center_guess, args=(edge_points[:, 1], edge_points[:, 0])) + mean_radius = np.mean(calc_distance(center, edge_points[:, 1], edge_points[:, 0])) + np.sqrt((edge_points[:, 1] - center[0])**2 + (edge_points[:, 0] - center[1])**2).mean() - if histogram_factor is not None: - hist, bins = np.histogram(np.ravel(diff), bins=256, range=(0, 1), density=True) - threshold = threshold_otsu(diff, hist = hist * histogram_factor) - else: - threshold = threshold_otsu(diff) - binary = (diff > threshold).astype(float) - - # Find the edges using the Sobel operator - edges = sobel(binary) - edge_points = np.argwhere(edges) - - # Use DBSCAN to cluster the edge points - db = DBSCAN(eps=10, min_samples=100).fit(edge_points) - labels = db.labels_ - if len(set(labels)) == 1: - raise ValueError("DBSCAN clustering resulted in only one group, check the parameters.") - - # Get the largest group of edge points - unique, counts = np.unique(labels, return_counts=True) - counts = dict(zip(unique, counts)) - largest_group = max(counts, key=counts.get) - edge_points = edge_points[labels == largest_group] - - # Fit a circle to the diffraction ring - def calc_distance(c, x, y): - Ri = np.sqrt((x - c[0])**2 + (y - c[1])**2) - return Ri - Ri.mean() - x_m = np.mean(edge_points[:, 1]) - y_m = np.mean(edge_points[:, 0]) - center_guess = x_m, y_m - center, ier = leastsq(calc_distance, center_guess, args=(edge_points[:, 1], edge_points[:, 0])) - mean_radius = np.mean(calc_distance(center, edge_points[:, 1], edge_points[:, 0])) + np.sqrt((edge_points[:, 1] - center[0])**2 + (edge_points[:, 0] - center[1])**2).mean() - - if return_plot: - fig, ax = plt.subplots(1, 3, figsize=(10, 4)) - circle = plt.Circle(center, mean_radius, color='red', fill=False) - ax[0].set_title('Diffractogram') - ax[0].imshow(dset.T, cmap='viridis') - ax[1].set_title('Otsu Binary Image') - ax[1].imshow(binary, cmap='gray') - ax[2].set_title('Edge Detection and Fitting') - ax[2].imshow(edges, cmap='gray') - ax[2].scatter(center[0], center[1], c='r', s=10) - ax[2].add_artist(circle) - for axis in ax: - axis.axis('off') - fig.tight_layout() + finally: + if return_plot: + fig, ax = plt.subplots(1, 4, figsize=(10, 4)) + ax[0].set_title('Diffractogram') + ax[0].imshow(dset.T, cmap='viridis') + ax[1].set_title('Otsu Binary Image') + ax[1].imshow(binary, cmap='gray') + ax[2].set_title('Smoothed Binary Image') + ax[2].imshow(smooth_binary, cmap='gray') + ax[3].set_title('Edge Detection and Fitting') + ax[3].imshow(edges, cmap='gray') + ax[3].scatter(center[0], center[1], c='r', s=10) + circle = plt.Circle(center, mean_radius, color='red', fill=False) + ax[3].add_artist(circle) + for axis in ax: + axis.axis('off') + fig.tight_layout() return center