Big-graph generation

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.

Logic

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.

Advantages

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.

Disadvantages

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])


Executing: ndmg_demo-dwi
Getting test data...
Archive:  /tmp/small_demo.zip
  inflating: small_demo/desikan.nii.gz  
  inflating: small_demo/KKI2009_113_1_DTI_s4.bval  
  inflating: small_demo/KKI2009_113_1_DTI_s4.bvec  
  inflating: small_demo/KKI2009_113_1_DTI_s4.nii  
  inflating: small_demo/KKI2009_113_1_MPRAGE_s4.nii  
  inflating: small_demo/MNI152_T1_1mm_brain_mask_s4.nii.gz  
  inflating: small_demo/MNI152_T1_1mm_s4.nii.gz  
Creating output directory: /tmp/small_demo/outputs
Creating output temp directory: /tmp/small_demo/outputs/tmp
This pipeline will produce the following derivatives...
DTI volume registered to atlas: /tmp/small_demo/outputs/reg_dwi/KKI2009_113_1_DTI_s4_aligned.nii.gz
Diffusion tensors in atlas space: /tmp/small_demo/outputs/tensors/KKI2009_113_1_DTI_s4_tensors.npz
Fiber streamlines in atlas space: /tmp/small_demo/outputs/fibers/KKI2009_113_1_DTI_s4_fibers.npz
Graphs of streamlines downsampled to given labels: /tmp/small_demo/outputs/graphs/desikan/KKI2009_113_1_DTI_s4_desikan.gpickle
Generating gradient table...
B-values shape (17,)
         min 0.000000 
         max 700.000000 
B-vectors shape (17, 3)
         min -0.991788 
         max 1.000000 
None
Aligning volumes...
Executing: eddy_correct /tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1.nii.gz /tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1_t2.nii.gz 16
Executing: bet /tmp/small_demo/KKI2009_113_1_MPRAGE_s4.nii /tmp/small_demo/outputs/tmp/KKI2009_113_1_MPRAGE_s4_ss.nii.gz -B
Executing: epi_reg --epi=/tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1_t2.nii.gz --t1=/tmp/small_demo/KKI2009_113_1_MPRAGE_s4.nii --t1brain=/tmp/small_demo/outputs/tmp/KKI2009_113_1_MPRAGE_s4_ss.nii.gz --out=/tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1_ta.nii.gz
Executing: flirt -in /tmp/small_demo/KKI2009_113_1_MPRAGE_s4.nii -ref /tmp/small_demo/MNI152_T1_1mm_s4.nii.gz -omat /tmp/small_demo/outputs/tmp/KKI2009_113_1_MPRAGE_s4_MNI152_T1_1mm_s4_xfm.mat -cost mutualinfo -bins 256 -dof 12 -searchrx -180 180 -searchry -180 180 -searchrz -180 180
Executing: flirt -in /tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1_ta.nii.gz -ref /tmp/small_demo/MNI152_T1_1mm_s4.nii.gz -out /tmp/small_demo/outputs/tmp/KKI2009_113_1_DTI_s4_t1_ta2.nii.gz -init /tmp/small_demo/outputs/tmp/KKI2009_113_1_MPRAGE_s4_MNI152_T1_1mm_s4_xfm.mat -interp trilinear -applyxfm
Beginning tractography...
Making Big Graph...
# of Streamlines: 8906
0
445
890
1335
1780
2225
2670
3115
3560
4005
4450
4895
5340
5785
6230
6675
7120
7565
8010
8455
8900
Generating graph for desikan parcellation...
{'ecount': 0, 'vcount': 70, 'region': 'brain', 'source': 'http://m2g.io', 'version': '0.0.50', 'date': 'Thu Aug  3 17:00:05 2017', 'sensor': 'Diffusion MRI', 'name': "Generated by NeuroData's MRI Graphs (ndmg)"}
# of Streamlines: 8906
0
445
890
1335
1780
2225
2670
3115
3560
4005
4450
4895
5340
5785
6230
6675
7120
7565
8010
8455
8900

 Graph Summary:
Name: Generated by NeuroData's MRI Graphs (ndmg)
Type: Graph
Number of nodes: 70
Number of edges: 773
Average degree:  22.0857
Execution took: 0:05:02.660297
Complete!

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)


# of Streamlines: 2
0
1

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))


8

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()))


8

We check that the voxel ids are the same:


In [8]:
print(poss_vertices == set(gra.nodes()))


True

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


True

In [ ]: