Algorithms¶
Fragment Generation¶

class
brainlit.algorithms.generate_fragments.
state_generation
(image_path: Union[str, Path], new_layers_dir: Union[str, Path], ilastik_program_path: str, ilastik_project_path: str, fg_channel: int = 0, soma_coords: List[list] = [], resolution: List[float] = [0.3, 0.3, 1], parallel: int = 1, prob_path: Union[str, Path] = None, fragment_path: Union[str, Path] = None, tiered_path: Union[str, Path] = None, states_path: Union[str, Path] = None)[source]¶ This class encapsulates the processing that turns an image into a set of fragments with endpoints etc. needed to perform viterbrain tracing.
 Parameters
image_path (str or pathlib.Path)  Path to image zarr.
new_layers_dir (str or pathlib.Path)  Path to directory where new layers will be written.
ilastik_program_path (str)  Path to ilastik program.
ilastik_project_path (str)  Path to ilastik project for segmentation of image.
fg_channel (int)  Channel of image taken to be foreground.
soma_coords (List[list])  List of coordinates of soma centers. Defaults to [].
resolution (List[float)  Resolution of image in microns. Defaults to [0.3, 0.3, 1].
parallel (int)  Number of threads to use for parallel processing. Defaults to 1.
prob_path (str or pathlib.Path)  Path to alrerady computed probability image (ilastik output). Defaults to None.
fragment_path (str or pathlib.Path)  Path to alrerady computed fragment image. Defaults to None.
tiered_path (str or pathlib.Path)  Path to alrerady computed tiered image. Defaults to None.
states_path (str or pathlib.Path)  Path to alrerady computed states file. Defaults to None.
 Raises
ValueError  Image must be four dimensional (cxyz)
ValueError  Chunks must include all channels and be 4D.
ValueError  Already computed images must match image in spatial dimensions.

compute_edge_weights
(self, frag_path=None)[source]¶ Create viterbrain object and compute edge weights

compute_states
(self, alg: str = 'nb')[source]¶ Compute entire collection of states
 Parameters
alg (string, optional)  algorithm to use for endpoint estimation. "nb" for neighborhood method, "pc" for principal curves method. Defaults to "nb"
 Raises
ValueError  erroneously computed endpoints of soma state
Connect Fragments¶

class
brainlit.algorithms.connect_fragments.
ViterBrain
(G: nx.Graph, tiered_path: str, fragment_path: str, resolution: List[float], coef_curv: float, coef_dist: float, coef_int: float, parallel: int = 1)[source]¶ 
compute_all_costs_dist
(self)[source]¶ Splits up transition computation tasks then assembles them into networkx graph

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

frag_soma_dist
(self, point: List[float], orientation: List[float], soma_lbl: int, verbose: bool = 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]

shortest_path
(self, coord1: List[int], coord2: List[int])[source]¶ Compute coordinate path from one coordinate to another.
 Parameters
 Raises
ValueError  if state sequence contains a soma state that is not at the end
 Returns
list of voxel coordinates of path
 Return type

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.
CubicHermiteChain
(x: np.array, y: np.array, left_dydx: np.array, right_dydx: np.array, extrapolate=None)[source]¶ A third order spline class (continuous piecewise cubic representation), that is fit to a set of positions and onesided derivatives. This is not a standard spline class (e.g. bsplines), because the derivatives are not necessarily continuous at the knots. A subclass of PPoly, a piecewise polynomial class from scipy.
 Parameters
x (np.array)  Independent variable, shape n.
y (np.array)  Independent variable, shape n x d.
left_dydx (np.array)  Derivatives on left sides of cubic segments (i.e. right hand derivatives of knots), shape n1 x d.
right_dy_dx (np.array)  Derivatives on right sides of cubic segments (i.e. left hand derivatives of knots), shape n1 x d.
extrapolate (np.array, optional)  If bool, determines whether to extrapolate to outofbounds points based on first and last intervals, or to return NaNs. If ‘periodic’, periodic extrapolation is used.

axis
¶
 Raises
ValueError  Left derivatives (left_dydx) and right derivatives (right_dydx) must have same shape.

class
brainlit.algorithms.trace_analysis.
GeometricGraph
(df: pd.DataFrame = None, root=1, remove_duplicates=False)[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.
 Parameters
df (pd.DataFrame, optional)  Indices, coordinates, and parents of each node in the swc. Defaults to None.
root (int, optional)  Sample index (in the dataframe) that is the root node. Defaults to 1.
remove_duplicates (bool, optional)  Whether to automatically remove consecutive nodes with the same location. Defaults to False. Defaults to False.

root
¶ Sample index corresponding to the root node.

spline_type
¶ Whether splines are bsplines or cubic hermite splines.

spline_tree
¶ Tree of splines that model the trace branches.

removed_duplicates
¶ Whether connected nodes with the same loc attribute were merged.
 Raises
ValueError  If two nodes in a row in the SWC file have the same location and the user has not opted to remove duplicates.

fit_spline_tree_invariant
(self, spline_type: Union[BSpline, CubicHermiteSpline] = BSpline, k=3)[source]¶ Construct a spline tree based on the path lengths.
 Parameters
spline_type (Union[BSpline, CubicHermiteSpline], optional)  spline type that will be fit to the data. BSplines are typically used to fit position data only, and CubicHermiteSplines can only be used if derivative, and independent variable information is also known. Defaults to BSpline.
k (int, optional)  Order of spline for bsplines. Defaults to 3.
 Raises
KeyError  Make sure all nodes have loc attribute.
ValueError  loc attributes must be flat arrays
ValueError  loc attributes must not be empty
ValueError  loc attributes must contain 3 coordinates
ValueError  loc attributes must be unique for different nodes.
ValueError  graph must be connected.
ValueError  graph cannot contain cycles.
ValueError  graph must be connected.
 Returns
graph of splines that model the branches.
 Return type
nx.DiGraph

brainlit.algorithms.trace_analysis.
compute_parameterization
(positions: np.array)[source]¶ Computes list of parameter values to be used for a list of positions using piecewise linear arclength.
 Parameters
positions (np.array)  nxd array containing coordinates of the n points in d dimentional space.
 Returns
TotalDist  n array containing parameter values, first element is 0.
 Return type
np.array
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)