---
title: Advanced 3D Volume Patterns
date: '`r Sys.Date()`'
output: rmarkdown::html_vignette
vignette: |
  %\VignetteIndexEntry{Advanced 3D Volume Patterns}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
params:
  family: red
  preset: homage
css: albers.css
resource_files:
- albers.css
- albers.js
includes:
  in_header: |-
    <script src="albers.js"></script>
    <script>document.addEventListener('DOMContentLoaded',function(){document.body.classList.add('palette-red');});</script>

---

```{r setup, echo = FALSE, message = FALSE}
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)
suppressPackageStartupMessages(library(neuroim2))
```

This article is now the detailed follow-on to `vignette("VolumesAndVectors")`.
Read that vignette first if you want the shortest introduction to `read_vol()`,
`read_vec()`, and the basic `NeuroVol` / `NeuroVec` mental model.

Use this article when you specifically want deeper 3D-volume details: masks,
coordinate conversion, manual construction, and slice-level inspection.

## Read one volume to establish context

```{r}
file_name <- system.file("extdata", "global_mask2.nii.gz", package = "neuroim2")
vol <- read_vol(file_name)
```

This article assumes you already know the basic `NeuroVol` story from
`vignette("VolumesAndVectors")`. The remaining sections focus on patterns that
are specific to 3D work.

```{r}
class(vol)
is.array(vol)
dim(vol)
vol[1, 1, 1]
vol[64, 64, 24]
```

## Coordinate conversion and spatial metadata

```{r}
sp <- space(vol)
sp
dim(vol)
spacing(vol)
origin(vol)
```

You can convert between indices, voxel grid coordinates, and real-world coordinates:

```{r}
idx <- 1:5
g <- index_to_grid(vol, idx)
w <- index_to_coord(vol, idx)
idx2 <- coord_to_index(vol, w)
all.equal(idx, idx2)
```

A numeric image volume can be converted to a binary image as follows:

```{r}
vol2 <- as.logical(vol)
class(vol2)
print(vol2[1, 1, 1])
```

## Masks and LogicalNeuroVol

Create a mask from a threshold or an explicit set of indices. Masks are `LogicalNeuroVol` and align with the 3D space.

```{r}
mask1 <- as.mask(vol > 0.5)
mask1

idx_hi <- which(vol > 0.8)
mask2 <- as.mask(vol, idx_hi)
sum(mask2) == length(idx_hi)

mean_in_mask <- mean(vol[mask1@.Data])
mean_in_mask
```

## Constructing volumes manually

We can also create a `NeuroVol` instance from an `array` or `numeric` vector.
First we construct a standard R `array`:


```{r}
    x <- array(0, c(64,64,64))
```

Now we create a `NeuroSpace` instance that describes the geometry of the image,
including at minimum its dimensions and voxel spacing.

```{r}
    bspace <- NeuroSpace(dim=c(64,64,64), spacing=c(1,1,1))
    vol <- NeuroVol(x, bspace)
    vol
```

We do not usually have to create `NeuroSpace` objects by hand because real
image files carry this information in their headers. In practice you usually
copy an existing space:


```{r}
    vol2 <- NeuroVol((vol+1)*25, space(vol))
    max(vol2)
    space(vol2)
```

## Slicing and quick visualization

The easiest way to view a volume is with `plot()`, which shows a 3 x 3 montage
of evenly-spaced axial slices:

```{r plot_montage, fig.alt='3x3 montage of axial slices.', fig.cap='Default plot() montage', fig.width=7, fig.height=7}
plot(vol)
```

You can also extract a single 2D slice for display using standard array indexing:

```{r slice_mid, fig.alt='Mid-slice of example volume (grayscale image).', fig.cap='Mid-slice of example volume'}
    z <- ceiling(dim(vol)[3] / 2)
    image(vol[,,z], main = paste("Slice z=", z))
```

## Reorienting and resampling

You can change an image’s orientation and voxel spacing. Use `reorient()` to remap axes (e.g., to RAS) and `resample_to()` to match a target space.

```{r}
    # Reorient the space (LPI -> RAS) and compare coordinate mappings
    sp_lpi <- space(vol)
    sp_ras <- reorient(sp_lpi, c("R","A","S"))
    g     <- t(matrix(c(10, 10, 10)))
    world_lpi <- grid_to_coord(sp_lpi, g)
    world_ras <- grid_to_coord(sp_ras, g)
    # world_lpi and world_ras differ due to axis remapping
```

Resample to a new spacing or match a target `NeuroSpace`:

```{r eval=FALSE}
    # Create a target space with 2x finer resolution
    sp  <- space(vol)
    sp2 <- NeuroSpace(sp@dim * c(2,2,2), sp@spacing/2, origin=sp@origin, trans=trans(vol))

    # Resample (trilinear)
    vol_resamp <- resample_to(vol, sp2, method = "linear")
    dim(vol_resamp)
```

## Downsampling

Reduce spatial resolution to speed up downstream operations.

```{r}
    # Downsample by target spacing
    vol_ds1 <- downsample(vol, spacing = spacing(vol)[1:3] * 2)
    dim(vol_ds1)

    # Or by factor
    vol_ds2 <- downsample(vol, factor = 0.5)
    dim(vol_ds2)
```

## Writing a NIFTI formatted image volume

When we're ready to write an image volume to disk, we use `write_vol`

```{r eval=FALSE}
    write_vol(vol2, "output.nii")
    
    ## adding a '.gz' extension results ina gzipped file.
    write_vol(vol2, "output.nii.gz")
```


You can also write to a temporary file during workflows:

```{r}
    tmp <- tempfile(fileext = ".nii.gz")
    write_vol(vol2, tmp)
    file.exists(tmp)
    unlink(tmp)
```

For reorientation, resampling, and downsampling, use
`vignette("Resampling")`, which now owns that topic directly.