A ipython notebook. Takes an image, converts its skeleton to a networkx graph

SkeletonToGraph.ipynb

{
 "metadata": {
  "name": "From shape skeleton to graph"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import Image, ImageDraw, ImageFont\n\nimport mahotas as mh\nfrom mahotas import polygon\nimport pymorph as pm\n\nimport networkx as nx\n\nfrom scipy import ndimage as nd\nimport skimage.transform as transform\nimport skimage.io as sio\nimport scipy.misc as sm\n\nimport os\nimport math\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def branchedPoints(skel):\n    branch1=np.array([[2, 1, 2], [1, 1, 1], [2, 2, 2]])\n    branch2=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 1]])\n    branch3=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 2]])\n    branch4=np.array([[2, 1, 2], [1, 1, 2], [2, 1, 2]])\n    branch5=np.array([[1, 2, 2], [2, 1, 2], [1, 2, 1]])\n    branch6=np.array([[2, 2, 2], [1, 1, 1], [2, 1, 2]])\n    branch7=np.array([[2, 2, 1], [2, 1, 2], [1, 2, 1]])\n    branch8=np.array([[2, 1, 2], [2, 1, 1], [2, 1, 2]])\n    branch9=np.array([[1, 2, 1], [2, 1, 2], [2, 2, 1]])\n    br1=mh.morph.hitmiss(skel,branch1)\n    br2=mh.morph.hitmiss(skel,branch2)\n    br3=mh.morph.hitmiss(skel,branch3)\n    br4=mh.morph.hitmiss(skel,branch4)\n    br5=mh.morph.hitmiss(skel,branch5)\n    br6=mh.morph.hitmiss(skel,branch6)\n    br7=mh.morph.hitmiss(skel,branch7)\n    br8=mh.morph.hitmiss(skel,branch8)\n    br9=mh.morph.hitmiss(skel,branch9)\n    return br1+br2+br3+br4+br5+br6+br7+br8+br9\n\ndef endPoints(skel):\n    endpoint1=np.array([[0, 0, 0],\n                        [0, 1, 0],\n                        [2, 1, 2]])\n    \n    endpoint2=np.array([[0, 0, 0],\n                        [0, 1, 2],\n                        [0, 2, 1]])\n    \n    endpoint3=np.array([[0, 0, 2],\n                        [0, 1, 1],\n                        [0, 0, 2]])\n    \n    endpoint4=np.array([[0, 2, 1],\n                        [0, 1, 2],\n                        [0, 0, 0]])\n    \n    endpoint5=np.array([[2, 1, 2],\n                        [0, 1, 0],\n                        [0, 0, 0]])\n    \n    endpoint6=np.array([[1, 2, 0],\n                        [2, 1, 0],\n                        [0, 0, 0]])\n    \n    endpoint7=np.array([[2, 0, 0],\n                        [1, 1, 0],\n                        [2, 0, 0]])\n    \n    endpoint8=np.array([[0, 0, 0],\n                        [2, 1, 0],\n                        [1, 2, 0]])\n    \n    ep1=mh.morph.hitmiss(skel,endpoint1)\n    ep2=mh.morph.hitmiss(skel,endpoint2)\n    ep3=mh.morph.hitmiss(skel,endpoint3)\n    ep4=mh.morph.hitmiss(skel,endpoint4)\n    ep5=mh.morph.hitmiss(skel,endpoint5)\n    ep6=mh.morph.hitmiss(skel,endpoint6)\n    ep7=mh.morph.hitmiss(skel,endpoint7)\n    ep8=mh.morph.hitmiss(skel,endpoint8)\n    ep = ep1+ep2+ep3+ep4+ep5+ep6+ep7+ep8\n    return ep\n\ndef pruning(skeleton, size):\n    '''remove iteratively end points \"size\" \n       times from the skeleton\n    '''\n    for i in range(0, size):\n        endpoints = endPoints(skeleton)\n        endpoints = np.logical_not(endpoints)\n        skeleton = np.logical_and(skeleton,endpoints)\n    return skeleton",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "image = Image.new(\"RGBA\", (600,150), (255,255,255))\ndraw = ImageDraw.Draw(image)\nfontsize = 150\nfont = ImageFont.truetype(\"/usr/share/fonts/truetype/liberation/LiberationMono-Regular.ttf\", fontsize)\ntxt = 'A'\ndraw.text((30, 5), txt, (0,0,0), font=font)\nimg = image.resize((188,45), Image.ANTIALIAS)\nplt.imshow(img)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 61,
       "text": "<matplotlib.image.AxesImage at 0xaf026ac>"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAksAAACiCAYAAABRaYlhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WuIXOd9P/DvzJnbzuzMnpndnZm9SNpQWbGVOLGK2kBJ\nYIsqh1AQDgEnppglcd6k9EXaQqK8itQ3WdOWkrYulGKKoOBab2rLbSxik6xpElJDK0OIQ22KJc/u\n3HbmzP1+zpn/C/+fx+ecOTsrrXakvXw/sOxqNJdnx4P15Xd+z+/xDIfDIYiIiIjIlfdhL4CIiIjo\nIGNYIiIiIhqDYYmIiIhoDIYlIiIiojEYloiIiIjGYFgiIiIiGuO+wtLNmzfx6KOP4pFHHsHzzz+/\nX2siIiIiOjA8e52zZBgGPvnJT+LNN9/E0tISfud3fgcvvfQSHnvssf1eIxEREdFD49vrA99++22c\nPn0aKysrAICvfe1rePXVV21hyePx3PcCiYiIiB4UtxrSnsPS1tYWTpw4If+8vLyM//qv/xq53/e/\n/31cuXJlry9DNOLKlSv8TNG+4eeJ9hM/T4fbTkWePYelu60abWxsyA/O6uoqVldX9/qSRERERPtm\nY2MDGxsbu95vz2FpaWkJmUxG/jmTyWB5eXnkfqurq0zZREREdOA4izhXr151vd+ed8OdP38e77//\nPm7fvo1+v4+XX34Zly5dcl0I0X7iZ4r2Ez9PtJ/4eTqa9rwbDgBef/11fPvb34ZhGHjuuefwve99\nz/7kHo9roxQRERHRQbNTbrmvsLTXFyUiIiI6aHbKLZzgTURERDQGwxIRERHRGAxLRERERGMwLBER\nERGNwbBERERENAbDEhEREdEYDEtEREREYzAsEREREY3BsEREREQ0BsMSERER0RgMS0RERERjMCwR\nERERjeF72As4iu738GCPx7NPKyEiIqL7xbA0IXsNTOLEYwYmIiKig4FhaZ8Nh0P5tVcMTERERAcH\nw9IE3E/IYUAiIiI6WBiW9pG1msTQQ0REdDQwLO2T/QpK99scTkRERPuLYWkfDIdDmKYJ0zTRbrdR\nqVRQrVbR7Xbl3wnOICXC0XA4RDgcxszMDGZmZhAKheD3+xEIBORj2MdERET04DEs7YEIONYQYxgG\ndF1HrVbD7du38cEHH6BarcIwDBiGIYOOW9gxTRPD4RDz8/M4efIkTp06BVVVEQ6H4fP5oCiKLVRZ\nX5uIiIgmi2FpH1jDUr1ex507d/DOO+8gn89jMBhA13VbWLIGHVGRMk0Tp06dgq7rmJ6eht/vh8/n\nw9TUlO11GJKIiIgerF3D0je+8Q38x3/8B5LJJH71q18BADRNw1e/+lXcuXMHKysruH79OlRVnfhi\nDwJrT5EIL7quo9FooF6vI5vNYnNzE5lMBvl8HqZp2ipLAGzfrZfw/H4/EokEYrGYvP/09DQURRl5\nTSIiInowdj3u5Otf/zpu3rxpu219fR0XL17Ee++9hwsXLmB9fX1iCzyInHOUBoMBGo0GisUitra2\nsLW1hUwmg83NTWxtbSGXyyGfzyOXyyGbzSKbzSKXy9m+8vk8stksMpkM7ty5g3w+j3q9DsMw5GsS\nERHRg7drZekLX/gCbt++bbvtxo0beOuttwAAa2trWF1dPXaByUpUlra3t22VpWKxaLvsZg1ZHo8H\nXq/XNsRSVJKmpqYQjUaRSqVkWBKPZ1WJiIjowdpTz1KhUEAqlQIApFIpFAqFHe975coV+fPq6ipW\nV1f38pIHWr/fR7VaxdbWFrLZLGq1GgaDwdjHuO2K6/f7qNVqyOVySKVSqNVq6PV6mJqa2rE5nIiI\niPZmY2MDGxsbu97vvhu8d/tH3BqWjgrn7+sWlnRdl/cd9/6IviUA6PV6qFarUBQFCwsLMizpug5F\nUaAoCscIEBER7RNnEefq1auu99u1Z8lNKpVCPp8HAORyOSSTyb08zaFkDSji58FggEqlsmNlybkL\nzm1XnMfjkaErl8uhWCyiXq+j1+vtOHqAfUxERESTt6ewdOnSJVy7dg0AcO3aNTz11FP7uqiDyHlA\nrmma6Ha7aDQaqFar0DQNpVIJlUoF7Xbb1pgthlIqioJgMAi/32/b4SaYponBYIBOp4NGo4FyuYxC\noYByuYx2u20bbumcu0RERESTsWtYeuaZZ/B7v/d7+N///V+cOHEC//zP/4zLly/jjTfewJkzZ/CT\nn/wEly9ffhBrfWicgcTj8djCUqVScQ1LzoClKAoCgQACgYAMS9aGb9M00e/30el0UK/XoWka8vk8\nyuUyWq2WHF6507qIiIho/+3as/TSSy+53v7mm2/u+2IOMufkbMMw0Ov10Gw2RypL1mqSCEui5ygY\nDMrQ42wCF7f1+33U63WUy2Xk83nEYjFMT0+PPCd7loiIiCaPE7z3aDAYoFaryZlJtVoN/X5/xyAT\nDoeRSqWQSqUwGAxQKpWwvb2NbrcLwzBGqka9Xg+apmFrawszMzOIx+Njz5gjIiKiyWBY2qPBYCCb\nsXO5HKrVqgxLziDj8XgQDoexsLCARx55BL1eD4qioNVqyaNQrFUj4KOwVKlUsLm5idnZWSwtLck+\nKOsOOiIiIposhqV74Jza7awsOXfAWb+Hw2Gk02l88pOfRLvdRqvVQj6fH2kGF7rdLjRNQzAYRDqd\nRrPZhGmatqDE6hIREdHkMSztwhpgrGe5dTodVCoV5PN5eTSJqBIJiqLA5/PB7/dDVVWkUimcPHkS\nrVYLhUIBmUwG3W4XwEdTwK2X2QaDAZrNJkqlkuyFqtfr8Pl88su6RgYnIiKiyWBYGmOnoGQYBrrd\nrrwMJ8KStbIkAkwgEMDU1JQMSydOnECr1UImk0EikUCj0YCu62i327bXEufNGYZhC0tTU1MIhUK2\nsGRdK0MTERHR/mJYukvWsKTruqws5XI5FAoF2X8EfBxcvF4vAoEAwuEwVFVFMpnEyZMn0Wg0kEql\nkEgkoGka2u02vF6vfB3g48pSu91GqVRCtVpFo9FANBqFoigIhUKc5k1ERPQAMCztwK2BWtd1NJtN\nNBoNFItFaJqGRqOBTqcjg5L1cFyfz4doNIr5+XnMzs5ienpazlqamZlBKpVCpVJBs9mUR5k4D94V\nl/zEzji/3w+Px4NIJCIDlrgvAxMREdH+29ME7+PCOlAS+DgsbW9vo1AooFKpoNFooNvtyktw1sDj\n9/tlWJqbm0MkEoGiKPD7/YjFYkilUpifn5chyvl4a1gql8vY3NxEsViUzd5ijURERDQ5DEu7sE7h\n1nUdjUYDpVJpbGVJVHhEWEomk7awFAgEEIvFkEwmMTs7K6tEIiRZn8NZWSoWi3Kat3WNRERENBkM\nS7uwVnrEuAAxW8k5LgCwBxdRQUqn00gmk4hGo3J3XDQalZWlaDSKQCAgL6tZA9pwOJQzl0RYcu68\n4+U3IiKiyWHP0j2wzlbKZrO2sGS9ZCe+i94ka1gSvUmisjQ/P49YLIZgMAhFUWCappzmLZ6n2+2i\nUqkgEAhgYWFB7qAzTdNWhRKvzfBERES0fxiW7oGzslStVkfGBVj5/f6RsCS2/MdiMSiKgrm5OVtl\nyRqWBFFZ0nVd7qYbDAYyUFkrUuI7AxMREdH+YFhycM5WMk0Tg8FABiVN01AsFlEul9FqteT0bfFY\n0cDt8/kQj8cRj8eRSCQwMzMjt/uL+UvD4RCxWAyJRALz8/PyYN5mswnDMGx9S91uFx6PB9VqFeVy\nGdvb2wA+mgweDodtgUk0hzMwERER3T+GJQu3RmnTNNHr9dDpdFCr1WRQKZfLaLfbrrOVgsEgIpEI\nVFVFIpHA7OwsZmZmEAwG5WWzQCAARVFkWEomk7Jxu9PpyKNNRGDr9XowDEOGpWKxKMcITE1N2dbO\noERERLR/GJYcrE3TokLT7/fRarVQrVahaZoMS+KSmfVxXq8XoVAI0WhUVpVEWLIKBAKyd0mEpVqt\nJoddijWIS3O6rqPb7crAViwWEQ6HEQqFZO+SuCzHoERERLR/GJZ2IS6/5fN52afU6/VkSHI2V++0\nA078vZXH40EwGISqqlhcXES9Xke9Xpczl5yGw6E8YHdzcxPhcBjRaNQ2RoCIiIj2F0cHWDi37IvZ\nSvV6HblcDtlsFpVKBb1ez1ZJss5IsjZ1p1IpW1gCYGveHg6HCAQCUFUVS0tLSCaTtgGVbqxhSQyo\n5JwlIiKiyWFY2oXbDjhnZck5iFJM504mk4jFYrK3yBnEACAYDCIej2NxcVGGK+sxJk5ijIDbNG8i\nIiLaf7wM9/9Ze32sgabf79suwzlnK1mDimjcto4L2OkoE8F5JMrMzAwikQi63a7rGAGxnkKhgO3t\nbdRqNbRaLSiKAq/XO7YqRURERPeOYQmj4wKAjy+X9Xo9VKtVWVmq1+u2sGRlDUuismS9DGcNS87z\n4wBgbm4OqqoiEomg2Wyi3++PDKgUYWk4HKJUKqFaraLdbiMYDModdtbQx2ZvIiKi+8Ow5GANF2LL\nvghL+XwezWZTDoQEPg4k1vlJIiyJo0ys1R5ndUmEpWAwKCtL4XAYgUAAhmHIwCSIfqlOp2OrLIXD\nYRmUxCU/6/qIiIhob3btWcpkMvj93/99fOpTn8KnP/1p/O3f/i0AQNM0XLx4EWfOnMGTTz6JarU6\n8cVOgltztGEYaLVacou+ODC33W7LoGQNIIqiIBwOIx6PY25uDvF4HKqqYnp62nbmm/MSnHgeMchS\nPEc6nd7xgF2xvsFgINeYzWZRKpXknCbnsStsACciItq7XcOS3+/H3/zN3+DXv/41fvnLX+KFF17A\nb37zG6yvr+PixYt47733cOHCBayvrz+I9U6E9TIXAOi6LoPI9vY2KpUK6vW6DEsAbE3diqIgEokg\nkUhgbm5OTuyORCJyEKVgrfxYd9RZw1IqlcLc3JxtMrd1152YKt5sNlEul7G1tYXt7W20Wi1bQHL+\nXkRERHTvdg1L6XQaTzzxBABgenoajz32GLa2tnDjxg2sra0BANbW1vDKK69MdqUTZg0XbpWler2O\nTqfjWlny+XwIh8NyAGU8HsfMzIysLIn73k1lSVVVpFKpsZUlMaRSrHFra8tWWXL+TkRERLR399Sz\ndPv2bdy6dQuf+9znUCgUkEqlAACpVAqFQsH1MVeuXJE/r66uYnV1dc+LnSRrkNF1HdVqFdlsFltb\nW6hWq+j3+66PEYMoZ2ZmsLCwgIWFBTkuwGm3/iHrgMpKpYJCoTCyu80a6rrdLsrlMjY3N5FIJLC4\nuDjSrM5+JSIiIncbGxvY2NjY9X53HZaazSa+8pWv4Ic//KHcvSWM+0fZGpYOA4/HY5uttLW1hUql\nIi+/ifsAH/cC+f1+qKqKhYUFpNNpzMzM2AZRWo9QGScQCMiZS9vb23Lm0k6X1MSAykwmg3Q6LWcu\nWStZd/O6REREx5GziHP16lXX+93VUMrBYICvfOUrePbZZ/HUU08B+KialM/nAQC5XA7JZPI+l3xw\n6LqOWq22a2VJEDvgFhcXZWXJ5/Pdc9+QtbI0Pz9vm9Hk9lzjpnm7zXQiIiKie7drWBoOh3juuedw\n9uxZfPvb35a3X7p0CdeuXQMAXLt2TYaow8IaPKwH5rbbbdRqNWiahmKxiFKphGazCV3X5eNEX5C1\nMVtVVSSTSczNzWF6eho+n2/HwLLT7X6/XzaKz87OQlVVRKNRTE1Nwe/3j6xZ13W02215wK+maSiX\ny6jX67YjWdjsTUREtHe7Xob7+c9/jn/5l3/BZz7zGZw7dw4A8IMf/ACXL1/G008/jRdffBErKyu4\nfv36xBe7X9xCw3A4xGAwQLfbRb1el2FJ7DLTdd32uOFwCJ/PJ5u7xS42Z0UI2PkymPN2n88nH5tI\nJKCqKmKxGNrtNrrdrgxs4rEiLA0GA1QqFWiahlKpBI/Hg2g0auubEv1SnLtERER0b3YNS5///Od3\nPHvszTff3PcFPSjOy1XWypI1LJVKJRiGAcMw5OOGw6E8WiQQCNgqS/Pz8wgGg7aeJWG30OT3+6Eo\nCqampmRlKRaLoVar2SpbYs26rsvdeaKqVC6XEQqFEAgEZG+ZsyJFREREd48TvPFRmNB1HY1GA8Vi\nUR6Y2+12YRiGDCjWZmtRBRIhSRxTIo4cEc+7E+eUbfFdvE4wGJRnzHU6HZimKUcDuI0gsDZ7+/1+\nhEIhJBIJnhVHRER0n45tWHJusdd1HfV6XR6YW6lU5GG21h4jUWXz+XzyANxkMol4PC6PKXE7m82t\nouMcTmm9PRQKQVVVpNNpNBoNtFqtkeGW1t+j0+mgXC4jk8lgenoaiUSCPUpERET74K52wx01bseA\nGIYxUlnq9Xry+BAxQVsEFmtYSqVSUFUV4XBYXkqzuptLX9ZjSkRlSYQltwGVziZxMXMpk8mM7Iwj\nIiKivTt2Ycl58K24rd/vo1arycqSc1yAc0eZ3+9HLBZDOp1GKpWS4wL2ylmBEmFpYWEB8/PzMixZ\n12FdT7fblYM0t7e3Ua/XZVO4dfYSERER3ZtjexnO2fMzGAzkZbhsNotarSYHUbpdKrNO7U6n0ztO\n7b7XtQjWmUu5XG5k5pKoRIl19Xo9VCoVeDweGZZ6vZ7csWetdnFHHBER0d07VmHJ2ackbhNhyVpZ\nqtfrI5Ul63cRltLptAxLe6ksWRu9rQKBAFRVhWmachyBqCy5rafX68mm9GKxiFqthm63i2AwKIOS\ns6mcgYmIiGh3xyYsuQUS0zTRbrfR6XRQLBZRLpdRrVbRbDZlv5I1VHm9XlmliUajUFVVHpwbDodd\nz3G717WJn8X8JtM0kUgkEI/HEY/HYRgGBoMB+v2+7ZKi2NHX7XbRaDRQLpeRy+VgGAZmZmY4c4mI\niGiPjk1YAkbPaDMMA61WSx5aWy6XUalU5MRut11qgUAAwWAQsVhs17Bkfc27WZs1mCmKgnA4DJ/P\nh3g8Lqd69/t9NJtNOV9JNH2L6phoVBdhSVEUWQUTz82ZS0RERHfvWIUlYHQHXLvdRrlcHqksWZvA\nRcDweDzyeBNrZUlVVSiKAkVRRhqvd7rMttvaFEWR1aVEIiGrSyIotVotuT6v1wvTNGEYBkzTtIWl\ncDiMmZkZWzWJYYmIiOjuHbuwZCUOzM3n89ja2kKlUrH1KTmbwEWFJpVKYWlpCfF4HKFQCIqiyH4i\n58DIew0kbsEqFAphdnYWJ0+exGAwwGAwQLVa3XGyunXmUiwWQzKZHAlIDEpERER359iNDrASgyiz\n2exIWHI77Nbv98sdaktLS0gkEjIsOYdEioqP21yknb7E45zE8ScnT55EOp1GNBq1Dah0cpu5tFOw\nIiIiovGOdWXJMAxbZck5W8kZYERYWlpaslWWrPOPrI/dS/XGeslPPI+1slSpVLC1tTX2GJNOpwNN\n0+Dz+XDixAm0Wi1WloiIiPboyIclt8taYlBjs9lEpVJBsVjE9vY2Go2GnK0EQO6GE8+h6zqazSZK\npRI+/PBDeDwetNttxGIxOfsIcK9K3e1arWFJfGmahmw2i2w2a6sUiftap38DQL/fR6PRgNfrlU3r\ntVoNkUgEfr/fNuKAO+KIiIjGO9JhyW03m9hiL3aVibAkQoh1EKX47gxL29vbUBQFrVYLhUIBU1NT\ntvByP1Ul8d0allqtFqrVKqrVqlynYRg7rrPf78s5UaVSCZqmoVqtAoAMTM7XZWAiIiJyd6TDEmDv\nIXKbR2StLOm6LrfkO59DPK7RaKBUKqHX66FYLMqeJdM0bX1B9xOWrOsVh/z2ej30+3052sBaTXJ+\n7/f7GAwGaDabtspSIBCQu/msr8mZS0RERDs78mHJajgcwjAMeSktl8uhUqmg0+nIS3OAe9DxeDww\nTVOewdbr9WT4cJ7ZttNz3O0axePFlxhEqes6Wq0Wer2e62s4/zwcDuVohEwmA+DjyeOc5k1ERHR3\njnxYclZrxNDGQqFgC0uGYYwEHeclNdM05X1brZYcGXAvs5TulQhpojep3+/bduxZf0fnqAMAaLfb\nKJVKyGQyCIVCtplL1t+RiIiI3B3psORWORFhqVgsIpvN2ipLzkGU1seK27rdLjqdzsilr0lyzm1y\nW6d1l5t1bdaZS6qqYmFh4YGsmYiI6Kg4smHJrWIijgSp1+vI5/PIZrMj4wLE/dyeb6ceoUnbaRq4\nc01uFaJutwtN07C1tYV0Oo1GoyGDoXhOVpaIiIh2diTDklujtHVH27ip3bs937jbJs3tyBLr7W7r\nEjOXFEXByZMnZViyXkJ0ex4iIiL6yNgJ3t1uF5/73OfwxBNP4OzZs/je974HANA0DRcvXsSZM2fw\n5JNPym3pB5G1cuKsLLkNotyJtYH7YV7Gcq7DGZqs34GP/huKQZbb29tyPILo0XJWlniJjoiIyG5s\nZSkUCuGnP/0pwuEwdF3H5z//efzsZz/DjRs3cPHiRXznO9/B888/j/X1dayvrz+oNY/l9o+9YRiy\n12h7e1tup6/X6+h2uzAMw7Xh2ev1wu/3Y2pqCuFwGH6/33UC9qQDhttrDYdD9Ho9tNtttNtt25wn\n6+MMw5C752q1mjw0WFVVRCIRRCIR23NzjAAREZHdrpfhxEyefr8PwzAQj8dx48YNvPXWWwCAtbU1\nrK6uHpiwBIz2+JimiXa7jWq1KsOSpmmo1+uyyuJ8vNfrhdfrRTAYRDwex9zcHCKRiDzz7WESwahW\nq6FUKskhm85J3mInXa/Xk5cfy+Uy8vk8hsMhFEWRYYljBIiIiNztGpZM08Rv//Zv4//+7//wrW99\nC5/61KdQKBSQSqUAAKlUCoVCYeILvVfWwCS2/FvDkqgsOXeWWS9NWcPS4uIi4vG4DFEPI1A4xwgU\nCgV5WdE6+kDcV1SWdF3HcDhEvV6XlaVQKITp6Wnb8zMoERERjdo1LHm9Xrzzzjuo1Wr44he/iJ/+\n9Ke2v99tN9WVK1fkz6urq1hdXd3zYvfKMAxbn1KlUrENdnSuX1EUhMNhRCIRpFIpnD59Go8++iiS\nyeSBqCyJUHfnzh0EAgEYhoFKpYJWq4Vms+k6LwqAHFD54YcfYmpqCqqqMiAREdGxtbGxgY2NjV3v\nd9e74WZmZvCHf/iH+O///m+kUink83mk02nkcjkkk8kdH2cNSw+LOKbEuQPOuX1eVG3E5SlRUTpz\n5gyeeOIJLC0tyed0zj6aJLfjVwAgFovJuVFi/a1Wa6RxW1TNOp2OPARY/G5WDE1ERHScOIs4V69e\ndb3f2LBUKpXg8/mgqio6nQ7eeOMNfP/738elS5dw7do1fPe738W1a9fw1FNP7evi99tulSVRKRJh\nQVSWEokElpaWcPr0aTzxxBM4deqU6+G8D4LzEhsABINBNJtNFItF2ey909BKADIshUIhLC4ujgQr\nVpmIiIhGjQ1LuVwOa2tr8pDYZ599FhcuXMC5c+fw9NNP48UXX8TKygquX7/+oNZ7TwzDkNWWarWK\nYrGIQqEgG7sF55Z7v9+PWCyGdDqNhYUFqKoKv98/cv8HFTDcdrl5PB5MTU0hkUjgxIkTaLVaqFQq\n8Hq9sn/JGex6vR5qtRr8fj9KpRKq1SqazSaCwSAURYHP57O9JoMTERHRLmHp8ccfx//8z/+M3J5I\nJPDmm29ObFF75QwxhmGg3++j3W6jUqnYwpKYreQWKnw+H2ZmZrCwsIDFxUWoqopAIGC7v7Mi8yB+\nN+fW/lAohEQigeXlZTml29qn5JzF1O12UavVYBiGLSwNh0MEg0FbWBLPwcBERETH3ZGZ4O12mUrX\ndfR6vZHK0mAwwGAwsIUP63P4fD7EYjEZlmZmZmyVJdM0H+BvZg8+Xq9XrjMUCmF2dhbdbhf5fB7T\n09O25nNnEOz1evIQYGtYElUltzPxGJiIiOi4OxJhyW1GkGmaaDab0DQNuVwOmqah3W7b5iqJYCFC\niDgCRDR3z8/PY35+HtFoVFZdRI+T2/TrSXEbawAAgUAAsVgMg8EAc3NziMfjchyAGBkgHm99X3Rd\nl4Fpc3MTg8FA/t7i/XDuqCMiIjqujkRYAkarKIZhyObnbDYLTdPQarVsvUrOS2iKoiAQCGB6ehqq\nqiKZTMqwZJ3ebQ1L4nkm/btZqzzOsKQoCubm5qCqKqLRqKyoibBkXaO4PNlsNlEqlZDJZODz+UZm\nLhEREdFHjkRYcvYRAR+Hpe3tbVtY0nVdVoes1Rbgo8tvwWDQVlmam5tDIBCw9fM4t+VPmvNymDUs\niZ178/PzUFUV09PT8ggX8d6IQZqiUV80vYvKUjQaxfz8vHw9VpaIiIg+diTCkltwEfOHxLgATdPk\nuADr40QoUBQF09PTmJ2dRTqdhqqq8jw4RVEAjM47epBBwq23Stzu8/nkqIPFxUU54bvVau34HJ1O\nB5qmIZPJYH5+Xo4REPezficiIjrOjkRYAkaHRFpnK21ubo7MVnLy+XyywpJOpxGPx0cOz32Ql97c\nWHuXBNF3ZZ0LJZraxzVot9ttlEol+P1+nDx5Eu12e+QQYYFN3kREdJwdibDk9g+5tbLkDEtujxGV\npWQyiVQqhXg8jqmpKfj9/pEt+AchOJimaTvDTsxcWlxclAfsOgOkc0BluVyWYwRarRZM04SiKLbL\njDxgl4iIjruHe8jZfXKGGNM05SyhcrksD8yt1Wpot9u2Ph5xtIlzXEAqlcLCwgJmZmbkbKWDyDql\nG4CcubS4uIi5uTmEw2FbJUr8vuKxuq6j0+mgVquhUqmgXC6jVCqhXq/LUOm87Pcg5kkREREdNIe2\nsuRW8RDnnzUaDZRKJZTLZWiahmq1isFgYNsd5gxaIiyl02k5iNI6WwkYDSgPk3MNoVAI8Xgcy8vL\nyGaziEQittEI1u/ARzOoxCiFSqUCTdOwvb0NVVXlwEvnY4iIiI6jQxuWgNFxAdbKUqlUgqZp0DQN\ntVrNFi7cDqa9m8qS2+DLB82tfwr4OCwNh0PMzc0hEom4NoWLsCcGc4rp5qKy5PV6ZVASDtLlRyIi\nogftUIclJ2efknMHnPV4EhE2AoEA/H4/VFVFPB7H7OysbO52Hv9xkCtLPp8PkUgEw+EQiUQCiUQC\ns7Oz8rLcpbYrAAAKFklEQVRav98fOWBXjEwQzd537tyB1+vlzCUiIiKLIxWWTNNEo9FAoVCwhSVr\nI7Qg5g8Fg0GEw2HE43EZMFRVlYfLCgchIAnOZu3hcAi/3y8Dnvg9ZmdnZSgUwzidM6KGwyHa7Ta2\nt7dx584dTE9P22YuERERHXdHKiyN2wFnraqInWQiLEWj0ZHKkniM+H7Qenecu9V8Pp+sLlkrS/1+\nXx4aLIgBlYKoLH344YdIJpNyjAAREREdgbAkZiqZpin7bwqFAnK5HKrV6sggSuDjXidFURCNRmWf\nkqqqtj4lt8BwkCpMglv/0uzsLE6cOCF3vWmaZutfsv5u3W4XmqYhEAjg1KlTqNVq6Ha78qw8zlwi\nIqLj7FCGJWvDsei96ff7rmHJWlVxhh9rU7fYAbdTWDoovUpWzktxYn1TU1MyLInDhK29Ws6wJMKU\nYRjY3t5GrVZDp9NBIBCQR6o4X/egvRdERESTcujCkjPAiINhRViqVqsyLInGZvE4t4NzrZUl6w64\ngzaEcifOM/GAj2cuLS8vy/PfrNO8nTvkOp2O7Sy9er2OTqcjj1JxPmbcZHAiIqKj5lCFJbfZSuLy\nW7VaRbFYRLVaRavVkg3NPp9PVkac/7iLxu50Oo10Oo2ZmRl5vMlh6dlxCyxijMDi4iKy2aw8567f\n749UlcRzDIdDDAYDtFotaJqGQqGAeDzuOnOJQYmIiI6TQxWWgNFLTyIsaZomw5IYtqgoiq0yIh7r\n9Xrh8XgQiURcwxIwetbcQQwH1lDnNnPJNE0kk0moqirnLhmGAcMwRhreTdOErutotVqoVCrI5/Pw\neDyYmpqyvSaDEhERHTeHLiwBoz1LbmFJ13V4vV5ZWRKBAPjo8puiKLbK0sLCAqanpw/0ESc7cR7u\nGwwG5fiD+fl5WVkyTRODwQD9fl82biuKgn6/D8MwRipL4nBeIiKi4+zQhSVnMBBDFOfm5tDv9zEY\nDDA9PY3HHnsMPp8Pfr8fXq/XNoRR7PJKJBL49Kc/jXQ6jXA4jEAgYNv5dRgqKG79Q4qiyNC3vLyM\nz372s/B4POh0OtB1Hbquy9EJXq8Xuq7L42B+67d+C6dPn8bS0hISiYRtmrfzvSciIjoODlVYcgYD\n8Q9+JBLB/Pw8AoEAYrEYTp06hUajISsnbpfhvF4vwuEwFhYWZFhy9jcdhgZvYHR9Xq9XhsSlpSV4\nvV6k02kMBgN5GU48znoZzjAMzM7OyhlNsVgM4XBY3tf53hMRER0HnuFddDIbhoHz589jeXkZr732\nGjRNw1e/+lXcuXMHKysruH79OlRVHX3yB9AoPRwObZURUUFy660Rt4kv0dM0rgn8sBOVJF3XZSgS\nFTa38/JEP5e4hOn3+0eOfSEiIjqKdsotXpf7jvjhD3+Is2fPyn9c19fXcfHiRbz33nu4cOEC1tfX\n93e190gEH7/fD7/fj0AggGAwKOcEiS/nbSIkHbWAZCWqayL4iPfB+V6I28R7eBzeGyIioruxa1ja\n3NzEj370I3zzm9+UaevGjRtYW1sDAKytreGVV16Z7CrHcFZBRCCw/iy+xG0iDPh8PllJOaqs74/1\nPdrpvbH+LHq7iIiIjrNdr6/86Z/+Kf7yL/8S9Xpd3lYoFJBKpQAAqVQKhUJhx8dfuXJF/ry6uorV\n1dW9r3YH9/sP+mGZqbQXoupGREREdhsbG9jY2Nj1fmN7lv793/8dr7/+Ol544QVsbGzgr//6r/Ha\na68hHo+jUqnI+yUSCWiaNvrkh2S442Fp5CYiIqLJ2Sm3jK0s/eIXv8CNGzfwox/9CN1uF/V6Hc8+\n+yxSqRTy+TzS6TRyuRySyeTEFn4v9hrMjkNI2st7cxzeFyIiot3c1W44AHjrrbfwV3/1V3jttdfw\nne98B7Ozs/jud7+L9fV1VKtV1ybvB1lZup/XOeqhgO8NERHR7u5rN5z1SQDg8uXLeOONN3DmzBn8\n5Cc/weXLl/dnlQ/BcQgDe/0dj8N7Q0REtJu7rizt6ckPSc8SERER0b5UloiIiIiOG4YlIiIiojEY\nloiIiIjGYFgiIiIiGmPiYeluJmMS3Qt+pmg/8fNE+4mfp6OJYYkOHX6maD/x80T7iZ+no4mX4YiI\niIjGYFgiIiIiGmPiQymJiIiIDot7Pkh3Ei9IREREdJjwMhwRERHRGAxLRERERGNMNCzdvHkTjz76\nKB555BE8//zzk3wpOqJWVlbwmc98BufOncPv/u7vAgA0TcPFixdx5swZPPnkk6hWqw95lXRQfeMb\n30AqlcLjjz8ubxv3+fnBD36ARx55BI8++ih+/OMfP4wl0wHn9pm6cuUKlpeXce7cOZw7dw6vv/66\n/Dt+po6GiYUlwzDwJ3/yJ7h58ybeffddvPTSS/jNb34zqZejI8rj8WBjYwO3bt3C22+/DQBYX1/H\nxYsX8d577+HChQtYX19/yKukg+rrX/86bt68abttp8/Pu+++i5dffhnvvvsubt68iT/+4z+GaZoP\nY9l0gLl9pjweD/7sz/4Mt27dwq1bt/ClL30JAD9TR8nEwtLbb7+N06dPY2VlBX6/H1/72tfw6quv\nTurl6AhzbhS4ceMG1tbWAABra2t45ZVXHsay6BD4whe+gHg8brttp8/Pq6++imeeeQZ+vx8rKys4\nffq0DOhEgttnCnDf0MTP1NExsbC0tbWFEydOyD8vLy9ja2trUi9HR5TH48Ef/MEf4Pz58/inf/on\nAEChUEAqlQIApFIpFAqFh7lEOmR2+vxks1ksLy/L+/H/WXQv/u7v/g6f/exn8dxzz8lLu/xMHR0T\nC0ucsUT74ec//zlu3bqF119/HS+88AL+8z//0/b3Ho+HnzXas90+P/xs0d341re+hQ8++ADvvPMO\nFhYW8Od//uc73pefqcNpYmFpaWkJmUxG/jmTydgSNtHdWFhYAADMz8/jy1/+Mt5++22kUink83kA\nQC6XQzKZfJhLpENmp8+P8/9Zm5ubWFpaeihrpMMlmUzK4P3Nb35TXmrjZ+romFhYOn/+PN5//33c\nvn0b/X4fL7/8Mi5dujSpl6MjqN1uo9FoAABarRZ+/OMf4/HHH8elS5dw7do1AMC1a9fw1FNPPcxl\n0iGz0+fn0qVL+Nd//Vf0+3188MEHeP/99+UOTKJxcrmc/Pnf/u3f5E45fqaOjolN8Pb5fPj7v/97\nfPGLX4RhGHjuuefw2GOPTerl6AgqFAr48pe/DADQdR1/9Ed/hCeffBLnz5/H008/jRdffBErKyu4\nfv36Q14pHVTPPPMM3nrrLZRKJZw4cQJ/8Rd/gcuXL7t+fs6ePYunn34aZ8+ehc/nwz/8wz/wkgmN\ncH6mrl69io2NDbzzzjvweDz4xCc+gX/8x38EwM/UUTLRs+GIiIiIDjtO8CYiIiIag2GJiIiIaAyG\nJSIiIqIxGJaIiIiIxmBYIiIiIhqDYYmIiIhojP8HYFBoXrSC+ysAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "img = np.array(img)\nim = img[:,0:50,0]\nim = im < 128\nskel = mh.thin(im)\n# To fix the algorithm, prune the skeleton 1,2,3...\nskel = pruning(skel,3)\nbp = branchedPoints(skel)\nep = endPoints(skel)\nedge = skel-bp\n\nedge_lab,n = mh.label(edge)\nbp_lab,_ = mh.label(bp)\nep_lab,_ = mh.label(ep)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "figsize(10,20)\nsubplot(141,frameon = False, xticks = [], yticks = [])\ntitle('pruning 0')\nimshow(mh.thin(im),interpolation = 'nearest')\nsubplot(142,frameon = False, xticks = [], yticks = [])\ntitle('pruning 1')\nimshow(pruning(mh.thin(im),1),interpolation = 'nearest')\nsubplot(143,frameon = False, xticks = [], yticks = [])\ntitle('pruning 2')\nimshow(pruning(mh.thin(im),2),interpolation = 'nearest')\nsubplot(144,frameon = False, xticks = [], yticks = [])\ntitle('pruning 3')\nimshow(pruning(mh.thin(im),3),interpolation = 'nearest')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 63,
       "text": "<matplotlib.image.AxesImage at 0xcf0066c>"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAACKCAYAAACuA5DbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC09JREFUeJzt3V9o1fUfx/HXmc7a3H82144aA81pG+1mF8s1Kggl9SKM\nYDSItIjA6CJMwcX+sBVtCNGFXgST1sVSgjDwX4MNE2N0E0nlIEjcMklQ53LYVsr3d9Fvx22dzfPn\n+933+32f5+PKwzmdfc45rz6+fH+/37OI4ziOAAAADMvyewEAAABeo/AAAADzKDwAAMA8Cg8AADCP\nwgMAAMyj8AAAAPMoPPPk5+fr8uXLfi8DIUeO4AZyhHSRofsoPPPcvn1blZWVnjz3Rx99pIqKChUW\nFuq1117T33//7cnPgf+8ytFPP/2krVu3qqysTFlZ/O9rnVc56uvrU11dnQoLC7V27Vrt379f9+7d\nc/3nwH9eZejo0aPauHGjCgsLVVpaqp07d+rq1auu/xw3mdox79696/cSFvT111+ru7tbQ0NDGh0d\n1aVLl9TW1ub3shBHkHO0YsUKNTU1qbe31++l4AGCnKO//vpLH3/8sW7cuKHvvvtOg4ODOnjwoN/L\nwjxBzlBDQ4POnTuniYkJjY6OKjc3V++8847fy1pU4AtPZWWlPvzwQ1VXV6ukpES7d+/W9PS0JOns\n2bNas2aNenp6VFFRod27d6uvr0+NjY1zniMrK0uXLl2SJL366qvas2ePduzYoYKCAtXX18fuS/ax\nAwMDqqqqUlFRkfbs2aOnn356wb+I+vr69Prrr2vTpk0qKipSa2urPv30UzffKizCSo42bNigXbt2\n6fHHH3f1/UFirOTozTffVENDg5YvX65oNKrm5mZ9++23rr5XiM9KhtauXatVq1ZJkhzH0bJly1RR\nUeHeG+WBwBceServ79fAwIB+/fVX/fLLL+rq6ordd+3aNY2Pj2tsbEyffPKJEvlNGceOHVN7e7vG\nx8e1fv16tbS0JP3Y69ev66WXXlJ3d7du3rypqqoqDQ8PKxKJxH2eixcvqra2Nnb7iSeeiK0dS8NC\njuA/izn65ptvVFNTk9BjkT4rGTp//ryKiopUUFCgsbExdXd3J/EuLL3AF55IJKK33npLq1evVnFx\nsVpaWvT555/H7s/KylJHR4eys7P18MMPJ/R8O3fuVF1dnZYtW6bm5mb98MMPST/21KlTqqmp0Qsv\nvKCsrCy9/fbbeuSRRxb8uZOTkyosLIzdLigokPTv8VV4z0qO4C+LOTpy5Ii+//577d27N6HHIz2W\nMvTUU0/p1q1bunLlirKzs/Xuu+8m8U4svcAXHunf0dmMRx99dM6JUWVlZVqxYkVSz1deXh77c05O\njiYnJ5N+7NWrV7VmzZo5j51/e7a8vDz9+eefsdsTExOS/j2DHkvDQo7gP0s5On78uA4cOKDTp0+r\npKQkqXUjdZYyJEnRaFSdnZ367LPPkln2kgtF4RkbG5vz52g0Grs9f9y2cuVK3blzJ3b7jz/+8GRN\n0WhUV65cid12HGfO7fmqq6vntO4LFy6ovLxcxcXFnqwP/2UhR/CflRydOXNGb7zxhk6cOKHq6mpP\n1oX4rGRotn/++Ue5ubleLM01gS88juPo8OHD+v3333Xz5k29//77ampqWvDxtbW1+vnnn3XhwgVN\nTU2pvb39P8+XzM9eyLZt2/Tjjz/qq6++0t27d3Xo0KFFg/jKK6+ot7dXIyMjGh8fV2dnp3bt2pXw\nWpAeKzmSpKmpqdhXGkxPT8dOeIT3rORoaGhIzc3N+vLLL1VXV5fwGpA+Kxnq7+/Xb7/9JkkaHR1V\nS0uLXnzxxYTX4ofAF55IJKKXX35ZW7Zs0bp16/TYY4/pvffem3P/bBs2bFBra6uee+45VVVVqbGx\ncc5jIpHIf/6b+fcn8tjS0lJ98cUX2rdvn0pLSzUyMqK6ujo99NBDcV/H1q1btW/fPj377LOqrKzU\nunXr1NHRkeS7gVRZydHly5eVm5urmpoaRSIR5eTkaNOmTUm+G0iVlRx1dXXp9u3bev7555Wfn6/8\n/Hxt3749yXcDqbCSoYsXL2rz5s3Ky8vTM888oyeffFI9PT1JvhtLzAm4yspKZ3Bw0O9lPNC9e/ec\naDTqnD171u+lIA5yBDeQI6SLDPkn8BOeIBsYGNCtW7c0PT2tDz74QJJUX1/v86oQNuQIbiBHSJf1\nDFF40jA8PKz169errKxMJ0+e1PHjxxcc/wELIUdwAzlCuqxnKOI4SZzxBAAAEEJMeAAAgHnLF7sz\nEmlfomXAb47T7tlzk6PMQIbgBnKEdC2UISY8AADAPAoPAAAwj8IDAADMo/AAAADzKDwAAMA8Cg8A\nADCPwgMAAMyj8AAAAPMoPAAAwDwKDwAAMI/CAwAAzKPwAAAA8yg8AADAPAoPAAAwj8IDAADMo/AA\nAADzKDwAAMA8Cg8AADCPwgMAAMyj8AAAAPMyqvC0qUNt6vB7GQAAYIllVOEBAACZabnfCwia2ROg\nDrX5uBIEEfkAgHBiwgMAAMyj8AAAAPM4pPV/nMwMN8zPEYe9EM9MTuLlY7H7AKSOCQ8AADAvIyc8\ni01z+FcVEkWOkC4my0hVotlhL7qPCQ8AADCPwgMAAMwzf0gr3thvsRMF+Z4VxJNqjsgQZrAXwQ3s\nRaljwgMAAMwzP+GZbbGGO3Pf7PZMM8Z8D8rC/Bzxr3TEw16EdLEXJY8JDwAAMI/CAwAAzMuoQ1qJ\nmD3q4zsykKp4hyWAZLAXwQ3sRfcx4QEAAOaZnfDQZpEuMgQ3kCOkiwy5gwkPAAAwj8IDAADMM3tI\na4Yb3zfA9xfAzRyRoczEXgQ3sBeljgkPAAAwz9SEx+0Tu7icLzO5+XlzaXFmYi+CG9iL3MWEBwAA\nmGdqwjNbph2bhPvIENxAjpAuMuQOJjwAAMA8Cg8AADDP7CEtr2Tq5XxwD5cWww3sRUhXpu1FTHgA\nAIB5JiY8Xl9ix+V89i3F58qlxfaxFyFd7EXeYcIDAADMo/AAAADzTBzSmpEJJ13Be+QI6SJDcAM5\nchcTHgAAYJ6pCc9SyrTL+eANLi1GutiL4IZM2IuY8AAAAPMoPAAAwLzQHtLy6/sDMvX7C6zy43Pk\nu1RsYS+CG9iLvMeEBwAAmBfaCc9sfp9klQkne1nn92fHiac2+P3ZsReFn9+fneW9iAkPAAAwj8ID\nAADMM3FIyw+ZdrIXvMGJp0gXexHckAl7ERMeAABgHoUHAACYR+EBAADmUXgAAIB5FB4AAGAehQcA\nAJgXcRzHWfDOSPsSLiUxYblkLmzfUOk47Z49NzlKDRm6jwyljhzdR45SYyVDTHgAAIB5ofjiwXgN\nOIiNMwxNPZORI6SLDMEN5MgfTHgAAIB5FB4AAGBeoA9phWXsF8/M2sOyXsvCmqPZ6w7Dei0La4Yk\n9qIgCWuOrOxFTHgAAIB5gZ7wzAhLo+S3Fgdb2HJEhoInbBmSyFEQhS1HVjLEhAcAAJhH4QEAAOZR\neAAAgHkUHgAAYF4gT1q2coKUZOdyvrCxlCGJS4v9YilH7EX+sJQhKdx7ERMeAABgHoUHAACYF8hD\nWjPCODKbYe37C8IsrDniu1SCI6wZktiLgiSsObKyFzHhAQAA5lF4AACAeRQeAABgnu/n8Cx2PDDM\nl7/FE9bflBt0YT6mnCwy5L1MyBM5QrrC+PczEx4AAGAehQcAAJjn+yGtGWEai6Uq3uWhYRwLBk0m\nvHfzX2MmHHbxi9U8WX1dQWftfQ/z62HCAwAAzPN9whPmtpiIeK/P+muG98gQACSHCQ8AADCPwgMA\nAMyj8AAAAPMoPAAAwDwKDwAAMI/CAwAAzKPwAAAA8yg8AADAPAoPAAAwj8IDAADMo/AAAADzKDwA\nAMA8Cg8AADCPwgMAAMyLOI7j+L0IAAAALzHhAQAA5lF4AACAeRQeAABgHoUHAACYR+EBAADmUXgA\nAIB5/wM/2Wni1FZT4AAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "subplot(221,frameon = False, xticks = [], yticks = [])\ntitle('skeleton')\nimshow(skel,interpolation = 'nearest')\nsubplot(223,frameon = False, xticks = [], yticks = [])\ntitle('skel-bp -> label')\nimshow(edge_lab,interpolation = 'nearest')\nsubplot(224,frameon = False, xticks = [], yticks = [])\ntitle('skel -> end points')\nimshow(ep_lab,interpolation = 'nearest')\nsubplot(222,frameon = False, xticks = [], yticks = [])\ntitle('branched points (bp)')\nimshow(bp_lab,interpolation = 'nearest')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 64,
       "text": "<matplotlib.image.AxesImage at 0xd0ed9ec>"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAANiCAYAAABl/BOpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuQleV9wPHfYVHBuNwRYZcVVrQkGqMNirE6kIrXYKWo\nEUTrJRltGqOYjAENCF4Sq5GO1TpNTE1pRMWSMaMSHYlaHGK9ZkKMKKIioAGV5SKguCz49g/HU1Z2\ncQVh4cfnM8MMey7vPuc97sN3n332tVQURREAAIm1ae0BAABsa4IHAEhP8AAA6QkeACA9wQMApCd4\nAID0BE8iEydOjLPPPnu7PQ/Y8fXp0yceffTR1h5GnHvuuTF+/PjP5ViTJ0+Oo48++nM51sZmzZoV\n/fv3/9yP25T6+vo48MAD4+23346IrTs/DzzwQIwYMeLzHF5KgieRUqm0zZ83ePDguP3227fo8wDb\nX6lU2uK5IeM4Nufoo4+OuXPntuixM2fOjN69e2/x57rtttti0KBB0aNHj4jYuvNz8sknx5w5c+LP\nf/7zFo9nVyB4EtnSa0h+luft6BMWsOXWr1/f2kPYZfz85z/fZGV9a64DPHLkyLjtttu2dlipCZ6d\n1PXXXx/V1dXRoUOH6N+/fzz22GONYqShoSFGjhwZp512WjQ0NMTixYvj1FNPjb333jtqa2vjlltu\nafbYTz31VBx55JHRuXPnOOSQQ+Lxxx+PiIgf/ehHMWvWrLjooouisrIyLr744oiI+N///d847LDD\nolOnTnH44YfHk08+WT7W4MGD48orr4yjjjoqOnToEMcff3wsW7ZsG50VoCnPPPNMHHjggdGlS5c4\n//zzo76+PiI+WqWorq6OG264IXr27Bnf+ta3YuXKlTF06NDYe++9o0uXLnHyySfHX/7yl/KxPu1r\n+ve//315/qipqYlf/epX5fuWL18eQ4cOjQ4dOsQRRxwR8+fPL983d+7cOPbYY6Nr167Rv3//mDZt\nWvm+ZcuWxd/93d9Fx44dY+DAgfHaa681+1oXLFgQbdq0iV/84hdRVVUVvXr1ikmTJpXvr6+vj9Gj\nR0dVVVVUVVXFpZdeGuvWrSufj41Xbfr06ROTJk2Kr3zlK9GpU6cYMWJE1NfXx3vvvRcnnnhiLF68\nOCorK6NDhw7x1ltvxTPPPBMDBgyIjh07xj777BM/+MEPmhzjokWLYv78+TFw4MBGt9fV1cVxxx0X\nHTp0iMGDB8eiRYvK97Vp0yZuueWW2G+//aJ79+7xwx/+sFEgDR48OH772982e16IiIKdzty5c4ve\nvXsXS5YsKYqiKBYuXFi89tprxcSJE4uzzjqrWLt2bXHSSScV5513XvHhhx8WGzZsKP76r/+6uOaa\na4qGhoZi/vz5RW1tbfHwww8XRVEUEyZMKM4666yiKIrizTffLLp27Vo89NBDRVEUxe9+97uia9eu\nRV1dXVEURTF48ODi9ttvL49l2bJlRadOnYopU6YUGzZsKO6+++6ic+fOxfLly4uiKIpBgwYV/fr1\nK1555ZVi7dq1xeDBg4uxY8dut3MFu7p99923+PKXv1y8+eabxfLly4u/+Zu/KcaNG1cURVH8z//8\nT9G2bdti7Nixxbp164q1a9cWy5YtK+69995i7dq1xerVq4vTTz+9GDZsWPl4m/uaXrBgQVFZWVlM\nnTq1WL9+fbFs2bJi9uzZRVEUxTnnnFN07dq1ePbZZ4v169cXo0aNKkaMGFEURVGsWbOmqK6uLiZP\nnlxs2LCh+OMf/1h069atePHFF4uiKIozzjijOOOMM4r333+/eOGFF4qqqqri6KOPbvL1vv7660Wp\nVCrOPPPM4v333y/+/Oc/F927dy8eeeSRoiiKYvz48cXXvva1YunSpcXSpUuLI488shg/fnz5fFRX\nV5eP1adPn2LgwIHFkiVLiuXLlxdf/OIXi5/97GdFURTFzJkzGz22KIriiCOOKKZMmVIURVG89957\nxVNPPdXkGKdPn14ceOCBjW4755xzisrKymLWrFlFfX19cckllxRHHXVU+f5SqVT87d/+bbFixYpi\n0aJFxQEHHFD8x3/8R/n+ZcuWFaVSqVi9enWTn5OisMKzE6qoqIj6+vqYM2dONDQ0RE1NTdTW1kZE\nxKpVq+L444+P/fffP375y19GqVSKZ599Nurq6mLcuHHRtm3b6Nu3b3z729+OqVOnbnLsKVOmxEkn\nnRQnnHBCREQMGTIkBgwY0Og7h2Kj7yp++9vfxl/91V/FqFGjok2bNjFixIjo379/3H///RHx0Y/A\nzjvvvOjXr1+0a9cuvvnNb8bs2bO35ekBNlIqleKiiy6Kqqqq6Ny5c/zoRz+Ku+++u3x/mzZt4qqr\nrorddtst2rVrF126dIm///u/j3bt2sVee+0VV1xxRXmV9+PjNfc1fdddd8Wxxx4bZ5xxRlRUVESX\nLl3iK1/5Svl5w4cPjwEDBkRFRUWMGjWq/Lzp06dH375945xzzok2bdrEIYccEsOHD49p06bFhg0b\n4t57742rr7462rdvHwceeGCcc845n/rjnwkTJkT79u3joIMOivPOO6/8mu+888648soro1u3btGt\nW7eYMGFC3HHHHc0e5+KLL4599tknOnfuHCeffHJ5zE19/t133z1eeeWVqKuriz333HOTFZyPrVy5\nMiorKze5fejQoXHUUUfF7rvvHj/+8Y/jySefbLS6NmbMmOjUqVP07t07Ro8e3eh9/Ph4K1eu3Ox5\n2ZUJnp1Qv3794qabboqJEydGjx49YuTIkbFkyZIoiiKeeuqpeOGFF2LMmDHlxy9cuDAWL14cnTt3\nLv+57rrr4p133tnk2AsXLoxp06Y1euwTTzwRb731VvkxG//obPHixVFTU9PoGPvuu28sXry4/PE+\n++xT/nv79u1jzZo1n8t5AFpm4x/T1NTUNPr67N69e+y+++7lj99///248MILo0+fPtGxY8cYNGhQ\nvPvuu43+gW/ua/qNN94of/PVlI836H7yeQsXLoynn3660bxz1113xdtvvx11dXWxfv36TV7DZ33N\nS5YsiYiIJUuWxL777tvs+fikzzJ/3X777TFv3rz44he/GIcffnizP2Lq3LlzrF69utFtpVIpqqur\nyx9/4QtfiC5dujQa2+bex4+P16lTp2bHt6sTPDupkSNHxqxZs2LhwoVRKpVizJgxUSqV4rjjjoux\nY8fGMcccUw6ampqa6Nu3b6xYsaL8Z9WqVTF9+vRNjltTUxNnn312o8euXr06fvjDH0bEppuWq6qq\nYuHChY1uW7hwYVRVVW2jVw58VhvvBVm0aFH06tWr/PEnv6YnTZoU8+bNi2eeeSbefffdePzxx6Mo\nihZtqK2pqdns/prNPW/QoEGbzDu33nprdOvWLdq2bbvJa/g0zb3mXr16xYIFC5q877No6hc4+vXr\nF3fddVcsXbo0xowZE6eddlqsXbt2k8cdfPDB8frrr8eHH37Y6PY33nij/Pc1a9bE8uXLG43tk69p\n43n2pZdeij59+sRee+31mV/LrkLw7ITmzZsXjz32WNTX18cee+wR7dq1i4qKivL9l112WZx55plx\nzDHHxLJly+Kwww6LysrKuOGGG2Lt2rWxYcOGeOGFF+K5557b5NhnnXVWPPDAAzFjxozYsGFDfPDB\nBzFz5szysmqPHj0aTWgnnXRSzJs3L+6+++5Yv3593HPPPTF37twYOnRo+TEtmSiBbaMoirj11lvj\nL3/5Syxfvjx+/OMfb/aaLWvWrIn27dtHx44dY/ny5XHVVVc1ecymnHnmmfHII4/EtGnTYv369bFs\n2bL405/+tNnnRER84xvfiHnz5sWUKVOioaEhGhoa4tlnn425c+dGRUVFDB8+PCZOnBhr166NF198\nMf7rv/7rU39j9Nprr421a9fGnDlzYvLkyXHGGWdExEffLF577bVRV1cXdXV1cfXVV2/Rdch69OgR\ny5Yti1WrVpVvmzJlSixdujQiIjp27BilUinatNn0n9nq6uro169fPP300+XbiqKIBx98MJ544olY\nt25djB8/Pr72ta81ipobb7wxVq5cGW+88UbcfPPN5dcUEfH444/HSSed9Jlfx65E8OyE6uvr4/LL\nL4/u3btHz549o66uLq677rqI+P/vOsaNGxfDhg2LIUOGxOrVq2P69Okxe/bsqK2tje7du8cFF1xQ\n/kLd+PoP1dXVcd9998VPfvKT2HvvvaOmpiYmTZpUnqwuueSS+PWvfx1dunSJ0aNHR5cuXWL69Okx\nadKk6NatW9x4440xffr06NKlS3m8G09MO8O1OCCTUqkUo0aNiuOOOy7222+/2H///WPcuHGN7t/Y\n6NGjY+3atdGtW7c48sgj48QTT9zkMc19TdfU1MSDDz4YkyZNiq5du8ahhx4azz///CaP++RxKisr\nY8aMGTF16tSoqqqKnj17xuWXX17+7al/+7d/izVr1sQ+++wT559/fpx//vmf+roHDRoU/fr1iyFD\nhsRll10WQ4YMiYiP5sYBAwbEwQcfHAcffHAMGDBgs+fjk+P9+P7+/fvHyJEjo7a2Nrp06RJLliyJ\nhx9+OA466KCorKyMSy+9NKZOnRp77LFHk8e68MILG+0d+vh9uuqqq6Jr167xxz/+MaZMmdLoOaec\nckp89atfjUMPPTSGDh3a6DxMnTo1Lrzwwk89L7uyUuHbbwCSWLBgQdTW1sb69eubXF3ZUaxbty4O\nPfTQeOyxxxrtbWpOmzZt4tVXX21yj9QDDzwQd955Z5O/iML/EzwApLGzBM9ntbngoWXy/NcAAJHz\nivAZX9P2ZoUHAEjPCg8AkF7bzd1ZKk3cTsMAdiRFMbG1h/C5MIfBrqe5+csKDwCQnuABANITPABA\neoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQn\neACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIH\nAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA\n9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP\n8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHptW3sAfD4mxFWf6/Guigmf6/EAoDVZ4QEA\n0hM8AEB6ggcASE/wAADp2bS8E2rpBuWWbjxu6njNfQ6bmQHYGVnhAQDSEzwAQHqCBwBIT/AAAOnZ\ntJzE1mwmbuq5zW1abup2G5kB2NFZ4QEA0hM8AEB6ggcASE/wAADp2bRMk5rbiNzSqzwDwI7ECg8A\nkJ7gAQDSEzwAQHqCBwBIz6blHZxNwgCw9azwAADpCR4AID3BAwCkJ3gAgPQEDwCQnt/S2gk19799\naC1N/SbZjjZGAHZtVngAgPQEDwCQnuABANITPABAejYt70B2hv+NRFObkXeGcQOwa7PCAwCkJ3gA\ngPQEDwCQnuABANKzaXkH54rFALD1rPAAAOkJHgAgPcEDAKQneACA9GxaZpto6urLNmAD0Fqs8AAA\n6QkeACA9wQMApCd4AID0bFpuJU1t6t1ZNbUZOdPrA2DnZ4UHAEhP8AAA6QkeACA9wQMApGfT8g7E\nlYgBYNuwwgMApCd4AID0BA8AkJ7gAQDSs2mZ7aapqy/bqA3A9mCFBwBIT/AAAOkJHgAgPcEDAKRn\n0/J20NRm3eya2oy8K54HAHYMVngAgPQEDwCQnuABANITPABAejYttxJXGP6Iqy8DsD1Y4QEA0hM8\nAEB6ggcASE/wAADp2bTMduPqywC0Fis8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuAB\nANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkF7b1h5AJhPiqm3y2F1NU+fmqpjQCiMB\nxhXrNrnt2tLurTCSHY+5audihQcASE/wAADpCR4AID3BAwCkZ9PyFvosm45tYvtsbOiGHYcNys0z\nt+9crPAAAOkJHgAgPcEDAKQneACA9GxaboGWbqK1gW3bcUVTALaGFR4AID3BAwCkJ3gAgPQEDwCQ\nnk3LW8iG2W2nqXPr6ssAbA0rPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwA\nQHqutPwJrui782jqvXIFbACaYoUHAEhP8AAA6QkeACA9wQMApGfTcgvYCNv6mnoPbDAHoKWs8AAA\n6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOm50jKpNHX1ZVfKBsAKDwCQnuAB\nANITPABAeoIHAEjPpuVPsMF15+G9AqClrPAAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6\nggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4\nAID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcA\nSE/wAADpCR4AID3BAwCkVyqKomjtQQAAbEtWeACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4A\nID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDS\nEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3B\nAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwA\nQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCk\nJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqC\nBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gA\ngPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBI\nT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQE\nDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AA\nAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQ\nnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJ\nHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuAB\nANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAg\nPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANIT\nPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcED\nAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABA\neoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQn\neACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIH\nAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA\n9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP\n8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQP\nAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA\n6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe\n4AEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6Qke\nACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA\n0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9\nwQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8\nAEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMA\npCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6\nggcASE/wtJKJEyfG2Wef/bk/b/DgwXH77bdvzdA+d5MnT46jjz66RY/d0vOytc8FGttWc9SOaFvN\nm9/5znfi2muv/dyPy5Zp29oD2FWVSqVt8rxSqbTFx94RbM3Yd+bXDTuabTVH7Yi21bz57//+7y1+\n7Lnnnhu9e/eOa6655nMfBx+xwtNKiqLYrs/7PKxbty7efffdbfo5tub1tea5gWy21xz19ttvb9Hn\ngc9K8Gxj119/fVRXV0eHDh2if//+8dhjj23ymIaGhhg5cmScdtpp0dDQEIsXL45TTz019t5776it\nrY1bbrnlM33OV199NQYOHBgdO3aMYcOGxYoVKyIiYsGCBdGmTZv4xS9+EVVVVdGrV6+YNGlSi4+7\ndOnSqKmpibPOOiseffTR+PDDDz/TuD52ySWXRE1NTXTs2DEGDBgQv//978v3lUql+OCDD2LEiBHR\noUOH+OpXvxrPP/98+f6tPTdAY60xR23svPPOi4EDB8bPf/7zWLlyZYufN3fu3Dj22GOja9eu0b9/\n/5g2bVr5vnPPPTe++93vxtChQ6NDhw5xxBFHxPz588v3/+53v4v+/ftHp06d4nvf+14URdFsqE2c\nODFOO+20Zuekl156KQYPHhydO3eOgw46KB544IFG4xg/fnxERMycOTOqq6vjX/7lX6JHjx7Rq1ev\nmDx5ckRE3HbbbXHXXXfFDTfcEJWVlXHKKadERMveG1pO8GxDL7/8ctx6663x3HPPxapVq2LGjBnR\np0+fRo/54IMPYtiwYdG+ffuYNm1aVFRUxMknnxyHHnpoLF68OB599NG46aabYsaMGS36nEVRxK9+\n9av4z//8z1iyZEm0bds2Lr744kaPmTlzZrz66qsxY8aMuP766+PRRx9t0bGrqqpi3rx5ceihh8al\nl14atbW1MWHChHj99ddb9PyPHX744fGnP/0pVqxYEWeeeWacfvrpsW7duvL477vvvvjmN79Zvn/Y\nsGGxYcOG+PDDD7fq3ACNtcYc9Un3339/XHHFFfHwww9Hnz59YtSoUfHII49sdqXovffei2OPPTbO\nOuusWLp0aUydOjX+6Z/+KV566aXyY+65556YOHFirFixIvr16xc/+tGPIiKirq4uTj311PjJT34S\ny5Yti/322y+eeOKJzf5I6/77729yTmpoaIiTTz45TjjhhFi6dGnccsstMWrUqJg3b15EbPqjsrff\nfjtWrVoVixcvjttvvz2++93vxrvvvhsXXHBBjBo1KsaMGROrV6+O++67r0XvDZ+N4NmGKioqor6+\nPubMmRMNDQ1RU1MTtbW1EfHRF8KqVavi+OOPj/333z9++ctfRqlUimeffTbq6upi3Lhx0bZt2+jb\nt298+9vfjqlTp7boc5ZKpfiHf/iH+NKXvhR77rlnXHPNNfHf//3fjSaPCRMmRPv27eOggw6K8847\nL+6+++4Wv6YePXrED37wg3j++efj3nvvjZUrV8bAgQPj61//eqPvejZn1KhR0blz52jTpk18//vf\nj/r6+nj55ZfL9w8YMCCGDx8eFRUV8f3vfz8++OCDePLJJ7f63ACNtcYc9Ult27aNU045Je69996Y\nP39+HHHEETFmzJjo06dP3HrrrU0+Z/r06dG3b98455xzok2bNnHIIYfE8OHDG63yDB8+PAYMGBAV\nFRUxatSomD17dkREPPjgg3HQQQeV55jRo0fHPvvss9kxNjcnPfXUU/Hee+/F2LFjo23btvH1r389\nhg4d2mhO3Xju3W233eLKK6+MioqKOPHEE2OvvfZqNPdt/NjNvTdsGcGzDfXr1y9uuummmDhxYvTo\n0SNGjhwZS5YsiYiP/sN+6qmn4oUXXogxY8aUn7Nw4cJYvHhxdO7cufznuuuui3feeWeT4//jP/5j\nVFZWRmVlZfzzP/9z+fbevXuX/15TUxMNDQ1RV1fX7P2LFy/e5NiLFi0qH7tDhw7Nvr6DDz449t9/\n/3j55ZdbvL/nxhtvjC996UvRqVOn6Ny5c7z77ruNxlddXV3+e6lUiurq6li8eHEsWrSoxecG+HTb\neo76pBNPPLE8rzT1jVbnzp3jy1/+chxyyCGxcuXKWLBgQZPHWbhwYTz99NONxnDXXXeV9wOVSqXo\n0aNH+fFPqMwcAAAJ50lEQVTt27ePNWvWRMRHPxbfeI6JaDwnNqW5OWnJkiWbPHffffdtck6NiOja\ntWu0afP//+zuueee5XF90ubeG7aM4NnGRo4cGbNmzYqFCxdGqVRqNHEcd9xxMXbs2DjmmGPKk0VN\nTU307ds3VqxYUf6zatWqmD59ekQ0/g2In/3sZ7F69epYvXp1jB07tnz7okWLGv19t912i27dujV7\nf1VV1SbjrqmpKR971apV5ds3bNgQDz30UIwcOTL23XffeOihh+KKK66IN998s0W/ej5r1qz46U9/\nGtOmTYuVK1fGihUromPHjo2+s3njjTfKf//www/jzTffjKqqqujdu3eLzw3QMttyjvqkhx56qDyv\njBw5snz7K6+8EuPHj4/a2toYPXp0HHzwwTF//vz46U9/2uRxampqYtCgQY3GsHr16mZXhDbWq1ev\nRnNMURSNPm5Kc3PSx8faeP5auHBhozm1pfNSU4/b3HvDZyd4tqF58+bFY489FvX19bHHHntEu3bt\noqKiotFjLrvssjjzzDPjmGOOiWXLlsVhhx0WlZWVccMNN8TatWtjw4YN8cILL8Rzzz0XEZ/+GxBF\nUcSUKVPipZdeivfffz+uvPLKOP300xt9MV177bWxdu3amDNnTkyePDnOOOOMFr2ed955J6qrq2Pc\nuHFx5JFHxmuvvRa//vWv4xvf+Eaj71o2Z/Xq1dG2bdvo1q1brFu3Lq6++upGQRUR8Yc//CF+85vf\nxPr16+Omm26Kdu3axRFHHLHV5wZorDXmqE86//zz48gjj4xVq1bFb37zm5g9e3Zccskl0bVr12af\nM3To0Jg3b15MmTIlGhoaoqGhIZ599tmYO3fup47hpJNOijlz5pTnmJtvvjneeuutzY6xuTnp8MMP\njz333DNuuOGGaGhoiJkzZ8b06dNjxIgR5XG09Hz06NGj0cbqlrw3fDaCZxuqr6+Pyy+/PLp37x49\ne/aMurq6uO666yKi8Wa2cePGxbBhw2LIkCGxevXqmD59esyePTtqa2uje/fuccEFF5Sj4NOuF/Hx\nHp5zzz03evbsGevWrYubb7650WMGDRoU/fr1iyFDhsRll10WQ4YMadHr+cIXvhAzZsyIP/zhD/G9\n730vunTp0qLnbTzmE044IU444YQ44IADok+fPtG+ffuoqalp9Nhhw4bFPffcE126dIk777wz7r33\n3qioqIiKioqtOjdAY60xR33Sd77znViyZEn867/+axxyyCEtes5ee+0VM2bMiKlTp0ZVVVX07Nkz\nLr/88vIvPzQ1ho8/7tatW0ybNi3Gjh0b3bp1i1dffTWOOuqoZj9XqVSKU045pck5affdd48HHngg\nHnrooejevXtcdNFFcccdd8QBBxzQ5Dg2d16+9a1vxYsvvhidO3eO4cOHb/a9YcuUCt8W7zIWLFgQ\ntbW1sX79+havyADsyq666qp49dVX44477mjtobCV/KsHAM2wJpCH4NnF+JEPQMv5UXkefqQFAKRn\nhQcASG+z/7f0UmnidhoGsCMpiomtPYTPhTkMdj3NzV9WeACA9AQPAJCe4AEA0hM8AEB6ggcASE/w\nAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8A\nkJ7gAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADp\nCR4AID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7g\nAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4A\nID3BAwCkJ3gAgPQEDwCQnuABANJr29oD4HMyceKOfTwAaEVWeACA9AQPAJCe4AEA0hM8AEB6Ni3v\njFq6oXhrHtfcc21mBmAnZIUHAEhP8AAA6QkeACA9wQMApGfTchZbs5l4azct28gMwA7OCg8AkJ7g\nAQDSEzwAQHqCBwBIz6ZlmuZKywAkYoUHAEhP8AAA6QkeACA9wQMApGfT8o7OJmEA2GpWeACA9AQP\nAJCe4AEA0hM8AEB6ggcASM9vae2MdrTf3GpqPDvaGAHYpVnhAQDSEzwAQHqCBwBIT/AAAOnZtLwj\n2Rk2+tqgDMBOyAoPAJCe4AEA0hM8AEB6ggcASM+m5R2dDcEAsNWs8AAA6QkeACA9wQMApCd4AID0\nbFpm23BFZgB2IFZ4AID0BA8AkJ7gAQDSEzwAQHo2LbeWTBt4bVAGYAdnhQcASE/wAADpCR4AID3B\nAwCkZ9PyjsRGXwDYJqzwAADpCR4AID3BAwCkJ3gAgPRsWmb7cUVmAFqJFR4AID3BAwCkJ3gAgPQE\nDwCQnk3L28OuuDHXBmUAdiBWeACA9AQPAJCe4AEA0hM8AEB6Ni23Fht4P2JzMwDbgRUeACA9wQMA\npCd4AID0BA8AkJ5Ny2w/NigD0Eqs8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDSEzwAQHqCBwBI\nT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6bVt7AJmMK9Y1efu1pYmb3Fb821XbeDQ7\nh9JFm942ITY9N1fFhO0wGgCyssIDAKQneACA9AQPAJCe4AEA0rNpeQs1t0G5pUoX2YQLANuLFR4A\nID3BAwCkJ3gAgPQEDwCQnk3LLdDSDcrXlnbfxiPZdbn6MgBbwwoPAJCe4AEA0hM8AEB6ggcASM+m\n5S1kg/K209Rm5KY2LQNAS1nhAQDSEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA\n0nOl5U8YV6xr7SHQQk1dfbmpqzQDgBUeACA9wQMApCd4AID0BA8AkJ5Nyy1wbWn31h7CLq+pzchN\nbVoGgKZY4QEA0hM8AEB6ggcASE/wAADpCR4AID3BAwCkJ3gAgPQEDwCQnuABANJzpWVSaerqy01d\npRmyGles2+Q2V4vfeZjDth0rPABAeoIHAEhP8AAA6QkeACC9UlEURbN3liZux6EAO4qimNjaQ/hc\nmMNg19Pc/GWFBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4A\nID3BAwCkJ3gAgPQEDwCQnuABANITPABAeoIHAEhP8AAA6QkeACA9wQMApCd4AID0BA8AkJ7gAQDS\nEzwAQHqCBwBIT/AAAOkJHgAgPcEDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASE/wAADpCR4AIL1S\nURRFaw8CAGBbssIDAKQneACA9AQPAJCe4AEA0hM8AEB6ggcASO//AHco2L/yBwQyAAAAAElFTkSu\nQmCC\n"
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "endpoints = np.zeros(skel.shape)\nedges = np.zeros(skel.shape)\n\nfor l in range(1,n+1):\n    cur_edge = edge_lab == l\n    cur_end_points = endPoints(cur_edge)\n    pruned_edge = pruning(cur_edge,1)\n    edges = np.logical_or(pruned_edge, edges)\n    endpoints = np.logical_or(endpoints,cur_end_points)\n    \nlab_bp, nbp = mh.label(bp)\nlab_ep, nep = mh.label(endpoints+ep)\npruned_edges, ned = mh.label(edges)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "plt.figsize(13,8)\nsubplot(321,frameon = False, xticks = [], yticks = [])\ntitle(str(np.unique(skel)))\nimshow(skel, interpolation='nearest')\n\nsubplot(322,frameon = False, xticks = [], yticks = [])\ntitle(str(np.unique(lab_bp)))\nimshow(lab_bp, interpolation='nearest')\n\nsubplot(323,frameon = False, xticks = [], yticks = [])\ntitle(str(np.unique(edge_lab)))\nimshow(edge_lab, interpolation='nearest')\n\nsubplot(326,frameon = False, xticks = [], yticks = [])\nedge_mask = lab_ep>0\nep_edgelab = edge_lab*edge_mask\ntitle(str(np.unique(ep_edgelab)))\nimshow(ep_edgelab,interpolation='nearest')#lab_ep\n\nsubplot(324,frameon = False, xticks = [], yticks = [])\ntitle(str(np.unique(pruned_edges)))\nimshow(pruned_edges, interpolation='nearest')\n\nsubplot(325,frameon = False, xticks = [], yticks = [])\ntitle(str(np.unique(lab_ep)))\nimshow(lab_ep, interpolation='nearest')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 66,
       "text": "<matplotlib.image.AxesImage at 0xd4333ac>"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAHbCAYAAAAK1NqsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHgZJREFUeJzt3X9s3PV9+PHXOQ4hiX8kyxJmCiUJjI2Ghbb5oZUxRiot\nzdStzkQHbRGhWRh0SqhR+g/RYHEGRGWspPwQq5JlUKuNOq0TaGNUmtaE0ooxxkpw5xA8iqCkjCyl\nxk5GcTz78/1jX4xNz7FDnPjlu8dDQrrcve/ufW1076fe77tLqSiKIgAAkqqZ6AkAAByLWAEAUhMr\nAEBqYgUASE2sAACpiRUAIDWxAgCkNiljpaamJurq6uKWW2454cd58cUXx2lWE2fnzp1RX19fMa8H\nIIvxWm/GQzW/10/KWImIaG9vj1tvvTUiIl566aWoqamJ+vr6wf8+9KEPTej8vvvd7w7Opa6ubtj8\nGhoa4sCBA+P2XOvWrYvDhw+P2+MB8I6h601ExN69e2PJkiUxc+bMWLp0aTz77LMj3nfPnj2xYsWK\nmDVrVixYsOCYz9PZ2RnNzc0xb968mDNnTqxatSo6OzsHb6/m9/pJGyvldHd3x+HDh+Pw4cPxzDPP\nTOhcfvM3f3NwLh0dHcPm19PTE2edddbg2P7+/omaJgDH4ejRo9Hc3Bxr1qyJN954I6655ppobm6O\nvr6+suPr6uri2muvjTvvvHPUx+7u7o7Vq1dHZ2dnHDx4MJYvXx7Nzc3j/RImpYqKlXd76qmn4iMf\n+UjMnj07zjzzzLjhhhtG/Av16KOPxqJFi6KhoSHOOuus+NKXvjR42yOPPBIf/OAHY/bs2fEbv/Eb\n8YMf/OC45vHuf9GgtbU1PvnJT8bVV18djY2N8eCDD8ZnP/vZYduMjz32WJx99tmDf3711Vfj8ssv\nj3nz5sXChQvj3nvvPa45AHDiHnvssejv74+WlpaYOnVq3HDDDVEURezevbvs+GXLlsVVV1016q7K\n22PXrl0bs2bNitra2rjxxhvj+eefj66urvF+GZNORcXKu6OgtrY27r777nj99dfjX/7lX+Lb3/52\n3H///WXvu27duti+fXv09PRER0dHfPSjH42IiGeeeSbWrVsXO3bsiJ/+9Kdx/fXXxyc+8Yk4evTo\nCc317//+7+MP/uAPoru7O6666qoolUpRKpXKjh0YGIjf+73fiw996EPx6quvxre//e348pe/HP/0\nT/90QnMA4Ph0dHTE4sWLh1130UUXDe6gj6fHH388mpqaYvbs2eP+2JNNRcXKL/7iL8bs2bNj9uzZ\ncdddd8WHP/zhWL58edTU1MQ555wT1113XXznO98pe9/TTjstOjo6oqenJxobGwc/87J9+/a4/vrr\nY9myZVEqlWLNmjUxbdq0ePLJJ09orhdffHF84hOfiIiI008/PSJ+Prbe9m//9m/xk5/8JG6++eao\nra2NBQsWxLXXXhvf+MY3TmgOAByfI0eORGNj47DrGhoaxv2zJAcOHIgNGzbEXXfdNa6PO1nVTvQE\nxtPrr78eNTXv9FdnZ2ds3Lgx/v3f/z3efPPN+N///d9YunRp2fv+3d/9Xdx2221x0003xeLFi+OL\nX/xi/Pqv/3q8/PLL0dbWNuzYpa+vL/7rv/7rhOY69DMro3n55Zfj1VdfHVbX/f39cemll57QHAA4\nPvX19dHT0zPsuu7u7mhoaBi35zh06FCsXLky1q9fH1deeeW4Pe5kVlE7K+/2x3/8x/GBD3wgXnjh\nheju7o7bb789BgYGyo5dunRpPPzww3Ho0KFYvXp1XHHFFRER8f73vz/+5E/+JLq6ugb/O3LkyAn9\nBSp35DNz5sx48803B//82muvDV4+++yzY8GCBcPm0NPTE4888sh7ngMAx2/RokXR3t4+7Lr29vZY\ntGjRuDx+V1dXrFy5MlavXh2bNm0al8esBBUdK0eOHIn6+vqYMWNG7N+/P/7yL/+y7Li+vr74+te/\nHt3d3TFlypSor6+PKVOmRETEH/3RH8VXvvKVeOqpp6Ioivif//mf+Md//Mc4cuTIe55XueOeD37w\ng/Hoo49GV1dXvPbaa/HlL3958Lbly5dHfX19/Pmf/3n87Gc/i/7+/viP//iPePrpp9/zHAA4fpdd\ndllMmTIl7rnnnujt7Y177rknampqBj/n+G5FUcRbb70VfX19URRF9Pb2jviZx56envjYxz4Wl1xy\nSWzduvVkvoxJZ9LGyrsX/HIfTv2Lv/iL2LVrVzQ0NMR1110Xn/rUp4aNG3r5a1/7WixYsCAaGxtj\n+/bt8fWvfz0iIpYsWRI7duyIDRs2xC/8wi/EL//yL0dbW9txz/fdz/vu+V599dVx0UUXxfz582PV\nqlXD5jplypR45JFHYu/evbFw4cKYO3duXHfddcO2Ikf6vAsAJ2bo++vUqVPj4Ycfjra2tpg9e3a0\ntbXFww8/HLW15T9V8Z3vfCdmzJgRH//4x+OVV16J6dOnx6pVq8qOfeihh+Lpp5+OBx54YMTf5arW\n9/pSMQlf+fTp02PatGnR0tISW7ZsmejpTLgHHnggNm7cGL29vbFv376YP3/+RE8JoCJkWm+q+b1+\nUsYKAFA9Ju0xEABQHY751eVSqfUUTYNKVxStEz0F4CSyXjBeyq0XdlYAgNTECgCQmlgBAFITKwBA\namIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBSEysAQGpiBQBITawAAKmJFQAgNbECAKQm\nVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAamIF\nAEhNrAAAqYkVACA1sQIApCZWAIDUaid6Ahltji1jHrslNp/EmQAAdlYAgNTECgCQmmOg/2+ko5+R\njnneHj/0fo6EAGD82VkBAFITKwBAao6ByhjLcc7bY4YeAzkSAoDxZ2cFAEhNrAAAqTkGOkFDj3uO\n58fkAICxsbMCAKQmVgCA1Kr6GMixDQDkZ2cFAEhNrAAAqVX1MdBQ4/0jbn4gDgDGh50VACA1sQIA\npFZ1x0An8xtAfiAOAMafnRUAIDWxAgCkVnXHQEP5lg4A5GdnBQBITawAAKlV9THQqeIH4gDgvbOz\nAgCkJlYAgNSq4hhoIn6gzQ/EAcD4sLMCAKQmVgCA1KriGGgo38YBgMnFzgoAkJpYAQBSq7pjoInm\nB+IA4PjYWQEAUhMrAEBqFXsMlOmH2PxAHAC8d3ZWAIDUxAoAkFrFHgMNlfVbN74ZBACjs7MCAKQm\nVgCA1KriGCgT3wwCgONjZwUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAahX1o3Aj\n/cjaZPjxNf9OEMDEuLk4Onj5ttJpEziTkVX7GmFnBQBITawAAKlN+mOgkY54Jts22WQ4qgKoRFmP\nfoaabGvaeLOzAgCkJlYAgNQm5TFQpRz9jKTaP/UNAEPZWQEAUhMrAEBqk/IYaKhKOSYZ+jp8MwgA\n3mFnBQBITawAAKmJFQAgNbECAKQ2aT5gW60fOvWbKwBUOzsrAEBqYgUASG3SHAMNVenHIX5zBQDe\nYWcFAEhNrAAAqYkVACA1sQIApCZWAIDUJuW3gaqVH4gDoBrZWQEAUhMrAEBqk+YYqFqPPar1dQPA\n2+ysAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDU\nxAoAkJpYAQBSEysAQGpiBQBITawAAKmJFQAgNbECAKQmVgCA1EpFURQTPQkAgJHYWQEAUhMrAEBq\nYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqVVMrNTU1ERdXV3c\ncsstEz2Vk2rnzp1RX18fNTU18eKLL070dAAmHevF5FMxsRIR0d7eHrfeeuvgn/fu3RtLliyJmTNn\nxtKlS+PZZ58d8b579uyJFStWxKxZs2LBggXHfJ7Ozs5obm6OefPmxZw5c2LVqlXR2dk54vht27bF\nueeeGw0NDXHGGWfE2rVr4/Dhw6O+nra2tqipqYmdO3cOXrdu3box3ReAkVkvJpeKipWhjh49Gs3N\nzbFmzZp444034pprronm5ubo6+srO76uri6uvfbauPPOO0d97O7u7li9enV0dnbGwYMHY/ny5dHc\n3Dzi+Obm5nj66aejp6cn9u/fHz/60Y/i9ttvP+ZzdHV1xdatW+PCCy+MUqk06pwAeG+sF/lVbKw8\n9thj0d/fHy0tLTF16tS44YYboiiK2L17d9nxy5Yti6uuumrUSn577Nq1a2PWrFlRW1sbN954Yzz/\n/PPR1dVVdvzChQtj9uzZERExMDAQNTU10dTUdMzn2LRpU7S0tMScOXNGnQ8A7531Ir+KjZWOjo5Y\nvHjxsOsuuuii6OjoGPfnevzxx6OpqWnwL1g5u3btisbGxpg7d27MnTs3WlpaRhz71FNPxfe///34\n3Oc+N+5zBWA460V+FRsrR44cicbGxmHXNTQ0jPv53YEDB2LDhg1x1113HXPcZz7zmeju7o7Ozs54\n7rnnYtu2bWXH9ff3x/r16+O+++6r2O08gEysF/lVbKzU19dHT0/PsOu6u7ujoaFh3J7j0KFDsXLl\nyli/fn1ceeWVY7rPeeedFzfddFO0tbWVvf3++++PxYsXx/LlywevK4piXOYLwM+zXuRXsbGyaNGi\naG9vH3Zde3t7LFq0aFwev6urK1auXBmrV6+OTZs2Hdd9+/r6YsaMGWVv2717dzz00EPR1NQUTU1N\n8cQTT8QXvvCF+PznPz8e0wbgXawX+VVsrFx22WUxZcqUuOeee6K3tzfuueeeqKmpiY9+9KNlxxdF\nEW+99Vb09fVFURTR29sbR48eLTu2p6cnPvaxj8Ull1wSW7duHXUuf/VXfxWHDh2KiIh9+/bFF7/4\nxbj88svLjn3wwQdj//798eyzz8bevXtj6dKl0draOuqnwQF4b6wXk0BRIUqlUvHCCy8Mu+6ZZ54p\nlixZUkyfPr1YsmRJsXfv3hHvv2fPnqJUKhWlUqmoqakpSqVSsWLFirJjH3zwwaJUKhUzZ84s6urq\nirq6uqK+vr545ZVXyo5fu3ZtccYZZxR1dXXF+eefX9xxxx3FwMDAmF7XZZddVuzcuXPYdQMDA0Wp\nVCp++MMfjukxAHiH9WLyKRVFZRxwTZ8+PaZNmxYtLS2xZcuWiZ7OSfPAAw/Exo0bo7e3N/bt2xfz\n58+f6CkBTCrWi8mnYmIFAKhMFfuZFQCgMtQe68ZSqfUUTYNKVxStEz0F4CSyXjBeyq0XdlYAgNTE\nCgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBSEysAQGpiBQBITawA\nAKmJFQAgNbECAKQmVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQ\nmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDUaid6Aim1tp6csQBUFuvFKWFnBQBITawA\nAKk5BnrbSNtzo10/9HZbfACVz3pxytlZAQBSEysAQGqOgcoZy/bcaNt6tvgAKp/14pSwswIApCZW\nAIDUHAOdKFt5AIyF9eI9s7MCAKQmVgCA1Kr7GMg2HABjYb2YUHZWAIDUxAoAkFp1HwMNNd5bfD71\nDVCZrBennJ0VACA1sQIApFZ9x0Anc4vNVh5A5bBepGFnBQBITawAAKlV3zHQULbeABgL68WEsrMC\nAKQmVgCA1Kr7GOhU8alvAMbCelGWnRUAIDWxAgCkVh3HQBOxlWYrD2DysV6kZGcFAEhNrAAAqVXH\nMdBQttgAGAvrRRp2VgCA1KpvZ2Wi+SAVAGNhvRhkZwUASE2sAACpVe4xUKYtM1t5AHllel+2XpRl\nZwUASE2sAACpVe4x0FBZt9Js9wHkkvW9uMrXCzsrAEBqYgUASK06joEyqfKtPADGyHoxyM4KAJCa\nWAEAUhMrAEBqYgUASE2sAACp+TYQANVtMnzTZjLM8SSyswIApCZWAIDUKuoY6Obi6ODl20qtg5eL\n+7ZMwGxGV9rwzuXN8c4ct8TmCZgNAEMNXVMmwm2l0yb0+TOxswIApCZWAIDUJv0x0Fi26UobHKsA\nVLuRPhIw0hqR9RimGj82YGcFAEhNrAAAqU3KY6CRjn6ybtkdr2rc4gM4GY736GeyqZb1ws4KAJCa\nWAEAUpuUx0BDVcrRz9Dtu6HbegCMj0o5+qnG9cLOCgCQmlgBAFITKwBAamIFAEhNrAAAqU2abwNN\n9D/VPVGq5Qd/AMbLSD8EV+kqeb2wswIApCZWAIDUJs0x0FCV8kNwI6nGH/wBOBkq5YfgRlIt64Wd\nFQAgNbECAKQmVgCA1MQKAJCaWAEAUpuU3waqVpX8gz8A46XSvwFUjeysAACplYqiKEa8sdR6CqdC\nJSuK1omeAnASWS8YL+XWCzsrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBA\namIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBSEysAQGpiBQBITawAAKmJFQAgtVJRFMVE\nTwIAYCR2VgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFIT\nKwBAaqc0VmpqaqKuri5uueWWU/m0jNHOnTujvr4+ampq4sUXX5zo6QCcFNWyFlXSe/op31lpb2+P\nW2+9dfDPe/fujSVLlsTMmTNj6dKl8eyzz4543z179sSKFSti1qxZsWDBgmM+T2dnZzQ3N8e8efNi\nzpw5sWrVqujs7Bxx/LZt2+Lcc8+NhoaGOOOMM2Lt2rVx+PDhUV9PW1tb1NTUxM6dO0cc09raGlOn\nTo36+vqor6+PhoaGeOmll0Yc39/fHzfffHO8733vi4aGhvjwhz8c3d3dZcf+zu/8zuDj1tfXx7Rp\n02Lx4sUjPvb3vve9WLZsWTQ2Nsa5554bO3bsGLxt3bp1Y3rNAJNdNaxFlfSePqHHQEePHo3m5uZY\ns2ZNvPHGG3HNNddEc3Nz9PX1lR1fV1cX1157bdx5552jPnZ3d3esXr06Ojs74+DBg7F8+fJobm4e\ncXxzc3M8/fTT0dPTE/v3748f/ehHcfvttx/zObq6umLr1q1x4YUXRqlUGnFcqVSKT3/603H48OE4\nfPhw9PT0xPz580ccv3nz5njyySfjySefjJ6envja174Wp59+etmx3/rWtwYf9/Dhw3HxxRfHFVdc\nUXZsf39//P7v/35cd9110d3dHX/zN38TGzdujPb29mO+ToBKVi1r0WQ2obHy2GOPRX9/f7S0tMTU\nqVPjhhtuiKIoYvfu3WXHL1u2LK666qpRS/btsWvXro1Zs2ZFbW1t3HjjjfH8889HV1dX2fELFy6M\n2bNnR0TEwMBA1NTURFNT0zGfY9OmTdHS0hJz5sw55riiKKIoilHnHPF/f+nuvvvu2LFjR5x99tkR\nEfGBD3wgpk2bNup9X3rppfjud78ba9asKXv7wYMH4/XXX4+rr746IiKWLl0aF1xwQTz33HNjmhtA\nJaqWtWgym9BY6ejo+Lkji4suuig6OjrG/bkef/zxaGpqGvxLUM6uXbuisbEx5s6dG3Pnzo2WlpYR\nxz711FPx/e9/Pz73uc+N+tylUin+4R/+IebMmRMXXnhhfOUrXxlx7A9+8IOora2Nv/3bv42mpqb4\nlV/5lbj//vtHfY6I/9sGvPTSS+P9739/2dvPPPPMWLx4cfz1X/919Pf3xxNPPBEvv/xyXHLJJWN6\nfIBKVC1r0WQ2obFy5MiRaGxsHHZdQ0PDuJ+xHThwIDZs2BB33XXXMcd95jOfie7u7ujs7Iznnnsu\ntm3bVnZcf39/rF+/Pu67774xbbldccUVsX///vjJT34SO3bsiD/7sz+Lb3zjGyPOtbu7O/7zP/8z\nXnrppfjmN78Zra2t8c///M+jPk9bW1t89rOfPeaY7du3x+bNm+P000+P3/qt34qtW7fG+973vlEf\nG6BSVctaNJlNaKzU19dHT0/PsOu6u7ujoaFh3J7j0KFDsXLlyli/fn1ceeWVY7rPeeedFzfddFO0\ntbWVvf3++++PxYsXx/LlywevO9YxzwUXXBC/9Eu/FKVSKT7ykY9ES0tLfPOb3yw7dvr06RER8ad/\n+qcxbdq0+LVf+7X41Kc+FY8++ugx5/y9730vDh48GJ/85CdHHPPjH/84fvd3fzd27doVfX190dHR\nEXfccceojw1QyaplLZrMJjRWFi1a9HMf7mxvb49FixaNy+N3dXXFypUrY/Xq1bFp06bjum9fX1/M\nmDGj7G27d++Ohx56KJqamqKpqSmeeOKJ+MIXvhCf//znT3jOI32TZ7Rq/upXvxqXX375iHOOiHji\niSfirLPOit/+7d+OiIjzzz8/Pv7xj8e3vvWt9z5hgEnOWjQJFKdQqVQqfvjDHw7++ejRo8U555xT\n3H333cVbb71V3H333cX8+fOLvr6+svcfGBgofvaznxWPPvpocc455xRvvfVW0dvbW3Zsd3d3sWzZ\nsmLDhg1jmtuOHTuK//7v/y6Koig6OjqKRYsWFV/60pfKjn3jjTeKgwcPFgcPHixee+214uKLLy62\nbdtW9PT0lB3/8MMPFz/96U+LgYGB4l//9V+LM888s2hraxtxLpdeemlx/fXXF729vcW+ffuKefPm\nFbt37x5x/Jtvvlk0NjYWe/bsOeZr3LdvXzFjxoxi9+7dxcDAQPHCCy8U5513XrFjx45h4979/xNA\nJam2tagS3tNPeay88MILw6575plniiVLlhTTp08vlixZUuzdu3fE++/Zs6colUpFqVQqampqilKp\nVKxYsaLs2AcffLAolUrFzJkzi7q6uqKurq6or68vXnnllbLj165dW5xxxhlFXV1dcf755xd33HFH\nMTAwMKbXddlllxU7d+4c8fZPf/rTxZw5c4q6urriV3/1V4t77733mI/34x//uFi1alVRV1dXLFy4\nsNi+ffsxx+/atauYP3/+mOb61a9+tbjggguK+vr64qyzzipuuummYbcPDAxUxF9sgJFU01pUKe/p\npaI4dQdc06dPj2nTpkVLS0ts2bLlVD0tY/TAAw/Exo0bo7e3N/bt23fM34IBmKyqZS2qpPf0Uxor\nAADHyz9kCACkJlYAgNRqj3VjqdR6iqZBpSuK1omeAnASWS8YL+XWCzsrAEBqYgUASE2sAACpiRUA\nIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBS\nEysAQGpiBQBITawAAKmJFQAgNbECAKQmVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWx\nAgCkJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBSEysA\nQGpiBQBITawAAKmJFQAgNbECAKQmVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCk\nJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIApCZWAIDUxAoAkJpYAQBSEysAQGpi\nBQBITawAAKmJFQAgNbECAKQmVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYA\ngNTECgCQmlgBAFKrnegJnDytI1xOpLW1/GUATp3J8F48GeZ4EtlZAQBSEysAQGoVdQw0t/jDwcuH\nSq2Dl28ujg5evq102qmc0s85UNw7ePms0jvXF/dtGbxc2rD5VE4JoOoMXxdaBy9nei8ePpd3rt8c\n71y/JapjvbCzAgCkJlYAgNRKRVEUI944ZGsMTkRRtE70FICTyHrBeCm3XthZAQBSEysAQGpiBQBI\nTawAAKmJFQAgNbECAKQmVgCA1MQKAJCaWAEAUhMrAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTE\nCgCQmlgBAFITKwBAamIFAEhNrAAAqdVO9ATGqjh7y+Dl0iubJ3AmE2dzvPO/wZaozv8NAEZzc3F0\n8PJtpdMmcCYTp9LWCzsrAEBqpaIoihFvLLWewqlQyYqidaKnAJxE1gvGS7n1ws4KAJCaWAEAUhMr\nAEBqYgUASE2sAACpiRUAIDWxAgCkJlYAgNTECgCQmlgBAFITKwBAamIFAEhNrAAAqYkVACA1sQIA\npCZWAIDUxAoAkJpYAQBSEysAQGpiBQBIrVQURTHRkwAAGImdFQAgNbECAKQmVgCA1MQKAJCaWAEA\nUhMrAEBq/w+cV8A+3i4sNAAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": "First make vertex from branched points and then link to the neighbouring endpoints"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "BPlabel=np.unique(lab_bp)[1:]\nEPlabel=np.unique(lab_ep)[1:]\nEDlabel=np.unique(pruned_edges)[1:]\n\nG=nx.Graph()\nnode_index=BPlabel.max()\nselected_ep_labels = []\nfor bpl in BPlabel:\n    #branched point location\n    bp_row,bp_col= np.where(lab_bp==bpl)[0][0], np.where(lab_bp==bpl)[1][0]\n    bp_neighbor= lab_ep[bp_row-1:bp_row+2,bp_col-1:bp_col+2]\n    \n    G.add_node(bpl,kind='BP',row=bp_row,col=bp_col,label=bpl)\n    #branched point neighborhood\n    neig_in_EP=lab_ep[bp_row-1:bp_row+2,bp_col-1:bp_col+2]\n    neig_inEplabel=np.unique(neig_in_EP)[1:]\n    print 'bp label',bpl,' pos',(bp_row,bp_col), 'ep neighbor label:',neig_inEplabel\n    \n    for epl in neig_inEplabel:\n        selected_ep_labels.append(epl)\n        node_index = node_index+1\n        #print 'bp label',bpl,' pos',(bp_row,bp_col),neig_inEplabel, 'current ep',epl\n        ep_row,ep_col= np.where(lab_ep==epl)[0][0], np.where(lab_ep==epl)[1][0]\n        G.add_node(node_index,kind='EP',row=ep_row,col=ep_col,label=epl)\n        G.add_edge(bpl,node_index, weight=1)\n        print 'bp label',bpl,':',(bp_row,bp_col),neig_inEplabel,' -ep',epl,'(',node_index,')',':',(ep_row,ep_col),\n        #print 'edge',(bpl, node_index)\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "bp label 1  pos (29, 16) ep neighbor label: [1 3 4]\nbp label 1 : (29, 16) [1 3 4]  -ep 1 ( 3 ) : (28, 16) bp label 1 : (29, 16) [1 3 4]  -ep 3 ( 4 ) : (29, 15) bp label 1 : (29, 16) [1 3 4]  -ep 4 ( 5 ) : (29, 17) bp label 2  pos (29, 30) ep neighbor label: [2 5 6]\nbp label 2 : (29, 30) [2 5 6]  -ep 2 ( 6 ) : (28, 30) bp label 2 : (29, 30) [2 5 6]  -ep 5 ( 7 ) : (29, 29) bp label 2 : (29, 30) [2 5 6]  -ep 6 ( 8 ) : (29, 31)\n"
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "figsize(5,7)\nnx.draw(G)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAIWCAYAAACbYsDtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUVPWdsPGnqqubXtibTUXQIBoV11FQXDAxoigqGJJh\nXpeIG4JmXCPHSUycN45xQ1wbd6NR0Yj7BkImKiqBAKOokTGMUSCyKDZb78t9/2h8hzTN0virpaue\nzzkcpe6tqi94Tj/+7lIVi6IoQpKkgOLpHkCSlH2MiyQpOOMiSQrOuEiSgjMukqTgjIskKTjjIkkK\nzrhIkoIzLpKk4IyLJCk44yJJCs64SJKCMy6SpOCMiyQpOOMiSQrOuEiSgjMukqTgjIskKTjjIkkK\nzrhIkoIzLpKk4IyLJCk44yJJCs64SJKCMy6SpOCMiyQpOOMiSQrOuEiSgjMukqTgjIskKTjjIkkK\nzrhIkoIzLpKk4IyLJCk44yJJCs64SJKCMy6SpOCMiyQpOOMiSQrOuEiSgjMukqTgjIskKTjjIkkK\nzrhIkoIzLpKk4IyLJCk44yJJCs64SJKCMy6SpOCMiyQpOOMiSQrOuEiSgjMukqTgjIskKTjjIkkK\nzrhIkoIzLpKk4IyLJCk44yJJCs64SJKCMy6SpOCMiyQpOOMiSQouke4BJGlrNmzYwMcff8zatWtp\n164dO++8M/369Uv3WNoG4yIpI3300UdMnjSJJ554gt3y8+kSi1ENfFpby+79+jF+wgRGjRpFYWFh\nukdVC2JRFEXpHkKSvlFVVcU5o0fz5owZnF9by/kNDfTeZHs98CpQ1r497+XlMeW55/je976Xpmm1\nJcZFUsaoqKjg+COPZLdFi3iwupp229j/j8Do4mLue+IJTj311FSMqO1kXCRlhCiKGDF0KF3efpuH\nq6uJbefz5gMnFBfz6htvcOihhyZzRLWCV4tJyghvvfUW/z17Nve1IiwA/wTcUFnJNZdemqzRtAOM\ni6SMUHbzzVxUWUnBVvb5K1AInNns8dOBBQsWsHjx4qTNp9YxLpLSbvny5bw+cyZnbeMo/UXAQNhs\nZVMInNPQwOTbb0/ShGot4yIp7d5++22GtGtHp63s8yTQBTgWaClBP6yr4z9ffTUp86n1jIuktCsv\nL6dbff0Wt68DfgVMouWwAHQDytetCz+cdohxkZR2iUSChtiWT+NfA5wH7Mzmh8S+UQ8k8vLCD6cd\n4h36ktKue/fuLNlCGN4D/gD818bfb2nlshToXloafjjtEOMiKe2OPfZYzq6v5zNgt2bb3gQ+A/ps\n/P0GoAH4GJi3yX6/LSriR2PGJHdQbTdvopSUES4bP56iBx7g+rq6f3i8Cli/8d8j4BaaYnMP8M06\n5Sugf2Eh//P3v9O1a9fUDKyt8pyLpIww7tJLeSCRYEWzx4uAHht/9QTab3xs0wNgNycSjDj1VMOS\nQVy5SMoY//eaa3hx0iRmVlTQeTuf82AsxnXduzP7/ffp1atXUufT9jMukjJGFEVcftFFzHj0UaZU\nVLDfVvatpGnF8lDnzkyfNYvvfve7qRpT28HDYpIyRiwWY1JZGRfdcAPHd+7MkA4deIqm+1wioI6m\nE/kXx2L0Lihg/tFHM/v99w1LBnLlIikj1dXV8fzzz1N2443MXbiQmvp6iMXYqXNn9jv4YKJEgtde\ney3dY2oLjIukNqGuro5EIkEsFmPdunX07duXRYsW0bNnz3SPphZ4WExSm5Cfn09s4138HTt2ZOTI\nkTz66KNpnkpb4spFUpv07rvvMmbMGBYtWvT/o5MuVVVVPPnkkzz729/y5apVNDQ00LVrV34wciRj\nzj2Xbt26pXW+dDAuktqkKIoYMGAAZWVlDBkyJC0zlJeXc90vf8kjDz/MoFiMszZsoA+QB6wCphYV\n8XxjI6cMH841v/kN/fv3T8uc6WBcJLVZt912G/PmzeOxxx5L+Xt/9tlnDBsyhCNXrODq2lq+s4X9\nVgOT43HuLCnhmVdf5cgjj0zlmGljXCS1WatXr6Zfv358+umnKb07/8svv+TwAw7gopUruayxcbue\nMx04s6SE199+mwMPPDC5A2YAT+hLarNKS0s58cQTU75yufDMMxn11VfbHRaA44G7Kir44bBhNDQ0\nJG+4DGFcJLVp559/Pvfffz+pOgjz+eef88abb3JNsw/YBDiGps8967Dx197Ntv8Y6F5Rwas58I2Z\nxkVSmzZkyBCqqqqYO3duSt7v3rvu4swooqSFbTHgbpo+xXk9TZ8m0Nz49eu5+8YbkzliRjAuktq0\neDzOeeedx/3335+S93vkgQcYW1Ozxe3bWj/9GPjz/PksX7486FyZxrhIavPOPvtsnnnmGdavX7/t\nnb+FhoYGVq5dy15b2edqoDtwJE1fdNZcIbBbu3b8/e9/T8aIGcO4SGrzevXqxfe+9z2mTJmS1Pep\nqqqiIC9viz84bwT+BnwBXACcDHzawn7FQGVlZXKGzBDGRVJW+ObEfjKVlJRQ29hI7Ra2DwRKgHzg\nLOAIoKVT92uiiE6dOiVpysxgXCRlhaFDh7Jq1Sree++9pL1HLBbjgH79+OO3eI0vgGV1dXznO1u6\n7TI7GBdJWSEvL49zzjmHBx54IKnvM+6qqyhr336zx9fSdKNkNVAPPA7MAk5ott/9iQT/Mno0HTp0\nSOqc6eYd+pKyxpIlSzjooINYunQpxcXFSXmPiooK+vTowYLKSvpu8vhXwInAIpo+W2xv4NfAsZvs\nUwfsVlTEtDlz2G+/rX3PZtvnykVS1ujTpw+DBg1i6tSpSXuPkpISLrnsMv65oIDqTR7vBsyl6Vsz\ny4F3+cewRMBF7dpxyODBWR8WMC6SskyyT+xHUUSiuJjFeXmcXFjIuu14TgNwSUEB8/r25bHnnkva\nbJnEuEjKKsOHD2fx4sV8/HFL98d/O42NjVx66aU8+eSTLFi0iD1Gj+agkhJuj8VY08L+NcATwOD2\n7fnwwAP5zzlzsv5cyzc85yIp61x99dXU1tYyceLEYK9ZU1PDT37yE5YvX84LL7xA586diaKId955\nh7JbbuG16dMZlpdH75oa8hob+bJdO14C9t9/f8ZPmMDJJ59MIpEINk+mMy6Sss7ixYsZPHgwS5cu\npV27dt/69datW8dpp51Gp06dePzxxyksLNxsn5UrV/Lyyy/z1VdfUV9fT5cuXTj22GPZa6+t3c+f\nvYyLpKz0/e9/nwsvvJAf//jH3+p1VqxYwYknnsigQYO46667yMvLCzRhdvOci6SsFOLE/uLFizni\niCMYMWIEZWVlhqUVXLlIykrV1dXsuuuuzJkzZ4fuhp8/fz4nn3wy1157LRdccEESJsxurlwkZaXC\nwkLOOOMMHnzwwVY/d8aMGQwbNoyysjLDsoNcuUjKWh999BFDhgzhuMGDeffddynfsIGCRIJepaWM\nPucczhs7ll69ev3Dc6ZMmcKll17K1KlTOeqoo9I0edtnXCRlpZdeeomfX3IJaz7/nMujiFOiiK5A\nLU0fi/9gYSFPRxEnHH88k+69l169enHbbbcxceJEXnvtNQYMGJDmP0HbZlwkZZ07Jk3ipp//nAer\nqhhK09cPt2QNcFMiwZSuXfnBKafw9ttvM23aNPr27buFZ2h7GRdJWeWR3/6Way+6iDcrK+mznc+5\nF5iQSPDu+++zzz77JHO8nGFcJGWNL7/8kr369uWdqir2buVzr04k+Pykk3ji+eeTMluu8WoxSVnj\nofvvZyRsFpa7gENo+v76MVt47oT6el6bPp2VK1cmc8ScYVwkZYWGhgbuue02xldVbbZtF+Aa4Jyt\nPL8zMAp48L77kjNgjjEukrLCnDlz6FhdzT+1sG0kcCpQuo3XOLe6mik7cF+MNmdcJGWFFStWsHts\nS9eFNdnWCebdgZWrVwebKZcZF0lZoa6ujoJtXJ+09fRAAVBTXx9splxmXCRlhc6dO/N1fOs/0ra1\nclkNdGnfPthMucy4SMoKhxxyCPNratjaQa1trVxejMU40o98CcK4SMoKpaWlnDJ8OA+3sHppAKqB\n+o3/XrPxn5tqBCYXFzP+yiuTPWpOMC6Sssb4K69kclERtc0e/zVQDNwIPAYUAf/RbJ/pQEnPnhx+\n+OHJHzQHeIe+pKwRRREjjz+ezrNm8XB19TYPg33jM+CIoiLu/f3vGT58eBInzB2uXCRljVgsxuPP\nPcd/9+/PmMLCzVYwLfkIOKa4mKuvv96wBGRcJGWVkpISZs6ezdqjj+a7JSXcHI/zVbN9ImAucHZR\nEUcWFnLdPfdw8aWXpmHa7OVhMUlZ689//jNlt9zC8y++yEHt2pG/bh357dvzOVBRVMS4yy9nzLnn\n0q1bt3SPmnWMi6Sst3r1aubPn8+oUaO488472WuvvRg4cCDxbdwXox1nXCTljM6dO/PZZ5/RuXPn\ndI+S9cy2pJxRW1tLQUFBusfICcZFUs4wLqljXCTlhIaGBhoaGsjLy0v3KDnBuEjKCXV1dRQUFBDb\nxsfyKwzjIiknfBMXpYZxkZQTPN+SWsZFUk4wLqllXCTlBOOSWsZFUk4wLqllXCTlhNraWvLz89M9\nRs4wLpJygiuX1DIuknKClyKnlnGRlBNcuaSWcZGUE4xLahkXSTnBuKSWcZGUE4xLahkXSTnBS5FT\ny7hIygmuXFLLuEjKCV6KnFrGRVJOcOWSWsZFUk4wLqllXCTlBOOSWsZFUk4wLqllXCTlBC9FTi3j\nIiknuHJJLeMiKSd4KXJqGRdJOcGVS2oZF0k5wbiklnGRlBOMS2oZF0k5wbiklnGRlBOMS2oZF0k5\nwftcUsu4SMoJrlxSy7hIygne55JaxkVSTnDlklrGRVJOMC6plUj3AJKULPX19UybNo1PPvmEv/3t\nb0ydOpU1a9YwbNgwEgl//CVTLIqiKN1DSFJIK1as4P7Jk7nvrrvoU1fHwJoa2tfWsj4/n3mFhXyW\nSHDBxRdz/rhx7LTTTukeNysZF0lZZebMmfyfkSM5rb6ecdXVHNDCPguByYWFPJNI8NgzzzB06NBU\nj5n1jIukrDF9+nTOOu00nq6s5Ojt2P9t4IdFRTw8dSonnnhissfLKcZFUlZYvHgxgw88kOcqKjii\nFc+bDZxSUsI7Cxaw5557Jmu8nOPVYpKywh0338wFNTWtCgvA4cD46mpuv+mmZIyVs1y5SGrzNmzY\nQJ8ePXi/qopdW9h+BvAHoALoBpwL/HyT7cuA/YqKWLJyJR06dEj+wDnAlYukNm/KE08wJB5vMSwA\nVwN/A9YBrwF3AtM22d4b+H48zuOPPZbcQXOIcZHU5r03Zw7fr6jY4vZ9gcJNfp8AejTb5/sVFbw/\nZ04SpstNxkVSm7d29Wo6bWOf8UAJTaH5BXBws+2dN76OwjAuktq84vbtqdrGPmXABmAmTXGZ22x7\n5cbXURjGRVKb16d/fz7ajs8NiwHHAD8CpjTb9pf8fHbdY4/ww+Uo4yKpzTvz7LN5Ih6ncjv3r6Pp\nENk3qoDH8vI4c8yY8MPlKOMiqc3r27cvgw8/nCdb2PYl8CRNlyE3ANOBp4FTN9nnKWDgwIF85zvf\nSfqsucK4SMoKP736av6jpISvmj0eA+6h6XLjUuAa4HfAoRu3rwauKynhp1dfnbJZc4E3UUrKGv92\n5ZX8YfJkXquspOt27F8OnFRczFEXXMCNkyYle7yc4spFUta47qabOPq88xhcXMzLNB0Ga0kD8Aow\nuLiYw84+m99MnJi6IXOEKxdJWeepp55i4rXXsmrJEsZWVzOosZEOwHpgbjzOvYWFdNt1V6649lpG\njx6d7nGzknGRlLXmzZvH/XfcwScffsi6devo0KED/ffdlwsuuYRDDz102y+gHWZcJEnBec5FkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWSFJxxkSQFZ1wkScEZF0lScMZFkhSc\ncZEkBWdcJEnBGRdJUnDGRZIUnHGRJAWXSPcA2aK2tpb333+fr7/+mng8TmlpKQcccAB5eXnpHk2S\nUs64fEtLly7lvrIyHpg8me6NjfSMx4mALxobqWjXjgsvuYRzL7iAHj16pHtUSUqZWBRFUbqHaIui\nKOKaCRMou/NOTo8ixtXUsE+zfeYDkwsLeQb45a9/zaVXXEEsFkvDtJKUWsZlB0RRxPlnnslHzz/P\nixUVdN/G/kuAE0tKOGXsWK6fODEVI0pSWnlCfwdc98tf8sFzzzFzO8IC0Ad4s6KCZ++5h3vLypI9\nniSlnSuXVlq9ejV79O7NX6qr2amVz/0LcEz79iz58ksKCwuTMZ4kZQRXLq3024ce4pR4vMWwLANO\nBkqBnYCfAg2bbN8HOAiYOnVq0ueUpHRy5dIKjY2N7LnLLjy+YgWDWth+GtAJuBcoB44DzqcpMt94\nAbhh332Z/eGHSZ9XktLFlUsrLFu2jMp16xi4he0fAf8MFAA9gRM2Prapk4D5H39MTU1N8gaVpDQz\nLq1QXl5OaSLBli4mPh54AqgC/g68Bgxrtk8C6FxQQHl5efIGlaQ0My6tkEgkqN/KUcRrgQ+BjsCu\nwKHAqS3sVx9F5OfnJ2NEScoIxqUVevTowfKaGmpb2BbRtHL5EVAJfAV8DUxott8aoLKhgU6dOiV1\nVklKJ+PSCt27d+eg/fbj+Ra2fUXTHfkXA/lAV+Bs4NVm+z0aizHihBNIJPzkHUnZy7i00vgJEyjr\n0GGzx7vRdPnxZJouP14DPAIcsMk+EVBWXMz4n/0sBZNKUvoYl1YaMWIEf00keLPZ4zHgWeAlmkLT\nH2gHTNpkn6eBRPfuHHXUUSmZVZLSxftcdsDrr7/OWSNGMLOqigHb+Zx3gVOLi3nlj39k4MAtXcws\nSdnBlcsOGDp0KLc9+CDHFhXxPNC4lX3rgUdpCstjzz5rWCTlBFcu38Ibb7zBZeefz/rlyxlXVcXo\nxkZ60HRu5Qvgsbw87s7Lg/x8XnrjDQ455JA0TyxJqeHK5Vs45phjWPDJJzw2cybvjxzJge3bUxSP\nUxKPc1jHjiw9/XRefvdduvbty/Lly9M9riSljCuXwL7569z0S8FeffVVrrjiCj744AMvQZaUE1y5\nBBaLxTb7tslhw4ax884788ADD6RpKklKLVcuKbJgwQJOOukkPvnkEzq0cJ+MJGUTVy4pcvDBB3Pc\nccdx0003pXsUSUo6Vy4ptGTJEg466CAWLlzILrvsku5xJClpjEuKXX311axatYoHH3ww3aNIUtIY\nlxRbu3Yte+65JzNmzGD//fdP9zhEUcTSpUv56quvaGxspGvXruy2227E4x4xlbTjjEsa3Hnnnbzy\nyitMmzYtbTNUVlYyZcoUym68kWXLlrFTfj7xWIxVdXUUd+rEuCuu4OxzzqFLly5pm1FS22Vc0qC2\ntpYBAwZw1113MXTo0JS//xOPP86/jh3LEbEY4zds4Dj+98qOCJhN06c3v9LYyM9/9SuumDBhs8ur\nJWlrjEuaPPvss/z7v/87CxYsIC8vL2Xve8ettzLxmmt4qbKSbR2UWwIMLy7mBz/5CRPvvtvASNpu\nHlhPk5EjR9KhQwd+97vfpew9pz79NDf/4hfM2o6wAPQBZlVWMvORR5h4443JHk9SFnHlkkZ/+tOf\nGDVqFJ988gnFxcVJfa/6+nq+06sXv1+9msNa+dwlwAGFhXz6xReeg5G0XVy5pNFhhx3G4MGDmTRp\n0rZ3/pZeeeUVdq2t3Sws7YEOm/xKAP/abJ8+wInxOI88/HDS55SUHVy5pNn//M//MGjQID766CN6\n9uyZtPc5fvBgzpo9m9O3sk8F0At4DTiy2bZ3gDE778yipUu9TFnSNhmXDHDZZZdRU1NDWVlZUl6/\noaGBwoIC1jc2UriV/R4Bfg0sbmFbBPQuLubtDz9k9913T8qckrKH/wuaAX7xi1/w9NNPs2jRoqS8\n/tq1aylJJLYaFmiKy1lb2BYDeiQSlJeXhx1OUlYyLhmgtLSUCRMmMGHChKS8fjwe3+pXMQN8DrwF\n/GQr+zRsfC1J2hZ/UmSIiy++mIULF/Lmm28Gf+2OHTtS09DAuq3s8zvgKKDvFrY3ACtqayktLQ0+\nn6TsY1wyRGFhIddffz1XXnkljY2NRFHEwoULefnll3nqqaeYPn06y5Yt26HXjsfjnDJ0KFu7o+ZR\ntr5qmQ706dOH3r1779AMknKLJ/QzSGNjI4cccgiHHHII8/74R8qXL2efRILixkbWxuPMq67mqCOO\nYPxVV3Hcccdt9yGq+fPnc+WVV/LZG2/wKU3nTzb1LjAUWAmUbOE1hrdvz2m3384555yzw38+SbnD\nuGSQWbNmcdqwYexfWcnPooih/OPSsgJ4Aihr3578Pn14YeZMdtpppxZfK4oiZsyYwY033shf//pX\nLrvsMu679VZuWraMk5vteyFQRdMJ/Za8D3y/pISlq1Yl/WZPSVkiUkaYNm1a1L24OHodomgbvxoh\n+nUiEe3Wo0e0dOnSf3idurq66PHHH48OOOCAaMCAAdGjjz4a1dbWRlEURW+99VbUvagoWrAd7/HN\nr88h2rWoKHri8cfT8dciqY0yLhngww8/jLqXlERvt+KHfgTRDYlEtH+/flFlZWW0YcOG6I477oj6\n9u0bDRkyJHrllVeixsbGzd7rmalTo+5FRdELGyO1tdd/B6LexcXR7RMnpuFvRVJb5mGxDHD6iBEc\n/OKLXNHK/xQRcHxJCUXHHsvs2bM5+uij+dnPfsagQYO2+ry33nqLsWecQby8nPEVFZwRRXTauK0K\neBoo69CBLwoKmHTPPfxw1Kgd+WNJymHGJc1WrVrFXn378rfqajo321YLjAP+AHwN9AN+A5ywyT6v\nAOM6d2bmnDnsueee2/2+URTx5ptvUnbzzbwwfToF8TjxWIyK+np+cNhhjJ8wgZNOOimlXwcgKXsk\n0j1ArnvwvvsYFYttFhaAepo+NPKtjf98Bfgx8AH/ez/KCUBeXR1r165t1fvGYjGOOeYYjjnmGKIo\nYv369TQ0NNCpUydvlJT0rflTJM1mvfYap1RVtbitGPgVTWEBOAnYHViwyT55wPC6Ot55550dniEW\ni9GxY0e6dOliWCQF4U+SNFuzZg1dt3PflcAnwL7NHu9aW0v511+HHUySvgXjkmb5+fnUbsd+dcDp\nwNlA8zMrtfE4Be3ahR5NknaYcUmzXrvswt+2sU8jcCZQCNzVwvbPioqS+l0wktRaxiXNRp93Hg91\n6LDF7RFwLvAl8AxN51g2tRqY1tDAiBEjkjajJLWWcUmzk08+mc/z83l/C9vHAYuAF4GWDnw9HI9z\n8vDhdOvWLWkzSlJrGZc0SyQSjP3pT7m2sHCz71z5HLiPps/26sX/fs/9lI3bvwLuKCxk/BVXpGxe\nSdoe3kSZAaqqqvjB4Ydz2F/+wi11dZt9anFL1gHDios5auxYbrj11mSPKEmt4solAxQVFfHiH/7A\n23vuyZmFhSzfxv7/BRxVXMwBo0dz/S23pGJESWoV45IhSktL+ePcuXQ54wz2KSrixyUl/AFYAawH\nltH0bZGHd+jAqaWljL3hBu5+4AFvepSUkTwsloHWrVvHww89xK8uv5yCkhIqa2vpVFzM/vvuy7ir\nrvIzvyRlPD9bLAN17NiRwUccQZ8BA1i4cGG6x5GkVvOYSoaaO3fuNj86X5IylXHJUHPnzmXgwIHp\nHkOSdohxyVBz5swxLpLaLE/oZ6A1a9bQu3dv1qxZQyLhaTFJbY8rlww0b948Dj74YMMiqc0yLhnI\nQ2KS2jrjkoG8UkxSW2dcMkwURa5cJLV5xiXDLF26lCiK6NOnT7pHkaQdZlwyzDeHxGKx7flsZEnK\nTMYlw3jzpKRsYFwyjOdbJGUDb6LMIPX19XTp0oUlS5bQpUuXdI8jSTvMlUsG+fjjj9l5550Ni6Q2\nz7hkEA+JScoWxiWDePOkpGxhXDKIV4pJyhae0M8QFRUV9OjRg6+//pp27dqlexxJ+lZcuWSIBQsW\nMGDAAMMiKSsYlwzhITFJ2cS4ZIg5c+Z4Ml9S1jAuGcKVi6RsYlwywMqVK1m7di177LFHukeRpCCM\nSwb4ZtU1i4utAAAEyklEQVQSj/ufQ1J28KdZBvCQmKRsY1wygHGRlG28iTLNGhsbKS0tZdGiRfTs\n2TPd40hSEK5c0mzx4sV06tTJsEjKKsYlzTwkJikbGZc085OQJWUj45JmfoeLpGzkCf00qqmpoWvX\nrqxatYqSkpJ0jyNJwbhySaOFCxfSv39/wyIp6xiXNPKQmKRslUj3ALlk3bp1zJkzh/LycvLy8njp\npZcYOXJkuseSpOCMSwp88MEHTJ40iSeffJL9Cwro3thIQyzG0g0buHb2bL5YsoQLxo+nd+/e6R5V\nkoLwhH4S1dXVMW7MGF579lnG1tVxXn09Ozfb5y/APe3a8Xgsxr9dey2XX3UVsVgsHeNKUjDGJUnq\n6uo47YQTYPZsnqyqYlun7JcCJxYXM3zsWH5z662pGFGSksYT+kly6YUX0vCnP/HcdoQFYFfgjcpK\nnr33Xu6/555kjydJSeXKJQk+++wzDtl7bz6trqZjK5+7EDihc2c+X7WK/Pz8ZIwnSUnnyiUJ7rv7\nbs5qbNxiWJ4E9gbaA3sAb2+ybX9gz4YGnn/++SRPKUnJ48olsJqaGvp0786s9evZs4XtM4Dzgd8D\nA4HlQAT/cKL/90DZwQfzxvz5SZ9XkpLBuAT24Ycf8qPBg/l4/foWtw+mKS5jtvIaVUDHeJza+nqv\nHJPUJnlYLLDy8nK6xlv+a20A5gOrgP40ncT/KVDdbL8iID8ep6KiIomTSlLyGJfACgoKqN3CtpVA\nHfAMTedZ3gP+C7iu2X4RUNfYSEFBQdLmlKRkMi6B9erViyW1tdS3sK1o4z9/CvQESoHLgVeb7bcU\n6FBYaFwktVnGJbC+ffvSr18/XmlhWxdgez7g5cG8PP5l9OjAk0lS6hiXJBg/YQJl7du3uG0McCfw\nJVAOTAJO3mR7HXB/QQHjLrss2WNKUtIYlyQYNWoU78XjzGph2zXAocCewD7APwE/32T7vbEY/ffe\nmwEDBqRgUklKDi9FTpLXX3+ds0aMYEZVFftt53NeAMZ27MisefPo379/MseTpKRy5ZIkQ4cO5faH\nHuLYoiKeghZP8H+jArg5Hmdcp068NHOmYZHU5hmXJPrn0aN5Zvp07thvP3YvLua6vDz+G/iapnMu\nC4BLCgroU1jIrCFDmDV/Poceemh6h5akADwsliLvvfcek2+9lRnTplG+YQN58TjdOnXih6efzgUX\nXUTfvn3TPaIkBWNcJEnBeVhMkhSccZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWS\nFJxxkSQFZ1wkScEZF0lScMZFkhSccZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWS\nFJxxkSQFZ1wkScEZF0lScMZFkhSccZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWS\nFJxxkSQFZ1wkScEZF0lScMZFkhSccZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWS\nFJxxkSQFZ1wkScEZF0lScMZFkhSccZEkBWdcJEnBGRdJUnDGRZIUnHGRJAVnXCRJwRkXSVJwxkWS\nFJxxkSQFZ1wkScEZF0lScMZFkhSccZEkBWdcJEnBGRdJUnDGRZIU3P8DrhmeZR+ERFcAAAAASUVO\nRK5CYII=\n"
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#Add isolated endpoints in the Graph\n#label of isolated end points\nlonely_ep = list(set(EPlabel) - set(selected_ep_labels))\n#lep:lonely ep \nfor lep in lonely_ep:\n    lep_row,lep_col= np.where(lab_ep==lep)[0][0], np.where(lab_ep==lep)[1][0]\n    node_index = node_index+1\n    G.add_node(node_index,kind='EP',row=lep_row,col=lep_col,label=lep)\n\n#print G.node[1]   \n#edges dict\nedges={}\nfor ed in EDlabel:\n    edges[ed] = np.where(pruned_edges==ed),np.sum(pruned_edges==ed)\n    #print edges[ed]\n#get nodes of kind EP\n\nep_nodes=[node for node,data in G.nodes(data=True) if data['kind'] == 'EP']\nprint ep_nodes\n#print edges.keys()\nep_edges_lab = (lab_ep>0)*edge_lab\n#print G.nodes()\n#",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[3, 4, 5, 6, 7, 8, 9, 10]\n"
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "figsize(5,8)\ntitle('Isolated endpoints are added')\nnx.draw(G)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAHnCAYAAADEl/S7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FOXhx/HP5s5u7gRCOFNFEVAr9QSsICKCnIpFoCJF\nLSCiohRUlB8elQoKKlbAAzkVDxSLKKioIKJcolAUSkWuAAZIICTZ3Hl+fyzE5MnmJJtg+b5fL16Q\nnZmdZzfhk9mZ2VmHMcYgIiJF/Op6ACIipxuFUUTEojCKiFgURhERi8IoImJRGEVELArjaaxjx47M\nmjWrVtb16KOPMmjQoFpZl5+fHz///PMp38/111/P/Pnza2BEvx0VfZ8SExP57LPPqnXfp7Ls/xqF\n0Udq4ofM4XDgcDgqvb7PP//8lNb1W/PRRx9VOua1+UvGlyr6PlXlZ6Yml/1fozD6SG3/kDkcDnSu\nftlq63uRn59fK+sR31IYa8FPP/1Ehw4diIqKol69evTv379o2tdff82ll15KVFQUl112Gd98843X\n+9i5cyedOnUiLi6OevXqccstt5CWlgbAoEGD2Lt3Lz179iQ8PJxnnnkGgLVr19KuXTuio6O56KKL\nWLVqVdH97dq1iw4dOhAREUGXLl04cuRIuY9h6dKlXHTRRURHR9O+fXv+/e9/F01LTExkypQp/P73\nvycqKor+/fuTk5NTNP3pp5+mYcOGNG7cmNdee63E/f7lL39h+PDhdOnShYiICDp27MjevXsr9fwU\n3wqcM2cOV155JWPGjCEmJoazzjqL5cuXA/Dwww+zevVqRo4cSXh4OPfccw8A9913H/Hx8URGRnLh\nhRfyww8/eH3ss2fPplWrVkRERHD22Wfz8ssvF01buXIljRs3ZvLkySQkJHD77bdjjOGpp56iefPm\nxMXFcfPNN3P06FGv933s2DF69OhB/fr1iYmJoWfPnuzfv7/S36f58+fTrFkz4uLimDhxYolpFY2j\nvGXPeEZ8IjEx0Xz22WfGGGP69+9vJk6caIwxJicnx6xZs8YYY0xKSoqJiooyCxYsMAUFBWbhwoUm\nOjrapKamGmOM6dixo5k1a5YxxpiffvrJrFixwuTm5prDhw+bq666yowaNcrr+owxJikpycTGxppl\ny5YZY4z59NNPTWxsrDly5IgxxpgrrrjCjB492uTm5povv/zShIeHm0GDBnl9LJs2bTL169c369ev\nN4WFhWbu3LkmMTHR5ObmFq378ssvNwcPHjSpqammZcuWZubMmcYYY5YtW2bi4+PNDz/8YDIzM82A\nAQOMw+EwO3fuNMYYM3jwYBMeHm5Wr15tcnJyzL333muuvPLKKj8/s2fPNoGBgebVV181hYWFZsaM\nGaZhw4ZFj6H4vMYYs3z5cnPxxRebtLQ0Y4wx27dvNwcPHvT6+D/88EPz888/G2OMWbVqlXE6nWbT\npk3GGGO++OILExAQYB588EGTm5trsrKyzHPPPWfatm1r9u/fb3Jzc82wYcPMgAEDvN53SkqKee+9\n90xWVpZJT083f/rTn0yfPn2Kppf3ffrhhx9MWFhY0XN3//33m4CAgKKfg/LGUdGyZzqF0UeKh+rW\nW281Q4cONUlJSSXmmTdvnrn88stL3Na2bVszZ84cY0zp/8zFLV682LRp08br+owx5qmnnioVuuuu\nu87MnTvX7NmzxwQEBBi32100beDAgeaWW27xuq7hw4eb8ePHl7itRYsW5ssvvyxa9+uvv140bezY\nsWb48OHGGGOGDBliHnrooaJpO3bsKBXG4tHIyMgw/v7+Zt++fVV6fmbPnm2aN29eNF9mZqZxOBwm\nOTm5aN5XX321aPrnn39uzj33XLN27VpTUFDg9XGXpU+fPub55583xnjCGBQUZHJycoqmt2zZssT3\n4sCBAyYwMLBS6/nuu+9MdHS0McaU+X06+X197LHHSjx3mZmZJigoqGjdZY0jPz+/wmXPdHopXQsm\nT56MMYbLLruM888/n9mzZwNw4MABmjZtWmLeZs2aceDAgVL3kZycTP/+/WncuDGRkZEMGjSIlJSU\nMte5Z88e3nnnHaKjo4v+rFmzhl9++YUDBw4QHR1NaGhoifWWd19TpkwpcV9JSUklxtmgQYOif4eG\nhpKZmQnAwYMHadKkSdE0+/E6HA4aN25c9LXL5SImJoYDBw5w8ODBSj8/9hicTicAGRkZJdZ10tVX\nX83IkSO56667iI+PZ9iwYaSnp3u932XLlnHFFVcQGxtLdHQ0H330UYnnvl69egQFBRV9vXv3bm64\n4Yai56pVq1YEBASQnJxc6r7dbjfDhg0jMTGRyMhIOnToQFpaGsaYCr9PBw4cKPHcOZ1OYmNjKzWO\ngwcPlrvsmU5hrAXx8fG8/PLL7N+/n5deeokRI0awc+dOGjVqxJ49e0rMu2fPHho1alTqPsaNG4e/\nvz9bt24lLS2N+fPnU1hYWDTdPrjQtGlTBg0axNGjR4v+pKenM3bsWBISEjh69Chut7vEess6QNG0\naVMefvjhEveVkZHBzTffXOFjT0hIKLHPsPi/wbMfbN++fUVfZ2RkkJqaSqNGjWjYsGGln5+KeHts\nd999Nxs3buTHH39kx44dPP3006XmycnJoW/fvowdO5ZDhw5x9OhRrr/++hIHurw998uXLy/xfLnd\nbhISEkrd/5QpU9ixYwfr168nLS2NVatWYTyv5Mr8Pp3UsGHDEs+d2+0uEeyyxtGwYUMSEhLKXfZM\npzDWgnfeeYekpCQAoqKicDgc+Pv7061bN3bs2MHChQvJz8/nrbfeYvv27fTo0aPUfWRkZOByuYiI\niGD//v2l/hPHx8ezc+fOoq9vueUWPvjgAz755BMKCgrIzs5m5cqV7N+/n2bNmnHJJZcwYcIE8vLy\n+Oqrr1i6dGmZ4//rX//KzJkzWb9+PcYYMjMz+fDDD0tsjdlOhqNfv37MmTOHbdu24Xa7eeyxx0rN\n+9FHH7FmzRpyc3MZP348bdu2pVGjRlV6fipiPz8bN25k3bp15OXl4XQ6CQkJwd/fv9Ryubm55Obm\nEhcXh5+fH8uWLeOTTz4pd13Dhw9n3LhxRb8EDh8+zJIlS7zOm5GRQWhoKJGRkaSmppZ4fir6PvXt\n25elS5cWPXf/93//V+KXZXnjuOmmm8pd9kynMNaCjRs3csUVVxAeHk7v3r2ZNm0aiYmJxMbGsnTp\nUqZMmUJcXBzPPPMMS5cuJSYmptR9TJgwgU2bNhEZGUnPnj3p27dviS2Vhx56iL///e9ER0czdepU\nGjduzL/+9S8mTpxI/fr1adq0KVOmTCn64X/jjTdYt24dMTExPP744wwePLjM8V988cW88sorjBw5\nkpiYGM455xzmzZtX5hZm8VOVunbtyqhRo+jUqRPnnnsu11xzTYnlHA4HAwcO5LHHHiM2NpbvvvuO\nBQsWAFTp+fF2elTxr++9914WLVpETEwMo0aN4vjx4wwdOpSYmBgSExOJi4tjzJgxpe43PDycadOm\n0a9fP2JiYli4cCG9e/cucz0n19WrV6+iI+1t27Zl/fr1Xp+rUaNGkZWVRVxcHO3ataNbt24l7q+8\n71Pr1q158cUXGThwIA0bNiQmJqbEbovyxtGqVatylz3TOYzRyW9Sd4YMGULjxo154okn6nooIkW0\nxSh1Sr+X5XSkMEqd0tvQ5HSkl9IiIhZtMYqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaF\nUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhF\nRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTE\nojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwK\no4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUeqE2+1m\n1qxZXNGqFTEuF0EBAcS6XLRt3ZrZs2eTlZVV10OUM5jDGGPqehBy5sjLy2PCgw/y8ksv0d7h4M6M\nDC4BwoF0YB0wIyyMtcZw5113MeHJJwkICKjbQcsZR2GUWuN2u7mxa1f8N25kelYWzcqZdxcw3Okk\n8PLLeefDDwkNDa2tYYoojFI7CgoKuLFrV8K/+oo52dlUZhswD/hzSAjmmmt4a8kS/Py050dqh37S\npFbMnz+fw998w2uVjCJAIDA/O5s9K1fy1ltv+XJ4IiUojFIrpk+axLjMTF4GLgFCgCHWPJ8B5wEu\noBOwFwgGHszMZPqkSbU5XDnDKYzicxs2bODQ3r10AxoB44HbrHmOAH2BJ4GjeOJ584lpvYBd//0v\nW7ZsqaURy5lOYRSfe/v11/lLdjb+wA1AbyDWmuc94Hw8cQwCHgU2AzuAAODWnBzeeuONWhuznNkU\nRvG55L17SSwsLHGbfcTvB+D3xb52As2BrSe+Tiwo4NC+fT4bo0hxCqP4XF5uLoHWbQ7r60wgwrot\nAsg48e9AIDcnxwejEylNYRSfi65fnxTrNnuLMQw4bt2WhufEb4CUE/cjUhsURvG5q7p04V/h4SVu\ns7cYW+PZp3hSJrDzxO0A/woP56rOnX02RpHiFEbxuRtvvJEfHA62AQVANpB/4t85J/6+Ac/+xPdO\nTH8MuAg4F9gC7A4MpFevXnUwejkTKYzic0FBQdwxfDgvBAXxBJ4DK5OABUAonlN04oB3gYeBGGAj\n8OaJ5V8IDmboyJF6z7TUGr0lUGpFcnIyl7RuzdMpKfSvwnLzHQ7Gx8Wx8ccfiYuL89n4RIrTFqPU\nivj4eD78/HPuDQ9ntsPew+jdS35+jI2I4KOVKxVFqVXaYpRatX37dnp17ky948cZkZ7OTXje9ndS\nNvA2MD08nGNRUSz9/HOaN29eN4OVM5bCKLUuPz+fDz74gOmTJrF582YuCA4m3BjSHQ625ORw8R/+\nwIgHHqB79+74+/vX9XDlDKQwSp3atWsXO3fuJD09nfDwcJo3b05iYmJdD0vOcAqjiIhFB19ERCwK\no4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCK\niFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iI\nRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgU\nRhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEU\nEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhER\ni8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEo\njCIiFoVRRMSiMIqIWBRGERFLQF0PQETOXAUFBSxfvpz3Fizg8IEDFBQUEFOvHp379KFfv36EhobW\nybgcxhhTJ2sWkTOW2+3mheeeY+ZzzxGblcXgjAyaAv7AIWBRWBgbgMFDhjD6oYdISEio1fEpjCJS\nqw4dOkTPTp1I+PlnHs7K4tIy5tsJPB8YyPsREXz4xRdccMEFtTZGhVFEas3x48e56uKL6bZnDxPz\n8nBUYpk3gNGRkazeuJHmzZv7eoiAwlhl2dnZLFq0iOXvvkvqoUP4+fkRl5BA7wED6NmzJwEB2m0r\nv03GGL755hvWrVvHsWPHCA4OpmHDhvTp04eoqKgaWcetN91EyNKlvJSTU6konvSSnx8vJiay+aef\ncDiqsmT1KIyVdPjwYaY89RSvvfIKbYB+6enEAwbYDywID2dPYCBDR45k1N/+Rnh4eN0OWKSSMjMz\neX3BAqZPnow7OZlueXlE5+aS7e/PT6GhfJafz019+zJi9GjatGlT7fUcPHiQVr/7Hbtzcois4rIG\nuCgsjCmLF9O5c+dqj6HyK5QKbd++3ZzVoIG5MyjI7ABjyvjzPZj+wcHmwrPPNklJSXU9bJEK7dy5\n07Ro0sT0dLnMJ2AKvPxc/wLmSX9/k+B0mr9PmGAKCwurta7HJ0wwQ0NCyvz/Y8DsABMM5hYv02aC\nuaFLlxp+BrxTGCuwe/du0ygmxsxyOMr9hp78UwjmyYAA06JJE5OSklLXwxcp065du0yjmBjzTz+/\nSv1s7wdzkdNpHho9ulrrS6xXz2yqYB3XgvkjmEFepqWDiQoONocOHarhZ6I0neBdDmMM/bp35760\nNG6r5B4HBzAuP58eycncMWCAbwcoUk15eXl0v/pqxh47xl2FhZVapiHwidvNOzNm8Mbrr1dpfQUF\nBew9coQLy5nnTSAauAbPS2dbGJAYHMzevXurtO7qUBjLsXbtWlJ37+a+goIqL/t4bi6rv/ySXbt2\n+WBkIqfm/fffJy4lhXsqGcWT6gEvu91MfPhhTBUOT2RlZRHs749/GdOPAxOAZ/EexZPCHA4yMzMr\nvd7qUhjLMf2ZZ7jT7fb6JG0DOgFRwDnA+9Z0JzC4sJCXXnjBx6MUqbrpkyZxV3q612m3AAlABHAW\n8KQ1vSNQmJLC6tWrK70+p9NJbmEh+WVMHw/cgWertLxjzseNISIiotLrrTafv1j/jTp+/LiJCA42\nKV72deSBOQfMsyf2KX4OxnVix3Hx+f4LJi4srNo7q0V8Ydu2baZBaKjJKWM/31YwWSf+vR1MPJhl\n1jwvOBymf8+eVVpvq6ZNzUov6/sOTGswuSe+nlDGwZfDYKJCQszRo0d99Mz8SluMZThw4ADxQUHE\neJm2HTgIjMLz2+1qoD0w35qvOeDOziYjI8OnYxWpii1btnBlQABBZUxvDYQU+zoAqG/Nc7UxbN60\nqUrrHXb//cxwuUrdvgrYDTTFs6U6BXgXuMSa7zU/P/r07l1j51SWR2EsQ0ZGBq4qnEhaCGz1cntY\nYGCt7BMRqay0tDQiKthvPgJw4YnkI8AfrOmRQJr1c52RkcHWrVtZs2YN33//PUeOHCkx/dbBg/m4\noIBfrPsaCvwMbAa+B4YD3YGPi81TCMwMDWXE6NEVP8AaoDCWITIykuNl7Jhugec36NNAHvAJ8CWQ\n5WXe43l5tbNPRKSSnE4nbr/y/+tPBzKAFXjCuN6angk4g4MB2Lx5M8MGD6ZJ/fr0a9eOv3XvzuAO\nHWjeuDE3XncdK1aswBhDVFQUt912G7eFhJTY1xiK5/9TfSAez9HnUCC22DyPBAaS2Lo1l15a1jur\na5be+VKGrKwsGsXGsjkriyZepv8buBvPVuKlQByelx+vFJvnW6BvvXrsSk6ulbcxiVTG6tWruat7\ndzanp1fqbXl34vnZfrbYbe8B09q0wRURweYNGxiWk8MdBQUUvwZOOvA6MD0sDFOvHu8uX87Bgwfp\nfe21dAHm5+URXMG6DfBkQADz4uP5+vvviYuLq8IjrT5tMZYhNDSUP//5z7xSxnufLwBWAkeAZXiu\nBHKZNc+MkBCG33uvoiinlfbt2+N2uVhbyfnz8LysLm6my8VPSUmcs3Ytu91uxltRBAjH87J4c0YG\n9+zeTbs2bejTpw/zFi2i8NpraedysRDI9bLOQjxbqz2cTt476yxWbdhQa1EEdFS6PD/88EOZR++2\nnDhylwnmaTBnFTuqZsCknjiClpycXNcPQ6SUZyZPNreEhpb6uT4EZiGYDDD5YJaDiQCz3nrbnsvh\nMH8LCKjUO2ZO/nkTTIOICJOUlGQKCgrMu+++a6657DITHxpq7g0IME+DmQrmIYfDnBsWZi783e/M\nzBkzTGZmZq0/P3opXYFe11xDszVrmGZdDWQs8Cqe36ZXAS/gOecLPL/tbg4NJa5fP2bMmVOr4xWp\njNTUVM5t2pQlmZm0K3b7EeAmPAdCDHAunn2MvU5MLwB6+PvzY0EBu6j6S84xAQHk3n47z8+cWXTb\n9u3bWbx4MUcOHiQ/L4/o+HiuvfZa2rVrV2evthTGCqSlpXHlH/5A96QkJubmVviDkAcMDQnh59at\n+WTNGoKDK9qLIlI3li1bxm033cQyt5uLKjF/ATAiOJj3/fyYnpVFX2v6P4E5ePa7DwBme7mPvUAb\np5O9hw7h8nLqzulC+xgrEBkZyRfr1vFly5Zc63LxIZ4fEFsu8DZwpcvFkXbt+HDlSkVRTmvdunXj\nn3Pn0sXp5EU8B0vKshHoERrKlnPPJQjo42WeRnjewXJbOffTFPijw1Hl91rXNm0xVlJOTg4LFy7k\nxaee4si+ffwpO5v4wkLP9RgDAngzMJAWLVty59ix9O3bVxesld+Mb7/9lokPP8wXq1YxwBiuz8kh\nGsjGc1DxlfBwDoWEcOeoUeQbw5FHH+XZ/LLe3OeJYxLetxjBc0XuxdddxzvLl9f0Q6kxCmM1bNiw\ngY8//piUgwfx8/cnLiGBXr160bp167oemki1JSUl8cqMGaz/4guOHT1KcEgICU2a8OehQ+nWrRv+\n/v48MHo0UVOn8lA59/MInos3lxXGT4FJl1zCig0bavwx1BRt1lTDpZdeWmsnmorUlsaNG/PYk/Yl\nI0ryDwykouvxVHS4pADwq+AE87p2eo9ORE4r9eLjSQoq613WHhW9BN0PxDVoUGNj8gWFUUQqrVev\nXizy8/P69teTKtpinBcezo2DBtXksGqcwigilXb22WdzySWX8LaXaQV4Dtjkn/h3DqXP4NgK/OTv\nT+/evX070FOkMIpIlYwYO5YpLhfZ1u1P4LlA8yRgAZ4LQRTfY2mAfwQHc8eIEQQGBtbOYKtJR6VF\npEoKCwsZ0KcPhZ9+yhvZ2VQ2cU8GBPBOYiJffvvtaX/FKYVRRKosOzubG7t2xWzYwBy3m/hy5s0A\nxgcF8XGDBnz69dc0atSotoZZbXopLSJVFhISwpIVKzj/tts4LySEAU4nq6HoVB4D/AjcExxMs5AQ\n9nfuzJrvv/9NRBG0xSgip+jYsWPMmzuXmVOm8J+kJMICA3Hn5eEKCODu0aMZOmIETZp4u6rp6Uth\nFJEak5+fT0ZGBrt27aJ///785z//qeshVYvCKCI1Lu/ER3qkpKTgdDrrejhVpn2MIlLjAgMDOe+8\n89i61dtHxJ3+tMUoIjWusLCQrl27Uq9ePc4//3zCwsI4++yz6dKly2/iylMKo4jUmNTUVOa89hoz\npk4lMDWV9vn5xACZgYF8GxhIUmAgw+6+mzuGD6fBafx+aYVRRGrEN998w43dutE5L48RbjdXUPp9\n09/h+ZC4xf7+zF+0iK5du9bBSCumMIrIKVuzZg19unRhrtvN9ZWZH+jrdPLKm2/Ss2dPXw+vyhRG\nETkl+/fv5+JWrZh7/DjXVWG5DUB3l4uV69fTqlUrXw2vWnRUWkROyfTnn6dfdnaVoghwKXBfVhZT\nK7g4bl3QFqOIVFtOTg7N6tdn1fHjtPAyvSOwjl8/KqAxsK3Y9ENAi5AQfj5wgOjoaN8Otgq0xSgi\n1bZ48WLON8ZrFMFz8OXkJxCmUzKKAPWB7n5+zJ1d1ifE1A2FUUSqbdP69XROL++DVyv+qIPObjeb\nvvqq5gZVAxRGEam2tMOHiapgnoeAesCVwCov0yOBtKNHa3pop0RhFJFqCw0PL/fzXyYBu4ADwFCg\nJ/CzNU8WEHqavZ9aYRSRamty1ln8GBxc5vTLABcQCNwKtAc+subZ7u9Pk+bNfTbG6lAYRaTaBgwc\nyCLgeDWXzwNmBQVx6x131OCoTp3CKCLV1rBhQzp36sR8L9PSgI/59ZMDXwdWA8XfBLgEOLtFCy64\n4AKfj7UqFEYROSV3P/ggk51Okq3b84DxeE7JqYfntJ1/ASdfNKcDj7pc3D1uXK2NtbJ0greInLL/\ne+ghlk2bxnK3m9hKzJ8B3OB0ctZNNzFzzhwcDvtyE3VLW4wicsoemziRa4YNo53TyTJ+/VAsm8Fz\nyk4Hl4umffrw4qxZp10UQVuMIlKD3n77bSaPH0/q/v0Mz8qibWEhEXi2EL9zOJjhckFMDPc98gi3\n33HHaRlFUBhFxAc2bNjAS889x7bNmzmenk54WBhnnXsufx01iquuuuq0DeJJCqOIiEX7GEVELAqj\niIhFYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqI\nWBRGERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhF\nYRQRsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRG\nERGLwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQR\nsSiMIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGL\nwigiYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiM\nIiIWhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigi\nYlEYRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIW\nhVFExKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEY\nRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFE\nxKIwiohYFEaplPz8fFJTU8nOzq7roYj4nMIoZTp27BjTnn+eVk2bEhIUxFkJCUS4XDSIjGTcmDHs\n2bOnroe/i/GSAAAMQklEQVQo4hMKo5SSn5/P2Hvu4XcJCXwzbhwv7dtHrjEcy80lp7CQVcePkz1t\nGhefdx59u3YlJSWlrocsUqMcxhhT14OQ00dOTg59u3Ujf9065rjdNChnXjcwITCQJfHxrPj6a5o0\naVJbwxTxKYVRihhjGHjDDeR98gkLs7IIrORyk/z9WdCsGV9t2kRkZKRPxyhSG/RSWoosW7aMf69Y\nwfwqRBFgbEEBf9i/n6f/8Q+fjU2kNmmLUYr06NiRG1et4jbr9jDAUezrLGAEMK3YbduATpGR7Dl0\niKCgIB+PVMS3tMUoAOzatYu169bR38u0DCD9xJ9fgFCgnzVPS6BlQQGLFy/27UBFaoHCKAAsWbKE\nGwFnBfMtAuKBK71MuzUjg/fmzavxsYnUNoVRADhy+DCNKnHy9lzg1jKmNQKOJCfX5LBE6oTCKAAU\nFhTgX8E8e4AvgcFlTPcHCgsLa3RcInVBYRQAYuLiOFTBQZP5wB+BZmVMPwREx8TU8MhEap/CKAB0\n7tyZ9wICyC9nnnmUvbUI8I7LxXU33VTDIxOpfTpdR4pc+fvfM3rLFm7wMu1roAuQDLi8TE8Cfu90\nsic5mbCwMF8OU8TntMUoRUY88ADPulx420s4D+iL9ygCvBAQwMCBAxVF+Z+gLUYpkpubyzVt23LF\n1q1Mzs0tcVJ3eRYBo2JiWLt5M40bN/blEEVqhbYYpUhQUBDvf/IJHzdpwsjgYPIqmN8ALwMjIyL4\nYMUKRVH+ZyiMUkJsbCyrv/2Wfe3bkxgayqMBAey35kkHZgIXulyM8fdn0rRptGnTht27dzNuzBi6\ntWtHu9atueaSS/hLv358/vnn6IWJ/JbopbSUaevWrcx49lneWLiQBgEBRPr5kWkM+3JyuKZjR0aM\nGUNubi5DhgyhzTnnsH7jRm4tLKRzbi5ReN5T/SPwclgYBVFRjHzgAe4cMQI/P/0+ltObwigVyszM\nZM+ePaSlpeFyuWjUqBGxsbEAzHntNe4fOpRJBQUMxPvBGQOsBsY5nTTo2JH5ixYRGhrqdV2FhYX8\n8ssvHDt2jMDAQOrVq0dUVJSvHpqIVwqjVNubCxcy5vbb+TQri/MqMX8OcEtICIUdO/L20qX4+//6\nXpuUlBRmz5rFjKlTyTh+nJiAAHKN4VBuLn+8/HJGjB1Lt27dSiwj4isKo1TL3r17+UPLlnzhdnNB\nFZbLAa51OrnxiScYdf/9FBYW8ui4cUx7/nl6+/kxwu3mMn69zFkW8DbwYng4h0NDWfDee7Rv376m\nH45ICQqjVMvDY8eSOW0az+XkVHnZtcCf4+P5T1IStw0YwM5ly3g3M7Pcj1EA+BAY4nTy2ltv0aNH\nj+oMW6RSFEapspycHJrVr8+q48dpYU3LBe4EPgNSgbOBfwBdi81jgIvDwmjeqRPJn33G8sxMvO9x\nLG090N3p5KOVK7n00ktP8ZGIeKfDg1Jly5Yto6UxpaIIkA80xXMVnuPA3/Fc1Lb4B606gBszMvhk\n6VIWVyGKAJcBU9xu/jZsWDVHL1IxhVGqbM+ePZyfm+t1mhOYgCeOAN2B3wGbrPm2AUMKC/F2LZ4k\noCcQCyQAdwMFxaYPAH7avp2tW7dW9yGIlEthlCpzu90488u7Ds+vkoEdQOtit2UBy4CRZSxzDxAH\nHAS+B1YB04tNDwT+mpfHjGefrdrARSpJYZQqi4yM5HglPvAqD/gz8Bfg3GK3r8ez7/HsMpb7AbgZ\nCMLzMQpdT9xW3ID8fD5csqRK4xapLIVRquzCCy9kpb8/5R21KwQGASHAP61pKXi2CMtyHfAGni3L\n/Xi2LrtZ88QDqenpVRm2SKUpjFJl7du3JyA2lpVlTDfA7cBh4F0o9ZEJBXiOXpflUWArEAE0AS4F\nentZh8NR2ev/iFSNwihV5nA4GDF2LC+6vF+d8U5gO7AECPYyfRvwSxnvlzZ4thj/BLiBI3hO+3nA\nmu8gEBseXo3Ri1RMYZRquWXQIDY4nbxhbbXtwXMpss1AAyD8xJ+FJ6YnAa+EhvJLcDD/8XK/R4Bv\n8RyYCQRi8Oyj/Mia7/WAAHrrYxTER3SCt1Tb1q1buaZdO55PT6d/JebfBXRzOrn9kUc4mppK1rRp\nPGud9mOAxsC9wGg8lzgbgufiFAtOzJMDNAsNZdV339GihbezKUVOjbYYpdrOP/98VqxZw9i4OAY6\nnawBrwdk9gOPBgRwRWgoIydOZMxDDzFs5Ejm+flhfwq1A3gP+ADPAZpz8LwcL35izlzg/AsuUBTF\nZ7TFKKfs2LFjzHntNR5/8EESAgLomp9PZF4eWX5+bHc6WVVQwICBAxlx3320bv3rGY2PPvwwHz73\nHJ+73VR2b+FKoJ/LxYqvv+bCCy/0xcMRURilZuzcuZO2bdsyf/58tmzZQtqxY4Q6nTRp0oQbbriB\ncC8HSowxjLzjDr5+6y3ey8zkd+XcvwHeAu5xuXhzyRI6derkq4ciojBKzXj88cc5cuQI06ZNq9Jy\nxhimTp7MxMcf548OByMyM+nMr/t40oB5DgfTXS4C6tVj3rvv0qZNm5oevkgJCqOcMmMMLVq0YP78\n+Vx++eXVuo/MzEzeeP11Xpw0iR179xIdFERuYSGZBQX07tqVEWPGcOWVV+rcRakVCqNUmdvt5ujR\no/j5+RETE8OWLVsYOHAgO3bsqJFwZWRkcPToUYKCgoiOjiaoEm8/FKlJAXU9APltyM7OZtGiRUyf\nNIlN27YRExxMgTGk5eaS2KABV1x9NcaYGgljWFgYYWFhNTBqkerRFqNUaPasWTxw331cZAx3ZWTQ\nnV9/o2YCbwIvuFykh4fz6htvcPXVV9fdYEVqgMIo5Xpi/HjmTp3K4go+28UAHwODQ0P555w5/Klf\nv1oaoUjNUxilTC/NmMGUv/2Nr9xu6ldymS1A59BQ3lm2jA4dOvhyeCI+ozCKV2lpaSQmJLA+K4tz\nqrjsEmBcs2b8e9cuHUWW3yS9JVC8mjd3Ll0djjKj+CbQEggDmgNfFZvWEyhISWH16tU+HqWIb2iL\nUUoxxtCqaVNeTkrij16mfwr8Fc/nPV+G5xJgBmhYbJ4XHA7WdO/Omx984PsBi9QwhVFK+e9//8s1\nbdqwJzMTby+E2+EJ45By7iMVaBwYiLuMD80SOZ3ppbSUcuTIERICArxGsQDP9RIP4bnyTRM8n+KX\nbc0XDRQWFuJ2u306VhFfUBjFq7IOmSTj+ZCrd/HsV/we+A7P50fbyzscDvSCRH6LFEYpJSYmhuT8\nfK/XVgw98ffdeD6QKha4n9JX2E7Ds9/R6XT6bJwivqIwSinnnHMO/uHhrPcyLRrPFbYr8iZw/dVX\n63Qd+U1SGKUUPz8/7rz/fl4MDfU6fQjwAp5PATyK5+raPYtNN8D0sDBGjB3r66GK+ISOSotXKSkp\nNG/cmM3Z2TS1puXj+UyWN/B8bvTNwGTg5DVwPgVGNmzItn378Cvj0wBFTmf6qRWvYmNjGf/YY3R3\nOjlqTQsAXsSztXgQeI5fo7gDuDU0lOdnzVIU5TdLP7lSpvvGjKHr7bfzR6eTnysx/9dAh9BQ/vHP\nf9K1a1dfD0/EZ/RSWir0/NSpPPbII3Ty82NEZiZX8+vpPHnA+8D08HD+4+/Py/Pn06NHj7obrEgN\nUBilUtLT01kwfz7Tn36ag7/8QlxgIAXAoZwc/nDBBYx44AH69OlDYGBgXQ9V5JQpjFIlxhiSk5NJ\nTU3F39+fuLg4YmNj63pYIjVKYRQRsejgi4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEY\nRUQsCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFE\nxKIwiohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQs\nCqOIiEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIw\niohYFEYREYvCKCJiURhFRCwKo4iIRWEUEbEojCIiFoVRRMSiMIqIWBRGERGLwigiYlEYRUQsCqOI\niEVhFBGxKIwiIhaFUUTEojCKiFgURhERi8IoImJRGEVELAqjiIhFYRQRsSiMIiIWhVFExKIwiohY\nFEYREYvCKCJiURhFRCz/D62FPockhoIOAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 70
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": "each endpoints pairs are linked"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# get all nodes of kind EP\n#get their label\nnode_to_label={}\nlabel_to_node={}\nfor node,data in G.nodes(data=True):\n    if data['kind']=='EP':\n        #print node,data['label']\n        label_to_node[data['label']]=node\nprint label_to_node",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "{1: 3, 2: 6, 3: 4, 4: 5, 5: 7, 6: 8, 7: 10, 8: 9}\n"
      }
     ],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#ep_labels from edge lab -> ep labels directly!\n#\n# node_index is not endpoint label!\n#\nmap_epEdlab_epPair={}\nprint EDlabel\n##plt.figsize(14,8)\nfor elab in EDlabel:\n    #subplot(3,3,elab)\n    mask = (ep_edges_lab==elab)# ep lab=edges lab\n    nodes_pair = np.unique(mask*lab_ep)[1:]\n    pair_weight = np.sum(pruned_edges==elab)\n    edge_image = (pruned_edges==elab)*elab\n    print 'edge lab',elab, 'ep2 lab',nodes_pair,'weight',pair_weight,'node1:',nodes_pair[0]\n    ##title(str(elab)+str(np.where(ep_edgelab==elab))+str(nodes_pair))\n    ##imshow(mask, interpolation='nearest')\n    node1 = label_to_node[nodes_pair[0]]\n    node2 = label_to_node[nodes_pair[1]]\n    G.add_edge(node1, node2, weight = pair_weight, image=edge_image)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[1 2 3 4]\nedge lab 1 ep2 lab [1 2] weight 45 node1: 1\nedge lab 2 ep2 lab [3 7] weight 11 node1: 3\nedge lab 3 ep2 lab [4 5] weight 10 node1: 4\nedge lab 4 ep2 lab [6 8] weight 10 node1: 6\n"
      }
     ],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "plt.figsize(5,8)\npos=nx.spring_layout(G)\n#nx.draw_circular(G)\n#subplot(121)\n#nx.draw(G, width=range(1,4))\n#imshow(edge_lab,interpolation='nearest')\n\n#subplot(122)\nedge_labels=dict([((u,v,),d['weight']) for u,v,d in G.edges(data=True)])\nnx.draw(G)\n#nx.draw_networkx_edge_labels(G,pos,edge_labels=edge_labels)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAJeCAYAAAB4aUjsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FPXCxfHvpmfTiIEQeuiEUESqVBEBFUGqYC9gQ8Vr\nRe/F3sGuiKIURUCkSxFUUBEBC0VaKIIECSG0hIT0su8fAV8MARKY2dlyPs/jo+xMssf13j37Kztj\nczgcDkRERAzkY3UAERHxPCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGcykVERAynchEREcOp\nXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGc\nykVERAynchEREcOpXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHD\nqVxERMRwKhcRETGcykVERAynchEREcOpXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXEREx\nnMpFREQMp3IRERHDqVxERMRwKhcRETGcykVERAynchEREcOpXERExHAqFxERMZzKRUREDKdyERER\nw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGcn9UBPMX69etZsmQJR1NS8PH1Japy\nZXr37k1cXJzV0UREnM7mcDgcVodwV3l5ecyYMYMPXnuNpL/+YmBeHpULCnAA+/39+dLPj/gmTRg+\nciR9+/bF19fX6sgiIk6hcjlPR48epW+PHti2bePhzEx6cfowMA+YA7wRGkqVtm2ZPn8+ISEhzg8r\nIuJkKpfzkJ6eTudWrbg8MZHX8/LOuXCVB9wVFMSepk1Z+tNPBAYGOiOmiIhltKB/Hm6/7jou3buX\nN8pQLAABwMScHKI2b+ahe+4xO56IiOU0cimn7du306VFCxKzsyk5/kgA7gPWAZWAMUDfU46nAnWC\ngtiemEh0dLRzAouIWEAjl3Ia9847DM3PP61YCoBrgT4Ul8h44CZg5ynnRAL9gQnjxzsjqoiIZTRy\nKYecnByqRUWxLiuLWiWObQYuBTJOeawn0BZ4/pTH1gIDKlXir5QUbDabuYFFRCyikUs5JCUlEeHj\nc1qxnEkRxaVzqpbAwdRUsrKyjA0nIuJCVC7lkJ6eTrhP6S9ZQyCa4nWWfOAbYAWQXcq5Ef7+HDt2\nzKyYIiKWU7mUQ0hICJlnmEX0B+YBi4AqwFvAdUD1Us49XlBAaGioWTFFRCyny7+UQ5UqVTiYl8dR\n4KJSjjcFfjjlz+2B20ucsxsICAggLCzMnJAiIi5AI5dyCAsLo0+vXkw+w0L8JiAHyAJeB1KA20qc\n86G/P7ffcYcW80XEo6lcymn4o48yzm6nqJRjU4CqQGXge+BbiqfLTsoGJvn6cs+IEeYHFRGxkLYi\nl5PD4aBdkyYM2b6dhwoLy/WzT/j7s6VDBxZ8/71J6UREXIPK5TwkJibSoUULXkxL47YyvHwO4HVf\nXz6OiWHVhg1UrFjR/JAiIhbStNh5qFWrFt+tWsVzlSpxv78/u85y7hbg1sBAPq1Zk+9WrVKxiIhX\nULmcp0aNGvHLpk3Yhw+nXWgoV4WG8inwNbAY+AToEBhI9woVqP3II/y8fj01a9a0NrSIiJNoWswA\n2dnZzJw5k69nzuTooUP4+PiQmZtLQWAgP/74I/7+/uf+JSIiHkTlYpI9e/bQrl07kpOTte1YRLyO\npsVMEhsbS2BgIDt27LA6ioiI06lcTNSlSxd+/PFHq2OIiDidysVEKhcR8VYqFxOdLBcta4mIt1G5\nmKhu3bo4HA52795tdRQREadSuZjIZrNpakxEvJLKxWQqFxHxRioXk6lcRMQbqVxM1rBhQ3JyckhM\nTLQ6ioiI06hcTGaz2ejcubNGLyLiVVQuTqCpMRHxNioXJ9DIRUS8jcrFCeLj40lLSyMpKcnqKCIi\nTqFycQIfHx86deqk0YuIeA2Vi5No3UVEvInKxUlULiLiTXSzMCcpLCykYsWKJCQkEBMTY3UcERFT\naeTiJL6+vnTs2JEVK1ZYHUVExHQqFyfS1JiIeAuVixOpXETEW2jNxYkKCgqIiopi165dVKxY0eo4\nIiKm0cjFifz8/Gjfvj0//fST1VFEREylcnEyTY2JiDdQuTiZykVEvIHWXJwsPz+fqKgoEhMTiYyM\ntDqOiIgpNHJxMn9/f9q2bcvKlSutjiIiYhqViwU0NSYink7lYgGVi4h4Oq25WCAnJ4eKFSuyf/9+\nwsPDrY4jImI4jVwsEBQURKtWrfj555+tjiIiYgqVi0U0NSYinszP6gDeqn379jzwwANERESQcewY\noeHhNGjQgD59+hAQEGB1PBGRC6I1Fyf7+++/+WjsWD4ZN44a6el08/UlrLCQDF9ffrPb2WKzMfTu\nu7n7/vupWbOm1XFFRM6LysWJFi9ezK3XXccN+fnck5dHXCnnbAM+DAxkqp8fk2fMoFevXs6OKSJy\nwVQuTrJgwQLuHDyYudnZXFqG89cAfYODGf/FF/Tp08fseCIihlK5OMH27dvpeMklLMrKok05fu43\n4Gq7nZ/WrqVRo0ZmxRMRMZx2iznBO6+9xvDc3HIVC0Br4IHcXN557TUzYomImEYjF5Olp6dTKyaG\nzdnZVCvl+D7gXmAVEAAMBN4GfE8c3w/EBweTeOCAvnApIm5DIxeTTZs6lW4+PqUWC8AIoCKQDGwA\nfgQ+OOV4VaC7jw+fT5liblAREQOpXEy28ddfuSwz84zHtwCDKR61VAauPPHYqS7LzGTjr7+allFE\nxGgqF5MdO3qUiLMc7wlMA7KBJOBr4KoS50QAx44cMSegiIgJVC4mCwkLI+ssx58FNgPhQA2KF/Gv\nLXFOJhCi9RYRcSMqF5PVatCATWe4nIuD4pHLICALOAwcBUaWOG9zQAC1GjQwM6aIiKG0W8xke/fu\npUXDhuzNySGkxLFDFK+zHAPCTjw2D3gK2HTiz5lAraAg1m7bRq1atZySWUTkQmnkYrKaNWvSqUMH\nppVyrCJQBRgHFAJpwKdA81PO+QLo0L69ikVE3IrKxQlG/Pe/vGS3k1LicRswB1hAcdHUBwKBt04c\nPwi8FBLCA08+6bSsIiJGULk4weWXX86tDz5Ir5AQDpU41hb4CUileJrsC6DSiX/uZbdz0wMPcMUV\nVzg3sIjIBVK5OMmzL71Ez+HDudRuZzaQf4bzCigezbS327ninnt47uWXnRdSRMQgWtB3snnz5vHG\ns8+ye8cO7szLo01hIWFABvCLjw/jbDYaxMfzyLPP0q9fP6vjioicF5WLRTZu3MjH773Hn1u2kJ6e\nTlhYGDGxsXy1dCkHDx7Ez083CRUR96VycTGtWrVizJgxdO3a1eooIiLnTWsuLmbgwIHMmjXL6hgi\nIhdEIxcXs3PnTjp37kxSUhI+Pup+EXFPevdyMfXr1yc6OppVq1ZZHUVE5LypXFyQpsZExN1pWswF\nbd26lZ49e5KYmKipMRFxS3rnckGNGzcmLCyM3377zeooIiLnReXiogYMGKCpMRFxW5oWc1F//PEH\n/fr1Y9euXdhsNqvjiIiUi0YuLqpZs2b4+vqyfv16q6OIiJSbysVF2Ww2BgwYwOzZs62OIiJSbioX\nFzZw4EBmzpyJZi5FxN2oXFxYy5YtycvLY/PmzVZHEREpF5WLC9PUmIi4K5WLi9O39UXEHalcXFzb\ntm1JTU1l27ZtVkcRESkzlYuL8/Hx0dSYiLgdlYsb0Lf1RcTdqFzcQMeOHUlOTmbXrl1WRxERKROV\nixvw9fWlX79+mhoTEbehcnETWncREXeiC1e6ifz8fKpWrcrvv/9OrVq1rI4jInJWGrm4CX9/f/r0\n6cOcOXOsjiIick4qFzcycOBATY2JiFvQtJgbycvLIyYmhs2bN1O1alWr44iInJFGLm4kICCAa665\nhrlz51odRUTkrFQubkbXGhMRd6BpMTeTk5NDTEwMO3bsIDo62uo4IiKl0sjFzQQFBXHVVVcxb948\nq6OIiJyRysUNaWpMRFydpsXcUGZmJlWrVmX37t1ERUVZHUdE5DQaubihkJAQunfvzldffWV1FBGR\nUqlc3JQuwy8irkzTYm4qPT2dGjVqsHfvXiIiIqyOIyLyLxq5uKnw8HAuu+wyFixYYHUUEZHTqFzc\nmC7DLyKuStNibiw1NZXY2Fj27dtHWFiY1XFERP6hkYsbi4yMpH379ixevNjqKCIi/6JycXO6DL+I\nuCJNi7m5w4cPU7duXZKTk7Hb7VbHEREBNHJxexUrVqR169YsXbrU6igiIv9QuXgAXWtMRFyNpsU8\nwIEDB4iLi+PAgQMEBgZaHUdERCMXTxATE0OzZs349ttvrY4iIgKoXDyGpsZExJVoWsxDJCUl0axZ\nM5KTkwkICLA6joh4OY1cPES1atVo2LAh33//vdVRRERULp5El+EXEVehaTEPsmfPHtq0acP+/fvx\n8/OzOo6IeDGNXDxIbGwstWrVYsWKFVZHEREvp3LxMLoMv4i4Ak2LeZidO3fSuXNn9u3bh6+vr9Vx\nRMRLaeTiYerXr090dDSrVq2yOoqIeDGViwfSZfhFxGqaFvNACQkJ9OjRg8TERHx89PlBRJxP7zwe\nKC4ujrCwMH777Tero4iIl1K5eChda0xErKRpMQ/1xx9/0LdvX3bv3o3NZrM6joh4GY1cPFSzZs3w\n8/Nj/fr1VkcRES+kcvFQNptNU2MiYhmViwc7eSFLzXyKiLOpXDxYy5YtycvLY/PmzVZHEREvo3Lx\nYDabTZfhFxFLqFw8nL6tLyJWULl4uLZt25KWlsa2bdusjiIiXkTl4uF8fHzo37+/Ri8i4lQqFy+g\nLcki4mwqFy/QoUMHkpOT2bVrl9VRRMRLqFy8gK+vL/369dPUmIg4jcrFS2hqTEScSReu9BIFBQVU\nqVKF33//nVq1alkdR0Q8nEYuXsLPz49rr72WOXPmWB1FRLyAysWL6Nv6IuIsmhbzInl5ecTExLBp\n0yaqVatmdRwR8WAauXiRgIAAevfuzdy5c62OIiIeTuXiZQYMGKAtySJiOk2LeZmcnBxiYmLYvn07\nlStXtjqOiHgojVy8TFBQEFdddRXz5s2zOoqIeDCVixfSZfhFxGyaFvNCWVlZVKlShd27dxMVFWV1\nHBHxQBq5eCG73U737t2ZP3++1VFExEOpXLyUpsZExEyaFvNSGRkZVKtWjb1791KhQgWr44iIh9HI\nxUuFhYXRtWtXFi5caHUUEfFAGrl4sSlTpjB79mxtSy6ngwcPsmjRIg4dOkRhYSGRkZF07dqVhg0b\nWh1NxGWoXLxYWloaNWvWJCkpibCwMKvjuLxVq1bxweuvs+jrr+np60uN3Fx8i4o4HBjIAoeDps2a\nMXzkSPr06YOfn5/VcUUspXLxcldffTW33norgwcPtjqKyyosLOTBe+5h0fTpjMjK4laHg4tKnJML\nzAHeCQ0lKC6OuUuXEhkZaUFaEdegNRcvp8vwn53D4eC2wYNJmDaNDZmZPFRKsQAEAtcDPx8/zsUb\nN9K1TRvS09OdnFbEdWjk4uUOHz5M3bp1SU5Oxm63Wx3H5bzy/PN89dprLM/KIriMP+MA7gsM5O/2\n7VmwfLmZ8URclkYuXq5ixYq0adOGJUuWWB3F5WRnZ/PG6NF8Xo5iAbAB7+bm8scvv7Bhwwaz4om4\nNJWLaGrsDGbMmEEbm426JR4PBcJO+csPGFHiHD/grtxcxr31lvlBRVyQpsWElJQUGjVqRHJyMkFB\nQVbHcRltGzfmqYQErjnLOZlADPA10LHEsQNAXFAQiSkphIeHmxVTxCVp5CJUrlyZZs2a8e2331od\nxWU4HA427NhBt3OcNwuozOnFAsWlExsQwM6dOw3PJ+LqVC4C6FpjJeXm5gKcc63lU+CWsxyvYLNx\n7Ngxo2KJuA2ViwDQv39/FixYQF5entVRXEJAQACFDgcFZzknEVgB3HqWc7JAu/DEK6lcBIBq1arR\nsGFDlmvrLAA+Pj5Ui4wk4SznTAE6AbXOcDwX2J2bS/Xq1Q3PJ+LqVC7yD02N/dutd93F+MDAMx7/\njLOPWuYAFzdvrnIRr6TdYvKPPXv20Lp1a5KTk3VtLGDfvn00q1ePvbm5hJY4tgroAaQAIWf4+U5h\nYTz86af069fP1JwirkgjF/lHbGwssbGxrFixwuooLqF69eq0a9uWUaUc+wwYwJmLZS7wd3AwvXv3\nNi2fiCtTuci/DBw4UF+oPGHevHn8umULMytU4HVf338d+5DinWKlWQbcZbcze/FijQDFa2laTP7l\nzz//pGPHjiQlJeFb4g3VWxQUFPC///2PGTNmMHPmTCpXrsyVnTvTITmZJ/PyqHOGn0sFxgGj/f35\n6rvv6Ny5sxNTi7gWjVzkX+rVq0dMTAyrVq2yOoolDhw4wBVXXMGGDRv4/fffad26NTVr1uTn9euJ\nuPtu2oaG0is0lC+B1cBvwGJgaHAwdYKC2Ni7N6HR0Rw9etTafxERi2nkIqd54YUXOHz4MO+8847V\nUZxq5cqVDBkyhGHDhvHUU0+VOnLLzs7myy+/ZNbEiRw6eJCCggIiIyPpdu21DL3rLipVqsTq1avp\n168f69ato2rVqhb8m4hYT+Uip0lISKBHjx4kJibi4+P5g1uHw8Hbb7/Na6+9xqRJk7jqqqsu+Hc+\n99xzrFy5kqVLl3rFayhSkv5XL6eJi4sjPDycX3/91eoopsvIyOC6665j6tSprFmzxpBiAfjf//5H\nVlYWb7/9tiG/T8TdqFykVN5wGf4tW7bQunVroqKiWLlyJbGxsYb9bj8/Pz7//HNeeeUV3dNFvJLK\nRUp18tv6njprOn36dC677DKefPJJPvzwQ1NuNVC7dm3eeustbrjhBrKysgz//SKuTGsuUiqHw0GD\nBg344osvaNmypdVxDJOXl8cjjzzC119/zezZs2nevLmpz+dwOLjxxhuJjIxk7Nixpj6XiCvRyEVK\nZbPZPO5aY3///TddunRh7969/P7776YXCxS/jh988AGLFi1i4cKFpj+fiKtQucgZDRw4kJkzZ3rE\n1NiyZcto06YNffv2Ze7cuVSoUMFpz12hQgWmTJnCnXfeyYEDB5z2vCJWUrnIGV1yySUUFBSwadMm\nq6Oct6KiIl5++WVuvvlmpk2bxsiRIy3ZGtypUyeGDRvGbbfdRlFRkdOfX8TZVC5yRjabjQEDBrjt\n1FhqairXXnstixYt4rfffqNr166W5nn66adJS0vj/ffftzSHiDOoXOSs3HVL8vr162nVqhX16tXj\nhx9+oFq1alZHwt/fn6lTp/LCCy+49WhQpCxULnJWbdu25dixYyQknO2ejK5l4sSJ9OjRg1deeYW3\n3noLf39/qyP9o27duowZM4YbbriB7Oxsq+OImEZbkeWcHnzwQSpVqsSoUaXd2cR1ZGdn88ADD7Bq\n1Spmz55NXFyc1ZFK5XA4GDJkCDExMf9cvy09PZ0vv/yS3X/+SeaxY4RHRRHfpAn9+vUj8Cx3wxRx\nVSoXOacVK1YwYsQIl/6m+e7duxk4cCANGjTgk08+ITS05L0jXUtqairNmzfnySefZOOvvzJjxgy6\n+fjQPDOTECADWBkWxh/AHXfdxb0jRlCzZk2LU4uUncpFzqmwsJBq1aqxcuVK6tWrZ3Wc0yxcuJCh\nQ4cyatQo7r//fmw2m9WRyuTxRx/lwzfe4BFfX+4qLKRKKefsAMYFBDA1IICpc+bQvXt3Z8cUOS8q\nFymTe++9l9jYWEaOHGl1lH8UFhbyzDPP8OmnnzJjxgzat29vdaQy+2jcOF579FGWZGXRoAzn/wQM\ntNuZMncuPXr0MDueyAVTuUiZLFu2jCeffNJlrpR86NAhbrjhBoqKipg+fTrR0dFWRyqzVatWMaB7\nd1ZmZVG3HD/3EzAgNJTft2zRFJm4PO0WkzLp0qULf/31F4mJiVZHYc2aNbRs2ZI2bdrwzTffuFWx\nALz+7LM8e4ZiuQwIBsJO/HXqloROwE15eYx7913zQ4pcIJWLlImfnx/XXnutpV+odDgcjB07lj59\n+vD+++/z0ksvlXq3SFe2b98+fvzpJ248w3EbMJbiBf0MoOQG8Hvy8pg4fjy5ublmxhS5YCoXKTMr\nL2SZmZnJTTfdxMcff8zq1avp06ePJTku1CcffsgNDgdn28t2tnnqBkBzYM6cOcYGEzGYykXK7PLL\nLychIYGkpCSnPu/27dtp27YtAQEBrF69mrp1y7NS4Vq2/PYbnc4x6ngSqAR0BH4s5XjHjAy2btli\nQjoR46hcpMwCAgLo3bs3c+fOddpzzpo1i06dOvGf//yHiRMnEhwc7LTnNkP6sWOEneX4a8BfwH7g\nLqA3sLvEOeFA+uHD5gQUMYjKRcpl4MCBTrnWWH5+Po888giPPfYYX3/9NcOGDXOb76+cTUhoKJln\nOd4GCAH8gVuADsDiEuccB0KceMsAkfOhcpFy6d69Oxs2bCAlJcW050hOTqZbt25s27aNtWvXetSd\nMGvHxfHHBW5C+MNuJ7ZOHYMSiZhD5SLlEhQUxNVXX828efNM+f0//vgjrVq1onv37ixYsICLLrrI\nlOexyh333MOEgADySzl2DFgK5AAFwFSKv9ty5SnnJAPfORwMHjzY9KwiF0LlIuVmxmX4HQ4HY8aM\nYfDgwUyaNImnnnrKkpt6mS0+Pp4GcXGUVs35wFNANMUL+mOB+cCpF9z5xNeX6wYNIiIiwvywIhdA\n39CXcsvKyqJKlSrs3r2bqKioC/59x44d4/bbbycpKYmZM2d6/LfP582bx8gbb2R1VhblGZftBDra\n7Xy3Zg1NmzY1K56IITzvo6GYzm6306NHD+bPn3/Bv2vjxo20bt2aKlWqsGLFCo8vFoC+ffvSZ9gw\netntHC3jz+wGrrTbefHNN1Us4hZULnJejJgamzJlCt26dePpp59m7NixXnXfktfeeouOQ4dyqd3O\nl1DqGgwU7wwbD3QIDubR0aO58+67nRdS5AJoWkzOS0ZGBtWrVycxMZEK5dwWm5uby3/+8x+WLVvG\n7NmzvfqT+Ny5c3nnxRfZkZDA0Lw8Li4sJJTixf0f/fyY7utL544deeSZZ+jUqZPVcUXKTOUi56Ww\nsJAOHTpQeOwYySkpHMvKItjfnxoxMdxy333cetttpZbO3r17GThwIDVq1GDSpEmEh4dbkN71bNmy\nhYnjxvHXtm0cz8gg5eBBomrUYPLnn3vFVKF4HpWLlNvEjz/mxVGjCEtP58GcHK4AKgDZwFZgvN3O\nkqIiBl93HWPef5+wsOLvpC9dupRbb72Vxx57jIcfftgjvhRplqlTpzJ//ny+/PJLq6OInBeVi5SZ\nw+HgsREjWDRxIpOysmhL8VV8S3MA+F9gIOtq1mTxjz/y8ccf89FHHzFt2jS6dOnixNTuaePGjQwe\nPJiEhJLXRRZxDyoXKbMXnn6aeW+8wXdZWUSW4XwH8JSfHx8HBlKveXNmzZpFlSql3cxXSsrNzSUi\nIoK0tDSCgoKsjiNSbtotJmWyceNGPnj9dRaVsVigeFTzQkEBl+Xk0KltWxVLOQQGBlK3bl2NXMRt\nqVykTMa99Rb35uURU+LxPGAoEEvx1XpbAEtOOW4DXigsZNInn+gGV+XUtGlTNm/ebHUMkfOicpFz\nSk9P54sZMxhWWHjasQKgJrACSAdeBK4DTr0ZcgOgucNh6V0s3VHTpk3ZtGmT1TFEzovKRc5p0aJF\ndPLzo2opx+zAMxQXDEAvoDawrsR5dxw/zrSPPjIxpedp0qSJykXclspFzunAgQPULuOUVgqwA4gv\n8XhtICU52eBknk3TYuLOVC5yTvn5+QQUFZ37POBG4DaKp8JOFQDk5Z/pIidSmtjYWNLS0khNTbU6\niki5qVzknCpUqMDRgICznlME3AwEAe+XcvwoEKm7J5aLj48P8fHxGr2IW1K5yDl16NCBrx2OM15c\n0UHxjrFDwGygtPssfhUQQMcePcyK6LG07iLuSuUi5xQfH0/9Ro1KvcEVwL3ANuAroLTrGmcCn/v4\ncNd995kV0WNp3UXclcpFymT4yJG8FxpKycs5JFJ8Sfg/gBgg7MRf0085ZwrQqUMHXYDxPGg7srgr\nXf5FyiQvL48OLVrQa+dOni3Hwvxaim9ytWTFClq2bGleQA918OBBGjRoQGpqqi70KW5FIxcpk4CA\nABYuX8606GhG+fufNoIpzQrg6uBgPp46VcVynqKjowkMDGTfvn1WRxEpF5WLlFnlypX5ef16foiP\np3loKB8CGSXOcQDLgAEhIQwIDWXqV1/Rt29f54f1IFp3EXekcpFyqVSpEj+tW8db8+bxbc+e1AoK\nokdEBIPCwugRHExNPz8eio2l+5gx7N6/nyuuuMLqyG5P6y7ijvysDiDux2az0a1bN7p160ZSUhKb\nNm3i2LFj7Nu3j8Tx4/lj2zatDxioSZMm/PDDD1bHECkXLeiLYbKysoiKiiI9PR1/f3+r43iMX3/9\nlbvuuosNGzZYHUWkzDQtJoax2+1Ur16dnTt3Wh3Fo8THx7N9+3YKCgqsjiJSZioXMZTWB4wXEhJC\ntWrVVNriVlQuYijtbDKHLgMj7kblIobSyMUcKm1xNyoXMZTKxRx6XcXdqFzEUHXr1iU5OZnjx49b\nHcWjqFzE3ahcxFB+fn40atSILVu2WB3Fo9SrV4+kpCQyMzOtjiJSJioXMZw+ZRvP39+fhg0bsnXr\nVqujiJSJykUMp8Vnc6i0xZ2oXMRwehM0h7YjiztRuYjhVC7m0Osq7kTlIoarUqUKBQUFpKSkWB3F\no2i6UdyJykUMZ7PZ9CnbBNWrVycnJ4dDhw5ZHUXknFQuYgqVi/FsNpvWXcRtqFzEFCoXc+h1FXeh\nchFTaH3AHHpdxV2oXMQUTZo0YevWrRQVFVkdxaNo5CLuQuUipggPD6dixYrs3r3b6igepUmTJmzZ\nskWlLS5P5SKm0eKz8SIjI4mIiCAxMdHqKCJnpXIR02gKxxx6XcUdqFzENHoTNIdGhOIOVC5iGu1s\nModKW9yec+66AAAgAElEQVSBykVM07BhQ/bs2UNOTo7VUTyKSlvcgcpFTBMQEEC9evVISEiwOopH\niYuLY9euXeTl5VkdReSMVC5iKq0PGC8oKIhatWqxbds2q6OInJHKRUyl9QFz6HUVV6dyEVPpTdAc\nWncRV6dyEVOpXMyh11VcncpFTFWrVi0yMjI4evSo1VE8itayxNWpXMRUJ+9BoikcY9WpU4fDhw+T\nnp5udRSRUqlcxHT6lG08X19fGjdurNIWl6VyEdNpfcAcel3FlalcxHR6EzSHRoTiylQuYrqT22Yd\nDofVUTyKtiOLK1O5iOmioqIICQnh77//tjqKRzk5IlRpiytSuYhTaGrMeJUrV8Zms5GcnGx1FJHT\nqFzEKVQuxrPZbHpdxWWpXMQptPhsDq27iKtSuYhT6BO2OfS6iqtSuYhTNG7cmJ07d5Kfn291FI+i\nEaG4KpWLOEVwcDA1a9Zkx44dVkfxKE2aNCEhIYHCwkKro4j8i8pFnEZTOMYLCwujcuXK7Nq1y+oo\nIv+ichGnUbmYQ6+ruCKViziN1gfModdVXJHKRZxGn7DNoe3I4opULuI0devW5eDBg2RkZFgdxaOo\ntMUVqVzEaXx9fYmLi9OnbIM1aNCAvXv3kp2dbXUUkX+oXMSpNIVjvICAAOrVq8fWrVutjiLyD5WL\nOJWmcMyh0hZXo3IRp1K5mEOvq7galYs41clts7oHibG0HVlcjcpFnComJgaAAwcOWJzEs2jkIq5G\n5SJOpXuQmKNWrVocP36cI0eOWB1FBFC5iAW0+Gw8m81GkyZN9LqKy1C5iNNp5GIOrbuIK1G5iNOp\nXMyh11VcicpFnC4+Pl73IDGBphvFlahcxOnCwsKIjo7WPUgMdnLNRdu8xRWoXMQSmsIxXsWKFbHb\n7ezdu9fqKCIqF7GGysUcmhoTV6FyEUvoTdAcKm1xFSoXsYTeBM2h7cjiKlQuYgndg8QcKm1xFSoX\nsYS/vz/169fXPUgM1rhxY3bu3El+fr7VUcTLqVzEMvqUbTy73U6NGjXYsWOH1VHEy6lcxDIqF3No\n3UVcgcpFLKMdY+ZQaYsrULmIZfQmaA6VtrgClYtYpkaNGmRmZuoeJAZTaYsrULmIZU7eg0RvhMaq\nW7cuBw4cICMjw+oo4sVULmIpfco2np+fH40aNdI2b7GUykUspXIxh15XsZrKRSylxWdzqFzEaioX\nsdTJctE9SIyltSyxmspFLBUZGUl4eDiJiYlWR/EoGhGK1VQuYjl9yjZe1apVKSgoICUlxeoo4qVU\nLmI5rQ8YT9u8xWoqF7GcysUcel3FSioXsZzeBM2hdRexkspFLBcXF8euXbvIy8uzOopHUWmLlVQu\nYrmgoCBiY2PZvn271VE8Snx8PFu3bqWoqMjqKOKFVC7iEvQp23gVKlTgoosuYvfu3VZHES+kchGX\noJ1N5tC6i1hF5SIuQSMXc+h1FauoXMQl6E3QHBoRilVULuIS6tSpw5EjR0hPT7c6ikfRtJhYReUi\nLsHHx4fGjRvrjdBgjRo14q+//iI3N9fqKOJlVC7iMjQ1ZrzAwEDq1KlDQkKC1VHEy6hcxGVofcAc\nel3FCioXcRkauZhD6y5iBZWLuIyT5aIbhxlLpS1WULmIy6hcuTJ+fn7s37/f6igeRdNiYgWVi7gU\nTeEYr3bt2qSmppKammp1FPEiKhdxKZrCMZ6Pjw/x8fFs2bLF6ijiRVQu4lJULubQ6yrOpnIRl6L1\nAXPodRVnU7mIS4mPj2fbtm0UFBRYHcWjaOQizqZyEZcSGhpKlSpV+PPPP62O4lFObpTQNm9xFpWL\nuBztGDNedHQ0AQEBJCUlWR1FvITKRVyOpnDMoXUXcSaVi7gclYs5NCIUZ1K5iMvRJ2xzqLTFmVQu\n4nLq169PUlISmZmZVkfxKCptcSaVi7gcf39/GjRowNatW62O4lG0zVucSeUiLknrA8YLDQ2latWq\n2uYtTqFyEZek9QFz6HUVZ1G5iEvSm6A5tO4izqJyEZekN0FzqLTFWVQu4pKqV69OTk4Ohw4dsjqK\nR9FaljiLykVcks1m06dsE9SvX599+/Zpm7eYTuUiLkvlYjxt8xZnUbmIy9IUjjlU2uIMKhdxWXoT\nNIdKW5xB5SIuq0mTJmzZsoWioiKro3gUlbY4g8pFXFaFChWoUKECe/bssTqKR9E2b3EGlYu4NH3K\nNl6NGjW0zVtMp3IRl6ZyMZ7NZqNJkyZadxFTqVzEpWnx2RyaGhOzqVzEpWnkYg69rmI2lYu4tEaN\nGrF7925yc3OtjuJRNCIUs6lcxKUFBgZSu3Zttm3bZnUUj3JyzUXbvMUsKhdxeZrCMd5FF11EeHg4\niYmJVkcRD6VyEZencjGHXlcxk8pFXJ7eBM2hdRcxk8pFXJ7eBM2h0hYzqVzE5cXGxpKamkpaWprV\nUTyKvusiZlK5iMvz8fEhPj5eoxeDxcXFsWvXLvLy8qyOIh5I5SJuQZ+yjRccHEytWrXYunUrmZmZ\nOBwOqyOJB1G5iFvQ+oCx9u7dy/9GjuTA7t20uuQSKlaoQICfH60aNuSTTz4hKyvL6oji5lQu4hZU\nLsY4dOgQA668kosbNOD422+zKj+fXIeD7IICMouKeHHHDr566CFqVKrE86NG6UuWct5sDo2FxQ0c\nOnSIBg0acPToUWw2m9Vx3FJiYiJXtG9P/0OHeDo/n5CznLsHuNFup2b37kyZNQs/Pz8npRRPoZGL\nuIVKlSoRGBhIUlKS1VHcUmpqKld16cLwAwd47RzFAhALLMvK4si33zLirru0HiPlpnIRt6GpsfP3\n0jPP0CE5mYfKMc0VBMzKymLpl1/y888/mxdOPJLKRdyGdoydn+zsbD6dNIknS9ly/D7QiuIiub2U\nnw0HRmRl8cHrr5sbUjyOykXchkYu52fGjBm0AeqUcqwa8BRwx1l+/laHg6+XLiUlJcWUfOKZVC7i\nNlQu52fO5Mnccvx4qcf6AdcCUWf5+QrAVb6+LFy40IR04qlULuI24uPj2b59O/n5+VZHcSuHUlKo\ncY5zzrVcXz03l8OHDxsVSbyAykXcht1up3r16vz5559WR3ErhUVF+J7jnHNt7vYtKqKgoMCoSOIF\nVC7iVjQ1Vn6RkZEcOsc55xq5HA4MJDIy0qhI4gVULuJWtGOs/Lr17cus4OCznnO2kUsu8JXDweWX\nX25oLvFsKhdxKxq5lN8dw4Yxr6iII6UcKwRygIIT/5x74u+nmg00bdaMRo0amRtUPIrKRdyKyqX8\nKlasSJ9rrmGcz+n/d38BsAOvAZ8DwcBLpxwvBN4NDeW+J55wRlTxILq2mLiVgoICwsPDOXjwIKGh\noVbHcRurV6+mR8eOzC4qokcZf8YB/CcggI0XX8y3P/+s64tJuWjkIm7Fz8+PRo0asXXrVqujuI1N\nmzZx/fXXM/CWW7gpJISZZfiZfIq/WPlDjRrMXbpUxSLlpnIRt6OpsbJbtGgR3bp14+WXX2bSpEl8\ns3IlIytX5tKwMKZQvN5yqmTgeV9fatvtrK5enWoNGhAeHm5BcnF3KhdxOyqXc3M4HLz11lvceeed\nzJ8/nxtuuAGAiy++mJ1JSTz5+edM7dCBakFBtIqIoHNEBM3Cw4kPDib55ptZvHo1G3ftIjU1ldGj\nR1v8byPuSGsu4naWLFnCmDFjWLZsmdVRXFJ+fj733Xcfa9asYcGCBdSqVeuM5+7fv5+kpCSysrKI\niIigbt26hIWF/XN83759tG7dmqlTp2orspSLykXcTlJSEi1atODgwYNWR3E5R48eZeDAgdjtdqZP\nn/6vojhfy5Yt46abbuL333+nWrVqBqQUb6BpMXE7VatWpaCgQFfpLWHHjh20a9eOFi1aMH/+fEOK\nBaBbt2488MADDBo0iLxSLtsvUhqVi7gdm81G06ZN2bx5s9VRXMby5cvp1KkTjz/+OG+88Qa+vue6\nmlj5PPHEE0RFRfH4448b+nvFc6lcxC1pUf//jR8/nuuvv54vvviCYcOGmfIcPj4+fPbZZyxYsIAZ\nM2aY8hziWbR5XdxS06ZN+fXXX62OYanCwkIeffRRFi9ezMqVK6lfv76pzxcZGcmsWbPo0aMHzZo1\nIy4uztTnE/emkYu4JW+/gGV6ejp9+vRh48aNrFmzxvRiOalFixa8+uqr9O/fn4yMDKc8p7gn7RYT\nt3Ts2DGqVatGeno6PqVcM8uT7dmzh969e9OhQwfee+89/P39nZ5h6NChZGZmMn36dGy2c90NRryR\nd/2/UjxGREQEUVFR7N692+ooTrVq1Srat2/PsGHDGDdunCXFAvD++++zfft23nvvPUueX1yfykXc\nlrct6k+dOpW+ffsyYcIEHnzwQUtHDMHBwcyePZsXX3yR1atXW5ZDXJfKRdyWt2xHLioqYtSoUYwa\nNYrly5dz1VVXWR0JgDp16jBhwgSuu+46faFVTqNyEbflDSOXrKwsBg8ezPfff88vv/xCkyZNrI70\nL7179+aWW27h+uuvp7Cw5G3GxJupXMRteXq57N+/n86dOxMUFMSyZcuIjo62OlKpnn/+eWw2G08/\n/bTVUcSFqFzEbTVs2JA9e/aQk1PywvHub926dbRt25b+/fvz2WefERQUZHWkM/L19WXatGlMmTKF\nBQsWWB1HXITKRdxWQEAAdevWJSEhweoohpozZw49e/bk7bff5r///a9bbPWNjo5mxowZDB061Ot2\n8EnpVC7i1jxpaszhcPDKK6/w4IMPsmTJEgYMGGB1pHK59NJLeeqppxgwYADZ2dlWxxGLqVzErXnK\njrHc3Fxuu+02Zs2axZo1a2jZsqXVkc7L/fffT6NGjbj//vutjiIWU7mIW/OEkcuhQ4fo1q0bx48f\nZ8WKFW59zxSbzcbHH3/M6tWrmTBhgtVxxEIqF3Fr7l4uW7ZsoW3btnTp0oWZM2cSEhJidaQLFhoa\nypw5c3jiiSdYt26d1XHEIioXcWs1a9YkPT2d1NRUq6OU25IlS+jatSvPPfccL730kkddI61Ro0aM\nHTuWgQMHcvToUavjiAV04Upxa9nZ2TRv3pwePXoQExNDREQELVq0oEOHDi67y8rhcPDee+/xyiuv\nMGvWLDp06GB1JNM89NBD7Ny5k6+++sqjylPOTeUibmnXrl18+O67TJ44kca5uTQvKiLU4eBYQADL\n/fzwi4pi+OOPc9PNNxt2u18j5OfnM2LECH766ScWLlxIbGys1ZFMlZ+fT9euXbnqqqv43//+Z3Uc\ncSKVi7idyRMn8tj993NHQQF35+dTp8RxB/AD8H5ICOtCQlj0/fc0btzY+UFLSE1NZdCgQQQEBPDF\nF18QHh5udSSnSEpKonXr1nz66ad0797d6jjiJBqnilv54L33eP6BB1iZnc1rpRQLgA3oCszOzOS5\nQ4e4/NJLLf+i5c6dO7n00ktp2rQpX331ldcUC0C1atWYOnUqN998M3///bfVccRJNHIRt7F8+XJu\nuuYaVmZnl1oqZ/KZzcYzlSqx8c8/LZki++GHHxg8eDDPP/88d999t9Of31W8+uqrzJ8/nx9//JGA\ngACr44jJNHIRt/Hyk08yppzFAnCLw8ElmZl8PmWKKbnOZsKECQwePJhp06Z5dbEAPP7440RHR/Po\no49aHUWcQCMXcQvbtm2jS4sW7M3JIbCU418AzwF/AzHAZKDjKceXAw/WqsXGv/5yyi6ywsJCRo4c\nyfz581m4cCENGzY0/TndQVpaGq1ateKFF17g+uuvtzqOmEgjF3ELH737LkMLCkotlm+BJ4BPgePA\nT3Da6KYrkH/kCD///LO5QYGMjAz69u3L2rVrWbNmjYrlFBUqVGDWrFmMGDGCLVu2WB1HTKRyEbew\nYc0auhUUlHrsmRN/tTnx5ypA1RLn2IBuBQWsX7/etIwAiYmJdOzYkZiYGJYuXUpUVJSpz+eOLr74\nYsaMGcOAAQPIyMiwOo6YROUibuHYsWNElPJ4IbAWOAjUB2oADwCl3eGlQm4ux44dMy3jmjVruPTS\nS7ntttsYP368Fq3P4rbbbqNz584MHToUzcx7JpWLuAV7cDClXcQ9BcgHZgMrgQ3AeuDFUs7N8vfH\nbrebkm/69On06dOH8ePH89BDD7ns1QFcybvvvsuuXbt45513rI4iJlC5iFuoHhvL1lIeDz7x9weA\nykAU8DCwuJRzEwIDqV69uqG5HA4HzzzzDE8++STfffcd11xzjaG/35MFBQUxa9YsXnnlFVauXGl1\nHDGYykXcwq333cdHoaGUnECJBMpSF7uBtQ4HvXv3NixTdnY2Q4YM4ZtvvuGXX36hWbNmhv1ub1G7\ndm0mTZrEkCFDSElJsTqOGEjlIm6hZ8+eHAsJ4ddSjt0OvAccAlKBt4CSFfKRvz+33XEHwcHBJX/8\nvCQnJ9OlSxd8fX35/vvvqVy5siG/1xtdffXV3HHHHQwZMoSCM2zaEPejchG34OPjw/2PPcZIu528\nEseeAloDDYDGQEvg1Eskbgcm+vlxz4gRhmTZsGEDbdu2pXfv3kydOpWgoCBDfq83e+aZZ/D392fU\nqFFWRxGD6EuU4jYKCwvpf+WV2Feu5NOcHMqyF+sv4Aq7nVHvvsvtQ4decIb58+czbNgwxo4dy3XX\nXXfBv0/+3+HDh2nZsiXvvPMOffv2tTqOXCCNXMRt+Pr6Mn3+fLI6dqSn3c4vcNoazEm5wDSgpc3G\nvU8/fcHF4nA4GD16NPfddx+LFy9WsZigYsWKfPnll9x11138+eefVseRC6SRi7idwsJC3n37bd4f\nM4YKWVncnZFBHGAH0oDlfn5M8PenSdOmBEdHExYWxtSpU//ZHlxUVMS+fftITU3F19eXSpUqnXXN\nJC8vj3vuuYf169ezYMECw3ecyb998MEHfPTRR6xevdq0reNiPpWLuK2ioiK++eYbJr//Pvv27CEr\nO5uIiAgubteOu0eMoFGjRmRnZ9O6dWsee+wx+vTpw+SJExn3xhtkpKVRyd+fQiA5N5f4Ro0YPnIk\n/fv3JzDw/y8yc/jwYfr3709UVBRTpkwhNDTUun9hL+FwOLj55pvx8/Nj0qRJ+s6Qm1K5iMfbtGkT\nHdu1w1ZYyDW+vgzPyuJSii8JA1AALAA+CAtjk83Gx59/Tu/evUlISOCaa65h0KBBvPzyy7pNrxNl\nZmbStm1bHnzwQe68806r48h5ULmIR3M4HNw/bBg/ff45i/PyzvmdmNXAgOBgbrjnHj77/HNGjx7N\nbbfd5oSkUtL27dvp1KkTixcvplWrVlbHkXJSuYhHe+nZZ5kzZgw/ZGVR1tuE7aD4IpijxozRvUcs\nNmvWLB577DHWrl3LRRddZHUcKQeVi3islJQUGsXGsiUn57SrJJ/LN8CD1auzde9ezflb7JFHHiEh\nIYGFCxdqatKN6L+UeKwJ48cz0GYrtVhuovjS/OEU3/vlpRLHuwN+aWl8//33JqeUc3n11VfJyMjg\npZdK/lcSV6aRi3ikwsJCaleuzLwjR7iklONbgLpAEMXf4O9C8d0rrzzlnHHAsiuvZNbXX5sdV85h\n//79tGrVismTJ9OjRw+r40gZaOQiHikhIYHA3NxSiwUgnuJiOckPiC5xzmBgyfLlZsSTcqpatSrT\np0/nlltuYe/evVbHkTJQuYhHOnr0KJV9fc96znAghOKiGQWnFVEkkFtQQF5eyauZiRW6dOnCI488\nwqBBg8jNzbU6jpyDykU8ko+PzxkvDXPSB8Bx4DuKy6W0Ky47HA4t6LuQRx99lKpVq/Lwww9bHUXO\nQeUiHikqKor9BQXnLBgbcBkwCJhe4thBICQwEH9/fxMSyvmw2WxMnjyZb775hs8//9zqOHIWKhfx\nSA0bNiQgIoKfy3h+PsVTZKea4uPDtb16GZxMLlRERASzZ8/moYceYtOmTVbHkTNQuYhH8vHx4d5H\nHmFcKRc+PAR8AWQChcBSYCZw7SnnFAEfBgczXF+idEnNmjXjzTffZMCAAaSnp1sdR0qhrcjisVJT\nU6lTrRq/ZWdT75THDwMDgT8ovmR/A4rXXPqccs5M4JX69Vm7fbvWXFzYvffey8GDB5k1a5b+O7kY\njVzEY0VGRvLqG29wtd3OqXdnrwj8QPEtkdMoXsg/tVh+A4bb7Xzw2Wd6w3Jxb7/9Nnv37uXNN9+0\nOoqUoHIRj3b3vfdy8yOP0N5u549znOsA5gJXBwcz8YsvaNeunRMSyoUIDAxk1qxZjBkzhhUrVlgd\nR06haTHxCpMnTuSJhx6icVERw48f51rg5B6wdGCKzcYHISEURUYy4YsvaN++vYVppbyWLFnC0KFD\nWbt2LTExMVbHEVQu4kXy8vKYM2cOY199lTWbNhEREEChw0FWQQF9evRg+GOPcdlll2kqzE0999xz\nLF++nGXLluHn52d1HK+nchGvlJeXR1paGn5+fkREROB7jm/zi+srKiqiV69eNG3alNGjR1sdx+up\nXETEYxw5coSWLVvy5ptv0r9/f6vjeDWVi4h4lN9++41evXqxcuVKGjRoYHUcr6XdYiLiUVq3bs3z\nzz/PgAEDyMzMtDqO19LIRUQ8jsPh4NZbb8XhcPCZvq9kCY1cRMTj2Gw2PvzwQ/744w8++ugjq+N4\nJY1cRMRj7dy5kw4dOrBw4ULatGljdRyvopGLiHis+vXr89FHHzFo0CAOHz78r2MOh4Njx46xd+9e\nDh48SH5+vkUpPZPKRUQ8Wr9+/Rg8eDA33ngjhYWFHDlyhNdHj6ZelSpUr1SJDnFxNK5Vi4vCwhh2\n002sW7fO6sgeQdNiIuLxCgoKuPzyyyEnh42bNnGtjw/Ds7JoQ/EN4wBSgAm+vnwYGEiNevX4fN48\nateubWFq96aRi4h4vKKiIkJsNop+/52dOTl8mpVFW/6/WAAqA/8tLOSvrCyu27yZDpdcwsaNGy1K\n7P40chERj+ZwOLj9+us5+tVXzM7Opqw3rf4CeCwqitUbNlC9enUzI3okjVxExKOtXLmSlQsXMr0c\nxQIwBLg5LY1nR440K5pHU7mIiEf7YMwYHsjKIuQs5+wEgoCbSzz+n8JCZs+ZQ2pqqnkBPZTKRUQ8\n1oEDB/j6m2+49Ryz//fBvxb3T4oGrvbx4dNJk0xK6LlULiLisb755ht6+PtT4SznfAFEAt0ovhtp\nSTdmZfHV1Kmm5PNkKhcR8VhHjhyhal7eGY+nA88Ab1F6sQBUAY4cPWp8OA+nchERj3WuzbBPAcOA\nqpw+JXaSrQy/R06ne4GKiMeKioriV39/KGX0sgFYBqw/8ecz1UcKcFFkpDkBPZjKRUQ81hVXXMF/\nCgpIB8JLHPsR2APUPPHn40AhkAD8fsp504OD6TVkiNlRPY6+RCkiHm3gVVdx+ZIlDC/xeDaQceKf\nHcDrFJfNh0DUicePAPWCgtj5999UrFjRGXE9htZcRMSjDX/sMd4LCSGnxOPBFG81jqb40i+hJx6L\nOuWcsT4+9LnmGhXLedDIRUQ8msPhYHCfPtiWLWNadja+Zfy5+cC9FSrw87p1uoDledDIRUQ8ms1m\n47OZMzl68cUMsNtJP8f5DmA8MDQ4mPnffKNiOU8qFxHxeEFBQSz64QeqDBpEbFAQ9wUGsrnEOWnA\nuzYbjUNDeb1qVWwhIVStWtWKuB5B5SIiXiEgIIBxkyez6c8/qfToo/SsUIGqdjtxYWHEhoRQMzCQ\n1ddcw/jFi9m+bx8PPvggQ4YM0R0qz5PWXETEKxUUFJCSkkJaWhrBwcFER0cTGhr6z/GioiJ69epF\n06ZNGT16tIVJ3ZPKRUTkDA4fPkzLli1577336NOnj9Vx3IrKRUTkLNasWcO1117LmjVrtLhfDlpz\nERE5i3bt2vHf//6XgQMHkpNT8tsyciYauYiInIPD4WDQoEFUqlSJcePGWR3HLWjkIiJyDjabjQkT\nJvDdd98xbdo0q+O4BY1cRETK6I8//uCKK65gxYoVxMXFWR3HpWnkIiJSRs2bN+e1115j4MCBZGZm\nWh3HpWnkIiJSTrfffjsFBQV89tln2GzFtxnLy8tjzZo1HDp0iMLCQiIjI2nbti3h4SUv9u8dVC4i\nIuWUlZVFu3btuP/++7n66qsZ/8EHfDJuHNWKiqgB+DocHPLxYWNeHkOGDOHehx6iadOmVsd2KpWL\niMh52LZtG21btsSnoICbbDbuzc2lcYlz9gOf+Pnxkb8/V/brx4eTJ+Pv729FXKfTmouIyHn47JNP\nqFpYyMa8PN4rpVgAqgJPFxSwIzubg3Pn0q9nT6+5VpnKRUSknD7+8ENmjxvHytxcapTh/BBgbnY2\njl9+YcSdd5odzyVoWkxEpBzy8/OpFR3N12lpNC/nz2YAtYOC+D0hgdjYWBPSuQ6NXEREymHevHk0\nKCykOfA+0AoIAm4vcd4yoBHFo5bLgb1AGHBLUREfvf++ExNbQyMXEZFy6NqqFfeuXct1wFyKP6Ev\nBbKBSSfOOQzUAyYAvYFRwE/AamAH0DEsjL8PHSIwMNDZ8Z1GIxcRkTJyOBz8vGEDvU/8uR9wLRBV\n4rw5QBNgABAAPAv8QXGxNDhx/o4dO5wR2TIqFxGRMsrKysLXZiO4xOMlp3+2wL/WY+wUj2RO3lr5\nIh8f0tLSTErpGlQuIiJlFBAQQH5h4WllYivx50yg5Pfyw4HjJ/4578Tv8mQqFxGRMvL39yfCbufv\nEo+XLJtQIL3EY8coXtAvAPbm5VG5cmWTUroGlYuISDkMGTKECb6+/3qs5MglnuI1lpMygV0nHl8E\n1KlTR1uRRUTk/w1/6CE+DgggHygEcigejRQCuSf+3o/i9ZU5J44/B1xM8WL+B6GhDB850oroTqWt\nyCIi5dS1dWv6r1vHkaIini9x7FngaYq/53I/kAi0AyZT/F2XAeHhJKakEBQU5MTEzqdyEREpp507\nd9KpVSs+TE+nbxl/ZhPQPTiYT+fOpWfPnmbGcwmaFhMRKaf69euzcNky7o2IYLSPD2e7bVgBMAPo\nFl8HqlIAAAIWSURBVBzMOxMnekWxgMpFROS8tGrVipVr1/LzZZdRMyiIBwMCWAccAo4C24EXfX2p\nbbfzTpMmzF66lMFDhlgb2ok0LSYicoESExMZP3Ysc6ZN41BaGoVFRUSGhnJFz57c+/DDtGjRwuqI\nTqdyERERw2laTEREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGcykVE\nRAynchEREcOpXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxE\nRMRwKhcRETGcykVERAynchEREcOpXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXERExnMpF\nREQMp3IRERHDqVxERMRwKhcRETGcykVERAynchEREcOpXERExHAqFxERMZzKRUREDKdyERERw6lc\nRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGcykVERAynchEREcOpXERExHAqFxERMZzK\nRUREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGcykVERAynchEREcOp\nXERExHAqFxERMZzKRUREDKdyERERw6lcRETEcCoXERExnMpFREQMp3IRERHDqVxERMRwKhcRETGc\nykVERAynchEREcOpXERExHAqFxERMZzKRUREDPd/qE1KA4BPGYoAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 73
    }
   ],
   "metadata": {}
  }
 ]
}

SkeletonToGraph.py

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

import Image, ImageDraw, ImageFont

import mahotas as mh
from mahotas import polygon
import pymorph as pm

import networkx as nx

from scipy import ndimage as nd
import skimage.transform as transform
import skimage.io as sio
import scipy.misc as sm

import os
import math

# <codecell>

def branchedPoints(skel):
    branch1=np.array([[2, 1, 2], [1, 1, 1], [2, 2, 2]])
    branch2=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 1]])
    branch3=np.array([[1, 2, 1], [2, 1, 2], [1, 2, 2]])
    branch4=np.array([[2, 1, 2], [1, 1, 2], [2, 1, 2]])
    branch5=np.array([[1, 2, 2], [2, 1, 2], [1, 2, 1]])
    branch6=np.array([[2, 2, 2], [1, 1, 1], [2, 1, 2]])
    branch7=np.array([[2, 2, 1], [2, 1, 2], [1, 2, 1]])
    branch8=np.array([[2, 1, 2], [2, 1, 1], [2, 1, 2]])
    branch9=np.array([[1, 2, 1], [2, 1, 2], [2, 2, 1]])
    br1=mh.morph.hitmiss(skel,branch1)
    br2=mh.morph.hitmiss(skel,branch2)
    br3=mh.morph.hitmiss(skel,branch3)
    br4=mh.morph.hitmiss(skel,branch4)
    br5=mh.morph.hitmiss(skel,branch5)
    br6=mh.morph.hitmiss(skel,branch6)
    br7=mh.morph.hitmiss(skel,branch7)
    br8=mh.morph.hitmiss(skel,branch8)
    br9=mh.morph.hitmiss(skel,branch9)
    return br1+br2+br3+br4+br5+br6+br7+br8+br9

def endPoints(skel):
    endpoint1=np.array([[0, 0, 0],
                        [0, 1, 0],
                        [2, 1, 2]])
    
    endpoint2=np.array([[0, 0, 0],
                        [0, 1, 2],
                        [0, 2, 1]])
    
    endpoint3=np.array([[0, 0, 2],
                        [0, 1, 1],
                        [0, 0, 2]])
    
    endpoint4=np.array([[0, 2, 1],
                        [0, 1, 2],
                        [0, 0, 0]])
    
    endpoint5=np.array([[2, 1, 2],
                        [0, 1, 0],
                        [0, 0, 0]])
    
    endpoint6=np.array([[1, 2, 0],
                        [2, 1, 0],
                        [0, 0, 0]])
    
    endpoint7=np.array([[2, 0, 0],
                        [1, 1, 0],
                        [2, 0, 0]])
    
    endpoint8=np.array([[0, 0, 0],
                        [2, 1, 0],
                        [1, 2, 0]])
    
    ep1=mh.morph.hitmiss(skel,endpoint1)
    ep2=mh.morph.hitmiss(skel,endpoint2)
    ep3=mh.morph.hitmiss(skel,endpoint3)
    ep4=mh.morph.hitmiss(skel,endpoint4)
    ep5=mh.morph.hitmiss(skel,endpoint5)
    ep6=mh.morph.hitmiss(skel,endpoint6)
    ep7=mh.morph.hitmiss(skel,endpoint7)
    ep8=mh.morph.hitmiss(skel,endpoint8)
    ep = ep1+ep2+ep3+ep4+ep5+ep6+ep7+ep8
    return ep

def pruning(skeleton, size):
    '''remove iteratively end points "size" 
       times from the skeleton
    '''
    for i in range(0, size):
        endpoints = endPoints(skeleton)
        endpoints = np.logical_not(endpoints)
        skeleton = np.logical_and(skeleton,endpoints)
    return skeleton

# <codecell>

image = Image.new("RGBA", (600,150), (255,255,255))
draw = ImageDraw.Draw(image)
fontsize = 150
font = ImageFont.truetype("/usr/share/fonts/truetype/liberation/LiberationMono-Regular.ttf", fontsize)
txt = 'A'
draw.text((30, 5), txt, (0,0,0), font=font)
img = image.resize((188,45), Image.ANTIALIAS)
plt.imshow(img)

# <codecell>

img = np.array(img)
im = img[:,0:50,0]
im = im < 128
skel = mh.thin(im)
# To fix the algorithm, prune the skeleton 1,2,3...
skel = pruning(skel,3)
bp = branchedPoints(skel)
ep = endPoints(skel)
edge = skel-bp

edge_lab,n = mh.label(edge)
bp_lab,_ = mh.label(bp)
ep_lab,_ = mh.label(ep)

# <codecell>

figsize(10,20)
subplot(141,frameon = False, xticks = [], yticks = [])
title('pruning 0')
imshow(mh.thin(im),interpolation = 'nearest')
subplot(142,frameon = False, xticks = [], yticks = [])
title('pruning 1')
imshow(pruning(mh.thin(im),1),interpolation = 'nearest')
subplot(143,frameon = False, xticks = [], yticks = [])
title('pruning 2')
imshow(pruning(mh.thin(im),2),interpolation = 'nearest')
subplot(144,frameon = False, xticks = [], yticks = [])
title('pruning 3')
imshow(pruning(mh.thin(im),3),interpolation = 'nearest')

# <codecell>

subplot(221,frameon = False, xticks = [], yticks = [])
title('skeleton')
imshow(skel,interpolation = 'nearest')
subplot(223,frameon = False, xticks = [], yticks = [])
title('skel-bp -> label')
imshow(edge_lab,interpolation = 'nearest')
subplot(224,frameon = False, xticks = [], yticks = [])
title('skel -> end points')
imshow(ep_lab,interpolation = 'nearest')
subplot(222,frameon = False, xticks = [], yticks = [])
title('branched points (bp)')
imshow(bp_lab,interpolation = 'nearest')

# <codecell>

endpoints = np.zeros(skel.shape)
edges = np.zeros(skel.shape)

for l in range(1,n+1):
    cur_edge = edge_lab == l
    cur_end_points = endPoints(cur_edge)
    pruned_edge = pruning(cur_edge,1)
    edges = np.logical_or(pruned_edge, edges)
    endpoints = np.logical_or(endpoints,cur_end_points)
    
lab_bp, nbp = mh.label(bp)
lab_ep, nep = mh.label(endpoints+ep)
pruned_edges, ned = mh.label(edges)

# <codecell>

plt.figsize(13,8)
subplot(321,frameon = False, xticks = [], yticks = [])
title(str(np.unique(skel)))
imshow(skel, interpolation='nearest')

subplot(322,frameon = False, xticks = [], yticks = [])
title(str(np.unique(lab_bp)))
imshow(lab_bp, interpolation='nearest')

subplot(323,frameon = False, xticks = [], yticks = [])
title(str(np.unique(edge_lab)))
imshow(edge_lab, interpolation='nearest')

subplot(326,frameon = False, xticks = [], yticks = [])
edge_mask = lab_ep>0
ep_edgelab = edge_lab*edge_mask
title(str(np.unique(ep_edgelab)))
imshow(ep_edgelab,interpolation='nearest')#lab_ep

subplot(324,frameon = False, xticks = [], yticks = [])
title(str(np.unique(pruned_edges)))
imshow(pruned_edges, interpolation='nearest')

subplot(325,frameon = False, xticks = [], yticks = [])
title(str(np.unique(lab_ep)))
imshow(lab_ep, interpolation='nearest')

# <headingcell level=3>

# First make vertex from branched points and then link to the neighbouring endpoints

# <codecell>

BPlabel=np.unique(lab_bp)[1:]
EPlabel=np.unique(lab_ep)[1:]
EDlabel=np.unique(pruned_edges)[1:]

G=nx.Graph()
node_index=BPlabel.max()
selected_ep_labels = []
for bpl in BPlabel:
    #branched point location
    bp_row,bp_col= np.where(lab_bp==bpl)[0][0], np.where(lab_bp==bpl)[1][0]
    bp_neighbor= lab_ep[bp_row-1:bp_row+2,bp_col-1:bp_col+2]
    
    G.add_node(bpl,kind='BP',row=bp_row,col=bp_col,label=bpl)
    #branched point neighborhood
    neig_in_EP=lab_ep[bp_row-1:bp_row+2,bp_col-1:bp_col+2]
    neig_inEplabel=np.unique(neig_in_EP)[1:]
    print 'bp label',bpl,' pos',(bp_row,bp_col), 'ep neighbor label:',neig_inEplabel
    
    for epl in neig_inEplabel:
        selected_ep_labels.append(epl)
        node_index = node_index+1
        #print 'bp label',bpl,' pos',(bp_row,bp_col),neig_inEplabel, 'current ep',epl
        ep_row,ep_col= np.where(lab_ep==epl)[0][0], np.where(lab_ep==epl)[1][0]
        G.add_node(node_index,kind='EP',row=ep_row,col=ep_col,label=epl)
        G.add_edge(bpl,node_index, weight=1)
        print 'bp label',bpl,':',(bp_row,bp_col),neig_inEplabel,' -ep',epl,'(',node_index,')',':',(ep_row,ep_col),
        #print 'edge',(bpl, node_index)

# <codecell>

figsize(5,7)
nx.draw(G)

# <codecell>

#Add isolated endpoints in the Graph
#label of isolated end points
lonely_ep = list(set(EPlabel) - set(selected_ep_labels))
#lep:lonely ep 
for lep in lonely_ep:
    lep_row,lep_col= np.where(lab_ep==lep)[0][0], np.where(lab_ep==lep)[1][0]
    node_index = node_index+1
    G.add_node(node_index,kind='EP',row=lep_row,col=lep_col,label=lep)

#print G.node[1]   
#edges dict
edges={}
for ed in EDlabel:
    edges[ed] = np.where(pruned_edges==ed),np.sum(pruned_edges==ed)
    #print edges[ed]
#get nodes of kind EP

ep_nodes=[node for node,data in G.nodes(data=True) if data['kind'] == 'EP']
print ep_nodes
#print edges.keys()
ep_edges_lab = (lab_ep>0)*edge_lab
#print G.nodes()
#

# <codecell>

figsize(5,8)
title('Isolated endpoints are added')
nx.draw(G)

# <headingcell level=3>

# each endpoints pairs are linked

# <codecell>

# get all nodes of kind EP
#get their label
node_to_label={}
label_to_node={}
for node,data in G.nodes(data=True):
    if data['kind']=='EP':
        #print node,data['label']
        label_to_node[data['label']]=node
print label_to_node

# <codecell>

#ep_labels from edge lab -> ep labels directly!
#
# node_index is not endpoint label!
#
map_epEdlab_epPair={}
print EDlabel
##plt.figsize(14,8)
for elab in EDlabel:
    #subplot(3,3,elab)
    mask = (ep_edges_lab==elab)# ep lab=edges lab
    nodes_pair = np.unique(mask*lab_ep)[1:]
    pair_weight = np.sum(pruned_edges==elab)
    edge_image = (pruned_edges==elab)*elab
    print 'edge lab',elab, 'ep2 lab',nodes_pair,'weight',pair_weight,'node1:',nodes_pair[0]
    ##title(str(elab)+str(np.where(ep_edgelab==elab))+str(nodes_pair))
    ##imshow(mask, interpolation='nearest')
    node1 = label_to_node[nodes_pair[0]]
    node2 = label_to_node[nodes_pair[1]]
    G.add_edge(node1, node2, weight = pair_weight, image=edge_image)

# <codecell>

plt.figsize(5,8)
pos=nx.spring_layout(G)
#nx.draw_circular(G)
#subplot(121)
#nx.draw(G, width=range(1,4))
#imshow(edge_lab,interpolation='nearest')

#subplot(122)
edge_labels=dict([((u,v,),d['weight']) for u,v,d in G.edges(data=True)])
nx.draw(G)
#nx.draw_networkx_edge_labels(G,pos,edge_labels=edge_labels)

# <codecell>