params <- list(family = "red", preset = "homage") ## ----echo = FALSE, message = FALSE-------------------------------------------- knitr::opts_chunk$set(collapse = T, comment = "#>", message = FALSE, warning = FALSE) suppressPackageStartupMessages({ library(purrr) library(neuroim2) }) options(mc.cores=1) ## ----------------------------------------------------------------------------- file_name <- system.file("extdata", "global_mask2.nii.gz", package="neuroim2") vol <- read_vol(file_name) roi <- spherical_roi(vol, c(20, 20, 20), radius = 5) cat("ROI contains", length(roi), "voxels\n") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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 ) ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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) ## ----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) ## ----------------------------------------------------------------------------- # 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") } ## ----------------------------------------------------------------------------- # 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") ## ----eval=FALSE--------------------------------------------------------------- # help(package = "neuroim2", topic = "roi")