density_map_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6758a718",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pydicom\n",
    "from pydicom import dcmread\n",
    "\n",
    "import pandas as pd\n",
    "from pandas import DataFrame\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from numba import jit\n",
    "import cv2\n",
    "\n",
    "import os\n",
    "\n",
    "import scipy as st\n",
    "from scipy.constants import e, h, c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2e920321",
   "metadata": {},
   "outputs": [],
   "source": [
    "def import_dicom(in_dir_path):\n",
    "    #names the path to the input model directory, and lists the files in this directory\n",
    "\n",
    "    in_dir_list = sorted(os.listdir(in_dir_path))\n",
    "    model_depth = len(in_dir_list)\n",
    "    model_data = []\n",
    "\n",
    "\n",
    "    i = 0\n",
    "    while i < model_depth:\n",
    "        path = in_dir_path + \"/\" + in_dir_list[i]\n",
    "        current_slice = pydicom.dcmread(path)\n",
    "        model_slice = current_slice.pixel_array\n",
    "        model_slice = np.asarray(model_slice)\n",
    "        model_data.append(model_slice)\n",
    "        i = i + 1\n",
    "    model_data = np.asarray(model_data, dtype = np.int16)\n",
    "    \n",
    "        \n",
    "    return model_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1a032728",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_model_kVp(in_dir_path):   \n",
    "    in_dir_list = sorted(os.listdir(in_dir_path))\n",
    "    path = in_dir_path + \"/\" + in_dir_list[0]\n",
    "    first_slice = pydicom.dcmread(path)\n",
    "    model_kVp = first_slice[0x0018, 0x0060].value\n",
    "    return model_kVp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d4a23c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_known_KeV(mat_dir_path):    \n",
    "    mat_dir_list = os.listdir(mat_dir_path)\n",
    "    first_file_path = mat_dir_path + \"/\" + mat_dir_list[0]\n",
    "    first_file = pd.read_csv(first_file_path)    \n",
    "    known_KeV = first_file[first_file.columns[0]].values.tolist()\n",
    "    known_KeV = np.asarray(known_KeV)\n",
    "    return known_KeV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b7c2500b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_known_CT_numbers(mat_dir_path):\n",
    "    mat_dir_path = mat_dir_path + \"/CT_numbers\"\n",
    "    mat_dir_list = sorted(os.listdir(mat_dir_path))\n",
    "            \n",
    "    tissue_types = []\n",
    "    for x in mat_dir_list:\n",
    "        tissue_types.append(x.replace(\".csv\", \"\"))\n",
    "        #x.replace() removes the file type, to leave just the tissue name\n",
    "       \n",
    "    tissue_data = []\n",
    "    for x in mat_dir_list:\n",
    "        mat_path = mat_dir_path + \"/\" + x\n",
    "        curent_file = pd.read_csv(mat_path)\n",
    "        current_tissue_data = curent_file[curent_file.columns[1]].values.tolist()\n",
    "        current_tissue_data = np.asarray(current_tissue_data)\n",
    "        tissue_data.append(current_tissue_data)\n",
    "    \n",
    "    tissue_types = np.asarray(tissue_types)\n",
    "    tissue_data = np.asarray(tissue_data)\n",
    "    \n",
    "    known_CT_numbers=[tissue_types, tissue_data]  \n",
    "    \n",
    "    return known_CT_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "be80c187",
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_CT_numbers(target_KeV, known_CT_numbers, known_KeV):\n",
    "    \n",
    "    #finds the closest two known KeVs\n",
    "    \n",
    "    #calculates the index of the nearest value more than the kVp\n",
    "    upper_i = 0\n",
    "        \n",
    "    while known_KeV[upper_i] < target_KeV:\n",
    "            upper_i = upper_i + 1\n",
    "\n",
    "\n",
    "    #calculates the index of the nearest value less than the kVp\n",
    "    lower_i = 0      \n",
    "    if target_KeV >= known_KeV[0]:\n",
    "        while known_KeV[lower_i] > target_KeV:\n",
    "            lower_i = lower_i - 1\n",
    "\n",
    "        \n",
    "        \n",
    "    upper_KeV = known_KeV[upper_i]\n",
    "    lower_KeV = known_KeV[lower_i]\n",
    "        \n",
    "    if (upper_KeV == target_KeV):\n",
    "        upper_KeV_weight = 1\n",
    "        lower_KeV_weight = 0\n",
    "    elif (lower_KeV == target_KeV):\n",
    "        upper_KeV_weight = 0\n",
    "        lower_KeV_weight = 1\n",
    "    else:\n",
    "        #solution to the simultaneous equations\n",
    "        #target_KeV = upper_KeV_weight*upper_KeV + lower_KeV_weight*lower_KeV\n",
    "        #AND\n",
    "        #upper_KeV_weight + lower_KeV_weight = 1\n",
    "        #which express the KeV as a weighted average of two known values\n",
    "        upper_KeV_weight = (target_KeV - lower_KeV)/(upper_KeV-lower_KeV + 0.00000000001)\n",
    "        #adding a very small number to to denominator to avoid ever dividing by 0\n",
    "        lower_KeV_weight = 1- upper_KeV_weight\n",
    "    \n",
    "    #finds the expected CT number of each material at the keV that the scan was taken at\n",
    "    predicted_CT_numbers = []\n",
    "  \n",
    "    material_ids = known_CT_numbers[0]\n",
    "    \n",
    "    material_number = []\n",
    "    for x in known_CT_numbers[1]:\n",
    "        value_to_append = upper_KeV_weight * x[upper_i] + lower_KeV_weight * x[lower_i]\n",
    "        material_number.append(value_to_append)\n",
    "        \n",
    "    \n",
    "    material_number_sorted = sorted(material_number)\n",
    "    sorting_indexs = np.argsort(material_number)\n",
    "\n",
    "    material_ids_sorted = [0] * len(material_ids)\n",
    "    \n",
    "    i = 0\n",
    "    for x in sorting_indexs:\n",
    "        material_ids_sorted[i] = material_ids[x]\n",
    "        i = i + 1\n",
    "        \n",
    "    material_ids_sorted = np.asarray(material_ids_sorted)\n",
    "    material_number_sorted = np.asarray(material_number_sorted)\n",
    "    \n",
    "    \n",
    "\n",
    "    predicted_CT_numbers = [material_ids_sorted, material_number_sorted]\n",
    "    return predicted_CT_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e8e91cda",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_filtered_spectrum(kVp, interval, anode_atomic_number, filter_thickness):\n",
    "    \n",
    "    K = anode_atomic_number\n",
    "    c = st.constants.c\n",
    "    h_joules = st.constants.h\n",
    "    \n",
    "    bins = np.arange(0, kVp, interval)  \n",
    "    unfiltered_binned_spectrum = np.zeros((2, len(bins)))\n",
    "    e = st.constants.elementary_charge\n",
    "    h = 1000 * e * h_joules\n",
    "    \n",
    "    total_integral = ((-K * np.square(kVp)) * (-kVp)) / (6 * np.square(c * h))\n",
    "    \n",
    "    i = 0\n",
    "    while i < len(bins):\n",
    "        E = bins[i]\n",
    "        start = ((-K * np.square(E)) * (2*E - 3*kVp)) / (6 * np.square(c * h))        \n",
    "        end = ((-K * np.square(E + interval)) * (2*(E + interval) - 3*kVp)) / (6 * np.square(c * h))\n",
    "        integral = end - start\n",
    "        integral_proportion = integral/total_integral\n",
    "        unfiltered_binned_spectrum[0][i] = E + 0.5 * interval\n",
    "        unfiltered_binned_spectrum[1][i] = integral_proportion\n",
    "        i = i + 1\n",
    "        \n",
    "    Al_file = pd.read_csv(\"./material_info/Al_filter.csv\")\n",
    "    \n",
    "    Al_lin_atten = Al_file[Al_file.columns[1]].values.tolist()\n",
    "    Al_lin_atten = np.asarray(Al_lin_atten)\n",
    "    \n",
    "    Al_known_KeV = Al_file[Al_file.columns[0]].values.tolist()\n",
    "    Al_known_KeV = np.asarray(Al_known_KeV)\n",
    "    \n",
    "    Al_filter = np.zeros((2, len(unfiltered_binned_spectrum[0])))\n",
    "    \n",
    "  \n",
    "    \n",
    "    Al_filter[0] = unfiltered_binned_spectrum[0]\n",
    "   \n",
    "    \n",
    "    i = 0\n",
    "    while i < len(Al_filter[0]):\n",
    "        dif = abs(Al_filter[0][i] - Al_known_KeV)\n",
    "        closest_index = np.argmin(dif)\n",
    "        \n",
    "        if Al_known_KeV[closest_index] == Al_filter[0][i]:\n",
    "            upper_KeV = Al_known_KeV[closest_index]\n",
    "            lower_KeV = Al_known_KeV[closest_index]\n",
    "            upper_weight = 1\n",
    "            lower_weight = 0\n",
    "            upper_u =  Al_lin_atten[closest_index]\n",
    "            lower_u = Al_lin_atten[closest_index]\n",
    "\n",
    "        elif Al_known_KeV[closest_index] < Al_filter[0][i]:\n",
    "            upper_KeV = Al_known_KeV[closest_index + 1]\n",
    "            lower_KeV = Al_known_KeV[closest_index]\n",
    "            upper_weight = (Al_filter[0][i] - lower_KeV)/(upper_KeV - lower_KeV)\n",
    "            lower_weight = 1 - upper_weight\n",
    "            upper_u =  Al_lin_atten[closest_index + 1]\n",
    "            lower_u = Al_lin_atten[closest_index]\n",
    "            \n",
    "        else:\n",
    "            upper_KeV = Al_known_KeV[closest_index]\n",
    "            lower_KeV = Al_known_KeV[closest_index - 1]\n",
    "            upper_weight = (Al_filter[0][i] - lower_KeV)/(upper_KeV - lower_KeV)\n",
    "            lower_weight = 1 - upper_weight\n",
    "            upper_u =  Al_lin_atten[closest_index]\n",
    "            lower_u = Al_lin_atten[closest_index - 1]\n",
    "            \n",
    "        Al_filter[1][i] = upper_weight*upper_u + lower_weight*lower_u\n",
    "        i = i + 1\n",
    "        \n",
    "        \n",
    "        \n",
    "    unfiltered_spectrum = unfiltered_binned_spectrum[1]\n",
    "    product = -Al_filter[1] * filter_thickness\n",
    "\n",
    "    filter_ratio = np.exp(product)\n",
    "    \n",
    "    filtered_spectrum_values = filter_ratio * unfiltered_spectrum\n",
    "    \n",
    "    normal_sum = np.sum(filtered_spectrum_values)\n",
    "    \n",
    "    normal_filtered_spectrum_values = filtered_spectrum_values/normal_sum\n",
    "    \n",
    "    normal_filtered_spectrum = np.zeros((2, len(normal_filtered_spectrum_values)))\n",
    "    \n",
    "    normal_filtered_spectrum[0] = unfiltered_binned_spectrum[0]\n",
    "    normal_filtered_spectrum[1] = normal_filtered_spectrum_values\n",
    "    \n",
    "    return(normal_filtered_spectrum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "392e2697",
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_CT_numbers_poly(mat_dir_path, kVp, interval, anode_atomic_number, filter_thickness):\n",
    "    #mat_dir_path is the path to material_info file, or equivalent\n",
    "    #interval refers to size of the invervals used when calculating x-ray spectra. Smaller interval is faster and lest accurate.\n",
    "    #filter_thickness is in centimetres\n",
    "    known_KeV = find_known_KeV(mat_dir_path)\n",
    "    known_CT_numbers = find_known_CT_numbers(mat_dir_path)\n",
    "    filtered_spectrum = generate_filtered_spectrum(kVp, interval, anode_atomic_number, filter_thickness)\n",
    "    \n",
    "    poly_CT_numbers_unsum = np.zeros((len(filtered_spectrum[0]), len(known_CT_numbers[0])))\n",
    "    \n",
    "    i = 0\n",
    "    while i < len(filtered_spectrum[0]):\n",
    "        line = predict_CT_numbers(filtered_spectrum[0][i], known_CT_numbers, known_KeV)\n",
    "        line_val = line[1]\n",
    "        line_val  = line_val  * filtered_spectrum[1][i]\n",
    "        poly_CT_numbers_unsum[i] = line_val\n",
    "   \n",
    "        i = i + 1\n",
    "    \n",
    "    poly_CT_numbers_val = np.sum(poly_CT_numbers_unsum, axis = 0)\n",
    "    \n",
    "    poly_CT_numbers = [line[0], poly_CT_numbers_val]\n",
    "    \n",
    "    return poly_CT_numbers\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "df9cda0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def identify_tissue(val, expected_CT_numbers):\n",
    "    #indentifies the closest match for a given CT number   \n",
    "    dif = abs(expected_CT_numbers - val)\n",
    "    tissue_index = np.argmin(dif)\n",
    "    tissue_type = expected_CT_numbers[0][tissue_index]\n",
    "    return tissue_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5d28235d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mass_density_dictionary(mat_dir_path):\n",
    "    #creates a dictionary of tissue types and their mass densities\n",
    "    \n",
    "    path = mat_dir_path + \"/mass_densities.csv\"\n",
    "    density_file = pd.read_csv(path)\n",
    "    density_file_material_list = density_file[density_file.columns[0]].values.tolist()\n",
    "    density_file_value_list = density_file[density_file.columns[1]].values.tolist()\n",
    "    mass_densities = {}\n",
    "\n",
    "    i = 0\n",
    "    while i < (len(density_file_material_list)):\n",
    "        mass_densities.update({density_file_material_list[i] : density_file_value_list[i]})\n",
    "        i = i + 1\n",
    "        \n",
    "    return mass_densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8273ce44",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_tissue_ranges(predicted_CT_numbers):\n",
    "    \n",
    "    #creates two lists of the predicted CT numbers, offset by one index, for efficient computation using vectorisation\n",
    "    #the principle is the cut-off between values is the mid point between them, which is also the average\n",
    "    x = predicted_CT_numbers[1]\n",
    "    y = predicted_CT_numbers[1]\n",
    "    \n",
    "    #offseting by one index\n",
    "    x = np.insert(x, 0, x[0])\n",
    "    y = np.insert(y, len(y)-1, y[-1])\n",
    "    \n",
    "    #finding average in vectorised way\n",
    "    z = x + y  \n",
    "    z = z/2\n",
    "    \n",
    "    material_count = len(predicted_CT_numbers[1])\n",
    "    tissue_ranges = np.zeros((material_count, 2))\n",
    "    \n",
    "    i = 0\n",
    "    while i < material_count:\n",
    "        tissue_ranges[i][0] = z[i]\n",
    "        tissue_ranges[i][1] = z[i + 1]\n",
    "        i = i + 1\n",
    "    \n",
    "    \n",
    "    #expands the top and bottom range\n",
    "    tissue_ranges[0][0] = tissue_ranges[0][0] * 2\n",
    "    tissue_ranges[-1][1] = tissue_ranges[-1][1] * 2\n",
    " \n",
    "\n",
    "    \n",
    "\n",
    "    return tissue_ranges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e5fa861d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_density_map(model, tissue_ranges, predicted_CT_numbers, mass_densities):\n",
    "    seg_density_map = np.zeros((len(tissue_ranges), len(model), len(model[0]), len(model[0][0])))\n",
    "    \n",
    "    i = 0\n",
    "    while i < len(tissue_ranges):        \n",
    "        lower = tissue_ranges[i][0]\n",
    "        upper = tissue_ranges[i][1]\n",
    "\n",
    "        mask = np.where((model > lower) & (model < upper), 1, 0)\n",
    "        \n",
    "        tissue_type = predicted_CT_numbers[0][i]\n",
    "        tissue_density = mass_densities[tissue_type]\n",
    "        \n",
    "        seg_density_map[i] = tissue_density * mask\n",
    "\n",
    "        i = i + 1\n",
    "        \n",
    "    density_map = np.sum(seg_density_map, axis = 0)\n",
    "\n",
    "    \n",
    "    return density_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "ebe0b9b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x28b30a66370>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABdmElEQVR4nO1dd5wU5f1+3qm7exy9KEVABBUsWNEoNkSKBaMhYkUskMRYEqPR+DNGE6MpGo2JMdYYJRILKvZCRLFgQyyIBQURBUGpd7dlyvv74533nXd2Zvd293bvjmMeP364nZ3y7sy83/dbny+hlCJGjBgxZChtPYAYMWK0P8SCIUaMGCHEgiFGjBghxIIhRowYIcSCIUaMGCHEgiFGjBgh1EwwEELGE0I+JoQsJYRcUqvrxIgRo/ogtchjIISoAD4BMBbASgBvAjiRUvph1S8WI0aMqqNWGsO+AJZSSj+nlOYAzAIwqUbXihEjRpWh1ei8/QB8KX1eCWBUoZ0JIXH6ZYwYtce3lNJepexYK8FAIrYFJj8hZDqA6TW6fowYMcL4otQdayUYVgIYIH3uD+BreQdK6a0AbgVijSFGjPaGWvkY3gQwlBAymBBiAJgCYE6NrhUjRowqoyYaA6XUJoT8FMAzAFQAd1JKF9fiWjFixKg+ahKuLHsQsSkRI0Zr4G1K6d6l7BhnPsaIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFgwxIgRI4RYMMSIESOEWDDEiBEjhFp1u26X0DQDqqI2u5/jOrBtC6aRENtc6kIhCihYNz1KKRSiIGdlAQCGbsKyc6CUgrf9MwwTBAS2Y8FxHKiqCsexAQCqqkFTdQCAZeegaToISGAclp0DIUTs548tF/g9LnVheePgMIxE4HwudQEAhJDAGNm5CDRNB6UuNFWH4zpQFRWO60BRFLiuC9u2oOsGFKKIsbmuAwBQFBWapov9+G8HEBiD7VggRBH78HPL583mMuxeEALLykLTDFDqglIKVdVACDtfLpcFwLexc7L7G7zPAGAaSVBQ5HIZAAS6boAQAgIifoeum6CUwrZzUFUNiqKK50upC00zxP3j51YUFYSQwH5sH0X8PvnZ8efHj8nfriqqGI+iqOL+cui6CcvKAah9R8etSjAAgKazB0wpFS8Ze6FzMDxBoLgum5DevgCQyTRCMwy4rgPXdWAYCbiuC5rLwDQS0HQDLnWh6yay2SbvpVFg2zk4jgNN06EoChzHAUCheAKKvQguFEWDoihiTJRSJjBUHURRxXduLiPGRKkLEA0ERLxIhCjQNN2bUEwAUNfxXigG3ftdum6y/VwHOSvLhJVugDhMIBBXhWVlxMTXNAO5XBq6noCum8jlMmIyaJrhCRkKQhSoqi6uDUKgqhpsx4JtW54QcaBrBrK5NBOymoZstgkAFffYHy97Lo50fDLZCbY3ibigNIyk95tdOI4NRVGZQNF0MZkJIeLeaxqbtFnvnrquA1VlU8Kyckil6mF4z0hVNWSzaXGOVKoeuVwGlpVDMtkJBvwFRdN0WFYOhChQvPM5jgVV1T0BQpFK1Ytnzd85WSC4rht4bxVFharqgedYS2x1giEK8stCKYXr2DCMBBzH8VYLNinYCudAIWylY6suBQiB47DJp6ps5dQ0JkQCGoKmQ1HYBHBdB5aTQyJRB9e14TgWMhn2WVEUcQwFRS7bBEM3oemG93Iw7cBxHDhOE3TdhGEkkc02wTRTcF2bCRYrC1034VAXiUQdEyyODdVbkbPZNFtxvd/OhR2boOwzF6CWlRMTznEaoKqa+J8LB/7SUuoinW4AAHHPTCMJzbs3qqLCti043uRnk8mAqmpCKPiaQA6OY7Oxe/s7DoGq6t7q7QCeVpJON0BV2W8xjaRYxdnx/ioOSpHNMYFHFBWqqop7xcfDtUXHdZjg8vbTNAPZbFrSZigymUZQ6jItzXuXuPCzrKy4r6oKITTyoSiKp7FY3pagVmAYSfFetAa2Kh+DrF7K0poQIlYPvk82m4bjWJ6K7ZsPmqZD0w1ks03I5dLgLyVf7WzPnCBEgaKoMM0UTCMJQhSxmnHzgE0aG7ZteZOdfZ/LZZDNNokVglJXTN78cfMXXlXZtZhWYiOXS3uTlnqrkA3bzjGNgBBf7eerOgDLyorx8wml6yYURYNhmEKd5veJvcxMS3Ic2xMQREwQANA1AwBBNpeB7VhIJOpgO5bQWjgymUZPvfa1OL6CclPC+9VwXSZ4mGnAtvNJw+898SanYSQASr1nxcbNBVI6vRkAvPvGhIPFhYm3avNrE0Kg6wmxgPAJzAUjWzjYOpvLZZiJ52kQfL9stgkKUYSQ5uPlaG7is3ezdRrDb1UaQ/7LKINrAHyF4mooExqGmCiapsO2ckINVhTVe0EhVjwOTdOZJmFnYGo6crmMUHspKEwz5a2swYdNCIFt21BVW7x4WUlrMIwEcrkMiKeiW1YOiqLCtnPeamV4gkUTY6OUwjBMz2xRvGNVZD0BoiiKJ5yYDU4pU+lVRQUFm/y2px7rugHbtoT5w1ZDjQkRonirLBMSjuvb/Px+snFlmQ2uGcIMsTxTRNdNEEKQy2Wg66aw+X3ThQhNBN7dZGaEK34z96Woqi6EgqYZ3gTngoyp+NRl98M0UyBEEfv7/gTm73Bdx9NsVCQSdQDY8ez+ukLou64DxzM3bdsSv5mbKbIJyxYfX1jzd0FVtcBClsk0en4bEnpfaoGtSmPId+4BvuRnqnyj93CYBmHoprD1KKVwHBu2lWMOR09l5A+bO8cAwDSToJKd6Loumpo2w/FMFK6iM0cXW4GYqWF4fgLuYMsINZtSKq1iFGwCew4+Qpg2II5XxHXZBGEvXC6XQSbTCNvKIZNp9PYjniZkiEnPzwlwB6inlXj3hE9w1+VOOPa9phnCOUqIAtNgvghdT0BVNZhmCpRSMXFNM+WZFb4DkENVdSQSncQKzTQwrhW4IRucn8OxLU9YZr1tvqOQCxjHtmB4wse2LWQ9QUCp6x3jCsHDfTa2bcFxbFhW1jMbs979Vbz3h2lkuqdVqaom3jfZ2cvBr8GFYHPQNB1qniOzliD5A24LEEJaZRCmmRKmAbc7maOQeY+znmqqKCp07yXnE5UQpqInEnVwHAuWlfNWhJxQOfN+lbAbVVUT3n8uCLgTUgZ/SQiIF2nIifEACHmp86/HhQBfgYSD01vtKgcbN39X5HPxiVOatzy40nNnb05yqFYC/nu5xuMLQ3bNQuMKe/7Z+BRFha6bzLxRFM8MdIVmxLW+0K8T9yIb0D4MIwlVVWFbOaFJcN8DBxce3DfjbRV/8fNlMvL3ZeNtSunepey4VZkSXLrzSWNZWRi6iWwuLaltTBXMSvarfDw7VgMhlnDS5e/je+p1sZIwUGHz67oacE4CCIUcOfjKn/87gsIoPwRZCQpNorBAEN9EhErl8/FoDNfMuFrtOA5yuax4FnIY0DRTnoMw3cxY2djYpKWeKecKDYhds7CwkhdFQhQYRkL4mgAE/ibEjzJxKIoKw0h6C0U2sBCAUqiaimSyk38+3QA800IWChEjC5gSsmbabkwJQsidhJA1hJAPpG3dCSHPEUI+9f7tJn13KSFkKSHkY0LIuFoNvBJwJ5tt5ZDLpeG6jpeXIKtyNPAvEwb+bcpkGoVTkr34Qb+FP1mZyixPfADC685V0xJHnrcC5guFaLAX138BuT2PCJOKXwdgK3kyWQ9dN8VKBSBwH+RzAmwyJ5P1MIwETDMlzALqaT7cZ+E4NjRVh6qqMAxT7Ms0KIZsNu0JZiL9L18/OH7HsYX9zx173LZnGgQzmfjvZ2ZOUmgDiUQnIbC4qSMLF8exoaq6MJGYc9U3A/h94iaTYSQDoe78ZxIFppX6v1Pej+dvyPvUGqX4GP4FYHzetksAzKWUDgUw1/sMQshwAFMAjPCOuZkQ0nxGUSuBJwJpuiEeIvEksD9h8n0H/gTkL7dlZYX6y+1e9pI0fzsJYaq0fF72EhARopNt7fDx8svBXnZmIhmeZ9y/3a7rSGo63zeJVKqe2feqhighQTx1XFFUgFKxL59o8m/hv1sRuRY8LKjCNJPCt6DrhrePKmL7uVwWuVzac5oqIvfDNJMiUsT/5ZOQCa1OSCTqPLtbg2kkheNXUdjE1XWD5RcYJgzdhKIowkGaStVD1XQmOHUTqspWdn5vWaiahUjZsey3mWbK0yqMgFnKx8ijNMVQLPIg/2bZnyA/09YKWTZrSlBKXyKEDMrbPAnAId7fdwOYB+CX3vZZlNIsgGWEkKUA9gXwWpXGWzZ03YRtW8JjzG84zz6zHV/FNQwz4O3nK7qsBQAQqjD7jmkc1JtAPILAJ0G+Wh/li/Dj88zv4Yftwiqj78giwtsvcgXMFFh+QlOkD4Ndw/GiAbqImth2ThoniyQQl01SVePOUTWwsltWRjjjDMmBxtVvbv5ks01CYwDghVMt4YjMf9H5ftwMYwJKkQSfv7+qJsX9MKUkJq7K89WcUiryFeTrcN9EFPjv0HRDTBLZtPCdxv455fBxS8AFkfxZ/rt0TbNyVOpj6EMpXQUAlNJVhJDe3vZ+ABZI+630toVACJkOYHqF1y8ZcjIJ8ZxLHHIWpLwvS+JhYSru/OOTlxDmqOSTmHnd2cRgIaWk8H4XcvjxnAbAf9mEx9xhiT+JRB2oMHUg8hnkye44jkg9lhNgEok64WB1vaxH5pzLiImdSNRBVVXv/6RIFVe14p5vPhEoNYRaLd9DPtl5opim6qAqDayAfJWNmpR8m5xcxf0GUaulfB5+Lwkhgd/BQ7MyarXycn9JtcD9Wq2Najsfo35BpKeEUnorgFuBWkcl2OrKVj9ZPfNCbLoB4iUAmWZKxKX5pGa5/Gzysbi6AXgqtZzIA7DUXUKIiHHzyQlA5ODzjERCFLbS5q1kiiKZEdJ2TTPYubysPdkHokVMZj9pSxcrI7e9QSlUVRUhM34fyoHmqeIAyzjkv9vP2lQD468E/N5WMonbYjIBLRc4+T6rtvodlQqGbwgh23rawrYA1njbVwIYIO3XH8DXLRlgS8GLiVzqFrzJqurXIuSvLLbFC5Z0kfSkKAqSyU6hlyD/WDmjkv8rq7DlPHT5XHIWJPf4RwkHjoCWJKn6jmMFBVGFMIxEoMApCnJtSozCcBwrcK/y8zVaC5WKtzkApnp/TwXwqLR9CiHEJIQMBjAUwBstG2IVUMIkLGhr6kZg0nFfRaUrAz+2JZOEmwB+YlH5iS88uQaA52ewxHfFclschyVZ8dyI/KrHQuC+AJbq7AauwUKXGfF9c2PIh3zOciZSoX35ueQx8zGVen4+pnLA7ivz9/gmcPD67aaIihByH5ijsSchZCWAKwBcC+B+QsiZAFYAmAwAlNLFhJD7AXwIwAZwDqW0JZk1LQZT53MwjGSLz8U85x0vWZT7HpjQU7wCId3LpHRDDjbXsQPRDq79cEEjV3ry8/MEMV54xp2zXLiwrEPmJDaMBBRF8zQa1csJ0IUpxIqsdJE0lJVqVvi5ediShx/zzT7bYinY3Ozj+7LJmZHyE1jEhPtteHWpnKgWtVDw31TO+8LMVyqKsPITv1ozGbHDZz7yCIJpJitaWbdUUEpZerBkwsiVfvyFlX0gYW+3H8rl+yuKhkymEcWTbAol4fiRlOLH8onq525w041PaNd1kUrVi4nM/R0KUcSEl9XxfNPNtnKgoKJKkwsG+d9sNu1no7aCGcQrNXnmZTq9GYoSTJJieTQVRyXizEcOtiq5njTeMgWDZWVZObdjg3pFRoCv2gquBtcVZdVc1Vc1HY5UDxCs9W9OHlNRfyGHM5s7TteNyBRpOc/BB4FpJARZDXe+sVJmtofppRAH4/u+is1CmVIo0dMuACYANN0Q5o8o0PKEJPFSqeWaBf4vLy4Td4PSwPcAPIIYTQgU/r28vVRwAeS6DtwCk78a4dBS0OEFQ+F03dLBowvFHjQXPvILytmEKjmfDFHOa2VFkVMulxHZeNxfkM02gVIKgzD+BJ4opGo6iHcfuCrP8vXZZNd1I6S28u/l0mqArcjFU5ULlwfnE7B4d0MQpbAUaV6XwHgWNE0XhVk8hs/PwQlmeN4J/32cz4HnL7B7qAihwnNOeEWpzzqV8AqkMsL05MKX52QArJScCx+xentFbrwSl3icHewrVxSJcSFf+L1g5oTlCWKeUAdwMhwrdFwt0OEFQzXAswCLTeIoh2Ih+5IQApRhwhFCkMulxYrK+BvYi2OaUpKPmQIg5xS4IvrAP3NVnVOjBSPMvjbAv+cCoxz1tfyCLf9eBOtPaGAihJx50grLK1H5caw0Pevt5jFI2ZYoEwd8arlkspO4h0zIJsTvYD6XjMj+tKwsKyX3tE+eJZmvkcnvgkyPx30pPNMycBc8oe3/zZDLZcSzbS1sNYJB5gMoF36JdDjzjZ+TO6Hyzy87p/jqQ6kr0oK57a+qutiHv0hsf1sUGzHGJn/FZCXAdt6LzhmN2KRWFCtvotLQv762QCP2Kw1y8Vh1UXhMjEJOC5g63HTkBVR8XFyLke8FJ4IBmIPUTwyzvG3M78JLqZlTVBUJatwxqygqTCMJy9NquFPVTxf3/TlRhXfi10nCID9lnheItUYBFbAVCYZikCdWFPhqImcOcmcYJ29xXTZBVUUNJAuxl5K9gDmP6JQTqfB6AoUooZwC7gEXtGi6ESg0YnRwtkh8ktV0Zu8C+ROhlqiNUCgOy+L1KqqYNH5iGk9Zd0U6Nedm4IVcWe/+AvCo8dKB1HIudLO5DAzPt8NJb13XQTabFu+BovhcFJxqj1HsuVBVP1oDQNRD5IO9f8XKxIt/X01sNYKhmGe5ucQcgKns/LOs0voJTEwYOJLqK2cf8iQknurLGYe8HaGphnBgMZXXT6nO5jKgniOQ25uKonocin7hECOacQLe9Y4Ev36FTQ5u3rB6CAuUMkclBYv3c4IVRlzLwozc/lektHXuG+JclNyfwTUIIPiOcPMFQIBEWK5k5Xkm8vugKKyqs7m086jf3drPcysSDMWr2gAIr7684ju2JQqtuEOKO/l0j2bNtq08qjXfNGDqvyUqH7n2wa7pJyhRSkUIjcfdOV04W5E427Ad0gI4fbpMdebXYhionJ+hfUHwNXhRDG42KUQR956zXXOb3BJsVUQ4FTmPI6fn43kbisq4LbmJwpG/qPCcjaj0ZV60xT/np5qXknqez7dRqFS7ltiKBENhrYAn9wCA4uU9yM5Ev8rS91HwMmIO5rlPCBo2P7kmLfgXDN1EziP0YEQxOVGzwLUL7jDL5bKMlt5jG855/JK+JkFgmknBg1hoMdlShIKsIeXzHTLwHAhHOPr45OT3hocc5WeQStUHCrDklHSqKDADY1Dgurz8m2tvOUSRqhTSQFuaACfXhzhOUDDwd6Q1nmnHS+MrE9lsWtj7msZyBTKZBqTTDcxkIKwAi3unAQitQOYeAAhyuaxH5uEJFGkVp9T1eii44FTzgO8UBZgwYOFIpipnPa4CruqG2aBcaXXyGZZrh+bPz6tN8xFcXf3zcE4IzrUJ+IKYE9sC8Aq0VI+HwfAIURKib4fucS7wf3n5e35Vpjyh+Yoe0BAdWzgf5WzOUmG3IGWZaz1y+BXIp9JrnXqTDq8xiKrEAuCrjuvYcF0ibE/AD1MaugnVczpxe56FAm0RS2dgtq1l5WDoJkBIIGzGveSct4GvjD7LNFs1DcMMmAacdUgWAmxF9AVTpU5GPgb//LLq4af6sl4Zuvgsp0QLshrNn2SMrckRGpHPRM0atHB2K35tQohw2PH0adUjdWGU+mnPfDOFv4f/dn7fuMOQc01UUrIsnxsofyIqasunFAsxB/uRcLQGFwOwFQiG5nj/+EttuTZUjw49kajz+kYwW992LK/ZC8+Pt4XHmk9w13WZD4JSKJ7UZyQqGgA/a09VNVCPXjyfTlxWYbnZwK+naYYfbxeclHKIsjIEXzTeRYqNi1Ol8/uYzabBmYySSV2sxlGTh5HmOiIHgCUI+aQuPP8iv46BN39hGhZr/JJI1ImELsBnWGb787HqQnhxRCUR1RrVqqVpLQFQCB1eMJQa2pFVRvbipgIqrQy+nTcz0bzMQTlDkHMFGkYCppEI8S4Iak+JZZk7JrmpQikVrdio60DVWRckn/Le9vwSfu4C8s7LOmPZEtELQvtqmhGooYia6LzsnEdBSgH7vQnBMK1peohAhQsI28p52ppPRa8oTIBq4BqILZ5FlN3fUQrcTG8xihIO+b1LaoWtQDBUBj/UGH4JOZehTMGlSz0o5P2AaE80n4Q+i5IeeLH9tmyuyLsH/Ea2tshm9JN4gv9SMZEURUMyaYhoCU/v5Y46n0Wq+ApbSXUpP2cxe513h1I8zYu38ePX5Ild7L6xtGNW/5GBYZji3udHlLZURDGUtzZiwVAh/Ew0BmYqBJl9i6EQ1yGHrF7LlYKKosDywp/ZbJNgRuZ9Ky0rC1VRkeO9JCgVvI2ygJMr9toDeL8G5mNJit8E5KdJe5WgHmksD8e4ji3qIkqBbI5UCzz/pOUsToUjD61VK9ExdK82QL4giPq+HESxPxX6mzMmmUYy4I0XDWGlzDreXbk9Qy6A4lEITlgrd97i0SFBZ0ep+H2yD6gU1MLsyK/GrBSyUAimSYc7WtUKscbQziETovJEKVGZ6ZHHEmIETAgeBWDqO8u0dGyramo2N0c4oqsmSwc32/L5MnjoU+46LsN1XahepifPBWlLP0O1Qoly3Qkvymrt37XVCAallerYqw2eZEO8PAmfdETxhAIRBCzpdAN03RDMQwCgKF7vS9CCZeCF0KVLL3Tu3BOjxx4tXkxFIZhw6lj06dxF7Pfi/IX4bNFn+PzDj7F48cvYuHFtVe8BG3OYMYn3ish5ncQIUURCE1B4okYRt7QX+C0RKw9BVwNbjWDYEmFbOU8rsEUok5sNbJUkIlUagOiRCPhRFt83UVqijqKoGDFiNI4+fTKOnDAafbt2Ra/OnUP7vfrpp9jQ2IiJI0di3yFDAACN2QyWfrMGTzw1H0/d+wiWLl2IzZu/a9E94JB5JnlyE+dekDNTc7mM11c0mGglr7oyS1MUuOO2tRm//DocP0cFQNVIe8tBh6d242AVb8X6BbY/yDUXQPDlrubLq2kGBg4cgSk/nY4ddtse43bbDZ2TxTkyFy5fjsZsFqN33DF67I6DBUuX4o4/34eFr8/DypUfi/yEcsE7jMvh1lSq3ousyP0k/Ya2ptQmjvFXZAW/I9B84VylpldL1H5KKZqaNiM/IpFM1ouS/Kjvy0DJ1G6xYGincF0XmUxjoFcDL7vO54KoBIQo2Gab7THmyB9g/KljMWrIEGzbtWtVxp6PVRs24Ov16/HArGcw8+YbsX79qmaPkfktHNsSvT94ViSvmFS8ikiewUmpK5ricnYn7mspZEblT+ZSTa4oIlh+bL45w7W/5s7L8lSCjXN5BCkWDDXAliAYRJ49YUSkABWrBeCnIXM1mu1anmDo3r0vdtppFMZPPRrfP+JADO7Vu/mDqojH3lmIl59YgIf+dTtWr/68pGN4G738dGUOfi942TlzWPLU6MJt6DjyBUOhVT9fEJTjqyhViOcLBlXVAg2MYsFQZbS1YChFxeQvHuchNI1EQJ2VX8RCnvoodOrUDfvsMwHjpk7A6H13xy79+0Nt4yzB9778Er8+98+YP/+BoqxGADOpeI1FMBri82xGTRozr+t0uRqWTPcmX7OWjst8wUCIIn53LBhqgGoIhuY6PpW6PzcJ5PRjmTY8k2kUbEH5kDWFYi+6oqjo3n1bHHjwJEy/9FQcMGxYmwuDfGQtC0+99x7uuupuvPbqI2hKb4rcL5fLCGYlJgxs0fPCdVmL+my2SUo44w10g9Tr/N5xlqUotNS/wK5b+X2OBYM8iHYmGFozbhzFGF3qyibzDsro1KkbJk46HVddcy56d+4MQ2vfwaesZWHmvPm4/mdXYPny9wvuFyVQ87cDfoIQQZgopRBkIV6p74bXxxQyeUoBp5eTO6XHzscaohzBwHsR1BqcNTiXy0b2wiyGKOGVSNRh/FHTMMPTEJR2GKcvhtUbN+Ku+5/ES7NfwKuvPix4EUoFF5S5XDoQui0F5WqDxcbQEjMjyvkop8QznpCK8xtiwcDBvdUtMSWq9dKEzmvlRPk1Jy1tNoxWIEQ5bOg+OPs3P8f0o8e1O5OhXGQtC3+d+QjmP/Qi5r14X9kTgZdoy4KhUmdta6OYKQHAy9OouFdKLBi8MyORqEMm01CyYIh6gVq6ChQ6PptNQ1VUWHZOlFiXX2OhY6+9xuHPd/4Oew4aVPEY2yPWNzbiHzMfxX03/RPLlr3bpmNpLfOyWFQCiAVDFc/td1IqR2OQ4868kKdS4o/msuzyPd3FUnnll1NVdUz76a9wzRXnIGls+eXGhfDld99h2gk/w+uvP15RghTgRzCqUf1YS/iCgXFmtJVgaL93qEoQJK9l1kqoWpAfgZG6VvZSRk10ZhY4QuBwh1oxjSFKKPzql2dAawOmotbEgB49cPt912Hqj3+F7t37VnQOTsDbHhbC5iCqRyPQWvUTHV4wcLgVrDQyB4KcgVgN8Cy9Ss95zHE/woU/Pw1zXn0Tlt22NGCtgUG9euHGa3+O391zC/r2HVrROTgxTntGfj/T/FLr1mpqu9UIhnKRv3JXy2nFGppYkWXGpdbb77jjKFz4m7Pw/pdf4ocHH4CU2boFNm2FNZs24Xs77YhZzz6AYUP3qegcrclpUAl4CbvcUSsoGFrHebrVCAa576SM/M5StR+HVpTdqbmXdtjQfXDH7FsxcuBATBw5stlip46EjJXDNb/5J1KGgdtm34qhQ/cq+xzttdyagz//tiy5BrYiwVAIrU3uUezFbG4siUQdLrj+/zBy4MBaDa9dY3Cv3jjv0tPxp9/dga51dbh99u0Vaw7tFZxwp1APj9Zij252RhBCBhBCXiCELCGELCaEnO9t704IeY4Q8qn3bzfpmEsJIUsJIR8TQsbV8geUinKTZWoFxsRU/mpAiIITpl2Akw8eXYNRbTnYfbvt8OurfozOiQT2GjwY1828EXV1Xcs6h+M4XpFabcAS1yo3V+T6jLaKopRyRRvAhZTSnQHsB+AcQshwAJcAmEspHQpgrvcZ3ndTAIwAMB7AzYSQNvf4tAe7kj/wShxg3bptg+nnnrDFJy9VAwN79hTkMQfuuCMOH3dSwX155Afw/Qu85wdLPy7NjHQcJ3LfqG2FulmXAqYt+F27w87HduJjoJSuopQu9P7eDGAJgH4AJgG429vtbgDHen9PAjCLUpqllC4DsBTAvlUed9mQW8G1FVgsvZLoiIKLrrsWu2+3XQ1GtWVDV1X84LzjAt56+TlbVkainmccF5bH02DbVmiiFXo+PKScj/xt1fJX8f4Zcl0ID3G3BspafgghgwDsAeB1AH0opasAJjwA8ML+fgC+lA5b6W1rU7REinO0VLBU0pcBAEaMOACnHnl4i67dkTFmxAiMGzdNfPaZpf32fjJ4L40oza3U51NIAOS/I6wytDxhIbfWa6vFrOS3lBDSCcBDAC6glEbXx3q7RmwL/TpCyHRCyFuEkLdKHUNL0FIyWLnpCf9c7kOr5CGrqo5TLpqBrqlU2cduLeicTOKI048QkSdVVaUsUjkpTGVtBOG3GUynGyJX4eaeVT5Hg3yNwH55/TNLeQcM3RSRq0q1zJaipNlCCNHBhMJMSulsb/M3hJBtve+3BbDG274SwADp8P4Avs4/J6X0Vkrp3qWmaLYUpSQ48YfAcw1kqKoq+hywfW043j5RDsWoh2nbOakZbGkYvvP+OGXCYWUdszVi8qEHYvjwAwLbOJmNoijimVpWRnTIBlgUIJNpEN/LjkPG0O13/Mp/ptxHkculQ+8Lh6azwrhyJrfcC6RQy8Bao5SoBAFwB4AllNLrpa/mAJjq/T0VwKPS9imEEJMQMhjAUABvVG/ItUMul0Y63YBstinAPCxLecexPCeiFuh85LrBMFLUw8xfPZpDMlmPUy/5UawtlIBudXX4/R2/Ryrl09rLFYm2nWOdy42kaGoDSG3z8hoJ87Ayq3ylgVAyY6tmHaldx4bmNSuWUchp2DzNHBNmcku+tjAnSmHwOADAqQDeJ4Qs8rb9CsC1AO4nhJwJYAWAyQBAKV1MCLkfwIdgEY1zKKVtm60B35QoVH7LOiWzRqrsJWC3hqucnLKdN1uRVUbmOwhmH0a9AOVyPIw7cipmHDO+rGO2Zhy0007Ye+/xeOml/wLwfTr5XcUB3gbO8lriJQLdu2Tw1HWZ8Ym/Q7wjGGdXKlYsVyrxTlSCU1sIhg5fXcmRSHQS7L1RnAa5XAa616JdBqVUUK2B0kBH5ihUqzy3vr4H7p/3ZEF69hjRuOmBObjszGklZQ5y9b9QwVI5FP3yuWQuynJAKUU63eCZP37DGbkyuLWIWraaoLisyvEHJtuTchPZ/ONMMwWFN4qVHJCyXenbrKWHlIrZnWPGnYgDhg0r6TwxfEw9aizGTTgjsE32FciIqleRwU2JUiCnultWtqLwuHw9f/ITYd4ACPxdS3R4wcC90m5EKqmqasKBWPwczN5MJjsFzIGoIijZIw4UckJaXsJMtADp1Kkbzrz4pC2Omq09oHMyiTEnjylYG1MuSvUJ2XZOvEummSy40FRyfVk4tZsEpy0dxSS+49hwmlHLZG+0/FCiohZymTY/NopUg6e5UuqGVhVFUXHaORfhwNiEqBgnjBmNnXfeT3yWtcRqIhim1IG8ruTVQLuNSmzpEIzBESqYqmoFacSD+/mT3k+c0QuqqByEkEjWKP6wozgeBw3aFT879yTo7Zw3oD2jW10dTr7wbKhKbdmx5WenKIXJVaqJmKilyijkP2jOUZhvCojeBFbOo2RTIvcrdt1i3+2x78E1axW3NeHE8Ydi2I5+Jr6cDVkN8JyXajvvGeMYKehLiH0M7QShTDZPkGi6AVX1exFWIxKhKComnjWhxeeJAfSsr8dpF/8okJdQzQzCUovhCl0zarvvuKatphkUQocXDNzHIDPvtgTVTk+Vzzdw4C4Yt9tuVT3/1oxTJo7BgAE7A2CCwbIyUp5AdZ5jc/Y/bzeYj2ILCeN8bFty3w4vGHjosJrOoGpCPt/ocUduVYxMtUaXVAqjxx4lPhtGsirvQSHzwXVdZLPpwPeGkShb/WeJTtzhHawaLebTqiY6vGDg9Vvlhq8K1d/XCt2798WM86bEIcoqQiEEZ59/Ajp37sk+S0K4UgHPe1tGgYe1s9mmQIatHKWShUZU2jQ3T/iC1lZ0AR1eMFTKqivX39f6wRCi4OgTpmGX/v1rep2tEf26d0d9ffeKj89/9oRE98Pkk9wwEjDNFFvdrbD5mi8IwqTDCkwjAd27BqWuH1mrUdg1CluBYCAoxJ9XCoqtENWCqmo46exJMTtTDbBNly44cMzR4jOPJhUDjzaUsyDImgGf8BT+OeRCrEIExDwFGwhquPLi1m44H7d0uK5TkH2nFBRaIaoJQpQ4b6GGGLrnUD86UQKHouup7y1hlCaERNbeAIBjW5EcpDzL0Qo4yv3Mx9akvu/wggGojinQHorNYlSGCUd8D6kU44gsJcSo6UbVncyBZChVi4w6ME2C+RYsKxM6LjYlqghepdZSrjzbzkXajNVC7HSsHXYdMAC77npQxcdX2wldLM2ZV2b6fCBuyGHZGujggoGUlPJcCjTNqJlJsdtuB2N4vzanxeywUAiB0vZE5SXB91OEp6bsg6j5OFrlKm0I3ryjpaphLQtZRh89dqtpM9dW2POwURUf2/p9HeSKymAX9FLLwFuKDi8YGNpvv8Jt+gzGtNOOwaerV7f1UDo0DpqwH5LJ+rYeRtGWiL6zk4oQZZiePi6iqhpImVyL+ahlxtluux+Kbbp0xiNPz6/J+WMwHLLzzthmm+3behhMMBQxB/gCJtPRlcsyXQ10cMHAJK9pplosGFDiAynXUbXHmD2gqxq69OxcydBilIDVGzcCADp37lHR8cVazpUbQlQUBapXsh91Lv7+RPkS4nBlFWHbuRYXUCmKUtTxKNPHV2KPGpqGs448ouLxxSiOex95Fqqi4PDJRze/cwSKNStyHLvsRaeQE5EQgkSiLuBHaE1hIKPDCwag9vnmpfA6REFRVJgp5nSMw5W1w/Ddh0JVFChq9V/3SpyBhUhdGONXBpqqg2frtkWzGWCrEAyVZ6+VfIUKM+Tq67vj+xMqj6/HKA0TR46sWbp5td8t13WhSiS1bUHrBmwFgkHXjYKpqW0NQggSeutkssUAth28TZvzHBQDD0e6ris0EUN6d+PMxyqCEWXULmOxJdh++5Goj/kXWg3j9tsLdXVdmt/RQ6n2fXMtAxzHKbk1oa6boFJI0s6rqYhZoqsIx7HbzFYrhmG77ha3n2vHkE3EYpEJx7aQyTQWFA6OY0FVddE2oBhc10E2lxFRMLl5cmsStbROGlUbgtFkFfYqx4hRErxQYlQRlqYbRTuUcfOFELXZ91DTdCiKKtoO6LrRJnkMHV4wmGaqpIq6UmBbORBFbfZ8juO0WT+AGLVBoXC1bVveZA5SsAG+2l9uxIr33GR/++8aZ3hqjXqJDm5KkKqyMCleA9PmkN+NKnIfRcOooyrP34/RPhAVrqxG7gE3fS0rGzhXTB9fZVTDxyBL8haDEIzYYVB1zhWjzRC1ALT0PZEFgVyCDYQ7oNUKHd6U4KiWOVEt9O83DAO6V85FuCViQ1MTvl6/HgDwzrLlWPS/Rei9XW+cd/KxMNtZ2NZ1XbiOXXP2rnxwKkFVUeEQJZAM1ZpZkFuNYGhv6NtvBwzoUVnu/paIB15/HTde+Ed8/PEbANjKZ1kZGEYSB43eA6OG7FDzMfSqr8duux2K+fPvb3bfbLYJrusWdSrWAjKVoNpKOQtR2GpMifaG7YfvVLVzbWhqwuPvvFPw+03pNNy8lcZyHNz25LN45ZNPYLWQ3ao5zFrwGi6cPA2LFs1FOr0Z6fRmWFYGppnCMcfPaDV27JRpYtvtSiPEYYzPZhtwMfhoCedkSxELBg8tpX4rF/sfs3/VzlVnmthr8ODI73K2jYsvvREbGhvztlu45fLr8MNDj8SVf7kLazZtrNp4ZPx3wQJcMuVsrFv3dWD7/vtPwk1zHsTNf7sUdWa48W9bgucOFMuSjFLp22OuTKXo4IKhdJvMdW2p5LV1m820FLqqFmyEa2ga/vLHn6N7p04AgDuf+x9mLXgNdWYCU849Cw0NG3DTVZfgB0edjaffe68q43EpxeZMBrc+/gwuPuFMfPvtSgAsp2TAgJ1x4W+vx78fuBEnHfC9dslc1Vzn6kLvVKnaRaXvVqXFepWgWR8DISQB4CUAprf/g5TSKwgh3QH8F8AgAMsB/JBSut475lIAZwJwAJxHKX2mJqMvAa7rlHQzdd0UDzw/xFlInePls6Xkr7uuK87bt+9QHLHn7iWNvxpIGmzla8pm8b975+KTJYuw6+x/4ryTj8XXn32N22+4EgsXPovTxizAGT//FQ6asB/23n4wutV1KovWvimbxSNvvY0Fjy3Ai089hmXL3hOMQ337DsUp552LaSce2aa+FSPZvCAihMA0C6eq8/eBU8w3By4IKKWhRkaFkqaau3atQZpbUQkbSR2ltIEQogN4GcD5AI4DsI5Sei0h5BIA3SilvySEDAdwH4B9AfQF8DyAYZTSgro6IaRGrlaCVKo+8mbatgU1Ly9BnryloBzBIL9EAwfugjcXzm311fKel+bjp0d9H45rY5ddRuP2B29GJzOByROnYsmSV8V+qVQX1NV1wWHjf4CzLjwJo4YMabYsfO2mTfjxGVdi3rxZyGSCZku/fsPwt4fuxNhddqnJ7yoHTy5ahBMPOhyO64cA5edY6jvAtEpbkA07jgPbzsEwEkUXEv4dP57AdzY6jlNUK6CUIp1uEH0tK8DblNK9S9mx2TtAGRq8j7r3PwUwCcDd3va7ARzr/T0JwCxKaZZSugzAUjAh0W7A+fvzHyB/IKWYEvwhllrtJl9r7ZoVePGjj8ocdcvRpVMdDG8l/OCD+bj7zjkY1KsXzr/+VzBNv2ajqWkj1q5dgf/ecz2OP2g8zj7nd9jQ1FTwvOsaGjB96hV45pm7QkJBVXWcdelFQijMXbwYM86/BuvzfB6thTrTBEJt4fznWOrCoKqq0DIppVAUJVTFW2zRpdSFquqByAN1Hbhu8XevXTE4EUJUQsgiAGsAPEcpfR1AH0rpKgDw/u3t7d4PwJfS4Su9bfnnnE4IeYsQ8lYLxt8som4k7xJUaH/mb3AKPgTG2GQFPpdjNzalN2HFitYnfz18xAgM6O9HQxybKXHHj9oXe+81PvKYTZu+xQP33oD/PDsv8nuXUlx+5T/w/PN3R65khx9+Gs47+VgALHryi1POw6y7rsP7X34Z2rc18MjMZ6pKqMojB6XY/7LQ4GnU8jbVq5NoDyhJMFBKHUrpSAD9AexLCCmmE0bpUaEZRim9lVK6d6mqTWWgot1YPqLUPUopcrkMdN0UZkZUtILlrOtCGFTiFHr81kfL2r8aSJkmfnbjr5FI1AEAPljwDpqyWSQNA7uPLvwYXNfB28++Hfnd1+vX48kH74kUCslkPc769WkieWnekiX4auWn6D9gJ+zQp08VflH52Lh2Y0tU8bJQij8gv8ltVbNrW4CyRkAp3QBgHoDxAL4hhGwLAN6/a7zdVgIYIB3WH0AwVtVqIEJVK0bmyRudZrNNUHmPQ4kcI39/gKmcMtuv7GAqBe++Ow/vtcGq+YP9RmH//ScBAN5++xmcMe2KFoUqHddFLpeO/G7YsH1wxK67AgC+/O47/O2XNyKdacCUn8xA327dKr5mS+C6pfErtFbqsW3nQouP67ol8zfUCs0KBkJIL0JIV+/vJIDDAXwEYA6Aqd5uUwHwJXAOgCmEEJMQMhjAUABvVHncZUFOFMlm0wG1n2easRAVK5+VoShKaH/5OP53uVJ+/frVuPKCG7A5w16AzZkMGrO1fxmShoHTfz0NhpFENtuEJx67FS98uARDRg4pGrfPpXOhJCkAUBUFhhHtwT/ytONBCMG/5s7DuNFHY8GCOdhuu+E486SjqvZ7ysGqDRvw6otPNLsf9zmUStLSUrs/ytelRHSiak2UcvVtAbxACHkPwJtgPobHAVwLYCwh5FMAY73PoJQuBnA/gA8BPA3gnGIRidZGfiRCRn75LEf+tmIqIiGkpJeFUhfz5t2Hfz32LADgllmPYcaM38FphfyJY/bcE1N/8ksxjsf/+QR+eNiBOOCA4woe8+7b8/Hd5s2h7X27dcNFf7kWqVSQGalz5544+NB98MW3a3HNTy7GihUfQlHUNtUWnl24SORUlIJizWH4945twSlC4gIUT57T9ejsSrlGg5HEZkscdXVQSlTiPUrpHpTS3Silu1BKr/K2f0cpHUMpHer9u0465mpK6RBK6Y6U0qdq+QPKRT5pC/crFHuwnHlHpokvBq6h8Ca4hc6dzTbhxkt/g8/WrMGYQ/ZB9226I2fXvtZeV1X89LyTsNtuhwIAMg1pdKurwx3/+QN22GHPyGP23v9Q9KivR9byVWwuxGYcMx7HTpkR2H/IkD0wcuB2eHL+m/jqq0+gaQZ2GLJnm2kLn61Zg1suu6Gg2ROF5ux9QdyaR6aSDznvofkFo7A/rDVp5Nvey1FjBEtY/c7BcvKSrptwbIvxQ3qTma8U7BgXhBBkMo2iuEYWEPkPjAsbztdn27lQXT3H6tWf4467HsGegwbhkl9Ow4/P+X3N0pNlDOndG4dMmgAAeOONp/DZmjXYpksXTP7RGaF9FUXFzvsPx53PzMUxR87Aqg0bAABX33wv7nj6efx11qN4+L5/iP0JUXDwpHGoMxMsqUfR8Ku//A0PP3dfm2gLqzZswJWX/g3vvz+vKufj/odivBuypiHC4M1oFraVi/QtuK4LXTe9aEo7ClduuaDIZBojfQS27T8EQggoPLUxj3WHCw4WvnSFP4Fnq0U1D3FdF5aVFTa7phlMmHgvFA9vcuHy8hPPoymbxfJvv8WTj9yJlevW1+qGBLDHIbvDNFNIN22G5Wkq9d2jO2J99/W3eP3x1/Hagkdx94NPs3HPeQE3XPgbXHPBBUinfTND103sN2Yv8bl3n0E44wcTMLBnzxr+mmhYjoOb/nk/Hnng74HtXHiXm55MS+xKJmsa/Brcl1XwmGaIgBgpTEwGWxUYRkKiyfJ/rqrqUFU/uUTXTZhmMpSwRCkTLjyMKScCAcwcyA9/sfi0ImLShBA2DlWDbVvI5dKBDlnvvTcPDy54HT06dUK37ttU78c3g5EDB4Ycjr236yXCmRyu6+CuG/+AT95nKc7L3l+GxStXYsmSBfhixWI0NQU1nFwujat/fIVIitJ1o834Fl795BPc+ZdrQrkLlVYu5nIZuGWEOx3HidQWo6IerutEOh3bInzZwQVD4bbh8qovI8osYJ2sXBa1kB6SbeVAiBIQMIAnCPKy4HxNQ/P8FQ5URfVCYzncfdVt6JJKYejQvdBWSOdycJ3oqsLGxg14+21W8rL0gw/x1fr12LTpW/H9gAE7o3NnXyNY8eUSLPmaRakPO+o4lnHYBnjuiVfQ0BDWwCilIIiOJhXzBWiaEXrexaCqKkwzKZik/fOEz6FpeZmQbdihvYMLhmjHou35E/jfst8hnwZcVVWYRhKmF5LjTkhKKRzXgWkmAy9XLpeBbeUCDqlcLiNWCPlalud7AIA33ngCN/5jFrYbNqjNJtGsl17BeZMmR04kVdUxduw0qKqOJUtewwP/nBPQlAYOHIGuXXuLz8lkPXbadlsAQKeudW3Sgm/1xo149O57Ir8jhBQkQhHO44hVnZuT5YIoakjbDJ6XmbLC1PVya+TvWxMdXDCwnhJM3Q8SaioK0yRc10E2y1ReXkMRMg1UDZpuwHUdZDINSKcb4Dh2qAKPJ8aQvLRWVjzDbjVLiqJeTrwDy8p5wsnFHdddjaOnjseO3oRqbTz/7+dCtQ4c9fXdcMqvToKiqNi06VvMuvu6gHrOf097gUspbrt3DpYtexcAm2i2lQtlGhZDoVW9EsHQHEEw00x9QUQUNVKbiRvOVAUEum5AU3WRXwD4ZgQredVhmilRVWcaSY+A03/J+QNSVc03TaSVX1Y9k8lOIROFF9dw+vmgA8mPUTelN+GlpxbU4D5EY1M6HZgoE86eKARYc8gXngsXPotVqz6P3HfJ60sCYc7WwMsff4y///bXALxnBIpsLhPpKOaaQTabbjYcXcjMqCTPwLZ8jZEn2HHkC5JS6zGqhQ4tGFgGmyGVtdoB9ZALCH7DHccJOAjzYVlZ6HqCtSrPazoqF9Pkg8fDNd2AqqqhlUiOUDx05+1YvnZtdW5AETRls7jsx1cHognphsIx/t12OxQ9PbIXGZpmoGvX3p4JFj3533tvHlZ8913LB10iLMfBbX+YKUwiHlnSdQPIayfH6mFsj8ot0Sw3gm3nIp9xqT0x5aQpzeurKo+lOZQquFuKDi0YOFitvOWZEKpQ+QvZbYZuQtP0kH+CZ6m1pJmMvEIZBqc086+xdu2XWLWxtnkMdz3/AsYcPBlvvPFEYOV/7u7nIguMCFGw/5EHBnwf228/EppmoHPnnhg16ujwRSiF5a3Oa9euwNknXCDyH2qNFz78EM89dW9ou21byFnZwOSXSVmcEuojHMcW+8mLTHPvg1xLY9u5in0GrWWudWjB4GeMuXAcG4QQv9pSMgH4Q+PsOmxlZ70GuaNRJtmg3qpT7sP1U1spALZSsRRtX9A4joVnn3qlavcgH5bj4Ml/PoH3P3gppBUV+j3JZD0mTjwwsG3yjDNx4IHHY8OGb9Czf0+kksH8h7Vrv8SjL7+Oxo3MZ/H2wmdx8c+uw8Nv1bTKHpbj4Lc/uSrkQOXPTy9AB1/MGSnDNFNiv0IRL/ma/jvGJrSiKGJBqEQ4NHfNaqFDCwZ+4zVNh2kmRSQB8JNNHNtCOt2AbDYdOE5RFCSTnUApq7rk6iY3O1w37NQsBh7N4IIAoEIwmGZKCAZV1TF+4gHVvRESvmtowJIlrwe2denaG0nDwPgzxmPw4N2hKCp69x4ovt9114Owx6BBgWM0Q8OxPz0erutg1bKvsfPwILmt49rINGbw1N0PA2DC+YNFr6JnfX1tfpiHNz77DEuXLgxs48lMmqYX5OEAUFLrN24uyhGEYuDnzDcfHdsSTu9y0Brt6YCtRDBwjUBRlFB4UdMNQRPOTQyZZ0HTdBhGMqT+aZohchyKvSC2lUM63eBxODD/QiJR52kJipdo5TuadM1AooZNTrbp0gV7jjo0sG2PvQ7DwJ49cea4MXju5Tl4dMF8XHzDHwCwCM7RZx0PhRA05XKB44Zs1xf19T3w4YevYJ/DRovt++8/CcedcB52GDYQY354pNj+5ZcfYc2mTTX7bQCwct06bN4c9GdYVhau6xQVCuWiFCeg0AIjTBRNN2CaqbK0huaKuqqJDi0YgLC65jhOgLqNUgpVZSuJ49iAR9gpQ1EUJBJ1ocalxeLSAoTANFOBY/n55KxMjpF7jMHwfqX1PigXLqX48Kuv0LC+QWwjRMHIw0aKz9t06YKd+/XD3dfcAoBxKpx29OEAgIdnPiOyNQFgvx12QM+e/bF69ef45K2PUFfXFQCrrLz9n7/GhN13R7dt/NqITKYRC19YVJPfVgjcv1TK6l4qTV850LwiqygQQmBZ2VavnCwFHV4w5Id4HIdFHDKZxoBpADD7zXYsUUhVDCLS0MzLVKyUO+rYATtsXxYzc6lY19CA6//9II4YdSjmPv9vsX3gwBGY9sOJ4nPOtvGXm2biw8UvQ9MMnHLhj9C9Uyesb2zEK888G3JO8gn34kv/RX0nJgReeul+PPv++wCA0fvshu7d+4r95z36FJ5+772alZf37tw5UAJOXQeGbkLX21fvCg6WTdv+GsJ1eMFg28EMMvYgVOE7yI9LE6KEnFD5WWi1gqpomHD6uKqf99vNm/Hrq27BlT+dgQ0b1giG5C5deuGCP18R6Enxm+vuxD///Bs4ro3xE8/ATyazMum3li3DF18sDpzX1HWMmfR9AED//jvimv/cglSqC9LpzbjirEuwOZPBLv374+gTpokw27vv/g+njTkKU6ddLkhqqokDd9wRQ4bsIT7zsuha9i7NZtOBCIVcmRu1yGSzwfyRfE00l8u0uRbR4QVDvjRmuQgmksl6GAbLYedVdrzYidG8+UxPjuuEYvS1EBSTJv8EE3evbr+JtZs24eixJ2PWnX8JrPadO/fEVXfegjPGjRHbFq9cif/846+w7RxSqS6Y/n+nwtDY/Vv4xmI0Nm4InFshBF17dwUADB26F/YcPFh4zTduXAvbcaAqCi655AxBJwewRK4nH7sTt89+qij7dCXQVRU/OOdUqN5z54ltUc8rvyS/UvBcCC4EeNGUoiggiirOzROoZF9HfuibVd46AU5RmS6gtdChBYOqqqEGH6riJzSJ0KSiCv5G9iK5gTwH00yGHFfVzkAbMmQPnP9/00ruNZGzbeRsu9m+kxnLgqrpItVZUVTst98xuPvZR3DmuDGBGoYX33oP333HCp923nk/HLLz8JLH77ouLNsWwkdRfIfqgB49cOes6zDxqOkirJnLpfGbH52NP94UzjdoKX48+SicMuOXorqVrcLpyNoHjlwuLXJM8p3QhSBqZhwbjmNDUTVQSgM08nJtBU+CkoVB/oRnUQ/FS79vRDrd0CbaQ4cWDPm5BpynMd8RpelGwEGkaXogtZmvONQjjpUTn/JToivBkCF74LbZtxfsP5mPrGXhol/fhAO/NwmTj/85Hnn7baRz0X6RAT16YMbvL0CPHv3QtWtvXHTNDbhv9t8im7+89uhrYmITohQsfDLNFAYOHxjY9vrrj+EHR5yApiYWdTh0wnGoT/h2ff/u3XHff67FzBeexJFHz0C3rn1Q37k7xk86qKTfXA5SpolLLjkDQ4aMBMD9QWpAhc/XInjImGs8rusEzFBeYMfTmJn5kBPvBV8o0umGgINWFpB8kXIcR2gXXGMF4Akjh2Xrajp0rz6nFCFVbbQ/r0eNwE0Fn6jFEg9N3gaEHYaMeIUJg/wiGsvKshr6MprPyNh++5G4dfZt2HfIkJKPMXUdv738R3jw0N2xYskK3Hnlv7DTLZdip759I/c/efSB2P3lp5HJ5bDX4MFQC2g7R804Eo88+HdQ6mLllx/h09WrMXSbbeC4LhY+5+cGpFJdcOiIoDbR2LhBmBq6nsABxx0Quo6qKBi7yy44+J6r8eW6C/G/txZh/6FDm/29D7z+Ol544EXscdhInD3xiGb3B5hAvPjm3+GcI49DLpcWviU5hV1+xvmLhWEkvLJ6nnjGNEpO6MPJe7jmqSiq59B2habCu53xLlUsCpEDQFkWrqqBKCp0hTOZux6PhyIWG8NIeNW5OSlTtvbYagSDbedESjRXG/kD4OFEQpTIlGAAIhmJ5S14lZIuy6jUE3XgPSj4i8Ez3fgKFBUuGzJkD9w6+zaMGrJD2b+nczKJM8YeBowFrHNOLRrJUAjBbgMYo7/junhxyRI8/+SruPickwOr+uBevVBf3x2bNn2LNWtXYN7C99D/8O548aOP8MLcWWK/3r23g+n5HvoM7ANNMwKr5F57HYHj9ws2H3MpxZKvvsIHK1fi6D32wJDevTGkxEneo74Txk45DOM8KvpS8cNRo7D0t3/E9f93ESwr0yypr+tR4XNzIBBmlEwGHmbm74+qasJ8IBIHCKvsTQvfFQ+F83Z4XGABitRhWxMl3zL/B2kl5iaOZntXtsogata7EjCNJDTdEKFJzsVgGEmoqipUOvklcBxWim2aKaiqGkiJ5tI7kegUoH8DfAeWXMlZqD8FNx/K0RRais2ZDC68+Ho89sAdsKwM5i1agF36+y1AXEox5rBT8OabjGK9vr4Hevbsj3R6M1av9isnzzzvCtz4+58BYAlFo/c6BGvXrgDA7Oi/PTYbp4wOplDPWbgQ5x9/OtavX42RI8dg2q9n4JRDDiqovVQLm9JpnHDc+ViwYA4sq3gUhFKKpqbNAFiZfDLZSTxLXuPA7H0ivstkGsW7xM/Bn7VlZZHLZZFMdpIS6JgpahpJEEWF41ii61nCW2Dy352o8VWI6vWu3NKR8zzEvp2niWQlXgYdCk/aORi6KR6249hwvQxHJgwILCsTqrTj/TBl52ZUxWWPHv1wy0PlmQ8txfyPP8aJP/gZZv3rOmze/B2GDt0HfbsGiVkVQnDk6ceKz5s3f4dly94NCIXOnXvipDP9oqmErqNbN5+Obu+9x2PyfqNC1//vDQ9i7doVsO0c3nrrKVw0ZSr+dv9jkb0qqonOySTuf/iv+Pnv/lSw/wUHz3TVNENMfABiMdB1UyLssUTVJvcb5Jds52e18mpfRVGR9fwKrKJTEyYq92W0tk8hHx1eMFDqimQmDiEkPHUuP2XVNJMhZyT37Nu2BdNMwjRTIlLBw1OlkngcfcI0jBoyBOlcrtmoQkthOQ6uuOEuTJs4BfPm3Se1pR+C7hFl1OPH7IcePQpnXp7/u6sDps+7K1Zg5cqPATAH3mmXnxXJ75htCnrWm5o24qbLfoPPvvmmot9VDuoTCVw8fQpGjTq62bJlXTdhe3UM/J1hIWwHjs0iPDK7E8+H4f4Fts3n/WARhiAZcSJRx0r3vYWHlYSb4lhN0yPzLmI+hioikegUKFKSJ65jW0x7KLFizTCSnnmhBSITPAxVilm2yy6jceFFp2HJV1/hsIN+gPN+8aeaZAG6lOJfc+fhx+ddgxt/80us/mZZSccN79cfl/39egwfHl3Ite1gXzuwHAd3XDsTTU0bQYiCAw44Hj+I0BYKYfU3y3DrrQ+VvH8l2JRO46E334TtOvjTbZejR49oBy0Hz0jlaj/3I+U3m1W90CT3OwAskiMXRnFNNX+x4HUzha4vh8aba1VQK3R4weA4ViSVNy+z5VTvzalusnnAHU6lkLQEz6Fg/8PHYHCv3njw4bl4/4OX8Mrcp2HXQGv4ZuNGXHX2+Zh1958DjsHmoBCC6UeNw78evhXbblvcKWo7Dt5fxErEu3ffFr+/+eKCeRjd+/aI3L5ofm26FzZls7jkmltwzIQzMGPCsfjHrMcwuFcvTDrpzEitQZ5whpHwzE0Ntp2DqmowjIQIdRtGAoaR8JyIhqiWlTNpS6WQi0q+KtQ9vaVh8XLQ4QVDPqMvj0gAnDQlJyolOc9C1EMJxJ4di/kdvP0sKxtInpEfoPyQVVXDsSezlOdhew8DwMJ8S2ugTr+zfDkaNpfWn+IPd/4X1/37wcC24f364YdnzwhNItv2hdgTixZh1aqlAIBJJ54ZcGTm46zzT0CXLr1C2z9asgALly9H1rKwOZOpms/h6hv/jZuvuQxvvfUUMplGXH3BuZj9+puYcc4P0bNn/5LOoaqqiDZYVlas3DwsTSlbVAzdRCLRCYlEHShlAiKXywR8Bvmp0XI+hSy4bY9NOp3e7IXUo9nMa40OLhjCNpmqalLRFKP2ZiGk4t2FbDsnHJm8BwU/D2sow16afOeR3OZs+PADhMORX+Kbb5bjNz+7IXCt9Y2N+HT16oJJS6VgSJ8+MMzCce9stkn4N+bOehoP/uPeUFPdaVOPCTgWAeCJW/ymsP+5dhYymUb06NEPZ50zueh49ho8GP377xjavmHjWvzy7Ctw3LHnYtSeh+GiK//WIqanpmwW97w0HzNvvjHAdmRZGdz7u7swsGcP7LpLOKkqYGJKnAdyVCK/poEQViUr12JonhbK/Q1Muyjs9Mw3HRRFFYKmtUhZIsfVZlduI0TV5Luu66mE0aFF13VhGAkkk52Eo1I2P0ScGjKFV1jrMM0UkgbzR7z6yGti++bNGwL73XzPwzhw9/1w+TW3VvIT4bguLjv3z9i48duC+7z55lP4dPVqAMC1t16Bq++4GnV5gqRLMtyAp6GBjXXpN9/gvXdfAAAc9cPTsUv/0lbhfFDq4rXXHsWLL87CihUf4tbrrsDMR58r+zyW4+Ca2+7D94/5CS74/gnYvHkdjjvhvEDznA8Wz8fSb9bgyOkRVHQS8p3I/PnK27gvIX8/XTeRStWLknzBw+GFzLnmkF8fIUPV9IILVGthqxMMUaCUilqJqO8Y0akXrvSkODc/8vPvVVVFMtkpMgty3bqvsXrjRqxvbMQb814oOJ7e2/VuES8gIQTf+/4BSCbCUQcOx7FER6U9Bw3CITsPD2kMXevqcMDBxwS2DdqJ+R2+WrcOa79diS5demH6uSdUrW8EpS6euffxsqM1b3z2GW664nK88spDyGabsN9+x+C2Wy7HGedfJsyhjRvX4tn/vY7Bg/qGunOXiwBjuAfZ11RoUvP3h5ut7Fy+qZBfDZx//jgqUQWUKnG5RJchawPyRJeTmliRllawDVk+li97H0u++grL1q7FJ58U5j6cfMiB+NWNN+GqS88uafz5UAjB+Sceiwv/8CdRZdgc7nnhRYwbc3LA36GrKrr1CeY6HHDs9wAAs+95Gq7roHPnnuhXhUa1mmagb9+h6Nt3aNn8BI3ZDP7x+3+LzlipVBeMOekImLqO4086Ap06Bcd34LBh6FWin6EQ9LxOY6VALobiwiP/nWGdrgr7FGLnYxXAcxhKmbRAsCmJoigiS7LYC8CluEzoWi4aGtYHuAm6plK48JTjQ6p9OVAVBWdOnohhO+7b7L6O6+KlB+fj3XdfwKPPvSy2b0qnsfDlMDFtzrbx9adfgVIXB4+bFOJxjHIgvrN8Ob766pPI6/fvvyP+eO89eP2deXjr3Zfw4KN/L5mspimbxc8uuh6PPfxP9Os3DOPHn4l/PfsofnoC03Tmvfi2IIbVNANde7dMU6gEfJHJZ+ySqyx5B7NiGkFrRiU6fK2ErhlQVM0jYo1+2SilcGwLtmNBU1kiEwtXqbBtJ5A5GQVelZfJNAZyJkrFBx/Mx6IvvsDoHcPOuZagZ309Trv4R7jszNdDtON9+gxG1xSzv9c1NODFZx8GIQQjdvOLmrKWha+++lR8VhQVqq7hxY8+wosvPgBNMzDxtHC9w5kzrkKnbp1w3dXnCz6HF19eiA0b1oT27dNnEG6efScOGz6iot+4rrERdtbCTXMewujhO2Fwr96B7x3LDvh7XMfLPakhr2Y+mlP/OYs0UZsLd3Mi4tpTyHdojQEAbMcP+cjOHxmObSGby0DXE56HWRMx6XynUzGU0gxkXWMj+nbrhoED/epEx7Fx9fl/wpctaMryxbffRkYxTpk4JjIasNNOo9C/e3cAwKOvvoG1a1itg+wr2JhOBzz0A/rvhAkjd8ebr7yLpqaNSKU6Y/ve4RDk5x9/iKcemonGrK86f/Dy4tB+mmbgpJ+cW7FQAFg5953/vAKnHXxQSChsaGoKVIXadg6bvtuElGli7HHHV3zNaoNnQ5ZCVttaDsmSBQMhRCWEvEMIedz73J0Q8hwh5FPv327SvpcSQpYSQj4mhFSfq6wMyKoXpVTECXlNveM40PRgbjyPMPBKu1JACGtcUuzBOa6NJ+98Ctt06RKYrJS6eOWV2TjrxIvw9frScg9krN20CSdNmo7Jx52Ptz4PtonrkkrhoCOOKXAkQ7YpK+jeZMx5/hWsX79afD5o3NHoWud7+QcP3g3b9egZOMYtoO6edsFk1Nf3gKow+nxNM/Cji6/E5edPLek3VoJ/PfosnnnmLvG5vr4H9tt/N2QtC+++/GbB41pSp2DbVrNt7uTrRO1bqBmS7LCsNcrRGM4HsET6fAmAuZTSoQDmep9BCBkOYAqAEQDGA7iZENL6GRoeZCZnmYBVVVmxDGf2kfMN+HG82g2oHpWbbbEHO2yvnQOOQUpdvPLqbEw74UKsXLeurHNuaGrCl19+hHnz7sNZx52Ft5f56c8KITj7/BPQrWufwDG6Gc0dwSnil61dg1k33iHU8G7dtsXZ508JaBRLlrwmQp4cb33+OT7+OJzNeOjw4Xj81efxyIKXcN7lf8BF116Pqy4+W5ga1cbaTZtw73W3BMyIpqZNWPzxMjiui6+/Whp5HC/HL4Zi70KhOgcOmc8xKnmJ+RHatoAKKFEwEEL6AzgSwO3S5kkA7vb+vhvAsdL2WZTSLKV0GYClAJr3gNUEvuc3XwLzh1Io+SQ/5FStMFHTxkZW2PTLszBwUJhF6dXXHsGPTrssFDoshq6pFHr1YlmHS5cuxNnHTccX3/o5DCMHDsSFf/qDoFUzjCQm/yysSruug+suuBYvLlmCv974H3z4IXM8ppKd8Ys/XYvdt9uOfe5cB1XR4NgWnnnuNVHrsWbTRjzy4Fw0NW3CiBEHImUE7fi9Bg/GocOH4+qLz8ZlPzqlZkIBAN5evjzUeMZxLGxet7nAEaWjJe+CrptF+1zyCsy2zGEASnc+3gDgYgCy+7kPpXQVAFBKVxFCuIHXD4Dcsnmlty0AQsh0ANPLHXA54IkqjuOEJjoHk9itp9C89tocrFx3Mfp16x5pU1Lq4tVXH8Ef/7EvZkw9Fn1LCAWmTAP19d3F50+Xvo0775mDK392BgAWoThvyiQMHNof9/z+P9h34r44dq+9xP6j99kNPXr0w3fffYV33nkOx35vvsgEJUTBERNPxbknHCO0heOPGI3ruvTC+vWrcP2vLsJLs/8HXTewcuXH+Pyzd9GnzyCc87uzI6ssWwOLvvgCl027KFQj0rlzTxx44B6g3n9RqHWeQCkTvq2FAlCCxkAIOQrAGkrp2yWeM+pXhZ4CpfRWSunepRJHtAQyISyH67o1rXsvZAtu3rQOs5+ZD0PTcNQpJ4Sq9gDm27ju/36GsQcejX/Nndds/YBCFBx87BHCXKDUxQO33Y41mzZK+xB8f++9MXv29bjkrCkBgpQR/frhp1f+Bskkk/uWlRGe72222R5/vOEXARNim65dccBoxvqcyTTi5ZcfxAsv/Aeffvo2evXeDjc9dBeOKJNtqVpoymZxwWm/xMcfvx76rq6uK3bo3RuPv7MIq1Z9HnF0y9HWPArVQini8QAAxxBClgOYBeAwQsi9AL4hhGwLAN6/PBa1EoBcTdMfwNdVG3EZ4Mw4+Q/Lz2Z0KiJzLVRoFbxGNOW349qY/+BLsBwHZ5x6DPr1G1bwHF988QEunzoD97zwYtFrJQ0DV5x/OsZNOlVs699/x0B36mJQFQU/P+14XPTH6wIpxADr/N01Fey41ZjNYuWKTwPbUqkumHjUdPzz0XswocoU+KVg7aZNeOrdd3HGtCvwzqK5oe8JUXDUlFPRo74ejZubmmVzqhR8tW8PzGgtQbOCgVJ6KaW0P6V0EJhT8X+U0lMAzAHAXcpTATzq/T0HwBRCiEkIGQxgKIDa1NY2A87gywUEB4sgpARZRikl0zKiNJB8FCNtef/9F7Hiu+8wsGdP/HDG9EitgWP9hm/wq1PPxo2zHsX6xsbQ95vSaZx48qV4ctEidO7RGYQo6N17IH5yzU/LSpBSFQW/OH0yLr3hrxi4nR8+XLN2BZ5fHAw1PjDvZXzwwXwArBP2eZf/Efe98CRmzrwGY0ZUHnqsBOsaGnDjrEcx4bAf4sSDx+LxObeEeoAALNJ04ulHQSEEL/wnOh29Gqt9ezADqoGWGFTXAhhLCPkUwFjvMyiliwHcD+BDAE8DOIdS2joxlhBYCSwLSwbDcbznhIxCvAxcqORymRBpRqHQUtTxHN98sxy33c4ISr7//cOQMOuiDhPYuHEtLp9+BsYefLxo/cbx1rLP8cLcWZh140OYM/MujB17Oh576YmAD6FUKITgwlOOx33PzMKuux4CVdWRTm/G3VfPDOzXqVs9evboh4lHTcfM/z2J3150FsaMGFGT1nrFsHrjRhwzfioun34GPvpoQVHeiX33PQq7DhiAr9evxztvzovcp5r+hZYKCN4MuVgvjFqiw5PBAiRAo8UanLKIhEwp77ou0ukGmB4hB8929LkbHEbiaaZAXUfkQzBa7yQcxwo5Ex3HEfvmF1XtvvtheO5/98HUdYwfNw2vvfYoSkG/fsMwe+5sjOjfHy6lOO30/8MjD/0dANC371A89fITGNK7dzNnaR7rGhrw78eex1N3zcHkn5+Es8YfLr5zKcW3mzeha6quppGF5nDNbffhmgt/WlIm4E2PPIxphx+K+R9/jKP2PTBSq2hPYM1uc+L91HWWTh2TwVYJhuGTuvIa+nyhAPgdqFWP1iufT0FVNdGEhng09LylPdsnPEGo60DxiD7zsXjxy3h7+XKoioIjTj6q5N/z1Vef4LZ/PCB+z6cfvQsA6NmzP66ZeXNVhAIAdO/UCReceCwef+q2gFAAmGbRu3OXNhEKjdkMbn/6eUyYMB03XXF5SUIhleqCgX3ZfXny4XmtklJcKXh9juXV7HBy2taqquTo4IKBhBicqPdS8E5DMmR/Q76ZIQsR+W/uuIxKatH0wg/UcWzMvvdpAMBBo/fANn1K60IFAI/MvAOLV64EACQTnZBKdcFF1/8Bx++zT8nnKBWtbR5EQS4wu/ZvM/GLKSdh/vz7RTVlc9h++91w4LAdsSmdxoJnXqpJApHM8AU073zk6fn5mbmu68Bx/PoO3pWqtX0XHVwwMB+DXELNi2dYC7DSvPYcvPpSNj/4ecseGXWx4Lm5WNfQgFFDdsDJPz235GO/++5r/P36mVjX0ICrb7sKdz83B2cdXVrzli0Rv/vznXh72TI898EHePXJeWVxWI4adTSuufMPMDQND72yAO+++7+qjo1T92WzTYLfAvD9UcW4FfLbJbKmRkwLLZYE1Rro8NWVrtepmpCW32gKKpqXBrZX6Kf5cMlruOfxuTjj2HHotV1vaJpRUgNTSh3cc9sfMPyAEZh6JFPzc46DXCvl0bc2pk6bhCeffxW3XPV7rP5mWVn3e9Cwodh7++2xoakJrzz8SoDFuRrgEzu/WYzqRZqiOBcKQZVIXHiPEv65tbEVOB/FVWp/iQrQvfs26NSpG1zXxapVS8sqkunVsz+Sqfrmd+wAcB0bX6/6rOyQYuf67ujajSV+rV79OXK51u8cXRn4lMh/b1s0VUp2Pm5FgiFGjK0ecVQiRowYlSMWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECCEWDDFixAghFgwxYsQIIRYMMWLECKEkwUAIWU4IeZ8QsogQ8pa3rTsh5DlCyKfev92k/S8lhCwlhHxMCBlXq8HHiBGjNihHYziUUjpS6pZ7CYC5lNKhAOZ6n0EIGQ5gCoARAMYDuJkQolZxzDFixKgxWmJKTAJwt/f33QCOlbbPopRmKaXLACwFsG8LrhMjRoxWRqmCgQJ4lhDyNiFkuretD6V0FQB4//b2tvcD8KV07EpvWwCEkOmEkLe4aRIjRoz2A63E/Q6glH5NCOkN4DlCyEdF9iUR22hoA6W3ArgVAAghoe9jxIjRdihJY6CUfu39uwbAw2CmwTeEkG0BwPt3jbf7SgADpMP7A/i6WgOOESNG7dGsYCCE1BFC6vnfAI4A8AGAOQCmertNBfCo9/ccAFMIISYhZDCAoQDeqPbAY8SIUTuUYkr0AfAwIYTv/x9K6dOEkDcB3E8IORPACgCTAYBSupgQcj+ADwHYAM6hlDo1GX2MGDFqAkJp25v3hJC1ABoBfNvWYykBPRGPs9rYUsa6pYwTiB7rQEppr1IObheCAQAIIW9JORLtFvE4q48tZaxbyjiBlo81TomOESNGCLFgiBEjRgjtSTDc2tYDKBHxOKuPLWWsW8o4gRaOtd34GGLEiNF+0J40hhgxYrQTtLlgIISM98qzlxJCLmkH47mTELKGEPKBtK3dlZgTQgYQQl4ghCwhhCwmhJzfHsdKCEkQQt4ghLzrjfPK9jhO6doqIeQdQsjj7XyctaVCoJS22f8AVACfAdgegAHgXQDD23hMBwHYE8AH0rY/ArjE+/sSAH/w/h7ujdkEMNj7LWorjXNbAHt6f9cD+MQbT7saK1jtTCfvbx3A6wD2a2/jlMb7cwD/AfB4e3323vWXA+iZt61qY21rjWFfAEsppZ9TSnMAZoGVbbcZKKUvAViXt7ndlZhTSldRShd6f28GsASsirVdjZUyNHgfde9/2t7GCQCEkP4AjgRwu7S53Y2zCKo21rYWDCWVaLcDtKjEvNYghAwCsAfYatzuxuqp54vACu2eo5S2y3ECuAHAxQBcaVt7HCdQAyoEGaWWXdcKJZVot2O0+fgJIZ0APATgAkrpJq+mJXLXiG2tMlbKamVGEkK6gtXd7FJk9zYZJyHkKABrKKVvE0IOKeWQiG2t+eyrToUgo601hi2lRLtdlpgTQnQwoTCTUjq7PY8VACilGwDMA6P8a2/jPADAMYSQ5WAm7WGEkHvb4TgB1J4Koa0Fw5sAhhJCBhNCDDCuyDltPKYotLsSc8JUgzsALKGUXt9ex0oI6eVpCiCEJAEcDuCj9jZOSumllNL+lNJBYO/h/yilp7S3cQKtRIXQWl7UIt7ViWAe9c8AXNYOxnMfgFUALDBJeyaAHmCEt596/3aX9r/MG/vHACa04jgPBFMH3wOwyPt/YnsbK4DdALzjjfMDAL/2trerceaN+RD4UYl2N06wKN673v+L+byp5ljjzMcYMWKE0NamRIwYMdohYsEQI0aMEGLBECNGjBBiwRAjRowQYsEQI0aMEGLBECNGjBBiwRAjRowQYsEQI0aMEP4f3J7qt+mRj2kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "in_dir_path = \"./infile\"\n",
    "mat_dir_path = \"./material_info\"\n",
    "model = import_dicom(in_dir_path)\n",
    "model_kVp = find_model_kVp(in_dir_path)\n",
    "\n",
    "interval = 5\n",
    "anode_atomic_number = 74\n",
    "filter_thickness = 0.5\n",
    "\n",
    "\n",
    "predicted_CT_numbers = predict_CT_numbers_poly(mat_dir_path, model_kVp, interval, anode_atomic_number, filter_thickness)\n",
    "tissue_ranges = find_tissue_ranges(predicted_CT_numbers)\n",
    "\n",
    "mass_densities = mass_density_dictionary(mat_dir_path)\n",
    "\n",
    "density_map = find_density_map(model, tissue_ranges, predicted_CT_numbers, mass_densities)\n",
    "\n",
    "plt.imshow(density_map[23], cmap=plt.cm.bone)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}