From 5fd226b1b6b7c0db72b324e0ac776166896f4f38 Mon Sep 17 00:00:00 2001 From: Mathieu Boudreau Date: Wed, 7 Feb 2024 12:24:19 -0400 Subject: [PATCH] File overwrite didn't work well in previous copy; fixed it now --- binder/requirements.txt | 5 + content/index.ipynb | 1206 +++++++++++++++++++++++++++++++-------- 2 files changed, 982 insertions(+), 229 deletions(-) diff --git a/binder/requirements.txt b/binder/requirements.txt index cd3607c..5e41401 100644 --- a/binder/requirements.txt +++ b/binder/requirements.txt @@ -1,3 +1,7 @@ +jupyter-book==0.13.0 +plotly +seaborn +markdown datalad matplotlib nibabel @@ -6,3 +10,4 @@ numpy pandas scipy statsmodels +repo2data \ No newline at end of file diff --git a/content/index.ipynb b/content/index.ipynb index d176545..d4b3e82 100644 --- a/content/index.ipynb +++ b/content/index.ipynb @@ -1,17 +1,126 @@ { "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [ + "remove_output", + "remove_input" + ] + }, + "outputs": [], + "source": [ + "build = 'archive' # 'archive' uses the neurolibre archive of the data., 'latest' would download the latest versions of the data.\n", + "notebook = 'neurolibre-figures' # neurolibre-figures: no processing, use downloaded processed data and do figures & stats,\n", + " # neurolibre-clean: if some processed data exists, cleanup; Reprocess data\n", + " # colab: download and process data from scratch in Google Colab\n", + "\n", + "try:\n", + " import google.colab\n", + " notebook = 'colab'\n", + "except:\n", + " pass" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "

\n", + "

\n", + "\n", + "

\n", + "Analysis code for the paper \"RF shimming in the cervical spinal cord at 7T\"\n", + "

\n", + "\n", + "\n", + "

\n", + "Daniel Papp1, Kyle M. Gilbert2,3, Gaspard Cereza1, Alexandre D’Astous1, Mathieu Boudreau1, Marcus Couch4, Pedram Yazdanbakhsh5, Robert L. Barry6,7,8, Eva Alonso Ortiz1, Julien Cohen-Adad1,9,10\n", + "

\n", + "\n", + "\n", + "
\n", + "\n", + "

\n", + "1NeuroPoly Lab, Institute of Biomedical Engineering, Polytechnique Montreal, Montreal, QC, Canada,\n", + "2 Centre for Functional and Metabolic Mapping, The University of Western Ontario, London, ON, Canada,\n", + "3 Department of Medical Biophysics, The University of Western Ontario, London, ON, Canada,\n", + "4 Siemens Healthcare Limited, Montreal, QC, Canada,\n", + "5 McConnell Brain Imaging Centre, Montreal Neurological Institute, McGill University, Montreal, QC, Canada,\n", + "6 Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, MA, USA,\n", + "7 Harvard Medical School, Boston, MA, USA,\n", + "8 Harvard-Massachusetts Institute of Technology Health Sciences & Technology, Cambridge, MA, USA,\n", + "9 Mila - Quebec AI Institute, Montreal, QC, Canada,\n", + "10 Functional Neuroimaging Unit, Centre de recherche de l'Institut universitaire de gériatrie de Montréal QC, Canada\n", + "\n", + "

\n", + "

\n", + "

\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "

Summary

\n", + "\n", + "Data was collected from five participants between two 7T sites with a custom 8Tx/20Rx parallel transmission (pTx) coil. We explored two RF shimming approaches from an MRI vendor and four from an open-source toolbox, showcasing their ability to enhance transmit field and signal homogeneity along the cervical spinal cord.\n", + "\n", + "The results indicate significant improvements in B1+ efficiency and cerebrospinal fluid / spinal cord signal ratio across various RF shimming conditions compared to the vendor based methods.\n", + "\n", + "The study's findings highlight the potential of RF shimming to advance 7T MRI's clinical utility for central nervous system imaging by enabling more homogenous and efficient spinal cord imaging. Additionally, the research incorporates a reproducible Jupyter Notebook, enhancing the study's transparency and facilitating peer verification. By also making the data openly accessible on OpenNeuro, we ensure that the scientific community can further explore, validate, and build upon our findings, contributing to the collective advancement in high-resolution imaging techniques.\n", + "\n", + "```{figure} ../featured.png\n", + "---\n", + "width: 900px\n", + "name: fig1\n", + "align: center\n", + "---\n", + "```" + ] + }, { "cell_type": "markdown", - "id": "459f60d5", "metadata": {}, "source": [ - "# Analysis code for the paper \"RF shimming in the cervical spinal cord at 7T\"\n", + "

\n", + "Figure 1. Overview of the RF shimming procedure. The top panel shows the RF coil used for the experiments, alongside the Tx coil geometry and the electromagnetic simulation results (on Gustav model) yielding the CP mode used for this coil. The bottom panel shows the RF shimming procedure (with approximate duration). First, GRE and tfl_rfmap scans are acquired (4min30s). Second, these images are transferred via ethernet socket from the MRI console onto a separate laptop running Shimming Toolbox and SCT (əs). Third, the spinal cord is automatically segmented to produce a mask that is resampled into the space of the individual coil magnitude and phase images of the tfl_rfmap scan (~5s). Fourth, the RF shim weights are calculated according to the defined constraints for each shim scenario (1min total).\n", + "

" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "

Statement of Need

\n", "\n", - "## Data\n", + "Advancing the development of 7T MRI for spinal cord imaging is crucial for the enhanced diagnosis and monitoring of various neurodegenerative diseases {cite:p}`Kearney2015-py` and traumas {cite:p}`David2019-jy`. However, a significant challenge at this field strength is the transmit field inhomogeneity {cite:p}`Ibrahim2001-xt,Collins2005-za,Roschmann1987-om,Yang2002-ui`. Such inhomogeneity is particularly problematic for imaging the small, deep anatomical structures of the cervical spinal cord , as it can cause uneven signal intensity and elevate the local specific absorption ratio, compromising image quality. This multi-site study explores several radiofrequency (RF) shimming techniques in the cervical spinal cord at 7T.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 1     |     Data\n", "\n", "The data can be downloaded from: https://openneuro.org/datasets/ds004906\n", "\n", "The structure of the input dataset is as follows (JSON sidecars are not listed for clarity):\n", + "\n", + "```{toggle}\n", "~~~\n", "ds004906\n", "├── CHANGES\n", @@ -50,9 +159,18 @@ "├── sub-04\n", "└── sub-05\n", "~~~\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2     |     Overview of processing pipeline\n", "\n", + "During the data acquisition stage, RF shimming was done using the [Shimming Toolbox](https://github.com/shimming-toolbox/shimming-toolbox) {cite:p}`DAstous2023` during the acquisition stage.\n", "\n", - "## Overview of processing pipeline\n", + "The post-processing pipeline uses the [Spinal Cord Toolbox](https://spinalcordtoolbox.com) {cite:p}`DeLeener201724`.\n", "\n", "For each subject:\n", "\n", @@ -67,18 +185,77 @@ " - Convert the B1 map to nT/V units\n", " - Extract the B1 map value within the SC\n", "\n", - ">Slow processes are indicated with the emoji ⏳" + ">Slow processes are indicated with the emoji ⏳\n", + "\n" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "4f78c5ea", + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 3     |     Requirements\n", + "\n", + "* Install [Spinal Cord Toolbox](https://spinalcordtoolbox.com/user_section/installation.html) {cite:p}`DeLeener201724`, eg\n", + "\n", + "```shell\n", + "# Install Python libaries\n", + "!wget -O requirements.txt https://raw.githubusercontent.com/shimming-toolbox/rf-shimming-7t/main/requirements.txt\n", + "!pip install -r requirements.txt\n", + "\n", + "# Install SCT ⏳\n", + "!git clone --depth 1 https://github.com/spinalcordtoolbox/spinalcordtoolbox.git\n", + "!yes | spinalcordtoolbox/install_sct\n", + "os.environ['PATH'] += f\":/content/spinalcordtoolbox/bin\"\n", + "```\n", + "\n", + "* Download the project repository\n", + "\n", + "```shell\n", + "git clone https://github.com/shimming-toolbox/rf-shimming-7t-neurolibre.git\n", + "cd rf-shimming-7t-neurolibre\n", + "```\n", + "\n", + "* Install requirements\n", + "\n", + "```shell\n", + "pip install -e binder/requirements.txt\n", + "```\n", + "\n", + "* Download data\n", + "\n", + "```shell\n", + "cd content\n", + "repo2data -r ../binder/data_requirement.json\n", + "```\n" + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# Necessary imports\n", + "# 4     |     Environment setup\n", "\n", + "In a Python shell:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Import the necessary modules." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [ + "remove_output", + "hide_input" + ] + }, + "outputs": [], + "source": [ "from datetime import datetime, timedelta\n", "import os\n", "import re\n", @@ -94,6 +271,8 @@ "from scipy.ndimage import uniform_filter1d\n", "from scipy.stats import f_oneway, ttest_rel\n", "import shutil\n", + "from shutil import rmtree\n", + "from pathlib import Path\n", "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n", "from statsmodels.stats.anova import AnovaRM\n", "import statsmodels.api as sm" @@ -101,10 +280,22 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "bc5a2ec0", - "metadata": {}, - "outputs": [], + "execution_count": 3, + "metadata": { + "tags": [ + "remove_input", + "remove_output" + ] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/Users/mathieuboudreau/neuropoly/github/rf-shimming-7t-neurolibre/content\n" + ] + } + ], "source": [ "# Start timer\n", "start_time = datetime.now()\n", @@ -115,60 +306,111 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "5efab0b8", - "metadata": {}, + "execution_count": 4, + "metadata": { + "tags": [ + "remove_input", + "remove_output" + ] + }, "outputs": [], "source": [ + "# Google Colab-Only\n", "# ⚠️ No need to run this cell if you run this notebook locally and already have these dependencies installed.\n", "\n", - "# Install packages on system\n", - "!sudo apt-get update\n", - "!sudo apt-get install git-annex\n", + "if notebook == 'colab':\n", + " # Install packages on system\n", + " !sudo apt-get update\n", + " !sudo apt-get install git-annex\n", "\n", - "# Install Python libaries\n", - "!wget -O requirements.txt https://raw.githubusercontent.com/shimming-toolbox/rf-shimming-7t/main/requirements.txt\n", - "!pip install -r requirements.txt\n", + " # Install Python libaries\n", + " !wget -O requirements.txt https://raw.githubusercontent.com/shimming-toolbox/rf-shimming-7t/main/requirements.txt\n", + " !pip install -r requirements.txt\n", "\n", - "# Install SCT ⏳\n", - "!git clone --depth 1 https://github.com/spinalcordtoolbox/spinalcordtoolbox.git\n", - "!yes | spinalcordtoolbox/install_sct\n", - "os.environ['PATH'] += f\":/content/spinalcordtoolbox/bin\"" + " # Install SCT ⏳\n", + " !git clone --depth 1 https://github.com/spinalcordtoolbox/spinalcordtoolbox.git\n", + " !yes | spinalcordtoolbox/install_sct\n", + " os.environ['PATH'] += f\":/content/spinalcordtoolbox/bin\"" ] }, { "cell_type": "code", - "execution_count": null, - "id": "3393a694", - "metadata": {}, + "execution_count": 5, + "metadata": { + "tags": [ + "remove_input", + "remove_output" + ] + }, "outputs": [], "source": [ "# Download data and define path variables\n", "\n", - "!datalad install https://github.com/OpenNeuroDatasets/ds004906.git\n", - "os.chdir(\"ds004906\")\n", - "!datalad get . # uncomment for production\n", - "# !datalad get sub-01/ # debugging\n", - "# Get derivatives containing manual labels\n", - "!datalad get derivatives" + "# Google Colab-Only\n", + "if notebook=='colab':\n", + " !datalad install https://github.com/OpenNeuroDatasets/ds004906.git\n", + " os.chdir(\"ds004906\")\n", + " !datalad get . # uncomment for production\n", + " # !datalad get sub-01/ # debugging\n", + " # Get derivatives containing manual labels\n", + " !datalad get derivatives\n", + "\n", + "if notebook!='colab':\n", + " os.chdir(\"../data/rf-shimming-7t/ds004906\")\n", + " if notebook=='neurolibre-clean':\n", + " flist = open('cleanup.txt', 'r')\n", + " for f in flist:\n", + " fname = f.rstrip('\\n')\n", + " fname = Path(os.getcwd()) / Path(fname)\n", + "\n", + " # or, if you get rid of os.chdir(path) above,\n", + " # fname = os.path.join(path, f.rstrip())\n", + " if os.path.isfile(fname): # this makes the code more robust\n", + " print('Removing file: ' + str(fname))\n", + " os.remove(fname)\n", + " elif os.path.isdir(fname):\n", + " print('Removing directory: ' + str(fname))\n", + " rmtree(fname)\n", + "\n", + " # also, don't forget to close the text file:\n", + " flist.close()\n" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "6fe120e8", + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# Define useful variables\n", - "\n", + "Define variables." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "tags": [ + "hide_input", + "remove_output" + ] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "path_data: /Users/mathieuboudreau/neuropoly/github/rf-shimming-7t-neurolibre/data/rf-shimming-7t/ds004906\n", + "shim_modes: ['CP', 'patient', 'volume', 'phase', 'CoV', 'target', 'SAReff']\n", + "subjects: ['sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05']\n" + ] + } + ], + "source": [ "path_data = os.getcwd()\n", "print(f\"path_data: {path_data}\")\n", "path_labels = os.path.join(path_data, \"derivatives\", \"labels\")\n", "path_qc = os.path.join(path_data, \"qc\")\n", "shim_modes = [\"CP\", \"patient\", \"volume\", \"phase\", \"CoV\", \"target\", \"SAReff\"]\n", - "shim_modes_MPRAGE = [\"CP\", \"CoV\"] # TODO: change variable name PEP8\n", - "# shim_modes = [\"CP\", \"CoV\"] # debugging\n", + "shim_modes_MPRAGE = [\"CP\", \"CoV\"] \n", + "\n", "print(f\"shim_modes: {shim_modes}\")\n", "subjects = sorted(glob.glob(\"sub-*\"))\n", "print(f\"subjects: {subjects}\")\n", @@ -180,155 +422,236 @@ }, { "cell_type": "markdown", - "id": "b0c4e225", "metadata": {}, "source": [ - "## Process anat/T2starw (GRE)" + "# 5     |     Process anat/T2starw (GRE)\n" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "e95d87a4-3194-4369-b6cb-933ad1d5cfa3", + "cell_type": "markdown", "metadata": {}, + "source": [ + "Run segmentation on GRE scan\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Run segmentation on GRE scan\n", "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " # The \"CoV\" RF shimming scenario was chosen as the segmentation baseline due to the more homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance in the C7-T2 region\n", - " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-SC_seg.nii.gz\")\n", - " if os.path.exists(fname_manual_seg):\n", - " # Manual segmentation already exists. Copy it to local folder\n", - " print(f\"{subject}: Manual segmentation found\\n\")\n", - " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_seg.nii.gz\")\n", - " # Generate QC report to make sure the manual segmentation is correct\n", - " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", - " else:\n", - " # Manual segmentation does not exist. Run automatic segmentation.\n", - " print(f\"{subject}: Manual segmentation not found\")\n", - " !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " # The \"CoV\" RF shimming scenario was chosen as the segmentation baseline due to the more homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance in the C7-T2 region\n", + " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-SC_seg.nii.gz\")\n", + " if os.path.exists(fname_manual_seg):\n", + " # Manual segmentation already exists. Copy it to local folder\n", + " print(f\"{subject}: Manual segmentation found\\n\")\n", + " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_seg.nii.gz\")\n", + " # Generate QC report to make sure the manual segmentation is correct\n", + " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", + " else:\n", + " # Manual segmentation does not exist. Run automatic segmentation.\n", + " print(f\"{subject}: Manual segmentation not found\")\n", + " !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "eb9a9267", + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# Copy CSF masks from the derivatives folder\n", + "Copy CSF masks from the derivatives folder.\n", "\n", - "# For more details about how these masks were created, see: https://github.com/shimming-toolbox/rf-shimming-7t/issues/67\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", - " if os.path.exists(fname_manual_seg):\n", - " # Manual segmentation already exists. Copy it to local folder\n", - " print(f\"{subject}: Manual segmentation found\\n\")\n", - " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", - " # Generate QC report to make sure the manual segmentation is correct\n", - " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", - " else:\n", - " # Manual segmentation does not exist. Panic!\n", - " print(f\"{subject}: Manual segmentation not found 😱\")" + "For more details about how these masks were created, see: [https://github.com/shimming-toolbox/rf-shimming-7t/issues/67](https://github.com/shimming-toolbox/rf-shimming-7t/issues/67)." ] }, { "cell_type": "code", - "execution_count": null, - "id": "862dc562", - "metadata": {}, + "execution_count": 8, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Crop GRE scan for faster processing and better registration results\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz\n", - " !sct_crop_image -i {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", + " if os.path.exists(fname_manual_seg):\n", + " # Manual segmentation already exists. Copy it to local folder\n", + " print(f\"{subject}: Manual segmentation found\\n\")\n", + " shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n", + " # Generate QC report to make sure the manual segmentation is correct\n", + " !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n", + " else:\n", + " # Manual segmentation does not exist. Panic!\n", + " print(f\"{subject}: Manual segmentation not found 😱\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Crop GRE scan for faster processing and better registration results\n" ] }, { "cell_type": "code", - "execution_count": null, - "id": "87030540", + "execution_count": 9, + "metadata": { + "tags": [ + "hide_input" + ] + }, + "outputs": [], + "source": [ + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n", + " !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz\n", + " !sct_crop_image -i {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "Label vertebrae on GRE scan" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Label vertebrae on GRE scan\n", "\n", "# Given the low resolution of the GRE scan, the automatic detection of C2-C3 disc is unreliable. Therefore we need to use the manual disc labels that are part of the dataset.\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n", - " !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n", - " # Generate QC report to assess labeled segmentation\n", - " !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n", + " !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n", + " # Generate QC report to assess labeled segmentation\n", + " !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "c46e8d9a", + "cell_type": "markdown", "metadata": {}, + "source": [ + "Register *_T2starw to CoV_T2starw ⏳" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Register *_T2starw to CoV_T2starw ⏳\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " for shim_mode in shim_modes:\n", - " # Don't do it for CoV_T2starw\n", - " if shim_mode != 'CoV':\n", - " !sct_register_multimodal -i {subject}_acq-{shim_mode}_T2starw.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=dl -qc {path_qc}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " for shim_mode in shim_modes:\n", + " # Don't do it for CoV_T2starw\n", + " if shim_mode != 'CoV':\n", + " !sct_register_multimodal -i {subject}_acq-{shim_mode}_T2starw.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=dl -qc {path_qc}" ] }, { "cell_type": "markdown", - "id": "aac9d889", "metadata": {}, "source": [ - "### Verify QC report (GRE segmentation)\n", - "\n", - "Open the quality control (QC) report located under `ds004906/qc/index.html`. Make sure the spinal cord segmentations are correct before resuming the analysis.\n", - "\n", - ">If you run this notebook on Google Colab, you can skip as the QC report cannot easily be viewed from Google Colab. If you *really* want to see the QC report, you can create a cell where you zip the 'qc/' folder, then you can download it on your local station and open the index.html file. \n" + "Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed." ] }, { "cell_type": "code", - "execution_count": null, - "id": "90deb5e5", + "execution_count": 12, "metadata": { - "scrolled": true + "tags": [ + "hide_input" + ] }, "outputs": [], "source": [ - "# Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed\n", "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", - " for shim_mode in shim_modes:\n", - " # Shim methods are registered to the CoV T2starw scan, so we need to use the added suffix to identify them\n", - " if shim_mode == 'CoV':\n", - " file_suffix = 'crop'\n", - " else:\n", - " file_suffix = 'reg'\n", - " fname_result_sc = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_sc}\n", - " fname_result_csf = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_csf}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"anat\"))\n", + " for shim_mode in shim_modes:\n", + " # Shim methods are registered to the CoV T2starw scan, so we need to use the added suffix to identify them\n", + " if shim_mode == 'CoV':\n", + " file_suffix = 'crop'\n", + " else:\n", + " file_suffix = 'reg'\n", + " fname_result_sc = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n", + " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_sc}\n", + " fname_result_csf = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n", + " !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_csf}" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "a00909b8", + "cell_type": "markdown", "metadata": {}, - "outputs": [], + "source": [ + "Make figure of CSF/SC signal ratio from T2starw scan" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "tags": [ + "hide_input", + "remove_output" + ] + }, + "outputs": [ + { + "ename": "FileNotFoundError", + "evalue": "[Errno 2] No such file or directory: '/Users/mathieuboudreau/neuropoly/github/rf-shimming-7t-neurolibre/data/rf-shimming-7t/ds004906/derivatives/results/sub-01_acq-CP_T2starw_label-SC.csv'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/var/folders/l1/1xswjffd73l8yp7dd7pq9lyw0000gn/T/ipykernel_3003/2567902214.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 54\u001b[0m \u001b[0;31m# Get signal in SC\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 55\u001b[0m \u001b[0mfile_csv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_results\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mf\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 56\u001b[0;31m \u001b[0mdf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile_csv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 57\u001b[0m \u001b[0mdata_sc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'WA()'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 58\u001b[0m \u001b[0mdata_sc_smoothed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msmooth_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata_sc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/pandas/io/parsers/readers.py\u001b[0m in \u001b[0;36mread_csv\u001b[0;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)\u001b[0m\n\u001b[1;32m 910\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwds_defaults\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 911\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 912\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0m_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 913\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 914\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/pandas/io/parsers/readers.py\u001b[0m in \u001b[0;36m_read\u001b[0;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[1;32m 575\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 576\u001b[0m \u001b[0;31m# Create the parser.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 577\u001b[0;31m \u001b[0mparser\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTextFileReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 578\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 579\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mchunksize\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0miterator\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/pandas/io/parsers/readers.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[1;32m 1405\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1406\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhandles\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mIOHandles\u001b[0m \u001b[0;34m|\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1407\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_make_engine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mengine\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1408\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1409\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/pandas/io/parsers/readers.py\u001b[0m in \u001b[0;36m_make_engine\u001b[0;34m(self, f, engine)\u001b[0m\n\u001b[1;32m 1659\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;34m\"b\"\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1660\u001b[0m \u001b[0mmode\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m\"b\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1661\u001b[0;31m self.handles = get_handle(\n\u001b[0m\u001b[1;32m 1662\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1663\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/pandas/io/common.py\u001b[0m in \u001b[0;36mget_handle\u001b[0;34m(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)\u001b[0m\n\u001b[1;32m 857\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mioargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mencoding\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m\"b\"\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mioargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 858\u001b[0m \u001b[0;31m# Encoding\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 859\u001b[0;31m handle = open(\n\u001b[0m\u001b[1;32m 860\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 861\u001b[0m \u001b[0mioargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/Users/mathieuboudreau/neuropoly/github/rf-shimming-7t-neurolibre/data/rf-shimming-7t/ds004906/derivatives/results/sub-01_acq-CP_T2starw_label-SC.csv'" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0UAAAk3CAYAAABWHaQxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACdn0lEQVR4nOzdf2zV9b348Vdb7ClmtuK4tMDquOp1blNBQXqrM8abziYj7PLHvevQACEyr5MZtdmd4A8650a5u2pIrnVE5q77xwubGWYZpM51ErNrb8iAJpoLGIYMYtYC10vLra6F9vP9Y7H7dhTllP6wvB+P5PzR997v83mf5Q369HN6TkGWZVkAAAAkqnC8NwAAADCeRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQtLyj6LXXXouFCxfGjBkzoqCgIF566aWPXLN9+/a4/vrrI5fLxRVXXBHPP//8MLYKAAAw8vKOou7u7pg9e3Y0NTWd1fy33347FixYELfeemu0tbXF/fffHytWrIiXX345780CAACMtIIsy7JhLy4oiC1btsSiRYvOOOfBBx+MrVu3xptvvjkw9tWvfjWOHz8ezc3Nw700AADAiJg02hdobW2NmpqaQWO1tbVx//33n3FNT09P9PT0DPzc398f7777bnzyk5+MgoKC0doqAADwMZdlWZw4cSJmzJgRhYUj8xEJox5F7e3tUV5ePmisvLw8urq64v3334/JkyeftqaxsTEee+yx0d4aAAAwQR0+fDg+9alPjchzjXoUDcfq1aujvr5+4OfOzs649NJL4/Dhw1FaWjqOOwMAAMZTV1dXVFZWxkUXXTRizznqUVRRUREdHR2Dxjo6OqK0tHTIu0QREblcLnK53GnjpaWloggAABjRX6sZ9e8pqq6ujpaWlkFjr7zySlRXV4/2pQEAAD5S3lH0f//3f9HW1hZtbW0R8aeP3G5ra4tDhw5FxJ/e+rZ06dKB+XfffXccOHAgvvWtb8XevXvjmWeeiZ/85CfxwAMPjMwrAAAAOAd5R9Fvf/vbuO666+K6666LiIj6+vq47rrrYs2aNRER8Yc//GEgkCIi/vqv/zq2bt0ar7zySsyePTuefPLJ+OEPfxi1tbUj9BIAAACG75y+p2isdHV1RVlZWXR2dvqdIgAASNhotMGo/04RAADAx5koAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEjasKKoqakpZs2aFSUlJVFVVRU7duz40Pnr16+Pz3zmMzF58uSorKyMBx54IP74xz8Oa8MAAAAjKe8o2rx5c9TX10dDQ0Ps2rUrZs+eHbW1tXHkyJEh57/wwguxatWqaGhoiD179sRzzz0XmzdvjoceeuicNw8AAHCu8o6ip556Kr72ta/F8uXL43Of+1xs2LAhLrzwwvjRj3405PzXX389brrpprj99ttj1qxZcdttt8XixYs/8u4SAADAWMgrinp7e2Pnzp1RU1Pz5ycoLIyamppobW0dcs2NN94YO3fuHIigAwcOxLZt2+JLX/rSGa/T09MTXV1dgx4AAACjYVI+k48dOxZ9fX1RXl4+aLy8vDz27t075Jrbb789jh07Fl/4whciy7I4depU3H333R/69rnGxsZ47LHH8tkaAADAsIz6p89t37491q5dG88880zs2rUrfvazn8XWrVvj8ccfP+Oa1atXR2dn58Dj8OHDo71NAAAgUXndKZo6dWoUFRVFR0fHoPGOjo6oqKgYcs2jjz4aS5YsiRUrVkRExDXXXBPd3d1x1113xcMPPxyFhad3WS6Xi1wul8/WAAAAhiWvO0XFxcUxd+7caGlpGRjr7++PlpaWqK6uHnLNe++9d1r4FBUVRURElmX57hcAAGBE5XWnKCKivr4+li1bFvPmzYv58+fH+vXro7u7O5YvXx4REUuXLo2ZM2dGY2NjREQsXLgwnnrqqbjuuuuiqqoq9u/fH48++mgsXLhwII4AAADGS95RVFdXF0ePHo01a9ZEe3t7zJkzJ5qbmwc+fOHQoUOD7gw98sgjUVBQEI888ki888478Vd/9VexcOHC+N73vjdyrwIAAGCYCrIJ8B62rq6uKCsri87OzigtLR3v7QAAAONkNNpg1D99DgAA4ONMFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQtGFFUVNTU8yaNStKSkqiqqoqduzY8aHzjx8/HitXrozp06dHLpeLK6+8MrZt2zasDQMAAIykSfku2Lx5c9TX18eGDRuiqqoq1q9fH7W1tbFv376YNm3aafN7e3vji1/8YkybNi1efPHFmDlzZvz+97+Piy++eCT2DwAAcE4KsizL8llQVVUVN9xwQzz99NMREdHf3x+VlZVx7733xqpVq06bv2HDhvjXf/3X2Lt3b1xwwQXD2mRXV1eUlZVFZ2dnlJaWDus5AACAiW802iCvt8/19vbGzp07o6am5s9PUFgYNTU10draOuSan//851FdXR0rV66M8vLyuPrqq2Pt2rXR19d3bjsHAAAYAXm9fe7YsWPR19cX5eXlg8bLy8tj7969Q645cOBA/PrXv4477rgjtm3bFvv374977rknTp48GQ0NDUOu6enpiZ6enoGfu7q68tkmAADAWRv1T5/r7++PadOmxbPPPhtz586Nurq6ePjhh2PDhg1nXNPY2BhlZWUDj8rKytHeJgAAkKi8omjq1KlRVFQUHR0dg8Y7OjqioqJiyDXTp0+PK6+8MoqKigbGPvvZz0Z7e3v09vYOuWb16tXR2dk58Dh8+HA+2wQAADhreUVRcXFxzJ07N1paWgbG+vv7o6WlJaqrq4dcc9NNN8X+/fujv79/YOytt96K6dOnR3Fx8ZBrcrlclJaWDnoAAACMhrzfPldfXx8bN26MH//4x7Fnz574+te/Ht3d3bF8+fKIiFi6dGmsXr16YP7Xv/71ePfdd+O+++6Lt956K7Zu3Rpr166NlStXjtyrAAAAGKa8v6eorq4ujh49GmvWrIn29vaYM2dONDc3D3z4wqFDh6Kw8M+tVVlZGS+//HI88MADce2118bMmTPjvvvuiwcffHDkXgUAAMAw5f09RePB9xQBAAARH4PvKQIAADjfiCIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApA0ripqammLWrFlRUlISVVVVsWPHjrNat2nTpigoKIhFixYN57IAAAAjLu8o2rx5c9TX10dDQ0Ps2rUrZs+eHbW1tXHkyJEPXXfw4MH45je/GTfffPOwNwsAADDS8o6ip556Kr72ta/F8uXL43Of+1xs2LAhLrzwwvjRj350xjV9fX1xxx13xGOPPRaXXXbZOW0YAABgJOUVRb29vbFz586oqan58xMUFkZNTU20traecd13vvOdmDZtWtx5551ndZ2enp7o6uoa9AAAABgNeUXRsWPHoq+vL8rLyweNl5eXR3t7+5BrfvOb38Rzzz0XGzduPOvrNDY2RllZ2cCjsrIyn20CAACctVH99LkTJ07EkiVLYuPGjTF16tSzXrd69ero7OwceBw+fHgUdwkAAKRsUj6Tp06dGkVFRdHR0TFovKOjIyoqKk6b/7vf/S4OHjwYCxcuHBjr7+//04UnTYp9+/bF5Zdfftq6XC4XuVwun60BAAAMS153ioqLi2Pu3LnR0tIyMNbf3x8tLS1RXV192vyrrroq3njjjWhraxt4fPnLX45bb7012travC0OAAAYd3ndKYqIqK+vj2XLlsW8efNi/vz5sX79+uju7o7ly5dHRMTSpUtj5syZ0djYGCUlJXH11VcPWn/xxRdHRJw2DgAAMB7yjqK6uro4evRorFmzJtrb22POnDnR3Nw88OELhw4disLCUf1VJQAAgBFTkGVZNt6b+ChdXV1RVlYWnZ2dUVpaOt7bAQAAxslotIFbOgAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDShhVFTU1NMWvWrCgpKYmqqqrYsWPHGedu3Lgxbr755pgyZUpMmTIlampqPnQ+AADAWMo7ijZv3hz19fXR0NAQu3btitmzZ0dtbW0cOXJkyPnbt2+PxYsXx6uvvhqtra1RWVkZt912W7zzzjvnvHkAAIBzVZBlWZbPgqqqqrjhhhvi6aefjoiI/v7+qKysjHvvvTdWrVr1kev7+vpiypQp8fTTT8fSpUvP6ppdXV1RVlYWnZ2dUVpams92AQCA88hotEFed4p6e3tj586dUVNT8+cnKCyMmpqaaG1tPavneO+99+LkyZNxySWXnHFOT09PdHV1DXoAAACMhryi6NixY9HX1xfl5eWDxsvLy6O9vf2snuPBBx+MGTNmDAqrv9TY2BhlZWUDj8rKyny2CQAAcNbG9NPn1q1bF5s2bYotW7ZESUnJGeetXr06Ojs7Bx6HDx8ew10CAAApmZTP5KlTp0ZRUVF0dHQMGu/o6IiKiooPXfvEE0/EunXr4le/+lVce+21Hzo3l8tFLpfLZ2sAAADDktedouLi4pg7d260tLQMjPX390dLS0tUV1efcd33v//9ePzxx6O5uTnmzZs3/N0CAACMsLzuFEVE1NfXx7Jly2LevHkxf/78WL9+fXR3d8fy5csjImLp0qUxc+bMaGxsjIiIf/mXf4k1a9bECy+8ELNmzRr43aNPfOIT8YlPfGIEXwoAAED+8o6iurq6OHr0aKxZsyba29tjzpw50dzcPPDhC4cOHYrCwj/fgPrBD34Qvb298Q//8A+DnqehoSG+/e1vn9vuAQAAzlHe31M0HnxPEQAAEPEx+J4iAACA840oAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpw4qipqammDVrVpSUlERVVVXs2LHjQ+f/9Kc/jauuuipKSkrimmuuiW3btg1rswAAACMt7yjavHlz1NfXR0NDQ+zatStmz54dtbW1ceTIkSHnv/7667F48eK48847Y/fu3bFo0aJYtGhRvPnmm+e8eQAAgHNVkGVZls+CqqqquOGGG+Lpp5+OiIj+/v6orKyMe++9N1atWnXa/Lq6uuju7o5f/OIXA2N/+7d/G3PmzIkNGzac1TW7urqirKwsOjs7o7S0NJ/tAgAA55HRaIO87hT19vbGzp07o6am5s9PUFgYNTU10draOuSa1tbWQfMjImpra884HwAAYCxNymfysWPHoq+vL8rLyweNl5eXx969e4dc097ePuT89vb2M16np6cnenp6Bn7u7OyMiD9VIQAAkK4PmiDPN7x9qLyiaKw0NjbGY489dtp4ZWXlOOwGAAD4uPmf//mfKCsrG5HnyiuKpk6dGkVFRdHR0TFovKOjIyoqKoZcU1FRkdf8iIjVq1dHfX39wM/Hjx+PT3/603Ho0KERe+EwlK6urqisrIzDhw/7/TVGlbPGWHHWGCvOGmOls7MzLr300rjkkktG7DnziqLi4uKYO3dutLS0xKJFiyLiTx+00NLSEt/4xjeGXFNdXR0tLS1x//33D4y98sorUV1dfcbr5HK5yOVyp42XlZX5Q8aYKC0tddYYE84aY8VZY6w4a4yVwsKR+8rVvN8+V19fH8uWLYt58+bF/PnzY/369dHd3R3Lly+PiIilS5fGzJkzo7GxMSIi7rvvvrjlllviySefjAULFsSmTZvit7/9bTz77LMj9iIAAACGK+8oqquri6NHj8aaNWuivb095syZE83NzQMfpnDo0KFB1XbjjTfGCy+8EI888kg89NBD8Td/8zfx0ksvxdVXXz1yrwIAAGCYhvVBC9/4xjfO+Ha57du3nzb2j//4j/GP//iPw7lURPzp7XQNDQ1DvqUORpKzxlhx1hgrzhpjxVljrIzGWcv7y1sBAADOJyP320kAAAATkCgCAACSJooAAICkfWyiqKmpKWbNmhUlJSVRVVUVO3bs+ND5P/3pT+Oqq66KkpKSuOaaa2Lbtm1jtFMmunzO2saNG+Pmm2+OKVOmxJQpU6KmpuYjzyZ8IN+/1z6wadOmKCgoGPg+OPgo+Z6148ePx8qVK2P69OmRy+Xiyiuv9M9Rzkq+Z239+vXxmc98JiZPnhyVlZXxwAMPxB//+Mcx2i0T0WuvvRYLFy6MGTNmREFBQbz00ksfuWb79u1x/fXXRy6XiyuuuCKef/75vK/7sYiizZs3R319fTQ0NMSuXbti9uzZUVtbG0eOHBly/uuvvx6LFy+OO++8M3bv3h2LFi2KRYsWxZtvvjnGO2eiyfesbd++PRYvXhyvvvpqtLa2RmVlZdx2223xzjvvjPHOmWjyPWsfOHjwYHzzm9+Mm2++eYx2ykSX71nr7e2NL37xi3Hw4MF48cUXY9++fbFx48aYOXPmGO+ciSbfs/bCCy/EqlWroqGhIfbs2RPPPfdcbN68OR566KEx3jkTSXd3d8yePTuamprOav7bb78dCxYsiFtvvTXa2tri/vvvjxUrVsTLL7+c34Wzj4H58+dnK1euHPi5r68vmzFjRtbY2Djk/K985SvZggULBo1VVVVl//RP/zSq+2Tiy/es/aVTp05lF110UfbjH/94tLbIeWI4Z+3UqVPZjTfemP3whz/Mli1blv393//9GOyUiS7fs/aDH/wgu+yyy7Le3t6x2iLniXzP2sqVK7O/+7u/GzRWX1+f3XTTTaO6T84fEZFt2bLlQ+d861vfyj7/+c8PGqurq8tqa2vzuta43ynq7e2NnTt3Rk1NzcBYYWFh1NTURGtr65BrWltbB82PiKitrT3jfIgY3ln7S++9916cPHkyLrnkktHaJueB4Z6173znOzFt2rS48847x2KbnAeGc9Z+/vOfR3V1daxcuTLKy8vj6quvjrVr10ZfX99YbZsJaDhn7cYbb4ydO3cOvMXuwIEDsW3btvjSl740JnsmDSPVBcP68taRdOzYsejr64vy8vJB4+Xl5bF3794h17S3tw85v729fdT2ycQ3nLP2lx588MGYMWPGaX/44P83nLP2m9/8Jp577rloa2sbgx1yvhjOWTtw4ED8+te/jjvuuCO2bdsW+/fvj3vuuSdOnjwZDQ0NY7FtJqDhnLXbb789jh07Fl/4whciy7I4depU3H333d4+x4g6Uxd0dXXF+++/H5MnTz6r5xn3O0UwUaxbty42bdoUW7ZsiZKSkvHeDueREydOxJIlS2Ljxo0xderU8d4O57n+/v6YNm1aPPvsszF37tyoq6uLhx9+ODZs2DDeW+M8s3379li7dm0888wzsWvXrvjZz34WW7dujccff3y8twanGfc7RVOnTo2ioqLo6OgYNN7R0REVFRVDrqmoqMhrPkQM76x94Iknnoh169bFr371q7j22mtHc5ucB/I9a7/73e/i4MGDsXDhwoGx/v7+iIiYNGlS7Nu3Ly6//PLR3TQT0nD+Xps+fXpccMEFUVRUNDD22c9+Ntrb26O3tzeKi4tHdc9MTMM5a48++mgsWbIkVqxYERER11xzTXR3d8ddd90VDz/8cBQW+m/znLszdUFpaelZ3yWK+BjcKSouLo65c+dGS0vLwFh/f3+0tLREdXX1kGuqq6sHzY+IeOWVV844HyKGd9YiIr7//e/H448/Hs3NzTFv3ryx2CoTXL5n7aqrroo33ngj2traBh5f/vKXBz5Jp7Kyciy3zwQynL/Xbrrppti/f/9AeEdEvPXWWzF9+nRBxBkN56y99957p4XPBzH+p9+hh3M3Yl2Q32dAjI5NmzZluVwue/7557P//u//zu66667s4osvztrb27Msy7IlS5Zkq1atGpj/n//5n9mkSZOyJ554ItuzZ0/W0NCQXXDBBdkbb7wxXi+BCSLfs7Zu3bqsuLg4e/HFF7M//OEPA48TJ06M10tggsj3rP0lnz7H2cr3rB06dCi76KKLsm984xvZvn37sl/84hfZtGnTsu9+97vj9RKYIPI9aw0NDdlFF12U/cd//Ed24MCB7Je//GV2+eWXZ1/5ylfG6yUwAZw4cSLbvXt3tnv37iwisqeeeirbvXt39vvf/z7LsixbtWpVtmTJkoH5Bw4cyC688MLsn//5n7M9e/ZkTU1NWVFRUdbc3JzXdT8WUZRlWfZv//Zv2aWXXpoVFxdn8+fPz/7rv/5r4H+75ZZbsmXLlg2a/5Of/CS78sors+Li4uzzn/98tnXr1jHeMRNVPmft05/+dBYRpz0aGhrGfuNMOPn+vfb/E0XkI9+z9vrrr2dVVVVZLpfLLrvssux73/tedurUqTHeNRNRPmft5MmT2be//e3s8ssvz0pKSrLKysrsnnvuyf73f/937DfOhPHqq68O+e9eH5ytZcuWZbfccstpa+bMmZMVFxdnl112Wfbv//7veV+3IMvcvwQAANI17r9TBAAAMJ5EEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQtLyj6LXXXouFCxfGjBkzoqCgIF566aWPXLN9+/a4/vrrI5fLxRVXXBHPP//8MLYKAAAw8vKOou7u7pg9e3Y0NTWd1fy33347FixYELfeemu0tbXF/fffHytWrIiXX345780CAACMtIIsy7JhLy4oiC1btsSiRYvOOOfBBx+MrVu3xptvvjkw9tWvfjWOHz8ezc3Nw700AADAiJg02hdobW2NmpqaQWO1tbVx//33n3FNT09P9PT0DPzc398f7777bnzyk5+MgoKC0doqAADwMZdlWZw4cSJmzJgRhYUj8xEJox5F7e3tUV5ePmisvLw8urq64v3334/JkyeftqaxsTEee+yx0d4aAAAwQR0+fDg+9alPjchzjXoUDcfq1aujvr5+4OfOzs649NJL4/Dhw1FaWjqOOwMAAMZTV1dXVFZWxkUXXTRizznqUVRRUREdHR2Dxjo6OqK0tHTIu0QREblcLnK53GnjpaWloggAABjRX6sZ9e8pqq6ujpaWlkFjr7zySlRXV4/2pQEAAD5S3lH0f//3f9HW1hZtbW0R8aeP3G5ra4tDhw5FxJ/e+rZ06dKB+XfffXccOHAgvvWtb8XevXvjmWeeiZ/85CfxwAMPjMwrAAAAOAd5R9Fvf/vbuO666+K6666LiIj6+vq47rrrYs2aNRER8Yc//GEgkCIi/vqv/zq2bt0ar7zySsyePTuefPLJ+OEPfxi1tbUj9BIAAACG75y+p2isdHV1RVlZWXR2dvqdIgAASNhotMGo/04RAADAx5koAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEjasKKoqakpZs2aFSUlJVFVVRU7duz40Pnr16+Pz3zmMzF58uSorKyMBx54IP74xz8Oa8MAAAAjKe8o2rx5c9TX10dDQ0Ps2rUrZs+eHbW1tXHkyJEh57/wwguxatWqaGhoiD179sRzzz0XmzdvjoceeuicNw8AAHCu8o6ip556Kr72ta/F8uXL43Of+1xs2LAhLrzwwvjRj3405PzXX389brrpprj99ttj1qxZcdttt8XixYs/8u4SAADAWMgrinp7e2Pnzp1RU1Pz5ycoLIyamppobW0dcs2NN94YO3fuHIigAwcOxLZt2+JLX/rSGa/T09MTXV1dgx4AAACjYVI+k48dOxZ9fX1RXl4+aLy8vDz27t075Jrbb789jh07Fl/4whciy7I4depU3H333R/69rnGxsZ47LHH8tkaAADAsIz6p89t37491q5dG88880zs2rUrfvazn8XWrVvj8ccfP+Oa1atXR2dn58Dj8OHDo71NAAAgUXndKZo6dWoUFRVFR0fHoPGOjo6oqKgYcs2jjz4aS5YsiRUrVkRExDXXXBPd3d1x1113xcMPPxyFhad3WS6Xi1wul8/WAAAAhiWvO0XFxcUxd+7caGlpGRjr7++PlpaWqK6uHnLNe++9d1r4FBUVRURElmX57hcAAGBE5XWnKCKivr4+li1bFvPmzYv58+fH+vXro7u7O5YvXx4REUuXLo2ZM2dGY2NjREQsXLgwnnrqqbjuuuuiqqoq9u/fH48++mgsXLhwII4AAADGS95RVFdXF0ePHo01a9ZEe3t7zJkzJ5qbmwc+fOHQoUOD7gw98sgjUVBQEI888ki888478Vd/9VexcOHC+N73vjdyrwIAAGCYCrIJ8B62rq6uKCsri87OzigtLR3v7QAAAONkNNpg1D99DgAA4ONMFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQtGFFUVNTU8yaNStKSkqiqqoqduzY8aHzjx8/HitXrozp06dHLpeLK6+8MrZt2zasDQMAAIykSfku2Lx5c9TX18eGDRuiqqoq1q9fH7W1tbFv376YNm3aafN7e3vji1/8YkybNi1efPHFmDlzZvz+97+Piy++eCT2DwAAcE4KsizL8llQVVUVN9xwQzz99NMREdHf3x+VlZVx7733xqpVq06bv2HDhvjXf/3X2Lt3b1xwwQXD2mRXV1eUlZVFZ2dnlJaWDus5AACAiW802iCvt8/19vbGzp07o6am5s9PUFgYNTU10draOuSan//851FdXR0rV66M8vLyuPrqq2Pt2rXR19d3bjsHAAAYAXm9fe7YsWPR19cX5eXlg8bLy8tj7969Q645cOBA/PrXv4477rgjtm3bFvv374977rknTp48GQ0NDUOu6enpiZ6enoGfu7q68tkmAADAWRv1T5/r7++PadOmxbPPPhtz586Nurq6ePjhh2PDhg1nXNPY2BhlZWUDj8rKytHeJgAAkKi8omjq1KlRVFQUHR0dg8Y7OjqioqJiyDXTp0+PK6+8MoqKigbGPvvZz0Z7e3v09vYOuWb16tXR2dk58Dh8+HA+2wQAADhreUVRcXFxzJ07N1paWgbG+vv7o6WlJaqrq4dcc9NNN8X+/fujv79/YOytt96K6dOnR3Fx8ZBrcrlclJaWDnoAAACMhrzfPldfXx8bN26MH//4x7Fnz574+te/Ht3d3bF8+fKIiFi6dGmsXr16YP7Xv/71ePfdd+O+++6Lt956K7Zu3Rpr166NlStXjtyrAAAAGKa8v6eorq4ujh49GmvWrIn29vaYM2dONDc3D3z4wqFDh6Kw8M+tVVlZGS+//HI88MADce2118bMmTPjvvvuiwcffHDkXgUAAMAw5f09RePB9xQBAAARH4PvKQIAADjfiCIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApA0ripqammLWrFlRUlISVVVVsWPHjrNat2nTpigoKIhFixYN57IAAAAjLu8o2rx5c9TX10dDQ0Ps2rUrZs+eHbW1tXHkyJEPXXfw4MH45je/GTfffPOwNwsAADDS8o6ip556Kr72ta/F8uXL43Of+1xs2LAhLrzwwvjRj350xjV9fX1xxx13xGOPPRaXXXbZOW0YAABgJOUVRb29vbFz586oqan58xMUFkZNTU20traecd13vvOdmDZtWtx5551ndZ2enp7o6uoa9AAAABgNeUXRsWPHoq+vL8rLyweNl5eXR3t7+5BrfvOb38Rzzz0XGzduPOvrNDY2RllZ2cCjsrIyn20CAACctVH99LkTJ07EkiVLYuPGjTF16tSzXrd69ero7OwceBw+fHgUdwkAAKRsUj6Tp06dGkVFRdHR0TFovKOjIyoqKk6b/7vf/S4OHjwYCxcuHBjr7+//04UnTYp9+/bF5Zdfftq6XC4XuVwun60BAAAMS153ioqLi2Pu3LnR0tIyMNbf3x8tLS1RXV192vyrrroq3njjjWhraxt4fPnLX45bb7012travC0OAAAYd3ndKYqIqK+vj2XLlsW8efNi/vz5sX79+uju7o7ly5dHRMTSpUtj5syZ0djYGCUlJXH11VcPWn/xxRdHRJw2DgAAMB7yjqK6uro4evRorFmzJtrb22POnDnR3Nw88OELhw4disLCUf1VJQAAgBFTkGVZNt6b+ChdXV1RVlYWnZ2dUVpaOt7bAQAAxslotIFbOgAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDShhVFTU1NMWvWrCgpKYmqqqrYsWPHGedu3Lgxbr755pgyZUpMmTIlampqPnQ+AADAWMo7ijZv3hz19fXR0NAQu3btitmzZ0dtbW0cOXJkyPnbt2+PxYsXx6uvvhqtra1RWVkZt912W7zzzjvnvHkAAIBzVZBlWZbPgqqqqrjhhhvi6aefjoiI/v7+qKysjHvvvTdWrVr1kev7+vpiypQp8fTTT8fSpUvP6ppdXV1RVlYWnZ2dUVpams92AQCA88hotEFed4p6e3tj586dUVNT8+cnKCyMmpqaaG1tPavneO+99+LkyZNxySWXnHFOT09PdHV1DXoAAACMhryi6NixY9HX1xfl5eWDxsvLy6O9vf2snuPBBx+MGTNmDAqrv9TY2BhlZWUDj8rKyny2CQAAcNbG9NPn1q1bF5s2bYotW7ZESUnJGeetXr06Ojs7Bx6HDx8ew10CAAApmZTP5KlTp0ZRUVF0dHQMGu/o6IiKiooPXfvEE0/EunXr4le/+lVce+21Hzo3l8tFLpfLZ2sAAADDktedouLi4pg7d260tLQMjPX390dLS0tUV1efcd33v//9ePzxx6O5uTnmzZs3/N0CAACMsLzuFEVE1NfXx7Jly2LevHkxf/78WL9+fXR3d8fy5csjImLp0qUxc+bMaGxsjIiIf/mXf4k1a9bECy+8ELNmzRr43aNPfOIT8YlPfGIEXwoAAED+8o6iurq6OHr0aKxZsyba29tjzpw50dzcPPDhC4cOHYrCwj/fgPrBD34Qvb298Q//8A+DnqehoSG+/e1vn9vuAQAAzlHe31M0HnxPEQAAEPEx+J4iAACA840oAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpw4qipqammDVrVpSUlERVVVXs2LHjQ+f/9Kc/jauuuipKSkrimmuuiW3btg1rswAAACMt7yjavHlz1NfXR0NDQ+zatStmz54dtbW1ceTIkSHnv/7667F48eK48847Y/fu3bFo0aJYtGhRvPnmm+e8eQAAgHNVkGVZls+CqqqquOGGG+Lpp5+OiIj+/v6orKyMe++9N1atWnXa/Lq6uuju7o5f/OIXA2N/+7d/G3PmzIkNGzac1TW7urqirKwsOjs7o7S0NJ/tAgAA55HRaINJ+Uzu7e2NnTt3xurVqwfGCgsLo6amJlpbW4dc09raGvX19YPGamtr46WXXjrjdXp6eqKnp2fg587Ozoj40/8BAABAuj5ogjzv7XyovKLo2LFj0dfXF+Xl5YPGy8vLY+/evUOuaW9vH3J+e3v7Ga/T2NgYjz322GnjlZWV+WwXAAA4T/3P//xPlJWVjchz5RVFY2X16tWD7i4dP348Pv3pT8ehQ4dG7IXDULq6uqKysjIOHz7srZqMKmeNseKsMVacNcZKZ2dnXHrppXHJJZeM2HPmFUVTp06NoqKi6OjoGDTe0dERFRUVQ66pqKjIa35ERC6Xi1wud9p4WVmZP2SMidLSUmeNMeGsMVacNcaKs8ZYKSwcuW8XyuuZiouLY+7cudHS0jIw1t/fHy0tLVFdXT3kmurq6kHzIyJeeeWVM84HAAAYS3m/fa6+vj6WLVsW8+bNi/nz58f69euju7s7li9fHhERS5cujZkzZ0ZjY2NERNx3331xyy23xJNPPhkLFiyITZs2xW9/+9t49tlnR/aVAAAADEPeUVRXVxdHjx6NNWvWRHt7e8yZMyeam5sHPkzh0KFDg25l3XjjjfHCCy/EI488Eg899FD8zd/8Tbz00ktx9dVXn/U1c7lcNDQ0DPmWOhhJzhpjxVljrDhrjBVnjbEyGmct7+8pAgAAOJ+M3G8nAQAATECiCAAASJooAgAAkiaKAACApH1soqipqSlmzZoVJSUlUVVVFTt27PjQ+T/96U/jqquuipKSkrjmmmti27ZtY7RTJrp8ztrGjRvj5ptvjilTpsSUKVOipqbmI88mfCDfv9c+sGnTpigoKIhFixaN7gY5b+R71o4fPx4rV66M6dOnRy6XiyuvvNI/Rzkr+Z619evXx2c+85mYPHlyVFZWxgMPPBB//OMfx2i3TESvvfZaLFy4MGbMmBEFBQXx0ksvfeSa7du3x/XXXx+5XC6uuOKKeP755/O+7sciijZv3hz19fXR0NAQu3btitmzZ0dtbW0cOXJkyPmvv/56LF68OO68887YvXt3LFq0KBYtWhRvvvnmGO+ciSbfs7Z9+/ZYvHhxvPrqq9Ha2hqVlZVx2223xTvvvDPGO2eiyfesfeDgwYPxzW9+M26++eYx2ikTXb5nrbe3N774xS/GwYMH48UXX4x9+/bFxo0bY+bMmWO8cyaafM/aCy+8EKtWrYqGhobYs2dPPPfcc7F58+Z46KGHxnjnTCTd3d0xe/bsaGpqOqv5b7/9dixYsCBuvfXWaGtri/vvvz9WrFgRL7/8cn4Xzj4G5s+fn61cuXLg576+vmzGjBlZY2PjkPO/8pWvZAsWLBg0VlVVlf3TP/3TqO6TiS/fs/aXTp06lV100UXZj3/849HaIueJ4Zy1U6dOZTfeeGP2wx/+MFu2bFn293//92OwUya6fM/aD37wg+yyyy7Lent7x2qLnCfyPWsrV67M/u7v/m7QWH19fXbTTTeN6j45f0REtmXLlg+d861vfSv7/Oc/P2isrq4uq62tzeta436nqLe3N3bu3Bk1NTUDY4WFhVFTUxOtra1DrmltbR00PyKitrb2jPMhYnhn7S+99957cfLkybjkkktGa5ucB4Z71r7zne/EtGnT4s477xyLbXIeGM5Z+/nPfx7V1dWxcuXKKC8vj6uvvjrWrl0bfX19Y7VtJqDhnLUbb7wxdu7cOfAWuwMHDsS2bdviS1/60pjsmTSMVBdMGslNDcexY8eir68vysvLB42Xl5fH3r17h1zT3t4+5Pz29vZR2ycT33DO2l968MEHY8aMGaf94YP/33DO2m9+85t47rnnoq2tbQx2yPliOGftwIED8etf/zruuOOO2LZtW+zfvz/uueeeOHnyZDQ0NIzFtpmAhnPWbr/99jh27Fh84QtfiCzL4tSpU3H33Xd7+xwj6kxd0NXVFe+//35Mnjz5rJ5n3O8UwUSxbt262LRpU2zZsiVKSkrGezucR06cOBFLliyJjRs3xtSpU8d7O5zn+vv7Y9q0afHss8/G3Llzo66uLh5++OHYsGHDeG+N88z27dtj7dq18cwzz8SuXbviZz/7WWzdujUef/zx8d4anGbc7xRNnTo1ioqKoqOjY9B4R0dHVFRUDLmmoqIir/kQMbyz9oEnnngi1q1bF7/61a/i2muvHc1tch7I96z97ne/i4MHD8bChQsHxvr7+yMiYtKkSbFv3764/PLLR3fTTEjD+Xtt+vTpccEFF0RRUdHA2Gc/+9lob2+P3t7eKC4uHtU9MzEN56w9+uijsWTJklixYkVERFxzzTXR3d0dd911Vzz88MNRWOi/zXPuztQFpaWlZ32XKOJjcKeouLg45s6dGy0tLQNj/f390dLSEtXV1UOuqa6uHjQ/IuKVV14543yIGN5Zi4j4/ve/H48//ng0NzfHvHnzxmKrTHD5nrWrrroq3njjjWhraxt4fPnLXx74JJ3Kysqx3D4TyHD+Xrvpppti//79A+EdEfHWW2/F9OnTBRFnNJyz9t57750WPh/E+J9+hx7O3Yh1QX6fATE6Nm3alOVyuez555/P/vu//zu76667sosvvjhrb2/PsizLlixZkq1atWpg/n/+539mkyZNyp544olsz549WUNDQ3bBBRdkb7zxxni9BCaIfM/aunXrsuLi4uzFF1/M/vCHPww8Tpw4MV4vgQki37P2l3z6HGcr37N26NCh7KKLLsq+8Y1vZPv27ct+8YtfZNOmTcu++93vjtdLYILI96w1NDRkF110UfYf//Ef2YEDB7Jf/vKX2eWXX5595StfGa+XwARw4sSJbPfu3dnu3buziMieeuqpbPfu3dnvf//7LMuybNWqVdmSJUsG5h84cCC78MILs3/+53/O9uzZkzU1NWVFRUVZc3NzXtf9WERRlmXZv/3bv2WXXnppVlxcnM2fPz/7r//6r4H/7ZZbbsmWLVs2aP5PfvKT7Morr8yKi4uzz3/+89nWrVvHeMdMVPmctU9/+tNZRJz2aGhoGPuNM+Hk+/fa/08UkY98z9rrr7+eVVVVZblcLrvsssuy733ve9mpU6fGeNdMRPmctZMnT2bf/va3s8svvzwrKSnJKisrs3vuuSf73//937HfOBPGq6++OuS/e31wtpYtW5bdcsstp62ZM2dOVlxcnF122WXZv//7v+d93YIsc/8SAABI17j/ThEAAMB4EkUAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNLyjqLXXnstFi5cGDNmzIiCgoJ46aWXPnLN9u3b4/rrr49cLhdXXHFFPP/888PYKgAAwMjLO4q6u7tj9uzZ0dTUdFbz33777ViwYEHceuut0dbWFvfff3+sWLEiXn755bw3CwAAMNIKsizLhr24oCC2bNkSixYtOuOcBx98MLZu3RpvvvnmwNhXv/rVOH78eDQ3Nw/30gAAACNi0mhfoLW1NWpqagaN1dbWxv3333/GNT09PdHT0zPwc39/f7z77rvxyU9+MgoKCkZrqwAAwMdclmVx4sSJmDFjRhQWjsxHJIx6FLW3t0d5efmgsfLy8ujq6or3338/Jk+efNqaxsbGeOyxx0Z7awAAwAR1+PDh+NSnPjUizzXqUTQcq1evjvr6+oGfOzs749JLL43Dhw9HaWnpOO4MAAAYT11dXVFZWRkXXXTRiD3nqEdRRUVFdHR0DBrr6OiI0tLSIe8SRUTkcrnI5XKnjZeWlooiAABgRH+tZtS/p6i6ujpaWloGjb3yyitRXV092pcGAAD4SHlH0f/93/9FW1tbtLW1RcSfPnK7ra0tDh06FBF/euvb0qVLB+bffffdceDAgfjWt74Ve/fujWeeeSZ+8pOfxAMPPDAyrwAAAOAc5B1Fv/3tb+O6666L6667LiIi6uvr47rrros1a9ZERMQf/vCHgUCKiPjrv/7r2Lp1a7zyyisxe/bsePLJJ+OHP/xh1NbWjtBLAAAAGL5z+p6isdLV1RVlZWXR2dnpd4oAACBho9EGo/47RQAAAB9noggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgacOKoqamppg1a1aUlJREVVVV7Nix40Pnr1+/Pj7zmc/E5MmTo7KyMh544IH44x//OKwNAwAAjKS8o2jz5s1RX18fDQ0NsWvXrpg9e3bU1tbGkSNHhpz/wgsvxKpVq6KhoSH27NkTzz33XGzevDkeeuihc948AADAuco7ip566qn42te+FsuXL4/Pfe5zsWHDhrjwwgvjRz/60ZDzX3/99bjpppvi9ttvj1mzZsVtt90Wixcv/si7SwAAAGMhryjq7e2NnTt3Rk1NzZ+foLAwampqorW1dcg1N954Y+zcuXMggg4cOBDbtm2LL33pS2e8Tk9PT3R1dQ16AAAAjIZJ+Uw+duxY9PX1RXl5+aDx8vLy2Lt375Brbr/99jh27Fh84QtfiCzL4tSpU3H33Xd/6NvnGhsb47HHHstnawAAAMMy6p8+t3379li7dm0888wzsWvXrvjZz34WW7dujccff/yMa1avXh2dnZ0Dj8OHD4/2NgEAgETldado6tSpUVRUFB0dHYPGOzo6oqKiYsg1jz76aCxZsiRWrFgRERHXXHNNdHd3x1133RUPP/xwFBae3mW5XC5yuVw+WwMAABiWvO4UFRcXx9y5c6OlpWVgrL+/P1paWqK6unrINe+9995p4VNUVBQREVmW5btfAACAEZXXnaKIiPr6+li2bFnMmzcv5s+fH+vXr4/u7u5Yvnx5REQsXbo0Zs6cGY2NjRERsXDhwnjqqafiuuuui6qqqti/f388+uijsXDhwoE4AgAAGC95R1FdXV0cPXo01qxZE+3t7TFnzpxobm4e+PCFQ4cODboz9Mgjj0RBQUE88sgj8c4778Rf/dVfxcKFC+N73/veyL0KAACAYSrIJsB72Lq6uqKsrCw6OzujtLR0vLcDAACMk9Fog1H/9DkAAICPM1EEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNKGFUVNTU0xa9asKCkpiaqqqtixY8eHzj9+/HisXLkypk+fHrlcLq688srYtm3bsDYMAAAwkiblu2Dz5s1RX18fGzZsiKqqqli/fn3U1tbGvn37Ytq0aafN7+3tjS9+8Ysxbdq0ePHFF2PmzJnx+9//Pi6++OKR2D8AAMA5KciyLMtnQVVVVdxwww3x9NNPR0REf39/VFZWxr333hurVq06bf6GDRviX//1X2Pv3r1xwQUXDGuTXV1dUVZWFp2dnVFaWjqs5wAAACa+0WiDvN4+19vbGzt37oyampo/P0FhYdTU1ERra+uQa37+859HdXV1rFy5MsrLy+Pqq6+OtWvXRl9f3xmv09PTE11dXYMeAAAAoyGvKDp27Fj09fVFeXn5oPHy8vJob28fcs2BAwfixRdfjL6+vti2bVs8+uij8eSTT8Z3v/vdM16nsbExysrKBh6VlZX5bBMAAOCsjfqnz/X398e0adPi2Wefjblz50ZdXV08/PDDsWHDhjOuWb16dXR2dg48Dh8+PNrbBAAAEpXXBy1MnTo1ioqKoqOjY9B4R0dHVFRUDLlm+vTpccEFF0RRUdHA2Gc/+9lob2+P3t7eKC4uPm1NLpeLXC6Xz9YAAACGJa87RcXFxTF37txoaWkZGOvv74+Wlpaorq4ecs1NN90U+/fvj/7+/oGxt956K6ZPnz5kEAEAAIylvN8+V19fHxs3bowf//jHsWfPnvj6178e3d3dsXz58oiIWLp0aaxevXpg/te//vV4991347777ou33nortm7dGmvXro2VK1eO3KsAAAAYpry/p6iuri6OHj0aa9asifb29pgzZ040NzcPfPjCoUOHorDwz61VWVkZL7/8cjzwwANx7bXXxsyZM+O+++6LBx98cOReBQAAwDDl/T1F48H3FAEAABEfg+8pAgAAON+IIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkDSuKmpqaYtasWVFSUhJVVVWxY8eOs1q3adOmKCgoiEWLFg3nsgAAACMu7yjavHlz1NfXR0NDQ+zatStmz54dtbW1ceTIkQ9dd/DgwfjmN78ZN99887A3CwAAMNLyjqKnnnoqvva1r8Xy5cvjc5/7XGzYsCEuvPDC+NGPfnTGNX19fXHHHXfEY489Fpdddtk5bRgAAGAk5RVFvb29sXPnzqipqfnzExQWRk1NTbS2tp5x3Xe+852YNm1a3HnnnWd1nZ6enujq6hr0AAAAGA15RdGxY8eir68vysvLB42Xl5dHe3v7kGt+85vfxHPPPRcbN2486+s0NjZGWVnZwKOysjKfbQIAAJy1Uf30uRMnTsSSJUti48aNMXXq1LNet3r16ujs7Bx4HD58eBR3CQAApGxSPpOnTp0aRUVF0dHRMWi8o6MjKioqTpv/u9/9Lg4ePBgLFy4cGOvv7//ThSdNin379sXll19+2rpcLhe5XC6frQEAAAxLXneKiouLY+7cudHS0jIw1t/fHy0tLVFdXX3a/KuuuireeOONaGtrG3h8+ctfjltvvTXa2tq8LQ4AABh3ed0pioior6+PZcuWxbx582L+/Pmxfv366O7ujuXLl0dExNKlS2PmzJnR2NgYJSUlcfXVVw9af/HFF0dEnDYOAAAwHvKOorq6ujh69GisWbMm2tvbY86cOdHc3Dzw4QuHDh2KwsJR/VUlAACAEVOQZVk23pv4KF1dXVFWVhadnZ1RWlo63tsBAADGyWi0gVs6AABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkbVhQ1NTXFrFmzoqSkJKqqqmLHjh1nnLtx48a4+eabY8qUKTFlypSoqan50PkAAABjKe8o2rx5c9TX10dDQ0Ps2rUrZs+eHbW1tXHkyJEh52/fvj0WL14cr776arS2tkZlZWXcdttt8c4775zz5gEAAM5VQZZlWT4Lqqqq4oYbboinn346IiL6+/ujsrIy7r333li1atVHru/r64spU6bE008/HUuXLj2ra3Z1dUVZWVl0dnZGaWlpPtsFAADOI6PRBnndKert7Y2dO3dGTU3Nn5+gsDBqamqitbX1rJ7jvffei5MnT8Yll1yS304BAABGwaR8Jh87diz6+vqivLx80Hh5eXns3bv3rJ7jwQcfjBkzZgwKq7/U09MTPT09Az93dXXls00AAICzNqafPrdu3brYtGlTbNmyJUpKSs44r7GxMcrKygYelZWVY7hLAAAgJXlF0dSpU6OoqCg6OjoGjXd0dERFRcWHrn3iiSdi3bp18ctf/jKuvfbaD527evXq6OzsHHgcPnw4n20CAACctbyiqLi4OObOnRstLS0DY/39/dHS0hLV1dVnXPf9738/Hn/88Whubo558+Z95HVyuVyUlpYOegAAAIyGvH6nKCKivr4+li1bFvPmzYv58+fH+vXro7u7O5YvXx4REUuXLo2ZM2dGY2NjRET8y7/8S6xZsyZeeOGFmDVrVrS3t0dExCc+8Yn4xCc+MYIvBQAAIH95R1FdXV0cPXo01qxZE+3t7TFnzpxobm4e+PCFQ4cORWHhn29A/eAHP4je3t74h3/4h0HP09DQEN/+9rfPbfcAAADnKO/vKRoPvqcIAACI+Bh8TxEAAMD5RhQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJG1YUdTU1BSzZs2KkpKSqKqqih07dnzo/J/+9Kdx1VVXRUlJSVxzzTWxbdu2YW0WAABgpOUdRZs3b476+vpoaGiIXbt2xezZs6O2tjaOHDky5PzXX389Fi9eHHfeeWfs3r07Fi1aFIsWLYo333zznDcPAABwrgqyLMvyWVBVVRU33HBDPP300xER0d/fH5WVlXHvvffGqlWrTptfV1cX3d3d8Ytf/GJg7G//9m9jzpw5sWHDhrO6ZldXV5SVlUVnZ2eUlpbms10AAOA8MhptMCmfyb29vbFz585YvXr1wFhhYWHU1NREa2vrkGtaW1ujvr5+0FhtbW289NJLZ7xOT09P9PT0DPzc2dkZEX/6PwAAAEjXB02Q572dD5VXFB07diz6+vqivLx80Hh5eXns3bt3yDXt7e1Dzm9vbz/jdRobG+Oxxx47bbyysjKf7QIAAOep//mf/4mysrIRea68omisrF69etDdpePHj8enP/3pOHTo0Ii9cBhKV1dXVFZWxuHDh71Vk1HlrDFWnDXGirPGWOns7IxLL700LrnkkhF7zryiaOrUqVFUVBQdHR2Dxjs6OqKiomLINRUVFXnNj4jI5XKRy+VOGy8rK/OHjDFRWlrqrDEmnDXGirPGWHHWGCuFhSP37UJ5PVNxcXHMnTs3WlpaBsb6+/ujpaUlqqurh1xTXV09aH5ExCuvvHLG+QAAAGMp77fP1dfXx7Jly2LevHkxf/78WL9+fXR3d8fy5csjImLp0qUxc+bMaGxsjIiI++67L2655ZZ48sknY8GCBbFp06b47W9/G88+++zIvhIAAIBhyDuK6urq4ujRo7FmzZpob2+POXPmRHNz88CHKRw6dGjQrawbb7wxXnjhhXjkkUfioYceir/5m7+Jl156Ka6++uqzvmYul4uGhoYh31IHI8lZY6w4a4wVZ42x4qwxVkbjrOX9PUUAAADnk5H77SQAAIAJSBQBAABJE0UAAEDSRBEAAJC0j00UNTU1xaxZs6KkpCSqqqpix44dHzr/pz/9aVx11VVRUlIS11xzTWzbtm2MdspEl89Z27hxY9x8880xZcqUmDJlStTU1Hzk2YQP5Pv32gc2bdoUBQUFsWjRotHdIOeNfM/a8ePHY+XKlTF9+vTI5XJx5ZVX+ucoZyXfs7Z+/fr4zGc+E5MnT47Kysp44IEH4o9//OMY7ZaJ6LXXXouFCxfGjBkzoqCgIF566aWPXLN9+/a4/vrrI5fLxRVXXBHPP/983tf9WETR5s2bo76+PhoaGmLXrl0xe/bsqK2tjSNHjgw5//XXX4/FixfHnXfeGbt3745FixbFokWL4s033xzjnTPR5HvWtm/fHosXL45XX301Wltbo7KyMm677bZ45513xnjnTDT5nrUPHDx4ML75zW/GzTffPEY7ZaLL96z19vbGF7/4xTh48GC8+OKLsW/fvti4cWPMnDlzjHfORJPvWXvhhRdi1apV0dDQEHv27InnnnsuNm/eHA899NAY75yJpLu7O2bPnh1NTU1nNf/tt9+OBQsWxK233hptbW1x//33x4oVK+Lll1/O78LZx8D8+fOzlStXDvzc19eXzZgxI2tsbBxy/le+8pVswYIFg8aqqqqyf/qnfxrVfTLx5XvW/tKpU6eyiy66KPvxj388WlvkPDGcs3bq1KnsxhtvzH74wx9my5Yty/7+7/9+DHbKRJfvWfvBD36QXXbZZVlvb+9YbZHzRL5nbeXKldnf/d3fDRqrr6/PbrrpplHdJ+ePiMi2bNnyoXO+9a1vZZ///OcHjdXV1WW1tbV5XWvc7xT19vbGzp07o6amZmCssLAwampqorW1dcg1ra2tg+ZHRNTW1p5xPkQM76z9pffeey9OnjwZl1xyyWhtk/PAcM/ad77znZg2bVrceeedY7FNzgPDOWs///nPo7q6OlauXBnl5eVx9dVXx9q1a6Ovr2+sts0ENJyzduONN8bOnTsH3mJ34MCB2LZtW3zpS18akz2ThpHqgkkjuanhOHbsWPT19UV5efmg8fLy8ti7d++Qa9rb24ec397ePmr7ZOIbzln7Sw8++GDMmDHjtD988P8bzln7zW9+E88991y0tbWNwQ45XwznrB04cCB+/etfxx133BHbtm2L/fv3xz333BMnT56MhoaGsdg2E9Bwztrtt98ex44diy984QuRZVmcOnUq7r77bm+fY0SdqQu6urri/fffj8mTJ5/V84z7nSKYKNatWxebNm2KLVu2RElJyXhvh/PIiRMnYsmSJbFx48aYOnXqeG+H81x/f39MmzYtnn322Zg7d27U1dXFww8/HBs2bBjvrXGe2b59e6xduzaeeeaZ2LVrV/zsZz+LrVu3xuOPPz7eW4PTjPudoqlTp0ZRUVF0dHQMGu/o6IiKiooh11RUVOQ1HyKGd9Y+8MQTT8S6deviV7/6VVx77bWjuU3OA/metd/97ndx8ODBWLhw4cBYf39/RERMmjQp9u3bF5dffvnobpoJaTh/r02fPj0uuOCCKCoqGhj77Gc/G+3t7dHb2xvFxcWjumcmpuGctUcffTSWLFkSK1asiIiIa665Jrq7u+Ouu+6Khx9+OAoL/bd5zt2ZuqC0tPSs7xJFfAzuFBUXF8fcuXOjpaVlYKy/vz9aWlqiurp6yDXV1dWD5kdEvPLKK2ecDxHDO2sREd///vfj8ccfj+bm5pg3b95YbJUJLt+zdtVVV8Ubb7wRbW1tA48vf/nLA5+kU1lZOZbbZwIZzt9rN910U+zfv38gvCMi3nrrrZg+fbog4oyGc9bee++908Lngxj/0+/Qw7kbsS7I7zMgRsemTZuyXC6XPf/889l///d/Z3fddVd28cUXZ+3t7VmWZdmSJUuyVatWDcz/z//8z2zSpEnZE088ke3ZsydraGjILrjgguyNN94Yr5fABJHvWVu3bl1WXFycvfjii9kf/vCHgceJEyfG6yUwQeR71v6ST5/jbOV71g4dOpRddNFF2Te+8Y1s37592S9+8Yts2rRp2Xe/+93xeglMEPmetYaGhuyiiy7K/uM//iM7cOBA9stf/jK7/PLLs6985Svj9RKYAE6cOJHt3r072717dxYR2VNPPZXt3r07+/3vf59lWZatWrUqW7JkycD8AwcOZBdeeGH2z//8z9mePXuypqamrKioKGtubs7ruh+LKMqyLPu3f/u37NJLL82Ki4uz+fPnZ//1X/818L/dcsst2bJlywbN/8lPfpJdeeWVWXFxcfb5z38+27p16xjvmIkqn7P26U9/OouI0x4NDQ1jv3EmnHz/Xvv/iSLyke9Ze/3117Oqqqosl8tll112Wfa9730vO3Xq1Bjvmokon7N28uTJ7Nvf/nZ2+eWXZyUlJVllZWV2zz33ZP/7v/879htnwnj11VeH/HevD87WsmXLsltuueW0NXPmzMmKi4uzyy67LPv3f//3vK9bkGXuXwIAAOka998pAgAAGE+iCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABIWt5R9Nprr8XChQtjxowZUVBQEC+99NJHrtm+fXtcf/31kcvl4oorrojnn39+GFsFAAAYeXlHUXd3d8yePTuamprOav7bb78dCxYsiFtvvTXa2tri/vvvjxUrVsTLL7+c92YBAABGWkGWZdmwFxcUxJYtW2LRokVnnPPggw/G1q1b48033xwY++pXvxrHjx+P5ubm4V4aAABgREwa7Qu0trZGTU3NoLHa2tq4//77z7imp6cnenp6Bn7u7++Pd999Nz75yU9GQUHBaG0VAAD4mMuyLE6cOBEzZsyIwsKR+YiEUY+i9vb2KC8vHzRWXl4eXV1d8f7778fkyZNPW9PY2BiPPfbYaG8NAACYoA4fPhyf+tSnRuS5Rj2KhmP16tVRX18/8HNnZ2dceumlcfjw4SgtLR3HnQEAAP+vvbuPrbusGz/+WTt6CpGWcc91D7+DExRBgQ03VgsuBFNpApnuD0PvQbZlARGdBGlUNh5WEV2nAlkihYUJYmJw0wXQuKWIlcUoNYvbmkDcRnDOLcR2m7p2Fm1Z+/39Yah3XQc7XR/ortcrOX/s8rrO9zrmYvj2ex7GUmdnZ+Tz+Tj77LOH7TlHPIqmTp0a7e3tA8ba29ujrKxs0LtEERG5XC5yudxx42VlZaIIAAAY1o/VjPjvFFVVVUVzc/OAsRdeeCGqqqpG+tIAAADvqOAo+sc//hGtra3R2toaEf/+yu3W1tbYv39/RPz7rW9Llizpn3/bbbfF3r1746tf/Wrs3r07Hn300fjxj38cd9555/C8AgAAgFNQcBT9/ve/j8svvzwuv/zyiIioq6uLyy+/PFatWhUREX/5y1/6Ayki4v3vf39s3rw5XnjhhZg1a1Y89NBD8b3vfS9qamqG6SUAAAAM3Sn9TtFo6ezsjPLy8ujo6PCZIgAASNhItMGIf6YIAADg3UwUAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJC0IUVRY2NjzJw5M0pLS6OysjK2bdv2tvPXrl0bH/rQh+LMM8+MfD4fd955Z/zrX/8a0oYBAACGU8FRtHHjxqirq4v6+vrYsWNHzJo1K2pqauLgwYODzn/66adjxYoVUV9fH7t27YonnngiNm7cGHffffcpbx4AAOBUFRxFDz/8cHz2s5+NZcuWxYc//OFYt25dnHXWWfHkk08OOv+ll16Kq666Km688caYOXNmXHvttbFo0aJ3vLsEAAAwGgqKop6enti+fXtUV1f/5wmKiqK6ujpaWloGXXPllVfG9u3b+yNo7969sWXLlrjuuutOYdsAAADDY2Ihkw8fPhy9vb1RUVExYLyioiJ279496Jobb7wxDh8+HB//+Mcjy7I4duxY3HbbbW/79rnu7u7o7u7u/3NnZ2ch2wQAADhpI/7tc1u3bo3Vq1fHo48+Gjt27IhnnnkmNm/eHA888MAJ1zQ0NER5eXn/I5/Pj/Q2AQCARE3Isiw72ck9PT1x1llnxaZNm2LhwoX940uXLo0jR47ET3/60+PWzJ8/Pz72sY/Fd77znf6xH/7wh3HrrbfGP/7xjygqOr7LBrtTlM/no6OjI8rKyk52uwAAwGmms7MzysvLh7UNCrpTVFJSEnPmzInm5ub+sb6+vmhubo6qqqpB17zxxhvHhU9xcXFERJyox3K5XJSVlQ14AAAAjISCPlMUEVFXVxdLly6NuXPnxrx582Lt2rXR1dUVy5Yti4iIJUuWxIwZM6KhoSEiIhYsWBAPP/xwXH755VFZWRmvvfZa3HfffbFgwYL+OAIAABgrBUdRbW1tHDp0KFatWhVtbW0xe/bsaGpq6v/yhf379w+4M3TvvffGhAkT4t57743XX3893vve98aCBQvim9/85vC9CgAAgCEq6DNFY2Uk3jcIAACMP2P+mSIAAIDTjSgCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASNqQoqixsTFmzpwZpaWlUVlZGdu2bXvb+UeOHInly5fHtGnTIpfLxYUXXhhbtmwZ0oYBAACG08RCF2zcuDHq6upi3bp1UVlZGWvXro2amprYs2dPTJky5bj5PT098clPfjKmTJkSmzZtihkzZsSf//znOOecc4Zj/wAAAKdkQpZlWSELKisr44orrohHHnkkIiL6+voin8/H7bffHitWrDhu/rp16+I73/lO7N69O84444whbbKzszPKy8ujo6MjysrKhvQcAADA+DcSbVDQ2+d6enpi+/btUV1d/Z8nKCqK6urqaGlpGXTNz372s6iqqorly5dHRUVFXHLJJbF69ero7e094XW6u7ujs7NzwAMAAGAkFBRFhw8fjt7e3qioqBgwXlFREW1tbYOu2bt3b2zatCl6e3tjy5Ytcd9998VDDz0U3/jGN054nYaGhigvL+9/5PP5QrYJAABw0kb82+f6+vpiypQp8fjjj8ecOXOitrY27rnnnli3bt0J16xcuTI6Ojr6HwcOHBjpbQIAAIkq6IsWJk+eHMXFxdHe3j5gvL29PaZOnTrommnTpsUZZ5wRxcXF/WMXX3xxtLW1RU9PT5SUlBy3JpfLRS6XK2RrAAAAQ1LQnaKSkpKYM2dONDc394/19fVFc3NzVFVVDbrmqquuitdeey36+vr6x1599dWYNm3aoEEEAAAwmgp++1xdXV2sX78+fvCDH8SuXbvi85//fHR1dcWyZcsiImLJkiWxcuXK/vmf//zn429/+1vccccd8eqrr8bmzZtj9erVsXz58uF7FQAAAENU8O8U1dbWxqFDh2LVqlXR1tYWs2fPjqampv4vX9i/f38UFf2ntfL5fDz//PNx5513xmWXXRYzZsyIO+64I+66667hexUAAABDVPDvFI0Fv1MEAABEvAt+pwgAAOB0I4oAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkjakKGpsbIyZM2dGaWlpVFZWxrZt205q3YYNG2LChAmxcOHCoVwWAABg2BUcRRs3boy6urqor6+PHTt2xKxZs6KmpiYOHjz4tuv27dsXX/7yl2P+/PlD3iwAAMBwKziKHn744fjsZz8by5Ytiw9/+MOxbt26OOuss+LJJ5884Zre3t646aab4v7774/zzz//lDYMAAAwnAqKop6enti+fXtUV1f/5wmKiqK6ujpaWlpOuO7rX/96TJkyJW6++eaTuk53d3d0dnYOeAAAAIyEgqLo8OHD0dvbGxUVFQPGKyoqoq2tbdA1v/nNb+KJJ56I9evXn/R1Ghoaory8vP+Rz+cL2SYAAMBJG9Fvnzt69GgsXrw41q9fH5MnTz7pdStXroyOjo7+x4EDB0ZwlwAAQMomFjJ58uTJUVxcHO3t7QPG29vbY+rUqcfN/+Mf/xj79u2LBQsW9I/19fX9+8ITJ8aePXviggsuOG5dLpeLXC5XyNYAAACGpKA7RSUlJTFnzpxobm7uH+vr64vm5uaoqqo6bv5FF10UL7/8crS2tvY/PvWpT8U111wTra2t3hYHAACMuYLuFEVE1NXVxdKlS2Pu3Lkxb968WLt2bXR1dcWyZcsiImLJkiUxY8aMaGhoiNLS0rjkkksGrD/nnHMiIo4bBwAAGAsFR1FtbW0cOnQoVq1aFW1tbTF79uxoamrq//KF/fv3R1HRiH5UCQAAYNhMyLIsG+tNvJPOzs4oLy+Pjo6OKCsrG+vtAAAAY2Qk2sAtHQAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkDSmKGhsbY+bMmVFaWhqVlZWxbdu2E85dv359zJ8/PyZNmhSTJk2K6urqt50PAAAwmgqOoo0bN0ZdXV3U19fHjh07YtasWVFTUxMHDx4cdP7WrVtj0aJF8eKLL0ZLS0vk8/m49tpr4/XXXz/lzQMAAJyqCVmWZYUsqKysjCuuuCIeeeSRiIjo6+uLfD4ft99+e6xYseId1/f29sakSZPikUceiSVLlpzUNTs7O6O8vDw6OjqirKyskO0CAACnkZFog4LuFPX09MT27dujurr6P09QVBTV1dXR0tJyUs/xxhtvxJtvvhnnnntuYTsFAAAYARMLmXz48OHo7e2NioqKAeMVFRWxe/fuk3qOu+66K6ZPnz4grP5bd3d3dHd39/+5s7OzkG0CAACctFH99rk1a9bEhg0b4tlnn43S0tITzmtoaIjy8vL+Rz6fH8VdAgAAKSkoiiZPnhzFxcXR3t4+YLy9vT2mTp36tmsffPDBWLNmTfziF7+Iyy677G3nrly5Mjo6OvofBw4cKGSbAAAAJ62gKCopKYk5c+ZEc3Nz/1hfX180NzdHVVXVCdd9+9vfjgceeCCamppi7ty573idXC4XZWVlAx4AAAAjoaDPFEVE1NXVxdKlS2Pu3Lkxb968WLt2bXR1dcWyZcsiImLJkiUxY8aMaGhoiIiIb33rW7Fq1ap4+umnY+bMmdHW1hYREe95z3viPe95zzC+FAAAgMIVHEW1tbVx6NChWLVqVbS1tcXs2bOjqamp/8sX9u/fH0VF/7kB9dhjj0VPT0985jOfGfA89fX18bWvfe3Udg8AAHCKCv6dorHgd4oAAICId8HvFAEAAJxuRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0oYURY2NjTFz5swoLS2NysrK2LZt29vO/8lPfhIXXXRRlJaWxqWXXhpbtmwZ0mYBAACGW8FRtHHjxqirq4v6+vrYsWNHzJo1K2pqauLgwYODzn/ppZdi0aJFcfPNN8fOnTtj4cKFsXDhwnjllVdOefMAAACnakKWZVkhCyorK+OKK66IRx55JCIi+vr6Ip/Px+233x4rVqw4bn5tbW10dXXFz3/+8/6xj33sYzF79uxYt27dSV2zs7MzysvLo6OjI8rKygrZLgAAcBoZiTaYWMjknp6e2L59e6xcubJ/rKioKKqrq6OlpWXQNS0tLVFXVzdgrKamJp577rkTXqe7uzu6u7v7/9zR0RER//4vAAAASNdbTVDgvZ23VVAUHT58OHp7e6OiomLAeEVFRezevXvQNW1tbYPOb2trO+F1Ghoa4v777z9uPJ/PF7JdAADgNPXXv/41ysvLh+W5Coqi0bJy5coBd5eOHDkS73vf+2L//v3D9sJhMJ2dnZHP5+PAgQPeqsmIctYYLc4ao8VZY7R0dHTEeeedF+eee+6wPWdBUTR58uQoLi6O9vb2AePt7e0xderUQddMnTq1oPkREblcLnK53HHj5eXl/iFjVJSVlTlrjApnjdHirDFanDVGS1HR8P26UEHPVFJSEnPmzInm5ub+sb6+vmhubo6qqqpB11RVVQ2YHxHxwgsvnHA+AADAaCr47XN1dXWxdOnSmDt3bsybNy/Wrl0bXV1dsWzZsoiIWLJkScyYMSMaGhoiIuKOO+6Iq6++Oh566KG4/vrrY8OGDfH73/8+Hn/88eF9JQAAAENQcBTV1tbGoUOHYtWqVdHW1hazZ8+Opqam/i9T2L9//4BbWVdeeWU8/fTTce+998bdd98dH/zgB+O5556LSy655KSvmcvlor6+ftC31MFwctYYLc4ao8VZY7Q4a4yWkThrBf9OEQAAwOlk+D6dBAAAMA6JIgAAIGmiCAAASJooAgAAkvauiaLGxsaYOXNmlJaWRmVlZWzbtu1t5//kJz+Jiy66KEpLS+PSSy+NLVu2jNJOGe8KOWvr16+P+fPnx6RJk2LSpElRXV39jmcT3lLo32tv2bBhQ0yYMCEWLlw4shvktFHoWTty5EgsX748pk2bFrlcLi688EL/HuWkFHrW1q5dGx/60IfizDPPjHw+H3feeWf861//GqXdMh79+te/jgULFsT06dNjwoQJ8dxzz73jmq1bt8ZHP/rRyOVy8YEPfCCeeuqpgq/7roiijRs3Rl1dXdTX18eOHTti1qxZUVNTEwcPHhx0/ksvvRSLFi2Km2++OXbu3BkLFy6MhQsXxiuvvDLKO2e8KfSsbd26NRYtWhQvvvhitLS0RD6fj2uvvTZef/31Ud45402hZ+0t+/btiy9/+csxf/78Udop412hZ62npyc++clPxr59+2LTpk2xZ8+eWL9+fcyYMWOUd854U+hZe/rpp2PFihVRX18fu3btiieeeCI2btwYd9999yjvnPGkq6srZs2aFY2NjSc1/09/+lNcf/31cc0110Rra2t86UtfiltuuSWef/75wi6cvQvMmzcvW758ef+fe3t7s+nTp2cNDQ2Dzr/hhhuy66+/fsBYZWVl9rnPfW5E98n4V+hZ+2/Hjh3Lzj777OwHP/jBSG2R08RQztqxY8eyK6+8Mvve976XLV26NPv0pz89CjtlvCv0rD322GPZ+eefn/X09IzWFjlNFHrWli9fnn3iE58YMFZXV5ddddVVI7pPTh8RkT377LNvO+erX/1q9pGPfGTAWG1tbVZTU1PQtcb8TlFPT09s3749qqur+8eKioqiuro6WlpaBl3T0tIyYH5ERE1NzQnnQ8TQztp/e+ONN+LNN9+Mc889d6S2yWlgqGft61//ekyZMiVuvvnm0dgmp4GhnLWf/exnUVVVFcuXL4+Kioq45JJLYvXq1dHb2zta22YcGspZu/LKK2P79u39b7Hbu3dvbNmyJa677rpR2TNpGK4umDicmxqKw4cPR29vb1RUVAwYr6ioiN27dw+6pq2tbdD5bW1tI7ZPxr+hnLX/dtddd8X06dOP+4cP/q+hnLXf/OY38cQTT0Rra+so7JDTxVDO2t69e+NXv/pV3HTTTbFly5Z47bXX4gtf+EK8+eabUV9fPxrbZhwaylm78cYb4/Dhw/Hxj388siyLY8eOxW233ebtcwyrE3VBZ2dn/POf/4wzzzzzpJ5nzO8UwXixZs2a2LBhQzz77LNRWlo61tvhNHL06NFYvHhxrF+/PiZPnjzW2+E019fXF1OmTInHH3885syZE7W1tXHPPffEunXrxnprnGa2bt0aq1evjkcffTR27NgRzzzzTGzevDkeeOCBsd4aHGfM7xRNnjw5iouLo729fcB4e3t7TJ06ddA1U6dOLWg+RAztrL3lwQcfjDVr1sQvf/nLuOyyy0Zym5wGCj1rf/zjH2Pfvn2xYMGC/rG+vr6IiJg4cWLs2bMnLrjggpHdNOPSUP5emzZtWpxxxhlRXFzcP3bxxRdHW1tb9PT0RElJyYjumfFpKGftvvvui8WLF8ctt9wSERGXXnppdHV1xa233hr33HNPFBX5/+Y5dSfqgrKyspO+SxTxLrhTVFJSEnPmzInm5ub+sb6+vmhubo6qqqpB11RVVQ2YHxHxwgsvnHA+RAztrEVEfPvb344HHnggmpqaYu7cuaOxVca5Qs/aRRddFC+//HK0trb2Pz71qU/1f5NOPp8fze0zjgzl77WrrroqXnvttf7wjoh49dVXY9q0aYKIExrKWXvjjTeOC5+3Yvzfn6GHUzdsXVDYd0CMjA0bNmS5XC576qmnsj/84Q/Zrbfemp1zzjlZW1tblmVZtnjx4mzFihX983/7299mEydOzB588MFs165dWX19fXbGGWdkL7/88li9BMaJQs/amjVrspKSkmzTpk3ZX/7yl/7H0aNHx+olME4Uetb+m2+f42QVetb279+fnX322dkXv/jFbM+ePdnPf/7zbMqUKdk3vvGNsXoJjBOFnrX6+vrs7LPPzn70ox9le/fuzX7xi19kF1xwQXbDDTeM1UtgHDh69Gi2c+fObOfOnVlEZA8//HC2c+fO7M9//nOWZVm2YsWKbPHixf3z9+7dm5111lnZV77ylWzXrl1ZY2NjVlxcnDU1NRV03XdFFGVZln33u9/NzjvvvKykpCSbN29e9rvf/a7/P7v66quzpUuXDpj/4x//OLvwwguzkpKS7CMf+Ui2efPmUd4x41UhZ+1973tfFhHHPerr60d/44w7hf699n+JIgpR6Fl76aWXssrKyiyXy2Xnn39+9s1vfjM7duzYKO+a8aiQs/bmm29mX/va17ILLrggKy0tzfL5fPaFL3wh+/vf/z76G2fcePHFFwf9315vna2lS5dmV1999XFrZs+enZWUlGTnn39+9v3vf7/g607IMvcvAQCAdI35Z4oAAADGkigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJJWcBT9+te/jgULFsT06dNjwoQJ8dxzz73jmq1bt8ZHP/rRyOVy8YEPfCCeeuqpIWwVAABg+BUcRV1dXTFr1qxobGw8qfl/+tOf4vrrr49rrrkmWltb40tf+lLccsst8fzzzxe8WQAAgOE2IcuybMiLJ0yIZ599NhYuXHjCOXfddVds3rw5Xnnllf6x//3f/40jR45EU1PTUC8NAAAwLCaO9AVaWlqiurp6wFhNTU186UtfOuGa7u7u6O7u7v9zX19f/O1vf4v/+Z//iQkTJozUVgEAgHe5LMvi6NGjMX369CgqGp6vSBjxKGpra4uKiooBYxUVFdHZ2Rn//Oc/48wzzzxuTUNDQ9x///0jvTUAAGCcOnDgQPy///f/huW5RjyKhmLlypVRV1fX/+eOjo4477zz4sCBA1FWVjaGOwMAAMZSZ2dn5PP5OPvss4ftOUc8iqZOnRrt7e0Dxtrb26OsrGzQu0QREblcLnK53HHjZWVloggAABjWj9WM+O8UVVVVRXNz84CxF154Iaqqqkb60gAAAO+o4Cj6xz/+Ea2trdHa2hoR//7K7dbW1ti/f39E/Putb0uWLOmff9ttt8XevXvjq1/9auzevTseffTR+PGPfxx33nnn8LwCAACAU1BwFP3+97+Pyy+/PC6//PKIiKirq4vLL788Vq1aFRERf/nLX/oDKSLi/e9/f2zevDleeOGFmDVrVjz00EPxve99L2pqaobpJQAAAAzdKf1O0Wjp7OyM8vLy6Ojo8JkiAABI2Ei0wYh/pggAAODdTBQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkLQhRVFjY2PMnDkzSktLo7KyMrZt2/a289euXRsf+tCH4swzz4x8Ph933nln/Otf/xrShgEAAIZTwVG0cePGqKuri/r6+tixY0fMmjUrampq4uDBg4POf/rpp2PFihVRX18fu3btiieeeCI2btwYd9999ylvHgAA4FQVHEUPP/xwfPazn41ly5bFhz/84Vi3bl2cddZZ8eSTTw46/6WXXoqrrroqbrzxxpg5c2Zce+21sWjRone8uwQAADAaCoqinp6e2L59e1RXV//nCYqKorq6OlpaWgZdc+WVV8b27dv7I2jv3r2xZcuWuO66605h2wAAAMNjYiGTDx8+HL29vVFRUTFgvKKiInbv3j3omhtvvDEOHz4cH//4xyPLsjh27Fjcdtttb/v2ue7u7uju7u7/c2dnZyHbBAAAOGkj/u1zW7dujdWrV8ejjz4aO3bsiGeeeSY2b94cDzzwwAnXNDQ0RHl5ef8jn8+P9DYBAIBETciyLDvZyT09PXHWWWfFpk2bYuHChf3jS5cujSNHjsRPf/rT49bMnz8/Pvaxj8V3vvOd/rEf/vCHceutt8Y//vGPKCo6vssGu1OUz+ejo6MjysrKTna7AADAaaazszPKy8uHtQ0KulNUUlISc+bMiebm5v6xvr6+aG5ujqqqqkHXvPHGG8eFT3FxcUREnKjHcrlclJWVDXgAAACMhII+UxQRUVdXF0uXLo25c+fGvHnzYu3atdHV1RXLli2LiIglS5bEjBkzoqGhISIiFixYEA8//HBcfvnlUVlZGa+99lrcd999sWDBgv44AgAAGCsFR1FtbW0cOnQoVq1aFW1tbTF79uxoamrq//KF/fv3D7gzdO+998aECRPi3nvvjddffz3e+973xoIFC+Kb3/zm8L0KAACAISroM0VjZSTeNwgAAIw/Y/6ZIgAAgNONKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABI2pCiqLGxMWbOnBmlpaVRWVkZ27Zte9v5R44cieXLl8e0adMil8vFhRdeGFu2bBnShgEAAIbTxEIXbNy4Merq6mLdunVRWVkZa9eujZqamtizZ09MmTLluPk9PT3xyU9+MqZMmRKbNm2KGTNmxJ///Oc455xzhmP/AAAAp2RClmVZIQsqKyvjiiuuiEceeSQiIvr6+iKfz8ftt98eK1asOG7+unXr4jvf+U7s3r07zjjjjCFtsrOzM8rLy6OjoyPKysqG9BwAAMD4NxJtUNDb53p6emL79u1RXV39nycoKorq6upoaWkZdM3PfvazqKqqiuXLl0dFRUVccsklsXr16ujt7T3hdbq7u6Ozs3PAAwAAYCQUFEWHDx+O3t7eqKioGDBeUVERbW1tg67Zu3dvbNq0KXp7e2PLli1x3333xUMPPRTf+MY3TnidhoaGKC8v73/k8/lCtgkAAHDSRvzb5/r6+mLKlCnx+OOPx5w5c6K2tjbuueeeWLdu3QnXrFy5Mjo6OvofBw4cGOltAgAAiSroixYmT54cxcXF0d7ePmC8vb09pk6dOuiaadOmxRlnnBHFxcX9YxdffHG0tbVFT09PlJSUHLcml8tFLpcrZGsAAABDUtCdopKSkpgzZ040Nzf3j/X19UVzc3NUVVUNuuaqq66K1157Lfr6+vrHXn311Zg2bdqgQQQAADCaCn77XF1dXaxfvz5+8IMfxK5du+Lzn/98dHV1xbJlyyIiYsmSJbFy5cr++Z///Ofjb3/7W9xxxx3x6quvxubNm2P16tWxfPny4XsVAAAAQ1Tw7xTV1tbGoUOHYtWqVdHW1hazZ8+Opqam/i9f2L9/fxQV/ae18vl8PP/883HnnXfGZZddFjNmzIg77rgj7rrrruF7FQAAAENU8O8UjQW/UwQAAES8C36nCAAA4HQjigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSNqQoamxsjJkzZ0ZpaWlUVlbGtm3bTmrdhg0bYsKECbFw4cKhXBYAAGDYFRxFGzdujLq6uqivr48dO3bErFmzoqamJg4ePPi26/bt2xdf/vKXY/78+UPeLAAAwHArOIoefvjh+OxnPxvLli2LD3/4w7Fu3bo466yz4sknnzzhmt7e3rjpppvi/vvvj/PPP/+UNgwAADCcCoqinp6e2L59e1RXV//nCYqKorq6OlpaWk647utf/3pMmTIlbr755pO6Tnd3d3R2dg54AAAAjISCoujw4cPR29sbFRUVA8YrKiqira1t0DW/+c1v4oknnoj169ef9HUaGhqivLy8/5HP5wvZJgAAwEkb0W+fO3r0aCxevDjWr18fkydPPul1K1eujI6Ojv7HgQMHRnCXAABAyiYWMnny5MlRXFwc7e3tA8bb29tj6tSpx83/4x//GPv27YsFCxb0j/X19f37whMnxp49e+KCCy44bl0ul4tcLlfI1gAAAIakoDtFJSUlMWfOnGhubu4f6+vri+bm5qiqqjpu/kUXXRQvv/xytLa29j8+9alPxTXXXBOtra3eFgcAAIy5gu4URUTU1dXF0qVLY+7cuTFv3rxYu3ZtdHV1xbJlyyIiYsmSJTFjxoxoaGiI0tLSuOSSSwasP+eccyIijhsHAAAYCwVHUW1tbRw6dChWrVoVbW1tMXv27Ghqaur/8oX9+/dHUdGIflQJAABg2EzIsiwb6028k87OzigvL4+Ojo4oKysb6+0AAABjZCTawC0dAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBpoggAAEiaKAIAAJImigAAgKQNKYoaGxtj5syZUVpaGpWVlbFt27YTzl2/fn3Mnz8/Jk2aFJMmTYrq6uq3nQ8AADCaCo6ijRs3Rl1dXdTX18eOHTti1qxZUVNTEwcPHhx0/tatW2PRokXx4osvRktLS+Tz+bj22mvj9ddfP+XNAwAAnKoJWZZlhSyorKyMK664Ih555JGIiOjr64t8Ph+33357rFix4h3X9/b2xqRJk+KRRx6JJUuWnNQ1Ozs7o7y8PDo6OqKsrKyQ7QIAAKeRkWiDgu4U9fT0xPbt26O6uvo/T1BUFNXV1dHS0nJSz/HGG2/Em2++Geeee25hOwUAABgBEwuZfPjw4ejt7Y2KiooB4xUVFbF79+6Teo677rorpk+fPiCs/lt3d3d0d3f3/7mzs7OQbQIAAJy0Uf32uTVr1sSGDRvi2WefjdLS0hPOa2hoiPLy8v5HPp8fxV0CAAApKSiKJk+eHMXFxdHe3j5gvL29PaZOnfq2ax988MFYs2ZN/OIXv4jLLrvsbeeuXLkyOjo6+h8HDhwoZJsAAAAnraAoKikpiTlz5kRzc3P/WF9fXzQ3N0dVVdUJ133729+OBx54IJqammLu3LnveJ1cLhdlZWUDHgAAACOhoM8URUTU1dXF0qVLY+7cuTFv3rxYu3ZtdHV1xbJlyyIiYsmSJTFjxoxoaGiIiIhvfetbsWrVqnj66adj5syZ0dbWFhER73nPe+I973nPML4UAACAwhUcRbW1tXHo0KFYtWpVtLW1xezZs6Opqan/yxf2798fRUX/uQH12GOPRU9PT3zmM58Z8Dz19fXxta997dR2DwAAcIoK/p2iseB3igAAgIh3we8UAQAAnG5EEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDSRBEAAJA0UQQAACRNFAEAAEkTRQAAQNJEEQAAkDRRBAAAJE0UAQAASRNFAABA0kQRAACQNFEEAAAkTRQBAABJE0UAAEDShhRFjY2NMXPmzCgtLY3KysrYtm3b287/yU9+EhdddFGUlpbGpZdeGlu2bBnSZgEAAIZbwVG0cePGqKuri/r6+tixY0fMmjUrampq4uDBg4POf+mll2LRokVx8803x86dO2PhwoWxcOHCeOWVV0558wAAAKdqQpZlWSELKisr44orrohHHnkkIiL6+voin8/H7bffHitWrDhufm1tbXR1dcXPf/7z/rGPfexjMXv27Fi3bt1JXbOzszPKy8ujo6MjysrKCtkuAABwGhmJNphYyOSenp7Yvn17rFy5sn+sqKgoqquro6WlZdA1LS0tUVdXN2CspqYmnnvuuRNep7u7O7q7u/v/3NHRERH//i8AAABI11tNUOC9nbdVUBQdPnw4ent7o6KiYsB4RUVF7N69e9A1bW1tg85va2s74XUaGhri/vvvP248n88Xsl0AAOA09de//jXKy8uH5bkKiqLRsnLlygF3l44cORLve9/7Yv/+/cP2wmEwnZ2dkc/n48CBA96qyYhy1hgtzhqjxVljtHR0dMR5550X55577rA9Z0FRNHny5CguLo729vYB4+3t7TF16tRB10ydOrWg+RERuVwucrnccePl5eX+IWNUlJWVOWuMCmeN0eKsMVqcNUZLUdHw/bpQQc9UUlISc+bMiebm5v6xvr6+aG5ujqqqqkHXVFVVDZgfEfHCCy+ccD4AAMBoKvjtc3V1dbF06dKYO3duzJs3L9auXRtdXV2xbNmyiIhYsmRJzJgxIxoaGiIi4o477oirr746Hnroobj++utjw4YN8fvf/z4ef/zx4X0lAAAAQ1BwFNXW1sahQ4di1apV0dbWFrNnz46mpqb+L1PYv3//gFtZV155ZTz99NNx7733xt133x0f/OAH47nnnotLLrnkpK+Zy+Wivr5+0LfUwXBy1hgtzhqjxVljtDhrjJaROGsF/04RAADA6WT4Pp0EAAAwDokiAAAgaaIIAABImigCAACS9q6JosbGxpg5c2aUlpZGZWVlbNu27W3n/+QnP4mLLrooSktL49JLL40tW7aM0k4Z7wo5a+vXr4/58+fHpEmTYtKkSVFdXf2OZxPeUujfa2/ZsGFDTJgwIRYuXDiyG+S0UehZO3LkSCxfvjymTZsWuVwuLrzwQv8e5aQUetbWrl0bH/rQh+LMM8+MfD4fd955Z/zrX/8apd0yHv3617+OBQsWxPTp02PChAnx3HPPveOarVu3xkc/+tHI5XLxgQ98IJ566qmCr/uuiKKNGzdGXV1d1NfXx44dO2LWrFlRU1MTBw8eHHT+Sy+9FIsWLYqbb745du7cGQsXLoyFCxfGK6+8Mso7Z7wp9Kxt3bo1Fi1aFC+++GK0tLREPp+Pa6+9Nl5//fVR3jnjTaFn7S379u2LL3/5yzF//vxR2injXaFnraenJz75yU/Gvn37YtOmTbFnz55Yv359zJgxY5R3znhT6Fl7+umnY8WKFVFfXx+7du2KJ554IjZu3Bh33333KO+c8aSrqytmzZoVjY2NJzX/T3/6U1x//fVxzTXXRGtra3zpS1+KW265JZ5//vnCLpy9C8ybNy9bvnx5/597e3uz6dOnZw0NDYPOv+GGG7Lrr79+wFhlZWX2uc99bkT3yfhX6Fn7b8eOHcvOPvvs7Ac/+MFIbZHTxFDO2rFjx7Irr7wy+973vpctXbo0+/SnPz0KO2W8K/SsPfbYY9n555+f9fT0jNYWOU0UetaWL1+efeITnxgwVldXl1111VUjuk9OHxGRPfvss28756tf/Wr2kY98ZMBYbW1tVlNTU9C1xvxOUU9PT2zfvj2qq6v7x4qKiqK6ujpaWloGXdPS0jJgfkRETU3NCedDxNDO2n9744034s0334xzzz13pLbJaWCoZ+3rX/96TJkyJW6++ebR2CangaGctZ/97GdRVVUVy5cvj4qKirjkkkti9erV0dvbO1rbZhwaylm78sorY/v27f1vsdu7d29s2bIlrrvuulHZM2kYri6YOJybGorDhw9Hb29vVFRUDBivqKiI3bt3D7qmra1t0PltbW0jtk/Gv6Gctf921113xfTp04/7hw/+r6Gctd/85jfxxBNPRGtr6yjskNPFUM7a3r1741e/+lXcdNNNsWXLlnjttdfiC1/4Qrz55ptRX18/GttmHBrKWbvxxhvj8OHD8fGPfzyyLItjx47Fbbfd5u1zDKsTdUFnZ2f885//jDPPPPOknmfM7xTBeLFmzZrYsGFDPPvss1FaWjrW2+E0cvTo0Vi8eHGsX78+Jk+ePNbb4TTX19cXU6ZMiccffzzmzJkTtbW1cc8998S6devGemucZrZu3RqrV6+ORx99NHbs2BHPPPNMbN68OR544IGx3hocZ8zvFE2ePDmKi4ujvb19wHh7e3tMnTp10DVTp04taD5EDO2sveXBBx+MNWvWxC9/+cu47LLLRnKbnAYKPWt//OMfY9++fbFgwYL+sb6+voiImDhxYuzZsycuuOCCkd0049JQ/l6bNm1anHHGGVFcXNw/dvHFF0dbW1v09PRESUnJiO6Z8WkoZ+2+++6LxYsXxy233BIREZdeeml0dXXFrbfeGvfcc08UFfn/5jl1J+qCsrKyk75LFPEuuFNUUlISc+bMiebm5v6xvr6+aG5ujqqqqkHXVFVVDZgfEfHCCy+ccD5EDO2sRUR8+9vfjgceeCCamppi7ty5o7FVxrlCz9pFF10UL7/8crS2tvY/PvWpT/V/k04+nx/N7TOODOXvtauuuipee+21/vCOiHj11Vdj2rRpgogTGspZe+ONN44Ln7di/N+foYdTN2xdUNh3QIyMDRs2ZLlcLnvqqaeyP/zhD9mtt96anXPOOVlbW1uWZVm2ePHibMWKFf3zf/vb32YTJ07MHnzwwWzXrl1ZfX19dsYZZ2Qvv/zyWL0ExolCz9qaNWuykpKSbNOmTdlf/vKX/sfRo0fH6iUwThR61v6bb5/jZBV61vbv35+dffbZ2Re/+MVsz5492c9//vNsypQp2Te+8Y2xegmME4Wetfr6+uzss8/OfvSjH2V79+7NfvGLX2QXXHBBdsMNN4zVS2AcOHr0aLZz585s586dWURkDz/8cLZz587sz3/+c5ZlWbZixYps8eLF/fP37t2bnXXWWdlXvvKVbNeuXVljY2NWXFycNTU1FXTdd0UUZVmWffe7383OO++8rKSkJJs3b172u9/9rv8/u/rqq7OlS5cOmP/jH/84u/DCC7OSkpLsIx/5SLZ58+ZR3jHjVSFn7X3ve18WEcc96uvrR3/jjDuF/r32f4kiClHoWXvppZeyysrKLJfLZeeff372zW9+Mzt27Ngo75rxqJCz9uabb2Zf+9rXsgsuuCArLS3N8vl89oUvfCH7+9//PvobZ9x48cUXB/3fXm+draVLl2ZXX331cWtmz56dlZSUZOeff372/e9/v+DrTsgy9y8BAIB0jflnigAAAMaSKAIAAJImigAAgKSJIgAAIGmiCAAASJooAgAAkiaKAACApIkiAAAgaaIIAABImigCAACSJooAAICkiSIAACBp/x9+2u0j/aOThAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "# Make figure of CSF/SC signal ratio from T2starw scan\n", "\n", @@ -372,10 +695,13 @@ "# Data storage for statistics\n", "data_stats = []\n", "\n", + "# Data storage for Plotly\n", + "t2_data_plotly = {}\n", + "\n", "# Iterate over each subject and create a subplot\n", "for i, subject in enumerate(subjects):\n", " ax = axes[i]\n", - "\n", + " t2_data_plotly[subject]={}\n", " for shim_mode in shim_modes:\n", " # Initialize list to collect data for this shim method\n", " method_data = []\n", @@ -406,14 +732,19 @@ "\n", " # If there's data for this shim method, plot it and compute stats\n", " if method_data:\n", + " t2_data_plotly[subject][shim_mode]=method_data\n", + "\n", " # Plotting each file's data separately\n", " for resampled_data in method_data:\n", " ax.plot(x_grid, resampled_data, label=f\"{shim_mode}\", linewidth=line_width)\n", + "\n", " \n", " # Compute stats on the non-resampled data (to avoid interpolation errors)\n", " mean_data = np.mean(data_sc_csf_ratio)\n", " sd_data = np.std(data_sc_csf_ratio)\n", " data_stats.append([subject, shim_mode, mean_data, sd_data])\n", + " else:\n", + " t2_data_plotly[subject][shim_mode]=None\n", "\n", " # Set x-ticks and labels for the bottom subplot\n", " if i == n_rows - 1: # Check if it's the last subplot\n", @@ -443,11 +774,22 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Perform statistics" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "adbc3046", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input", + "report_output" + ] + }, "outputs": [], "source": [ "# Perform statistics\n", @@ -515,124 +857,159 @@ }, { "cell_type": "markdown", - "id": "bda7c3ac", "metadata": {}, "source": [ - "## Process fmap/TFL (flip angle maps)" + "# 6     |     Process fmap/TFL (flip angle maps)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Register TFL flip angle maps to the GRE scan ⏳\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "e0cde907", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Register TFL flip angle maps to the GRE scan ⏳\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " !sct_register_multimodal -i {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -d ../anat/{subject}_acq-CoV_T2starw_crop.nii.gz -dseg ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", + " for shim_mode in shim_modes:\n", + " !sct_register_multimodal -i {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -d ../anat/{subject}_acq-CoV_T2starw_crop.nii.gz -dseg ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=slicereg,metric=CC -qc {path_qc}" ] }, { "cell_type": "markdown", - "id": "39de4623", "metadata": {}, "source": [ - "### Verify QC report (B1maps to GRE registration)\n", - "\n", - "Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis." + "Warping spinal cord segmentation and vertebral level to each flip angle map" ] }, { "cell_type": "code", "execution_count": null, - "id": "6a6d7670", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Warping spinal cord segmentation and vertebral level to each flip angle map\n", + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", + " for shim_mode in shim_modes:\n", + " # Use linear interpolation to preserve partial volume information (needed when extracting values along the cord)\n", + " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x linear -o {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz\n", + " # Use nearest neighbour (nn) interpolation because we are dealing with non-binary discrete labels\n", + " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert).\n", "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " # Use linear interpolation to preserve partial volume information (needed when extracting values along the cord)\n", - " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x linear -o {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz\n", - " # Use nearest neighbour (nn) interpolation because we are dealing with non-binary discrete labels\n", - " !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz" + "The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage, then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "a835cdee", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert)\n", - "# The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage,\n", - "# then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n", "\n", "GAMMA = 2.675e8; # [rad / (s T)]\n", "requested_fa = 90 # saturation flip angle -- hard-coded in sequence\n", "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " # Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n", - " with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n", - " metadata = json.load(f)\n", - " ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n", - " print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n", + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", + " for shim_mode in shim_modes:\n", + " # Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n", + " with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n", + " metadata = json.load(f)\n", + " ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n", + " print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n", "\n", - " # Open flip angle map with nibabel\n", - " nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n", - " acquired_fa = nii.get_fdata()\n", + " # Open flip angle map with nibabel\n", + " nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n", + " acquired_fa = nii.get_fdata()\n", "\n", - " # Siemens maps are in units of flip angle * 10 (in degrees)\n", - " acquired_fa = acquired_fa / 10\n", + " # Siemens maps are in units of flip angle * 10 (in degrees)\n", + " acquired_fa = acquired_fa / 10\n", "\n", - " # Account for the power loss between the coil and the socket. That number was given by Siemens.\n", - " voltage_at_socket = ref_voltage * 10 ** -0.095\n", + " # Account for the power loss between the coil and the socket. That number was given by Siemens.\n", + " voltage_at_socket = ref_voltage * 10 ** -0.095\n", "\n", - " # Compute B1 map in [T/V]\n", - " b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n", + " # Compute B1 map in [T/V]\n", + " b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n", "\n", - " # Convert to [nT/V]\n", - " b1_map = b1_map * 1e9\n", + " # Convert to [nT/V]\n", + " b1_map = b1_map * 1e9\n", "\n", - " # Save as NIfTI file\n", - " nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n", - " nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")" + " # Save as NIfTI file\n", + " nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n", + " nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "740fda91", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n", - "\n", - "for subject in subjects:\n", - " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - " for shim_mode in shim_modes:\n", - " fname_result_b1plus = os.path.join(path_results, f\"{subject}_TB1map_{shim_mode}.csv\")\n", - " !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o {fname_result_b1plus}" + "if notebook != 'neurolibre-figures':\n", + " for subject in subjects:\n", + " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", + " for shim_mode in shim_modes:\n", + " fname_result_b1plus = os.path.join(path_results, f\"{subject}_TB1map_{shim_mode}.csv\")\n", + " !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o {fname_result_b1plus}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make figure of B1+ values along the spinal cord across shim methods" ] }, { "cell_type": "code", "execution_count": null, - "id": "b0185bb0", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input", + "remove_output" + ] + }, "outputs": [], "source": [ - "# Make figure of B1+ values along the spinal cord across shim methods\n", - "\n", "# Go back to root data folder\n", "os.chdir(os.path.join(path_data))\n", "\n", @@ -674,12 +1051,17 @@ "# Data storage for statistics\n", "data_stats = []\n", "\n", + "# Data storage for Plotly\n", + "b1_data_plotly = {}\n", + "\n", "# Iterate over each subject and create a subplot\n", "for i, subject in enumerate(subjects):\n", + " \n", " ax = axes[i]\n", " \n", " os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", - "\n", + " b1_data_plotly[subject]={}\n", + " \n", " for shim_mode in shim_modes:\n", " # Initialize list to collect data for this shim method\n", " method_data = []\n", @@ -704,12 +1086,16 @@ " if method_data:\n", " # Plotting each file's data separately\n", " for resampled_data in method_data:\n", + " b1_data_plotly[subject][shim_mode]=resampled_data\n", " ax.plot(x_grid, resampled_data, label=f\"{shim_mode}\", linewidth=line_width)\n", "\n", " # Compute stats on the non-resampled data (to avoid interpolation errors)\n", " mean_data = np.mean(wa_data)\n", " sd_data = np.std(wa_data)\n", " data_stats.append([subject, shim_mode, mean_data, sd_data])\n", + " else:\n", + " b1_data_plotly[subject][shim_mode]=None\n", + "\n", "\n", " # Set x-ticks and labels for the bottom subplot\n", " if i == n_rows - 1: # Check if it's the last subplot\n", @@ -739,15 +1125,24 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Perform statistics" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "c09f3c05", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input", + "report_output" + ] + }, "outputs": [], "source": [ - "# Perform statistics\n", - "\n", "# Convert to DataFrame and save to CSV\n", "df_stats = pd.DataFrame(data_stats, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n", "df_stats.to_csv(os.path.join(path_results, 'stats_b1plus.csv'), index=False)\n", @@ -809,40 +1204,66 @@ "print(f\"Paired t-tests:\\n{comparison_results}\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 7     |     Paper figures" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Prepare RF shimming mask for figure\n" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "1dde4c00", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input" + ] + }, "outputs": [], "source": [ - "# Prepare RF shimming mask for figure\n", - "\n", "# Select subject to show\n", "subject = 'sub-05'\n", "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", "\n", "# Create RF shimming mask from the segmentation (re-create the live RF shimming procedure)\n", "file_mask = f\"{subject}_acq-anatCP_TB1TFL_mask-shimming.nii.gz\"\n", - "!sct_maths -i {subject}_acq-anatCP_TB1TFL_seg_labeled.nii.gz -thr 3 -uthr 9 -o {file_mask}\n", - "!sct_maths -i {file_mask} -bin 1 -o {file_mask}\n", - "!sct_create_mask -i {subject}_acq-anatCP_TB1TFL.nii.gz -p centerline,{file_mask} -size 28mm -f cylinder -o {file_mask}\n", + "if notebook != 'figures':\n", + "\n", + " !sct_maths -i {subject}_acq-anatCP_TB1TFL_seg_labeled.nii.gz -thr 3 -uthr 9 -o {file_mask}\n", + " !sct_maths -i {file_mask} -bin 1 -o {file_mask}\n", + " !sct_create_mask -i {subject}_acq-anatCP_TB1TFL.nii.gz -p centerline,{file_mask} -size 28mm -f cylinder -o {file_mask}\n", "\n", - "# Dilate mask\n", - "!sct_maths -i {file_mask} -dilate 2 -shape disk -dim 2 -o {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz\n", - "# Subtract dilated mask with mask to get the edge\n", - "!sct_maths -i {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz -sub {file_mask} -o {file_mask}" + " # Dilate mask\n", + " !sct_maths -i {file_mask} -dilate 2 -shape disk -dim 2 -o {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz\n", + " # Subtract dilated mask with mask to get the edge\n", + " !sct_maths -i {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz -sub {file_mask} -o {file_mask}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create figure of B1+ maps" ] }, { "cell_type": "code", "execution_count": null, - "id": "93d72ec3", - "metadata": {}, + "metadata": { + "tags": [ + "hide_input", + "report_output" + ] + }, "outputs": [], "source": [ - "# Create figure of B1+ maps\n", - "\n", "# Select subject to show\n", "subject = 'sub-05'\n", "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n", @@ -937,14 +1358,340 @@ "# cbar = plt.colorbar(im, cax=cbar_ax)\n", "cbar_ax.set_title('nT/V', size=12)\n", "plt.savefig(os.path.join(path_results, 'fig_b1plus_map.png'), dpi=300, format='png')\n", - "plt.show\n" + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Figure 2. B1+ efficiency for one participant (sub-05) across all seven RF shimming conditions. The top left panel shows the tfl_b1map magnitude image with an overlay of the mask that was used to perform RF shimming. Text inserts show the mean (in nT/V) and CoV (in %) of B1+ efficiency along the spinal cord between C3 and T2. " ] }, { "cell_type": "code", "execution_count": null, - "id": "9128d799", + "metadata": { + "tags": [ + "hide_input", + "remove_output" + ] + }, + "outputs": [], + "source": [ + "# PYTHON CODE\n", + "# Module imports\n", + "\n", + "# Base python\n", + "import os\n", + "from os import path\n", + "from pathlib import Path\n", + "\n", + "# Graphical\n", + "\n", + "import plotly.graph_objs as go\n", + "from IPython.display import display, HTML\n", + "from plotly import __version__\n", + "from plotly.offline import init_notebook_mode, iplot, plot\n", + "config={\n", + " 'showLink': False,\n", + " 'displayModeBar': False,\n", + " 'toImageButtonOptions': {\n", + " 'format': 'png', # one of png, svg, jpeg, webp\n", + " 'filename': 'custom_image',\n", + " 'height': 2500,\n", + " 'width': 500,\n", + " 'scale': 2 # Multiply title/legend/axis/canvas sizes by this factor\n", + " }\n", + " }\n", + "\n", + "init_notebook_mode(connected=True)\n", + "\n", + "from plotly.subplots import make_subplots\n", + "import plotly.graph_objects as go\n", + "\n", + "import seaborn as sns\n", + "\n", + "# Set the color palette\n", + "pal=sns.color_palette()\n", + "\n", + "# Imports\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "## Setup for plots\n", + "fig = make_subplots(rows=5, cols=2, vertical_spacing = 0.025,\n", + " subplot_titles=(\n", + " 'sub-01',\n", + " 'sub-01',\n", + " 'sub-02',\n", + " 'sub-02',\n", + " 'sub-03',\n", + " 'sub-03',\n", + " 'sub-04',\n", + " 'sub-04',\n", + " 'sub-05',\n", + " 'sub-05',))\n", + "\n", + "t2_datasets={}\n", + "b1_datasets={}\n", + "\n", + "t2_data = []\n", + "b1_data = []\n", + "\n", + "legend_bool = True\n", + "for subject in subjects:\n", + " index = 0\n", + " t2_datasets[subject]={}\n", + " b1_datasets[subject]={}\n", + "\n", + " for shim_mode in shim_modes:\n", + " t2_datasets[subject][shim_mode]={}\n", + " b1_datasets[subject][shim_mode]={}\n", + "\n", + " t2_data=go.Line(\n", + " x=x_grid,\n", + " y=t2_data_plotly[subject][shim_mode][0],\n", + " name=shim_mode,\n", + " legendgroup=shim_mode,\n", + " line=dict(color='rgb'+str(pal[index]), width=3),\n", + " showlegend=False\n", + " )\n", + "\n", + " b1_data=go.Line(\n", + " x=x_grid,\n", + " y=b1_data_plotly[subject][shim_mode],\n", + " name=shim_mode,\n", + " legendgroup=shim_mode,\n", + " line=dict(color='rgb'+str(pal[index]), width=3),\n", + " showlegend=legend_bool\n", + " ) \n", + "\n", + " \n", + " t2_datasets[subject][shim_mode]=t2_data\n", + " b1_datasets[subject][shim_mode]=b1_data\n", + "\n", + " index += 1\n", + " legend_bool=False\n", + " \n", + "\n", + "index = 1\n", + "# For z-ordering \n", + "for subject in subjects:\n", + " for shim_mode in shim_modes:\n", + " fig.add_trace(\n", + " t2_datasets[subject][shim_mode],\n", + " row=index, col=1\n", + " )\n", + " fig.add_trace(\n", + " b1_datasets[subject][shim_mode],\n", + " row=index, col=2\n", + " )\n", + " index+=1\n", + "\n", + "\n", + "index = 1\n", + "for subject in subjects:\n", + " if index == 5:\n", + " x_title = 'Vertebral Levels'\n", + " showticklabels=True\n", + " else:\n", + " x_title = None\n", + " showticklabels=False\n", + "\n", + " fig.update_xaxes(\n", + " type=\"linear\",\n", + " autorange=True,\n", + " title=x_title,\n", + " showgrid=True,\n", + " gridcolor='rgb(169,169,169)',\n", + " tickvals=label_positions,\n", + " ticktext=vertebral_levels,\n", + " showticklabels=showticklabels,\n", + " title_font_family=\"Times New Roman\",\n", + " title_font_size = 20,\n", + " linecolor='black',\n", + " linewidth=2,\n", + " tickfont=dict(\n", + " family='Times New Roman',\n", + " size=16,\n", + " ),\n", + " row=index, col=1\n", + " )\n", + " if index == 1:\n", + " fig.update_yaxes(\n", + " type=\"linear\",\n", + " title={\n", + " 'text':'CSF/Cord T2*w signal ratio',\n", + " 'standoff':0\n", + " },\n", + " range=[1.05, 1.4],\n", + " showgrid=True,\n", + " gridcolor='rgb(169,169,169)',\n", + " title_font_family=\"Times New Roman\",\n", + " title_font_size = 20,\n", + " linecolor='black',\n", + " linewidth=2,\n", + " tickfont=dict(\n", + " family='Times New Roman',\n", + " size=16,\n", + " ),\n", + " row=index, col=1\n", + " )\n", + " else:\n", + " fig.update_yaxes(\n", + " type=\"linear\",\n", + " title={\n", + " 'text':'CSF/Cord T2*w signal ratio',\n", + " 'standoff':0\n", + " },\n", + " showgrid=True,\n", + " gridcolor='rgb(169,169,169)',\n", + " title_font_family=\"Times New Roman\",\n", + " title_font_size = 20,\n", + " linecolor='black',\n", + " linewidth=2,\n", + " tickfont=dict(\n", + " family='Times New Roman',\n", + " size=16,\n", + " ),\n", + " row=index, col=1\n", + " )\n", + "\n", + " fig.update_xaxes(\n", + " type=\"linear\",\n", + " autorange=True,\n", + " title=x_title,\n", + " showgrid=True,\n", + " gridcolor='rgb(169,169,169)',\n", + " tickvals=label_positions,\n", + " ticktext=vertebral_levels,\n", + " showticklabels=showticklabels,\n", + " title_font_family=\"Times New Roman\",\n", + " title_font_size = 20,\n", + " linecolor='black',\n", + " linewidth=2,\n", + " tickfont=dict(\n", + " family='Times New Roman',\n", + " size=16,\n", + " ),\n", + " row=index, col=2\n", + " )\n", + " fig.update_yaxes(\n", + " type=\"linear\",\n", + " title={\n", + " 'text':'B1+ efficiency [nT/V]',\n", + " 'standoff':0\n", + " },\n", + " showgrid=True,\n", + " gridcolor='rgb(169,169,169)',\n", + " title_font_family=\"Times New Roman\",\n", + " title_font_size = 20,\n", + " linecolor='black',\n", + " linewidth=2,\n", + " tickfont=dict(\n", + " family='Times New Roman',\n", + " size=16,\n", + " ),\n", + " row=index, col=2\n", + " )\n", + "\n", + " index+=1\n", + "\n", + "fig.update_layout(height=1800, width=900)\n", + "\n", + "for i in fig['layout']['annotations']:\n", + " i['font'] = dict(size=22)\n", + "\n", + "fig.update_layout(legend=dict(\n", + " yanchor=\"top\",\n", + " y=0.999,\n", + " xanchor=\"left\",\n", + " x=0.01,\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=12\n", + " ),\n", + " bordercolor=\"Black\",\n", + " borderwidth=1.5\n", + " ),\n", + " legend_tracegroupgap=0,\n", + " paper_bgcolor='rgb(255, 255, 255)',\n", + " plot_bgcolor='rgb(255, 255, 255)',\n", + ")\n", + "\n", + "fig.add_annotation(\n", + " dict(\n", + " x=-0.07,\n", + " y=-0.04,\n", + " showarrow=False,\n", + " text='A',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=48\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ))\n", + "fig.add_annotation(\n", + " dict(\n", + " x=0.5,\n", + " y=-0.04,\n", + " showarrow=False,\n", + " text='B',\n", + " font=dict(\n", + " family='Times New Roman',\n", + " size=48\n", + " ),\n", + " xref='paper',\n", + " yref='paper'\n", + " ),\n", + ")\n", + "#iplot(fig, filename = 'figure4a', config = config)\n", + "plot(fig, filename = 'figure1.html', config = config)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "remove_input", + "report_output" + ] + }, + "outputs": [], + "source": [ + "display(HTML('figure1.html'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Figure 3. B1+ efficiency (A) and CSF/Cord signal ratio from the GRE scan (B) across subjects and across different RF shimming conditions. Data were measured in the spinal cord from C3 to T2 vertebral levels. To match the x-ticks across subjects, the C2-C3 and the T2-T3 intervertebral discs of each subject were aligned with that of the PAM50 template {cite:p}`DELEENER2018170`, and the curves were linearly scaled.\n" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "# References \n", + "\n", + "```{bibliography}\n", + ":filter: docname in docnames\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "remove_input", + "remove_output" + ] + }, "outputs": [], "source": [ "# Indicate duration of data processing\n", @@ -966,6 +1713,7 @@ } ], "metadata": { + "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", @@ -981,9 +1729,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.8" + "version": "3.9.13" } }, "nbformat": 4, - "nbformat_minor": 5 + "nbformat_minor": 2 }