R/automate_particle_filespecs.R

In OpenSpecy: Analyze, Process, Identify, and Share Raman and (FT)IR Spectra

#' @rdname automate_particle_analysis
#' @export
automate_particle_analysis.FileSpecs <- function(
    x, library, output_dir = NULL, images = NULL, bottom_left = NULL,
    top_right = NULL, origins = NULL, material_col = "material_class",
    library_id_col = "sample_name",
    particle_id_strategy = c("collapse", "partial_collapse",
                             "nonspatial_collapse", "all_cell_id", "raw"),
    spectral_smooth = FALSE, sigma1 = c(1, 1, 1),
    sigma2 = c(3, 3), close = FALSE,
    close_kernel = c(4, 4), sn_threshold_min = 0.04,
    sn_threshold_max = Inf, cor_threshold = 0.7, area_threshold = 1,
    label_unknown = FALSE, remove_materials = NULL, remove_unknown = FALSE,
    pixel_length = 25, metric = "sig_times_noise", abs = FALSE,
    collapse_function = stats::median,
    outputs = c("details", "summary"),
    process_args = list(), specs_steps = c("pca", "kmeans"),
    specs_centers = NULL, file_processing = c("stream", "memory"), ...) {
  .reject_removed_particle_args(list(...))
  file_processing <- match.arg(file_processing)
  .validate_particle_sn_thresholds(sn_threshold_min, sn_threshold_max)
  .filespec_validate_object(x)
  .filespec_validate_source(x, strong = FALSE)
  if (identical(file_processing, "memory")) {
    views <- split_spec(x, by = "region")
    if (!length(views)) views <- list(source = x)
    if (is.null(names(views)) || any(!nzchar(names(views)))) {
      names(views) <- paste0("region_", seq_along(views))
    }
    output_names <- .particle_filespec_output_names(x, names(views))
    counts <- vapply(views, .filespec_n_spectra, integer(1))
    materialize_started <- proc.time()[["elapsed"]]
    .particle_progress(
      "file-backed source", "read into memory",
      sprintf(
        "%s regions; %s spectra; %s bands",
        format(length(views), big.mark = ","),
        format(sum(counts), big.mark = ","),
        format(length(.filespec_axis(x)), big.mark = ",")
      )
    )
    maps <- lapply(seq_along(views), function(i) {
      sample_name <- names(views)[[i]]
      .particle_progress(
        sample_name, "memory materialization",
        sprintf("region %s/%s; %s spectra", i, length(views),
                format(counts[[i]], big.mark = ","))
      )
      map <- decompress_spec(views[[i]], index = seq_len(counts[[i]]))
      attr(map, "particle_output_name") <- output_names[[i]]
      .particle_progress(
        sample_name, "memory materialization complete",
        sprintf("elapsed %.1f s",
                proc.time()[["elapsed"]] - materialize_started)
      )
      map
    })
    names(maps) <- names(views)
    return(automate_particle_analysis.default(
      x = maps, library = library, output_dir = output_dir, images = images,
      bottom_left = bottom_left, top_right = top_right, origins = origins,
      material_col = material_col, library_id_col = library_id_col,
      particle_id_strategy = particle_id_strategy,
      spectral_smooth = spectral_smooth, sigma1 = sigma1, sigma2 = sigma2,
      close = close, close_kernel = close_kernel,
      sn_threshold_min = sn_threshold_min, sn_threshold_max = sn_threshold_max,
      cor_threshold = cor_threshold, area_threshold = area_threshold,
      label_unknown = label_unknown, remove_materials = remove_materials,
      remove_unknown = remove_unknown, pixel_length = pixel_length,
      metric = metric, abs = abs, collapse_function = collapse_function,
      outputs = outputs, process_args = process_args,
      specs_steps = specs_steps, specs_centers = specs_centers,
      file_processing = "memory"
    ))
  }
  strategy <- .normalize_particle_strategy(particle_id_strategy)
  if (!strategy %in% c("collapse", "all_cell_id")) {
    stop("FileSpecs particle analysis currently supports ",
         "'particle_id_strategy = \"collapse\"' or \"all_cell_id\"",
         call. = FALSE)
  }
  if (!identical(match.fun(collapse_function), base::mean)) {
    stop("FileSpecs particle analysis currently requires ",
         "'collapse_function = mean'", call. = FALSE)
  }
  if (isTRUE(spectral_smooth) &&
      (!is.numeric(sigma1) || length(sigma1) != 3L || anyNA(sigma1) ||
       any(sigma1 < 0))) {
    stop("'sigma1' must be a nonnegative numeric vector of length 3 when ",
         "spectral_smooth = TRUE", call. = FALSE)
  }
  if (identical(metric, "entropy")) {
    stop("FileSpecs entropy S/N requires explicit global breaks and is not ",
         "part of the initial particle pipeline", call. = FALSE)
  }

  outputs <- .normalize_particle_outputs(outputs)
  if (!is.null(output_dir)) {
    dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  }
  views <- split_spec(x, by = "region")
  if (!length(views)) views <- list(source = x)
  if (is.null(names(views)) || any(!nzchar(names(views)))) {
    names(views) <- paste0("region_", seq_along(views))
  }
  output_names <- .particle_filespec_output_names(x, names(views))

  sample_results <- lapply(seq_along(views), function(i) {
    .particle_progress(names(views)[[i]], "region", sprintf("%d of %d", i,
                                                            length(views)))
    .automate_particle_filespec_region(
      x = views[[i]], library = library, sample_name = names(views)[[i]],
      output_name = output_names[[i]],
      output_dir = output_dir, image = .indexed_argument(images, i),
      bottom_left = .indexed_argument(bottom_left, i),
      top_right = .indexed_argument(top_right, i),
      origin = .particle_origin(origins, i), material_col = material_col,
      library_id_col = library_id_col,
      spectral_smooth = spectral_smooth, sigma1 = sigma1,
      sigma2 = sigma2, close = close,
      close_kernel = close_kernel,
      particle_id_strategy = strategy,
      sn_threshold_min = sn_threshold_min,
      sn_threshold_max = sn_threshold_max, cor_threshold = cor_threshold,
      area_threshold = area_threshold, label_unknown = label_unknown,
      remove_materials = remove_materials, remove_unknown = remove_unknown,
      pixel_length = pixel_length, metric = metric, abs = abs,
      outputs = outputs, process_args = process_args
    )
  })
  names(sample_results) <- names(views)

  details_all <- .bind_particle_tables(
    lapply(sample_results, .sample_particle_item, "particle_details_csv")
  )
  summary_all <- .bind_particle_tables(
    lapply(sample_results, .sample_particle_item, "particle_summary_csv")
  )
  if (!is.null(output_dir)) {
    .write_particle_all_outputs(output_dir, details_all, summary_all, outputs)
  }
  structure(
    list(samples = sample_results, particle_details_all_csv = details_all,
         particle_summary_all_csv = summary_all),
    class = c("OpenSpecyParticleAnalysis", "list")
  )
}

.particle_filespec_output_names <- function(x, region_names) {
  paths <- as.character(x$source$members$path)
  extensions <- tolower(tools::file_ext(paths))
  preferred <- if (identical(x$source$backend, "h5")) {
    which(extensions %in% c("h5", "hdf5"))
  } else {
    which(extensions %in% c("dat", "img"))
  }
  if (!length(preferred)) preferred <- seq_along(paths)
  stem <- if (length(preferred)) {
    tools::file_path_sans_ext(basename(paths[[preferred[[1L]]]]))
  } else {
    "source"
  }
  source_region_count <- if (identical(x$source$backend, "h5")) {
    length(x$source$layout$regions)
  } else {
    1L
  }
  names <- if (length(region_names) > 1L || source_region_count > 1L) {
    paste(stem, region_names, sep = "_")
  } else {
    rep(stem, length(region_names))
  }
  make.unique(names)
}

.automate_particle_filespec_region <- function(
    x, library, sample_name, output_name, output_dir, image, bottom_left, top_right,
    origin, material_col, library_id_col, spectral_smooth, sigma1, sigma2,
    close, close_kernel, particle_id_strategy, sn_threshold_min,
    sn_threshold_max, cor_threshold, area_threshold, label_unknown,
    remove_materials, remove_unknown, pixel_length, metric, abs, outputs,
    process_args,
    chunk_size = getOption("OpenSpecy.filespec.chunk_size", 8192L)) {
  time_start <- Sys.time()
  .particle_progress(sample_name, "index")
  index <- data.table::copy(.filespec_index(x))
  if (!nrow(index)) {
    stop("the FileSpecs region view contains no spectra", call. = FALSE)
  }
  if (is.null(image) && "particle_image" %in% outputs) {
    visual <- .filespec_materialize_visual(x)
    if (!is.null(visual$image)) {
      image <- visual$image
      bottom_left <- visual$bottom_left
      top_right <- visual$top_right
    }
  }
  source_axis <- .filespec_axis(x)
  bands <- which(
    (source_axis >= 750 & source_axis <= 2200) |
      (source_axis >= 2420 & source_axis <= 4000)
  )
  if (!length(bands)) {
    stop("the FileSpecs axis does not overlap the particle S/N ranges",
         call. = FALSE)
  }
  cache_key <- digest::digest(list(
    schema = "filespec-particle-collapse-4", source = x$source$id,
    view = x$view, metric = metric, abs = abs,
    strategy = particle_id_strategy,
    spectral_smooth = isTRUE(spectral_smooth), sigma1 = sigma1,
    sigma2 = sigma2, close = close,
    close_kernel = close_kernel, sn_min = sn_threshold_min,
    sn_max = sn_threshold_max, area = area_threshold,
    image = .filespec_image_identity(image, bottom_left, top_right),
    library = if (identical(particle_id_strategy, "all_cell_id")) {
      digest::digest(library, algo = "sha256")
    } else NULL,
    process_args = if (identical(particle_id_strategy, "all_cell_id")) {
      process_args
    } else NULL,
    material_col = material_col, library_id_col = library_id_col
  ))
  cache_file <- .filespec_cache_path(x, "particle-collapse", cache_key)
  cached <- if (file.exists(cache_file)) {
    tryCatch(readRDS(cache_file), error = function(e) NULL)
  } else {
    NULL
  }

  if (is.null(cached)) {
    .particle_progress(
      sample_name, "streaming signal/noise",
      sprintf("%s spectra", format(nrow(index), big.mark = ","))
    )
    snr <- .filespec_particle_snr(
      x, index = index, bands = bands, metric = metric, abs = abs,
      spectral_smooth = spectral_smooth, sigma1 = sigma1,
      chunk_size = chunk_size, sample = sample_name
    )
    threshold <- snr > sn_threshold_min & snr < sn_threshold_max
    threshold[is.na(threshold)] <- FALSE
    display <- .filespec_particle_display(index, snr, threshold)
    display <- .attach_particle_image(display, list(image), list(bottom_left),
                                      list(top_right), 1L)
    threshold_state <- .particle_threshold_state(threshold)

    if (identical(threshold_state, "none")) {
      cached <- list(snr = snr, threshold = threshold,
                     feature_metadata = NULL, collapsed = NULL)
    } else {
      material <- max_cor_val <- NULL
      if (identical(particle_id_strategy, "all_cell_id")) {
        .particle_progress(sample_name, "streaming pixel identification")
        pixel_matches <- .filespec_particle_matches(
          x, eligible = threshold, library = library,
          process_args = process_args, material_col = material_col,
          library_id_col = library_id_col,
          spectral_smooth = spectral_smooth, sigma1 = sigma1,
          chunk_size = min(chunk_size, 1000L), sample = sample_name
        )
        match_index <- match(index$col_id, pixel_matches$object_id)
        material <- pixel_matches[[material_col]][match_index]
        max_cor_val <- pixel_matches$max_cor_val[match_index]
        material[!threshold] <- "background"
      }
      .particle_progress(sample_name, "streaming particle means")
      partition <- .filespec_collapse_connected_mean(
        x, eligible = threshold, material = material, snr = snr,
        max_cor_val = max_cor_val, area_threshold = area_threshold,
        spectral_smooth = spectral_smooth, sigma = sigma1,
        shape_kernel = sigma2, close = close, close_kernel = close_kernel,
        chunk_size = chunk_size, sample = sample_name
      )
      cached <- list(
        snr = snr, threshold = threshold,
        feature_metadata = partition$display$metadata,
        collapsed = partition$analysis_units
      )
    }
    .filespec_atomic_save_rds(cached, cache_file)
  } else {
    .particle_progress(sample_name, "reuse cached particle means")
  }

  .particle_threshold_state_message(cached$threshold, sample_name,
                                    particle_id_strategy)

  display <- .filespec_particle_display(index, cached$snr, cached$threshold)
  if (!is.null(cached$feature_metadata)) {
    display$metadata <- data.table::copy(cached$feature_metadata)
  }
  display <- .attach_particle_image(display, list(image), list(bottom_left),
                                    list(top_right), 1L)
  plot_outputs <- .particle_pre_match_plots(
    display, output_name, output_dir, outputs, pixel_length, origin,
    sn_threshold_min, sn_threshold_max
  )
  if (is.null(cached$collapsed)) {
    out <- .empty_particle_result(sample_name, x, time_start, outputs,
                                  plot_outputs, output_dir,
                                  output_name = output_name,
                                  pixel_length = pixel_length)
    return(out)
  }

  .particle_progress(sample_name, "library matching")
  proc_map <- .process_for_particle_match(cached$collapsed, library,
                                          process_args)
  proc_map <- .append_particle_matches(
    proc_map, library = library, material_col = material_col,
    library_id_col = library_id_col
  )
  proc_map <- .filter_particle_matches(
    proc_map, material_col = material_col, cor_threshold = cor_threshold,
    label_unknown = label_unknown, remove_materials = remove_materials,
    remove_unknown = remove_unknown
  )
  display <- .join_particle_display_matches(display, proc_map, material_col)
  details <- if ("details" %in% outputs) {
    .particle_details_table(proc_map, output_name, material_col,
                            cor_threshold, pixel_length, origin)
  } else NULL
  summary <- if ("summary" %in% outputs) {
    .particle_summary_table(
      proc_map, output_name, material_col, pixel_length, display
    )
  } else NULL
  plot_outputs <- utils::modifyList(
    plot_outputs,
    .particle_post_match_plots(
      display, proc_map, output_name, output_dir, outputs, material_col,
      pixel_length, origin, cor_threshold
    )
  )
  .particle_progress(sample_name, "outputs")
  elapsed <- Sys.time() - time_start
  if (!is.null(output_dir)) {
    .write_particle_outputs(
      output_dir, output_name, x, proc_map, details, summary, outputs,
      material_col, pixel_length, origin, elapsed
    )
  }
  result <- list(
    sample_id = output_name,
    particle_details_csv = details,
    particle_summary_csv = summary,
    particles_raw_rds = if ("raw" %in% outputs) x else NULL,
    particles_rds = if ("processed" %in% outputs) proc_map else NULL,
    particle_image = plot_outputs$particle_image,
    particle_heatmap = plot_outputs$particle_heatmap,
    particle_heatmap_thresholded = plot_outputs$particle_heatmap_thresholded,
    cor_heatmap = plot_outputs$cor_heatmap,
    sn_histogram = plot_outputs$sn_histogram,
    cor_histogram = plot_outputs$cor_histogram,
    time_rds = if ("time" %in% outputs) elapsed else NULL
  )
  .particle_progress(sample_name, "complete")
  result
}

.filespec_particle_snr <- function(x, index, bands, metric, abs,
                                   spectral_smooth, sigma1, chunk_size,
                                   process = NULL, sample = NULL) {
  if(!is.null(process) && !is.function(process)) {
    stop("'process' must be NULL or a function", call. = FALSE)
  }
  chunk_size <- .filespec_bounded_chunk_size(length(bands), chunk_size)
  chunks <- if(isTRUE(spectral_smooth)) {
    col_chunk <- .filespec_column_chunk_id(index, chunk_size)
    if(is.null(col_chunk)) {
      stop("spectral_smooth requires a complete rectangular row/col grid ",
           "for this FileSpecs region", call. = FALSE)
    }
    split(seq_len(nrow(index)), col_chunk)
  } else {
    .filespec_particle_chunks(x, index, chunk_size)
  }
  out <- rep(NA_real_, nrow(index))
  completed <- cumsum(lengths(chunks))
  started <- proc.time()[["elapsed"]]
  for (i in seq_along(chunks)) {
    rows <- chunks[[i]]
    values <- if (isTRUE(spectral_smooth)) {
      .filespec_smoothed_values(x, index, rows, bands = bands,
                                sigma1 = sigma1)
    } else {
      .filespec_read_values(x, index = rows, bands = bands)
    }
    block <- as_OpenSpecy(
      values$wavenumber, spectra = values$spectra,
      metadata = data.frame(col_id = colnames(values$spectra)),
      coords = "gen_grid", session_id = FALSE, compute_file_id = FALSE
    )
    if(!is.null(process)) {
      block <- process(block)
      if(!inherits(block, "OpenSpecy")) {
        stop("the streamed S/N process must return OpenSpecy", call. = FALSE)
      }
    }
    out[rows] <- sig_noise(block, metric = metric, spatial_smooth = FALSE,
                          abs = abs)
    .particle_stream_progress(
      sample, "streaming signal/noise", i, length(chunks), completed[[i]],
      nrow(index), started
    )
  }
  out
}

# Purely geometric column grouping from row/col grid coordinates, independent
# of backend. `.filespec_particle_chunks()` only applies it for the h5
# backend (matching its original chunking heuristic); halo-based
# spectral_smooth uses it directly for any backend, since correctness of the
# padding math depends only on a complete rectangular grid, not on how the
# backend physically reads pixels.
.filespec_column_chunk_id <- function(index, chunk_size) {
  if (!all(c("row", "col") %in% names(index))) return(NULL)
  rows <- sort(unique(index$row))
  columns <- sort(unique(index$col))
  complete <- nrow(index) == length(rows) * length(columns) &&
    !anyDuplicated(index[, c("row", "col"), with = FALSE])
  if (!isTRUE(complete)) return(NULL)
  columns_per_chunk <- max(1L, floor(as.integer(chunk_size) / length(rows)))
  ceiling(match(index$col, columns) / columns_per_chunk)
}

.filespec_particle_chunks <- function(x, index, chunk_size) {
  positions <- seq_len(nrow(index))
  if (!identical(x$source$backend, "h5")) {
    return(split(positions, ceiling(positions / as.integer(chunk_size))))
  }
  col_chunk <- .filespec_column_chunk_id(index, chunk_size)
  if (is.null(col_chunk)) {
    return(split(positions, ceiling(positions / as.integer(chunk_size))))
  }
  split(positions, col_chunk)
}

# Read a halo-padded column block around `rows` (positions into `index`, a
# single region's complete row x col grid), 3-D Gaussian-smooth it with the
# same mmand::gaussianSmooth() call the eager reader uses, then trim back to
# exactly the requested pixels. The halo equals mmand's own kernel radius for
# `sigma1`, so trimmed values are numerically identical to smoothing the full
# region at once, without ever reading more than one padded column slab.
.filespec_smoothed_values <- function(x, index, rows, bands, sigma1) {
  target <- index[rows]
  if (!all(c("row", "col") %in% names(index))) {
    stop("spectral_smooth requires row/col grid coordinates for this ",
         "FileSpecs region", call. = FALSE)
  }
  rows_all <- sort(unique(index$row))
  cols_all <- sort(unique(index$col))
  complete <- nrow(index) == length(rows_all) * length(cols_all) &&
    !anyDuplicated(index[, c("row", "col"), with = FALSE])
  if (!isTRUE(complete)) {
    stop("spectral_smooth requires a complete rectangular row/col grid for ",
         "this FileSpecs region; irregular regions are not supported",
         call. = FALSE)
  }

  target_cols <- sort(unique(target$col))
  halo <- .gaussian_kernel_half_width(sigma1[[3L]])
  col_lo_i <- max(1L, match(min(target_cols), cols_all) - halo)
  col_hi_i <- min(length(cols_all), match(max(target_cols), cols_all) + halo)
  padded_cols <- cols_all[col_lo_i:col_hi_i]

  block_rows <- which(index$col %in% padded_cols)
  block <- .filespec_read_values(x, index = block_rows, bands = NULL)
  sel <- block$index
  nband <- length(block$wavenumber)

  row_map <- match(sel$row, rows_all)
  col_map <- match(sel$col, padded_cols)
  lin <- (col_map - 1L) * length(rows_all) + row_map
  arr <- matrix(NA_real_, nrow = nband,
                ncol = length(rows_all) * length(padded_cols))
  arr[, lin] <- block$spectra
  dim(arr) <- c(nband, length(rows_all), length(padded_cols))

  smoothed <- mmand::gaussianSmooth(arr, sigma = sigma1)

  band_keep <- if (is.null(bands)) {
    seq_len(nband)
  } else {
    match(x$source$axis[bands], block$wavenumber)
  }
  out_row <- match(target$row, rows_all)
  out_col <- match(target$col, padded_cols)
  out_lin <- (out_col - 1L) * length(rows_all) + out_row
  spectra <- matrix(smoothed, nrow = nband)[band_keep, out_lin, drop = FALSE]
  colnames(spectra) <- target$col_id

  list(wavenumber = block$wavenumber[band_keep], spectra = spectra,
       index = target)
}

.gaussian_kernel_half_width <- function(sigma) {
  if (!is.finite(sigma) || sigma <= 0) return(0L)
  size <- ceiling(6 * sigma)
  if (size %% 2L == 0L) size <- size + 1L
  as.integer((size - 1L) / 2L)
}

.filespec_particle_display <- function(index, snr, threshold) {
  md <- data.table::copy(index)
  md$snr <- as.numeric(snr)
  md$threshold <- as.logical(threshold)
  ids <- if ("col_id" %in% names(md)) as.character(md$col_id) else
    as.character(md$source_id)
  spectra <- matrix(0, nrow = 1L, ncol = nrow(md),
                    dimnames = list("preview", ids))
  as_OpenSpecy(0, spectra = spectra, metadata = md, coords = NULL,
               compute_file_id = FALSE)
}

# Match every eligible file-backed pixel while retaining only its best result.
# Query and reference score matrices are both bounded so all-cell grouping does
# not recreate the complete H5 spectra matrix or a library-by-map matrix.
.filespec_particle_matches <- function(
    x, eligible, library, process_args, material_col, library_id_col,
    spectral_smooth, sigma1, chunk_size, sample = NULL) {
  index <- .filespec_index(x)
  if (!is.logical(eligible) || length(eligible) != nrow(index)) {
    stop("'eligible' must have one logical value per file-backed spectrum",
         call. = FALSE)
  }
  eligible[is.na(eligible)] <- FALSE
  positions <- which(eligible)
  empty <- data.table::data.table(
    object_id = character(), max_cor_val = numeric()
  )
  empty[[material_col]] <- character()
  if (!length(positions)) return(empty)

  chunk_size <- .filespec_bounded_chunk_size(
    length(.filespec_axis(x)), chunk_size
  )
  chunks <- if (isTRUE(spectral_smooth)) {
    col_chunk <- .filespec_column_chunk_id(index, chunk_size)
    if (is.null(col_chunk)) {
      stop("spectral_smooth requires a complete rectangular row/col grid ",
           "for this FileSpecs region", call. = FALSE)
    }
    split(seq_along(positions), col_chunk[positions])
  } else {
    split(seq_along(positions), ceiling(seq_along(positions) / chunk_size))
  }

  prepared <- NULL
  result <- vector("list", length(chunks))
  completed <- cumsum(lengths(chunks))
  started <- proc.time()[["elapsed"]]
  for (i in seq_along(chunks)) {
    rows <- chunks[[i]]
    source_rows <- positions[rows]
    values <- if (isTRUE(spectral_smooth)) {
      .filespec_smoothed_values(
        x, index, source_rows, bands = NULL, sigma1 = sigma1
      )
    } else {
      .filespec_read_values(x, index = source_rows)
    }
    query <- .filespec_values_to_OpenSpecy(x, values)
    query <- .process_for_particle_match(query, library, process_args)

    if (is_OpenSpecy(library)) {
      if (is.null(prepared)) {
        reference <- library
        if (!identical(reference$wavenumber, query$wavenumber)) {
          reference <- conform_spec(
            reference, range = query$wavenumber, res = NULL, allow_na = FALSE
          )
        }
        prepared <- .filespec_prepare_correlation_reference(reference)
      }
      matches <- .filespec_match_prepared_best(query, prepared)
      lib_md <- data.table::as.data.table(library$metadata)
      material_index <- match(matches$library_id, lib_md[[library_id_col]])
      material <- if (all(c(library_id_col, material_col) %in% names(lib_md))) {
        as.character(lib_md[[material_col]][material_index])
      } else {
        as.character(matches$library_id)
      }
    } else {
      prediction <- data.table::as.data.table(match_spec(query, library))
      if (!all(c("x", "name", "value") %in% names(prediction))) {
        stop("the model library returned an unsupported match table",
             call. = FALSE)
      }
      prediction <- prediction[order(x, -value)][!duplicated(x)]
      object_index <- suppressWarnings(as.integer(prediction$x))
      matches <- data.table::data.table(
        object_id = colnames(query$spectra)[object_index],
        library_id = as.character(prediction$name),
        match_val = as.numeric(prediction$value)
      )
      material <- matches$library_id
    }
    aligned <- match(colnames(query$spectra), matches$object_id)
    if (anyNA(aligned)) {
      stop("streamed pixel matches do not align with their source spectra",
           call. = FALSE)
    }
    block <- data.table::data.table(
      object_id = colnames(query$spectra),
      max_cor_val = as.numeric(matches$match_val[aligned])
    )
    block[[material_col]] <- material[aligned]
    result[[i]] <- block
    rm(values, query, matches, block)
    if (i %% 5L == 0L || i == length(chunks)) invisible(gc(verbose = FALSE))
    .particle_stream_progress(
      sample, "streaming pixel identification", i, length(chunks),
      completed[[i]], length(positions), started
    )
  }
  data.table::rbindlist(result, use.names = TRUE)
}

.filespec_prepare_correlation_reference <- function(reference) {
  values <- make_rel(reference$spectra, na.rm = TRUE)
  values <- .matrix_mean_replace(values)
  list(
    wavenumber = reference$wavenumber,
    library_id = colnames(reference$spectra),
    scaled = .scale_correlation_spectra(values)
  )
}

.filespec_match_prepared_best <- function(query, prepared,
                                          library_block_size = 1000L) {
  if (!identical(query$wavenumber, prepared$wavenumber)) {
    stop("processed query and reference axes do not match", call. = FALSE)
  }
  query_values <- make_rel(query$spectra, na.rm = TRUE)
  query_values <- .matrix_mean_replace(query_values)
  scaled_query <- .scale_correlation_spectra(query_values)
  query_count <- nrow(scaled_query)
  best_value <- rep(NA_real_, query_count)
  best_index <- rep(NA_integer_, query_count)
  query_columns <- seq_len(query_count)
  starts <- seq.int(1L, nrow(prepared$scaled), by = library_block_size)

  for (start in starts) {
    rows <- seq.int(
      start, min(nrow(prepared$scaled), start + library_block_size - 1L)
    )
    scores <- tcrossprod(prepared$scaled[rows, , drop = FALSE], scaled_query)
    ranked <- scores
    ranked[!is.finite(ranked)] <- -Inf
    local_index <- max.col(t(ranked), ties.method = "first")
    candidate_value <- scores[cbind(local_index, query_columns)]
    candidate_index <- rows[local_index]
    update <- is.na(best_index) |
      (!is.na(candidate_value) &
         (is.na(best_value) | candidate_value > best_value))
    best_value[update] <- candidate_value[update]
    best_index[update] <- candidate_index[update]
  }
  data.table::data.table(
    object_id = colnames(query$spectra),
    library_id = prepared$library_id[best_index],
    match_val = best_value
  )
}

.filespec_mean_features <- function(x, index, feature_metadata, feature_ids,
                                    axis, spectral_smooth, sigma1,
                                    chunk_size, unit_metadata = NULL,
                                    sample = NULL) {
  chunk_size <- .filespec_bounded_chunk_size(length(axis), chunk_size)
  .filespec_retained_mean_capacity(length(axis), length(feature_ids))
  ids <- as.character(feature_metadata$feature_id)
  keep <- ids %in% feature_ids
  selected <- which(keep)
  sums <- matrix(0, nrow = length(axis), ncol = length(feature_ids),
                 dimnames = list(as.character(axis), feature_ids))
  counts <- integer(length(feature_ids))
  if (isTRUE(spectral_smooth)) {
    col_chunk <- .filespec_column_chunk_id(index, chunk_size)
    if (is.null(col_chunk)) {
      stop("spectral_smooth requires a complete rectangular row/col grid ",
           "for this FileSpecs region", call. = FALSE)
    }
    chunks <- split(selected, col_chunk[selected])
  } else {
    chunks <- split(selected, ceiling(seq_along(selected) /
                                        as.integer(chunk_size)))
  }
  completed <- cumsum(lengths(chunks))
  started <- proc.time()[["elapsed"]]
  for (i in seq_along(chunks)) {
    rows <- chunks[[i]]
    block <- if (isTRUE(spectral_smooth)) {
      .filespec_smoothed_values(x, index, rows, bands = NULL,
                                sigma1 = sigma1)
    } else {
      .filespec_read_values(x, index = rows)
    }
    groups <- match(ids[rows], feature_ids)
    for (group in unique(groups)) {
      cols <- which(groups == group)
      sums[, group] <- sums[, group] + rowSums(block$spectra[, cols,
                                                             drop = FALSE])
      counts[[group]] <- counts[[group]] + length(cols)
    }
    .particle_stream_progress(
      sample, "streaming particle means", i, length(chunks), completed[[i]],
      length(selected), started
    )
  }
  spectra <- sweep(sums, 2L, counts, "/")
  md <- if (is.null(unit_metadata)) {
    data.table::copy(feature_metadata[match(feature_ids, ids)])
  } else {
    unit_metadata <- data.table::as.data.table(unit_metadata)
    id_column <- intersect(c("col_id", "unit_id", "feature_id"),
                           names(unit_metadata))
    if (!length(id_column)) {
      stop("collapsed unit metadata requires an identifier column",
           call. = FALSE)
    }
    unit_index <- match(feature_ids,
                        as.character(unit_metadata[[id_column[[1L]]]]))
    if (anyNA(unit_index)) {
      stop("collapsed unit metadata does not align with streamed features",
           call. = FALSE)
    }
    data.table::copy(unit_metadata[unit_index])
  }
  md$col_id <- feature_ids
  as_OpenSpecy(axis, spectra = spectra, metadata = md, coords = NULL,
               compute_file_id = FALSE)
}

# Collapse a file-backed map through the same connected-region contract as the
# in-memory app path. Only a one-row geometric display and bounded spectral
# blocks are materialized; one running sum/count pair is retained per particle.
.filespec_collapse_connected_mean <- function(
    x, eligible, material = NULL, snr = NULL, max_cor_val = NULL,
    area_threshold = 1, spectral_smooth = FALSE, sigma = c(1, 1, 1),
    shape_kernel = c(3, 3), close = FALSE, close_kernel = c(4, 4),
    chunk_size = 8192L, sample = NULL) {
  started <- proc.time()[["elapsed"]]
  .filespec_validate_object(x)
  index <- data.table::copy(.filespec_index(x))
  if (!is.logical(eligible) || length(eligible) != nrow(index)) {
    stop("'eligible' must have one logical value per file-backed spectrum",
         call. = FALSE)
  }
  eligible[is.na(eligible)] <- FALSE
  if(!is.null(material) && length(material) != nrow(index)) {
    stop("'material' must have one value per file-backed spectrum",
         call. = FALSE)
  }
  if (is.null(snr)) snr <- rep(NA_real_, nrow(index))
  if (length(snr) != nrow(index)) {
    stop("'snr' must have one value per file-backed spectrum", call. = FALSE)
  }
  if (!is.null(max_cor_val) && length(max_cor_val) != nrow(index)) {
    stop("'max_cor_val' must have one value per file-backed spectrum",
         call. = FALSE)
  }
  display <- .filespec_particle_display(
    index, snr = snr, threshold = eligible
  )
  if (!is.null(max_cor_val)) display$metadata$max_cor_val <- max_cor_val
  partition <- .partition_particle_map(
    display, eligible = eligible, strategy = "collapse",
    material = material, collapse_function = base::mean,
    area_threshold = area_threshold, shape_kernel = shape_kernel,
    close = close, close_kernel = close_kernel
  )
  message(
    "FileSpecs collapse: ", nrow(index), " source spectra; ",
    sum(eligible), " retained; ",
    sum(unique(stats::na.omit(partition$pixel_to_unit$unit_index)) > 0L),
    " connected units."
  )
  if (is.null(partition$analysis_units)) return(partition)

  feature_ids <- colnames(partition$analysis_units$spectra)
  partition$analysis_units <- .filespec_mean_features(
    x, index = index, feature_metadata = partition$display$metadata,
    feature_ids = feature_ids, axis = .filespec_axis(x),
    spectral_smooth = spectral_smooth, sigma1 = sigma,
    chunk_size = chunk_size,
    unit_metadata = partition$analysis_units$metadata,
    sample = sample
  )
  partition$settings$file_backed <- TRUE
  partition$settings$chunk_size <- .filespec_bounded_chunk_size(
    length(.filespec_axis(x)), chunk_size
  )
  partition$settings$elapsed_seconds <-
    proc.time()[["elapsed"]] - started
  message(
    "FileSpecs collapse complete in ",
    sprintf("%.2f", partition$settings$elapsed_seconds), " seconds."
  )
  partition
}

.filespec_bounded_chunk_size <- function(n_bands, requested,
                                         max_bytes = getOption(
                                           "OpenSpecy.filespec.max_block_bytes",
                                           64 * 1024^2)) {
  n_bands <- suppressWarnings(as.integer(n_bands))
  requested <- suppressWarnings(as.integer(requested))
  max_bytes <- suppressWarnings(as.numeric(max_bytes))
  if(length(n_bands) != 1L || is.na(n_bands) || n_bands < 1L ||
     length(requested) != 1L || is.na(requested) || requested < 1L ||
     length(max_bytes) != 1L || is.na(max_bytes) || !is.finite(max_bytes) ||
     max_bytes < n_bands * 8) {
    stop("FileSpecs block settings cannot fit one spectrum in the memory bound",
         call. = FALSE)
  }
  as.integer(min(requested, floor(max_bytes / (n_bands * 8))))
}

.filespec_retained_mean_capacity <- function(
    n_bands, n_features,
    max_bytes = getOption("OpenSpecy.filespec.max_live_bytes", 768 * 1024^2)) {
  bytes <- as.double(n_bands) * as.double(n_features) * 8 * 2
  max_bytes <- suppressWarnings(as.numeric(max_bytes))
  if(length(max_bytes) != 1L || is.na(max_bytes) || !is.finite(max_bytes) ||
     max_bytes <= 0) max_bytes <- 768 * 1024^2
  if(!is.finite(bytes) || bytes > max_bytes) {
    stop(
      "Collapsed particle means would exceed the file-backed live-memory ",
      "bound. Increase the minimum particle area or split the map.",
      call. = FALSE
    )
  }
  list(bytes = bytes, max_bytes = max_bytes)
}

.filespec_image_identity <- function(image, bottom_left, top_right) {
  if (is.character(image) && length(image) == 1L && file.exists(image)) {
    info <- file.info(image)
    image <- list(path = normalizePath(image, winslash = "/"),
                  size = info$size, mtime = info$mtime,
                  sha256 = digest::digest(image, algo = "sha256", file = TRUE))
  } else if (!is.null(image)) {
    image <- list(class = class(image), dim = dim(image),
                  sha256 = digest::digest(image, algo = "sha256",
                                          serialize = TRUE))
  }
  list(image = image, bottom_left = bottom_left, top_right = top_right)
}