---
title: 'End-to-end: From NIfTI to clusters'
name: end-to-end-export
description: Read 4D data and a mask, fit a checked partition, and round-trip its NIfTI label map.
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  base_size: 13
  content_width: 72
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts
vignette: |
  %\VignetteIndexEntry{End-to-end: From NIfTI to clusters}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
if (requireNamespace("albersdown", quietly = TRUE)) albersdown::albers_register_fonts()
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.align = "center", fig.retina = 2,
  out.width = "100%", fig.width = 8.5, fig.asp = 0.56,
  message = FALSE, warning = FALSE
)
options(width = 88)
suppressPackageStartupMessages({
  library(neurocluster)
  library(neuroim2)
})
parcel_colors <- c("#0072B2", "#E69F00", "#009E73", "#CC79A7",
                   "#56B4E9", "#D55E00", "#F0E442", "#000000")
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(paste0(
  '<script>document.addEventListener("DOMContentLoaded",function(){',
  'document.body.classList.add("palette-%s","preset-%s");',
  '});</script><style>:root{--content:%sch}</style>'
), params$family, params$preset, params$content_width))
```

# What does the complete workflow do? {#workflow}

The production path has five explicit boundaries: read the 4D image and mask,
check that their geometry matches, fit a partition, validate the returned object
against those inputs, and write then reread the label map. This article executes
every boundary using temporary NIfTI files created from a deterministic fixture.

For your analysis, replace only `vec_file` and `mask_file` with your own paths.
The remaining calls are the same.

# How do you create a reproducible input boundary? {#fixture}

The generated data contain four spatial regions with distinct time series. We
write them first so the example exercises the same I/O path as external files.

```{r write-inputs}
base_source <- generate_synthetic_volume(
  scenario = "gaussian_blobs", dims = c(16, 16, 8),
  n_clusters = 4, n_time = 40, noise_sd = 0.08, seed = 42
)
dims <- base_source$dims
ijk <- arrayInd(seq_len(prod(dims)), dims)
inside <- ((ijk[, 1] - 8.5) / 7.5)^2 +
  ((ijk[, 2] - 8.5) / 7.0)^2 + ((ijk[, 3] - 4.5) / 3.8)^2 <= 1
mask_array <- array(as.integer(inside), dims)
spatial_space <- neuroim2::NeuroSpace(
  dims, spacing = c(2, 2.5, 3), origin = c(-18, 7, 11)
)
series_space <- neuroim2::NeuroSpace(
  c(dims, 40L), spacing = c(2, 2.5, 3), origin = c(-18, 7, 11)
)
series_values <- array(
  as.numeric(base_source$vec), dim = dim(base_source$vec)
)
source_data <- list(
  vec = neuroim2::NeuroVec(series_values, series_space),
  mask = neuroim2::NeuroVol(mask_array, spatial_space),
  truth = base_source$truth[inside], dims = dims
)
stopifnot(all(sort(unique(source_data$truth)) == 1:4), any(!inside), any(inside))
vec_file <- file.path(tempdir(), "neurocluster-input-4d.nii.gz")
mask_file <- file.path(tempdir(), "neurocluster-input-mask.nii.gz")
neuroim2::write_vec(source_data$vec, vec_file)
neuroim2::write_vol(source_data$mask, mask_file)
```

# How do you read and check the data? {#read-inputs}

Use `read_vec()` for a 4D image and `read_vol()` for a 3D mask. The mask controls
both which voxels enter clustering and the order of returned labels.

```{r read-inputs}
vec <- neuroim2::read_vec(vec_file)
mask <- neuroim2::read_vol(mask_file)
input_max_abs_error <- max(abs(as.array(vec) - as.array(source_data$vec)))
c(vec_dimensions = paste(dim(vec), collapse = " x "),
  mask_dimensions = paste(dim(mask), collapse = " x "),
  included_voxels = sum(as.array(mask) > 0),
  max_signal_roundtrip_error = signif(input_max_abs_error, 3))
```

```{r check-inputs, include=FALSE}
stopifnot(
  identical(dim(vec)[1:3], dim(mask)),
  identical(dim(vec), c(16L, 16L, 8L, 40L)),
  isTRUE(all.equal(neuroim2::space(vec), neuroim2::space(source_data$vec))),
  isTRUE(all.equal(neuroim2::space(mask), neuroim2::space(source_data$mask))),
  isTRUE(all.equal(neuroim2::trans(neuroim2::space(vec)),
                   neuroim2::trans(neuroim2::space(mask)))),
  input_max_abs_error < 2e-6,
  identical(as.integer(as.array(mask)),
            as.integer(as.array(source_data$mask))),
  all(is.finite(as.array(mask))),
  sum(as.array(mask) > 0) == length(source_data$truth)
)
```

# How do you fit and validate the partition? {#fit}

ReNA is used here to keep the first production-style path concise. Method choice
and tuning should be justified separately; see [Compare methods](compare-methods.html)
and [Choose parameters](choose-parameters.html).

```{r fit-partition}
fit <- cluster4d(
  vec, mask, n_clusters = 4, method = "rena",
  connectivity = 26, parallel = FALSE
)
validation <- validate_cluster4d(fit, vec, mask)
c(valid = validation$valid, requested_k = 4, actual_k = fit$actual_k)
```

```{r measure-partition}
metrics <- cluster_metrics(fit, truth = source_data$truth)
metrics[c("ari_truth", "size_summary", "spatial_rms_distance_mm")]
```

```{r check-fit, include=FALSE}
stopifnot(
  validation$valid,
  length(validation$errors) == 0L,
  identical(fit$actual_k, 4L),
  is.finite(metrics$ari_truth), metrics$ari_truth > 0.75,
  all(metrics$size_summary > 0),
  is.finite(metrics$spatial_rms_distance_mm)
)
voxel_series <- as.matrix(neuroim2::series(
  vec, which(as.array(mask) > 0)
))
temporal_sd <- apply(voxel_series, 2, sd)
sd_array <- array(NA_real_, dim(mask))
sd_array[which(as.array(mask) > 0)] <- temporal_sd
label_array <- as.array(fit$clusvol)
```

```{r input-output-figure, echo=FALSE, fig.cap='The input temporal variability and fitted parcel map at the same middle axial slice. The quantitative color bar reports temporal standard deviation; the parcel colors are categorical and intentionally use a different scale.', fig.alt='An axial heat map of temporal standard deviation with a numbered color bar beside a categorical four-parcel map, all within an ellipsoidal mask.'}
old_par <- par(no.readonly = TRUE)
layout(matrix(c(1, 2, 3), nrow = 1), widths = c(1, 0.16, 1))
par(mar = c(1, 1, 3, 1))
sd_slice <- sd_array[, , 4]
sd_colors <- hcl.colors(64, "YlOrRd", rev = TRUE)
sd_limits <- range(sd_slice, na.rm = TRUE)
image(t(sd_slice[, ncol(sd_slice):1]), col = sd_colors, zlim = sd_limits,
      axes = FALSE, useRaster = TRUE, asp = 1, main = "Input temporal SD")
box(col = "#4A4A4A")
par(mar = c(1, 0, 3, 3))
legend_values <- seq(sd_limits[1], sd_limits[2], length.out = length(sd_colors))
image(1, legend_values, matrix(legend_values, nrow = 1), col = sd_colors,
      axes = FALSE, xlab = "", ylab = "", useRaster = TRUE)
axis(4, las = 1, cex.axis = 0.8)
par(mar = c(1, 1, 3, 1))
label_slice <- label_array[, , 4]
image(t(label_slice[, ncol(label_slice):1]),
      col = parcel_colors[seq_len(fit$actual_k)],
      zlim = c(0.5, fit$actual_k + 0.5), axes = FALSE,
      useRaster = TRUE, asp = 1, main = "Fitted parcels")
box(col = "#4A4A4A")
layout(1)
par(old_par)
```

# How do you export and prove the round trip? {#export}

Write the `ClusteredNeuroVol` directly. It carries the source geometry and final
integer labels. Then read the file back and compare every voxel with the object
that produced it.

```{r export-result}
output_file <- file.path(tempdir(), "neurocluster-parcels.nii.gz")
neuroim2::write_vol(fit$clusvol, output_file)
exported <- neuroim2::read_vol(output_file)
c(path = output_file, bytes = file.info(output_file)$size,
  labels = paste(sort(unique(as.integer(as.array(exported)))), collapse = ", "))
```

```{r check-export, include=FALSE}
stopifnot(
  file.exists(output_file),
  identical(dim(exported), dim(mask)),
  isTRUE(all.equal(neuroim2::space(exported), neuroim2::space(mask))),
  identical(as.integer(as.array(exported)),
            as.integer(as.array(fit$clusvol))),
  identical(sort(unique(as.integer(as.array(exported)))), c(0L, seq_len(fit$actual_k))),
  all(as.array(exported)[as.array(mask) == 0] == 0)
)
unlink(c(vec_file, mask_file, output_file))
```

# What should you preserve with the result? {#provenance}

Keep the input mask, image geometry, preprocessing record, package version,
method, complete parameters, requested and actual K, validation outcome, and
selection diagnostics. The NIfTI label map is an exchange artifact; the full
`cluster4d_result` retains richer feature-center and provenance information for
auditing and downstream R workflows.