extract rois supports percentage or proportion inputs; prompt_input allows for explicit empty return types; brain_mask

is prompted for during roi extraction

Author: michaelhallquist

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BrainGnomes
 Title: An R package for preprocessing and analyzing fMRI data
-Version: 0.7-1
-Date: 2025-09-24
+Version: 0.7-2
+Date: 2025-09-28
 Authors@R: 
     person("Michael", "Hallquist", email="michael.hallquist@gmail.com", role = c("aut", "cre"))
 Description: This R package handles preprocessing of fMRI data on HPC clusters.

NEWS.md

@@ -1,5 +1,15 @@
+# BrainGnomes 0.7-2
+
+* Allow `min_vox_per_roi` to be specified as a percentage or proportion of atlas voxels during ROI extraction
+* ROI extraction allows for an optional brain mask that is applied to the atlas and time series data
+* Add ability to set an optional field to NA (empty) in `prompt_input` when a default is provided
+* Support notch filtering of motion parameters prior to calculation of framewise_displacement
+* Support different head sizes for calculation of framewise displacement
+
 # BrainGnomes 0.7-1
 
+Released 2025-09-24
+
 * UI/UX improvements to ask query enable/disable for postprocessing and ROI extraction during edit_project
 * Do not ask about ROI extraction details if no postprocessing streams are defined
 * Corrected validation of band-pass cutoffs and clarified temporal filtering prompts/documentation

R/extract_rois.R

@@ -9,8 +9,7 @@
 #' the input timeseries.
 #'
 #' @param bold_file Path to a 4D NIfTI file containing postprocessed BOLD data.
-#' @param atlas_files Character vector of atlas NIfTI files with integer
-#'   ROI labels.
+#' @param atlas_files Character vector of atlas NIfTI files with integer ROI labels.
 #' @param out_dir Directory where output files should be written.
 #' @param log_file If not `NULL`, the log file to which details should be written.
 #' @param cor_method Correlation method(s) to use when computing functional
@@ -19,12 +18,14 @@
 #'   computation. Multiple methods may be supplied.
 #' @param roi_reduce Method used to summarize voxel time series within each
 #'   ROI. Options are "mean" (default), "median", "pca", or "huber".
-#' @param brain_mask Optional brain mask NIfTI file. If \code{NULL}, a mask
-#'   is generated by excluding voxels with zero variance across time.
-#' @param min_vox_per_roi The minimum number of voxels required for an ROI to
-#'   be extracted and entered into correlations. If the ROI is smaller than this,
-#'   it show up as NA in the outputs, keeping the dimensionality of the connectivity
-#'   matrix consistent across inputs. Default: `5`
+#' @param mask_file Optional path to a mask NIfTI file. Voxels outside of this mask
+#'   are excluded from ROI extraction and connectivity calculation. Note that
+#'   constant and zero voxels are always automatically removed by extract_rois.
+#' @param min_vox_per_roi Minimum ROI size requirement. Supply a positive integer
+#'   to require at least that many ROI voxels survive masking and are non-zero, or provide 
+#'   a proportion (e.g., `0.8`) or percentage string (e.g., `80%`) to require that fraction of
+#'   the ROI voxels to remain. ROIs failing this check are set to `NA`, preserving
+#'   consistent ROI matrix size. Default: `5`.
 #' @param save_ts If `TRUE`, save the ROI time series (aggregated using `roi_reduce` method)
 #'   to `_timeseries.tsv`. files. Useful for running external analyses on the ROIs. Default: `TRUE`.
 #' @param rtoz If `TRUE`, using Fisher's z (aka atanh) transformation on correlations to make them
@@ -35,24 +36,28 @@
 #' @return A named list. Each element corresponds to an atlas and contains
 #'   paths to the written timeseries (\code{timeseries}) and correlation
 #'   matrix (\code{correlation}, or \code{NULL} if not computed).
-#' @importFrom checkmate assert_file_exists assert_character assert_directory_exists assert_integerish assert_flag
+#' @importFrom checkmate assert_file_exists assert_character assert_directory_exists assert_flag
 #' @export
 extract_rois <- function(bold_file, atlas_files, out_dir, log_file = NULL,
                          cor_method = c("pearson", "spearman", "kendall", "cor.shrink"),
                          roi_reduce = c("mean", "median", "pca", "huber"),
-                         brain_mask = NULL, min_vox_per_roi = 5, save_ts = TRUE, rtoz = FALSE,
+                         mask_file = NULL, min_vox_per_roi = 5, save_ts = TRUE, rtoz = FALSE,
                          overwrite = FALSE) {
   checkmate::assert_file_exists(bold_file)
   checkmate::assert_character(atlas_files, any.missing = FALSE, min.len = 1)
   checkmate::assert_directory_exists(out_dir, access = "w")
   cor_method <- match.arg(cor_method, several.ok = TRUE)
   roi_reduce <- match.arg(roi_reduce)
-  checkmate::assert_integerish(min_vox_per_roi, len = 1L, lower = 1L)
+  checkmate::assert_string(mask_file, null.ok = TRUE, na.ok = TRUE)
+  if (isTRUE(is.na(mask_file[1L]))) mask_file <- NULL
+
+  min_vox_spec <- parse_min_vox_per_roi(min_vox_per_roi)
   checkmate::assert_flag(save_ts)
   checkmate::assert_flag(rtoz)
+  checkmate::assert_flag(overwrite)
 
   lg <- lgr::get_logger_glue("extract_rois")
-  lg$config(NULL)
+  lg$config(NULL) # reset logger object to clear any appender files
   if (!is.null(log_file)) lg$add_appender(lgr::AppenderFile$new(log_file), name = "extract_logger")
 
   # Read 4D NIfTI
@@ -64,16 +69,41 @@ extract_rois <- function(bold_file, atlas_files, out_dir, log_file = NULL,
   n_time <- dim_img[4]
   mat <- matrix(bold_img, prod(dim_img[1:3]), n_time)
 
-  # Determine brain mask
-  if (is.null(brain_mask)) {
-    # if no mask provided, remove zero voxels and constant. Check all(v==0) because it's faster than var()
-    brain_mask_vec <- apply(mat, 1, function(v) !all(v == 0) && stats::var(v) > 0)
-  } else if (checkmate::test_file_exists(brain_mask)) {
-    brain_mask_vec <- as.vector(RNifti::readNifti(brain_mask)) > 0
-  } else {
-    to_log(lg, "fatal", "brain_mask must be a valid NIfTI file or NULL")
+  # Start with BOLD-derived mask that drops constant voxels or those with NAs.
+  # Use the !all(zero) to screen out 0 voxels because is it faster than computing the variance
+  mask_vec <- apply(mat, 1L, function(v) {
+    var_ts <- stats::var(v)
+    !anyNA(v) && # no NAs
+      !all(abs(v) < 2 * .Machine$double.eps) && # not all zero
+      !is.na(var_ts) && # variance is defined
+      var_ts > 2 * .Machine$double.eps # variance is positive
+  })
+
+  # Handle user-specified mask, if provided
+  if (!is.null(mask_file)) {
+    if (!checkmate::test_file_exists(mask_file)) {
+      to_log(lg, "fatal", "mask_file must be a valid NIfTI file or NULL")
+    }
+
+    mask_img <- RNifti::readNifti(mask_file)
+    mask_dims <- dim(mask_img)
+    if (length(mask_dims) > 3L) mask_dims <- mask_dims[1:3]
+    if (!identical(mask_dims, dim_img[1:3])) {
+      to_log(lg, "fatal", "Mask dimensions {paste(mask_dims, collapse = 'x')} do not match BOLD grid {paste(dim_img[1:3], collapse = 'x')}")
+    }
+
+    provided_mask_vec <- as.vector(mask_img > 0)
+    if (length(provided_mask_vec) != length(mask_vec)) {
+      to_log(lg, "fatal", "Mask voxel count ({length(provided_mask_vec)}) does not match BOLD grid ({length(mask_vec)})")
+    }
+
+    # intersect mask file with internal automask (for 0/constant voxels)
+    mask_vec <- mask_vec & provided_mask_vec
   }
 
+  mask_vec[is.na(mask_vec)] <- FALSE
+  mask_vec <- as.logical(mask_vec)
+
   compute_correlation <- !is.null(cor_method) && length(cor_method) > 0L
   bids_info <- as.list(extract_bids_info(bold_file))
   sub_id <- bids_info$subject
@@ -100,35 +130,32 @@ extract_rois <- function(bold_file, atlas_files, out_dir, log_file = NULL,
     }
 
     atlas_vec <- as.vector(atlas_img)
+
     if (!checkmate::test_integerish(atlas_vec, tol = 1e-6)) stop("Atlas ", atlas, " contains non-integer labels (outside tolerance).")
-    roi_vals <- sort(unique(atlas_vec[atlas_vec > 0]))
+    roi_vals <- sort(unique(atlas_vec[atlas_vec > 0 & mask_vec]))
 
     # ROI reduction
     ts_mat <- sapply(roi_vals, function(lbl) {
-      vox <- which(atlas_vec == lbl & brain_mask_vec)
-      if (length(vox) < min_vox_per_roi) {
-        to_log(lg, "info", "Fewer than {min_vox_per_roi} voxels in ROI {lbl}. Dropping")
+      roi_voxels <- sum(atlas_vec == lbl)
+      required_vox <- compute_min_vox_required(min_vox_spec, roi_voxels)
+      roi_idx <- which((atlas_vec == lbl) & mask_vec) # get voxel positions for this ROI in the brain mask
+      if (length(roi_idx) < required_vox) {
+        req_txt <- as.character(format_min_vox_requirement(min_vox_spec, roi_voxels))
+        to_log(lg, "info", "ROI {lbl} has {length(roi_idx)} usable voxels but requires {req_txt}. Dropping")
         rep(NA_real_, n_time)
       } else {
-        vals <- mat[vox, , drop = FALSE]
-        bad <- apply(vals, 1, function(ts) any(is.na(ts)) || all(ts == 0) || stats::var(ts) == 0)
-        if (sum(!bad) < min_vox_per_roi) {
-          to_log(lg, "info", "Fewer than {min_vox_per_roi} good time series in ROI {lbl}. Dropping")
-          rep(NA_real_, n_time)
-        } else {
-          roivox <- t(vals[!bad, , drop = FALSE]) # time x voxels
-          if (roi_reduce == "pca") {
-            pc <- stats::prcomp(roivox, scale. = TRUE)$x[, 1]
-            mn <- rowMeans(roivox) # because direction of eigenvector is arbitrary, ensure it scales positively with the mean
-            if (stats::cor(pc, mn) < 0) pc <- -pc
-            pc
-          } else if (roi_reduce == "median") {
-            apply(roivox, 1, median)
-          } else if (roi_reduce == "huber") {
-            apply(roivox, 1, function(x) huber(x)$mu)
-          } else { # mean
-            rowMeans(roivox)
-          }
+        roi_vox <- t(mat[roi_idx, , drop = FALSE]) # time x voxels
+        if (roi_reduce == "pca") {
+          pc <- stats::prcomp(roi_vox, scale. = TRUE)$x[, 1]
+          mn <- rowMeans(roi_vox) # because direction of eigenvector is arbitrary, ensure it scales positively with the mean
+          if (stats::cor(pc, mn) < 0) pc <- -pc
+          pc
+        } else if (roi_reduce == "median") {
+          apply(roi_vox, 1, median)
+        } else if (roi_reduce == "huber") {
+          apply(roi_vox, 1, function(x) huber(x)$mu)
+        } else { # mean
+          rowMeans(roi_vox)
         }
       }
     })
@@ -275,4 +302,92 @@ huber <- function(y, k = 1.5, tol = 1.0e-6) {
     mu <- mu1
   }
   list(mu = mu, s = s)
-}
+}
+
+
+
+#' Compute ROI voxel number requirement given a parsed specification
+#' @param roi_voxels the number of voxels in an ROI to be potentially extracted
+#' @keywords internal
+#' @noRd
+compute_min_vox_required <- function(spec, roi_voxels) {
+  checkmate::assert_list(spec, any.missing = FALSE)
+  checkmate::assert_number(roi_voxels, lower = 0, finite = TRUE)
+
+  if (spec$type == "count") {
+    req <- spec$value
+  } else if (spec$type == "fraction") {
+    req <- ceiling(spec$value * roi_voxels)
+  } else {
+    stop("Unknown specification type for min_vox_per_roi", call. = FALSE)
+  }
+
+  req <- as.integer(req)
+  if (is.na(req) || req < 1L) req <- 1L
+  return(req)
+}
+
+#' Internal utilities shared across ROI extraction functions.
+#' @keywords internal
+#' @noRd
+parse_min_vox_per_roi <- function(spec) {
+  if (length(spec) != 1L) {
+    stop("min_vox_per_roi must be a single value", call. = FALSE)
+  }
+
+  if (is.na(spec)) {
+    stop("min_vox_per_roi cannot be NA", call. = FALSE)
+  }
+
+  # Numeric input: allow integer counts or proportions in (0, 1]
+  if (is.numeric(spec)) {
+    if (spec >= 1 && checkmate::test_integerish(spec, tol = 1e-8)) {
+      return(list(type = "count", value = as.integer(round(spec))))
+    } else if (spec > 0 && spec <= 1) {
+      return(list(type = "fraction", value = as.numeric(spec)))
+    }
+  }
+
+  if (is.character(spec)) {
+    trimmed <- trimws(spec)
+    if (trimmed == "") stop("min_vox_per_roi cannot be an empty string", call. = FALSE)
+
+    if (grepl("%$", trimmed)) {
+      pct <- suppressWarnings(as.numeric(sub("%$", "", trimmed)))
+      if (is.na(pct)) stop("Percentage min_vox_per_roi must contain a valid number before '%'", call. = FALSE)
+      frac <- pct / 100
+      if (frac <= 0 || frac > 1) stop("Percentage min_vox_per_roi must be between 0% and 100%", call. = FALSE)
+      return(list(type = "fraction", value = frac))
+    }
+
+    # Handle numeric strings recursively
+    num <- suppressWarnings(as.numeric(trimmed))
+    if (!is.na(num)) {
+      return(parse_min_vox_per_roi(num))
+    }
+  }
+
+  stop("min_vox_per_roi must be a positive integer, a proportion in (0, 1], or a percentage string such as '80%'", call. = FALSE)
+}
+
+
+
+#' Format min_vox_per_roi requirement for display or storage
+#' @keywords internal
+#' @noRd
+format_min_vox_requirement <- function(spec, roi_voxels = NULL, digits = 1) {
+  checkmate::assert_list(spec, any.missing = FALSE)
+  checkmate::assert_number(digits, lower = 0, finite = TRUE)
+
+  if (spec$type == "count") {
+    return(glue::glue("{spec$value} voxels"))
+  }
+
+  pct <- round(spec$value * 100, digits)
+  if (is.null(roi_voxels)) {
+    return(glue::glue("{pct}% of ROI voxels"))
+  }
+
+  req <- compute_min_vox_required(spec, roi_voxels)
+  return(glue::glue("{pct}% of {roi_voxels} voxels (>= {req})"))
+}

R/setup_extract.R

@@ -1,7 +1,7 @@
 #' @keywords internal
 manage_extract_streams <- function(scfg, allow_empty = FALSE) {
   extract_field_list <- function() {
-    c("input_streams", "atlases", "roi_reduce", "correlation/method", "rtoz")
+    c("input_streams", "atlases", "mask_file", "roi_reduce", "correlation/method", "rtoz")
   }
 
   show_val <- function(val) {
@@ -191,7 +191,7 @@ setup_extract_stream <- function(scfg, fields = NULL, stream_name = NULL) {
 
   if (is.null(fields)) {
     fields <- c(
-      "extract_rois/input_streams", "extract_rois/atlases", "extract_rois/roi_reduce", "extract_rois/save_ts",
+      "extract_rois/input_streams", "extract_rois/atlases", "extract_rois/mask_file", "extract_rois/roi_reduce", "extract_rois/save_ts",
       "extract_rois/correlation/method", "extract_rois/rtoz", "extract_rois/min_vox_per_roi"
     )
   }
@@ -206,11 +206,25 @@ setup_extract_stream <- function(scfg, fields = NULL, stream_name = NULL) {
   }
 
   if ("extract_rois/atlases" %in% fields) {
-    atlas <- prompt_input("Atlas NIfTI file(s) (comma separated):", type = "character")
-    atlas <- trimws(strsplit(atlas, ",")[[1]])
+    atlas <- prompt_input("Atlas NIfTI file(s) (separate using commas for multiple):", type = "character", split="\\s*,\\s*")
     excfg$atlases <- atlas
   }
 
+  if ("extract_rois/mask_file" %in% fields) {
+    mask_file <- prompt_input(
+      instruct = glue("\n
+        ROI extraction always removes constant and zero voxels from its calculations.
+        In addition, you may provide a mask file that only retains voxels in the mask
+        for ROI extraction and connectivity calculation. This mask is applied in addition
+        to masking done in postprocessing, if requested, as well as the automatic removal
+        of constant/zero voxels.
+      "),
+      prompt = "Mask file for ROI extraction",
+      type = "character", required = FALSE,
+      default = excfg$mask_file
+    )
+  }
+
   if ("extract_rois/roi_reduce" %in% fields) {
     reduce <- select_list_safe(c("mean", "median", "pca", "huber"),
                                multiple = FALSE, title = "ROI reduction method")
@@ -258,16 +272,52 @@ setup_extract_stream <- function(scfg, fields = NULL, stream_name = NULL) {
   }
 
   if ("extract_rois/min_vox_per_roi" %in% fields) {
-    excfg$min_vox_per_roi <- prompt_input(instruct = glue("
-      Sometimes smaller ROIs, especially near the edge of the brain, may not have enough voxels to yield
-      a reasonable average time series. In this case, it's often best to exclude these ROIs, rather than
-      include a noisy estimate. How many voxels must an ROI have to be extracted (and entered into correlation)?
-      Note that ROIs less than this number will be set to NA in the output, but will still be labeled correctly rather
-      than being dropped altogether."),
-
-    "Enter the minimum number of voxels for valid ROIs:", type = "integer", lower = 1L, default=5L)
+    default_value <- if (is.null(excfg$min_vox_per_roi)) "5" else as.character(excfg$min_vox_per_roi)
+    repeat {
+      response <- prompt_input(
+        instruct = glue("
+          Sometimes smaller ROIs, especially near the edge of the brain, may not have enough voxels to yield
+          a reasonable average time series. In this case, it's often best to exclude these ROIs, rather than
+          include a noisy estimate. ROIs failing this threshold will be set to NA in outputs, preserving labeling.
+          You can enter an absolute voxel count (e.g., 5) or a percentage/proportion (e.g., '80%' or 0.8) of atlas voxels."),
+        prompt = "Enter the minimum voxel requirement for valid ROIs:", type = "character", len = 1L,
+        default = default_value
+      )
+
+      response <- trimws(response)
+      parsed <- tryCatch(parse_min_vox_per_roi(response), error = function(e) {
+        cat(conditionMessage(e), "\n")
+        NULL
+      })
+
+      if (!is.null(parsed)) {
+        excfg$min_vox_per_roi <- canonicalize_min_vox_per_roi(parsed)
+        break
+      }
+
+      default_value <- response
+    }
   }
 
   scfg$extract_rois[[stream_name]] <- excfg
   return(scfg)
 }
+
+
+#' Convert a parsed specification to a canonical configuration value
+#' @keywords internal
+#' @noRd
+canonicalize_min_vox_per_roi <- function(spec, digits = 1) {
+  checkmate::assert_list(spec, any.missing = FALSE)
+
+  if (spec$type == "count") return(spec$value)
+
+  pct <- round(spec$value * 100, digits)
+  fmt <- if (abs(pct - round(pct)) < .Machine$double.eps^0.5) {
+    as.integer(round(pct))
+  } else {
+    pct
+  }
+
+  return(paste0(fmt, "%"))
+}