---
title: "Publication-quality activation tables"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Publication-quality activation tables}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  dpi = 96
)
```

## Overview

`cluster_table()` is fmrireport's standalone function for building
publication-quality activation tables from any statistical brain volume.
It follows COBIDAS reporting guidelines: hierarchical cluster/sub-peak
structure, world coordinates, p-values, and optional anatomical labeling.

This vignette walks through creating realistic synthetic activations,
detecting clusters, formatting tables, and adding atlas labels.

## Creating synthetic activation maps

Real fMRI analyses produce `NeuroVol` objects from `fmrireg` or similar
packages. For demonstration, we'll build a synthetic t-statistic volume
with Gaussian-shaped activation blobs placed at known locations.

```{r make-synthetic}
library(fmrireport)
library(neuroim2)

# Define a brain-like volume in MNI-ish space
# 2mm isotropic, 91x109x91 grid, origin at (-90, -126, -72)
sp <- NeuroSpace(c(91, 109, 91),
                 spacing = c(2, 2, 2),
                 origin = c(-90, -126, -72))

# Start with Gaussian noise (simulating a t-stat map under the null)
set.seed(2024)
arr <- array(rnorm(prod(dim(sp))), dim(sp))
```

Now add Gaussian-shaped activation blobs. We'll place them at grid
coordinates and use a smooth falloff to simulate realistic spatial
extent:

```{r add-blobs}
# Helper: add a 3D Gaussian blob centered at grid coordinate (cx, cy, cz)
add_gaussian_blob <- function(arr, cx, cy, cz, peak_value, sigma = 4) {
  dims <- dim(arr)
  # Work in a local window for efficiency
  radius <- ceiling(3 * sigma)
  x_range <- max(1, cx - radius):min(dims[1], cx + radius)
  y_range <- max(1, cy - radius):min(dims[2], cy + radius)
  z_range <- max(1, cz - radius):min(dims[3], cz + radius)

  for (xi in x_range) {
    for (yi in y_range) {
      for (zi in z_range) {
        d2 <- (xi - cx)^2 + (yi - cy)^2 + (zi - cz)^2
        arr[xi, yi, zi] <- arr[xi, yi, zi] +
          peak_value * exp(-d2 / (2 * sigma^2))
      }
    }
  }
  arr
}

# Place blobs at grid coordinates that map to known cortical MNI regions.
# With origin=(-90,-126,-72) and spacing=2mm:
#   grid (27,51,65) -> MNI (-38,-26,56) = Left SomMot
#   grid (65,51,65) -> MNI ( 38,-26,56) = Right SomMot
#   grid (46,55,67) -> MNI (  0,-18,62) = Medial motor
#   grid (41,21,39) -> MNI (-10,-86, 4) = Left V1
#   grid (51,21,39) -> MNI ( 10,-86, 4) = Right V1

# Blob 1: Left somatomotor cortex (large, strong activation)
arr <- add_gaussian_blob(arr, cx = 27, cy = 51, cz = 65, peak_value = 8, sigma = 5)

# Blob 2: Right somatomotor cortex (moderate)
arr <- add_gaussian_blob(arr, cx = 65, cy = 51, cz = 65, peak_value = 6, sigma = 4)

# Blob 3: Medial motor / SMA (moderate)
arr <- add_gaussian_blob(arr, cx = 46, cy = 55, cz = 67, peak_value = 5.5, sigma = 3)

# Blob 4: Left primary visual cortex (small, strong)
arr <- add_gaussian_blob(arr, cx = 41, cy = 21, cz = 39, peak_value = 7, sigma = 3)

# Blob 5: Right primary visual cortex (small, moderate)
arr <- add_gaussian_blob(arr, cx = 51, cy = 21, cz = 39, peak_value = 5, sigma = 3)

vol <- NeuroVol(arr, sp)
```

Let's visualize the synthetic activation map:

```{r plot-montage, fig.height = 3.5}
# Show axial slices at the levels where we placed blobs
slices <- c(39, 55, 60, 65, 67)
plot_montage(vol, zlevels = slices, along = 3L,
             cmap = "inferno", ncol = 5L,
             title = "Synthetic t-statistic map")
```

## Basic cluster detection

Detect clusters at t > 3.0 with default settings:

```{r basic-cluster}
ct <- cluster_table(vol,
                    threshold = 3.0,
                    stat_type = "t",
                    df = 100,
                    min_cluster_size = 20)
print(ct)
```

The output shows:

- **Cluster ID** and **size** (voxels and mm^3)
- **Peak coordinates** in world space (mm)
- **Peak statistic** and **p-value**
- **Sub-peaks** (local maxima) indented below each cluster

## Controlling sub-peak detection

The `local_maxima_dist` parameter sets the minimum distance (mm) between
reported sub-peaks. A larger value gives fewer, more distinct peaks:

```{r subpeak-dist}
# Tight distance: many sub-peaks
ct_tight <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                          local_maxima_dist = 8, max_peaks = 5)
cat("Sub-peaks with 8mm distance:\n")
nrow(ct_tight$peaks)

# Wide distance: fewer sub-peaks
ct_wide <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                         local_maxima_dist = 25, max_peaks = 5)
cat("Sub-peaks with 25mm distance:\n")
nrow(ct_wide$peaks)
```

## Sorting options

By default, clusters are sorted by size (descending). You can sort by
peak statistic instead:

```{r sort-options}
ct_by_stat <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                            sort_by = "peak_stat")
# First cluster has the highest absolute peak stat
ct_by_stat$clusters[, c("cluster_id", "k", "peak_stat")]
```

## P-value computation

`cluster_table()` computes voxel-level p-values from the peak statistic
and degrees of freedom. Different statistic types are supported:

```{r pvalues}
# t-statistic
ct_t <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100)
head(ct_t$clusters[, c("peak_stat", "peak_p")])

# z-statistic (provide df to enable p-value computation)
ct_z <- cluster_table(vol, threshold = 3.0, stat_type = "z", df = 1)
head(ct_z$clusters[, c("peak_stat", "peak_p")])

# Without df, p-values are NA
ct_nodf <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = NULL)
head(ct_nodf$clusters[, c("peak_stat", "peak_p")])
```

## Flat data.frame output

`as.data.frame()` produces a table where cluster-level rows have all
columns populated and sub-peak rows have `NA` for `k` and `Vol_mm3`:

```{r flat-df}
df <- as.data.frame(ct)
df
```

## MNI coordinate column names

When you specify a known coordinate space, column headers adapt:

```{r mni-columns}
ct_mni <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                        coord_space = "MNI152")
df_mni <- as.data.frame(ct_mni)
names(df_mni)
```

## Formatted tables with tinytable

For inclusion in Quarto/R Markdown reports, `format_cluster_tt()` produces
a `tinytable` object with an informative caption:

```{r tinytable-output, eval = requireNamespace("tinytable", quietly = TRUE)}
tt <- format_cluster_tt(ct_mni, max_sub_peaks = 2, digits = 2)
tt
```

## Atlas labeling

When you provide an atlas, `cluster_table()` automatically labels each peak
with anatomical region names, hemisphere, and (for parcellation atlases like
Schaefer) network assignment.

```{r atlas-label, eval = requireNamespace("neuroatlas", quietly = TRUE)}
library(neuroatlas)

# Load the Schaefer 200-parcel, 7-network atlas
atlas <- get_schaefer_atlas(parcels = "200", networks = "7")

ct_labeled <- cluster_table(vol, threshold = 3.0,
                            stat_type = "t", df = 100,
                            atlas = atlas,
                            coord_space = "MNI152")
print(ct_labeled)
```

The labeled table now includes region name, hemisphere, and network for
each peak coordinate:

```{r atlas-df, eval = requireNamespace("neuroatlas", quietly = TRUE)}
as.data.frame(ct_labeled)
```

## Combining with overlay plots

You can pair the cluster table with a brain overlay visualization:

```{r overlay-plot, fig.height = 4}
# Create a thresholded overlay for visualization
arr_thresh <- as.array(vol)
arr_thresh[abs(arr_thresh) < 3.0] <- NA
vol_thresh <- NeuroVol(arr_thresh, sp)

# Create a background volume (uniform gray for this demo)
bg_arr <- array(1000, dim(sp))
bg <- NeuroVol(bg_arr, sp)

plot_overlay(bg, vol_thresh,
             zlevels = c(39, 55, 60, 65),
             along = 3L,
             ov_cmap = "inferno",
             ov_thresh = 0,
             ncol = 4L,
             title = "Thresholded activation (t > 3.0)")
```

## Connectivity options

By default, clusters use 26-connectivity (face, edge, and vertex neighbors).
You can switch to more conservative connectivity:

```{r connectivity}
ct_6 <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                      connectivity = "6-connect")
ct_26 <- cluster_table(vol, threshold = 3.0, stat_type = "t", df = 100,
                       connectivity = "26-connect")
cat("6-connect clusters:", ct_6$n_clusters, "\n")
cat("26-connect clusters:", ct_26$n_clusters, "\n")
```

Stricter connectivity (6-connect) tends to break large clusters into
smaller pieces.

## Edge cases

The function handles edge cases gracefully:

```{r edge-cases}
# No suprathreshold voxels
ct_empty <- cluster_table(vol, threshold = 100)
ct_empty$n_clusters
print(ct_empty)

# Very liberal threshold (many small clusters filtered by min_cluster_size)
ct_liberal <- cluster_table(vol, threshold = 1.5, min_cluster_size = 50)
cat("Clusters with k >= 50:", ct_liberal$n_clusters, "\n")
```

## The cluster_table object

A `cluster_table` is an S3 list with these elements:

```{r structure}
str(ct, max.level = 1)
```

| Element | Type | Description |
|---|---|---|
| `clusters` | data.frame | One row per cluster |
| `peaks` | data.frame | Sub-peaks with rank within cluster |
| `threshold` | numeric | Detection threshold used |
| `stat_type` | character | `"t"`, `"z"`, `"F"`, or `"other"` |
| `df` | numeric | Degrees of freedom (NULL if not provided) |
| `n_clusters` | integer | Number of clusters |
| `coord_space` | character | Coordinate space label |
| `atlas_name` | character | Atlas name or NA |

## Summary

`cluster_table()` provides a complete pipeline from statistical volume to
publication-ready activation table:
1. Cluster detection via `neuroim2::conn_comp()`
2. Size filtering
3. World coordinate transformation
4. Sub-peak (local maxima) identification
5. Optional atlas labeling via `neuroatlas`
6. Parametric p-value computation

Output can be printed to the console, converted to a flat data.frame, or
formatted as a tinytable for direct inclusion in reports.
