{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from scipy.interpolate import RectBivariateSpline\n", "from skimage.data import shepp_logan_phantom\n", "from skimage.transform import radon, rescale, rotate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Analytic Reconstruction Methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the following setup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An object with unknown attenuation coefficient $\\mu(x,y)$ can be scanned in multiple positions. During scanning, photons are sent through the object and the resulting intensity is picked up by a scanner. It is known that\n", "\n", "$$I = I_0 e^{-\\int \\mu(x,y)ds}$$\n", "\n", "where $ds$ is a particular straight line through the object. If is more useful to write\n", "\n", "$$\\ln(I_0/I) = \\int \\mu(x,y)ds$$ \n", "\n", "The quantity on the left can be measured for a variety of scanner positions. Shown above are two possible scanner positions with corresponding values of $\\ln(I_0/I)$. \n", "\n", "For the rest of this we will call \n", "\n", "* $\\mu(x,y)=f(x,y)$ \n", "* $\\ln(I_0/I)=p(r,\\theta)$.\n", "\n", "The parameters are as follows:\n", "\n", "* The parameter $\\theta$ species the tilt of the scanner (i.e. scanner 1 at $\\theta=90^{\\circ}$ and scanner 2 at $\\theta=-45^{\\circ})$. \n", "* The parameter $r$ specifies the horizontal location on the the scanner (i.e. scanner 1 measures largest value at $r=0$ and decays as $r\\to \\pm$)\n", "\n", "**The goal is to reconstruct $f(x,y)$ given $p(r, \\theta)$**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll start with $f(x,y)$, use this to get $p(r,\\theta)$, and then look at techniques to reobtain $f(x,y)$ from $p(r,\\theta)$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#image = shepp_logan_phantom()\n", "image = np.ones([100,100])\n", "# Resize Image\n", "diag = len(np.diag(image)//2)\n", "image = np.pad(image, pad_width=diag+10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a meshgrid of $x$ and $y$ values that correspond to coordinates of the square. $x=0$ and $y=0$ correspond to the point of rotation" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "_ = np.linspace(-1, 1, image.shape[0])\n", "xv, yv = np.meshgrid(_,_)\n", "image[(xv-0.1)**2+(yv-0.2)**2<0.01] = 2\n", "\n", "# Create a rotated image\n", "image_rot = rotate(image, 45)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the square and its rotation" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAEyCAYAAACh0Ed7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdDUlEQVR4nO3df5Ac5Z3f8c9nVwgbLAUEBgTINsgb50QSy7IMcvCdwbI40B2RSSUXqBhRV3bJJKcqi6RSpSpXOfen48v5KFcwWE7IQSoHh7Exsi0sIcUp3XHRRYsiQAvoJGTuWFZBsUxAgK0f7Dd/bC83WmZ2fnTPzDPd71fV1k53P93zNC2++uh5pqcdEQIAAEBahvrdAQAAALwbIQ0AACBBhDQAAIAEEdIAAAASREgDAABIECENAAAgQYWENNv32j5ie1+D7bb9TdsHbT9te1nNtutt78+2bSyiPwDQKuoXgFQVNZL2x5Kun2X7DZJGsp91ku6WJNvDku7Kti+RdIvtJQX1CQBa8ceifgFIUCEhLSJ2SvrFLE3WSLo/puySdI7thZKulHQwIg5FxAlJD2ZtAaAnqF8AUjWnR+9ziaSXapbHs3X11l9V7wC212nqX7Ea1vDHz9L87vQUQJKO6dWfR8T7+/DW1C8AuXRav3oV0lxnXcyy/t0rIzZJ2iRJ870grvLK4noHIHnb4+G/7tNbU78A5NJp/epVSBuXtKhm+VJJE5LmNlgPAKmgfgHoi159BcdmSWuzu6RWSHotIg5L2i1pxPZltudKujlrCwCpoH4B6ItCRtJsPyDpGknn2x6X9O8knSFJEXGPpC2SVks6KOktSb+bbTtle72krZKGJd0bEWNF9AkAWkH9ApCqQkJaRNzSZHtI+r0G27ZoqggCQM9RvwCkiicOAAAAJIiQBgAAkCBCGgAAQIIIaQAAAAkipAEAACSIkAYAAJAgQhoAAECCCGkAAAAJIqQBAAAkiJAGAACQIEIaAABAgghpAAAACSKkAQAAJIiQBgAAkCBCGgAAQIIIaQAAAAkipAEAACSIkAYAAJAgQhoAAECCCGkAAAAJIqQBAAAkiJAGAACQIEIaAABAgghpAAAACSKkAQAAJKiQkGb7etv7bR+0vbHO9n9re2/2s8/227YXZNtetP1Mtm20iP4AQDuoYQBSNCfvAWwPS7pL0ipJ45J2294cEc9Ot4mIP5D0B1n7GyXdERG/qDnMtRHx87x9AYB2UcMApKqIkbQrJR2MiEMRcULSg5LWzNL+FkkPFPC+AFAEahiAJBUR0i6R9FLN8ni27l1snyXpeknfq1kdkrbZftL2ukZvYnud7VHboyd1vIBuA4CkHtQw6heATuSe7pTkOuuiQdsbJT0xY5rg6oiYsH2BpMdtPx8RO991wIhNkjZJ0nwvaHR8AGhX12sY9QtAJ4oYSRuXtKhm+VJJEw3a3qwZ0wQRMZH9PiLpEU1NPQBAr1DDACSpiJC2W9KI7ctsz9VUEds8s5HtvyPp05IerVl3tu15068lXSdpXwF9AoBWUcMAJCn3dGdEnLK9XtJWScOS7o2IMdu3Z9vvyZreJGlbRLxZs/uFkh6xPd2XP4mIn+TtEwC0ihoGIFWOGLyPR8z3grjKK/vdDQA9tD0efjIilve7H3lRv4Dq6bR+8cQBAACABBHSAAAAEkRIAwAASBAhDQAAIEGENAAAgAQR0gAAABJESAMAAEgQIQ0AACBBhDQAAIAEEdIAAAASREgDAABIECENAAAgQYQ0AACABBHSAAAAEkRIAwAASBAhDQAAIEGENAAAgAQR0gAAABJESAMAAEgQIQ0AACBBhDQAAIAEEdIAAAASREgDAABIECENAAAgQYQ0AACABBUS0mxfb3u/7YO2N9bZfo3t12zvzX6+2uq+ANBt1DAAKZqT9wC2hyXdJWmVpHFJu21vjohnZzT9s4j47Q73BYCuoIYBSFURI2lXSjoYEYci4oSkByWt6cG+AFAEahiAJBUR0i6R9FLN8ni2bqZP2n7K9mO2r2hzX9leZ3vU9uhJHS+g2wAgqQc1jPoFoBO5pzsluc66mLG8R9IHI+IN26sl/UDSSIv7Tq2M2CRpkyTN94K6bQCgA12vYdQvAJ0oYiRtXNKimuVLJU3UNoiI1yPijez1Fkln2D6/lX0BoMuoYQCSVERI2y1pxPZltudKulnS5toGti+y7ez1ldn7Hm1lXwDoMmoYgCTlnu6MiFO210vaKmlY0r0RMWb79mz7PZL+qaR/afuUpF9KujkiQlLdffP2CQBaRQ0DkCpP1ZnBMt8L4iqv7Hc3APTQ9nj4yYhY3u9+5EX9Aqqn0/rFEwcAAAASREgDAABIECENAIAuG1p2RfNGwAyENAAAumho2RU6sHYeQQ1tI6QBANAl0wFNEkENbSviiQMAAGCGF+5c8a51B9bOk9au0Mj9xzS5h29rwewYSQMAoGDNRsymR9eA2RDSAAAoUO0UZ7N2wGwIaQAAFKTVgCbxGTU0R0gDAKAg7U5jEtQwG0IaAAAF6DRsEdTQCHd3AgCQQztTnI0cWDtPI7qCOz5xGkbSAADoUBEBbRojapiJkAYAQIeK/ioNghpqEdIAAOhAt8IUQQ3TCGkAALSpyGnOeviyW0iENAAA2tLtgFb7Pqg27u4EAKBF9Z7H2S085xOMpAEA0IJ+jWwx9VldhDQAAJro1RTnbO+P6mG6E6fp5VB+2SzesKvfXQDQBf0OaBJfdltVjKQBANBACgFtGl/NUT2MpKE0Dv7Ot1tq9+GHvtTlngAog5QC2rTpmwkYua8GRtJQCq0GtHbbAqimFANaLUbUqoGRNAy0TgPX9H6MqgGYKfWAJvEZtaooZCTN9vW299s+aHtjne3/wvbT2c9f2P5ozbYXbT9je6/t0SL6g2ooYkSMUTVI1DD8rUEIaNP4jFr55R5Jsz0s6S5JqySNS9pte3NEPFvT7GeSPh0Rr9q+QdImSVfVbL82In6ety+oDsIVikINQ61BCWjTGFErtyJG0q6UdDAiDkXECUkPSlpT2yAi/iIiXs0Wd0m6tID3BQpB4Ks8ahgkDe7nvAYtWKJ1RYS0SyS9VLM8nq1r5AuSHqtZDknbbD9pe12jnWyvsz1qe/SkjufqMAZbN0IVQa3Sul7DqF9pG1p2hV64c8VAh50X7lwxsCETjRVx44DrrIu6De1rNVXgPlWz+uqImLB9gaTHbT8fETvfdcCITZqaYtB8L6h7fJQfYQpd0PUaRv1K1yB9Bq0Zpj7Lp4iRtHFJi2qWL5U0MbOR7X8o6T9JWhMRR6fXR8RE9vuIpEc0NfUA9BwBsLKoYRVWloA2jZsJyqWIkLZb0ojty2zPlXSzpM21DWx/QNL3Jd0aEX9Vs/5s2/OmX0u6TtK+AvoEAK2ihlVUWcNM2YJnleWe7oyIU7bXS9oqaVjSvRExZvv2bPs9kr4q6TxJ37ItSaciYrmkCyU9kq2bI+lPIuInefuEcmKkC91ADaumMk1z1jO0jGnPMijky2wjYoukLTPW3VPz+ouSvlhnv0OSPjpzPdAvB3/n23zBbQVRw6ql7AFN4vNpZcFjoQAAlVGFgDbtwNp53PU54AhpAIBKqFJAq1XFcy4LQhoAoPSqGtCmMZo2mAhpAIBSq3pAk/hqjkFFSAMAlFrVA9o0gtrgKeTuTqAXPvzQl7r+NRzc2QmURxVG0IYu+qV+a2RM4788V9+9fLtenXxL5w6dpTcmf6VJhb7w4o06b+6b2rb7H2jopKf+e6xdocUbdvW762gBIQ0AUDplD2hzLn5Lk5PWzqu/pYuG3/fO+nOHzpIkvW/oPZKk716+XZK0+/079c9/eruGXp/6a5/vURsMTHdioHRzpItRNKAcyh7QJOn5T/1XjX7q26cFtNl84swz9KfX3qPzRo7qg39/gqnPAUFIAwCUStkD2q+veFaSNH/ovW3t94kzz9D//NhD2rHkh5L4jNogIKRh4HRjxItRNKAcyh46fn3Fs/ovH/izjvcfyv7a/8jSv5FU/kA76AhpqDwCGlAOZZ/mzBvQav347z72TlAre7AdZNw4gIE0Hazy3u1JQAMGX9nDmSQ9sOY/6hNnnlHoMe++/CFdc2jDO3d8jtx/jJsJEsNIGgZanpBFQAPKoewBTZIWzTle+DE/MGee/vTae95Z5jNq6SGkYeB1ErYIaEA5VCFU+MJftXwXZ7tmjs5VIfAOEqY7UQq1oavRFCjBDCiXKkxzStLwnMmevh/foZYOQhpKhzAGlF9VApokff4ju7t6/LmXvqkT42e/s3xg7TyNiKCWAqY7AQADpUoBbfK9k7rlnNGevy+fT0sDIQ0AMDCqFNAkyceHdGyyu5Ne3/jYQ3XXH1g7Ty/cuaKr743ZEdIAAAOhagFNkjwpvd3lv6pPxuwhkBG1/iGkAQCSV8WAJkmT55zSx+fO7ep7XHXmEU2+t/HNCUx99g8hDQCQtKoGNEmKE0P6m1PHuvoeh069RzoxexwgqPUHIQ0AkLSqBjRJGn5rSBNvt/cg9XZ958inNfR283YEtd4jpAEAkjS07Ao+uC7p8098savH/+neJS23Jaj1FiENAJCcKk9xznTeuW909fhDv2wvChDUeoeQBgBIDgHtbx09cF6/u/AuBLXeKCSk2b7e9n7bB21vrLPdtr+ZbX/a9rJW9wWAbqOGpYW//N/te292J7SufPbGjvclqHVf7pBme1jSXZJukLRE0i22Z05w3yBpJPtZJ+nuNvYFgK6hhqWFac765rqFT/a36eW3j+nQzy7MdQyuVXcVMZJ2paSDEXEoIk5IelDSmhlt1ki6P6bsknSO7YUt7gsA3UQNSwQBrbEvb/u8fvjWWYUe84dv/D0Nvzmc+ziMpnVPESHtEkkv1SyPZ+taadPKvpIk2+tsj9oePanjuTsNAJmu1zDqV3MEtNkNnbDu+NGtuvwH6woJax9+6Ev6D1t/u4CeMe3ZTUWENNdZFy22aWXfqZURmyJieUQsP0NnttlFAGio6zWM+tXc5J4xjdzf3S9tLYOhE9aXt30+1zHuea3uWEjHRu4/psk9Y4UeE1OKeGrruKRFNcuXSpposc3cFvYFgG6ihiVics+YRsSIWjNDJ6x/9fJVevXE2Xrgsv/e8n5rDvymLjnr/2nbE0sL7Q8BrXuKGEnbLWnE9mW250q6WdLmGW02S1qb3SG1QtJrEXG4xX0BoJuoYQlhRK01255Yqt27R3TLzz6jF5s8NmrviV9p5bM3aux/f6jwgMa16q7cI2kRccr2eklbJQ1Lujcixmzfnm2/R9IWSaslHZT0lqTfnW3fvH0CgFZRw9IzuWdMWsuTBlqxe/eIPvPMv9Ydv7FVr546W//n+Hwtes+revn4OTpr+KQWnXlU33jiukJuEJiJac7uc0Tdj4Albb4XxFVe2e9ulBKPYOnc4g27+t2FUtseDz8ZEcv73Y+8qF+t42aCNBHO2tdp/eKJAwCAJDH1mSYCWu8Q0gAAySIQpIXQ3FuENABA0ggGaWCas/cIaQCApDHt2X8EtP4gpAEAkkdQ6x8CWv8Q0gAAA2Fyzxh3UfcYAa2/CGkAgIHCiFpvEND6j5AGABgoTH12HwEtDYQ0AMDAIah1DwEtHYQ0AMBAIqh1BwEtHYQ0AMDAIlAUi9CbltwPWAcAoJ8Wb9jFcz5zmp7inOx3R3AaRtIAAAOPqc/O8Rm0dBHSAAClQFDrDAEtXYQ0AEBpEDjaQ6hNGyENAFAqBI/WMM2ZPm4cAACUyuSeMS3eI24mmMXiDbu4SWAAMJIGACglPqNWH/9NBgchDQBQWkznnY4pzsFCSAMAlBojR1MIaIOHkAYAKDWmPQlog4qQBgAovaoHNQLaYOLuTgBAJUzf9fnCnSv63ZWeYQRtsDGSBgColKqMqBHQBh8hDQBQKVWY+iSglUOukGZ7ge3HbR/Ifp9bp80i2z+1/ZztMdtfrtn2+7Zftr03+1mdpz8A0A5qWHWVPagR0Moh70jaRkk7ImJE0o5seaZTkv5NRPyapBWSfs/2kprtfxQRS7OfLTn7AwDtoIZVWFmDTJnDZ9XkDWlrJN2Xvb5P0udmNoiIwxGxJ3t9TNJzki7J+b4AUARqWMWVLdAwzVkueUPahRFxWJoqZJIumK2x7Q9J+pikv6xZvd7207bvrTfVULPvOtujtkdP6njObgOApB7VMOpXuib3jGnxhl0DH9ZG7j829TxOAlqpNA1ptrfb3lfnZ007b2T7fZK+J2lDRLyerb5b0mJJSyUdlvSHjfaPiE0RsTwilp+hM9t5awAVlkINo36lb9A/o0Y4K6em35MWEZ9ttM32K7YXRsRh2wslHWnQ7gxNFbf/FhHfrzn2KzVtviPpR+10HgCaoYahVZN7xqS1g/cdaiP3H9NkvzuBrsg73blZ0m3Z69skPTqzgW1L+s+SnouIb8zYtrBm8SZJ+3L2BwDaQQ3DaQZtNI3PoJVb3pD2NUmrbB+QtCpblu2LbU/f5XS1pFslfabObepft/2M7aclXSvpjpz9AYB2UMNwmkGa9iSglV+ux0JFxFFJK+usn5C0Onv955LcYP9b87w/AORBDUM904+PGlp2hQ6sndfv7tS1eMMupjgrgCcOAABQR6qjVIMy0of8CGkAADSQWiBiirNaCGkAADSQ0mfUCGjVk+szaSifxRt29bsLAJCUyT1jGlF/P59GQKsmRtIAAGii3yNqBLRqYiQNAIAW9GNEjRG0amMkDQCAFvVyRI2ABkIaAABt6FVQI6CBkAYAQJu6HdRSuaMU/UVIAwCgA90a6WKaE9MIaQAAdKjoES8CGmpxdycAAB0q8jmfPI8TMzGSBgBATnk/o8Zn0FAPIQ0AgAJ0Ok3JFCcaIaQBAFCQdkfECGiYDSENAICCtDPtSUBDM4Q0AAAKNLlnTIs37Jq1DQENrSCkAQDQBY1G1AhoaBUhDQCALqg39UlAQzsIaQAAdEltUCOgoV18mS0AAF00uWdMI7qCgIa2MZIGAECXEdDQCUIaAABAgghpAAAACcoV0mwvsP247QPZ73MbtHvR9jO299oebXd/AOgGahiAlOUdSdsoaUdEjEjakS03cm1ELI2I5R3uDwBFo4YBSFbekLZG0n3Z6/skfa7H+wNAHtQwAMnKG9IujIjDkpT9vqBBu5C0zfaTttd1sD8AdAM1DECymn5Pmu3tki6qs+krbbzP1RExYfsCSY/bfj4idraxv7LCuE6S3qOz2tkVQIWlUMOoXwA60TSkRcRnG22z/YrthRFx2PZCSUcaHGMi+33E9iOSrpS0U1JL+2f7bpK0SZLme0E06zcASGnUMOoXgE7kne7cLOm27PVtkh6d2cD22bbnTb+WdJ2kfa3uDwBdRA0DkKy8Ie1rklbZPiBpVbYs2xfb3pK1uVDSn9t+StL/kvTjiPjJbPsDQI9QwwAkK9ezOyPiqKSVddZPSFqdvT4k6aPt7A8AvUANA5AynjgAAACQIEIaAABAgghpAAAACSKkAQAAJIiQBgAAkCBCGgAAQIIIaQAAAAkipAEAACSIkAYAAJAgQhoAAECCCGkAAAAJIqQBAAAkiJAGAACQIEIaAABAgghpAAAACSKkAQAAJIiQBgAAkCBCGgAAQIIIaQAAAAkipAEAACSIkAYAAJAgQhoAAECCCGkAAAAJIqQBAAAkiJAGAACQoFwhzfYC24/bPpD9PrdOm4/Y3lvz87rtDdm237f9cs221Xn6AwDtoIYBSFnekbSNknZExIikHdnyaSJif0QsjYilkj4u6S1Jj9Q0+aPp7RGxJWd/AKAd1DAAycob0tZIui97fZ+kzzVpv1LSCxHx1znfFwCKQA0DkKy8Ie3CiDgsSdnvC5q0v1nSAzPWrbf9tO176001AEAXUcMAJKtpSLO93fa+Oj9r2nkj23Ml/WNJ361ZfbekxZKWSjos6Q9n2X+d7VHboyd1vJ23BlBhKdQw6heATsxp1iAiPttom+1XbC+MiMO2F0o6MsuhbpC0JyJeqTn2O69tf0fSj2bpxyZJmyRpvhdEs34DgJRGDaN+AehE3unOzZJuy17fJunRWdreohnTBFlRnHaTpH05+wMA7aCGAUhW3pD2NUmrbB+QtCpblu2Lbb9zl5Pts7Lt35+x/9dtP2P7aUnXSrojZ38AoB3UMADJajrdOZuIOKqpu51mrp+QtLpm+S1J59Vpd2ue9weAPKhhAFLGEwcAAAASREgDAABIECENAAAgQYQ0AACABBHSAAAAEkRIAwAASBAhDQAAIEGENAAAgAQR0gAAABJESAMAAEgQIQ0AACBBhDQAAIAEEdIAAAASREgDAABIECENAAAgQYQ0AACABBHSAAAAEkRIAwAASBAhDQAAIEGENAAAgAQR0gAAABJESAMAAEgQIQ0AACBBhDQAAIAEEdIAAAASlCuk2f5ntsdsT9pePku7623vt33Q9saa9QtsP277QPb73Dz9AYB2UMMApCzvSNo+Sf9E0s5GDWwPS7pL0g2Slki6xfaSbPNGSTsiYkTSjmwZAHqFGgYgWblCWkQ8FxH7mzS7UtLBiDgUESckPShpTbZtjaT7stf3Sfpcnv4AQDuoYQBS1ovPpF0i6aWa5fFsnSRdGBGHJSn7fUGjg9heZ3vU9uhJHe9aZwFghtw1jPoFoBNzmjWwvV3SRXU2fSUiHm3hPVxnXbSw3+k7RGyStCnr07Ht8XCzf/2W2fmSft7vTvQJ515N50v6YCc7plDDZtSv/7s9Hn5T1b6WnHv1VPncJekjnezUNKRFxGc7OXCNcUmLapYvlTSRvX7F9sKIOGx7oaQjLR5zf0Q0/JBv2dkerer5c+6VPvcPdbJvajUsIt7PteTcq6bK5y5NnX8n+/ViunO3pBHbl9meK+lmSZuzbZsl3Za9vk1SK/+qBYBeooYB6Iu8X8Fxk+1xSZ+U9GPbW7P1F9veIkkRcUrSeklbJT0n6aGIGMsO8TVJq2wfkLQqWwaAnqCGAUhZ0+nO2UTEI5IeqbN+QtLqmuUtkrbUaXdU0soO3npTB/uUSZXPn3Ovpq6cOzWsLzj3aqryuUsdnr8j2v4MPwAAALqMx0IBAAAkiJAGAACQoIEIaXmfrzfIWn02oO0XbT9je2+nt/qmotl19JRvZtuftr2sH/3slhbO/xrbr2XXeq/tr/ajn0Wzfa/tI7b3Ndg+sNedGkYNm7F9YP8sN1PV+iV1qYZFRPI/kn5NU18E9z8kLW/QZljSC5IulzRX0lOSlvS77wWc+9clbcxeb5T07xu0e1HS+f3ubwHn2/Q6auoD3Y9p6ktGV0j6y373u8fnf42kH/W7r10499+QtEzSvgbbB/a6U8OoYTPaDOyf5QLOvZT1Kzu3wmvYQIykRf7n6w2yqj0bsJXruEbS/TFll6Rzsi8SLYOy/jluKiJ2SvrFLE0G9rpTw6hhM9oM7J/lJsr6Z7gl3ahhAxHSWjTb8/UGWavPNw1J22w/aXtdz3pXvFauY1mvtdT6uX3S9lO2H7N9RW+61ndlvu5Sec+PGladGkb9ml3b1z3X96QVyQk8X69fZjv3Ng5zdURM2L5A0uO2n89S/aBp5ToO7LVuQSvntkfSByPiDdurJf1A0ki3O5aApK87NYwalqlyDaN+za7t655MSIvuPl8vabOdu+2Wng0YU1++qYg4YvsRTQ07D2KBa+U6Duy1bkHTc4uI12teb7H9LdvnR0TZH16c9HWnhtVHDatUDaN+za7t616m6c7Znq83yJo+G9D22bbnTb+WdJ2kuneXDIBWruNmSWuzO2VWSHptejqlBJqev+2LbDt7faWm/j8+2vOe9l6Zr7tEDaOGDT7q1+zav+79vhuixTsmbtJUAj0u6RVJW7P1F0vaMuPOib/S1N0lX+l3vws69/Mk7ZB0IPu9YOa5a+pOmqeyn7FBP/d611HS7ZJuz15b0l3Z9mfU4G65Qf1p4fzXZ9f5KUm7JP2jfve5oPN+QNJhSSez/9+/UJbrTg2jhpXlz3IB517K+pWdW+E1jMdCAQAAJKhM050AAAClQUgDAABIECENAAAgQYQ0AACABBHSAAAAEkRIAwAASBAhDQAAIEH/H+4GSh8yw9JBAAAAAElFTkSuQmCC\n", "text/plain": [ "