Refactor sharpness filters and improve single_sided_correction (#573)

Author: faridyagubbayli

kwave/utils/filters.py

@@ -4,13 +4,13 @@
 import numpy as np
 import scipy
 from scipy.fftpack import fft, fftshift, ifft, ifftshift
-from scipy.signal import convolve, lfilter
+from scipy.signal import lfilter
 
 from kwave.utils.conversion import create_index_at_dim
 
 from .checks import is_number
 from .data import scale_SI
-from .math import find_closest, gaussian, next_pow2, norm_var, sinc
+from .math import find_closest, gaussian, next_pow2, sinc
 from .matrix import num_dim, num_dim2
 from .signals import get_win
 
@@ -27,30 +27,25 @@ def single_sided_correction(func_fft: np.ndarray, fft_len: int, dim: int) -> np.
         dim: The dimension along which to apply the correction.
 
     Returns:
-        The corrected FFT of the function.
+        None, modifies the input array in place to have the corrected FFT of the function.
     """
-    # Determine the slice to use based on FFT length
-    if fft_len % 2:
-        # odd FFT length - multiply all elements except the first one by 2
-        dim_slice = slice(1, None)
-    else:
-        # even FFT length - multiply all elements except the first and last ones by 2
-        dim_slice = slice(1, -1)
-
-    # Create a slice tuple with the appropriate slice at the specified dimension
-    idx_all = [slice(None)] * func_fft.ndim
-    idx_all[dim] = dim_slice
-    idx_tuple = tuple(idx_all)
+    # Create a slice object for each dimension
+    slices = [slice(None)] * func_fft.ndim
 
-    # Apply the correction
-    func_fft[idx_tuple] = func_fft[idx_tuple] * 2
+    if fft_len % 2:  # odd FFT length
+        # Set slice for the specified dimension to select all elements except the first
+        slices[dim] = slice(1, None)
+    else:  # even FFT length
+        # Set slice for the specified dimension to select all elements except first and last
+        slices[dim] = slice(1, -1)
 
-    return func_fft
+    # Apply the slicing and multiply by 2
+    func_fft[tuple(slices)] *= 2
 
 
 def spect(
     func: np.ndarray,
-    Fs: float,
+    fs: float,
     dim: Optional[Union[int, str]] = "auto",
     fft_len: Optional[int] = 0,
     power_two: Optional[bool] = False,
@@ -62,7 +57,7 @@ def spect(
 
     Args:
         func: The signal to analyse.
-        Fs: The sampling frequency in Hz.
+        fs: The sampling frequency in Hz.
         dim: The dimension over which the spectrum is calculated. Defaults to 'auto'.
         fft_len: The length of the FFT. If the set length is smaller than the signal length, the default value is used
                  instead (default = signal length).
@@ -133,10 +128,10 @@ def spect(
     slicing[dim] = slice(0, num_unique_pts)
     func_fft = func_fft[tuple(slicing)]
 
-    func_fft = single_sided_correction(func_fft, fft_len, dim)
+    single_sided_correction(func_fft, fft_len, dim)
 
     # create the frequency axis variable
-    f = np.arange(0, num_unique_pts) * Fs / fft_len
+    f = np.arange(0, num_unique_pts) * fs / fft_len
 
     # calculate the amplitude spectrum
     func_as = np.abs(func_fft)
@@ -152,7 +147,7 @@ def spect(
 
 
 def extract_amp_phase(
-    data: np.ndarray, Fs: float, source_freq: float, dim: Tuple[str, int] = "auto", fft_padding: int = 3, window: str = "Hanning"
+    data: np.ndarray, fs: float, source_freq: float, dim: Tuple[str, int] = "auto", fft_padding: int = 3, window: str = "Hanning"
 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
     """
     Extract the amplitude and phase information at a specified frequency from a vector or matrix of time series data.
@@ -163,7 +158,7 @@ def extract_amp_phase(
 
     Args:
         data: Matrix of time signals [s]
-        Fs: Sampling frequency [Hz]
+        fs: Sampling frequency [Hz]
         source_freq: Frequency at which the amplitude and phase should be extracted [Hz]
         dim: The time dimension of the input data. If 'auto', the highest non-singleton dimension is used.
         fft_padding: The amount of zero padding to apply to the FFT.
@@ -191,7 +186,8 @@ def extract_amp_phase(
     data = win * data
 
     # compute amplitude and phase spectra
-    f, func_as, func_ps = spect(data, Fs, fft_len=fft_padding * data.shape[dim], dim=dim)
+    f, func_as, func_ps = spect(data, fs, fft_len=fft_padding * data.shape[dim], dim=dim)
+
     # correct for coherent gain
     func_as = func_as / coherent_gain
 
@@ -211,99 +207,6 @@ def extract_amp_phase(
     return amp.squeeze(), phase.squeeze(), f[f_index]
 
 
-def brenner_sharpness(im):
-    num_dim = im.ndim
-    if num_dim == 2:
-        # compute metric
-        bren_x = (im[:-2, :] - im[2:, :]) ** 2
-        bren_y = (im[:, :-2] - im[:, 2:]) ** 2
-        s = np.sum(bren_x) + np.sum(bren_y)
-    elif num_dim == 3:
-        # compute metric
-        bren_x = (im[:-2, :, :] - im[2:, :, :]) ** 2
-        bren_y = (im[:, :-2, :] - im[:, 2:, :]) ** 2
-        bren_z = (im[:, :, :-2] - im[:, :, 2:]) ** 2
-        s = np.sum(bren_x) + np.sum(bren_y) + np.sum(bren_z)
-    return s
-
-
-def tenenbaum_sharpness(im):
-    num_dim = im.ndim
-    if num_dim == 2:
-        # define the 2D sobel gradient operator
-        sobel = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
-
-        # compute metric
-        s = (convolve(sobel, im) ** 2 + convolve(sobel.T, im) ** 2).sum()
-    elif num_dim == 3:
-        # define the 3D sobel gradient operator
-        sobel3D = np.zeros((3, 3, 3))
-        sobel3D[:, :, 0] = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]])
-        sobel3D[:, :, 2] = -sobel3D[:, :, 0]
-
-        # compute metric
-        s = (
-            convolve(im, sobel3D) ** 2
-            + convolve(im, np.transpose(sobel3D, (2, 0, 1))) ** 2
-            + convolve(im, np.transpose(sobel3D, (1, 2, 0))) ** 2
-        ).sum()
-    return s
-
-    # TODO: get this passing the tests
-    # NOTE: Walter thinks this is the proper way to do this, but it doesn't match the MATLAB version
-    # num_dim = im.ndim
-    # if num_dim == 2:
-    #     # compute metric
-    #     sx = sobel(im, axis=0, mode='constant')
-    #     sy = sobel(im, axis=1, mode='constant')
-    #     s = (sx ** 2) + (sy ** 2)
-    #     s = np.sum(s)
-    #
-    # elif num_dim == 3:
-    #     # compute metric
-    #     sx = sobel(im, axis=0, mode='constant')
-    #     sy = sobel(im, axis=1, mode='constant')
-    #     sz = sobel(im, axis=2, mode='constant')
-    #     s = (sx ** 2) + (sy ** 2) + (sz ** 2)
-    #     s = np.sum(s)
-    # else:
-    #     raise ValueError("Invalid number of dimensions in im")
-
-
-def sharpness(im: np.ndarray, mode: Optional[str] = "Brenner") -> float:
-    """
-    Returns a scalar metric related to the sharpness of a 2D or 3D image matrix.
-
-    Args:
-        im: The image matrix.
-        metric: The metric to use. Defaults to "Brenner".
-
-    Returns:
-        A scalar sharpness metric.
-
-    Raises:
-        AssertionError: If `im` is not a NumPy array.
-
-    References:
-        B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard, "Automatic sound speed selection in
-        photoacoustic image reconstruction using an autofocus approach," J. Biomed. Opt., vol. 16, no. 9, p. 090501, 2011.
-
-    """
-
-    assert isinstance(im, np.ndarray), "Argument im must be of type numpy array"
-
-    if mode == "Brenner":
-        metric = brenner_sharpness(im)
-    elif mode == "Tenenbaum":
-        metric = tenenbaum_sharpness(im)
-    elif mode == "NormVariance":
-        metric = norm_var(im)
-    else:
-        raise ValueError("Unrecognized sharpness metric passed. Valid values are ['Brenner', 'Tanenbaum', 'NormVariance']")
-
-    return metric
-
-
 def fwhm(f, x):
     """
     fwhm calculates the Full Width at Half Maximum (FWHM) of a positive
@@ -341,7 +244,7 @@ def half_max_x(x, y):
 
 
 def gaussian_filter(
-    signal: Union[np.ndarray, List[float]], Fs: float, frequency: float, bandwidth: float
+    signal: Union[np.ndarray, List[float]], fs: float, frequency: float, bandwidth: float
 ) -> Union[np.ndarray, List[float]]:
     """
     Applies a frequency domain Gaussian filter with the
@@ -351,7 +254,7 @@ def gaussian_filter(
 
     Args:
         signal:         Signal to filter [channel, samples]
-        Fs:             Sampling frequency [Hz]
+        fs:             Sampling frequency [Hz]
         frequency:      Center frequency of filter [Hz]
         bandwidth:      Bandwidth of filter in percentage
 
@@ -362,9 +265,9 @@ def gaussian_filter(
 
     N = signal.shape[-1]
     if N % 2 == 0:
-        f = np.arange(-N / 2, N / 2) * Fs / N
+        f = np.arange(-N / 2, N / 2) * fs / N
     else:
-        f = np.arange(-(N - 1) / 2, (N - 1) / 2 + 1) * Fs / N
+        f = np.arange(-(N - 1) / 2, (N - 1) / 2 + 1) * fs / N
 
     mean = frequency
     variance = (bandwidth / 100 * frequency / (2 * np.sqrt(2 * np.log(2)))) ** 2
@@ -456,7 +359,7 @@ def filter_time_series(
     assert not isinstance(kgrid.t_array, str) or kgrid.t_array != "auto", "kgrid.t_array must be explicitly defined."
 
     # compute the sampling frequency
-    Fs = 1 / kgrid.dt
+    fs = 1 / kgrid.dt
 
     # extract the minimum sound speed
     if medium.sound_speed is not None:
@@ -490,7 +393,7 @@ def filter_time_series(
     if ppw != 0:
         filtered_signal = apply_filter(
             signal,
-            Fs,
+            fs,
             float(filter_cutoff_f),
             "LowPass",
             zero_phase=zerophase,
@@ -532,7 +435,7 @@ def filter_time_series(
 
 def apply_filter(
     signal: np.ndarray,
-    Fs: float,
+    fs: float,
     cutoff_f: float,
     filter_type: str,
     zero_phase: Optional[bool] = False,
@@ -545,7 +448,7 @@ def apply_filter(
 
     Args:
         signal: The input signal.
-        Fs: The sampling frequency of the signal.
+        fs: The sampling frequency of the signal.
         cutoff_f: The cut-off frequency of the filter.
         filter_type: The type of filter to apply, either 'HighPass', 'LowPass' or 'BandPass'.
         zero_phase: Whether to apply a zero-phase filter. Defaults to False.
@@ -564,13 +467,13 @@ def apply_filter(
 
         # apply the low pass filter
         func_filt_lp = apply_filter(
-            signal, Fs, cutoff_f[1], "LowPass", stop_band_atten=stop_band_atten, transition_width=transition_width, zero_phase=zero_phase
+            signal, fs, cutoff_f[1], "LowPass", stop_band_atten=stop_band_atten, transition_width=transition_width, zero_phase=zero_phase
         )
 
         # apply the high pass filter
         filtered_signal = apply_filter(
             func_filt_lp,
-            Fs,
+            fs,
             cutoff_f[0],
             "HighPass",
             stop_band_atten=stop_band_atten,
@@ -584,7 +487,7 @@ def apply_filter(
             high_pass = False
         elif filter_type == "HighPass":
             high_pass = True
-            cutoff_f = Fs / 2 - cutoff_f
+            cutoff_f = fs / 2 - cutoff_f
         else:
             raise ValueError(f'Unknown filter type {filter_type}. Options are "LowPass, HighPass, BandPass"')
 
@@ -602,7 +505,7 @@ def apply_filter(
         N = int(N)
 
         # construct impulse response of ideal bandpass filter h(n), a sinc function
-        fc = cutoff_f / Fs  # normalised cut-off
+        fc = cutoff_f / fs  # normalised cut-off
         n = np.arange(-N / 2, N / 2)
         h = 2 * fc * sinc(2 * np.pi * fc * n)
 

kwave/utils/math.py

@@ -357,22 +357,6 @@ def sinc(x: Union[int, float, np.ndarray]) -> Union[int, float, np.ndarray]:
     return np.sinc(x / np.pi)
 
 
-def norm_var(im: np.ndarray) -> float:
-    """
-    Calculates the normalized variance of an array of values.
-
-    Args:
-        im: The input array.
-
-    Returns:
-        The normalized variance of im.
-
-    """
-    mu = np.mean(im)
-    s = np.sum((im - mu) ** 2) / mu
-    return s
-
-
 def gaussian(
     x: Union[int, float, np.ndarray],
     magnitude: Optional[Union[int, float]] = None,