R/current_source_density.R

In craddm/eegUtils: Utilities for Electroencephalographic (EEG) Analysis

#' Channel interpolation
#'
#' Interpolate EEG channels using a spherical spline (Perrin et al., 1989;
#' 1990). The data must have channel locations attached.
#'
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#'
#' @param data Data as a `eeg_data` or `eeg_epochs` object.
#' @param bad_elecs Name(s) of electrode(s) to interpolate.
#' @param ... Other parameters passed to the functions.
#' @references * Perrin, F., Pernier, J., Bertrand, O., & Echallier, J. F.
#'   (1989). Spherical splines for scalp potential and current density mapping.
#'   Electroencephalography and Clinical Neurophysiology, 72, 184-187
#'  * Perrin, F., Pernier, J., Bertrand, O., & Echallier, J. F.
#'   (1990). Corrigenda EEG 02274. Electroencephalography and Clinical
#'   Neurophysiology, 76, 565
#' @export

interp_elecs <- function(data, bad_elecs, ...) {
  UseMethod("interp_elecs", data)
}

#' @export
interp_elecs.default <- function(data, bad_elecs, ...) {
  stop("Not implemented for objects of class", class(data))
}

#' @describeIn interp_elecs Interpolate EEG channel(s)
#' @export
interp_elecs.eeg_data <- function(data,
                                  bad_elecs,
                                  ...) {

  if (is.null(data$chan_info)) {
    stop("No channel locations found.")
  }

  bads_check <- bad_elecs %in% data$chan_info$electrode

  if (!all(bads_check)) {
    stop("No specified channels found.")
  }

  if (any(!bads_check)) {
    warning("Electrode(s) not found: ",
            paste0(bad_elecs[!bads_check],
                   collapse = " "))
    bad_elecs <- bad_elecs[bads_check]
  }

  data$chan_info <- validate_channels(data$chan_info,
                                      channel_names(data))
  missing_coords <- apply(is.na(data$chan_info),
                          1,
                          any)
  missing_chans <- names(data$signals)[missing_coords]

  if (any(missing_coords)) {
    message("Coords missing for electrodes ",
            paste0(missing_chans,
                   collapse = " "))
  }

  chan_info <- data$chan_info[!missing_coords, ]

  # Note - this checks for old-style electrode locations first, then for new
  # style
  check_ci_str(chan_info)

  if (all(c("cart_x", "cart_y", "cart_z") %in% names(chan_info))) {
    xyz_coords <- chan_info[, c("cart_x", "cart_y", "cart_z")]
    #normalise to unit sphere
    xyz_coords <- norm_sphere(xyz_coords)
  } else {
    xyz_coords <- sph_to_cart(chan_info$theta,
                              chan_info$phi,
                              1)
  }

  bad_select <- toupper(chan_info$electrode) %in% toupper(bad_elecs)
  xyz_bad <- xyz_coords[bad_select, ]
  xyz_good <- xyz_coords[!bad_select, ]

  if (nrow(xyz_bad) == 0) {
    return(data)
  }

  sigs_select <- toupper(names(data$signals)) %in%
    toupper(chan_info$electrode)

  bad_cols <- toupper(names(data$signals)) %in% toupper(bad_elecs)
  final_cols <- sigs_select & !bad_cols

  weights <- spheric_spline(xyz_good,
                            xyz_coords)

  data$signals <- interp_chans(data$signals,
                               bad_elecs,
                               missing_coords = missing_coords,
                               weights)
  data
}


#' Calculate a spherical spline smooth for interpolation of electrodes
#'
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#'
#' @param good_elecs Electrode positions that do not need interpolation
#' @param all_elecs Electrode positions including those to be interpolated
#' @keywords internal

spheric_spline <- function(good_elecs,
                           all_elecs) {

  lec <- compute_g(good_elecs,
                   good_elecs)

  diag(lec) <- diag(lec) + 1e-05

  ph <- compute_g(all_elecs,
                  good_elecs)

  g_dims <- dim(lec)
  lec <- cbind(lec, 1)
  lec <- rbind(lec, 1)
  lec[g_dims[1] + 1, g_dims[2] + 1] <- 0
  invC <- MASS::ginv(lec, tol = 0)
  W <- cbind(ph, 1)
  W <- W %*% invC[, 1:g_dims[1]]
  W
}

#' Interpolate channels
#' @param .data Channel data containing all data
#' @param bad_chans Vector of names of bad channels
#' @param missing_coords Logical vector indicating any channels in the data that
#'   had no associated coordinates
#' @param weights Spherical spline weights for interpolation
#' @keywords internal
interp_chans <- function(.data,
                         bad_chans,
                         missing_coords = FALSE,
                         weights) {

  bad_cols <- (toupper(names(.data)) %in% toupper(bad_chans)) | missing_coords
  weight_rows <- names(.data[, !missing_coords]) %in% bad_chans
  new_chans <- weights[weight_rows, , drop = TRUE] %*% t(.data[, !bad_cols])
  .data[, bad_chans] <- t(new_chans)
  .data
}

#' Convert to Current Source Density
#'
#' Convert an `eeg_data` or `eeg_epochs` object to using Current
#' Source Densities. This command uses a spherical spline algorithm (Perrin et
#' al., 1989) to compute scalp surface Laplacian/current source density
#' estimates of scalp potentials, a reference-free measure of electrical
#' activity that emphasises more local spatial features
#'
#' @param data `eeg_data` or `eeg_epochs` object
#' @param ... Other parameters
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#' @references * Perrin, F., Pernier, J., Bertrand, O., Echallier, J.F.
#'   (1989). Spherical splines for scalp potential and current density mapping.
#'   Electroencephalography and Clinical Neurophysiology, 72(2), 184-187. PMID:
#'   2464490
#'  * Kayser, J., Tenke, C.E. (2006). Principal components analysis
#'   of Laplacian waveforms as a generic method for identifying ERP generator
#'   patterns: I. Evaluation with auditory oddball tasks. Clinical
#'   Neurophysiology, 117(2), 348-368.
#'  * Kayser, J., Tenke, C.E. (2015).
#'   Issues and considerations for using the scalp surface Laplacian in EEG/ERP
#'   research: A tutorial review. International Journal of Psycholphysiology,
#'   97(3), 189-209
#' @examples
#' csd_epochs <- compute_csd(demo_epochs)
#' plot_butterfly(csd_epochs)
#'
#' # Compare the topographies of the CSD vs average referenced data
#' topoplot(demo_epochs, c(.2, .21))
#' topoplot(csd_epochs, c(.2, .21))
#' @export
compute_csd <- function(data,
                        ...) {
  UseMethod("compute_csd", data)
}

#' @describeIn compute_csd Default method to detect unknown classes.
#' @export
compute_csd.default <- function(data,
                                ...) {
  warning("compute_csd not implemented for objects of class", class(data))
}

#' @param m smoothing constraint (higher = more rigid splines)
#' @param smoothing lambda constant. Added to the Defaults to 1e-05
#' @param scaling Default scaling (1) is uV / m^2. Note that this depends on the
#'   units of the electrode co-ordinates.
#' @describeIn compute_csd Transform `eeg_data` to CSD
#' @export
compute_csd.eeg_data <- function(data,
                                 m = 4,
                                 smoothing = 1e-05,
                                 scaling = 1,
                                 ...) {
  convert_to_csd(data,
                 m,
                 smoothing,
                 scaling)
}

#' @describeIn compute_csd Transform `eeg_data` to CSD
#' @export
compute_csd.eeg_epochs <- function(data,
                                   m = 4,
                                   smoothing = 1e-05,
                                   scaling = 1,
                                   ...) {
  convert_to_csd(data,
                 m,
                 smoothing,
                 scaling)
}

#' Calculate current source densities
#'
#' @param data `eeg_data` object
#' @param m smoothing constraint (higher = more rigid)
#' @param smoothing lambda constant
#' @param scaling Default scaling (1) is uV / m^2.
#' @keywords internal

convert_to_csd <- function(data,
                           m = 4,
                           smoothing = 1e-05,
                           scaling = 1) {

  orig_elecs <- names(data$signals)
  data_chans <- toupper(names(data$signals)) %in% toupper(data$chan_info$electrode)
  scaling <- scaling * scaling

  if (any(!data_chans)) {
    stop("No channel information found for ",
         paste0(names(data$signals)[!data_chans],
                collapse = " "),
         ". Either remove channel or add channel info.")
  }

  missing_coords <- apply(data$chan_info,
                          1,
                          function(x) any(is.na(x)))

  if (any(missing_coords)) {
    stop("No coordinates for ",
      paste0(data$chan_info$electrode[missing_coords],
             collapse = " "),
      ". Remove channels before applying CSD."
    )
  }

  # Convert data to average reference
  data <- eeg_reference(data,
                        verbose = FALSE)

  if (all(c("cart_x", "cart_y", "cart_z") %in% names(data$chan_info))) {
    xyz_coords <- data$chan_info[, c("cart_x", "cart_y", "cart_z")]
    #normalise to unit sphere
    xyz_coords <- norm_sphere(xyz_coords)
  } else {
    xyz_coords <- sph_to_cart(data$chan_info$theta,
                              data$chan_info$phi,
                              1)
  }

# Optimize matrix operations
  g_matrix <- compute_g(xyz_coords, xyz_coords, m = m, iter = 50)
  h_matrix <- compute_h(xyz_coords, xyz_coords, m = m, iter = 50)

  diag(g_matrix) <- diag(g_matrix) + smoothing
  g_inverse <- solve(g_matrix)
  g_inverse_col_sums <- colSums(g_inverse)
  g_inverse_total_sum <- sum(g_inverse_col_sums)

  signal_matrix <- as.matrix(data$signals[, data_chans])

  # Vectorize operations
  intermediate_result <- signal_matrix %*% t(g_inverse)
  row_averages <- rowSums(intermediate_result) / g_inverse_total_sum
  centered_result <- intermediate_result - outer(row_averages, g_inverse_col_sums)
  csd_values <- (centered_result %*% h_matrix) / scaling

  data$signals[, data_chans] <- as.data.frame(csd_values)
  names(data$signals) <- orig_elecs
  data$reference$ref_chans <- "CSD"
  data
}

#' Compute the g function for two sets of channel locations on the
#' unit sphere.
#'
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#'
#' @param xyz_coords A set of electrode locations on a unit sphere.
#' @param xyz_elecs A set of electrode locations on a unit sphere.
#' @param m Interpolation constant (higher = less flexible)
#' @param iter Number of iterations for calculations
#' @keywords internal

compute_g <- function(xyz_coords,
                      xyz_elecs,
                      m = 4,
                      iter = 7) {

  # Convert inputs to matrices
  xyz_coords <- as.matrix(xyz_coords)
  xyz_elecs <- as.matrix(xyz_elecs)

  # Calculate cosine similarity
  EI <- xyz_coords %*% t(xyz_elecs)

  # Pre-allocate matrix
  g <- matrix(0, ncol = ncol(EI), nrow = nrow(EI))

  # Pre-compute coefficients
  coeff <- (2 * seq_len(iter) + 1) / (seq_len(iter)^m * (seq_len(iter) + 1)^m)

  # Compute all Legendre polynomials up to degree 'iter'
  P <- legendre_polynomials(iter, EI)

  # Sum the polynomials
  for (i in seq_len(iter)) {
    g <- g + coeff[i] * P[i + 1, ]  # Use the i-th degree term
  }

  # Final scaling
  g / (4 * pi)
}

#' Compute the h function for two sets of channel locations on the
#' unit sphere.
#'
#' @author Matt Craddock \email{matt@@mattcraddock.com}
#'
#' @param xyz_coords A set of electrode locations on a unit sphere.
#' @param xyz_elecs A set of electrode locations on a unit sphere.
#' @param m Interpolation constant (higher = less flexible)
#' @param iter iterations for calculations
#' @keywords internal

compute_h <- function(xyz_coords,
                      xyz_elecs,
                      m = 4,
                      iter = 50) {

  # Convert inputs to matrices
  xyz_coords <- as.matrix(xyz_coords)
  xyz_elecs <- as.matrix(xyz_elecs)

  # Calculate cosine similarity
  EI <- xyz_coords %*% t(xyz_elecs)

  # Pre-allocate matrix for results
  h <- matrix(0, ncol = ncol(EI), nrow = nrow(EI))

  # Pre-compute coefficients
  coeff <- (-2 * seq_len(iter) - 1) / (seq_len(iter)^(m - 1) * (seq_len(iter) + 1)^(m - 1))

  # Compute all Legendre polynomials up to degree 'iter'
  P <- legendre_polynomials(iter, EI)

  # Sum the polynomials
  for (i in seq_len(iter)) {
    h <- h + coeff[i] * P[i + 1, ]  # Use the i-th degree term
  }

  # Final scaling
  -h / (4 * pi)
}

#' Calculate Legendre polynomials up to degree n
#'
#' @param n Maximum degree of Legendre polynomials
#' @param x Input values
#' @return Matrix of Legendre polynomial values
#' @keywords internal
legendre_polynomials <- function(n, x) {
  P <- matrix(0, nrow = n + 1, ncol = length(x))
  P[1, ] <- 1
  if (n > 0) P[2, ] <- x

  for (k in 2:n) {
    P[k + 1, ] <- ((2 * k - 1) * x * P[k, ] - (k - 1) * P[k - 1, ]) / k
  }

  P
}