Untitled4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "050fc33a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SimpleGVXR 2.0.6 (2023-05-23T20:52:16) [Compiler: Microsoft Visual Studio] on Windows\n",
      "gVirtualXRay core library (gvxr) 2.0.6 (2023-05-23T20:53:11) [Compiler: Microsoft Visual Studio] on Windows\n"
     ]
    }
   ],
   "source": [
    "import os, sys, time, math\n",
    "from natsort import natsorted, ns\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "import cma\n",
    "from PIL import Image\n",
    "import glob\n",
    "import re\n",
    "\n",
    "\n",
    "from gvxrPython3 import gvxr\n",
    "from gvxrPython3 import json2gvxr\n",
    "\n",
    "from gvxrPython3.utils import visualise # Visualise the 3D environment if k3D is supported\n",
    "from gvxrPython3.utils import plotScreenshot # Visualise the 3D environment using Matplotlib\n",
    "\n",
    "from gvxrPython3.utils import loadSpekpySpectrum # Generate and load an X-ray spectrum using Spekpy\n",
    "from gvxrPython3.utils import loadXpecgenSpectrum # Generate and load an X-ray spectrum using xpecgen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f7aef2c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#data_path = \"4thOCtober_image/\"\n",
    "data_path = \"C:/Users/snn23kfl/project/4thOCtober_image/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7be97d94",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C:\\Users\\snn23kfl\\project\n"
     ]
    }
   ],
   "source": [
    "current_folder = str(globals()['_dh'][0])\n",
    "print(current_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ab9ea012",
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_images(image_paths):\n",
    "    \n",
    "    \"\"\"Average a list of images.\"\"\"\n",
    "    # Load the first image to get the shape\n",
    "    sample_image = cv2.imread(image_paths[0], cv2.IMREAD_GRAYSCALE)\n",
    "    if sample_image is None:\n",
    "        raise ValueError(f\"Failed to load image: {image_paths[0]}\")\n",
    "    \n",
    "    avg_image = np.zeros_like(sample_image, dtype=float)\n",
    "    \n",
    "    for path in image_paths:\n",
    "        image = cv2.imread(path, cv2.IMREAD_GRAYSCALE)\n",
    "        if image is None:\n",
    "            raise ValueError(f\"Failed to load image: {path}\")\n",
    "        avg_image += image.astype(float)\n",
    "    \n",
    "    avg_image /= len(image_paths)\n",
    "    \n",
    "    return cv2.medianBlur(avg_image.astype(np.single), 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "355e14c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def flatField(img, white, dark, epsilon=0.0):\n",
    "\n",
    "    temp_white = np.copy(white)\n",
    "    temp_img = np.copy(img)\n",
    "    \n",
    "    test = white - dark == 0\n",
    "    temp_white[test] += 1\n",
    "\n",
    "    if len(img.shape) == 2:\n",
    "        temp_img[test] += 1\n",
    "        return (temp_img - dark + epsilon) / (temp_white - dark + epsilon)\n",
    "    elif len(img.shape) == 3:\n",
    "        flat = np.zeros(img.shape, dtype=np.single)\n",
    "        for i, proj in enumerate(temp_img):\n",
    "            proj[test] += 1\n",
    "            flat[i] = (proj - dark + epsilon) / (temp_white - dark + epsilon)\n",
    "        return flat\n",
    "    else:\n",
    "        raise IOError(\"Bad image dimension: \" + str(img.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "85196c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "#path for my recently acquired images stored in PhD file\n",
    "\n",
    "\n",
    "\n",
    "dark_field_paths = glob.glob(data_path + '/darkfd/darkfd_*.tiff')\n",
    "dark_field_paths = natsorted(dark_field_paths, key=lambda y: y.lower())\n",
    "\n",
    "white_field_paths = glob.glob(data_path + '/whitefd/whitefd_*.tiff')\n",
    "white_field_paths = natsorted(white_field_paths, key=lambda y: y.lower())\n",
    "\n",
    "raw_image_paths = glob.glob(data_path + '/raw_images/raw_image_*.tiff')\n",
    "raw_image_paths = natsorted(raw_image_paths, key=lambda y: y.lower())\n",
    "\n",
    "if len(raw_image_paths) == 0:\n",
    "    raw_image_paths = glob.glob(data_path + '/rawimages/raw_images_*.jpg')\n",
    "    raw_image_paths = natsorted(raw_image_paths, key=lambda y: y.lower())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c0c283d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Average the dark field and white field images\n",
    "I_dark = average_images(dark_field_paths)\n",
    "I_white = average_images(white_field_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1347007b",
   "metadata": {},
   "outputs": [],
   "source": [
    "I_raw = []\n",
    "for fname in raw_image_paths:\n",
    "    I_raw.append(cv2.imread(fname, \n",
    "                           cv2.IMREAD_GRAYSCALE))\n",
    "I_raw = np.array(I_raw, dtype=np.single)\n",
    "# I_raw = average_images(raw_image_paths)\n",
    "\n",
    "# for image in I_raw:\n",
    "#     image /= np.mean(image[:,0:200])\n",
    "\n",
    "# I_white /= np.mean(I_white[:,0:200])\n",
    "# I_dark /= np.mean(I_dark[:,0:200])\n",
    "\n",
    "I_flat1 = flatField(I_raw, I_white, I_dark)\n",
    "I_flat2 = flatField(I_raw, I_white, np.zeros(I_dark.shape))\n",
    "\n",
    "# Clamp pixel values\n",
    "I_flat1[I_flat1<0] = 0\n",
    "I_flat2[I_flat2<0] = 0\n",
    "\n",
    "I_flat1[I_flat1>1] = 1\n",
    "I_flat2[I_flat2>1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "07154898",
   "metadata": {},
   "outputs": [],
   "source": [
    "I_flat = I_flat2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d19fecde",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sod = 43\n",
    "#sod =47\n",
    "sod = 48\n",
    "sdd = 61"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3dcff2fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ti90Al10 = 5.68\n",
    "x_src = 0\n",
    "y_src = 0\n",
    "z_src = sod\n",
    "\n",
    "x_det = 0\n",
    "y_det = 0\n",
    "z_det = -(sdd - sod)\n",
    "\n",
    "x_obj = 5\n",
    "y_obj = 0\n",
    "z_obj = 0\n",
    "\n",
    "alpha_x = 0 \n",
    "alpha_y = 0\n",
    "alpha_z = 0\n",
    "\n",
    "x_rot = 0\n",
    "y_rot = 0\n",
    "z_rot = 0\n",
    "\n",
    "x_default = [\n",
    "    x_src,\n",
    "    y_src,\n",
    "    z_src,\n",
    "\n",
    "    x_det,\n",
    "    y_det,\n",
    "    z_det,\n",
    "\n",
    "    x_obj,\n",
    "    y_obj,\n",
    "    z_obj,\n",
    "\n",
    "    alpha_x, \n",
    "    alpha_y,\n",
    "    alpha_z,\n",
    "\n",
    "#     x_rot,\n",
    "#     y_rot,\n",
    "#     z_rot,\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "eda086c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Create an OpenGL context: 512x512\n"
     ]
    }
   ],
   "source": [
    "# Initialise gVXR using our JSON file\n",
    "json2gvxr.initGVXR(\"simulation1.json\", \"OPENGL\", 3, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "19579e84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set up the beam\n",
      "\tSource position: [0.0, 0.0, 0.0, 'mm']\n",
      "\tSource shape: PointSource\n"
     ]
    }
   ],
   "source": [
    "# Load our source properties\n",
    "json2gvxr.initSourceGeometry(\"simulation1.json\")\n",
    "json2gvxr.initSpectrum(verbose=0)\n",
    "\n",
    "gvxr.setSourcePosition(x_src, y_src, z_src, \"mm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "54551648",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAverageEnergy(k, f):\n",
    "    \n",
    "    avg = 0\n",
    "    for energy, count in zip(k, f):\n",
    "        avg += energy * count\n",
    "        \n",
    "    return avg / np.sum(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0dea486d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set up the detector\n",
      "\tDetector position: [0.0, 0.0, 0.0, 'mm']\n",
      "\tDetector up vector: [0, -1, 0]\n",
      "\tNumber of pixels: [640, 480]\n",
      "\tDetector number of pixels: [640, 480]\n",
      "\tPixel spacing: [0.23, 0.23, 'mm']\n"
     ]
    }
   ],
   "source": [
    "# Load our detector\n",
    "json2gvxr.initDetector(\"simulation1.json\")\n",
    "    # \"NumberOfPixels\": [2880, 2880],\n",
    "\n",
    "gvxr.setDetectorPosition(x_det, y_det, z_det, \"mm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9a6852a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load our samples\n",
    "json2gvxr.initSamples(\"simulation1.json\", verbose=0)\n",
    "\n",
    "gvxr.setDensity(\"cuboid\", Ti90Al10, \"g/cm3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "79c51d65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGFCAYAAABQVuEuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0vUlEQVR4nO3df1RVdb7/8deR3xCQgnIg0cHJWxrWbQ79QHO0q+Ky1GlqTTb9vCttIs1EdDWpNf5YN2l1G2NMhathVlY4d8y5toab4JoiG+2WBCs1l9WShIxzSVNAU0DY3z/8cq6bX54D+3Dg8HystdeKvT9nn8/nGGe/eH/2D5thGIYAAAAsNsDXHQAAAP6JkAEAALyCkAEAALyCkAEAALyCkAEAALyCkAEAALyCkAEAALwi0NcdAADAH5w/f14NDQ3d3k9wcLBCQ0Mt6JHvETIAAOim8+fPKykpSU6ns9v7stvtKi8v94ugQcgAAKCbGhoa5HQ6VVlZqaioqC7vp7a2VomJiWpoaCBkAACA/xMZGanIyMguv97fnvRByAAAwCKGYXQrKPhbyODqEgAA4BVUMgAAsAiVDDNCBgAAFiFkmDFdAgAAvIJKBgAAFqGSYUbIAADAIoQMM6ZLAACAV1DJAADAIlQyzAgZAABYhJBhRsgAAMAihAwzzskAAABeQSUDAACLUMkwI2QAAGARQoYZ0yUAAMArqGQAAGARKhlmhAwAACxCyDBjugQAAHgFlQwAACxCJcOMkAEAgEUIGWZMlwAAAK+gkgEAgEWoZJgRMgAAsAghw4yQAQCARQgZZpyTAQAAvIJKBgAAFvK3akR3EDIAALAI0yVmTJcAAACvIGQAAGCRlkpGd5au2LBhg5KSkhQaGiqHw6E9e/Z02r64uFgOh0OhoaEaMWKEcnNz27TJzs7WNddco7CwMCUmJmrhwoU6f/68R/0iZAAAYBFfhIxt27YpIyNDy5YtU2lpqcaPH69p06apoqKi3fbl5eW64447NH78eJWWlmrp0qV66qmntH37dlebt956S88884yWL1+uw4cPKy8vT9u2bdOSJUs86ptPQ4anyQsAAJitWbNGs2fP1pw5czRq1ChlZ2crMTFROTk57bbPzc3VsGHDlJ2drVGjRmnOnDl69NFH9dJLL7na7Nu3T+PGjdP999+vn/3sZ0pLS9Nvf/tb7d+/36O++SxkeJq8AADo7ayqZNTW1pqW+vr6dt+voaFBJSUlSktLM61PS0vT3r17233Nvn372rSfOnWq9u/fr8bGRknSbbfdppKSEn366aeSpKNHj6qgoEB33nmnR5+Hz64uuTR5SRfnfnbt2qWcnBxlZWV1+trm5mZ9//33ioyMlM1m64nuAgD6KMMwVFdXp4SEBA0Y4N2/ra26uiQxMdG0fvny5VqxYkWb9idOnFBTU5Pi4uJM6+Pi4uR0Ott9D6fT2W77Cxcu6MSJE4qPj9d9992nH374QbfddpsMw9CFCxf0xBNP6JlnnvFoPD4JGS3Jq3VnO0pe9fX1phR3/PhxjR492uv9BAD4j8rKSg0dOtSr72FVyKisrFRUVJRrfUhISKeva/0Ht2EYnf4R3l77S9d/+OGHev7557Vhwwbdcsst+uabb7RgwQLFx8frueeec3s8PgkZniavrKwsrVy5ss361v8IAAC0Vltbq8TEREVGRvq6K26Liopy6/gWGxurgICANsfO6urqNsfYFna7vd32gYGBiomJkSQ999xzeuihh1yzDWPGjNHZs2f1u9/9TsuWLXO7IuTTm3G5m7yWLFmizMxM188t/8O4+48AAEBPTK/39M24goOD5XA4VFRUpF//+teu9UVFRfrVr37V7mtSU1P13nvvmdYVFhYqJSVFQUFBkqSffvqpTZAICAjweHw+CRmeJq+QkJDLlooAAPA1X9zxMzMzUw899JBSUlKUmpqqjRs3qqKiQunp6ZIu/qF+/PhxvfHGG5Kk9PR0rVu3TpmZmXrssce0b98+5eXl6Z133nHtc8aMGVqzZo1uvPFG13TJc889p5kzZyogIMDtvvkkZHQleQEAgLZmzZqlkydPatWqVaqqqlJycrIKCgo0fPhwSVJVVZXpys2kpCQVFBRo4cKFWr9+vRISErR27Vrdc889rjbPPvusbDabnn32WR0/flyDBw/WjBkz9Pzzz3vUN5vhoxulb9u2TQ899JByc3NdyWvTpk06dOiQ64PpSG1traKjo1VTU8N0CQCgUz1xzGh5jy+//LJb537U1dVp9OjRfnN889k5GZdLXgAA9DU8IM3Mpyd+zp07V3PnzvVlFwAAgJfwqHcAACxCJcOMkAEAgEUIGWY8hRUAAHgFlQwAACxCJcOMkAEAgEUIGWaEDAAALELIMOOcDAAA4BVUMgAAsAiVDDNCBgAAFiFkmDFdAgAAvIJKBgAAFvK3akR3EDIAALAI0yVmTJcAAACvoJIBAIBFqGSYETIAALAIIcOM6RIAAOAVVDIAALAIlQwzQgYAABYhZJgRMgAAsAghw4xzMgAAgFdQyQAAwCJUMswIGQAAWISQYcZ0CQAA8AoqGQAAWIRKhhkhAwAAixAyzJguAQAAXkElAwAAi1DJMCNkAABgEUKGGdMlAADAK6hkAABgESoZZoQMAAAsQsgwI2QAAGARQoYZ52QAAACvoJIBAIBFqGSYETIAALAIIcOM6RIAAOAVVDIAALCQv1UjuoOQAQCARZguMWO6BAAAeAWVDAAALEIlw4yQAQCARQgZZkyXAAAAr6CSAQCARahkmBEyAACwCCHDjJABAIBFCBlmnJMBAAC8gkoGAAAWoZJhRsgAAMAihAwzpksAAIBXUMkAAMAiVDLMCBkAAFiEkGHm8XTJRx99pBkzZighIUE2m01//etfTdsNw9CKFSuUkJCgsLAwTZw4UYcOHTK1qa+v1/z58xUbG6uIiAjNnDlT3333XbcGAgAAehePQ8bZs2d1ww03aN26de1uf/HFF7VmzRqtW7dOn332mex2u6ZMmaK6ujpXm4yMDO3YsUP5+fn6+OOPdebMGU2fPl1NTU1dHwkAAD7WUsnozuJPPJ4umTZtmqZNm9buNsMwlJ2drWXLlunuu++WJL3++uuKi4vT22+/rccff1w1NTXKy8vTm2++qcmTJ0uStm7dqsTERO3evVtTp07txnAAAPAdpkvMLL26pLy8XE6nU2lpaa51ISEhmjBhgvbu3StJKikpUWNjo6lNQkKCkpOTXW1aq6+vV21trWkBAAC9m6Uhw+l0SpLi4uJM6+Pi4lzbnE6ngoODNXDgwA7btJaVlaXo6GjXkpiYaGW3AQCwhK+mSzZs2KCkpCSFhobK4XBoz549nbYvLi6Ww+FQaGioRowYodzcXNP2iRMnymaztVnuvPNOj/rllftk2Gw208+GYbRZ11pnbZYsWaKamhrXUllZaVlfAQCwii9CxrZt25SRkaFly5aptLRU48eP17Rp01RRUdFu+/Lyct1xxx0aP368SktLtXTpUj311FPavn27q827776rqqoq13Lw4EEFBAToN7/5jUd9szRk2O12SWpTkaiurnZVN+x2uxoaGnTq1KkO27QWEhKiqKgo0wIAQG/ji5CxZs0azZ49W3PmzNGoUaOUnZ2txMRE5eTktNs+NzdXw4YNU3Z2tkaNGqU5c+bo0Ucf1UsvveRqM2jQINntdtdSVFSk8PBw34aMpKQkV2daNDQ0qLi4WGPHjpUkORwOBQUFmdq0pKSWNgAA9Getz0Osr69vt11DQ4NKSkpM5zlKUlpaWofnOe7bt69N+6lTp2r//v1qbGxs9zV5eXm67777FBER4dE4PL665MyZM/rmm29cP5eXl6usrEyDBg3SsGHDlJGRodWrV2vkyJEaOXKkVq9erfDwcN1///2SpOjoaM2ePVuLFi1STEyMBg0apMWLF2vMmDGuq00AAOiLrLq6pPW5h8uXL9eKFSvatD9x4oSampo6PReyNafT2W77Cxcu6MSJE4qPjzdt+/TTT3Xw4EHl5eV5OhzPQ8b+/ft1++23u37OzMyUJD3yyCPasmWLnn76aZ07d05z587VqVOndMstt6iwsFCRkZGu17z88ssKDAzUvffeq3PnzmnSpEnasmWLAgICPB4AAAC9hVUho7Ky0nRqQEhISKev8/RcyPbat7deuljFSE5O1s0339x559vhcciYOHFipx+gzWbTihUr2k1cLUJDQ/XKK6/olVde8fTtAQDwe+6efxgbG6uAgIBOz4VszW63t9s+MDBQMTExpvU//fST8vPztWrVKg9HcBFPYQUAwEI9edJncHCwHA6H6TxHSSoqKurwPMfU1NQ27QsLC5WSkqKgoCDT+j//+c+qr6/Xgw8+6HHfJEIGAACW8cXVJZmZmXr11Ve1efNmHT58WAsXLlRFRYXS09MlXbwNxMMPP+xqn56ermPHjikzM1OHDx/W5s2blZeXp8WLF7fZd15enu666642FQ538RRWAAD6sFmzZunkyZNatWqVqqqqlJycrIKCAg0fPlzSxSs4L71nRlJSkgoKCrRw4UKtX79eCQkJWrt2re655x7Tfr/66it9/PHHKiws7HLfbEYfvFF6bW2toqOjVVNTwz0zAACd6oljRst7/PnPf1Z4eHiX9/PTTz/p3nvv9ZvjG5UMAAAswgPSzDgnAwAAeAWVDAAALEIlw4yQAQCARQgZZoQMAAAsQsgw45wMAADgFVQyAACwCJUMM0IGAAAWIWSYMV0CAAC8gkoGAAAWoZJhRsgAAMAihAwzpksAAIBXUMkAAMAiVDLMCBkAAFiEkGHGdAkAAPAKKhkAAFiESoYZIQMAAIsQMswIGQAAWISQYcY5GQAAwCuoZAAAYBEqGWaEDAAALELIMGO6BAAAeAWVDAAALEIlw4yQAQCAhfwtKHQH0yUAAMArqGQAAGARpkvMCBkAAFiEkGHGdAkAAPAKKhkAAFiESoYZIQMAAIsQMswIGQAAWISQYcY5GQAAwCuoZAAAYBEqGWaEDAAALELIMGO6BAAAeAWVDAAALEIlw4yQgX6jubnZ736B+yObzaYBAyjConciZJgRMtBvnD17Vk6nU/X19T55/+4cGG02m4U9sV53vxibm5vdahcSEqL4+HhdccUV3Xo/AD2DkIF+wzAM1dbW6vvvv1ddXZ3bBzbJ84DQWShwJzB4Giq8GUK6EiAu95qOtnf2usjISF111VWKi4vzuD9AT6GSYUbIQL/T2NioH374QadOnVJDQ4Mk90OElQHhcu2s2k93uPuF11k7TwJF63XNzc0KDg7WwIEDFRoa6lZfAF8iZJgRMtBvnTt3Tj/88IPOnz/f7i92Vw/elwss3goPVoUNbwYLd6pHLa+12WwKDQ3V4MGDNXDgQLf6BKB3IWSg3zIMQ01NTTp9+rTq6up04cKFDgNCRwdwbwURK97DapcLH+5OP7XeT3vVi8DAQEVGRmrw4MGedRLwMSoZZoQM9FstB2/DMHT+/HnV1taqoaHB9Evu7gG+s9DgjXMwfMkbFYxLqxfBwcG68sor2+yrL31G6L8IGWaEDPR7LQev5uZmnTlzRufOnVNzc3O7BzV313W23t3tVr2mI55+mVl5Mmd71YwBAwYoLCxMV155paSLYyVYoK8hZJgRMgCZKxGNjY2uE0JbWBUuentVw50vOG+d5BkcHKywsLDLvj+AvoOQgX6tN/6lbBiGT/rlb39BAb5AJcOMkAFcRn+pYlx6jsrl2nXUpqNt7a3vbD9AX0XIMPPoDkNZWVm66aabFBkZqSFDhuiuu+7SkSNHTG0Mw9CKFSuUkJCgsLAwTZw4UYcOHTK1qa+v1/z58xUbG6uIiAjNnDlT3333XfdHA/hQdwJGbzr/wJ2+dDdMWd0fAL2TRyGjuLhY8+bN0yeffKKioiJduHBBaWlpOnv2rKvNiy++qDVr1mjdunX67LPPZLfbNWXKFNXV1bnaZGRkaMeOHcrPz9fHH3+sM2fOaPr06WpqarJuZIAFunvy5+W2ubPdHS0H4ksXK/bZ1e1d+dwIE/AHLZWM7iz+xKPpkvfff9/082uvvaYhQ4aopKREv/zlL2UYhrKzs7Vs2TLdfffdkqTXX39dcXFxevvtt/X444+rpqZGeXl5evPNNzV58mRJ0tatW5WYmKjdu3dr6tSpFg0N8NzlDnRW3yDLm1eYtNfO0y+wy01peDrl4W77lr4TOtAX+VtQ6I5uPcqwpqZGkjRo0CBJUnl5uZxOp9LS0lxtQkJCNGHCBO3du1eSVFJSosbGRlObhIQEJScnu9q0Vl9fr9raWtMCeJMnf1V39fwMTw6gVlUourKfrgYvTz4/wgTgn7ocMgzDUGZmpm677TYlJydLkpxOpyS1eYBRXFyca5vT6XQ9i6CjNq1lZWUpOjratSQmJna120CX9eQJoJe29dYBuCthw9Ntnk4tAX0d0yVmXQ4ZTz75pL744gu98847bba1/hJx55K8ztosWbJENTU1rqWysrKr3UY/Z/XtwT3dV2/8696Kqs3ltgH9BSHDrEshY/78+dq5c6c++OADDR061LXebrdLUpuKRHV1tau6Ybfb1dDQoFOnTnXYprWQkBBFRUWZFqAnWfEXeW94qqoV72vFuAkkQP/gUcgwDENPPvmk3n33Xf39739XUlKSaXtSUpLsdruKiopc6xoaGlRcXKyxY8dKkhwOh4KCgkxtqqqqdPDgQVcboLfzVkXE1wdfb/XfqnH5+vMBLodKhplHV5fMmzdPb7/9tv7rv/5LkZGRropFdHS0wsLCZLPZlJGRodWrV2vkyJEaOXKkVq9erfDwcN1///2utrNnz9aiRYsUExOjQYMGafHixRozZozrahOgN+mJqYTedPB056ZcXbkhl7vtPH0SLtCbcDMuM48qGTk5OaqpqdHEiRMVHx/vWrZt2+Zq8/TTTysjI0Nz585VSkqKjh8/rsLCQkVGRrravPzyy7rrrrt07733aty4cQoPD9d7772ngIAA60YGXEZXD1q94TLWnmB1v3vrOAEr+aqSsWHDBiUlJSk0NFQOh0N79uzptH1xcbEcDodCQ0M1YsQI5ebmtmlz+vRpzZs3T/Hx8QoNDdWoUaNUUFDgUb88qmS4+9fJihUrtGLFig7bhIaG6pVXXtErr7ziydsDlmh9Ayh321rRzqrXefL67vxl1NVbf3elmtHRCaiEE6Bz27ZtU0ZGhjZs2KBx48bpP/7jPzRt2jR9+eWXGjZsWJv25eXluuOOO/TYY49p69at+sc//qG5c+dq8ODBuueeeyRdPNVhypQpGjJkiP7yl79o6NChqqysNBUM3MGzS9Av9dRVEt44QHangmB1KbarIcSd/XY0bQL0Zr6YLlmzZo1mz56tOXPmSJKys7O1a9cu5eTkKCsrq0373NxcDRs2TNnZ2ZKkUaNGaf/+/XrppZdcIWPz5s368ccftXfvXgUFBUmShg8f7nHf+C0GLnG5KodV98pwZ3t77a26GZenr+nK9u5+fp62AXoDq6ZLWt+Asr6+vt33a2hoUElJiekGl5KUlpbW4Q0u9+3b16b91KlTtX//fjU2NkqSdu7cqdTUVM2bN09xcXFKTk7W6tWrPX78ByED+P8u/cvZ2wc1X5/P4Ov3d3f/hAv0V4mJiaabULZXkZCkEydOqKmpqdObYLbmdDrbbX/hwgWdOHFCknT06FH95S9/UVNTkwoKCvTss8/qj3/8o55//nmPxsF0CeAmq+6O6a337Oq+rZju8GTaxFtTLEBvYNV0SWVlpemeUCEhIZ2+ztObYLbX/tL1zc3NGjJkiDZu3KiAgAA5HA59//33+vd//3f94Q9/cHs8hAygG7oSAnrbTbm6e8mp1a8B+jKrQoa7N56MjY1VQEBApzfBbM1ut7fbPjAwUDExMZKk+Ph4BQUFma76HDVqlJxOpxoaGhQcHOzWeJguAXqhnp4mYFoC6JuCg4PlcDhMN7iUpKKiog5vcJmamtqmfWFhoVJSUlwneY4bN07ffPONmpubXW2++uorxcfHux0wJEIGYDl/fr6HP48NsIIv7pORmZmpV199VZs3b9bhw4e1cOFCVVRUKD09XdLF5389/PDDrvbp6ek6duyYMjMzdfjwYW3evFl5eXlavHixq80TTzyhkydPasGCBfrqq6/0t7/9TatXr9a8efM86hvTJYC6dgWFN5+O6gvevBy19X6ZRoG/8sUlrLNmzdLJkye1atUqVVVVKTk5WQUFBa5LTquqqlRRUeFqn5SUpIKCAi1cuFDr169XQkKC1q5d67p8Vbp44mlhYaEWLlyo66+/XldddZUWLFig3//+9x71jZAB9JC+cKmmOwd/AgLQ+8ydO1dz585td9uWLVvarJswYYI+//zzTveZmpqqTz75pFv9ImSg37Pypk++Dgm+RgBBf8ezS8wIGQAAWISQYUbIALqgN9wu3Fu8UY2gwoH+gpBhxtUl6Ld68qDeWwKEVfxtPAC8g0oG+iVf3L0TZlQ34I+oZJgRMtCvtA4MBAgAViJkmDFdgn6nNwaL3tan3tYfAH0TlQwAACxCJcOMkAEAgIX8LSh0B9MlAADAK6hkAL2AYRi96jwI/hIDuobpEjNCBgAAFiFkmDFdAlyGv/3S9xaefK78GwB9E5UMoAf0tumQ7uqpgz7hAn0NlQwzQgbQBd4IDb0liHjjS87fvjiBjhAyzAgZAHo1f/vShX8jZJhxTgZgIX/7gvCUleNvbm5us+/+/vkCfQ2VDKCHuDMd4uspE3cO4hzogY5RyTAjZABuaO/g761A4Kug4a0vt/b2629fpEALQoYZ0yWAxTr7kujrXyA9Nba+/jkBuIhKBtAL9XQ1g4M6YA0qGWaEDPRr3f2F7koYcPc1PRU03P0MuvJZ+dsXJnA5hAwzpkvQ73T1l9gXd6j05hUVVu6bu3cCaA+VDKAd3q4ieLp/q/vj6YHe28GA4AF/QSXDjJCBfqXlF/jSL4Lm5uYOr4Bw54qSjgLA5YJBV4KGpG6FDW9MeXS03d2rSi5dZxhGm/tjAH0JIcOMkAG0YvVVEt64/XgLd8/t8JaeuqLE3754gf6CkAFcovXBzNsnaXY3hPTWaQxvnkwK9GZUMswIGcD/526Z3urzKXx9l8+OdHWaxKr2XX0N4EuEDDNCBnAZVlQz3AkaUvfOt7BKd28t7q0qhr99+cI/ETLMuIQV6AJPTna0Yr89xVv950sX6J+oZABusOLZJb3tJlztva+V7TprT3CAv6KSYUbIACzWnWmTS9tJPTN9YtWNtHh2CUDIaI3pEsBNnvxFbtXBuDfd8bMrY6KKAfRvVDLQr116c6721ndneuNyFQ1399+6f335Zlzd7RMBBb0dlQwzQgbggd5w+Wp7X0IdhZzu6unLWP3tCxb9E/8f/x+mS9BvdVTKb32ba3df58n7eLLdHS39vnSxYp9d3d7VW4oD8C9UMtAveeOAxn0y3H894K+YLjEjZACX4cnlq929ssSXYcOKcyg4ART9HSHDjOkS9Du++iV2Nzj4qprR2/sHoO+hkgG0EhQUpMBA869GewfWjg62lzsIe3KQ9vYBvSfukeFuFcNbJ68CPYlKhhkhA/3apV8IAwYM0BVXXKGIiIjLHtwHDGhbBOzsNX35r//unADqyUPnpIufU8tn629ftugfCBlmhAz0ezabTcHBwYqNjTWta6+dO9oLIJ68vqvte+KKEsn94ODuvttbP2DAgD4dzNB/ETLMPAoZOTk5ysnJ0bfffitJuu666/SHP/xB06ZNk3Txw1m5cqU2btyoU6dO6ZZbbtH69et13XXXufZRX1+vxYsX65133tG5c+c0adIkbdiwQUOHDrVuVICbAgICdOWVVyo6Otqtg5onFQx39tdTB1Krb3h1uTZWBBGbzaaAgIAu7QdA7+BRyBg6dKheeOEFXX311ZKk119/Xb/61a9UWlqq6667Ti+++KLWrFmjLVu26J/+6Z/0b//2b5oyZYqOHDmiyMhISVJGRobee+895efnKyYmRosWLdL06dNVUlLCFwp6VFhYmBISEkzrWh/0PQ0QvSlYtGbVHTitCCGetGt9fgzQm1HJMPPot3fGjBmmn59//nnl5OTok08+0ejRo5Wdna1ly5bp7rvvlnQxhMTFxentt9/W448/rpqaGuXl5enNN9/U5MmTJUlbt25VYmKidu/eralTp1o0LKB9hmEoKChIQ4YM0eDBg13ruxIa+vLVGN68XNXT93CHzWZTUFCQZfsDvIWQYdblPxGampr0n//5nzp79qxSU1NVXl4up9OptLQ0V5uQkBBNmDBBe/fu1eOPP66SkhI1Njaa2iQkJCg5OVl79+7tMGTU19ervr7e9XNtbW1Xu41+bMCAAYqKilJYWJirnN/R+RMd8UZgsHKf3nyYmie6Ol1yOSEhIR7/mwHwHY9DxoEDB5Samqrz58/riiuu0I4dOzR69Gjt3btXkhQXF2dqHxcXp2PHjkmSnE6ngoODNXDgwDZtnE5nh++ZlZWllStXetpVwCQiIkIjRozwdTfQTYQM9GZUMsw8DhnXXHONysrKdPr0aW3fvl2PPPKIiouLXdvdvTOiJ22WLFmizMxM18+1tbVKTEz0tOvo5ziREIC3ETLMPP6TIDg4WFdffbVSUlKUlZWlG264QX/6059kt9slqU1Forq62lXdsNvtamho0KlTpzps056QkBBFRUWZFgAA0Lt1u+5oGIbq6+uVlJQku92uoqIi17aGhgYVFxdr7NixkiSHw6GgoCBTm6qqKh08eNDVBgCAvqq9JyJ7uvgTj6ZLli5dqmnTpikxMVF1dXXKz8/Xhx9+qPfff182m00ZGRlavXq1Ro4cqZEjR2r16tUKDw/X/fffL0mKjo7W7NmztWjRIsXExGjQoEFavHixxowZ47raBACAvorpEjOPQsb//u//6qGHHlJVVZWio6N1/fXX6/3339eUKVMkSU8//bTOnTunuXPnum7GVVhY6LpHhiS9/PLLCgwM1L333uu6GdeWLVuYKwcAwM/YjD4Ym2praxUdHa2amhrOzwAAdKonjhkt7/HAAw8oODi4y/tpaGjQW2+95TfHN26lBwCARZguMSNkAABgEUKGGXe1AQAAXkElAwAAi1DJMKOSAQCARXx1n4wNGzYoKSlJoaGhcjgc2rNnT6fti4uL5XA4FBoaqhEjRig3N9e0fcuWLbLZbG2W8+fPe9QvQgYAAH3Ytm3blJGRoWXLlqm0tFTjx4/XtGnTVFFR0W778vJy3XHHHRo/frxKS0u1dOlSPfXUU9q+fbupXVRUlKqqqkxLaGioR31jugQAAIv4YrpkzZo1mj17tubMmSNJys7O1q5du5STk6OsrKw27XNzczVs2DBlZ2dLkkaNGqX9+/frpZde0j333ONqZ7PZXI8M6SoqGQAAWMiKqZLa2lrTUl9f3+57NTQ0qKSkRGlpaab1aWlprqejt7Zv37427adOnar9+/ersbHRte7MmTMaPny4hg4dqunTp6u0tNTjz4KQAQBAL5OYmKjo6GjX0l5FQpJOnDihpqamNg8ZjYuLa/PA0hZOp7Pd9hcuXNCJEyckSddee622bNminTt36p133lFoaKjGjRunr7/+2qNxMF0CAIBFrJouqaysNN3xMyQkpNPX2Wy2Nvtpve5y7S9df+utt+rWW291bR83bpx+8Ytf6JVXXtHatWvdGMlFhAwAACxiVciIiopy67bisbGxCggIaFO1qK6ublOtaGG329ttHxgYqJiYmHZfM2DAAN10000eVzKYLgEAoI8KDg6Ww+FQUVGRaX1RUZHGjh3b7mtSU1PbtC8sLFRKSoqCgoLafY1hGCorK1N8fLxH/SNkAABgEV/cJyMzM1OvvvqqNm/erMOHD2vhwoWqqKhQenq6JGnJkiV6+OGHXe3T09N17NgxZWZm6vDhw9q8ebPy8vK0ePFiV5uVK1dq165dOnr0qMrKyjR79myVlZW59ukupksAALCILy5hnTVrlk6ePKlVq1apqqpKycnJKigo0PDhwyVJVVVVpntmJCUlqaCgQAsXLtT69euVkJCgtWvXmi5fPX36tH73u9/J6XQqOjpaN954oz766CPdfPPNHvWNR70DAPxaTz7q/de//nWHUw7uaGxs1I4dO/zm+MZ0CQAA8AqmSwAAsAgPSDMjZAAAYBFChhnTJQAAwCuoZAAAYBEqGWaEDAAALELIMGO6BAAAeAWVDAAALEIlw4yQAQCARQgZZkyXAAAAr6CSAQCARahkmBEyAACwCCHDjJABAIBFCBlmnJMBAAC8gkoGAAAWoZJhRsgAAMAihAwzpksAAIBXUMkAAMAiVDLMCBkAAFiEkGHGdAkAAPAKKhkAAFjI36oR3UHIAADAIkyXmBEyAACwCCHDjHMyAACAV1DJAADAIlQyzAgZAABYhJBhxnQJAADwCioZAABYhEqGGSEDAACLEDLMmC4BAABeQSUDAACLUMkwI2QAAGARQoYZ0yUAAMArqGQAAGARKhlmhAwAACxCyDAjZAAAYBFChlm3zsnIysqSzWZTRkaGa51hGFqxYoUSEhIUFhamiRMn6tChQ6bX1dfXa/78+YqNjVVERIRmzpyp7777rjtdAQAAvUyXQ8Znn32mjRs36vrrrzetf/HFF7VmzRqtW7dOn332mex2u6ZMmaK6ujpXm4yMDO3YsUP5+fn6+OOPdebMGU2fPl1NTU1dHwkAAD7WUsnozuJPuhQyzpw5owceeECbNm3SwIEDXesNw1B2draWLVumu+++W8nJyXr99df1008/6e2335Yk1dTUKC8vT3/84x81efJk3Xjjjdq6dasOHDig3bt3WzMqAAB8gJBh1qWQMW/ePN15552aPHmyaX15ebmcTqfS0tJc60JCQjRhwgTt3btXklRSUqLGxkZTm4SEBCUnJ7vatFZfX6/a2lrTAgAAejePT/zMz8/X559/rs8++6zNNqfTKUmKi4szrY+Li9OxY8dcbYKDg00VkJY2La9vLSsrSytXrvS0qwAA9ChO/DTzqJJRWVmpBQsWaOvWrQoNDe2wnc1mM/1sGEabda111mbJkiWqqalxLZWVlZ50GwCAHsF0iZlHIaOkpETV1dVyOBwKDAxUYGCgiouLtXbtWgUGBroqGK0rEtXV1a5tdrtdDQ0NOnXqVIdtWgsJCVFUVJRpAQAAvZtHIWPSpEk6cOCAysrKXEtKSooeeOABlZWVacSIEbLb7SoqKnK9pqGhQcXFxRo7dqwkyeFwKCgoyNSmqqpKBw8edLUBAKAvopJh5tE5GZGRkUpOTjati4iIUExMjGt9RkaGVq9erZEjR2rkyJFavXq1wsPDdf/990uSoqOjNXv2bC1atEgxMTEaNGiQFi9erDFjxrQ5kRQAgL6EczLMLL/j59NPP61z585p7ty5OnXqlG655RYVFhYqMjLS1ebll19WYGCg7r33Xp07d06TJk3Sli1bFBAQYHV3AACAj9iMPhibamtrFR0drZqaGs7PAAB0qieOGS3v4XA4uvUHc1NTk0pKSvzm+MazSwAAsAjTJWaEDAAALELIMOvWA9IAAAA6QiUDAACLUMkwI2QAAGARQoYZ0yUAAMArqGQAAGARKhlmhAwAACxCyDBjugQAAHgFIQMAAIv46gFpGzZsUFJSkkJDQ+VwOLRnz55O2xcXF8vhcCg0NFQjRoxQbm5uh23z8/Nls9l01113edwvQgYAABbxRcjYtm2bMjIytGzZMpWWlmr8+PGaNm2aKioq2m1fXl6uO+64Q+PHj1dpaamWLl2qp556Stu3b2/T9tixY1q8eLHGjx/vcb8kQgYAAH3amjVrNHv2bM2ZM0ejRo1Sdna2EhMTlZOT02773NxcDRs2TNnZ2Ro1apTmzJmjRx99VC+99JKpXVNTkx544AGtXLlSI0aM6FLfCBkAAFjEqkpGbW2taamvr2/3/RoaGlRSUqK0tDTT+rS0NO3du7fd1+zbt69N+6lTp2r//v1qbGx0rVu1apUGDx6s2bNnd/nzIGQAAGARq0JGYmKioqOjXUtWVla773fixAk1NTUpLi7OtD4uLk5Op7Pd1zidznbbX7hwQSdOnJAk/eMf/1BeXp42bdrUrc+DS1gBALCIVZewVlZWmh71HhIS0unrbDZbm/20Xne59i3r6+rq9OCDD2rTpk2KjY31qP+tETIAAOhloqKiTCGjI7GxsQoICGhTtaiurm5TrWhht9vbbR8YGKiYmBgdOnRI3377rWbMmOHa3tzcLEkKDAzUkSNH9POf/9ytcTBdAgCARXr66pLg4GA5HA4VFRWZ1hcVFWns2LHtviY1NbVN+8LCQqWkpCgoKEjXXnutDhw4oLKyMtcyc+ZM3X777SorK1NiYqLb/aOSAQCARXxxx8/MzEw99NBDSklJUWpqqjZu3KiKigqlp6dLkpYsWaLjx4/rjTfekCSlp6dr3bp1yszM1GOPPaZ9+/YpLy9P77zzjiQpNDRUycnJpve48sorJanN+sshZAAA0IfNmjVLJ0+e1KpVq1RVVaXk5GQVFBRo+PDhkqSqqirTPTOSkpJUUFCghQsXav369UpISNDatWt1zz33WN43m9EHb5ReW1ur6Oho1dTUuDVnBQDov3rimNHyHiNHjlRAQECX99PU1KSvv/7ab45vVDIAALAID0gz48RPAADgFVQyAACwCJUMM0IGAAAWIWSYMV0CAAC8gkoGAAAW8rdqRHcQMgAAsAjTJWaEDAAALELIMOOcDAAA4BVUMgAAsAiVDDNCBgAAFiFkmDFdAgAAvIJKBgAAFqGSYUbIAADAIoQMM6ZLAACAV1DJAADAIlQyzAgZAABYhJBhxnQJAADwCioZAABYhEqGGSEDAACLEDLMCBkAAFiEkGHWJ0NGyz9CbW2tj3sCAOjtWo4V/nYA7wv6ZMioq6uTJCUmJvq4JwCAvqKurk7R0dFefQ8qGWZ9MmQkJCToyy+/1OjRo1VZWamoqChfd6nH1NbWKjExkXH3E4ybcfcH3h63YRiqq6tTQkKC5ftu770IGf+nT4aMAQMG6KqrrpIkRUVF9atfxhaMu39h3P0L47aetysYaF+fDBkAAPRGVDLMCBkAAFiEkGHWZ+/4GRISouXLlyskJMTXXelRjJtx9weMm3HDP9gMf4tNAAD0sNraWkVHR2vQoEEaMKDrf783Nzfrxx9/VE1NjV+cl8N0CQAAFmG6xKzPTpcAAIDejUoGAAAW8rdqRHcQMgAAsEh3A4a/BZQ+OV2yYcMGJSUlKTQ0VA6HQ3v27PF1l7rlo48+0owZM5SQkCCbzaa//vWvpu2GYWjFihVKSEhQWFiYJk6cqEOHDpna1NfXa/78+YqNjVVERIRmzpyp7777rgdH4ZmsrCzddNNNioyM1JAhQ3TXXXfpyJEjpjb+OO6cnBxdf/31rpsOpaam6r//+79d2/1xzO3JysqSzWZTRkaGa50/jn3FihWy2WymxW63u7b745hbHD9+XA8++KBiYmIUHh6uf/7nf1ZJSYlru7+OveWcjO4sfsXoY/Lz842goCBj06ZNxpdffmksWLDAiIiIMI4dO+brrnVZQUGBsWzZMmP79u2GJGPHjh2m7S+88IIRGRlpbN++3Thw4IAxa9YsIz4+3qitrXW1SU9PN6666iqjqKjI+Pzzz43bb7/duOGGG4wLFy708GjcM3XqVOO1114zDh48aJSVlRl33nmnMWzYMOPMmTOuNv447p07dxp/+9vfjCNHjhhHjhwxli5dagQFBRkHDx40DMM/x9zap59+avzsZz8zrr/+emPBggWu9f449uXLlxvXXXedUVVV5Vqqq6td2/1xzIZhGD/++KMxfPhw41//9V+N//mf/zHKy8uN3bt3G998842rjb+NvaamxpBkREdHG1deeWWXl+joaEOSUVNT4+shWaLPhYybb77ZSE9PN6279tprjWeeecZHPbJW65DR3Nxs2O1244UXXnCtO3/+vBEdHW3k5uYahmEYp0+fNoKCgoz8/HxXm+PHjxsDBgww3n///R7re3dUV1cbkozi4mLDMPrPuA3DMAYOHGi8+uqr/WLMdXV1xsiRI42ioiJjwoQJrpDhr2Nfvny5ccMNN7S7zV/HbBiG8fvf/9647bbbOtzuj2NvCRlRUVFGdHR0l5eoqCi/Chl9arqkoaFBJSUlSktLM61PS0vT3r17fdQr7yovL5fT6TSNOSQkRBMmTHCNuaSkRI2NjaY2CQkJSk5O7jOfS01NjSRp0KBBkvrHuJuampSfn6+zZ88qNTW1X4x53rx5uvPOOzV58mTTen8e+9dff62EhAQlJSXpvvvu09GjRyX595h37typlJQU/eY3v9GQIUN04403atOmTa7t/jx2g+kSkz4VMk6cOKGmpibFxcWZ1sfFxcnpdPqoV97VMq7Oxux0OhUcHKyBAwd22KY3MwxDmZmZuu2225ScnCzJv8d94MABXXHFFQoJCVF6erp27Nih0aNH+/WYJSk/P1+ff/65srKy2mzz17HfcssteuONN7Rr1y5t2rRJTqdTY8eO1cmTJ/12zJJ09OhR5eTkaOTIkdq1a5fS09P11FNP6Y033pDkv//eaKtPXl1is9lMPxuG0Wadv+nKmPvK5/Lkk0/qiy++0Mcff9xmmz+O+5prrlFZWZlOnz6t7du365FHHlFxcbFruz+OubKyUgsWLFBhYaFCQ0M7bOdvY582bZrrv8eMGaPU1FT9/Oc/1+uvv65bb71Vkv+NWbp418qUlBStXr1aknTjjTfq0KFDysnJ0cMPP+xq549j724lgkqGD8XGxiogIKBNiq2urm6TiP1Fy5nonY3ZbreroaFBp06d6rBNbzV//nzt3LlTH3zwgYYOHepa78/jDg4O1tVXX62UlBRlZWXphhtu0J/+9Ce/HnNJSYmqq6vlcDgUGBiowMBAFRcXa+3atQoMDHT13R/HfqmIiAiNGTNGX3/9tV//e8fHx2v06NGmdaNGjVJFRYUk//79ZrrErE+FjODgYDkcDhUVFZnWFxUVaezYsT7qlXclJSXJbrebxtzQ0KDi4mLXmB0Oh4KCgkxtqqqqdPDgwV77uRiGoSeffFLvvvuu/v73vyspKcm03V/H3R7DMFRfX+/XY540aZIOHDigsrIy15KSkqIHHnhAZWVlGjFihN+O/VL19fU6fPiw4uPj/frfe9y4cW0uSf/qq680fPhwSf3r97vfs/Y8Uu9ruYQ1Ly/P+PLLL42MjAwjIiLC+Pbbb33dtS6rq6szSktLjdLSUkOSsWbNGqO0tNR1We4LL7xgREdHG++++65x4MAB47e//W27l3oNHTrU2L17t/H5558b//Iv/9JrL/UyDMN44oknjOjoaOPDDz80Xd73008/udr447iXLFlifPTRR0Z5ebnxxRdfGEuXLjUGDBhgFBYWGobhn2PuyKVXlxiGf4590aJFxocffmgcPXrU+OSTT4zp06cbkZGRru8rfxyzYVy8TDkwMNB4/vnnja+//tp46623jPDwcGPr1q2uNv429parS8LDw42IiIguL+Hh4X51dUmfCxmGYRjr1683hg8fbgQHBxu/+MUvXJc99lUffPCBIanN8sgjjxiGcfFyr+XLlxt2u90ICQkxfvnLXxoHDhww7ePcuXPGk08+aQwaNMgICwszpk+fblRUVPhgNO5pb7ySjNdee83Vxh/H/eijj7r+3x08eLAxadIkV8AwDP8cc0dahwx/HHvLvR+CgoKMhIQE4+677zYOHTrk2u6PY27x3nvvGcnJyUZISIhx7bXXGhs3bjRt97ext4SMsLAwIzw8vMtLWFiYX4UMHvUOAEA3tTzqPSwsrFsnphqGoXPnzvGodwAAYNbdv9v97e9+QgYAABYhZJgRMgAAsAghw6xPXcIKAAD6DioZAABYhEqGGSEDAACLEDLMmC4BAABeQSUDAACLUMkwI2QAAGARQoYZ0yUAAMArqGQAAGARKhlmhAwAACxCyDBjugQAAHgFlQwAACzkb9WI7qCSAQBANwUHB8tut1uyL7vdruDgYEv25Ws2g8gFAEC3nT9/Xg0NDd3eT3BwsEJDQy3oke8RMgAAgFcwXQIAALyCkAEAALyCkAEAALyCkAEAALyCkAEAALyCkAEAALyCkAEAALzi/wFPUl8eH9o9owAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's get an x-ray image\n",
    "fig = plt.figure()\n",
    "plt.imshow(gvxr.computeXRayImage(), cmap=\"gray\")\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6ca730a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getXrayImage(x, take_screenshot=False):\n",
    "    \n",
    "    backup = gvxr.getLocalTransformationMatrix(\"root\")\n",
    "\n",
    "    # Move source, det, object using x\n",
    "    x_src = x[0]\n",
    "    y_src = x[1]\n",
    "    z_src = x[2]\n",
    "    gvxr.setSourcePosition(x_src, y_src, z_src, \"mm\")\n",
    "    \n",
    "    x_det = x[3]\n",
    "    y_det = x[4]\n",
    "    z_det = x[5]\n",
    "    gvxr.setDetectorPosition(x_det, y_det, z_det, \"mm\")\n",
    "\n",
    "    x_obj = x[6]\n",
    "    y_obj = x[7]\n",
    "    z_obj = x[8]\n",
    "\n",
    "    alpha_x = x[9]\n",
    "    alpha_y = x[10]\n",
    "    alpha_z = x[11]\n",
    "\n",
    "    gvxr.resetSceneTransformation();\n",
    "\n",
    "    \n",
    "    \n",
    "#     gvxr.rotateNode(\"root\", rot_angle, 1, 0, 0)\n",
    "    \n",
    "#     gvxr.translateNode(\"root\", x_rot_axis_pos, y_rot_axis_pos, z_rot_axis_pos, \"mm\")\n",
    "\n",
    "    \n",
    "    \n",
    "    gvxr.translateNode(\"root\", x_obj, y_obj, z_obj, \"mm\")\n",
    "    \n",
    "    gvxr.rotateNode(\"root\", alpha_x, 1, 0, 0)\n",
    "    gvxr.rotateNode(\"root\", alpha_y, 0, 1, 0)\n",
    "    gvxr.rotateNode(\"root\", alpha_z, 0, 0, 1)\n",
    "    \n",
    "    gvxr.translateNode(\"root\", -x_obj, -y_obj, -z_obj, \"mm\")\n",
    "    \n",
    "    test_image = np.array(gvxr.computeXRayImage(), dtype=np.single)\n",
    "\n",
    "    if take_screenshot:\n",
    "        global screenshot\n",
    "        gvxr.displayScene()        \n",
    "        screenshot = gvxr.takeScreenshot()\n",
    "    \n",
    "    gvxr.setLocalTransformationMatrix(\"root\", backup)\n",
    "    \n",
    "    return test_image / gvxr.getTotalEnergyWithDetectorResponse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bc3c5c5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cma\n",
    "from tifffile import imread, imwrite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "a81ce9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_nb_projs = 200\n",
    "angular_span = 360\n",
    "angular_step = 360 / 103\n",
    "\n",
    "fnames_images = []\n",
    "angles = []\n",
    "\n",
    "for i in range(total_nb_projs):\n",
    "    angles.append(i*angular_step)\n",
    "    fnames_images.append(current_folder + \"/4thOctober_image/raw_images/raw_image_\" + str(i) + \".tiff\")    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "34a45d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(\"angles (first 5):\", angles[:5])\n",
    "# print(\"fnames_images (first 5):\", fnames_images[:5])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "9c5149d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setNumberOfProjections(n):\n",
    "    global nb_projs_in_fitness, ref_images, crop_ref_images, fitness_angles, fitness_fnames, test_images\n",
    "\n",
    "    nb_projs_in_fitness = n\n",
    "    \n",
    "    ref_images = []\n",
    "    fitness_angles = []\n",
    "    fitness_fnames = []\n",
    "\n",
    "    for i in range(nb_projs_in_fitness):\n",
    "        index = int(i / nb_projs_in_fitness * total_nb_projs / 2)\n",
    "        print(i, angles[index], fnames_images[index])\n",
    "        ref_images.append(imread(fnames_images[index]))\n",
    "        fitness_angles.append(angles[index])\n",
    "        fitness_fnames.append(fnames_images[index])\n",
    "\n",
    "    ref_images = np.array(ref_images, dtype=np.single)\n",
    "    ref_images -= ref_images.mean()\n",
    "    ref_images /= ref_images.std()\n",
    "\n",
    "    fitness_angles = np.array(fitness_angles, dtype=np.single)\n",
    "\n",
    "    test_images = np.zeros(ref_images.shape, dtype=np.single)\n",
    "    \n",
    "    crop_ref_images = ref_images[:, 211:470, 100:600]\n",
    "    crop_ref_images -= crop_ref_images.mean()\n",
    "    crop_ref_images /= crop_ref_images.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "2bb968d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def displayResult(x):\n",
    "    transform(x)\n",
    "    \n",
    "    ZNCC = 100 * (crop_ref_images * crop_test_images).mean()\n",
    "    print(\"Overall ZNCC:\", \"{:.4f}\".format(ZNCC) + \"%\")\n",
    "    \n",
    "    for i, (angle, fname) in enumerate(zip(fitness_angles, fitness_fnames)):\n",
    "        \n",
    "        ref = np.copy(crop_ref_images[i])\n",
    "        test = np.copy(crop_test_images[i])\n",
    "\n",
    "        ref -= ref.mean()\n",
    "        ref /= ref.std()\n",
    "        \n",
    "        test -= test.mean()\n",
    "        test /= test.std()\n",
    "        \n",
    "        ZNCC = 100 * (ref * test).mean()\n",
    "\n",
    "        offset=1.5\n",
    "        # imwrite(current_folder + \"/test_\" + str(i) + \"tif\", test_images[i])\n",
    "        fig, axs = plt.subplots(1, 3, figsize=(10, 4))\n",
    "        plt.suptitle(fname + \", ZNCC=\" + \"{:.4f}\".format(ZNCC) + \"%\")\n",
    "        axs[0].imshow(ref_images[i],cmap=\"gray\", vmin=-offset, vmax=offset)\n",
    "        axs[1].imshow(test_images[i],cmap=\"gray\", vmin=-offset, vmax=offset)\n",
    "        # axs[2].imshow(np.abs(ref_images[i] - test_images[i]),cmap=\"gray\")\n",
    "        im = axs[2].imshow((ref_images[i] - test_images[i]),cmap=\"gray\", vmin=-offset, vmax=offset)\n",
    "        cbar = fig.colorbar(im)\n",
    "        \n",
    "        for ax in axs:\n",
    "            ax.set_xlim([100, 600])\n",
    "            ax.set_ylim([211, 470])\n",
    "\n",
    "        fig = plt.figure(figsize=(10, 4))\n",
    "        plt.plot(crop_ref_images[i][crop_ref_images.shape[1] // 2, :], label=\"Experimental\")\n",
    "        plt.plot(crop_test_images[i][crop_ref_images.shape[1] // 2, :], label=\"Simulated\")\n",
    "        plt.legend()\n",
    "        plt.show()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "19f3d233",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.0 C:\\Users\\snn23kfl\\project/4thOctober_image/raw_images/raw_image_0.tiff\n",
      "1 174.75728155339806 C:\\Users\\snn23kfl\\project/4thOctober_image/raw_images/raw_image_50.tiff\n"
     ]
    }
   ],
   "source": [
    "setNumberOfProjections(2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "1b2da784",
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform(x):\n",
    "    \n",
    "    global test_images, fitness_angles, crop_test_images\n",
    "    backup = gvxr.getLocalTransformationMatrix(\"root\")\n",
    "    \n",
    "    for i, angle in enumerate(fitness_angles):\n",
    "    \n",
    "        # Rotate\n",
    "        gvxr.resetSceneTransformation();\n",
    "        gvxr.translateNode(\"root\", x[0], x[1], x[2], \"mm\")\n",
    "        gvxr.rotateNode(\"root\", angle, 0, 0, 1)        \n",
    "        gvxr.translateNode(\"root\", x[3], x[4], x[5], \"mm\")\n",
    "\n",
    "\n",
    "        # Simulate\n",
    "        test_images[i] = np.array(gvxr.computeXRayImage(), dtype=np.single)  \n",
    "        \n",
    "        gvxr.setLocalTransformationMatrix(\"root\", backup)\n",
    "\n",
    "    test_images -= test_images.mean()\n",
    "    test_images /= test_images.std()\n",
    "\n",
    "    crop_test_images = test_images[:, 211:470, 100:600]\n",
    "    crop_test_images -= crop_test_images.mean()\n",
    "    crop_test_images /= crop_test_images.std()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "37214161",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fitness(x):\n",
    "    \n",
    "    global ref_images, test_image, fitness_angles, nb_projs_in_fitness, crop_ref_images\n",
    "    \n",
    "    transform(x)\n",
    "        \n",
    "    return np.abs(crop_ref_images - crop_test_images).mean()\n",
    "    # return np.abs(ref_images - test_images).sum() / nb_projs_in_fitness\n",
    "    # return 1.0 - (ref_images * test_images).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c1a82779",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "could not broadcast input array from shape (480,640) into shape (480,640,3)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[27], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m displayResult(\u001b[38;5;241m6\u001b[39m\u001b[38;5;241m*\u001b[39m[\u001b[38;5;241m0.0\u001b[39m])\n",
      "Cell \u001b[1;32mIn[23], line 2\u001b[0m, in \u001b[0;36mdisplayResult\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mdisplayResult\u001b[39m(x):\n\u001b[1;32m----> 2\u001b[0m     transform(x)\n\u001b[0;32m      4\u001b[0m     ZNCC \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m100\u001b[39m \u001b[38;5;241m*\u001b[39m (crop_ref_images \u001b[38;5;241m*\u001b[39m crop_test_images)\u001b[38;5;241m.\u001b[39mmean()\n\u001b[0;32m      5\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mOverall ZNCC:\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{:.4f}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(ZNCC) \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m%\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "Cell \u001b[1;32mIn[25], line 16\u001b[0m, in \u001b[0;36mtransform\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m     12\u001b[0m     gvxr\u001b[38;5;241m.\u001b[39mtranslateNode(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mroot\u001b[39m\u001b[38;5;124m\"\u001b[39m, x[\u001b[38;5;241m3\u001b[39m], x[\u001b[38;5;241m4\u001b[39m], x[\u001b[38;5;241m5\u001b[39m], \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmm\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m     15\u001b[0m     \u001b[38;5;66;03m# Simulate\u001b[39;00m\n\u001b[1;32m---> 16\u001b[0m     test_images[i] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(gvxr\u001b[38;5;241m.\u001b[39mcomputeXRayImage(), dtype\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39msingle)  \n\u001b[0;32m     18\u001b[0m     gvxr\u001b[38;5;241m.\u001b[39msetLocalTransformationMatrix(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mroot\u001b[39m\u001b[38;5;124m\"\u001b[39m, backup)\n\u001b[0;32m     20\u001b[0m test_images \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m=\u001b[39m test_images\u001b[38;5;241m.\u001b[39mmean()\n",
      "\u001b[1;31mValueError\u001b[0m: could not broadcast input array from shape (480,640) into shape (480,640,3)"
     ]
    }
   ],
   "source": [
    "displayResult(6*[0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ba4cf2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}