params <-
list(family = "red")

## ----echo = FALSE, message = FALSE--------------------------------------------
knitr::opts_chunk$set(collapse = T, comment = "#>")
library(assertthat)
library(purrr)
library(neuroim2)
options(mc.cores=1)

## -----------------------------------------------------------------------------
# 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")

## -----------------------------------------------------------------------------
# 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")

## -----------------------------------------------------------------------------
# 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")

## -----------------------------------------------------------------------------
# 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")

## -----------------------------------------------------------------------------
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 - 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")

## -----------------------------------------------------------------------------
# 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")

## -----------------------------------------------------------------------------
# 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")

