include target snr

Author: RafaelSdeSouza

R/segment.R

@@ -17,6 +17,19 @@
 #' @param redshift Numeric, the redshift to apply for wavelength correction. Defaults to 0 (no correction).
 #' @param scale_fn A function used to scale each row of the 2D representation of the data cube.
 #'   Defaults to \code{\link[base]{scale}}. If you have a custom scaling function, pass it here.
+#' @param target_snr Optional minimum accepted SNR per cluster. When supplied,
+#'   Capivara chooses the largest number of clusters whose minimum cluster SNR
+#'   remains above this threshold.
+#' @param var_cube Optional variance cube matching the input cube. Used only
+#'   when \code{target_snr} is supplied.
+#' @param k_values Optional candidate cluster counts tested when
+#'   \code{target_snr} is supplied.
+#' @param wavelength_range Optional wavelength interval used to compute SNR when
+#'   \code{target_snr} is supplied.
+#' @param snr_stat Either integrated SNR or median per-wavelength SNR when
+#'   \code{target_snr} is supplied.
+#' @param variance_inflation Multiplicative factor applied to propagated
+#'   variances when \code{target_snr} is supplied.
 #'
 #' @details
 #' Steps performed by the function:
@@ -26,7 +39,9 @@
 #'   \item Scales the data row-wise using \code{scale_fn}.
 #'   \item Computes pairwise distances between rows using \code{\link{torch_dist}}.
 #'   \item Performs hierarchical clustering using Ward's D2 method via \code{\link[fastcluster]{hclust}}.
-#'   \item Cuts the dendrogram into \code{Ncomp} clusters and reshapes the results into a 2D cluster map.
+#'   \item Cuts the dendrogram into \code{Ncomp} clusters and reshapes the results into a 2D cluster map,
+#'   or, when \code{target_snr} is supplied, chooses the largest cut whose minimum cluster SNR
+#'   remains above the requested threshold.
 #'   \item Calculates the signal-to-noise ratio (SNR) for each cluster.
 #' }
 #'
@@ -56,8 +71,35 @@
 #' }
 #'
 #' @export
-# Cluster the IFU cube data
-segment <- function(input, Ncomp = 5, redshift = 0, scale_fn = median_scale) {
+segment <- function(input,
+                    Ncomp = 5,
+                    redshift = 0,
+                    scale_fn = median_scale,
+                    target_snr = NULL,
+                    var_cube = NULL,
+                    k_values = NULL,
+                    wavelength_range = NULL,
+                    snr_stat = c("integrated", "median_per_wavelength"),
+                    variance_inflation = 1) {
+  if (!is.null(target_snr)) {
+    if (!missing(Ncomp)) {
+      stop("Specify either `Ncomp` or `target_snr`, not both.")
+    }
+
+    return(choose_ncomp_by_snr(
+      input = input,
+      target_snr = target_snr,
+      var_cube = var_cube,
+      k_values = k_values,
+      wavelength_range = wavelength_range,
+      redshift = redshift,
+      scale_fn = scale_fn,
+      snr_stat = snr_stat,
+      variance_inflation = variance_inflation,
+      na_safe = TRUE
+    ))
+  }
+
   .segment_core(
     input = input,
     Ncomp = Ncomp,

R/segment_big_cube.R

@@ -183,8 +183,8 @@
 #' @param refine_max_pixels Maximum number of boundary pixels to refine.
 #' @param verbose Logical.
 #' @return A list containing the full-resolution \code{cluster_map}, block
-#'   metadata, medoid centers, axis/header data, the original cube, and
-#'   \code{backend = "medoid"}.
+#'   metadata, medoid centers, a per-cluster SNR summary, axis/header data, the
+#'   original cube, and \code{backend = "medoid"}.
 #' @seealso \code{\link{segment}}, \code{\link{segment_starlet}}
 #' @export
 segment_big_cube <- function(input,
@@ -391,8 +391,16 @@ segment_big_cube <- function(input,
     }
   }
 
+  sn <- .compute_signal_noise(IFU2D)
+  cluster_snr <- .compute_cluster_snr(
+    clusters = cluster_map[valid_pix],
+    signal_valid = sn$signal[valid_pix],
+    noise_valid = sn$noise[valid_pix]
+  )
+
   list(
     cluster_map = cluster_map,
+    cluster_snr = cluster_snr,
     block_map = block_map,
     block_size = block_size,
     n_blocks = n_blocks,

R/segment_starlet.R

@@ -14,6 +14,19 @@
 #' @param engine Segmentation backend after masking: \code{"standard"} uses
 #'   \code{\link{segment}}, while \code{"big_cube"} uses
 #'   \code{\link{segment_big_cube}}.
+#' @param target_snr Optional minimum accepted SNR per cluster for the exact
+#'   backend. When supplied, the masked cube is cut at the largest
+#'   \code{Ncomp} whose minimum cluster SNR remains above this threshold.
+#' @param var_cube Optional variance cube used only when \code{target_snr} is
+#'   supplied with \code{engine = "standard"}.
+#' @param k_values Optional candidate cluster counts tested only when
+#'   \code{target_snr} is supplied with \code{engine = "standard"}.
+#' @param wavelength_range Optional wavelength interval used to compute SNR only
+#'   when \code{target_snr} is supplied with \code{engine = "standard"}.
+#' @param snr_stat Either integrated SNR or median per-wavelength SNR when
+#'   \code{target_snr} is supplied with \code{engine = "standard"}.
+#' @param variance_inflation Multiplicative factor applied to propagated
+#'   variances when \code{target_snr} is supplied with \code{engine = "standard"}.
 #' @param collapse_fn Function used to build the white-light image.
 #' @param starlet_J Number of starlet scales.
 #' @param starlet_scales Integer vector of scales kept in the reconstruction.
@@ -33,6 +46,12 @@ segment_starlet <- function(input,
                             redshift = 0,
                             scale_fn = median_scale,
                             engine = c("standard", "big_cube"),
+                            target_snr = NULL,
+                            var_cube = NULL,
+                            k_values = NULL,
+                            wavelength_range = NULL,
+                            snr_stat = c("integrated", "median_per_wavelength"),
+                            variance_inflation = 1,
                             collapse_fn = collapse_white_light,
                             starlet_J = 5,
                             starlet_scales = 2:5,
@@ -45,6 +64,11 @@ segment_starlet <- function(input,
   engine <- match.arg(engine)
   mode <- match.arg(mode)
   mask_mode <- match.arg(mask_mode)
+  snr_stat <- match.arg(snr_stat)
+
+  if (engine == "big_cube" && !is.null(target_snr)) {
+    stop("`target_snr` is currently supported only with `engine = \"standard\"`.")
+  }
 
   cubedat <- .as_cubedat(input)
   mask_info <- build_starlet_mask(
@@ -74,12 +98,26 @@ segment_starlet <- function(input,
       ...
     )
   } else {
-    segment(
-      input = masked_input,
-      Ncomp = Ncomp,
-      redshift = redshift,
-      scale_fn = scale_fn
-    )
+    if (is.null(target_snr)) {
+      segment(
+        input = masked_input,
+        Ncomp = Ncomp,
+        redshift = redshift,
+        scale_fn = scale_fn
+      )
+    } else {
+      segment(
+        input = masked_input,
+        redshift = redshift,
+        scale_fn = scale_fn,
+        target_snr = target_snr,
+        var_cube = var_cube,
+        k_values = k_values,
+        wavelength_range = wavelength_range,
+        snr_stat = snr_stat,
+        variance_inflation = variance_inflation
+      )
+    }
   }
 
   c(out, list(