#
# <<< channel positions >>> --------
#
#' Transformation between spherical, cartesian and geographical coordinates
#'
#' This package contains several functions which require a data.frame of the
#' exact channel positions (e.g., \code{\link{plot2dview}},
#' \code{\link{chanNb}}). To facilitate the transformation between various
#' coordinate systems, \code{sph2cart}, \code{sph2geo}, \code{cart2sph},
#' \code{cart2geo}, \code{geo2sph}, \code{geo2cart} transform spherical (sph),
#' cartesian (cart), or geographical (geo) coordinates into the respective
#' coordinate system.
#' @name coordinates
#' @param ch_pos a data frame or matrix containing the spherical, cartesian,
#' or geographical coordinates of the electrode positions. It should contain
#' at least the following (named) columns, respectively: theta and phi;
#' x, y, and z; or lat and long. The rownames attribute of \code{ch_pos} is
#' interpreted as the names of the electrodes.
#' @param r radius (default = 1)
#' @param deg logical variable indicating whether spherical or geographical
#' coordinates are or should be given in degrees (TRUE, default)
#' @param long360 logical variable; if TRUE (default), longitudes range from
#' 0 to 360 (otherwise from -180 to 180)
#' @param orient a character value specifying the orientation in the
#' geographical coordinate system; "northpole" (default) or "equatorial"
#' @return A data.frame with converted coordinates
NULL
#' @rdname coordinates
#' @export
# Transform cartesian to spherical coordinates
sph2cart <- function(ch_pos, r = 1, deg = TRUE) {
ch_pos <- as.data.frame(ch_pos)
if (!all(c("theta", "phi") %in% colnames(ch_pos)))
stop("Either theta or phi angles are missing!")
if (!is.null(ch_pos$r)) r <- ch_pos$r
if (deg)
ch_pos[, c("theta", "phi")] <- ch_pos[, c("theta", "phi")]/180*pi
theta <- ch_pos$theta
phi <- ch_pos$phi
out <- data.frame(
x = r * sin(theta) * cos(phi),
y = r * sin(theta) * sin(phi),
z = r * cos(theta)
)
rownames(out) <- rownames(ch_pos)
# return
out
}
#' @rdname coordinates
#' @export
# Transform cartesian to spherical coordinates
cart2sph <- function(ch_pos, deg = TRUE) {
ch_pos <- as.data.frame(ch_pos)
if (!all(c("x", "y", "z") %in% colnames(ch_pos)))
stop("Either x, y, or z coordinates are missing!")
x <- ch_pos$x
y <- ch_pos$y
z <- ch_pos$z
r <- sqrt(x^2 + y^2 + z^2)
theta <- acos(z/r)
phi <- atan2(y, x)
out <- data.frame(theta, phi, r)
rownames(out) <- rownames(ch_pos)
if (deg) out[, 1:2] <- out[, 1:2] * 180/pi
# return
out
}
#' @rdname coordinates
#' @export
# Transform spherical to geographical coordinates
sph2geo <- function(ch_pos, r = 1, deg = TRUE, long360 = TRUE,
orient = c("northpole", "equatorial")) {
ch_pos <- as.data.frame(ch_pos)
if (!all(c("theta", "phi") %in% colnames(ch_pos)))
stop("Either theta or phi angles are missing!")
if (!is.null(ch_pos$r)) r <- ch_pos$r
if (!deg) ch_pos <- ch_pos * 180/pi
orient <- match.arg(orient)
if (orient == "equatorial") {
temp <- sph2cart(ch_pos)[, c("z", "x", "y")]
colnames(temp) <- c("x", "y", "z")
temp <- cart2sph(temp)
ch_pos$theta <- temp$theta
ch_pos$phi <- temp$phi
}
theta <- ch_pos$theta
phi <- ch_pos$phi
long <- ifelse(theta < 0, 180, 0) + phi
long[long < 0] <- 360 + long[long < 0]
if (!long360) long[long > 180] <- long[long > 180] - 360
lat <- 90 - abs(theta)
out <- data.frame(long, lat, r)
rownames(out) <- rownames(ch_pos)
# return
out
}
#' @rdname coordinates
#' @export
# Transform geographical to spherical coordinates
geo2sph <- function(ch_pos, r = 1, deg = TRUE) {
ch_pos <- as.data.frame(ch_pos)
if (!all(c("lat", "long") %in% colnames(ch_pos)))
stop("Either latitude [lat] or longitude [long] coordinates are missing!")
if (!is.null(ch_pos$r)) r <- ch_pos$r
neglong <- ch_pos$long < 0
long360 <- if (any(neglong)) F else T
theta <- (90 - ch_pos$lat) * sign(90 - ch_pos$long) * sign(270 - ch_pos$long)
theta <- ifelse(ch_pos$lat == 45 & theta == 0, 45, theta)
phi <- -(ifelse(theta < 0, 180, 0) - ch_pos$long)
phi <- ifelse(phi > 180, phi - 360, phi)
out <- data.frame(r = r, theta = theta, phi = phi)
backcheck <- sum(abs(as.matrix(
sph2geo(out, long360 = long360)[, c("long", "lat")] -
ch_pos[, c("long", "lat")]))) < 1e-8
if (backcheck) {
if (!deg) out <- out/180 * pi
return(out)
} else {
stop("Oops, something went wrong with the transformation and I don't know why.")
}
}
#' @rdname coordinates
#' @export
# Transform geographical to cartesian coordinates
geo2cart <- function(ch_pos, r = 1, deg = TRUE) {
#sph2cart( geo2sph(ch_pos, r, deg) )
ch_pos <- as.data.frame(ch_pos)
if (!all(c("long", "lat") %in% colnames(ch_pos)))
stop("Either longitudes (long) or latitudes (lat) are missing!")
if (!is.null(ch_pos$r)) r <- ch_pos$r
if (deg) {
ch_pos[, c("long", "lat")] <- ch_pos[, c("long", "lat")] * pi/180
}
long <- ch_pos$long
lat <- ch_pos$lat
x <- r * cos(long) * cos(lat)
y <- r * sin(long) * cos(lat)
z <- r * sin(lat)
out <- data.frame(x = x, y = y, z = z)
rownames(out) <- rownames(ch_pos)
# return
out
}
#' @rdname coordinates
#' @export
# Transform cartesian to geographical coordinates
cart2geo <- function(ch_pos, deg = TRUE) {
ch_pos <- as.data.frame(ch_pos)
if (!all(c("x", "y", "z") %in% colnames(ch_pos)))
stop("Either x, y, or z coordinates are missing!")
x <- ch_pos$x
y <- ch_pos$y
z <- ch_pos$z
r <- sqrt(x^2 + y^2 + z^2)
long <- atan2(y, x)
lat <- asin(z/r)
lat[r == 0] <- 0
out <- data.frame(long, lat, r)
rownames(out) <- rownames(ch_pos)
if (deg) out[, 1:2] <- out[, 1:2] * 180/pi
# return
out
}
#' Find channel neighbours
#'
#' \code{chanNb} finds neighbouring channels
#' @param ch_pos a data.frame of electrode positions,
#' see \code{\link{coordinates}}
#' @param check_alpha a two-element numeric vector defining the range which is
#' supposed to contain the optimal value of alpha
#' @param alpha a numeric value which influences the allowed distance between
#' neighbouring electrodes; if other than NULL (the default), check_alpha is
#' ignored
#' @param ... parameters to \code{\link{sph2cart}}
#' @export
#' @return An electrode neighbourhood matrix
chanNb <- function(ch_pos, check_alpha = c(0.1, 10), alpha = NULL, ...) {
options(rgl.useNULL = TRUE)
reqFn(c("alphashape3d", "geometry"))
if (!all(c("x", "y", "z") %in% colnames(ch_pos))) {
if (all(c("theta", "phi") %in% colnames(ch_pos))) {
ch_pos <- sph2cart(ch_pos, ...)
} else {
stop(paste0("Channel coordinates should be spherical (polar)",
"or cartesian coordinates!"))
}
}
ch_pos <- ch_pos[, c("x", "y", "z")]
channames <- paste(1:nrow(ch_pos), rownames(ch_pos), sep = ". ")
if (is.null(alpha)) {
if (!(requireNamespace("shiny") & requireNamespace("shinyRGL"))) {
stop(paste0(
"If you wish to determine alpha by an interactive plot, ",
"install the 'shinyRGL' package first."), call. = FALSE)
}
alpha <- shiny::runApp(list(
ui = shiny::pageWithSidebar(
# Application title
shiny::headerPanel("Find channel neighbours"),
# Sidebar with a slider input for number of points
shiny::sidebarPanel(
shiny::sliderInput("alpha",
"Alpha value: ",
min = min(check_alpha),
max = max(check_alpha),
value = 1, step = 0.1),
shiny::actionButton("submit", "Use selected alpha")
),
# Show the generated 3d scatterplot
shiny::mainPanel(
shinyRGL::webGLOutput("chanPlot", height = "700px")
)
),
server = function(input, output) {
output$chanPlot <- shinyRGL::renderWebGL({
#bg3d("grey40")
a <- suppressWarnings(alphashape3d::ashape3d(as.matrix(ch_pos),
input$alpha, pert = TRUE))
plot(a, walpha = TRUE, transparency = 0.95,
col = c("red", "red2", "red"), shininess = 100)
rgl::text3d(ch_pos$x, ch_pos$y, ch_pos$z, cex = 0.6,
channames, col = "black")
})
shiny::observe({
if (input$submit == 0)
return()
shiny::stopApp(input$alpha)
})
}
))
}
a <- suppressWarnings(alphashape3d::ashape3d(as.matrix(ch_pos),
alpha, pert = TRUE))$edge
a <- a[,c(1, 2, ncol(a))]
out <- matrix_(0, nrow(ch_pos), nrow(ch_pos))
rownames(out) <- rownames(ch_pos)
mirror <- out > 0
for (i in 1:nrow(ch_pos)) {
nb <- sort(unique(c(
which(mirror[i,]),
i,
a[a[,1]==i & a[,3]>0, 2])))
out[i, 1:length(nb)] <- nb
mirror[nb, i] <- T
}
out <- out[, apply(out > 0, 2, any)]
# return
out
}
#' Cosine of the angles between electrodes
#'
#' \code{cosAngle} computes the cosine of angles between electrodes or
#' between any N-dimensional vectors, or between two sets of electrodes or
#' any N-dimensional vectors.
#' @param x a matrix or data.frame with at least two columns and rows
#' @param y If not NULL (default), a matrix or data.frame with at least two
#' columns and rows
#' @param coords a logical variable; if TRUE (default), data in x are
#' coordinates. Coordinates should be named as "x", "y", "z" or "theta", "phi",
#' and if names are present, only coordinate columns are used in the
#' computations.
#' @param units_in_rows a logical variable indicating whether the units (vectors)
#' are in the rows of x (TRUE, default) or in the columns.
#' @param check_params logical; if TRUE (default), the appropriateness of
#' input data are checked before computing the cosine of angles. Set it to
#' FALSE if you really know what you are doing.
#' @export
#' @return A symmetric matrix
cosAngle <- function(x, y = NULL, coords = TRUE, units_in_rows = TRUE,
check_params = TRUE) {
checkFn <- function(xx) {
if (!units_in_rows) xx <- t(xx)
if (length(dim(xx)) <= 1 || ncol(xx)<2) {
stop("Provide a matrix or data.frame with at least 2 columns.")
}
if (coords) {
if (!is.null(colnames(xx))) {
if (!all(c("x", "y", "z") %in% colnames(xx)) &&
all(c("theta", "phi") %in% colnames(xx))) {
xx <- sph2cart(xx)
}
ind <- na.omit(match(c("x", "y", "z"), colnames(xx)))
if (length(ind) < 2) {
stop("Provide a matrix or data.frame with x and y coordinates.")
} else {
xx <- xx[, ind]
}
}
}
if (is.data.frame(xx)) xx <- as.matrix(xx)
return( xx )
}
if (is.null(y)) {
if (check_params) {
x <- checkFn(x)
}
len <- sqrt(rowSums(x^2))
out <- (x %*% t(x)) / outer(len, len)
rownames(out) <- colnames(out) <- rownames(x)
} else {
if (check_params) {
x <- checkFn(x)
y <- checkFn(y)
}
lenx <- sqrt(rowSums(x^2))
leny <- sqrt(rowSums(y^2))
out <- (x %*% t(y)) / outer(lenx, leny)
rownames(out) <- rownames(x)
colnames(out) <- rownames(y)
}
# return
out
}
#' Spherical spline interpolation
#'
#' \code{chanInterp} performs spherical spline interpolation. It can impute
#' bad channels (i.e. channels with missing values) or interpolate to
#' artbitrary positions on the scalp.
#' @param dat a numeric vector, matrix, data.frame or array with named dimnames
#' @param ch_pos channel (electrode) positions; should be a matrix or data.frame
#' with the following column names (order is not important): "theta", "phi"
#' (spherical coordinates), or "x", "y", "z" (cartesian coordinates)
#' @param interp_pos interpolation positions, if they differ from ch_pos;
#' should be in the same format as ch_pos (default: NULL)
#' @param maxNA numeric value given as ratio (0 < maxNA < 1) or integer
#' (maxNA >= 1 ); it determines the maximum ratio or number of non-missing
#' channels in a time sample. If exceeded, no interpolation occurs in the given
#' time sample.
#' @param m integer value (default: 4); the flexibility of the spline
#' @param N integer value (default: 7); number of terms in the Legendre
#' polynomial. You should probably increase it (even up to 100L) for
#' realistic (i.e. non-spherical) or high-density electrode arrangements.
#' @param lambda numeric value (default: 1e-10); smoothing factor
#' @param type character value, do not use yet
#' @param alarm_tolerance numeric value (default: 1e-2); if the maximal absolute
#' interpolation error at any time sample exceeds this limit, a message is
#' shown or an error is thrown depending on \code{error_on_alarm}. If set to
#' NULL, no check is performed.
#' @param error_on_alarm defaults to TRUE, see \code{alarm_tolerance}
#' @export
#' @return An object having the same attributes as dat
chanInterp <- function(dat, ch_pos, interp_pos = NULL, maxNA = 0.3,
m = 4L, N = 7L, lambda = 1e-10,
type = c("voltage", "laplacian", "scd"),
alarm_tolerance = 1e-2, error_on_alarm = TRUE) {
message("\n****\nStart interpolation...\n")
# Function to compute G matrices:
# G = g(cos(ch_pos)) / interpG = g(cos(ch_pos, interp_pos))
compGmat <- function(cos_angles, m, N) {
n <- seq.int(N)
series <- (2*n + 1) / (n^m * (n+1)^m)
k <- 1/4 / pi
P <- polynomial.values(
legendre.polynomials(N, normalized = FALSE)[-1],
cos_angles)
G <- k * Reduce("+", mapply("*", P, series, SIMPLIFY = FALSE))
return( G )
}
# Function to compute C coefficients
compCmat <- function(x, G, tol) {
Gx <- array_(1, dim(G) + 1)
Gx[-1, -1] <- G
Gx[1, 1] <- 0
diag(Gx) <- diag(Gx) + tol
invG <- corpcor::pseudoinverse(Gx, tol)[, -1]
return( invG %*% x )
}
# Function to perform spherical spline interpolation
interpFn <- function(y, m, N, lambda, type) {
na_ind <- is.na(y)
message("...Looking for missing value patterns - ", appendLF = FALSE)
if (is.null(interp_pos)) {
temp <- colMeans(na_ind)
keep_columns <- temp <= maxNA & temp > 0
if (!any(keep_columns)) return(y)
yorig <- y
y <- y[, keep_columns, drop = FALSE]
na_ind <- na_ind[, keep_columns, drop = FALSE]
missing_patterns <- fastUnique(na_ind, margin = 2L)
mischan <- which(rowAnys(missing_patterns))
mischan <-
if (!is.null(rownames(ch_pos))) {
paste(rownames(ch_pos)[mischan], collapse = ", ")
} else {
paste(seq_along(mischan), collapse = ", ")
}
message("Done")
message("...There are missing values in the following channels: ",
mischan)
} else {
missing_patterns <- fastUnique(na_ind, margin = 2L)
message("Done")
yy <- array_(0, c(nrow(interp_pos), ncol(y)),
list(rownames(interp_pos), colnames(y)))
names(dimnames(yy)) <- names(dimnames(y))
}
if (!is.null(alarm_tolerance)) {
dev <- rep(FALSE, ncol(y))
names(dev) <- colnames(y)
maxdev <- 0
}
message("...Perform interpolation - ", appendLF = FALSE)
G0 <- compGmat(cosAngle(ch_pos, check_params = FALSE), m, N)
for (i in 1:ncol(missing_patterns)) {
na_vec <- missing_patterns[, i]
y_ind <- colSums(na_ind==na_vec) == length(na_vec)
ch_good <- ch_pos[!na_vec, , drop = F]
ch_interp <-
if (is.null(interp_pos)) ch_pos[na_vec,,drop = F] else interp_pos
#
G <- G0[!na_vec, !na_vec]
Coef <- compCmat(y[!na_vec, y_ind, drop = F], G, lambda)
interpG <- compGmat(cosAngle(ch_good, ch_interp, check_params = FALSE),
m, N)
if (is.null(interp_pos)) {
res <- crossprod(Coef[-1, , drop = F], interpG) + Coef[1, ]
y[na_vec, y_ind] <- t(res)
} else {
yy[, y_ind] <- t(crossprod(Coef[-1, , drop = F], interpG) +
Coef[1, ])
}
if (!is.null(alarm_tolerance)) {
ch_interp <- ch_pos
interpG <- compGmat(cosAngle(ch_good,
ch_interp, check_params = FALSE),
m, N)
res <- crossprod(Coef[-1, , drop = F], interpG) + Coef[1, ]
mdevs <- colMaxs(abs(y[, y_ind] - t(res)), na.rm = TRUE)
maxdev <- max(c(maxdev, mdevs))
dev[y_ind] <- temp <- ( maxdev > alarm_tolerance )
if (error_on_alarm && temp)
stop(sprintf("Deviance %f exceeds threshold %f",
maxdev, alarm_tolerance))
}
}
if (!is.null(alarm_tolerance)) {
if (any(dev)) {
devind <- which(dev)
devind <- devind[1:min(c(20, length(devind)))]
message(
"****\nThe maximum absolute interpolation error exceeds the limit
at the following time samples (only the first 20 are shown):\n-----\n",
paste(names(dev)[devind], collapse = " ")
)
message(paste("Maximal deviation = ", maxdev, sep = ""))
}
}
message("Done")
if (is.null(interp_pos)) {
yorig[, keep_columns] <- y
return(yorig)
} else {
return(yy)
}
}
#
type <- match.arg(type)
if (is.null(lambda)) {
lambda <- if (type == "voltage") 1e-7 else 1e-5
}
#
# check if ch_pos has an appropriate format
if (!all(c("x", "y", "z") %in% colnames(ch_pos))) {
if (all(c("theta", "phi") %in% colnames(ch_pos))) {
ch_pos <- sph2cart(ch_pos)
} else {
stop("Provide channel locations in spherical or cartesian coordinates!")
}
}
ch_pos <- as.matrix(ch_pos[, c("x", "y", "z")])
#
# do the same for the interpolation locations
if (!is.null(interp_pos)) {
if (!all(c("x", "y", "z") %in% colnames(interp_pos))) {
if (all(c("theta", "phi") %in% colnames(interp_pos))) {
interp_pos <- sph2cart(interp_pos)
} else {
stop("Provide interpolation locations in spherical or cartesian coordinates!")
}
}
interp_pos <- as.matrix(interp_pos[, c("x", "y", "z")])
}
#
# set up interpolation infos
argnames <- setdiff(names(as.list(args(chanInterp))), "dat")
argnames <- argnames[argnames != ""]
procstep <- list(
what = "interpolation", call = match.call(),
options = mget(argnames), nr_of_missings = sum(is.na(dat)))
#
# check input data and run interpolation
dim_names <- dimnames(dat)
if (is.list(dat) || !is.numeric(dat)) {
stop("Provide a numeric vector, matrix, data.frame or array as input!")
} else if (is.null(interp_pos) &&
procstep$nr_of_missings == 0) {
setattr(dat, "processing_steps",
c(attr(dat, "processing_steps"), list(procstep)))
message("...No missing data to interpolate.\nInterpolation finished.")
return(dat)
} else if (is.vector(dat)) {
if (length(dat) != nrow(ch_pos)) {
stop("Number of channels and number of data points do not match!")
}
out <- interpFn(as.matrix(dat), m, N, lambda, type)
} else {
if (!any(dim(dat) == nrow(ch_pos))) {
stop("Number of channels and data size do not match!")
}
target_dim <- if (!is.null(dim_names$chan)) "chan" else 1
arg_list <- list(m = m, N = N, lambda = lambda, type = type)
if (length(arg_list) == 0) arg_list <- NULL
newdims <-
if (is.null(interp_pos)) list() else list(chan = rownames(interp_pos))
out <- fnDims(dat, "chan", interpFn, arg_list = arg_list,
newdims = newdims, vectorized = TRUE)
out <- apermArray(out, names(dim_names),
keep_attributes. = TRUE)
}
#
if (is.null(interp_pos)) {
NAs_left <- sum(is.na(out))
procstep$nr_of_interp <-
procstep$nr_of_missings - NAs_left
message("Number of NAs after bad channel interpolation: ",
NAs_left)
}
message("Interpolation finished.")
setattr(out, "processing_steps",
c(attr(dat, "processing_steps"), list(procstep)))
# return
out
}
#' Project channel positions onto 2D plane
#'
#' \code{project3dMap} projects 3D coordinates onto a 2D plane or vice versa
#' @param pos electrode positions (matrix or data.frame)
#' @param r radius
#' @param projection character string (default = "laea"). See
#' \url{http://www.remotesensing.org/geotiff/proj_list/} for common projections
#' @param projref projection reference (pole [ = default] or equator)
#' @param origo a named character vector of lat ( = latitude) and long (longitude)
#' @param inverse if set to TRUE, back-projection is performed (default = FALSE)
#' @export
#' @return A data.frame containing the projected coordinates
project3dMap <- function(pos, r = 1,
projection = "laea", projref = c("pole", "equator"),
origo = c(lat = ifelse(projref == "pole", 90, 0),
long = ifelse(projref == "pole", 270, 0)),
inverse = FALSE) {
#
xyfn <- function(dat) data.frame(x = dat[,1], y = dat[,2])
projfn <- function(p, projcall, ...) {
res <- lapply(unique(projcall), function(x) {
ind <- projcall == x
out <- project(as.data.frame(p[ind, , drop = F]), proj = x, ...)
out <- as.data.frame(out)
rownames(out) <- which(ind)
return( out )
})
res <- Reduce("rbind", res)
res <- res[order(as.numeric(rownames(res))), , drop = F]
return( xyfn(res) )
}
#
if (!is.data.frame(pos) & is.list(pos)) {
stop("Pos can be a vector, matrix, or data.frame.")
}
if (is.atomic(pos) && is.vector(pos)) {
if (length(pos) != 2) {
stop("If pos is a vector, it must contain exactly two elements!")
}
pos <- cbind(x = pos[1], y = pos[2])
}
pos <- as.data.frame(pos)
if (is.null(pos$r)) pos$r <- 1
projref <- match.arg(projref)
proj4call <- paste0("+proj=", projection,
" +lat_0=", origo[1],
" +lon_0=", origo[2],
paste0(" +R=", pos$r))
#
if (inverse) {
if (is.null(colnames(pos))) colnames(pos)[1:2] <- c("x", "y")
if (!all(c("x", "y") %in% colnames(pos))) {
stop("Provide a matrix or data.frame with x and y coordinates!")
}
mapxy <- projfn(as.matrix(pos[, c("x", "y")]), proj4call, inverse = TRUE)
colnames(mapxy) <- c("long", "lat")
} else {
if (!all(c("long", "lat") %in% colnames(pos))) {
if (!all(c("theta", "phi") %in% colnames(pos)))
pos <- cart2sph(pos)
if (projref == "equator") {
temp <- sph2cart(pos)
colnames(temp) <- c("y", "z", "x")
pos <- cart2sph(temp)
}
pos <- sph2geo(pos, long360 = FALSE)
}
mapxy <- projfn(pos[, c("long", "lat")], proj4call)
}
rownames(mapxy) <- rownames(pos)
# return
mapxy
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.