Graph Analysis

def contacts(
    labels: NDArray[Any],
    connectivity: Literal[4, 6, 8, 18, 26] = 26,
    surface_area: bool = True,
    anisotropy: tuple[Union[int, float], Union[int, float], Union[int, float]] = (1, 1, 1),
) -> dict[tuple[int, int], Union[int, float]]:
    """
    Get the N-connected region adjacency graph and contact area between regions.

    Parameters:
    - labels: 3D numpy array of integer segmentation labels
    - connectivity: 6, 18, or 26 (default) for 3D; 4 or 8 for 2D
    - surface_area: Return contact surface area (True) or simple voxel count (False)
    - anisotropy: Weights for x, y, and z dimensions for computing surface area

    Returns:
    - dict: Mapping {(label_1, label_2): contact_value, ...}
            contact_value is surface area (float) or voxel count (int)
    """

import cc3d
import numpy as np

# Create test segmentation
labels = np.zeros((100, 100, 100), dtype=np.int32)
labels[20:40, 20:40, 20:40] = 1  # Component 1
labels[35:55, 35:55, 35:55] = 2  # Component 2 (overlapping)
labels[60:80, 60:80, 60:80] = 3  # Component 3 (separate)

# Get contact surface areas between all pairs
surface_contacts = cc3d.contacts(labels, connectivity=26, surface_area=True)

# Print all contacts
for (label1, label2), area in surface_contacts.items():
    print(f"Contact between regions {label1} and {label2}: {area:.2f} surface area")

# Get simple voxel contact counts instead of surface area
voxel_contacts = cc3d.contacts(labels, connectivity=6, surface_area=False)
for (label1, label2), count in voxel_contacts.items():
    print(f"Regions {label1} and {label2} share {count} face-adjacent voxels")

# Handle anisotropic voxel spacing (e.g., microscopy data)
# Voxel dimensions: 0.5μm x 0.5μm x 2μm (z is thicker)
anisotropic_contacts = cc3d.contacts(
    labels, 
    connectivity=26,
    surface_area=True,
    anisotropy=(0.5, 0.5, 2.0)
)

# Different connectivity patterns
face_contacts = cc3d.contacts(labels, connectivity=6)   # Face-adjacent only  
edge_contacts = cc3d.contacts(labels, connectivity=18)  # Faces + edges
corner_contacts = cc3d.contacts(labels, connectivity=26) # Faces + edges + corners

# Check if specific regions are adjacent
if (1, 2) in surface_contacts:
    contact_area = surface_contacts[(1, 2)]
    print(f"Regions 1 and 2 are adjacent with area {contact_area}")
else:
    print("Regions 1 and 2 are not adjacent")

# Find all neighbors of a specific region
region_of_interest = 2
neighbors = []
for (label1, label2) in surface_contacts.keys():
    if label1 == region_of_interest:
        neighbors.append(label2)
    elif label2 == region_of_interest:
        neighbors.append(label1)
print(f"Region {region_of_interest} neighbors: {neighbors}")

def region_graph(
    labels: NDArray[np.integer],
    connectivity: Literal[4, 6, 8, 18, 26] = 26,
) -> set[tuple[int, int]]:
    """
    Get the N-connected region adjacency graph of a 3D image.

    Parameters:
    - labels: 3D numpy array of integer segmentation labels
    - connectivity: 6, 18, or 26 (default) for 3D; 4 or 8 for 2D

    Returns:
    - set: Set of edges between labels as (label1, label2) tuples
    """

import cc3d
import numpy as np

# Create complex segmentation
labels = create_complex_segmentation()  # Your segmentation function

# Get adjacency graph
edges = cc3d.region_graph(labels, connectivity=26)

# Print graph structure
print(f"Found {len(edges)} adjacency relationships")
for label1, label2 in sorted(edges):
    print(f"Regions {label1} and {label2} are adjacent")

# Convert to other graph representations
# NetworkX graph
import networkx as nx
G = nx.Graph()
G.add_edges_from(edges)

# Adjacency list
from collections import defaultdict
adjacency_list = defaultdict(list)
for label1, label2 in edges:
    adjacency_list[label1].append(label2)
    adjacency_list[label2].append(label1)

# Find connected components in the adjacency graph
connected_groups = list(nx.connected_components(G))
print(f"Found {len(connected_groups)} connected groups of regions")

# Find regions with most neighbors
degree_centrality = nx.degree_centrality(G)
most_connected = max(degree_centrality.items(), key=lambda x: x[1])
print(f"Most connected region: {most_connected[0]} with {most_connected[1]:.2f} centrality")

# Compare different connectivity patterns
edges_6 = cc3d.region_graph(labels, connectivity=6)
edges_26 = cc3d.region_graph(labels, connectivity=26)
print(f"6-connected: {len(edges_6)} edges")
print(f"26-connected: {len(edges_26)} edges")

def voxel_connectivity_graph(
    data: NDArray[IntegerT],
    connectivity: Literal[4, 6, 8, 18, 26] = 26,
) -> NDArray[IntegerT]:
    """
    Extracts the voxel connectivity graph from a multi-label image.

    A voxel is considered connected if the adjacent voxel has the same label.
    The output is a bitfield representing allowed directions for transit between voxels.

    For 2D connectivity (4,8): returns uint8 with bits 1-8 encoding directions
    For 3D connectivity (6): returns uint8 with bits 1-6 encoding face directions  
    For 3D connectivity (18,26): returns uint32 with bits 1-26 encoding all directions

    Bit encoding for 3D (26-connected):
    Bits 1-6: faces (+x,-x,+y,-y,+z,-z)
    Bits 7-18: edges 
    Bits 19-26: corners
    Bits 27-32: unused (zero)

    Parameters:
    - data: Input labeled array
    - connectivity: Connectivity pattern to analyze

    Returns:
    - NDArray: Same shape as input, with connectivity bitfields
    """

import cc3d
import numpy as np

# Create test segmentation
labels = np.random.randint(0, 5, (50, 50, 50), dtype=np.int32)

# Generate voxel connectivity graph
vcg = cc3d.voxel_connectivity_graph(labels, connectivity=26)

# Analyze connectivity patterns
print(f"VCG dtype: {vcg.dtype}")  # uint32 for 26-connected
print(f"VCG shape: {vcg.shape}")  # Same as input

# Check connectivity at specific voxel
x, y, z = 25, 25, 25
connectivity_bits = vcg[x, y, z]
print(f"Voxel ({x},{y},{z}) connectivity: {connectivity_bits:032b}")

# Count total connections per voxel
total_connections = np.zeros_like(vcg, dtype=np.int32)
for bit in range(26):  # For 26-connected
    bit_mask = 1 << bit
    total_connections += ((vcg & bit_mask) > 0).astype(np.int32)

# Find highly connected voxels
highly_connected = np.where(total_connections > 20)
print(f"Found {len(highly_connected[0])} highly connected voxels")

# Different connectivity patterns
vcg_6 = cc3d.voxel_connectivity_graph(labels, connectivity=6)   # uint8
vcg_18 = cc3d.voxel_connectivity_graph(labels, connectivity=18) # uint32
vcg_26 = cc3d.voxel_connectivity_graph(labels, connectivity=26) # uint32

# Analyze 2D connectivity (using 2D slice)
labels_2d = labels[:, :, 25]
vcg_2d_4 = cc3d.voxel_connectivity_graph(labels_2d, connectivity=4)  # uint8
vcg_2d_8 = cc3d.voxel_connectivity_graph(labels_2d, connectivity=8)  # uint8

# Extract specific directional connections
# For 3D 26-connected, bit positions:
# 1: +x, 2: -x, 3: +y, 4: -y, 5: +z, 6: -z
plus_x_connections = (vcg & 1) > 0  # Can move in +x direction
minus_z_connections = (vcg & (1 << 5)) > 0  # Can move in -z direction

print(f"Voxels with +x connectivity: {np.sum(plus_x_connections)}")
print(f"Voxels with -z connectivity: {np.sum(minus_z_connections)}")

def color_connectivity_graph(
    vcg: NDArray[VcgT],
    connectivity: Literal[4, 6, 8, 18, 26] = 26,
    return_N: bool = False,
) -> Union[NDArray[VcgT], tuple[NDArray[VcgT], int]]:
    """
    Color the connectivity graph back into connected components.

    Given a voxel connectivity graph, this function treats it as an undirected
    graph and returns connected components based on the connectivity patterns
    encoded in the bitfield.

    Parameters:
    - vcg: Voxel connectivity graph from voxel_connectivity_graph()
    - connectivity: Must match the connectivity used to generate vcg
    - return_N: If True, also return the number of components

    Returns:
    - NDArray: Labeled image with connected components
    - tuple[NDArray, int]: If return_N=True, includes component count
    """

import cc3d
import numpy as np

# Start with original labels
original_labels = np.random.randint(0, 10, (100, 100, 100), dtype=np.int32)

# Extract voxel connectivity graph
vcg = cc3d.voxel_connectivity_graph(original_labels, connectivity=26)

# Reconstruct connected components from connectivity graph
reconstructed_labels = cc3d.color_connectivity_graph(vcg, connectivity=26)

# Get component count
reconstructed_labels, N = cc3d.color_connectivity_graph(
    vcg, connectivity=26, return_N=True
)
print(f"Reconstructed {N} components")

# Compare original vs reconstructed
print(f"Original max label: {np.max(original_labels)}")
print(f"Reconstructed max label: {np.max(reconstructed_labels)}")

# Note: Labels will be different but connectivity should be preserved
# Verify connectivity is preserved by checking region adjacency
original_graph = cc3d.region_graph(original_labels, connectivity=26)
reconstructed_graph = cc3d.region_graph(reconstructed_labels, connectivity=26)

print(f"Original graph edges: {len(original_graph)}")
print(f"Reconstructed graph edges: {len(reconstructed_graph)}")

# Modify connectivity graph before reconstruction
modified_vcg = vcg.copy()

# Remove some connections (set bits to 0)
# For example, remove all +x connections (bit 0)
modified_vcg = modified_vcg & ~1  # Clear bit 0

# Reconstruct with modified connectivity
modified_labels = cc3d.color_connectivity_graph(modified_vcg, connectivity=26)

# This should result in more components since we removed connections
original_N = len(np.unique(reconstructed_labels)) - 1  # Subtract background
modified_N = len(np.unique(modified_labels)) - 1
print(f"Original: {original_N} components")
print(f"After removing +x connections: {modified_N} components")

# Use with different connectivity patterns
vcg_6 = cc3d.voxel_connectivity_graph(original_labels, connectivity=6)
labels_6 = cc3d.color_connectivity_graph(vcg_6, connectivity=6)

# 2D example
labels_2d = original_labels[:, :, 50]
vcg_2d = cc3d.voxel_connectivity_graph(labels_2d, connectivity=8)
reconstructed_2d = cc3d.color_connectivity_graph(vcg_2d, connectivity=8)

# Complete region analysis pipeline
labels = your_segmentation_image()

# Get adjacency relationships
edges = cc3d.region_graph(labels)
contacts = cc3d.contacts(labels, surface_area=True)

# Analyze voxel-level connectivity
vcg = cc3d.voxel_connectivity_graph(labels)

# Create connectivity-based alternative labeling
alternative_labels = cc3d.color_connectivity_graph(vcg)

# Compare region structures
print(f"Original regions: {len(np.unique(labels))}")
print(f"Alternative regions: {len(np.unique(alternative_labels))}")

# Convert to NetworkX for advanced analysis
import networkx as nx

contacts = cc3d.contacts(labels, surface_area=True)
G = nx.Graph()
for (u, v), weight in contacts.items():
    G.add_edge(u, v, weight=weight)

# Network metrics
centrality = nx.betweenness_centrality(G, weight='weight')
clustering = nx.clustering(G, weight='weight')
communities = nx.community.greedy_modularity_communities(G)

print(f"Found {len(communities)} communities in the region graph")