257 add smoothing functions

eis_toolkit/raster_processing/derivatives/parameters.py

@@ -11,7 +11,8 @@
     NonSquarePixelSizeException,
 )
 from eis_toolkit.raster_processing.derivatives.partial_derivatives import coefficients
-from eis_toolkit.raster_processing.derivatives.utilities import reduce_ndim, scale_raster, set_flat_pixels
+from eis_toolkit.raster_processing.derivatives.utilities import scale_raster, set_flat_pixels
+from eis_toolkit.utilities.miscellaneous import reduce_ndim
 from eis_toolkit.utilities.checks.raster import check_quadratic_pixels
 from eis_toolkit.utilities.conversions import convert_rad_to_deg, convert_rad_to_rise
 from eis_toolkit.utilities.nodata import nan_to_nodata, nodata_to_nan

eis_toolkit/raster_processing/derivatives/utilities.py

@@ -5,22 +5,6 @@
 from beartype.typing import Union
 
 
-@beartype
-def reduce_ndim(
-    data: np.ndarray,
-) -> np.ndarray:
-    """
-    Reduce the number of dimensions of a numpy array.
-
-    Args:
-        data: The input raster data as a numpy array.
-
-    Returns:
-        The reduced array.
-    """
-    return np.squeeze(data) if data.ndim >= 3 else data
-
-
 @beartype
 def scale_raster(
     data: np.ndarray,

eis_toolkit/raster_processing/filters/focal.py

@@ -0,0 +1,149 @@
+from numbers import Number
+
+import numpy as np
+import rasterio
+from beartype import beartype
+from beartype.typing import Literal, Optional
+
+from eis_toolkit.raster_processing.filters.kernels import basic_kernel, gaussian_kernel, mexican_hat_kernel
+from eis_toolkit.raster_processing.filters.utilities import apply_correlated_filter, apply_generic_filter, check_inputs
+from eis_toolkit.utilities.miscellaneous import cast_array_to_float, reduce_ndim
+from eis_toolkit.utilities.nodata import nan_to_nodata, nodata_to_nan
+
+
+@beartype
+def _focal_median(window: np.ndarray) -> Number:
+    weighted_value = np.nanmedian(window) if sum(np.isnan(window)) != len(window) else np.nan
+    return weighted_value
+
+
+@beartype
+def focal_filter(
+    raster: rasterio.io.DatasetReader,
+    method: Literal["mean", "median"] = "mean",
+    size: int = 3,
+    shape: Literal["square", "circle"] = "circle",
+) -> np.ndarray:
+    """
+    Apply a basic focal filter to the input raster.
+
+    Parameters:
+      raster: The input raster dataset.
+      method: The method to use for filtering. Can be either "mean" or "median". Default to "mean".
+      size: The size of the filter window. E.g., 3 means a 3x3 window. Default to 3.
+      shape: The shape of the filter window. Can be either "square" or "circle". Default to "circle".
+
+    Returns:
+      np.ndarray: The filtered raster array.
+
+    Raises:
+      InvalidRasterBandException: If the input raster has more than one band.
+      InvalidParameterValueException: If the filter size is smaller than 3.
+                                      If the filter size is not an odd number.
+                                      If the shape is not "square" or "circle".
+    """
+    check_inputs(raster, size)
+
+    kernel = basic_kernel(size, shape)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    if method == "mean":
+        out_array = apply_correlated_filter(raster_array, kernel)
+    elif method == "median":
+        out_array = apply_generic_filter(raster_array, _focal_median, kernel)
+
+    out_array = nan_to_nodata(out_array, raster.nodata)
+    return cast_array_to_float(out_array, cast_float=True)
+
+
+@beartype
+def gaussian_filter(
+    raster: rasterio.io.DatasetReader,
+    sigma: Number = 1,
+    truncate: Number = 4,
+    size: Optional[int] = None,
+) -> np.ndarray:
+    """
+    Apply a gaussian filter to the input raster.
+
+    Parameters:
+      raster: The input raster dataset.
+      sigma: The standard deviation of the gaussian kernel.
+      truncate: The truncation factor for the gaussian kernel based on the sigma value.
+                Only if size is not given. Default to 4.0.
+                E.g., for sigma = 1 and truncate = 4.0, the kernel size is 9x9.
+      size: The size of the filter window. E.g., 3 means a 3x3 window.
+            If size is not None, it overrides the dynamic size calculation based on sigma and truncate.
+            Default to None.
+
+    Returns:
+      np.ndarray: The filtered raster array.
+
+    Raises:
+      InvalidRasterBandException: If the input raster has more than one band.
+      InvalidParameterValueException: If the filter size is smaller than 3.
+                                      If the filter size is not an odd number.
+                                      If the resulting radius is smaller than 1.
+    """
+    check_inputs(raster, size, sigma, truncate)
+
+    kernel = gaussian_kernel(sigma, truncate, size)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    out_array = apply_correlated_filter(raster_array, kernel)
+    out_array = nan_to_nodata(out_array, raster.nodata)
+
+    return cast_array_to_float(out_array, cast_float=True)
+
+
+@beartype
+def mexican_hat_filter(
+    raster: rasterio.io.DatasetReader,
+    sigma: Number = 1,
+    truncate: Number = 4,
+    size: Optional[int] = None,
+    direction: Literal["rectangular", "circular"] = "circular",
+) -> np.ndarray:
+    """
+    Apply a mexican hat filter to the input raster.
+
+    Circular: Lowpass filter for smoothing.
+    Rectangular: Highpass filter for edge detection. Results may need further normalization.
+
+    Parameters:
+      raster: The input raster dataset.
+      sigma: The standard deviation.
+      truncate: The truncation factor.
+                E.g., for sigma = 1 and truncate = 4.0, the kernel size is 9x9.
+                Default to 4.0.
+      size: The size of the filter window. E.g., 3 means a 3x3 window. Default to None.
+      direction: The direction of calculating the kernel values.
+                 Can be either "rectangular" or "circular". Default to "circular".
+
+    Returns:
+      np.ndarray: The filtered raster array.
+
+    Raises:
+      InvalidRasterBandException: If the input raster has more than one band.
+      InvalidParameterValueException: If the filter size is smaller than 3.
+                                      If the filter size is not an odd number.
+                                      If the resulting radius is smaller than 1.
+    """
+    check_inputs(raster, size, sigma, truncate)
+
+    kernel = mexican_hat_kernel(sigma, truncate, size, direction)
+
+    raster_array = raster.read()
+    raster_array = reduce_ndim(raster_array)
+    raster_array = nodata_to_nan(raster_array, raster.nodata)
+
+    out_array = apply_correlated_filter(raster_array, kernel)
+    out_array = nan_to_nodata(out_array, raster.nodata)
+
+    return cast_array_to_float(out_array, cast_float=True)

eis_toolkit/raster_processing/filters/kernels.py

@@ -0,0 +1,121 @@
+import numpy as np
+from beartype import beartype
+from beartype.typing import Literal
+from scipy.signal import ricker
+
+
+@beartype
+def get_kernel_size(sigma, truncate, size):
+    """
+    Calculate the kernel size and radius based on the given parameters.
+
+    Args:
+        sigma: The standard deviation.
+        truncate: The truncation factor.
+        size: The optional size of the kernel.
+
+    Returns:
+        A tuple containing the calculated size and radius of the kernel.
+    """
+    if size is not None:
+        radius = int(size // 2)
+    else:
+        radius = int(float(truncate) * float(sigma) + 0.5)
+        size = int(2 * radius + 1)
+
+    return size, radius
+
+
+@beartype
+def create_grid(radius, size):
+    """
+    Create a grid of coordinates.
+
+    Args:
+        radius: The radius of the grid.
+        size: The size of the grid.
+
+    Returns:
+        Tuple with x and y coordinates of the grid.
+    """
+    y, x = np.ogrid[-radius : (size - radius), -radius : (size - radius)]
+    return x, y
+
+
+@beartype
+def basic_kernel(size: int, shape: Literal["square", "circle"]) -> np.ndarray:
+    """
+    Generate a basic kernel of a specified size and shape.
+
+    Args:
+        size: The size of the kernel.
+        shape: The shape of the kernel. Can be either "square" or "circle".
+
+    Returns:
+        np.ndarray: The generated kernel.
+    """
+    if shape == "square":
+        kernel = np.ones((size, size))
+    elif shape == "circle":
+        radius = int(size // 2)
+        x, y = create_grid(radius, size)
+        mask = x**2 + y**2 <= radius**2
+        kernel = np.zeros((size, size))
+        kernel[mask] = 1
+
+    return kernel
+
+
+@beartype
+def gaussian_kernel(sigma, truncate, size):
+    """
+    Generate a Gaussian kernel for image denoising.
+
+    Args:
+        sigma: Standard deviation.
+        truncate: Truncation value.
+        size: Size of the kernel. If not provided, it will be calculated based on sigma and truncate.
+
+    Returns:
+        The Gaussian kernel.
+    """
+    size, radius = get_kernel_size(sigma, truncate, size)
+
+    x, y = create_grid(radius, size)
+    kernel = 1 / (2 * np.pi * sigma**2) * np.exp((x**2 + y**2) / (2 * sigma**2) * -1)
+    kernel /= np.max(kernel)
+
+    return kernel
+
+
+@beartype
+def mexican_hat_kernel(sigma, truncate, size, direction):
+    """
+    Generate a Mexican Hat kernel for denoising (circular) or edge detection (rectangular).
+
+    Args:
+        sigma: The standard deviation of the Gaussian function.
+        truncate: The truncation factor for the Gaussian function.
+        size: The size of the kernel.
+        direction: The direction of the kernel, either "rectangular" or "circular".
+
+    Returns:
+        The Mexican Hat kernel.
+    """
+    size, radius = get_kernel_size(sigma, truncate, size)
+
+    if direction == "rectangular":
+        ricker_wavelet = ricker(size, sigma)
+        kernel = np.outer(ricker_wavelet, ricker_wavelet)
+    elif direction == "circular":
+        x, y = create_grid(radius, size)
+        kernel = (
+            2
+            / np.sqrt(3 * sigma)
+            * np.pi ** (1 / 4)
+            * (1 - (x**2 + y**2) / sigma**2)
+            * np.exp(-(x**2 + y**2) / (2 * sigma**2))
+        )
+    kernel /= np.max(kernel)
+
+    return kernel