Added get_atlas_values for continuous atlases such as binary or pro…

…bablistic `ROI`. The radius refers to maximum `RAS` distance instead of voxel indexing distance, hence more accurate when the atlas volume has imbalanced slice count.
dipterix · Oct 6, 2024 · 0c3164e · 0c3164e
1 parent 161545b
commit 0c3164e
Show file tree

Hide file tree

Showing 4 changed files with 180 additions and 11 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,7 +1,8 @@
 ## Changes since last CRAN release
-* `11c7a6e5 (HEAD -> master)` [_`dipterix`_]: Renamed `active-voxel` to `column-row-slice`
+* `7370d864 (HEAD -> master)` [_`dipterix`_]: Added `get_atlas_values` for continuous atlases such as binary or probablistic `ROI`. The radius refers to maximum `RAS` distance instead of voxel indexing distance, hence more accurate when the atlas volume has imbalanced slice count.
+* `161545b2 (origin/master, origin/HEAD)` [_`dipterix`_]: Renamed `active-voxel` to `column-row-slice`
 * `ad4de70a` [_`dipterix`_]: `Voxel` filter is linear now when displayed at side slices only and when the slice mode is not `active-voxel`
-* `8f62ef5f (origin/master, origin/HEAD)` [_`dipterix`_]: Fixed drifting issue when visualizing via active `voxel` mode; Added direction arrow helper to `DBS` (or electrodes with non-zero model up vectors)
+* `8f62ef5f` [_`dipterix`_]: Fixed drifting issue when visualizing via active `voxel` mode; Added direction arrow helper to `DBS` (or electrodes with non-zero model up vectors)
 * `28fe8ae6` [_`dipterix`_]: Support `WebAssembly`
 * `b07713dd` [_`dipterix`_]: Avoid using pandoc to save the whole page self-contained
 * `bb7ba967` [_`dipterix`_]: Fixed the depth issue in electrode material shader

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: threeBrain
 Type: Package
 Title: Your Advanced 3D Brain Visualization
-Version: 1.1.1.9022
+Version: 1.1.1.9023
 Authors@R: c(
     person("Zhengjia", "Wang", email = "[email protected]", role = c("aut", "cre", "cph")),
     person("John", "Magnotti", email = "[email protected]", role = c("ctb", "res")),

diff --git a/R/class_brainelectrodes.R b/R/class_brainelectrodes.R
@@ -427,11 +427,96 @@ BrainElectrodes <- R6::R6Class(
       return(positions[, seq_len(3), drop = FALSE])
     },
 
-    get_atlas_labels = function(atlas, radius = 1, lut = NULL) {
+    get_atlas_values = function(atlas, radius = 1.5, ...) {
+      # DIPSAUS DEBUG START
+      # self <- raveio::rave_brain("demo/DemoSubject")$electrodes
+      # atlas <- "~/rave_data/raw_dir/DemoSubject/rave-imaging/fs/mri/aparc.a2009s+aseg.mgz"
+      # radius <- 2
+
+      if(!is.data.frame(self$raw_table)) {
+        return(NULL)
+      }
+
+      if(!inherits(atlas, "threeBrain.volume")) {
+        atlas <- read_volume(atlas)
+      }
+
+      sub_tbl <- as.matrix(self$raw_table[, c("Coord_x", "Coord_y", "Coord_z")])
+      dist <- rowSums(sub_tbl^2)
+      valids <- !is.na(dist) & dist > 0
+      scanner_ras <- self$apply_transform_points(sub_tbl, from = "tkrRAS", to = "scannerRAS")
+      ras_to_ijk <- solve(atlas$Norig)
+
+      # IJK starts from 0, 0, 0
+      ijk <- round((ras_to_ijk %*% rbind(t(scanner_ras), 1))[seq_len(3), , drop = FALSE])
+
+      atlas_shape <- dim(atlas$data)[seq_len(3)]
+      atlas_cumprod <- cumprod(c(1, atlas_shape))
+      atlas_n <- atlas_cumprod[[4]]
+      atlas_cumprod <- atlas_cumprod[seq_len(3)]
+
+      # construct neighboring indices
+      if( radius > 0 ) {
+        # columns of ras_to_ijk are incremental steps along voxel-index space
+        max_index_radius <- max(abs(ras_to_ijk[, 1:3])) * radius
+        # IJK offsets
+        deltas <- t(as.matrix(expand.grid(
+          seq.int(-max_index_radius, max_index_radius),
+          seq.int(-max_index_radius, max_index_radius),
+          seq.int(-max_index_radius, max_index_radius)
+        )))
+        # actual offsets in RAS
+        ras_delta <- atlas$Norig[1:3, 1:3] %*% deltas
+        # distance
+        distance <- sqrt(colSums(ras_delta^2))
+        sel <- distance <= radius
+        distance <- distance[sel]
+        deltas <- deltas[, sel, drop = FALSE]
+        deltas <- colSums(deltas * atlas_cumprod)
+      } else {
+        deltas <- 0
+        distance <- 0
+      }
+
+      voxel_count <- length(distance)
+      sel <- distance == 0
+
+      unknown_labels <- data.frame(
+        CenterValue = NA_real_,
+        MeanValue = NA_real_,
+        VoxelCount = NA_integer_
+      )
+
+      value_rows <- lapply(seq_len(ncol(ijk)), function(ii) {
+        if(!valids[[ii]]) {
+          return(unknown_labels)
+        }
+        ijk0 <- sum(ijk[, ii] * atlas_cumprod) + 1
+        if(is.na(ijk0) || ijk0 <= 0 || ijk0 > atlas_n) {
+          return(unknown_labels)
+        }
+        values <- atlas$data[ijk0 + deltas]
+        center_value <- values[ sel ]
+        values <- values[!is.na(values)]
+        if(!length(values)) {
+          return(unknown_labels)
+        }
+
+        data.frame(
+          CenterValue = center_value,
+          MeanValue = mean(values),
+          VoxelCount = length(values)
+        )
+      })
+      return(do.call("rbind", value_rows))
+    },
+
+    get_atlas_labels = function(atlas, radius = 1.5, lut = NULL) {
       # DIPSAUS DEBUG START
       # self <- raveio::rave_brain("demo/DemoSubject")$electrodes
       # atlas <- "~/rave_data/raw_dir/DemoSubject/rave-imaging/fs/mri/aparc.a2009s+aseg.mgz"
       # lut = NULL
+      # radius <- 2
 
       if(!is.data.frame(self$raw_table)) {
         return(NULL)
@@ -461,15 +546,28 @@ BrainElectrodes <- R6::R6Class(
       atlas_cumprod <- atlas_cumprod[seq_len(3)]
 
       lut_keys <- names(lut$map)
-      if(radius >= 1) {
-        deltas <- as.matrix(expand.grid(
-          seq.int(-radius, radius),
-          seq.int(-radius, radius),
-          seq.int(-radius, radius)
-        ))
-        deltas <- colSums(t(deltas) * atlas_cumprod)
+
+      # construct neighboring indices
+      if( radius > 0 ) {
+        # columns of ras_to_ijk are incremental steps along voxel-index space
+        max_index_radius <- max(abs(ras_to_ijk[, 1:3])) * radius
+        # IJK offsets
+        deltas <- t(as.matrix(expand.grid(
+          seq.int(-max_index_radius, max_index_radius),
+          seq.int(-max_index_radius, max_index_radius),
+          seq.int(-max_index_radius, max_index_radius)
+        )))
+        # actual offsets in RAS
+        ras_delta <- atlas$Norig[1:3, 1:3] %*% deltas
+        # distance
+        distance <- sqrt(colSums(ras_delta^2))
+        sel <- distance <= radius
+        distance <- distance[sel]
+        deltas <- deltas[, sel, drop = FALSE]
+        deltas <- colSums(deltas * atlas_cumprod)
       } else {
         deltas <- 0
+        distance <- 0
       }
 
       unknown_labels <- data.frame(

diff --git a/inst/prototypes/GeometryMaker/sEEG.R b/inst/prototypes/GeometryMaker/sEEG.R
@@ -336,3 +336,73 @@ proto <- seeg_prototype(
   overwrite = TRUE
 )
 
+# ---- sEEG-PMT-2102-??-SP350 ---------------------------------------------
+
+# 8 contact - 26.5mm recording length 2mm contacts 3.5mm spacing and .8mm diameter
+# 10 contact- 33.5mm recording length 2mm contacts 3.5mm spacing and .8mm diameter
+# 12 contact- 40.5mm recording length 2mm contacts 3.5mm spacing and .8mm diameter
+# 14 contact- 47.5mm recording length 2mm contacts 3.5mm spacing and .8mm diameter
+# 16 contact- 56.5mm recording length 2mm contacts 3.5mm spacing and .8mm diameter
+
+probe_head <- 0
+width <- 2
+contact_spacing <- 3.5
+# overall_length <- probe_head + width + contact_spacing * (n_contacts-1)
+diameter <- 0.8
+
+for( n_contacts in c(8, 10, 12, 14, 16) ) {
+  contacts <- probe_head + width / 2 + 0:(n_contacts-1) * contact_spacing
+  overall_length <- ceiling(max(contacts) + width / 2 + 300)
+  proto <- seeg_prototype(
+    type = sprintf("sEEG-PMT-2102-%02d-SP350", n_contacts),
+    description = c(
+      sprintf("PMT sEEG - %d contacts", n_contacts),
+      "Contact length   : 2   mm",
+      "Central spacing  : 3.5 mm",
+      "Tip size         : 0   mm",
+      "Diameter         : 0.8 mm"
+    ),
+    center_position = contacts,
+    contact_widths = width,
+    diameter = diameter,
+    overall_length = overall_length,
+    default_interpolation = sprintf("%.1fx%d", contact_spacing, n_contacts - 1L),
+    overwrite = TRUE
+  )
+}
+
+# ---- sEEG-PMT-2102-16-SP??? ---------------------------------------------
+
+# 16 contact- 61.5mm recording length 2mm contacts 3.97mm spacing and .8mm diameter
+# 16 contact- 68.5mm recording length 2mm contacts 4.43mm spacing and .8mm diameter
+
+probe_head <- 0
+width <- 2
+contact_spacing <- c(3.97, 4.43)
+# overall_length <- probe_head + width + contact_spacing * (n_contacts-1)
+diameter <- 0.8
+n_contacts <- 16
+
+for( spacing in contact_spacing ) {
+  contacts <- probe_head + width / 2 + 0:(n_contacts-1) * spacing
+  overall_length <- ceiling(max(contacts) + width / 2 + 300)
+  proto <- seeg_prototype(
+    type = sprintf("sEEG-PMT-2102-16-SP%.0f", spacing * 100),
+    description = c(
+      sprintf("PMT sEEG - %d contacts", n_contacts),
+      "Contact length   : 2    mm",
+      sprintf("Central spacing  : %.2f mm", spacing),
+      "Tip size         : 0    mm",
+      "Diameter         : 0.8  mm"
+    ),
+    center_position = contacts,
+    contact_widths = width,
+    diameter = diameter,
+    overall_length = overall_length,
+    default_interpolation = sprintf("%.2fx%d", spacing, n_contacts - 1L),
+    overwrite = TRUE
+  )
+}
+
+
+