In this demo, we will verify that our big-graph generation code is functioning properly on a small portion of a real DWI dataset that we can manually verify very easily.
The logic of this function is essentially the same as that of the downsampling code written originally by Disa and Greg. The primary difference here is that, instead of looking for the ROI indices a particular voxel is part of, we instead this time define each voxel as its own index, and independently count streamlines for that particular voxel based on the Morton index of a given position.
This approach has the advantage that it is purely data-derived, and the graph we end up with will be totally invertible since the only way a voxel ends up as part of the graph is by a streamline existing. By definition of a streamline, a streamline must be between two or more voxels, so then each point will be connected to some other point.
This approach has the disadvantage that our ultimate graph will only have vertices if there exists a streamline linking the corresponding voxel of a given vertex to some other vertex. This means that small registration differences between subjects may lead to different vertex counts and different Morton indices corresponding to the same anatomical region due to that anatomical region being shifted by a voxel or two.
We begin by running Greg's small demo:
In [1]:
import ndmg
import ndmg.utils as mgu
# run small demo for experiments
print(mgu.execute_cmd('ndmg_demo-dwi', verb=True)[0])
The approach we will take is to take 2 fibers from our graph and verify that we end up with the appropriate voxels in our streamlines being connected:
In [2]:
import numpy as np
fibs = np.load('/tmp/small_demo/outputs/fibers/KKI2009_113_1_DTI_s4_fibers.npz')['arr_0']
small_fibs = fibs[1:3]
In [3]:
from ndmg.graph import biggraph as mgg
from ndmg.graph.zindex import XYZMorton
In [4]:
g1 = mgg()
g1.make_graph(small_fibs)
In [5]:
import networkx as nx
gra = nx.Graph()
gra.add_weighted_edges_from(g1.edge_list)
First, we should check to see that our graph ends up with the right number of vertices. We begin by looking at the floored values of the above voxel positions, since our image resolution is at 1mm scale:
In [6]:
poss_vertices = set() # use a set since we want unique elements
streamlines = []
for stream in small_fibs:
vertices = set()
for vertex in stream:
mid = str(XYZMorton(tuple(np.round(vertex)))) # morton index for vertex
vertices.add(mid)
poss_vertices.add(mid)
streamlines.append(vertices)
print(len(poss_vertices))
and we see that there are 8 unique possible vertices, defining a vertex as a unique point in 3-dimensional space at 1mm resolution. We then can check out the number of unique vertices in our corresponding graph:
In [7]:
print(len(gra.nodes()))
We check that the voxel ids are the same:
In [8]:
print(poss_vertices == set(gra.nodes()))
Indicating that our vertex indices appear to be correct. Let's check our streamlines to verify that the vertices each streamline is incident to are fully connected (and consequently have nonzero edge weight) in our resulting graph:
In [9]:
from itertools import combinations
edgect = 0 # count the number of edges we should have
for stream in streamlines:
combns = combinations(stream, 2) # stream is a list of vertices
for comb in combns:
edgect += 1
if gra.get_edge_data(*comb) == 0: # check the particular combination
raise ValueError('Edge should exist that isnt in the graph!')
Since we don't get any errors here, it is clear that every element that is in our graph should, in fact, be there. Using set notation, what we have shown is that:
\begin{align*} A \subseteq B \end{align*}where $A$ is the set of edges that we expect to have, and $B$ is the set of edges that actually exist in our resulting graph. However, we also want to show that:
\begin{align*} B \subseteq A \end{align*}so that we can conclude that $B = A$, or that our graph exactly matches the result we expect to end up with. To do this, we can simply check that the edges of $A$ are the only edges in $B$:
In [15]:
print(edgect == .5*nx.to_numpy_matrix(gra).sum()) # multiply by 2 the expected count since the graph is directed
# whereas the edgecount is undirected
In [ ]: