Algorithms¶
Fragment Generation¶

class
brainlit.algorithms.generate_fragments.
state_generation
(image_path, ilastik_program_path, ilastik_project_path, chunk_size=[375, 375, 125], soma_coords=[], resolution=[0.3, 0.3, 1], parallel=1, prob_path=None, fragment_path=None, tiered_path=None, states_path=None)[source]¶ 

compute_states
(self)[source]¶ Compute entire collection of states
 Raises
ValueError  erroneously computed endpoints of soma state


brainlit.algorithms.generate_fragments.
get_seed
(voxel)[source]¶ Get a seed point for the center of a brain volume.
 Parameters
voxel (tuple:)  The seed coordinates in x y z.
 Returns
A tuple containing the (x, y, z)coordinates of the seed.
 Return type

brainlit.algorithms.generate_fragments.
get_img_T1
(img)[source]¶ Converts a volume cutout to a SimpleITK image, as wel as a SimpleITK image with scaled intensity values to 0255.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to convert to a SimpleITK image.
 Returns
img_T1 (SimpleITK.SimpleITK.Image)  A SimpleITK image.
img_T1_255 (SimpleITK.SimpleITK.Image)  A SimpleITK image with intensity values between 0 and 255 inclusive.

brainlit.algorithms.generate_fragments.
thres_from_gmm
(img, random_seed=2)[source]¶ Computes a numerical threshold for segmentation based on a 2Component Gaussian mixture model.
The threshold is the minimum value included in the Gaussian mixture modelcomponent containning the highest intensity value.

brainlit.algorithms.generate_fragments.
fast_marching_seg
(img, seed, stopping_value=150, sigma=0.5)[source]¶ Computes a fastmarching segmentation.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
seed (tuple)  The seed containing a coordinate within a known segment.
stopping_value (float)  The algorithm stops when the value of the smallest trial point is greater than this stopping value.
sigma (float)  Sigma used in computing the feature image.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
level_set_seg
(img, seed, lower_threshold=None, upper_threshold=None, factor=2, max_rms_error=0.02, num_iter=1000, curvature_scaling=0.5, propagation_scaling=1)[source]¶ Computes a levelset segmentation.
When root mean squared change in the level set function for an iteration is below the threshold, or the maximum number of iteration have elapsed, the algorithm is said to have converged.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
seed (tuple)  The seed containing a coordinate within a known segment.
lower_threshold (float)  The lower threshold for segmentation. Set based on image statistics if None.
upper_threshold (float)  The upper threshold for segmentation. Set based on image statistics if None.
factor (float)  The scaling factor on the standard deviation used in computing thresholds from image statistics.
max_rms_error (float)  Root mean squared convergence criterion threshold.
num_iter (int)  Maximum number of iterations.
curvature_scaling (float)  Curvature scaling for the segmentation.
propagation_scaling (float)  Propagation scaling for the segmentation.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
connected_threshold
(img, seed, lower_threshold=None, upper_threshold=255)[source]¶ Compute a thresholdbased segmentation via connected region growing.
Labelled pixels are connected to a seed and lie within a range of values.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
seed (tuple)  The seed containing a coordinate within a known segment.
lower_threshold (float)  The lower threshold for the region growth. Set by a 2component Gaussian mixture model if None.
upper_threshold (float)  The upper threshold for the region growth.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
confidence_connected_threshold
(img, seed, num_iter=1, multiplier=1, initial_neighborhood_radius=1, replace_value=1)[source]¶ Compute a thresholdbased segmentation via confidenceconnected region growing.
The segmentation is based on pixels with intensities that are consistent with pixel statistics of a seed point. Pixels connected to the seed point with values within a confidence interval are grouped. The confidence interval is the mean plus of minus the "multiplier" times the standard deviation. After an initial segmentation is completed, the mean and standard deviation are calculated again at each iteration using pixels in the previous segmentation.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
seed (tuple)  The seed containing a coordinate within a known segment.
num_iter (int)  The number of iterations to run the algorithm.
multiplier (float)  Multiplier for the confidence interval.
initial_neighborhood_radius (int)  The initial neighborhood radius for computing statistics on the seed pixel.
replace_value (int)  The value to replace thresholded pixels.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
neighborhood_connected_threshold
(img, seed, lower_threshold=None, upper_threshold=255)[source]¶ Compute a thresholdbased segmentation via neighborhoodconnected region growing.
Labelled pixels are connected to a seed and lie within a neighborhood.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
seed (tuple)  The seed containing a coordinate within a known segment.
lower_threshold (float)  The lower threshold for the region growth. Set by a 2component Gaussian mixture model if None.
upper_threshold (float)  The upper threshold for the region growth.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
otsu
(img, seed)[source]¶ Compute a thresholdbased segmentation via Otsu's method.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type

brainlit.algorithms.generate_fragments.
gmm_seg
(img, seed, random_seed=3)[source]¶ Compute a thresholdbased segmentation via a 2component Gaussian mixture model.
 Parameters
img (cloudvolume.volumecutout.VolumeCutout)  The volume to segment.
random_seed (int)  The random seed for the Gaussian mixture model.
 Returns
labels  An array consisting of the pixelwise segmentation.
 Return type
Connect Fragments¶

class
brainlit.algorithms.connect_fragments.
ViterBrain
(G, tiered_path, fragment_path, resolution, coef_curv, coef_dist, coef_int, parallel=1)[source]¶ 
compute_all_costs_dist
(self, frag_frag_func, frag_soma_func)[source]¶ Splits up transition computation tasks then assembles them into networkx graph
 Parameters
frag_frag_func (function)  function that computes transition cost between fragments
frag_soma_func (function)  function that computes transition cost between fragments

compute_all_costs_int
(self)[source]¶ Splits up transition computation tasks then assembles them into networkx graph

frag_frag_dist
(self, pt1, orientation1, pt2, orientation2, verbose=False)[source]¶ Compute cost of transition between two fragment states
 Parameters
pt1 (list of floats)  first coordinate
orientation1 (list of floats)  orientation at first coordinate
pt2 (list of floats)  second coordinate
orientation2 (list of floats)  orientation at second coordinate
verbose (bool, optional)  Print transition cost information. Defaults to False.
 Raises
ValueError  if an orientation is not unit length
ValueError  if distance or curvature cost is nan
 Returns
cost of transition
 Return type
[float]

frag_soma_dist
(self, point, orientation, soma_lbl, verbose=False)[source]¶ Compute cost of transition from fragment state to soma state
 Parameters
 Raises
ValueError  if either distance or curvature cost is nan
ValueError  if the computed closest soma coordinate is not on the soma
 Returns
cost of transition [list of floats]: closest soma coordinate
 Return type
[float]

Trace Analysis¶

brainlit.algorithms.trace_analysis.
speed
(x: np.ndarray, t: np.ndarray, c: np.ndarray, k: np.integer, aux_outputs: bool = False)[source]¶ Compute the speed of a BSpline.
The speed is the norm of the first derivative of the BSpline.
 Parameters
x  A 1xL array of parameter values where to evaluate the curve. It contains the parameter values where the speed of the BSpline will be evaluated. It is required to be nonempty, onedimensional, and realvalued.
t  A 1xm array representing the knots of the Bspline. It is required to be a nonempty, nondecreasing, and onedimensional sequence of realvalued elements. For a BSpline of degree k, at least 2k + 1 knots are required.
c  A dxn array representing the coefficients/control points of the Bspline. Given n realvalued, ddimensional points ::math::x_k = (x_k(1),...,x_k(d)), c is the nonempty matrix which columns are ::math::x_1^T,...,x_N^T. For a BSpline of order k, n cannot be less than mk1.
k  A nonnegative integer representing the degree of the Bspline.
 Returns
A 1xL array containing the speed of the BSpline evaluated at x
 Return type
speed
References: .. [Rbbccb9c002d51] Kouba, Parametric Equations.

brainlit.algorithms.trace_analysis.
curvature
(x: np.ndarray, t: np.ndarray, c: np.ndarray, k: np.integer, aux_outputs: bool = False)[source]¶ Compute the curvature of a BSpline.
The curvature measures the failure of a curve, r(u), to be a straight line. It is defined as
\[k = \lVert \frac{dT}{ds} \rVert,\]where T is the unit tangent vector, and s is the arc length:
\[T = \frac{dr}{ds},\quad s = \int_0^t \lVert r'(u) \rVert du,\]where r(u) is the position vector as a function of time.
The curvature can also be computed as
\[k = \lVert r'(t) \times r''(t)\rVert / \lVert r'(t) \rVert^3.\] Parameters
x  A 1xL array of parameter values where to evaluate the curve. It contains the parameter values where the curvature of the BSpline will be evaluated. It is required to be nonempty, onedimensional, and realvalued.
t  A 1xm array representing the knots of the Bspline. It is required to be a nonempty, nondecreasing, and onedimensional sequence of realvalued elements. For a BSpline of degree k, at least 2k + 1 knots are required.
c  A dxn array representing the coefficients/control points of the Bspline. Given n realvalued, ddimensional points ::math::x_k = (x_k(1),...,x_k(d)), c is the nonempty matrix which columns are ::math::x_1^T,...,x_N^T. For a BSpline of order k, n cannot be less than mk1.
k  A nonnegative integer representing the degree of the Bspline.
 Returns
A 1xL array containing the curvature of the BSpline evaluated at x
 Return type
curvature
References: .. [Rce97e449f49f1] Máté Attila, The Frenet–Serret formulas.

brainlit.algorithms.trace_analysis.
torsion
(x: np.ndarray, t: np.ndarray, c: np.ndarray, k: np.integer, aux_outputs: bool = False)[source]¶ Compute the torsion of a BSpline.
The torsion measures the failure of a curve, r(u), to be planar. If the curvature k of a curve is not zero, then the torsion is defined as
\[\tau = n \cdot b',\]where n is the principal normal vector, and b' the derivative w.r.t. the arc length s of the binormal vector.
The torsion can also be computed as
\[\tau = \lvert r'(t), r''(t), r'''(t) \rvert / \lVert r'(t) \times r''(t) \rVert^2,\]where r(u) is the position vector as a function of time.
 Parameters
x  A 1xL array of parameter values where to evaluate the curve. It contains the parameter values where the torsion of the BSpline will be evaluated. It is required to be nonempty, onedimensional, and realvalued.
t  A 1xm array representing the knots of the Bspline. It is required to be a nonempty, nondecreasing, and onedimensional sequence of realvalued elements. For a BSpline of degree k, at least 2k + 1 knots are required.
c  A dxn array representing the coefficients/control points of the Bspline. Given n realvalued, ddimensional points ::math::x_k = (x_k(1),...,x_k(d)), c is the nonempty matrix which columns are ::math::x_1^T,...,x_N^T. For a BSpline of order k, n cannot be less than mk1.
k  A nonnegative integer representing the degree of the Bspline.
 Returns
A 1xL array containing the torsion of the BSpline evaluated at x
 Return type
torsion
References: .. [R8b689f5f8f911] Máté Attila, The Frenet–Serret formulas.

class
brainlit.algorithms.trace_analysis.
GeometricGraph
(df=None)[source]¶ The shape of the neurons are expressed and fitted with splines in this undirected graph class.
The geometry of the neurons are projected on undirected graphs, based on which the trees of neurons consisted for splines is constructed. It is required that each node has a loc attribute identifying that points location in space, and the location should be defined in 3dimensional cartesian coordinates. It extends nx.Graph and rejects duplicate node input.

fit_spline_tree_invariant
(self)[source]¶ Construct a spline tree based on the path lengths.
 Raises
ValueError  check if every node is unigue in location
ValueError  check if every node is assigned to at least one edge
ValueError  check if the graph contains undirected cycle(s)
ValueErorr  check if the graph has disconnected segment(s)
 Returns
nx.DiGraph a parent tree with the longest path in the directed graph
 Return type
spline_tree

Soma Detection¶

brainlit.algorithms.detect_somas.
find_somas
(volume: np.ndarray, res: list)[source]¶ Find bright neuron somas in an input volume.
This simple soma detector assumes that somas are brighter than the rest of the objects contained in the input volume.
To detect somas, these steps are performed:
Check input volume shape. This detector requires the x and y dimensions of the input volumes to be larger than 20 pixels.
Zoom volume. We found that this simple soma detector works best when then input volume has size 160 x 160 x 50. We use ndimage.zoom to scale the input volume size to the desired shape.
Binarize volume. We use Otsu thresholding to binarize the image.
Erode the binarized image. We erode the binarized image with a structuring element which size is directly proportional to the maximum zoom factor applied to the input volume.
Remove unreasonable connected components. After erosion, we compute the equivalent diameter d of each connected component, and only keep those ones such that 5mu m leq d < 21 mu m
Find relative centroids. Finally, we compute the centroids of the remaining connected components. The centroids are in voxel units, relative to the input volume.
 Parameters
volume (numpy.ndarray)  The 3D image array to run the detector on.
res (list)  A 1 x 3 list containing the resolution of each voxel in nm.
 Returns
label (bool)  A boolean value indicating whether the detector found any somas in the input volume.
rel_centroids (numpy.ndarray)  A N x 3 array containing the relative voxel positions of the detected somas.
out (numpy.ndarray)  A 160 x 160 x 50 array containing the detection mask.
Examples
>>> # download a volume >>> dir = "s3://openneurodata/brainlit/brain1" >>> dir_segments = "s3://openneurodata/brainlit/brain1_segments" >>> volume_keys = "4807349.0_3827990.0_2922565.75_4907349.0_3927990.0_3022565.75" >>> mip = 1 >>> ngl_sess = NeuroglancerSession( >>> mip=mip, url=dir, url_segments=dir_segments, use_https=False >>> ) >>> res = ngl_sess.cv_segments.scales[ngl_sess.mip]["resolution"] >>> volume_coords = np.array(os.path.basename(volume_keys).split("_")).astype(float) >>> volume_vox_min = np.round(np.divide(volume_coords[:3], res)).astype(int) >>> volume_vox_max = np.round(np.divide(volume_coords[3:], res)).astype(int) >>> bbox = Bbox(volume_vox_min, volume_vox_max) >>> img = ngl_sess.pull_bounds_img(bbox) >>> # apply soma detector >>> label, rel_centroids, out = find_somas(img, res)