# Install from local source
devtools::install()

# Or build and install
R CMD build .
R CMD INSTALL neurocluster_*.tar.gz

library(neurocluster)
library(neuroim2)

# Load your data
mask <- read_vol("brain_mask.nii")
vec <- read_vec("fmri_data.nii.gz")

# Cluster with default method (supervoxels)
result <- cluster4d(vec, mask, n_clusters = 100)

# 3D convenience: you can pass a single-volume NeuroVol directly;
# it will be wrapped to a one-frame NeuroVec internally
vol <- read_vol("anatomical.nii.gz")
vol_mask <- vol > 0
result_3d <- cluster4d(vol, vol_mask, n_clusters = 200, method = "slic")

# Try different methods
result_snic <- cluster4d(vec, mask, n_clusters = 100, method = "snic")
result_msf <- cluster4d(vec, mask, n_clusters = 100, method = "slice_msf")

# Load test data
mask <- NeuroVol(array(1, c(64,64,32)), NeuroSpace(c(64,64,32)))
vec <- replicate(100, NeuroVol(array(rnorm(64*64*32), c(64,64,32)),
                 NeuroSpace(c(64,64,32))), simplify=FALSE)
vec <- do.call(concat, vec)

# Cluster with balanced spatial/feature weighting
result <- cluster4d(vec, mask, 
                   n_clusters = 50,
                   spatial_weight = 0.5)

# Access results
clusters <- result$cluster          # Cluster assignments
centers <- result$centers            # Feature centers
spatial_centers <- result$coord_centers  # Spatial centroids

# Get intelligent parameter suggestions
n_voxels <- sum(mask > 0)
n_timepoints <- dim(vec)[4]

# Get balanced recommendations
params <- suggest_cluster4d_params(n_voxels, n_timepoints, 
                                   priority = "balanced")

# Use recommended settings
result <- cluster4d(vec, mask, 
                   n_clusters = params$n_clusters,
                   method = params$recommended_method)

# Supervoxels with custom kernel bandwidths
result <- cluster4d(vec, mask, n_clusters = 100,
                   method = "supervoxels",
                   sigma1 = 2.0,    # Feature kernel
                   sigma2 = 3.0,    # Spatial kernel
                   use_gradient = TRUE)

# Slice-MSF with consensus clustering
result <- cluster4d(vec, mask, n_clusters = 100,
                   method = "slice_msf",
                   num_runs = 5,      # Multiple runs
                   consensus = TRUE,  # Consensus clustering
                   stitch_z = TRUE)   # Reduce z-plane artifacts

# FLASH-3D with temporal smoothing
result <- cluster4d(vec, mask, n_clusters = 100,
                   method = "flash3d",
                   lambda_t = 1.5,    # Temporal regularization
                   bits = 128)        # Quantization precision

# Run multiple methods
methods <- c("supervoxels", "snic", "slice_msf")
results <- lapply(methods, function(m) {
  cluster4d(vec, mask, n_clusters = 50, method = m)
})
names(results) <- methods

# Compare results
comparison <- compare_cluster4d(results)
print(comparison)

# Use gradient-based initialization for better seeds
result <- cluster4d(vec, mask, n_clusters = 100,
                   method = "supervoxels",
                   use_gradient = TRUE)

# Extract ClusteredNeuroVol
clusvol <- result$clusvol

# Get cluster sizes
cluster_sizes <- table(result$cluster)

# Find cluster for specific voxel
voxel_coord <- c(32, 32, 16)
voxel_idx <- coord_to_index(mask, voxel_coord)
voxel_cluster <- result$cluster[voxel_idx]

# Extract time series for a cluster
cluster_id <- 5
cluster_mask <- result$cluster == cluster_id
cluster_timeseries <- series(vec, which(cluster_mask))

# Plot cluster centers
plot(result)

# Visualize spatial distribution
library(neuroim2)
writeVol(result$clusvol, "clusters.nii.gz")

# View specific cluster
cluster_vol <- NeuroVol(array(0, dim(mask)), space(mask))
cluster_vol[result$cluster == 5] <- 1
writeVol(cluster_vol, "cluster_5.nii.gz")

# Original supervoxels syntax
result <- supervoxels(vec, mask, K = 100, alpha = 0.3)

# Original SNIC syntax  
result <- snic(vec, mask, K = 100, compactness = 5)

# Original slice_msf syntax
result <- slice_msf(vec, mask, target_k_global = 100)