{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dijkstra and Fast Marching Algorithms\n",
    "=====================================\n",
    "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n",
    "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n",
    "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n",
    "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
    "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n",
    "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n",
    "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n",
    "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n",
    "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n",
    "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n",
    "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n",
    "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n",
    "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n",
    "$\\newcommand{\\CC}{\\mathbb{C}}$\n",
    "$\\newcommand{\\RR}{\\mathbb{R}}$\n",
    "$\\newcommand{\\EE}{\\mathbb{E}}$\n",
    "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n",
    "$\\newcommand{\\Ww}{\\mathcal{W}}$\n",
    "$\\newcommand{\\Vv}{\\mathcal{V}}$\n",
    "$\\newcommand{\\Nn}{\\mathcal{N}}$\n",
    "$\\newcommand{\\NN}{\\mathcal{N}}$\n",
    "$\\newcommand{\\Hh}{\\mathcal{H}}$\n",
    "$\\newcommand{\\Bb}{\\mathcal{B}}$\n",
    "$\\newcommand{\\Ee}{\\mathcal{E}}$\n",
    "$\\newcommand{\\Cc}{\\mathcal{C}}$\n",
    "$\\newcommand{\\Gg}{\\mathcal{G}}$\n",
    "$\\newcommand{\\Ss}{\\mathcal{S}}$\n",
    "$\\newcommand{\\Pp}{\\mathcal{P}}$\n",
    "$\\newcommand{\\Ff}{\\mathcal{F}}$\n",
    "$\\newcommand{\\Xx}{\\mathcal{X}}$\n",
    "$\\newcommand{\\Mm}{\\mathcal{M}}$\n",
    "$\\newcommand{\\Ii}{\\mathcal{I}}$\n",
    "$\\newcommand{\\Dd}{\\mathcal{D}}$\n",
    "$\\newcommand{\\Ll}{\\mathcal{L}}$\n",
    "$\\newcommand{\\Tt}{\\mathcal{T}}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\al}{\\alpha}$\n",
    "$\\newcommand{\\la}{\\lambda}$\n",
    "$\\newcommand{\\ga}{\\gamma}$\n",
    "$\\newcommand{\\Ga}{\\Gamma}$\n",
    "$\\newcommand{\\La}{\\Lambda}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\Si}{\\Sigma}$\n",
    "$\\newcommand{\\be}{\\beta}$\n",
    "$\\newcommand{\\de}{\\delta}$\n",
    "$\\newcommand{\\De}{\\Delta}$\n",
    "$\\newcommand{\\phi}{\\varphi}$\n",
    "$\\newcommand{\\th}{\\theta}$\n",
    "$\\newcommand{\\om}{\\omega}$\n",
    "$\\newcommand{\\Om}{\\Omega}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This numerical tours details the implementations\n",
    "of Dijkstra and Fast Marching algorithms in 2-D.\n",
    "\n",
    "\n",
    "The implementation are performed in Matlab, and are hence quite slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:3 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:3.\n",
      "WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:6 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:6.\n",
      "WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:13 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:13.\n",
      "WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:19 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:19.\n",
      "WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:33 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:33.\n",
      "WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:36 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:36.\n",
      "WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:44 overwritten at /Users/quentin/.julia/v0.5/NtToolBox/src/ndgrid.jl:44.\n"
     ]
    }
   ],
   "source": [
    "using PyPlot\n",
    "using NtToolBox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Installation\n",
    "------------\n",
    "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Navigating on the Grid\n",
    "----------------------\n",
    "We use a cartesian grid of size $n \\times n$, and defines operators to navigate in the grid.\n",
    "\n",
    "\n",
    "We use a singe index $i \\in \\{1,\\ldots,n^2\\}$ to index a position on\n",
    "the 2-D grid.\n",
    "\n",
    "\n",
    "Size of the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 40;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The four displacement vector to go to the four neightbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "neigh = [1 -1 0 0; 0 0  1 -1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity of implementation, we use periodic boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(::#1) (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boundary = x -> mod(x-1,n)+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given grid index |k|, and a given neighboring index k in ${1,2,3,4}$,\n",
    "|Neigh(k,i)| gives the corresponding grid neighboring index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind2sub1 = k -> Int.([rem(k-1, n)+1; (k - rem(k-1, n) - 1)/n + 1])\n",
    "sub2ind1 = u -> Int((u[2]-1)*n + u[1])\n",
    "Neigh = (k,i) -> sub2ind1(boundary(ind2sub1(k) + neigh[:,i]));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that these functions are indeed bijections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[13,27]\n",
      "134\n"
     ]
    }
   ],
   "source": [
    "println( ind2sub1( sub2ind1([13, 27]) ) )\n",
    "println( sub2ind1( ind2sub1(134) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dikstra Algorithm\n",
    "-----------------\n",
    "The Dijkstra algorithm compute the geodesic distance on a graph.\n",
    "We use here a graph whose nodes are the pixels, and whose edges defines the\n",
    "usual 4-connectity relationship.\n",
    "\n",
    "\n",
    "In the following, we use the notation $i \\sim j$ to indicate that an\n",
    "index $j$ is a neighbor of $i$ on the graph defined by the discrete\n",
    "grid.\n",
    "\n",
    "\n",
    "The metric $W(x)$. We use here a constant metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = ones( (n,n) );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set $\\Ss = \\{x_0\\}$ of initial points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = [Base.div(n,2), Base.div(n,2)];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the stack of available indexes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = [sub2ind1(x0)];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the distance to $+\\infty$, excepted for the boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = ones( n,n ) + Inf\n",
    "u = ind2sub1(I)\n",
    "D[u[1],u[2]] = 0;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the state to 0 (unexplored), excepted for the boundary point $\\Ss$ (indexed by |I|)\n",
    "to $1$ (front)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = zeros(n,n)\n",
    "S[u[1],u[2]] = 1;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a callbabk to use a 1-D indexing on a 2-D array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "extract   = (x,I) -> x[I]\n",
    "extract1d = (x,I) -> extract(vec(x),I);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step of each iteration of the method is\n",
    "to pop the from stack the element $i$ with smallest current distance $D_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "j = sortperm( extract1d(D,I)  )\n",
    "if ndims(j)==0\n",
    "    j = [j]   # make sure that j is a list a not a singleton\n",
    "end\n",
    "j = j[1]\n",
    "i = I[j]         \n",
    "a = deleteat!(I,j);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We update its state $S$ to be dead (-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = ind2sub1(i);\n",
    "S[u[1],u[2]] = -1;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve the list of the neighbors that are not dead and add to I those that are not yet in it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "J = [] \n",
    "for k in 1:4\n",
    "    j = Neigh(i,k)\n",
    "    if extract1d(S,j)!= -1\n",
    "        # add to the list of point to update\n",
    "        append!(J,j)           \n",
    "        if extract1d(S,j)==0\n",
    "            # add to the front\n",
    "            u = ind2sub1(j)\n",
    "            S[u[1],u[2]] = 1\n",
    "            append!(I,j)\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAABoFJREFUeJzt3EFOxDAQAEEH7b8tv9x8AGnDgnADVedEmltrDuNr770HABz2dnoAABhDkACIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgITH6QE+stY6PQIAXzDn/PQ/NiQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIECQAEgQJgARBAiBBkABIeJweAP6aOefTb9ZaPzAJ/C42JAASBAmABEECIEGQAEgQJAASBAmABEECIEGQAEgQJAASvNQAN915gWGMMa7revrN3vvpN15z4L+xIQGQIEgAJAgSAAmCBECCIAGQIEgAJAgSAAmCBECCw1i46e6hqqNXeI0NCYAEQQIgQZAASBAkABIECYAEQQIgQZAASBAkABIcxsI3c/QKr7EhAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACdfee58eAgBsSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkCBIACYIEQIIgAZAgSAAkvAMfaB7HM/N9XAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31190eb90>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imageplot(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Update neighbor values.\n",
    "For each neightbo $j$ of $i$, perform the update,\n",
    "assuming the length of the edge between $j$ and $k$ is $W_j$.\n",
    "$$ D_j \\leftarrow \\umin{k \\sim j} D_k + W_j. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "DNeigh = (D,k) -> extract1d(D,Neigh(j,k))\n",
    "dx = 0; dy = 0; w = 0\n",
    "for j in J\n",
    "    dx = min(DNeigh(D,1), DNeigh(D,2))\n",
    "    dy = min(DNeigh(D,3), DNeigh(D,4))\n",
    "    u = ind2sub1(j)\n",
    "    w = extract1d(W,j);\n",
    "    D[u[1],u[2]] = min(dx + w, dy + w)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 1__\n",
    "\n",
    "Implement the Dijkstra algorithm by iterating these step while the\n",
    "stack |I| is non empty.\n",
    "Display from time to time the front that propagates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAGgCAYAAACOmZ8zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAFoZJREFUeJzt3c+r/Nddx/FXpKBQ6EIQpKmCTdU2rRBsqIsWErBUcFEpuHDXPyXJn+LOhSC6UCwVWmgXLWkJqG21TQraiiC4KBRaEK6LWJN5f07mM+eeuXfeM/fx2M2Pe7+Tb5J58+F5Puc8c3d3dxcAoJ1fuvQHAADGDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaOp9l/4AI8888+qlPwK8p7u7Vy/9EXgPvjvo7D7fHa6kAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmmp5wAYAj++VvHbw+LW8cqFPwi+4kgaApgxpAGjKkAaApjRpgCeqNujPbt6hUV+aK2kAaMqQBoCmDGkAaEqTBngCan9Otg36M79Rnvj3+hMa9WNzJQ0ATRnSANCUIQ0ATWnSADdo/x7oQYP+THn4tfL6TqN++xmd+pxcSQNAU4Y0ADRlSANAU5o0wA3Ya9Cb/pxsGnQ+ffzP2G/UiXupz8uVNAA0ZUgDQFOGNAA0ZUgDQFMWjt2gt8rCjQ9buAE3Z3qhWF0klmwXiu0sHKs2C8kSh3KcmStpAGjKkAaApgxpAGhKk74BtUH/1ufK61/SqOHaLTfoUW8uz/34hV89ePxs/vvUj/fOHzt5KIdGfZwraQBoypAGgKYMaQBoSpO+MrU/J9sGnc+X1+vv0KihtdqfkzM06EGTrg369bx4+IYXXj94qFE/PlfSANCUIQ0ATRnSANCUJt3c3j3QSTYNOl84/jv3GnWiU8Nj2rsHOllv0LU/J9sG/Xo++Z6fMclFGvXbzzzd7yNX0gDQlCENAE0Z0gDQlCbdzG6Drv052TTor37wUwePX/rCN4/+mbVRJ+6lhoc0vQ93stygN/dAZ9ugvzV4z1GlUSfznXq/USdP+V5qV9IA0JQhDQBNGdIA0JQhDQBNWTh2YdMLxQYbldSFYl/OHx6+4YOHD/cWkiUO5YBzml4oVheJJcsLxUYbldSFYqPFZdMWNzzZLCRLnvShHK6kAaApQxoAmjKkAaApTfqRrTbo2p+TbYP+Uv7o+IfQqOFBLTfo2p8Hz8026NFGJfVn/uurv3n4+kuDzzHrAody3FKjdiUNAE0Z0gDQlCENAE1p0g+o9udkvUFv7oHOtkF/859LSPr4sU+ZTaNO5g/l0Kh5qmp/Ts7QoAdNerVBj+6Brg06XymvR6O+NFfSANCUIQ0ATRnSANCUJn1Gu/dAJ8sNenQP9KZB/2V5PZONOpm+l3qvUSc6Nbdh7x7oZL1B1/6crDfoTX9ONg1687i4lUb99jPX8X3kShoAmjKkAaApQxoAmtKkF0zvw50sN+hNf042DTp/Mfhz3/07LtCoE/dSc52m9+FOlhv06J7m5Qb9lcHnrM+N3nNEbdTJGTp1adTJGc6k3jTq5FrupXYlDQBNGdIA0JQhDQBNGdIA0JSFYxOmF4qVRWLJGRaK1UViyXah2Pd+WF4fLeN6159RF5Il04dy7C0kSxzKwXWYXihWF4klywvF6iKx5AwLxerj4XPfKK//weCHjuu44clmIVlyNYdyuJIGgKYMaQBoypAGgKY06SNWG3Ttz8kZGvRoo5LaoPPn5fUvlt9xvFEn99jwRKPmRkw36NqfB8/NNujan0c/M92g6+Mkmwadvxu96V2/4zYadTJ/KMelGrUraQBoypAGgKYMaQBoSpM+ojbR2kz3yu6wy35w+9SB0n6H9zBXtTHXBv3R8vqflZ//0+2v/NTHv3rw+HP5+4PHn80/HDx+6T/KP+tfDT7n3xw+/OGXDh9r0HTw5fpEaZXDe253bJrp4BCJWbXtjg672LXbmMvrL+f448Fzv/bSvx08fjGH/+yfLI9fzLc2v7L+zLNvlL/Pr+f448G/s6+Vf6/137v7pAGAowxpAGjKkAaApjTpCeXu43yxNNX9u48HnfoxGvVOg679OTlDgy79Odk26Pr3CR1sW+ThWpQn1ahfztzjrDfo+v5kvUHX/pz0bdCVK2kAaMqQBoCmDGkAaEqTnrDXqi7RqJMTOvVOg679OVlv0LU/J9sG3bUBwbvNNupkvlO3bdQvH39c+3Oy3qA3/TlZbtCbe99zPd8/rqQBoClDGgCaMqQBoClDGgCasnBswexCsuQeh3LsLSRLdjc8mT0sI1lfKDbaqORaFmrAMbsLyZLlDU82C8mS5cVkdSFZcsJispcPH+5tVJKcYaFYXRQ2em5yodg1f/e4kgaApgxpAGjKkAaApjTpMzqlVc1ueHKORj19WEay3KCvuQHBjPF/67dxKMfsYRlvv2exQY+a9BNq0JUraQBoypAGgKYMaQBoSpN+QKe0quVGnex26unDMhINGhbMHsrRtVHPHpYx+pnpBj34u3hKDbpyJQ0ATRnSANCUIQ0ATWnSj2x2v++9Rp3s30s9uw93okHDOV1ro57dhztZb9C1PydPq0FXrqQBoClDGgCaMqQBoClN+sIepVFP7sOdaNDwkGYbdXKGM6nv0ain9+FOlht07c/J0/7+cSUNAE0Z0gDQlCENAE0Z0gDQlIVjzcwuJEtOWEw2eVjG+HMAD2V3IVmyvOHJZiFZsruYbPqwjNFzkwvFfPccciUNAE0Z0gDQlCENAE1p0s2d0qr2NjxxWAZcl/H/k49/KMf0YRmJBn1mrqQBoClDGgCaMqQBoClN+sqc0qpqo9ag4frNHspxlkY9eVhGokGfmytpAGjKkAaApgxpAGhKk74Be61KA4Lb8xiNenYf7kSDPjdX0gDQlCENAE0Z0gDQlCZ9gzQgeHpmG3VyQqee3Id7/DlY4UoaAJoypAGgKUMaAJoypAGgKQvHAG7Q7kKyZHfDE4dlXJ4raQBoypAGgKYMaQBoSpMGeALG/fj4hica9OW5kgaApgxpAGjKkAaApjRpgCdq715qDfryXEkDQFOGNAA0ZUgDQFOaNABJNOiOXEkDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA0ZUgDQFOGNAA09czd3d3dpT8EALDlShoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaApQxoAmjKkAaCp9136A4w888yrl/4IV+1v89rB4z/OKxf6JLfp7u7VS38E3sNT/u54pfx///vl9fo4ST702+WJj5bH3zt8+KPvHz7+9uB31ude8/3z/+7z3eFKGgCaMqQBoClDGgCaatmkmVMb9Kd/ubz+c40abs1sg97052TboD92/M/8UH3i+6N3HaqfU6Oe40oaAJoypAGgKUMaAJrSpK9M7c/JtkF/4IXy+hvld2jUcHWWG3Ttz8m2Qe806WrTqJPdTq1Rz3ElDQBNGdIA0JQhDQBNGdIA0JSFY83tbVSSbBeKpTz+QHl5byFZYjEZXFJdXJWcYaHYaFHY4sKxkdkNTywkO86VNAA0ZUgDQFOGNAA0pUk3s9egN/052TTo/N7xP2OvUSc2PIFLqr159Nx0gz6hSf/sE4ePf2X04Sbd51AO3uFKGgCaMqQBoClDGgCa0qQvbLpBj5p0bdCj9xxRG3XiUA54TH+9c3hGcoYGPWjStUG/+f4PHzx+7hNvHTx+jEZd/y6S5E+e8PeNK2kAaMqQBoCmDGkAaEqTfmTLDXp0D3R5z8/K4/t0pNn9vjVqON1eg97052S5Qdf+nGwb9A/ykcM3vP/w4SUadbL9+3pKjdqVNAA0ZUgDQFOGNAA0pUk/oNqfkzM06ME90LVBf+f9hzHq+Re+e/BYo4bHNd2ga39Olht07c/JtkG/mecGf/C77DTqZL1Tbxp1snsv9S03alfSANCUIQ0ATRnSANCUIQ0ATVk4dkZ7G5Uk6wvF6iKxZLtQ7Dt5/vANZbHHJRaSJRaT8XQsLxQbHIaxulBss1FJtgvFRu856v3bpzocynFLC8lcSQNAU4Y0ADRlSANAU5r0gunDMpLlBl37c7Jt0P+S3xn8we+y06iT+Y6016gTG55wm2oPTc7QoE9o0rMNerRRyfRmJqdocCjHLTVqV9IA0JQhDQBNGdIA0JQmPWG6QY+a9GKD3twDnW2D/u7gPUcN7nVcvZe6NurEoRzchr17oJMzNOhBk15t0A9yn/QpGjbq5Ho6tStpAGjKkAaApgxpAGhKkz5iuUHX/jx4z2yDHt0DXRv0qFtPa7Dft0ZNB9P7cCfLDbr252S9QZ9yn/RbPy7veXb7OZY1aNTJ9dxL7UoaAJoypAGgKUMaAJrSpP9P7c/JGRr04D7p1QY9uge6/sx3f1zec46upFHzRCyfBZ0sN+jan5P1Bj26B3rToN84/D/5rdqxL9Cok/VOvWnUydXs9+1KGgCaMqQBoClDGgCaMqQBoKknu3Bsb6OSZH2hWF0klqwvFBttVLJZKPb64TKLzWKzG1lIllhMxrrlhWKDwzBWF4rd5zCMvdc3i8SSzUKx/FN9w+UXkiU9DuW41EIyV9IA0JQhDQBNGdIA0NSTadLTh2Ukyw269udkvUFv+nOyadB5o77h8o16+yn27TXqxIYnzKmdMTlDgz6hSc826FMOw5hu0LU/J9sGvWnS1U6jTm72UI5LNWpX0gDQlCENAE0Z0gDQ1M026ekGPWrSiw16dE/zcoOu/TnZNujXt285tNOok/WuNLjXcfVe6tqoE4dycNzePdDJGRr0oEmvNuj73Ce926BHvXm6SVfb/4s7HMpxiUadPEyndiUNAE0Z0gDQlCENAE3dTJNebtC1Pw/eM9uga39OztCgB/cKbxr0bpOutgXnVvf71qhv2/Q+3Mlyg679OVlv0KfcJz3doB+kSY9cfr/vSzTq5GHupXYlDQBNGdIA0JQhDQBNXWWTrv05OUODHtwnvdqgR/cfLzfoUW+uz/3nT8rrozuM91x+v2+NmlnT+3Anyw269udkvUGP7pNebtCnNOn/+VF5fVNi76Ffo07WO/Xwb2b3XupXp/8cV9IA0JQhDQBNGdIA0JQhDQBNXcXCsb2NSpL1hWJ1kViyvlBsdMDG8kKx0cKxulAsXy+vf7r8DgvJfqEuJEssJrtm04dljJ6bXCh2n8Mwpg/LSNYXio0WjtWFYvl2eb3+jttYSJZcbsOTWa6kAaApQxoAmjKkAaCplk16+rCMZLlB1/6crDfoTX9O1hv0pj8nmwadbwze8+7fcRuNevsp9u016sSGJ1dt9rCMwXOzDfqUwzCmG3Ttz8l6g97052TToDeP6++of8YjNOrkdg7luAdX0gDQlCENAE0Z0gDQVMsmXRtgbYSjjjhbVWtvGPXO0b11UwYdZXvoxmT5GPXj2pg3yuu/Xn7Hi4Mfqc/V5v/izw4efuzZ7xw8fj6Hj5PkY+W5382/Hv2Z539a7pMe/HvfNP1/PP76T8rjr/98+ys16Cv2vfVfUf+P3BzOsPq9MFK+K4ZddrWSjvpxbcwb5ciS95XfUfr98Lm9xy8cfpd8+Nk3N7/yI/nBwePn8ubS60ny3E9Lg64Nv46EvcfJ5r+/H5X7pO/TqF1JA0BThjQANGVIA0BTLZt0tbnrd9ARa6debdTJoFOfo0WV9rTcqJNtp66Neq9Bj5r0YoOu/Tk5Q4MeNenFBr1zRzlX5hwNsOrQqJNRpz7Dnby1U9dGvdeg79Okdxr0sB8vNujan5MzNOjB+of631+961yTBoAbYkgDQFOGNAA0dRVN+rVy3+or2Z4BXDv1aqNOTriX+loa9V6DHu2Fvtiga38e/cx0g679efCe2QZd/9vium12nr5Eo04udC/1AzTq2XueR89NNujal095z16D3vTnZLlB1/6c7O+E/vnBx9jjShoAmjKkAaApQxoAmjKkAaCpq1g4Vo0W+2wWk1lI9o7JwzKS9YViowM2lheKDTYzsVCMd9ssHKsGi31WF5ON/o/scSjHGRaS3Wfh2OJCsXttZjJ7WMboucmFYqP/1upz9fvm1cHP7HElDQBNGdIA0JQhDQBNXWWTHtnd8ORmG/XokxSTh2Uk6w269udkvUHX/pxo0Bw6aeOj6kYP5ThLo548LCNZb9AnbWYy26BPadKTDXrUpB/i+8aVNAA0ZUgDQFOGNAA0dTNNulpt1Ml8p95t1Ml6ixocBr97L/XkYRnJeoPe9OdkuUHX/pxo0BynUb9j26iT3U49eVhGst6gh/dJrzboUZNebNCP9V3jShoAmjKkAaApQxoAmrrZJl3NNupk/V7qUe25xL3Us/twJ2do0KMmvdiga39ONGjmnLTvf3WjjTo54V7qyX24k/UGXftzcoYGXfpz0rdBV66kAaApQxoAmjKkAaCpJ9Okq5Pun7yR/b5n9+Ee/czyWdCD9zgLmg6m76W+RKNOLnIv9ew+3Ke8Z3of7mS5Qdf+nPRt0JUraQBoypAGgKYMaQBoypAGgKae7MKx6qRNDq50IdnsYRnJGRaKDTYzsVCMa7C6kCxZX0w22gjpEhue3OswjNWFYqPDMBYXitVFYaPnun7fuJIGgKYMaQBoypAGgKY06SNmD+Xo2qhnD8tI1ht07c+JBs11mm7Uyc0cyjF7WMbwPbMN+pQmPdmgR036Wr5vXEkDQFOGNAA0ZUgDQFOa9ITVRp3Md+rdRp3stqjpwzKS5QZd+3OiQXMbnlKjvtd90qsNetSkFxv0NX/XuJIGgKYMaQBoypAGgKY06QWzjTpZv5d6tKfv3r3U0/twJ8sNuvbn5Lq7ELyXk/b9r66kUc/uw52coUGX/pw8rQZduZIGgKYMaQBoypAGgKY06TM66f7JC+z3Pb0P9+A9zoKG062eSf0ojTrZ7dTT+3Anyw269ufkaTXoypU0ADRlSANAU4Y0ADRlSANAUxaOPaCTNjl4hIVks4dlJBaKwTmtLiRL1heTjTZC2tvwZPqwjNFzkwvF6qKw0XNP6fvGlTQANGVIA0BThjQANKVJP7LZQznO0ahnD8tINGh4SNONOrnIhifTh2WMnpts0KMm/ZS/b1xJA0BThjQANGVIA0BTmvSFrTbq5IROPXlYRqJBw2Pq2qinD8tIlhu075pDrqQBoClDGgCaMqQBoClNupnZRp3s30s9uw/36HMAj+ekff+rB2jUs/twJxr0ubmSBoCmDGkAaMqQBoCmNOnmTrp/cudeavtww/VbPZP6Xo16ch/uRIM+N1fSANCUIQ0ATRnSANCUIQ0ATVk4dmVO2uTAQjG4easLyZITFpNNHpYxes73zRpX0gDQlCENAE0Z0gDQlCZ9A/balCYEt2+6USe7G57MHpYx+hyscSUNAE0Z0gDQlCENAE09c3d3d3fpDwEAbLmSBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCmDGkAaMqQBoCm/hcui15MccQleQAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x318d40810>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "include(\"NtSolutions/fastmarching_0_implementing/exo1.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the geodesic distance map using a cosine modulation to make the\n",
    "level set appears more clearly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAADy9JREFUeJzt3VGKXNkRBNAY4yV4duJZnAZrEJbRLK69E3kP7V+1MVRKJKmo5pzv4qlIvergfsTNX15fX18DAD/ZX372FwCARCABUEIgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVDhrz/7C/w/H/P7w8/869cvjx/0n0+zf/Bvjz/34esfDz/z5bfPDz/z6d+TL5R8+vvjz/z+8vHhZ/789Z+PH7Q4p398ffx/9/m3Px9+ZnNOH18+PPzM9fs0etbSvJO9mW/NO1mc+dLvN9n7DW/9fpPF3/Dx+5SX77+32wkJgAoCCYAKAgmACgIJgAoCCYAKAgmACgIJgAqVPaTLfkJy2zGa9BOS447RcE6XHaPpnE47Rotz2vpOo+ckydfHH9n6vxt1VJLRd9qa0+i3koy+09bfgslzkoy+09bfgtX3afakN5yQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoUFmMtZTreZdyveclb6vLyw5LoeNnHZZnp8+6ntPod3dYnp0+67I8m8xm/nnwnf6XExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVOm9qWLqFYboi+PIWhvEq5cPbBaYrpy9vYRivUn7SOZ3eUpDsra9eus0hed45Xd7mkCyuQ9+6zSEZ3ugw+cxbTkgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFSoLMZulV6nK4JPS6/DFcGXZc7NOa2VXt/5nC5LocnezLfKs8nezBvntFWeTfb+1m2VZ5PvKNB+JyckACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqVBZjt4pgk4Jiclx6HRQ5k9sy53ROp6XXJ57T5H3aesen21nXtqpulWeTte2zzzqncbl0afvsWnk2GRZov78864QEQAWBBEAFgQRABYEEQAWBBEAFgQRABYEEQAWBBECFymLsVklxUlBMbkuvkyJnclvmnBRek9vSa+OcVt+nw7JjsrdVdXPb79b22c1tv41zanyfRgXaF8VYAJ6UQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoELnTQ1LrfnxiuDD2wWmq5Qvb2EYr5x+x3O6fp8u2/fJ4prvrVsKkr116Eu3OSTPO6fG9+nT6ElvOSEBUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVChshi7VlIcrgi+LHNOVwSfll7f+Zwa36fLsmOyt756qxSa7L2bW+XZZO/dbJzT9fv0I5yQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoUFmMXSspTgqKuS1zToqcyXHp9YnndFp6XZzTadkxWdsWulYKTfa2qm6VZ5O17bPPOqfV92n2pDeckACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqFBZjN0qKU4KisltmXNSeE1uS6+Nc5oUXpPb0uvmnC7LjsnettCt8mzSOaet7bObW5Eb5zR5n74Mt89+ywkJgAoCCYAKAgmACgIJgAoCCYAKAgmACgIJgAoCCYAKAgmACp03NSy15qcrgi9vYRivUj68XeB6Tmtrx5OnndNp+z7ZW1+9dJtD8rxzurzNIXneOSWTz7zlhARABYEEQAWBBEAFgQRABYEEQAWBBEAFgQRABYEEQIXKYuxWSXG6Ivi09DpcEXxZ5tyc01rp9Z3P6bLsmOy9m1vl2WRvHXrjnLbKs8neu3k+px/ghARABYEEQAWBBEAFgQRABYEEQAWBBEAFgQRABYEEQIXKYuxWEWxSUEyOS6+DImdyW+aczum09PrO53RaCk32toVulWeTte2zzzqncbl0afvs9ZyS7y/POiEBUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVChshi7VVKcFF6T29LrpMiZ3JY5J0XO5Lb0+u7ndFgKTfa2hW5uRTanT48/M33WYXl2+qy8KMYC8KQEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFTpvali6hWG8IvjwdoHpiuDL2wXGK6fNafCNsrfme+mWgmRxffXWLQWJOS3O6fI2h2T4uxs96S0nJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKlQWY9dKr8MVwZdlzumK4NMypzmdz+myFJrsrfneKoUme79hc8ppeTaZF2i/lxMSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVKouxa6XXSUExt2XOSZEzOS5zmtP5nE5LocnaVtW1Umiyt1XVnG7Ls8ls++zsSW84IQFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUKGyGLtVUpwUFJPbMuekyJncljnNKedzuiw7JntbVbdKoYk5PfOcJs/6PNw++y0nJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACp03tSw1Jqfrgi+vF1gvEr58HYBc7qf02n7Ptlb8710S0FiTu99TsnkM285IQFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUKGzGDsoXY0KXMMVuluruUdFt+F32lrNbU6pnNPHlw8PP7O5cnpSDv7w9Y+Hn9lczT0pUZvT887pRzghAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQobMYOyldbZUdk1Hh8bQUmoy+02kpNDGnxTmtlRSHBcXLMuekyJkclznNaXVOkw3LyfD38g0nJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKlQWYyelq83th1vbQrdKocneVtWtUmhiTqtzWiopToqcyW2Zc1LkTMypcU6zwuvsN5wXxVgAnpRAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKgQuVNDaMW8Naa6GRvffXSLQXJ4prvrVsKEnNanNPW7QLTVe+Xtwts/u7MKae3MGz+7j6NnvSWExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUqi7GT0tVWeTbZWxG8VQpN9tZ8b5VCk7013+a0V+acrno/LXMOV71fll7f+5y2Sq+bc/oRTkgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFToLMYOSldr5dlkb/vsVik0WduqulYKTfa2qprTWplzWlA8LXMOipzJben1med0WXrdnNPw1/KGExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUqi7Gj0tXi9sOt7bObWyK3tqpuzmlrq+pWeTZ53jltlRQnv5Xktsw5Kbwm5jQpvCa3pdfNOX0eblj+lhMSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFSpvarhs3yeL69C3bnNI9tZ8H8/p9DaH5GnntNWa33yftm5hmK56f89z2lo7njzxnDL5zFtOSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVKgsxl6WHZO9FcFb5dlkb5Vy45y2yrPJ3jr06zmtlRSHq94vS6/TVe+nZc7FOW2VXt/7nH6EExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABU6i7GT0tVW2TFZ2xa6Vp5N9rbPPumcpltH17bPHs9praQ4LChell4nRc7kuMw5nNNl6fWZ5zR5n5Lh1udvOCEBUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVChshg7KV1tlR2TvW2hm1sit7bPbm6JbJzT1vbZ8/dpqaQ4Kyjell4nRc7kdk6TwmtyW3ptnNPm+5QXxVgAnpRAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKgQuVNDaft+2RvffXSLQXJ4jr0rdsckqedU+P7tNWan656v7xdYPN92rqFYfN9etY5nb9Poye95YQEQAWBBEAFgQRABYEEQAWBBEAFgQRABYEEQAWBBECFymLsZdkx2VtfvVUKTfZWKW+VZ5O9VcqNc7p+n7ZKitNV76dlzuGq98vS63ufU+P79COckACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqNBZjB2UrtbKjsnettCtUmiytlV1rTyb7G2ffdI5bb5PWyXFaUHxtMw5KHImt6XXZ57TZel1OqffXz4+/MzwF/yGExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUqi7Gj0tXi9sOtbaFb5dlkb6vq5jbNre2zm9s0G+d0WVKcFDmT2zLnpPCa3JZeG+c0Kbwmt+/T5G9vMvtb92W4YflbTkgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUqLypYXTjwdKa6GRxffXWbQ7J3prvpVsKksV16Fu3OSRPO6et2wU257R1C4M57a0dT25vYdj8+5RMPvOWExIAFQQSABUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUqi7GT0tVWeTbZWxG8VZ5N9tahb5VCk72V01vl2WRv5fT1nNbKnMNV75el1+mq99PS6+Kctkqvm3NaK70uzulHOCEBUEEgAVBBIAFQQSABUEEgAVBBIAFQQSABUEEgAVChsxg7KV1tlWeTte2za+XZZG/77FYpNFnbqrpWnk32ts8ez2mtzDksKF6WXidFzuS49Dqc02XpdTqn09Lr4vuUDH8v33BCAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKgQmUxdlK62twSubV9dnNL5Nb22c0tkVtbVTe3jm5tn918ny7LnLOC4m3pdVJ4TW7nNCm8Jrel10nhNbktvW6+T3lRjAXgSQkkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqVN7UcNq+T/bWoS/d5pAsrkPfus0h2VvzvXSbQ7K4Dv34fdpqzW/OaesWhus5ba0dT25vYdj83VW+T6MnveWEBEAFgQRABYEEQAWBBEAFgQRABYEEQAWBBEAFgQRAhV9eX19ff/aXAAAnJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACoIJAAqCCQAKggkACr8Fz9NsOLJYx2BAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31925bf50>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "displ = D -> cos(2*pi*5*D/maximum(D) )\n",
    "imageplot(displ(D))\n",
    "set_cmap(\"jet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fast Marching\n",
    "-------------\n",
    "The Dijstra algorithm suffers from a strong metrization problem, and it\n",
    "actually computes the $\\ell^1$ distance on the grid.\n",
    "\n",
    "\n",
    "The Fast Marching algorithm replace the graph update by a local\n",
    "resolution of the Eikonal equation. This reduces significantly the grid\n",
    "bias, and can be shown to converge to the underlying geodesic distance\n",
    "when the grid step size tends to zero.\n",
    "\n",
    "\n",
    "Over a continuous domain, the distance map $D(x)$ to a set of seed\n",
    "points $ \\Ss $ is the unique solution in the viscosity sense\n",
    "$$ \\forall x \\notin \\Ss, \\quad \\norm{\\nabla D(x)} = W(x)\n",
    "      \\qandq\n",
    "    \\forall y \\in \\Ss, \\quad D(y) = 0. $$\n",
    "\n",
    "\n",
    "The equation is then discretized on a grid of $n \\times n$ pixel, and a\n",
    "solution $ (D_{k,\\ell})_{k,\\ell=1}^n \\in \\RR^{n \\times n} $ is found by using\n",
    "an upwind finite difference approximation, that is faithful to the\n",
    "viscosity solution\n",
    "$$ \\forall (k,\\ell) \\notin \\tilde \\Ss,  \\quad\n",
    "     \\norm{ (\\nabla D)_{k,\\ell} } = W_{k,\\ell}$\n",
    "  \\qandq\n",
    "  \\forall (k,\\ell) \\notin \\tilde \\Ss, \\quad D_{k,\\ell}=0,\n",
    "$$\n",
    "where $\\tilde \\Ss$ is the set of discrete starting points (defined here\n",
    "by |x0|).\n",
    "\n",
    "\n",
    "To be consisten with the viscosity solution, one needs to use a\n",
    "non-linear upwind gradient derivative. This corresponds to computing\n",
    "the norm of the gradient as\n",
    "$$ \\norm{ (\\nabla D)_{k,\\ell} }^2 =\n",
    "      \\max( D_{k+1,\\ell}-D_{k,\\ell}, D_{k-1,\\ell}-D_{k,\\ell}, 0 )^2 +\n",
    "      \\max( D_{k,\\ell+1}-D_{k,\\ell}, D_{k,\\ell-1}-D_{k,\\ell}, 0 )^2.\n",
    "$$\n",
    "\n",
    "\n",
    "A each step of the FM propagation, one update $ D_{k,\\ell} \\leftarrow d $\n",
    "by solving the eikonal equation with respect to $D_{k,\\ell}$ alone.\n",
    "This is equivalent to solving the quadratic equation\n",
    "$$ (d-d_x)^2 + (d-d_y)^2 = w^2 \\qwhereq w=W_{k,\\ell}. $$\n",
    "and where\n",
    "$$ d_x = \\min(D_{k+1,\\ell},D_{k-1,\\ell}) \\qandq\n",
    "   d_y = \\min(D_{k,\\ell+1},D_{k,\\ell-1}). $$\n",
    "\n",
    "\n",
    "The update is thus defined as\n",
    "$$\n",
    "      d = \\choice{\n",
    "          \\frac{d_x+d_y+ \\sqrt{\\De}}{2} \\quad\\text{when}\\quad \\De \\geq 0, \\\\\n",
    "          \\min(d_x,d_y)+w \\quad \\text{otherwise.}$\n",
    "      }$\n",
    "      \\qwhereq\n",
    "      \\De = 2 w^2 - (d_x-d_y)^2.\n",
    "$$\n",
    "\n",
    "\n",
    "Note that in the case where $\\De<0$, one has to use the Dijkstra\n",
    "update.\n",
    "\n",
    "\n",
    "Once the Dijstra algorithm is implemented, the implementation of the Fast\n",
    "Marching is trivial. It just corresponds to replacing the graph udpate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "D[u[1],u[2]] = min(dx + w, dy + w);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "by the eikonal update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "Delta = 2*w - (dx-dy)^2\n",
    "if (Delta>=0)\n",
    "    D[u[1],u[2]] = (dx + dy + sqrt(Delta))/ 2\n",
    "else\n",
    "    D[u[1],u[2]] = min(dx + w, dy + w)\n",
    "    end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 2__\n",
    "\n",
    "Implement the Fast Marching algorithm.\n",
    "Display from time to time the front that propagates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAGgCAYAAACOmZ8zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3U1sXeed3/H/jV4oiRYpXyq8MinZpEJasi2pqW3ERdICg+lgkFVmmQCDbLPIxkCz8SKAbSALb1LAmyyyzcaLosDMIiiMzqAoJkESuNPEcSd2pNhMJMsSbcp6sax33y7SYvx8n8c8l7Zk/kl+P7v/vYdH555z73l08HteesPhcBiSJCmdL6z3AUiSpDYbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpLavt4H0NLrPbfehyB9ouHwufU+BH0C7x3K7NPcO3ySliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSktq+3gcgSbr7no3ni3rQ2KaPegL1ZdQXUJ9v7PP5eLbz2DQ6n6QlSUrKRlqSpKRspCVJSspMWpI2oBeQOS/g/TnUBxv7mEYo3Rsv6+HVsl5GKH2msc/jOK5TeP8ZM+s18UlakqSkbKQlSUrKRlqSpKRspCVJSsqOY5K0AbyEDlnH8f7iZFnvOIQNphs7nUKNffQulfVgBfVyvcsTp8v6JPYxh8/xLTuSrconaUmSkrKRliQpKRtpSZKSMpPehF5E5vO0mY+U2g/wm/1KY5unxsp6YhEbzKM+jHqmsVOusIHJTAKTmVQrbJytd7njzbJ+9K2yPniyrF++UX72X9W7jO9v4XuYT9KSJCVlIy1JUlI20pIkJWUmndyPkVWdaGyzsK2spzD28W9XMOH9nfL9Vxv7/M4WzoCke42LY/wF3v8Ks+KI6HFg9FHUzKhHyaQxdnrYMU46OC66kUkHx2fj353A/emvfov3mXtHxHacr620SIdP0pIkJWUjLUlSUjbSkiQlZSadzCvIXp7gfLuPNP6I4yOR+Uxhvt0pjFt86nf1Lh9fLo/jyS2UAUn32l+gfop58Zcbf3QMdUdGPUQmfbZfB93nY1DUV2JvUe/tXynqwfz5op65UAfIPWbSA9TMvTE2+6lfV7ussm+OK9/M46h9kpYkKSkbaUmSkrKRliQpKTPpzxnHPX8D7w+ewgtPomYOFVFn0sy3OJYRmXRgnGJExBOvlPW5X5bH/ffY3nHV0id7mXNzMx5mBs3ffWsbZNTnsKD0UswV9elqAPMImXQgk44ykz7Ux+LRETHXXyrqA5MYbM35wXdVu6h85XpZX26Mpd6sfJKWJCkpG2lJkpKykZYkKSkbaUmSkrLj2D3ETmIREd/kwu3/ERt8DfW/K8thY5KDU/2DRX0es+YPjpWz4i9cOFPUPU7MH1F1PhtggpRv/gO2x8LtdiTTVvYSfvtP4XdfLZbBiUpak5mgM9mbhw4U9e/jSFG/gbrVcewsfugXY19R74uLRT2DXqj8+9Y+Hl58o6gP7zpX/U3hev1S72pZP/WLsn4J959vbaL7j0/SkiQlZSMtSVJSNtKSJCVlJn0PcaKSiEYG/deo8f4fj32xqH/bmM2EWRRzIuZID/fLjOj4X9azmTw0/W75AiYgmMD23/hpWX+n2qO0Ob3Q6HvCX+kE+31gMYzqD5hRR51B/zZOoC53wkyak5tE1Dn1ytWy88nUeLk6z6EoJy/hZCgR9YQoNwKBPKLxw9eRUSN/jogIzIcygUWDjr9W1rwmz2zgjNonaUmSkrKRliQpKRtpSZKSMpO+i15BDlItlhFRj4NGBv36sYeK+ufx1aL+340BlMymmDMxRzoerxb1ciNX+uqxnxf10fhjuQFyowEyold+Wed0T27gXEj6JAuN17DWRb0ITkdGzcUyIuq+J8ygWf9LPFrUb6yUfx8Rcet19C4pp1CIDw6WfWLOHsW46qlyTHRExM3YWb32cWNxo6j3LF4r6gOXEEBHBNb1qBYNWsQ6HwuNXWxUPklLkpSUjbQkSUnZSEuSlJSZ9GfAubmfmMYGrYXbMRc3x0Ezg/6f8R9WfT8i4uQb/6Z8AWMG3zz2WFGfPfJAUXNcY8vuYx8W9UPLGEeNjOiJt+p9/HjZ+b218T2L3/1cY5sdnCb78Or1EHVrTDPHPbNmBv3aH9F/5ZUd9YHiXhFLqHEYt86VGfZrTzYmGS+71cTOuFnUe+NKUXN+8MHh31S77CFzRjeb2IH7zRwyaV6ziIjnN8j9xydpSZKSspGWJCkpG2lJkpIyk/4MTvCFR1DX02xX60FzbCPHQTODPvkz5M8REf8N9SuokY2f/Dr2wbHbUedG0xio+OCX/7GoeyexA2ZdEXFiuX5N2mg4q8DB1kbsn8Jll1Gf7feLurX2M19jbl2Ng2YG/U+N4+S94nXUnGP8Iuqoc+437iuPY98UMmfcS7i2wKE+A+iI2ZkL5Qs8nzjfvCb1TBAbh0/SkiQlZSMtSVJSNtKSJCVlIy1JUlJ2HPsMFrbhBU6izzoiTvXLLg31pPlld7RqohJ2EouI+C+oX79c1kuYRB9O7q87o80ceaeo5zHLwZH+74t6cR4z8zc++wI7rtxZ9bCklPqop/lCRMRU1x+V5Xl0bWIdEXEWvaXYkaxaLIOdN9lJLKLRmexVvF91jy3V62vErQPlcZz+2uoL/ozy2Wen0XGM5xPnm9ekjz/fSHySliQpKRtpSZKSspGWJCkpM+k1eBGTtE8xd2LNAfcRcR5hVFfONFKuxAw6/jPe/0/YB7Krxhz5p4+Ux8Hj5OdYnEEmzXMR9fl6EQtuPL1BJrzX1sYeHr3xxkaTqLHNEO9zkZvWojcXEQCvXMUPCj/BarEMTlQSEVUGHf+1tdHH9oGMeqGxDY6Dx3lxvPwco3x2nq/qnHe8P2EmLUmS7jYbaUmSkrKRliQpKTPpNWBm+rcryKhX8AdnozI4Vq4yUU0ujzGEbx57rNwBFsuIiHocNDPoo3if+zhW75LHweMcBFbL4GfluYiIFby2gWMibWHsATK8Wm/Tu4QXrq7+/t5+uaANF7iJiNiH1S2mxssf1AcHv1j+wRx2wMUyIrrHQXMZIe6D/0ZEtboFj5Ofg5+19dm7zmfgfV4TXrONxCdpSZKSspGWJCkpG2lJkpIyk/4MTmHu6am3sAHriFi4UA4ifLj/RlEfx7jFs0ceKOqTX6/n2a5wHDQz6K+X5eKR31S74HE8HOVx8nNUn7Xx2Xm+nndctDYg9qVYbnSuGLBPRvVH2H7+fFlHWUd09185e7Scy+DWOdwHyii4jeOgmUGP0J9lx9EyAe7u39L92Xm+qvOJ881rspH7v/gkLUlSUjbSkiQlZSMtSVJSZtKfAWe9fep3eOG39d/0Fsv6+F+WGy1jLdVqHtuv1fus1oPmXNzIjZhBfzV+Xu3z38avy+PEh+mVb9efleci6vMlbURMTDlldkTEgBkq5xFAPXOhTE0P9csc989/wvnzy3vFxalyTuzXnuSNYEd9oFwPmnNxz6FmBv3krWqXR6bK/itzmEScGXWVWV9oJMgd54+ZNa9JI+XeMHySliQpKRtpSZKSspGWJCkpG2lJkpKy49hn8B1MxvH4crngxhOvNP6o7PsRD02/W9RfPVZ34vq41uTzM0feKerTRw4VNTtmcKISdhKLqDuTPfRaeZzxC/wBPuv/YseZqM+XtBFxEp7j8Xy1zQn0+9rxJjYof6LRQz3XX6r2eRG9vNip9GbsLP/gobJ8474j1T5vHcCEJ+xxhcUyOFEJO4lFRDwa/1Jug4mQWLNjWY/nKiKCr6G+hfNd7nFjT5zkk7QkSUnZSEuSlJSNtCRJSZlJ30VPIvc498s6qxpM4YXxsjwafyzq3cc+LOrpxrD8eSQwnPSAE9pzsQxOVBLRyKD/ARv8rCzP/7KseS6kzepU47WTl8r6US44g74pmJckDkxiBxHx8GL5u70RY6se1864WdT7puoVNk5/rQzDV66WN6ip8XLlCvZvYZ4cUWfOvL/w/nOAJ+v1apcRJ1HjfHIXrWuyUfkkLUlSUjbSkiQlZSMtSVJSZtL30N83Xvsmst0JbnC1LB9aLrPhB7/8j9U+j/R/X9TnY7qoB5h9fuFCORiyWiwjoh4HjQz6Mj5H67NKW8Ezjf4Xcxg7fRCZ6gT7pkyiRl+ViIjDu86VL2Bs9VjcKGrOqTBo9GdhxnxxvByLvS/KHJv9W/j3EXUmzQz68Gl8jtewg8bCRMypL+N88k9a12Sj8klakqSkbKQlSUrKRlqSpKTMpO+h5lzVN8qs6hs/Ld8elMMSq8XNexwvGBGL82XGvDiDCXi5QDrHbLYyIMzFzXHQzKCdl1v6V9/C7+Fl/O7/Cr+5HjPoXd3/xuHrZba7Z/FaUXflyRER5zFAm/OBd+XarUyaY6ercdDMoNknhu9HxBDn65dl/F6d783EJ2lJkpKykZYkKSkbaUmSkjKT/pwxu/0O3n8F830/wfy4kdfEPGqOwWTOzX3+rt4l14N2Lm7p0/sV6okLZf1Ua64Cuo4acyocuFRmv4PDvynqQ/06P/6smfTMBXyQaKwHzbm42QemK6OOiF/hn+H53Mx8kpYkKSkbaUmSkrKRliQpKRtpSZKSsuNYMuyg9ePlsiPZCXToiohY+KeynkLHsRV0HDt1p6xfbRyHk5NId8/38XvajgU4OM/IV9hJLCJ66CgWmCOE62f00E9sdqbu5DU7Xb42xEIfPf4bvP/U86NEsOMYJ2BCRzJOVMJOYhER/wM1z+dm5pO0JElJ2UhLkpSUjbQkSUmZSSc3Sjb87J0y3+ojN2LE8/wWynOkjJ7Bb/AHyKgvN3LZp35R1hMdi/FUa1/MNA6kX5bVQh/MwXlco2TSmDzpMjJqLpbRmqhkK2XQ5JO0JElJ2UhLkpSUjbQkSUn1hsPhcL0Pgnq959b7EKRPNBw+t96HoE+wme8dLyG3Po73FzHGecchbDDd2CkX48E+qrHYzMEb8zbcQhZ+Evvg+hrf2kJ586e5d/gkLUlSUjbSkiQlZSMtSVJSjpOWpA2A2e0LyKgXkP3OoT7Y2Od0xzjpIcZJL2Oc9JnGPpdQn0LNMeJanU/SkiQlZSMtSVJSNtKSJCVlJi1JG1BXtvssMutBY5s+MuYJ1JexPafuxhLWEeHaAHebT9KSJCVlIy1JUlI20pIkJWUjLUlSUnYck6RNyA5cm4NP0pIkJWUjLUlSUjbSkiQlZSMtSVJSNtKSJCVlIy1JUlI20pIkJWUjLUlSUjbSkiQlZSMtSVJSNtKSJCVlIy1JUlI20pIkJWUjLUlSUjbSkiQlZSMtSVJSNtKSJCVlIy1JUlI20pIkJWUjLUlSUjbSkiQlZSMtSVJSNtKSJCVlIy1JUlI20pIkJWUjLUlSUjbSkiQlZSMtSVJSveFwOFzvg5AkSTWfpCVJSspGWpKkpGykJUlKykZakqSkbKQlSUrKRlqSpKRspCVJSspGWpKkpGykJUlKykZakqSkbKQlSUrKRlqSpKRspCVJSspGWpKkpGykJUlKykZakqSkbKQlSUrKRlqSpKRspCVJSspGWpKkpGykJUlKykZakqSkbKQlSUrKRlqSpKRspCVJSspGWpKkpGykJUlKykZakqSkbKQlSUrKRlqSpKRspCVJSspGWpKkpLav9wG09HrPrfchSJ9oOHxuvQ9Bn8B7x2fzbDxf1M/Hs+t0JJvTp7l3+CQtSVJSNtKSJCVlIy1JUlIpM2lJ0tq8gDy5j/cnGn/D19gg/BT7vIz3WUdEXED9jLn2Z+KTtCRJSdlIS5KUlI20JElJmUlL0gbwI+TDs3ifGfRImfRYWW9Hi3D7dllfvoG6sU9m0n+H434b73/XzHpVPklLkpSUjbQkSUnZSEuSlJSNtCRJSdlxLDlOeN/q/LEX9W7U11BfQd3q/OHE+tLnh53C5hrbsKPYAPU0eor1JrEB64gIdByLXaivl+UEOo7FpXqXQ7y2jJ5k57E9J0xZqne5pTuX+SQtSVJSNtKSJCVlIy1JUlJm0uvsReQxzJmq3KmxD+bUu7eV9bU7Zc0Mermxz0dwXMyRnt7CGZG0Vlz8Yg7vL6Cew284ImKq62bA2UtGyaTHUbNFwGQmcRV1I5Pu4bUBMukBbjizuLlM434VEfFSR269mRfx8ElakqSkbKQlSUrKRlqSpKTMpO+hHyJHiWhkT6hnkUVNTWED5k4REfeh5thHjm38oCwf5Yz4EbGyUtZvIyfipPmn8Pff28QZkdSFfU0ewfvVfYBjnB9s7JQZNGtm1uysMkom3TFOepRMuur0wg4tyKSn8Dn6f6p32cc9irdBnu/N1GfGJ2lJkpKykZYkKSkbaUmSkjKTvot+0pFDRUQsIi+emMEGXbkTM+qIOldizRyJNfLniIgp5kao586W9Sxyb56LiIhvb6KcSPr/OO92RP3bZz2YxwvMoHlfaL32AOqu/iufVybNPi48zndQ415SzTkeEfPIqfe8VdY7sD2vyUae+9snaUmSkrKRliQpKRtpSZKSMpP+DDhW+Djen2/lSodQd2VRo2TSHA/ZlStxHGMjk64m9EZuNIGM6InTZd3H9hH1+fqbDZwTaeti3snffUTECfY9YQZ9GHXXfSGivjes8V5xq7EY/ZXJ8mZxE5Ms7MQkC3svlTeTHa3F6Hk/4b2kKztvHCez8wHO725k1NvRR6bVb2Cj5NQ+SUuSlJSNtCRJSdlIS5KUlI20JElJ2XFsDV5G54MTeH+wiBdYR0SsdRKDrslNIupJCromM+EEBOzY0XqNHcE6OqnM8xgiYs/Jsub5/OsN0pFDW0vXYhnsJBbR6CjGewE7jnXdFyLqzmW4V1welFN6vLet/FFeib3VLvlaZ8exySur1hER+2fKnmMT52+VG3RNstK4d1SdYbEQEfuanUBHsttcZCg2zqIcPklLkpSUjbQkSUnZSEuSlJSZ9CqYmT6OHGSqaxb9T5NJo76FyelXJuvZ55krfRi7i3pPXCvqvVHmSFOX6lnyd3ASfC7E/ilyJU5A8PjvyvrlO2bUyqfrZ17lzxH1b/9h1F2ZdSOTvo5M+uz4gaJejkFRX4x9Rf0+6oiID6p7x56i3hMfFvV9uHfcHxerfS5vK1+bnjlf1DOT54p6133YQSuTZkvV0XIxo37k9XobJOXxAu73zyS5//gkLUlSUjbSkiQlZSMtSVJSZtL/Dxd/iKjHQVcZNGfWP4q6lUnjtSHGS57tl2HveeRMy42B0nUmvXquxEx6erIeKD2YRI40Xa7k3mMmzRypMXaU3zbOs3/itbJ2QQ6th5fwvVvA+4OuxTJar3Vl0KgvzHNgcMQ7GBh9tqNewS+slUkzt+4aJ70PGXQrk57CChu8P10bL+9PDxwtJ2Hoj3FFoBHcRn2nLAeNcdILGEt9od4kBZ+kJUlKykZakqSkbKQlSUpqy2bSPxlh4fZqLm5m0syguwZURsTVxfL/RUtjc0V9GhP0Mmd6r0pyIy7G/UV9JcqBh3vjg6LeF+8X9f5qlfaIGUzWfbF/uqjnxpeKenzXR+UORvlmIUdibnQcc33zmkVEfNucWp/Rjzoy6Dn2v+AYZs6pHbHmDPrd+fI3y/tA6zXW59FfheOmmVFH1Dn1NfRn2Y3+LMygmT9HRExH2Z+lq8/MDeTgd+bLe01ExBdxD6swk2aszfULImIO00NcQCjN78V31+le45O0JElJ2UhLkpSUjbQkSUltmUz6hx1rws5zHeeIzhypqrHTy4/sCDq1rUy8/hBfKuquTJp1RJ01dY2TZo7E/DmiHj/JXOnGWJkjLTxyqqgnqplxo86NOHYRudE86gv1YVbX9Xtm1FqjOdaYo7/XteZ7a+3njjn5OQ6av/ul6qjq17ruFaPMsVBl0jeQSY+tnklPNxaj572D84PfjJ1FfYeLQzdsmy8HNfdvIHRm5sz6cr3PHjJpZtTLGGu9XnySliQpKRtpSZKSspGWJCkpG2lJkpLaMh3HOEHBIheAGGVCgo4JCjhRCTuJRUS8gdXf/4AjW2vnkIh6gpOLl8qOG/smy84enLykNckBO3+wM9ptdvZAeWTxjWqf49cx4Qk7d6DjBudJWKznTYi3GxPnS5+EE1RERMyinhrgBfa34k+w1ekU95PrqLlYxigdx97Ca0u4IZ2NB1CX/8byeX6wiI/OYWUcrJeBW0ksHyh/tO8P6kU7OJlS570DtnF1jIgYQy/TPYfeLOpdnOuE95LW6hm4n/C6z6Kj6npNbuKTtCRJSdlIS5KUlI20JElJbdpM+kXkB3N4f+IuTEgwxMLuXCyDE5X8+bUygz6Fbboy6jNvN8LzM1gg/r2yPLe/XCHg3MHyw1+crXOlteZI25EjjY3VYfGjh8scqceciDXmSZio502IOSzczuv+tJOb6GOYP0dEVEktM2jWD6BuZdJ47ez4gbLuyKRbC2wwg676r9wo/+bSqfLfjKXGcZ5DjUwaXVPiowNlhn1uDjfBiLi2UN47bo6xA1CJ946dcbPahhMy7R4v68Mz+CC8lzT6s2AdkOo6D5BJt747nwefpCVJSspGWpKkpGykJUlKatNm0syZZhmpdo19bL2GTPpsv8x6R8uV5oqaGfQbcaSo3/0/+EdfbxznEmrmTIimYq7MsM8c5UohETceW1uOxMxob1yp/mZfvwy8Zh9EcMTMmee/Xh8gZv9U1oMkk+IrhxfQR6Hf2GaaL7LmzYTTCjS+l5cH5eI6y9hJvRjG9Krv//m1MgyvMujX8UPnveJU1M6gRn+W2I/6IGqOT46IS7dxHEfLcif6q/DecV/j3sH7Cev9gzJ0nljGAj/1VBD1dUUGze9FH7crfrciIp65B31gfJKWJCkpG2lJkpKykZYkKalNk0k/i3ygipG6cqRGrsTXbmF8JBdVZ47UypW65uitMuhfYwevNY6TWRNzJuZIzJ2wfnpExLsI4MceWz1H2ocBlqwj6gXipx8oQ54dXdekcY14XQfItfm9eN5x01tKV9wcEdGbxAusJzp20sg739uG+fQx4Jjz5TOz5r0lor6fVOOgmUF31RFr78/Cn/Xtxj7h0naMEX+M/VfKYLt177gfr/F88nxPTOGDtC48ryuuO78XzKRbu7wXfJKWJCkpG2lJkpKykZYkKalNk0kzXqjiy64cqTWODjtZmSxDimVswHWd22Mdy9equbiZGzGDZkbd2qYrk26MbaxgOvAz+8rj3DdbZkRTmBx3upoYtz5fK5PlwMQD01gEdpRrhOs6zfm+G3+irYPXv/l96MqkO+pbjZ1eib1F/T4yVNZ1Rl13wKjWg17CBuyb0nUvaf3N9TIvjnPlPNwj3TvYqpTLS8fy/vJz3D/gevfsNNN9/ni+b02UmfQOXsOINV/nkb5L94BP0pIkJWUjLUlSUjbSkiQlZSMtSVJSm6bj2F7UVaiPzgsxjrrVCwAdB9g5gfXFuL+o2Rkkou5cFmfQQ2sJf8COHa3OH2fQ2SP+Ge8/jvfRGYTnJqIxsX55nO/Nlp+Dn5XnIqL7/B2YRMcxXhNes4jq2Pkn/F5oa6k6+7TWjeFr/J511Fcm8RuO+rv9QXWvWL0jFOuIiI/O4R/mxCPsMLqEurXAxnV28HwV75/APtB5rXXv4KFjQhR+jvcH3Z+d54vns7q34Jr0xxszNnVdZ3wv+N2ZKOd3umd8kpYkKSkbaUmSkrKRliQpqU2TSe9mvQ0vdOVOdaxUbfMh/pUPke1eQUDD9yMiLl5C3sJx+105E+uIqDLo+O+tjT62j3/fvU8eB46Tn+PDydXPRUR9Png+O69JK5PGdeV1332n8TfatLigCm9w21t3PH7PuA3fR32zurnUr/G7z/ev4f1rN+p7R7XuBOuuewknKomIKoOOnzW2+fg+voZ/o3GcPI6O4+ZnvTZW73Ot57O6JrsamfQarzu/O9sbmfS9WNDHJ2lJkpKykZYkKSkbaUmSkto0mfQ11swimR9cRd2ILLjNHvwre2L1xcv5fkTEvskykDm3HytEcJF1Lo7BOqIxDprw/ij75HFg3DQ/R9e5aG3D89l5Tfh+RHVded35vdDmxgzwp8gIb99u/BG/Z9yG76PeWd1c6tf43ef7u/H+7rH63sHuLNVwYs5twN9wKz/mOOgK3t+FffDfaB1Hx3Hzs/JcRKz9fFbXpHV/X+N153en9VW6Gxk0+SQtSVJSNtKSJCVlIy1JUlKbJpO+gvoyN2BEynyz+oOIwFTSe/GvsN4X7xf1VKxUu9yP184dnCk3mMPgPI45bC66jpyI46CZOR9DvdDY5Rzqg2VAw8/Bz8pzEdF9/ni+q2vSyqRxPvgn/F5oa+H34XJjbGs1BzO/Zx313kt14Ll3svzm3VfdK8o+Hfd31BERywfKf/ijA5g4gL9z7qJ17+Bc3BwHzQya94q5xj55HMitv4DPMcpn5/ni+azuLbwmrXtH13XG94LfnVaTcS/4JC1JUlI20pIkJWUjLUlSUpsmk2Y+sIz60Qt4gXFxHR9XO5m6VIam05PlBsxpZ+Jstctq3eXZctDgmaOL5R+0xvcRp8nmXNzMiJgrMaOOiDiKXcyeLmp+NtY8FxER0zihPJ/VRRvlGuG6chefV26knKpMurUR+0Kssd7R2Ckz6a7clX06+FuJqNddPjd3uNyAmXNrIC/x3sGx1BwHPYca94mI6MytpwflGtbVfaHxQ+86f8ykq2vCa9h6raMe6bt0D/gkLUlSUjbSkiQlZSMtSVJSNtKSJCW1aTqOcWLzRzCx/gr6IkyxX0bdT6N6bcc7ZT2YLDtAsPPUxWpm+fo1Ll5+47FysfJ348FyB1yoPKKe0J6LvX+Kzh9ffOxP+JOloj4Ua+tIFhExiPJ88XxW12CEa8Treh7v34sJ77VxsL8o64iIIToI9bom1enqhBoR+2fKF5e3dXVs3QSOAAAJ8ElEQVQUK7+5rXvHFfTyurZQ3jsu3W6tdvExrbs9/xlOnsR7S1cn1IjqfjK5UN6QeG/gfYHnIqI1WVJ5PvffwUXgNWld+I7OZfxejPJduhd8kpYkKSkbaUmSkrKRliQpqU2TSRNTjbfvlHWVSdcRav1aGdPGzHSZSlzslzntldhb7ZIZ9O3Y1viH/9XYY+Ws7mf2Hao3OtixKEeVK5UzpHCikog6g/5SnFr1fWbUrCMiZi4gxcH5rM4360Ymzetap1nayp5Bn4S/Q1+ViIhlfC0HDBv5pXqAO6j/3Ynzt4p6eqbcCe8NrD8Y4d5xc6zsv8Is+NJ2ZNScuCSi7q/CtS2YWXf1b4k6gz40xv4r76BmRl2fUG7D3Jrnu3NipIj6unJipAurvl19t+4Vn6QlSUrKRlqSpKRspCVJSmrTZtJPd2RRc8g7J5iPRkRMo+6XZQ/13PhSUd9gZhTdGfT2KEPWPfFhUe+brRdEf28Wi3ZcKoOkfZMYUzjCQiDMlLsy6i/FH8rtb5TbR0T03sQLPOcd9eVGvwH+K7zu0se93XiN0eSga3w+x/dPRQ33hpnJMqe9Nl7my1XeHDurXXbdO3aOlf1Xzj5W3juW9w+qv/no3Hj5Qkcm/YUDV4uai2VEtOZMKE/YXLxV1CP1Z+E+r2IyiK7+LLxmEZ3zMFT9mhq7+Dz4JC1JUlI20pIkJWUjLUlSUps2k6ZTqGfL+CaeqGOQOmuaRI04Z3zXR0W98Aj/1QjGSp0ZdMfi8BERKzjQDyfLfIv75D5amTRfY07EDHrhTvlZx0+W5yIiIk6ifqujxjU5iWsWUV9XaTXfbfRZ+Cn6q8wijJxi3xT+XJA/R0R1r9iFMcoPHC13ciPK/it3OvLniO57x974oKjvH9T9Wd4flKHztRvlvWP3WLnP+3E/mm6MaeZc3F33kq45FyIiHsA+dnET1l0ZdUSVQa/gujODbn13Pg8+SUuSlJSNtCRJSdlIS5KU1JbJpL+HPOEnyKH6jcxiHpkzM+jgMGiczYnAfLIRcWTxjXIXGNu4N64UNTPp1lqrF+P+oua6s8ym9sX7Rc1x0xEj5EgYB11l0L+rdlln0h31W7gmrV3yukprtYR6GvPB9zFev8e+KRONnXbcO/pj5fz5d+ZbnWJK25BB74ybRX1fx71jfzWpf8T7GAh9bQyZdKyeSbf6yPAexbm4u8ZFtzLp/lvl+VrrHAutTHqIbZZw3ZfqP1kXPklLkpSUjbQkSUnZSEuSlJSNtCRJSW2ZjmP07REWg9+DjkyDjo5ildv1S+PXyw5Wjx4uV53Y1199soDlatWPesF4TtZfT3JQdjAZaUKCC+WS59ViGV2dwiIiXkeNnmDn8Te/xea8ZtLdwEkqXmKn0vKrH/PslMROYhERu1B33Cu+iM6d2+Y5s0/EWJSdTLt+1+zkxU5iEREX8dpN9IbdiX+TndH4b0R0T5bEmhOVVJ3EIrrvLx0TIVUdySJiCdeVEyOt1+Ql5JO0JElJ2UhLkpSUjbQkSUlt2Uya/qaRP7yMbOpxZKjVWu/MoBsLQkS5Znr0kIvMPli+MP1AWa9M1qPy60x6d1HviWtFzexq6tKlap87uEh612QBzIRamTTO3wrqV7E5M2np87CEmutn7MF3veqrElEtpLPWO23/Rp3L7jlUdgTZPb56Js28uZVJf7DG/iycMKWVSXdNwDRz9VxRV4tlNPLjzgyafWRQn68j/iqDXmr8sxn4JC1JUlI20pIkJWUjLUlSUmbSq/hnvoAJ2E+8VtYDZtDInyMigvEvMmkOWd6BYdEHpuv8+MAkXuO4TR4Hd1EPk65fYxT+KTJpjoNmBs3z/f0k4xS1tTyD792L6JuyA9vvbuSdrTU3Cuy/wrpx79hVDqWOwzNltrt/UI5Pfm9b2WuGfVdar3WNk2buzToiYv+d8jgmzmOhId5LRsmkuzJo3FsuY/vW4jx8jdc9C5+kJUlKykZakqSkbKQlSUrKTHoVzCh+gGyqHH0ccRy5yHwrk+Ya6cx+Z1Bzqu5qcHbUARjnDeaQy8sdx9Q6rq5MGrnSW41F1jnumbUZtDJ6Gt/LH+E+sL0xH8IJZKJVRo3+LdVv9FP0Z5lYLrPfiakys741UdYREVcmy5tFZyZ9qTzQHbyXRHTf43hv6Lq3RNS5NTJpZtCv4pq0Mmle16x8kpYkKSkbaUmSkrKRliQpKTPpNejKTH+CrOpCI5ddRF4zwbyGGfQomTTHRXeNk2Y9SiaN+jI+28kRMiDXg9ZmwHWGmVFHRNzG7+ERrKXeOadCK+vlnAr83fLegEnHd0zWu+yPIwzfhborK6+nbeg+Tq4L0JVRR1Q5Nefi5v2GdZa1oT8Nn6QlSUrKRlqSpKRspCVJSspGWpKkpOw4dhexY9QPGx1K3kaHkTl0gJhFB4mpjs4gERFxH2ouQs9OKpiov+roEREr6OzxNiZfWML2XED9exu4o4a0Fq1OSVyUA0tMxAJ+93PogNVrdchiB6zzqAeoOYNKo+NY1cm0ayKkUTqOsdMbj7OjU+qwMZnJEu5RvN+wo9hGmahkFD5JS5KUlI20JElJ2UhLkpRUbzgcDtf7IKjXe269D+Fzw+yKsRJrzm0SUUdPu7eV9TXkyYyMGBFF1DES682U+azVcPjceh+CPkHWe8cL+J3P4f0F1HP4DUdETHXdDNhfhRn0KJk0eyndRj1KJt2xEAhvOCu4uSxx8ZGoM+gl1FwMKatPc+/wSVqSpKRspCVJSspGWpKkpMykk3sWWVa1eHxE7EW9G/U11FdQt+byf36DZDzrwUw6r4167+AiHXONbWZRVxE1MuneKJk051ToGifNORcamfQQry0jk2b/lrdRL9W73NALZHycmbQkSZuIjbQkSUnZSEuSlJRzdydnNixtfqNkrsytmVH3L6xet/qzTCCT3o4W4TbGSV9GJt3qz8Jh0ayZQW+WvPle8UlakqSkbKQlSUrKRlqSpKTMpCVpA+jKbjk/OKfybmbSyJi3o+bU3cygR8mkN8q82ln5JC1JUlI20pIkJWUjLUlSUjbSkiQlZccxSdoE7kYHLS7o42RK688naUmSkrKRliQpKRtpSZKSMpOWJEWEGXRGPklLkpSUjbQkSUnZSEuSlFRvOBwO1/sgJElSzSdpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpKykZYkKSkbaUmSkrKRliQpKRtpSZKSspGWJCkpG2lJkpL6vwHeRwzL5T0fAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x319220f90>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "include(\"NtSolutions/fastmarching_0_implementing/exo2.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the geodesic distance map using a cosine modulation to make the\n",
    "level set appears more clearly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2w1nW57/HPLy01UckdpltSwIcDCAbSmHuiovIBzhFrq6nYcqCEhB6UqB3pIlvbWBrtIjQLFCwcV4Ik7hLOwYceqGhSJwQFWR4lwDZuSdqGiqml/s4fnZmYppnPxz3nnK7mvF9/f+b6frnve92X94zX72ratm0FAMBf2Wv+2hcAAECiIQEAiqAhAQBKoCEBAEqgIQEASqAhAQBKoCEBAEqgIQEASqAhAQBKoCEBAEqgIQEASqAhAQBK2PuvfYG/5Ev6hM3MHP41X2hjV3jiNJs46AX/DNpdvzjMZtaMji6kHweZzoE+s3NLH5s55sVHg9Okpwcc6kM7FgeVtgWZMUFG0oAg91kf6X+Rfw1O052+kKTh2mAzGzTcZu7UaTaz/bpjojvpi0Fm2+oglGQGBBlJh06ykYO27bCZR/fxr0G/QbuTG6l7q8+8K6gzeo3P9H3rE0El6el9myA130eGddnInA3+u1eSPjP8Wh/a8Oqf280vJABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACU3btq/+fxb/v6zxIxrhjNHno/NGt9+3mZ9eearNzO/0Zx2QXEhSxxyfmfSZb9jMjeP8jJXu6AluJEmPBZlgjuHaA22k42MLg7OkuZphM/2uCGZQlvjIpoeDCyl7lY4MMkMHB6EJQUbSzsv9PNoMzbWZnq9P8Yd9/JnkSpKCWcLklRrbYSMTVwVzOpIWf+mjNtMz09d5NjhrWncQkvSOy+6ymTXNyUGlf/aRYFZJkq7fcIHNTNFNUa098QsJAFACDQkAUAINCQBQAg0JAFACDQkAUAINCQBQAg0JAFACDQkAUELNwdimK0j5odfT21uj81ZMOcdm5izydZLRtFG/CkKSmt+/4kNHPxlUSgYCPxhkJC32i9DmTbzIZi657nqbeeaS6EZa/KLPPBXUOTzIDAkykjRkL5/pfTnIBGc9HmQk6eAgM2kfnznwap+5+qKPBKdJ02+8zocmJcsjvx1kggFxSdp8iI20r/P/Hb/2CH+UH8f/o5mTfWb8wmU2s7I5OzgtGJ6VpLFdNtKuykrtiV9IAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEoo+qeHXNjO6fdBmfjrFrx2XpO7gKQzBsLR+2/p1y0O+si2oJOnTtwQh/zppzMU2MvFH4XrnHr/e+UG/2Vi3BWcdHWQkqWNYELrMRx6Y4J9CcauSSXfpPr3NZk7UvTZztvyTRt6yJHmSgaQrfaRno89sDo46M8hI0vHBhutJHd+wmRvfHTyFYfU1wY0k6U0+8uVzbaT3UwNs5g1NsuxeCr6e1Bl8Qb1jYbIK/fjgNCl5Okbbzghr/Qm/kAAAJdCQAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJdCQAAAllByM7fviDpvZ9ZXDbGZOZ3bepCCztX2LzfzDtPW+0IJ0QM8PamrpOBvZeO5RNnPcoC3JhdS91WcGB3XOmuQz938rWxh+npbYzKOf9O+d5iWn+WHW/31ikAne32DAVtODMpKO+eoDNrNUE2zmhA/5xerLFyc3kh4OMp0DfeahLYNsZtgtvwxOk3Resnc7eH+n+oH0n88fEZwlDWz8e7c4qDOz22f6fuqJoJL09L4bbKZtT4lq7YlfSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASSg7G6meNjcwf7cucGB63f7Lpddo2X2jBnOA0P8wqSa/Z4QdaX97cx2ZWB6/Tz5ILSeoM5vgeWhcMKT7+kC80Yt/gRpJ+szgIbQsyhweZtwcZSRoaZDYFmeSdeTzISNIAH3njJJ9Z/4KNbDz8OF9H0nEj/UB2dzBrnrwrY9YEIUl7Hb3bZl45NBmyDQZsp84M6ki98wfYzHPB9tn7grOmha9Ts9m3jXZiVmtP/EICAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlFByMHZN4wdjtwV1On6VndcsC16CTyebXsfYRL+2b1BHenKRH9ZdOMXX2S84q2NBEJI0/qJlNrNy5Ad8ofXzg9OeCTKS9AkfmfV6GxnyhfttZrIWJRfSGK22mdXBZ2WRJttM7+dOCG4kafbvgtDXgsyBPjJiWlBHOn3dd2xmxXXn2EzPVH/W88mFJE1Z6DOHTPZDqDubXcFpq4OMpC/77bPtOf47s+cIf9SA4DqStHewQfskBVPNf4ZfSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJKPqmhO3hSQ2ewLbw565XswKP9EwiSCfXX7Hinzby8wq8dl6RrgqcwHB/UGfOIzzS/Dz8Cw5K128lrOcNH7gieCCDphtPOt5kP9yzxhS73kbu3BheStDnIHB1kThkYhK4IMpK+2THBZi6882ZfaGzyBI25QUaS/FMYtNGvg29f578vVh+b3Ed6MMhcHDzNYa/xySr0nwSnSdFTSzb717Jd7n9/dGdb1dU5Igite/WthV9IAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJKDsZqkB90m7TlGzZzY3NmeOAtPrI0WCPc39974ejkPtIxQeYkP3un/ZYGb+/k4N8vSXrKR87z66tnLbnMZr5w5VXJhbSm02e+H9Q5PMicHGQkaWAwNLg12O6c3PvxICNldx/d7TOfu+xSm5k94crgNElLk1X2B/vIonNt5Pnz/N+mJN0TzK0/GtSZssZnmu3hV+951wQh/xpMbG+zmcWDPhqcJc0JhsRn/idaC7+QAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJdCQAAAl0JAAACWUHIzdqQNs5pBxz/pCd3RlB47xuY0/OspmdjZbbGZ7ch9JHcmm158kQ6/Beku9KchIuvUMG3nirL42c2Cfp23mS89FN9KYJHOJzzw870ibuSJZKytpxXPjbWb8/its5vJgHezg6Y9Fd1p9dZAJ6nxmf595ZvdBQSXpsOW7fOjs24NKv/aRRcEKZkntO/0AbU+wfbZ/cFa/dlCQkoa9+5c+tLrLZ8b6zJOr/HevJN3Z+Kn8DgZjAQB/q2hIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASSg7G9n1xh808vW+yT/Nt0XkTW18r2aTYHWxR7FyQ3EhqRgdvy7Bk0+t+PvILP/AqSe1Lfmhw1Um+Tm9w1gw/WypJWnn7e2xm/M9+4AudHhy2a1UQkrJ/4RAf6TvOZ1YGR0la8fb32szpZ/zQZub6ed7kXyZJGnePzzR7B38Hb02GZ58PMpI2+s2r7Rr/d9A91R/VOTC5ULodO9kJfK9NHPRCthd51wWH+dAyBmMBAH+jaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEko+qaEJhoC1o9tnFndG57V7BZPXF/g6nSN8Zvy6ZcGNpJXNcUHqxz5y6zQbafv7f78k9QRPYXhtUOfc23xm6D+uDSpJvQec4EO75wSVkpt/MMhIGhashN8YrN3Wt4PMH4KMpD4zbWTIs/fbzKZ/HWUzt5wZ3Si6eUfyNIftwVfY2fOD0yTpXTZxevuQzawYeY7NdK+PLqTOm3ymeTl4DSYF35mHZt+Za58YajMnaFNUa0/8QgIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUUHMwtlkcpPz03bz2n6Lz3t1cbzOPBnUGt4NsZtjIXwaVJK3v8pnzfOaJJX1tZl3ztD9L0jNB5tx1PtO8GHzkTloTnCZFw8H7BsN+D/sxzRuOnBjcR/rwvy2xmW++eYLNXPjYjf6wwclAr6QXgqHIYChU94y2kXafbND6lpE+c2BQZ2R7kM0cNmFXUEnS0i6fGeEzG9cdZTMPN1v8WZKOCTI/aj9iM9ObfwkqBVPrkvq177GZJ3VEVGtP/EICAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlFB0MLbLh671mXbvbECva2qQmeQzzeznfaj/t3xGUrKddFb7RZu5tM9VNrPguehCmhHMzDV/nwy9Lg9OC15LSZrdYSPf7TzNZt53/l02s8rPu0qS/N5VKdhzq3F+dlbfu/nUoJL0/u47fWhWT1BpPx+556ygjtT+u//7nBtsn526v89ctfvS4EbS7OazQSrY5Lv9QzbSzgpeS0ldi4PMAp9pXgr+Nj/e5TOSdLLPtXdnpfbELyQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACUUHY1fbTEfrd7h+fV+/RVGSbn/RZ4a2Q2xmVL9NvtBv5gQ3knTHTBtp1/rBwq5gWWrX+ORC0tDb19pMb/O7oNI2H1npB14lqR3qX4O1fpGvVgRnZftipYGTfWbrIp8J9sUqfOs0KlhO2mwKvgpOT4ZnBwQZaUj7epvZdMYom+kK3ryuZGGupGZU8BqMDf6G3+j/ftfuHBrcSNrU9NrMGfv4Oh97wW/G7mmS/bSSdLBNtO3xYa0/4RcSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoAQaEgCghJpPahjoM09uPcBmljS7o/MuHuYzx25YbzOPNuuC084JMtINrR/3P7bxO7VfCs7a3b4nSEnjD/hBUCwYiZ/tHx/Rnp+tn18ePIUhWYbeMc9nPnnJlUElad5yvy57+ll+tfxXr77MZnqmR1dKFo/rrORpDjcHXxezwsci9PGfgxXPvteXaX5oM3tHF5Ieaf3e+Aub4DEbWmYTx7QjgzrSI8NH2Mw1G32dCW0fmzlk4LPJlaRtc22kbWdktfbALyQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTUHY6/zmfbXwfruz2fndd3sM819wcs0r8tnZgUZSe1/Cf59F/g6XZf4TPOB8CMwOljdvK9f3fzd50+zmf7NXcmN5Jc7Sx3BvHLzXPAajH4wOE2SbgsyZ/rIGr8Cut0/GyDuCWYwhwR1tren2sz797szqCTpheDztMZ/ntrvBH8rVycXkrpu8pnmfwafldldPjM9yEhqTwz+fef7Ol3/7DPNm8LvgqldNtK2PvPn+IUEACiBhgQAKIGGBAAogYYEACiBhgQAKIGGBAAogYYEACiBhgQAKCFdpPj/VP+LHvWhYIrv8PC8ByYc40PB4Fly4pAv3J8UkoJNqMm/7+F5R/rQG4JCkqTXBgf+wUbed74feu0KbiNJXcGm12zodXFwmt9SLCkaGtSC5T4z2n9WmjXZIGM7LxiuDLbPdgXvXfIZkCQNCD5PpwfH/dZ/xg+/+rHgQpIu95EhW/z70js7+OsMPruS9EDrv58OPz/4zvQLptW/N6gjafvUg6Pcq8UvJABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJJTfGTta1NjOj+YTNPBWed2d7qc3Mbt4XVPKDk19p/b9NkoY3823m6KBOZ3uDzSxpDgsqSdIJNnFD+0mbOazxE3qDo/tI17RX2sy85r8FlfxA4MEvnRTUkf7jpv4283cXbLeZp/a+JzgtGOqWNL397zZzcXOZzTwcnPVEOyFISRc2Xw1Sfgh1QvuEzXQ3FwZnSZuDzIZ2ms18qvl4UOnZICPNar9nM6c1V9lMMso6t/1akJJuaI6ymbYdF9XaE7+QAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJZR8UsPVushmjm2ut5kT98rOO/+l79rMXU0yVd1hE2vboUEd6d6m12amjfB1DljzpM3s7nNTciVp2Awbaf+HX5XdfYQ/qnNyciGpGRt8fM/u8plg7Xj7Nv9vk6SuDwWZb/lMc2/wb1vQ5TOSdKvPtXcE790if1Tnr4L7SGr+a/Dv2zjXRvrsvsBmnh19SHIlzV/vM29rh9jMqGZTcFpPkJFObf0TYG7e+/02c9/L/qxH2o8kV9L05hKbacPvuj3xCwkAUAINCQBQAg0JAFACDQkAUAINCQBQAg0JAFACDQkAUAINCQBQwt5/7Qv8JRs03GZOCer0BoNgknSi7rWZu5SsMPfDcKs1JqgjDZcfjN0aDPGN33+FzSyRH/STJG38tY18881+ffUJ8ivMtwYDmJI0faFf3TxPZ/pCC5bbyN9d69eOS9J/6P/MCnN9yN9Jyb9N0vSz/Ou09Wxfxy+xzz4DkqLPk4LPZvIZT/5WJOnoIJP9DSeDsdn6+RPlV5gn33VHBmctD757/+ixIMNgLADgbxQNCQBQAg0JAFACDQkAUAINCQBQAg0JAFACDQkAUAINCQBQQsmNsW/Wozbzb0OOtZmFD2fnndj6AbURzSNBpYU2MaQdFdSRNg3yuYVbfZ13tH4cbsgbtgU3krTLb+/Utk/YSHvp62ymy8/O/jE3z2eatwYf8dGLg9P85k5J0tSzfCYYxJWCLcVrJgV1pPYXfhts13RfpyuYeW2u+n1wI0kDvuYzff2W4t7fDrCZnzbJIKc0ZaDPDN2y1mZ6G5+RpgQZaX3rv+vua/x35pTB/qw39ybfc9L2ZpXNtO3FUa098QsJAFACDQkAUAINCQBQAg0JAFACDQkAUAINCQBQAg0JAFACDQkAUELJjbHbrws2KQYDeo9/PjvvLUv8UJmCoUHNe9xGej+XDcPpCh95/AKfGTw9GAhc6SOSpNF/CA58rY187/lTbWb8kruSG6kneF/adX4otFmTDM8+GNxI0oKuIBRsel1zvI20+/t/m5S9TuODOt+72b932s9/Bv4o+DwFn83kM740uI2k6O+u93PJ3tzbfST5TlH2/fSvSaHgOzP67pUkfTvMvTr8QgIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRQcoV5E6wRfnKrXye9pNkdnXfxMJ85dsN6m3m0WRecdk6QkW5oJ9vMsY3f8/1ScNbu9j1BShp/wA+CYt0+M7vTRtrzsycQLB/kM88HdTqCVeifvOTKoJI0b/mlNjP9rKts5qtXX2YzyRMYJGm/IHPWFp9pbg6+LmYFnwFJ6uM/Byuefa8v0/zQZtJH0jzS+scZXNgsCiots4lj2pFBHemR4SNs5pqNvs6Eto/NHDLw2eRK0ra5NtK2fv38n+MXEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKCEmoOxzWqb6Wj9Wt+v7/uR6LzbX/SZoe0QmxnVb5Mv9Js5wY0k3THTRtq1fni0y88eqivZXS1p6O1rbaa3+V1QaZuPrOwI6kjtUP8arA2GZ1cEZ00MMpI00M80a2swW3ljcFb41mlUMvS6KfgqOL0nOG1AkJGGtK+3mU1njLKZruDN6wpndZtRwWswNvgbfqP/+127c2hwI2lT02szZ+zj63zshettpqdJV5gfbBNte3xY60/4hQQAKIGGBAAogYYEACiBhgQAKIGGBAAogYYEACiBhgQAKIGGBAAooehgbJcPXesz7d7Z1tGuqUFmks80s4PdpP2/5TOSpA/axKz2izZzaR+/mXTBc9GFNOM2n2n+Pvg4nbQ8OC3Z8yppth+g/W7naTbzvvPvsplVfkGvJOn+IHNCkBnnl5fqezefGlSS3t99pw/NSoZeg92z95wV1JHaf/d/n3PP9HWm7u8zV+32W3wlaXbz2SD1bR/Z/iEbaWcle3ylrsVBZoHPNC8Ff5sf7/IZSTrZ59q7s1J74hcSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoISig7GLg5SfmJvX/lN03rsbv0nR76eVBrd+Nemwkb8MKkla3+Uz5/nME0v62sy65ml/lqRngsy563ymeTEZnl0TnCZJP/aRfYO1uQ//wUZuODLbGfvhf/MTtN98s596vfCxYGfs4NcmV5JeSFamvstH7hltI+0+2UD6LSN95sCgzsj2IJs5bMKuoJKkpV0+M8JnNq47ymYeboI1vpKSHa4/av127OnNvwSVgul3Sf3a99jMkzoiqrUnfiEBAEqgIQEASqAhAQBKoCEBAEqgIQEASqAhAQBKoCEBAEqgIQEASqAhAQBKqPmkhsOC0I5g8nxxMKEvqd3LT5Z3X+DrdI7wmfHrlgU3klY2xwWp4CkFt06zkbZ/Nlnfc5LPJM8NODcYBh/6j2uDSlLvAcEy8N1zgkrJzf1aeUnSsDf5zMZfB4WCVdnyT5iQJPWZaSNDnvXL1zf96yibuSVYOy5lN++4x2ea7cFX2Nnzg9Ok5GkVp7cP2cyKkefYTPf66ELqvMlnmpeD12BS8J15aPadufaJoTZzgjZFtfbELyQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACSUHY/u+uMNmnt73+0Glt0XnTWx9rcWDPmoz3Vv9WZ0LkhtJzejgbRl2S1BpPx/5xRlBHal9yQ/QrgqGZ3uDs2aMD0KSVt7uVymP/9kPfKHTg8N2rQpCUvYvHOIjfcf5zMrgKEkr3v5emzn9jB/azNwV/qzgXyZJGpcMve4d/B289fbgtOeDjKSN59pIuyYYpJ/qj+ocmFxImrTlGzZzY3NyUOlemzjohaSOtOuC4OkFy159a+EXEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKAEGhIAoAQaEgCgBBoSAKCEkoOxO3WAzRwy7llf6I6u7MAxPrfxR0fZzM5mi81sT+4jqeMRn2l+Erx1kxcGpwUbTiXpVj9A+8RZfW3mwD5P28yXnotupDFJ5hKfeXjekTZzhS4PTpNWPOenesfv7ydML9cVNjN4+mPRnVZfHWSCOp/Z32ee2X1QUEk6bPkuHzo7GXoNtu8umhLUkdp3+qHXnmN9nf7BWf3aQUFKGvbuX/rQ6i6fGeszT67y372SdGez22Y6/hOthV9IAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJKDsZqkB9Oy7YonhkeGGxeXXqxjbT9/b0Xjk7uIx0TZE7ys2nab2kyPJtsnpWkp3zkvGk2MmvJZTbzhSuvSi6kNZ0+k+wWPjzIZLs0pYEjfGbrep9J7v14kJGyu4/u9pnPXXapzcyecGVwmqSl84PQwT6yyG95ff48/7cpSff08ZlHgzpT1vhMsz386j3vmiDkX4OJ7W02k2zGlqQ5wXbsmQzGAgD+VtGQAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJdCQAAAl0JAAACWUfFJDd+Onqjvn+DrNWa9kBx69LAgdaBOv2fFOm3l5RTAKLumaYOPy8UGdMckq9N+HH4Fhm4JQ8lrO8JE7/OstSTecdr7NfLhniS8UbCe/O5hOl6TNQeboIHPKwCDkt5xLkr7ZMcFmLrzzZl9o7DPBaXODjCSd4yMbh9pI+zr/fbE6WDsuSQ8GmYsX+sxe4/1jVF459CfBaZIUvOab/WvZLve/P7pnJveROoOnkWgdT2oAAPyNoiEBAEqgIQEASqAhAQBKoCEBAEqgIQEASqAhAQBKoCEBAEooORi7JhiM3RbU6fhVdl6zLHgJPp2sER5jE/3avkEd6clFR9rMwmB4dr/grI4FQUjS+Iv80OvKkR/whdYnq6uTAUxJ+oSPzHq9jQz5wv02M1mLkgtpjFbbzOrgs7JIk22m93MnBDeSNPt3QehrQSYYWB7h19hL0unrvmMzK67zA589U/1ZzycXkjQlGHo9ZPJjNrOz2RWctjrISPryxTbSnuO/M3uO8EcNCK4jSXu3b7GZk7Q+rPYn/EICAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlFByMFY/80Ne80f7MieGx+3f+iHUIdO2+UILgjW2GhdkpNfsOMpmXt7st8+uDl6nnyUXUrYl8qF1g2xm2OMP+UIj9g1uJOk3i4PQtiBzeJB5e5CRJL/lVEq27ybvzONBRopGHt84yWfWv2AjGw8/zteRdNzILTbTHcxWJu/KmDVBSNJeRyebXn8ZVFrlI1Oz9ay98wfYzHONH9a9LzhrWvg6NZt922gnZrX2xC8kAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAl7/7Uv8Jf0fesTNrOr+zCbmdOZnTcpGCr7eeunQv8h2ZC4INk8K71yqB94bJb64bSNrR+w7RzkBxSlbEhxcONrtZP8Htv7dw5JrqTztMRmHv3kJF9oXnLavUlIUk+QOSbIBCuBpwdlJB3z1QdsZmkw0HvCh3ptZvni5EbSd4NM50CfeWiLH8ZubkmGWSWNDgZa9aiPBEOvP58fTJpLekPw/fS9oM7Mbp9JvnslSaPv9pmJp2S19sAvJABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAklV5g3za9tZnT7oM38dMqp0Xndi3xmclDnt8kq9K9sCypJ+vQtQci/ThpzsY1M/NH84Cxpcc9HbebBC3yd24Kzjg4yktQxLAhd5iMPTPBPTrhVZweHSffpbTZzYvDUh7N1q828ZUnw1ABJutJHejb6zObgqDODjCQdf5PPTOr4hs3c+O5pvtDq7Akp0pt85Mvn2kjvpwbYTPIEBkkKvp7UGXxBvWPhXTazpjk+OE2Svm0TbTsjrPUn/EICAJRAQwIAlEBDAgCUQEMCAJRAQwIAlEBDAgCUQEMCAJRAQwIAlFB0MLYrSH3eJk5v/WChJK2Yco7NzAmm004Ozhr1qyAkqfn9Kz509JNBpWTo9YNBRtJiPzw6b+JFNnPJddfbzDOY2CzZAAAD3ElEQVSXRDfS4hd95qmgzuFBJluqLg3Zy2d6Xw4ywVl+0f0fHRxkJu3jMwde7TNXX/SR4DRp+o3X+dCkZPDXD2lKwfCsJG0+xEba1/n/jl97hD/q+8l9JM0Mhl7HL1xmMyubZLD7n4OMpLFdNtIm2+D/DL+QAAAl0JAAACXQkAAAJdCQAAAl0JAAACXQkAAAJdCQAAAl0JAAACXUHIwdHoQ2dgUhPzwrSaNbP6L20yv99tn5nf6sA5ILSeqY4zOTPhNs0xwXDATe0RPcSJKSDZef8JFrD7SRjo8tDM6S5spvpex3xW5faImPbHo4uJCyV8nvFpaGDg5CE4KMpJ2X97GZGZprMz1fn+IP+/gzyZUkfS3IBK/U2A4bmbgq3Ir8Jb8VuWemr/NscNa07iAk6R2XJZtek7H8YOh1WFdQR7p+g18NPUXBSuA/wy8kAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAklB2O/FAxXzhweDNVFw7NSsk3yoBf8y7TrF4fZzJrR0YX04yDTOdBndm7xA5HHvJhs5ZSeHnCoD+1YHFTaFmTGBBlJA4LcZ32k/0X+NThNd/pCkoZrg81skJ/+vlOn2cz26/wWX0nSF4PMttVBKMkMCDKSDp1kIwdt22Ezj+7jX4N+g4LhaEndW33mXUGd0Wt8pu9bnwgqSU/v2wSpYPA3GHqdsyEYbJf0meHX+tCGV99a+IUEACiBhgQAKIGGBAAogYYEACiBhgQAKIGGBAAogYYEACiBhgQAKIGGBAAooeSTGjTcTyYv3ODXFn9kXLhC9w6/ulnBZL0Wn2IjP584IjhLOmnkAzYzZ72vc3hwVscHgpCk+5cNsZmxusNmdp5yhD/s+w8mV5Lk189LyUrtg4NM+FSEaEF5sug8eYLGU0FGkvzaeClYg33y8TbS7+5fBWdJd2iszZxwTq/N9HzHn/V4ciFJM4M/z3vWvcVm/uHG4I9z0t3BjSQFT/7Q2Bk2cv2qYO348J7kQuraGGT+E62FX0gAgBJoSACAEmhIAIASaEgAgBJoSACAEmhIAIASaEgAgBJoSACAEmoOxgIA/r/DLyQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJNCQAQAk0JABACTQkAEAJ/wvHQ34Y31sjKgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31dbdf550>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imageplot(displ(D))\n",
    "set_cmap(\"jet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computation of Geodesic Paths\n",
    "-----------------------------\n",
    "We use a more complicated, non-constant metric, with a bump in the\n",
    "middle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "x = collect(linspace(-1, 1, n))\n",
    "(Y, X) = meshgrid(x, x)\n",
    "sigma = .2\n",
    "W = 1 + 8 * exp(-(X.^2 + Y.^2) / (2*sigma^2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAFJFJREFUeJzt3duS28ixBVBqrr79/5fa4xmPbZ2n46jOkbI6CYDcTa31pAqAJNjdoQzkLlR9+vz58+cbADzZd8++AAC43RQkAEIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiPDDsy/gSz59+vTsSwDggHvW7XaHBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQQUECIIKCBEAEBQmACAoSABEUJAAiKEgARPjh2RcAyT59+vTsSxj5/Pnzsy8B7uYOCYAIChIAERQkACLIkHh5Hy0HOuLId5U/8WzukACIoCABEEHLjg/pGW24xNbfmW22yffT3uMK7pAAiKAgARBBQQIgggyJSFfmNYn5072ZzLOmeXefK1/iXu6QAIigIAEQQUECIIIMiac5M8u5971Ssqr13CszmPW9r3ru6Kq8jNfnDgmACAoSABG07HiYI+2xe9tfZ5575muPvO+jWl6T9t6Rlp4WHv/PHRIAERQkACIoSABEkCFxqoSc6MwM6RkZ0y5TSVi2p37OkWtK+D5kcIcEQAQFCYAIChIAEWRIHHJVxnIk27nqtVfmTZ0uR9llOd3xq/KZ6TV1PLP0bXGHBEAEBQmACFp2jD1iGZ9pe2zSajvyOe89dqVJC687Xs+9qr131RRx7bvX4w4JgAgKEgARFCQAIsiQ2HrUkj5dtnPma+89t/vMe45/zSQH2h3rxpNzqzOnY3fbXDzqGsjgDgmACAoSABEUJAAiyJD4g6u2b7gqB6rj77777t3nnpk3VY/IkKY50OS1//3vf9997lmOLDskU/r43CEBEEFBAiCCggRABBkS46zjrHXjznw+aM2NumP1+FV505fG73Xls0RdLrQeu93eXn89dtXzTdWR9fWse/fxuEMCIIKCBEAELbtv1FXLAZ3V8tq10urxdTx5bfc+02s6sszQvcsB7Vpp9fg6rsfq91uP1+vvPrceO5Nlhl6bOyQAIihIAERQkACIIEP6RjxqOaBJhjTJaybj77///t3nTt53cv1fGr/XmVO3JxlSHf/nP//53793363LjXbXeC/LDL0ed0gARFCQAIigIAEQQYbEeBuFs54tOpLtdDlRPdadW9/3hx9+ePe5CRnSLgf697///e7XrpnR7dYvHdSdO8mX3nP8vY4sM0QGd0gARFCQAIigZfei7m0Xfem1R1p2906x3rXduvFu2vfalqstuu59p1PRzzKZul1bafX7rS282s7rvk/XovvS8XuduezQpE1nGngGd0gARFCQAIigIAEQQYb0Qs5cHui9rz2yLUQ3lbtmObusp8uFuvGRDKm+tqrv9V67PGbNfnYZUh13GVI3PmtK+9RkSaLONCMyZfw53CEBEEFBAiCCggRABBnSN+Ksbcjr+MhzSF0ONM2Qfvrpp6+e++OPP351PMmbJksS3W79sknVZEuJNRfqMqLd+Pfff39zbPJc1ZHt2ju7LSXu/ZwjW1XwOO6QAIigIAEQQcvuA7tq6vZVK3hPlv+pbbZJG25t3+3O7Y7VcbcS+JeOdz/jam0Z1fZRt2J3PVbbcN149/v417/+9b9/1+uv79t9v0k7bHfu+t3PbA12tPcexx0SABEUJAAiKEgARJAhvahnbT8xyZAmS/rssp6ff/75Xcfq+EjetMuU1u87yZDqVO5uF9hJZnS7vc2FfvvttzfHrtrxdnL+brfce3O5HTlRBndIAERQkACIoCABEEGG9IFclQvt3vfIc0jdlhKT55Dq+E9/+tOb8ZoFdZlRfW3NkOq56/FdhtQ9h3Rk6aDJ8j9rRnS7/TEnWq9x8gzZ5G+vZl67XGiSIa3jXe5zZGmhjrzpOu6QAIigIAEQQUECIIIM6YXcmzFN17KbbEs+eQ6py2u6bKcer/lSN96d2z3ftHtW6qznkGpO1GVIXWa0u6azMpddptJtr1H/fiZ501nr6dX3khE9jjskACIoSABE0LJ7UdM2XHdsshNqN5V7snTQbkmfbir3X/7ylzfH/vznP391PGnv1Wuq4649Npn2vdsFdp3aXad5T6am73a8fe/13m5v21r12GS8m/Z9b3tv2s7WpnsOd0gARFCQAIigIAEQQYYU7shyQfd+xmTa9y5v6nKUbtr0ZAuJ2+1t1lMzo5opreMuX6rv2y1XdLv1ec1k2ne3VNDt9jY3+vXXX98c22VIk6noXbbTZUiT7TPqNdVzu7+v3bTvR03dNkX8PO6QAIigIAEQQUECIIIM6QM7azuKaYY02X5iHe+2a5hkSN3zQrsMaT3+t7/97d3n7jKkLr+ZPIdUM5Zui4ndNdRnjTo1/1izoF2GtJ5bM6Q6rr/37tmiOl7PnfydHtkywnYTj+MOCYAIChIAEbTsXsikDdcdu2q1713Lbh3X1tNu2Z5uenZt4a1tur/+9a9vjtXx2sLbtexqW/Gsll1dwXtt2e1adJO2brdk0W45o7WtWI/tfu/r+ZO/p90U8cnffEeL7nHcIQEQQUECIIKCBEAEGdI36qxMadfz77Y76LajqNnIZBr4LkNac6FdhrSOu2WFbrd+KaEjGVLdYmJdLqjmMUd2pu2WLNpNRV9/P/XY5Pder6nb+uRIZiQXyuQOCYAIChIAERQkACLIkMJctd3Emc8hdX38LlM6M0PqngHabUu+ZkF1qaAuQ6rn1vd91HNIk/ftlgOqWU833m2Vvo67LS++NO6ea5v87SVkSrKqY9whARBBQQIggpbdB3JmO2/yXpNp30eWGVpbOd1K4Ltxt6xQHU+miHcrgd9uf2zZdd+n6pbpqa/t2nS73Vm7NlxtDXarr3crm3ctud141wLuWnadM1tp2nLXcYcEQAQFCYAIChIAEWRIL+qsjGg3nkz77pYVut3eZg+7c7vxbquKdTzJm2pmtMuQ1pxlMu27Tr8+khnVXKj77t3SR7uffze9f5eBPWPadyUXyuAOCYAIChIAERQkACLIkL5RR55purePv8tRuu3Od5lSt3RNd27Nfbolieq5u/H6OTVXqdZnj87cqqJeU7d1xeRn2v0+Jr/n2+26XOiM1/FY7pAAiKAgARBBQQIgggzphUyW25+8z5Fl/btza67SZQndufX47jmYyfM13Rptk/X1dhnS5FmjNSeaXH8dT7aJmPw+dudWk/Xp7t3CvHuf2+2Pzx2txz2T9DjukACIoCABEEHL7gNLmMp65jV02xBMtiU4c0uMyVT0ruW1mwq9qi26yRTrq7Z6mCzpM/muO4l/41p413GHBEAEBQmACAoSABEUJAAiKEgARFCQAIigIAEQwXNIH1i33MmzruGI9fmb+r51XJ/VWY93x+rxem43rls91HPXLSRut9nzOOtr6/t025RPrr+Oz/yZducekfDMT8I1fCvcIQEQQUECIIKW3QtZWwtHljvZtcu695qcW1tT67n1dd259Xg9t7ba1nF3rI7r+/7+++9vxpMVu6v1vev71vF67uT667geq99vHU9+H7tzq+61k3OP/I1Pj3MNd0gARFCQAIigIAEQQYb0jTrSI5/08SfTgSdTrOvxNWfZ5Sjrub/99tubY+turPX4bvuJar3mybk1M6rXuI4n11/fe/Jzmvw+Jr/n2+3+v6cJmdDH4A4JgAgKEgARFCQAIsiQXtRkWaFp377r+XdL1exyoPV5ld253bjmKN24Hvv111/fjLvtwqt6zet71/yp6p53qjnQP//5z69ebx1PvnsdT57XWr/77jmwyd/I5G/vyueSeAx3SABEUJAAiKBl94Gcubr3WUsJnbmydrckzm45nXW8a8P99NNP//v3jz/++ObYZGr3rq14b7uvvk9t2a3f55dffnlzbG3n1XPreNfeW3+mu59/13K8agXyZ7XktPeu4w4JgAgKEgARFCQAIsiQwly1C2z3vo+a9t1tb3C7vc0hdplFt0ROzYW68ffff//mWPfz3mVG9ZrOypC6TKxmRv/4xz/ejLuMaTdlfP2ZdvlSHU+2tbjdnjPt+6ocSL50jDskACIoSABEUJAAiCBD+kYd2Ya8ew6py1nqMz5HMqT6bM763vVz6njNjXYZ3fpdd88H/fzzz1/93DOfQ1rznJr71Myoy5Qmzyx121jU8TRD6rauOLJVxXuPkcMdEgARFCQAIihIAESQIb2QSZ/8yHNIk/XpumykG++eO+pyoe7YTv2ua95Rr2ldE+92O2/rit06fl2GVHOhmin9/e9//+q5XYa026pisjV6tx3FkXXvjmRKZ22VzjHukACIoCABEEHL7gM7sszQkSVYJtO+Jy27te1T22x12nFtga3ffbIcUFWnJHfbWtSWXTfd/Mj08m4pod3yP7Utt44nW1d0ywrdbv1WFWe27O79uz3SgtO+exx3SABEUJAAiKAgARBBhhRu7V+ftRVF9xnvGa99/ZrldNuS13NrlrAe3+VAXYa0m2K9fp9uGZvb7W0eUnOTXYa0fofJtO/dNa0Z0m679m7c5Uv13Prdj2RI3d/IbvuJs6Z9n0nGdB53SABEUJAAiKAgARBBhvSi7l1G6EuvrX399fwuD6jnTjKkmo3sMqROd42TZ6PqVug1M6rH793monsWql7jbluOLkOa5E277SfW39fuZ9ptR9Ftb3679c/AWSro43OHBEAEBQmACFp2L+TeKeK7dkXXwqufU8fddObuGnfXf+T7TaZ9r62oXcuujidT0bvpzEdWRZ/sNttN7d69b7fE0q5l94zVvne08J7DHRIAERQkACIoSABEkCF9IGdtN7GzmwbeTb2tr13zgl3etI5rNrJz1pYY9XPX3GiXIdWc6Kxp3/X6uwxpN17znclyQLtp3+u4m9a9Gx/ZfqKabFUxIV+6jjskACIoSABEUJAAiCBDelFX5k1db77LlGp20Dmy1cYuL1ivo8uM6niXGXXPIU0ypN1zSJMMrBtP8qbduet1TJcO6p5DmmxLLhf6+NwhARBBQQIggpbdB/aoaeDda2uLpZpcU3furjXYXVM3zXjX8lrbcHUX2zo+sotttyRO9312Sx9Nlh3q2nB1OaDuc6fTvidLB01adpMp4h3tvMdxhwRABAUJgAgKEgARZEjfiC5vmvbIJztz7jKmez5jd3yXQ6wZRp2q3eVEu2neVX2v99pNj1/zmkleVl87yZt2U7cnGdJZW0oc2X5CLpTJHRIAERQkACIoSABEkCG9kMkW5vdud15fW92bGU1NtsSo+cea/RzJkHbjs3SZy5kZ0uR5p248yYzqeLJ00MSZOSnXcYcEQAQFCYAIChIAEWRIL+pZ69xVZ2VKuyyhyyFqLtRt39BtMbHLjHZbtL/Xke+6y8/uzZ92506u6VHPIU3IjDK4QwIggoIEQAQtO7btvbPaGbv23ZFlYdb22W75oq4NV1tR6/H6c0lo2e2+62SKddeWO3Pq9qO2lNCG+3jcIQEQQUECIIKCBEAEGdI3YjIN/FmZ0uSaJtlCl/1MpnIfyYwmSzntjnXLJE0ypclrJ1uLT5f/ecaWEvKlTO6QAIigIAEQQUECIIIM6Rt171YV9fyrevG7fGmSq9TsZz0+yYV2mdGjnkM6krlMnmG6NxeaXtO9ywHJjF6POyQAIihIAETQsmO8MnjX7jvSRjnSHlvfq7bhus/ZtQYn11Sd1bLrjl/ZHnvUKtxXLQekTffxuEMCIIKCBEAEBQmACDIk/iBxmaHJ5+6u6UhW9d5j7zn+NUdylKumjJ+Za3WvnRw7ci6Z3CEBEEFBAiCCggRABBkSW2dlSkd6/GfmWl2GVB3JkK5yJHM5a5mes87dvfasc/kY3CEBEEFBAiCClh1jR1YKX505RXyynNG9bbirpnnvnNWie9Rrz1z+58zXks8dEgARFCQAIihIAET49DmwKfus6bQcd+R3N3ntWdOxr8yFrlo66Mhrz5pW/ajp2YH/PfFO9/zu3CEBEEFBAiCCggRABM8hcarJc0dHXvuobOGq73PEozIZORGP5g4JgAgKEgARFCQAIsiQeJjJmnK711ZnbXMxMbmmM9/3iHvf+8xrkhnxNe6QAIigIAEQwdJBRLryb+AZf1/PmiL+jM8M/C+FJ7B0EAAfloIEQAQFCYAIpn0T6aop1e9576TPfNb7pn4ur80dEgARFCQAIihIAESQIfEhTTKMs547evXc5NW/H/ncIQEQQUECIIKWHS/vSCvqoy1jpe3GR+YOCYAIChIAERQkACLIkKAhk4HHcYcEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiPDDsy/gSz5//vzsSwDgwdwhARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARBBQQIggoIEQAQFCYAIChIAERQkACIoSABEUJAAiKAgARDh/wCAHYukN2JikAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31dd46d90>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imageplot(W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = Int.([round(.1*n); round(.1*n)]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 3:__ \n",
    "\n",
    "Compute the distance map to these starting point using the FM algorithm. _Important:_ use symetric boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzsnXu4VXW1/scOEARk0wbcyEU2yhZUbgq6TTB24hWvpXnFxKJHLT1ZmdrdU1lZeelkpSc6YuGF8lJKioq4SfCIQYIgYogs4iKbWyAXNyLu3x/nedYc72ey52Kh53dmPeP96zsYa68113fNtSbzfcd4R0Vzc3OzBQKBQCDwf4wP/V8fQCAQCAQCZnFBCgQCgUBOEBekQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AJxQQoEAoFALhAXpEAgEAjkAnFBCgQCgUAuEBekQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AJxQQoEAoFALtD6//oAdoev2bcknmDji+t1Pz1QHzwJfzxnuwteRHIJ4rWId+7hEZqZtXHr/ZGrRXx0shzeXlNjNez2hb9LPN4mFNeft59LruddG/WP70yWM+dpahaO6G237ofcxVUaV1yh8cab2hXX38VndfsrX9UHX+vWU1filR5G7N8PjwobdY2GfW5bXFz/wG6Q3IXT/6gPdofcgI3hPvkvyAnIDTsD//B9DW8deGVx/ZXGH0vuvSs66IP/4Pcma1/MzAYly3bnaOonGh7/+SkS/9iuK66P/Omr+mAc/0T39SjgCGoQX+rOmYrrNTf/Ov0+XGe6F0/9/KwkuNYUTQ/hHxa4NU5UO1/DsdXFZeWENZK6te2XJP70jPv1b917eHq2pniOdHJrnhK14zVuuj1Zf6XDzZK7Y8Z1+mB872zxTBfMQJI/5acmy2MGa2qCht84/GvF9XdX/ECTX9Bw7iMaT3Prd3EEX98L3+64QwoEAoFALhAXpEAgEAjkArmk7DxFZ2a27qYDfVJRWIZ/8DfUy5EjJdcGsacAOiH3FuItbr0KOVKBjo6ZM0JT6/tKuG6rUpITvp7sRXvbLrl/G/8fEndqTN7fSBxC42qNPfHBo58Fhmgk2IGqx5uK6/oxDZJ79vCPSTx/5DFJMK2XPtG7PfHK/oW5h4s0nHOYhMvnDyiu5w05QnIXDFXKruLQZN0H/Mtf8ar+U+c+DeOpBwbskIGvFdcHV78uuSX9huCP/d7w3CNl5869JpzTr+s5XTA9v1ZY7+L6yANxwGCeqzIoO34b1rpDrMa51sPelPgAwwP8W++OJy5wL/z7I0m0QcM1CWW3eWW1pBoP1tgQ+r3IOgIz/VXhvvCQ2rk97doXya6WHcuR8Kebe+GOZBNSyl7ahsOTF2oCC9quS8tHYKZ78ba9f8QdUiAQCARygbggBQKBQCAXyCVll6qk8zRdAdRNqubFEyu8uSZFxGo4x+Wk7uEbEXu6o1T1nqdrcGNbAIU3QamodR2TvfjNFy6RXI9WSn18+jxXKYRDqkM1oj9CvrPXER+Jir32jsI7bsyfJXesPS/x/HpH2Q3FE8/h/vtXJvkB+nWx7pPNSZYLhgyS1KKqgyQ+vPaN4roGFEUV2DF/FCTOdq7QuA0ovL6O6OpnSyWXouw6uvVWVo8VEPujwqdXUFp09eYDJF5RmVB21gNPC8rO0zP8JpGe8UdRjdO/29qtElfvjwf4QyZNVeBe+J8sHsUWDde79ZoKSa09WN9sM957RWWy3i/jCMyULMMRmG1G7Fi6zn2VS2vTXc/5nZ1JkPkjKfWJuCMhZbeJYedk3aFSct0r9Q3wiD7oC0jcIQUCgUAgF4gLUiAQCARygbggBQKBQCAXyKWGlHJfkNLuLM3IzGxftz4UOeg1vVDTONKtB2jKFkNTmunilapZpI/R602pImsNC/tqPCkp211So7rDg2edK/FhAxJ97ZiT50uu10J92hqnC1Eboaa0YJvGdU6vqZrVJLnhI+ZI3G1E4jyxbii0wTl0Y/B6AfcJR7UeXeALE43gNesvKcZeQ6qAjkINqeDWVLVWQh/oqyYbVrMj+evebSE48a17HWUxdZMsvQD7skY1pCZoMKuHuDfc2xT4Ovi/LKPIOFXqzCrv6v1xzN3dOdS1nebKKjTGJ7S+hbWZ/cPpJmZm66s6StytS6J78QjwDZXvT6r0mSeNe3BniDn7dVYFamNKQ8oq+yZci8hWpKAhbXHa1BYoZt076Em+b1v92zY7ShxGmYg7pEAgEAjkAnFBCgQCgUAuEBekQCAQCOQC+dSQ5mzHP3idJUszMjM70q1P1NTpeOgFGrY7PSF4B1cukNzLm1UnapriGPYHQL5POTPjGGlOU0JTmuO0qgfUKfyJ4WMkruuZ2BLXjVENqWKmhDbIaUjsoqKmxHzdYhegR2noCP2HQfZycT2dGlJHOJ9v9Tod94VkfEHD1xOtrdBYo4+s1ljaz3BIVdDashQL7lNfaCUd/v5ecd27FhoSXJTEMielIZXR+QKthDYxa4ckDTdvVas21Wl/tSHay66XdO8NNoraSWXXJN7cmd5B5XS+4Ki8doJ92QQNaYtBQ6p0GlIZugl/uVL6jdsb7kPnVhpv7My+yX1bWO8O7hx5HxoSt39fSHytQ0MKBAKBwL8i4oIUCAQCgVwgn5RdarCet40hecDSbkfTYaYbh7odP0yHl51pjxXXw0zLl+dWDpf40YuTUVzTB4ALVDbAbJKnDkl2LEBMh3K3F9PqNTVF75+fvTxx2h5dNU1yI4crVdhrerLuCU6OVBTjlY6a6oXD77/jNY3b/q24nj4QT1SDeKGn7Pg505QF9jMrE8ruvdd1+N2Kaq1v3unCNrCMoWGUJ0YyKnh3e0g+rq7VUueOvdZJvLVrNxeVU2icUepslqLsPFX1j1Z6onaq1GPq1ModwS59Hu6FUFUZdjlmu6Gq2mZRdjTu4XnhkUHZpexyPixxiqrqkGxcOTRVqux7W8vxfjin6eYPFtH0PMjaB7PMYaOg8N52z7ud5xreexvSlyWOolzEHVIgEAgEcoG4IAUCgUAgF4gLUiAQCARygZxqSCw09nwoSyFhB+TlHGhGFw77L4nHY/zs8Yv/Owme0b8dOVo1mEEDEvGkepiKB/df82n9Y89fT8Hxpox6WO7s9mI9LIoatNx89gVHF9fUvKgh+dLnXtjuxRqm1JuCW/dCNXOHV9+TuN/QZKTEhwYoof5ejWo9ttALOtRNSig4K1tYm9lqzFlYW5mUVffcX5+nqpWE1hraSdYRpbQTd1p0gZDSpYPGqiGVsg7yyNBNzFLaibfMYelzH2pITsNozfeWdRQZuolZWjuRmPpr6jzIKneGbuKr4zN0E7PdaUjJshzdhJZKRr3J7QU1o9Iakh+hUeqn2x8JbLaadBTHdkvaL942tGLgK2rYixg/EQgEAoF/ScQFKRAIBAK5QFyQAoFAIJAL5FRDYkOHZ20x9pojJJwdEPuMUprRU/8tsU10a2hIdPQ5fpz725M01zhMm1umX+CErXk43pUc453RzCJjLMxs3kgJm+Yk2sOi0Tri+61BsInpm/DtVOXY+ZHZlwS7HMMIht5DE5GpR7U+eGV3vne/N6U0JMReK0HvzXrMVfBxzy767jqBt983QztJ9ZzwsS7O1E3MoJ18QLqJWUpD2pplEwO9wGsnpfpN5GUzdBMzs32ztBNOn7AKxFk/WXzzTjvJ0E3MzN6hONK2hXWJI2D3TzMOqcJN2tjH3pHcvjyjUnvhUU4fEg6iSf92h+2z27WZpd8s+5JKHEW5iDukQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AI51ZDIxPq+DHjXqYwiIyS8N50Z+ozMVDMys2X3J+sZOIJR92vc162PP1Cfd8EAfd3nTz+2uJaxFWZmD9CL72XEXuNAz1IBD3UNRK+PPlhSy1rVSDzkwKT5qLpSnyZLNzErMa4aEpjvv2Evzsqu1JC8h1spe33w7Z4mT/mWcdSA007w3ttQR8nYi3J6TkpqSKJdfUC6iVmmdpKpmyAu9UPhv7FZuomZWXt8dqKlZOompcDfDXcgGbqJmdk71E78ceDNl6ObvIs+Nj+6ohU+O8apTfdx6uRLvbJbU2fUd7DLPfEuvijPif9lxB1SIBAIBHKBuCAFAoFAIBfIKWVHeCoHAwIGaOgnvXKERKqUG7Gn6Qoljqiv/1u4AQ0boK/rj+nFAaPwTBx4wNEDnrIDzQMqxFvmbLCukmJsVQllR5qqFGUnpiQo6a3IsIkpbY3iUQ5JZCZU1ValqUhNiT0KKQlaB5VzBGWU+LYycDl7/U3MoKnMdkPPtNrterfHsJfHlKKpMuyXzLAXJV/zgy40/h+8y73w4fv4ldyJj8MffWucA4z/14BjyjwniDK+H3uDuEMKBAKBQC4QF6RAIBAI5AJxQQoEAoFALvBPoiH5QmOUPi9WDeblzcmIhtQIBoyQoB0QS7slx38Y3cLazOaavq4/ptRsh9T4idRQAwdYvbBEtley7IJZ1oy9NLUTuk/KEgfwJ00F7ekR+xJr2rWkRiUIStW1UktwuhG0qX1Qjy3WNSzVBo2fdRSpI8iwWWFZcarMmHrgHoNHgRjH5LWslI7FN1uytHj3r9qaMgRi6hQSl3zNjNHcWWjNsIR+48My9oFok6HLUbdK6VjvCxk/7Uj586DkOVHG92NvEHdIgUAgEMgF4oIUCAQCgVwgLkiBQCAQyAVyqiGRF/f9NxjBMFM1JG/N8+jFZ0jOjx03wwgJUzugvuxZgk5k45Ll9AEfkdSjpq8rdkEz8Tx8P6mB4X4v0LNUg4e6nqx+tlRSfXcV9LFuTEQj+o5KaUjSKYVpGqaTN6zR/cMGPhiylupnpY4C1kL+TMYY7M7wEhLbHrx36mlZikXqy5Ohp3HUQ7aeBvufzL0o8RWGnuZ7wVJ9YRw97vS1crSCCmqbOAba9mz3n2VJLa0cVc/FOCb2hTGW48BLZp0TPIIsPY39cSkrJ+6FHEc5Wlq2rtjWfdBtKaqyNeqDFo2AuEMKBAKBQC4QF6RAIBAI5AI5pezA+9gqt16iqZWDNH4goYWmDzhdUtXDYEWNSa/i2g07IFJ2nqabYOM1N1df1x7wx0t7bLyf1K24n+cKZ/ChGrYbntRyH2aLJNdpAZ53WbJcpZlM0tBMvdetB5IHafime8DqRjwYk13VOrws4tDEGak7UxtajvFxbCyDvqTJU4qycxtFx3HGyiryE8jiSUBdki7Dy3i6MuU4Dspuu4tLEURyFLRjwr5sBX0pVk5wak+/clkkagLsC+nKTPoSLNYHRV9mUpdmJejLctoi8A3OoC8zqcvdvOxeFuG3iLhDCgQCgUAuEBekQCAQCOQCcUEKBAKBQC6QUw2Jk0TXtrA2S/n/TDkzWYM/v/+aT0vcOEy1Kj/plSMkaAfkS7tTmtHtOMQpGcebej9UbNxedEXZdL2GdZUvFtep0RsIvXS1Eqmseb1mqDbvq7lttfp/nNetX3H93mKICQU8sewF1ZtMJUt1o16a6mGrNd7ofJOw/VT4PGXOI0hpSCyBd1X6jSjZb9yAEn4pgedRZClZKB/HOc+pI1kTfGXSiZltcTpKWUX4mMLLj+ofOEiJU3ZSfOW91NMytLTdxV5D2l5GK0BqznGGnpappf3PAwDfDlBKvcn4af+AWgHMwjooEAgEAv+iiAtSIBAIBHKBuCAFAoFAIBfIqYZ0NGKvcryOHK13HIs76URNocdh+gWq/Tx/+rHFtR87boYREgY7oAdMMQWxPe3WPF5ywf0Qu704AanTtUngY/ZscT1yLUZtQENa6TQk9iERHLLey7cToTXqtbb9NbZDkmAhnqjAV/KjOEopWT01rEmWbQa8hZS+UIWzTYK8RBlFFAx+WXhEqfY5t0+r0bC1cyUUKOnJ4lFwL7yahZEk0Iyop1W7Pe6yGU1XG1oOqRVk6mnUkKCtbcBBbtjg4pSd1PvQ07xWgn2hnRRjbym1pYzxLCkNib1pbqNoJ8U43ZPlda5S6o07Ep64H1Bvmln0IQUCgUDgXxRxQQoEAoFALpBPym44br3neB8f3jCTcPJUFR47BX5A85RL8DTciwMwI5aTXr1rd8oOiKXdnqbj8YN6omeR34sLNHVqz8clPtmeTIKn8LSzNfSEJI1qSDvUIPau4rQvmod/WGCDfVKxFSWmqem5HiyyBhflmM6aLsskdTBpXu/W9HdN8Qg8JQFCJUVl2oEaruudcEYrrLcmyTxL7X0WcWime4GjgG2S9VJS5QDHUbZZgceiBN4fBY+APxyyNxnUpZk6wJuZ7Vzj3k+KsuPZ6d8PicMM+jKj/N3MrMs27LkLOcO5jMLzTPoyZSe1LctOikdSiixze8NWgAzKrmMJyu5tWAlF2XcgEAgE/iURF6RAIBAI5AJxQQoEAoFALpBPDWks4vXOn6bAuRDUa7ymtAA5KAQrYVH0gK9hpkJAdcHrQhwhQTsgz/eW0Ixq4MXj9qL2rPmSOtcelPiYxS7/pKRsJfSbQgtHt7sjHMTSVeeitG6E+pDMtjqJ181ywgo1pJSQQu0k46iolTjpiqM3+tvf9LHu49r55p4fAcu8qzl6Ax/dUjt4t2szK1NDIrxSkT1FuFefgsS9zQlH0M+ySuB5jmTqadCQ3uqhWs+bFJX8e0+NJOFeeNWCP18ZI0kyyt/NzNrxK+skJmpIWUX4KTupDA1pPYStrWsgdGVqSKWstdyRwCooa6Lyh3e1XP5uZvYWJsiWspQqF3GHFAgEAoFcIC5IgUAgEMgF4oIUCAQCgVwglxpSty8oub1uq9MhJhymDy6w8t9rSsuRY88SieOX3ZpsMJlkX6+fxSqbqR1QCc1Ip6HLXnzKfiu5c3ephmS/c+vpmkIbkjDo3EGaF7VHr5G5Fq3n7VhJMbYGt05pSNTePCPN/a/RcCDSTtcaBO2wdjUGbDj5bwk4cna9eKT6jtBalLJRssRG6XXuKvvapC+mlIbk1KzW2KcBGtI2qZ8tTQJufwkbJQ+eM9VeK0E/1opWulGZPVmchZKyDiqjM8zpRh17rZMUR5Iw9D8NWSqWme4Fj4h6WrOL1zK5pkLjzJ4sHkVGZxjtpKC/+p6sThvwW8aRJHiqrK6wvUHcIQUCgUAgF4gLUiAQCARygbggBQKBQCAXyKWGNN4mSDzh64mwsq4jCOpJ0GDmeB75Rc2V7Bfa2MJ6d/CMKTt3OILdjZCgTx96rqif+b24zO6WXKcJ4HsfSZYzwYmz7cUj1RnFhhvY+m0ck8yGbsAc9YXzj9IHN7j1u+oxlz34osSeDke2f9KDdQTFKo5vd2MwqDJSKfQKTQ0PEZpRM7S2BZaMLFk3F+ctR3GI1sOjIDvv9oaCH7S1VE/WrteSAF+Hlfg68Cg8qKe18bIQPirqZ8u4k/7kbGrWXKaqR50RmozTkHp3UOO+AygasSerjD4k/3WpbpV9SKurkkdzrH2qByvVk5XlZUdVz+1NSkNSQ7ouXqzK6MfiEZhld4XtDeIOKRAIBAK5QFyQAoFAIJAL5JKy+7z9XOL2lowp+M0XLpHckpoh+scPOEpsWr3m1uvU1/T0Vl8MTaogq8QU3E1XjMj0k14xQoJ2QCzt9jRdz4mgEe/X8GXHVGFebMriwx99HXIVdGc6U8Mn7eTdrs0sPS3Xj+lIWTmRFvW0A7golnnXa1jn6Nk6Frkj3O6oqlSVMeAZl1qOVQBF91KVngcv+Qe8gL/NpOz4aZFDdZwY9qXdMbqnLIHv9JKjekDZFTSUb0CqzBux2Calyt8PkXjJKp0qrHwyjyKLOMS+dETZtB9Jgufty9fBKA5vKZVFGprp3nTBV58uSX5yMKcIp956quy7nEJ8tzco8+7WU3m5HubeLMvfyyj7To3e2AvEHVIgEAgEcoG4IAUCgUAgF4gLUiAQCARygVxqSD3vUuLy38b/R3Hdo5WSnA+eda7ETwwfkwRT2knOGkDwzhupccGtMarX8FRSuUprnXrEpydPxrHjHCFBOyAp7YZmtGSGxl4qIdvMAlmvptXy+KEZvVJ3kMSPuQe8OuNIffBUPFeTFyoKlg1f6o3nxUfV7WSt0/2o/Tl5ltl49yj7XuDGMrPKlQXWNT6gjnWMhrN9eb+ZvbjNqXMz9bHWxHEmPBIPlMC3dgoOjmF45VyN+eZ9GhJqAa/q9QGqWP1Y3ux0ozUDdOaCjLE3M5uHL5PYKLEVgKpFxugNlsA7GyWOsa/ZUdDHvqHhSmcplaVimWFvqDPSRsnZJqUslChovsvS7qwjocbtfusweoO2STKKgxpSibJvj9CQAoFAIPAvg7ggBQKBQCAXiAtSIBAIBHKBXGpIdqeGnRoTLvXT56mQctgAtUap65koKc9e/jHJzb5AOf6mOWDGPZdNPhc8rOen2w1XzaKuUi2LPmbPFtcnY7a4jB030xESZmIH9DIccThSwrPv5HPRGmIjPX39cc1tu0T/n/KgqU73yOazk+ABPDG1EhnpkdV3hKMcgJ6SUzT0e2pmVu89iqCtNUNG8apWqfHtg9u6AHrNqjo9f56zj0q8dWq3JEjtC3vg/N5wXzJsk+o1daw9L3HdRpxfrh9qCXpvylCxrAtPKHdM8yCqMk5ZOYm8Qw0py6inRlMYvdFuaLKnHGPfYcF7+mC4WvmjoIpFnVGULGhG/OgK7pgLjTVMAtQZswaGQyX2Y8vxMtSQZKw9NKSd+MqWo2LtDeIOKRAIBAK5QFyQAoFAIJAL5JKymwlqaqTnEmB3cszJSknUjUni0VXTJDe3Ui2iF43W6bOvjz64uN4Ai9wu8PHwkzfpqDwMnMTItc7I5ylTPIkYk169azftgLIIMDIqp7Ic9RNurW5MNrnt+Rqbxk0THW3yBzxvipvK8hlnne6wZAmKrvtZWpdL6rN2tuNYQdm9iI3ydAzpF/q0d/KnzGjNTRNPqHQsJfArSXbQed5TU32QQ9m0e5mBw/4iKVKZFTy/HGVH0jDL6ZznU6rVwdlNzYYZ1cKlsGanjZLwZaSpCHcit0bZN45pcGVim8TvaGpiLyg7fxRZpKEZ6EwMH9hWq//n987n7y3uYEgCWSXwPHNxVF5iqNEUbZR672iZsluJicplkIZ7hbhDCgQCgUAuEBekQCAQCOQCcUEKBAKBQC6QSw1pFuJGx2vWTdJcL9j4VzgJY+RwVV0YvzVIedhlrWqK61IaUt9dheK60wKwzCxr9TFqtVeWKOX2tDL5W3K2nudPaUafRJwM4bWH+54qqd/YpyR+9SHY+PjPYA2nwHLEhD9qDi3A89a7NcZ0nGmPSTzG1ILJhztxAlEr8UfEcuajQet73Wj58d0k9biNkXjdZNT8io0SP1kWWXsFEJoRbJN8FT734eSNENCgSS5z2knBsuGVrFpoI5wi/EKPZAzMs6btFjYNJfwpDckLOqVaA2qSZQkrp6FucvBQe0mT+N5tgHVQVgk8NaRefooERMjX2uqoDRnFwREkBb4S9TT/O8NvP85kL88OUB+0GghmHZa5EnhMzuU+8DfI/4JyX/YGcYcUCAQCgVwgLkiBQCAQyAXighQIBAKBXCCXGhJ5Sq9KkNOsARc8yMW9wJ+T3+3UV7WfIQe63pAq9ImQ2vZcK2UU/OlKF1NhKSDO6sKgAsOB7GIH9Akkx2s4ZeDxxfUEJGfMQBPQBDzXHN+xklL8EHsNAEfcHcKE00aOrlMt5OPeQ8nMuj+KBgnXb9OAFLs5MtQaaw9txNsq/QEeS49sgOcSbZRWev8pKlnsbvFeVDgqdW6yUUOmupSOK6n4PZ4WfUje1IqnNDUAUfi4LyqfySj7GavqNcmRJJt4jhSsZfCsd3sDzajbKBVA/OiNqtmYJ4Mv4uu7NPZnODt+6CImvUc4xRfZYS3H1JC28pzI+jXgp1WjodOQDuq5FCmNxVUJdlI8R7J07NRY+71A3CEFAoFAIBeIC1IgEAgEcoFcUnY0lPGUC29ieUvp2bKeoM56IWbJb7UbdNkG5b87t2nc6GghUkI0Cs9yDuZNOotc/THWIZea9OoZJNgBsbTb03RPzAK/dweed2oz/uFpt16OHAkOX4yO+mVQUZXj1xTX59tkyZ2yDOXMoKaWuarqUuSYP7/qWM4MFm7G0MQh/j67SJ/3TpTepmyUPJ3JM5UEx4nJEiXvVVfpGfYp+01xfdQs8D5wi58JCsa3EWR9UmZmff0Jh32Z2kM5vEftjCR4ABNh1cHLsk2wSvjU93NHDaemY0Efi/M5Ha2wbQWkPTWVmpaL2Jef7xymKTqdr5nvpi+TskvZSfHXwgPnT2eU1rtj6gdPIk7PlZeFdVApIyfv8F1T2eLD9hhxhxQIBAKBXCAuSIFAIBDIBeKCFAgEAoFcIJca0sUgbWc5ipkO7VmaEll7Os5zwuG+m3e/NkuXO/q4lC4kr4GYOhb56RFuLypGIHmmhn7SK0dI0A5ISrupGT2IOFW363cyS6ExE6L/bKSu0PCitvcla7tPk7/VsBnOQeWoNfU+gHy2cbzqH3fbuOL6xT+h9pnl8H5qrZnp2cpPHh/mSKdHXaupy1vdJfGn17qpyTiGZWh1oFrjPy1qRif2wD+4vXnjzO6Sutcu1td5xumDLH/fitEPZY3e0CnP/nSqOnsVUs9IfPhi5wcE56ZF0EqoA3vwO9qLuuMRyXJepfoZvcSZGN42KaUhUY/lQBB/DuGoMC3Xa0gcvdF/G/bfhcvQV5M1IdZM9bU2B5R48B4g7pACgUAgkAvEBSkQCAQCuUBckAKBQCCQC+RSQ6qAtjDStaAcCaugBegP8uwotQRqPczvLUoMExYNowa5Qeh3as/eIi9bQDN6pe4giR90jT0cO54aIeG1h1SfETWjlxF7zr8GuRM1HNk+WeNzPevw+yW+zO4urrtPgoiHvqOp+PA8+061BiqEVbuWGVoq3YWDvGepi2/HExWojXDEhAe6yPrBHuiGZPmZw1XU+6Ldpo+9JVlugLbGlh+e4155OLUtkhdq2PT5ZD3BPiu5ScvH6YPvdOs52/HELyLOMi3CPg1Hf43rXTu5lfoicay9hBh5Uc74do61L2d8+9wdGePbN3Gf2MFI+H2CVotRHN3qEhv31+aSAAAgAElEQVSlQfBJaoffUP+jmTU03ayEjdJB9r4Rd0iBQCAQyAXighQIBAKBXCAuSIFAIBDIBXKpIW28SXtBqh5PrOPbw9KsDuPC61yLzEr0GhT4Oog9j/wuctwozzGXMILXEcfsFwDFTJv/jWOSvfAW/2Zmj0FUemRz0ujTNBFHhdHvOkLiaSTZscVeI69E4ICHouvHSTAjT9bXGY8mmqOeco0Zd0vKZqNnI0vVglpmw7jH7pgmDThHUrfZF/Wx33MaxrQNeCIqNuxWc0fStV5TN2p4zmnJB/RN+47kun17q8Tbf56sH8bYBGoA7KH5RKtk3WkcktdpeFuHa4rrn227WpM/hJogvWtP4IkLiKnyOT2tKxQb+PoNHP2X4pojSWpnQYN5Klm+DE8/dvwQft8Gsz8LYy+WDEyUlOftWMltnqb9Wya/V1Sy+IuUodh0bq8pnOOD3DeEGhLHt2dpSPzmZ/7WpcS28hF3SIFAIBDIBeKCFAgEAoFcIJeU3XftWxLXj2koro8b82fJVc3CJEg/MRZ3qr1w206rdXOMTDPKyStQnm1d3Jq39LQW8R4tKBldN6KjxLzlb3BGN6TsXp0BcspbtnAUwhqOtfVmOyQweKNeg9jRdAPxZq/S8MiLE9//z9nPJXf6HPjcOAZvCVIsqCY55rf4VO7/5RpOH/OR4vpmu15y675yoD54ot+Lh/HEpFgwLrSjo1R/qKmzLtaS9x/ZV4rrPt9eJ7ntt0ho97pzsxRFdz5Kuzs5urJZv2Z2U9WXJb55W7I3W2/opg++04BH3ZqUL6knzumtT5ZjNVN51RqJvaXUuRun6IPRGrDd0fukeFnmTSpKPknOfTlJw+fsuN2uzSztJiXUM43QeFbzqBzfTxoak12G29zieuhm8N34XVzpfhdLjZugDVe1n1JNL6q9QNwhBQKBQCAXiAtSIBAIBHKBuCAFAoFAIBfIpYZ0+ytflfjZwz9WXMtYYjMbPkLrvoeOSESk/jtek1yHV9/TF/q7huas1yugIRk1pP3dGpYZ22r1Ov9a2/7FNUca02qEGtLC+UclASjzlMPPzBYDSxHHwhaT4+cICdgB+dJuaEYDP/MXiT9vvyiuL1z8R32wSkrW6DSAUhY4NYjP83Q7LIrmj9d61O84jXLhzUfpg3+CJ7Z73ZqKDUjzdlpC7q2GzvvMPZL6qX1B4u7XJlZJjdCMdJi77gU/qbH74x8+r+GabyVzpqmf3b70Bn3wN1zJ+wO0l6Ke5s8vnk/Q1uxUDb1907WqCY9vq60Bfnx7BVoDDCNJZrjvMFVSHiErlg/z/zBac4uH6oiMaW4mxson8Uyp8e1ezcoaemGWVgSdOIrS89rD50s83NWXt6EAC0nJK1mlrIJ4RFn6+N4g7pACgUAgkAvEBSkQCAQCuUBckAKBQCCQC+RSQ+II5/kjE8J0fr2Sp91GqBDkLTP6t/2b5PoN1br/3kO1MamLa0TaD2zqFgw8b3Qi0ptoRHodzP5rdkhxvQA9GOtmoe+lwVqOKQs1cRS056dL9Th4+xY2EJyg4UjYlDiNxvcZmalmZGb26WWu3wbjGzg6wXeysB+C/Q+fQH+NH1my/DrtmfmG3STxjJ+78e2QTdIeSwW3xj511BEfHCd+6fm/LK7v3PY5ybWD9rZoYrJ+DEdAGys//PxE8vaw/5l5ofaqfd++Xlw/8SfMb/8enusFb8VDOyDqH/58QuNO63qNofF96BuJ2PO5ahUW2bvWc4JT0DAqfTa+Dt6Yp1Rn3QiO4jjerT+uqcdtjMY7XEydlzY9clRZQy/MUsqW7z2q1xS1dYkxemMnNKSCz+EI+L07tBX/IVm+NVwVJ76bPUHcIQUCgUAgF4gLUiAQCARygXxSdlPh2jvNudzSemeoUl7TXTwdUxQ/NEBruXtUq3eQp+zam05z3G5KW21w3kGrG5Wye28xasT9LTJv4UvF73rLH5ZuFxBnzcDlzbcvxYXvyNkakmLxrt2kVFKl3Y6me2uiph7LcKpOUXSIO+GY1t2UWDB9FT49U+75pD5YKGFSdKQ63T51RVn3RA2/fJpyXj9Z8s0kAEXX8BRitybVwfde6+11vq+5W3tfKfFNu74u8cYbXOGuDqY1a2rAP/h6YVK+LAB2RGLNYZrCZ9XteqXZL3c+RFfYXfoqd+GcdqXec+H0z+pmT4jxfBqBuBNKu+28ZDm1hzraPwIOb/MDztGblJ2z8PkfZJV690EMi6X6ZNlrtPKTH7XnJO45x+0b9ullDGPOsgvip9yF9fGORnyp1RGSwhyAPULcIQUCgUAgF4gLUiAQCARygbggBQKBQCAXyKeGRFuSdx2TOQck5hyYp3R0Wk+Npt6rUW1nZXd9rpV+WqVOhTDbini9W69BrpARb92OJDULlnJ7zjlLIzLT0luaymBURXdnQ3IuHgrO/6zDdVSCn/SaGiEBOyBf2k3NqICX9e4/ZyBXjRESb92iJaZfsR8V1/ff+2l9MN6PNXndKEMzMjPr5XSjBzV1c51OUb1uOkQZZ9szCRMZ+Kr+0xrbG0mVgeyFy4cU19+2f5fcU5PP0gezlHuhP7/oa0M1wZ9PVF1gJ+XPIez3qNHqcXUZxgFfuvZ3ScCxFnrq2Uy3j3/FQ/nt8OfT0cjV8u1cqOErxyd+YA+YlvfPnIv37svPC5wqTN3XK1scLwHNiCMmXLdCPfpDjoOGZH6yNmSsrIYQ6pepiRJHIHYV/rOxy6EhBQKBQOCfFnFBCgQCgUAuEBekQCAQCOQCOdWQyAb7mAwoeNitruNgIboPFtKbvwvicswuPBdM3ngt4sYW1mbp98p+D499EbNLwDO+wzRVj4c6zr9yPMZEt71PYnL+Rz3lGqtgl9OIMdLeDogdGGTQvW7Ua7zmtv1U/+/0pVa3SnzPZNd/k6kZmek5BG1twJkaT03MVH7bR4WGsRMfkng7eo0muLY32vpjAIPV+VYpWCx9rcc3Jf7Ba99JgmvwRFN5LtKIqODWPJ+wF/4o6zGIAJ/PQRe/Ulz7ERFmOnbczKz2cfQZuo9nJ0ZIPIGeGW+8w28Ke428ojGMY8gv0XDVWD0b77bLiuv7NlysD8Y5r2NgZiHJs97vI/TwdohP1/Cg0ckenwD9r3YO9tSPb9d2y8xOKP6i1PbAP2AfF9cmvVMvpma9l4+4QwoEAoFALhAXpEAgEAjkAjml7Fiy7CkwOuTyBtTHnHdIiiIr5tbQc/ntFta7i+mhmwXShv4mukQp9wA34fMUTdkFGh5dl9zTn4+ZpKRYuk8Cb+IYvCWo+s4qJCalwtJuT9M1gbb6Ulul6H79EPgxTyFtLWUH5KiF4UqedWxYJ/HDHRLjnhO/rc7mi74joT2CV/H+8NdVaq4NJtNOGZ/YS1+F2vnlVw/QB9/hp7feq7nUe+U5Ll7hmqJllNvTgafpJGCeM+e6mvgB0zGflcNlQcvNde5YLOXOopdqELPMuK//h7GaWzVeKbo7TfsK7t6VUHY7b8d3Ei7jasPP/c/yGUftOb+zaMcYY38qrk+2JzUJKyrvo/QqUlke46kybzD/POTn7LjiehamXe8N4g4pEAgEArlAXJACgUAgkAvEBSkQCAQCuUBONSQQvrbIrcFPp8qoPUPKYluyp4w/KFC78nx1lkZklmbGXSkoxmlwaoTnoLuf9YakzkT578ed4nHKshmSM0xyNZRyz3ZV37T8ZxG71404RoF2QL60m5rRfz70BX3wODzZ1iw7IBDf9Yl20utZtWqaYqdJPOSyJP/0RH0aFvjCdMjOcVNHm7GHl1T9p8STfv3ZJIA8Zk0P4R+8HQ3PJxSUd0Up7ji3vkr1jfP6qHbodaGzN+tchTYUzLwuhNNpNrogXtZQFGIqLmwN8BpHPXU5HeQquuniM3W0wwTUrU/YofHmn3T3D1Zs4jvwyhd/U/gOXDH6AHx20IxGHa6WS6e5Te7+OHTdZzSc6zY1S4cz01+gwSzzxm/MG8O7S/xn+2hxveaZg/TBHOmxB4g7pEAgEAjkAnFBCgQCgUAuEBekQCAQCOQC+dSQaIcyx41EXozxyOubNRZrFFr4UOEg3+v7h8hmZ/U0kccnb+xZ2l6aUko2rRN5C/p6TXU7WUdBf8yeLa7ZpzAGzR/dH3UcNPSNZvSJTMW2eQadHVc1iD/RNllz7DhHSHg7oFSf0Tg8cWavUb2mTtH40CcSzv9p9OL0PE3f7GS3F+znoCY2+FsaT/n3RET6xAYVXXYehXNmToMLGjSXOr/cfO3h+D7coOHR56igM841kbGXqOr+Jv3jDF1o0QqNvapFzaKUUuu/LZyQfXQrjbt4TYOaESbVT3WNSPea2v9MWgXvoDvaaTzRrddQM6J66M8ZflboFezsPi/0Bva6WPXMj6Oz7ZTV7kPAd/QtHJL/PEqpWqJ90v3nJA0b7GOI65OATYihIQUCgUDgnxVxQQoEAoFALhAXpEAgEAjkArnUkPrcpvOel893Xl5z8OCFFRq/7kZzr+yrOTi02ybEYldHbQqv43euKx5KXajGrWlHNxQxxhbX9p9fXNfZi5L7qP1ZYj/WuHY23iw4Zy8xLUMzERlydn55dY3eV+eBoK5wutG6m3QuvB87boYREhhvUNqfrj5Znl4vmSGPvSDxn3clvROdRqlW+Cu8ea9CXtdBc+3h0Xb1STdLfMct1yXBtRwLcbO1DIw7OAXKivP5u6b/D/QY7D8kPuh+HS3iR4Lv1NPHGtDa4hUNqrFZCut+yLE/i+fMYb73hd5pNKg7OVn+daA+02NwR3zEPl5cz3/mGH2eiXhejKe3Jv+lKNVt5/Vkvls08jjd6ENXbJNUlj+gmek4d3jXzdqhcZajJ3W6w/w/QPdZPFT7t6bZCRKvfKbWJxU/sLIRd0iBQCAQyAXighQIBAKBXCCXlN0PULs6b8gRxfWCIXpL/Jr1l7jQWFNcv/c6OBZSdmAzhMLbCopO2Sazzm5Nig6V3W0GJIWXNV2WSe4wsUUyGyQFm2ZH2Lziug7UQc/ZoA5mtLA2s52gojw9w3LmLELCTAtZTwUrymmty6/rVlx/1X4oufvv/XTLf1sORWcmNN2Rj+mYiFnbjpO4nWNR/mOepFJDRr7u2OKdyvzZ0Mr/lnj+caCFZnrLH/1cU8dfn8RtHtRC3Z92GSfxlfPuSYKL9Gk2/E7jybs0Lrg1S/azDK84S7YP4kPdgyvICZGjy6CpN47Q8uvnMdLgOWdV04A9fPEV8Hue8foDXnPedvwDvjDyrcga2GCmu6PHlGpXcN0Ml1RnT9bt+Xt8E53711ytEE99hz2lSnMyltILLYqByU96jtTMntylsXlHKcope4G4QwoEAoFALhAXpEAgEAjkAnFBCgQCgUAukEsN6cLpf5T4gqFJvKhKLc5TGlJ1TXG9orq35Fabequvty4Sb3LC0DvWVnL7mNZVdnaCU1fTkt4etlriGsfcHwwtpL/9TeLa1RC6PC/L6lNwts0ufhH0MzlmXxbKEl6OGj8a8TBfmo4REvPHq4DwDbupuJ5yD7xdoDdZUxkjJGAH5Eu7U5oRpJ1b3fgM6mP/BhlicUOilhy6StsRrCvsZt7VkRny9Wp3o6YaNPxF3bji+spH79Hk9Ro+6g6DylTWwGwzk6LdYfsjSduqQRk56ELNLl5SpSLqUjtYYn5nF1lip7MAZdMvbVbBqanBiVUNOCbGog/ORZI7R8Mjv5P8RkBRa+38djg5BzZo5x2efLaXofb8yOn4lt6v4Uonc3HUe5bKRf2vC8fWuGr56b0/IqnH4c+08UEoUjIhg3vMGv7SiDukQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AK51JAMNv4Vjp8+vFZHczP2vhg7VUKytZXqa0MNaYszPXnb2ktuX9O+hf3cePSUhrRRBZwKPyUC/QOpmGKP0zu247EL1HlEnoqMOHtOvHZCN6N6xNVn4B+c9jN9jHLO38GHN+Pnbq76tXiepqxeI/jg1+uYCD9CwkztgNqBI/eakZm+9ysxwmD6n/T9jH7y+SQ4pRHH+0PE9RpekcQX/vK/JHXfvM/oYw9JlhPxORfwKv7zulqlTuvEmRiQ7baNSf4P+nRb7fF5CQ1Cf3NazzKoUW9Cj129I4k3L4XmUkBPXwHH+HoLazOzxRnxu/w8+OXxL5SlEZml1UTfaYVzsTsarfwYiavUcuy8g7XX6HK7q7geNVutwNxkEDMzewt2Xw1uzXfDHjIv8Q1jryDO+TVnJrPgab/01FI0JrGfa7H/XeT+h4YUCAQCgX9SxAUpEAgEArlALim7Btjc9HFxDd2klTkwOzBZtkFZa8/9lUrr2QW10ZVuDSoEVd9m3hmZRs60RvZV4Drk1Xa+qfESOC57p206H2U5MPMWnvYhg926jrf0pH3gvD1pwDnF9c2oSV5481H6YHGBKmUH5ApUh58qmV7PKo/FSa/etZt2QCRjPE039U9a533qvQ364LF+Wijsve1KDacpVfXM6IQSO/4itRl6CCW9vgi5Hq8yjpNRk0p6u3/oWZKaaJdJnKJc7nT0Ga1e+HF4a631yG21DJAOY1Ey4y1uTeMqfrk2trA2SxPTHjwL+I2gv5Gj6dA2wEmvHcevK64v6qD2P5eYUnYjZzmqeYI+z3YdEGtP4zeHrvsepN1P8L9lPH8u1PARO9utP67JB0C3TjHA/1hnHeGeIe6QAoFAIJALxAUpEAgEArlAXJACgUAgkAvkUkPixFJf4FsF2jgVuxJfGn5UwXa9E0ZKtPHTKmjRDhv/na7keiN0nyzWm4WqZMG3ICbb7pGlE3ECwNGcduqlE9DGG8erJc5d8Pi5zb5YXK/7yoGSs5/wKLPsgDBdc0Cid3RsWCepKXaaxD1P053zk145QoJ2QL60O60Z0Z/pCbe+UTLdmlUQXLta+fbt7vy6ESX6nCt6oyvPXvO7SsmNMS0Zf+LXTuSDHmCLqTQ+htgXDPPs4s5RC8obqAvxG+9FZH4jBms4AGnXrWDnamrICJ1D4ie7cupr7aP4PB5Iltsf1dQTOEdY8e4/jRrk6hG38boRzpEpvY+XeLITxZZPxkZwku7WrD6V93++xB1SIBAIBHKBuCAFAoFAIBeIC1IgEAgEcoFcakg8qLdaWJulXUi8rkKGuTV0oH2h/bRxMY8hi11n9wMf6/PlsqzeRp7TAmoQD3a9U52GIzkasdONZgzVARN3Y+7yPUsxJ+J7TiuZyHd0L+KCW0M56XWOxlOT53q4gzZDDblMuevJsFXxPVlfBw3uR0iYwQ5I+ozMVDMys9Y3JseAGebzHlWboZnaEmTT3Prb7J9DD9DovonWM/2y0zU5Ucfem/kxF/xG4IVSZ4kf48HHcjS3R6neordbWO8u5jfEg988/q1XSqkZqRWY+SkYHJ/B78cJGvYZtdilpknuDOhyZ258qriu+D2eF61rjclDU9MyqLBmjRKBLGrV7DVyvVIvjBgiqfvsYolnzHWCGVsFU6PeqbF6LZe/uOUj7pACgUAgkAvEBSkQCAQCuUAuKTvcPUuhKsuks4iDrJLp/02wHNtv8n7IkTQhCVHj1rWlJnx6ixNQdMuP7ybxHxxnd59dJLkXYadjt+N1pvnCdtrp0IfYWbJ0BUWHktLf9knqU0/89kzJPT1RH0tf4etcWTuYtfSkV3Ht5vHfKJGn6eZNUoruoUv0L2mccqOjRSc9rO/9krvw5q/wtJweQ/qscF5Op4BmG6dhu1P0G9O/MplQ3BmTjYktltStr8UxrG5Uz673Cu4DKOCJ1pSIvQ3RJuQY17h1d+R6IXZ+Oh2HahvBoA46MXYYfJQ+as8V1/Ug17o9Bd8kTx8/qamXcer5FpaCZaMGsf9W9j0JyXEaLj4zoal/Y5+S3P2r9PsuFkYpa6AZiMsxMCofcYcUCAQCgVwgLkiBQCAQyAXighQIBAKBXCCXGtIwTCgd5uj1nSs0txKl21nm9Fl6k5kWo7LkMksXYrEsY68TUQ2o5vgMTLkVV/yhyMEWf1Vd8krToMQ9Dg/6RzYkAsfOO3HEsMW3wiL8gy+D5S7Dxr/j+cl6oqZurrta4rETHyquF31HH0s7KU7IaO+koKGVOurBuqoVkk561REStAPypd3UjDj+40tf03jcTb8oru85CqMq5vyHtYiON2qMKZ3XjP5BcX0tvJp6PoLP47t4bi8bsfq6ErE/N+m8A/3yreHJN2RR3WGSex3aQgHqyGr3QmvR3LAak2n95OZqfAJdMCOjhyWzXfqhqPoQe03iwRu1raDCS5g8+RBvcJrli2gtodGO/3T4m1JycrP/Co/F63xSBbQJTmf87TacuHfg+yCl3qrdptVa/jL6MnyMat4LxB1SIBAIBHKBuCAFAoFAIBeIC1IgEAgEcoFcakj2fcSOxmwDF5W+GAne13PkJPk3l4jdyOBm8OsV3KkOLazNUg4mQotjWoNxfDgkmGanG71UpcnZppY/z9lHi2tqSOsm44WdDT41irSpCe1CvPoGO6B26DVyetSXT/uepK6bfofE269K1pjmnBrXMPhbGl990s3F9fzjIK69e6vGnp3H2HGOkPB2QOzAoGY05qaHJH5iX6d0Nd2Iv4YKNiEZh/Czz+jM+Ksm/Fof6w55Ms7xAl4la6h3KXhlkf1y1DsObZVoC8ccOl9yxwzSODXqobaFtZltHK56x4c3NhXXFaW+z16wwe+E0Y2JUomLNyC3IEMn4k8OFRe/pxyaPoojYmgH5HQj32dkppqRmdmEHUm89Sfag0gt17bOdcFfkSxlTeVGvddQIS8fcYcUCAQCgVwgLkiBQCAQyAXighQIBAKBXCCXGtKtA7Vn45CBSc9AX7DkNTs07vD395KAhG4pTcmNEK5oQo6tLJ7vLWVI51op1vXWuelL7WCJX7P+Ei9w6slLaER6cVudxFunOq54Ko6B8Uo/WpmNFjTCJ45Mlh3P1BR87y49/5fF9U+WfFOTn9dwgtt/ev6do1OXbcq/6z/ccct1STBTtZzUaX5FfXH5zOhjJbUdY+19x9WNGPXu+4zMoBmZQTf6tqR6Neser1iYaFdv4Vy7aYfG3uSfYwjOZ78QYzeiJPV9QNPMyy7PIR08Y6Y5XaXNQs1VIa7B33o15FB8l6oOxRfR59/UlPE7686nZnzX16Jdq4A/9d8ObhM777xOxN6iGsTum2OD+dnQnw6jx/8yImn+utsukxx7jUQ3Yl/hGvYVvujWfHfsqjxSw45uFPwF9r4Rd0iBQCAQyAXighQIBAKBXKCiubm5+f/6IIhWjdskPrg6oTf62VLJ9bYVLcbV1ii5LrZB4v1sS4vxPvaO5N6xfSTe4kilTdZZco206nec3Qp4A5Gyo83KurmuXBtjFVIuHz5eyXJNlm77WlbepnPyo1KD1rU+Wf9QU+d95h6J79k2rrhuB0Zr0lMa+7Lq62Bj0/oNjduCg9nZ1RMnP9MHt7tRwgvf/q/i+r6LPiO5G+/XP/WTXu/dgBESR2GExBx9HU/TjWzWqaPP3aL8zORrk3UBz3I97aQcBfPwSafq89j5Ei8ytfFpZQm3Nsh0BMPH7FmJz3bF91X3gw/DZNTt7rOcpV/fVEU1zzYPnnlZM21LlbT7MyJrivPu4qzJzlnHyFLuI1nK7Qf2sqz7kxr+sQfOEffZpkZI0A5oolunKDoSrn5kDN8dKDo7UUNfbf4N3bXmPiQwSyPukAKBQCCQC8QFKRAIBAK5QFyQAoFAIJAL5FJDqkB5rcgq9Cxh7FzYO/bSscVdOuy5huS5drO0hrTd2hfXKQ1pg2pIO1c65psV1Ywx8tgWtrA2M2tqxD94tp7G91mmJqxTH6Fhv8Ea35gsz7pYRZc7Mc6h+2WJ1tMwUZ+mAa/q1ZC6X2nukvH/KfGkoz6rD5hzswugeb1QL2Fz26TE+qEjLBPnOO2q4il8Va7gCAl9nV7NCR+/4pZDJPfLayUUd6lzdKKEffPLX5X4e/c4b61rcAibqBdwpLxXYaAPDAfn7/SBgy5/RVKX2514aCJsVU2A3gQfqO2Yiu01J34dskqs3w+oblA58fH+yNUgHuz0zjYYy8GvkteN3hilM9gfM22heNDOlXjmbKffsJT7AcRiB/Qikvzd8O8W33VTjRKHZHZDsjxvmOrHk+1SKxdxhxQIBAKBXCAuSIFAIBDIBeKCFAgEAoFcIJ8aUsVK/IsXhlpOmZmZp2W7Ise4M2L/3DRVoi3JVrfehNx6xGvcmm+NMXqltCulgFyWqQnZdrLkXnxDb8FI2IXcoOE5pyUzj39kX5HcQdeukXjRLcn6dxlHYGY21vVhTPmdWgOd8etn9MHjG/DXrs+q/nrJ/OLZcRJfeUjCdd8Iqe1G9IKM/t1jxfX0isPxmo9pOOHfJGyuS7SqyZifwdPrnCeTde+T/ia5lfvCY6bJa1dUWbirnIXi9QPqSzxnfFcNmsgugArzjWQ59nAVAMdD8Bg1B5qG72lCb1rjPI39x8Uzeou1DO53qX6nXl4XYh8YLX+8oxemeG88XvuDnnFjYZ60kyX3KDSkdfdiZIwfNU4rsFRToh8jUarP0OtGp2jqbB3HQs3y1FEPF9efs59L7nTDd3YPEHdIgUAgEMgF4oIUCAQCgVwgp5Qdy2mzZldmxXSq5WN564rb00z4bSNZQNptYwvrPYm9BVApsxT/fnoih3LOdhkuvShJ/szhOtn1m/ad4rrPt7W0vvE7EoqDCUttv0QqxFkj1fRQw5nl+2LMaGoC68XJ66zXMvZ3ZqkP0UQ3BbYGzzKgWR97wGWOj52I1+yo8c+26NTOT7VLJr3+Eo7d16O0u/eXE5puZQXPU1gh2XnJcrxaA9G+pU8ftdpaPsPtI5zZ7Q+0m5rs1qT3eH45Su900ISgfE8d8bDEl9ndxfW5G6dIruJxvMwctyYPx2pmTHYVcMozt9w59KemOoN+faNvohPMRsuBn5XDeXgAACAASURBVOJsZtbgWgNenY2yezhR2RTEi/3n04AkC+b9Y/k7yPnLjrI/G6mrNDx+tB7UFXZXcf3JxTjgAeVfWuIOKRAIBAK5QFyQAoFAIJALxAUpEAgEArlALifGZusqBeSoTPi3RI2olGFIOdvh9ZxS5vY7M3Kl4I+ZRDd5fF+PCs0I5ajeAqTqKtUHLm91l8RftNsk7vbtpOZ9+y2SEtXBTHcGhcNmX9fwaz2SibLLr6ZmxCmw4MVPSd77T7uM05xWgcsZNA4jAMbYf+k/TFzmAlgs/UHDqyb8WmI/6fXr0MtoB6Sl3TfqgxF3a/57cT0L40xqf4Q+AjpIuSmkf3xExxt8yW6V+I0rnT3TndRFWcTvdK4pGMIw5TwJnzhXz4QnrkniUSO0nnnM2D9JXDc2KRlvayrMHWCrJfZWYK13qaC0vVV7ibdgRnGjMwwqQER6zdQGyk91nmfqRbV8Ns5j//Z0Ikm6cjs1p9eXcpcq2fe/FRwhgR8DbwdEzWiU6kIs4f/kMpennRF00j1B3CEFAoFAIBeIC1IgEAgEcoG4IAUCgUAgF8hpHxL1At/nw14JajJetfigzOrfL7wOlGV0b5atE8GzpDU0jeFufYKmaBs/akhCZn/KfiO5T6/FHG/oRNudQ8i9GFdNZtu77584VnMv/HaIxB95zfnEpHoY/h2x6hK2OOnHaX5b+8kexYgJfwaNfUlzFXPxuuO9rqJ9Rtc0q1XKbdVfk/hW5+z0pSclZRVv4nXG+d47/fC6Natf1toZfYrrZfX6NFR2+O3wxkJj2Qf2Uw1/+fFkfMDnXpuoSfauzfMCCOZLpHA0YtcHg3OE5227+kRP7lypnl19oS93dp5eHDXjx8f8T6zfww3OZ2zFLt2ojYuh3Xp7oxc0JX1TjN+lwEfNaDlif+ZSD6ee7Pa4M3rVOELiimR56rCWe8TMoBmZmUwhwc+G/T36kAKBQCDwT4q4IAUCgUAgF8hn2Xe7czRu8tQb/UEYl0PvkdLLovvKKS9vj9iXlHI6K+LWKGf2HAunUR6DuN49dNhfJDXG1IPlXOdTctQsjKJF+eaG32r8sKugJUVHr+kTvRPy9zX3bdJw4iR8L54J+zJcaYhr+v8gCS7Shy7AM13dNlnfP/QsTV5ogDuHTtFjuBZ1rZNhvj7KrR8+CZM3P8zXca0NsANiaben6XRGZ7q0fnAPjSe5yugbV2huBP74yrHJs5/8W+UcL3xJR5S++EX3bm9nmTE/y1mInU3UJEz7Rdw0MKG013RXentNv4P0b73zP5392alBN3/v2L8GuQJi79qznjQVzz5P0/Hbw3YXwp9/KK0nDTrA2TeRXh2vb/bCnvcV16TvT1kI+nUinstxxE/jfMIMgT1C3CEFAoFAIBeIC1IgEAgEcoG4IAUCgUAgF8inhkTLidedflPAiNg1iNe3sDbTKa9maR5ZRkowSQ3JoR1i8tV+Um135GoQw2nE60btjlGOeXjlXImPteeL64/Zs5I7eaNywRV+SidqhZdN15gOJ575ZrHp2P3xD9cly1t7XymppyZDv5nq7Wlopw8NBiMNrrakbHoD3g/VwE5OK5lol2lyMUf4Op1inGZ6PqKfRwF/eb6r0v+Rna/JTYvwaKe+YYQE7YBuduuUZgT7mYdH6L6N/fkTxfUy2MRQj1rlJpSOW6JCyuzf10s87rZfJM8zVD9nu+azGm96Gq/kJ8g+gdxsDRe6M24h9NdpPBt9nrpuKWx3a2o71K3XZuSyxskQ1KKzRsjAGgyDXr1u1P3SN5BSg6/zXXzMrPn6PNCP7RENn3ZvnWpZaEiBQCAQ+KdFXJACgUAgkAvEBSkQCAQCuUAuNaTjP6/2FN7+ffXmAyTXVIDVzpoW1mZmm0rETc5y5l1oRtwp7+aSpRmZmXmZq5fqA736FCSugRJxmCVawyCwtMPhS1K3MeF/K2BVY9CFzOVnon/grxqmWHDPbJ/fFsnPazjzwqQn5aZdmDfxPfytPebW4NO7aj/K0eeoJnbQ/cmHPRmjq2v4Mp9Mlk8tPTPjGPSv252CnfiuhqnBIk5DWmSwb0n1oCR9Ixw7zhES/nXYZ0TN6Jxfa//ZeZ9PlKLJw8dJ7usn63PdujlZ/xJSzpU6ucIm/vZzxXXvS/WE+t5ANKBdA3Vhpu+padBcyj5nQQtrs2xbLuZKwWvIpfoXyxkp43uJqBGxi2+Yhr6nL0MzMjP7yJDkC38uZqOfg7jPI+uS4PeSsp0YIf/EZo1dB1nZg3V2h7hDCgQCgUAuEBekQCAQCOQCcUEKBAKBQC6QSw3px755xcxWOC+vFZXq67V6iJLoa4ckjTCbIO78A/FWjC32lvS7rJXk9rF3JG7v+hRobd/FdNxztetN4Jjl3qZ8ez9T/aD/rteK604vgbvWNiS1vocN/rLFGvvOD3b8kCEn0/0JtzWdrtDcmm9VSvx9N6d84w14poW03y+49QhNjWOotvje+r6gmdQkjm1j3P/D7qxAtuUBGv0r/6ap1ZYNp6+1Mghb9OZzbPzyGei5gr9eP+c1OAnH4PuMzFQzMjP7nRspsf0R7c157Ekd6XH16GR9G8aM/BLn05XuT7975w8kd+BJeo7f0PBDiTfe4ETWiZg/sX6ZxqIbler5KeUNt7dgT6L/LDk+ht6VXifCOBl6VWLSuJ2eLHudpt8delWeYY8mf7YWAjLHRCQPtUY8tAEPzfqt4DvdG8QdUiAQCARygbggBQKBQCAXyCVld+RPX9X4QBejzBXO/PZWdXI7/Y9WStGRwtsCyu4dx7GQsiPlUg5l12VzUivZBiXW9nfEZLGWZOR0m2zJihZTKRrLkxkkIGhsfypKuzuNS9bN39LczXa9xE/8yZnb3IEnTpkS+dJclAZfpUTi+bA/2fnnZM3y02GwM3q67bFJwImeKWuXhILpTI4uZT0FOFsVluzPJwfp9+J2zfzxEa2xHtv7qeKaIyRoB8TSbk/TTfnqJyV39Q9ulvhndyaf5UWX6PPSZuhXjln7LGjcz940SeKaCwsS//DHyetMP/10ydmkvho3uJj8kZHe8/TrFuS2Wzaypjxn0XKgpbuCEva03HBN+fExZmZVpyh9/LFWiR3YyaZ9HSch7jPdlXI/agq0hLzs6Fe2fBQsG/7djmjxUXuOuEMKBAKBQC4QF6RAIBAI5AJxQQoEAoFALpBLDYmjrm3/FtZm3nHFzMw67Z9oDZ0q10muD2LrgOfyWgl3hnrBthbWZulqUy8pYcx1qnQY8Ur3+AIeytg/dZbJvZmy4NSMTqROx7Herir/pqovS+r2pZgL4e2BmhrwRCzbdaOvz9bMeX3uk7jqfh3D3OAsTVKDQlBO+5L3YEnpEPygWZ7tUNlyysxE8+M4kEnDMZJhjtPI/qCf3pfsVonP+mnyhjh2nNoO7YB8aTc1oztmabtF3dikOWDsgockd8KP9Hm9GjgRUs64azQ+cbXOyBj05URfmzBqvOQeHHWuxPPnH5MEDfq8Ngd60+supo1YqVE0fqRMR+RoDeZHytD9h6XcTjc66PBXJFWHURvH2Z8l/qg9V1wfvlBHSqDq2+yZZLkdE+Nn4ffK68386eJ3qQ9irxv1HWXvG3GHFAgEAoFcIC5IgUAgEMgF4oIUCAQCgVwglxrSROgsVS4mo8+OAN9Z1ElbiawTuOA2HJ1Qjoa0I1luBye7BbGXkEqZm2QNPGYnBS1+PLhP5H6dWmN965CELtGEkRK3dUhEgZu3ad+RfQN9Fy/48duYYZDq73CWOSolpCz0yZn7Fi2eEzZIw79Z/ySgtpCxq1soJkBr456/7M7bszH7+TN4fzbHK3naY/XGlao3/fKXif3PlWNVNVqlLT8yQsJM7YB8n5GZakZm2lN27M0qRIx8VTduo5vawV6WX+H7fD6mkHRfnBzkN8bdIrkzRmgTzTNDEiOo54ccK7k5aOxZvvzgJChADUmNnkFchobUrib51vatLEjuYFiB+X60I+wlyQ1DU9xBC3FyzmhhbZayCnvZ9ae9jIdyoIc/43kOUxKrZ0/iGBdoW9teIe6QAoFAIJALxAUpEAgEArlARXNzc/P/9UEQN1bQgTkByxDJrGUZfmQ9dnf5LHgGr9T8yHJmT2aB74e3197Gg6XctaiINV+i+XFNvXFmd4knmFJGP9t2dXG99YZu+sd38HRy1tQZTtpmZlaf2AUNfPYvkvrr5qMlbgMa7iZHURypKTv1VxqPHp/wS9P3g1XN1hvx18l77dWshOqKSYdIPBH2Ov58GqtV63bwhQslfqPCP9dNLR6DmZktTnodlvbXCcoHHaM0Dye9+ndwEV6lr1Z92xs3J+fBbfZFyf1sCaha92STYcdEGytSqvVuPZh2OmMQOxeldSOUS5snI1V1Su9q8Ku0Edth+0jc2lmF7Quboa6wBtvftS/0RTMGKbvatY7CnmcK2lhlxEtgGZVFy7EFhL97vpNmMHJ1sGZLfR7O5X3+8epePsTgjr8HiDukQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AK5LPuuQew50Cx9hvlS9jn/v+A5W274foipE1W3sDYz64ey9i5eOFI6XTUjM+GCp/bQ5L12scSTlo/Tv/2he0d34nntYcReN+LsWYyYcKXQHC/RRqumbRE4dK/FscSdgtqbXk+ghUwKyVm0uhF13rCFYYmsFEr/XnOXX6gbd/0FP0uCB7hPv9PwgiuLywtfekBSs39fL/GVOrlCJr3SZoh2QL60+2e3qGZ0T61Ol730m8kxnn+tPs9kjE2hpuRn3C6BbjIScbUr9+82TD+8E49QS6ITB7j4IH2enRBg3+6oykqrdxMNqf229yRXgVJ6sfviOBlOxPBvHvvSjHgBekB8mmps1mzcUqXcvuujmjMkqBmhtPvp2mSs7b1QJSdmHFNLiDukQCAQCOQCcUEKBAKBQC4QF6RAIBAI5AK57ENq7qJ9SGsdQcqBBdSJvL0OhxSX0p+yeoJSIw0cqPtk9QtRM6IuVI2RBm18HwB7idhs5Hs4wAW/0GOIxE9aMpfgUTtDcn99ZqTEKZ1IXHw4H5nGMX439HXs3MMkPOj3iR3/n+w0yQ04Tw1PJkOT8VMkvsZGF3DzH+7wZnG9uV13PPhGxE4reUGPd/Nw7V3Z2VbPoJ+5qffXYdRJ09Z2End5xZ2dA7XPxexniBP7HLtGP6tLb/ulxBPnfE7/1L2dX0HfoC7h+7nOZH/QNzWcfuZHiuvjJ/63JjFO5gl8Hv6M4XeQH2WNW6d67Tia5kC3PgA5jg5pZy2DtkLUkLyAA5uk7Rgns8rZihXwNKV0oSxNvJyexMFZPYnQjDZ+UjfmIdNxIA+6+KlXzpJc8+G7PdRMxB1SIBAIBHKBuCAFAoFAIBfIZdl3BVxJqt1tbzUnrpLd2NzC2iw92XWHhs2Ow3NVn2Zm1hol1hX+Tpau4aQDfExaIYtmMDPzbhy4914zQF/IW6fMNrXwftY+JvGMVfVJ8AD4Cq0kNptD8tMX6i5GjuSmO44apbzsCg0/Zb8prgdMhycx3I1Jb3gqtEIdTGxJVS+JNy/1RGkp8yZHpRU0s6hO388xh86XuI1zB+KUzhMnKA80dnzibzTpdEyTnULSxW3G7WqUdM/QKyXufanWx3/3zh8U15/F/nPSq6fS3kb5NUu7j9+Y0HQvjFN6+BjTfTn1p/q3+zsLHfUbT1P0nsZKWRLht8H/VpD6I8WV9UOY1VpiplRaKZrNSwrl2oj5Y2ZjAEu5h/nfFbr5swXkzGQ5s1bPp8d80swewSjnJX90nzWc5tnqsCeIO6RAIBAI5AJxQQoEAoFALhAXpEAgEAjkArnUkOZfpyJAD0vKdLuthdcLyiqFxKW+RA0JcYWj9dtAQzJoSDJFEiW9KcK6i1vDfeatHqq5rGilfu+vO3b4NdNxBwtgFu81pIVLUac7DSM9pvocjnfrIvwDmf2CW1MzwlyI1vXJGprFqNFTJb7I3IwGOBDNhj5Abl5eFZLLUjtY/6Hg96KUwZRj+jHA83Uw98cMUq2kymlI1DtOhBXS+PHJmI5JN1BDUpseraO+V1PX6N9+b6DWXB94UqIpffYmJf3HXSOhTHrl8dMO6Hz3MtSM/jpOP5Ajq/TZhrkJJbVP6fPOgs5bcGvqNdQVffx+Rs0Q5bSLZIHtIZSTU1Zhbn0YfkdSVmG+7WO0pt6o01YH3wJCzeiJVagDnwS92beAzOFOZDXL7B5xhxQIBAKBXCAuSIFAIBDIBeKCFAgEAoFcIJca0nX2Y4kPcEJR9f4qJlTvr50KnW3TbtdmZvtJF0B6NHH7VGdAgl0QkfzI460wBPoHxiNvsK7FdSOY4jchKq0w1ZCWObOUJav660HNA5/re0Ve0FQq3uT3jXY/EAhSbL1nvzn0+FQNnW7U7Xr15r/M7pa49nE33vlxSaVGNBOiUoBPf82wbwUflNKQnGIADanAQSkDNPRZTBK37eirGjUn0elOHaEC2hPnfkIf/KAf5z5Lc5ue1vgaHfFxQ8MPk+O7sCC5E1fr+Ibzv56sJ0LLSY2QcKcM+4yoGc08U3tdRh6SnH+djtC/PfVJjTe4c/xV6LxZ1jv8lKkDlQP+aO7bwtosbRXm5WT2EtWgn7ET+ulEKB2GHNy+VtUlQjZ7EKdBVPIa0pqHMKfjD3idKYg3+eY1as3nW7mIO6RAIBAI5AJxQQoEAoFALhAXpEAgEAjkArnUkJ76udqYW68W1mZm3dUTrLKr05DaZmtI7aEh7WPvFNetTAlqakjbHVv8trWXXEpD2pBoSDvXwEVrpYYyR4Exc7SRkzxnJ/PBBbemRlRqCIDXjeo1NV7DD30jafa6HHMsLl2L0dyuLWYuDp8Whjwi6ctAC9Yig4ee7JOeE5lItcChGQScvx+lTg0p5W3nfL8uG67a2hPXUEPyuhDVHPD4M7UHaOMNyRfohz9W08hBX14gcffFiRlk/QRJiZuhmaqQ3pvOTPuMzFQzMjObNOCc4vq4f/+z5PqMWidxl2fc88Bfj9LnW65HsREaGDWlltXjdDdN1riZaupAXTSWU4ZjIKgZsbfIndeL+/aR1PN2bIsxNaQ3ZmAuhNeFtDXQbCF/C6BZyvlH58HyEXdIgUAgEMgF4oIUCAQCgVwgl5SdwdrevNNFV+S6aunz5s7dd7s2MwOTlp4S6WPuDOtEPVO4CTlQO7a+hbVZqpQ4ReE1+YG+BSSzzFKyjPsJkhB9EMO/vqvjFsbiodcqhfq56p8X11fYXfpYTKLd6Uq9WYhO4oDshi+D3ThCP9gFtDMSyi5rX8zkRMDnvJZmLzioQx2vuC9ehuzric4y59zrtbZ21AjlUWaMPSUJJnG2AMm0Bg0nJh/Y9NNPl9SEUcq3fmPcLcX1YNBwS0CXebKPxb+0A2Jpt6fpfmGfl9yg45VGPO745LF9FiqdZzjGTo7C60SLsVKjaTyybMPMdC4EKTqOk/FdHbC4WtVXiWiety9ZsnFzUffNcTMrZ7iTkdZgjKUlhA0WjPmb48lO/o6Uj7hDCgQCgUAuEBekQCAQCOQCcUEKBAKBQC6QTw2p6SGNC46kLZQaRuxjGneQ48RIhr0GFQ4Wkfo5GCwzpobBmRn+8VkDkXd3HB58797cniOyj9ZwOPbpgmRZeZWKYOPbao3v5yzRkHrehfd6v4ZPOF6fTDU/9aPJ67uxzCyBfWkz6melAp77TbjPEhoSy743DlftqurQRE+rQrUsy9gbnf5RDdukMWP/JPGMc7M0JBaYYxT8eldPP0nrjh8cda7EZ4x4tLgeMkZrqllynaVecoQE7YB8aTc1o1/Y5yR+zM4orocN1IMYPFD/9mBbWlz32KEiUocN7+lBZDlI8VzDuJltXZL/1ze2VV2R54i3BluK8SVsT2C8oDHRlN6biYNQ1yeNWR6fsgbzOhHOl5LWWlm/I+Uj7pACgUAgkAvEBSkQCAQCuUBckAKBQCCQC+RTQ7IFiL15Bw85azhxKdOPrLfPv6U+8+4e5sxUUyr12HIGIvP9eKWFA5FrEDv7n354ryfgoSot2MDRfymuZey4mX3KfiNxzwlON1JHHJsJNyOa4HiQne4Cu33noG/P2Ucl1dQABUpet1QfkhsSUKMZji/58EbtwfIfB5VOamSe1a8G5183Vjt72tUnx9w0EO9tIYca8Lvk4gbVkObPP0biZ4YkJ8KQk1R3oM5V446ZO1pAvAHvz9sB+T4jM9WMzMx+N//SZN1Vxxsc2nORxKIhtVUNqXMP2Ir1aNlCKm0bplZhm1yD4wY0ImVpSCuXqoZkC6HVoq9KYuYKiM3vBb9ZWYM6CJ657FF0o0Q60wupfMQdUiAQCARygbggBQKBQCAXyCllxyJfT2uxpDrLpzePyKIYzdKl6v6WmftSjbjGPS1yA/FQz86Aoqs6W2/pT26ldboft0eK63M3qs1NBWg5eyBZzgVVQ3sg/0nWIFdfiX8Yg+camJB6DXQgb8DfvusLk7Om4ZrJHndnRou3K2hH82ayLHWWyqvilGhrWjfduTKhm9Z0J2XHc4KUnXvv9C9q0PD5IUn5/CUjfiu5bsPUH+tQP8kVT8sd5qRXX0JOOyCWdgtN9xMts3+1tU6ifbWfizklgBZktBXzIKtOazDP/pWyAvNxlrO/2W5oOH/esnSbD/bnZqnSbX/2kfLlRGjEnjo/u8TL7AHiDikQCAQCuUBckAKBQCCQC8QFKRAIBAK5QE41pPMRZ1nvZM1+JHP/QZZcZ5WiU4fwZaJZGpFZWhdy5dsdURaKqlEb4NacNqkVvdZt1N+L62MxBfIEe0bik001pNpZjgj/vSlQDjzbUd00taG24N/5KOTaQDOyT2roy4NffAV/3YC/FZWD5wj335Xx9mIGs0SoIaEKPAtyVsN75wDTkuW+Ti9Y0+8gffA0agA8F/2uYyzvHC3bneNGlM7DCXXiEepVU+tO0yr4ImUNSTEzlUNQzkw7IF/aTc3IMJlWvt6lNCOOlGjpeczSn6vXlDhehqNppFWg1M7QhMnH5diG8TeGLSG+qQJ7OhC/OWwJcbrRoaOoCuO59gBxhxQIBAKBXCAuSIFAIBDIBeKCFAgEAoFcIJ8a0ljw+GtcTI6W8dYW1mZpLtiaSz0gA+wncuBodM9Pp0awI2a/hI+zNCMzazc00QcGVyr3PhTk/HDnSX+sPS+5wxe/oU+McQHmRlJvn6GpGRgF7dUaqn2pkRJu3Zci0gUaTsUDHrGPJ8GD+FvarKQbPBzAr2fsfw/faGSWFsXcXpSjThr6dPaDXtDZCxPojUprYFkaEjSL11VDWr784OJ6UR8dhXDiAMw7cKO6q0toSNymt5xE1gntNd7+h7H0GZntpl/IqZZbefZRV8kav12ONRhfJ0vjpg5USvP24O9PhvZstcihl8j/jtTjoado2O2sv0vs9eUxpmNSzP5g5SLukAKBQCCQC8QFKRAIBAK5QC4pu8oJ6r+xeaW7HV2DMsQsCo8ll6TwmvBcTRk0HHeqXQtrs3RJadcW1mYpiq5jL7VO6d1hRXFdA6rpYHiN9Le/FdeHmTofD7WXJK6a7WpXOW2S9dkvaPjyCrfGQzlv0pMOJBUwl9aG+eGnYzW3+Ex1Gb7XLpZ4/jOurj3FFMxF7EkkUjWgN5zlUseh+tn0o9eLshnW7MrAS5HBWV/E1ruUw9uvlaN6UpY37RFnUVGgjGh7U0i+D6v7qGu1odrcDkiWpGL5rSKJ1eickTpphXtq0qu4dpeyAxKajhZKpMfK+SnM+jTLImeBLPd+xizvx3nrrcPYAkKn/PpkWXmKngT1bZ+VOKslpHY2fJI4zHgPEHdIgUAgEMgF4oIUCAQCgVwgLkiBQCAQyAVyqSHd2vZLEjcenPChaw/Wstx/gET30xs32Ycl9zY4Wk5+3GH7tHhMrVGLu4+9U1y3x+TQrDLdLmKDZFYNe5AesInxtjF9oSHV7NC4w4L3kgDTWFOlz55SX6ipReDxOU7A60TUA6gX1Lj1CORq+Q+XJMtV45U/n2DjJZ606hKJbaJbz9PPI60feJ6/hN1+4p5jgzro8xxir+lj4cSz1tU3U7HInGXcQXPbW+l5KudtOZ0KKWCfMsYqbML3bCeqptu48SAsqOaPDA9ZziHYL3XY8J7EMum15AgJfyTllFSXQtbPJq3BqAv5uAtytPShSOZ0IbaAZI2XgWbUfYS2dYxwbR/H2XOSqzfVkIYsRl2+bwlRecnsUSsbcYcUCAQCgVwgLkiBQCAQyAXighQIBAKBXCCXGtKnZ9yv/+Co02bQrOur1Dd+i/Pp2QI+l/E71hZxoiG9a60k94FpSNvUOKUdbFYgIWlvywrk4PAjGgb0jA147Ovu7RTwNLR6yRqAzE4JmpSMcFvcaTSSF+J1xybPdqddLrkJO1RDsjvQ/CV2QfAzSr0jz+MfqinYMXm7/WGm47QHbwSfDrGt4NZUMDI7TrCpPG+3+78uY8RFSWSMWaC++nZHVcHatEs0mXJ/VGRvYD3Fk2+/Hhk9WKkREn6XSx0VP3jfNUfFL2vcDBU0fkOc/kcNrAZxGeNlPnSMbtzQ6qTvcDj68OrQaOitwwYsRichv0qMXY/iIvzmqNnUniHukAKBQCCQC8QFKRAIBAK5QFyQAoFAIJAL5FJDsusRO92oolJT3bpo80S3Shd3gDkX+jsgIaknXSvkMBJAuHvy3ox9bwW99zcgpqbk8jsx7WAleja8UsLhx3xaT81T3yDIivvOnUHIHUYR6Xi3Pk9Trxyvhmh322XJetdlktv8E8xZmIjXafK8ODun2HPiffFguAW7/T6jkoauj6JHo4IegHhZ7+zFI6Cy0Muf17CNa0R/ygYvPtDLkb1FmY1K0Eboyehiaqit3uUX4p8dcFps7Zp3KC9Rq/Jxlo+lmbYW1SAHzajjQPVOPKRD0vd2GE62QXCVPMI1HnL0TLfZef3jHwAAGJFJREFUaDjz5zF9LNG/uCRDNi3gT0NDCgQCgcA/LeKCFAgEAoFcIJeU3dO4bfSUEY05SCd1cjTcvqAg2pCiY9y6hbVZmvnw8Q5NbQdlt8XFpeZJpqZp+ucp8bc+X8oYxZM1pI9KGNvbYE8p0WKepd1ukOvUHjrl9QE7X+L7NiQjJXbejk92Ap53DQdf+JOGO8PBF+6gj0HqXA1PsGnFdb01aHKWhhtA2ZEm9aBJTJveLtDBrVbAP6zY5R7MkRGpMyiLkEXxeQYVtS+owPbb1NLHU9jluhkJcUiqHPEu/w9Z30kzK88eCPSlp+nQcZA6Z7onr9Oxq8686dKhZauw3ujj4HgZxv3NU3Y6XqbnEnzunmrTboW0jZizDqNtGAaspGg5nm3vF3GHFAgEAoFcIC5IgUAgEMgF4oIUCAQCgVwglxoSqHlhd7NMO8zM2jg9pzW0nVIGIBkDzFPw7DSpazLXnsUv9dhS+Sxk6UKpMmO3pkNJL2gYqRHInkM/SVOLh+qo8cdtTHH9iBeUzGzm3BP1j71O9ABecxM1I54lns2msnikht2dKnaBpoaM0HntZ9hjxXW3p1Aui0NYgEpof0Q8T2sQi1CHWvrX7BB93sVO5SvwiVjwzzMo4yxhiXJXv1QtpAItB761gaoVz+nM7zB1LI7i8NY7HJeRslF6t4X17o4CsT8OaEan1j0ssbfi2R/7Xw0l0ec5aqbvRhUEK2gN5kfKsLOBsSvPfoul2vhdLLh1ObZhZvrZUXveG8QdUiAQCARygbggBQKBQCAXiAtSIBAIBHKBXGpIVACy9JoPug7+/zcyR1mbsvyZIwtMu23I5/aCHY20tnD88RGIMWp8ycBEgXrOjpPcND+vwcwe35FoSJsfgP0PdaKpPqAvz18R85P3uwPNqDWapZxu1HG82rOcq3Ms7MyNTyXB4/o0G1RuMlD1ct6y72gwLLC8TvdGX92nBRSVfB8JG0VS3U9UdPy3C/1Z+Hja1SR7TG0kNSbFfRzUHahiZQ315pd/Wxf9P7OMUteWn7SmVNIUK+OovIbUXd8BxzeMd+Jnz3k4L/+uoXn7r6xRM7uLXdvSW8gVMnQhqorvo1Mts0fxaP4g7QXiDikQCAQCuUBckAKBQCCQC+SSsjsDcZZ9Dm83t2fk3k9JNZFVIp5FSTBHepLWSNXOKaVLFyTJAx3o1izdpv+PY4F2DtPUvErl8GbDH+h5O7a4JmW38km80JQW1mZmBVqd+zpqclEkgrhzntaq19RYPPSq5uLyog73Sep8myxxxe9d8KQ+zYso8yZZ5s+RQ5kjTerMpbnf88iheqpwfbPmUgQN4XkVEDCo/+9bWUjWrC8nneTefCkand+Bam/hhXO8sS2dzt0DaJuUcj735wy/7SVMyBwzSDsg0pdC092Jp6VNj6Pp3sIJ0wjajZ+kj7nH5dBw/O3yZwRNttgSciRouQrfEjLS3jfiDikQCAQCuUBckAKBQCCQC8QFKRAIBAK5QC41pFravXupgZYllBZ86Scnt4KjbYao5Idg7kSuDXaqtdN2Kjhpk2MtvP0Jy30ZUxfyMUu3D0Ts5Jtttfp/jdfa9pd4kZvnOA/eQC8hnrtjuMSbp7n64AYcwzTEwqHPRXIBYm9cUmrGKkq7PYE9DqlrNDzv4N8U15fYbyRX++hKfbBziXl5saayyrzN1B7oSE4rRin9xuOTk+g5+6jkls/GyFIZJ8A9LKXgOJWga4WmoGsdbEt3uzYzs2UabvfaSIkjSKk3XjfCOb0aJ73E+KhSZeCZR1KiicLZJqVHSED88XoaNKOXMUrHO/yU0n0Yl6N5e0Ws1HgZrxP17Y0kbcOGI3YTZdaNUt+nblkH2ALiDikQCAQCuUBckAKBQCCQC8QFKRAIBAK5QC41pKbbNW7nKVu2rpCI9RoTNSTEFbCr96MrUn1G3Ck/WrmEZb4QuuwlQtwMDWl1VcIAk09fYUr4Fpxq8To6CDjCwGtIa+YfpC8KS5xU7DWMhcgZx0R41ryUub3f9Rrkjtaw82Ea+zESV2nqvMPvkfhyu6u4HjkLlkSwM2p0zkGlzIuyhl60h2bkpnKYmdkzznKpgX1UUzVUDYlKVql+LXdesBcK+sAgp0/VroVgg3EHq9x3i72CBL8CclpDw+A5LjE1JIxZL88Up72Gbj5LNTqCUjZKGXZAnArhu+tKaUI8Qq+9ZdmGmalOVAudug17Ev15QI0Irlurhusr/9n1IVL7/IWVj7hDCgQCgUAuEBekQCAQCOQCuaTsvtLhZom79k14us59tbazs7Uc7wfyoD1u6fexdyRu5cyFWpv6wrwrHJ3ZO662e4ftI7mtKGzd4mJxKzaz9RjTuRZ1343uZrwkZddYU1y/txi8Iam1hXuYMzPbRCrEExG0+Cln5iSJB88lgOMaAOoJk14/dEXCGV1SraXcl9lEiUfNfjEJJkjKtj+qcYNbFzSVonVpDzTYvx1QdG+MUmvtJ+3k4vrV2ShpZyn9u56m434TGf7MoGcOOvwViY+wl5KAFjhgCgtuTaKslGO02FxhE5eCel651MUpp3PuRRZ5COKQ03JrkmVvb7Nt6UmvWXZApHU9TVeqHLtXRr4X20OyrMJ4YmaUcv91f33wXFNfsVnONsxMbcSWzB0iuV/AkmxPEHdIgUAgEMgF4oIUCAQCgVwgLkiBQCAQyAVyqSHdMeM6/QfH77bprprEfp2VJ+7cKtGQqBkx3hds995qSNvBkr+NElLRkLaphrR1DcjrNbBz8Rb7BU2lyl59nvw6/3arZ7NZOry8xAt5ZrzUVE6v/fRBbrCG7RzxfQoeeq6GvS7WY/ZjIy4yHSlx5HQU3979/9o7lyCrqisML6P4wgeRIAQkNISOBnygUIASpVVEfKFWQKOJGisOHKQq81RmGWSSYQYZWNHy/QzxhU8QFKxgaFFE09oFdgyYVoKlRFGCSAap6vuv79j7cjVVblP/N9qr9+nb5+67T++661+P1nDXsjz1GFIDGjKFwNL88xnuv1DGS/PUw7E42Q+pnZvWNpvnptD6Uufcz7lLDfHtyVPshDpT48vX52v3YcuUij5R7uhiaS352LdNzsqKpidERMQmeT7aaki6N7kuuKsuTE/VqYE0Nfk99L2Q0kFsIVF6OqgZXYD9czhCrpMWxNBt6kTyOfePz2oUOxBrqbBeCIvrYW9fgdpOq2TMffpMdIy/IRljjKkCH0jGGGOqwAeSMcaYKqhSQ4obYYvMsmdUzkd5r2GLZ5YlfWizbUQnq/HJMOOI3AIjIpfFZ4l82mzDPFiY+5Teei1pQn96qSEy8zXaNRDQDJx22RTq7IZmxDIlF8sYmtH86bl+zuWRxZ8lIrxMuA+6yl3Z3Lm8NX4KPn90mEh6SBfmemAfjlyjuKo1fHD8wjR1P97g9jvEN89W741yTNT4FK4/kkGkS8cxi/IeOTOeTfaUTbL5oCG9giUuNb1gWZujqH9IXgz1jYaGpPlQA/xL3OOlrB9k+UBqO+LE7UNjakgHbMFLqYaEqVJOFvOMGprRNbBlf718bF5ErtMb0Wo3wzXdCLt/szyXa6BhtysjlvLTqEXzg26PvyEZY4ypAh9IxhhjqsAHkjHGmCqoU0PqY0C76kSNBsiwDxtmHBEB/yjpZDU+LU3ug60aDfUZ2iU9h556/q56rPk6nTRAbtPeOXm/0V6bRbVUJ+rBpcg1mnJuq5bahfFomrsolid70dur8y+rTvRwntqKS1fJmGpMqQ35fMyNpWb0k2z+ZV4rGeSeuDLNrVl3Xr74dhn38XNl4wud5/6HTse6ZaLTnX1gThQ5K57L1+q6QUMqNb0oNLz4LyfBlj2yIU5NU6+8g4uTZlHSRUmbhg3Yxt8b+frQ+Ph4PU9SaJRSdyUVKyL/92rU9GMuEfbXH45tiZJa+zCimS+0ZbNoSuvxf491CV8qzA3yHbChRqmuojUkY4wxX1N8IBljjKmCOl12AR9Luk0W/edbUHdTu2sxX3TDEf0qy1/k19xPC3MMDG33WvsL3zsdKerCKPWbjGg4XUZJaSSGbs+F3dMaHnduP6ZWJXuB9Fk4P55Ic+OWf5BsePAipLNrL/xJdHiVGjZ0wVY33eSFmISLrm9xLo10c1w/NL5r29X5YrS9yF1hV2GSd6yfbaPpRTbhFj3uotbicI2nb0I8szyG/bkDQ3ENuXumjccPEIneN7m1bmx38Nka1NMphhnTTa3rhLuiHxGuzWnimpoWr+VJeK12Sth3Kfw9Ij91jRYS8HAxtFvddPe+el2+mOWmNDyb7WQanXbflHG7djLtAtu/HP6GZIwxpgp8IBljjKkCH0jGGGOqoFINibelusr/1mf51dOJJsZw7FIIfLuSPl2t4SiEhTKS+0TYqhv9IE91T3852WfE80NjhhWfCbt7vTi3n4zMimzuXJvttVICiIGpDKLWFaeU0AM7hXb/NM/1L83FX26KG5J920dS++V3qFN1N/5Qqt1PP34pGH12nmLbDrR6v1DEt4XQkBq6nOgQLF5ErUR3XiPMm6Hn2DPaBntdoH4OM0AGhjWiuU6qjUKgwZ7+xtzcd+QkeccT+vFuscEGZO+105CSWsu245ADWQ4ohXZTM6ImuXWHGHwiBmDrXXeaLqKffCOQvWP8DckYY0wV+EAyxhhTBT6QjDHGVEGlGtIFsFUFoI9zF+xSzk8n+UKktFRfJjeK+UGd6EKjs6ntNVjbno59takRwR4z561kq399VvSmuVmoMaMa0oT18LAz3Uzt3E07et/N9iv41VILba6wuuoXHJ3nRrAckGgwzDOiZnTT7mx/+NsxLeMWvO6HvfiBZktR9eLnLrrRCdgD0IxOP2Vlsi+Jh4bGk1ZuzxdDp9souUelhhcRWT2YyfyaedncNie/H9WQtq6G1tNo3645Qe9yEsiNHIRcO+TLzRi7IdmnasITy+kg/WlAxlS4ufeSysLKOnjutIVEBMoBsQ1E0owiQnL6mppRSYdv0+q9kbOo/0jQKuQL4G9IxhhjqsAHkjHGmCqo02U3F+VPSh1X2Z2VdoJVuL9MmZ4Rw4w/B13lUZhjF9tvFexxmKNbrmuYcUTECbmt7ZQJm4fGUxFmzFIpJ8FBpvaMD3JdkhFwtSXXAipGw9sXG99ujRmoygImpVBuBp+eBnumhtsWurxGRPx53ilD41vj2jSXwroDLrqIHIo7iPIz8QJsdWfS0YN3MEpcI3DRjbsul/9Zgvjgi98VF95DaSp2IZReQ7253iU3KCO349xsPhNnD28/na9t7Jm0M9rdlfjE2oSe0/U8Q/10uIedcNl1UmO8W13ECPPuH58faHZ6TVW7WQ6oGNrdiSOxC3NspQsXsa4ry4Z9AfwNyRhjTBX4QDLGGFMFPpCMMcZUQZ0aEstgDMqYGlLJpp70CUrkfML2E2JTXuJKqY2qMA1daNQw48+zqRONa2k/YybkMNfx8fawdhdCPbtSifmIqdHSkL4LDen4j7KT/NBSh0nGX8O3vUfsjeggUSp0T3WANHzzMp59YJ4bDb0g6UbQjB6ZeE6y74wfD40bLSRYDugW/J2kG0GgaSgPGm7Llqp4A0tkfEPWBn8U9yT7h6wxo511UTloba6ek0K9qZIyi+Bk1eXQWnfLnLypn4aotGX1dJ0EbDFRanzBEGXRorGE4+ZlrW0O8gzGrJN/Htj/f92dbVX/uE7UM0foRoWuRc1oI/eB3kejhcQAbNWN2miS+oFNxf9ItpehLes6aQ5b6bIGWXv8DckYY0wV+EAyxhhTBT6QjDHGVEGVGtKvpv8y2Tumt5Jx3ofo8i+U2lH7Y5TB2BWHJ3t3HJzsvbIceyMLEQfG3mQfEi1H8sHx7zR3OMoZHSnljo5E6aNREMFGxw7Y/xwaj49/pLmx0CEmRqvWy8Tduef0yDc/S3a8IWO66Tuwt6K1NXWhARl30vy4WHIlmirLNPXNQ8OIS7I5uLiVDLIsLktz9yCxZ3Wv9HOgtnk77EY5IM01KmlGEbn1+Hl5akk248bW8KoJd6apK6EhTVqG8kCSe7QRLn9msmhWHtefqUZpzRfnKW29/Xl2PCJjlsQpNr7gGiKxRyWMnjw1T0paReQSVxGRSxZBFx3AX9V9TG2z0YpDywNBj9kAUal/M/Ixk5aVNeFy4wt+enhAFoluxL22KGuUsydkrW22aG9zGrl17JHRHn9DMsYYUwU+kIwxxlRBlS67X//9N8n+RL4Hvz8yl2cuuex24Sv9xxW47I6Ay+6be7PL7qgdKF+kkd45yrtsc+4t2Opqw7VvoogyA23VpiOqk36TJbccnC9x8nj8oFSeBi6jlRNPT/bD4sNbFpenub/dg1BVdcs9EoClqF+EXXIvwR2jrV4vw9TPs3nBzD8Oja+NW9Pc3LW5Y2/cl813pHJQ6W4j8ufDux2LCt4aSr+mO4cVP4wPZPCBKfl3H1eDLrpSnfE2RaJ6WsOjFw2mKXYrPqEPf0c8U6/h+SgFnrMW9uSJ+IG46V48Nu/yXvrw1iAEO7ns6Bznk6f7rStPMbRb3HTdP8v75/L4U7LPjmeSrS67Y9Zm9x6rvO8P/oZkjDGmCnwgGWOMqQIfSMYYY6qgSg0pfpHNQ6Xi+bijc/2ZcSNRj0Yd3yzpMxI23/0hhXtiKSGVlOA6jY8KNudw+w1H/o79nIuIPTK/Fa/L3pr6UtSBWLanE12IpVPUp06PP3WibtWJZmKS5X8WZrNvRqubK8OKl6PHxJObRdO4G/50Rqq+pHogW9wyULrUDoGB6gjtVt0ImtE587N4dX3cPDRetAn3dFs29yzP9ioZD+COSuWB5lALQduO95a2HjZqRo9tw8VZlojYpDuKGlJpTaFsnYjPUmS5nkOy9tEDLaTx0YpeQ7WGd6RqTSPMm20vRPvsxSZfTw2JIfCDuk5UsvhUahkl3BXL/0hoNzWjG5Dr0L0SNYu0yzBbhaA01f7gb0jGGGOqwAeSMcaYKvCBZIwxpgqq1JB6l2X7qGHGERGHQfc5THSjEdSEaFNjUtDCAGlIWVOivoTy9LtEN/oYetNOvC71mp37OUebZXna2SWoLWh5FOZdNHQiWcfR3Zg8FbbmFiGHYcus3MJgFdpgPx0LhsZP7M0a0nv3467UTc7cog9ZJ0lLpTAnhn57Fo7RvJhCC4mIpBtRM6Iff+mbMn8LXgfPzmNtWn4ok2D36PPCVu9Ls/mAvCGWY4rb8aA18rm0NUcpyyci3yU0pAXZHHPpWzK1Is2d0ofPGRpSv0wPtLkj3V2ncQtAr9k2q3XB2jgjzW1f8Z18caOMkmqWpcbpEenJZNtx3JOWA2KeUUMzYvksWdZeCNWUgfcHf0MyxhhTBT6QjDHGVIEPJGOMMVVQpYbELsaqYfCGR0CvOUhsah+N3y3cA6+lTKRQSeC1ewpz1HJKr1XK/2kH36tqcUdijm7w42B3yXgs81OYXKQ28x9Qj66vu6UPPBdnprln46xkr0I/ga0rRKCiRvE47D7NLWJrceYWdVLWH2/oCNE4cleL1EIiIten0zyjCGhGERG/l/G9eeop+PE7aSnB0mNHqW50Bf5Od9bE7hcNqf/BUwKTmffZOkHvkk8EVUrR5ZibBunqfEmEOZ9JMcyRgV6jd8QdwKqEKo0ewLwjdHp4Vvb189CQUpJYRKN1eu4Dw3Wiui65R7wnrJvWo5uNVu6Q3hr2Y7LfuNesIRljjPna4gPJGGNMFVTpsmvn1vp/puRmpKuAtrre+AWebjh1hHTljh4x4tu4GN0Cko+CLjq4B3bOar2jDQfmOO91MTvZL4jLiyGxgytwE/Trqs0SJsFOrupcaBfKXSoMA99HF9xL6qa7Mb/uFZNyp9drpI3ExX0r01wj1FbcdE+hY+8ruLTkAKOLbjI77Upo98vn5Jj9O+LqZD/56qUtg51013NN2VlUQ5jbdIEdNbk1hovu+/NzQ40L49Ghcfe6QsmbiHgNXsSBGB66Omfrw4UtsX3+Ecl+TlzP/b1wbbKbSaNVcykknnc1rTWcm2cmzcmtgrXTa6OFBJ4lhnbrk8Q0lC+CvyEZY4ypAh9IxhhjqsAHkjHGmCo4YN++ffu+6pswxhhj/A3JGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfhAMsYYUwU+kIwxxlSBDyRjjDFV4APJGGNMFfwHbLZ943NG6LAAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x319244610>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "include(\"NtSolutions/fastmarching_0_implementing/exo3.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the geodesic distance map to $\\Ss$\n",
    "has been computed, the geodesic curve between any point $x_1$ and $\\Ss$\n",
    "extracted through gradient descent\n",
    "$$ \\ga'(t) = - \\eta_t \\nabla D(\\ga(t)) \\qandq \\ga(0)=x_1 $$\n",
    "where $\\eta_t>0$ controls the parameterization speed of the resulting\n",
    "curve.\n",
    "\n",
    "\n",
    "To obtain unit speed parameterization, one can use $\\eta_t =\n",
    "\\norm{\\nabla D(\\ga(t))}^{-1}$ (one need to be careful when\n",
    "$\\ga$ approaches $\\Ss$ since $D$ is not smooth on $\\Ss$).\n",
    "\n",
    "\n",
    "Compute the gradient $G_0(x) = \\nabla D(x) \\in \\RR^2$ of the distance map.\n",
    "Use centered differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "G0 = grad(D);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalize the gradient to obtained $G(x) = G_0(x)/\\norm{G_0(x)}$, in order to have unit speed geodesic curve (parameterized\n",
    "by arc length)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = sqrt(sum(G0.^2, 3))\n",
    "U = zeros(n,n,2)\n",
    "U[:,:,1] = d\n",
    "U[:,:,2] = d\n",
    "G = G0 ./ U;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geodesic is then numerically computed using a discretized gradient\n",
    "descent, which defines a discret curve $ (\\ga_k)_k $ using\n",
    "$$ \\ga_{k+1} = \\ga_k - \\tau G(\\ga_k) $$\n",
    "where $\\ga_k \\in \\RR^2$ is an approximation of $\\ga(t)$ at time\n",
    "$t=k\\tau$, and the step size $\\tau>0$ should be small enough.\n",
    "\n",
    "\n",
    "Step size $\\tau$ for the gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = .8;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the path with the ending point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = Int(round(.9*n)) + im*Int(round(.88*n))\n",
    "gamma = [x1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a shortcut to interpolate $G$ at a 2-D points.\n",
    "_Warning:_ the |interp2| switches the role of the axis ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "Geval = (G,x) -> bilinear_interpolate(G[:,:,1], imag(x), real(x) ) + im * bilinear_interpolate(G[:,:,2],imag(x), real(x));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the gradient at the last point in the path, using interpolation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the descent and add the new point to the path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Complex{Float64},2}:\n",
       " 0.32013+0.947373im"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = Geval(G, gamma[end] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "append!(float(gamma), gamma[end] - tau*g);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 4__\n",
    "\n",
    "Perform the full geodesic path extraction by iterating the gradient\n",
    "descent. You must be very careful when the path become close to\n",
    "$x_0$, because the distance function is not differentiable at this\n",
    "point. You must stop the iteration when the path is close to $x_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"NtSolutions/fastmarching_0_implementing/exo4.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the geodesic curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3U+sLEd1x/Ff99z3nv8RFLOLiYT8R8kOBVtxFCBhxSLZhUUihYBhhSFYxsiyTXaRsEmIDDY2LP0gQkAisssGFokBC5B4TmCVYGMlkh8r4gThP+/Pne4sZubdmrozp7q6q7ure74fKTJ9u3um7r3z7kmfU3WqqOu6FgAAIyvHHgAAABIBCQCQCQISACALBCQAQBYISACALBCQAABZICABALJAQAIAZIGABADIAgEJAJAFAhIAIAsEJABAFo7GHsAuRVGMPQQAQAdt+nbzhAQAyAIBCQCQBQISACALBCQAQBYISACALBCQAABZICABALJAQAIAZIGABADIAgEJAJAFAhIAIAsEJABAFghIAIAsEJAAAFkgIAEAskBAAgBkgYAEAMgCAQkAkAUCEgAgCwQkAEAWCEgAgCwQkAAAWSAgAQCyQEACAGSBgAQAyAIBCQCQBQISACALBCQAQBYISACALByNPYApeYukmyS9Kul/Rh4LAMwNT0gBb5Z0n6QXJP1C0n+t//vC+utvHm1kADAvRV3X9diD8BVFMfYQJEnvlfRNSTesj93oXa3/+7qk90n61oDjAoDctQktBKQ93ivpnyUVkhbGdUtJtaQ/FkEJADYISIm8WdLLkq6XHYw2lpLekPRWSb/scVwAMBVtQkvGNaS7JN2//u+wPqhVmq5JMNL6uhskfaC3EQHA/GX6hPSMpHu0SpjVks5L+vBg7/+CpFsVF60rSS9JuqOXEQHAtMwoZVdrFYw2akm/K+lHvb/3W7SaRdfl/lcSjQUApmpGKTu/hlRIeucg73xTx/vflGQUAHB4Mg1IfmStJT03yDu/2vH+XyUZBQAcnkwD0nmdBKVNDan/dJ206sDwok7WGTVVre8jXQcA7WQakD6sVc3o/vV/h5vQIElfaHnfk0lHAQCHJdNJDaxDAoApm9GkhnH9Uqt2QLVWwcay6dTwJyIYAUAXBKQ9vqVVO6A3dFaVClXezL9q/X9vSPojSd8efIQAMC8EJMO3dIveqp/rfn1eL+nWrXMvaaH7Jd0ighEApEANyfQeSf+y/t+1btYrepN+pV/pGb2ivx5xXACQtzahhQ36TC9oVSVaSCr0it6iV3SzpNfGHRYAzBApO9NFSQ9pe6FuIekzWiXrAACpEJCCLuh0K6MjSbePMBYAmC8CUtALOt23YalVXwYAQCoEJABAFghIQXfo9I9pIVJ2AJAWASmIlB0ADIGABADIAgEpiJQdAAyBgBT0qnZvGMjiWABIiYAUdJN2b6l+4whjAYD5IiAFMakBAIZAQAIAZIGAFMSkBgAYAgEpaNPx21VLumuEsQDAfLEfUiOflPRZbU9uOJb0Nq06gmOu8vss2jL854wD1eazyBNSI3T8BoC+EZAaYaYdAPSNgAQAyAJbmDdizbSjhpS7qdWBuujyvVJ/wth4QmqEmXYA0DcCUiMXJT2k7Z52haTPSLpllBEBwNyQsmvMmmlH2m5oY6Thckz9pUyzxXx/pPfQB56QGiNtBwB9IiA1RtoOAPpEQIrCAlkA6As1pCibBbJuHGeBbB/6rNfkWH9qW5MZa5q39b7Ul9AWT0gAgCwQkKLsWyB73whjAYB5odt3lFsk/bdWQchF5+/USNn1r6/3zPBPCkbQ5nNADSnKRUmPS3rQ+zrrkdpIGRjavlYugc+9ts8/6O5r97XuKMfgi2kgZRftCbEeCQDSIyBFYz0SAPSBlF0rtBFqo0t6rG36K+W1Ke/t8rpDpbxi0ntdUnqk8LDBE1IrtBECgNQISK2QtgOA1AhIrdFGCABSoobUGm2EdsmhTpSyhjRGjSlUU8mhbY//Pl3GlMP3gzzwhAQAyAIBqbV9bYRI2QFAGwSk1phpBwAp0cuuk09K+qy2JzccVl+7vmosXWo7fd3bZ73JYv0TDf3zHePeLq8bkuGfK+zR5nfFE1InzLQDgFQISJ2QtgOAVEjZdXZ4absh2vjEpsdiUm1d3qfpuT6lSpf518a8bqprQ/emuhbDI2U3CtJ2AJACAakz0nYAkAIBqTP62gFACrQOSmLe21EM1dLHqu2kvLfttdZ7tjm/T8rp2NZxzLW+lFtIWNtcDDUG5IEnpCRe1fYTktbHr40wFgCYJgJSEjfp9BNSIenGEcYCANNEQEqCiQ0A0BXrkJKZz3qkvrZv6KsO5B+XZdn42pT1Jt8QNaTYOlDMvVVV9fK6Pute69qQDP+0HRTWIY2K9UgA0AUBKRnSdgDQBQEpGdYjAUAXrENKaprrkWJrHan6xqVcH+TWjaxz/vm+6k27jpvqcy2RVRdyz0nb4/fP9bW+yWdtlR6zNTr1pGngCSkp0nYA0BYBKSnSdgDQFim75KaRtuurHVCqlFcoleafd49j7rVeJ3ZMXdoMtW0HFEql+efdY/+c//255/3xW+/rn0uJNkPzxhNSci9I8v9BLiW9OMJYAGA6CEgAgCwQkJK7Q6d/rAuxQBYAbNSQktvMtFs4X9vMtHt2lBFJw7UDiqkhxdRrYo4Xi0Xja2NeN2b8u46bSjl1O6aG5B8vlyczRkPfm1U3Co2xLWtK+K7zba/FcHhCSo6ZdgDQBgGpF/S1A4BYBKResEAWAGKx/URv8tqOIuVW3KnWFnWp7Vh1Iv+cda3/ukdHR42vzaGGFKoDHR8fN77XrRn5x6FrY8ZkHfe1zUXKbS3QDNtPZGVf2u6+EcYCAPkjIPVmV9pOkh4QkxsA4DSmfffmoqTHJT3ofX2YNkJd0p4p2wG1nWIdSrtZx6Fp325azk/RWa8bOxU9lZip234qzf/+3BSen86zvh//df3fs3++rZRth2JSRkwDzwNPSL16QkxuAIBmCEi9Yk0SADRFQOoda5IAoAlqSL3bdP92Y38/3b9Ttgdqem+XbSGsqdx+LSdU67HqQtZxlxqSf6/Pf62mQvUYt/YTM3Xbv9evIVnHqaa0x4ppSWSJrRGx2+w4eEICAGSBgNQ7un8DQBMEpN7RRggAmqB10CCGaSPUtj1QzLojKW5LBqsGY9WBYmtIZ8+e3XvtmTNn9h7H1JtiWhJJdtskX8yWEm5dyKoRhY6vXr26dc4/dq+9cuXK3jH414bGZLUdss5Jcdu3W22GfF3aDmE3Wgdli5l2ABBCQBoEaTsACCFlN5j0absu07zbduyW0k3lttJyfpotJg3npu9C11rn/GOrE/iu89bP2Gell6yO3f65UBrOPQ5d66bpYtJ71jnJnpoeStlZbZO6dAaP+TOY4Z/MLJGyyxppOwCwEJAGQ9oOACwEpMHQ1w4ALLQOGpS1ad9DSd9prO0nYraUiGnpE6r1nDt3rtE5/7hLvSlUU3K/35gaUqiOEjN126oLXb58eetcXzvexlwfqu20rcuFsP1EHnhCGhSb9gHAPgSkQW027fMxuQEACEiDY9M+ANiFGtLgNpMb3DVJm8kNX5e1JqmvulDoda17U65Dco9DNaPrrrtu69itBVk1I/9ev4bkX+ueD9WQrHVIXVoHxbT/8Vv8+HUid4wxW7LHfPZi1gf5xzHXhuo+MWOm3pQHnpBGwZokAPARkEax2bTP1c+mfQAwFQQkAEAWqCGNwtq0r/12FG1rTH32sotZh2TVa6zajn/ery9Zx6FrrfVNobVSqdYhxfSNs2pGoTGlqrnEbvXg1pz8z09Mvckaf2ydhy3Mx8ET0ihoIwQAPgLSKGgjBAA+UnajsWbadd9FNjYNZ52L2QnVmsod0zoo1NLHmsp9ww03bJ27/vrr9x7HpPf8MfnHVnosZtp3aMdVd2q3P807Zmp6aMfbpuOV7J1cY45D077bpvdi09mk6cbBE9JoSNsBgIuANBrSdgDgIiCNigWyALBBDWlUm7Sdm8vfpO2elTTMdu5dpn2H6k1WHcWaNh2zhYS0Xevxa0Z+Tck9tupL/uta7Yoku14TM+3bahUkbdeNLl26tHUuVEOKmYpu1XasGlLM9hn+mPxrrc9XaNr3UFO3mSKeDk9IoyJtBwAbBKTRWZv2AcDhICCNjk37AECihpSBzaZ9D3pf30xu+PneO1NtRxFbQ4rZfsI9Dm3XEFNDstYLhWpI7vmbbrqp8bWhGpJVv4lZh+TXWKwtJkJj8NcaWfz6h1sLCtWQ3Gv9GpJ/7P/erbVF/rF7bczntMuWEWw3MRyekLLApn0AQEDKwv7JDXVN2g7AYSBll43urYRi0nDWub66fYdSdu6xn3oKte2xpmf7KTw3TXfjjTdunfOP3RReKGXnpxVTpez8Dt5uyi6UootJ61oti0LtjNy0on8u9Ht3r4/5PIWmiMd85i2k6IbDE1I22LQPwGEjIAEAskBAyoa1aR8AzB81pGyE2willKqmFMr5W9sdWNtR+LWRmGngoRqSWxcK1ZDcY6utkGS3EupSQ/K3mHDbBfn1mC4701oti0JT0d3fj38u5vfuj8na+qRLzYi6UJ54QsrGvpl2jzHTDsBBICBlhe7fAA4XASkrbNoH4HBRQ8rKJm33WZ08KW3Sdl9XUbTf2jzlOiQrj2/VlFLWkKw1QKFtyd1akN8qyKoh+df6rzvUOqSY17XaAfm1Hus4tFW6e2xtebHr2FrXFvPZy6GmRK2qG56QsvO8SNsBOEQEpOyQtgNwmEjZZaYoLqqud6ftpG90Stttv0/z1ikx0767tBlyUzlWJ/DQsdVWyD+OmSJudQKXTqfsrO/HZ7Xp8e+10nSh3VmtNJyfGrS6r1udza2UXOg4lAK2UnaWlKk00nL94QkpS7vTdnX98TEGAwCDICBlad+mfZ9gTRKA2SIgZWiVlnt8xxkmNwCYL2pI2XpSq23Mt1sJ1fVdKopwK6FUNaLQccy0b6utkLRdewhdax2Htqpwj2PqTX7NKFRDcussMdO+/enXXWpGfl3I+t6t1kehn781vT9UAxtj2rePulAeeELK1OopaVcroUdJ2wGYJQJS1liTBOBwEJCytntNUl2zJgnA/FBDytj+NUmP9tpKKObeLltV+KztzkM1Jat1jXWtX/exWhL514aO3ffx6yo+d+1Ryq0q/DFZW1fE/Eyt30fM71nqry6U4j4Miyek7JG2A3AYCEjZI20H4DAQkDLHbDsAh4Ia0iRYabuTOlJMu31Lyrb+1rV+XcWqJVjX+udD62Bi1tdYPdpi+uuFakgxa43cOlHM+P3jmG0iYn4foWt9Mf3p2m5hbr2OdHrdkXueNUnD4QlpEkjbAZg/AtIEkLYDcAhI2U3G/g7gZfnIGAOSlHY6rbUNQcy2BCm3xIiZim6lvEJToV1+ii5minVfWz3EtPSJ+V5DcpiuTVuh4fCENBl0AAcwbwSkiVil7T634wxrkgDMAwFpQoriSe2a3FBVd44xHABIioA0IaunpIfF5AYAc0RAmpiioJUQgHkiIE3O7jVJpO0ATB0BaWJI2wGYK9YhTVBRPK+63rUm6bZOW1K0kXJNhrv+xn9d/9hfq+Oet8755/1rrWN/qwf/WncLCSluPY57r/861jblMeP3j1P+TK1ru8hhzU8OYzgUPCFN0r5WQqTtAEwXAWmC9qXt6pq0HYDpImU3UfvSdsvlx7RYPNKp3UkoXWa9Vsy1fmrKvda/z7rWP+9f66fa3GPrnH/sv+7Vq1e3jmM6dvvc1/Zf1z92r40Zv3/sn/O/P/c45vcRutZn3RtzbZfPeOx59IMnpMmilRCAeSEgTRSthADMDQFpwmglBGBOqCFNWFFcVF0/LOlvddK9YbMm6R/NKeBdcuQxefyY6cAxU6z9826dJVRHca+9fPny1jl3N1b/fGj7CZ875phr/ZqRP0b3OGb8/mvH/Jxifh8xv2ep/ecpBjWhaeAJaeL2tRKq69vGGA4AtEZAmrx9a5LeMcZgAKA1AtLE7VuTVFWsSQIwLdSQZmDfmqSqulVl+fL6mv1bQcfm7a2cv9WqJlQHcterhK61jv06inXsn7t06dLWsbVduM8fs/vafv3JZ6138utAb7zxxt7x+scx37t/HLNey/3eQ+vAYj4jMZ+9PtclYRg8Ic0CrYQATB8BaQb2txJ6jCngACaDlN2E+GkFNw23O223UFU9J+leleWXzdeKed++OmtbLXFC7XTc41Aa7uzZs9f+95kzZ7bOxUztDqUV26b7/NfxU3bu9/P6669vnXPTef61/nEovef+TEM/fyvl2FcH8rFScqT3+sMT0mzsayW0UFV9kQkOALJHQJqJoriooviIdgcl1iUByB8BaUbK8hkVxe9rdzsh1iUByBs1pMxYdaImyvKCquoR1fXfyG0nVFWfVlX9w7V2Qu7rDjXt29reQNquQ4RqFlaLHL8uZB0vFoutc9bPO1Qz8seUqoZk1cT8mtFrr722dWzVmEJTxt2fqVVf8o9jtrWQxpn23VcdiPpSNzwhzRDthABMEQFplmgnBGB6CEgztJrg8IhoJwRgSqghzdT+dkIf02LxqU7bkFvrkKw6i7/Gp0sNyV+b4762/z7+sVs3CtXo3O81tD7o3Llze9835Tokt57j1338mpFVU4pZs2RtY+Efx9aQrK0rumxV0fQc8sET0kwVxYvaNQW8qu7nKQlAlghIM7VK2+3e4ryq7h58PAAQQkCasbJ8Srufkr6q5fKewccDABZqSDNyOk/+soriI6rrL2n7V73Qcvm0yvLb623Q49ZzxPSns2oj1nFo3ZFVF7LOhfjfq1vv8Mfk9sST0m1dEerjZ9WQ/LqQX1N69dVX915r1ZBCW1XEbI1ubUfRpe9dl5pSqq3S0Q1PSDO36t7w/h1nWJcEIC8EpANQFN8X7YQA5I6U3YQ1bTO0WZd0up3Qo1ouv6Oi+DfzdbtM+45J2blpHz/N5k879lNg7vce0w7I509Jtra18FN21nTzLtPLrVZCofY/flrOPY7ZusJqKyTZW1WkTNm1bR3UJQVH+m44PCEdiN3thBZaLr/HBAcAWSAgHYh965KkhY6Pn2JtEoDREZAORFFcVFneK/ZLApArakiZc/PXsVtR+MryGUk/WW9r7tZZai2Xv6OiePbUe+46dvP6fi3H2pbcv9avJbjnQ3Ugq4YUmmLtfj9WGxtpux7i101CNST3e4iZ9h0ak1tDCm3Xbh1b9SX/Wv9771JDsj4joe0nUk37TokaUzo8IR2Ysrywt/FqVd051rAAgIB0iPZNcDg+/i4THACMhoB0gKwJDsvl00xwADAKakgzZee1X1ZZ3quq+pK2a0nSqvnqrSrLl699xc/ru/Uaqx7gXxtTQ/JrI6EaksUaY8zaKH8rdL9m5J9vu82FtRbKH2NoWw6rhhRTbwptP+H+vkI/U2s7Cmt7c8leA0eroOnjCelAleUzKst3ig4OAHJBQDpg+yY4LJefZoIDgMGRspuRNlPEd+8su9DVq9/RYvFRLRbnzam3/vv4x9Z0ZmuMofHHTIG30jMxU6xDKTv/OGYqujWduUtX9JjdZq2p3aHXtVoshVJ2Y3T7DiGFNw6ekA4cExwA5IKAdOCK4qIWCzo4ABgfAQkqy/NaLN4lJjgAGBM1pAlput1Ek3t9RfEjY4uK76ksL+x8LWtKuLRdLwjVm9xjvzYSkmpLDP993bpRqIbk14lSTfv2x2/VkELHbn0nph1QaNq3e2xN6w4dd9l+whezVUUM6kv94QkJ1+zr4HD16nfo4ACgdwQkXGNvUfE0U8EB9IqAhGvsLSp4UgLQr6LOMCHadZuFQxHzc7LqNf65ur5Ly+X3dLqtkCQd69y531ZRXJR0um5i1VH8LSX8Y7cmE9umxz32t4WwrrXO+cehmpG1DimmhhRahxRTA7OOY+pNoWutulbK1kHWmqWYLcy71JQy/JOZpTY/J56QcEpZXjCnglfV3UMPCcABICBhp/1TwaWrV7+i4+MPDj8oALNGym5G2rbisdJ5VfUhVdUXdXqFwFJnzvyBFovnt75qpeysc9J2yss6t+u8m6aLSe+F2v9YY/KPu+xia7XEsbqVh1ofxbQdslJtfjsg631DY7LScim7fbedIh5zDvuRskNyq67gf7HjzGqSA09KAFIhICGoLH8gazo4/e4ApEBAQpDd726hq1c/OvSQAMwQNaQZ6TIN3Dq3Oa6qO1VVz+n0dPBjnT37W+t1TPungXeZIh6qIbnnY+pCode1amD+vT7/tZry6yg+tyYT2rHXqufE1JtCdSCrhtRlKneqGpKPad79o4aEXq029PvcjjNMBQfQHQEJUcryKe1K3R0f/z1dHAB0QkBClFU96aOSjr0z9LsD0A01pJlK1VZo3/Fy+T5V1dd2vNpSR0cf02JxXlJcDck9Dq35sY5D64NS1ZBCx6nE1Fy61JBi1julav/jH/fVDij2z1yGfxYnhxoSBhOaCs6TEoBYBCS0Ep4KTmdwAHEISGitLM/r6Ojd4kkJQArUkA5E2z53u8772ypU1T1aLr+k3dtVnNSUYmowMWuW/PMxPedCa4tixhSqvTUVqo3ErNux6kQx9afQtTFjso671IW6rDvK8M/g5FFDwiiszuCrJ6WntFy+jxZDAEwEJCQR2kPp+PirunTpP2jGCmAvUnYHIlVbIf983G6z0mbbirK8ELVdQ5cp411Sg+6x/73mkLILbVXRpU2PlbLrMnXbOp9yKjftgcZFyg6jO3lS8hfObjADD8BuBCQkt5p9d4fK8s/EDDwATRGQ0IvVOqV/WrcZ2h2ULl/+V2pKAK6hhnSg+t6qwj2uqjuNutJJTcmqIYXqN1ZdyLo3Zip3l5pR6OcdM2U5ZguGVLWePreFaFtD8lEzygs1JGTJnoFHTQnACgEJg6CrA4AQAhIGs3pS2l9TunLlWWpKwAGjhoTon3fbmpK7Ffr+mtKxjo4+qMXihyqKi41eV2q/XiimLtRmm442+qy5xKxhalsXih1TX+2A2l6LNKghYRKadHW4fPk/eVoCDgwBCaOwa0oSdSXg8JCywylDtRmSpKr6kKrqi5KOTp1bWaosn9DR0dONU3h9peFiU3KpUnbW+T7TY0N04Q7da10bkuGftoNCyg6TU5bPaLG43ezqUFUP6MqVnzI1HJi5ff9vKTCYorioovimiuLXjH2VVttYSL9SWf5ARfHzgUcJoG88ISEb9r5K0mbCw5UrP2XCAzBD1JAQ1FeboX3nw3UlabUT7btVlhdOvXbKulDKnXab6lJH6WvKeMq6lnVvzLku16J/1JAwC2X5zLpb+N/J2sbi+Pi7Wi4/oar6Q3ajBWaAJyQEDf2E5P7vur5FVXW3quqr2l1bqiUVWj0xfUyLxXmekBoe84SEPvGEhNkJb2Ox+aPPuiVg6nhCQrSU65T2Xbt7zdKdqqrntH97dElaarH4KxXF8yqKF9cz+Jo/IVljijnXpy5PFKna9KS6NnRvqmsxvDa/HwISoo0VkCR/wsMmXec7SeMtFh/V0dGXo9+nyXkCUvdrQ/emuhbDIyBhEGMGJEnrCQy3q67foap6TKEnps0GgLHvEzpPQOp+bejeVNdieAQkDK7L7yrNFPFmabyy/NS1NF5Z7l9UO8YkhpBUAWioe1O2/0l5L4ZFQMLgxg5IUrs03mJxvtcxpURASnMvhkVAwuByCEjSKo1X17epru9UXYfSeMdaLD6wbkF0cesMAan7vQQkSAQkjCCXgOSq67uMDQBdq07ii8VT1wITAan7vQQkSAQkZKCvABVzbVEUDdN4G9s1Jmubixh9LYztcm+qSQNDTT7I8M8TGiIgYXS5BCTJT+M9Krs3nltjuldleT56TDFjtBCQ0tyLcRGQMLqcApJrFZx+z2hB5FpqsXiXyvICAann1+nzXoyLgITR5RqQNpp1Epc2abyy/LdTabwUY7QQkNLci3ERkJC1HIKVtHpaqqq/lPQJNZ0qXhSPnKox5bgwNuW9BB10QUBC1nIJSBsnHR+aTBU/HZzK8metnpy6IiBhCghIyFpuAcm9tlnHhw03OH1+a9r4EAhImII2v1e2nwAkleUFleW9OtkQ0PrHdLLlRV1/UsfHL6qq7ul1fMAh4AkJo0n5e277Wrs6Pki3r/dV2kwVD61jkqRjFcX7VRTf71xj6vOfZNvXTjmmDP/koAek7DApOQYkV7vgtHsCRAwCEuaAgIRJyT0guXYHp71Xy5qdF34vAhKmj4CE2ejzM9D1tVeLbD+uk2nj5tU6CU6fU1F8YWdgGuOfYV/vmeGfFIyAgITZyDkgbWy6P0hfU7PZeZL71CS9cC04EZAwNwQkzMYUAtJGXX9Ydf0lNZ8AcTqlV9c/XZ+7Q26g6hMBCX0iIGE2phSQJHeR7V2Km523uaZaH5eSlpIePvUUlX7MBCT0h4CEgzHGZyS+A4QbnGKdPEVJn1NRPJlN7WnM98V0EJBwMHIOSK7TEyCaPDXtsvupiYCEXBGQcDCmEpA22qf0tl5FJ09ND1+rPU2tnx4OAwEJB2NqAcm1Ozgt12cXip0YIT0kqd9606l3z+/PBjJDQAJ2yPnztAlO0ovrr9wu6S5Jj6l57Wl3cFppPmsvwz8FmDACErDDFD9Pq0B1n+JrT/tm7T0uaffEiJP3zO5PASaMgATsMOXP08kTlPvU1GVixP70XoZ/CjBhBCQgsZw+i2mCk5/eu6BVem/4iRGYNwISkFiun8V+g5O0qT0RqNAWAQlIbAqfxd3ByZ211+hVtLv2xFMU2iEgAYlN7bNY17+h7Vl790l6QO1rT6T40A4BCcAO7hPUZ5QmOD0u6Yn110nv4TQCEoCAVMFJIr0HCwEJQIRdwSm2Y4SL9B5OEJAAtLSvY0Tbpyhm8B06AhKAxHY9RcViBt8hIiAB6NGmndFm1l6X9J60e5IEgWkuCEgABpA6vbdBem9OCEgARpRiBh/pvbkgIAHIBDP4Dh0BCUCGQim+WMzgmwICEoAJ8SdJdJkYsWv/JyZJjImABGCCUqf3NkjvjYmABGDihlqgS3DqGwEJwEylnMFH7WlCCuWlAAADKElEQVQIBCQAByCU4muC6eV9IyABODB+io/9n3JBQAIA9n/KAgEJALaw/9NYCEgAsBfdI4ZEQAKARtj/qW8EJADoJPP9n26QdFbSFUmvt3uJoRCQACCZTPZ/uk7S2yXdLelm5+uvSPqhpB9LuhQ5lAEQkAAguRH3f7pN0p9KOrN+m9K5vVq/7VVJ35D0sxZD6BEBCQAG0/P+T7d9W/rzn5x8eZ/NrV9VVkGJgAQAo0g8g++6/5Ue+E3p6A2prMLXV5KOtcoEZpK+axNa2lTsAABbLq7/71lJX1fn/Z/e/hXpzOtS0fCPeqlVWu/tWtWVJoonJAAYRNP9n2rpvjukX3+peUCSVk9J/yfpye4jTaFNaLEykwCAZC5qNYnhbZLeI+lBrfJs0iq9t07x3fAL6eafxQUjafXX/GZJ1ycY6khI2QHAoALpvbN3dHv5c5Le6PYSYyEgAcBoNsHJOb7ybLeXvNzt9jGRsgOAnLyu1aLXBpPrtlTr+yb6dCQRkAAgPz9U/FrbQpOeYScRkAAgPz/WqgND06ekan39j3sb0SAISACQm0tatQOSwkFpc/4bymZRbFsEJADI0c+0agd0rNWSJT8wVeuvHyu7tkFtsTAWAHIW6vb978pyZh297ABgzq7Xap3RZWU/m46ABADIAq2DAACTRUACAGSBgAQAyAIBCQCQBQISACALBCQAQBYISACALBCQAABZICABALJAQAIAZIGABADIAgEJAJAFAhIAIAsEJABAFghIAIAsEJAAAFkgIAEAskBAAgBkgYAEAMgCAQkAkAUCEgAgCwQkAEAWCEgAgCwQkAAAWSAgAQCycDT2AHap63rsIQAABsYTEgAgCwQkAEAWCEgAgCwQkAAAWSAgAQCyQEACAGSBgAQAyAIBCQCQBQISACALBCQAQBYISACALBCQAABZICABALJAQAIAZIGABADIAgEJAJAFAhIAIAsEJABAFghIAIAsEJAAAFkgIAEAskBAAgBkgYAEAMgCAQkAkAUCEgAgCwQkAEAWCEgAgCwQkAAAWSAgAQCyQEACAGSBgAQAyAIBCQCQBQISACAL/w9vnsZ76YcrSgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31f364950>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imageplot(W) \n",
    "set_cmap(\"gray\")\n",
    "h = plot(imag(gamma), real(gamma), \".b\", linewidth=2)\n",
    "h = plot(x0[2], x0[1], \".r\", markersize=20)\n",
    "h = plot(imag(x1), real(x1), \".g\", markersize=20);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}