R/surface_dist.R

In freesurferformats: Read and Write 'FreeSurfer' Neuroimaging File Formats

# Functions for computing Euclidean distance between vertex coordinates.
# These functions are not really specific for meshes as they simple work with vertex coords (=points) without
# any knowledge or consideration of mesh topology (graph distance) or the surface (geodesic distance).

#' @title Find vertex index closest to given query coordinate using Euclidean distance.
#'
#' @param surface an fs.surface instance or a nx3 numerical matrix representing mesh points.
#'
#' @param point_coords nx3 matrix of query coords. If a vector, will be transformed \code{byrow} to such a matrix.
#'
#' @return named list with entries: 'vertex_id' integer vector, the index of the closest vertex, and 'dist': double vector, the Euclidean distance to that vertex.
#'
#' @family Euclidean distance util functions
#'
#' @export
#' @importFrom stats dist
closest.vert.to.point <- function(surface, point_coords) {

  # Check params and extract vertex coords if needed.
  if(is.fs.surface(surface)) {
    vertices = surface$vertices;
  } else {
    if(is.matrix(surface)) {
      if(ncol(surface) == 3L) {
        vertices = surface;
      } else {
        stop("Matrix in parameter 'surface' must have 3 columns.");
      }
    } else {
      stop("Parameter 'surface' must be a numerical matrix or an fs.surface instance.");
    }
  }

  if(is.vector(point_coords)) {
    point_coords = matrix(point_coords, ncol = 3L, byrow = TRUE);
  }
  if(ncol(point_coords) != 3L) {
    stop("Parameter coords must have 3 columns.");
  }

  num_query_points = nrow(point_coords);
  if(num_query_points < 1L) {
    stop("Parameter 'point_coords' must not be empty.");
  }

  results = list('vertex_id'=rep(NA, num_query_points), 'dist'=rep(NA, num_query_points));

  for(point_idx in seq(num_query_points)) {
    point = point_coords[point_idx, ];
    dists = apply(vertices, 1, euclidian.dist, point);
    results$vertex_id[point_idx] = which.min(dists);
    results$dist[point_idx] = min(dists);
  }

  return(results);
}


#' @title Compute Euclidean distance from all mesh vertices to given point.
#'
#' @param surface an fs.surface instance
#'
#' @param point_coords double vector of length 3, the xyz coords of a single point.
#'
#' @family Euclidean distance util functions
#'
#' @return double vector of distances
#'
#' @export
vertexdists.to.point <- function(surface, point_coords) {
  assert.surface(surface);
  return(apply(surface$vertices, 1, euclidian.dist, point_coords));
}


#' @title Compute Euclidean distance.
#'
#' @param x1 numerical vector, coords of first point
#'
#' @param x2 numerical vector, coords of second point
#'
#' @return the Euclidean distance between x1 and x2.
#'
#' @keywords internal
euclidian.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))

#' @title Compute Euclidean distance between two vertices v1 and v2.
#'
#' @param surface an fs.surface instance
#'
#' @param v1 positive integer, vertex index of v1
#'
#' @param v2 positive integer, vertex index of v2
#'
#' @return double, the Euclidean distance between v1 and v2
#'
#' @family Euclidean distance util functions
#'
#' @export
vertex.euclid.dist <- function(surface, v1, v2) {
  assert.surface(surface);
  euclidian.dist(surface$vertices[v1, ], surface$vertices[v2, ])
}