Viterbi Algorithm for Connecting Fragments in Long-Range Paths

The Viterbi Algorithm is a dynamic programming approach to finding an optimal path between a starting state and a goal state. For our neuron reconstruction problem, the states are defined as fragments, while the traversals are the connections made between endpoints using evidence observed from image intensity data. This notebook illustrates a simple example of the Viterbi Algorithm connecting synthetic fragments in a 10x10 grid example.

[1]:
import numpy as np
from brainlit.algorithms.connect_fragments.tests.grid_generator import grid_gen
from brainlit.algorithms.connect_fragments.dynamic_programming_viterbi import viterbi_algorithm
import matplotlib.pyplot as plt
/opt/buildhome/python3.7/lib/python3.7/site-packages/nilearn/datasets/__init__.py:96: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
  "Numpy arrays.", FutureWarning)

Running a 10x10 Grid Example

[2]:
# Load the 10x10 example from grid_generator
img, lbls, _ , somas = grid_gen(10)
plt.figure()
plt.title("Image data")
plt.imshow(img[:,:,0])
plt.figure()
plt.title("Pre-generated Labels")
plt.imshow(lbls[:,:,0])

# Initiate the algorithm class
alg = viterbi_algorithm(img, lbls, somas, [1,1,1])
print("Soma label: ", somas)
Soma label:  {5: [(9, 9, 0)]}
../../_images/notebooks_algorithms_viterbi_tutorial_3_1.png
../../_images/notebooks_algorithms_viterbi_tutorial_3_2.png
[3]:
# Manually identify the endpoints. Note that Labels 2 and 5 are "blobs" and do not have endpoints
endpoints = {}
endpoints[1] = ((0,0,0),(0,2,0))
endpoints[3] = ((3,3,0),(6,5,0))
endpoints[4] = ((7,7,0),(7,8,0))
# Assign the endpoints
alg.end_points = endpoints

Running Viterbi Algorithm

Within the Viterbi class object, we first need to compute the distance matrices. After that, we can run the Viterbi algorithm to find the best path.

[4]:
def print_path(alg, path):
    c = alg.connection_mat
    path_lbls = path[1]
    for i in range(len(path_lbls)-1):
        from_lbl = path_lbls[i]
        to_lbl = path_lbls[i+1]

        print(f"From {from_lbl} to {to_lbl}: {c[0][from_lbl][to_lbl]}, {c[1][from_lbl][to_lbl]}")
[5]:
alg.compute_all_dists()
top_path, sorted_paths = alg.viterbi_frag(1, K=4, somas=alg.somas)
print(top_path)
print_path(alg, top_path)
(1.1241229611526466, [1, 3, 4, 5, 5])
From 1 to 3: [0 2 0], [3 3 0]
From 3 to 4: [6 5 0], [7 7 0]
From 4 to 5: [7 8 0], [9 9 0]
From 5 to 5: [0 0 0], [0 0 0]

From the path found, we can observe that the algorithm traversed from the closest endpoints on each fragment to the goal state soma with label 5, starting from fragment label 1. We can also observe that the algorithm chose to ignore the blob fragment label 2, as it is not a soma state nor is it biologically sound to have a blob as an intermediate state.