R/spat_filter.R

In neuroim2: Data Structures for Brain Imaging Data

#' @include all_class.R
NULL
#' @include all_generic.R
NULL

#' Spatial Filtering Methods for Neuroimaging Data
#'
#' @name spatial-filter
#' @description Methods for applying spatial filters to neuroimaging data
NULL

#' @importFrom methods new
#' @importFrom stats dnorm
NULL

#' Gaussian Blur for Volumetric Images
#'
#' @description
#' This function applies an isotropic discrete Gaussian kernel to smooth a volumetric image (3D brain MRI data).
#' The blurring is performed within a specified image mask, with customizable kernel parameters.
#'
#' @param vol A \code{\linkS4class{NeuroVol}} object representing the image volume to be smoothed.
#' @param mask An optional \code{\linkS4class{LogicalNeuroVol}} object representing the image mask.
#'   This mask defines the region where the blurring is applied. If not provided, the entire volume is processed.
#' @param sigma A numeric value specifying the standard deviation of the Gaussian kernel. Default is 2.
#' @param window An integer specifying the kernel size. It represents the number of voxels to include
#'   on each side of the center voxel. For example, window=1 results in a 3x3x3 kernel. Default is 1.
#'
#' @return A \code{\linkS4class{NeuroVol}} object representing the smoothed image.
#'
#' @details
#' The function uses a C++ implementation for efficient Gaussian blurring. The blurring is applied
#' only to voxels within the specified mask (or the entire volume if no mask is provided).
#' The kernel size is determined by the 'window' parameter, and its shape by the 'sigma' parameter.
#'
#' @examples
#' # Load a sample brain mask
#' brain_mask <- read_vol(system.file("extdata", "global_mask_v4.nii", package = "neuroim2"))
#'
#' # Apply Gaussian blurring to the brain volume
#' blurred_vol <- gaussian_blur(brain_mask, brain_mask, sigma = 2, window = 1)
#'
#' # View a slice of the original and blurred volumes
#' image(brain_mask[,,12])
#' image(blurred_vol[,,12])
#'
#' @seealso
#' \code{\link{NeuroVol-class}}, \code{\link{LogicalNeuroVol-class}}, \code{\link{bilateral_filter}}
#'
#' @references
#' Gaussian blur: https://en.wikipedia.org/wiki/Gaussian_blur
#'
#' @export
gaussian_blur <- function(vol, mask, sigma = 2, window = 1) {
  if (!inherits(vol, "NeuroVol")) {
    cli::cli_abort("{.arg vol} must be a {.cls NeuroVol} object.")
  }
  if (window < 1) {
    cli::cli_abort("{.arg window} must be >= 1, not {.val {window}}.")
  }
  if (sigma <= 0) {
    cli::cli_abort("{.arg sigma} must be positive, not {.val {sigma}}.")
  }
  if (!missing(mask)) {
    if (!inherits(mask, "NeuroVol")) {
      cli::cli_abort("{.arg mask} must be a {.cls NeuroVol} object.")
    }
  }

  if (missing(mask)) {
    mask.idx <- seq_len(prod(dim(vol)))
    target_space <- space(vol)
  } else {
    mask.idx <- which(mask != 0)
    target_space <- space(mask)
  }

  arr <- as.array(vol)
  farr <- gaussian_blur_cpp(arr, as.integer(mask.idx), as.integer(window), sigma, spacing(vol))

  out <- NeuroVol(farr, target_space)
  out
}

#' Edge-Preserving Guided Filter for Volumetric Images
#'
#' @description
#' This function applies a guided filter to a volumetric image (3D brain MRI data)
#' to perform edge-preserving smoothing. The guided filter smooths the image while
#' preserving edges, providing a balance between noise reduction and structural preservation.
#'
#' @param vol A \code{\linkS4class{NeuroVol}} object representing the image volume to be filtered.
#' @param radius An integer specifying the spatial radius of the filter. Default is 4.
#' @param epsilon A numeric value specifying the regularization parameter. It controls
#'   the degree of smoothing and edge preservation. Default is 0.49 (0.7^2).
#'
#' @return A \code{\linkS4class{NeuroVol}} object representing the filtered image.
#'
#' @details
#' The guided filter operates by computing local linear models between the guidance
#' image (which is the same as the input image in this implementation) and the output.
#' The 'radius' parameter determines the size of the local neighborhood, while 'epsilon'
#' controls the smoothness of the filter.
#'
#' The implementation uses box blur operations for efficiency, which approximates
#' the behavior of the original guided filter algorithm.
#'
#' @examples
#' # Load an example brain volume
#' brain_vol <- read_vol(system.file("extdata", "global_mask_v4.nii", package = "neuroim2"))
#'
#' # Apply guided filtering to the brain volume
#' \donttest{
#' filtered_vol <- guided_filter(brain_vol, radius = 4, epsilon = 0.49)
#'
#' # Visualize a slice of the original and filtered volumes
#' oldpar <- par(mfrow = c(1, 2))
#' image(brain_vol[,,12], main = "Original")
#' image(filtered_vol[,,12], main = "Filtered")
#' par(oldpar)
#' }
#'
#' @references
#' He, K., Sun, J., & Tang, X. (2013). Guided Image Filtering. IEEE Transactions
#' on Pattern Analysis and Machine Intelligence, 35(6), 1397-1409.
#'
#' @seealso
#' \code{\link{gaussian_blur}}, \code{\link{bilateral_filter}}, \code{\link{NeuroVol-class}}
#'
#' @export
guided_filter <- function(vol, radius = 4, epsilon = 0.7^2) {
  if (!inherits(vol, "NeuroVol")) {
    cli::cli_abort("{.arg vol} must be a {.cls NeuroVol} object.")
  }
  if (radius < 1) {
    cli::cli_abort("{.arg radius} must be >= 1, not {.val {radius}}.")
  }
  if (epsilon <= 0) {
    cli::cli_abort("{.arg epsilon} must be positive, not {.val {epsilon}}.")
  }

  mask_idx <- which(vol !=0)
  mean_I = box_blur(vol, mask_idx, radius)
  mean_II = box_blur(vol*vol, mask_idx, radius)
  var_I = mean_II - mean_I * mean_I
  mean_p = box_blur(vol, mask_idx, radius)
  mean_Ip = box_blur(vol*vol, mask_idx, radius)

  cov_Ip = mean_Ip - mean_I * mean_p
  a = cov_Ip / (var_I + epsilon)
  b = mean_p - a * mean_I
  mean_a = box_blur(a, mask_idx, radius)
  mean_b = box_blur(b, mask_idx, radius)
  out = mean_a * vol + mean_b
  ovol = NeuroVol(out, space(vol))
  ovol
}

#' Apply a bilateral filter to a volumetric image
#'
#' This function smooths a volumetric image (3D brain MRI data) using a bilateral filter.
#' The bilateral filter considers both spatial closeness and intensity similarity for smoothing.
#'
#' @param vol A \code{\linkS4class{NeuroVol}} object representing the image volume to be smoothed.
#' @param mask An optional \code{\linkS4class{LogicalNeuroVol}} object representing the image mask that defines the region where the filtering is applied. If not provided, the entire volume is considered.
#' @param spatial_sigma A numeric value specifying the standard deviation of the spatial Gaussian kernel (default is 2).
#' @param intensity_sigma A numeric value specifying the standard deviation of the intensity Gaussian kernel (default is 25).
#' @param window An integer specifying the number of voxels around the center voxel to include on each side. For example, window=1 for a 3x3x3 kernel (default is 1).
#'
#' @return A smoothed image of class \code{\linkS4class{NeuroVol}}.
#'
#' @examples
#' brain_mask <- read_vol(system.file("extdata", "global_mask_v4.nii", package="neuroim2"))
#'
#' # Apply bilateral filtering to the brain volume
#' filtered_vol <- bilateral_filter(brain_mask, brain_mask, spatial_sigma = 2,
#' intensity_sigma = 25, window = 1)
#'
#' @export
bilateral_filter <- function(vol, mask, spatial_sigma=2, intensity_sigma=1, window=1) {
  if (window < 1) {
    cli::cli_abort("{.arg window} must be >= 1, not {.val {window}}.")
  }
  if (spatial_sigma <= 0) {
    cli::cli_abort("{.arg spatial_sigma} must be positive, not {.val {spatial_sigma}}.")
  }
  if (intensity_sigma <= 0) {
    cli::cli_abort("{.arg intensity_sigma} must be positive, not {.val {intensity_sigma}}.")
  }
  if (!missing(mask)) {
    if (!inherits(mask, "NeuroVol")) {
      cli::cli_abort("{.arg mask} must be a {.cls NeuroVol} object.")
    }
  }

  if (missing(mask)) {
    mask.idx <- seq_len(prod(dim(vol)))
    target_space <- space(vol)
  } else {
    mask.idx <- which(mask!=0)
    target_space <- space(mask)
  }

  arr <- as.array(vol)
  farr <- bilateral_filter_cpp(arr, as.integer(mask.idx), as.integer(window), spatial_sigma, intensity_sigma, spacing(vol))

  out <- NeuroVol(farr, target_space)
  out
}

#' Apply a bilateral filter to each volume of a NeuroVec
#'
#' This function applies a bilateral filter to each volume of a NeuroVec object.
#' The filter is applied using a specified spatial and intensity sigma, and a given window size.
#'
#' @param vec A NeuroVec object containing the volumes to be filtered.
#' @param mask A binary mask specifying the region of interest. If not provided, the whole volume is considered.
#' @param spatial_sigma The spatial sigma for the bilateral filter (default = 2).
#' @param intensity_sigma The intensity sigma for the bilateral filter (default = 1).
#' @param window The size of the window for the bilateral filter (default = 1).
#' @return A NeuroVec object with the filtered volumes.
#' @examples
#' brain_mask <- read_vol(system.file("extdata", "global_mask_v4.nii", package="neuroim2"))
#' vec <- read_vec(system.file("extdata", "global_mask_v4.nii", package="neuroim2"))
#' out <- bilateral_filter_vec(vec,brain_mask)
#' @noRd
bilateral_filter_vec <- function(vec, mask, spatial_sigma=2, intensity_sigma=1, window=1) {
  if (!inherits(vec, "NeuroVec")) {
    cli::cli_abort("{.arg vec} must be a {.cls NeuroVec} object.")
  }
  if (window < 1) {
    cli::cli_abort("{.arg window} must be >= 1, not {.val {window}}.")
  }
  if (spatial_sigma <= 0) {
    cli::cli_abort("{.arg spatial_sigma} must be positive, not {.val {spatial_sigma}}.")
  }
  if (intensity_sigma <= 0) {
    cli::cli_abort("{.arg intensity_sigma} must be positive, not {.val {intensity_sigma}}.")
  }

  if (missing(mask)) {
    mask.idx <- seq_len(prod(dim(vec)[1:3]))
    target_space <- space(vec[[1]])
  } else {
    if (!inherits(mask, "NeuroVol")) {
      cli::cli_abort("{.arg mask} must be a {.cls NeuroVol} object.")
    }
    mask.idx <- which(mask!=0)
    target_space <- space(mask)
  }

  res<- lapply(seq_len(dim(vec)[4]), function(i) {
    vol_i <- vec[[i]]
    arr <- as.array(vol_i)
    farr <- bilateral_filter_cpp(arr, as.integer(mask.idx), as.integer(window), spatial_sigma, intensity_sigma, spacing(vec)[1:3])
    NeuroVol(farr, target_space)
  })

  do.call(concat,res)

}

#' Apply a 4D bilateral filter to a NeuroVec
#'
#' This function applies a full 4D bilateral filter to a \code{NeuroVec},
#' smoothing jointly across space (x, y, z) and time (t). The filter uses
#' spatial, temporal, and intensity kernels to preserve edges while reducing
#' noise, leveraging a parallel C++ backend for performance.
#'
#' @param vec A \code{\linkS4class{NeuroVec}} object (4D image).
#' @param mask An optional \code{\linkS4class{LogicalNeuroVol}} or \code{\linkS4class{NeuroVol}}
#'   specifying the spatial region to process. If omitted, the entire spatial
#'   extent is processed.
#' @param spatial_sigma Numeric; standard deviation of the spatial Gaussian (default 2).
#' @param intensity_sigma Numeric; standard deviation of the intensity Gaussian (default 1).
#' @param temporal_sigma Numeric; standard deviation of the temporal Gaussian (default 1).
#' @param spatial_window Integer; half-width of the spatial window in voxels (default 1),
#'   e.g., 1 => 3x3x3 spatial neighborhood.
#' @param temporal_window Integer; half-width of the temporal window in frames (default 1),
#'   e.g., 1 => 3 timepoints (t-1, t, t+1).
#' @param temporal_spacing Numeric; spacing of the temporal dimension (e.g., TR in seconds).
#'   Default is 1. This sets the temporal scale used for the temporal kernel.
#'
#' @details
#' Parameter guidance and units:
#' - spatial_sigma: Measured in physical units (millimeters). Distances are
#'   computed using \code{spacing(vec)[1:3]}, so choose \code{spatial_sigma}
#'   relative to voxel size. As a rule of thumb, set it to about 1-2 voxel sizes
#'   (e.g., 2-4 mm for 2 mm isotropic data) for moderate smoothing.
#' - intensity_sigma: Dimensionless multiplier of the global intensity standard
#'   deviation. Internally, the filter uses exp(-(dI)^2 / (2 * (intensity_sigma * sigma_I)^2)),
#'   where sigma_I is the standard deviation of all finite voxel intensities within
#'   the mask across time. Start with 1.0 for moderate smoothing; use 0.5-0.8 to
#'   preserve more edges, or 1.5-2.0 for stronger smoothing.
#' - temporal_sigma: Measured in \code{temporal_spacing} units (e.g., seconds).
#'   Typical values are 0.5-2 x TR. Larger values blend more across time.
#'
#' Choosing the neighborhood window sizes:
#' - spatial_window controls the discrete spatial support. A common choice is
#'   \code{ceiling(2 * spatial_sigma / min(spacing(vec)[1:3]))}, which covers
#'   ~95% of a Gaussian's mass.
#' - temporal_window similarly can be set to \code{ceiling(2 * temporal_sigma / temporal_spacing)}.
#'
#' Quick presets (typical fMRI with 2-3 mm voxels and TR~2s):
#' - Light: spatial_sigma = 1 x min(spacing), intensity_sigma = 0.8,
#'   temporal_sigma = 0.5 x TR, windows = 1
#' - Moderate (default-ish): spatial_sigma = 1.5 x min(spacing), intensity_sigma = 1.0,
#'   temporal_sigma = 1 x TR, windows = 1-2
#' - Strong: spatial_sigma = 2 x min(spacing), intensity_sigma = 1.5,
#'   temporal_sigma = 1.5 x TR, windows = 2
#'
#' Tip: If your time axis has known TR, pass it via \code{temporal_spacing}.
#' For NIfTI inputs, you can get TR via:
#' \preformatted{
#'   hdr <- read_header(nifti_path)
#'   tr  <- hdr@header$pixdim[5]
#'   out <- bilateral_filter_4d(vec, mask, temporal_spacing = tr)
#' }
#'
#' @return A \code{\linkS4class{NeuroVec}} with filtered data.
#'
#' @examples
#' \donttest{
#' vec <- read_vec(system.file("extdata", "global_mask_v4.nii", package = "neuroim2"))
#' mask <- read_vol(system.file("extdata", "global_mask_v4.nii", package = "neuroim2"))
#' out  <- bilateral_filter_4d(vec, mask,
#'                             spatial_sigma = 2, intensity_sigma = 1,
#'                             temporal_sigma = 1, spatial_window = 1,
#'                             temporal_window = 1, temporal_spacing = 1)
#' }
#'
#' @seealso \code{\link{bilateral_filter}}, \code{\link{NeuroVec-class}}, \code{\link{NeuroVol-class}}
#' @export
bilateral_filter_4d <- function(vec,
                                mask,
                                spatial_sigma = 2,
                                intensity_sigma = 1,
                                temporal_sigma = 1,
                                spatial_window = 1,
                                temporal_window = 1,
                                temporal_spacing = 1) {

  if (!inherits(vec, "NeuroVec")) {
    cli::cli_abort("{.arg vec} must be a {.cls NeuroVec} object.")
  }
  if (!is.numeric(spatial_window) || spatial_window < 1) {
    cli::cli_abort("{.arg spatial_window} must be a numeric value >= 1, not {.val {spatial_window}}.")
  }
  if (!is.numeric(temporal_window) || temporal_window < 0) {
    cli::cli_abort("{.arg temporal_window} must be a numeric value >= 0, not {.val {temporal_window}}.")
  }
  if (spatial_sigma <= 0) {
    cli::cli_abort("{.arg spatial_sigma} must be positive, not {.val {spatial_sigma}}.")
  }
  if (intensity_sigma <= 0) {
    cli::cli_abort("{.arg intensity_sigma} must be positive, not {.val {intensity_sigma}}.")
  }
  if (temporal_sigma <= 0) {
    cli::cli_abort("{.arg temporal_sigma} must be positive, not {.val {temporal_sigma}}.")
  }
  if (temporal_spacing <= 0) {
    cli::cli_abort("{.arg temporal_spacing} must be positive, not {.val {temporal_spacing}}.")
  }

  # Determine mask and target space
  if (missing(mask)) {
    mask.idx <- seq_len(prod(dim(vec)[1:3]))
    target_space <- space(vec)
  } else {
    if (!inherits(mask, "NeuroVol") && !inherits(mask, "LogicalNeuroVol")) {
      cli::cli_abort("{.arg mask} must be a {.cls NeuroVol} or {.cls LogicalNeuroVol} object.")
    }
    if (!all(dim(mask) == dim(vec)[1:3])) {
      cli::cli_abort("{.arg mask} spatial dimensions {.val {dim(mask)}} must match spatial dims of {.arg vec} {.val {dim(vec)[1:3]}}.")
    }
    if (!all(spacing(mask) == spacing(vec)[1:3])) {
      cli::cli_abort("{.arg mask} and {.arg vec} must have identical spatial spacing.")
    }
    mask.idx <- which(mask != 0)
    # Keep original 4D space
    target_space <- space(vec)
  }

  # Assemble spacing with temporal component for the 4D kernel
  sp4 <- c(spacing(vec)[1:3], temporal_spacing)

  arr <- as.array(vec)
  farr <- bilateral_filter_4d_cpp_par(arr,
                                       as.integer(mask.idx),
                                       as.integer(spatial_window),
                                       as.integer(temporal_window),
                                       spatial_sigma,
                                       intensity_sigma,
                                       temporal_sigma,
                                       sp4)

  DenseNeuroVec(farr, target_space)
}

#' Laplacian Enhancement Filter for Volumetric Images
#'
#' @description
#' This function applies a multi-layer Laplacian enhancement filter to a volumetric image (3D brain MRI data).
#' The filter enhances details while preserving edges using a non-local means approach with multiple scales.
#'
#' @param vol A \code{\linkS4class{NeuroVol}} object representing the image volume to be enhanced.
#' @param mask A \code{\linkS4class{LogicalNeuroVol}} object specifying the region to process. If not provided,
#'   the entire volume will be processed.
#' @param k An integer specifying the number of layers in the decomposition (default is 2).
#' @param patch_size An integer specifying the size of patches for non-local means. Must be odd (default is 3).
#' @param search_radius An integer specifying the radius of the search window (default is 2).
#' @param h A numeric value controlling the filtering strength. Higher values mean more smoothing (default is 0.7).
#' @param mapping_params An optional list of parameters for the enhancement mappings.
#' @param use_normalization_free Logical indicating whether to use normalization-free weights (default is TRUE).
#'
#' @return A \code{\linkS4class{NeuroVol}} object representing the enhanced image.
#'
#' @export
laplace_enhance <- function(vol, mask, k = 2, patch_size = 3, search_radius = 2,
                          h = 0.7, mapping_params = NULL,
                          use_normalization_free = TRUE) {

  if (!inherits(vol, "NeuroVol")) {
    cli::cli_abort("{.arg vol} must be a {.cls NeuroVol} object.")
  }
  if (k < 1) {
    cli::cli_abort("{.arg k} must be >= 1, not {.val {k}}.")
  }
  if (patch_size < 3 || patch_size %% 2 != 1) {
    cli::cli_abort("{.arg patch_size} must be an odd integer >= 3, not {.val {patch_size}}.")
  }
  if (search_radius < 1) {
    cli::cli_abort("{.arg search_radius} must be >= 1, not {.val {search_radius}}.")
  }
  if (h <= 0) {
    cli::cli_abort("{.arg h} must be positive, not {.val {h}}.")
  }

  # Create default mask if not provided
  if (missing(mask)) {
    mask <- LogicalNeuroVol(array(TRUE, dim(vol)), space(vol))
  } else {
    if (!inherits(mask, "LogicalNeuroVol")) {
      cli::cli_abort("{.arg mask} must be a {.cls LogicalNeuroVol} object.")
    }
  }

  # Call C++ implementation
  farr <- fast_multilayer_laplacian_enhancement_masked(
    as.array(vol),
    as.logical(mask),
    as.integer(k),
    as.integer(patch_size),
    as.integer(search_radius),
    h,
    mapping_params,
    use_normalization_free
  )

  # Return enhanced volume
  out <- NeuroVol(farr, space(vol))
  out
}