0
# Core Skeletonization
1

2
The core skeletonization functionality implements the TEASAR-derived algorithm for converting densely labeled image volumes into skeleton representations. This module provides the primary `skeletonize` function along with specialized point connection capabilities.
3

4
## Capabilities
5

6
### Main Skeletonization Function
7

8
Skeletonizes all non-zero labels in a 2D or 3D image using the TEASAR algorithm with support for parallel processing, soma detection, and various quality improvements.
9

10
```python { .api }
11
def skeletonize(
12
    all_labels,
13
    teasar_params=DEFAULT_TEASAR_PARAMS,
14
    anisotropy=(1,1,1),
15
    object_ids=None,
16
    dust_threshold=1000,
17
    progress=True,
18
    fix_branching=True,
19
    in_place=False,
20
    fix_borders=True,
21
    parallel=1,
22
    parallel_chunk_size=100,
23
    extra_targets_before=[],
24
    extra_targets_after=[],
25
    fill_holes=False,
26
    fix_avocados=False,
27
    voxel_graph=None
28
):
29
    """
30
    Skeletonize all non-zero labels in a given 2D or 3D image.
31

32
    Parameters:
33
    - all_labels (numpy.ndarray): 2D or 3D array of integer labels
34
    - teasar_params (dict): Algorithm parameters controlling path tracing
35
        - scale (float): Rolling ball invalidation scale multiplier (default: 1.5)
36
        - const (float): Minimum invalidation radius in physical units (default: 300)
37
        - pdrf_scale (float): Penalty distance field scale factor (default: 100000)
38
        - pdrf_exponent (int): Penalty field exponent, powers of 2 are faster (default: 8)
39
        - soma_detection_threshold (float): DBF threshold for soma detection (default: 750)
40
        - soma_acceptance_threshold (float): DBF threshold for soma acceptance (default: 3500)
41
        - soma_invalidation_scale (float): Soma-specific invalidation scale (default: 2)
42
        - soma_invalidation_const (float): Soma-specific invalidation constant (default: 300)
43
        - max_paths (int, optional): Maximum paths to trace per object
44
    - anisotropy (tuple): Physical dimensions of voxels (x,y,z) (default: (1,1,1))
45
    - object_ids (list, optional): Process only specified label IDs
46
    - dust_threshold (int): Skip connected components smaller than this (default: 1000)
47
    - progress (bool): Show progress bar (default: True)
48
    - fix_branching (bool): Improve quality of branched shapes (default: True)
49
    - in_place (bool): Modify input array in place (default: False)
50
    - fix_borders (bool): Center skeleton where shape contacts border (default: True)
51
    - parallel (int): Number of processes (1=single, <=0=all CPUs) (default: 1)
52
    - parallel_chunk_size (int): Skeletons per progress update (default: 100)
53
    - extra_targets_before (list): Extra target points before main algorithm
54
    - extra_targets_after (list): Extra target points after main algorithm
55
    - fill_holes (bool): Fill holes in connected components (default: False)
56
    - fix_avocados (bool): Use heuristics to combine nuclei with cell bodies (default: False)
57
    - voxel_graph (numpy.ndarray, optional): Custom connectivity graph
58

59
    Returns:
60
    dict: Dictionary mapping label IDs to Skeleton objects
61

62
    Raises:
63
    DimensionError: If input dimensions are incompatible
64
    """
65
```
66

67
#### Usage Example
68

69
```python
70
import kimimaro
71
import numpy as np
72

73
# Load segmented volume
74
labels = np.load("neurons.npy")
75

76
# Basic skeletonization
77
skeletons = kimimaro.skeletonize(labels, anisotropy=(4, 4, 40))
78

79
# Advanced skeletonization with custom parameters
80
skeletons = kimimaro.skeletonize(
81
    labels,
82
    teasar_params={
83
        "scale": 2.0,        # More aggressive invalidation
84
        "const": 500,        # Larger minimum radius
85
        "pdrf_scale": 50000, # Lower penalty field influence
86
        "soma_detection_threshold": 1000,  # Higher soma threshold
87
        "max_paths": 200,    # Limit paths per object
88
    },
89
    anisotropy=(16, 16, 40),        # 16x16x40 nm voxels
90
    dust_threshold=2000,            # Skip small objects
91
    parallel=4,                     # Use 4 processes
92
    fix_avocados=True,              # Combine soma artifacts
93
    fill_holes=True,                # Fill holes in shapes
94
    object_ids=[1, 5, 10, 15],      # Process specific labels only
95
    extra_targets_after=[(100, 200, 50)],  # Add extra endpoint
96
)
97

98
# Access individual skeletons
99
neuron_1 = skeletons[1]
100
print(f"Skeleton has {len(neuron_1.vertices)} vertices")
101
print(f"Skeleton has {len(neuron_1.edges)} edges")
102
```
103

104
### Point Connection
105

106
Creates a direct skeleton path between two specified points on a binary image, useful for tracing specific connections or when exact endpoints are known.
107

108
```python { .api }
109
def connect_points(
110
    labels,
111
    start,
112
    end,
113
    anisotropy=(1,1,1),
114
    fill_holes=False,
115
    in_place=False,
116
    pdrf_scale=100000,
117
    pdrf_exponent=4
118
):
119
    """
120
    Draw a centerline between preselected points on a binary image.
121

122
    Parameters:
123
    - labels (numpy.ndarray): Binary image (2D or 3D)
124
    - start (tuple): Starting point coordinates (x, y, z) in voxel space
125
    - end (tuple): Ending point coordinates (x, y, z) in voxel space
126
    - anisotropy (tuple): Physical voxel dimensions (default: (1,1,1))
127
    - fill_holes (bool): Fill holes in shapes before tracing (default: False)
128
    - in_place (bool): Modify input array in place (default: False)
129
    - pdrf_scale (float): Penalty field scale factor (default: 100000)
130
    - pdrf_exponent (int): Penalty field exponent (default: 4)
131

132
    Returns:
133
    Skeleton: Single skeleton connecting the two points
134

135
    Raises:
136
    ValueError: If start and end points are in disconnected components
137
    """
138
```
139

140
#### Usage Example
141

142
```python
143
import kimimaro
144
import numpy as np
145

146
# Create binary image of single object
147
binary_neuron = (labels == 12345)
148

149
# Connect two specific points
150
skeleton = kimimaro.connect_points(
151
    binary_neuron,
152
    start=(50, 100, 25),    # Starting coordinates
153
    end=(200, 300, 75),     # Ending coordinates  
154
    anisotropy=(32, 32, 40) # Voxel dimensions in nm
155
)
156

157
# The skeleton traces the shortest path through the object
158
print(f"Path length: {skeleton.cable_length()} nm")
159
```
160

161
### Synapse Target Conversion
162

163
Converts synapse centroid locations to skeleton target points for enhanced skeletonization around synaptic sites.
164

165
```python { .api }
166
def synapses_to_targets(labels, synapses, progress=False):
167
    """
168
    Convert synapse centroids to skeleton targets.
169

170
    Given synapse centroids and desired SWC labels, finds the nearest 
171
    voxel to each centroid for the corresponding label.
172

173
    Parameters:
174
    - labels (numpy.ndarray): Labeled volume
175
    - synapses (dict): Mapping of {label_id: [((x,y,z), swc_label), ...]}
176
    - progress (bool): Show progress bar (default: False)
177

178
    Returns:
179
    dict: Mapping of {(x,y,z): swc_label, ...} for use as extra_targets
180
    """
181
```
182

183
#### Usage Example
184

185
```python
186
import kimimaro
187

188
# Define synapse locations for different labels
189
synapses = {
190
    1: [((120, 200, 30), 3), ((150, 220, 35), 3)],  # Two synapses on neuron 1
191
    5: [((80, 180, 25), 4)],                         # One synapse on neuron 5
192
}
193

194
# Convert to target format
195
targets = kimimaro.synapses_to_targets(labels, synapses, progress=True)
196

197
# Use targets in skeletonization
198
skeletons = kimimaro.skeletonize(
199
    labels,
200
    extra_targets_after=list(targets.keys()),
201
    anisotropy=(16, 16, 40)
202
)
203
```
204

205
## Algorithm Details
206

207
The TEASAR algorithm works through several key phases:
208

209
1. **Distance Transform**: Computes distance-to-boundary (DBF) for each connected component
210
2. **Root Finding**: Identifies starting points using soma detection or maximum DBF values
211
3. **Penalty Field**: Combines DBF and euclidean distance to guide path selection
212
4. **Path Tracing**: Uses Dijkstra's algorithm to find optimal paths through penalty field
213
5. **Invalidation**: Applies rolling-ball invalidation around traced paths
214
6. **Iteration**: Repeats until all significant regions are covered
215

216
### TEASAR Parameters Explained
217

218
- **scale & const**: Control the rolling-ball invalidation radius = `scale * DBF + const`
219
- **pdrf_scale & pdrf_exponent**: Shape the penalty field = `pdrf_scale * DBF^pdrf_exponent + euclidean_distance`
220
- **soma_detection_threshold**: DBF values above this trigger soma detection analysis
221
- **soma_acceptance_threshold**: DBF values above this are treated as confirmed somas
222
- **soma_invalidation_scale & soma_invalidation_const**: Special invalidation for soma regions
223

224
### Performance Considerations
225

226
- **Parallel Processing**: Use `parallel > 1` for large datasets with many objects
227
- **Dust Threshold**: Higher values skip small artifacts but may miss fine structures
228
- **Memory Usage**: Large volumes may require careful memory management
229
- **DBF Computation**: Most expensive phase, scales with object complexity
230
- **Path Tracing**: Scales with skeleton complexity and path count