diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 44ddd5e1b..f98a807cd 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -2,6 +2,28 @@ Changelog ========= +v 3.21.0 +-------- + +CT +^^ + +.. warning:: + + In the last release, the noise power spectrum was not being calculated correctly. + We recommend re-running analyses that were using NPS values. + +* The noise power spectrum introduced last version was not working correctly. + The NPS was not subtracting the mean value from the ROI. This has been fixed. + However, as a result of reworking the calculation, the NPS now has its own module: + ``pylinac.core.nps``. This contains several modules for independent calculation of the NPS + and associated metrics like average power, etc. See :ref:`nps`. +* The NPS is now calculated over square ROIs approximately the same size as the circular + uniformity ROIs rather than one central ROI. This is because the resulting spectra + is smoother when averaged using multiple, separate ROIs. +* The ``power_spectrum`` property of the CTP486 module has been renamed to ``power_spectrum_2d`` + and another property, ``power_spectrum_1d`` has been added. + v 3.20.0 -------- diff --git a/docs/source/index.rst b/docs/source/index.rst index 0a2899df6..b36e87a1f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -46,6 +46,7 @@ topics/images topics/xim topics/contrast + topics/nps topics/mtf topics/scale topics/profiles diff --git a/docs/source/topics/nps.rst b/docs/source/topics/nps.rst new file mode 100644 index 000000000..83a7474b7 --- /dev/null +++ b/docs/source/topics/nps.rst @@ -0,0 +1,181 @@ +.. _nps: + +=========== +Noise Power +=========== + +The noise power or noise power spectrum is a metric for evaluating +the noise in an image, and has more characteristics than simply using +the standard deviation. + +The noise power spectrum is calculated by taking the 2D Fourier +transform of :math:`N` images and extracting a radially-average 1D profile. + +The 2D noise power is calculated using ICRU Report 87, equation 11.1 [1]_.: + +.. math:: + + NPS(f_x, f_y) = \frac{\Delta_x \Delta_y}{N_x N_y} \frac{1}{N} \sum_{i=1}^{N} | DFT_{2D} [ I_{i}(x,y) - \overline{I_i} ] | ^2 + +where :math:`\overline{I_i}` is the mean of the image and :math:`I_i(x,y)` is the :math:`i`-th image, :math:`\Delta_x` and :math:`\Delta_y` are the pixel sizes in mm, :math:`N_x` and :math:`N_y` are the number of pixels in the x and y directions, and :math:`N` is the number of images or ROIs evaluated over. +The :math:`DFT_{2D}` is the 2D discrete Fourier transform. + +To calculate the 1D noise power spectrum, the 2D noise power is radially averaged. + +.. math:: + + f_r = \sqrt{f_x^2 + f_y^2} + +Example +------- + +Let's create a synthetic image and calculate the noise power spectrum, showing +how images even with the same standard deviation can have different noise power spectra. + +.. plot:: + + import matplotlib.pyplot as plt + import numpy as np + + from pylinac.core.nps import noise_power_spectrum_2d, noise_power_spectrum_1d, average_power + + np.random.seed(123) + + + def generate_gaussian_noise_map(shape: (int, int), scale: int=10, intensity: float=0.5): + """Generate a Gaussian noise map with varying intensities (clumps). + + Parameters: + - shape: Shape of the noise map (height, width). + - scale: Scale of the Gaussian clumps. + - intensity: Intensity of the noise. + + Returns: + - Gaussian noise map as a NumPy array. + """ + # Create low-resolution noise + low_res_shape = (shape[0] // scale, shape[1] // scale) + low_res_noise = np.random.normal(loc=0, scale=intensity, size=low_res_shape) + + # Upscale the noise to the original resolution + noise_map = np.kron(low_res_noise, np.ones((scale, scale))) + + # Ensure the noise map matches the exact original shape, trimming if necessary + noise_map = noise_map[:shape[0], :shape[1]] + + return noise_map + + + def apply_noise_clumps(image: np.ndarray, noise_map: np.ndarray) -> np.ndarray: + """Apply a Gaussian noise map to an image. + + Parameters: + - image: Original image as a NumPy array. + - noise_map: Gaussian noise map. + + Returns: + - Noisy image as a NumPy array. + """ + noisy_image = image + noise_map + dtype_min, dtype_max = np.iinfo(image.dtype).min, np.iinfo(image.dtype).max + noisy_image = np.clip(noisy_image, dtype_min, dtype_max) + + return noisy_image + + def generate_noisy_image(shape: (int, int), scale: int=10, intensity: float=0.5, dtype: np.dtype = np.uint16) -> np.ndarray: + """Generate a noisy image with Gaussian clumps. + + Parameters: + - image: Original image as a NumPy array. + - scale: Scale of the Gaussian clumps. + - intensity: Intensity of the noise. + + Returns: + - Noisy image as a NumPy array. + """ + latent = np.zeros(shape, dtype=dtype) + noise_map = generate_gaussian_noise_map(latent.shape, scale, intensity) + noisy_image = apply_noise_clumps(latent, noise_map) + return noisy_image + + + img1 = generate_noisy_image(shape=(200, 200), scale=20, intensity=50, dtype=np.uint8) + img2 = generate_noisy_image(shape=(200, 200), scale=2, intensity=50, dtype=np.uint8) + + # 2D power spectrum + # we use pixel size of 1 for simplicity + nps2d_1 = noise_power_spectrum_2d(pixel_size=1, rois=[img1]) + nps2d_2 = noise_power_spectrum_2d(pixel_size=1, rois=[img2]) + + # 1D power spectrum + nps1d_1 = noise_power_spectrum_1d(spectrum_2d=nps2d_1) + nps1d_2 = noise_power_spectrum_1d(spectrum_2d=nps2d_2) + + # plot the two images, their 2D power spectra, and a combined 1d power spectra plot + fig = plt.figure(figsize=(10, 10)) + img1_ax = plt.subplot2grid((3, 2), (0, 0), fig=fig) + img1_ax.imshow(img1, cmap='gray') + img1_ax.set_title('Image 1') + img1_ax.set_xlabel(f'Std: {img1.std():.2f}') + img1_ax.set_ylabel(f'Mean: {img1.mean():.2f}') + + img2_ax = plt.subplot2grid((3, 2), (0, 1)) + img2_ax.imshow(img2, cmap='gray') + img2_ax.set_title('Image 2') + img2_ax.set_xlabel(f'Std: {img2.std():.2f}') + img2_ax.set_ylabel(f'Mean: {img2.mean():.2f}') + + nps2d_1_ax = plt.subplot2grid((3, 2), (1, 0)) + nps2d_1_ax.imshow(nps2d_1, cmap='viridis') + nps2d_1_ax.set_title('NPS 2D') + + nps2d_2_ax = plt.subplot2grid((3, 2), (1, 1)) + nps2d_2_ax.imshow(nps2d_2, cmap='viridis') + nps2d_2_ax.set_title('NPS 2D') + + # turn off ticks + for ax in [img1_ax, img2_ax, nps2d_1_ax, nps2d_2_ax]: + ax.set_xticks([]) + ax.set_yticks([]) + + x_vals = np.arange(len(nps1d_1))/len(nps1d_1) + nps1d_ax = plt.subplot2grid((3, 2), (2, 0), colspan=2) + nps1d_ax.plot(x_vals, nps1d_1, label=f'Image 1; avg power: {average_power(nps1d_1):.2f}', color='b') + nps1d_ax2 = nps1d_ax.twinx() + nps1d_ax2.plot(x_vals, nps1d_2, label=f'Image 2; avg power: {average_power(nps1d_2):.2f}', color='g') + nps1d_ax.set_title('NPS 1D') + nps1d_ax.set_xlabel('Frequency ($mm^{-1}$)') + nps1d_ax.set_ylabel('NPS Image 1') + nps1d_ax.grid(True) + nps1d_ax2.set_ylabel('NPS Image 2') + nps1d_ax.legend(loc='upper right') + nps1d_ax2.legend(loc='lower right') + + plt.tight_layout() + plt.show() + +Even with roughly the same mean and standard deviation, the noise power spectrum is vastly different. + +.. note:: + + The images are the same size in pixels. This image generator is a simplistic approach to generating + synthetic CT images, but is useful for demonstrating the noise power spectrum. + +References +---------- + +.. [1] `International Commission on Radiation Units and Measurements. (2017). ICRU Report 87 `__ + + +API +--- + +.. autofunction:: pylinac.core.nps.noise_power_spectrum_2d + +.. autofunction:: pylinac.core.nps.noise_power_spectrum_1d + +.. autofunction:: pylinac.core.nps.average_power + +.. autofunction:: pylinac.core.nps.max_frequency + +.. autofunction:: pylinac.core.nps.plot_nps1d diff --git a/pylinac/core/contrast.py b/pylinac/core/contrast.py index 2ab0b6382..0fce5396b 100644 --- a/pylinac/core/contrast.py +++ b/pylinac/core/contrast.py @@ -3,7 +3,6 @@ import numpy as np from ..core.utilities import OptionListMixin -from . import validators class Contrast(OptionListMixin): @@ -136,41 +135,3 @@ def weber(feature: float, background: float) -> float: def ratio(feature: float, reference: float) -> float: """The ratio of luminescence""" return feature / reference - - -def power_spectrum_1d(array: np.ndarray) -> np.ndarray: - """Get the 1D noise power spectrum from a 2D numpy array. ChatGPT-made. - - Parameters - ---------- - array: np.ndarray - The 2D numpy array. - - Returns - ------- - array: np.ndarray - The 1D power spectrum as an array. This spectrum is radial averaged. - - References - ---------- - https://bertvandenbroucke.netlify.app/2019/05/24/computing-a-power-spectrum-in-python/ - https://chat.openai.com/share/7d138eb7-bb20-428f-ad61-7c63f2ec435e - """ - validators.double_dimension(array) - # Apply Fourier Transform - f_transform = np.fft.fft2(array) - f_shift = np.fft.fftshift(f_transform) - - # Calculate 2D power spectrum - power_spectrum_2d = np.abs(f_shift) ** 2 - - # Convert to 1D power spectrum by radial averaging - y, x = np.indices(power_spectrum_2d.shape) - center = np.array([(x.max() - x.min()) / 2.0, (y.max() - y.min()) / 2.0]) - r = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2) - r = r.astype(int) - - radial_spectrum = np.bincount(r.ravel(), power_spectrum_2d.ravel()) / np.bincount( - r.ravel() - ) - return radial_spectrum diff --git a/pylinac/core/nps.py b/pylinac/core/nps.py new file mode 100644 index 000000000..0ef69875f --- /dev/null +++ b/pylinac/core/nps.py @@ -0,0 +1,148 @@ +from typing import Iterable + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.axis import Axis + +from . import validators + + +def radial_average(arr: np.ndarray) -> np.ndarray: + """Calculate the radial average of a 2D array about the center pixel""" + # Determine the center of the array + center = np.floor(np.array(arr.shape) / 2) + + # Calculate the Euclidean distance from the center for each element + y, x = np.indices(arr.shape) + r = np.sqrt((x - center[1]) ** 2 + (y - center[0]) ** 2) + + # Bin the elements based on their distance from the center + r = r.astype(int) + + # Create an array to hold the sum of elements in each bin and their counts + tbin = np.bincount(r.ravel(), arr.ravel()) + nr = np.bincount(r.ravel()) + + # Avoid division by zero + nonzero = nr != 0 + radial_mean = np.zeros(nr.shape) + radial_mean[nonzero] = tbin[nonzero] / nr[nonzero] + + return radial_mean + + +def noise_power_spectrum_2d( + pixel_size: float, rois: Iterable[np.ndarray] +) -> np.ndarray: + """Calculate the noise power spectrum in 2D and 1D for a set of square ROIs. + + Based on ICRU 87, equation 11.1 and 11.2. Calculates the 2D FFT of each ROI, then averages the FFTs. + FFTs are shifted so that the zero frequency is in the center and then radially averaged + + Notes + ----- + The ROIs should be 2D and square (i.e. not a rectangle, circle, etc). + The smallest dimension of the ROIs will be used. I.e. the ROIs can + vary in size, but only the smallest dimension will be used. In practice, + ROIs can sometimes vary depending on how they are extracted from an image (at least in pylinac). + Sometimes the ROI may be 32 pixels wide and 30 pixels tall, for instance. In this case, the + 30x30 subset of the ROI will be used. + + + Parameters + ---------- + pixel_size : float + The size of each pixel in mm. + rois : list of ndarray + A list of the ROIs to calculate the NPS over. Each ROI should be a 2D array. + + Returns + ------- + nps2d : ndarray + The 2D noise power spectrum. + + References + ---------- + [ICRU 87] https://www.aapm.org/pubs/protected_files/ICRU/ICRU_Report_87_Radiation_Dose_and_Image-Quality_Assessment_in_Computed_Tomography_AAPM.pdf + """ + length = min(min(roi.shape) for roi in rois) + ffts = np.zeros((length, length, len(rois))) + for idx, roi in enumerate(rois): + rroi = roi[0:length, 0:length] + b = np.abs(np.fft.fft2(rroi - np.mean(rroi))) ** 2 + s = np.fft.fftshift(b) + ffts[:, :, idx] = s + s = np.mean(ffts, axis=-1) + nps2d = pixel_size**2 / length**2 * s + return nps2d + + +def noise_power_spectrum_1d(spectrum_2d: np.ndarray) -> np.ndarray: + """Calculate the 1D noise power spectrum from a 2D spectrum. + + Parameters + ---------- + spectrum_2d : ndarray + The 2D noise power spectrum. + + Returns + ------- + ndarray + The 1D noise power spectrum. + """ + validators.double_dimension(spectrum_2d) + return radial_average(spectrum_2d) + + +def average_power(nps1d: np.ndarray) -> float: + """Average power given a 1D spectra. + + Parameters + ---------- + nps1d : sequence of ndarray + The 1D noise power spectra to average. + + Returns + ------- + float + The average noise power spectrum value + """ + validators.single_dimension(nps1d) + # Calculate the weighted average x position + x_positions = np.linspace(0, 1, len(nps1d)) + return np.average(x_positions, weights=nps1d) + + +def max_frequency(nps1d: np.ndarray) -> float: + """Determine the maximum frequency of the NPS given a 1D spectrum.""" + validators.single_dimension(nps1d) + return np.argmax(nps1d) / len(nps1d) + + +def plot_nps1d(nps1d: np.ndarray, ax: Axis | None = None) -> Axis: + """Plot the 1D noise power spectrum. + + This will plot the power spectrum as a function of frequency. + We also need to remove frequencies above the Nyquist frequency. + + Parameters + ---------- + nps1d : ndarray + The 1D noise power spectrum. + ax : matplotlib.axis.Axis, optional + The axis to plot on. If None, a new figure will be created. + + Returns + ------- + matplotlib.axis.Axis + The axis the plot was created on. + """ + validators.single_dimension(nps1d) + if ax is None: + _, ax = plt.subplots() + ax.plot(np.linspace(0, 1, len(nps1d)), nps1d) + ax.set_title("1D Noise Power Spectrum") + ax.set_xlabel("Frequency ($mm^{-1}$)") + ax.set_ylabel("NPS / ($HU^2 mm^2$)") + ax.grid(True) + return ax diff --git a/pylinac/ct.py b/pylinac/ct.py index 7417946ec..2458189bf 100644 --- a/pylinac/ct.py +++ b/pylinac/ct.py @@ -32,19 +32,23 @@ import matplotlib.pyplot as plt import numpy as np from matplotlib.axes import Axes -from matplotlib.axis import Axis -from matplotlib.figure import Figure from py_linq import Enumerable from scipy import ndimage from skimage import draw, filters, measure, segmentation from skimage.measure._regionprops import RegionProperties from .core import image, pdf -from .core.contrast import Contrast, power_spectrum_1d +from .core.contrast import Contrast from .core.geometry import Line, Point from .core.image import ArrayImage, DicomImageStack, ImageLike, z_position from .core.io import TemporaryZipDirectory, get_url, retrieve_demo_file from .core.mtf import MTF +from .core.nps import ( + average_power, + max_frequency, + noise_power_spectrum_1d, + noise_power_spectrum_2d, +) from .core.profile import CollapsedCircleProfile, FWXMProfile from .core.roi import DiskROI, LowContrastDiskROI, RectangleROI from .core.utilities import ResultBase @@ -928,6 +932,7 @@ class CTP486(CatPhanModule): roi_dist_mm = 53 roi_radius_mm = 10 nominal_value = 0 + nps_rois: dict[str, RectangleROI] roi_settings = { "Top": { "value": nominal_value, @@ -984,33 +989,25 @@ def plot_profiles(self, axis: plt.Axes | None = None) -> None: axis.set_title("Uniformity Profiles") def _setup_rois(self) -> None: + """Generate our NPS ROIs. They are just square versions of the existing ROIs.""" super()._setup_rois() - self.noise_roi = RectangleROI( - array=self.image, - width=self.roi_dist_mm * 1.5 / self.mm_per_pixel, # convert to pixels - height=self.roi_dist_mm * 1.5 / self.mm_per_pixel, # convert to pixels - angle=0, - dist_from_center=0, - phantom_center=self.phan_center, - ) + self.nps_rois = {} + for name, setting in self.roi_settings.items(): + self.nps_rois[name] = RectangleROI( + array=self.image, + width=setting["radius_pixels"] * 2, + height=setting["radius_pixels"] * 2, + angle=setting["angle_corrected"], + dist_from_center=setting["distance_pixels"], + phantom_center=self.phan_center, + ) def plot(self, axis: plt.Axes): - """Plot the ROIs but also the noise power spectrum ROI""" - self.noise_roi.plot2axes(axis, edgecolor="green", linestyle="-.") + """Plot the ROIs but also the noise power spectrum ROIs""" + for nps_roi in self.nps_rois.values(): + nps_roi.plot2axes(axis, edgecolor="green", linestyle="-.") super().plot(axis) - def plot_noise_power_spectrum(self, show: bool = True) -> (Figure, Axis): - """Plot the noise power spectrum of the Uniformity slice.""" - fig, axis = plt.subplots() - axis.semilogy(self.power_spectrum) - axis.set_title("Noise Power Spectrum") - axis.set_xlabel("Spatial Frequency") - axis.set_ylabel("Power Intensity") - axis.grid(which="both") - if show: - plt.show() - return fig, axis - @property def overall_passed(self) -> bool: """Boolean specifying whether all the ROIs passed within tolerance.""" @@ -1035,21 +1032,27 @@ def integral_non_uniformity(self) -> float: return (maxhu - minhu) / (maxhu + minhu + 2000) @cached_property - def power_spectrum(self) -> np.ndarray: + def power_spectrum_2d(self) -> np.ndarray: """The power spectrum of the uniformity ROI.""" - return power_spectrum_1d(self.noise_roi.pixel_array) + return noise_power_spectrum_2d( + pixel_size=self.mm_per_pixel, + rois=[r.pixel_array for r in self.nps_rois.values()], + ) + + @cached_property + def power_spectrum_1d(self) -> np.ndarray: + """The 1D power spectrum of the uniformity ROI.""" + return noise_power_spectrum_1d(self.power_spectrum_2d) @property def avg_noise_power(self) -> float: """The average noise power of the uniformity ROI.""" - spectrum = self.power_spectrum - frequencies = np.arange(len(spectrum)) - return np.sum(frequencies * spectrum) / np.sum(spectrum) + return average_power(self.power_spectrum_1d) @property - def max_noise_power_frequency(self) -> int: + def max_noise_power_frequency(self) -> float: """The frequency of the maximum noise power. 0 means no pattern.""" - return int(np.argmax(self.power_spectrum)) + return max_frequency(self.power_spectrum_1d) class CTP528CP504(CatPhanModule): diff --git a/tests_basic/core/test_contrast.py b/tests_basic/core/test_contrast.py index 0a5b18abc..9b1d1fbbe 100644 --- a/tests_basic/core/test_contrast.py +++ b/tests_basic/core/test_contrast.py @@ -3,7 +3,7 @@ import numpy as np from pylinac.core import contrast -from pylinac.core.contrast import Contrast, power_spectrum_1d +from pylinac.core.contrast import Contrast class TestContrastAlgorithms(TestCase): @@ -99,48 +99,3 @@ def test_contrast_ratio_bad_array(self): arr = np.array((0.5, 1, 1.5)) with self.assertRaises(ValueError): contrast.contrast(arr, Contrast.RATIO) - - -class TestCalculate1DPowerSpectrum(TestCase): - def test_random_noise(self): - # Test with a random noise image - test_image = np.random.rand(256, 256) # noqa: NPY002 - output = power_spectrum_1d(test_image) - self.assertEqual(output.ndim, 1) - self.assertTrue(len(output) > 0) - self.assertEqual(np.argmax(output), 0) - - def test_average_power_uniform_image(self): - # Test with an image with a constant value - test_image = np.ones((256, 256)) - output = power_spectrum_1d(test_image) - self.assertEqual(output.ndim, 1) - self.assertTrue(len(output) > 0) - self.assertEqual(np.argmax(output), 0) - frequencies = np.arange(len(output)) - avg_power = np.sum(frequencies * output) / np.sum(output) - self.assertAlmostEqual(avg_power, 0, places=5) - - def test_average_power_sinusoidal_image(self): - # Test with an image with a sinusoidal pattern - # this oscillates in the y-direction - test_image = np.ones((256, 256)) - test_image = test_image * np.sin(np.arange(256) * 2 * np.pi / 16)[:, np.newaxis] - output = power_spectrum_1d(test_image) - self.assertEqual(output.ndim, 1) - self.assertEqual(np.argmax(output), 15) - frequencies = np.arange(len(output)) - avg_power = np.sum(frequencies * output) / np.sum(output) - self.assertAlmostEqual(avg_power, 15, delta=1) - - def test_edge_cases(self): - # Test with an empty array - test_image = np.array([]) - with self.assertRaises(ValueError): - power_spectrum_1d(test_image) - - # Test with a very small array - test_image = np.random.rand(2, 2) # noqa: NPY002 - output = power_spectrum_1d(test_image) - self.assertEqual(output.ndim, 1) - self.assertTrue(len(output) > 0) diff --git a/tests_basic/core/test_nps.py b/tests_basic/core/test_nps.py new file mode 100644 index 000000000..fd62ca768 --- /dev/null +++ b/tests_basic/core/test_nps.py @@ -0,0 +1,152 @@ +import math +from unittest import TestCase + +import numpy as np + +from pylinac.core.nps import ( + average_power, + max_frequency, + noise_power_spectrum_1d, + noise_power_spectrum_2d, +) + +np.random.seed(123) # noqa: NPY002 + + +def generate_gaussian_noise_map( + shape: (int, int), scale: int = 10, intensity: float = 0.5 +): + """Generate a Gaussian noise map with varying intensities (clumps). + + Parameters: + - shape: Shape of the noise map (height, width). + - scale: Scale of the Gaussian clumps. + - intensity: Intensity of the noise. + + Returns: + - Gaussian noise map as a NumPy array. + """ + # Create low-resolution noise + low_res_shape = (shape[0] // scale, shape[1] // scale) + low_res_noise = np.random.normal( # noqa: NPY002 + loc=0, scale=intensity, size=low_res_shape + ) + + # Upscale the noise to the original resolution + noise_map = np.kron(low_res_noise, np.ones((scale, scale))) + + # Ensure the noise map matches the exact original shape, trimming if necessary + noise_map = noise_map[: shape[0], : shape[1]] + + return noise_map + + +def apply_noise_clumps(image: np.ndarray, noise_map: np.ndarray) -> np.ndarray: + """ + Apply a Gaussian noise map to an image. + + Parameters: + - image: Original image as a NumPy array. + - noise_map: Gaussian noise map. + + Returns: + - Noisy image as a NumPy array. + """ + noisy_image = image + noise_map + dtype_min, dtype_max = np.iinfo(image.dtype).min, np.iinfo(image.dtype).max + noisy_image = np.clip(noisy_image, dtype_min, dtype_max) + + return noisy_image + + +def generate_noisy_image( + shape: (int, int), + scale: int = 10, + intensity: float = 0.5, + dtype: np.dtype = np.uint16, +) -> np.ndarray: + """ + Generate a noisy image with Gaussian clumps. + + Parameters: + - image: Original image as a NumPy array. + - scale: Scale of the Gaussian clumps. + - intensity: Intensity of the noise. + + Returns: + - Noisy image as a NumPy array. + """ + latent = np.zeros(shape, dtype=dtype) + noise_map = generate_gaussian_noise_map(latent.shape, scale, intensity) + noisy_image = apply_noise_clumps(latent, noise_map) + return noisy_image + + +class Test2DSpectrum(TestCase): + def test_single_roi(self): + roi = generate_noisy_image((300, 300), scale=30, intensity=500, dtype=np.uint16) + nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[roi]) + self.assertEqual(nps2d.shape, roi.shape) + + def test_multiple_rois(self): + roi1 = generate_noisy_image( + (300, 300), scale=30, intensity=500, dtype=np.uint16 + ) + roi2 = generate_noisy_image( + (300, 300), scale=10, intensity=100, dtype=np.uint16 + ) + nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[roi1, roi2]) + self.assertEqual(nps2d.shape, roi1.shape) + + def test_take_smallest_shape(self): + roi1 = generate_noisy_image( + (300, 300), scale=30, intensity=500, dtype=np.uint16 + ) + roi2 = generate_noisy_image( + (200, 200), scale=10, intensity=100, dtype=np.uint16 + ) + nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[roi1, roi2]) + self.assertEqual(nps2d.shape, roi2.shape) + + +class Test1DSpectrum(TestCase): + def test_1d_spectrum_uniform(self): + nps2d = np.ones((300, 300)) + nps1d = noise_power_spectrum_1d(nps2d) + self.assertAlmostEqual(nps1d[0], 1, delta=0.0001) + + def test_1d_spectrum(self): + roi = generate_noisy_image((300, 300), scale=30, intensity=500, dtype=np.uint16) + nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[roi]) + nps1d = noise_power_spectrum_1d(nps2d) + # shape is same as diagonal distance from center to corner + self.assertEqual(len(nps1d), math.ceil(300 * math.sqrt(2) / 2)) + + +class TestAvgPower(TestCase): + def setUp(self) -> None: + self.roi = generate_noisy_image( + (300, 300), scale=30, intensity=500, dtype=np.uint16 + ) + self.nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[self.roi]) + + def test_avg_power(self): + nps1d = noise_power_spectrum_1d(self.nps2d) + avg_power = average_power(nps1d) + self.assertAlmostEqual(avg_power, 0.0207, delta=0.005) + + def test_odd_roi_size_same_as_even(self): + nps1d = noise_power_spectrum_1d(self.nps2d) + avg_power_even = average_power(nps1d) + nps2d_odd = noise_power_spectrum_2d(pixel_size=1, rois=[self.roi[:-1, :-1]]) + avg_power_odd = average_power(noise_power_spectrum_1d(nps2d_odd)) + self.assertAlmostEqual(avg_power_even, avg_power_odd, delta=0.0005) + + +class TestFrequency(TestCase): + def test_frequency(self): + roi = generate_noisy_image((300, 300), scale=30, intensity=500, dtype=np.uint16) + nps2d = noise_power_spectrum_2d(pixel_size=1, rois=[roi]) + nps1d = noise_power_spectrum_1d(nps2d) + f = max_frequency(nps1d) + self.assertAlmostEqual(f, 0.0094, delta=0.0001) diff --git a/tests_basic/test_cbct.py b/tests_basic/test_cbct.py index 3c862a8a4..4ff2f454a 100644 --- a/tests_basic/test_cbct.py +++ b/tests_basic/test_cbct.py @@ -584,7 +584,7 @@ class CatPhan604Test(CatPhanMixin, TestCase): unif_values = {"Center": -3, "Left": 0, "Right": 0, "Top": 0, "Bottom": 0} mtf_values = {50: 0.43} lowcon_visible = 1 # changed w/ visibility refactor in v3.0 - avg_noise_power = 17.65 + avg_noise_power = 0.252 class CatPhan504Mixin(CatPhanMixin): @@ -812,7 +812,7 @@ class CBCT2(CatPhan504Mixin, TestCase): lowcon_visible = 2 avg_line_length = 49.9 slice_thickness = 2.4 - avg_noise_power = 0.05 + avg_noise_power = 0.395 class CBCT3(CatPhan504Mixin, TestCase): @@ -991,7 +991,7 @@ class CBCT11(CatPhan504Mixin, TestCase): avg_line_length = 49.94 lowcon_visible = 4 slice_thickness = 2.35 - avg_noise_power = 14.6 + avg_noise_power = 0.226 class CBCT12(CatPhan504Mixin, TestCase): @@ -1498,7 +1498,7 @@ class CatPhan604wJig(CatPhan604Mixin, TestCase): unif_values = {"Center": 6, "Left": 4, "Right": 13, "Top": 10, "Bottom": 7} mtf_values = {50: 0.28} lowcon_visible = 1 - avg_noise_power = 0.96 + avg_noise_power = 0.287 class CatPhan604wJig2(CatPhan604Mixin, TestCase):