## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 5
)

## ----setup--------------------------------------------------------------------
library(neuromapr)

## ----mesh---------------------------------------------------------------------
vertices <- matrix(c(
  0, 0, 0,
  1, 0, 0,
  0.5, 1, 0,
  1.5, 1, 0,
  0, 2, 0,
  1, 2, 0
), ncol = 3, byrow = TRUE)

faces <- matrix(c(
  1, 2, 3,
  2, 3, 4,
  3, 4, 6,
  3, 5, 6
), ncol = 3, byrow = TRUE)

## ----graph--------------------------------------------------------------------
g <- make_surf_graph(vertices, faces)
g

## ----igraph-plot, fig.cap = "Surface graph of the toy mesh. Vertices are positioned at their x-y coordinates, edges are weighted by Euclidean distance."----
layout_xy <- vertices[, 1:2]
edge_weights <- round(
  igraph::E(g)$weight, 2
)

plot(
  g,
  layout = layout_xy,
  vertex.size = 20,
  vertex.color = "steelblue",
  vertex.label.color = "white",
  vertex.label.cex = 0.9,
  edge.label = edge_weights,
  edge.label.cex = 0.8,
  edge.color = "grey40",
  main = "Surface mesh graph"
)

## ----geodesic-full------------------------------------------------------------
dmat <- get_surface_distance(vertices, faces)
dmat

## ----geodesic-partial---------------------------------------------------------
dmat_partial <- get_surface_distance(
  vertices, faces, source_vertices = c(1, 6)
)
dmat_partial

## ----compare-distances--------------------------------------------------------
euclid <- as.matrix(dist(vertices))
geodesic <- get_surface_distance(vertices, faces)

max(abs(geodesic - euclid))

## ----vertex-areas-------------------------------------------------------------
areas <- vertex_areas(vertices, faces)
areas

## ----total-area---------------------------------------------------------------
sum(areas)

## ----triangle-check-----------------------------------------------------------
tri_area <- function(v, f) {
  a <- v[f[1], ]
  b <- v[f[2], ]
  cc <- v[f[3], ]
  ab <- b - a
  ac <- cc - a
  cross <- c(
    ab[2] * ac[3] - ab[3] * ac[2],
    ab[3] * ac[1] - ab[1] * ac[3],
    ab[1] * ac[2] - ab[2] * ac[1]
  )
  0.5 * sqrt(sum(cross^2))
}

total_tri <- sum(vapply(
  seq_len(nrow(faces)),
  function(i) tri_area(vertices, faces[i, ]),
  numeric(1)
))

all.equal(sum(areas), total_tri)