Package 'neurocluster'

Description

Clusters fMRI voxels into spatially-coherent groups based on temporal correlation and spatial proximity. Includes optional refinement for boundary corrections. The algorithm supports parallel processing via the future package for improved performance.

Usage

acsc(
  bvec,
  mask,
  block_size = 2,
  ann_k = 10,
  alpha = 0.5,
  correlation_metric = c("pearson", "spearman", "robust"),
  spatial_weighting = c("gaussian", "binary"),
  refine = TRUE,
  max_refine_iter = 5,
  K = NULL,
  knn_proj_dim = 0L,
  knn_proj_method = c("none", "dct", "rp"),
  knn_proj_seed = 1L
)

Arguments

Details

ACSC builds feature-neighbor edges from the declared correlation metric and unions them with six-neighbor block-grid edges. Edge correlations are mapped to non-negative Louvain similarities as (correlation + 1) / 2. At alpha = 0, feature-neighbor selection is bypassed entirely, so the graph topology and weights depend only on spatial adjacency. Optional DCT or random projection changes feature-neighbor candidate selection, but edge weights are always recomputed using the declared metric. Louvain calls use a fixed internal seed and restore the caller's RNG state, so identical ACSC inputs produce the same partition independently of the surrounding random stream.

Boundary refinement uses the same metric representation as graph building. Pearson uses centered unit vectors, Spearman uses centered average ranks, and robust mode uses median-centered Tukey-bisquare residuals. The C++ path accelerates centroid dot products and recomputes the active boundary after every synchronous pass; failures fall back to the equivalent R path.

After refinement, labels are flood-filled into connected components. If K is supplied, the shared adjacency-preserving exact-K engine merges or splits those components and fails closed when the target is topologically infeasible.

Value

A standardized cluster4d_result (also classed acsc_result) with final contiguous labels, original-space centers, physical coordinate centers, provenance, and these ACSC-specific elements:

cluster_map

3D array with cluster labels per voxel.

graph

An igraph object used for clustering.

init_block_label

Initial coarse partition (3D array) matching mask dimensions.

Description

Diagnostic function to detect whether cluster IDs are spatially ordered along any axis. Useful for identifying clustering algorithms that may produce ordered IDs.

Usage

check_cluster_ordering(cluster_result)

Arguments

Value

A list with correlation statistics:

• cor_x: Spearman correlation between cluster ID and x-coordinate

• cor_y: Spearman correlation between cluster ID and y-coordinate

• cor_z: Spearman correlation between cluster ID and z-coordinate

• max_abs_cor: Maximum absolute correlation across all axes

• is_ordered: Logical indicating if any axis shows strong ordering (|r| > 0.3)

Examples

## Not run: 
result <- cluster4d(vec, mask, K = 100, method = "flash3d")

# Check for ordering
ordering <- check_cluster_ordering(result)
print(ordering)

# If strongly ordered, randomize
if (ordering$is_ordered) {
  result <- randomize_cluster_ids(result)
}

## End(Not run)

Description

This function extracts the cluster class identifiers from a cluster result object. It is a method for the cl_class_ids generic from the clue package.

Usage

## S3 method for class 'cluster_result'
cl_class_ids(x)

Arguments

Value

An integer vector of cluster assignments, one for each data point.

See Also

cl_class_ids for the generic function.

Description

Computes partition, within-cluster temporal, and physical spatial estimands from final voxel labels. Temporal coherence is the mean Pearson correlation over all unordered within-cluster voxel pairs. Spatial dispersion is the root mean squared Euclidean distance from each voxel coordinate to its final cluster centroid; lower values are more compact.

Usage

cluster_metrics(result, feature_mat = NULL, coords = NULL, truth = NULL)

Arguments

Value

A named list. Always contains n_clusters and size_summary. With truth, it adds adjusted Rand, variation of information in bits, and pairwise Dice. With feature_mat, it adds temporal_pairwise_correlation. With coordinates, it adds spatial_rms_distance_mm (or spatial_rms_distance when explicit coordinates are not declared to be millimetres by a canonical result).

Description

Performs spatially-constrained clustering on 4D neuroimaging data using various algorithms. This is the main entry point for all clustering methods in the neurocluster package.

Usage

cluster4d(
  vec,
  mask,
  n_clusters = 100,
  method = c("supervoxels", "snic", "slic", "corr_slic", "brs_slic", "slice_msf",
    "flash3d", "g3s", "rena", "rena_plus", "mcl", "acsc", "commute"),
  spatial_weight = 0.5,
  max_iterations = NULL,
  connectivity = NULL,
  parallel = TRUE,
  verbose = FALSE,
  ...
)

Arguments

Value

A cluster4d_result object (also inherits from cluster_result) containing:

Algorithm Comparison

The Speed column reports single-threaded wall-clock tiers measured with cluster4d() defaults on the neuroim2 bundled mask (global_mask_v4.nii, 29,532 voxels), 20 timepoints of Gaussian noise and n_clusters = 200: Fast is under 1 s, Moderate a few seconds, and Slow tens of seconds. Timings depend on hardware, mask size, and parameters; treat the tiers as relative guidance, not a ranking of quality. Parallel, Iterations, and Spatial weight state whether the common parallel, max_iterations, and spatial_weight arguments are active.

Parameter Guidelines

For whole-brain parcellation:

• n_clusters: 100-1000 depending on desired granularity

• spatial_weight: 0.4-0.6 for balanced clustering

• connectivity: 26 for smoother boundaries

For ROI analysis:

• n_clusters: 10-100 depending on ROI size

• spatial_weight: 0.2-0.4 to emphasize functional similarity

• connectivity: 6 for more discrete parcels

For high-resolution data (< 2mm):

• method: "flash3d" for speed; Slice-MSF remains experimental

• n_clusters: Scale with voxel count (roughly n_voxels/200)

See Also

Method-specific functions: cluster4d_supervoxels, cluster4d_snic, cluster4d_slic, cluster4d_slice_msf, cluster4d_flash3d, cluster4d_g3s, cluster4d_rena, cluster4d_rena_plus, cluster4d_mcl, cluster4d_acsc, cluster4d_commute

Legacy functions (deprecated): supervoxels, snic, slic4d_supervoxels, slice_msf, supervoxels_flash3d

Examples

# Simple synthetic example (runs quickly)
library(neuroim2)
set.seed(1)
mask <- NeuroVol(array(1, c(4,4,4)), NeuroSpace(c(4,4,4)))
vec <- NeuroVec(array(rnorm(4*4*4*20), c(4,4,4,20)),
                NeuroSpace(c(4,4,4,20)))
result <- cluster4d(vec, mask, n_clusters = 3, method = "g3s",
                   max_iterations = 1)
print(result$n_clusters)


# More realistic examples with larger data
mask <- NeuroVol(array(1, c(20,20,20)), NeuroSpace(c(20,20,20)))
vec <- replicate(50, NeuroVol(array(runif(20*20*20), c(20,20,20)),
                              NeuroSpace(c(20,20,20))), simplify=FALSE)
vec <- do.call(concat, vec)

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

# Fast clustering with FLASH-3D (hash-based)
result <- cluster4d(vec, mask, n_clusters = 100, method = "flash3d")

# Emphasize spatial compactness
result <- cluster4d(vec, mask, n_clusters = 100, spatial_weight = 0.8)

# Explicitly evaluate the experimental Slice-MSF method
result <- cluster4d(vec, mask, n_clusters = 100, 
                   method = "slice_msf",
                   num_runs = 3,  # slice_msf-specific parameter
                   consensus = TRUE)

# Get parameter suggestions for your data
n_vox <- sum(mask > 0)
n_time <- dim(vec)[4]
params <- suggest_cluster4d_params(n_vox, n_time, priority = "quality")
result <- cluster4d(vec, mask, 
                   n_clusters = params$n_clusters,
                   method = params$recommended_method)

Description

Wrapper for Adaptive Correlation Superclustering algorithm with standardized interface. ACSC uses graph-based clustering with Louvain community detection and optional boundary refinement.

Usage

cluster4d_acsc(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  max_iterations = 5,
  verbose = FALSE,
  block_size = 2,
  correlation_metric = c("pearson", "spearman", "robust"),
  spatial_weighting = c("gaussian", "binary"),
  refine = TRUE,
  max_refine_iter = NULL,
  ...
)

Arguments

Value

A cluster4d_result object

Description

Hybrid method: fast sketch-SLIC coarse assignment, then exact-correlation refinement only on boundary voxels.

Usage

cluster4d_brsslic(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.05,
  max_iterations = 2,
  connectivity = 6,
  parallel = FALSE,
  verbose = FALSE,
  embedding_dim = 24L,
  sketch_repeats = 1L,
  coarse_alpha = NULL,
  boundary_passes = 1L,
  global_passes = 0L,
  refine_spatial_weight = NULL,
  refine_l2_weight = 0,
  refine_stride = NULL,
  seed = 1L,
  min_size = NULL,
  ...
)

Arguments

Details

Both the coarse and exact-refinement spatial weights use the convex blend (1 - w) * correlation_distance + w * scaled_spatial_distance. refine_l2_weight is a separate optional feature penalty, so a strictly spatial-only refinement uses refine_spatial_weight = 1 and refine_l2_weight = 0. The hash-sketch artifact is explicitly marked as a pre-refinement embedding under metadata$approximation$coarse_embedding. Public centers and coordinates are final-label means in original feature and physical coordinate spaces, respectively.

Value

A cluster4d_result object

Description

Wrapper for commute_cluster with standardized interface.

Usage

cluster4d_commute(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  verbose = FALSE,
  ncomp = NULL,
  ...
)

Arguments

Value

A cluster4d_result object

Note

Masks with more than 3,000 voxels use a sparse shift-invert eigensolver (the required package eigencore); smaller masks use the dense eigen() path. Pass eigen_solver (see commute_cluster) to choose explicitly.

Description

Wrapper for a correlation-first SLIC variant that builds a compact random projection embedding of voxel time series, then runs local 3D SLIC updates.

Usage

cluster4d_corrslic(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  max_iterations = 6,
  connectivity = 6,
  parallel = FALSE,
  verbose = FALSE,
  embedding_dim = 64L,
  adaptive_embedding = FALSE,
  embedding_basis = c("hash", "dct"),
  embedding_whiten = FALSE,
  sketch_repeats = 1L,
  alpha = NULL,
  assign_stride = 1L,
  quantize_assign = FALSE,
  refine_exact_iters = 0L,
  refine_boundary_only = TRUE,
  refine_stride = 1L,
  refine_alpha = NULL,
  seed = 1L,
  min_size = NULL,
  ...,
  .input_contract = NULL
)

Arguments

Details

spatial_weight (or its alias alpha) is a convex blend: (1 - w) * correlation_distance + w * scaled_spatial_distance. Coarse correlation is approximated by a hash or DCT sketch. Assignment may additionally use int8 quantization or z-plane subsampling; exact refinement may use a temporal stride. These choices and their stage are recorded under metadata$approximation. The public centers are always recomputed from the final labels in the original time-series space; sketch centers are separately typed metadata and are never substituted for them.

Value

A cluster4d_result object

Description

Wrapper for FLASH-3D algorithm with standardized interface.

Usage

cluster4d_flash3d(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.25,
  max_iterations = 4,
  verbose = FALSE,
  lambda_t = NULL,
  bits = 64,
  dctM = 16,
  ...
)

Arguments

Value

A cluster4d_result object

Description

Compresses voxel features, selects low-gradient seeds, and assigns voxels by deterministic multi-source Dijkstra propagation. Edges are exactly the requested masked-grid neighbors; excluded gaps are never bridged. The edge objective is alpha * feature_distance + (1 - alpha) * physical_distance / compactness.

Usage

cluster4d_g3s(
  vec,
  mask,
  K = 100,
  n_components = 15,
  variance_threshold = 0.95,
  alpha = 0.5,
  compactness = NULL,
  max_refinement_iter = 3,
  verbose = FALSE,
  use_irlba = TRUE,
  use_rsvd = TRUE,
  connectivity = 26
)

Arguments

Value

A typed g3s_result and cluster4d_result with K by T centers, physical coordinate centers, graph provenance, and physical-scale metadata.

References

Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1, 269-271.

Description

Builds a sparse voxel graph, runs an MCL flow process, and maps attractors to cluster labels. The implementation is optimized for sparse neuroimaging graphs.

Usage

cluster4d_mcl(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.2,
  max_iterations = 8,
  connectivity = 6,
  verbose = FALSE,
  inflation = 1.6,
  expansion = 2L,
  loop_weight = 1,
  prune_k = NULL,
  prune_threshold = 1e-06,
  tol = 1e-04,
  feature_metric = c("correlation", "euclidean"),
  feature_sigma = NULL,
  spatial_sigma = NULL,
  exact_k = FALSE
)

Arguments

Details

Every expansion, inflation, and pruning step is checked to remain finite, non-negative, and column stochastic. The parallel argument was removed because this sparse Matrix implementation has no parallel backend. Supplying it is therefore an ordinary unused-argument error rather than a silently sequential execution.

Value

A cluster4d_result. Metadata records natural and final K, convergence, selected inflation, candidate search results, flow sparsity, final normalization error, and exact-K policy.

Description

Wrapper for ReNA algorithm with standardized cluster4d interface.

Usage

cluster4d_rena(
  vec,
  mask,
  n_clusters = 100,
  max_iterations = 50,
  connectivity = 26,
  verbose = FALSE,
  exact_k = TRUE,
  ...
)

Arguments

Value

A cluster4d_result object

Description

Cluster4d wrapper for ReNA++

Usage

cluster4d_rena_plus(
  vec,
  mask,
  n_clusters = 100,
  r = 2,
  lambda = 1,
  grad_img = NULL,
  connectivity = 26,
  max_iterations = 50,
  verbose = FALSE,
  ...
)

Arguments

Value

A cluster4d_result.

Description

Wrapper for SLIC algorithm with standardized interface.

Usage

cluster4d_slic(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  max_iterations = 10,
  connectivity = 26,
  parallel = TRUE,
  verbose = FALSE,
  preserve_k = FALSE,
  seed_relocate = "none",
  ...,
  .input_contract = NULL
)

Arguments

Value

A cluster4d_result object

Description

Wrapper for slice_msf algorithm with standardized interface. Slice-MSF is experimental and is deliberately excluded from suggest_cluster4d_params() recommendations pending release recertification.

Usage

cluster4d_slice_msf(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.2,
  connectivity = 26,
  parallel = TRUE,
  verbose = FALSE,
  num_runs = 1,
  consensus = FALSE,
  stitch_z = TRUE,
  r = 12,
  gamma = 0,
  z_mult = 0,
  min_size = NULL,
  seed = 1L,
  ensemble_fraction = 0.8,
  ...
)

Arguments

Details

Unlike direct slice_msf(), this wrapper always requests exactly n_clusters, defaults to one non-consensus run, and chooses a conservative automatic min_size. spatial_weight maps to compactness = spatial_weight * 10 and then to FH scale 2 / (compactness + 1); it is not a convex spatial/feature blend. Positive gamma values fail with class slice_msf_unsupported_gamma.

Value

A cluster4d_result object with natural/exact-K and effective-control provenance in metadata$exact_k_repair; see slice_msf().

Description

Wrapper for SNIC algorithm with standardized interface.

Usage

cluster4d_snic(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  max_iterations = 100,
  verbose = FALSE,
  ...
)

Arguments

Value

A cluster4d_result object

Description

Wrapper for supervoxels algorithm with standardized interface.

Usage

cluster4d_supervoxels(
  vec,
  mask,
  n_clusters = 100,
  spatial_weight = 0.5,
  max_iterations = 50,
  connectivity = 27,
  parallel = TRUE,
  verbose = FALSE,
  sigma1 = 1,
  sigma2 = 2.5,
  use_gradient = TRUE,
  converge_thresh = 0.001,
  ...
)

Arguments

Value

A cluster4d_result object

Description

Compares predicted cluster assignments against ground truth using multiple metrics that are invariant to cluster label permutation.

Usage

clustering_accuracy(predicted, truth)

Arguments

Details

The metrics compare the two induced equivalence relations; cluster names do not need to match. All observations are included, including a label of zero. Missing or non-finite labels are rejected rather than silently removed. Computation uses only observed contingency cells and is linear in the number of observations for singleton-heavy partitions.

Adjusted Rand Index (ARI): Measures agreement between two clusterings, adjusted for chance. ARI = 1 means perfect agreement, ARI = 0 means random, ARI < 0 means worse than random.

Variation of information (VI): The sum of the two conditional entropies, measured with base-2 logarithms. VI is zero exactly for identical partitions and increases with disagreement.

Pairwise Dice: 2 * TP / (2 * TP + FP + FN) over unordered pairs of observations, where a positive pair is co-clustered. If both partitions contain only singleton clusters, their pairwise Dice is defined as one.

Matched accuracy: The greatest fraction of matched observations under a one-to-one mapping between predicted and truth labels. It is solved exactly on the sparse observed contingency graph. accuracy is retained as an alias for existing callers; it is not the former greedy approximation.

Value

A list with adjusted Rand index (ari), variation of information in bits (variation_of_information_bits), pairwise Dice (pairwise_dice), arithmetic normalized mutual information (nmi), exact maximum one-to-one label-matched accuracy (matched_accuracy), its compatibility alias accuracy, partition counts, and a permutation-invariant perfect flag.

Examples

# Perfect clustering
truth <- rep(1:3, each = 10)
pred <- rep(c(2, 3, 1), each = 10)  # Same clustering, different labels
clustering_accuracy(pred, truth)  # ARI = 1

Description

Builds an explicitly declared symmetric k-nearest-neighbor graph in physical coordinate space, computes a Laplacian commute-time embedding, and applies deterministically initialized k-means.

Usage

commute_cluster(
  bvec,
  mask,
  K = 100,
  ncomp = ceiling(sqrt(K * 2)),
  alpha = 0.5,
  sigma1 = 0.73,
  sigma2 = 5,
  connectivity = 6L,
  weight_mode = c("binary", "heat"),
  bad_voxel_policy = c("error", "zero_constant"),
  variance_tol = 1e-09,
  eigen_tol = 1e-10,
  eigen_solver = c("auto", "dense", "sparse"),
  verbose = FALSE
)

Arguments

Details

The implementation is deterministic. RNG initialization for the sparse eigensolver and stats::kmeans() is scoped to a fixed local seed, after which the caller's prior RNG state (including absence of .Random.seed) is restored exactly. With the dense solver it forms a full Laplacian eigendecomposition, requiring quadratic memory and cubic time, so eigen_solver = "auto" switches to the sparse solver above 3,000 voxels. Sparse solver memory and time depend on factorization fill-in, the requested embedding dimension, and convergence. Sparse eigenpairs must have original-Laplacian relative residuals at most 1e-10 and orthogonality error at most 1e-8. Residuals are scaled by (norm(L, "I") + abs(lambda)) * norm(v, "2"). These checks and eigencore's separate certificate at the solver tolerance are retained in metadata$spectral$diagnostics. The sparse solver never silently falls back to a dense solve. A full-spectrum request (ncomp = N - 1) uses the dense solver even when eigen_solver = "sparse".

Value

A commute_time_cluster_result with contiguous labels, K-by-T original-space feature means, K-by-3 physical coordinate means, the N-by-ncomp embedding, graph and spectral provenance, and explicit data-policy provenance.

Description

Compares clustering results over exactly the same included mask voxels and physical coordinate system. Ambiguous center-based scores are not used.

Usage

compare_cluster4d(
  ...,
  metrics = c("summary", "spatial_dispersion"),
  feature_mat = NULL
)

Arguments

Details

Every result must carry canonical mask and affine provenance. The exact included voxel indices and physical geometry are checked before any cross-result comparison. Partition metrics use only observed contingency cells and never construct an N-by-N co-membership matrix. Pairwise agreement columns are repeated on the two result rows so the return value remains a data frame.

Value

A data frame with one row per result and explicitly named metrics.

Description

Reduces the dimensionality of fMRI time series data using Singular Value Decomposition (SVD) while preserving the majority of variance. This is a key component of the G3S (Gradient-Guided Geodesic Supervoxels) algorithm, enabling 10-20x speedup in similarity computations.

Usage

compress_features_svd(
  feature_mat,
  n_components = 15,
  variance_threshold = 0.95,
  use_irlba = TRUE,
  use_rsvd = TRUE
)

Arguments

Details

Algorithm

The compression follows these steps:

1. Center and scale: Each timepoint is centered (mean=0) and scaled (sd=1) to ensure equal contribution from all timepoints.

2. SVD decomposition: Computes X = U * D * t(V) where:

   • U (N x k): Left singular vectors (voxel loadings)

   • D (k x k): Diagonal matrix of singular values

   • V (T x k): Right singular vectors (time loadings)

3. Compression: Retains only the first M components where M is chosen to balance speed (fewer components) and accuracy (more variance explained).

4. Normalization: Each compressed feature vector is normalized to unit length, enabling fast cosine similarity via dot products: cor(x, y) approx x * y

Performance Characteristics

• Memory: Compressed data uses M/T of original size (e.g., 15/300 = 5%)

• Speed: Similarity calculations are M/T times faster (e.g., 20x for 15/300)

• Accuracy: Typically preserves 95%+ of signal with M=10-20 for fMRI data

Choosing n_components

• Small (M=10): Maximum speed, good for noisy or heavily smoothed data

• Medium (M=15): Balanced, recommended default for most fMRI data

• Large (M=25): High accuracy, use for high SNR data or critical applications

The function will automatically increase n_components if needed to meet variance_threshold. If that requires the full available rank, the exact full-rank decomposition is used.

Value

A list with components:

Reproducibility

Both randomized backends draw random projections. They are run under a fixed internal seed and the caller's RNG state is restored afterwards, so repeated calls on the same data return the same decomposition and the function never perturbs the surrounding stream. Without that, two identical cluster4d_g3s runs could disagree on half their voxel labels.

See Also

cluster4d_g3s for the full G3S clustering algorithm. transform_new_data_svd for applying the compression to new data.

Examples

# Simulate fMRI time series data
n_voxels <- 1000
n_timepoints <- 300
latent <- matrix(rnorm(n_voxels * 8), n_voxels, 8)
loadings <- matrix(rnorm(n_timepoints * 8), n_timepoints, 8)
feature_mat <- latent %*% t(loadings) +
  matrix(rnorm(n_voxels * n_timepoints, sd = 0.05), n_voxels, n_timepoints)

# Compress to 15 dimensions (typical for G3S)
compressed <- compress_features_svd(feature_mat, n_components = 15)
print(compressed$variance_explained)  # At least 0.95
print(ncol(compressed$features))      # 15 here; may increase when necessary

# More aggressive compression
compressed_fast <- compress_features_svd(feature_mat, n_components = 10)

# Conservative compression (preserve 98% variance)
compressed_accurate <- compress_features_svd(
  feature_mat,
  n_components = 15,
  variance_threshold = 0.98
)

Description

The compute_centroids function calculates the center and centroid of each cluster given a feature matrix, grid, and cluster assignments.

Usage

compute_centroids(feature_mat, grid, assignment, medoid = FALSE)

Arguments

Value

A list containing two elements:

Examples

# Simple synthetic example
feature_mat <- matrix(rnorm(100), nrow = 10, ncol = 10)
grid <- cbind(x = runif(10), y = runif(10), z = runif(10))
assignment <- rep(1:2, each = 5)
centroids <- compute_centroids(feature_mat, grid, assignment)
print(names(centroids))

## Not run: 
  # Larger example with real neuroimaging data
  # Assuming `feature_mat`, `grid`, and `assignment` are available
  centroids <- compute_centroids(feature_mat, grid, assignment)
  # To compute medoids instead of means
  medoids <- compute_centroids(feature_mat, grid, assignment, medoid=TRUE)

## End(Not run)

Description

Calculates the functional gradient for each voxel as the average dissimilarity to its spatial neighbors in feature space. Low gradient values indicate stable, homogeneous regions (good seed locations); high values indicate boundaries.

Usage

compute_functional_gradient(feature_mat, coords, k_neighbors = 26)

Arguments

Details

For each voxel i, the functional gradient is:

$grad(i) = \frac{1}{k} \sum_{j \in N_k(i)} (1 - f_i \cdot f_j)$

where $N_k(i)$ are the k nearest spatial neighbors and $f_i \cdot f_j$ is the dot product (cosine similarity for unit-length vectors).

Value

Numeric vector of length N containing gradient values for each voxel.

Examples

## Not run: 
features <- matrix(rnorm(500 * 15), 500, 15)
features <- t(apply(features, 1, function(x) x / sqrt(sum(x^2))))
coords <- matrix(rnorm(500 * 3), 500, 3)

grad <- compute_functional_gradient(features, coords, k_neighbors = 26)
hist(grad, breaks = 50, main = "Functional Gradient Distribution")

## End(Not run)

Description

Removes trends from each voxel's timeseries by projecting out a basis. Supports polynomial or DCT (discrete cosine transform) basis functions.

Usage

detrend_basis_cpp(data, basis)

Arguments

Details

For orthonormal basis B, computes residuals as: y_detrend = y - B * (B^T * y)

This is equivalent to regressing out the basis and keeping residuals.

Value

Numeric matrix with basis components removed from each row

Examples

## Not run: 
# Remove linear + quadratic trend (polynomial degree 2)
basis <- make_poly_basis(n_time = 100, degree = 2)
data_detrend <- detrend_basis_cpp(data, basis)

# Remove low-frequency drift with DCT (first 5 components)
basis <- make_dct_basis(n_time = 100, n_basis = 5)
data_detrend <- detrend_basis_cpp(data, basis)

## End(Not run)

Description

Removes low-frequency components using discrete cosine transform.

Usage

detrend_dct_cpp(data, n_basis = 4L)

Arguments

Details

DCT detrending is commonly used in fMRI preprocessing. The number of basis functions controls the high-pass filter cutoff:

• n_basis = 1: mean removal only

• n_basis = 2: ~ linear detrend

• n_basis = 4-6: typical for fMRI (removes drift < ~0.01 Hz for TR=2s)

Higher n_basis = more aggressive high-pass filtering.

Value

Numeric matrix with low-frequency components removed

Description

Removes polynomial trend up to specified degree from each voxel.

Usage

detrend_poly_cpp(data, degree = 1L)

Arguments

Details

Convenience wrapper that generates polynomial basis and calls detrend_basis_cpp. For degree=1, this is equivalent to detrend_time_cpp but slightly slower.

Value

Numeric matrix with polynomial trend removed

Description

Fast removal of trends from each voxel's timeseries in 4D neuroimaging data. Supports linear, polynomial, and DCT (discrete cosine transform) detrending.

Usage

detrend_time(
  vec,
  mask = NULL,
  method = c("linear", "poly", "dct"),
  degree = 2,
  n_basis = 4
)

Arguments

Details

Three detrending methods are available:

Linear (method="linear"): Fastest option. Removes linear drift by fitting y = a + b*t.

Polynomial (method="poly"): Removes polynomial trend up to specified degree. degree=1 is linear, degree=2 adds quadratic, etc. Uses orthonormalized basis for stability.

DCT (method="dct"): Removes low-frequency components using discrete cosine transform basis. This is a standard high-pass filter for fMRI data. The n_basis parameter controls the cutoff:

• n_basis=1: mean removal only

• n_basis=2: approximately linear detrend

• n_basis=4-6: typical for fMRI (removes ~0.01 Hz drift for TR=2s)

All methods use fast C++ parallel implementation via RcppParallel.

Value

A NeuroVec with trend removed from each voxel

Examples

## Not run: 
vec <- neuroim2::read_vec("functional.nii.gz")

# Linear detrend (fastest)
vec_lin <- detrend_time(vec)

# Quadratic detrend
vec_quad <- detrend_time(vec, method = "poly", degree = 2)

# DCT high-pass filter (common for fMRI)
vec_dct <- detrend_time(vec, method = "dct", n_basis = 5)

## End(Not run)

Description

Removes linear trend from each voxel's timeseries using fast parallel linear regression.

Usage

detrend_time_cpp(data)

Arguments

Details

For each voxel, fits a linear model y = intercept + slope * t and subtracts the fitted values: residuals = y - fitted.

Uses RcppParallel for fast parallel computation across voxels. The linear regression is computed efficiently using precomputed time statistics.

Value

Numeric matrix with linear trend removed from each row (voxel)

Description

Backwards-compatible alias for rena_plus

Usage

er_ml_rena(bvec, mask, ...)

Arguments

Value

Same as rena_plus.

Description

Uses the elbow method or cumulative variance to suggest an appropriate number of components for SVD compression.

Usage

estimate_optimal_components(
  feature_mat,
  max_components = 50,
  method = c("variance", "elbow"),
  variance_target = 0.95
)

Arguments

Value

Integer; suggested number of components

Examples

## Not run: 
data <- matrix(rnorm(1000 * 300), 1000, 300)
n_opt <- estimate_optimal_components(data)
print(n_opt)

## End(Not run)

Description

Identifies optimal seed locations for G3S clustering by finding local minima in the functional gradient field. Unlike uniform grid seeding, this approach places seeds in the centers of stable functional regions rather than on boundaries.

Usage

find_gradient_seeds_g3s(
  feature_mat,
  coords,
  K,
  k_neighbors = 26,
  oversample_ratio = 3,
  min_separation_factor = 1,
  distance = c("cosine", "euclidean"),
  knn = NULL,
  spatial_scale = NULL
)

Arguments

Details

Algorithm

The seeding process follows these steps:

1. Functional Gradient Computation: For each voxel, compute the average dissimilarity to its k nearest spatial neighbors in feature space:

   $grad(i) = \frac{1}{k} \sum_{j \in N(i)} dist(f_i, f_j)$

   where $dist(f_i, f_j) = 1 - f_i \cdot f_j$ (cosine distance).

2. Candidate Selection: Rank voxels by gradient value and select the top K * oversample_ratio voxels with lowest gradient (most stable regions).

3. Spatial Separation: Starting with the lowest-gradient candidate, greedily select seeds that maintain minimum spatial separation from all previously selected seeds. If that pool cannot provide K seeds at the physical inhibition radius, deterministic farthest-point sampling over the full mask fills the shortfall.

Why Gradient-Based Seeding?

Traditional uniform grid seeding often places seeds on functional boundaries (high gradient regions), leading to:

• Poor initial cluster centroids

• Slower convergence

• Less biologically plausible parcellations

Gradient-based seeding ensures:

• Seeds are in functional cores (low gradient = stable, homogeneous regions)

• Better initial centroids = faster, more accurate clustering

• Biologically plausible: aligns with cortical organization

Comparison with Other Methods

• vs. Uniform Grid: G3S gradient seeding reduces boundary artifacts and improves convergence speed by 2-3x.

• vs. K-means Initialization: G3S respects spatial structure, preventing disconnected clusters.

• vs. Random Seeding: G3S is deterministic given the data, improving reproducibility.

Value

Integer vector of length K containing the row indices of selected seed voxels in the feature_mat/coords matrices.

See Also

compute_functional_gradient for gradient computation details. cluster4d_g3s for the full G3S clustering algorithm. spatial_gradient for spatial-only gradient (used in SNIC).

Examples

## Not run: 
# Simulate compressed fMRI features
n_voxels <- 1000
n_components <- 15
features <- matrix(rnorm(n_voxels * n_components), n_voxels, n_components)
features <- t(apply(features, 1, function(x) x / sqrt(sum(x^2))))  # Normalize

# Spatial coordinates (10x10x10 grid)
coords <- as.matrix(expand.grid(x = 1:10, y = 1:10, z = 1:10))

# Find 20 seeds
seeds <- find_gradient_seeds_g3s(features, coords, K = 20)
print(length(seeds))  # Should be 20 or close to it

# Visualize gradient values
grad_vals <- compute_functional_gradient(features, coords)
print(summary(grad_vals[seeds]))  # Should be low values

## End(Not run)

Description

Deprecated: This function is now available as find_gradient_seeds() in cluster4d_init.R. This version is maintained for backward compatibility.

This function finds the initial cluster centers for a supervoxel algorithm. Supervoxels are used to partition 3D image data into volumetric regions, grouping similar voxels together. The initial cluster centers are crucial for the performance and quality of the final supervoxels.

Usage

find_initial_points(cds, grad, K = 100)

Arguments

Value

A list containing two elements: selected - a vector of the selected indices corresponding to the initial cluster centers, coords - a matrix or data frame with the spatial coordinates of the initial cluster centers.

See Also

find_gradient_seeds, initialize_clusters

Description

Convenient synthetic data for examples, vignettes, and tests. Each scenario produces a NeuroVec (time series volume), a NeuroVol mask, and the ground-truth cluster labels so you can benchmark algorithms or build illustrations quickly.

Usage

generate_synthetic_volume(
  scenario = c("gaussian_blobs", "z_layers", "checkerboard"),
  dims = c(16, 16, 8),
  n_clusters = 4,
  n_time = 30,
  spread = c(3.5, 3.5, 2),
  noise_sd = 0.05,
  amplitude_range = c(0.8, 1.2),
  spacing_mm = c(3, 3, 3),
  pattern_type = c("orthogonal", "harmonic", "random"),
  seed = NULL
)

Arguments

Value

A list with elements vec, mask, truth, coords, patterns, weights, dims, n_clusters, and scenario.

Examples

syn <- generate_synthetic_volume(
  scenario = "gaussian_blobs",
  dims = c(12, 12, 6),
  n_clusters = 4,
  seed = 1
)
str(syn)
truth_vol <- array(syn$truth, syn$dims)
clusvol <- neuroim2::ClusteredNeuroVol(syn$mask > 0, clusters = syn$truth)
# Visualise a single slice (requires neuroim2 plot method):
## Not run: 
plot(clusvol)

## End(Not run)

Description

This function tests whether a cluster result object represents a partition. It is a method for the is.cl_partition generic from the clue package. For cluster_result objects, this always returns TRUE since cluster results represent valid partitions where each data point belongs to exactly one cluster.

Usage

## S3 method for class 'cluster_result'
is.cl_partition(x)

Arguments

Value

TRUE, indicating that cluster results are always valid partitions.

See Also

is.cl_partition for the generic function.

Description

Replace each voxel by the mean of its k nearest neighbors in its local spatial neighborhood.

Usage

knn_shrink(bvec, mask, k = 5, connectivity = 27)

Arguments

Value

A SparseNeuroVec or similar object with the smoothed data.

Examples

## Not run: 
mask <- neuroim2::NeuroVol(array(1, c(20,20,20)), neuroim2::NeuroSpace(c(20,20,20)))
bvec <- replicate(10,
                  neuroim2::NeuroVol(array(runif(20*20*20), c(20,20,20)),
                           neuroim2::NeuroSpace(c(20,20,20))),
                  simplify=FALSE)
bvec <- do.call(neuroim2::concat, bvec)
sbvec <- knn_shrink(bvec, mask, k=3)

## End(Not run)

Description

A scikit-learn style make_blobs for neuroimaging demos. Builds a small 2D/3D mask with three spatial bands (left/mid/right) and clearly separated time-series patterns plus mild Gaussian noise. Designed to be fast, deterministic, and visually intuitive in vignettes and tests.

Usage

make_block_synthetic(
  dims = c(12, 12, 1),
  ntime = 80,
  noise = 0.15,
  seed = NULL
)

Arguments

Value

List with components:

vec

NeuroVec 4D volume

mask

NeuroVol mask

truth

Integer ground-truth labels (length = prod(dims))

Description

Creates a discrete cosine transform (type II) basis matrix for detrending.

Usage

make_dct_basis(n_time, n_basis)

Arguments

Details

DCT basis functions are excellent for removing low-frequency drift in fMRI. The first basis is constant (mean), subsequent bases capture increasingly higher frequency components.

Value

Numeric matrix (n_time x n_basis) with orthonormal DCT basis

Description

Creates a polynomial basis matrix for detrending.

Usage

make_poly_basis(n_time, degree)

Arguments

Details

Uses Gram-Schmidt orthonormalization for numerical stability. Degree 0 = mean removal, 1 = linear detrend, 2 = quadratic, etc.

Value

Numeric matrix (n_time x (degree+1)) with orthonormalized polynomial basis

Description

Generates a volume with perfectly separable clusters for testing clustering algorithm correctness. Each cluster is a spatially contiguous cubic sub-region with identical timeseries within the cluster and orthogonal signals between clusters.

Usage

make_synthetic_clusters(
  n_time = 100,
  noise_sd = 0,
  n_clusters = 27,
  dim_per_cluster = 3,
  seed = NULL
)

Arguments

Details

For the default n_clusters=27, creates a 9x9x9 volume divided into 27 cubic 3x3x3 clusters arranged in a 3x3x3 grid:

z=1-3:           z=4-6:           z=7-9:
  y                y                y
  ^                ^                ^
  | 7 | 8 | 9 |    | 16| 17| 18|    | 25| 26| 27|
  | 4 | 5 | 6 |    | 13| 14| 15|    | 22| 23| 24|
  | 1 | 2 | 3 |    | 10| 11| 12|    | 19| 20| 21|
  +-----------> x  +-----------> x  +-----------> x

Each cluster has a unique orthogonal signal generated from QR decomposition of a random matrix, ensuring zero correlation between clusters.

This creates an "easy" test case where any reasonable clustering algorithm should achieve perfect or near-perfect recovery with noise_sd=0.

Value

A list with:

Examples

# Default: 27 clusters, perfect separability
data <- make_synthetic_clusters()
table(data$true_labels)  # 27 voxels per cluster

# 8 clusters (2x2x2 grid of 3x3x3 cubes = 6x6x6 volume)
data <- make_synthetic_clusters(n_clusters = 8)

# With noise for more realistic testing
data <- make_synthetic_clusters(noise_sd = 0.5)

Description

The merge_clus function combines a set of clustering results using a consensus clustering algorithm.

Usage

merge_clus(x, ...)

Arguments

Value

A ClusteredNeuroVol instance.

See Also

cl_consensus, as.cl_hard_partition, cl_ensemble

Examples

## Not run: 
# Assuming clustering1, clustering2, and clustering3 are objects of class "cluster_result"
merged_clustering <- merge_clus(clustering1, clustering2, clustering3, method="SE")

## End(Not run)

Description

This method of merge_clus is specifically designed to merge clustering results for ClusteredNeuroVol objects by performing consensus clustering across multiple clustering results.

Usage

## S3 method for class 'cluster_result'
merge_clus(x, method = "SE", ...)

Arguments

Value

A ClusteredNeuroVol instance representing the consensus partition.

See Also

cl_consensus, as.cl_hard_partition, cl_ensemble

Description

The meta_clust generic function is used to define methods for different classes of input objects.

Usage

meta_clust(x, cuts, ...)

Arguments

Value

Depends on the method called.

Description

The meta_clust function performs meta clustering on a given clustering result by applying hierarchical clustering on the cluster centers.

Usage

## S3 method for class 'cluster_result'
meta_clust(x, cuts = NULL, ...)

Arguments

Details

Parallelization Status

Currently NOT parallelized. Meta-clustering performs hierarchical clustering on cluster centers using sequential R functions.

Sequential Operations:

1. Distance Matrix Computation:

   • Computes correlation distance (1 - cor) between all cluster centers

   • O(K^2) pairwise correlations where K = number of input clusters

   • Sequential double loop in base R

2. Hierarchical Clustering: hclust()

   • Standard hierarchical agglomeration (Ward.D by default)

   • Sequential merging of closest clusters

   • O(K^2 log K) complexity

3. Tree Cutting: cutree()

   • Cuts dendrogram at multiple heights

   • Creates nested cluster assignments

   • Linear in number of clusters

Why Not Parallelized:

• Small scale: Usually operates on hundreds of clusters, not thousands of voxels

• R built-ins: Uses standard hclust() which is sequential

• Fast enough: Typically completes in seconds even for large K

• Memory efficient: Only stores K x K distance matrix

Potential for Parallelization:

• Distance matrix computation could use parallel correlation

• Some hierarchical clustering variants support parallelization

• Multiple cuts could be computed in parallel

Performance Characteristics:

• Input size: K clusters from initial clustering (typically 100-1000)

• Complexity: O(K^2) for distances, O(K^2 log K) for clustering

• Memory: O(K^2) for distance matrix

• Speed: Usually < 1 second for K < 500

Performance Tips:

• Pre-reduce clusters: Start with fewer initial clusters if meta-clustering is slow

• Use appropriate linkage: Ward.D is slower but often gives better results

• Limit cuts: Fewer cut levels = faster processing

• Alternative: Use consensus clustering (merge_clus()) for different approach

Use Cases:

Meta-clustering is useful for:

• Creating multi-resolution parcellations

• Hierarchical organization of functional regions

• Reducing large numbers of clusters to interpretable groups

• Finding stable cluster boundaries across scales

Value

A list containing:

See Also

hclust, cutree

Description

Reports build-time SIMD feature macros and (when available) a small set of runtime CPU feature flags, plus the currently-selected int8 dot-product kernel used by neurocluster's correlation-SLIC code paths.

Usage

neurocluster_simd_info()

Details

This is primarily a development/benchmarking aid. For definitive evidence that an ISA-specific kernel is in use, disassemble the package shared library and look for the corresponding instructions.

Value

A named list with fields compiler, arch, compile, runtime, dot_i8_kernel, shared_library, and notes.

Description

Combined fast preprocessing: removes volume mean offsets and linear trends in a single efficient operation.

Usage

normalize_detrend(vec, mask = NULL)

Arguments

Details

Combines normalize_volumes and detrend_time in a single pass for efficiency. The operations are applied in order:

1. Remove volume-wise mean offset

2. Remove linear trend from each voxel's timeseries

Uses fast C++ parallel implementation via RcppParallel.

Value

A NeuroVec with volume offsets and linear trends removed

Examples

## Not run: 
vec <- neuroim2::read_vec("functional.nii.gz")
vec_clean <- normalize_detrend(vec)

## End(Not run)

Description

Combined operation: first removes volume mean offsets, then removes linear trend from each voxel's timeseries. More efficient than calling both functions separately.

Usage

normalize_detrend_cpp(data)

Arguments

Details

Performs two operations in a single parallel pass:

1. Centers each volume to zero mean

2. Removes linear trend from each voxel's (now centered) timeseries

Value

Numeric matrix with volume offsets and linear trends removed

Description

Fast removal of volume-wise mean offset from 4D neuroimaging data. Each volume/timepoint is centered to have zero mean across voxels.

Usage

normalize_volumes(vec, mask = NULL)

Arguments

Details

This function removes global signal offset that can vary across timepoints due to scanner drift, physiological noise, or other factors. For each timepoint t:

$output[,,,t] = input[,,,t] - mean(input[,,,t])$

Uses fast C++ parallel implementation via RcppParallel.

Value

A NeuroVec with each volume centered to zero mean

Examples

## Not run: 
vec <- neuroim2::read_vec("functional.nii.gz")
vec_normalized <- normalize_volumes(vec)

## End(Not run)

Description

Removes the mean offset from each volume/timepoint in a matrix. This centers each volume to have zero mean across voxels.

Usage

normalize_volumes_cpp(data)

Arguments

Details

For each timepoint t, computes mean_t = mean(data[, t]) and subtracts it from all voxels: output[, t] = data[, t] - mean_t

Uses RcppParallel for fast parallel computation across voxels.

Value

Numeric matrix with each column (volume) centered to mean zero

Description

Creates visualization of clustering results. Shows axial, sagittal, and coronal slices through the clustered volume.

Usage

## S3 method for class 'cluster4d_result'
plot(x, slice = NULL, view = "all", colors = NULL, zlim = NULL, ...)

Arguments

Value

Invisibly returns the plotted data

Description

Methods for printing, summarizing, and plotting cluster4d results. Print cluster4d result

Usage

## S3 method for class 'cluster4d_result'
print(x, ...)

Arguments

Value

Invisibly returns x

Description

Print Method for G3S Results

Usage

## S3 method for class 'g3s_result'
print(x, ...)

Arguments

Description

Randomly permute cluster IDs to break any spatial ordering. This is useful when visualizing clusters with continuous color scales, as spatially-ordered cluster IDs create artificial gradients that can be misleading.

Usage

randomize_cluster_ids(cluster_result, seed = NULL)

Arguments

Details

Many clustering algorithms assign cluster IDs sequentially as clusters are encountered during processing. When voxels are processed in memory order (which is typically z-major in 3D arrays), this creates systematic ordering where cluster ID correlates with z-position (inferior to superior in brain imaging).

This spatial ordering creates two problems:

1. Visualization artifacts: When viewed with continuous color scales (e.g., viridis, rainbow), ordered cluster IDs create artificial gradients that suggest spatial structure where none exists. Randomized IDs show each cluster as a distinct color patch.

2. Interpretation bias: Ordered IDs may mislead users into thinking cluster numbering has anatomical significance (e.g., "cluster 1 is always inferior, cluster 100 is always superior").

This function breaks the ordering by randomly permuting cluster IDs while maintaining all cluster memberships and properties.

Value

A modified cluster_result object with randomized cluster IDs. All components (cluster vector, clusvol, centers, coord_centers, etc.) are updated consistently.

Examples

## Not run: 
# Standard workflow with ordered IDs
result <- cluster4d(vec, mask, K = 100, method = "flash3d")

# Randomize for better visualization
result_random <- randomize_cluster_ids(result, seed = 42)

# Compare visualizations
# Ordered IDs show gradient from inferior to superior
plot(result$clusvol)

# Randomized IDs show distinct color patches
plot(result_random$clusvol)

# Works with other supported clustering methods
result_snic <- cluster4d(vec, mask, n_clusters = 100, method = "snic")
result_snic_random <- randomize_cluster_ids(result_snic)

## End(Not run)

Description

Performs spatially-constrained hierarchical clustering using recursive 1-nearest neighbor aggregation. ReNA is a fast, linear-time agglomerative algorithm that avoids the "percolation" problem of standard hierarchical methods by using 1-NN graphs to ensure balanced cluster sizes.

Usage

rena(
  bvec,
  mask,
  K = 100,
  connectivity = 26,
  max_iterations = 50,
  verbose = FALSE,
  exact_k = TRUE
)

Arguments

Details

Algorithm Description

ReNA recursively aggregates features using a 1-Nearest Neighbor (1-NN) graph to ensure clusters remain roughly balanced in size and spatially compact. The algorithm proceeds in iterations:

1. Connectivity Constraint: Clustering respects the spatial topology graph G. Distance is only calculated between spatially adjacent voxels.

2. Similarity Calculation: For all connected edges, compute squared Euclidean distance in feature space.

3. 1-NN Subgraph: For every voxel, identify its single nearest neighbor. This forms a directed subgraph Q. The 1-NN graph is less likely to "percolate" (create giant components) compared to k-NN graphs where k >= 2.

4. Connected Components: Extract weakly connected components of Q. These components become the clusters for the current iteration.

5. Reduction: Features within each component are averaged to form a single "super-feature" for the next iteration. The graph G is contracted so edges between new clusters exist if any constituent voxels were connected.

6. Stopping Condition: Repeat recursively until the number of clusters drops to the desired K. If exact_k = TRUE, edges are pruned by distance to ensure exactly K clusters.

Performance Characteristics

• Complexity: O(N log N) per iteration where N = number of voxels

• Memory: O(N) for graph and assignments

• Speed: Generally faster than iterative methods (supervoxels)

• Iterations: Typically 5-15 iterations to reach K clusters from N voxels

Advantages over Other Methods

• No percolation: Unlike single-linkage, doesn't form one giant cluster

• Balanced sizes: 1-NN graph encourages roughly equal cluster sizes

• Topology-aware: Respects spatial structure via connectivity graph

• Deterministic: No random initialization

• Fast convergence: Linear-time operations per iteration

Parallelization Status

Currently NOT parallelized. ReNA uses a sequential recursive algorithm that processes the graph in a hierarchical manner. Each iteration depends on the previous iteration's clustering result.

Sequential Operations:

1. Distance Computation: Pairwise distances for connected voxels (C++)

2. 1-NN Finding: Each voxel finds its nearest neighbor (C++)

3. Component Detection: Union-Find algorithm for connected components (C++)

4. Feature Aggregation: Mean pooling within components (C++)

5. Graph Contraction: Building reduced graph for next iteration (C++)

Why Not Parallelized:

• Sequential dependency: Each iteration requires previous iteration's result

• Fast per-iteration: Linear-time operations make parallelization overhead high

• Small iteration count: Typically 5-15 iterations total

• C++ optimization: Core operations already optimized in C++

Performance Tips:

• Reduce connectivity: Use connectivity=6 instead of 26 for larger graphs

• Use smaller K: Fewer target clusters means fewer iterations

• Pre-smooth data: Reduces noise, improves cluster coherence

• Parallel alternatives: In the unified API, supervoxels and slic support parallel = TRUE. The other methods do not advertise a parallel clustering contract.

Value

A list of class rena_cluster_result (inheriting from cluster_result) with the following elements:

clusvol

An instance of type ClusteredNeuroVol.

cluster

A vector of cluster indices equal to the number of voxels in the mask.

centers

A matrix of cluster centers with each row representing the feature vector for a cluster.

coord_centers

A matrix of spatial coordinates with each row corresponding to a cluster.

n_clusters

Actual number of clusters found.

method

Clustering method used ("rena").

parameters

List of all parameters used.

metadata

List containing iteration count and convergence information.

Note

Consider using cluster4d with method = "rena" for a standardized interface across all clustering methods.

References

Hoyos-Idrobo, A., Varoquaux, G., Kahn, J., & Thirion, B. (2019). Recursive Nearest Agglomeration (ReNA): Fast clustering for approximation of structured signals. Pattern Recognition, 94, 17-28.

See Also

cluster4d for unified interface, supervoxels for iterative heat kernel clustering, snic for non-iterative priority queue clustering

Examples

## Not run: 
# Small synthetic example
library(neuroim2)
mask <- NeuroVol(array(1, c(10,10,10)), NeuroSpace(c(10,10,10)))
vec <- replicate(20,
                 NeuroVol(array(runif(10*10*10), c(10,10,10)),
                          NeuroSpace(c(10,10,10))),
                 simplify=FALSE)
vec <- do.call(concat, vec)

# Run ReNA clustering
rena_res <- rena(vec, mask, K=50, connectivity=26)
print(rena_res$n_clusters)

# With verbose output
rena_res <- rena(vec, mask, K=50, verbose=TRUE)

## End(Not run)

Description

Two-stage pipeline:

1. Edge-aware reciprocal-nearest-neighbor ReNA coarsening toward K' = ceiling(r * K) super-voxels.

2. Spatially-constrained Ward refinement on the super-graph to reach exactly K clusters.

Usage

rena_plus(
  bvec,
  mask,
  K = 100,
  r = 2,
  lambda = 1,
  grad_img = NULL,
  connectivity = 26,
  max_iterations = 50,
  verbose = FALSE,
  ...
)

Arguments

Details

A reciprocal-neighbor level is retained only when applying the complete level would leave at least K' super-voxels. This prevents the coarsening stage from overshooting below K and leaves the remaining one-at-a-time adjacent merges to Ward refinement. Ward costs are recomputed from the current size-weighted means before every accepted merge; equal-cost edges are resolved by ascending super-voxel identifiers. Because the adjacency constraint can expose a new edge after a merge, individual Ward deltas need not be ordered, although every delta is non-negative and the cumulative objective is non-decreasing. When Ward refinement is needed, its complete trace and queue diagnostics are available in result$metadata$ward.

Value

A cluster4d_result with additional class rena_plus_cluster_result.

Description

Compute and return the gradient volume that would be used for seed relocation in slic4d_supervoxels. Useful for debugging and visualization.

Usage

slic4d_grad_summary(
  bvec,
  mask,
  method = c("correlation", "intensity", "spatial")
)

Arguments

Value

A 3D array containing the gradient values.

Description

Cluster a 4D NeuroVec (x,y,z,time) into compact, spatially contiguous 3D supervoxels using an enhanced SLIC-style algorithm with mask-aware seeding, gradient-based seed relocation, and topology-aware K preservation.

Usage

slic4d_supervoxels(
  bvec,
  mask,
  K,
  compactness = 10,
  max_iter = 10,
  n_threads = 0,
  step_mm = NULL,
  n_components = 0,
  feature_norm = c("zscale", "l2", "none"),
  seed_method = c("mask_poisson", "mask_grid", "grid", "farthest"),
  seed_relocate = c("none", "correlation", "intensity", "spatial"),
  seed_relocate_radius = 1L,
  connectivity = c(26L, 6L),
  strict_connectivity = TRUE,
  enforce_connectivity = TRUE,
  preserve_k = FALSE,
  topup_iters = 2L,
  min_size = 0L,
  verbose = FALSE,
  .input_contract = NULL
)

Arguments

Details

Connectivity is the final partition constraint. With strict_connectivity = TRUE, every returned label has exactly one connected component under the requested neighborhood. If preserve_k = TRUE, an adjacency-preserving repair reaches exactly K whenever K is at least the number of connected mask components. When K is below that topological minimum, strict connectivity takes precedence, the minimum feasible number of labels is returned, and a warning is emitted. Returned feature and coordinate centers are recomputed from the final public labels in the original feature space and physical coordinate system.

Value

A cluster4d_result (also inheriting from cluster_result) with elements:

• clusvol: ClusteredNeuroVol with the final labels

• cluster: integer vector of length = #masked voxels

• centers: matrix (actual_k x d_feat) center features

• coord_centers: matrix (actual_k x 3) spatial centers in mm

Note

Consider using cluster4d with method = "slic" for a standardized interface across all clustering methods.

Examples

## Not run: 
  library(neuroim2)
  # Basic usage
  res <- slic4d_supervoxels(bvec, mask, K = 1000, compactness = 15)
  
  # With mask-aware seeding and gradient relocation
  res <- slic4d_supervoxels(bvec, mask, K = 1000, 
                           seed_method = "mask_poisson",
                           seed_relocate = "correlation",
                           preserve_k = TRUE)

## End(Not run)

Description

Compresses each voxel time series into a non-DC DCT sketch, uses cosine dissimilarity on a masked 3-D neighborhood graph, and applies Felzenszwalb-Huttenlocher segmentation. Multi-run fits use seeded DCT subspace and mode-weight resampling followed by uniformly weighted consensus. Nonzero sketches are normalized to unit length. A zero-information sketch remains zero and therefore has edge dissimilarity one to every neighbor.

Usage

slice_msf(
  vec,
  mask,
  target_k_global = -1,
  r = 12,
  compactness = 5,
  min_size = 80,
  num_runs = 3,
  consensus = TRUE,
  stitch_z = TRUE,
  nbhd = 8,
  gamma = 0,
  k_fuse = NULL,
  min_size_fuse = NULL,
  use_features = FALSE,
  lambda = 0.7,
  z_mult = 0,
  seed = 1L,
  ensemble_fraction = 0.8
)

Arguments

Value

A slice_msf_cluster_result and cluster4d_result. cluster is in mask order, centers is K by T, and coord_centers is K by 3. Multi-run results also contain native runs. metadata$exact_k_repair records the natural and requested K, pre/post sizes, whether and how repair ran, connectivity, minimum size, and effective gamma, z smoothing, seed, run count, and consensus controls for both natural and exact-K calls.

Experimental status

slice_msf() is available for explicit evaluation but is not a recommended production method. A previous reliability-weighting path collapsed distinct feature vectors to zero-distance edges, and the former exact-K repair could create singleton-dominated partitions. Those paths are now rejected or replaced and covered by regression tests, but release-level recertification is still pending. See vignette("slice-msf") for the evaluation contract.

References

Felzenszwalb, P. F., & Huttenlocher, D. P. (2004). Efficient graph-based image segmentation. IJCV, 59(2), 167-181.

Description

Builds neighbor-edge co-association scores from at least two run results and optionally mixes those scores with DCT-sketch similarity.

Usage

slice_msf_consensus(
  run_results,
  mask,
  nbhd = 8,
  k_fuse = 0.3,
  min_size_fuse = 80,
  use_features = FALSE,
  lambda = 0.7,
  stitch_z = TRUE,
  gamma = 0
)

Arguments

Value

A list with full-volume consensus labels and complete params.

Description

Low-level diagnostic interface returning native full-volume labels, a signed adjacent-pair temporal-smoothness diagnostic, and an r_used by voxel DCT sketch. Temporal smoothness does not affect segmentation.

Usage

slice_msf_single(
  vec,
  mask,
  r = 12,
  k = 0.32,
  min_size = 80,
  nbhd = 8,
  stitch_z = TRUE,
  gamma = 0,
  z_mult = 0,
  dct_frequencies = NULL,
  dct_weights = NULL
)

Arguments

Value

A list with full-volume labels, temporal_smoothness, sketch, and params.

Description

The SNIC function performs a spatially constrained clustering on a NeuroVec instance using the Simple Non-Iterative Clustering (SNIC) algorithm.

Usage

snic(vec, mask, compactness = 5, K = 500, max_iter = 100)

Arguments

Details

Performance Optimization (2025)

The SNIC implementation has been highly optimized using lightweight C++ structs and in-place operations, providing 10x-50x speedup over the original implementation. Key optimizations include:

• Elimination of Rcpp::List overhead in priority queue (uses lightweight struct)

• In-place centroid updates with no memory allocations

• Inline 26-connectivity neighbor iteration

• Direct pointer access for array operations

• Efficient scalar math (eliminates temporary vector allocations)

Parallelization Status

Currently NOT parallelized. SNIC uses a sequential priority queue-based algorithm that processes voxels in order of their distance from cluster centers.

Sequential Operations:

1. Initialization: Gradient-based seed selection using find_initial_points()

   • Finds K seed points with high gradient and spatial separation

   • Sequential search through candidate voxels

2. Priority Queue Processing (C++ implementation):

   • Maintains a global priority queue of voxels to be assigned

   • Each voxel assignment depends on previously processed neighbors

   • Voxels are processed in order of their combined distance metric

3. Distance Computation: For each voxel, calculates:

   • Feature distance to nearest cluster center

   • Spatial distance weighted by compactness parameter

   • Combined into single priority score

Why Not Parallelized:

• Sequential dependency: Priority queue enforces strict processing order

• Global state: Each voxel assignment affects subsequent assignments

• Algorithm design: SNIC's key innovation is its non-iterative, sequential nature

• Coherent clusters: Sequential processing ensures connected components

Performance Characteristics:

• Complexity: O(N log N) where N = number of voxels

• Memory: O(N) for priority queue and assignments

• Speed: Generally faster than iterative methods (supervoxels)

• Single pass: Processes each voxel exactly once

Performance Tips:

• Reduce K: Fewer clusters means less competition for voxels

• Adjust compactness: Higher values create more local clusters, faster processing

• Pre-smooth data: Reduce noise to improve gradient-based initialization

• Use smaller masks: Process ROIs separately if possible

• Parallel alternative: In the unified API, supervoxels and slic support parallel = TRUE. The other methods do not advertise a parallel clustering contract.

Comparison with Other Methods:

• Faster than: supervoxels() due to non-iterative nature

• Parallel choices: Use one of the two unified methods listed above when a tested sequential/parallel clustering contract is required.

• More coherent than: Methods without spatial priority (ensures connectivity)

Value

A list of class snic_cluster_result with the following elements:

clusvol

An instance of type ClusteredNeuroVol.

gradvol

A NeuroVol instance representing the spatial gradient of the reference volume.

cluster

A vector of contiguous positive cluster labels, one per included mask voxel.

centers

A K-by-T matrix of mean feature values in the original input space.

coord_centers

A matrix of spatial coordinates with each row corresponding to a cluster.

Note

Consider using cluster4d with method = "snic" for a standardized interface across all clustering methods.

See Also

supervoxels

Examples

## Not run: 
mask <- NeuroVol(array(1, c(20,20,20)), NeuroSpace(c(20,20,20)))
vec <- replicate(10, NeuroVol(array(runif(202020), c(20,20,20)),
NeuroSpace(c(20,20,20))), simplify=FALSE)
vec <- do.call(concat, vec)

snic_res <- snic(vec, mask, compactness=5, K=100)

## End(Not run)

Description

The spatial_gradient function calculates the spatial gradient of a NeuroVol instance within the specified mask.

Usage

spatial_gradient(vol, mask, sigma = 0.5)

Arguments

Value

A NeuroVol instance containing the spatial gradient values for the input vol.

See Also

spatial_laplacian, weighted_spatial_adjacency

Examples

## Not run: 
mask <- neuroim2::NeuroVol(array(1, c(20,20,20)), neuroim2::NeuroSpace(c(20,20,20)))
input_vol <- neuroim2::NeuroVol(array(runif(20*20*20), c(20,20,20)),
  neuroim2::NeuroSpace(c(20,20,20)))
gradient_vol <- spatial_gradient(input_vol, mask)

## End(Not run)

Description

Provides parameter recommendations based on data size and user priorities. Experimental methods are reported separately and are never returned as the recommended method.

Usage

suggest_cluster4d_params(
  n_voxels,
  n_timepoints,
  priority = c("balanced", "speed", "quality", "memory")
)

Arguments

Value

A list with suggested parameters for supported methods plus experimental_methods and an experimental_note. Experimental methods are not recommended or assigned suggested parameter sets.

Examples

# Get parameter suggestions for a typical fMRI dataset
params <- suggest_cluster4d_params(50000, 200, priority = "balanced")
print(params$recommended_method)
print(params$n_clusters)

# Speed-optimized parameters for large dataset
speed_params <- suggest_cluster4d_params(100000, 300, priority = "speed")
print(speed_params$recommended_method)

# Quality-optimized parameters for smaller dataset
quality_params <- suggest_cluster4d_params(10000, 150, priority = "quality")
print(quality_params$n_clusters)

# Memory-efficient parameters
memory_params <- suggest_cluster4d_params(200000, 100, priority = "memory")
print(memory_params$recommended_method)

Description

Summarize cluster4d result

Usage

## S3 method for class 'cluster4d_result'
summary(object, ...)

Arguments

Value

A summary list (invisibly)

Description

Cluster feature data on a cortical surface or mesh using a supervoxel-like approach.

Usage

supervoxel_cluster_surface(
  bsurf,
  K = 500,
  sigma1 = 1,
  sigma2 = 5,
  iterations = 50,
  connectivity = 6,
  use_medoid = FALSE
)

Arguments

Value

A list with:

clusvol

A NeuroSurface storing the final clustering result.

clusters

Integer vector of cluster assignments (one per vertex).

centers

Matrix of cluster centers.

coord_centers

Matrix of spatial centroid coordinates.

index_sets

List of vertex indices for each cluster.

Description

Cluster feature matrix (rows = time points) in a "supervoxel" style but over temporal dimension.

Usage

supervoxel_cluster_time(
  feature_mat,
  K = min(nrow(feature_mat), 100),
  sigma1 = 1,
  sigma2 = 3,
  iterations = 50,
  TR = 2,
  filter = list(lp = 0, hp = 0),
  use_medoid = FALSE,
  nreps = 5
)

Arguments

Value

A list of cluster results (one per repetition), each of which has the same structure as supervoxel_cluster_fit().

Examples

## Not run: 
feature_mat <- matrix(rnorm(100 * 10), 100, 10)
library(future)
plan(multicore)
cres <- supervoxel_cluster_time(t(feature_mat), K=5)

## End(Not run)

Description

Cluster a NeuroVec instance into a set of spatially constrained clusters.

Usage

supervoxels(
  bvec,
  mask,
  K = 500,
  sigma1 = NULL,
  sigma2 = 2.5,
  iterations = 50,
  connectivity = 27,
  use_medoid = FALSE,
  use_gradient = TRUE,
  alpha = 0.5,
  parallel = TRUE,
  grain_size = 100,
  num_threads = NULL,
  verbose = FALSE,
  converge_thresh = 0.001
)