---
title: Dilating Atlas Parcels Safely
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: interaction
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Dilating Atlas Parcels Safely}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) {
  knitr::opts_chunk$set(dev = "ragg_png")
}
if (requireNamespace("systemfonts", quietly = TRUE) &&
      requireNamespace("albersdown", quietly = TRUE)) {
  albersdown::albers_register_fonts()
}
if (requireNamespace("ggplot2", quietly = TRUE) &&
      requireNamespace("albersdown", quietly = TRUE)) {
  ggplot2::theme_set(
    albersdown::theme_albers(
      family = params$family,
      preset = params$preset
    )
  )
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 8,
  fig.height = 4.5
)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    "});</script>"
  ),
  params$family,
  params$preset
))
```


A parcellation can leave small gaps between labelled regions. If those gaps
are valid analysis voxels, regional summaries silently omit them. Dilation can
close the gaps, but only after you state two contracts: which voxels are
eligible, and how far a parcel may grow.

This article uses a tiny synthetic atlas so every assignment can be inspected
and every claim is executed offline. The same rules apply to a whole-brain
atlas.

## What does dilation do?

`dilate_atlas()` considers voxels that are inside a supplied mask but not
already labelled. A candidate is assigned only when a labelled voxel lies
within `radius`, measured in **voxel units**. When several parcels are nearby,
inverse-distance-weighted voting over at most `maxn` neighbours chooses the
label.

The function never decides what anatomy is appropriate. That is the mask's
job.

## How do you build a transparent example?

Our atlas has two one-voxel parcels on a straight line. IDs are deliberately
non-contiguous so a label-mapping error cannot pass unnoticed.

```{r make-atlas}
library(neuroatlas)

dims <- c(9L, 9L, 5L)
space <- neuroim2::NeuroSpace(
  dim = dims,
  spacing = c(2, 2, 2),
  origin = c(0, 0, 0)
)

labels <- array(0L, dim = dims)
labels[2, 5, 3] <- 1L
labels[8, 5, 3] <- 5L

toy_atlas <- structure(
  list(
    name = "two-parcel example",
    atlas = neuroim2::NeuroVol(labels, space),
    ids = c(1L, 5L),
    labels = c("A", "B"),
    orig_labels = c("A", "B"),
    hemi = c(NA_character_, NA_character_)
  ),
  class = c("toy_atlas", "atlas")
)
```

The eligible mask is a one-voxel-wide bridge between A and B. One additional
voxel sits far away and demonstrates the radius guard.

```{r make-mask}
eligible <- array(FALSE, dim = dims)
eligible[2:8, 5, 3] <- TRUE
eligible[9, 9, 5] <- TRUE

mask <- neuroim2::LogicalNeuroVol(eligible, space)
```

```{r validate-inputs, include = FALSE}
stopifnot(
  identical(dim(toy_atlas$atlas), dim(mask)),
  identical(neuroim2::space(toy_atlas$atlas), neuroim2::space(mask)),
  sum(labels > 0) == 2L,
  sum(eligible) == 8L
)
```

## What changes with a small radius?

A radius of one voxel grows each parcel by one grid step. The middle of the
bridge and the distant candidate remain unassigned.

```{r dilate-one}
radius_one <- dilate_atlas(
  toy_atlas,
  mask,
  radius = 1,
  maxn = 2
)

before <- neuroim2::as.dense(toy_atlas$atlas)
after_one <- neuroim2::as.dense(radius_one$atlas)

c(
  labelled_before = sum(before > 0),
  labelled_after = sum(after_one > 0),
  eligible = sum(eligible)
)
```

```{r validate-one, include = FALSE}
stopifnot(
  after_one[3, 5, 3] == 1L,
  after_one[7, 5, 3] == 5L,
  after_one[5, 5, 3] == 0L,
  after_one[9, 9, 5] == 0L,
  setequal(unique(after_one[after_one > 0]), toy_atlas$ids)
)
```

## What changes with a larger radius?

At three voxels, every point along the bridge can reach a parcel. The distant
voxel is still protected because it lies more than three grid steps from the
nearest label.

```{r dilate-three}
radius_three <- dilate_atlas(
  toy_atlas,
  mask,
  radius = 3,
  maxn = 2
)

after_three <- neuroim2::as.dense(radius_three$atlas)

c(
  labelled_before = sum(before > 0),
  labelled_after = sum(after_three > 0),
  eligible = sum(eligible)
)
```

```{r validate-three, include = FALSE}
stopifnot(
  all(after_three[2:8, 5, 3] > 0),
  after_three[3, 5, 3] == 1L,
  after_three[7, 5, 3] == 5L,
  after_three[9, 9, 5] == 0L,
  setequal(unique(after_three[after_three > 0]), toy_atlas$ids)
)
```

The radius changed proximity, not anatomy. If the far voxel had been close
enough, the algorithm would have labelled it because the mask declared it
eligible.

## How do you see the assignment?

The table below shows the bridge from left to right. Zero means unassigned.

```{r assignment-table}
data.frame(
  x = 2:8,
  before = before[2:8, 5, 3],
  radius_1 = after_one[2:8, 5, 3],
  radius_3 = after_three[2:8, 5, 3]
)
```

The midpoint is equidistant. Its label is a deterministic consequence of the
grouped label ordering used to break an equal vote, so do not interpret a tie
as anatomical evidence. In real work, inspect parcel boundaries and quantify
sensitivity to radius.

## What must be true for real data?

Before dilating a real atlas, verify all of the following:

1. **Same template space.** The atlas and mask describe the same anatomy.
2. **Same voxel grid.** Dimensions and spatial geometry match exactly.
3. **Appropriate mask.** A cortical atlas should use a cortical ribbon or an
   equivalently justified mask, not an unrestricted whole-grey-matter mask.
4. **Interpretable radius.** Convert voxel units to millimetres using the grid
   spacing, and report both.
5. **Stable identities.** The output contains only the input parcel IDs.
6. **Sensitivity.** Report how many voxels are added and where they are added
   for plausible radii.

Grid resampling does not transform one anatomical template into another. For
example, the nonlinear route between `MNI152NLin6Asym` and
`MNI152NLin2009cAsym` is currently recorded as planned:

```{r transform-status}
atlas_transform_plan(
  "MNI152NLin6Asym",
  "MNI152NLin2009cAsym",
  data_type = "voxel"
)
```

Do not build a mask in one of those spaces, reslice it onto the other's grid,
and treat the result as registered anatomy unless an actual validated transform
has been applied.

## How do you use the result?

The returned object preserves the atlas class, IDs, and labels, so it can enter
the same volume workflows as its input:

```{r downstream}
image <- neuroim2::NeuroVol(
  array(seq_len(prod(dims)), dim = dims),
  space
)

reduce_atlas(radius_three, image, mean)
```

```{r validate-downstream, include = FALSE}
summary <- reduce_atlas(radius_three, image, mean)
stopifnot(
  nrow(summary) == length(toy_atlas$ids),
  all(is.finite(summary$value))
)
```

Return to `vignette("neuroatlas-overview")` for the general reduction workflow,
or read `vignette("working-with-templateflow")` before combining external
templates and masks.