---
title: Displaying Surfaces with RGL
author: Bradley Buchsbaum
date: '`r Sys.Date()`'
output:
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 5
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
vignette: |
  %\VignetteIndexEntry{Displaying Surfaces with RGL}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
params:
  family: teal
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html

---

This vignette demonstrates how to display 3D brain surface meshes using the `rgl` plotting tools provided by the `neurosurf` package, primarily through the `plot()` method which utilizes the `view_surface()` function internally.

For interactive HTML widgets, see `vignette("interactive-surfaces")`.
For high-level, multi-view layouts with colourbars and atlas outlines, see `vignette("surface-figures")`.

## Setup and Loading Data

First, we set up `knitr` options to embed `rgl` plots directly into the HTML output using WebGL and prevent standalone `rgl` windows from popping up during knitting. We then load example left and right hemisphere white matter surfaces included with the package and prepare some data (smoothed geometry, curvature, random values) for the examples.

```{r setup, include = FALSE}
# Check for required packages - skip vignette build if rgl unavailable
has_rgl <- requireNamespace("rgl", quietly = TRUE)
if (!has_rgl) {
  knitr::opts_chunk$set(eval = FALSE)
  message("Skipping vignette: rgl package not available")
}

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 = "#>",
  warning = FALSE,
  message = FALSE,
  fig.width = 7,        # Default figure width
  fig.height = 5        # Default figure height
)

if (has_rgl) {
  # Use the NULL device to suppress external windows and embed widgets
  options(rgl.useNULL = TRUE, rgl.printRglwidget = TRUE)

  # Let rgl autoprint widgets into the document (>= 0.102.21)
  rgl::setupKnitr(autoprint = TRUE)

  # Additional knitr defaults for figure sizing in vignettes
  knitr::opts_chunk$set(
    dev = 'png',
    out.width = "100%",
    out.height = "400px"
  )
  options(rgl.max.lights = 2)
  library(rgl)
  options(rgl.useNULL = TRUE)
}
library(neuroim2)
library(neurosurf)

# Counter for unique figure names
.render_counter <- new.env(parent = emptyenv())
.render_counter$n <- 0

# Treat uniform black/white PNGs as failed snapshots so the vignette can fall
# back to interactive widgets instead of embedding bad files.
snapshot_is_usable <- function(path,
                               uniform_fraction = 0.995,
                               near_black = 0.02,
                               near_white = 0.98,
                               low_variation_sd = 0.005,
                               low_variation_span = 0.025) {
  if (!is.character(path) || length(path) != 1L || !nzchar(path) ||
      !file.exists(path)) {
    return(FALSE)
  }

  img <- try(png::readPNG(path), silent = TRUE)
  if (inherits(img, "try-error")) {
    return(TRUE)
  }

  rgb <- if (length(dim(img)) >= 3L) {
    img[, , seq_len(min(3L, dim(img)[3L])), drop = FALSE]
  } else {
    img
  }

  vals <- as.numeric(rgb)
  vals <- vals[is.finite(vals)]
  if (!length(vals)) {
    return(FALSE)
  }

  if (mean(vals <= near_black) >= uniform_fraction ||
      mean(vals >= near_white) >= uniform_fraction) {
    return(FALSE)
  }

  !(stats::sd(vals) < low_variation_sd &&
      diff(range(vals)) < low_variation_span)
}

snapshot_current_scene <- function(file) {
  if (rgl::rgl.useNULL() && !requireNamespace("webshot2", quietly = TRUE)) {
    return(character())
  }

  if (rgl::rgl.useNULL()) {
    rgl::snapshot3d(file, webshot = TRUE)
  } else {
    rgl::rgl.snapshot(file)
  }

  if (snapshot_is_usable(file)) file else character()
}

# Helper function: try snapshot_surface first, fall back to rglwidget
render_surface <- function(surfgeom, ..., width = 800, height = 600) {
  # Use knitr's figure path to ensure proper asset tracking
  .render_counter$n <- .render_counter$n + 1
  fig_file <- knitr::fig_path(paste0("-render-", .render_counter$n, ".png"))
  dir.create(dirname(fig_file), recursive = TRUE, showWarnings = FALSE)

  img_path <- try(snapshot_surface(surfgeom, file = fig_file, width = width, height = height, ...), silent = TRUE)
  if (!inherits(img_path, "try-error") && snapshot_is_usable(img_path)) {
    knitr::include_graphics(img_path)
  } else {
    # Fallback to rglwidget
    rgl::open3d()
    view_surface(surfgeom, new_window = FALSE, ...)
    widget <- rgl::rglwidget()
    rgl::close3d()
    widget
  }
}

# Helper function: render multiple viewpoints in a grid
render_multi_view <- function(surfgeom, viewpoints = c("lateral", "medial", "ventral", "posterior"),
                              ..., width = 800, height = 600) {
  paths <- lapply(seq_along(viewpoints), function(i) {
    vp <- viewpoints[i]
    .render_counter$n <- .render_counter$n + 1
    fig_file <- knitr::fig_path(paste0("-multiview-", .render_counter$n, ".png"))
    dir.create(dirname(fig_file), recursive = TRUE, showWarnings = FALSE)
    try(snapshot_surface(surfgeom, file = fig_file, viewpoint = vp, width = width, height = height, ...), silent = TRUE)
  })

  valid <- sapply(paths, function(p) !inherits(p, "try-error") && snapshot_is_usable(p))

  if (all(valid)) {
    knitr::include_graphics(unlist(paths))
  } else {
    # Fallback to rglwidget with mfrow3d
    n <- length(viewpoints)
    ncol <- min(n, 2)
    nrow <- ceiling(n / ncol)
    rgl::open3d()
    rgl::mfrow3d(nrow, ncol, sharedMouse = TRUE)
    for (i in seq_along(viewpoints)) {
      if (i > 1) rgl::next3d()
      view_surface(surfgeom, viewpoint = viewpoints[i], new_window = FALSE, ...)
    }
    widget <- rgl::rglwidget()
    rgl::close3d()
    widget
  }
}

# Helper function: render the left and right hemispheres as a single figure.
# Each hemisphere can take its own background colors (e.g. per-hemisphere
# curvature). When static snapshots succeed we emit two side-by-side images;
# otherwise we fall back to ONE combined rglwidget (a single mfrow3d scene).
# Emitting two independent rglwidgets does not lay out reliably in HTML and
# can leave only the first hemisphere visible, so we always collapse the
# fallback into a single widget.
render_hemi_pair <- function(lh, rh, bgcol_lh, bgcol_rh, ...,
                             viewpoint = "lateral", width = 900, height = 700) {
  hemis <- list(list(geom = lh, bg = bgcol_lh),
                list(geom = rh, bg = bgcol_rh))

  paths <- lapply(seq_along(hemis), function(i) {
    .render_counter$n <- .render_counter$n + 1
    fig_file <- knitr::fig_path(paste0("-hemipair-", .render_counter$n, ".png"))
    dir.create(dirname(fig_file), recursive = TRUE, showWarnings = FALSE)
    try(snapshot_surface(hemis[[i]]$geom, file = fig_file, bgcol = hemis[[i]]$bg,
                         viewpoint = viewpoint, width = width, height = height, ...),
        silent = TRUE)
  })

  valid <- sapply(paths, function(p) !inherits(p, "try-error") && snapshot_is_usable(p))

  if (all(valid)) {
    knitr::include_graphics(unlist(paths))
  } else {
    # Fallback: both hemispheres in a single mfrow3d widget.
    rgl::open3d()
    rgl::mfrow3d(1, 2, sharedMouse = TRUE)
    view_surface(lh, bgcol = bgcol_lh, viewpoint = viewpoint, new_window = FALSE, ...)
    rgl::next3d()
    view_surface(rh, bgcol = bgcol_rh, viewpoint = viewpoint, new_window = FALSE, ...)
    widget <- rgl::rglwidget()
    rgl::close3d()
    widget
  }
}

# Load example surfaces
white_lh_asc <- system.file("extdata", "std.8_lh.smoothwm.asc", package="neurosurf")
white_rh_asc <- system.file("extdata", "std.8_rh.smoothwm.asc", package="neurosurf")
white_lh <- read_surf(white_lh_asc)
white_rh <- read_surf(white_rh_asc)

# Prepare data for examples
white_lh_base <- smooth(white_lh, type = "HCLaplace", delta = 0.2, iteration = 5)
white_rh_base <- smooth(white_rh, type = "HCLaplace", delta = 0.2, iteration = 5)
set.seed(123) # for reproducibility
random_vals <- rnorm(length(nodes(white_lh_base)))

# Keep coarse meshes for algorithms; use display copies we may refine
white_lh_display <- white_lh_base
white_rh_display <- white_rh_base
curv_lh_display <- curvature(white_lh_display)
curv_rh_display <- curvature(white_rh_display)
random_vals_display <- random_vals
random_vals_display_smooth <- random_vals_display

```

```{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-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

```{r upsample-hidden, echo=FALSE}
# Gently refine the std8 mesh to reduce faceting in the vignette visuals.
# Apply the same display pipeline to both hemispheres so side-by-side panels
# are comparable.
refine_display_surface <- function(surfgeom) {
  surf_display <- surfgeom

  if (requireNamespace("Rvcg", quietly = TRUE)) {
    try({
      if (is.function(Rvcg::vcgSubdivide)) {
        surf_display@mesh <- Rvcg::vcgSubdivide(
          surf_display@mesh,
          type = "loop",
          iterations = 2,
          silent = TRUE
        )
      }
      surf_display <- smooth(
        surf_display,
        type = "HCLaplace",
        delta = 0.15,
        iteration = 8
      )
    }, silent = TRUE)
  }

  surf_display
}

smooth_surface_values <- function(surfgeom, vals, passes = 1) {
  tryCatch({
    surf_vals <- NeuroSurface(
      surfgeom,
      indices = seq_len(length(vals)),
      data = vals
    )
    for (i in seq_len(passes)) {
      surf_vals <- smooth(surf_vals)
    }

    out <- values(surf_vals)
    if (length(out) != length(vals) || any(!is.finite(out))) {
      vals
    } else {
      out
    }
  }, error = function(e) vals)
}

white_lh_display <- refine_display_surface(white_lh_base)
white_rh_display <- refine_display_surface(white_rh_base)

# Final safety: ensure data vectors match the current display vertex count
nv_display <- length(nodes(white_lh_display))
if (length(random_vals_display) != nv_display) {
  random_vals_display <- rnorm(nv_display)
}

# Smooth the overlay data and curvature on the display mesh to avoid speckles
# in the static vignette snapshots.
random_vals_display_smooth <- smooth_surface_values(
  white_lh_display,
  random_vals_display,
  passes = 1
)
curv_lh_display <- smooth_surface_values(
  white_lh_display,
  curvature(white_lh_display),
  passes = 5
)
curv_rh_display <- smooth_surface_values(
  white_rh_display,
  curvature(white_rh_display),
  passes = 5
)
```

## Basic Surface Plotting

The simplest way to display a `SurfaceGeometry` object is using the `plot()` method. By default, it renders the surface with a light gray background. We can specify a `viewpoint`.



```{r basic-plot, results='asis'}
# Plot the smoothed left hemisphere from a lateral viewpoint
render_surface(white_lh_display, viewpoint = "lateral", lit = TRUE)
```

## Coloring Based on Curvature

Surface curvature helps distinguish gyri (outward folds) from sulci (inward folds). The `curvature()` function calculates this, and `curv_cols_smooth()` maps the values to a continuous grayscale gradient (dark in sulci, light on gyri) for natural-looking shading. For a simpler binary split, see `curv_cols()`. Either way, pass the resulting colors to the `bgcol` argument of `plot()`.

```{r curvature-plot, results='asis'}
# Calculate curvature colors
curv_colors <- curv_cols_smooth(curv_lh_display)

# Plot with curvature background from a medial viewpoint
render_surface(white_lh_display, bgcol = curv_colors, viewpoint = "medial", specular = "black")
```

## Overlaying Data Values

Often, we want to visualize data mapped onto the surface vertices (e.g., activation values, thickness). We can pass a vector of values to the `vals` argument. The `cmap` argument specifies the color map, and `irange` defines the data range to map onto the colormap. Values outside `irange` are clamped to the minimum or maximum color.

```{r data-overlay, results='asis'}
# Overlay random data using a rainbow colormap
# Map data range from -2 to 2 onto the colormap
render_surface(white_lh_display, vals = random_vals_display_smooth, cmap = rainbow(256),
               irange = c(-2, 2), thresh = NULL, viewpoint = "lateral", specular = "gray")
```

## Thresholding Data Visualization

The `thresh` argument (a vector of two values, `c(lower, upper)`) can be used with `vals` to make parts of the surface transparent. Vertices where the corresponding value in `vals` is *inside* this range (between `lower` and `upper`) are rendered transparently; values outside remain opaque. This is useful for masking out a band of values.

```{r threshold-plot, results='asis'}
# Same data overlay as above, but make values between -1 and 1 transparent
render_surface(white_lh_display, vals = random_vals_display_smooth, cmap = rainbow(256),
               irange = c(-2, 2), thresh = c(-1, 1), viewpoint = "lateral", lit = TRUE)
```

## Direct Vertex Coloring

Instead of mapping data values to a colormap, you can provide a vector of specific hex color codes directly to the `vert_clrs` argument. This overrides `vals` and `cmap`. The vector length must match the number of vertices.

```{r vertex-color-plot, results='asis'}
# Color vertices based on their x-coordinate (e.g., red for positive x, blue for negative)
x_coords <- coords(white_lh_display)[, 1]
vertex_colors <- ifelse(x_coords > median(x_coords), "#FF0000", "#0000FF") # Red/Blue

render_surface(white_lh_display, vert_clrs = vertex_colors, viewpoint = "ventral", lit = TRUE)
```

## Controlling Transparency

The `alpha` argument controls the overall transparency of the surface, ranging from 0 (fully transparent) to 1 (fully opaque).

```{r alpha-plot, results='asis'}
# Plot the surface with 60% opacity (40% transparent)
render_surface(white_lh_display, vals = random_vals_display_smooth, cmap = heat.colors(256),
               irange = c(-2, 2), alpha = 0.6, viewpoint = "posterior")
```

## Adjusting Lighting and Material

The appearance of the surface is affected by lighting. The `specular` argument controls the color of specular highlights (shininess). Setting it to `"black"` creates a matte appearance.

```{r lighting-plot, results='asis'}
# Plot with a matte finish (no specular highlights)
render_surface(white_lh_display, vals = random_vals_display_smooth, cmap = topo.colors(256),
               irange = c(-2, 2), specular = "black", viewpoint = "lateral", lit = TRUE)
```

## Snapshotting to an image (for knitr/CI)

Use `snapshot_surface()` to render an off-screen PNG and include it directly:

```{r snapshot-example}
.render_counter$n <- .render_counter$n + 1
snapshot_file <- knitr::fig_path(paste0("-snapshot-", .render_counter$n, ".png"))
dir.create(dirname(snapshot_file), recursive = TRUE, showWarnings = FALSE)

img_path <- try(snapshot_surface(white_lh_display,
                                 file = snapshot_file,
                                 vals = random_vals_display_smooth,
                                 cmap = viridis::viridis(256),
                                 viewpoint = "lateral",
                                 specular = "black",
                                 width = 1200, height = 900),
                silent = TRUE)

if (!inherits(img_path, "try-error") && snapshot_is_usable(img_path)) {
  knitr::include_graphics(img_path)
} else {
  rgl::open3d()
  view_surface(white_lh_display,
               vals = random_vals_display_smooth,
               cmap = viridis::viridis(256),
               viewpoint = "lateral",
               specular = "black",
               new_window = FALSE)
  widget <- rgl::rglwidget()
  rgl::close3d()
  widget
}
```

## Changing Viewpoints

The `viewpoint` argument can be set to common anatomical views like `"lateral"`, `"medial"`, `"ventral"`, or `"posterior"`. The function automatically selects the correct left/right version based on the surface's hemisphere information (`surf@hemi`).

```{r viewpoints-plot, results='asis', fig.height=8, out.width="48%", fig.show='hold'}
# Display multiple viewpoints with curvature shading
render_multi_view(white_lh_display,
                  viewpoints = c("lateral", "medial", "ventral", "posterior"),
                  bgcol = curv_cols_smooth(curv_lh_display), specular = "black")
```

## Displaying Two Hemispheres

For lateral views, each hemisphere is best rendered separately since the camera can only face one direction. We render the left and right lateral views side by side.

```{r two-hemispheres-plot, results='asis', out.width="49%", out.height=NULL, fig.show='hold'}
# Render both hemispheres as a single figure so they always appear together
# (two side-by-side images, or one combined widget when snapshots are
# unavailable). Leave some extra margin so the static PNGs do not feel cramped.
render_hemi_pair(
  white_lh_display,
  white_rh_display,
  bgcol_lh = curv_cols_smooth(curv_lh_display, quantiles = c(0.02, 0.98)),
  bgcol_rh = curv_cols_smooth(curv_rh_display, quantiles = c(0.02, 0.98)),
  viewpoint = "lateral",
  specular = "black",
  zoom = 0.92,
  width = 900,
  height = 700
)
```

## Adding Spheres to the Surface

The `spheres` argument allows you to draw spherical markers at specified coordinates. It requires a data frame with columns `x`, `y`, `z`, and `radius`. An optional `color` column can specify colors for each sphere.

```{r spheres-plot, results='asis'}
# Define coordinates for some spherical markers
# Sample some vertex indices safely from available vertices
n_vertices <- nrow(coords(white_lh_display))
sample_indices <- sample(1:n_vertices, size = min(3, n_vertices))

peak_coords <- data.frame(
  x = coords(white_lh_display)[sample_indices, 1],
  y = coords(white_lh_display)[sample_indices, 2],
  z = coords(white_lh_display)[sample_indices, 3],
  radius = c(3, 4, 2.5)[1:length(sample_indices)],
  color = c("yellow", "cyan", "magenta")[1:length(sample_indices)]
)

# Plot the surface with curvature shading and add the spheres
render_surface(white_lh_display, bgcol = curv_cols_smooth(curv_lh_display),
               viewpoint = "lateral", specular = "black", spheres = peak_coords)
```

## Plotting Other NeuroSurface Objects

The `plot()` method also works for other classes like `NeuroSurface`, `LabeledNeuroSurface`, and `ColorMappedNeuroSurface`. These objects already contain data and potentially color mapping information. The `plot` method extracts this information and passes the appropriate arguments (like `vals`, `cmap`, `irange`, `thresh`, `vert_clrs`) to the underlying `view_surface` function.

```{r neurosurface-plot, results='asis'}
# Create a NeuroSurface object with the random data
nsurf <- NeuroSurface(white_lh_display, indices = 1:length(random_vals_display), data = random_vals_display)

# Plot the NeuroSurface - uses data stored within the object
# We can still override or add parameters like cmap, irange, thresh, alpha etc.
render_surface(geometry(nsurf), vals = values(nsurf), cmap = heat.colors(128),
               irange = c(-2.5, 2.5), viewpoint = "lateral")
```


## Showing an activation map overlaid on a surface mesh

We will plot surface in a row of 3. We generate a set of random values and then smooth those values along the surface to approximate a realistic activation pattern.

In the first column we display all the values in the map. Next we make values between (-0.2, 0.2) transparent. In the last panel we additionally add a cluster size threshold of 30 nodes.

`surface_montage()` handles the per-panel rendering and layout for us: it captures each panel as a static image and tiles them into one figure (falling back to a single interactive widget when static snapshots are unavailable). Each panel is either a surface or a `list(surface, ...overrides)`, and arguments shared by every panel (here `cmap` and `irange`) are passed once.

```{r activation-map, results='asis', fig.width=9, fig.height=3.5, out.width="100%"}
vals <- rnorm(length(nodes(white_lh_base)))
ssurf <- smooth(NeuroSurface(white_lh_base, indices = seq_along(vals), data = vals))
csurf <- cluster_threshold(ssurf, size = 30, threshold = c(-0.2, 0.2))

surface_montage(
  list(
    ssurf,                                # all values
    list(ssurf, thresh = c(-0.2, 0.2)),   # band around zero made transparent
    list(csurf, thresh = c(-0.2, 0.2))    # + cluster-size threshold (>= 30 nodes)
  ),
  cmap = rainbow(100), irange = c(-2, 2), ncol = 3
)
```

## Showing two hemispheres in same scene

For views where the left-right axis maps to the screen (posterior, anterior, dorsal), both hemispheres can share a single scene since their coordinates naturally separate (LH at x < 0, RH at x > 0).

```{r two-hemi-posterior, results='asis', fig.width=8, fig.height=5}
# Two hemispheres shown from posterior viewpoint
.render_counter$n <- .render_counter$n + 1
posterior_file <- knitr::fig_path(paste0("-posterior-", .render_counter$n, ".png"))
dir.create(dirname(posterior_file), recursive = TRUE, showWarnings = FALSE)

img_path <- try({
  file <- posterior_file
  rgl::open3d()
  rgl::par3d(windowRect = c(0, 0, 1200, 600))
  rgl::bg3d(color = "white")

  # LH and RH sit naturally at x<0 and x>0; small offset adds breathing room
  view_surface(white_lh_display, bgcol = curv_cols_smooth(curv_lh_display),
               viewpoint = "posterior", new_window = FALSE, offset = c(-5, 0, 0))
  view_surface(white_rh_display, bgcol = curv_cols_smooth(curv_rh_display),
               viewpoint = "posterior", new_window = FALSE, offset = c(5, 0, 0))
  rgl::view3d(fov = 0, zoom = 0.55,
              userMatrix = rbind(c(1,0,0,0), c(0,0,1,0), c(0,-1,0,0), c(0,0,0,1)))
  snapshot_current_scene(file)
}, silent = TRUE)
try(rgl::close3d(), silent = TRUE)

if (!inherits(img_path, "try-error") && snapshot_is_usable(img_path)) {
  knitr::include_graphics(img_path)
} else {
  # Fallback to rglwidget
  rgl::open3d()
  view_surface(white_lh_display, bgcol = curv_cols_smooth(curv_lh_display),
               viewpoint = "posterior", new_window = FALSE, offset = c(-5, 0, 0))
  view_surface(white_rh_display, bgcol = curv_cols_smooth(curv_rh_display),
               viewpoint = "posterior", new_window = FALSE, offset = c(5, 0, 0))
  rgl::view3d(fov = 0, zoom = 0.55,
              userMatrix = rbind(c(1,0,0,0), c(0,0,1,0), c(0,-1,0,0), c(0,0,0,1)))
  rgl::rglwidget()
}
```


## Next Steps

For **interactive 3D visualization** with `surfwidget()`, see `vignette("interactive-surfaces")`.

For **publication-quality multi-view figures**, see `vignette("surface-figures")`.
