StructureTensor2D_Examples-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D Structure tensor -  a few examples\n",
    "A few examples of using 2D structure tensor. The examples show different ways of visualizing the orientations, visualizing the distribution of orientations, and how parameter choices influences the orientationa analysis.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import skimage.io\n",
    "import scipy.ndimage\n",
    "import matplotlib.pyplot as plt\n",
    "import st2d    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1\n",
    "Orientations computed from simple testing image, and visualized using orientation arrows or color. \n",
    "\n",
    "When visualizing predominant orientation information as color without the overlay of intensity image, it may seem as orientation changes abruptly in some areas. However while dominant orientation changes, the anisotropy is low in these areas. This means that transition from one orientation (ellipse elonagted in one direction) to another orientation (ellipse elongated in another directrion) is over an isotropy (a circle). Therefore, the full orientation information (dominant orientation and anisotropy) is smoothly changing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The image has a shape (50, 50), i.e. 2500 pixels.\n",
      "Structure tensor information is carried in a (3, 2500) array.\n",
      "Orientation information is carried in a (2, 2500) array.\n",
      "Anisotropy information is carried in a (2, 2500) array.\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'plot_orientations' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-a3b5350fc3be>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msharex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msharey\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgray\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0mplot_orientations\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvec\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     18\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Orientation as arrows'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0morientation_st_rgba\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhsv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marctan2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvec\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvec\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'plot_orientations' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABH4AAAEvCAYAAAAzXwbsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dW2xd133v+98QxZskUhRFXSiSEnWXLUu2bPkSO0biNGnTpEj64AM0DTZcIIVf9gZS7AJteg5wcAqcAu1Lb0CxD4ymqB+CpN1NigRpnNRxrDiJY8WSLauSdb9QokSJIkWKoiTex3nQ8sz4/0kuUhcuLk9/P4Dh+eNYl7HmWrNBh+f8zRBjFAAAAAAAAPJnwXxPAAAAAAAAAHODhR8AAAAAAICcYuEHAAAAAAAgp1j4AQAAAAAAyCkWfgAAAAAAAHKKhR8AAAAAAICcWjjfEwAw/5qammJ7e/t8TwOYF/v37++JMa6Y73lMhWMTH2Ucm0B5Ktdjk+MSH3XFjk0WfgCovb1d+/btm+9pAPMihNAx33OYDscmPso4NoHyVK7HJsclPuqKHZtc6gUAAAAAAJBTLPwAH1EhhBdDCPtCCPuuXLky39MBUMCxCZQnjk2g/HBcArPDwg/wERVjfCnGuDvGuHvFirK7TBv4yOLYBMoTxyZQfjgugdlh4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCn7mnhJ4Tw2RDCsRDCyRDC1+7XpAAAAAAAAHDv7nrhJ4RQIekfJP22pAclfSmE8OD9mhgAAAAAAADuzcJ7eO4Tkk7GGE9LUgjhW5K+KOn96Z4QQogu28kstNMZGxvLtmM0T1VFRYXJfnymvGDB9GteExMTJldWVprs5+0fX11dXfTxixYtyrZra2uLvtbo6KjJ6T6RpJs3b5o8Pj5ucvq5R0ZGis7TP9e/l/8cfh8W20/Dw8NmzH9//rl+Lv7xXvo5Z3qsf22/z4vtQ8l+rmK/Iz+XkZERjY2NhSIPBwAAAADgvrqXhZ8WSeeT3CnpyTt5Af//7K9atcrk7u7ubNsvHDQ0NJjsFzX84/14TU2Nyen/g37jxg0z1tTUZLJfMBkcHDR548aNJvvFnV27dmXb27dvL/paly5dMvnq1asm79+/3+Rr166ZnC7enD171oytW7eu6Gv39PSYXFVVZfKSJUtMXrFihcnpPj558qQZW7ZsmcnNzc0m9/f3m+y/by/9vv1j/UKO38d+8cy/t198S79/v0/84lg6l+PHj085dwAAAAAA5sq9LPxMdeZCnPSgEF6U9OI9vA8AAAAAAADuwr0s/HRKaktyq6SL/kExxpckvSRNvtTLn4XjpWfK+DNCOjs7TfaXJXn+7CJ/eY8/66MYf+aMv7TInzGycuVKk9OzQPxj/Vk0/vK3jo4Ok69fv26y36dDQ0PZtj8T6cSJEyb7M648f1mT32e9vb0m19XVZduLFy82Y36fnTt3zuT6+vqi75V+LsnutzNnzhSd99KlS4uOp5fiTTX39KyeW7dumbHVq1eb3NXVlW37M48AAAAAAJhr93JXr7clbQ4hrA8hVEn6PUnfuz/TAgAAAAAAwL266zN+YoxjIYT/IelHkiok/VOM8fB9mxkAAAAAAADuyb1c6qUY4w8k/eA+zQUAAAAAAAD30T0t/NypyspKc4estP9Emnw3rbT/xnfh+LtC9fX1FX1v3+PiO2JSvpPH36mpsbGx6Gv7PiE/1/Tx/r38Xbn8PvHZfw7fdTQwMJBtz3Sbev/enn++7wzydwVLe3v85/SdPf6W6f779p1Aa9asmXae69evN9l36/jfne/88Xdx89JeH38ns7Vr1077Wr4DCQAAAACAuXYvHT8AAAAAAAAoYyz8AAAAAAAA5BQLPwAAAAAAADlV0o6f0dFRdXd3Z3nTpk1mPO1O8RYtWmTyTJ0+3vDwsMm1tbUmL1z4612RzlGa3NHju24GBwdN3rZtm8m+QybtCKqpqTFjvsMn7ejx85xqbkeOHDF5dHQ02/Z9QKtWrTLZdxX5LiPfw+Pn6qW9Pr57yPP7KJ23JLW1tZm8efNmk9PvzD/X9wPt2LHDZP8dXLhwwWS/35YvX55tf/zjHzdjHR0dJqffzy9/+UsBAAAAAFBKnPEDAAAAAACQUyz8AAAAAAAA5FRJL/Wqqqoyt+H2t+z20ltlX7lypehj/aVCM92a3Lt+/fq0Y/61/O3dZ7pkyl9WVl9fn237y6f8ZUUXL14smv2lRXV1dSant1hvaGgwYz09PSYvWbLEZH95lZ+bv+TNf+708jx/OZV/rL+1vP++/e3f/eV4VVVV2ba/xfrDDz9c9L137dpl8rFjx0w+evSoyekt23/zN3/TjPnLyl555ZVs2/8OAAAAAACYa5zxAwAAAAAAkFMs/AAAAAAAAOQUCz/AR1QI4cUQwr4Qwr6ZLqUEUDocm0B54tgEyg/HJTA7Je34GRkZ0dmzZ7PsO2X8bdHT24n73hZ/YC9evNjkmW5dfu7cuWnnOdN7tbe3m+xvQ+97dvzj0+6b3t5eM1ZdXW1y2l0j2c4eaXLPjn9+S0tLtu17dvzt2v28/S3Y/T4eHh4u+vj0luvpZ5ak8fFxk9PeI0l67733TE5voS5JO3fuNDndL9u2bTNjvg/I9/L4fejnkt6WXrL7cc+ePWbM9wel+99/l/MtxviSpJckaffu3XGGhwMoEY5NoDxxbALlh+MSmB3O+AEAAAAAAMgpFn4AAAAAAAByioUfAAAAAACAnCppx4+3ceNGky9cuGBy2iEzMjJixtatW2fy+fPnTfa9Lb6np6KiwuSFC3+9K0ZHR83Ypk2bTO7p6THZP/4//uM/TH7qqadMTrt1Qghm7OjRo0XnXVtba/LNmzdN9j096eN9F9HSpUtN9v1Afp/7jqC2traic03ndu3aNTP2xBNPmOz7gT796U+bnH4/kvTuu+9O+/hjx46ZMf9b+fu//3uTGxoaiuZly5aZ3NfXl23/4he/MGP++0k7f/r7+wUAAAAAQClxxg8AAAAAAEBOsfADAAAAAACQUyz8AAAAAAAA5FRJO35qa2u1devWLJ86dcqML1q0yOTm5uZs2/fHdHR0mPzQQw+ZfOjQIZN9v42Xdu2sXr3ajPkOn8HBQZN9r8769etNfuedd0zetm1btu27b1atWmXys88+a/Jbb71lsu8qampqMjntlfF9Qn6ftLa2mlxVVWWy7/y5ceOGyb4LKd0v77//vhn7+c9/bvIjjzxi8tDQkMn+t+Ln2tnZmW3738b169dN9vvs7NmzJscYTfYdP2nX0Uz7ZHx8PNv2vyMAAAAAAOYaZ/wAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6VtOOnurra9N9UV1eb8cbGRpPPnDmTbT/99NNm7MKFCyY3NDSY7LtyfA9PS0vLtK/ne3IWLrS7qb293WTfP7Ry5UqTfR9N2hG0ceNGM7Zu3TqT+/r6TF6woPhaXW1trcnpPr1586YZ6+3tLZrTfiBp8j7bsWOHyZcuXZr29Xy/kN+n+/btM9n36ixevNhk30+UPt939Fy+fNlk/zvzHUC+t+f06dMmr127NtseGxszY/630NbWNu1jAQAAAACYa5zxAwAAAAAAkFMs/AAAAAAAAOQUCz8AAAAAAAA5VdKOn8rKSq1evTrLaf+JNLlTZvv27dn2wMCAGauqqjK5pqbG5PR9PnjvVLHens7OzqLvNTo6avLmzZuLjvtul/r6+mx7y5YtRed57do1k32vztKlS032XTh79uzJtn03Udo1NNVz/T6amJgw+fjx4yann0uy3Tnj4+NmzHfhbNu2zWTfy+N7ePzzi/H71Hf4eP61/W/r7Nmz0z7X77P0N+33AQAAAAAAc40zfgAAAAAAAHJqxoWfEMI/hRC6QwiHkr81hhBeDSGcKPx7WbHXAAAAAAAAQOnN5oyff5b0Wfe3r0l6Lca4WdJrhQwAAAAAAIAyMmPHT4zxjRBCu/vzFyV9srD9sqQ9kv50pteqqqrSunXrsrxx40YzfvHiRZP7+vqy7Vu3bpmxNWvWmOz7Zfzj0/eVJve8dHd3Z9u1tbVmrKWlxWTf1XLp0iWTN23aZLLvCGpqatJ0fL/M8uXLTb569arJMUaTb968aXLanePnuWLFCpN7e3tN3rBhg8l+n/nOIN8BlPbdnD592oz5z3Xs2DGTh4eHTfbfp3+vdB/7z+F7kPw+q6urM3nZsuInsPX09GTbCxbYtVP/2unvync9AQAAAAAw1+6242dVjLFLkgr/Xnn/pgQAAAAAAID7Yc7LnUMIL4YQ9oUQ9t24cWOu3w4AAAAAAAAFd3s798shhOYYY1cIoVlS93QPjDG+JOklSdq2bVt85plnsjF/2VJDQ4PJ6WU0/tIefwt1r7Gx0eQjR46Y3NXVZXJ66ZhfoPKXKflbxftbkftLu/ylQ+nlWP7W4P6W6n4uM93Gvq2tzeT0duLV1dVmzF8S5V/bX5rkL7datGiRyYcOHTI5vfzOX5rnLzMbGBgomisqKkz2l+Oll2v5y/7SSwalybex95edHT161GR/6V56OZf/nS1ZssTk9LKwEIIAAAAAACiluz3j53uSXihsvyDpu/dnOgAAAAAAALhfZnM7929K+qWkrSGEzhDCVyT9paTPhBBOSPpMIQMAAAAAAKCMzOauXl+aZug37vNcAAAAAAAAcB/dbcfPXRkbGzO9Pjt27DDjv/jFL0xOe122bt1qxnxHzOjoqMm+K2fLli0m+16XtPPn6aefNmO+o8c/19+G3t8G3XcEpZ0zvi/Idxldu3bNZN+z43t6it3W3t9+vaOjw+SdO3ea7G/fnvYFSZNvHe9vi16s3+bcuXMmz3SLdf98332U7lM/L29oaMjk48ePm+z7gy5cuGByuk99p49/77QnyX9GAAAAAADm2pzf1QsAAAAAAADzg4UfAAAAAACAnGLhBwAAAAAAIKdK2vGzYMECLVq0yOTUxz/+cZPTvhvfs5J2BUlSS0uLyZWVlSb7Xh7flfPcc89l2763xfcLnT171uTHHnvM5AMHDphc7HP5eZ04ccLkpUuXmuzntnnz5qLvnX7OkydPmjG/z86cOWOy7/Txz/d8J9Dly5ez7eHhYTPme5C8hoYGk/fv32/ywoX2p9vU1JRt+14k3wc1k6qqKpN9J1DK9yT5HiXf/wQAAAAAQClxxg8AAAAAAEBOsfADAAAAAACQUyz8AAAAAAAA5FRJO36uXr2qb3zjG1nesWOHGfe9LrW1tdm274jxz+3t7TV5zZo1JvuulYqKCpPHxsay7cbGxknzTp07d67oez/xxBMmp103ktTV1ZVtv/nmm2bMd/j4zp50n0hSdXW1yX4f7t27N9uuq6szY93d3Sb7XiTfZTQwMGCy7/SJMZqc7mM/T79Pjx07pmLWrl1r8vXr101O5+57k3yX1COPPGLywYMHTfa9SelvQ5JGR0ez7cOHD5uxmzdvmpzug/HxcQEAAAAAUEqc8QMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOVXSjp/x8XH19/dn+T//8z/N+KpVq0xevHhxtu27a3wvy7Jly0z2fTS+98WPp302ixYtMmO+j8Z3APlel+985zsmNzc3m3zjxg1Nx/fk+M/t+4L8XE+cOGFy2id09OhRMzYxMWGy70Hy436f+e6b+vr6acdramrMWF9fn8l+n/quI7/P/D5NfzttbW1mrKWlxWTfdfTlL3/5jt473cd+Hr7/KX2tzs5OAQAAAABQSpzxAwAAAAAAkFMs/AAAAAAAAOQUCz8AAAAAAAA5VdKOn9HRUXV3d2c57fuRpMrKSpPTPhvfs1JbW1v0vXxXzvj4uMkVFRUm37p1a9r38v00vkcn7SKSpN7eXpN9V87q1auzbd8P1N7ebvL169dN9p0+PT09Jvs+otdffz3b9p08Pvt5NjU1Fc1DQ0Mmj4yMaDqjo6Mmz/T9+Nd+4IEHTF6xYoXJmzZtyrYfe+wxM+b7g3zf08qVK00+deqUyf77XLp0abbt98mTTz5p8rp167Ltv/iLvxAAAAAAAKXEGT8AAAAAAAA5xcIPAAAAAABATrHwAwAAAAAAkFMl7fipqKjQkiVLslxfX2/GfedP+ljf6eM7enxfTfpcyXb4SJP7aEII2bbvo/EWLrS7zffTVFVVmZx2FUnSmTNnsu1ly5aZMd9t09raarLvE/KdP/690/6azs5OM+Y7fRoaGkz234fv5fGPv3btmsnpPq6rqzNjfh96vjsn7dWRJnchPfPMM9m27+zxHU3+d+f5fd7X12dy+ltpa2szY/77efzxx7Ptv/u7vyv6vgAAAAAA3G+c8QMAAAAAAJBTLPwAAAAAAADkVEkv9WpsbNTv//7vZ/n06dNmfGBgwOTNmzdn2/6yJH+L7fTyKWnyJU/+UjB/6Vd6W/SWlpai81qwwK6X+cvG/C3a/Xull6n5W8d3dXWZ7C9Z85e8+du9p7eKl6Rdu3Zl2/7W7/6SNv85ZroUz38Hfm7pJVZ+n/h5+surtmzZYnL6W5AmXyKXXgrmL7fyl6gdP37cZH/ZWXNzs8lr1qwxOf3+/KV21dXVJvvLAAEAAAAAKCXO+AEAAAAAAMgpFn4AAAAAAAByioUf4CMqhPBiCGFfCGHflStX5ns6AAo4NoHyxLEJlB+OS2B2StrxU1dXp09+8pNZ9j0uviMm7U+prKw0Y4ODg0Wz76/xHT/+Ft9vv/12tu07eQ4cOGDy8PBw0ffyfUQdHR0mp5/Lz8vfzt3fCn7FihUm+w4gn9POmQcffNCMHTp0yOTGxkaTZ+rl8XNbvny5yel+WrVqVdF5fuITnzDZ9yz5XiV/u/f01vLp7dalyd/14cOHTb548aLJ/nMXuxV9eht5aXI/0KlTp7Jt/7uZbzHGlyS9JEm7d++OMzwcQIlwbALliWMTKD8cl8DscMYPAAAAAABATrHwAwAAAAAAkFMzLvyEENpCCK+HEI6EEA6HEL5a+HtjCOHVEMKJwr+XzfRaAAAAAAAAKJ3ZdPyMSfrjGOM7IYQ6SftDCK9K+gNJr8UY/zKE8DVJX5P0p8VeaGRkRBcuXMjy+Pi4GV+6dKnJZ86cybZ37Nhhxnzvju++8T09AwMDJj/33HMmd3Z2Zttr1qwxY9u2bTM57XiRpOPHj0/7WpK0bJldE/vVr36Vbfvuohs3bph87do1k2O0l676rhy/T9O+Gt/x4/f3z372M5Pr6+tNTr87afI+r6ioMDntDPL9P/61ff+N7/Txn3PXrl3Tzu1HP/qRGfO/hTfffNNk38nU399vcltbm8np9+9fy++DdB/39PQIAAAAAIBSmvGMnxhjV4zxncL2dUlHJLVI+qKklwsPe1nS787VJAEAAAAAAHDn7qjjJ4TQLmmXpL2SVsUYu6Tbi0OSVt7vyQEAAAAAAODuzXrhJ4SwRNK3Jf1RjHFgpscnz3sxhLAvhLDPX0IDAAAAAACAuTObjh+FECp1e9HnGzHG7xT+fDmE0Bxj7AohNEvqnuq5McaXJL0kSRs3bowXL17MxqqqqsxjfZ/Nu+++m22///77Zsz37qRdNpJ05coVk32Py9tvv21yTU3NtK+1atWqou/1+OOPm7xz506T/+AP/sDk73znO9n2wYMHzZj/nKOjoyYfPnzY5Lq6OpN7e3tNXrx4cbbte3ZaW1tN/tjHPmbyG2+8YfLatWtN9l1HtbW1JqcdQL4PyO8z3+EzNjZmcgjB5CNHjpg8ODiYbaefWZr8OS5dumTy0NCQyb5v6OrVqyY3NTVl26dOnTJjvpso7Tny/U0AAAAAAMy12dzVK0j6uqQjMca/Toa+J+mFwvYLkr57/6cHAAAAAACAuzWbM36ekfTfJP1XCOGDW2n9n5L+UtK/hhC+IumcpP9jbqYIAAAAAACAuzHjwk+M8eeSwjTDv3F/pwMAAAAAAID7ZVYdP/fLjRs3tHfv3mnHBwZsZ3TarXPo0CEzdvLkSZMXLVpksu9aSXtZJNvpI0n19fXZtu8Dunz5ssl+3L/3ypX2Bme+t2fr1q3Z9kMPPWTGOjs7Tfb9NO+9957JvuvGd+mk3Th9fX1mLMZYdN6+88fPbfXq1Sb7np6NGzdm21u2bDFjDzzwgMkjIyMmp509ktTdbSuk3nrrLZNPnDiRbft95HuP/D7y8/bfp+90SjuBfJ+Qf6329vZs2+9/AAAAAADm2h3dzh0AAAAAAAAfHiz8AAAAAAAA5BQLPwAAAAAAADlV0o6foaEhHT9+PMuVlZWTxlNdXV3Zdm1trRnzHT3p60qT+4La2tpMfvLJJ6d9/Pj4uBlLu2okqaqqymT/+DNnzpjs+20uXryYbftum+rqapO3b99usn+8743x75120Ph95ve/77ZJ+2mkyT1JaX+QJG3YsMHkHTt2ZNvNzc1mzO+z4eFhk/1cb9y4YfLrr79uctp1tHTpUjPme3guXLhgcktLS9G5XblyxeS022jhQnsI+ee2trZm2wcPHhQAAAAAAKXEGT8AAAAAAAA5xcIPAAAAAABATrHwAwAAAAAAkFMl7fjxjh49arLvS0l7XSoqKszY6OioyQsW2DWsW7duFc3++f39/dm27/Tp7u42Oe1tkSZ34/geHt8xs2nTpmy7s7PTjI2NjZl87do1k9evX2+y7+nxXThpd47vA/Lzvnr1qsnr1q0z+Qtf+ILJfp9OTEyYnPYq+Q6fdH9LUmNjo8mXL1822e+n+vp6k9Pfjv++/G9n5cqVRcf951ixYoXJN2/ezLb99+V/G+n35fubAAAAAACYa5zxAwAAAAAAkFMs/AAAAAAAAORUSS/1mpiYMJcirV271oz7S27S2577y3f8JTb+MjF/CZS/hfeSJUumfS9/a3h/+/bnn3/eZP85/CVQ/jKldK7+Eqf0VuGSdPjwYZO3bt1qsr9VvL+dezoXPzY4OGjyzp07i+b00i1JGhgYKJr95V0pf8v1s2fPmnzq1CmTz507Z/KJEydMHhoayrZra2vNmL+k7fr16yb7S9q2bdtm8oEDB0xOLwXzt7T3v7OOjo5s218CCAAAAADAXOOMHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIqZJ2/CxcuFDLly/PcnNzsxk/ePCgyTU1Ndm272nxvSzvvPOOyf622r29vSa/9957Jvvbpqd8143vAPLdLQ0NDSYfOXJk2tf2XTe+X+bzn//8tM+VpPHxcZMfeughk9Pbpnd1dZkx34Xjb1vue5XSHiRJGh0dLTqXdNzvf387d99ldPLkyaLv7W/nnnY+pb8baXJHU0tLi8nbt283efHixSanv1n/en6f+s4m3wEEAAAAAEApccYPAAAAAABATrHwAwAAAAAAkFMs/AAAAAAAAORUSTt+Fi1apEcfffTXb77Qvr3vXkl7d3zfzP79+4s+d3h42GTfZ3Pz5k2TN23alG37fpmzZ88Wfa6fm3+vuro6k5988sls23fC+O4a3xHj89WrV02+fPnytHNtamoyY743yXf2+H14/vz5ou/d19dncvr9+p4jv097enqKvvfg4KDJaaePZPuJ/Of0PUpDQ0NFX3um/qf0+12wwK6dxhhNTveJnzMAAAAAAHONM34AAAAAAAByioUfAAAAAACAnGLhBwAAAAAAIKdK2vGzZMkSfexjH8vy6dOnzXhbW5vJv/zlL7Nt3wf02GOPmdzR0WHywMCAyWfOnDG5oaHB5LTXx3fC+M4e3+nj38uPe2kPjO+28a/1zDPPmFxdXV308b63p6urK9v2/UDLli0zeWJiwmTf6eM7ajo7O032fUX9/f3TPtbP238fftx/J+Pj4yaPjIxM+b7S5B4e3y/kO3z8fvLdR5WVldPO0+/D9vb2bPutt94SAAAAAAClxBk/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTJe34iTGafpu1a9ea8aGhIZMffPDBbHvlypVmzPe4+H6Zffv2mbxjxw6T+/r6TL58+XK27TtffPaGh4eLjns//elPs23/OdKuIWlyr86tW7dMbm5uNnnFihUmb9y4Mdv2PUeHDx822XcTHT9+3OTFixebXFNTY3JFRYXJ169fz7YfffRRM/bjH//Y5Js3b5rs96nvLtq0aZPJly5dmva1fB+Q/61s377dZL+f/OdK+4QeeughM+b7oNLvy/f/AAAAAAAw1zjjBwAAAAAAIKdmXPgJIdSEEH4VQngvhHA4hPDnhb+vDyHsDSGcCCH8Swihau6nCwAAAAAAgNmazRk/w5I+FWN8WNIjkj4bQnhK0l9J+psY42ZJfZK+MnfTBAAAAAAAwJ2aseMn3i7lGSzEysI/UdKnJP1+4e8vS/p/JP2vYq9VV1enT3ziE1k+deqUGb969eq02ffLtLa2muy7cdJuG2lyZ8z+/ftNbm9vz7Z9r47vAwohmLxmzRqTfcdMZ2fntK+3Z88eM+Z7j3zfzLJly0yuq6szeenSpSanvUm+62ZgYMDkEydOmFxfX2/y4OCgybt27TLZv37ay5P24kjSww8/bHJVlT1hrKenx+S0L2iquaxbt27K95VkeqWkyd1EV65cMdn3QRV7r4sXL5qxyspKk9Pf4SuvvCIAAAAAAEppVh0/IYSKEMIBSd2SXpV0SlJ/jPGDFZJOSS1zM0UAAAAAAADcjVkt/MQYx2OMj0hqlfSEpAemerS33r4AABvQSURBVNhUzw0hvBhC2BdC2OfPrAAAAAAAAMDcuaO7esUY+yXtkfSUpIYQwgeXirVKujjNc16KMe6OMe72txoHAAAAAADA3Jmx4yeEsELSaIyxP4RQK+nTul3s/Lqk5yV9S9ILkr4702uNjY2Z7hbfxeI7ZdKunBs3bpgx3xnT399vsu+f8d03aU+LJH3rW9/KticmJsyYf2/v2rVrJldXV5u8ZMmSabPvtjl37lzR9/ZnTfkuHN8/1NjYmG13dXVNOyZN7ujxXUdtbW0m+9dbsMCuI6Z9RA0NDWbMdy7578+P+94d33XU3NycbR89etSM+e/61q1bJvvOH//b2rJli8np53z22WfN2Pvvv29y2h3l3wcAAAAAgLk248KPpGZJL4cQKnT7DKF/jTF+P4TwvqRvhRD+X0nvSvr6HM4TAAAAAAAAd2g2d/U6KGnXFH8/rdt9PwAAAAAAAChDd9TxAwAAAAAAgA+P2Vzqdd/09/fru9/9dRWQ73F5+umnTU67VDZu3GjGfF+Kf63Vq1cXnYvvp0nf+9ChQ2asu7vbZN9F5Dt+fM+O75QZHBwsOrdir3Xxou3QrqioMNn38KTvtXz5cjPmO3nq6upM9j06ly5dMnn37t0m+26k9PUXL15sxjZs2GCy7zo6e/asyf5z+T6itC9q586dZsyXis90dzn/+G3btpmcdgA98MBUN7j7tbSnCgAAAACAUuOMHwAAAAAAgJxi4QcAAAAAACCnSnqp18jIiDo7O7PsL1P6wQ9+YPKjjz6abZ8+fdqM+UuD/O3bjx8/brK/NMhferRq1aps29+O/eGHHzb51VdfVTH+c/lLvSorK7Ntf6vxkydPmuxv1+5vix5jNHl0dHTax/t9ls5jqrx582aT/e3f/SVuPqeXZ/nLyNLbr0vS0NCQyf5yOn/5lb/0b+HCX/+U/aVc/lIt/9vwt2v3/GVm6e3if/jDH5qxYpe7+d/gfAshvCjpRUlau3btPM8GwAc4NoHyxLEJlB+OS2B2OOMH+IiKMb4UY9wdY9ztF9YAzB+OTaA8cWwC5YfjEpgdFn4AAAAAAAByioUfAAAAAACAnCp5x8+ZM2eyPFNfzXvvvZdtp7frlqSmpqai2ffyXLhwwWR/KmD6fN8B42/nvnXrVpP7+vpMHhgYMNl355w7dy7b9p0+aa/RVPwt2P1t6devX29yut/87dn9/va9Ob6bqLa21uRnn33W5FOnTpmcdgr5/Z128kiTb8Hu+3D844eHh01Ovz//W/G9R+3t7SannT3S5P3k5/LOO+9k2ytXrjRjvsuoo6Mj2/b9PwAAAAAAzDXO+AEAAAAAAMgpFn4AAAAAAAByioUfAAAAAACAnCppx8/o6Kjpy/G9O75jJu1HeeKJJ8yY74x56623TH7uuedMbmlpMfn99983+ROf+ES27btYbt68aXLaXTPV4/247wBKe2FqamrM2NjYmMm+B8n3zfieHd85c/78+Wzb9//4PqHe3l6Tly5danLauSRN/pwPPvjgtM9/7LHHzFjacyRN/u79fqivry861yVLlky5LUlXrlwxee/evSb73qTDhw+b7Puj0u4jv098n1M6bzp+AAAAAAClxhk/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTJe34GRkZMb09vsdl4UI7nf7+/mzbd91cvHjRZN9t4/tnrl27ZnJ7e7vJp0+fzrZ938ypU6dMXrx4scmNjY0mDwwMmOz7hbq6urJt3y+zbt06kwcHB01ua2sz+amnnjI5xmjy2rVrs+3q6mozdvnyZZN9p4/vyvnDP/xDk/0+379/v8m7d++e9r3S/S3JdD9Jkz/3pUuXTE73oX/+hg0bzJjv4QkhmNzT02Py1atXi85l2bJl2fbRo0fNWG1trcnp/ve/dwAAAAAA5hpn/AAAAAAAAOQUCz8AAAAAAAA5xcIPAAAAAABATpW042diYsL0pfg+FN8Rk3a1HDp0yIw99NBDJvseFv/4RYsWmZz2tEjS5s2bs+29e/cWfS/fRXTy5EmTfWfMgQMHTP7kJz+ZbfvuIZ9958/27dtNHh8fN9l3Fy1fvjzbnpiYKDpP3ze0c+dOk33X0blz50xevXq1yen3WVdXZ8Z8vnHjhsk1NTUmV1RUmOz7h9L94j/X+fPnTd6zZ4/J/rfgf5e+G+ngwYPZdlNTkxnzXUUnTpzItq9fvy4AAAAAAEqJM34AAAAAAAByioUfAAAAAACAnGLhBwAAAAAAIKdK2vETY9TY2FiWQwhmfNWqVSan48PDw2bs0qVLJq9cudLkN99802Tf09Pb22ty2hnje1taW1tNjjGaPDo6avLIyIjJvnenr68v296yZYsZ27p1q4r5zGc+Y/KxY8dM9p0/V65cybYbGhrM2NmzZ032XTdvvPGGyb77xncd+S6djo6ObLuqqsqM+X169epVk323ke9o8n1CaQfQ0qVLzZjvIvLv5b9P//2nn0Oyv52uri4zNjQ0ZHK6j9LfPgAAAAAApcAZPwAAAAAAADnFwg8AAAAAAEBOsfADAAAAAACQUyXt+JFsF8vJkyfN2IIFdh2qra0t2+7p6TFjN2/eNHnNmjVFx31vS319/bR5+fLlRR/ru27Wr19vsu+fOXLkiMlpl87GjRvN2Nq1a00eGBgw2fcH+ff2Dhw4kG37fdTf329y2gckSTdu3DDZ9/T43h3/HaWdQm+//XbR1/6t3/otkxcvXmxy+ruRJvdBpfvp29/+thl76623TPb9Tr7zx3f6+O8//W35niTfo5T2VNHxAwAAAAAoNc74AQAAAAAAyKlZL/yEECpCCO+GEL5fyOtDCHtDCCdCCP8SQqia6TUAAAAAAABQOndyqddXJR2R9MF1L38l6W9ijN8KIfx/kr4i6X8Ve4Hx8XFdu3Yty/7SF3/pUHo5kL98yt9K3F++U1NTY/KhQ4dM9pd+pZcK+Vus+9t9t7S0FB33lwpt2LDB5PRSsccee8yM+Uug/Dz97dqPHj1qsr9c67XXXsu2X3311WnnIU3eh/778e/tb5vu55JKb2E/lfR3MdVrPfDAAyb7z5nO/Sc/+UnRx/rP4W9D7y8zu3Tpksnpb9Hflt6/V3t7e7Z9/fp1AQAAAABQSrM64yeE0Crp85L+sZCDpE9J+rfCQ16W9LtzMUEAAAAAAADcndle6vW3kv5E0kQhL5fUH2P84JSQTkktUz0RAAAAAAAA82PGhZ8Qwu9I6o4x7k//PMVD4xR/UwjhxRDCvhDCPn9JFAAAAAAAAObObDp+npH0hRDC5yTV6HbHz99KagghLCyc9dMq6eJUT44xviTpJUmqqqqKS5YsycZ8D093d7fJExMT2fbKlSvN2KlTp0z2t0x/6qmnTPa32fa9Lunrb9261Yz52337W6j7W8en85ak6upqk9Pb1vtb2Ptbrvvbha9YscJk3ynj+2jSPqKLF+1X5PeBz/69/Fz93NLvVrL7ze9D/336Tp+6ujqTff+Tv537u+++m20PDg6aMf99+B4lz/9WLly4YHLaB+XnOTQ0ZPKtW7emnQcAAAAAAHNtxjN+Yox/FmNsjTG2S/o9ST+JMX5Z0uuSni887AVJ352zWQIAAAAAAOCOzfp27lP4U0n/M4RwUrc7f75+f6YEAAAAAACA++FObueuGOMeSXsK26clPXH/pwQAAAAAAID74Y4Wfu7VxMSEhoeHs9zb21v08cuWLcu2ff+P77LZtm1b0dd67733TN61a5fJq1evzraPHTtmxpqbm00+fvy4yb5v5rHHHjP5+vXrJqcdNL6j55VXXjH58uXLJvv+mb1795pcW1tr8vj4eLbte5GWL19u8vbt201+8803Ta6vrzfZf3/+s4yMjGTb6fcuSf39/Sb/7Gc/M9n39Pgund27d2s6/vvq6Ogw2X8fnt/nxfjuIf/a6T5KvwsAAAAAAErhXi71AgAAAAAAQBlj4QcAAAAAACCnWPgBAAAAAADIqZJ2/IyPj8/Y65NKO2Ju3bpV9LHV1dUmnz9/3uQ1a9aY3NfXN+14TU2NGfP9Qg8//LDJTU1NRedWUVFh8o9//ONsu7Ky0oz5HpiZuo0WLrRf4dGjR01OO3/S/SlN7rI5c+aMyTFGk8+dO2fy0qVLTe7s7DR55cqV2bb/XNeuXSv6Xr4rZ/HixSb39PSYnPZB+ddqbW0t+loXLlwweWxszGTfGZR2Gfl9uGnTJpNPnjwpAAAAAADmC2f8AAAAAAAA5BQLPwAAAAAAADnFwg8AAAAAAEBOlbTjZ+HChWpsbMyy76/x0p4X39PiO398x8/mzZtN9p0y6TwkqaurK9tOu2mmyjN1+hw6dMjk9vZ2kz/96U9n26+88ooZGx0dNXnFihUm+8/9zjvvmFxVVWVy2kdTX19vxjo6OkyemJgw2ffsDA4OFh33nT/pe/uOJc/v4wceeMBkv8+feOIJk998881s23fy+D4gv09979Tw8LDJvlfJdwilfKdP2j00MDAw7fMAAAAAAJgLnPEDAAAAAACQUyz8AAAAAAAA5BQLPwAAAAAAADlV0o6f8fHxSV07xSxevDjbDiFMeq3pHitJY2NjJj/++OMmL1hg17yeffbZbPuhhx4qOi//2r4DZtWqVSafOHHC5LT3xXff+B4e33Xz9a9/3WTffeP7adK5+q4b39Hj+4Xq6upMrqysNNn3DfnvNv2cfX19ZqympsZk39G0du1ak/1+On/+/LRz3bhxoxnzXUYLF9qffdpFJE3uMvI5/e35/X/jxg2T+/v7s+1i3UAAAAAAAMwFzvgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJwqacfPggULtGTJkiw3Nzeb8bNnz5qcdsjU1taasZ07d5rsO2F8j8u2bdtM9n01aR+Nf67vwnnrrbdMXrlypckPP/ywyb7bpbe3N9v2fULd3d0mHzp0yOTnn3/e5J/+9Kcmd3Z2mpx26XR1dZmxgYGBaR871Vy8hoYGk4eGhkxOe338a1dVVZnse5N8p4/vcPKdQWm3kZ+H//78e1+4cKHouO+TWrduXbbt+4PSMcl2U128eFEAAAAAAJQSZ/wAAAAAAADkFAs/AAAAAAAAOVXy27mnlzml25LU2Nho8ubNm7Ntf7tvf8v0Bx980GR/mdLTTz9t8n/913+ZvGXLlmz7+PHjZszf5txflnTq1CmTjxw5YrKf+/bt27NtfwmUz48++qjJ/nbtfq5+PL1VeXo5mzT5cjd/mZjn97l/L3+J1NKlS7Nt/7muXr067WOl25cFptLbokvS5cuXTU4vC/SX/bW2tpr8xhtvmNzW1mayv1X8ihUrTPbff8pf+lXOQggvSnpRmrzPAMwfjk2gPHFsAuWH4xKYHc74AT6iYowvxRh3xxh3+8UtAPOHYxMoTxybQPnhuARmh4UfAAAAAACAnGLhBwAAAAAAIKdK2vETQjBdL01NTWbc34Y7vYV7S0vLtGNTPdd34/jxTZs2mXzz5s1s+xe/+IUZ8x0/586dM9l3+Pjbt/t+m+XLl2fbFRUVZsx326xfv97kS5cumexvH+5v/57emt53KPnbu/uenWvXrpnse3Vmku63iYmJoo9tbm422Xc0pV1FklRfX29y+ls6efKkGfOfw98a3s/N90X5fZ7+9vxt5/17pb8NPwYAAAAAwFzjjB8AAAAAAICcYuEHAAAAAAAgp1j4AQAAAAAAyKmSdvxUVlZq1apV046nfTSS7eVJO3gkacOGDSZv377dZN+Nc/36dZN9v9D58+ez7bVr15qx/v5+k3/4wx+a3NPTY7LvyvGfeXx8fNp5P/fccyYfPnxYxTz++OMm+x6ZAwcOZNsnTpwwY21tbSb7fqDW1laTOzs7i87FS/ttent7zZi/3aLv+PHjg4ODJg8PD5vc0dGRbe/du9eM+e/edzL5cf9b831SaYeQ/x3533A6fuPGDQEAAAAAUEqc8QMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOVXSjp+RkZGes2fPdkhqktQz0+O7urqmHfv3f//3+zizzKzmNU9KPrdZdvpMO6/Lly9P+yT/3X7zm98smu/SnO8z3+/kuc6ldXM5FwAAAAAAvJIu/MQYV0hSCGFfjHF3Kd97Nsp1XlL5zq1c5yWV99wAAAAAACgFLvUCAAAAAADIKRZ+AAAAAAAAcmq+Fn5emqf3nUm5zksq37mV67yk8p4bAAAAAABzbl4WfmKMZfn/kJfrvKTynVu5zksq77kBAAAAAFAKXOoFAAAAAACQUyVd+AkhfDaEcCyEcDKE8LVSvvcUc/mnEEJ3COFQ8rfGEMKrIYQThX8vm4d5tYUQXg8hHAkhHA4hfLWM5lYTQvhVCOG9wtz+vPD39SGEvYW5/UsIoarUcyvMoyKE8G4I4fvlNC8AAAAAAOZLyRZ+QggVkv5B0m9LelDSl0IID5bq/afwz5I+6/72NUmvxRg3S3qtkEttTNIfxxgfkPSUpP9e2E/lMLdhSZ+KMT4s6RFJnw0hPCXpryT9TWFufZK+Mg9zk6SvSjqS5HKZFwAAAAAA86KUZ/w8IelkjPF0jHFE0rckfbGE72/EGN+QdNX9+YuSXi5svyzpd0s6KUkxxq4Y4zuF7eu6vZDRUiZzizHGwUKsLPwTJX1K0r/N59xCCK2SPi/pHws5lMO8AAAAAACYT6Vc+GmRdD7JnYW/lZNVMcYu6fYCjKSV8zmZEEK7pF2S9qpM5la4nOqApG5Jr0o6Jak/xjhWeMh8fa9/K+lPJE0U8vIymRcAAAAAAPOmlAs/YYq/xRK+/4dKCGGJpG9L+qMY48B8z+cDMcbxGOMjklp1+yyuB6Z6WCnnFEL4HUndMcb96Z+neCi/NwAAAADAR8rCEr5Xp6S2JLdKuljC95+NyyGE5hhjVwihWbfPaim5EEKlbi/6fCPG+J1ymtsHYoz9IYQ9ut1D1BBCWFg4u2Y+vtdnJH0hhPA5STWS6nX7DKD5nhcAAAAAAPOqlGf8vC1pc+FOS1WSfk/S90r4/rPxPUkvFLZfkPTdUk+g0E3zdUlHYox/XWZzWxFCaChs10r6tG53EL0u6fn5mluM8c9ijK0xxnbd/l39JMb45fmeFwAAAAAA861kCz+Fsy7+h6Qf6fZiwb/GGA+X6v29EMI3Jf1S0tYQQmcI4SuS/lLSZ0IIJyR9ppBL7RlJ/03Sp0IIBwr/fK5M5tYs6fUQwkHdXsh7Ncb4fUl/Kul/hhBO6na3ztfnYW5TKdd5AQAAAABQEqW81Esxxh9I+kEp33M6McYvTTP0GyWdiBNj/Lmm7qeR5n9uB3W7bNr//bRu9/3MuxjjHkl7CttlMy8AAAAAAOZDKS/1AgAAAAAAQAmx8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAOcXCDwAAAAAAQE6x8AMAAAAAAJBTLPwAAAAAAADkFAs/AAAAAAAAORVijPM9BwDzLIRwRVKHpCZJPfM8namU67yk8p1buc5LKr+5rYsxrpjvSUyFY/OelOvcmNfslfuxeUPlt88+UI7fp1S+85LKd27lOK+yPDY/BP+bKZXv3Mp1XlL5zq0c5zXtscnCD4BMCGFfjHH3fM/DK9d5SeU7t3Kdl1TecytX5brPynVeUvnOjXnlRznvs3KdW7nOSyrfuZXrvMpZOe+zcp1buc5LKt+5leu8psOlXgAAAAAAADnFwg8AAAAAAEBOsfADIPXSfE9gGuU6L6l851au85LKe27lqlz3WbnOSyrfuTGv/CjnfVaucyvXeUnlO7dynVc5K+d9Vq5zK9d5SeU7t3Kd15To+AEAAAAAAMgpzvgBAAAAAADIKRZ+ACiE8NkQwrEQwskQwtfmeS7/FELoDiEcSv7WGEJ4NYRwovDvZfMwr7YQwushhCMhhMMhhK+W0dxqQgi/CiG8V5jbnxf+vj6EsLcwt38JIVSVem6FeVSEEN4NIXy/nOb1YcCxOat5leWxWe7HZWEuHJt3iWNzVvPi2Lz7OXJs3iWOzVnNi2Pz7ub3oT4uWfgBPuJCCBWS/kHSb0t6UNKXQggPzuOU/lnSZ93fvibptRjjZkmvFXKpjUn64xjjA5KekvTfC/upHOY2LOlTMcaHJT0i6bMhhKck/ZWkvynMrU/SV+ZhbpL0VUlHklwu8yprHJuzVq7HZrkflxLH5l3h2Jw1js27x7F5Fzg2Z41j8+58qI9LFn4APCHpZIzxdIxxRNK3JH1xviYTY3xD0lX35y9Kermw/bKk3y3ppCTFGLtijO8Utq/r9v/hbymTucUY42AhVhb+iZI+Jenf5nNuIYRWSZ+X9I+FHMphXh8SHJuzUK7HZjkflxLH5j3i2JwFjs27w7F5Tzg2Z4Fj887l4bhk4QdAi6TzSe4s/K2crIoxdkm3/8dK0sr5nEwIoV3SLkl7VSZzK5x+ekBSt6RXJZ2S1B9jHCs8ZL6+17+V9CeSJgp5eZnM68OAY/MOlduxWcbHpcSxeS84Nu8Qx+Yd4di8exybd4hjc9Y+9MclCz8AwhR/43Z/0wghLJH0bUl/FGMcmO/5fCDGOB5jfERSq27/F68HpnpYKecUQvgdSd0xxv3pn6d4KL+3qbGv7kA5HpvleFxKHJv3AfvqDnBszh7H5j1jX90Bjs3ZyctxuXC+JwBg3nVKaktyq6SL8zSX6VwOITTHGLtCCM26/V8BSi6EUKnb/wP5jRjjd8ppbh+IMfaHEPbo9jXbDSGEhYX/GjEf3+szkr4QQvicpBpJ9br9X0zme14fFhybs1Tux2aZHZcSx+a94ticJY7NO8axeW84NmeJY/OO5OK45IwfAG9L2lxopq+S9HuSvjfPc/K+J+mFwvYLkr5b6gkUruX9uqQjMca/LrO5rQghNBS2ayV9Wrev135d0vPzNbcY45/FGFtjjO26/bv6SYzxy/M9rw8Rjs1ZKNdjs1yPS4lj8z7g2JwFjs07x7F5zzg2Z4Fj887k5riMMfIP//DPR/wfSZ+TdFy3r6P9v+Z5Lt+U1CVpVLf/y81XdPs62tcknSj8u3Ee5vVx3T6F86CkA4V/Plcmc9sp6d3C3A5J+r8Lf98g6VeSTkr635Kq5/F7/aSk75fbvMr9H47NWc2rLI/ND8NxWZgPx+bd7TeOzZnnxbF5b/Pk2Ly7/caxOfO8ODbvfo4f2uMyFCYNAAAAAACAnOFSLwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcoqFHwAAAAAAgJxi4QcAAAAAACCnWPgBAAAAAADIKRZ+AAAAAAAAcur/BzAn303GwscNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x360 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "filename = 'example_data_2D/drawn_fibres_B.png';\n",
    "sigma = 0.5\n",
    "rho = 2\n",
    "\n",
    "image = skimage.io.imread(filename)\n",
    "S = st2d.structure_tensor(image, sigma, rho)\n",
    "val,vec = st2d.eig_special(S)\n",
    "print(f'The image has a shape {image.shape}, i.e. {image.size} pixels.')\n",
    "print(f'Structure tensor information is carried in a {S.shape} array.')\n",
    "print(f'Orientation information is carried in a {vec.shape} array.')\n",
    "print(f'Anisotropy information is carried in a {val.shape} array.')\n",
    "\n",
    "# visualization\n",
    "figsize = (20,5)\n",
    "fig, ax = plt.subplots(1, 5, figsize=figsize, sharex=True, sharey=True)\n",
    "ax[0].imshow(image, cmap=plt.cm.gray)\n",
    "plot_orientations(ax[0], image.shape, vec)\n",
    "ax[0].set_title('Orientation as arrows')\n",
    "orientation_st_rgba = plt.cm.hsv((np.arctan2(vec[1], vec[0])/np.pi).reshape(image.shape))\n",
    "ax[1].imshow(plt.cm.gray(image)*orientation_st_rgba)\n",
    "ax[1].set_title('Orientation as color on image')\n",
    "ax[2].imshow(orientation_st_rgba)\n",
    "ax[2].set_title('Orientation as color')\n",
    "anisotropy = (1-val[0]/val[1]).reshape(image.shape)\n",
    "ax[3].imshow(anisotropy)\n",
    "ax[3].set_title('Degree of anisotropy')\n",
    "ax[4].imshow(plt.cm.gray(anisotropy)*orientation_st_rgba)\n",
    "ax[4].set_title('Orientation and anisotropy')\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2\n",
    "Orientation analysis of a slice showing carbord fibres. This example shows how orientation may be collected and shown as a histogram over angles or as a polar histogram. \n",
    "\n",
    "Here I collect the orientation information on the whole image. However, as shown above, the orientation information is reliable only in areas of hight anisotropy. For this reason, and depending on problem at hand, it may be desirable to weight the orientation using the anisotropy (remove parts with low anisotropy) or the intensity (remove the background and keep only information of fibres). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'example_data_2D/10X.png';\n",
    "sigma = 0.5\n",
    "rho = 15\n",
    "N = 180 # number of angle bins for orientation histogram\n",
    "\n",
    "# computation\n",
    "image = skimage.io.imread(filename)\n",
    "image = np.mean(image[:,:,0:3],axis=2).astype(np.uint8)\n",
    "S = st2d.structure_tensor(image, sigma, rho)\n",
    "val,vec = st2d.eig_special(S)\n",
    "angles = np.arctan2(vec[1], vec[0]) # angles from 0 to pi\n",
    "distribution = np.histogram(angles, bins=N, range=(0.0, np.pi))[0]\n",
    "\n",
    "# visualization\n",
    "figsize = (10,5)\n",
    "fig, ax = plt.subplots(1, 2, figsize=figsize, sharex=True, sharey=True)\n",
    "ax[0].imshow(image,cmap=plt.cm.gray)\n",
    "ax[0].set_title('Input image')\n",
    "orientation_st_rgba = plt.cm.hsv((angles/np.pi).reshape(image.shape))\n",
    "ax[1].imshow(plt.cm.gray(image)*orientation_st_rgba)\n",
    "ax[1].set_title('Orientation as color on image')\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=figsize)\n",
    "bin_centers = (np.arange(N)+0.5)*np.pi/N # halp circle (180 deg)\n",
    "colors = plt.cm.hsv(bin_centers/np.pi)\n",
    "ax[0].bar(bin_centers, distribution, width = np.pi/N, color = colors)\n",
    "ax[0].set_xlabel('angle')\n",
    "ax[0].set_xlim([0,np.pi])\n",
    "ax[0].set_aspect(np.pi/ax[0].get_ylim()[1])\n",
    "ax[0].set_xticks([0,np.pi/2,np.pi])\n",
    "ax[0].set_xticklabels(['0','pi/2','pi'])\n",
    "ax[0].set_ylabel('count')\n",
    "ax[0].set_title('Histogram over angles')\n",
    "polar_histogram(ax[1], distribution)\n",
    "ax[1].set_title('Polar histogram')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OCT example\n",
    "A example of of computing orientation information from the OCT image of retina."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'example_data_2D/OCT_im_org.png';\n",
    "sigma = 0.5\n",
    "rho = 5\n",
    "N = 180 # number of angle bins for orientation histogram\n",
    "\n",
    "# computation\n",
    "image = skimage.io.imread(filename)\n",
    "S = st2d.structure_tensor(image, sigma, rho)\n",
    "val,vec = st2d.eig_special(S)\n",
    "angles = np.arctan2(vec[1], vec[0]) # angles from 0 to pi\n",
    "distribution = np.histogram(angles, bins=N, range=(0.0, np.pi))[0]\n",
    "\n",
    "# visualization\n",
    "figsize = (10,5)\n",
    "fig, ax = plt.subplots(1, 2, figsize=figsize, sharex=True, sharey=True)\n",
    "ax[0].imshow(image, cmap=plt.cm.gray)\n",
    "ax[0].set_title('Input image')\n",
    "orientation_st_rgba = plt.cm.hsv((angles/np.pi).reshape(image.shape))\n",
    "ax[1].imshow(plt.cm.gray(image)*orientation_st_rgba)\n",
    "ax[1].set_title('Orientation as color on image')\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=figsize)\n",
    "bin_centers = (np.arange(N)+0.5)*np.pi/N # halp circle (180 deg)\n",
    "colors = plt.cm.hsv(bin_centers/np.pi)\n",
    "ax[0].bar(bin_centers, distribution, width = np.pi/N, color = colors)\n",
    "ax[0].set_xlabel('angle')\n",
    "ax[0].set_xlim([0,np.pi])\n",
    "ax[0].set_aspect(np.pi/ax[0].get_ylim()[1])\n",
    "ax[0].set_xticks([0,np.pi/2,np.pi])\n",
    "ax[0].set_xticklabels(['0','pi/2','pi'])\n",
    "ax[0].set_ylabel('count')\n",
    "ax[0].set_title('Histogram over angles')\n",
    "polar_histogram(ax[1], distribution)\n",
    "ax[1].set_title('Polar histogram')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating the effect of integration size rho\n",
    "Integration size rho is a crucial parameter when working with structure tensor, and especially if ofientation information is to be sampled om the image. This example shows how small rho captures the randomness of the individual fibre, and equidistantly placed orientation arrows will appear random. Larger rho captures the orientation on a larger scale, and flow appears smoother. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'example_data_2D/short_fibres.png'\n",
    "image = skimage.io.imread(filename)\n",
    "image = np.mean(image[:,:,0:3],axis=2)\n",
    "image -= np.min(image)\n",
    "image /= np.max(image)\n",
    "s = 128 # quiver arrow spacing\n",
    "sigma = 0.5\n",
    "figsize = (20,10)\n",
    "\n",
    "rhos = [2,10,20,50]\n",
    "\n",
    "fig, ax = plt.subplots(2, len(rhos), figsize=figsize, sharex=True, sharey=True)\n",
    "for k in range(4):\n",
    "\n",
    "    # computation\n",
    "    rho = rhos[k] # changing integration radius\n",
    "    S = st2d.structure_tensor(image,sigma,rho)\n",
    "    val,vec = st2d.eig_special(S)\n",
    "\n",
    "    # visualization\n",
    "    ax[0][k].imshow(image,cmap=plt.cm.gray)\n",
    "    plot_orientations(ax[0][k], image.shape, vec, s = s)\n",
    "    ax[0][k].set_title(f'Rho = {rho}, arrows')\n",
    "    intensity_rgba = plt.cm.gray(image)\n",
    "    orientation_st_rgba = plt.cm.hsv((np.arctan2(vec[1], vec[0])/np.pi).reshape(image.shape))\n",
    "    ax[1][k].imshow((0.5+0.5*intensity_rgba)*orientation_st_rgba)\n",
    "    ax[1][k].set_title(f'Rho = {rho}, color')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating the effect of downsampling\n",
    "This example demonstrates the robustness of structure tensor to downsampling. In highest resolution, individual fibres can be distinguished, and large rho is used to capture a smooth flow. As smaller resolutions, individual fibres cannot be distinguished, but the orientation information, now using a smaller rho, is still rather stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'example_data_2D/short_fibres.png'\n",
    "downsampling_range = 4\n",
    "figsize = (20,10)    \n",
    "fig, ax = plt.subplots(2, downsampling_range, figsize=figsize)\n",
    "\n",
    "for k in range(downsampling_range):\n",
    "\n",
    "    # downsampling and computation\n",
    "    scale = 2**k\n",
    "    f = 1/scale # image scale factor\n",
    "    s = 128//scale # quiver arrow spacing\n",
    "    sigma = 0.5 # would it make sense to scale this too?\n",
    "    rho = 50/scale # scaling the integration radius\n",
    "    image = skimage.io.imread(filename)\n",
    "    image = np.mean(image[:,:,0:3],axis=2)\n",
    "    image = skimage.transform.rescale(image,f,multichannel=False)\n",
    "    image -= np.min(image)\n",
    "    image /= np.max(image)\n",
    "    S = st2d.structure_tensor(image,sigma,rho)\n",
    "    val,vec = st2d.eig_special(S)\n",
    "\n",
    "    # visualization\n",
    "    ax[0][k].imshow(image,cmap=plt.cm.gray)\n",
    "    plot_orientations(ax[0][k], image.shape, vec, s = s)\n",
    "    ax[0][k].set_title(f'Downsample = {scale}, arrows')\n",
    "    intensity_rgba = plt.cm.gray(image)\n",
    "    orientation_st_rgba = plt.cm.hsv((np.arctan2(vec[1], vec[0])/np.pi).reshape(image.shape))\n",
    "    ax[1][k].imshow((0.5+0.5*intensity_rgba)*orientation_st_rgba)       \n",
    "    ax[1][k].set_title(f'Downsample = {scale}, color')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing the dominant orientation with the optical flow\n",
    "When imaged structures are mostly unidirectional, and we are interested in computing the deviation from this given direction, it may be favorable to compute optical flow. While orientation from structure tensor yields a $(x,y)$ vector per pixel, optical flow in $y$ direction yields a $(x,1)$ vector. However, where assumption on uni-directionality is broken, optical flow is still returning a flow in $y$ direction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = skimage.io.imread('example_data_2D/drawn_fibres_B.png');\n",
    "\n",
    "# computing structure tensor, orientation and optical flow\n",
    "sigma = 0.5\n",
    "rho = 5\n",
    "S = st2d.structure_tensor(image,sigma,rho)\n",
    "val,vec = st2d.eig_special(S) # dominant orientation\n",
    "fx = st2d.solve_flow(S) # optical flow\n",
    "\n",
    "# visualization\n",
    "figsize = (10,10)\n",
    "fy = np.ones(image.shape)\n",
    "fig, ax = plt.subplots(2,2,figsize=figsize)\n",
    "\n",
    "ax[0][0].imshow(image,cmap=plt.cm.gray)\n",
    "plot_orientations(ax[0][0], image.shape, vec)\n",
    "ax[0][0].set_title('Orientation from structure tensor, arrows')\n",
    "ax[0][1].imshow(image,cmap=plt.cm.gray)\n",
    "plot_orientations(ax[0][1], image.shape, np.r_[fx,np.ones((1,image.size))])\n",
    "ax[0][1].set_title('Orientation from optical flow, arrows')\n",
    "intensity_rgba = plt.cm.gray(image)\n",
    "orientation_st_rgba = plt.cm.hsv((np.arctan2(vec[1],vec[0])/np.pi).reshape(image.shape))\n",
    "orientation_of_rgba = plt.cm.hsv((np.arctan2(1,fx)/np.pi).reshape(image.shape))\n",
    "ax[1][0].imshow(intensity_rgba*orientation_st_rgba)\n",
    "ax[1][0].set_title('Dominant orientation from structure tensor, color')\n",
    "ax[1][1].imshow(intensity_rgba*orientation_of_rgba)\n",
    "ax[1][1].set_title('Orientation from optical flow, color')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}