Commit a113197b authored by Guillaume HEUSCH's avatar Guillaume HEUSCH

rppg-pad: up to pulse

parent 63607465
Pipeline #25126 passed with stage
in 1 minute and 3 seconds
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Let's start with a real video sequence\n",
"\n",
"![](data/frame1.png)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"import bob.io.base\n",
"import bob.io.image\n",
"from rppg_ltss import detect_and_plot_face\n",
"\n",
"image = bob.io.base.load('data/frame1.png')\n",
"face = detect_and_plot_face(image)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADYCAYAAAD7yhhyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsfXe4bFV5/rt2m5lzzm0UaxKjP0uMiWLvih2p0otSBClSRCwIKqgoAgIiIgKXSBXpKCqo0URjEiuWGGNL1GhU5Aq3njIzu6zfH9/3rbXXmplzztx77pzLcb3Pc5999569V9vf2etd3/qK0lojICAgIODBj2ixGxAQEBAQsDAIH/SAgICAJYLwQQ8ICAhYIggf9ICAgIAlgvBBDwgICFgiCB/0gICAgCWC8EEPCAgIWCLYog+6UmoXpdTPlVL/o5Q6baEaFRAQEBAwPNTmOhYppWIAvwDwCgC/A/BdAAdrrX+ycM0LCAgICJgvki149lkA/kdr/SsAUErdBGAvAAM/6K0s1cvHGkgiWhgopeiIqPZ/7fzWA+3/Z/Aio7eEwZOXf6/yrtjmqPpNXKwafI//vOnn8PBb31vd5pRaK50r0EqbS3SQmqX8is6kUzryh6JWamXLByAEgg6Vc60suwCAoigBAO28cs61Sk07SpRueVDIixJFVW3JIGwWglwjyPVWlGt6prhfa73jXL3dkg/6IwH8X+38dwCe7d+klDoGwDEAsKyV4ZAX/j22H2sBABqNBjdiDBH/MWQRdSRJqGmJdgdM8ciWJXU8RjawgZEnKcq8BKd9XI77B5co9w8qjmPzu9wjbYZOnHP/j1YpZX5L0sh91oce/IdcRW77zcfCa49Sqrcc1dt3wI5rVRX2/8oVyG6Re+V3nfO4apk2aG/Q86gNwL6voqKyiqJApbvm/wCwfv1vAABr164FAPz0/7p8fT0AYCZ+hOnf+nLSebaAwq/+tL5vH7cGglwHuQZGI9cA8JM/3P+bvh31sCUf9H7TZg9V0FqvBrAaAB66YlxXVWVnSDNB2qKqigZdBgqVV2RlygUARKqnylrdg5tnnpeXGblCJBCBrwuV+WORPwZvGH2BdO7leraEyljW5wr87ILvnfK5jLdScc+YyDtIEvtHT/CZqKq9U/kDU075pgu1lyJjIpAPoT3CrV9okrZ99xnnqBDk2i0nyPW2Iddbsin6OwB/WTv/CwB/2ILyAgICAgK2AFvC0L8L4HFKqUcD+D2AgwAcMtsDSgFJCmgzIynzg7nmzbBRRfolM8tqZjg8IZaz6A/trCpF2/NBS0PDUiJhK3LsZTJRJMsxf6nqzfR1pmao2MBGD+yPMstwGSteTiOuX4WCgoo9piZLR1O+MBg+q+x4yNJUnqmqjtOOSDecsmLdsm1Uk1wvPZtXoh+UpX9uyvbfQRIvBwBkKdXfTH8OABhvTgEAulP0e1EUyCKRD1qaalUsGlsPch3kemvK9TDY7A+61rpQSp0I4EsAYgBXaa3/a3PLCwgICAjYMmwJQ4fW+m4Adw/3UO9JXf/kz+N2t1c7v9v96dnMLntKq/2vf7k+++i3MeNf0z2bRn3KkI0Z3dvnOpS/nV4v19OrRmY3XvfcH0Uu2/KZTGX0inFtHIQZcb+YCcaRqxOU38uSGU+kDPP034/RW3o6R2qPtI03oVj3KMcsI+aWJJHT3yiKbJ+NfFSzssCtjiDXTr98BLneArkeAlv0QR8WCjFSPYFMdsSVCEFp+lEofyk145wDLSkMAJDofHB93s50XaY0aKe5kk0OXvaJ2VCmaPkVRQnXn5p2yBJXVr6ad7jtH08s/+EbrPmaLr3NFPNMf0F1b6KddaNaqKjNkWzCyPJMRaa9VpibTpNkT67iTlTRNKpKLAOkYllquxYXpVglsEDmSds2lcfNbCFJhbxrn6oWnxYoa5YBABBH1Obly1YCALZrUZt1k8pct26azjsdRClbJGh+Nq42b0NuARDkOsj1VpXrIRBc/wMCAgKWCEbK0BcC2ltuzrYw7SFstZuV2RySjRdvA8gsB1VvPV6l/mrS+C/UltXRAPrYY/vr/uj0Q3ttMUezTJPz3qWv7jGD8/trl5++iZssUU2ZkXdf3ZRMRc69sgS2zLSnl3b56p23WmMAgGaTGEyjQYyqKPqxV7u0fjAiyHWQ68FyPX8Ehh4QEBCwRDBahq6IOfgzY6ya9p6IZykzS3p6O2+ToG6WFOnZ9U31+TmJMr7WuzkEAJpnZKOLFFMqBcRl4tzr16p07F2oUJoZ3duI8RmGx6AAO0fnvtkYXy+47NR4vSuUKJx7UyU8gcv02qiRGkcSj8igqkqnPp+tRLomRto1sxLmKX4npZYWRsat2bAqedX8Lsaa5OnczOihsYwYTdHehGkQy1Ex1x1vWDx2EuQ6yPVWlOthEBh6QEBAwBLBaBl6Hx2Y+aHPf/ue+xecaVd5Pw02e/B1e8rTMVo9nlevrlsWsS7RkA9mNr4pk1OVmcrpTBxNTHsG68zm0qZJu6pKGwYhQ5Dnef+b5bQWT6P3HXn39pjcDX5/vllZ/Z3470finIhr9vg4sZVWi2OkNIklJdMxEmFirCvWcbJoVi5BrmvlBLleeLkeAoGhBwQEBCwRjNYOPQLSTPXoslDTg6meI+vXRB/lTZVOsFQ/ktxsjYl8XZs4IlAZhbFnrW2xMxrlgJJFn+bpD1WtkWJDbPpjosIJo7HV9e7Ks7u46YPLVnRVi5wXlVye6AVdJlNGru4RVVHTYXbdn7yxqrSvt6z31wvoxHpZw9NqOlvl6UMjZieR2GwnpFtsTZBeeNkycr9utycxNVXTTwMo02hW5ro1EeQ6yPXWlOthEBh6QEBAwBLBItih21nPBlOav62lf2dc0xlG3vykZnEF93/qDXjk6sjE80wphbjsHx/alNXvOl+KexS92rlBvNr66eQqeOxjgN2wUgqo3F1531Ki9JMx1JhN5FkOVN60v8s5lwAAvnjGSXyDHjgWFXsQRrFbXxzHKLlgo5sVNaxYMKTE3JpNYi0TE6R7nJpuoVmK7W5pylskgs4Icu3d7dzwYJDrnq/hNiLXwyAw9ICAgIAlgtHq0BWQpBqxFxYzqmUdMSEzzczo2ev6gXyGmpPsvYOSBMhO9Gxkz+jpBhClfva6wozK0i3Z9yKLKzvj+8nUlPL1g56NMWzsDT/eR7v0Ghu7usd6IgDlZcspmv3HeIqbE+u45l2XOvdIv+WdN1ISOV1ViDh8aSw64i7/xu1IMjpvtCjmxfgEycnYZIGxgv6f58TYuqpp9NKjRpDrpSHXFevYDdveRuR6GCyq679dTuheM6q+Z73QevCyyGxq+NHbVC2Wnb+8E7OoWZa1pbdE9B0itMnFaJfeplwTGMh9tjJRhay5l/U8lmWyH9nNrQ/1dtWW0gCQ5+75vpdeBgC49Y1v4Ov2WbOc5E2jjrd8FrQ7nMKr0ogl+4r34UpM1CQ6vPx95wMAvvyeU2E6aDaUeKmqXZO3NKU/xFaTzLyajaY1jzNmbPHsX6sRIsj14sl1Vfkf9PnLtU5sTHNgOLkWZyIVxQsv10MgqFwCAgIClghGHD5XIYliqIqT6mkxqSpqJkLeMihyZ9w4knCZPCPnpWFEfu7ALp8ffe2tAIArDz2Af1dmqVRWLiuRWMgNLuPw624AAFx3+OvMPTPFtNOGniN4ecZdSqK4Zp7mOSSY6ulc5uOq0mZpashA5S4JK895o4rb5lx769puucypV7CpzZnIy3XWsaKge0+88VrMhtdeeHnPtYuOfhMAu+HTjGksxLlC0FHj0JrGUV551OR3ws4iYsqX8rPNmB0yknGAw4qWvBxPsGzxVC6LJNeRx4KHkWvT9toqYCnIdUe3nd/7ybWM6ybvnRQmZC39XlRdG7ucMyUNkmudaPPsQsv1MAgMPSAgIGCJYKQMXUNDaw3NjMKomgBzTcWii6r9CDtr5xw0XmbvZpIYJ4mCN34qngl9t+AOn8dJAhRuICUJpi98Ki9dhtMtrLmUbMxEno7R6Ce9IE260nZzyGchXhZzw3CiCCp2gwX1ZB435fOsLo4RqoTWrulUuyMJFdz6Nk6SU0NXT2NmhphFd5pyHb77Za8GAEx52dA//JXPAwDe/IrduV4g5U2h9L77AABZRmy1lRDDkYzngsnJGbQy1kcKyyk5uD+7O0cSaIrb3DSu0k00MnqX3Q4xMmJVi6NEXyy5NrpeReO1OXJdx1KQ63ZF8iDj2U+u5R34cl14Om5HrlOS30Fy3R2jZ8bHx7eCXM8fgaEHBAQELBGMXIceQyESF2OemKMsMSeyW2yOzBrefA2xwosP38/5fUORG5OsmS7NbnJeip8zYx1fT6MIReHOxkkkpkV0b+Yx2amarjGR3X5pfySmW+wMwIzCsCSUhgU0Im9H37MKkDlWVaonD2PEZlVS7iGs37/6sNcDAPKa0wX7JeD4G2/GbHjzrbf2XDvk73Zy+qMjm6asjj+uIXZf6MowiSile278zj/PWu9BF5zbc+2z734z1ctWDZwJDAmzoqxB9TWaTYwl1NdCREkNSrWw9bElcl0Z65PI+X0+ci0MMI6Hl2sb4tf2Y1uR645nXthPrgWbpjn0LDPyyYLCzXa7VOZMtR4zM8TiO9MNpz++XFccutbo0PvItejUG2xOKOcrWLe+alWEhz2E/j82RuUlPBibK9fDIDD0gICAgCWCkduhR0oZ11dhEQW6Pema/KzlMNcJYlM61Z5BwUb47/vMP85a9xm3frbn2tt3eTkA4PwvfmXWZ0/+5I09187fdy+nH4b1mOQFdDjltttxyQG0sogS13JAdI+W0Yi+0o6BQPSGxraXIQyk4MTDGkAnp5E8Z/ddAADrNijn3kvu+RoA4JC/+TsAwKbu/baejZu4Htbdso2Cb6kiusmq5hQjZrP7PPm59Awn4hWmc+MPibmf8tL98Ki/XA4A2H77HQAAE5qZmrAS7+VL/Y1GhkaD2WldJb2IduibK9eizxU1cT+5nu66DFazE0/MbK5u+93D0BN3/FPD0Ok+kb8kjhFxI5ws9Bgs1zFgUt4tpFzPFGypInsIfeRaWLYv1xs7xNBln8GR62Kc6+kv1yK7Mlb95LrD985wm+TemYTe/fTMNJKY5ZoHaoLzf2+2XA+BwNADAgIClghGytAjpZBFWU9uK61qu/9yTVyWvTlqusN6RZ6R/zS5yczwh7/wWW4ZOZVx3be+CwA45OnPBED6NhM45/d/AgAcuNPTANgZVxjG9fd8BwBw1HOeD4B0ZpliG9T1xAJS9jwTRtNIxcvLtr3NrsNityvHVLtBeBJmTrTT7w5UmxVrvvXOxkli1DPVOnq2qjDFrE50iff9kV51u+3a6d57L43jpmKmxlTIGkCYWZPZh5zv9mhi9anYytbDxIqylptYRhP0rMcdHli3DnGDmSdbEBSJ6CXpPC1nnGdERLIsxViLnp3k9F1Jt1o0gr4Qct2dRa47nuJY5NrYpXfoPkeuGUaePbm2Ona67sg1y8Fccp2lMWLlJqNeCLmeaZP8SV+GkesZHgth8MPIdao6fB+PwxByPSUWSGWx4HI9DAJDDwgICFgiWJRYLiYekTkqM+ufcs1tsz77zps+3XNt3yc/AYDVUaXi3eUFE+p27OwuM7p40MkuuehBY8+i44H71wIAskaGLCKlmLCPZkbPZByboWj0ekd2u+yFJpYEPC1HsRcnpaZf0xwo/w2fuqWnz3Wc9oW7eq4d8fePAQB02sQ62m1iFF2P3Rcle69FUc2+2NXDRszCYraYaPGOvmV/NQZp4nCIrph1t6kb3KjSOTZNkr5z3QYqL22tBAAkKT2zxzmrAQCff+ex9Lux/W2iwWwnFhap1GKq0AHMLtfSTmGdInfCyCfX0/kkr7bun5qspYLrL9diNSHsfhi5lneXMUvdLLnWGnEs+vXh5dpYtfDfpTD0dpvOO8zCpwveI8pzbOoQcx0k1zmvZqz37GjkWnF9VTW3XPuev3PJ9TAY8QddIY5TKLP84ghwKoUEPb749eSeL4K4YZpMes66nTY8j3gaqVVEcO+ffABYP8nl0CC0WlR+plxnlmKaXnrRzZFPs+OBCA0b8pvlmTeQa9fSH1qj0UDWoE0PeZmdgl7MTT8evCn7rjt7zQMB4MP77Q/ALruSWoQ+zYGGProfbb5O8ljIcvKML30NAPCWZ72Mfp9i86xOB5MbH+B7eVNOsXt1QuW/dIdH0PX2pOmXjHkja5lrALBqbDsA1kRrbNwVzDSZqX186Bn5cGm4pmHHvegFAIBxXaLMeHlZUFvHcxr7pHAnYllIpuIinWlkEanKInYF10U9ENaoMbdcy/jIu5uapr6vX09ytZG0CuZjtp4/7FROf7k2m68zw8u1nNsNucFynfKHKC/svQBtLspvGqJioTZrL1NRP7mWMZnquHI93WlxH1idMUXnnU6EyVzu7S/XMetEcv4+jEqukdg+zSXX4zzZZhx1cS65HgZB5RIQEBCwRDAnQ1dK/SWA6wA8DLSbsVprfbFSajsANwP4awD/C+AArfW6ucrTsOErZaNBR3ZWlCWTLMMmJ91NvMmpKee82Rw3YSmNm7Xk7uNl0QF/SxuemG5zvQoTnHm7ywyj4o2amN1zJbbRi1f9FQAgYfOsNElQ8YwvIaCZLODgJ72E2sTuwrJUveK7d+GtL9oDADAxwa7SshnFpmmJuGgrzkcIBV25Jlp57jp2COR3YSJFUfQs61q8lDNspcFLRmZUYyvGTXnLlxFzEaY2wcxN2ijPioldHMcoCnZHj1tOPWXlqhoUbxApKHTZNivi+NcFh+OVo0C6K44uSRIjkUwxfE9VFFhEir7Zcr1xA6kQJqfcts9Hrk0qHya/w8h1yht/8TzkWvGKo8PBqzQ7+qhSI+eVgOTLrNi5SrHKZXPk2qqllqZcV/L65ivXQ2A+DL0A8Fat9RMBPAfACUqpvwVwGoB/0lo/DsA/8XlAQEBAwCJhToautb4XwL38/01KqZ8CeCSAvQDszLddC+BrAN4xZ3nQqFj/JCE+q6hjHBSmWU+3gXWMa9d4ZjuRuzGTVGPW/Tj2zKrKcekDAGBsYsIUo9RGAEC3y7Nz4urPdNfeCwAqZjOiShm3bVXyRmfOYUeZ0hRCcQo7vB0ODCRsQPR24LHIeSyymBlPFBn9mWQvqdgl2g+NOsWrlqJm3pZlbqYTxbrFsTFePXC4z8YOtGEzviy2oUEb7rMpn+Y5scqWXs71URlZAygK0SVy15lYCPuoWCfYmaJ30u12kbKO9rwvfQmzYY9zLuu5dvUxLwUArGnQmHTu3WgcRhYDmyvX3S4xakR0/cEk1ypv2vYuoFx3OrTh+WCX64yd6vQU763wJmuxnB3QeHM74wxHWQqsWkbhcutyPQyG0qErpf4awFMBfBvAQ/ljLx/9hwx45hil1D1KqXsmeWc6IODBjiDXAdsi5m3lopSaAHA7gDdrrTcOSo/lQ2u9GsBqAPirHVZpoNed9ZRrPoMLD6NQrB1JuDBDs2ZR0Mx+0BPZ8ceE+KRnW40ECZtejTE5EF1mQ4v+jPWEkbgnAxUHQxI93arteKZldpAomq2nOfhPHDfN/SXveJscgX4AHdnRV3X9F+veePYvTPIA2VF3kwdoXXOjZv2nkqBCrFM9+Wmks4/Y9C1ObBjS2DAlcQJxEyhkiegcG3y9RMR61YL1gsK+xHxMdLjdwk1I0Cm6yFgvqWM3KFSkXOcNlHyOEiXrTs/gMLzZMnpm2TJikW++9pMAgNtPPQEAMNbYaOpN/vRH7hf1M0oKDGnhtUWYj1wDlTEfHCTXFeujRZf94JJrK98LKdeyN2DMCh+kci1Kci2mlGxaWbESXfor70xFMZLY3X+LkoXXoUMplYI+5jdore/gy/cppR7Ovz8cwJqhag4ICAgIWFDMx8pFAfgEgJ9qrT9c++mzAA4HcC4f75y7OgWtYhN4RteCTFUlMZduO+WjzGbiJkuza6o42BAzmWamIZP2eJPDU3KgIEnqLcwtqzl8FGg5LStllu4yS4loZpyZ4dm2tLvpM6WwKwnmT+fCiszON1d8/NN2QgS2n81bTvuRyO6/sBCemVVcs+2V0KNcrrFRleu8r5BytvEEAO+6S1vGOIi+TTDMfWCLiW6uUE5y+9mO1rhtV9QmsQ54911kb3/Gq2mFoJNpNJkpJRIeVomVg+x5sJMKUzUKfCQB/+mZcWZmLU8VHjEr6oLriBMsmyCmOT7GLtjJH9Djez8y9JdrjRwVO24tRbmOUEE44VKQ646XYHpryXUcSRC3+cr1/DEflcvzARwK4D+VUj/ka+8EfchvUUodBeC3APYfquaAgICAgAXFfKxc/g391IOElw1fpapZI9QYlZB1k+6KDtZO1p0RRdc0NpbVdthdt+bYhL90rQSqqjLefHEi9rnMKJj9dAvZ6GKb2UqS71ZIeIYVtzhlQsS69rWCSMeIvVCkxp7WpBFzj2+749O4aJ+9nXZLeqrKlC8eiWIhYS0atKQ84w6VhetuXZSeHjSqsPr79wAAjnzW3zr1xhCXaFdc1k9RvemYwnSXVlhZ6tpTjzdpbJqRuEqzd5xumGS5ckhTscd1xU2sH4xzv4rQYFa1nK0CxsdbiOLhLAIWFr1yrVQU5No/QkGxnI5KrgtOfs1iPFCuUw51ICusrSXX0sa4cgObDZLrYRA8RQMCAgKWCEYay+UhD6zDyVf1j2ny9k/1TzBx3DOfCqAWECiSwDk0F401bPyIOOaZL5YQpRQkx8zePH9VqKw9KXtzzXTcJALnf4XC5p74gsfRMznZtZZlaViJYUbK7sID1ovN6hpjG7SosEwFAFKeUhs8e2exDWLVbLlMpsn6z5SZy8wMMdKymuL+rqA+QBvrCRMAiVW2xhuvGEMdeWx30zducuf5yS4nPohclvLHnMoY6zZMuRnrjIVNdlkv22I2Mp5xWWmC1NN5NxrCyIhF3vSWY6hMbk7BZUQ6wTgrmJePSdqwxQufCwUgihExK6yMR2UE+RNLWPYkZGvc9JKGjEiuZ9quM/eo5TqOY6ikdMrZ2nItY1Lk7grDl+tonNrYam1duY7Z6zSap1wPg5F+0O/bfhUu2vOVaGpZ7tHxxOs+jfMPpgzzG6eoSWJW1eqIeRU3NXaFSkc5NOczLErtlCtmSSagDn+PLvjK9/DG5z+Zy6FrM235I3Q/Wpsm2ZSplKWqQhx1nDaoODe/AfZF2ahuupaFXVyXeTMslYBHHOGuFpVwYsINwpSOkbmaLP8mZ0jgO9x26S8AVJH7eTMmUmYsCue8XdgQC+vX0Zh/5Y+/wmy49Ztf7rm2904UK/3TP/wxAOD1z3kmt5H/mDj0Qas1hiyWAFH8MeN+RYkdA6C21OdxLSsNzWM/MUFL04nxMTNOo4bWQFXWg4OxbFYwJnjSx4yX9Yn21FgLINeV1kaOB8l1kbuZdkYt11EcI20kTjlbW64lUNiGTRJAi8435hLRkNrRWU9yP8HqjjSeQos3XRs8NpIntMv1jkKuh0FQuQQEBAQsEYyUoSsoJCo1M6LWtfmkYicDCbfJ0XVU4eZp1ImE5eQNhngDKp7VxI1aAuhsnJLgP5w1Xdvlyy//tJ6qlUBAncypR7B2E7WrKTkY4xjN1DP2j73NImEwYl6W2OzoJZsv6ogY2qmf/Q4G4dArbuh7/bMnvhIAsInD5W7aQHkTy02WAUQtN3TwVT/9KQDgsCc+FgBw3c+/N7DeQcz8b1p/AwD42czPAACPGXsGAKCVbcAyZjXNB9yNyfXr2GyPAxNFsmmlgUwClcXsIs3vtpTAZV7s7lJkADGShPrXaBK7G88iLBJBHyjXVV6NXK7XTrnZfny5zrw8mltLrmOJ8d3k99TgEABZDCa5JvRuwgS8ybq1YeS6itgpidUcpWYWLurBqo2c1SPtaTIJlNDbU9xGE8a3oLI3cqC0VpZj2TIOG4A2l8FqInZ6GoVcD4PA0AMCAgKWCEaesUhrySvoRjsVXZiwHGEqkm3bZA2XAEJKdGapCb8pSQM6HNR/wyY32E9HW8eBPz1AjOIHD6ydtb2f++Uve64d+nhiqtf/gpjqYU96LLdFXHlZN1dJ0P3EmqsxCxJd40f2exEAYOVKDiY0Rkx3/4/egK+852gA1vEhLWhDqzXF+wtNojpKS5jONpcdIWq4+SAF4j5+2N/SZrOwhCjr4hPf+xEA4FUPIz34A2tpbGbYPVxVbi7EX02TmeOLdnwyxtg/PdOuzq/LYQnA+sIGM8RGmqHioEmKN4LEfMtsrMVWV0tgRqcSm8WHAy4tHx83m3qLgX5yrSs9UrkuyxIbpt3kEToXpyQOq9CgZ1qKN+Z4M0/FGmUhm4TMGJUX1nYIuZZ3kcpKgLMdNZsZWk0320/Km4VSzzByXZocqTKgbg7VVqbMflyHdeZtjr0TaTH35NUDv6OM39H2K1o9ci2rom7Oy4oRyPUwCAw9ICAgYIlgtAxdKURpiqotLIIuf+TgfRCxflBchxMxnUpdPV5sdtNpLsq7yzE5SWz7yn8jvfAej3sMAOBz/z3YSmMQM3/ixJMAAD+d/C8AwL5PJja+rFnL0bjRTbqhIMxG9IkE0WwW3QIFm7Q1U2KsEnynJTpGtn5oZlYfmsX0mwQESiRxhrCdhNmwhDflFHVpmkKNuQGOBKXkfuRxHuNd9OU72LCqK1e6+S6RE6MRUy0fX//Nj/peB4C7fvvzvtff8oLdUEj5zNQSZl9ZxD7vlZh7sR6YHUGiKkK7kFRg1PdHrNoeaTzyBSdhgFwnUYKi2nK53sDyNjND72HTJgm4RWWsKyldXVEUqMoVXB4zWZDeOGK2O9G6DwAwtpwYeqspSSQ60BslhR7LDovOQsp1M9MLKtel9sLp9pHrTmec+7WBx7G/XDdYH7+8SWP2yO0yjI2JR9a482zJ1jQQK5tOYY4LLdfDIDD0gICAgCWC0VMabdlDVcvCbXbQJXRnIvTADdBfQfR69Oz6tQXWcxovQZd1i6/4y0cBAHJOI/VAh5jMf67dhMc2HwkARm+YJTQrNxOX0T7kIQ8FAKSa6piamkYRu/WJPamx7TWpwkRHp0yy5+WOAAAgAElEQVRfBdLfmGftRI6x3cWPFVtIaMlOzsxGMq2bUJ6iV6R2XLfmAWDNA+iH63/ys77XT33VC83/x8eIWaxcQSwn5QVJnlObX57R2H3lPirryGc+AxkH5I/YeuPj3/oXAMAbn/Yc6oNYWaTstNGIDJsQ93STvoyvC9uKjd0z22jrCq98zxU9ffh4356NCFtRrqemyCJLxlBsvCWcLceuonRlLAtFV+RarEvovTz0oQ8DAGy/PTG/fnItbSjZ/n0h5TqO0s2Wa10V1jmIq5Uk0ZLSUfT+Beu4yzxFyiEN5pLrLKNVqqwuWkkM9tpHxBYrSUv2nCSkgivXzS2Ua5uWjp6ZGAs69ICAgIA/S4zcDj1C1BPIJ44VZKIXXaPs/8qMI4kCqtxNIDu1cTkm17u63Ypn58ZyYhzjvNMetZn9rt2EJrvfiv4s4Sl/hQTDIVKE7ViHdfat/zKwX5/82W/7Xj/q759CfYkagBY9u7gQE1toxMQK0micjzWGXnGQfi0zOOvaeFd8rEEWMcvHVwEAmhnZ7R6+Ywvx2PZcH43gP/ya2njo42nVcv0vfuO09UNf+lfz/6u+/x8D+0r4k3N21XfvGXjnw3ak/plUZRz+NEkSk5ghTZh58svnU2MNIGFpTSIErfGFM4+j35g95usfwIYbZk9lt7WwteU6L+8FYK1CBsn19HSJyY2zy/XKFcQ0Ra6raQ5IVcSYiiedNkjbDevmtHLG3n4z5DpGutlyXRRdROUyHki25Za9CEn7x4knOhwTYHJt18heGlH/tltBz45xfIKK25Ml4t1K15fFkbHoSRKxChKPUC9l4ALJtVjtyGpoPA3BuQICAgL+LDFyHbqKlGFVAq1LxBIW0+ij6beY9WwVJ5Qtmcmc/aXvDqzjy/cyg7zXvX7A0/8f/+9+7LiKGESbd6QbrANescwdkpxThr3phS8AQLa/+RQxhqv+69cAgNc94VHcZravjUkfH1uXOsNqFOsHtXZjUpggR8rWr0xMEOJ1eSlhWSXcJx1bY8QOxtlmtd1um4TFYnUgiFkvefgTnwjAJveNIuDKH1H8laOf8v/4Gtv65u6YxBGN3RU/onfw9ue/wLDSZkZWFh/42t0AgImJFrdVbI4lSFNiM5rZ6KF8lAQRzKAkXgjrmSttU7vJGCVJZljNYmCh5LrDwbSKzhQiti+XwFZynGiybTfbO2ec8WKmkaEVzyHXvGckct1mC5pupzQvwFi3MPlVEEbp+TZsllzHCyvXqVifuYG2YpaVvFOi7EoyakKD2XbS8uWarjc5fV0a29VWxvRawhJPTPDf+gjkehiM9IOuQYMqpluyXCnKNpTEcRZzIDnnXIXgzRC2+MFbnkfqjN//njJtA8Cn7yX39kMe+3gAwPiqtlPPstqCZMcxKveW33lmdV6CkA99/u6efhz+pEc75ylv/IjAJ5B40iLcsRF8s3RjZ51Dr765p3zBs99zbt/rXz+dcmyKE0qrRfWPc/b3OElQVCv6PptJzGuJNc1/iKqW6aaBCaf9ccP9UKWJuwxs1ZaFslw959V7AgCWtWQzUD74HFkvVUDJjhcsD2nG93L8b4mTVpqlvziGVRBv95KVGEmU2b+gEWMh5VqOaTWDMYmymHLwKnZUGc9YrnlpH7MaoIEULS43Z7f6rEFj1mqRmiZqc/C7+8lJrTsj4QU0TCQOcW5h2dgcuZb3ImEJqlI2OhOz6WpUVKwmyQsJLDZ/uZYIjXaS6ZVrE298nnItE3OW2s+jyLVU02CRH4lcD4GgcgkICAhYIhjtpqgiZb+YJdXNu+T/ksmlYpMpSQKjOd+guOuK6dMO28dot3nKYxXLyuW0ZPr4d38xsC23/KL/xt+xz6CAU1fcQxt9x+xEYXbB+QfzPEfimS1asy6XIZr40YhqS0/J+E33XPm6fQDYZeV247RZtcuHP4J/P/0tVLWEBmX2IWaZRc6BgngTabvtqY3T09OYZiYm7uFHP47iukfGaYSXiKY9tu2p6AUko4zkbWTW0ErdAEnNRmqz6SSS5YbHhCmnkuwsTE8SFZudQhkbyegS8bOamZUxkxOmFcV2WSsxmFRp2jdqbG25brbcDPORF0+8YAecZtyESiX+uGSPp3rknUa8Uuiy+zu0uKfnPRmplNp8uc5ZlSRu9mMc6CtCjDJyQ9xqxe70vNIeRq7bEhxL+tdHriMvGNZccm3MQ5N4m5DrYRAYekBAQMASwUgZ+o5/uh/HX74aq488AoB1kVaxssHu5WaeaipINnTO8ydG+bwJsWyHaXS7XNBP6LCS9iTx9hdRbkxhAkVBM/LF3/oBjt+JdPDCPoTR+OHkMzB74IzgWUIZ3gHgmKdQWAAzS4sOlcuyOSAz078GZxEXJiabOcJoio6dkatC3LLFiYGZTS5ux8Tqm6ybW76KXJvTZoSYoqgi4lCg+UbWxYnjgmTS4XZopXDKU57JfeR7RGEoYU65P5KF5rxdd6HzSBndYpq5Or+Ug0DFEoAInLVFK8C4u3NgKh434dnK8A12Aa9tuAkjko2mKh6OySwktKbkBEXhBudaKLmO2VxRAkPlbGZnVgIVjXFZlUhTCfLFDi8e6y4165x5i7CfXJt+qYWX6xiVeVcLKdeX/oxCdZz0d/R3XZdrYcKD5Prcr5JJ8jm7vJJ/FzPDB59cB4YeEBAQsEQwUoa+Zoftccnee6FpMnMTlFaIxVlByc46mzZFbhbzrCEsiPV4cWacGd6/KzNMZi5KkglKONPM6n5bLTfBgMTD9zObx56TQKQiFIWbqsvmg2RTppoJEwCkqZ03I7EyYbO1Muf+cLVFZUP8SnZwST0m2bdkl1zGKuKg+K2I1hdJkiHJ2LmKuyxmV6J7zLmfElhMx5lNmdUQMzUeY3HMSN0s8Els7zN9z1x9q6Qgk98TYZCx1U8anay4fItZm7jF8zgUJjyuRmx0/nzvIlm4EDSgSygvl+RCybXoUZXoZgfIdbfThRY5HSDXvv521HKd6GSryPVpT3saAOCS738fAPC2nXbi5gyW6/czM3//rpT+srEE5Dow9ICAgIAlghG7/pPBfwwxsDdTsmUmWkJK8uzN9p0ya6YJzfSiP6zQrO2Wc0JpsRiIXOsAXSPfraYNiFMvDxyg/9Tnki4uk6BKbEmgtUYLrj6tilzdIgyDYUbTyG0KsIT0g5ZB8TOlhAOl/t715hOx24c/BgD40snH8DjJjE/llvz6JFVZmootLJA1KXbBxATpNseYXbXbVM8Up/kSPWY3jgyTKVlHro2jlLBLbjPnDItTCRYVm538JHGTJIgpr7C5mO2UoyhC4jmqKNZHWmWjm+G9Uta1XrLai5PGYjL0bUWu0ygxAbUGyXVkgoKxBcaI5brUOUo25ShlRbOAcn3Wc54OALjgWxRK+7RnP6dHrt//VQpzcdYur6IxWEJyHRh6QEBAwBLBomQE8ENuks5O97nWezSznsxcWkH0TT5zkTK1YT307Fmveg6iQnv3uvWYEiSFmNF/kX2q0x8J8ymppcyszsdUIa7p5RxIGqxSlJ1ytLeI3tXvuzCLStKaMcvTACK216liDrG6kthXp8M6W3YfF9vf6ZJc1QEgVw3us7BL7hesdQOdW0/ChPW9YgljdIveUdzK4zgy+mWxC5bdfz/IlXlHRvdYoTKiK+8lNvrhxcLWluvTP/0NAMAH96CQxL5cK9VnzLiMc75OuuXTX/B0fnbx5FosOLamXF/4iucBAM798jdw1kueCwA486vfBAC85+UvpXoWUK73vuhKAMDn33Hs0HK9z9kfBQDc+o4T+sr1MAgMPSAgIGCJYFEZumEpUKgFFHXu7UlyDJ+1RdbeU7l6QwXRUUl98kgEeDvOElRIrAEMU++jf48T11pdK+08k6QScpOOWapt8lxvxrX6NWYLpRwt20xMYuHU6V8aNbl/rBeNuqbNSpIAc9viMW6T6D85c1dRMLPJZ9DtEhNqd5Y5bRTmkkB264W528D9JhmBeOUZD0Lpp7AymKPY+wrq3qrSj3p/NetcdVkiryTkKhW47/mfwPuxuNjacn3ePpRQ/J13fB0AcO6ez+ff5ZFeuX7PPxKrP+OlxOrTavHlOhKLkC2Q64NW34r54G0AAGbm5uo/AwDOefVe1J8FkOvPn/oGAMDu512BL73rBKcNg+R6t/MoScvt76D7+8l1FLle2XNh3h90RX7A9wD4vdZ6d6XUowHcBGA7AN8HcKjWun/SSekIOOIZ90/MllRkOynLE1kKmtULlyHLI+MGXVtmm5yHXK4WU63I/WhXZQXlL394M0Jchq0FGv8e26VzkvoOGO4fMluiIZZ4yElks75XHDSIm53yKzjypusxCDtfeHnf619+01u47dTPMrZ/rGUpf1DSPzqNOBBSxstrcbZIGgm6XXaOkM2p0p0YFW/mGMMq+fiUgC4kHyVfEpM0z93avgtj4WbfoPGHl81FueyqzpQGKnb4Umq4JenWwKjl+sP7vgQAcNrtXwUAXLA3ndfl+p1300f/A696MQC7uT9quT7qxsFyvSW4/ZgDqY0SB12cB3Xl3HfwNZ/BlQfT+Bx9I43XhbuTmeLbvnDnvOu79qiDAcwt13e//VjsevalAIAvncYfdk+udzufPuR3voUmAfnb6SfXQkbni2FULicD+Gnt/DwAF2mtHwdgHYCjhqo5ICAgIGBBMS+GrpT6CwC7ATgbwFsUTdkvBXAI33ItgPcCuGyusjSUdX+WTZ2qxmQ0N0nODUsQ0yZ3kwmI7BJXMnAL24FkYReGx2oVXRr3XBvC02VIfmhtVVubRv6oebOoLENtZpLYMJdMnAzE3IsZ8y2vPxYAsJxdqLMsw0s/dj4A4Nun0sJRZW7FzQ75QUfMnAtmdUWhTW5PUSGVMujMYCRIkzEzSyokwrbYUV1MGmXcjEohctkYSqAySR75Gq8WZENNwoEaCh9r41xj1VuuWsK8aQlgJaoFDZQxO2D1yMPiYDHk+uL9Xw4AeNutXwEAXPial+LUu0mdcN7uOwMAkmp2uT79C/+8Gb2dP24+sleumw0OCyyBrjy57rBcm2Bd4GNRmPDCc8n1gdfeTvW/YU+MsRx/4uDXAADeeuNnAADj+5HKxch1LHItm7W2jYffcOO8+/z5grJp7XYuMfW7TuHzi2il/blTjuVOcJtnkesSrhPmXJgvQ/8IgFNhF2zbA1ivtQR5xu8APLLfg0qpY5RS9yil7plst4dqXEDAtoog1wHbIuZk6Eqp3QGs0Vp/Tym1s1zuc2tfqqS1Xg1gNQD81Y47aECZTB2CssbMhI3Y+k07+Iq3+aKstlGcJuzGkOh1JaC+NEoZnaaSOJXicMHd8AmTqgfQ8adBYV/eRomwMqWtiZ5YhqWsr2yZPIaci9HoJ2v1iYlU4josmGBK3HbFrBwV0GXWLlO2BE0yczBvVomlnYoiJKw7b3HMVZ0mzj2SgABGb2113f4mkdmjkPyU8qy830RD9wykK1ayl2TM8ziZwH4f7b+nMEpsS3KtTNH/zJuAAPC1efXjgt1eZuuMvffhWj6aVYQk8IjjyGwSNjhLlzgdtTjcbGOsV64lAJlsii6kXB94DTPzw/enZ3XeI9c3HLYfAOC4624DAFx10GtkAJyO1+X66kMPoPokHIIkAUklhAhv+icxdr+U5fPDXBzo/M7u0dQNzk7ly3UURdASisP8HQ1niDgflcvzAeyplNoVQBPAchBjX6mUSpil/wV6cv0EBAQEBIwSc37QtdanAzgdAJihv01r/Vql1K0A9gNZuhwOYF5bxlrb2Uemv7pDhgT6t6ZfcM59JlNHBJfJWAbjlYXeQDkwv7H+zCtb1eqNvCBMWqfOs7IUiJihR5U1f0rZIYIJOZpsldTkUASZcbyxdYhur4fJJGLWxbr6btP0t5AxFXbdlf4KM3fHGYUy+tyEKZ/ov8VcsJGKG7frKFFVtfflhVyNzFjwuxCzubwyLs8GHhtRspyxudEAALeffDw67KYtDGn/c+fcvtmq6CfXb7z2Dlx6KLG/YeX66OtvwZXMCn25Pv6mz7l1cxEf3uuVeMud/wgAuGD3V/ZtZ8/Ks/7biORaKTdF20LI9X5Xk477Rh4zs6KeRa5vPZKY+pFXEVO/8YgDuJtbJtefPvb1AIDXrL6a20JN2fNScj66841uKA+78rKmotbatZ8yZDC2xLHoHaAN0v8B6dQ/sQVlBQQEBARsIYZyLNJafw2snNNa/wrAszanUpk8K89Gl38dVLf8x7len79kcjNHYaFyNJNfZGsZMAFGA35QSll7aOMYIx2SU2GwwkAio3cUW15z5OA+chS7XrF3pYJETyjli+ODZEB3GY2KEpTMmDq8YdeZIkZblZVzTE0i3ASVcQeX9yHMSQbUbaNBbMfKvh7t3KuMfbD0r6ol73UZoHGZ9py7ZEiSOELEKQGRL74dusCX60sP3QcnXH8HAOCyQw7s/4wn18d8khxmrnzd/ubacVyG4NL9dwcAHH/b5wEAH917VwBAVFW4aA9KOvK2z30RAHDRHq9ynt0m5JrlV/TvCyHXN7yW2LXop4eR69uPpGdTn99uplzvd9X1XO7h/Nu1fE7Mfd/LVtP50WTpPZtca2/VNBeC639AQEDAEsHiuP57noB1bz9L/uZpX6zrtr6ehYq4Ior+suaFV1na0bfYyJ/rareV2r3H6qWZ/fK9sYQoTWK0WpKiiz3BImHgHFa07HB3vLCnAKan6J6crSYkuH9pbJrFTZjDmzZipBxyNWcdnIThFE+6MhcbcPYYTRvGesIOvfRPvNf4KttQ90uObW2wxcZWfhNbai6zqmxo2f4E3bTDkFgesgpAxGnYjMfeNoB+cn3pYWRtccJ1NwMALn9df6YuzHz16+j+Y6+3ru2rD9mbyhUWfxvZUX987z0AAHEfuf7oHuQNecrnvuCcbwtyHRUi23RcSnK931XX49bDD6XnPcsmCYh2+9HkIbrvlf8AALj5uDea0fDluohyDIPA0AMCAgKWCBaFoVv9Nc+ytZlM+8Rc92fQzgOGkfc/Ko+FK6V6bck9RD311mZricvBs2hu4khwmi2Jk9LgcJxxYlKNVWxGm0tigRmODdFmG18lHpyWeU5Nkedcwso2CaKVa/G447J4NlcActl1jyiYUXNc6mMd5DR72rXpmQiJsQ8WpmICSHnxUiSolrUAsIMpdv02uH9/ayJC/1WS9n8213mcSw2diPcoHW886XisvXl+AZu2FmaT68teS8z8jZ8kpr76tWRhccwNZGFhDX5uNb9bvwl6V8fcRMz88n324md82/Zeub50r90AAG+68y5qx567e60evVyrgvTf4u25FOR6n09cBwC49YjXGq/VXrkWXww63HbUEQCAAy8nC62bjjuuR66rYjiGPtIP+kPufwBv+oere65fwstMoDYWcj5HjGula3d4G0xm6ehnkvH+3w+93/u6o484QNB5zv7I3YLzHbJTRXOMboijBAVvQuY5Cwu/25j/aA6+9pqBbdn1yo/0vX7r0WT+JEJWKDG3qhCzYMjyMWuyySG/8qn1VG+bc4wmqkJD4m+LUMmk4gm+8qLw1Te1lR9MaLZhlnt9HYucSUAxiY9tHMS0eE3bzba804cNjBYmJya3oxb3yvTl4wfTh/0E/rAf+0m6fvkh9m9AHvDlevUB9CGPiuHl+rK96EPeGylh68q1fDvjWqiJhM0hBd0ZksGORPss3Q/vg0Gu7zjqMPpPWfaRawE7DXly/aljj+brRX+5HgJB5RIQEBCwRDBShr5mh+1x8Wt2Q9mZdhvhzHa+s4lXiD/p6dpNcjBmSv6StDbL9hB/94Jv3qWNE4CGAm1cdNrEXGbazD7EHEtc57mImU4OcEClvONuuIip1HWvJZOmLG5LY3HQtZ8CANxxBP1W1BgqAMxMSW5CXm5qaldZVsZNm/1usGoltVVYmDgu5QW9i3Y7tqFbOVa6hCewm3zuUtQca7GwLSkxTzltrLnUmGdMgDQvEJHZpJPfCw4joDUqztCe87VC5wPUOqOBRonSy1jkyhC17fgbbwEAnPApvsq3XHFwT4FYCnItMpRJ0CmtoDHp3Ft0PbkuxEz3wS3X+19FL/nm1+/PT0peVfp9vnI9DAJDDwgICFgiGClDJ91nYc27ZOasbz74etSe4DS+LkvbcrS45UoYWL7F7I7W3Ol9t3NvbvM3Uk1+gEqjyGnYWE2HXPSGzBIqrq+dE3uY3NQ2oWjzdovvpdl/LCX2IdnDW80ZADYAEgDM5BsAwNQrjGJyRvSXNIuvnaJZfXp6GnEiTIazwWhiRa0W1V9q3piq2EGjk5lxlMBgmUla4OsWxQFExsy21Y6b2ZH2npVNJW0CEVm454VQwVKu2/5XnEulKHN+Us+5L7K1MB+5Pv5G2uy87GDaDL3sIGGBxNi1ot8vP2g/KXTJynUFNmXkyoNczy7XwyAw9ICAgIAlgtEydGiUZWlmNXFvj2ssvKrc2UxVnuOIt9tcP7UqVHc2NQ4CtV3nng1o37Sxx+VWCi+Rz3D7O1y+CcVL1zsdZi3s5PCnBzZhcpKYxPQkDfn42DgAYLtlKwHUGAazlzSx47BhitgNW3shy8jaYHKG6lm7di0A4Nf/R3rDje1NJkjSypUrAABjHMZ02XLuZ5vqMWnmStvnbs6OI6zbZKsvk41dmKIfTAuw6yftXZGgQ8YQCb0WTH6Gc/E/sWnZiLXEcWRCrFbMgotqbouorYXZ5PotnyKzxEsPJOcg8fP25frjB+wDADj+Jrr/iv32seUvMbkWqRF2H+R6drkeBoGhBwQEBCwRjNaxSLPrr7jaitt24d3kYG7W1ctg+pdVd+XVvsF7T71V318rXQA5MYmKmYyMouyOl6JXZKuHtes3YM2aNQCAqU3U54mJCQBAd4Zm5x122AEAEGVsNVBpXLjrK+j/zNCTiihFl5WbD2wgpnTvH/8IAPjDGrqvgw4mmtSWzz7wEwDAaQ97LrWRnTTSXNJ6sQNIBVT8YjpdSYRL9bRiCT/qBmfSJiRsfZxc2+gexugMs/+jn6TYLVMCISkoIKJ3cOSlVwEAPnbU6xdJg46Bcv2W62/Dh1ln3qgK76H+cm2Y+i134PJ99x5w74NbrsWqQ3TkmyPXjTFq66ZJbuMSk2upIK+Go+iBoQcEBAQsEYxWh15p5FMllo8vA2B1aIp3dgEg93bpE882Oa68Oagq+6TzYninTsozecaEbvVvdi0W6tYAJXvOybHSEoCI7WeJUKBgBrL6xz/FvttvDwBoz8jOPt2UF0QxioR0j1WT+js2FpsQt11mGxm3sd2mZ+9dQ/fev5bOVUTHZVmKlcu5H5Lukj352puorbm3w15EudVtsr3zTIdXBukYHTmcqHWw67WLNhYYPJy+5Yfs+CtV2THuse0VHS7rOk2iC9EvKpOIW1C0O46r/SgxSK4vOmh/49I/l1wrT66v2GfvPkkw5IJX/yLJNfIKbDY9crleNs4Vb8NybTG8XOdmfyEw9ICAgIA/S4yUoVelxszGDsYVJ5RNeAdcWW+owptiEs/qJarcOBAKpQ02PweTqV9QyrV9VV6qp5J1nhIc386uCpqDJcmxYlZSciKK7ib2gCssC1u1jHbl252GU97n15GO+8DtHg4AaHZY1xh1zWzb5a3uTkl0u90lNlJpYhhfnPlfAMCeqx5NZTSbWL6Cqfka7i8PbFFKElpqo9i767iE5hgXYm+csy1sx3gM8oqK42eoSMTHMpLYBDii+rqlN34myJFNPq3NGLvedpKzWI55LY1b5Bs/FVVPsKpRYUvk+hK2fjn+VrJuuXw/TqH2IJDrVpINlOucg2dNdqkPmyPXcUr6+O04jV1drpct4/5tw3It8Pf45iXXvHpRxXAyPeJNUY2qU6KT03IsYsGvVM2ER6K18WCWSoLTSKcT5zwGTIaV3r/nwRuqsXZfkPwBSLzjngRJlWxgaKS8b1GIDy//3ValuEG7SywAmJighx6Wtpx68Es6NBu8mcObOtNViYgzt3zsW98AAJz0VEkQRWNgTLUYn133a3tyr9vfc77/b/4QAADOft5LAABpopFEkh9RdsPo0O3wZpIIvqg7ZOxi+w5sDGt2xTaxqGksqvpLMqoC/uPwv1ucS7XkMo7gTDD9cPLNN2Hwr1sZCyDXgoJPHwxynSTVQLnuxhnfM1iuxbSwNLZ5rlw3W6SuyRo8SUURWhx+QEwcoxmWN+6v5DaNWY2yuHLNt8wh1yWkrQlQSvlUXiubwTAIKpeAgICAJYKRmy2qQgMRzYgn30yc6oKD9zezZiRLRJ7wyliyirDTg6oVBpqRZAbs5S2zmTwKc2GGwhOuXYrCOcIwHG1m8qwh7ICG8V3f+e7A2j7y45/N0hag0eBAW7wUL6rSCQ4EAHkuzIgzuGRu4J63Pu9p1BelEHM/zvvm9wEAZz37Fc6zWcqbVLHkV8xrQZ44eBGXIZtJ4lgim0fCWnRqXb4l+JOJLS2e0n1WjtYZg5mfZ2JnN57ofPWRFKK0IaltABx+GWVSv/J1e+H+u7/WW8ko4Ml1Bd5UhJ63XAtKte3JtcQt91UJShVIU6qPibpRm3Ul5HG1+XLdbNF9SWZXGzGzbIkX3uAQA75cx9uAXJtz73dfrkU0kig2JtxVLuWH8LkBAQEBf5YYKUNXCkgbCt4EjQrKTFe17R3+kTdbTBCc0vkZVWRcdX2G5wcicuv02ublDKwq5RwloaLWMMGFTLB93kT50AueDcBu9hTsKPHe7/1gYDsEq3/wzZ5rxz39Oc755T8evAIAgAu/8f2BvzUmmAWldBzLmHEwxSgLoBBztZIZsDAJZiXdkkOVdnksxOGkVGCVpkkqImpEyYIjG1LWXE6ba5LKYaBLjHJ/L6oSFdzVyaaZ/0VZdbEY8OW6kEQGiOYt14JKSeLUbUeuKzap9OVaoayFPrJsWhUAACAASURBVKbNStE5R5EYL4iLvNW/2xyefC2S9ytliTEB95s3PpMkgRYHLd493JblWlB572s+cr1xegO36TcYBoGhBwQEBCwRjJShR7FCazxF5eXB0trqk2Jfn6ppilSGxRXOOZDCzks9/GRgW6wKkXfFPUeBypgluUVrDShmUbLjLXq1yuQ7lF15W9/pz30cAGA6d9t08T2/AAC88RlPAQBcds9/mN8u/963Bra/H970TGL0cRyjwXrDc79JzD9pshNHRsym1XAtFvKONkGldO4mHBBdrgQvKkTpxynDdKmsu7RYIzGDsvpE16IF0LVQr/1Ns0R3a/I4xqJUVigqV7fYGGtbK4URw5frSnS2Opq3XF/I6eVOueUOAMDH99sf27pcZyo29eS56ygF3keQdyz7JEopxMr+HwCqyOWVmsex5JVOxCnrqipGaqxOWL++Tcu1wNsnmYdcz3Qpj/CqsTaGQWDoAQEBAUsEI2XoD9uwAafd9bme6++46Zaeax/anxPiiq5PgvzLZCeMA5WZPaMBgYcEzjwpW/aVlOMdfRfpWpmSRdw4b0iaKCX6Q7GJtcObsp1uJqpLz2RAynrjs3cy12JFLO5j3/oOAOCEZz3baZNiq4pLv026c5PdPI6NO/OZL3k5AGBlylYVKdcbN536oyIzQaWEMJXGAsKzlRaSabIiK2jJqK7EBZpZSCK6RwlNWjezqMzzdHDfWFkwC2drC0nuq1CgyOn/FxxCwawe9SiNtHEfFgtRVKHN6m8J9pSoCiI5pejGxTprkFwzHgxynalxE7DLl+uicjtUJ+EiG3KHylxWHZt75bxXro0sbMtyLb3w9gznI9fNBjlsPeRREsPj55gPAkMPCAgIWCIYKUP/44rlOH/nF2CSmcD7br8TAHDugfuae5SxXhA7WuWc98xACrOb5dZRozZ2FvaYi/+Iua7M72ZGNVHBbKhOus6uxtz2D7ziOajapAsTu1Y/LZVc1zU9cOQxVh17SYi9wVBcRpwkhjmYI3vsgcuQNFhivxsnETJJhSVJGFgt2vUyaokVgBwrrVAwvckNk6GH05jedWIiH9U5oaw0uCN+8pLYtXqI+HelgCIX3SatYia22w5R7IaFGB00VKTtePjpygDD9EzwrCUg17HOkOfimenKdaxE1ywWJbPItXcu4mB0zPrBKde2Q+7psHI9DObF0JVSK5VStymlfqaU+qlS6rlKqe2UUl9WSv03H1cNVXNAQEBAwIJivgz9YgBf1Frvp5TKAIwBeCeAf9Jan6uUOg3AaQDeMVshCgXicg3SzJ11lI6Mni7xbFC1koBBoovk2daEB62s7k3Kk7gLXv2Vrsez0M41o3MUe125Fy7jAYCu0e1x3AUOSh9zQtxSwo1KcK6iMHo8iWkhROrk5z2Tx0BL4+sj4/XAZVCym/7mFz7L6z+M25swsApkg1tw6rOS26w48BJiZdWv7G1XGeYk4T3hlKltLFGjA+6KRyLbLGuxfpD2GDairG2yoa3cVhn7VsZtJhlosJVFUZYmBKnYJT+wZiPKwk/OOxr4cm28CLFiaLm+6AAK1nXirbfjkgNpf+DBKNe9AcX6ybV79OXapsCrFfEgkOsDr7gaAPCpY46gtokH8WbK9TCY84OulFoO4EUAjgAArXUXQFcptReAnfm2awF8DXN+0CskagoRR1Mz0Kq2AeSZWcFdqpY9373KvAAbnU7K9aqRJ3Q9FrdTfE3w3XaIVERRZCNZe5tFiGVJyi7T5o/LVlDxRpOex+JIqn7T82kzNNLuBz1m8y5zbgIIKcD8MfAfo0StlHP5TJTiKaHNcjXnZXMpagIWuNKERzR/AfQorJOIqBtMppVcXpiMEX804th8mVTsBpISdUTBTYzZu0NyXZZaI0v548L9mp7soCy9Fz4i+HJt5GGL5NoGb1qqcm0djB58cv36624AANxw+BF0c12uGYPk+piPkxHIFW8+mOqfQ66HwXxULo8B8CcAVyulfqCU+gel1DiAh2qt7wUAPj6k38NKqWOUUvcope6Z7Pq2qgEBD04EuQ7YFjEflUsC4GkATtJaf1spdTFIvTIvaK1XA1gNAI/aYUw3VjXRWbcOAHD+Hn8PAKja65DwzIRUllSyNHSbWEauc0AW2RkwqcRVmuuWCbe+XwEiFeaS2TxiBjOgH5KFW0OhNBshsnlCGySSTUmSpnTZUaGoKrO8UsbbSDYEvbYp29/IrDXpYJfedHPCm1ixHHn5GSllWYHEXC5kiU/X2wkHMZLlrdbGkaMwpnZSHzfR2Na5o1QpZYJKmX4ZBsgsRViQeX81B3ZjweeyoYhdviMuY5q/m0nSwAmrrwIAfPLMEwEAy1YkSBuj2xSdTa4FBdLNlmsAKHmjL8j1tifXqw88CABwzLXXAACuef3rnHZQPW7Yg8OuIFZ/+Rvp2TjvlevWdtTn5dlKACTXw2A+DP13AH6ntf42n98G+sDfp5R6OADwcc1QNQcEBAQELCjm/Pxrrf+olPo/pdQTtNY/B/AyAD/hf4cDOJePd85ZWZJgu4dtj3ITZV9Yv349AKCjGli2jPIxVik5BhQlTVuiExN9nnHQqPll6DpFqf/mKRKNClLXgxV5zGiAqZiZu7W2Zmn8SMkMRnR1mgP3S6KAsiqNnk55+lAtemRpqrZ6N9+F2D8X3WIi5oq69jojl1lIMgEJJSo6OvE6jnRl2VxtbOtH0b8aQy3Z0IO2D5m2it6QyxeGLjRPKaC0+ltTTq1cYeay2VlyYgRV4yFPfNoTAQBjrRjNlutUMipsDbkGaiEEglzbH7cxuf6H1xIzf8PVnwQAXHPkoeZe0e8fdgWFCb/qmNdSm3iDtZ9cL38oORQ9+tGPBkByPQzmy+dPAnADW7j8CsDrQez+FqXUUQB+C2D/oWoOCAgICFhQzOuDrrX+IYBn9PnpZUNVlo1jh796FqpN3wMAdFiZtWnDGhSSFKAk0y9h3amYIYk1QOXvJtsNKS0zsQT39+pXJjxnBRPW0zPjsqEvfR2kmDRpFOJ6b46eqy+zUMlzWOjSsC3labkSbc2dACBiEyc6Z1ZjqJgkRWAGw2zlfV/+VwDAB171ElNuBXdml/NKu0H+teFoGsbSTdiPtE3CwSoZE3GI4HboyPTZkkceL7FYKF13b2JBzMQNI3LHsTD1E7rMxhpNm+Bi2Q5/CQAYb04gThaJoW8FuT7/kF3x9k99HgBwyYGvoTL+zORarF0iTJpyt1W5vvy1ZLFy3FXX4/WfoHuFma9+wyHcAHpmNrlurXoEAFeuh0Fw/Q8ICAhYIhip63+SNbHjI56AsRkK3t4aJ/1iY82vDTNZ36Y5Js9p1pZAPSYjubAUYxtrWURldtxdRwWrk7U6WqtjFKbi2Qn3OGQoc57z7C9hN6FdhxZxMTaOBErbUJ2wukQ6J5zzL5TE+cyXPdf8rvzXY9iOa9ssiGOb9CHy5upckZ7OtxMua7pWozuU9huPDNcsTxi1uHerquamzeMm56KfzFl3LGOmY2UUkWLuXCj3/UmY0ZOuvRWD8PgX7Gn+v1iO/wsh1++65W4AwHkHU3+GkeuTbvssAOBje+++zcm1SaocWcY+X7m2Jt0PHrm+7JBDAHwKAHDFEcTMZUXjy7Xf1qrS0BxcTGWUILuE57MzBwJDDwgICFgiGG0KuihBOrYdlj2G7M+x418DACY2rMTMzAwA4NjTyFbznL2fBwDIDUuQ3WW2M02EoVv2Iam/jHcee7ZZ7aGwk7LHo64yR1fnKHpd8eQrq9LoRU2AJZNmTBiUlMWMB5G1lZcASMxG0sidrZsN6zWn2Q74nXd/BQDwoVfv7NyrlJuGLYvEHVn1pMHSHOznXZ8lY6SzXrMLN90yNpsMuDDX6OjqLVO2KTZBlXSFWKw2RKdp6pVbRC/LVh5a11gcHQuxjBDVakptvugIstud5tCiK1etwgkf+xgA4L+/RWOj4hXoHHYYFgPzkevol6QHvv/+PwHoleuz96Yk3u+6kdj2h/bfxZQ/SK5PZmb+0X12BwBU5bYn19bZtFeupZWJbqEOZXTqXL+xABos17JUKHWbm744cn3izTfjePqE4dhriKlffNiB0nw6slwXEqpC2L+KjIWPXOuq4axcAkMPCAgIWCIYKUMHFLRKES/bHgCw3YodAQDb71hg40Y3CE3JaaByTlybcTwPMwPJLFtaXaNhMhIDQvS1wlpqOjRJgyWoz+iAnYlF/yXVlGWJQkt6LWYu2rUptrv3dDj3q1/H6S95Mdcj+jvxePP0lLUQsn58o1O/8DXMhrd+vjdl3bmvfgGVlbgzvRamLP0ryxqTcW2Zy5o+l9os/6FDrGoX/aTHJtEC2weLgbCujPWBCcolbEe5z5b8btKUmF2raa1ZGiwXnWo4JrOwmFuup9f8HwDgvj/8AcBgub7ggF0BAKfecrcJjFR4OvS3334XAODCvYmZR+XiyLWCqunz+8u1OGH2l2t+NvesaZihWxv6rmlH5THWbUWuT/3Mp+k2DVx08H4AgEteS207+bqbAQAfOfRA59l+ci2yvblyPfIPeqUyxEbRT4PVaPw1Gg33g17FtAyrctf8qeAXKs4Ncf0FS4wdGWzZvPGWRfSgVdkAdtsl58VU15gYyTn9XmqY5VAcS9AgWRLKUlQ+TFYIEm7c+7/2r5gNp9/1jYG/nfX0V1IbeKmWl7Sc/+CP6ZlTn/RUqitJAN4sisV/t7nWKauc2sht5KZrBcVR9MqChOnCfyPn4BNfTBu14rJdipMIH8tYWbOuikZSNrgU/D84+/E2pm7yMTCONHRPp5BcknxbizaKWiseYcorCpKTpOrWPG5GjbnlOm7Shul85fqD+74aWn2B7tFfAmAdYc7bazc658h/iyXXSc3hR7QiJsdom8vnFtTlwP/oqC59xKxci2qE/445BjnJdeW0UeRaJh3dmOQ2ctl95NqMTSqqjs2X63d8htReF+yzF995p8nbKnJ90UGU7+GU6+nDfs7Br3PaUZfrRmMHAJ5cD4GgcgkICAhYIhgpQycSoMwM+IRXvnrgvWd+9ut9r79vP1qSmqVeZc2SzF6GbPSU3pKqHvjIC14kjEWu59plMh2ZTbXNpG43bVzzJ+0tkQEg79BmzenPpPjnBbMG2Ru74Ec/6NvfOs783j/O+vuH/ovKOOmvnopLfvvdWe896yvfGfjbSU93fci6YmonG1vikCGbdBHMUt+6U7PaS7lqr8q4UldQlWRpYpWYMWN0Xc7FKTtOiQE3xlbgni9+FQCwkjcIU53PO8HPQsOXaxPlQDURcd8yZmGNMXIU0VNuWNTKuM5buf7APvL3QUz9g3vxRuk2ItdKV4aRy/68L9dGbcJsXmuNTJyApL6cHMVMPtJSzPtYNVJvj3JVlNEYG0nIyiSjlZCoMQqtjDx1c9fcs8sSszly/c7PkNrr3New2qumvpGVhi/XHz6A4tuffiOFCXj//hQmoC7XaTrOY2HlehgEhh4QEBCwRDBiHTqZDEo4yf/84hcBAE09Y1jgk/Y8AABw6dG0gfDAGlICGwcF0cUZNm5ZROkx9NiEKpVNUduOKLKshp4pvWPNxLFeqNYoS2GdfOT+fPQbg1nxuQN+O+GxT+17/cTHPh0Jm3d95JfuZuebH/V87he17eLfuL938zaOffjjAdiNnw53/vo1vwIAHPzwRwEAsozqaJczuPl3vwcArF13v1Ne3iE2qSVPYyzZym1YVWFtOYd8NeZX6G86qgGTP9bsgyj3Ht+dXPSmjUZmNhMN/NTqI0ZdrgV5BVScRGHZCgq6NMHBujpsziiJIYz+vSbX776DmPkZXN4776S/lw/u/iq+dTRybdZUJev5mYVPdzumwbyfWnPiE9ZrV2TgK5V28wbbJChylTc0+W9bsjvleQHJWWrkekZ06lRfh3/PMupDu5wxslJhmqoR00rFq6XNkOsP7r2H09/TP/05AMDZe++OQs0u1x/Y//D65QWV68DQAwICApYIRszQNWie8kyrVAQVuTPRylWUc3rTekoaYBw0NAe7l13tuumRF+pSXH0NazBkRNs8g8aEyZo50XX53Q05UFUVrvgO6aqP2unJAIA0dufF4578FABAh5nt1T//GY54zBOoDaWXbqvr6kMF3W6B3HMRFkgIVlmdnPDIpwMALv3997hsoNKuY4fx2pHTXIIZUTtuvu/35rcbfvW/zr1Xf+8/+7bjlBfvDIDGxGaRF70huI1yN1u31ANwKdHviomb18bC1dmuHCcWOzEx7qT6MvUslpHLPOS61SKrhfnK9Rl3fAFnvYaY+Pv2FJ31F7lcsno569Uv42e5vi2U67wjAbToN1+uhUmLXFfdyrr6+3It5ote3jyttZFrYy3j1uKsuqnt0q7Bci0LDsktKnLdqSq7GuHQGAkzclWIAxPXw+xe9O/DyLVpO/TCy/UQCAw9ICAgYIlg5HboURT7kxoqHSFKXd3R8lUUbrRx7+8AAFPTbF9akrVIFNEsWmpVc0CQgEPMyOEFwanNssbW1KRKo4PszlfMCmoX+GgZZmeSdHK5N4tuWk9trQfh0ayfFP3jdb//BWbD6t/+x8DfPva7bw/8DQAu/8NPe64d/Rc7Oee33v+rWcvohyOfRckkWjHpgxuGbVteYBiaMA1JjmAsKGDORVcqAZfkXVjLGPo9YcaUNKi+rGV1xWXJTieV7xQ+Sswt1xOsQ59Lrs/4zD8DAM7aaxfDnqHcVdz7dn0JAOA9d/8Tne+yM4Atl+uyTWMpcpvXmD9gY1nJ74nKjP5X5Nowcmm6YfDW2qni7N8RO5mZVG3ChnkzLEkkbR9bOcWxoeYR19s0dvcsX7HrPBTHY6YcJU5cYkUVc98lKJhJyjF/uX73p8kO/QN772b6sNByPQwCQw8ICAhYIhi5lYtSUQ/jkNkVAH79+RsBANvfS+m8/sD6pfvW3Ec38AyZpqJrqnXBN0Q31/lyvU4/j1ePFbPV+dF/JCi9LcVU59kFW08xW6ZYbKho6/PItz72eUDieiKOsY7zpMdS0LNL/oe8S497/LMBAGmrxCX/cQ8A4MRnEJv/2D0/BAAc/dy/pXuYwiReurF6PcZDz7jz07l9Jbb/hpGJnbG5w9VbCiaWjZmjsClhOyRDi8jR55Dr1hjZF2+/A3kC+nL9/i/R+zh7HwrSZTIZA0bQztzz5QCA995FAcmwmzB1ssl/76t33iK5jmSfw5Nr4zJfyqALo4175NowcnMP/Z5KmFutUXmfHct+2ZItcZ9tZIm9z5PrvPTCaDfcMMVpywb+0qmrm++kNqQAsHlybcrWW1Ou54/A0AMCAgKWCEbO0EutEWm32qqIbIwECVTD9rorV5LOsdkkJtNmvbWJ8RSP1wryGFJNX1uHhraB/82RZ0TRu/NO+yXfHOzBef0v/6fv9Rt/36ufvuY3/wUAOO4xpId+42OJ9Zp1BtvAfuTnPwIAnPo3T0VScvJYZg5n/4LszU/7G/LkTGOyjIh51/69//EvAIAVyzMglmQE4jnItr8eYXvIQzloVmbHbsUyL6RvKklsxSpBQoryeaTMykUsJMTiQ3t2tEWdoYtHHp9LUmjxxhM29O5brgMA3Pa886m+pISKJE6PMLa4t3MjxJbK9WkvIw/iJnsUzybXpk4eyzM5rPIwcm2O7AugUWGMdbkSkVb6Y+KzcGwdE1ZXK8vw5Z1xf8UHJGHWnYheWSkkJhous2m4TNnKNclOkrDVidZzyrV4joqnbJxVyNmaRfJ2mMTSWFy5FpvzZovHaKBczx+jd/2vKsAL2KSgUIqagpdb4xNk9L+CN5PGxmhZkk+TwBfGHb1WjpfJpQd1C6o5/vblRb7pOX9H9XJwnqqqcMU9/w0AeDtvEsYsGOd+hz7Gp+zkOgvpQuMjPyb1RaMhH2kRfN7AjdzFUtaIkXmmYIJmkwUgFoF3x7PZSFBAlsvu8s+vZ9kyXgrXsh2lLBXvfgn1L1Num1G4gqmUMn/k8ocEMTUD/XH2yz5vnY3kvcEp1zfhynjpnTVS6EiWvu5SfzEwark+c3dStZx1F6lazthtZ9uQecq1UR3W1JQZt1HM+kSuS5NL1DVc0IXN/KX8d+fJdVJzs/flOjK/9ZfrupzNJddxizcgYyvXSc4Th/mQ8ybkFsj1e+6kMBzv25Nj0deips5XrqXNCynXQeUSEBAQsEQwWpWL1kBVofD3LEsgZfOuLi8FC27aqoc+FACw/cMfDgDYuG49gFpuxqoyjKKbd7kambUlHKtbXxxFZuNSlnVKdDjisCSBesz0Ks4Jdg5ssJoi0m4FYy1vmaTtM82mmGJxdSY7uvtIlFRQzJqVx6qj1DPT9IIo5VEBBbr27m8ODscLAG+/+5s919KEwplKHsiicBli1RMbukDELu6Vya4jkCxKfG+fNtiMU8xctJjCuaqfFatoQzFOm+h4TldVtYjcZLHkWsAbklsq1y2WIylf5FonojLolWth7xLKd5Bcy4dGRb1yLbl+o7i/XGtrT2jkWgKIxYnLsiVQWaRthqGUx1H6ujByLR0SlVYv5pLrR6yi+PkLKdeBoQcEBAQsEYx8UxToNS6LVGQ2KmRCqngubLGOcYcdKQvMUV+8CgBw6stoYzDP8x69sIG3KWp0kdCGOsikHLFeS0KIxoUbiCjx9F6A1XN61mo4+5v39G8PgPN+0H+T9cxnPd05T5IYMVydnumW2RGWQ++mr2SKOfeFHMgrcss6jRNtXMgbalVideixbChx/2JxI/dctuvZcKwutVdX7ja+9l+vX2exk8YgvPjgQ3uu/ezfv85twWJaLVIbvPNh5PqBFZTJ6P77KTDafOT63Xu8CADw/s/RZviZe+28RXKdmIxfXI2Ra2bSvpt7FRm9d+Vl3oq0+0dRl2VfroX5+05J1t3ejqzItQn9Ye0J6WDaaLYkTX8kW9aCynUf69Bevbesjtx3v2LlSgDAxARtlCdJitxrU3D9DwgICPgzxaKYLfrBd+rWDF1DFHnnd4zDjq7c0X0mo2fyzjRiNmVSJj0UQSa3wrgBS+Aem83bzMA8e/tJFkSPKYiiCO95KQXfErMrVZEe770706qh8qZxVSV439dJV33mC57n/DaIVVZJatKG+TP+Gd8anJwCAN797//Sc+2sF1Ju0cgLuFRy8KhIWdd15dECpab4yE2uJDiY6B4LY74mbttGl+rp90uTPCAyikfp35l7v8Y5Hxt/GADgHddfBgD44RfvBgCsWrUd2qJrlHp78sGPFlsq1+MryPpl3UYK2jWMXKN2viVyrUBOOL5ca+3qxwV1uYzYYsRA9z/VQI9cl8omvwBgcpnaiuz9msdPVh5lyanhJG2eXZrw+cLK9Qc+928AgDP2JAewgnX6s8m1j//f3pUGy1Fd5+90z5t5epIQT2ADBmLjikNMiM0WAlhYu4TE6tgVQ/kHMVTJGAirCwljJAMOBQSDFQLEKruIs5QhYIwxuxbEYmLFQMCAMbG8BAQCyWUjgfSWmembH+ec233vzHtvxk+vezLcr0rqmZ7uvqdvf6/vueeepVTidZVKhfs7lnWWuqGMh6paE0Lof0BAQMB7Erlr6MaYBj9xoshqFMbzGNHI1ylTd3P2903mwIvtQ+/Y9JqqnagW4o9tajfkNl07ntrrbGFp61GQ2uL0+DguO+eirhqFagVwrkHZyt2i/Vz5xJMYDV/90ejeKaPh0qOORk/k2fxtVXbf3p7aIG2fwH8+kGvp/Wr4eGK3NsRf/a5t6S7Xtmoo1SAbbMRev/l+6JRJQ6rH2DJltWrDWkKeGC+v95SUAL/7HRc9bofXSxfzOsnVDzyOy46fKW3nw2v7HLz7s/ZvDYPPpOmFN1tIIrcknam7mrLeMRGlvPYftfptW812gngt0LTE9n7/AF7bNL5j8LodtKShE9GFRPQSEb1IRN8lol4iOoCINhDRL4joDiIqj32lgICAgICJwpgaOhHtC+A8AAcZYwaI6D8AnApgMYAbjTG3E9E/ATgTwK1jN5nYFXdFNTMIJTK6HX7iiaNe5fK71jXuW3gUAFjf0br4farukGTGr6r6s6rtyyuOYQsxyCkxpdpL5EWYmcgf+V3txKo2AEzExQEuO/ZwOUbsiDLiX/sUe8FcdMRh6sFtNYmVP3kWAHD2kWyrJw0pFjXi1mfZu2ZbdSd6ZfQviebQE5UdUa5ayDb1UqYfNBG/8WgRkav1aOk4DR9Pkjj1jPCi42oN+qT0N0X23rVUF3nH1j2fXC0ePFxNbDi8JlQysb9ykTda43WqdbK0lclsO+8Xf/T+d94BAAxv35HRzlSj7UBea8RmIkUvdEYi7WuRh3omiVbKa0nTq371VvN3eU2UduRYvFY/dWXSrua1ombGz2vrpz4Gr9tBqzb0EoBJRFQC0AdgM4A5AO6S378D4JS2Wg4ICAgI2KUYU0M3xrxORNcDeBXAAIBHADwD4G2TLutvArBvs/OJaAmAJQCw9/vez6PSCOltgTSK6sd3fx8AMCkedg79+MlcPPp711wCANi44Wls374dAJDUNA+Gjsje5TOrzr4Jzrd/+t4gUcaOqRqEjtrkD4u+YzqAry08lmWrSdSfRtZZm7Zv204T4ydeet5BjSaEa3tU1E1iNTG1bfaoppGqMu42c53Gx+N6WRC5MvN5zR10Uy2LGs6x+8iTSaDJix647hYAwF6yP4pi1BPVorRv8nVE/0N5rX3cEyt3+BlO2Y1t6fvsy39GOzZttryutcjrS47/JK67n/3yLzlhhtOeYlfzWicAsRaWHoHX2W8+r4e8QjQ+r1XrjqKoMF5f/zAnxtP1iroZP699L5iRed06WjG59AM4GcABAN4GcCeARU0ObdqyMWYVgFUA8NE//oiBMYBxDf21JF2sTJDmTQaAgWGeymmdP8WU3dnta/oe0zAwwNPUgSG3fqIG09gHGqeETSxZ5D79k9baiQAAFgdJREFUKZT8YF9IUTb4wWVNusCq19ckQPpAqzYAQsOqdb1PTS3ZBVtuL7Gh10niy6bZ6SbJ/fHuL/7FnwEAeiJCT5mP0UyMuiXbfsWROXKmhTvc5nTaZ9ww57qtzp4uGmlSJJ37adUb/XtzCKd9brPQuX8kanGpSUy9UAG1aoQk0QUlvVYdjYPKxGFX87p3Ci/yj5fXaXvqAtgZvNbgIYNGXidmSK6hg4vLax2D4gJ5rdCv+fK6dbRicpkH4NfGmK3GmCqAuwEcA2B3McEAwH4A3mir5YCAgICAXYpW3BZfBXAUEfWBTS5zATwN4FEAnwFwO4DTAfxgrAtldYAshoYGrStPTcOM7SISHzM87Go/UySvdH//7ti6dQsAYGDIdXfyp5uqLURRZDUIO+2x7mXi1qUy+wl7TIJMqiE5JnKOJS/vcfam48itimJnefLh0mOP5ParVZBXOf2Mww4U+eVasVdpPUpd0nrKMiXVmY31JnSDOpz7s/3lJ3CSTeQLnT0145YGwE44rOth45PXRSn9JbK1HSPn/qyWZ6fxZLU537WtCOwKXlckde14eX3BcVyB6vqHfgQA+NLCI51ziuK1LooakzTwOrJavlwrJ16vuJ9r8y5fdIQndIqbHnsGAHDeXO7HmqQ4+MYa7t+LJHVGFkXyekwN3RizAbz4+SyAF+ScVQCWAriIiDYC2APAt9tqOSAgICBgl6KlwCJjzAoAK7zdvwJwZHvNGaBWRd1LCdkT92asd6I5J1ohRGXg7Y/v/GcAacrN6R+YimnbONnN79/makY7pVhAVOHkNzpa18RlsIQIRqpq2+oi9V5HpmrsBiqoaxgZgIxb6STS8GMdgTXbpy4YUSUT3KDah9o4h+QYCQkvizZChB7RZBKpODGpprZHsVeTW2m9FLFccRyjpIn/E9UMWRZN5xqL/VArr3MxBrFD2uAMRqK2RtU0jIZZ6wJYjy2fWiMNQhFNVK/p1Ws0AK6590EAwOUnLGDZKqrBSDtVDoOn4aly0pC0N4gErhZXpqjBtSw/jJ/Xw1KdKhKb8Hh5bSWrun/igddp3684mV9fV/6AU2ksPU5mxxleK4aF17eu4RQeZ8/jRWeb8TfD69hWSGLhxsvrdhBC/wMCAgK6BIWkz/UtjhFRxsbX3E5oF31VFZQRbPKUPrx/L0lB+gaPdJqCNJYK4ZWKlEFLWMMZGh6yKWIVfnXtxNrktN3U5l4y7m/ZMl78VbQE1VJNZA+2VdJtRXBXc40k2IIHd7VdSjtqX9Pk/p6t0dZkLJUQqyZmbYpyTuRqC6mpkZzPWVjbrU2q5NpWs4FTqRHV9ZhQ2HqOGc0j8uyvsfXaUE1JbY264p8+K9eunJ/bYnMUz+ub13JgWjKXE8h9fR2nkLhw1tG8fxReX7f6YQDAZfMXiizvDV5ffhLb0K+9lzX1i+cdpkLbY29dy5r5OQtmyMl58rp1BA09ICAgoEuQb5Fowy6apolvZepRoaHDzYMBEpswiDFp6r6YNI0P6p/GI94729iO+PfrOQHWl0Q7iXrS8Ou6519qzKD7Xc3GkWc4NEBdfIojOx6qXVS1Md8LoG7P14lAuoadOFsHnkk46nF3WBuj2F1V1iiObHpRP1FQVNKgBtF01OsiSex6QuKFXquvcRryLdfUjaHM83NX50s2GEaLbIvdt5ber9o7S6Kpaci02t9teLlR7wCgLu1ogYY66oXp53nzemCQ/alj4V2W1+fM/BgA4OZ1zwMAbpLrVYXfI/H6H1Y/jaXzZvFvGgL/HuP1suNZ+/76/fzeIJMq4l+cL++QFnhdiTXF767hdTsIGnpAQEBAl6CA9LnZz402I3i2Rt+27Ws2BhH6p3Pq0S2TtwIABgeHnGPUdlXVtJhJM12ueTtp6HpqN05HeNEGYtfGnFi7qY6Xo+iOnqnZ+t7GsS34a+9DNRVrm3O1lWxazvSY5j6w5PncGqKM/c+1lTfYHnWvtuH0p9sH6llwk/hFN8OK7z/YdP/Kz3HJOS2cnMYipPJUa5qMql6oCT1XXsdqax6Z16qpG/opAIAM29bPnnWItM/n/OM63n/evCNQCbzmvZlL/O38v3R+9Xnd1K9K1w2i9J6zMita5nUbCBp6QEBAQJegoCLR7qiTZEp3JXU3p4nRpD921HTTkNZKZQyDbYuT+vi3yVPd26Jh9g6o9Eox1jjBTi/vhvFsmH7X2IT5Bki8FKWabIjkHE1YVLPHlTI+0ibzP2C0/JbaWCUqLkGamB8274ab3pQizWMROVtDEYxqIQ2JgmLn2DTdaZIpuOtpYJH47KsGqAU7rEcD2aIA5GuicqVz5nHpPfXsjaIINzzCtsqrP7sYADCp4hY+qFXZ/7qOAb6Wke/JDlDUhywSJIUWuADy4/W7b28GAFAP98tovD575sHy6UWRkXHzetbMz5Fi63nx+qurN+Ar846WTmmP18t+uBrXniJppMbB62X3PITRoHbz8xfMwMrVzNGzRWaf1zbZmq4zGEBM6HZLNCRb4bU861Z53Q6KqVjk/eFlk081uutEI+xnVGs1zD7jQmffaV6b1zzxSoMcX5z7UXfHSDEp9p0qDyyKRpxOjrgfZInnh9On7k/+NSjNBqdTNrgkjvypqp2ORg3CWKJr+8b7A2xyrF7BVpbxaksmdlGHGkxUlvga1OFNO8uZZGuRVzMz0wUAgKpkqKxWNWgmzWhXtzkG2pua7mpMBK8rFQ52K1e4r0ol7sNJvVpLhq8/LNmdSqVSI491UW82v9hvXf8cb+0BT9tPF8mgO5G8Xr7gaHztEXYBvHz+MfLT6Ly+7P7VAIBrTl6UGUBG5/UFt4+cieS6T7MC8W70LoA0BP/v7mEXz4sWslwmAc6dy59vWcO/fWEWBx+Nxms1hRXB62ByCQgICOgSFKKhW+tFkzzF6WfyjtER2N1fjfrw2L99BwCw8zesffz0hRcAAEvvWMPbT3yQz+np5y0AeBXa6+SrNlYiaVbyMCNBr1RJsdpBJMmTdDXFVkjRW0nNJ8ZrRxUZTYlqExdFBBh5PDr1VStK5roAQGl+UduGyqbH1Pzcy87dyXQ6bj6+x8RVdZKkKu27FWbYfU1mMNKM5os2NXGXE9esnl5OE1vKVMOxbl4qjSglfmUXQKauGLIdly4aGRSZoGsieK3om8qa+tR+fg7VHW/ysRKQUo65T03d2CAdhfL6m4/y30R9ziHynf9Wlsw51Ip1w5r269he8am5In7rvF6+cBYA4KqH1wMAVhw3U0UQWfjTV0Qzv/pTi2R3yuuld943qlzXnXoygHRWmUWiJhHh9RV3PcDXPJ7vBVZjTnl99hxeHP3mOk7otWQGz3ia8Vr5vOt43TqChh4QEBDQJShmUbQhnDUTmDKSBmPPhXdcar8qSbKf3t6Kc47auVRbaepw5NsCrYuT7LdW5kbtt8FIZi+Z+d03oLm30SCHtpVK0uRU1X4aFojI/ZzBiCHFI85Q0OAqZvtElctMtaB0cZRRzaR2zTYTxxGu/WvWvCoanFJ3F0UVqtEkTo1Rdf+T4Bq4GnERyIvXlQrPENX2OxqvVTNfMvvPtSUAwBdEU18lbotL5h6Ki4/jqlrp4mTV+a51LvXvKaYYK76/tklPjI7li2bzVjT1Kx9aDwD4qqSiXf4Az6xT98FGt9ZrP3uS892mWPC4Mxqvl9/BmvnyT3NyOBrS9R+5ZhNenzWLF5FXree1h7M++TGnmTiOGv5OxsvrdhA09ICAgIAuQc4auoFBFX5FtThO7X61qtiX9CA/6bwk+YEdkUtWUzGaSnNS2TnFehtkEuck5Hpj1MRenS7Wu6vXSckO26hZe7GOxK7dUhMtVWsaiGEySYrcMZSM6+2gQRsUUSqDXlfSpMbi7aD2tsT2UbrSX/OTjUm3WY+VhuCqkTUZlLUdscdLO7qtmwS3rGHb4pnHssai2seUPrb/ahh3jw2HJpsKFZLy1UaAa19INfS+EpckozrbKas7S4in7CY3pHwYdp5vvsiX17Z02yi8/ub6lwAAZ4gGaWo6w3R5/fk5nIjqtrXP4oLjOPS9HV4vO3GOnOMeS+I+aXlt0lnelQ+uc469yn5aDwC4Qvf/lSYJS+9PbfVD4jo5Hl5fcdqJfI7YzJvx2s6CJJWw8vpcmeH05sHrNhA09ICAgIAuQUHpc13U6/WMNuMbpH1bo19tOy1HJXnybVrRqxezdgKpnj6qFiq/3fYEh0qffgyPwCqX1dwzKVGt07/nkzpaOw1mVtd9dow1bd/+6ntKpDImDTZa8W9Vu11Dua2Rx/akqgmIPF9j9WipNfrK+omVGlOiUuOagLejUubnuGMHJ6My2/g5Tu+fjkQK/vb19aUyFG1E9zBRvNZ73t6E16qZL5l5kFzY7dO6x1Xtss/PORy3PcRBNKqpTxSvLz9hvnPoleLNcsUpswAAK+5Zz/tPUW8qab0AXvtrHYXwug0EDT0gICCgS5Bz+lyDer1uR1VFvWZSQ+EILiN+Mh4bUYWq/dzXzwUBot++DSAtG4VeHu2M9biIrWdFzfqE+mOb+J1r5GNdR+bUxzqrOQCwofGJ531CmXHTb4W8iDrKqDax10+qbWn2TRP5PraZ755N07fvNmoyqSbgN4uqm1pYtSQNba9nVum15JhqlyW5YS0jpqlES1FsyRfb9QW3cHefhIBH0t6W16SQw7bNmPZH7D89bdo0AMDk3abAII0+zROdxOvz5rJ/dKmuhaUl3Nw+o5F5fe4CDcl3n8NE8frqB9lD5suLP8nHSNaCL5/I/uBX37Pa+c4Hvfd43Q6Chh4QEBDQJegIG3o2Ja2OMWrXoqi5XTTr471zJye3eef11wEAAwMDcq7nQSLfS6VY6w3Y1J21xB2+bV4O2aRpdGO78u1rVfZcTbGZ8XKw3jKeStF4e+kO/9jUM0Lshd6KflYOP6+Irz02ajKZY/2JgXoUiDwrHxw5onDV48823f/lkzgaUG2QRJTmusjYH7NQDUn7efMbnJTq5Z+9jMqv2O544IEHAgD2P2Cq9XnvFHQSr9VzxZ9dFsFr1cyXLZwtQrh5S7SPLlnEnLnuh3z80sWzJpTXeq1aEx4pb8tl9jTqidy+mEhet4OCXujuw3e8zYRwibg9JeK109OjGdjEdUtyKg8OTsKLL7wKAHjzNe6Md3fwSX2Q6UrM2wreAQAMmWGsXPPiqBL+y1MvNN3/NzMOR9WaWkTkSBePeIfSXae1lDnWR+wl34kyfUMNObO9jJCemShrgmkkthcoRW67owXN17VSvfxBL5nNYdC2FqKp4VuPcyj5BQv4NyV+r4hcFlL3RBqkQtYtzuaa90SOeniBaFBDsatM7rc2vYFfPrcNALDldX7mMxYfjuFqsQm6OoHXJdLMiNKeE3SUTv81JL4qcfAE5MLrGx5+EhfPZxMLLK9G5/WF82cBAG58YD3OXzDba2nX8bpe06Ce9Nr6fMrSBz0iaiV2X9YTyet2EEwuAQEBAV2CjjC5AKZh2uNPVcvlivP93Xc59eV/bngWjz32GJ8zzFPUwUFeZDj4w6zB7L333gCA+s7t0obBlxYc6rRTrfHixk0S0rtkxkH2WACoQhagEoOahuVGrtugDfgQfehbj/0EAHDWnJn2Tv3wXyJfh8h8b6h7ytuVjz4OADhv5kz39ygNOElG0Jy0fX/6ecvaJ3DW3GObnlOru8ENRi5uc0GX0ntSjUZTiJZK/v3qdBRpeuBMrvnsMbpD3fUU9XrdTos3btwIANjz55MxODjQVP7ikD+vyTOXRN5WGaX1QTWQrR1e2/S2mQXKVnl9wfxjgLqrRaeZbz1fR/1deH3OzE8gGW4eaDMSr0fT0FvhtQ2Ik2C+InjdDoKGHhAQENAlyNltEeLe1WjrbFyI4FGsVGJtJI5Zg97yFtuYnn+eq5o//uRT+M2r/wsAGBjmY3rLHG7+J2KX3CnVvvvKGpKbAOKwr+5U5KUdLccadKAaVkWkMhgQ8UmDQRJXc/G146Fa6oJGXqBAzXNpooy20mhDdy88ZFPTmoZ2E+/YVINKMv+7GBzWZ+CO84lUztG+0GcRxZHznT+LLVE0GJvKQO5P3b/IRLbCj60qb+OlxKYpyaGGRaMrEaeNndq3DZXKTjmYXfl++epkDA25tWTzwv8HXouiiaoZP69tOlgzMbz2i3+Mh9fu32L7vNZ9cazPLX9et4OgoQcEBAR0Cagx5ecENka0FcAOAL/NrdHWsCeCTK2g02X6oDHmfXkLEHjdNjpRrk6XqSVu5/pCBwAietoYc0SujY6BIFNrCDJ1vhxZdKJMQGfK1S0yBZNLQEBAQJcgvNADAgICugRFvNBXFdDmWAgytYYg08joFDmy6ESZgM6Uqytkyt2GHhAQEBAwMQgml4CAgIAuQXihBwQEBHQJcnuhE9FxRPQKEW0komV5tevJsD8RPUpELxPRS0R0vuyfTkSriegXsu0vQLaYiP6biO6T7wcQ0QaR6Q4iKo91jQmQaXciuouIfi59dnTRfUVEF8qze5GIvktEvUX3VeD2mLJ1FLe7mde5vNCJ80neDGARgIMAnEZEB+XRtocagIuNMR8FcBSAc0SOZQDWGmM+AmCtfM8b5wN4OfP9WgA3iky/B3BmATKtBPCQMeZPAXxc5Cusr4hoXwDnATjCGHMwOKPrqSiwrwK3W0Kncbt7eW2MmfB/AI4G8HDm+6UALs2j7THk+gGA+QBeAbCP7NsHwCs5y7EfmERzANwHTjX9WwClZv2Xk0y7Afg1ZOE8s7+wvgKwL4DXAEwH5yG6D8DCIvsqcHtMOTqK293O67xMLiqwYpPsKwxE9CEAhwLYAGAvY8xmAJDt+3MW5xsALkGaW2gPAG8bYzSzUxH99WEAWwHcJtPlbxHRZBTYV8aY1wFcD+BVAJsBbAPwDIrtq8Dt0dFp3O5qXuf1Qm+Wnbswf0kimgLgewAuMMZsL0oOkeUEAFuMMc9kdzc5NO/+KgE4DMCtxphDwblKCrEPK8SueTKAAwB8AMBksKnDR5591QnPyiJwe0x0Na/zeqFvArB/5vt+AN7IqW0HRNQDJvy/G2Pult1vEdE+8vs+ALbkKNInAJxERL8BcDt4avoNALsTkaY3LqK/NgHYZIzZIN/vAv8hFNlX8wD82hiz1RhTBXA3gGNQbF8Fbo+MTuR2V/M6rxf6TwB8RFZty2CD/705tW1BRATg2wBeNsbckPnpXgCny+fTwfbHXGCMudQYs58x5kPgfllnjPkcgEcBfKYImUSuNwG8RkQHyq65AH6GAvsKPCU9ioj65FmqTEX2VeD2COhEbnc9r3M0/C8G8D8Afgngsrza9WSYAZ62/BTAc/JvMdiutxbAL2Q7vSD5ZgG4Tz5/GMB/AdgI4E4AlQLkOQTA09Jf9wDoL7qvAFwB4OcAXgTwrwAqRfdV4HZL8nUMt7uZ1yH0PyAgIKBLECJFAwICAroE4YUeEBAQ0CUIL/SAgICALkF4oQcEBAR0CcILPSAgIKBLEF7oAQEBAV2C8EIPCAgI6BL8H78Og39hWkMhAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from rppg_ltss import build_and_plot_mask\n",
"build_and_plot_mask(face)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Now the mean color is computed within the green area.\n",
"Note that eyes have been omitted so that this region is more \"stable\" (blinking can interfere)\n",
"\n",
"Now let's process the whole sequence"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import bob.io.video\n",
"video = bob.io.base.load('data/andre-real.avi')\n",
"\n",
"from rppg_ltss import process_one_sequence\n",
"real_pulse = process_one_sequence(video)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's do the same thing with an attack\n",
"\n",
"![attack with a mask](data/frame1-attack.png)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"video = bob.io.base.load('data/maskattack.avi')\n",
"attack_pulse = process_one_sequence(video)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extracting features from pulse signals"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from rppg_ltss import get_ltss\n",
"\n",
"real_green_pulse = real_pulse[:, 1]\n",
"\n",
"get_ltss(real_green_pulse)"
]
},
{
"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.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
#!/usr/bin/env python
# encoding: utf-8
import numpy
from matplotlib import pyplot
import bob.io.base
import bob.io.image
import bob.ip.facedetect
import bob.ip.dlib
from bob.rppg.base.utils import crop_face
from bob.rppg.base.utils import build_bandpass_filter
from bob.rppg.cvpr14.extract_utils import kp66_to_mask
from bob.rppg.cvpr14.extract_utils import compute_average_colors_mask
from bob.rppg.cvpr14.filter_utils import detrend
from bob.rppg.cvpr14.filter_utils import average
def load_and_plot_frame(image):
bob.io.image.imshow(image)
def detect_and_plot_face(image):
bbox, quality = bob.ip.facedetect.detect_single_face(image)
face = crop_face(image, bbox, bbox.size[1])
bob.io.image.imshow(face)
return face
def build_and_plot_mask(face):
# get the dlib keypoints
detector = bob.ip.dlib.DlibLandmarkExtraction()
ldms = detector(face)
# image with all keypoints
display_ldms = numpy.copy(face)
for p in ldms:
bob.ip.draw.plus(display_ldms, p, radius=3, color=(255, 0, 0))
# get the mask
ldms = numpy.array(ldms)
mask_points, mask = kp66_to_mask(face, ldms, 10, False)
# image with mask keypoints and lines
face_mask = numpy.copy(face)
for k in list(range(len(mask_points))):
bob.ip.draw.cross(face_mask, (int(mask_points[k][0]), int(mask_points[k][1])), 4, (255,0,0))
for k in list(range(len(mask_points)-1)):
bob.ip.draw.line(face_mask, (mask_points[k][0], mask_points[k][1]), (mask_points[k+1][0], mask_points[k+1][1]), (0,255,0))
bob.ip.draw.line(face_mask, (mask_points[0][0], mask_points[0][1]), (mask_points[8][0], mask_points[8][1]), (0,255,0))
# display
f, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True)
ax1.imshow(numpy.rollaxis(numpy.rollaxis(display_ldms, 2),2))
ax2.imshow(numpy.rollaxis(numpy.rollaxis(face_mask, 2),2))
pyplot.show()
def process_one_sequence(video):
nb_frames = video.shape[0]
face_color = numpy.zeros((nb_frames, 3), dtype='float64')
bandpass_filter = build_bandpass_filter(30, 32, plot=False)
detector = bob.ip.dlib.DlibLandmarkExtraction()
for i, frame in enumerate(video):
ldms = detector(frame)
ldms = numpy.array(ldms)
mask_points, mask = kp66_to_mask(frame, ldms, 10, False)
face_color[i] = compute_average_colors_mask(frame, mask, False)
pulse = numpy.zeros((nb_frames, 3), dtype='float64')
for i in range(3):
detrended = detrend(face_color[:, i], 300)
averaged = average(detrended, 5)
from scipy.signal import filtfilt
pulse[:, i] = filtfilt(bandpass_filter, numpy.array([1]), averaged)
colors = ['r', 'g', 'b']
for i in range(3):
f, ax = pyplot.subplots(1, 2, figsize=(20, 5))
ax[0].plot(range(face_color.shape[0]), face_color[:, i], colors[i])
ax[0].set_title('Original color signal')
ax[1].plot(range(face_color.shape[0]), pulse[:, i], colors[i])
ax[1].set_title('Pulse signal')
pyplot.show()
return pulse
def get_ltss(signal):
from scipy.fftpack import rfft
window_size = 30
nfft = 128
window_stride = int(window_size / 2)
# log-magnitude of DFT coefficients
log_mags = []
# go through windows
for w in range(0, (signal.shape[0] - window_size), window_stride):
# n is even, as a consequence the fft is as follows [y(0), Re(y(1)), Im(y(1)), ..., Re(y(n/2))]
# i.e. each coefficient, except first and last, is represented by two numbers (real + imaginary)
fft = rfft(signal[w:w+window_size], n=nfft)
print(fft)
print(fft.shape)
# the magnitude is the norm of the complex numbers, so its size is n/2 + 1
mags = numpy.zeros((int(nfft/2) + 1), dtype=numpy.float64)
# first coeff is real
if abs(fft[0]) < 1:
mags[0] = 1
else:
mags[0] = abs(fft[0])
# go through coeffs 2 to n/2
index = 1
for i in range(1, (fft.shape[0]-1), 2):
mags[index] = numpy.sqrt(fft[i]**2 + fft[i+1]**2)
if mags[index] < 1:
mags[index] = 1
index += 1
# last coeff is real too
if abs(fft[-1]) < 1:
mags[index] = 1
else:
mags[index] = abs(fft[-1])
log_mags.append(numpy.log(mags))
# build final feature
log_mags = numpy.array(log_mags)
mean = numpy.mean(log_mags, axis=0)
print(mean)
#std = numpy.std(log_mags, axis=0)
#ltss = numpy.concatenate([mean, std])
#return ltss
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment