Merged in feature/RAM-3294_nps (pull request #345)

Feature/RAM-3294 Fix NPS

Approved-by: Randy Taylor

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
 --------
 

docs/source/index.rst

@@ -46,6 +46,7 @@
    topics/images
    topics/xim
    topics/contrast
+   topics/nps
    topics/mtf
    topics/scale
    topics/profiles

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 <https://www.aapm.org/pubs/protected_files/ICRU/ICRU_Report_87_Radiation_Dose_and_Image-Quality_Assessment_in_Computed_Tomography_AAPM.pdf>`__
+
+
+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

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