From 9e22b14874dc76255ee06a3ebacb89013226880d Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Mon, 25 Mar 2019 11:06:22 -0400 Subject: [PATCH] add updated notebook --- .gitignore | 4 + notebooks/graph_test.ipynb | 682 ++++++++++++++++++++++++++++++++++++ testing/create_test_data.sh | 93 +---- 3 files changed, 703 insertions(+), 76 deletions(-) create mode 100644 notebooks/graph_test.ipynb diff --git a/.gitignore b/.gitignore index cb86687..171dd0f 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,10 @@ *.nii.gz *.fib *.fib.gz +*.mif +*.mif.gz +*.txt +*.tck # vim *.swp diff --git a/notebooks/graph_test.ipynb b/notebooks/graph_test.ipynb new file mode 100644 index 0000000..01d3a1d --- /dev/null +++ b/notebooks/graph_test.ipynb @@ -0,0 +1,682 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tour of MITTENS features\n", + "\n", + "Here we download a ground-truth FOD dataset (in DSI Studio format) that was created from the ISMRM 2015 tractometer streamline set. This is the best-possible scenario where FODs are directly from fibers, not from a reconstruction.\n", + "\n", + "We also download an atlas and a seed region to build connectivity matrices and run shortest-path tractography." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\t\t\t\n", + "\t\t" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens._mittens:\n", + "Using\n", + "------\n", + " Step Size:\t\t0.8660 Voxels \n", + " ODF Resolution:\todf8\n", + " Max Angle:\t\t35.00 Degrees\n", + " Angle Weights:\tflat\n", + " Angle weight power:\t1.0\n", + "INFO:mittens.external.dsi_studio:Loading DSI Studio ODF data\n", + "WARNING:mittens.external.dsi_studio:Unable to load real affine image \n", + "INFO:mittens._mittens:Loaded ODF data: (69248, 321)\n", + "INFO:mittens._mittens:Assuming ODF vertices are LPS+\n", + "INFO:mittens._mittens:Assuming ODF vertices are LPS+\n" + ] + } + ], + "source": [ + "# Only necessary on mac os 10.14+\n", + "import os\n", + "os.environ['KMP_DUPLICATE_LIB_OK']='TRUE'\n", + "\n", + "# Remove data from previous runs\n", + "import glob\n", + "for f in glob.glob(\"*.nii*\") + glob.glob(\"*.fib*\") + glob.glob(\"*.mat\"):\n", + " os.remove(f)\n", + "\n", + "# Download testing data\n", + "from urllib.request import urlretrieve\n", + "urlretrieve(\n", + " \"https://upenn.box.com/shared/static/4t9x64evipkxaieq6fzy5fzjymcb3ont.gz\", \n", + " \"TOD_2mm.fib.gz\")\n", + "urlretrieve(\n", + " \"https://upenn.box.com/shared/static/n3vv3ob7oeyq8nt2fuzilla07wca1h07.gz\", \n", + " \"endpoint_atlas_min_500.nii.gz\")\n", + "urlretrieve(\n", + " \"https://upenn.box.com/shared/static/o4134g0o18ndhyp69ofmz2zoavdyvztg.gz\", \n", + " \"cst_seed.nii.gz\")\n", + "\n", + "\n", + "from mittens import MITTENS\n", + "mitn = MITTENS(reconstruction=\"TOD_2mm.fib.gz\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculating transition probabilities\n", + "\n", + "We see in the output above that the default parameters of $\\theta_{max}=35$ degrees and $s=\\sqrt{3}/2$ voxels will be used for calculating transition probabilities. Now we actually calculate the probabilities and save them in NIfTI format." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens._mittens:ODF 0/69248\n", + "INFO:mittens._mittens:ODF 10000/69248\n", + "INFO:mittens._mittens:ODF 20000/69248\n", + "INFO:mittens._mittens:ODF 30000/69248\n", + "INFO:mittens._mittens:ODF 40000/69248\n", + "INFO:mittens._mittens:ODF 50000/69248\n", + "INFO:mittens._mittens:ODF 60000/69248\n", + "INFO:mittens._mittens:Calculating None-Ahead CoDI\n", + "INFO:mittens._mittens:Calculating Order1 KL Distance\n", + "INFO:mittens._mittens:Writing singleODF results\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_a_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ai_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_as_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_i_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_l_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_la_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_lai_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_las_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_li_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_lp_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_lpi_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_lps_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ls_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_p_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_pi_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ps_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_r_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ra_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_rai_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ras_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_ri_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_rp_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_rpi_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_rps_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_rs_prob.nii.gz\n", + "INFO:mittens._mittens:Writing TOD_2mm_singleODF_s_prob.nii.gz\n", + "INFO:mittens._mittens:Pre-computing neighbor angle weights\n", + "INFO:mittens._mittens:ODF 0/69248\n", + "INFO:mittens._mittens:ODF 10000/69248\n", + "INFO:mittens._mittens:ODF 20000/69248\n", + "INFO:mittens._mittens:ODF 30000/69248\n", + "INFO:mittens._mittens:ODF 40000/69248\n", + "INFO:mittens._mittens:ODF 50000/69248\n", + "INFO:mittens._mittens:ODF 60000/69248\n", + "/Users/mcieslak/projects/MITTENS/mittens/_mittens.py:406: RuntimeWarning: divide by zero encountered in true_divide\n", + " self.doubleODF_results = outputs / np.nansum(outputs, 1)[:,np.newaxis]\n", + "/Users/mcieslak/projects/MITTENS/mittens/_mittens.py:406: RuntimeWarning: invalid value encountered in true_divide\n", + " self.doubleODF_results = outputs / np.nansum(outputs, 1)[:,np.newaxis]\n", + "INFO:mittens._mittens:Calculating Double ODF CoDI\n", + "INFO:mittens._mittens:Calculating CoAsy\n" + ] + } + ], + "source": [ + "import os\n", + "#os.makedirs(\"101915/data/mittens_output\")\n", + "mitn.calculate_transition_probabilities(output_prefix=\"TOD_2mm\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This writes out NIfTI images for the transition probabilities to each neighbor using both the singleODF and doubleODF methods. Additionally, you will find volumes written out for CoDI and CoAsy. Writing to NIfTI is useful to quickly assess the quality of the output. the CoDI volumes should look a lot like GFA or FA images.\n", + "\n", + "One can load theses images directly to create a ``MITTENS`` object. This bypasses the need to re-calculate transition probabilities and is the recommended way to access transition probabilities before building and saving Voxel Graphs. \n", + "\n", + "### Loading from NIfTI outputs\n", + "\n", + "Here we demonstrate loading transition probabilites from NIfTI files. We verify that their contents are identical to those generated by calculating transition probabilities from a ``fib.gz`` file.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens._mittens:\n", + "Using\n", + "------\n", + " Step Size:\t\t0.8660 Voxels \n", + " ODF Resolution:\todf8\n", + " Max Angle:\t\t35.00 Degrees\n", + " Angle Weights:\tflat\n", + " Angle weight power:\t1.0\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_mask.nii.gz\n", + "INFO:mittens._mittens:Used to mask from 411516 to 69248 voxels\n", + "INFO:mittens._mittens:Loading singleODF results\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_a_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_a_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_as_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_as_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_i_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_i_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_l_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_l_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_la_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_la_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_lai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_lai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_las_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_las_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_li_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_li_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_lp_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_lp_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_lpi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_lpi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_lps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_lps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ls_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ls_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_p_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_p_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_pi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_pi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_r_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_r_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ra_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ra_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_rai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_rai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ras_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ras_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_ri_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_ri_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_rp_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_rp_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_rpi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_rpi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_rps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_rps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_rs_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_rs_prob.nii.gz\n", + "INFO:mittens._mittens:Loading NIfTI Image TOD_2mm_singleODF_s_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_s_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_singleODF_CoDI.nii.gz\n", + "INFO:mittens._mittens:Reading doubleODF results\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_a_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_a_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_as_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_as_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_i_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_i_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_l_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_l_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_la_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_la_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_lai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_lai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_las_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_las_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_li_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_li_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_lp_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_lp_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_lpi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_lpi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_lps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_lps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ls_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ls_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_p_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_p_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_pi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_pi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_r_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_r_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ra_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ra_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_rai_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_rai_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ras_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ras_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_ri_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_ri_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_rp_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_rp_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_rpi_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_rpi_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_rps_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_rps_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_rs_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_rs_prob.nii.gz\n", + "INFO:mittens._mittens:Loading TOD_2mm_doubleODF_s_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_s_prob.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_CoDI.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image TOD_2mm_doubleODF_CoAsy.nii.gz\n", + "INFO:mittens._mittens:Assuming ODF vertices are LPS+\n", + "INFO:mittens._mittens:Assuming ODF vertices are LPS+\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "singleODF outputs are identical: True\n", + "doubleODF outputs are identical: True\n" + ] + } + ], + "source": [ + "nifti_mitn = MITTENS(nifti_prefix=\"TOD_2mm\")\n", + "\n", + "import numpy as np\n", + "print(\"singleODF outputs are identical:\", \n", + " np.allclose(mitn.singleODF_results,nifti_mitn.singleODF_results))\n", + "print(\"doubleODF outputs are identical:\", \n", + " np.allclose(mitn.doubleODF_results,nifti_mitn.doubleODF_results,equal_nan=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You'll notice loading from NIfTI is much faster!\n", + "\n", + "## Building and saving a voxel graph\n", + "\n", + "Transition probabilities need to be converted to edge weights somehow. The ``MITTENS`` object accomplishes this through the ``build_graph`` function, which offers a number of options. Here is the relevant portion of the ``build_graph`` documentation\n", + "\n", + " Schemes for shortest paths:\n", + " ---------------------------\n", + "\n", + " ``\"negative_log_p\"``:\n", + " Transition probabilities are log transformed and made negative. This is similar to the Zalesky 2009 strategy.\n", + "\n", + " ``\"minus_iso_negative_log\"``:\n", + " Isotropic probabilities are subtracted from transition probabilities. Edges are not added when transition probabilities are less than the isotropic probability.\n", + "\n", + " ``\"minus_iso_scaled_negative_log\"``:\n", + " Same as ``\"minus_iso_negative_log\"`` except probabilities are re-scaled to sum to 1 *before* the log transform is applied. \n", + " \n", + " \n", + "You also have to pick whether to use singleODF or doubleODF probabilities. You can easily create graphs using either." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 69248/69248 [00:05<00:00, 11715.18it/s]\n", + "INFO:mittens.voxel_graph:Converting networkit Graph to csr matrix\n", + "INFO:mittens.voxel_graph:Saved matfile to TOD_2mm_doubleODF_nlp.mat\n" + ] + } + ], + "source": [ + "doubleODF_nlp = nifti_mitn.build_graph(doubleODF=True, weighting_scheme=\"negative_log_p\")\n", + "doubleODF_nlp.save(\"TOD_2mm_doubleODF_nlp.mat\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Notes\n", + "\n", + "#### Masking\n", + "If no spatial mask is specified when the ``MITTENS`` object is constructed, the generous brain mask estimated by DSI Studio will be used. This means the file sizes will be somewhat larger and the calculations will take longer. However, we recommend using the larger mask, as you can easily apply masks to the voxel graphs. It's better to have a voxel and not need it than to re-run transition probability calculation.\n", + "\n", + "#### Affines\n", + "``MITTENS`` mimics DSI Studio. Internally everything is done in LPS+ orientation and results are written out in RAS+ orientation. The affine used in the transition probability NIfTI images is the same as the images written by DSI Studio. This is a bizarre area of scanner coordinates and more closely matches the coordinates used by streamlines.\n", + "\n", + "### Verifying that the Voxel Graph is preserved\n", + "\n", + "Here we load a Voxel Graph directly from a matfile and check that it exactly matches the graph produced above." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens.voxel_graph:Loading voxel graph from matfile\n", + "INFO:mittens.voxel_graph:Loading graph from matfile\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "row indices match: True\n", + "column indices match: True\n", + "Weights all match: True\n" + ] + } + ], + "source": [ + "from mittens import VoxelGraph\n", + "mat_vox_graph = VoxelGraph(\"TOD_2mm_doubleODF_nlp.mat\")\n", + "\n", + "\n", + "from networkit.algebraic import adjacencyMatrix\n", + "matfile_adj = adjacencyMatrix(mat_vox_graph.graph,matrixType=\"sparse\")\n", + "mitn_adj = adjacencyMatrix(doubleODF_nlp.graph,matrixType=\"sparse\")\n", + "\n", + "\n", + "mat_indices = np.nonzero(matfile_adj)\n", + "indices = np.nonzero(mitn_adj)\n", + "\n", + "print(\"row indices match:\", np.all(mat_indices[0] == indices[0]))\n", + "print(\"column indices match:\",np.all(mat_indices[1] == indices[1]))\n", + "\n", + "print(\"Weights all match:\", np.allclose(matfile_adj[indices],mitn_adj[indices]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conclusions\n", + "\n", + "Here we've shown that transition probabilities and edge weights are preserved over the three ways they can be represented in ``MITTENS``. Next we show useful operations that can be performed using a VoxelGraph.\n", + "\n", + "## Connectomics \n", + "\n", + "Add an atlas to the graph:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens.spatial:Loading NIfTI Image ../testing/endpoint_atlas_min_500.nii.gz\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11420 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10470 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10828 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10880 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10803 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10755 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10859 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10486 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10529 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10748 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11100 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10534 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10715 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10706 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.12950 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10354 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10571 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11023 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11104 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10701 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10912 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10756 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10649 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11004 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10545 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11473 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10387 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11465 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10730 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10587 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11189 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11076 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11077 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10981 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10751 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11170 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.10189 sec\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.13730 sec\n", + "38it [00:00, 108.66it/s]\n" + ] + } + ], + "source": [ + "mat_vox_graph.add_atlas(\"../testing/endpoint_atlas_min_500.nii.gz\")\n", + "asym_raw_prob_graph, asym_mean_prob_graph, asym_path_length_graph, \\\n", + "conj_raw_prob_graph, conj_mean_prob_graph = mat_vox_graph.build_atlas_graph()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/mcieslak/miniconda3/envs/qsiprep/lib/python3.7/site-packages/IPython/core/magics/pylab.py:160: UserWarning: pylab import has clobbered these variables: ['f', 'indices']\n", + "`%matplotlib` prevents importing * from pylab and numpy\n", + " \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD0CAYAAAC7DZs3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXeYHMWZ/z89YfNKq4C0CosCSEW0OAQmGBvZcMeBAz7OZ99xcLYxP3wEnxBJBCMQBgzYQhYmGDDh7jA2xkbYZGyTDEeSCEKEEgKz0iqspEXSxtlJ/fujp2emu2pne1czyy5THz37PN2t6pqanq7uqrfe9/tatm1jMBjKi9An3QCDwTD0mI5vMJQhpuMbDGWI6fgGQxliOr7BUIaYjm8wlCGm4xsMIwQhxCFCiGc0x78qhHhVCPGiEOL/BanLdHyDYQQghLgA+CVQ5TseBZYC/wAcCZwmhGjsrz7T8Q2GkcEHwAma43sDa6WU26WUceB54PP9VRYpcuMMBkMG+fabdiodzDP2o/Ub316wYEEs79BtUsrbsnVJ+XshxHTNqaOAnXn7HcDo/j7PdHyDoUSkUmlmNdYFKtsTnxiTUh40iI9pB+rz9uuBHf2dZDq+wVBK7HSpP+FdYJYQYizQCXwB+Gl/J5mObzCUknRpOr4Q4kSgTkp5mxDiHOAJHJvdnVLKDf2db5noPIOhNLyz6jV7z/GVgcq+tSm2cu7cuYMZ6g8K88Y3GEqFTcne+LuK6fgGQ8mwh2KOPyiGpOMLIULAzcAcoBc4VUq5tsif8Tq5ZY2/SSm/W4Q6DwGulVLOE0LsCdyN8xxfDZwppdylX9VX/4HAQ8D7mf++RUp53yDqjAJ3AtOBSuBK4J1itb2P+luK0fZM/WHgdkAAKeC7gFWM9vdR9+hitV1LOlW0qorJUL3xvw5USSkPE0IcCiwBji9W5UKIKgAp5bwi1nkBcDLQlTl0PfBDKeUzQohf4LR/eRHrPxC4Xkq5ZPCtBuAkoE1KebIQYhzwOvAGxWu7rv4ritR2gK8CSCk/J4SYh3PdLYrTfl3dDxWx7V5sG1JBO3646B9fiKHy3DsCeBxASvkSUGwjxhygRgjxpBDiqczDZVfxe0rNBZ7NbD8GHF2C+r8shHhOCHGHEKK+j/P6437g0rz9JMVte1/1F6PtSCkfBE7L7E4DWilS+wvUXZS2+7Gxse10oL+hZqg6vt+7KCWEKOZooxtn7fIY4D+BX+1q/VLK3wOJvEOWlNJdAgnkHTXA+l8BzpdSfgH4ELhskPV2Sik7Mjfw74AfUsS291F/Udqe9xlJIcR/Az/PfEYx2++vu6ht9+Aa94L8DTFD1fH93kUhKWWyiPWvAe6RUtpSyjVAGzCpiPUD5P86gbyjBshyKeVKdxv4u8FWJIRoAp4G/ldKeS9Fbrum/qK13UVK+W1gNs6cvDrvv3a5/b66nyx22z3Y6WB/Q8xQdfwXgOMAMsPwt4pc/yk4dgOEEJNxRhibivwZr2fmhQDHAn8tcv1PCCE+m9k+ClhZqHBfCCEmAk8CC6WUd2YOF63tfdRflLZn6j9ZCHFRZrcb56G1ohjt76PuB4rVdhXbMe4F+Rtihsq4txz4eyHE/+EYanbZ4u7jDuBuIcTzOAOsU4o8ogA4F7hdCFGB4yb5uyLXfzpwoxAiDmwmNxcdKBcDY4BLhRDuXHw+cEOR2q6r/xzgZ0VoO8ADwF1CiOeAKHA2TpuLce11da+nONddzzBdzjOeewZDiXj7tZftPSs7A5VdHWswnnsGw6cDG1LFHngWB9PxDYZSYYNtl7cDj8FQhpS5y67BULYM0yCdIdfcE0IU12o6RHWXun7T9k+m/lK3vdzX8fMp5YUu7Y9o2v5J1F3q+ktXt23W8Q2G8mSYzvEHtY4/0DDbN15/za6IOIOLrW072G1cA2+sWuMpE7LUwYdlBWtPKDNwqRtdQ+fO7j5Kqd/Tf0R7LfIaMWp0De07uwmhNiztr01Tl+47pjPetKNG19K+s4t0wN8j7Ksrrlk2qopEs9u1o2roau8mpbkRLd/3SQW0RLvfx217UvPmiobUd0s6QGfIb9Oohlrad3T1cf2816smFFXKxNLqtbEz59WPrqEjc8/Yvrp2C1Ur543es3Hb3Llzd+v3CwBvv/qCvUdyXZCivF0xe0Ss4w8ozLYiEmJWY85Vf1ZjPYcdepGnTG1Flf805eYGtJ2uJuLIG12x7CIWzf+xtvP4f1RAKRdLxZUykVAuXPLaGxex8KwrqNTcXL3phGdf1wlqI+p3dD/zupsWccGZV9Cd7FXK6BhVUePZX9+xTSkjGqZkty9bdgGL51/HjkSXUq7C1zm39wZzOqmPOh3DbfuW7p1Kmcl1Y5Vjuu/o/y2iedd9yU2Xc+6Zl2uvn/+6z6ltUsq817NZORbPPAyuvfFSFp71IwCSvgfeafWfUc477oFzmpWDfWIXzbjX38tWCLEQ+DecuJjrpJQPF6pvsHP8UofZGgwjn+JG52VftsCFZGJTAIQQ+wMnAofiZNS5QghRo60lw2Df+Now27784+PJNO9v7gAglkjx/uYObrntx54yxRjqT9m9kSuWXdRHqV0f6k/dfRLX3rioJEP9qbtP5rqbFpVsqD+5qZHLll1QkqG+2/ZSDfWnTpvMkpsuL8lQf8ruk7j2xks9x1x0Q/2BYQd24NmxY8d4IcSKvEOehBr4XrZCiPyX7d7AM1LKGIAQ4n3gM8BLfX3eYDv+gMJs84f672/uMEN9zWeaof6IGeorxwoScKjf0NCwrZ+EGoVetm8BF2U0EiqAw4HbNHVkGWzHfwFHxui3QcJs31i1JtvRb7ntxxx26EVs++hJT5nOM7+nnHfWy2OUY4enapVjN8Zlvw2Oa5761eEKz35M83DOLxPCojpcobwhAWrCXhnl9rRqZPzPqr2UY5dsex5wHhRbe9qZNUqVEYj5bm5QO2tDlXpd8r+zbTv72genr+3xqOYNqXmYug9FC4tIKMyTDQcrZfY5ZKtybP/n1IdU1PeQv6di71z7QlXcU7E3J3S/oZx3cMMenv1XO/6mlGmIqtcmkbk2tm1nty3fEPPnO9QI3eOUI/1QPKt+ny9bKeW7QogbcdSJ1gIvA+pFzj95kI1YDsQyYbZLgQWDrMdg+PRiZ4J0gvz1T5+aFkKI3YDxUsojcEKwm3BESftkUG/8jMLpfw7mXIOhfCieVR+NpkUmg85aHMHQmUKIV4E4jpRYQeOCceAxGEqFTdGG+n28bN/L2/7+QOobko4fskJZ45277Z/T1910h3Le1IN+qBw78fAW5dgdz1Vk6raoDlVoLdf5xiKXiOU7ptrsqA/nLLthK0R9uFqZCwJUWt5L6Z+DA3wYVufqdXnXpa6iymNMzLZBU1eVz8DYHlZtCuOiuSlhxAoxLlpPd1g1rI2JeOfAurb3aAyfbhvCVojRkRp+Uam+ZJZ2q7aBOo2Rzm/IXV6d2z/Osni0OkR9QrWy7xX26m5+EN7SZzvz6Q07xyzLojKz7Tfa6to5YIZpkI554xsMJaOoQ/2iYjq+wVAqijjULzam4xsMJcO88Q2G8sOmvDX3LCtnvHG3/c45OkPelSuuVI51fP8U5djWeDvgeF5tjbcX9O4rhM6Ale+sk7LTdKZiioeXjoTGc29OWDUy3ZVwjG1p26Y70UtrLFiuiNm1kz37OnfWuJ3nwINN3E5qPff8xtDOZEwpo/N8dK9Xyk6zPdFFW1Q978VVU5RjLZ0vKsdGV3o9EddWdmS3e0mxNtWhdcJ6I9Hm2Q9pDK89afV37Uw4bU3b6ex2Muz9zXTuxgPDSG8ZDOWJGeobDGWGbeb4BkN5MkwT1piObzCUknJ/4+cb3EJYSpSdziNPZ8irv/VO5djY/U4CHE+8sdE6Ehrjm87o4/fc04Xb5keuhS2LmnCl1njo9zzbmVQ96T6IqO1yvcZCGQ+yCZVqBmgl1h+43Re4eHSPUoSKPG9CC4sKK6JE4gGMjdR59nWGvFhK9Tqsymt7baSSySFV++ELR6q5S3f78yjlmN9jcVo416YKQkwL17FK49U4Nzres/9YQg0NVjw0gYqwc20sy8puV/miNXX3zIAwQ32DoRwxKbQMhvLEzPENhjLD1dwbhgy64wshXicnBfQ3KWWxc94bDCOcT9kcXwhRBSClnBekfIhQ1nPO3fbLZbmhtfm4Hnn5uIa8fFasvgdw9PxWrL6HjtNUo+ARL6jWryReY9vHvR1KGa9sk0XYCtEYqVfKVfvCcldo2v5kfL16XiQTUoxFdaSCTbHtSpmOuNr2k0LTPPsJjVfbHtGcBbDSirBHdAxHWaon3fGHeQ2rX3xeNYbVRdRj54dmANBIJReF92RpXNWQ//Lz6i2mMxQmE16vv5ejG7PbR9kJXu7dqA2TvWe7V45rTv00pcx9B6sehaeudKTxa0OVHFw3HYAX2z/wlNFp/A2YT5nn3hygRgjxZKaOizMy2waDwcUGOz085/iDzaSzP46G9y+BWTgif6Ivpd1nnnrG3tr6MeBIYG9Yp6qe6pZOdD7xuqWZvfd13jyxRIqqaJjURx8pZd7r1MhK+z5S51+fv8TX2DSRzetbFWFIUNV/uzT+4TrVYFeNd+q0ybQ0b9QupaU0v1GNb+lJ58deH869scZOHcfHLW3U22obGuq8bV3TGeyemIjThsqmcfSub2MLhWMdXPzKuKBKneeLgUxqmsim9a1a9WT/PeK/LgAza9Tzmrud+2jc1PG0tTi6lJ0p78hAF/+wx+zpgTPevP3cE/aMdx4NUpR3D/6PgvUGSKhxHk5CjTRwtZRyeaHPG+wbfw2wVkppA2uEEG3AJEAdywJbWz9m0XxHR9+VwPZTHQo41I/WKcfyh/qzGuvpWLRUKfMtzVDf/7Bp7VEDZGbWNWa3F16/gGvPWRpsqN+lDntHR9V1bjcox5WQ1qkH6Yb6+zd4h7Tre1RR1S+Omp3d/rfrTuHXF9zJUUm1Df6h/snPqx1YhzvUn3b9iTSfcy9LbfU7V1jqLfZBl/rg98uRT6+bkN2+eOm5XL1gCd0pVT2ozTc9CzrUv2yl40twyk9O487zHSXqIEP9+564XTnWJzaQKlpCzD6zVwkhGoD/AvYEaoE3cDT6+mSwHf8UYH/gDCHEZBzNb9VTI4/8p7Xuya274XWOMjrnHHdOnzplAR2LllJ/m+rkw77/Vqh5APQk1Rs+X1IrhEWlFWGcpXGCsbzOP29qnE10uBrzadJ0J3vpTqg3d6Umqs/fCXSjieZkzpklbqdoTu5k/o43lXJdLx3m2e9MvhOo7T/ofQWAm9LH84OOV5heO0EpE0cdieh09f0jT93+xq6PlfN0UmVBOCbtOEqNIpzdfsr3+w+27hxFNe4VSqjRBTTjdPpaoN8PHWzHvwO4WwjxPM5z7ZRCCTUMhrJkAMt5ATLp9Je9aj3wDhAG1CG1j8HKa8dxcnUZDIZCBLShBcikUyh71bE4U+0Zmf0nhBAvSClf6auywSbUMBgM/WIXM2lmnwk1gO1AD9CbyZ+3A2goVJnx3DMYSoUNFG85r8+EGlLKPwohjgZeEkKkgeeBPxWqbIg6vp23TOVs+5efdLr3Orks3bKf65yz+Jtpx3qvMeS9+favlWNd80/17E9erv5Ia7tzFuhYOsHa7s2sC6kWdP8yY0dCtcRHqvo3FjXWqvkCdYY7fyTZhu2qRX1MRW4FJEWa9lQP9RWqNv0dSW++Od3n6ajKOB9ZhKiKVHADjUqZC0OqQ5JuCTnle+s9OTd3rTbVWDw5N0zTX1QzUk3Ue49sSqgrM4c9rxqEd8SdvHjLUsezqMPZjvmMe7p7csAMUUINKeVlwGVB6zNvfIOhVNg2drJoy3lFxXR8g6GUDFPPPdPxDYZS8inz1TcYDP1h2+X9xneMm7Znu9pnnNL54OvQlXOj7Cyrb7kkvyEPoHbZLz37dY98WSmT7z1oZfZ1xi//54ZDwQxkOV92CwtLcV0FqIxoJMF8Ls610f4TUYatkDbhp98nPajklHtt3Ouyz+Gq0fPBXvXGn/GMWr+/XXa+Hc929nVt939HXWyA7vvkH3O3/efqPm/AfJrCcg0GQ0DK+Y1vMJQltsmkYzCUHbZNeS/n2bZNLJNnzd2O+a+HOo3V5rLTSWC7yjmJdIrWnh3aKDudc45/Tr9h7SNKmal7fiX3PbBJpFPauV+FL89ae68qr61zXJlY7XhWRqwwE6sbaE+o5+2Iq/numju2ePZ1seNvfdyc3e5Jxnnr42bSmjfQ9linZ398tSp/3ZVQQ1t3xJx2JdIpWrt2MO0Jte0TqlW58NmjJivH/nKU12ax3yO577fk20mOfeVj6ir6t2PocgPq2u5GQabsdPb739NwhKfMhcn3lPMGjBnqGwzlRplb9Q2GssTGzPENhrLEvPENhvJjuIptBur4QohDgGullPOEEHsCd+MMZFYDZ2Yih/rGsnIyRpltvwNPfViNGtM5Y+hyv7kS2JWhKDPrGj1yWS75UXYufmmvfEOeS8vah7Pb72/uoGXtw3T9QHUG6mnx/sBHrVGNkP7vDNAWdwyTaWy6U71aw2RDZa1ybFylV/evK6kasO6Kiux2baSWxxsOZUKtaoBr3NurW3fMa+rPOaNaldU6rmEiAFOi9VzdOI9f9qjGsNERVeNvY0yV0JrzqNdINypPnzBkhRgVrdHqLfq1BvesUSMEl1jqvXVqog1wohz3aWgCYFF6rafMHprvPCBsG4apVb9f9zIhxAU4arquSfV64IdSys/jOG0dX7rmGQwjGDceP8jfEBPEr/QD4IS8/bnAs5ntx4Cji90og+FTwzDt+IF09YUQ04HfSCkPFUJslFJOzhz/Eo7QppreJo+nn3rG3tLqDK2m7j6JlnWblGG2zv9dp7wb1vpPO8dc3XudOm9Mo+XuL6XTtP/M/rNydWR0+9PrPlLK+WX0ZUyjGqzNHeCUc6+Lbp1dd23806C0Rlh1upVb9w43jSO1vo1oSC0XrfIeW9OtXgddG0ZnnC9qp46lq+VjtqXV6YbuPF3WH/+PEcp7J7m5GCKauuK2ty6dn0eT5v32UeY8V7MfVPVnv28GwNQ9pwTW1V/9+B/s3ZffHKQo7592deB6i8FgjHv5d0k9jr5XQba0trHwrCsAuPbGRSw864pAc3x/ggPQz/Hdm8vVvR/sHF+XUMM/x5/VWE/XVcuUcv45/r+uUefSheb47nXRaejr5vj+G7zfOf6yb9M1/78DzfG/o5njj9L8PseFnDn+QT85kRXn37tLc3z/AyJf037xsgu5bP41geb4M2smKmV0c/yrMnP8S352Hled/VNAfTk0VY5Vzrv6wauUY31SROmtQgk1hBAHAD/LK34o8HUp5eN91TeYjv+6EGKelPIZHHXPp/s7IYSVvVHd7SCRULpMOrq3uZvgIprJa6fTvdfJZSmRXZo25Bvy0ifPp+uqZdT+/JdKufT3vfn6HmlXPdb+aYvqVeY+DEJYVIcraEftmDr8o6G0ZrTyu+rc9/maBX+stnhohyrRdeeKfTz7leEtShkd9/Q6CShm2r3c0/uBdsSkS5SiyxPgvx/8ddnYTAqriUya2erZnx1WNSb3OFJ90FT9OXfdqzJt9Hs/jgkVIXde8YbxfSbUkFK+AcwDEEL8C7CxUKeHwXX8c4HbhRAVwLvA7wZRh8FQFhRxOa9QQg0AhBC1wGLgC/1VFqjjSyk/whk+IKVcAxwZvL0GQ5liA8lgHb8ICTUAvgfcL6VUh7c+jAOPwVAqbDvwG38XE2q4/DvwjSCfZxJqGAylpHjLeYUSaiCEGA1USim1iWv9DMkbP42dTY3sbvut8zpLvA7d8pCbqTaERbUVURJYgl6yy7+8plu+ybfWp+POvt+QB1B/qzdRp3Xm95QyyVbVWJkz0tna5UvQp8CuiHjbqvt+i0VuJWNzVZLFYjOPvKoaMB+s9t5427s7lTJjIqpFPT+pZCQUJp1Sb+CYrS6j6pKmhnzNyl+1CFmOQbhK8x39SUYPTKmG3Z1vK4e00lt+2bMm1NWAAVO8GJ2CCTWA2cBHQSszQ32DoUQ4AjzFMe4FSKjxKo7lPxCm4xsMpWR4RuWajm8wlIwivvGLzdB0fNvOzZ8y2+1pr6OKbn6t86TbmVQdXFbE2wE4IR1nRdc63tTkPNPlsvNLYOvksvKj7BbH0vzrmm6tc45/Tl930x1KmY9nf0055to6XNnxnoQanafLIdja43WY1EUyfu613HW4qjvFP7/WxcZO1Znl/tQqz36+15zLNrtdObaxy6krlkogd27QOubonHp0x/w+u1t6cytXiXSKLb072RxRHXj8dd2ealbK3LJBvY+6U45tIGmn2BJ3PsvvbNRsB3OmKoSt8U4eDpg3vsFQKmzMUN9gKEeGqfKW6fgGQ8kwb3yDoTwp6zd+yAplDUbu9n9W7eUp82FYdfSYozEWfRBRDTVPxh1npbAVYnRUDQMFiFT1n5tPp02QH0obspwIOl2Und85R2fIa17zR+VY7dR5AMRTSdZ1bOWMSZ9Tyvx251vKsYMb9vDsb4hvV8pMrciFlVaGosysmkCnJnx3atU4z74uHLpXo2fQVDcecAyzTXXj2b9K1cv/UXWvcuyQDWuVY37HrK+O3T+7PSpcxT+M2ou/dn2onLe8Zo5n/6x0i1ImqnH8cZ2l7Lztlu4273kag/NAGMYiu+aNbzCUjOGbQct0fIOhlNipImTcLQGm4xsMpcIGOz2CO75PXvtA4CHg/cx/3yKlvK9UDTQYRjIjdqifkdc+GXAtWgcC10splwT9kDTpbNJMd/uSbc97yugSIt6VUA1DOu+w6ohjgEumU7TGdtCdVM/T4fd2cxNY5uNq4oEjjNkW79Bq5/kj63TagK4hL5+ulmcAR8+vq+UZrba/Tpzyua3vevZ13nBbKnPeb12pXl5r/5tSBuBvXa3edkbV30KnBeh+50Q6yabu7XzUrkp2vVyjejk2f2O6cuzIx71eck/szH2/f0z18MTOdz1a+y7fjHlD7yrC6i2tM9rujDufl0yn2NLjXKc1B+7uKTPvPTVKcSDYgG2P3De+K6/9v5n9uYAQQhyP89Y/W0rZ0dfJBkPZMoyNe4OR1/4usEpKuVIIcQkwRkp5XqHzHXltRw1o6u6TaVm3UYl91qV51sdtq09QV4Bz6rTJtDRv1EpN6/HWpYtpzxex7Esa3MEvDqkST6lv7gPnOEq4rnT3qrfeV8ro3ua6a+Mn/5ruPm0K65o39HuO/7zc56nX1L13pk2bQnPzBm07I5q4iX1Hq6M22e6tP3+U0zRtCuubNxDWLst57yOdYKrux3BHK27bAfar8Y7kdBLpM2dPCyyD/eYDf7DH/viGIEXZ8ovrhr289nIppRshshz4eX8nbGndxgVnOvLa1920iAvOvIKtPd6gD91Q3y+yAIWH+ktuupxzz7y8qEN9N5gD+pYGB43qraZjruvYqhzLH+rPaqzni0csUsrohvr+YB5dpxtVmRsaL7vlR8w//VKlDKhD4YEO9W++9WrO+P7F9CbVtf4JmqH+qq+pqalO8g31W2O5IKSlNy9mwRmXaYf6O+Jen4qBDvVvve0avn/ahYA61P8XzVD/3sdvVY4VIsDz+RNhMB3/CSHED6SUrwBHASv7OyFt29nO6G7PGjXJ2xDNmyH/x3eZUKneSJtijvOKnVGx0T0wGmvHKMf8o472hBqNlZ/LLm2n6Yj3BJLA1kXZ6Zxz3Dn9sluu4ItHLPLo+LvsPuuryrEpNV6nG50Ueb6KTUUowuTqsWzsUaPz/Nde92CzKtQ3qfubhiyLynCUyTWqDr3OeWb2cnXk4Y+UHFORU/wJW2HGVNRpnYg+O2qmZ39NrFUp449kBK9Mu7t98GqvRqUuSnEg2Dakk8NT3W4wHf904EYhRBzYDJxW3CYZDJ8eRvQb3yev/RpweAnbZDB8OijiOn6hTDqZ/z8WuCyz+xpOFus+HzvDcxxiMHxKsG0r0F8Aspl0gAtxMukAIISoB34CfEVKeSiO6Ob4QpUZzz2DoVQMYDkvQEKNQpl0DseR214ihJgJ/FJKqVqS8/jEOr4/QWF9wEgoXY441+Kcsm064j1ay79Olrsy4i3ntxCDN2Fl2AppE1iCKoGtk8vSRdm5Fnsbm0Q6qTXkrXv/IeXY5D2O8+x3xtWIuoaqXFsT6SStsR3UR1XJaP+12RZTZbZ00t9jKh0DXMgKURetYluvet6EKtUY25tSjXThtLcNrcmcQS5pp2jt2aFdubg46r3OX9f8hrq2uytBlmVl74N8Q677f7uCjUU6oANPgIQahTLpjAe+CBwAdAJ/FUK8mMl6pcW88Q2GEpJOFW02XSiTThvwqpRyM4AQ4jmch0CfHd/M8Q2GUmFntPUD/AWgUCadlcB+QojxQogIjiH+nUKVmTe+wVAiHCGOovnqF8ykI4S4CHgiU/a3UsrVhSozHd9gKCFB5/j9ESCTzm+A3wStb0g6ftgKMaqixrPt19GvCqkGudm1qpTT7aoDHieFpgFQE65g/4ZpHjfbbP0ab7Qav456hxpdNq4yN62ycHK46YxF/lx2Om8xv1wW5KLsXE19v0ceqIY8gI0fPOrZf3a/i9TPOyYnJbW+oYL3vjaJ058dpZT7l16vwe/K2nVKmbBmVuh60oWtEHXhKuZX76uU+fbfq550M+5XvQf9rrb/MTonqTU+XMOpY+fy2853/adxcdgbH9ZUo65iVWjyMm6IOdfGlVMDqKzw3oN14V3z3IPAS3VDjnnjGwylIvj8fcgxHd9gKBFOdiTzxjcYyo70SJbe2lXiqSTrO7Z5tvOdSwDaw2rEmy4u/Gg1OjTrBBNPJ1nfs03rrLNhuzpv9Yef6j6vK0+OOk2armRM60Tkj+XX5bLTSWDnO6XY2NooO51zjn9Of+TqHytlZopc1uSfHpvg2Ie3kExvUsq96LN/6MKA/Q5K/v/bFNvOgZb64xz28E7lmC5y0S9lfW9Hzih9YLqHeztWs7Vbrcuvz6ALydbZd/zRouDkAMwn34loMAxEgWeoHw/mjW8wlArbdHyDoQwJ7rI71J50puMbDCXCRi/BNhwuqgeRAAAca0lEQVQwHd9gKCEj0qovhIgCdwLTgUrgShwf4LtxHmarcQL+CwYfVkWiiIYpnm2/sWhctF45L26rBiWdM8YeUcerpz5cxRdHzaY5qRqB8qWcXPxGwLc+blbK3BUV2e1aq4q7ooLfVas/5mKx2bP/uddUyan8XHYurgR2yAoxqrLGI5fl4jeEgtc5B7yGPJcP5YPZ7fc3d/ChfJBvzT1bKXf7HO/1siLq90t2qu+uE95xfouqUJTZdZO5S3PeXw9QowHnPKf+1v5IuH1rp2a3q0MV7Fs7ld6aRuU8P/ftrRr3TpeqluL6hGO4qwxFmVE7EVBzBi6pmqmcN1CGqwNPf1OLk4A2KeXngWOBG4HrgR9mjlnA8aVtosEwMrGBFFagv6Gmv45/P5AvzZrE0dV/NrP/GHB0CdplMHwqSNvB/oaaoLr69cAfgduBn0opJ2eOfwk4RUp5UqHzn3nqWXvbFmdoOrmpkY3rNyuujBHNGrpOeEG3Pl6ZGf6PnTqOj1vaiGvWwlMBtPb9QgwAsyJ5QhxN40itb2O75gE9pco7LXmvu7DqrUtXJq7A1b33xzCAfl19vwbv2vTbO1Vxi/332zO77er2f/jOeqXctBq1rX50SjIfZJbtJzZNpHV9K1WaadjUGrXt73RofAJ817Q6L45it6m7sbVlq9Z/ws/MKo2seUyddrn3SGPTRDavd+IJ/DEYTZYqphKd1RhY//7l3zxsd83/nyBFGf3owuGlqy+EaMIJCbxZSnmvEOK6vP+uB/r1cti2pY3F853TLlt2AYvnX1eSOf6/XXcKv77gTu0cvz2lOpcEmeM/3nBodrt22bfpmv/f/DHAHP+fX1OVYGZWqXryblorV/d+crVqB9DJjL/3Na88+bEPqwFG/jn+rMZ6fvjl25Ryg53jn5GZ45+75L9Ycu4NiIja9usOUBWgjn2uTTlWaI5/xk9P5+bzbqFXcz/40c3xFxeY41+49ByuWXA9oJnjW+ocf8Kj5/fbhhwW9icwjA9Cf8a9icCTwFlSyr9kDr8uhJgnpXwGZ97/dH8fkrLT7Eh0ebb92Wi6w+oP5p6Tjy4n3VGWYzist0Mclaxh/o43lTL1FaqRyX+z6bLFTKjNeRTGQ2km1Hbz0A7VC/CRV711bexUI9A6k6oHnh+d7r1OLssfZafzyMs35J3yk9P44Zdv476VP1PKbTvhe579I99TM6LpItXODjUBsBtRTrMn8au02qEP/z/1N9R51/k992R37vvE0glk9yaPF6WL/+H9+TfUpBshjUfhzkwOhXg6SUsmUs8/2joltUo5T816UJhhmkGr3zf+xcAY4FIhhDvXnw/cIISoAN4FflfC9hkMI5oR+caXUs7H6eh+jixNcwyGTw82I/eNbzAYdoFiLdUFSKhxA/A5wJ2nHS+lVI1dGUzHNxhKhA0UMSo3m1AjI7a5BK8PzYHAMVLKbdqzfQxJx7ewsoYTd9tvpBsTUb3TdBJXYyOqB97xh7UA0FIX5/jDWuh66TClzB3JvynH/GG422NqdtTGvXOGrk1VaRr37uDOFfso5R6s9lq979cYhqZWqbJaf+vKSVPZtq1NHqoLM/bLZflDa8Frrd9Yk+L2OTsVQx7A+Afu8OyLA/9LKRNHXfL7ykHOdd9Qm+Arh7TwmxXqb3hCzSzl2P8m1fwC/u84O89LrzIUZXZNIx0aSbW2hNcQOad6ilLm9nnq7/qHP88GYGqommuqDwDgmqQ3RfkxdbOV8wZKunhz/D4TamRGA7OA2zIG+TuklHcWqsy88Q2GEhLUNydAJp1CCTVqcdLVXw+EgaeFECuklOrbJ4Pp+AZDiRiIcS9AJp1CCTW6gWVSym4AIcRTOLaAPju+SahhMJSQtGUF+gtAoYQas4HnhRDhTGDdETgZc/tkSN74KTvF9t5Oz3Y86vXC0rmq6hxedG68X3zemRdf+g2bk5+P05lUk4jo5sl+2abx1ar09DGv5Z7ZF3fbfOe1NJVh1Utue7d3HlkbUR1e/J5hkJP/ClkhaqNVWannfHS57PwS2Dq3Xr8HnhWxtM45/jn9A6/doJRp/fL/U47940vOb3Hht+B7L9m0xFTPx/UVGo9MTVt7k16X4w3RnExZwk5qZcsAYimvm/UHCdWJ6Et/Us/rSDqS9FfZMa5MONuJtNeOsdnW6LwNkCK64feXUONXwEtAAvgfKeXbhSozQ32DoUTYQLJItr0ACTWuA64jIKbjGwwlwsYqplW/qJiObzCUECO9ZTCUIcNUVv+T6/h+HYCelBoLrzPk+bXPAeoiqtPLYOhKqMa3GdW5UNqwFWJUWI2UAxjjcyzaZqsGOTfXXD4dcceAlLbTdMR7sCrUO0XnyOTPZaczmOWH0tppZ18XZed3ztEZ8iY+crtyjM9819smjQFV1y6d7kG3JneAH10+Qr8T1pr2jUqZOQ3TlWO9Iee3sKycc1lLp9cwuDqYtb1PjK++wVCmmKG+wVBm2FbxrPrFxnR8g6GEjMihfh/y2i3AQ4Ab0XCLlPK+ErbRYBixDFN17X7f+K689slCiHHA68AVwPVSyiVBPyRkhbLyUe62PwqtSiNEqTP4VYXVcueHZgAwkQrOD83gB72vqOdFVI84v/zXjpgqE3Vcw8Ts9miiHBeayD29Hyjl/N9nY5cqodVUN1455hrubNsmZae1slRjKtWIRJ2h0I+rew9wbo+jkefKZeXjRtm5uB55HnyGPIAXVt0FOHp+L6y6i7/b70SlTIMm6lLniVhTM9qzn2/QtG1nf3SFWtfOuPc30xnydB6Ti0J7ADCJyuz2whrvtS+UKDQII9m4dz9eaS1XXlsIIY7HeeufLaVU/UANBsOw7fiDkdeuBFZJKVcKIS4Bxkgpzyt0/tNPPWNvaXX0AabuPpmWdRsVmWzdUpBuGcvvXw/QiBPbX9k0jt71baxPq29uSxOP5K/J76sNMCVP/bd26li6Wj7mY1t9K/vRLTsWks6eNm0Kzc0btN9Pl77bf710b6f8UZQrgb0b6ohpdK23re+rl0/LXvtOB3LS3e+8/aHaTs11j2lGK/5vHc5LOz5l90Y2rNusXd713yO6kaNOlrsRZwRY0TSO+HpnGW+D7R0Z6D5vz9kzAstgP3vfo/aq8+4NUpTDH1ww7OW1G6SU7oLqcpw44IJsad3GBWdeAcB1Ny3igjOvUIbGoyOqOup2jcpubURV2b0o7OjHT7v+RJrPuZcfdAxuqN/apa4TX904L7t90E9OZMX59wYa6sudG5QyuqH+pm4n+OTmW6/mjO9fTKVmKlMXVdfe/evxm2JqEMvsusnZbVcC+zR7klLuK4d4h/rf0w31NeQP9Wc11vPNo69VyuiG+u+2tyjH/NdvbN705oplF7Fo/o+1D2b/UF/Uq0IcuqH+xdZ0AKYs/Xc2LPgVAAvj3uCupCY/w4N/uls51hdFVuApKoOR135CCPEDKeUrwFHAyv4+JJlOsaV7p2f7yYaDPWV+Uale5Lao+oNNDqkPiKVxJ1LtYuIstdcxvVbVr78BNe/aPod7VYqmPdGtlPllTzYOgunpGL/seU/7JkinvMd0HXj/qsnKsY/anUg/G5veZILJNao2/bZedU48v3pfz/6BGgnp/Fx2VVYEERmrlcD2K+fooux0IzJ3Tn/5soV88+hreX21+nbrXuiPK4Fpv1VHcvkJNABeOCL3O2+oD/HCETXM+pMqIX7g6Bme/S0J9Vr5ZdQBzok5oepL0idwbkwftq77zgNl16wEpWMw8trnAD8TQsSBzcBpJWyfwTBiGbFpsgvIax9emuYYDJ8uRuRQ32Aw7BrD1apvOr7BUCJG7FC/WERDESbXjfVs73OIN5ni0m71Er24SrXQfuFI1cDz5efzpLutCHGNSeXCkGr1frDX+5kTqkcrZfJXG8JWiNGRGsUQBRCzvUtUOgPgj6rVZcCXM44rkVCYCTWjiVpqpOGEKrVd3/77Vs/+YQ+ruRP+ekAuknBjTZLrDtiqzWXnl8AOKpflWuzDhGiI1GoNeTXX/kI9tvxryjG/Aa7ltVwbkt8K0fJaPROrVQPmEWHvSslDKdVAq0u0mv08K7ftN+YdWjfDf9oAsQNl+A1Cfwk18so8AvxBSqle+DyM2KbBUELSAf8CkE2oAVyIk1DDz5WAuiykwXR8g6FE2EAq4F8APAk1AI+zjxDiGzjPkMeCVGbm+AZDCQlq1d+VhBpCiP2AE4FvAIuCfJ7p+AZDiXCCdILN8XcxocZ/AFOAp3AiaeNCiI+klI/3VdmQdPx0XtSZu73/c16vuTqNDn1L54vKsd3+rGrfu37xvekEH3Rt1ka46WISZjzjfRzPHqV61m2M5aLsEukkG2Mfa73y0r76dca9QzasVY41f2M6AOtHR1n1tQnMXq66+vZq/P5n3O+N/utJqJGMc57L3SdLvpXk2OfatNfGn8suqFyWG2UXSyd4t71F65GnM+Q1r/mjcuyg/U7y7B+zI+f3vzTVy4IdHzIqqnpt3rzDmzeiQRPB163Judfe6xgBU+l0dlseNNVT5itrVC/HgVJEq/4LwFeB3/oTakgpL3C3hRCXA5sLdXowb3yDoWQUOSy3YEKNgVZmOr7BUEKKtZzXX0KNvHKXB6nPdHyDoYSUtQOPwVCO2EBymHb9Iev4+cavtG0T9XlJ6UIgR1eqxhx/3DZAMqOHb+OE/eoMeam0Otvye4v95SjVI2/Oo3ntspx2+kVEAELKIbWM7jse+bhjWLr0H9Oc9Hg34ZBGdCOtHqsIe3+6qEbkw/P9Mh5qunL+dvkTWIJe996Vy7JwfhedR6MuJNZvyANYsfoez/5M8XVPHZFQmK0x1TvRL8fVkVS9+7o0bXeNlWnbzm0nvNehLa5+3kAZnt3evPENhpIxkjX3DAbDLqBb1h0OmI5vMJSQEfvGF0KEcUQ2BY5b8XdxpnV344xmVgNnZpYbtFhYRDNzc3f7noq9PWWWV6vz2LWVqnjvtLAqNf1y1MmXVhGKML1ugnaO/+Rc1TZg+/xU9ntki1Im32kkRIjaSJX2KV7pE3nc0qvOD786dn/l2BM73wUc56DW2A7GVKjfrzWpagH+x+g5nv17O1YrZfatzTmkVIcq2Ld2KrJbjW6cXeOVJdsQVSMZdbhCl2ErzNjKOo9clkt+lJ1LvnOOS/6cHuBD+WB2+/3NHXwoH2TKnl/ut+1bExrBZ1Wmkba4Uy4aDjOpdgwAB7ztvTbVYdVmMRDsIkbnFZsgQTpfBZBSfg7HD/j6zN8PpZSfx3kIHF+yFhoMIxg74N9Q02/Hl1I+SE5XbxrQiqOt/2zm2GPA0SVpncEwgnGX84L8DTWBdPUBhBD/DfwTTgTQ3VLKyZnjXwJOkVKqazQZnn7qGXvL5oyu/rTJtDRvZEbI65u/XbPs06sJWKzQPKu6MiIYk5omsml9q/L/ALNrNGFSvq/+Tpfqo56vae/qu+vw6+HrpKBHaVJUt6ec5aemaVNY37zBoyfvopN5Hh/2Dqs/TqvLWPnLa7tN3Y2tLVu1mvb+aUrCPwfqA/fWca/LXvXqb5PoUo+t0fjO+5f99t9vz+y2q9v/5lvv+0+jJuwdx+uulQ63nJvnAdR4C12Og4Ho6j9238P27xfeHag9p//+ouGlq+8ipfy2EGIh8DKQnyS+HlAnoXls2byNc8+8HIAlN13OuWderszxH9XN8VMB5/i9zg938dJzuXrBkkHP8Y99RU17lT/HX7zsQi6bf82g5/j/MGov5Zg7x19682IWnHGZfo6vyQt/6ti5nv3+5vhn/PR0bj7vlmBz/PjA5viu7n3QOf4CzRzf75/hn+PPaqxn3hGX+k/j70ZN9+xr5/ga3Dm+m+cB1JRtujn+A3+6K1D9LiPZuHcyMFVK+WOgG+e7rBBCzJNSPgMcCzxdqI6Q5RjF8rdP6H7DU6Y+Ua2cp4sSW6Vx4HEj+9K2TXeqV5+37i9qXf63TF2F+kYeG811xIgVYmy0jklh9Wau8r2pN0fUMn/tUm9498EStsKMitZoc+LpHjS/7XzXs7+1W33Q9OZ16DQ2vXaSrqTqzNKheQP70T18XOcZG5tEOqXVvdfJZemi7PzOOfmGvGW3/Ih5R1zKhrWPKOfttfc3PPsH105Tyjzepj4U3fwFrsEW1Pvh9aPGKOepv2BhRvJy3gPAXUKI54AocDbwLnC7EKIis/27AucbDGXJiHbgkVJ2Ad/U/NeRxW+OwfDpwm83GC4YBx6DoUTY2KRG8FDfYDAMkpE8x99l0thZo5W7fXDDHp4ye4VV7fg3Eqr00dyomnH2nu2OoTBpp2jr7dBG8NVEVfctf1SaLnpufU9OIixuJ1nfs41mtirluhNeA5nuB19eM0c59s3Y2wCk7BQ74l18dtRMpczFmrZfHPZar3VLTzp037HNZwmPpVSZLV2qbjdTbcpOszPepSSwBFX3HlS5LFCj7PJXGmrClfzdqOmKIQ/gvXe95qUD9/t3pcy02t2UY+2ZKD4bO2vN98uSHf+Ceh8tVY4UZsTO8Q0Gw+AYiNjmUGM6vsFQQoo11O8vk44Q4kzgOzjPmyuklA8Xqs8k1DAYSshQZNIRQowHzsDJYn0UcIsQouDcb0je+DWhKHNqmzzbr3b8zVPmg7AaGaebtz6WUB1V5tQ7Ths14Yrstp9NCdUBxa+ksyOh5pXbM2+uWRmKMrNmIrPDDUq5A1PeefjtqWalzFnpFuWYq6RjWRYV4QhrYqrL8dfjaruaarxzZ51s9n17545trrK5b+9ePv+G6jwzp9qbo/ADjW1lTftG9byG6QBUhaKI+ilsSbQrZXS57HQS2H7lnHwPvKSdYmuiQ+uc45/Tv7b6V0qZ7x90gXJsbSbiscKK0FQ1DoB1Ma/k+//spt5/AxHcttHLuusIkFDDk0lHCJF175VSbhNCzMkk15gO7JBSFvxgM9Q3GEpG8ACcAAk1+sykA5Dp9GcBi4Eb+vs8M9Q3GEqEE3Ib7F8ACmXSAUBKeSMwCfiCEOKLhSozb3yDoYQU0arfZyYdIYQAfgz8M5DAMf4VNB2Yjm8wlJCgc/wAFMykI4R4E3gRZ6DxmJTy2QJ1DU3Hj6WTvNez2bPdEPUaeKpCaj66nrTqSBLRxKvfd7ATcbahxs5u+znseTVO22887Eqo5y6xclGDNiGWWNXscaQa/bfzbe/+LRvUz4tq2p69MWxnWxcF54a/5lNheX+6Kk0I6ekyZ4T8bizMYtlAyFKj5W6f1+nZ/9KflCJZQ14+nSnneqWx6UzFtFLa/naCPpedIoGtkcvSRdn5nXN0hrxbV1ynHNv5ne8CsKnG5oEDHOeyg17yltnUquZpHIgYVzGDdPrLpCOlXIwzvw+EeeMbDCWkrF12DYbyZPiKbZqObzCUCBv9NG04MFh57dHAQ4ArgnaLlPK+UjXSYBipjOShflZeWwgxD0da+yHgeinlkkInutjYWRktdzvhk9XqDavGvU6Nsc2fMw7g1JWOgeeUfw1z2cpRHJNWI/12xFcqx/zGPX+EHcCpeV5sl9hJrkq0UfVn1cSj1KUxYOme/jvj3dn/2xnvJqTJuVcdUT9vQ8zrQ6bz3Fuf560Yt1OsT+xgZ0L1pPvDn2d79juSSvZlekOqJNiikBNh2UgFF1vTOSe2SimjM/i196ptcPPXubiaeJCJuox3ZOWyPHX5PP7WanIQuIa8fEbf7Wjnbdnckd3eNs0rFn1fvXrdT1aOFMAewUIcUsoHhRCuw3++vLYQQhyP89Y/W0oZTOXQYCgTPinN/CAMVl57CrBKSrlSCHEJMEZKeV5f5z791DP21lbnDTVl90lsWLdJWd/UvRnSmjekrlxtyFn7GTd1PG0t2xiFumy2MaX6u/vRvZHzl8lc+W7dW9mPTuZZd6WTGRnuadOm0Ny8QVuX7jv7Rxi6N0u+8m9j00Q2r2/VCphODXmFTjfZ6khL14ZJmTW3iqZxxNe30aKR+NZdKl3mYn/7o+Hcb+hKYIc0jqb+obRu+XCPGvXahGdMB3LS3QCvr1rjKbNbRI0pGLdnY2AZ7OW/edD+6bn9es8CcMMffzIi5LUPl1K6d+ly4OeFzt3a2sbCs34EwLU3XsrCs36kDPUrd2Gof3DddABO+clp3Hn+bdqh/qKO/of622OdSpl9Gpqy25f87DyuOvunVGnSQfvr2qJJsax7sGzpccrdets1fP+0C7UPlcqIem380s+6of6M2onZ7QuXnsM1C66nJaaGmVxTfYBn/8qEOtSv0KTXdof6U5b+OxsW/IpzizjUd9NaQU4C21XD9ZznEw1xA27ycdfp83GH9650N8BnD1noKXPqhEOU807+/fnKsUIMV6t+v776QoiThRAXZXZdee0HhBCfzRw7ClB7lcFQ5rjReUH+hprBymuvB24UQsSBzeRSbGmxsbNDX3fb/ybwJ6QASIbV4bLOQ+3F9g8A+GYqxovtH/BUUvX4i2mO+cNy72k4QimzKJ3VOiBt28TSCa0MVdKXOadaMypo6daEux64OwBbaypYc+DuHLx6m1LG/zYEqKzwXq9YSn2ruZ514Iw2OlMx7Zv7mqQ3Q40uC1BLp9r2hTXOKONaO8bC+DvK/4Ne6kseNFU5lk54y+UnsEzbjjyWbvTgH+n4Q2tB9ciDnCHv1tuvzb7ptzf/2VNmv32+pZw3IOMekBqm4lu7Iq99ePGbYzB8mvhk3uZBMA48BkOJMJp7BkOZUtZv/N1C1ZxW/xnP9s93eO2BdRqLbVQzH9XJcfnz8unktaOaY/4544Uax5U9qidktytCEZoqxzImpLa1Ce+SWLOtWq5132fee85KwpWxNP/yXqfWcq2b29b5Mu+2ahxXllTlpLrHWJUssWZySkq1vB9T53Xg2Wyry3KrNW3Id8pK2intfP7QOlVy+ytrVHtBm28VJH/VImRZVIcrtLns/BLYOrksXZSd65yzW6Q2a733z+lXv6M6o65aVzA/rIJ54xsMZYarwDMcMR3fYCgZ9sh12TUYDIPDtkdwdJ7BYBg8Q5hQYwHwr5ndRzOKPH0S2Fd/V1i5cuVWoBkc/fCGhgbVy6IIlLLuUtdv2v7J1D+IuqfNnTtXTcan4f5fP2Bf8l9XB6r014/fWtBXXwhxAvA1KeV3MmKbF0kpj8/830zgt8AhOKaFvwKnSylVS26GIXnj518oIcSKfvTDB00p6y51/abtn0z9pax7ANLZQegzoQaOJ+0/SilTAEKIKKAXn8xghvoGQwkJatwLkEmnz4QaUsoEsC2TNusnwOtSSm+ooQ/T8Q2GEhL0jR8gk07BhBpCiCrgTqADJ49eQT6Jjn9b/0WGZd2lrt+0/ZOpv2R129ikNLoMg6RQQg0L+APwlJTy2iCVDYlxz2AoR+779e/s888MJnW//E9392fcc636nyGTUAM4DlgLhIFfA/lxiBdJKV/sqz4z1DcYSkixXHb7S6gBqL7eBTAd32AoEbZd5kE6BkN5Ylx2DYayxATpGAxliBnqGwxlxohOoWUwGAaLmeMbDGWHseobDGWKkd4yGMoQ88Y3GMoMmxGcLddgMAyWosbjFxXT8Q2GEqLLDDwcMB3fYCgRRl7bYChHzHKewVCOmKSZBkPZsffes5546cU/jg9YvGQqxTqMAo/BUIaoWQ4NBsOnHtPxDYYyxHR8g6EMMR3fYChDTMc3GMoQ0/ENhjLEdHyDoQwxHd9gKENMxzcYypD/D+mR7/NFqhHiAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "%pylab inline\n", + "matshow(adjacencyMatrix(conj_mean_prob_graph).toarray())\n", + "colorbar()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing shortest-path mapping from a source region\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:mittens.spatial:Loading NIfTI Image cst_seed.nii.gz\n", + "INFO:mittens.spatial:Loading NIfTI Image cst_seed.nii.gz\n", + "INFO:mittens.voxel_graph:Computed shortest paths in 0.11079 sec\n" + ] + } + ], + "source": [ + "raw_scores, scores, path_lengths, backprop = mat_vox_graph.shortest_path_map(\"cst_seed.nii.gz\")\n", + "mat_vox_graph.save_nifti(raw_scores, \"raw_scores.nii.gz\")\n", + "mat_vox_graph.save_nifti(backprop, \"backprop.nii.gz\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualizing one of the backprop maps," + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAASkAAAD7CAYAAAA2GGISAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXl4HNWV6H/dWr3LOxIYzJYTJ4BJBMHsDju8EJgkJDwgz5gwDANfQnAySUhgEmf5MpDgyTJDwDAYJ5M9DklYA4PZwRg8DjsHTDDYWLZs2bJkW0uru94f1apbMlJ3S+5uVavPz19/PlV1q86t7tKpe+4999yY53kYhmFElfhwV8AwDCMTZqQMw4g0ZqQMw4g0ZqQMw4g0ZqQMw4g0ZqQMw4g0ZqQMw4g0ZqQMw4g0lUM5SUTiwE3AbKALuFRV1+SzYoZhGDBEIwWcC9Sq6tEiMge4EThnoMIPPbTca2pqpr5+Gk1NzUNUOXSGS2+56i7Hey6m7lmzDt7S2Ng4dajn68vPe8lUbjNNOrqTf21sbDxjqLrywVCN1HHA/QCqukJEjshUuKmpmXnzF7B0ySLmzV8wRJVDZ7j0lqvucrznYupeueKut/fk/GTK4+D68TmVfeHtbVP2RFc+iA1l7p6I3AYsU9X70tvvAAeoak9/5XtbUjNnzmDt2nV7VOGhMFx6y1V3Od5zMXXPmnXwqsbGxowNg0y88sJq7+C9xuZU9oV1bXukKx8MtSXVBowLbccHMlBgLaly012O91xM3StX3LXnF0kO+OcaOYY6uvckcBZAuk/qxbzVyDCMAuPheamcPlFgqC2pO4FTReQpIAbMz1+VDMMoKB6QioYByoUhGSlVTQGX57kuRkT43eS5gTzF6w7k+2urAvmGDY8Ws0pGXvEgIq2kXBhqS8owjFImlRzuGuSMGSnDKDc8z//kRKygVckFM1IjhF9O+Wggn3yYPwy+cVwtvw25bp9peSSna3XHQg9m6Fk+L7UzkMfVO33XNT08uMqGWNjgrnN4p3u7v1HtHs0vb1w+5Osb78XDw0vm2pKqyl6kwJiRMoxyoxw6zg3DKHGs49woNquq3UN36Bo/zjbRFWdN9eD7FHbGQ+cM8CyfW73NFQm5bN/ckN31C7t4R3W6oMKNFe5xTIWqcH7DUYH8mw3PZL2+kQ3POs4Nw4g41pIyDCOyeF5JTYsxI1XC3FB/UiDPO+SdQH7nhToAPC/GYV2Db9a3D3Ky1DGdLuDzE/VHUlc1hk/UH8kfm57tt/zbdAXygfGaQG4LuZnbYu5NfxCjAvnQyTP7XOvFlrWDq6yB7+5ZS8owjKjigedZn5RhGJElf9NismXpFZGrgfPTm/eq6kIRqQAWAUcANcC3VPXugXSYkSphzqhyI2xVU6sDeepeOwDYWpWkgsHnC3st7tyxvULumNDdX3FmTG4L5PO2TWaiF+e8xBj+OMD1m73OQN5cWRvIaytcP8lWLxHIH/Scu9dY09DnWi+ydgAtRkby5+4NmKVXRA4ALgSOwo/OelxE7gQ+DFSp6rEisjdwXiYFthCDYZQjXiq3T3b6ZOnFbx31sg44Q1WT6aQEVUAncDqwXkTuAW4FMibIspaUYZQbngepRPZyQGt76xQReS60a7GqLg5tjwe2h7aTIlKpqj2qmgC2iEgM+AGwWlVfF5EpwMHAx4ATgCXp//vFjFSJcXHDMYG836edK5dsci7arjbf9Usl45y19fFB6xgfeiz+FgoSndrpXLOGmh39njsxlaQy/f9AVIQmrbbE3fXfSLYHckfY3Ys7d29crO8je269e3H/qek5jBzwch/dq6ur26KqmdIHZ8zSKyK1wO1AO3BFencLcLeqesCjIvK+THUwd88wypH8uXsDZulNt6D+DDyvqv+kqr1vridC58wG3iED1pIyjLIjr3FS78nSKyILgDVABXAiUCMiZ6bLX4PfD/UzEVmRPidjAk0zUoZRbuQxC8IAWXpfC8m19M8lueowI1UC/GtoQu6XPuv6aio+ekogxze6FvPl/+NPwr0qxwUgd+fFVGsgz47XBfIr1S7Moa7DyROqOwJ5Wu0uemIpptXu4pJQ/9ntG54K5PaUC2VYX+Hkph4XyrBPldMbJrHbjGfrhxoKHl4yt47zKGBGyjDKEZtgbBhGZBnE6F4UMCNVApyR2BXI1ZcvDOTkuxrIqZcfCOTHNr0MwOcSzg0bDA9ufD6Q9w65bBNibjA4F0/yi5VO/zc++IFA/t1Wt3L3/+JSEu8bcvGe3bE2kB9LuAj19q6h3ZOxGyOtJSUiRwHXq+pcETkIuAO/++0l4Mp055lhGCVBabWkssZJichXgNtwvfSLgGtV9Xj84cNzClc9wzDyjkc+46QKTi4tqTeBTwC/SG83Ar0rQ94HnIYfK2HkkTsnnRjIh/369ECO1Yxx8ugJgfyjX7qJwPlk+Y5gQjv7jv1g1vJjxnaxoyLFmLFdVNe6qPPuzopA3hlzvuKMUK6oh7rdCGV3KClb2MVLReQPp7TxoKd0kt5lbUmp6jIgPF4ZS4ezgx/qPuG9ZxmGEVlGYEtqd8I1Hwe0DlSwl/r6aSxdsoiZM2ewdMmiIajcM4ZL757onlrppkP9fdT4QI5tdPPbvIR7PxxzwwWBvDRx9h7pDlNd4dZdmxx3cXnx0AO8o8LJsRgkG6axY+GV9Fm+L9TRfkIytOBCSNec0JLuiZR70/f06T/J3GNfir918SmtPqmhGKnVIjJXVR8BzgSyLg/S1NTMvPkLWLpkEfPmLxiCyj1juPTuie7WLxwZyLUn9//gJzc2B/JhF/1j3nQPxBV7Hx/IZ3Y4C3TwZJfXKhbz2LHwSsZ+8z/7uHsdO5yxe6BzYiAnQq7ffbveDOS32zcHcldP/3ms+qMUf+vBsnJFxswmuRGRVlIuDMVIfQm4VUSqgVeBP+S3SoZhFJSRuDioqq4F5qTl1/EnDRqGUZLkL31wMbBgzmHmxr3cii+XfdX1MVV9JrvbkFpxf0HqNBB/aXs1kCvGu+DMmd1u5K5mlN+XFItDLO5cuXD/VNjFS4T6mHb1uJxYg3HxjEHildbonhkpwyhHvKFNPh8OzEgZRjky0vqkjMJxUWhRz6rPLMla3utyc91u/3ZzhpL5Z337lkBuH+fchdYdLjRh+iiXVjhe4d7WdTPc/MMpOjmQn/DcyGB3qnRckJLGJhgbhhF5rOPcMIzo4kHSVjA2cmTUyQcPqnznt/8lkK/a+GKGkoXljlCmzUNCI5Rnp3x3b/cXdcyW/IgOIzFOyjCMkYT1SRmGEXWsT8rIROeGwS/Y2cv65RXZCxWZL29cHsgneEdCMs727aOoHeNG68zdixAeeENcpGM4MCNlGGWHuXuGYUQZDxvdM/JL4vc/CuRru6Pn7oX5tzjMA5bG4Qfdg/PxmndmTU1m5AVrSRmGEWUsBMEwjMhjE4yNfLLt5y8F8p+bNEPJ4eePTc9yTuL/8semZ1lYfXiwfyyWeiU6mLtnGEaU8chtddeIYEbKMMoRG90zdmefcVOyFwoRTsny7abwudF298K07BgdyNPYkaGkUVQ8D8/cPcMwIo25e4ZhRJo8zd0TkThwEzAb6AIuVdU1oeNXA+enN+9V1YWhY+8HngGmq2rnQDrMSBWJj41/fyAnm98K5Ipp+/dbPvGbHwfy9w7ZFMi3bShA5QrEGxUuY+es0P7HvK2B/ODG5wM5FlqtwSuhIfKSw/Py2ZI6F6hV1aNFZA5wI3AOgIgcAFwIHIXfXf+4iNypqi+IyPh02a4Brhtg0z4NoxzpSeb2yc5xwP0AqroCOCJ0bB1whqomVTUFVAGdIhIDFgNfB3aRhYwtKRGpAm4HZgI1wHeBV4A78C3jS8CV6QoYhlEKeLmvu9fa2jpFRJ4L7VqsqotD2+OB7aHtpIhUqmqPqiaALWmj9ANgtaq+LiLfAu5R1edFJGsdsrl7FwEtqvpZEZkMrAb+Blyrqo+IyM34Tbs7s2oqMw6dPBOAUZU1HDp5Jl+oaQuODeTihel46I1A/sYr00NH1ry3cJH40JQDA3n1ljf7LTN78v6Mqqxh9uT9eb4qEez/eKjMOZ5biOHg0NLtN73bfwqbQybtF8gvbX17sNU2+iNHd6+urm6Lqh6RoUgbMC60HVfVIEePiNTiN3TagSvSuy8C1ovI54C9gAeAEwZSkM1I/Z6+y6j3AI3Ao+nt+4DTMCNlGCWD35DKm/PzJHA28Lt0n1SQ0zrdgvozsFxVr+/dr6oHhcqsxbchAxLLpYNSRMYBfwFuBX6oqg3p/ScBl6jqRZnOf+ih5V5TUzMzZ85g7dp1WfXlm+HQO6qyBoD6GdNpWreJAypCq/y+f9+s5yffWhvI73a4d8mWRO7xRvm+79Hpe4K+qw2HGVVZE9zz2Fh1sL9hlGtVbetw+ztCKxhvTrQPeM2g/AB6exmuZ6yYumfNOnhVY2NjptZNRl568G5v5v/cnlNZ/fQ3MuoKje4dBsSA+cBZ+E3+CuDXwIrQKdeo6tOh89cC79+j0T0RmYHfUrpJVX8lIjeEDo8DsubXaGpqZt78BSxdsoh587MvH55viqV30ijX6r15dKMv/OhC+OLvOeTF7wzqWo+c8utA/km1++N9ptW5ez0p17G5o7sjkI+b5i+BfvWNn+eyf/xasD+fS5ePrR7Vr+5/3+skpvzwAt79l2WcEHMubvUop3tb68RAfsINAPK25/pQlzU9G8j7T9grkD80ekafejzfsT6Q32xtGrZnDIr3nK1ccdceXiF/o3vp/ujLd9v9WkiuJQOqOjObjoyjeyIyHd9f/Kqq9pre1SIyNy2fCQw9F65hGMWnN+ldLp8IkK0l9XVgInCdiFyX3ncV8BMRqQZepW+flWEYJcCIyXGuqlfhG6XdObEw1TEMo+DkN5iz4FjEeR7Z2uH6jj7+pt8H9cbGdo4cZH8UQFOF+2nOT7o+nPPHHdlv+XPPdVHcn7vb75iujMX3uB/qiKlu8dL/jLl6jB3jOrC/0VkVyMd5bXgkOY42JtS5vqqtW91k4/fXbQvkNR1u8vTGmBtcOGn6oYG8fJNbBHVydXi0G1o6++9sN7JgE4wNw4gslk/KMIzIY0aqPLlzUv666iYkB9ccf3hZXSAv2tefkLy1Ksl36j8a7L832RTITzeHR4kHZm5VQyB3dzkXb8Y/u1CAn/zWRcdX1KRorkky7YDtdLW5x2tXj3MJV+x0LtvHxmwJ5He6JgTyxlAdDp9yQCCvaW8iTGun5akaNJ6HN8jnazgxI2UY5Yi1pAzDiDIjJgTBGByf2PZYIHdkKJcLi6vdxPJ3ut3I3bZu596Ep6bMnjDTndzkj5JdlYjxYKo52F0VGj3bdFowfYobX9q7j+65nc4VOGCM0z3xA+6uuh5rCeSaie6BT+yIgwdeKsaRf3832N+y69VA7o2IBzgl9Ah+qtMFD95Y7abONCXcpIYDx7roc4BVncM34bpksY5zwzAiT+l0SZmRMoyyw/PwekrHSpmRyiP5THl7/8a/Dar8Y50vB3LvhNzuVA/rOlv6Ld+szm26ambfnMS1M5xbWDHZBWHGQpOKve4gZRAv/c65Zpck1/HtzhSf0Z207HITjMNMrxgT2uo/q8GWpFstpyI0xXTVFnPv8kLp2CgzUoZRbvj5pKxPyjCMKGMtKWNP6dzgMuDUNhyfoeR7eWu7HwrZnUyQSPX0W+bMrS4P04OVDX2OxSvdKN6Nj7vAywle/5l9rtvycJ/trmQ3r297t9+yANPj4RRD/bt71TH3aL687Z0Br2UMAWtJGYYReawlZRhGZPHA67+BHUnMSBWItksvASB56dW0XfvvjL8tt5zSvQzWxbumYW4gf3/DI4G8vn3L7kUBWDLVzen77W6Rp21r3Gv2wS63OsuUyrGBvLVnJ0NlgleRtcyYmJvrN6oqlFu9e8BU2MYgyNMCxkXBjJRhlBse5u4ZhhFtrCVlMO1eP33J0vM6OfzeN9jx5LLgWOWxn8x6/n9NOymQP9e8PGv5e7v6XzQzvNJK76gfwKNVblTtjg1PDXjd8Dy7J7a49C6J5NA7NWqIZS1zQigL6KpQEKm5e3uOhxkpwzCijAdeMvuLIiqYkTKMMsRaUsZ7qd9/UMXPPd4FQ35uWYaCaZ5veavf/WEXL0wmFy/ME82v5FRuMIxNDe4tXl/rXL9NO7ZlKGnkRDqVTqlgRsowypAR1ZISkQrgVkCAJP5a7zHgDvw+uJeAK9PLLRuGEXE8wPNGVkvqbABVPTa9vPoifCN1rao+IiI3A+cAdxasliOAcJqTXKi57EK3sexf81yb4WXiIF9nlbHswZ/GIPBKqyXV/4zREKr6J+Cy9OZ+wCagEXg0ve8+4JSC1M4wjAIQI5XM7RMFYrkmahORpcA/AJ8C7lDVhvT+k4BLVPWigc596KHlXlNTMzNnzmDt2nV5qPbgGC69Yd0fmrVfsC+XVpW3y63Mu/qNgTMK5KK72GTTu1/l+EAeW+HirRJJ987cGXPyNs+twryrJ3OcVBR+60Iza9bBqxobG48Y6vnPL/uzV/e9n+ZUdsut1++RrnyQc8e5qs4Tka8CzwDhv7JxQGv/Z/k0NTUzb/4Cli5ZxLz5C4ZW0z1guPSGde948LvBvsoP5rI+n0uRsu8vvh/IdbfknrEzqt93eH3CaaEl15u3u7mBz1W7dC5/TrrMoc+19M3MuXtQaRR+60KzcsVde3S+H8wZjVZSLuTScf5ZYB9V/T6wC3/Wz3MiMldVHwHOBB7OcAnDMCJGvjJdi0gcuAmYjZ8c7FJVXRM6fjVwfnrzXlVdKCITgP8GxgPVwAJVfXogHVn7pIA/Ah8SkceAvwJfBK4EForI02klfxjszRmGMUyk46Ry+eTAuUCtqh4NfA24sfeAiBwAXAgcAxwNnCYihwELgIdU9UTgYuA/MynI2pJS1Z3Ap/s5lL81xUcgFzTMAWBS1VguaJiDtzGUXfKDg7tW9Revcxu3ZJ/3F3UmxbqzFwoxKpS2JZbDvD8jO3kMQTgOuB9AVVeISLj/ah1whqomAUSkCugE/h2XkrUyvW9ALJjTMMoMz4uRzHHkrrW1dYqIPBfatVhVF4e2xwPbQ9tJEalU1R5VTQBbRCQG/ABYraqv9xYUkb3w3b4vZqqDGSnDKENybUnV1dVtUdVMo3tthEd5IK6qwWiGiNQCtwPtwBWh/YcCvwG+rKqPkgEzUnlkQq1bT+5wz1+vbjRxDvdG473V/9y6XIhPmBbIbd8/M5DHX3PfkK9ZbJZPOiaQR1cNzt2rDgVzdicTeatTOZPH0b0n8QO+fycic4AXew+kW1B/Bpar6vWh/R8Afg98RlWfz6bAjJRhlBte/kb38GeanCoiT+HPRJkvIguANUAFft91jYj0vl2vwe9grwV+LCIA21X1nIEUmJEyjDIjn3FS6Tm7l++2+7WQXMt7GdAg9YcZqTwyfVRdINek31SxtFz5qcv6P2mQxD94WGgr2u7eQ5OODeTaUGR5PPT3UVGZfRLZYJecN7IRI5nKJfooGpiRMoxyI7/uXsExI2UYZYYHpEZYqhYDOKjOLUW+pnVDv2U27NwayPO/4q8Vt3ZajPlX1RCv26vfcwZLbJ9Zgbz9Gy6edsL3Mo7iDgsnb30ykJ+Z9pFArqpKBnJ1bZL++FpT9sUnjKEz0vJJGYYxwjB3zzCMyDKYzJxRaG+ZkTKMcsPLfXQvCgYiCnWIFOHFNMOMqqgO5JOmHxrI3066KPPDLnchIVWf/jwAsW3JQO6lZ60bUq+cefjgKljhJtvG990nkDeedFAg77W8b86lKDCq1kWKV9e4cIRkwv2x7DXJJfp7bZL7jsO0tfUNu6mIO7/l7a4xTKscxx8nncgntkavjy4qeOlPqWBGyjDKEBvdMwwj0tjoXgmzvXtnIE+scelsK0I5t3862sn7LpjhypzmVniJjUpPDN/eTmzUOBK/XRQc6/hrMAeTqnqXibnqY6cGcuXR/9Bv/WKhSczEXT2q967p/4YiwmHvrA5kfd8hgVw7YXAThsdO7hrwWOOK51jacw7zzNXLiIefXrdUMCNlGGWIF4lxu9wwI2UYZYYH9Ji7V7rs6HaZTH87yk3m3W+KSz5Yf7VzV6rO3n0CuI/Xu4qJ5+Ele/jFv7nzN1S4UbnLQ9Hrdec59zLVutFdLORqpraE0hCH3L2KhomB/NrB/sjYjprBLUhaSD5Zf2Qgf7ndORudbdkX/nxwY9aUQ8agiFlLyjCMaGN9UoZhRBprSZUwP5pyfCBv73Hvm66O0FfV7dLfJu66OZBTa1yK4NikCQB4jZ8i8cAf+Ow3XQrgipPPD+TkA78MZE/dqF989snuum1bnO7QysakXP1SW93+ZE+vGzi8D+JnGo4K5IpQXZKhiWM9nptgvDnh7uGlrW8XuHbly4gc3RORacAq4FSgB7gD/15fAq5MZ+czDKNESJZQSyrrBJ70Wlm3AB3pXYuAa1X1ePxX9aBSgRqGMbx4QCqW2ycK5NKS+iFwM34CdYBGoDda7j7gNPxk7COCNyrdvLJdMff1HDHRjfrtWOYCE1vXjQ7kqYe4MrHKFgC8WV0knvs7zz88OTj29Hd+HMgX778+kCfeESz+2mdEj57Q6ipdLtiU7W7EMNXqdK/a4ee+2i9Z0Wee4fJNzp0sFBc2zGFS1VgubJhDVegdmAzNFtsv7r6zlz13bzewbyCfhbl7hSRVQi2pmJchsYyIXAzso6rfFZFH8BOuL1fVhvTxk4BLVPWiTEoeemi519TUzMyZM1i7dl3eKp8rg9E7rcotIVYd+iEnhpZhCq2wRLLb/SFWjnJebyx9avfk6VS3bGJXuzN4O2PuupNrXMR1xUwXmhCrcOW9PkZql5M7nWHydjq5td2fhFwzYzIb1m8K9rcnOig0k6rGMmmfyWxd39LnzyD8lFWGjnTg+qSme+6LXdMT6nsbBMP1jBVT96xZB69qbGzMtBZeRp7+zT3epi/+MntBYMY9X9ojXfkgW0vqEsATkVOAw4GfA9NCx8cBrdmUNDU1M2/+ApYuWcS8+QuGXNmhMhi9V+99QiDvnXRfz6f2cvFMNXWutTVwS8r/Q1z/2avY5xc/7tOSeqbGTWEZqCUVH+uMZarFlUmFYoa8118N5O5nNJCffyTdklp0Af/x5ZuC/cVqSZ1/wyX85iu3D9iSmhRaNv3lpHt8FnRPCOR5Wx8fkv7hesaKqXvlirv26PwR1XGuqsFfbKgl9QMRmauqjwBnAg8XsoLFZqvnWjZH9Lg3++c3uTlzx29yf0wfCQV/xl91P33vBM7ErhibX6jm8VpnmL65wX1l32gKKX/fxwPxlmknBfJF91zgrrvDtTBSLdsCuXuza5F0pFtqHjE6vMIvpnnF3m5ENI4/kjcxZIgAJnnuUdsQcvH+Z+MLTi5cFY3dSMVKx90bSgjCl4BbRaQaeBX4Q36rZBhGIfGA/jPLR5OcjZSqzg1tnjhQOcMwok9URu5ywYI56ZtpMxbq1L1ix8pA/sB4l5JlXZVz/a7veCmQW15ve8+1l3Z1MO+NwfcF/VOzWy3ltM85n3DK510fZmqL0+f1uHo39Ph9ZlV4nFFZH+x/us/CsnvGVxvmBnJNaLJqEo9KYIJXwQTP9Uk1xV0/3i+bn8tbPYzB4xErqdE9M1KGUYZY+mDDMCKNuXsRYtqYOirjFUwbU0dNaBGDydVuiL8qFPjUnHKxRNs7XeBkc60LnLy5xS10kEg6N6ZQLGt2Lts/dbnMlLEadz819c61OqRtMwDbqnroDL0zz65vDOS7mlblpHvKaDeSuWWX+w72SoXDCxxjUzGqvBgzeuK8VelGO7d57nvqSAycXdMoPCMqBMEwjBFIDJLWkjIMI6pYSyoCjK52a7NVxyuJxWJUxyupiTv3KDyK15p0Lt4zm13kdpg3W5v63V8MvrzRjfR99HtudO/Ac52jVTHJZeGsrfPd1Filx+ldLnByWpVz3eoajg7kX2x4uo++j0x9XyC3hNKnfH+scxdDWWIYm3Iu5ZpqODDmsanSYzNO9682PTvQ7RnDgBkpwzAiTb5SnItIHLgJmA10AZeq6prQ8auB3gRq96rqQhEZBfw3/hS7dmCeqm4eSEduay0bhjFi6HX3cvnkwLlAraoeDXwNCCagisgBwIXAMcDRwGkichjwz8CL6XRPPweuzaRgxLekEl4Sz/P/Hx2ar9SVcnPaunByTaVbTr0rnH0gIlRWDm5Cw8TRzpU9ibDsOGp6eAumhjKSTo2772BXT/8jmVrtHqOmWDcJPJpiCTaG3OjuZOHnEBq5k8dpMccB9wOo6goRCWdMWAecoapJCHLTdabPuSFd5j7gukwKRryRMgyjL71J73KhtbV1ioiEpwgsVtXFoe3xwPbQdlJEKlW1R1UTwBYRiQE/AFar6usiEj6nHZhABsxIGUYZkmvHeV1d3RZVzZRPqg0/ZVMvcVUNmtwiUgvcjm+MrujnnKzpnkaMkYqHMll2JpyL4nkeSS9FW9cuquPuditDAZxvbQ+tcRdxmne6/FX7pXZmKOkzPpRRdCD+DwOPXFbWuMe56Z3xgfxG0q0RGI65aU51kcCjOdXFKx0uB5cRLfI4uvckcDbwOxGZAwQTVdMtqD/jJ8q8frdzzgJW4qd7ypg8bMQYKcMwcsMjr3P37gROFZGn8Nc8mC8iC4A1QAV+xpQaETkzXf4a4GfAUhF5AugGLnjvZR1mpAyjDMnX3L30SlG7L+MdTrdRS/+cl6uOEWmkUp5rzPakkoBHTypJd8qNTvWUVNovx5RRbsSs660dgVyz/9j+irOr3Y1Wrm6fFMivVrt36d7JvpEoFaHX7DGjtg65rkmvlEIGy4cRm/TOMIyRgkeqhJK1mJEyjDKklNq4I8ZIpQZwLRLJHjzPI5HsYdOObf2WKSVWJuoC+cMvODft4P2zp4ypT7lRz58m3NLtB1VN7FPuhJ7RDJU2r4skHm1eF3PG7h/sbw6leYlikGw5keeO84IzYoyUYRi5Yy0pwzAiiwf0xEqnLWVGqgSoH+tG5baEht4a9g/PRhhDf1SEsmNWxpz8bmeLOzNeTV+G7u6t69pKd6qHdV1++r2LAAAMT0lEQVRbaat012kY4+6hlIJnRyqlY6LMSBlG2TEik96JyGrchMC3gFuAHwM9wAOqurAw1TMMoxCMqBCE9ATBPouDisjfgE8CfwfuEZEPq+r/FqqShmHkl9IxUbm1pGYDo0XkgXT5bwE1qvomgIj8FTgZMCNVICZWu2jyj4TSAf9dJweyjHchBRVj/TkPsd1SGlbFXSP/S2MOC+Svbu47v/Mzk6f2Ww9vgHSOx/e4yPeLj6plw9g4jx9VS+dmV/5P6z8QyF+wPqlhZSS6e7uAHwK3AQfjJ6kKp1ZoBw7IdIH6+mksXbKImTNnsHTJoqHWdcgMl9586a4NJeIbg5PDNmjDmPfGSSWmTKdn4RXBdkVo+svsUBaI23rO7nPevhUu80ZnaOXh8ISvGSHt8dAjv2FMksSU6Wy4bAGhvIIcknD1Xpr42Hvqmi9K/bcuFskSakvlYqReB9aoqge8LiLbgUmh41nzwTQ1NTNv/gKWLlnEvPkLhl7bITJcevOl+wOT9g3kH+OWe6+NuRlYBxzpWlLEfWOy4dKrqf3mz4LdW1vdaNvzFW40cPeW1E8nHx/I4bl7uzqcoVmXcucfEHMpYxoat7LhsgU0LF5E52ZnCJ9Yv3cgf2GTW1gi35T6b50LK1fctUfnj8SW1CXAocAVItKAPz69U0QOxO+TOh2wjvMC8srWdwJ5Zf2BgXxct8sV1b3NtWxi6TAFLwnxiv7fmDMS7jHdPbXvfj1uO+zihRaF6fMeroi7rS2vj6anM86W10ezvd1NgH+yKnvuK6N4eCOsJfVfwB3p3C8evtFKAb/EzxfzgKo+U7gqGoaRb0ZUS0pVB0pKNSf/1TEMo9B4lgXBKCTfaHo4kP8weW4gv/V31+cj+K5VvLOCmm7nBiZDHd/vVrn9FzT0fd/8IrR6zlVd7hHpTrk+pn3ocrp7XP9Uy85xzEpVsnLnFL7T8bdgf9PmoeelMvJP6ZgoM1KGUXZ4QE8JmSkzUoZRhoy0jnMjonyq5ZFAXtjw0UBuS6/kckisglfaXcBnW5Vz9/btduEL+1IzoI53c6jHppDr+K22Z/lJ8hyubXuWrR3tOZxtDAcjquPcMIyRh7WkDMOILCMxmNMoAb65wY363VB/EuA/jBtCv/B+Cff2XFflRuoeiPedMLBXbFQgT6cqkMeGAjurQvKve1yw6daOdpKppLl6ESfpWUvKMIyIYnFShmFEHuuTMoaVrzT5E3iXJj7GyzG3mOhb1c7F2+i5/TH6pmAJx9BsC01iXo/LiPBhz7mEb+/anIdaG8XE+qQMw4gsfse5taQMw4gw5u4ZkeFXG1YU5LpLC3JVoxh42OieYRiRxkb3DMOIONZxbhhGZPHIX5+UiMSBm/AXbOkCLlXVNbuVmQo8BRyqqp0iMgH4Df6Ktt3ARao64Ooc8YEOGIYxckmlXb5snxw4F6hV1aOBrwE3hg+KyOnAA8D00O6LgRdV9QTgt8C/ZFJgRsowyhDP83L65MBxwP0AqroCOGK34yngFCCc9fBF/AVcAMYDCTJg7p5hlBkeuS9p1draOkVEngvtWqyqi0Pb43GrmwMkRaRSVXsAVPVBABEJX7YFOE1EXsFfeep4MmBGyjDKjtxH9+rq6rao6u6tozBtuFYRQLzXQGXgm8ANqnqLiBwGLAMOG6iwuXuGUYbk0d17EjgLQETm4Lty2diGa30147fGBsRaUoZRZuR5WsydwKki8hQQA+aLyAL8BYX/MsA51wG3icgVQBXwj5kU5GSkROQa4ONANf5w46PAHfj3+xJwpaqWUuiFYZQ1+QpBSP/dX77b7tf6KTczJG8g3frKhazunojMBY4BjgVOBGYAi4BrVfV4fOt5Tq4KDcMYXjzPnxaTyycK5NIndTq+n3kncBdwN9CI35oCuA9/iNEwjJIgtxipqEydycXdmwLsB3wM2B/4C34Pfu8dtAMTClM9wzAKQVQMUC7kYqRagNfSy62riHTiu3y9jANa+z0zTX39NJYuWcTMmTNYumTR0Gs7RIZLb7nqLsd7Hm7dg8GDXEfuIkEuRuoJ4CoRWQTU48+3eUhE5qrqI8CZwMMZzqepqZl58xewdMki5s1fsKd1HjTDpbdcdZfjPRdT98oVd+3xNUZUS0pV7xaRE4CV+H1YVwJvAbeKSDXwKvCHgtbSMIy8MuKS3qnqV/rZfWKe62IYRhHw8Eh6pRMxZMGchlGGjLQ+KcMwRhgjqk/KMIyRRT6T3hUDM1KGUXZ4pMzdMwwjylhLyjCMyOLP3bPRPcMwIoy5e4ZhRBYv/a9UMCNlGGWItaQMw4g01pIyDCOy+NNiksNdjZwxI2UYZYhNizEMI9LYtBjDMCKL51lLyjCMSGPTYgzDiDg2umcYRmTxsGkxhmFEHOuTMgwjwliflGEYEcZG9wzDiDwWJ2UYRqSxlpRhGJFlxC1pJSIXAxenN2uBw4G5wI+BHuABVV1YmOoZhlEISqnjPJ6tgKreoapzVXUusAr4AnAzcAFwHHCUiHy4oLU0DCOveJ6X0ycK5OzuicgRwAeBa4CrVfXN9P6/AicD/1uQGhqGkVfyuaSViMSBm4DZQBdwqaqu2a3MVOAp4FBV7RSRCmARcARQA3xLVe8eSMdg+qS+DiwExgNtof3twAGZTqyvn8bSJYuYOXMGS5csGoTK/DBcestVdzne83DrHhT5DUE4F6hV1aNFZA5wI3BO70EROR34N2B66JzPAlWqeqyI7A2cl0lBTkZKROqA96vqwyIyHhgXOjwOaM10flNTM/PmL2DpkkXMm78gF5V5Zbj0lqvucrznYupeueKuPbxCXoM5jwPuB1DVFWmPK0wKOAW/q6iX04EXReQeIAZ8PpOCWC4WVUQ+Dpyiql9Ib/8N+CTwd+AeYKGqPjPQ+atWrdoMvJ1VkWEYubBfY2Pj1KGevGrVqvuBKbmU3bRpU+3VV1/dGdq1WFUX926IyG3AMlW9L739DnCAqvaEryMia/EbOp0i8iCwHrgEOAH4jqqeMFAdcnX3BN8g9XI58EugAn90b0ADBbAnX6hhGPmlsbHxjMGUP+usszIdbqOvZxXf3UD1Qwtwt6p6wKMi8r5MhXMyUqr6g922VwBzcjnXMIwRzZPA2cDv0n1SL+ZwzhPAWcAyEZkNvJOpsAVzGoaxJ9wJnCoiT+H3L80XkQXAGlX9ywDn3Ar8TERWpM+5PJOCnPqkDMMwhouswZyGYRjDiRkpwzAiTcH7pHKJSC2Q3qOA61V1rogcBNyBH2z7EnClquZ9hqWIVAG3AzPxI2m/C7xSJN0V+L6+AElgPr6/X3Ddaf3T8GNhTsWf01ksvauB7enNt4BbKMK8UhG5Bvg4UI3/fD9Kke653ChGSyqISAW+hh+RWlBE5CvAbfgTosEPwb9WVY/H/8M9Z6Bz95CLgJa0njOB/yii7rMBVPVY4F/TeouiO22cbwE60ruKpbcWoHduqarOpwjzSkVkLnAMcCxwIjCD4v3OZUcxjFSfiFT8+TqF5k3gE6HtRvw3HcB9+BGwheD3wHWh7Z5i6VbVPwGXpTf3AzYVSzfwQ3zjsCG9XSy9s4HRIvKAiCwXkROAGlV9Mx2D0zuvNN+cjj/UfidwF3A3xbvnsqMYRmo8rjkOkBSRgrqZqroMSIR2xdIPLfhzDScUSO8OVW0XkXHAH4Bri6U7rb9HRJYCP03rL7judCqfzar619DuYt3zLnwDeTr+MPaS9L5eCqV7Cv7L9jxcYHO8WL9zuVEMIzWUiNR8E+4byDrXcE8QkRnAw8AvVPVXxdQNoKrzgPfh90+NKoLuS/DjZB7BzzX2c2BaEfQCvA78t6p6qvo6/stwUhF0twB/VdVuVVWgk75GqeC/czlRDCP1JH50KYOISM03q9P9COD3FT1eCCUiMh14APiqqt5eZN2fTXfmgt+aSAHPFVq3qp6gqiem8439Dfh/wH3FuGd8A3kjgIg0AKOBnSJyoIjE8FtYhdD9BHCGiMTSescADxXpnsuOYkScvycitQg6d+dLwK0iUg28iu8KFYKvAxOB60Skt2/qKuAnRdD9R2CJiDwGVAFfTOsrxn3vTrG+7/8C7hCRJ/BH1S7BN845zysdCqp6d7r/ayX+i/5K/JHF4fiuRzwWcW4YRqSxYE7DMCKNGSnDMCKNGSnDMCKNGSnDMCKNGSnDMCKNGSnDMCKNGSnDMCKNGSnDMCLN/wdRDEyKD+0GFAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import nibabel as nb\n", + "bp_data = nb.load(\"backprop.nii.gz\").get_fdata()\n", + "imshow(bp_data[:,40,::-1].T, vmin=.16, vmax=0.27)\n", + "colorbar()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/testing/create_test_data.sh b/testing/create_test_data.sh index 5ce20a6..8e6c6f7 100644 --- a/testing/create_test_data.sh +++ b/testing/create_test_data.sh @@ -1,77 +1,18 @@ #!/bin/bash -root=`pwd` -test_dirs="test_data DSI-Studio MRTRIX3 Dipy" - -for test_dir in ${test_dirs} -do - - if [ -d ${test_dir} ]; then - rm -rf ${test_dir} - fi - mkdir ${test_dir} - -done - -wget 'http://tractometer.org/downloads/downloads/ismrm_challenge_2015/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2.zip' - -unzip ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2.zip -rm ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2.zip -bval=${root}/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2/NoArtifacts_Relaxation.bvals -bvec=${root}/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2/NoArtifacts_Relaxation.bvecs -dwi=${root}/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2/NoArtifacts_Relaxation.nii.gz - -# First use mrtrix to make some masks -dwi_mif=${root}/MRTRIX3/dwi.mif -mask_mif=${root}/MRTRIX3/dwi_mask.mif -mask=${root}/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2/dwi_mask.nii.gz -dwi_las=${root}/ISMRM_2015_Tracto_challenge_ground_truth_dwi_v2/dwi_las.nii.gz -mrconvert -strides -1,2,-3,4 ${dwi} ${dwi_las} -mrconvert ${dwi_las} -fslgrad ${bvec} ${bval} ${dwi_mif} -dwi2mask ${dwi_mif} ${mask_mif} -mrconvert ${mask_mif} ${mask} - -process_dsi_studio () { - - cd DSI-Studio - SRCGZ=`pwd`/dwi.src.gz - dsi_studio \ - --action=src \ - --mask=${mask} \ - --source=${dwi} \ - --bval=${bval} \ - --bvec=${bvec} \ - --output=${SRCGZ} - - dsi_studio \ - --action=rec \ - --method=4 \ - --mask=${mask} \ - --csf_calibration=1 \ - --record_odf=1 \ - --deconvolution=1 \ - --param2=0.5 \ - --source=${SRCGZ} - - cd ${root} -} -#FIBGZ=${root}/DSI-Studio/dwi.src.gz.odf8.f5.bal.csfc.de0.5.rdi.gqi.1.2.fib.gz - - -process_mrtrix () { - - cd MRTRIX3 - dwi2response tournier ${dwi_mif} dwi_response.txt - dwi2fod csd ${dwi_mif} dwi_response.txt dwi_fod.mif -mask ${mask_mif} - #tckgen dwi_fod.mif dwi_fod.tck -seed_image ${mask_mif} -mask ${mask_mif} -select 10000 - #mrview ${dwi_mif} -tractography.load dwi_fod.tck - -} - - - -process_dsi_studio -process_mrtrix - - - - +rm *nii* *mif *fib* +wget 'https://upenn.box.com/shared/static/25h8ozw929xhn1wmtj4wndrd29ym9xb7.tck' +mv 25h8ozw929xhn1wmtj4wndrd29ym9xb7.tck ismrm2015.tck +tckmap -vox 5 -tod 6 ismrm2015.tck TOD_5mm.nii.gz +3dresample -orient RAI -inset TOD_5mm.nii.gz -prefix TOD_5mm_LPS+.nii.gz +mrconvert TOD_5mm_LPS+.nii.gz TOD_5mm_LPS+.mif + + +3dcalc -a power_5mm.nii -expr 'step(a)' -prefix mask_5mm.nii.gz +tckmap -vox 2 -tod 6 ismrm2015.tck TOD_2mm.nii.gz +3dresample -orient RAI -inset TOD_2mm.nii.gz -prefix TOD_2mm_LPS+.nii.gz +mrconvert TOD_2mm_LPS+.nii.gz TOD_2mm_LPS+.mif + +CONVERT=/Users/mcieslak/projects/qsiprep/qsiprep/cli/convertODFs.py +python $CONVERT \ + --mif TOD_2mm_LPS+.mif \ + --fib TOD_2mm.fib