--- title: Regions of Interest (ROI) date: '`r Sys.Date()`' output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Regions of Interest (ROI)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} params: family: red css: albers.css resource_files: - albers.css - albers.js includes: in_header: |- --- ```{r, echo = FALSE, message = FALSE} knitr::opts_chunk$set(collapse = T, comment = "#>") library(assertthat) library(purrr) library(neuroim2) options(mc.cores=1) ``` # Introduction to Regions of Interest ## What are ROIs? Regions of Interest (ROIs) are fundamental tools in neuroimaging analysis that allow researchers to focus on specific brain areas or patterns. Rather than analyzing every voxel in the brain independently, ROIs enable: - **Targeted analysis** of anatomically or functionally defined brain regions - **Dimensionality reduction** by aggregating signals within regions - **Statistical power** improvement through spatial averaging - **Hypothesis-driven** investigations of specific brain areas ## ROI Types in neuroim2 The **neuroim2** package provides comprehensive support for creating and manipulating ROIs: | ROI Type | Function | Description | Use Case | |----------|----------|-------------|----------| | **Spherical** | `spherical_roi()` | Sphere around a center point | Searchlight analysis, seed-based connectivity | | **Cuboid** | `cuboid_roi()` | Rectangular box region | Grid-based parcellation | | **Square** | `square_roi()` | 2D square in one plane | Slice-based analysis | | **Clustered** | `ClusteredNeuroVol` | Parcellation-based regions | Atlas-based analysis | | **Custom** | `Kernel` | Weighted regions | Gaussian-weighted, gradient-based | ## Quick Start ```{r} # Load the package library(neuroim2) # Read a brain volume file_name <- system.file("extdata", "global_mask2.nii.gz", package="neuroim2") vol <- read_vol(file_name) # Create a simple spherical ROI roi <- spherical_roi(vol, c(20, 20, 20), radius = 5) cat("ROI contains", length(roi), "voxels\n") ``` --- # Basic ROI Operations ## Creating a Spherical ROI Spherical ROIs are the most commonly used type in neuroimaging. They're particularly useful for searchlight analyses and seed-based connectivity studies. ```{r} # Load an example brain volume file_name <- system.file("extdata", "global_mask2.nii.gz", package="neuroim2") vol <- read_vol(file_name) # Create a spherical ROI centered at voxel coordinates [20,20,20] # with a 5mm radius, filling all values with 100 sphere <- spherical_roi(vol, c(20, 20, 20), radius = 5, fill = 100) # Examine the ROI cat("Number of voxels in ROI:", length(sphere), "\n") cat("ROI dimensions:", dim(sphere), "\n") cat("All values are 100:", all(sphere == 100), "\n") ``` ### Performance Note: C++ vs Pure R Implementation ```{r} # The spherical_roi function can use either C++ (fast) or pure R (slower) # C++ is the default and recommended for performance sphere_cpp <- spherical_roi(vol, c(20, 20, 20), radius = 5, use_cpp = TRUE) sphere_r <- spherical_roi(vol, c(20, 20, 20), radius = 5, use_cpp = FALSE) # They produce the same voxel set (up to ordering and type) all.equal( coords(sphere_cpp)[order(coords(sphere_cpp)[,1], coords(sphere_cpp)[,2], coords(sphere_cpp)[,3]),], coords(sphere_r)[order(coords(sphere_r)[,1], coords(sphere_r)[,2], coords(sphere_r)[,3]),], check.attributes = FALSE ) ``` ## Creating a Spherical ROI Around Real-World Coordinates Often, you'll have coordinates in millimeter space (e.g., from published studies) rather than voxel indices. ```{r} # Define a real-world coordinate in mm rpoint <- c(-34, -28, 10) # Convert real-world coordinates to voxel coordinates vox <- coord_to_grid(vol, rpoint) cat("Real coordinate:", rpoint, "-> Voxel coordinate:", vox, "\n") # Create spherical ROI with 10mm radius sphere <- spherical_roi(vol, vox, radius = 10, fill = 1) cat("ROI contains", length(sphere), "voxels\n") # Verify the center of mass is close to our target coords <- index_to_coord(vol, indices(sphere)) center_of_mass <- colMeans(coords) cat("Center of mass:", center_of_mass, "\n") cat("Original point:", rpoint, "\n") ``` ## Creating Multiple Spherical ROIs Efficiently When creating many ROIs, use the vectorized `spherical_roi_set()` function for better performance: ```{r} library(neuroim2) vol <- read_vol(system.file("extdata", "global_mask_v4.nii", package = "neuroim2")) d <- dim(vol) # Define multiple ROI centers centers <- rbind( c(floor(d[1]/3), floor(d[2]/3), floor(d[3]/3)), c(floor(d[1]/2), floor(d[2]/2), floor(d[3]/2)), c(floor(2*d[1]/3), floor(2*d[2]/3), floor(2*d[3]/3)) ) # Efficient vectorized creation rois <- spherical_roi_set(bvol = vol, centroids = centers, radius = 5, fill = 1) cat("Created", length(rois), "ROIs\n") # Compare with loop approach (slower but equivalent) rois2 <- lapply(seq_len(nrow(centers)), function(i) { spherical_roi(vol, centers[i,], 5, fill = 1) }) cat("Loop method also created", length(rois2), "ROIs\n") ``` ## Creating Cuboid and Square ROIs ### Cuboid ROIs (3D boxes) ```{r} # Create a cuboid ROI - a 3D rectangular box sp1 <- NeuroSpace(c(20, 20, 20), c(1, 1, 1)) # Create a 7x7x7 cube (surround=3 means 3 voxels on each side of center) cube <- cuboid_roi(sp1, c(10, 10, 10), surround = 3, fill = 5) cat("Cuboid ROI contains", length(cube), "voxels\n") cat("All values are 5:", all(cube == 5), "\n") # Get the coordinates vox_coords <- coords(cube) cat("Coordinate ranges:\n") cat(" X:", range(vox_coords[,1]), "\n") cat(" Y:", range(vox_coords[,2]), "\n") cat(" Z:", range(vox_coords[,3]), "\n") ``` ### Square ROIs (2D in 3D space) ```{r} # Create a square ROI - a 2D square fixed in one dimension sp1 <- NeuroSpace(c(20, 20, 20), c(1, 1, 1)) # Create a 5x5 square in the z=10 plane (fixdim=3 fixes the z dimension) square <- square_roi(sp1, c(10, 10, 10), surround = 2, fill = 3, fixdim = 3) cat("Square ROI contains", length(square), "voxels (should be 25)\n") # Verify it's planar vox_coords <- coords(square) cat("All z-coordinates are the same:", length(unique(vox_coords[,3])) == 1, "\n") cat("Z-plane:", unique(vox_coords[,3]), "\n") ``` --- # Working with ROI Data ## Converting ROIs to Sparse Volumes Sparse volumes are memory-efficient representations that only store non-zero values: ```{r} # Create a spherical ROI sphere <- spherical_roi(vol, c(50, 10, 10), radius = 10, fill = 1) cat("ROI contains", length(sphere), "voxels\n") # Convert to SparseNeuroVol - memory efficient storage # Prefer the provided coercion helper sparsevol <- as.sparse(sphere) # Verify properties cat("Sum of values match:", sum(sparsevol) == sum(sphere), "\n") cat("Dimensions match:", all(dim(sparsevol) == dim(vol)), "\n") # Memory comparison roi_size <- object.size(sphere) sparse_size <- object.size(sparsevol) cat("Memory usage - ROI:", format(roi_size, units = "auto"), "\n") cat("Memory usage - Sparse:", format(sparse_size, units = "auto"), "\n") ``` ## Extracting Time-Series from ROIs (4D Data) ROIs are particularly useful for extracting time-series from 4D functional data: ```{r} # Load a 4D dataset (using mask file as example - normally this would be fMRI data) vec4d <- read_vec(system.file("extdata", "global_mask_v4.nii", package = "neuroim2")) cat("4D data dimensions:", dim(vec4d), "\n") # Create an ROI roi <- spherical_roi(vol, c(12, 12, 12), radius = 6) # Extract time-series from the ROI roi_series <- series_roi(vec4d, roi) cat("ROI time-series dimensions:", dim(roi_series), "\n") # Get mean time-series across ROI (average across voxels) mat_roi <- values(roi_series) # T x N matrix mean_series <- rowMeans(mat_roi) # length equals time dimension cat("Mean time-series length:", length(mean_series), "\n") ``` ## Using ROIVec for 4D ROI Data The `ROIVec` class is designed for efficient storage of 4D ROI data: ```{r} # Create a 4D NeuroSpace vspace <- NeuroSpace(dim = c(10, 10, 10, 20), spacing = c(1, 1, 1)) # Define ROI coordinates roi_coords <- matrix(c(5,5,5, 6,5,5, 5,6,5, 5,5,6), ncol = 3, byrow = TRUE) # Create synthetic time-series data for each voxel n_timepoints <- 20 n_voxels <- nrow(roi_coords) roi_data <- matrix(rnorm(n_timepoints * n_voxels), nrow = n_timepoints, ncol = n_voxels) # Create ROIVec object roi_vec <- ROIVec(vspace, roi_coords, roi_data) cat("ROIVec created with", nrow(roi_coords), "voxels and", nrow(roi_data), "timepoints\n") # Access as matrix of values (T x N) roi_matrix <- values(roi_vec) cat("Matrix dimensions (T x N):", dim(roi_matrix), "\n") ``` --- # Searchlight Analyses Searchlight analysis is a powerful technique for multivariate pattern analysis that examines local neighborhoods throughout the brain. ## Basic Searchlight ```{r} library(purrr) # Generate exhaustive searchlight covering all voxels slist <- searchlight(vol, eager = TRUE, radius = 8) cat("Number of searchlights:", length(slist), "\n") # Compute mean value in each searchlight searchlight_means <- slist %>% purrr::map_dbl(~ mean(vol[coords(.)])) cat("Computed means for", length(searchlight_means), "searchlights\n") cat("Mean range:", range(searchlight_means, na.rm = TRUE), "\n") ``` ## Randomized Searchlight Randomized searchlight samples voxels without replacement, ensuring coverage while reducing redundancy: ```{r} # Randomized searchlight - each voxel appears in at least one searchlight set.seed(42) # For reproducibility random_lights <- vol %>% random_searchlight(radius = 8) %>% purrr::map_dbl(~ mean(vol[coords(.)])) cat("Random searchlight count:", length(random_lights), "\n") ``` ## Clustered Searchlight Use anatomical or functional parcellations to define ROIs: ```{r} # Create a clustering over the voxel space grid <- index_to_coord(vol, which(vol > 0)) set.seed(123) kres <- kmeans(grid, centers = 50, iter.max = 500) # Create ClusteredNeuroVol kvol <- ClusteredNeuroVol(vol, kres$cluster) cat("Created", length(unique(kres$cluster)), "clusters\n") # Run clustered searchlight cluster_means <- vol %>% clustered_searchlight(cvol = kvol) %>% purrr::map_dbl(~ mean(vol[coords(.)])) cat("Cluster mean range:", range(cluster_means, na.rm = TRUE), "\n") ``` --- # Image Patches Patches are regular, fixed-size regions that tile the image space, useful for convolutional approaches: ```{r} # Create 3x3x1 patches covering the volume pset <- patch_set(vol, dims = c(3, 3, 1)) cat("Number of patches:", length(pset), "\n") # Compute statistics for each patch patch_means <- pset %>% purrr::map_dbl(~ mean(.)) cat("Patch mean range:", range(patch_means, na.rm = TRUE), "\n") # Restrict patches to masked regions only pset_masked <- patch_set(vol, dims = c(3, 3, 1), mask = as.logical(vol)) cat("Masked patches:", length(pset_masked), "\n") # Patches are guaranteed to be equal size (padded at edges if needed) patch_sizes <- pset_masked %>% purrr::map_int(~ length(.)) cat("All patches same size:", length(unique(patch_sizes)) == 1, "\n") cat("Patch size:", unique(patch_sizes), "voxels\n") ``` --- # Advanced ROI Techniques ## Custom Weighted ROIs with Kernels Create Gaussian-weighted or other custom-shaped ROIs: ```{r} # Create a Gaussian kernel for weighted ROI kern_dim <- c(5, 5, 5) # 5x5x5 kernel voxel_dim <- c(1, 1, 1) # 1mm isotropic voxels # Create Gaussian-weighted kernel gauss_kernel <- Kernel(kerndim = kern_dim, vdim = voxel_dim, FUN = dnorm, sd = 1.5) # Embed kernel at a specific location embedded <- embed_kernel(gauss_kernel, space(vol), center_voxel = c(20, 20, 20)) cat("Kernel weights sum to:", sum(embedded), "\n") ``` ## Working with ClusteredNeuroVec for Parcellated Data Efficiently handle parcellated 4D data where voxels are grouped into regions: ```{r} # Create synthetic 4D data sp4 <- NeuroSpace(c(10, 10, 10, 20), c(1, 1, 1)) arr <- array(rnorm(10*10*10*20), dim = c(10, 10, 10, 20)) vec <- NeuroVec(arr, sp4) # Create mask for central region sp3 <- NeuroSpace(c(10, 10, 10), c(1, 1, 1)) mask_arr <- array(FALSE, dim = c(10, 10, 10)) mask_arr[3:8, 3:8, 3:8] <- TRUE mask <- LogicalNeuroVol(mask_arr, sp3) # Assign voxels to 5 random clusters n_voxels <- sum(mask_arr) clusters <- sample(1:5, n_voxels, replace = TRUE) cvol <- ClusteredNeuroVol(mask, clusters) # Create clustered representation - one time-series per cluster cv <- ClusteredNeuroVec(vec, cvol) # Access cluster time-series efficiently cluster_matrix <- as.matrix(cv) # T x K matrix cat("Cluster matrix dimensions:", dim(cluster_matrix), "\n") cat("(", dim(cluster_matrix)[1], "timepoints x", dim(cluster_matrix)[2], "clusters)\n") ``` ## ROI Set Operations Combine and manipulate multiple ROIs: ```{r} # Create two overlapping spherical ROIs roi1 <- spherical_roi(vol, c(20, 20, 20), radius = 6, fill = 1) roi2 <- spherical_roi(vol, c(23, 20, 20), radius = 6, fill = 1) # Get indices for set operations idx1 <- indices(roi1) idx2 <- indices(roi2) # Intersection - voxels in both ROIs intersection_idx <- intersect(idx1, idx2) cat("Intersection:", length(intersection_idx), "voxels\n") # Union - voxels in either ROI union_idx <- union(idx1, idx2) cat("Union:", length(union_idx), "voxels\n") # Difference - voxels in roi1 but not roi2 diff_idx <- setdiff(idx1, idx2) cat("Difference (roi1 - roi2):", length(diff_idx), "voxels\n") # Calculate overlap percentage overlap_pct <- length(intersection_idx) / length(union_idx) * 100 cat("Overlap:", round(overlap_pct, 1), "%\n") ``` --- # Best Practices and Performance ## Memory Management When working with many ROIs, consider memory usage: ```{r} # Compare memory usage of different ROI storage methods n_rois <- 100 centers <- matrix(runif(n_rois * 3, 5, 15), ncol = 3) # Method 1: List of ROI objects (flexible but more memory) roi_list <- lapply(1:nrow(centers), function(i) { spherical_roi(vol, centers[i,], radius = 6, fill = 1) }) # Method 2: Single sparse volume with labeled regions (ensure increasing indices) all_indices <- list() all_labels <- list() for (i in 1:nrow(centers)) { roi <- spherical_roi(vol, centers[i,], radius = 6) idx <- indices(roi) idx <- sort(unique(idx)) all_indices[[i]] <- idx all_labels[[i]] <- rep(i, length(idx)) } idx_all <- unlist(all_indices) lab_all <- unlist(all_labels) ord <- order(idx_all) idx_all <- idx_all[ord] lab_all <- lab_all[ord] # Ensure strictly increasing indices by removing duplicates keep <- !duplicated(idx_all) idx_all <- idx_all[keep] lab_all <- lab_all[keep] combined_sparse <- SparseNeuroVol(lab_all, space(vol), indices = idx_all) cat("Memory usage:\n") cat(" ROI list:", format(object.size(roi_list), units = "auto"), "\n") cat(" Sparse combined:", format(object.size(combined_sparse), units = "auto"), "\n") ``` ## Choosing ROI Sizes Guidelines for selecting appropriate ROI sizes: ```{r} # Demonstrate effect of ROI size on coverage and overlap radii <- c(6, 8, 10, 12) coverage_stats <- data.frame( radius = radii, n_voxels = numeric(length(radii)), pct_overlap = numeric(length(radii)) ) center1 <- c(20, 20, 20) center2 <- c(25, 20, 20) # 5 voxels apart for (i in seq_along(radii)) { roi1 <- spherical_roi(vol, center1, radius = radii[i]) roi2 <- spherical_roi(vol, center2, radius = radii[i]) coverage_stats$n_voxels[i] <- length(roi1) overlap <- length(intersect(indices(roi1), indices(roi2))) total <- length(union(indices(roi1), indices(roi2))) coverage_stats$pct_overlap[i] <- overlap / total * 100 } print(coverage_stats) ``` ## Parallel Processing Tips For large-scale ROI analyses, consider parallel processing: ```{r, eval=FALSE} # Example of parallel searchlight analysis (not run in vignette) library(parallel) library(foreach) library(doParallel) # Setup parallel backend n_cores <- detectCores() - 1 cl <- makeCluster(n_cores) registerDoParallel(cl) # Parallel searchlight computation searchlight_results <- foreach(roi = searchlight_list, .packages = c("neuroim2")) %dopar% { # Your analysis function here mean(data[coords(roi)]) } # Clean up stopCluster(cl) ``` --- # Troubleshooting and Tips ## Common Issues and Solutions ### Issue: ROI extends outside brain mask ```{r} # Problem: ROI at edge of brain edge_vox <- c(2, 2, 2) # Near edge of volume # Solution 1: Use nonzero=TRUE to keep only mask voxels roi_filtered <- spherical_roi(vol, edge_vox, radius = 5, nonzero = TRUE) cat("Filtered ROI size:", length(roi_filtered), "voxels\n") # Solution 2: Check if ROI is valid before analysis if (length(roi_filtered) < 10) { cat("Warning: ROI too small for reliable analysis\n") } ``` ### Issue: Different coordinate systems ```{r} # Always verify your coordinate system test_coord_mm <- c(-34, -28, 10) # MNI coordinates test_coord_vox <- coord_to_grid(vol, test_coord_mm) # Verify round-trip conversion back_to_mm <- grid_to_coord(vol, matrix(test_coord_vox, nrow = 1)) cat("Original (mm):", test_coord_mm, "\n") cat("Voxel:", test_coord_vox, "\n") cat("Back to mm:", back_to_mm, "\n") ``` --- # See Also For related functionality, see these other vignettes and help pages: - `vignette("NeuroVector")` - Working with 4D neuroimaging data - `?read_vec` and `?read_vol` - Reading neuroimaging files - `?NeuroSpace` - Understanding coordinate systems and spaces - `?ClusteredNeuroVol` - Parcellation-based analyses - `?SparseNeuroVol` - Memory-efficient sparse representations - `?series` and `?series_roi` - Time-series extraction For the complete list of ROI-related functions: ```{r, eval=FALSE} help(package = "neuroim2", topic = "roi") ```