new documentation

Author: RafaelSdeSouza

NAMESPACE

@@ -9,6 +9,8 @@ export(cube_to_matrix)
 export(median_scale)
 export(plot_cluster)
 export(plot_cluster_spectra)
+export(reconstruct_cluster_cube)
+export(reconstruct_flux_preserving_cube)
 export(segment)
 export(segment_approx)
 export(segment_blockward)

NEWS.md

@@ -8,6 +8,7 @@ Released: 2026-03-18
 - `build_starlet_mask()` and `segment_starlet()` for Sagui-style white-light starlet masking before clustering.
 - `summarize_cluster_spectra()` for median, summed, and inverse-variance-weighted cluster spectra.
 - `choose_ncomp_by_snr()` for variance-aware component selection from an SNR threshold.
+- `reconstruct_cluster_cube()` and `reconstruct_flux_preserving_cube()` for representative and flux-preserving model cubes.
 - `scripts/run_manga_starlet_comparison.R` for a full-frame MaNGA comparison workflow.
 
 ## Changed

R/SpecMean.R

@@ -24,6 +24,10 @@
 #'   \item Generates a table linking each pixel's original spectrum to its assigned cluster.
 #' }
 #'
+#' This legacy helper is retained for compatibility. For new workflows, prefer
+#' \code{\link{reconstruct_cluster_cube}} or
+#' \code{\link{reconstruct_flux_preserving_cube}}.
+#'
 #' @return A list containing:
 #' \describe{
 #'   \item{\code{MeanSpec}}{A data frame with the median spectrum for each cluster (rows = clusters, columns = wavelength).}
@@ -58,6 +62,13 @@ SpecMean <- function(cluster_result) {
   cluster_map <- cluster_result$cluster_map
   summary <- summarize_cluster_spectra(cluster_result)
   wavelength <- summary$wavelength
+  recon <- reconstruct_cluster_cube(
+    cluster_result = cluster_result,
+    template = "median",
+    preserve_mask = FALSE,
+    fill_mode = "na",
+    return_residual = FALSE
+  )
 
   IFU2D <- cube_to_matrix(cubedat)
   class2D <- as.vector(cluster_map)
@@ -72,19 +83,7 @@ SpecMean <- function(cluster_result) {
     check.names = FALSE
   )
 
-  median_cube <- matrix(NA_real_, nrow = nrow(IFU2D), ncol = ncol(IFU2D))
-  for (i in seq_along(summary$cluster_ids)) {
-    group <- summary$cluster_ids[i]
-    group_indices <- which(class2D == group)
-    median_cube[group_indices, ] <- matrix(
-      summary$median_spectra[i, ],
-      nrow = length(group_indices),
-      ncol = ncol(IFU2D),
-      byrow = TRUE
-    )
-  }
-
-  reshaped_cube <- array(median_cube, dim = dim(cubedat$imDat))
+  reshaped_cube <- recon$model_cube$imDat
 
   post_process_table <- as.data.frame(IFU2D)
   post_process_table$class <- class2D

R/reconstruct_cluster_cube.R

@@ -0,0 +1,152 @@
+#' Reconstruct a Model Cube from Cluster Representative Spectra
+#'
+#' This function builds a model cube from a segmentation result by assigning one
+#' representative spectrum to every pixel in each cluster. It is useful for
+#' visualization, denoising, and downstream spectral analysis where a
+#' cluster-level spectral template is desired.
+#'
+#' @param cluster_result A list returned by a segmentation function.
+#' @param template Representative spectrum to assign to each cluster.
+#'   \code{"median"} is robust, \code{"mean"} is the arithmetic cluster mean,
+#'   and \code{"weighted_mean"} uses inverse-variance weighting.
+#' @param var_cube Optional variance cube. Required when
+#'   \code{template = "weighted_mean"}.
+#' @param variance_inflation Multiplicative factor applied to propagated
+#'   variances when \code{var_cube} is supplied.
+#' @param preserve_mask Logical; if \code{TRUE}, channels that were non-finite in
+#'   the original cube remain masked in the reconstructed cube.
+#' @param fill_mode Either \code{"na"} or \code{"zero"} for unassigned pixels and,
+#'   when \code{preserve_mask = TRUE}, for masked spectral channels.
+#' @param return_residual Logical; if \code{TRUE}, also return the residual cube.
+#'
+#' @return A list with the reconstructed cube, optional residual cube, the
+#'   template spectra, a cluster summary object, and a global flux comparison
+#'   between the original and reconstructed cubes.
+#' @export
+reconstruct_cluster_cube <- function(cluster_result,
+                                     template = c("median", "mean", "weighted_mean"),
+                                     var_cube = NULL,
+                                     variance_inflation = 1,
+                                     preserve_mask = TRUE,
+                                     fill_mode = c("na", "zero"),
+                                     return_residual = TRUE) {
+  template <- match.arg(template)
+  fill_mode <- match.arg(fill_mode)
+
+  if (identical(template, "weighted_mean") && is.null(var_cube)) {
+    stop("`var_cube` must be supplied when `template = \"weighted_mean\"`.")
+  }
+
+  cubedat <- .as_cubedat(cluster_result$original_cube)
+  flux_mat <- cube_to_matrix(cubedat)
+  cluster_vec <- as.vector(cluster_result$cluster_map)
+
+  if (length(cluster_vec) != nrow(flux_mat)) {
+    stop("Cluster map dimensions do not match the input cube dimensions.")
+  }
+
+  summary <- summarize_cluster_spectra(
+    cluster_result = cluster_result,
+    var_cube = var_cube,
+    variance_inflation = variance_inflation
+  )
+
+  template_spectra <- switch(
+    template,
+    median = summary$median_spectra,
+    mean = summary$mean_spectra,
+    weighted_mean = summary$weighted_mean_spectra
+  )
+
+  model_mat <- matrix(NA_real_, nrow = nrow(flux_mat), ncol = ncol(flux_mat))
+
+  for (i in seq_along(summary$cluster_ids)) {
+    idx <- which(cluster_vec == summary$cluster_ids[i])
+    if (!length(idx)) {
+      next
+    }
+
+    model_mat[idx, ] <- matrix(
+      template_spectra[i, ],
+      nrow = length(idx),
+      ncol = ncol(flux_mat),
+      byrow = TRUE
+    )
+  }
+
+  if (isTRUE(preserve_mask)) {
+    orig_finite <- is.finite(flux_mat)
+    if (fill_mode == "na") {
+      model_mat[!orig_finite & !is.na(cluster_vec)] <- NA_real_
+    } else {
+      model_mat[!orig_finite & !is.na(cluster_vec)] <- 0
+    }
+  }
+
+  unassigned <- is.na(cluster_vec)
+  if (fill_mode == "zero") {
+    model_mat[unassigned, ] <- 0
+  }
+
+  original_sum_spectrum <- colSums(flux_mat[!unassigned, , drop = FALSE], na.rm = TRUE)
+  model_sum_spectrum <- colSums(model_mat[!unassigned, , drop = FALSE], na.rm = TRUE)
+  flux_difference <- model_sum_spectrum - original_sum_spectrum
+  flux_check <- list(
+    original_sum_spectrum = original_sum_spectrum,
+    model_sum_spectrum = model_sum_spectrum,
+    difference = flux_difference,
+    max_abs_difference = max(abs(flux_difference), na.rm = TRUE)
+  )
+
+  model_cube <- cubedat
+  model_cube$imDat <- array(model_mat, dim = dim(cubedat$imDat))
+
+  template_df <- data.frame(
+    cluster = summary$cluster_ids,
+    template_spectra,
+    check.names = FALSE
+  )
+
+  out <- list(
+    model_cube = model_cube,
+    template = template,
+    template_spectra = template_df,
+    cluster_summary = summary,
+    flux_check = flux_check
+  )
+
+  if (isTRUE(return_residual)) {
+    residual_cube <- cubedat
+    residual_cube$imDat <- array(flux_mat - model_mat, dim = dim(cubedat$imDat))
+    out$residual_cube <- residual_cube
+  }
+
+  out
+}
+
+#' Reconstruct a Flux-preserving Model Cube from Cluster Spectra
+#'
+#' This is a convenience wrapper around \code{\link{reconstruct_cluster_cube}}
+#' that uses the arithmetic cluster mean while preserving the original spectral
+#' mask. In this mode, the reconstructed cube preserves the summed flux spectrum
+#' of the segmented cube on the observed support, which makes it a sensible
+#' default for later spectral fitting.
+#'
+#' @param cluster_result A list returned by a segmentation function.
+#' @param fill_mode Either \code{"na"} or \code{"zero"} for unassigned pixels
+#'   and masked channels.
+#' @param return_residual Logical; if \code{TRUE}, also return the residual cube.
+#'
+#' @return The same structure returned by \code{\link{reconstruct_cluster_cube}}.
+#' @export
+reconstruct_flux_preserving_cube <- function(cluster_result,
+                                             fill_mode = c("na", "zero"),
+                                             return_residual = TRUE) {
+  reconstruct_cluster_cube(
+    cluster_result = cluster_result,
+    template = "mean",
+    preserve_mask = TRUE,
+    fill_mode = fill_mode,
+    return_residual = return_residual
+  )
+}

R/summarize_cluster_spectra.R

@@ -2,7 +2,8 @@
 #'
 #' This is the recommended post-segmentation summary layer for Capivara. It keeps
 #' the segmentation step separate from downstream spectral products and can return
-#' median spectra, summed spectra, and inverse-variance-weighted mean spectra.
+#' median spectra, arithmetic mean spectra, summed spectra, and
+#' inverse-variance-weighted mean spectra.
 #'
 #' @param cluster_result A list returned by a segmentation function.
 #' @param var_cube Optional variance cube matching the dimensions of the original
@@ -11,8 +12,9 @@
 #'   variances. Use this to account for covariance if needed.
 #'
 #' @return A list with wavelength coordinates, cluster ids, cluster sizes,
-#'   median spectra, summed spectra, and, when \code{var_cube} is supplied,
-#'   propagated sum variances and inverse-variance-weighted means.
+#'   per-wavelength finite counts, median spectra, mean spectra, summed spectra,
+#'   and, when \code{var_cube} is supplied, propagated sum variances and
+#'   inverse-variance-weighted means.
 #' @export
 summarize_cluster_spectra <- function(cluster_result,
                                       var_cube = NULL,
@@ -39,7 +41,14 @@ summarize_cluster_spectra <- function(cluster_result,
     ncol = ncol(flux_mat),
     dimnames = list(as.character(cluster_ids), wave_names)
   )
+  mean_spectra <- median_spectra
   sum_spectra <- median_spectra
+  finite_counts <- matrix(
+    0L,
+    nrow = length(cluster_ids),
+    ncol = ncol(flux_mat),
+    dimnames = list(as.character(cluster_ids), wave_names)
+  )
   n_spaxels <- integer(length(cluster_ids))
 
   has_var <- !is.null(var_cube)
@@ -63,6 +72,9 @@ summarize_cluster_spectra <- function(cluster_result,
     n_spaxels[i] <- length(idx)
     median_spectra[i, ] <- apply(X, 2, stats::median, na.rm = TRUE)
     sum_spectra[i, ] <- colSums(X, na.rm = TRUE)
+    finite_counts[i, ] <- colSums(is.finite(X))
+    mean_spectra[i, ] <- sum_spectra[i, ] / finite_counts[i, ]
+    mean_spectra[i, finite_counts[i, ] == 0] <- NA_real_
 
     if (has_var) {
       V <- var_mat[idx, , drop = FALSE]
@@ -88,7 +100,9 @@ summarize_cluster_spectra <- function(cluster_result,
     wavelength = wavelengths,
     cluster_ids = cluster_ids,
     n_spaxels = stats::setNames(n_spaxels, cluster_ids),
+    finite_counts = finite_counts,
     median_spectra = median_spectra,
+    mean_spectra = mean_spectra,
     sum_spectra = sum_spectra
   )
 

README.md

@@ -128,6 +128,18 @@ spec_summary <- summarize_cluster_spectra(res, var_cube = var_cube$imDat)
 inspection. For flux and SNR calculations, `sum_spectra` or the
 inverse-variance-weighted summary are usually better choices.
 
+### Reconstructed Cubes
+
+Use `reconstruct_cluster_cube()` to build a representative cube from
+cluster templates, or `reconstruct_flux_preserving_cube()` when you want
+a model cube that preserves the summed flux spectrum of the segmented
+data for later spectral fitting.
+
+``` r
+rep_cube <- reconstruct_cluster_cube(res_star, template = "median")
+fit_cube <- reconstruct_flux_preserving_cube(res_star)
+```
+
 ## Release Notes
 
 See [NEWS.md](NEWS.md) for the `0.2.0` release summary.
@@ -167,7 +179,6 @@ research. A BibTeX entry for the paper is:
     **539**(4), 3166–3179. <https://doi.org/10.1093/mnras/staf688>
 3.  **Torch in R**: Paszke, Adam, et al. “PyTorch: An Imperative Style,
     High-Performance Deep Learning Library.” Advances in Neural
-    Information Processing Systems 2019.
-------------------------------------------------------------------------
-
-For more information, check the [Capivara GitHub webpage](https://rafaelsdesouza.github.io/capivara/).
+    Information Processing Systems 2019. For more information, check the
+    [Capivara GitHub
+    webpage](https://rafaelsdesouza.github.io/capivara/).

README.qmd

@@ -121,6 +121,18 @@ spec_summary <- summarize_cluster_spectra(res, var_cube = var_cube$imDat)
 For flux and SNR calculations, `sum_spectra` or the inverse-variance-weighted
 summary are usually better choices.
 
+### Reconstructed Cubes
+
+Use `reconstruct_cluster_cube()` to build a representative cube from cluster
+templates, or `reconstruct_flux_preserving_cube()` when you want a model cube
+that preserves the summed flux spectrum of the segmented data for later
+spectral fitting.
+
+```R
+rep_cube <- reconstruct_cluster_cube(res_star, template = "median")
+fit_cube <- reconstruct_flux_preserving_cube(res_star)
+```
+
 ## Release Notes
 
 See [NEWS.md](NEWS.md) for the `0.2.0` release summary.

index.qmd

@@ -98,6 +98,13 @@ res <- segment(cube, Ncomp = choice$Ncomp)
 spec_summary <- summarize_cluster_spectra(res, var_cube = var_cube$imDat)
 ```
 
+### Reconstructed Cubes
+
+```R
+rep_cube <- reconstruct_cluster_cube(res_star, template = "median")
+fit_cube <- reconstruct_flux_preserving_cube(res_star)
+```
+
 ## Attribution
 
 Please cite de Souza et al. if you find this code useful in your research.