From ac4ae83c3edf6aff00f34acc1a3b3e3b94806ff3 Mon Sep 17 00:00:00 2001 From: jennmald Date: Mon, 27 Apr 2026 16:30:26 -0400 Subject: [PATCH 1/5] cosmic analysis script rewritten in numpy --- src/cditools/analysis-scripts/cosmics.py | 91 ++++++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 src/cditools/analysis-scripts/cosmics.py diff --git a/src/cditools/analysis-scripts/cosmics.py b/src/cditools/analysis-scripts/cosmics.py new file mode 100644 index 0000000..a66c983 --- /dev/null +++ b/src/cditools/analysis-scripts/cosmics.py @@ -0,0 +1,91 @@ +""" +Cosmics for a counting detector will trigger values of 1 or 2, +which are hard to see in default color settings, so we mask out +high and zero values so that cosmics stand out more. + +Data is accessed like data[image][row][column], +rather than data[image][x][y] +""" +from typing import Any + +import h5py +from numpy.typing import NDArray +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +from skimage.measure import label +from scipy.stats import poisson + +def setup(det_name:str, date_dir:str, file_name:str, proposal_id:str) -> NDArray[np.floating[Any]]: + """ Setups the data to use later for analysis. + + Args: + det_name (str): name of directory for the detector, e.g. "merlines-1" + date_dir (str): data in format "YYYY/MM/`DD" + file_name (str): file name without the .h5 extension, e.g. "scan_0000" + proposal_id (str): current PASS proposal ID + + Returns: + data: dataset to be analyzed, in the form of a 3D array with dimensions (num_images, rows, cols) + """ + file_path = f'/nsls2/data/cdi/proposals/commissioning/pass-{proposal_id}/assets/{det_name}/{date_dir}/{file_name}.h5' + with h5py.File(file_path, 'r') as f: + return np.asarray(f['entry']['data']['data']) + +def movie(data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5) -> None: + # Set colormapping and mask + cmap = mpl.colormaps['gray'].with_extremes(under='aliceblue', over='r') + norm = mpl.colors.Normalize(vmin, vmax) + for n in range(data.shape[0]): + plt.imshow(data[n], cmap=cmap, norm=norm) + plt.title(f"Image {n}") + plt.draw() + plt.pause(0.1) + +def label_cosmics(data: NDArray[np.floating[Any]]) -> tuple[list[int], float]: + counts = [] + cosmic_sum = 0 + for n in range(data.shape[0]): + image = data[n].copy() + #find indices of top 3 pixels + flat_idx = np.argpartition(image, -3, axis=None)[-3:] + coords = np.unravel_index(flat_idx, image.shape) + # set them to 0 + image[coords] = 0 + _, num = label(image, return_num = True) + counts.append(num) + cosmic_sum+=num + + average = cosmic_sum/data.shape[0] + return counts, average + +def plot_counts(counts: NDArray[np.floating[Any]]): + counts = np.array(counts) + k_vals = np.arange(counts.max() + 1) + fit_vals = poisson.pmf(k_vals, mu=np.mean(counts)) + + plt.hist(counts, bins=np.arange(counts.max()+2)-0.5, density=True, alpha=0.6) + plt.plot(k_vals, fit_vals, 'o-', label="Poisson fit") + plt.xlabel("Counts per image") + plt.ylabel("Probability") + plt.title("Histogram of counts") + plt.legend() + plt.show() + +def check_empty(data: NDArray[np.floating[Any]]): + for n in range(data.shape[0]): + if not np.any(data[n]): + print(f"Image {n} is empty.") + +def find_cosmics(data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5): + cosmic_mask = (data >= vmin) & (data <= vmax) + cosmic_points = set(np.flatnonzero(np.any(cosmic_mask, axis=(1, 2)))) + return cosmic_points + +def show_image(data: NDArray[np.floating[Any]], image_index:int, vmin: int = 1, vmax: int = 5): + cmap = mpl.colormaps['gray'].with_extremes(under='aliceblue', over='r') + norm = mpl.colors.Normalize(vmin, vmax) + plt.imshow(data[image_index], cmap=cmap, norm=norm) + plt.title(f"Image {image_index}") + plt.colorbar() + plt.show() From ba1a04ba9836beb4124cc3beb5bb0b6fab20ac16 Mon Sep 17 00:00:00 2001 From: jennmald Date: Tue, 28 Apr 2026 12:14:43 -0400 Subject: [PATCH 2/5] docs and tests --- src/cditools/analysis-scripts/cosmics.py | 91 -------------- src/cditools/analysis_scripts/cosmics.py | 153 +++++++++++++++++++++++ tests/test_cosmics.py | 26 ++++ 3 files changed, 179 insertions(+), 91 deletions(-) delete mode 100644 src/cditools/analysis-scripts/cosmics.py create mode 100644 src/cditools/analysis_scripts/cosmics.py create mode 100644 tests/test_cosmics.py diff --git a/src/cditools/analysis-scripts/cosmics.py b/src/cditools/analysis-scripts/cosmics.py deleted file mode 100644 index a66c983..0000000 --- a/src/cditools/analysis-scripts/cosmics.py +++ /dev/null @@ -1,91 +0,0 @@ -""" -Cosmics for a counting detector will trigger values of 1 or 2, -which are hard to see in default color settings, so we mask out -high and zero values so that cosmics stand out more. - -Data is accessed like data[image][row][column], -rather than data[image][x][y] -""" -from typing import Any - -import h5py -from numpy.typing import NDArray -import numpy as np -import matplotlib.pyplot as plt -import matplotlib as mpl -from skimage.measure import label -from scipy.stats import poisson - -def setup(det_name:str, date_dir:str, file_name:str, proposal_id:str) -> NDArray[np.floating[Any]]: - """ Setups the data to use later for analysis. - - Args: - det_name (str): name of directory for the detector, e.g. "merlines-1" - date_dir (str): data in format "YYYY/MM/`DD" - file_name (str): file name without the .h5 extension, e.g. "scan_0000" - proposal_id (str): current PASS proposal ID - - Returns: - data: dataset to be analyzed, in the form of a 3D array with dimensions (num_images, rows, cols) - """ - file_path = f'/nsls2/data/cdi/proposals/commissioning/pass-{proposal_id}/assets/{det_name}/{date_dir}/{file_name}.h5' - with h5py.File(file_path, 'r') as f: - return np.asarray(f['entry']['data']['data']) - -def movie(data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5) -> None: - # Set colormapping and mask - cmap = mpl.colormaps['gray'].with_extremes(under='aliceblue', over='r') - norm = mpl.colors.Normalize(vmin, vmax) - for n in range(data.shape[0]): - plt.imshow(data[n], cmap=cmap, norm=norm) - plt.title(f"Image {n}") - plt.draw() - plt.pause(0.1) - -def label_cosmics(data: NDArray[np.floating[Any]]) -> tuple[list[int], float]: - counts = [] - cosmic_sum = 0 - for n in range(data.shape[0]): - image = data[n].copy() - #find indices of top 3 pixels - flat_idx = np.argpartition(image, -3, axis=None)[-3:] - coords = np.unravel_index(flat_idx, image.shape) - # set them to 0 - image[coords] = 0 - _, num = label(image, return_num = True) - counts.append(num) - cosmic_sum+=num - - average = cosmic_sum/data.shape[0] - return counts, average - -def plot_counts(counts: NDArray[np.floating[Any]]): - counts = np.array(counts) - k_vals = np.arange(counts.max() + 1) - fit_vals = poisson.pmf(k_vals, mu=np.mean(counts)) - - plt.hist(counts, bins=np.arange(counts.max()+2)-0.5, density=True, alpha=0.6) - plt.plot(k_vals, fit_vals, 'o-', label="Poisson fit") - plt.xlabel("Counts per image") - plt.ylabel("Probability") - plt.title("Histogram of counts") - plt.legend() - plt.show() - -def check_empty(data: NDArray[np.floating[Any]]): - for n in range(data.shape[0]): - if not np.any(data[n]): - print(f"Image {n} is empty.") - -def find_cosmics(data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5): - cosmic_mask = (data >= vmin) & (data <= vmax) - cosmic_points = set(np.flatnonzero(np.any(cosmic_mask, axis=(1, 2)))) - return cosmic_points - -def show_image(data: NDArray[np.floating[Any]], image_index:int, vmin: int = 1, vmax: int = 5): - cmap = mpl.colormaps['gray'].with_extremes(under='aliceblue', over='r') - norm = mpl.colors.Normalize(vmin, vmax) - plt.imshow(data[image_index], cmap=cmap, norm=norm) - plt.title(f"Image {image_index}") - plt.colorbar() - plt.show() diff --git a/src/cditools/analysis_scripts/cosmics.py b/src/cditools/analysis_scripts/cosmics.py new file mode 100644 index 0000000..84e2418 --- /dev/null +++ b/src/cditools/analysis_scripts/cosmics.py @@ -0,0 +1,153 @@ +""" +Cosmics for a counting detector will trigger values of 1 or 2, +which are hard to see in default color settings, so we mask out +high and zero values so that cosmics stand out more. + +Data is accessed like data[image][row][column], +rather than data[image][x][y] +""" + +from __future__ import annotations + +from typing import Any + +import h5py +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +from numpy.typing import NDArray +from scipy.stats import poisson +from skimage.measure import label + + +def setup( + det_name: str, date_dir: str, file_name: str, proposal_id: str +) -> NDArray[np.floating[Any]]: + """Sets up the data to use later for analysis. + + Args: + det_name (str): name of directory for the detector, e.g. "merlines-1" + date_dir (str): data in format "YYYY/MM/`DD" + file_name (str): file name without the .h5 extension, e.g. "scan_0000" + proposal_id (str): current PASS proposal ID + + Returns: + data: dataset to be analyzed, in the form of a 3D array with dimensions (num_images, rows, cols) + """ + file_path = f"/nsls2/data/cdi/proposals/commissioning/pass-{proposal_id}/assets/{det_name}/{date_dir}/{file_name}.h5" + with h5py.File(file_path, "r") as f: + return np.asarray(f["entry"]["data"]["data"]) + + +def movie(data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5) -> None: + """Creates a movie of the images in the dataset, with a colormap that highlights values between vmin and vmax. + + Args: + data (NDArray[np.floating[Any]]): data set with entries in the form of (num_images, rows, cols) + vmin (int, optional): minimum value for cosmics. Defaults to 1. + vmax (int, optional): maximum value for cosmics. Defaults to 5. + """ + cmap = mpl.colormaps["gray"].with_extremes(under="aliceblue", over="r") + norm = mpl.colors.Normalize(vmin, vmax) + for n in range(data.shape[0]): + plt.imshow(data[n], cmap=cmap, norm=norm) + plt.title(f"Image {n}") + plt.draw() + plt.pause(0.1) + + +def label_cosmics(data: NDArray[np.floating[Any]]) -> tuple[list[int], float]: + """Labels the cosmics in each image and counts the number of cosmics per image, + then calculates the average number of cosmics across all images. + + Args: + data (NDArray[np.floating[Any]]): data set with entries in the form of (num_images, rows, cols) + + Returns: + tuple[list[int], float]: A tuple containing a list of counts of cosmics per image and the average number of cosmics across all images. + """ + counts = [] + cosmic_sum = 0 + for n in range(data.shape[0]): + image = data[n].copy() + # find indices of top 3 pixels + flat_idx = np.argpartition(image, -3, axis=None)[-3:] + coords = np.unravel_index(flat_idx, image.shape) + # set them to 0 + image[coords] = 0 + _, num = label(image, return_num=True) + counts.append(num) + cosmic_sum += num + + average = cosmic_sum / data.shape[0] + return counts, average + + +def plot_counts(counts: NDArray[np.floating[Any]]): + """Plot the histogram of the cosmic counts and fit a poisson distribution + + Args: + counts (NDArray[np.floating[Any]]): array of counts of cosmics per image + """ + counts = np.array(counts) + k_vals = np.arange(counts.max() + 1) + fit_vals = poisson.pmf(k_vals, mu=np.mean(counts)) + + plt.hist(counts, bins=np.arange(counts.max() + 2) - 0.5, density=True, alpha=0.6) + plt.plot(k_vals, fit_vals, "o-", label="Poisson fit") + plt.xlabel("Counts per image") + plt.ylabel("Probability") + plt.title("Histogram of counts") + plt.legend() + plt.show() + + +def check_empty(data: NDArray[np.floating[Any]]) -> list[int]: + """Checks for images that have values of all zero + + Args: + data (NDArray[np.floating[Any]]): data set with entries in the form of (num_images, rows, cols) + Returns: + list[int]: list of indices of images that have all zero values + """ + zero_images = [] + for n in range(data.shape[0]): + if not np.any(data[n]): + zero_images.append(n) + return zero_images + + +def find_cosmics( + data: NDArray[np.floating[Any]], vmin: int = 1, vmax: int = 5 +) -> set[int]: + """Masks images to find the cosmics in + + Args: + data (NDArray[np.floating[Any]]): data set with entries in the form of (num_images, rows, cols) + vmin (int, optional): minimum value for cosmics. Defaults to 1. + vmax (int, optional): maximum value for cosmics. Defaults to 5. + + Returns: + set[int]: set of points that are identified as cosmics, based on the values in the data set being between vmin and vmax + """ + cosmic_mask = (data >= vmin) & (data <= vmax) + return set(np.flatnonzero(np.any(cosmic_mask, axis=(1, 2)))) + + +def show_image( + data: NDArray[np.floating[Any]], image_index: int, vmin: int = 1, vmax: int = 5 +): + """Display the desired image by index + + Args: + data (NDArray[np.floating[Any]]): data set with entries in the form of (num_images, rows, cols) + image_index (int): index of the image to display + vmin (int, optional): minimum value for cosmics. Defaults to 1. + vmax (int, optional): maximum value for cosmics. Defaults to 5. + """ + cmap = mpl.colormaps["gray"].with_extremes(under="aliceblue", over="r") + norm = mpl.colors.Normalize(vmin, vmax) + plt.imshow(data[image_index], cmap=cmap, norm=norm) + plt.title(f"Image {image_index}") + plt.colorbar() + plt.show() diff --git a/tests/test_cosmics.py b/tests/test_cosmics.py new file mode 100644 index 0000000..1f64355 --- /dev/null +++ b/tests/test_cosmics.py @@ -0,0 +1,26 @@ +from __future__ import annotations + +import numpy as np + +from cditools.analysis_scripts.cosmics import check_empty, find_cosmics, label_cosmics + + +def test_check_empty(): + data = np.array([[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [0, 0]]]) + empty_indices = check_empty(data) + assert empty_indices == [0, 2] + + +def test_find_cosmics(): + data = np.array([[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [5, 0]]]) + cosmic_points = find_cosmics(data, vmin=1, vmax=5) + assert cosmic_points == {1, 2} + + +def test_label_cosmics(): + data = np.array([[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [5, 0]]]) + labeled_data = label_cosmics(data) + expected_labeled_data = np.array( + [[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [5, 0]]] + ) + assert np.array_equal(labeled_data, expected_labeled_data) From e9c279e7e3ded812c1618c0d9dab3edd32b95da8 Mon Sep 17 00:00:00 2001 From: jennmald Date: Thu, 30 Apr 2026 11:44:04 -0400 Subject: [PATCH 3/5] clean up copilot suggestions --- src/cditools/analysis_scripts/cosmics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cditools/analysis_scripts/cosmics.py b/src/cditools/analysis_scripts/cosmics.py index 84e2418..7914f3d 100644 --- a/src/cditools/analysis_scripts/cosmics.py +++ b/src/cditools/analysis_scripts/cosmics.py @@ -27,7 +27,7 @@ def setup( Args: det_name (str): name of directory for the detector, e.g. "merlines-1" - date_dir (str): data in format "YYYY/MM/`DD" + date_dir (str): data in format "YYYY/MM/DD" file_name (str): file name without the .h5 extension, e.g. "scan_0000" proposal_id (str): current PASS proposal ID From 76c40e9245edbca5dbeb14615250d3e3305b64a4 Mon Sep 17 00:00:00 2001 From: jennmald Date: Fri, 1 May 2026 12:32:18 -0400 Subject: [PATCH 4/5] copilot suggestions --- src/cditools/analysis_scripts/cosmics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cditools/analysis_scripts/cosmics.py b/src/cditools/analysis_scripts/cosmics.py index 7914f3d..d27fb0b 100644 --- a/src/cditools/analysis_scripts/cosmics.py +++ b/src/cditools/analysis_scripts/cosmics.py @@ -83,7 +83,7 @@ def label_cosmics(data: NDArray[np.floating[Any]]) -> tuple[list[int], float]: return counts, average -def plot_counts(counts: NDArray[np.floating[Any]]): +def plot_counts(counts: NDArray[np.floating[Any]]) -> None: """Plot the histogram of the cosmic counts and fit a poisson distribution Args: From 379eed16eacaaeddc97dba93e3de3bf742d95e33 Mon Sep 17 00:00:00 2001 From: jennmald Date: Fri, 1 May 2026 12:41:27 -0400 Subject: [PATCH 5/5] fix test case for labeling cosmics --- tests/test_cosmics.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_cosmics.py b/tests/test_cosmics.py index 1f64355..0e9d346 100644 --- a/tests/test_cosmics.py +++ b/tests/test_cosmics.py @@ -20,7 +20,5 @@ def test_find_cosmics(): def test_label_cosmics(): data = np.array([[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [5, 0]]]) labeled_data = label_cosmics(data) - expected_labeled_data = np.array( - [[[0, 0], [0, 0]], [[1, 0], [0, 0]], [[0, 0], [5, 0]]] - ) - assert np.array_equal(labeled_data, expected_labeled_data) + expected_labeled_data = ([0, 0, 0], 0.0) + assert labeled_data == expected_labeled_data