From 8b814e080e1a2a49b438a2518ba55c6245922bf8 Mon Sep 17 00:00:00 2001 From: dipterix Date: Wed, 17 Apr 2024 16:33:13 -0400 Subject: [PATCH] `plot_slices` allows overlay images --- .Rbuildignore | 3 +- CHANGELOG.md | 5 +- DESCRIPTION | 2 +- R/class_brain.R | 27 ++++++- R/plot_volume-slices.R | 177 +++++++++++++++++++++++++++++++++-------- man/plot_slices.Rd | 11 ++- 6 files changed, 187 insertions(+), 38 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index 14fbe499..c2b04d4c 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,7 +8,8 @@ ^cran-comments\.md$ ^shortcuts\.md$ ^README\.md$ -^adhoc/ +^adhoc +.ipynb_checkpoints .git .Rprofile .Rhistory diff --git a/CHANGELOG.md b/CHANGELOG.md index 0604572b..9192dd61 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,7 @@ ## Changes since last CRAN release -* `e3c54181 (HEAD -> master)` [_`dipterix`_]: Fixed right-click not settings crosshairs correctly on template brain -* `8ccd7706 (origin/master, origin/HEAD)` [_`dipterix`_]: Drag and drop feature supports `STL` format now; Color map is remembered regardless of the file extension +* `e66ea42d (HEAD -> master)` [_`dipterix`_]: `plot_slices` allows overlay images +* `5ad19c55 (origin/master, origin/HEAD)` [_`dipterix`_]: Fixed right-click not settings crosshairs correctly on template brain +* `8ccd7706` [_`dipterix`_]: Drag and drop feature supports `STL` format now; Color map is remembered regardless of the file extension * `d6de9e10` [_`dipterix`_]: Added `conform_volume` to conform images that simulate `FreeSurfer` conform algorithm * `40d10a3b` [_`dipterix`_]: Exported `write.fs.surface` * `8f356c31` [_`dipterix`_]: `volume_to_surf` can take objects (from `read_volume`) as input diff --git a/DESCRIPTION b/DESCRIPTION index c030b014..1ea1bc9d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: threeBrain Type: Package Title: Your Advanced 3D Brain Visualization -Version: 1.1.0.9004 +Version: 1.1.0.9005 Authors@R: c( person("Zhengjia", "Wang", email = "dipterix.wang@gmail.com", role = c("aut", "cre", "cph")), person("John", "Magnotti", email = "John.Magnotti@Pennmedicine.upenn.edu", role = c("aut", "res")), diff --git a/R/class_brain.R b/R/class_brain.R index 43669cde..abc482e0 100644 --- a/R/class_brain.R +++ b/R/class_brain.R @@ -910,7 +910,12 @@ Brain2 <- R6::R6Class( }, plot_electrodes_on_slices = function( - electrodes_to_plot = "all", volume = NULL, elec_table = NULL, + electrodes_to_plot = "all", + volume = NULL, + overlays = NULL, + overlay_colors = DEFAULT_COLOR_DISCRETE, + overlay_alpha = 0.3, + elec_table = NULL, zoom = 1, adjust_brightness = NA, electrode_color = "green", electrode_size = 2, verbose = TRUE, ..., @@ -976,6 +981,24 @@ Brain2 <- R6::R6Class( volume$data[volume$data > 255] <- 255 } + if( length(overlays) ) { + if(length(overlay_colors) < length(overlays)) { + overlay_colors <- rep(overlay_colors, ceiling(length(overlays) / length(overlay_colors))) + } + overlays <- lapply(seq_along(overlays), function(ii) { + overlay_img <- overlays[[ii]] + if(is.character(overlay_img)) { + overlay_img <- read_volume(overlay_img) + } + list( + volume = overlay_img, + color = overlay_colors[[ii]] + ) + }) + } else { + overlays <- NULL + } + # img_height <- 480 # png(filename = file.path(subject$imaging_path, "snapshots%03d.png"), width = 3*img_height, height = img_height, bg = "black") @@ -1006,6 +1029,8 @@ Brain2 <- R6::R6Class( progress$inc(detail = sprintf("Generating graphs for electrode %s", dipsaus::deparse_svec(ii))) plot_slices( volume, + overlays = overlays, + overlay_alpha = overlay_alpha, positions = scanner_ras[ii,], main = sprintf('%s (Ch=%.0f,ScanRAS=%.1f,%.1f,%.1f)', elec_table$Label[ii], diff --git a/R/plot_volume-slices.R b/R/plot_volume-slices.R index 75d277df..5bf1a0f7 100644 --- a/R/plot_volume-slices.R +++ b/R/plot_volume-slices.R @@ -1,5 +1,10 @@ #' Plot slices of volume -#' @param volume path to volume +#' @param volume path to volume (underlay) +#' @param overlays images to overlay on top of the underlay, can be either +#' a vector of paths to the overlay volume images, or a sequence of named lists. +#' Each list item has \code{'volume'} (path to the volume) and \code{'color'} +#' (color of the overlay) +#' @param overlay_alpha transparency of the overlay; default is 0.3 #' @param transform rotation of the volume in scanner \code{'RAS'} space #' @param positions vector of length 3 or matrix of 3 columns, the \code{'RAS'} #' position of cross-hairs @@ -25,8 +30,8 @@ #' @returns Nothing #' @export plot_slices <- function( - volume, transform = NULL, positions = NULL, zoom = 1, pixel_width = 0.5, - col = c("black", "white"), normalize = NULL, zclip = NULL, + volume, overlays = NULL, transform = NULL, positions = NULL, zoom = 1, pixel_width = 0.5, + col = c("black", "white"), normalize = NULL, zclip = NULL, overlay_alpha = 0.3, zlim = normalize, main = "", title_position = c("left", "top"), fun = NULL, nc = NA, which = NULL, ...) { # DIPSAUS DEBUG START @@ -39,6 +44,15 @@ plot_slices <- function( # positions = rnorm(12) # nc <- 1 # title_position <- "top" + # overlays <- "~/rave_data/raw_dir/YAB/rave-imaging/fs/mri/aseg.mgz" + # + # # test with overlays + # brain <- raveio::rave_brain("YAEL/AnonSEEG") + # root_path <- "/Users/dipterix/Library/CloudStorage/Box-Box/RAVEExternal/YAEL-Andrew-Stanford/rave_data/raw_dir/AnonSEEG/rave-imaging/atlases/OCD Response Tract Atlas (Li 2020)/mixed/" + # pixel_width <- 1 + # volume <- file.path(brain$base_path, "mri", "T1.mgz") + # overlays <- list.files(root_path, pattern = "\\.nii.gz$", full.names = TRUE) + # overlay_alpha <- 0.3 title_position <- match.arg(title_position) @@ -48,6 +62,38 @@ plot_slices <- function( if(!inherits(volume, "threeBrain.volume")) { stop("`volume` must be character or threeBrain.volume") } + + overlays <- lapply(seq_along(overlays), function(ii) { + item <- overlays[[ ii ]] + default_color <- col2hexStr(DEFAULT_COLOR_DISCRETE[[(ii - 1L) %% length(DEFAULT_COLOR_DISCRETE) + 1L]], alpha = overlay_alpha) + if(is.character(item)) { + ovol <- read_volume(item) + return(list( + volume = ovol, + color = default_color, + world2ijk = solve(ovol$Norig) + )) + } + if(!is.list(item)) { return(NULL) } + if(inherits(item, "threeBrain.volume")) { + return(list( + volume = item, + color = default_color, + world2ijk = solve(item$Norig) + )) + } + if(!"volume" %in% names(item)) { + return(NULL) + } + if(length(item$color) != 1) { + item$color <- default_color + } else { + item$color <- col2hexStr(item$color, alpha = overlay_alpha) + } + item$world2ijk <- solve(item$volume$Norig) + return( item ) + }) + if(is.null(transform)) { transform <- diag(1, 4) } else { @@ -113,8 +159,11 @@ plot_slices <- function( more_args$zlim <- zlim more_args$useRaster <- TRUE more_args$main <- '' + more_args$xlab <- "" + more_args$ylab <- "" more_args$x <- x more_args$y <- x + default_add <- isTRUE(more_args$add) # more_args$add <- FALSE oldpar <- graphics::par(no.readonly = TRUE) @@ -193,6 +242,37 @@ plot_slices <- function( panel_last <- function(...) {} } + # create template positions to calculate world positions + # pre-apply inverse of the rotation to template position + wpos_axial <- transform_inv %*% pos # varying RA, S=0 -> axial + wpos_sagittal <- transform_inv %*% pos[c(3,1,2,4), , drop = FALSE] # varying AS, R=0 -> sagittal + wpos_coronal <- transform_inv %*% pos[c(1,3,2,4), , drop = FALSE] # varying RS, A=0 -> coronal + + get_ijk <- function(point_position, which_slice, world2ijk, volume_data) { + wpos <- switch( + which_slice, + "axial" = { wpos_axial + point_position }, + "sagittal" = { wpos_sagittal + point_position }, + "coronal" = { wpos_coronal + point_position }, + { stop("Invalid slice type") } + ) + # world to voxel index + IJK <- round( world2ijk[c(1, 2, 3), ] %*% wpos ) + + # Remove invalid indices + shape <- dim(volume_data) + cumshape <- cumprod(c(1, shape))[seq_len(3)] + sel <- IJK[1,] >= shape[[1]] | IJK[2,] >= shape[[2]] | IJK[3,] >= shape[[3]] + IJK[, sel] <- NA + IJK[!is.na(IJK) & IJK < 0] <- NA + + # actual indices + idx <- cumshape %*% IJK + 1L + slice <- volume_data[idx] + dim(slice) <- c(nx, nx) + slice + } + lapply(seq_len(npts), function(ii) { pos_pt <- c(positions[ii, ], 0) @@ -211,51 +291,84 @@ plot_slices <- function( if(!length(which) || 1 %in% which) { # Axial # translate x transform_inv x translate^-1 x Norig - IJK <- round(world2ijk[c(1, 2, 3), ] %*% (transform_inv %*% pos + pos_pt)) + 1L - sel <- IJK[1,] > shape[[1]] | IJK[2,] > shape[[2]] | IJK[3,] > shape[[3]] - IJK[,sel] <- NA - IJK[IJK < 1] <- NA - idx <- t(IJK - 1) %*% cumshape + 1 - slice <- nu(volume$data[idx]) - - dim(slice) <- c(nx, nx) - more_args$z <- slice + + underlay <- get_ijk(pos_pt, which_slice = "axial", world2ijk = world2ijk, volume$data) + more_args$add <- default_add + more_args$col <- pal + more_args$z <- nu(underlay) adjust_plt() do.call(graphics::image, more_args) + + more_args$z <- NULL + more_args$add <- TRUE + lapply(overlays, function(item) { + if(!is.list(item)) { return() } + overlay <- get_ijk(point_position = pos_pt, which_slice = "axial", world2ijk = item$world2ijk, volume_data = item$volume$data) + invalids <- is.na(overlay) | overlay < 0.5 + if(!all(invalids)) { + overlay[invalids] <- NA + more_args$col <- item$color + more_args$z <- overlay + do.call(graphics::image, more_args) + } + }) + panel_last( ii, 1 ) } if(!length(which) || 2 %in% which) { # Sagittal - IJK <- round(world2ijk[c(1, 2, 3), ] %*% (pos[c(3,1,2,4), , drop = FALSE] + pos_pt)) + 1L - sel <- IJK[1,] > shape[[1]] | IJK[2,] > shape[[2]] | IJK[3,] > shape[[3]] - IJK[,sel] <- NA - IJK[IJK < 1] <- NA - idx <- t(IJK - 1) %*% cumshape + 1 - slice <- nu(volume$data[idx]) - - dim(slice) <- c(nx, nx) - more_args$z <- slice + underlay <- get_ijk(pos_pt, which_slice = "sagittal", world2ijk = world2ijk, volume$data) + more_args$add <- default_add + more_args$col <- pal + more_args$z <- nu(underlay) adjust_plt() do.call(graphics::image, more_args) - panel_last( ii, 2 ) + + more_args$z <- NULL + more_args$add <- TRUE + lapply(overlays, function(item) { + if(!is.list(item)) { return() } + overlay <- get_ijk(pos_pt, which_slice = "sagittal", world2ijk = item$world2ijk, item$volume$data) + invalids <- is.na(overlay) | overlay < 0.5 + if(!all(invalids)) { + overlay[invalids] <- NA + overlay_color <- item$color + more_args$col <- overlay_color + more_args$z <- overlay + do.call(graphics::image, more_args) + } + }) + + panel_last( ii, 1 ) } if(!length(which) || 3 %in% which) { # Coronal - IJK <- round(world2ijk[c(1, 2, 3), ] %*% (pos[c(1,3,2,4), , drop = FALSE] + pos_pt)) + 1L - sel <- IJK[1,] > shape[[1]] | IJK[2,] > shape[[2]] | IJK[3,] > shape[[3]] - IJK[,sel] <- NA - IJK[IJK < 1] <- NA - idx <- t(IJK - 1) %*% cumshape + 1 - slice <- nu(volume$data[idx]) - - dim(slice) <- c(nx, nx) - more_args$z <- slice + underlay <- get_ijk(pos_pt, which_slice = "coronal", world2ijk = world2ijk, volume$data) + more_args$add <- default_add + more_args$col <- pal + more_args$z <- nu(underlay) adjust_plt() do.call(graphics::image, more_args) - panel_last( ii, 3 ) + + more_args$z <- NULL + more_args$add <- TRUE + lapply(overlays, function(item) { + if(!is.list(item)) { return() } + overlay <- get_ijk(pos_pt, which_slice = "coronal", world2ijk = item$world2ijk, item$volume$data) + invalids <- is.na(overlay) | overlay < 0.5 + if(!all(invalids)) { + overlay[invalids] <- NA + overlay_color <- item$color + more_args$col <- overlay_color + more_args$z <- overlay + do.call(graphics::image, more_args) + } + }) + + panel_last( ii, 1 ) } NULL diff --git a/man/plot_slices.Rd b/man/plot_slices.Rd index 6aa16cc3..35359a8d 100644 --- a/man/plot_slices.Rd +++ b/man/plot_slices.Rd @@ -6,6 +6,7 @@ \usage{ plot_slices( volume, + overlays = NULL, transform = NULL, positions = NULL, zoom = 1, @@ -13,6 +14,7 @@ plot_slices( col = c("black", "white"), normalize = NULL, zclip = NULL, + overlay_alpha = 0.3, zlim = normalize, main = "", title_position = c("left", "top"), @@ -23,7 +25,12 @@ plot_slices( ) } \arguments{ -\item{volume}{path to volume} +\item{volume}{path to volume (underlay)} + +\item{overlays}{images to overlay on top of the underlay, can be either +a vector of paths to the overlay volume images, or a sequence of named lists. +Each list item has \code{'volume'} (path to the volume) and \code{'color'} +(color of the overlay)} \item{transform}{rotation of the volume in scanner \code{'RAS'} space} @@ -43,6 +50,8 @@ one pixel is 0.5 millimeters wide} \item{zclip}{clip image densities; if specified, values outside of this range will be clipped into this range} +\item{overlay_alpha}{transparency of the overlay; default is 0.3} + \item{zlim}{image plot value range, default is identical to \code{normalize}} \item{main}{image titles}