ampullaImage.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<script>\n",
    "  function code_toggle() {\n",
    "    if (code_shown){\n",
    "      $('div.input').hide('500');\n",
    "      $('#toggleButton').val('Show Code')\n",
    "    } else {\n",
    "      $('div.input').show('500');\n",
    "      $('#toggleButton').val('Hide Code')\n",
    "    }\n",
    "    code_shown = !code_shown\n",
    "  }\n",
    "\n",
    "  $( document ).ready(function(){\n",
    "    code_shown=false;\n",
    "    $('div.input').hide()\n",
    "  });\n",
    "</script>\n",
    "<form action=\"javascript:code_toggle()\"><input type=\"submit\" id=\"toggleButton\" value=\"Show Code\"></form>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Avoid computing the boundaries again\n",
    "*note*: For the next 3 scenarios, it is necessary to click on \"cell>Run All above\" after browsing to the target cell and clicking on it, this will load the needed dependencies.  \n",
    "- 1D displacement test\n",
    "- 2D source point test\n",
    "- Generate EVs in a pre processed cross section  \n",
    "\n",
    "[load the dependencies](#Dependencies)\n",
    "\n",
    "*Otherwise:*  \n",
    "[Load the dependencies](#Dependencies) first and then proceed to [Process a new cross section](#The-full-pipeline) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from imageio import imread\n",
    "from skimage import measure\n",
    "from skimage.transform import rescale\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import collections as mc\n",
    "%matplotlib inline\n",
    "from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar\n",
    "import matplotlib.font_manager as fm\n",
    "from shapely.geometry import Polygon\n",
    "from shapely.geometry.polygon import LinearRing\n",
    "import os\n",
    "import sys\n",
    "import math\n",
    "import random\n",
    "from math import radians\n",
    "import copy\n",
    "import pickle\n",
    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\n",
    "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n",
    "from struct import pack\n",
    "\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "from utilities.geometry import Plane, middle_point, two_points_normals_2d, project_point_on_direction, distance\n",
    "\n",
    "\n",
    "class Segment(object):\n",
    "    def __init__(self, x1, y1, x2, y2, mx, my, nx, ny, l, s_id):\n",
    "        self.p1x = x1\n",
    "        self.p1y = y1\n",
    "        self.p2x = x2\n",
    "        self.p2y = y2\n",
    "        self.mpx = mx\n",
    "        self.mpy = my\n",
    "        mag = np.sqrt(nx * nx + ny * ny)\n",
    "        self.direction_x = nx / mag\n",
    "        self.direction_y = ny / mag\n",
    "        self.normalx = nx\n",
    "        self.normaly = ny\n",
    "        self.length = l\n",
    "        self.section_id = s_id\n",
    "        \n",
    "    def location(self):\n",
    "        return np.array([self.mpx, self.mpy])\n",
    "\n",
    "    def direction(self):\n",
    "        return np.array([self.direction_x, self.direction_y])\n",
    "\n",
    "class Cell(Segment):\n",
    "    def __init__(self, _id, p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id):\n",
    "        super(Cell, self).__init__(p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id)\n",
    "        self._id = _id\n",
    "        self.radius = diameter / 2.\n",
    "    \n",
    "    def getXmlDescription(self):\n",
    "        raise NotImplementedError(\"Method not yet implemented!!!\")\n",
    "    \n",
    "    def getCoords(self):\n",
    "        return [(self.p1x, self.p1y),(self.p2x, self.p2y)]\n",
    "    \n",
    "    def angle(self):\n",
    "        return np.arctan2(self.p2y - self.p1y, self.p2x - self.p1x) #* 180./np.pi\n",
    "    \n",
    "    def getXmlDescription(self):\n",
    "        return \"\"\"\n",
    "    <xagent>\n",
    "        <name>{}</name>\n",
    "        <id>{:d}</id>\n",
    "        <p1_x>{:f}</p1_x>\n",
    "        <p1_y>{:f}</p1_y>\n",
    "        <p2_x>{:f}</p2_x>\n",
    "        <p2_y>{:f}</p2_y>\n",
    "        <x>{:f}</x>\n",
    "        <y>{:f}</y>\n",
    "        <direction_x>{:f}</direction_x>\n",
    "        <direction_y>{:f}</direction_y>\n",
    "        <angle>{:f}</angle>\n",
    "        <time_since_last_secreted>0</time_since_last_secreted>\n",
    "        <section>{:d}</section>\n",
    "    </xagent>\"\"\".format(self.name, self._id, self.p1x, self.p1y, self.p2x, self.p2y, self.mpx, \n",
    "                            self.mpy, self.direction_x, self.direction_y, self.angle(), self.section_id)\n",
    "        \n",
    "class SecretoryCell(Cell):\n",
    "    def __init__(self, _id, p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id):\n",
    "        super(SecretoryCell, self).__init__(_id, p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id)\n",
    "        self.name = 'SecretoryCell'\n",
    "\n",
    "\n",
    "class CiliaryCell(Cell):\n",
    "    def __init__(self, _id, p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id):\n",
    "        super(CiliaryCell, self).__init__(_id, p1x, p1y, p2x, p2y, mpx, mpy, dirx, diry, diameter, section_id)\n",
    "        self.name = 'CiliaryCell'\n",
    "\n",
    "\n",
    "class EV:\n",
    "    # to convert:\n",
    "    # N square metres to square micrometre: multiply N * 1e+12\n",
    "    # N metres to micrometres: multiply N * 1e+6\n",
    "    # N micrometres to metres: divide N / 1e+6   \n",
    "    def random_diameter_nm(self):\n",
    "        return random.randint(3, 30) * 10\n",
    "    \n",
    "    def mass_from_diameter(self):\n",
    "        # update function to compute the mass as a function of the volume\n",
    "        # assume spherical shape and same material density \n",
    "        # with an estimated mass of 5.47990700763866E-20 for an EV of radius 100 nm\n",
    "        # Mass = density * volume, Volume = 4/3*Pi*radius^3\n",
    "        # return 5.47991E-20 * self.diameter_um\n",
    "        ev_mass_per_radius_um = {\n",
    "        30:1.47957E-21, 40:3.50714E-21, 50:6.84988E-21, 60:1.18366E-20, 70:1.87961E-20,\n",
    "        80:2.80571E-20, 90:3.99485E-20, 100:5.47991E-20, 110:7.29376E-20, 120:9.46928E-20,\n",
    "        130:1.20394E-19, 140:1.50369E-19, 150:1.84947E-19, 160:2.24457E-19, 170:2.69228E-19,\n",
    "        180:3.19588E-19, 190:3.75867E-19, 200:4.38393E-19, 210:5.07494E-19, 220:5.835E-19,\n",
    "        230:6.6674E-19, 240:7.57542E-19, 250:8.56235E-19, 260:9.63148E-19, 270:1.07861E-18,\n",
    "        280:1.20295E-18, 290:1.33649E-18, 300:1.47957E-18}\n",
    "        \n",
    "        return ev_mass_per_radius_um[self.diameter_nm]\n",
    "    \n",
    "    def __init__(self, *args, **kwargs):\n",
    "        self._id = kwargs.get('_id', 0)\n",
    "        self.x_1 = kwargs.get('x_1', 0) # holds the starting position for this EV\n",
    "        self.y_1 = kwargs.get('y_1', 0)\n",
    "        self.vx = kwargs.get('vx', 0)\n",
    "        self.vy = kwargs.get('vy', 0)\n",
    "        self.diameter_nm = self.random_diameter_nm()\n",
    "        self.diameter_um = kwargs.get('diameter_nm', self.diameter_nm) / 1000\n",
    "        self.diameter_m = self.diameter_um / 1e6\n",
    "        \n",
    "        self.colour = '()'\n",
    "        \n",
    "        self.radius_um = (self.diameter_um / 2)\n",
    "        self.radius_m = self.radius_um / 1e+6\n",
    "        \n",
    "        self.mass_Kg = self.mass_from_diameter()\n",
    "\n",
    "        # Compute the radius-dependent diffusion coefficient defined as\n",
    "        # D = (boltzmann tempK)/(6 pi fluidViscosity particleRadius)\n",
    "        self.diffusion_rate_m = ((const_Boltzmann * const_Temperature_K) \n",
    "                                 / (6. * math.pi * const_water_viscosity * self.radius_m))\n",
    "        self.diffusion_rate_um = self.diffusion_rate_m * 1e+12\n",
    "        \n",
    "        # compute the mean displacement per second which is diffusion-rate dependat\n",
    "        self.velocity_ums = math.sqrt(4 * self.diffusion_rate_um)\n",
    "        self.velocity_ms = math.sqrt(4 * self.diffusion_rate_m)\n",
    "        \n",
    "        # initially the EV is located at the origin and is static\n",
    "        self.direction_x = 0\n",
    "        self.direction_y = 0\n",
    "        self.x = self.x_1\n",
    "        self.y = self.y_1\n",
    "        \n",
    "        # Originally, the EV is assigned to the 0 visual section\n",
    "        self.section_id = 0\n",
    "        \n",
    "    def unit_normal_from_direction(self, dir_x, dir_y):\n",
    "        magnitude = math.sqrt(dir_x * dir_x + dir_y * dir_y)\n",
    "        self.direction_x = dir_x / magnitude\n",
    "        self.direction_y = dir_y / magnitude\n",
    "\n",
    "    def set_direction(self, dir_x, dir_y):\n",
    "        self.unit_normal_from_direction(self, dir_x, dir_y)\n",
    "        \n",
    "    def location(self):\n",
    "        return np.array([self.x, self.y])\n",
    "    \n",
    "    def getXmlDescription(self):\n",
    "        return \"\"\"\n",
    "    <xagent>\n",
    "        <name>EV</name>\n",
    "        <id>{:d}</id>\n",
    "        <x>{:f}</x>\n",
    "        <y>{:f}</y>\n",
    "        <x_1>{:f}</x_1>\n",
    "        <y_1>{:f}</y_1>\n",
    "        <vx>{:.10f}</vx>\n",
    "        <vy>{:.10f}</vy>\n",
    "        <mass_kg>{:.25f}f</mass_kg>\n",
    "        <colour/>\n",
    "        <radius_um>{:f}</radius_um>\n",
    "        <radius_m>{:.10f}</radius_m>\n",
    "        <diameter_um>{:f}</diameter_um>\n",
    "        <diameter_m>{:f}</diameter_m>\n",
    "        <diffusion_rate_um>{:.25f}</diffusion_rate_um>\n",
    "        <diffusion_rate_m>{:.25f}</diffusion_rate_m>\n",
    "        <direction_x>{:f}</direction_x>\n",
    "        <direction_y>{:f}</direction_y>\n",
    "        <velocity_ums>{:.10f}</velocity_ums>\n",
    "        <velocity_ms>{:.10f}</velocity_ms>\n",
    "        <closest_ev_id>-1</closest_ev_id>\n",
    "        <closest_ev_distance>-1.f</closest_ev_distance>\n",
    "        <closest_cell_id>-1</closest_cell_id>\n",
    "        <closest_cell_distance>-1.f</closest_cell_distance>\n",
    "        <section>{:d}</section>\n",
    "    </xagent>\"\"\".format(self._id, self.x, self.y,self.x_1, self.y_1, self.vx, self.vy, self.mass_Kg,\n",
    "                        self.radius_um, self.radius_m, self.diameter_um, self.diameter_m,\n",
    "                        self.diffusion_rate_um, self.diffusion_rate_m, \n",
    "                        self.direction_x, self.direction_y, self.velocity_ums, self.velocity_ms, self.section_id)\n",
    "\n",
    "def plot_normal(p1, p2, mp, symbol, normal):\n",
    "    new_p = project_point_on_direction(direction=normal, source=mp)\n",
    "    plt.plot([mp[0],new_p[0]], [mp[1],new_p[1]], symbol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Image processing code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def get_image_foreground(img):\n",
    "#    return (img[:,:,3] != 0) & (img[:,:,1] > 240)\n",
    "\n",
    "def show_image_foreground(img, foreground):\n",
    "    fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1)\n",
    "    ax1.axis([-50, 2600, -50, 1900])\n",
    "    ax2.axis([-50, 2600, -50, 1900])\n",
    "    ax1.imshow(img)\n",
    "    ax2.imshow(foreground)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_contours_in_foreground(img, img_upscaled, foreground, save=True, show=True, min_perimeter = 90):\n",
    "    contours = measure.find_contours(foreground, 0.8)\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1, 2,figsize=[20,10])\n",
    "    \n",
    "    greater_than_N = []\n",
    "\n",
    "    # show image\n",
    "    ax1.set_title('Sow ampulla cross section (orig. size)')\n",
    "    ax1.imshow(img)\n",
    "    ax1.set_aspect('equal')\n",
    "    \n",
    "    ax2.set_title('Identified contours (upscaled image)')\n",
    "    ax2.set_aspect('equal')\n",
    "    ax1.axis([0, img.shape[1], 0, img.shape[0]])\n",
    "    ax2.axis([0, img_upscaled.shape[1], 0, img_upscaled.shape[0]])\n",
    "    ax2.imshow(img_upscaled, interpolation='nearest', cmap=plt.cm.gray)\n",
    "\n",
    "    for n, contour in enumerate(contours):\n",
    "        ring = LinearRing(contour)\n",
    "        if ring.length >= min_perimeter:\n",
    "            greater_than_N.append(n)\n",
    "            ax2.plot(contour[:,1], contour[:,0], linewidth=2, color='orange')\n",
    "\n",
    "    fontprops = fm.FontProperties(size=12)\n",
    "    scalebar1 = AnchoredSizeBar(ax1.transData,\n",
    "                               45, '100 um', 'lower right', \n",
    "                               pad=0.1,\n",
    "                               frameon=False,\n",
    "                               size_vertical=1,\n",
    "                               fontproperties=fontprops)\n",
    "    scalebar2 = AnchoredSizeBar(ax2.transData,\n",
    "                               100, '100 um', 'lower right', \n",
    "                               pad=0.1,\n",
    "                               frameon=False,\n",
    "                               size_vertical=1,\n",
    "                               fontproperties=fontprops)\n",
    "\n",
    "    ax2.add_artist(scalebar2)\n",
    "    ax1.add_artist(scalebar1)\n",
    "    if save:\n",
    "        fig.savefig('01_ampulla_cross_section_and_identified_contours_'+base_filename+'.png', bbox_inches='tight', dpi=300)\n",
    "    elif show:        \n",
    "        fig.showfig()\n",
    "\n",
    "    return contours, greater_than_N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def simplify_contours(contours, greater_than_N, avg_cell_size, tolerance, scale = 1.0):\n",
    "    selected_simplified_contours = []\n",
    "    polygons = []\n",
    "\n",
    "    total_perimeter = 0.\n",
    "    secretory_cells_fitted = 0.\n",
    "    ciliary_cells_fitted = 0.\n",
    "    total_vertices = 0.\n",
    "\n",
    "    for i in range(len(greater_than_N)):\n",
    "        contour = measure.approximate_polygon(contours[greater_than_N[i]], tolerance)\n",
    "        # scale the contours here\n",
    "        contour *= scale\n",
    "        selected_simplified_contours.append(contour)\n",
    "        \n",
    "        polygon = Polygon(contour)\n",
    "        polygons.append(polygon)\n",
    "\n",
    "        total_perimeter += polygon.length\n",
    "        total_vertices += contour.shape[0]\n",
    "\n",
    "        for j in range(contour.shape[0]-1):\n",
    "            d = distance(contour[j], contour[j+1])\n",
    "            #print(f'edge d:{d} would fit {int(d / avg_cell_size)} secretory cells')\n",
    "            if d >= avg_cell_size:\n",
    "                # generate n secretory cells in this edge\n",
    "                # 1 compute the edge normal\n",
    "                # 2 subsection\n",
    "                secretory_cells_fitted += int(d / avg_cell_size)\n",
    "            else:\n",
    "                # generate a ciliary cell patch here\n",
    "                ciliary_cells_fitted += 1\n",
    "    print(f'|{tolerance}|{total_vertices}|{total_perimeter:.2f}|{secretory_cells_fitted}|{ciliary_cells_fitted}|-> total:{secretory_cells_fitted+ciliary_cells_fitted}')\n",
    "    return selected_simplified_contours, polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the specified edge/shape\n",
    "def display_edge_with_normals(polygons, avg_cell_size, i=0):\n",
    "    e_coords = np.array(polygons[i].exterior.coords)\n",
    "    fig, ax = plt.subplots(1, 1, sharex='col',figsize=[80,40])\n",
    "    # plt.figure()\n",
    "    for idx in range(len(e_coords) - 1):\n",
    "        p1 = e_coords[idx]\n",
    "        p2 = e_coords[idx + 1]\n",
    "        normal = two_points_normals_2d(p1, p2)[0] # [0 if i in inside_another else 1]\n",
    "\n",
    "        d = distance(p1, p2)\n",
    "        if d >= avg_cell_size:\n",
    "            source_p = p1\n",
    "            n = d // avg_cell_size\n",
    "            # print(f'edge {idx} d:{d:.3f} sub-sections:{n}')\n",
    "            p1p2 = p2 - p1\n",
    "            if n >= 2:\n",
    "                # fin n full size cells\n",
    "                step = 1./n\n",
    "                for st in range(int(n)):\n",
    "                    dest_p = project_point_on_direction(source=source_p, direction = p1p2, dist=step)\n",
    "                    mp = middle_point(source_p, dest_p)\n",
    "                    plot_normal(source_p, dest_p, mp, '-g', normal)\n",
    "                    # the edge\n",
    "                    plt.plot([source_p[0], dest_p[0]], [source_p[1], dest_p[1]], '-.')\n",
    "                    source_p = dest_p\n",
    "            else:\n",
    "                # fit one full size cell\n",
    "                extra = d - avg_cell_size\n",
    "                dest_p = project_point_on_direction(source=source_p, direction = p1p2, dist=1./d * avg_cell_size)\n",
    "                plt.plot([source_p[0], dest_p[0]], [source_p[1], dest_p[1]], '--')\n",
    "                mp = middle_point(source_p, dest_p)\n",
    "                plot_normal(source_p, dest_p, mp, '--b', normal)\n",
    "                source_p = dest_p\n",
    "                # fit one partial size cell\n",
    "                plt.plot([source_p[0],p2[0]], [source_p[1],p2[1]], '-k')  \n",
    "                mp = middle_point(source_p, p2)\n",
    "                plot_normal(source_p, p2, mp, '--b', normal)\n",
    "\n",
    "        else:\n",
    "            # the normal\n",
    "            mp = middle_point(p1, p2)\n",
    "            plot_normal(p1, p2, mp, '-b', normal)\n",
    "            # the edge\n",
    "            plt.plot([p1[0],p2[0]], [p1[1],p2[1]], '-y')\n",
    "\n",
    "    fontprops = fm.FontProperties(size=12)\n",
    "    scalebar = AnchoredSizeBar(ax.transData,\n",
    "                               22.5, '50 um', 'lower left', \n",
    "                               pad=0.1,\n",
    "                               frameon=False,\n",
    "                               size_vertical=1,\n",
    "                               fontproperties=fontprops)\n",
    "\n",
    "    ax.add_artist(scalebar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deprecated code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# identify the coordinates for the segments, the midpoints, the directions, the lengths \n",
    "def process_edge(segments, p1, p2, avg_cell_size, section_id):\n",
    "    normal = two_points_normals_2d(p1, p2)[0]  # [1 if i in inside_another else 0]\n",
    "\n",
    "    d = distance(p1, p2)\n",
    "    # is this section wider/larger than the avg cell size?\n",
    "    if d >= avg_cell_size:\n",
    "        source_p = p1\n",
    "        diameter = d / avg_cell_size\n",
    "\n",
    "        # how many cells can we fit?\n",
    "        n = int(diameter)\n",
    "        p1p2 = p2 - p1\n",
    "        if n >= 2:\n",
    "            # fit n cells\n",
    "            step = 1./n\n",
    "            for st in range(int(n)):\n",
    "                # compute the coordinates for this source point\n",
    "                dest_p = project_point_on_direction(source=source_p, direction = p1p2, dist=step)\n",
    "                mp = middle_point(source_p, dest_p)\n",
    "\n",
    "                segments.append(Segment(source_p[0], source_p[1], dest_p[0], dest_p[1], mp[0], mp[1], normal[0],normal[1], diameter, section_id))\n",
    "\n",
    "                source_p = dest_p\n",
    "        else:\n",
    "            # fit one cell\n",
    "            dest_p = project_point_on_direction(source=source_p, direction = p1p2, dist=1./d * avg_cell_size)\n",
    "            mp = middle_point(source_p, dest_p)\n",
    "\n",
    "            segments.append(Segment(source_p[0], source_p[1], dest_p[0], dest_p[1], mp[0], mp[1],normal[0],normal[1], avg_cell_size, section_id))\n",
    "\n",
    "            source_p = dest_p\n",
    "            mp = middle_point(source_p, p2)\n",
    "\n",
    "            segments.append(Segment(source_p[0], source_p[1], p2[0], p2[1], mp[0], mp[1], normal[0],normal[1], d - avg_cell_size, section_id))\n",
    "\n",
    "    else:\n",
    "        # only one cell fits in this section\n",
    "        mp = middle_point(p1, p2)\n",
    "        segments.append(Segment(p1[0], p1[1], p2[0], p2[1], mp[0], mp[1], normal[0],normal[1], d, section_id))\n",
    "    return segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def identify_segments_in_polygons(polygons, avg_cell_size):\n",
    "    segments = []\n",
    "    for i in range(len(polygons)):\n",
    "        e_coords = np.array(polygons[i].exterior.coords)\n",
    "\n",
    "        for idx in range(len(e_coords) - 1):\n",
    "            p1 = e_coords[idx]\n",
    "            p2 = e_coords[idx + 1]\n",
    "            segments = process_edge(segments, p1, p2, avg_cell_size, i)\n",
    "            #print(f'Computed {len(segments)} segments in this polygon')\n",
    "    return segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def produce_cells(segments, probabilistic = False, threshold_secretory = 0.8):\n",
    "    # above this threshold the cell will be secretory\n",
    "    \n",
    "    # use the identified segments to assign the cells\n",
    "    cells = []\n",
    "    secretory = []\n",
    "    ciliary = []\n",
    "\n",
    "    countCiliary = 0\n",
    "    countSecretory = 0\n",
    "    cell_id = 0\n",
    "    section = 0\n",
    "\n",
    "    for s in segments:\n",
    "        if probabilistic:\n",
    "            if random.random() >= threshold_secretory:\n",
    "                cell = SecretoryCell(cell_id, s.p1x, s.p1y, s.p2x, s.p2y, s.mpx, s.mpy, s.normalx, s.normaly, s.length, s.section_id)\n",
    "                secretory.append(cell)\n",
    "            else:\n",
    "                cell = CiliaryCell(cell_id, s.p1x, s.p1y, s.p2x, s.p2y, s.mpx, s.mpy, s.normalx, s.normaly, s.length, s.section_id)\n",
    "                ciliary.append(cell)\n",
    "        else:\n",
    "            # to contol the experimetation, we can use a deterministic aproach\n",
    "            cell = SecretoryCell(cell_id, s.p1x, s.p1y, s.p2x, s.p2y, s.mpx, s.mpy, s.normalx, s.normaly, s.length, s.section_id)\n",
    "            secretory.append(cell)\n",
    "\n",
    "        cells.append(cell)\n",
    "        cell_id += 1\n",
    "        cell = None\n",
    "\n",
    "    return cells, secretory, ciliary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# secretory, ciliary, collection_ciliaries, collection_secretories\n",
    "def export_computed_positions(target_filename, threshold_secretory, cells=None, secretory=None, ciliary=None):\n",
    "    s = int(threshold_secretory * 100)\n",
    "    c = str(100 - s)\n",
    "    s = str(s)\n",
    "    filename = f'{base_filename}__distribution_{s}-{c}_s-c_'\n",
    "    if secretory:\n",
    "        flnm = f'{filename}_secretory.pickle'\n",
    "        with open(flnm, 'bw') as pickle_target:\n",
    "            pickle.dump(secretory, pickle_target)\n",
    "            print(f'Exported file: {flnm}')\n",
    "    if ciliary:\n",
    "        flnm = f'{filename}_ciliary.pickle'\n",
    "        with open(flnm, 'bw') as pickle_target:\n",
    "            pickle.dump(ciliary, pickle_target)\n",
    "            print(f'Exported file: {flnm}')\n",
    "    if cells:\n",
    "        flnm = f'{filename}_allCells.pickle'\n",
    "        with open(flnm, 'bw') as pickle_target:\n",
    "            pickle.dump(cells, pickle_target)\n",
    "            print(f'Exported file: {flnm}')\n",
    "\n",
    "def load_distributions(base_filename, pct_sec, pct_cil):\n",
    "    s = str(pct_sec)\n",
    "    c = str(pct_cil)\n",
    "    filename = f'{base_filename}__distribution_{s}-{c}_s-c_'\n",
    "    with open(f'{filename}_allCells.pickle', 'rb') as pickle_source:\n",
    "        loaded_cells = pickle.load(pickle_source)\n",
    "\n",
    "    with open(f'{filename}_secretory.pickle', 'rb') as pickle_source:\n",
    "        loaded_secretory = pickle.load(pickle_source)\n",
    "\n",
    "    with open(f'{filename}_ciliary.pickle', 'rb') as pickle_source:\n",
    "        loaded_ciliary = pickle.load(pickle_source)\n",
    "\n",
    "    return {'cells': copy.deepcopy(loaded_cells), \n",
    "            'ciliary': copy.deepcopy(loaded_ciliary), 'secretory': copy.deepcopy(loaded_secretory)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# File format description:\n",
    "def export_v1_bin_files(polygons, source_points):\n",
    "    # the old format depended on having the sections as 'visual sections'\n",
    "    visual_sections = len(polygons)\n",
    "    \n",
    "    boundaries_filename = '0_boundaries.evbin'\n",
    "    source_points_filename = '0_starting_points.evbin'\n",
    "    \n",
    "    with open(boundaries_filename, 'bw') as f_out:\n",
    "        # write the number of sections\n",
    "        f_out.write(pack('<i', visual_sections))\n",
    "        \n",
    "        # write the number of vertices per section\n",
    "        vertices_per_visual_section = [len(polygon.exterior.coords)+1 for polygon in polygons]\n",
    "        #print(vertices_per_visual_section)\n",
    "        f_out.write(pack('<'+'i'*len(vertices_per_visual_section), *vertices_per_visual_section))\n",
    "       \n",
    "        #for section in range(visual_sections):\n",
    "        #    vertices_per_visual_section = len(polygons[section].exterior.coords)\n",
    "        #    f_out.write(pack('i', vertices_per_visual_section))\n",
    "        max_x = None\n",
    "        max_y = None\n",
    "        min_x = None\n",
    "        min_y = None\n",
    "        \n",
    "        # the file should contain N vertices, and N-1 angles.\n",
    "        # these vertices will provide the perimeter for each visual section\n",
    "        for polygon in polygons:\n",
    "            coords =  np.append(np.array(polygon.exterior.coords), np.array([polygon.exterior.coords[0]]), 0)\n",
    "            #coords = np.array(polygon.exterior.coords)\n",
    "            #print(coords.shape)\n",
    "            #print(coords[:,1])\n",
    "            n = len(coords)\n",
    "            # firstly, all the x coordinates are written\n",
    "            f_out.write(pack('<'+'f'*n, *coords[:,0]))\n",
    "            \n",
    "            mx_x = max(coords[:,0])\n",
    "            mn_x = min(coords[:,0])\n",
    "            if max_x is None or max_x < mx_x:\n",
    "                max_x = mx_x\n",
    "            if min_x is None or min_x > mn_x:\n",
    "                min_x = mn_x\n",
    "        #print(coords[-10:,0])\n",
    "        \n",
    "        for polygon in polygons:\n",
    "            coords =  np.append(np.array(polygon.exterior.coords), np.array([polygon.exterior.coords[0]]), 0)\n",
    "            n = len(coords)\n",
    "            # then, all the y coordinates are written\n",
    "            f_out.write(pack('<'+'f'*n, *coords[:,1]))\n",
    "            \n",
    "            #print(f'wrote {n} xs and ys')\n",
    "            # find max x,y and min x,y for this section\n",
    "            \n",
    "            mx_y = max(coords[:,1])\n",
    "            mn_y = min(coords[:,1])\n",
    "            if max_y is None or max_y < mx_y:\n",
    "                max_y = mx_y\n",
    "            if min_y is None or min_y > mn_y:\n",
    "                min_y = mn_y\n",
    "        # print(coords[-10:,1])\n",
    "\n",
    "        for polygon in polygons:\n",
    "            coords =  np.append(np.array(polygon.exterior.coords), np.array([polygon.exterior.coords[0]]), 0)\n",
    "            n = len(coords)\n",
    "            # finally, all the angles are written\n",
    "            for s in range(1, n):\n",
    "                f_out.write(pack('<f', np.arctan2(coords[s, 1] - coords[s - 1, 1], coords[s, 0] - coords[s - 1, 0]) ))\n",
    "            #f_out.write(pack('f', np.arctan2(coords[0, 1] - coords[n-1, 1], coords[0, 0] - coords[n-1, 0]) ))\n",
    "            # print(f'wrote {n-1} angles')\n",
    "        \n",
    "        # at the very last we need to write the max x,y and min x,y values\n",
    "        f_out.write(pack('ffff', max_x, max_y, min_x, min_y))\n",
    "        # print(max_x, max_y, min_x, min_y)\n",
    "        print(f'boundaries file written: {boundaries_filename}')\n",
    "\n",
    "    \n",
    "    # A separate file contains the starting points. \n",
    "    # These are the mid points for the secretory cells under version 2 of the model\n",
    "    number_of_starting_points = len(source_points)\n",
    "    \n",
    "    # for each starting point we need to export:\n",
    "    # x, y, angle (in degrees), the index of the visual section it belongs to\n",
    "    with open(source_points_filename, 'bw') as bout:\n",
    "        bout.write(pack('i', number_of_starting_points))\n",
    "        \n",
    "        for sp in source_points:\n",
    "            bout.write(pack('ffff', sp.mpx, sp.mpy, np.arctan2(sp.direction_y, sp.direction_x), sp.section_id))\n",
    "        \n",
    "        print(f'source_points file written: {source_points_filename}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_distributions(secretory, ciliary, threshold_secretory, save=False):\n",
    "    \n",
    "\n",
    "    collection_secretories = np.array([s.getCoords() for s in secretory])\n",
    "    collection_ciliaries = np.array([c.getCoords() for c in ciliary])\n",
    "\n",
    "    col_ciliaries = np.zeros((collection_ciliaries.shape[0], 2, 2))\n",
    "    col_secretories = np.zeros((collection_secretories.shape[0], 2, 2))\n",
    "    col_ciliaries[:,:,0] =collection_ciliaries[:,:,1]\n",
    "    col_ciliaries[:,:,1] =collection_ciliaries[:,:,0]\n",
    "    col_secretories[:,:,0] = collection_secretories[:,:,1]\n",
    "    col_secretories[:,:,1] = collection_secretories[:,:,0]\n",
    "\n",
    "    s = int(threshold_secretory * 100)\n",
    "    c = str(100 - s)\n",
    "    s = str(s)\n",
    "\n",
    "    zoom = 14\n",
    "    fig, ax1 = plt.subplots(1, 1,figsize=[20,10])\n",
    "\n",
    "    pct_secretory = 100 / len(cells) * len(secretory)\n",
    "    pct_ciliary = 100 / len(cells) * len(ciliary)\n",
    "    ax1.set_title('Hypothetical cell distribution at the ampulla')\n",
    "    ax1.add_collection(mc.LineCollection(col_secretories, color='red', linewidths=2))\n",
    "    ax1.add_collection(mc.LineCollection(col_ciliaries, color='blue', linewidths=2))\n",
    "    ax1.legend(ax1.collections[:2], [f'Secretory cells: {len(secretory):,d} ({pct_secretory:.1f}%)', \n",
    "                                     f'Ciliary cells: {len(ciliary):,d} ({pct_ciliary:.1f}%)']);\n",
    "    #ax1.autoscale()\n",
    "    ax1.axis([0, 2550, 0, 1850])\n",
    "    ax1.set_aspect('equal')\n",
    "    fontprops = fm.FontProperties(size=12)\n",
    "    scalebar = AnchoredSizeBar(ax1.transData,\n",
    "                               45, '100 um', 'lower right', \n",
    "                               pad=0.1,\n",
    "                               frameon=False,\n",
    "                               size_vertical=1,\n",
    "                               fontproperties=fontprops)\n",
    "    ax1.add_artist(scalebar)\n",
    "\n",
    "    axins = zoomed_inset_axes(ax1, zoom, loc=3) # zoom = 6\n",
    "    axins.add_collection(mc.LineCollection(col_secretories, color='red', linewidths=2))\n",
    "    axins.add_collection(mc.LineCollection(col_ciliaries, color='blue', linewidths=2))\n",
    "    axins.set_title(\"Zoomed in x14\")\n",
    "\n",
    "    # sub region of the original image\n",
    "    x1, x2, y1, y2 = 940, 990, 600, 640\n",
    "    axins.set_xlim(x1, x2)\n",
    "    axins.set_ylim(y1, y2)\n",
    "\n",
    "    plt.xticks(visible=False)\n",
    "    plt.yticks(visible=False)\n",
    "    fontprops = fm.FontProperties(size=12)\n",
    "    scalebar = AnchoredSizeBar(ax1.transData,\n",
    "                               22.5* zoom, '50 um', 'lower right', \n",
    "                               pad=0.1,\n",
    "                               frameon=False,\n",
    "                               size_vertical=1,\n",
    "                               fontproperties=fontprops)\n",
    "    axins.add_artist(scalebar)\n",
    "\n",
    "    # draw a bbox of the region of the inset axes in the parent axes and\n",
    "    # connecting lines between the bbox and the inset axes area\n",
    "    mark_inset(ax1, axins, loc1=2, loc2=4, fc=\"none\", ec=\"0.5\")\n",
    "    if save:\n",
    "        fig.savefig('03_ampulla_hypothetical_cell_distributions_'+s+'-'+c+'_s-c_secretory_'+base_filename+'.png', bbox_inches='tight', dpi=300)\n",
    "    else:\n",
    "        fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_initial_impulse(ev, scale):\n",
    "    init_angle = random.randrange(0,359,1)\n",
    "    direction_x = math.cos(radians(init_angle))\n",
    "    direction_y = math.sin(radians(init_angle))\n",
    "    return displace_in_direction(ev, direction_x, direction_y, scale)\n",
    "\n",
    "def displace_in_direction(ev, dir_x, dir_y, scale):\n",
    "    # compute the unit normal\n",
    "    ev.unit_normal_from_direction(dir_x, dir_y)\n",
    "\n",
    "    ev.vx = ev.velocity_ums * ev.direction_x\n",
    "    ev.vy = ev.velocity_ums * ev.direction_y\n",
    "\n",
    "    # displace the ev further by a sigle step\n",
    "    ev.x = ev.x_1 + ev.vx * scale\n",
    "    ev.y = ev.y_1 + ev.vy * scale\n",
    "    return ev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_evs_from_source_points(max_evs, source_points):\n",
    "    \"\"\"\n",
    "    Generates the N requested source points by randomly selecting from those provided.\n",
    "    The source points provided are copied to a stack, then an iterative process pops one element from the stack for each EV. \n",
    "    When the stack is empy a new stack is generated.\n",
    "    This process is repeated until the desired number of source points is generated\n",
    "    \"\"\"\n",
    "    evs_per_source_point = {}\n",
    "    \n",
    "    source_point_stack = copy.deepcopy(source_points)\n",
    "    random.shuffle(source_point_stack)\n",
    "    evs = []\n",
    "    source_point = source_point_stack.pop()\n",
    "    \n",
    "    for i in range(max_evs):\n",
    "        # generate a single EV\n",
    "        ev = EV(_id=i)\n",
    "        \n",
    "        # we need to keep track of the evs assigned to each cell\n",
    "        if source_point._id in evs_per_source_point:\n",
    "            # if we already have evs in the current source point we fetch the last cell and use its location\n",
    "            # as the initial location for the new EV. However, to avoid detecting false positive collisions\n",
    "            # we compute a new location from the radius of the current last element\n",
    "            \n",
    "            last_element_index = len(evs_per_source_point[source_point._id]['evs'])\n",
    "            last_ev = evs_per_source_point[source_point._id]['evs'][last_element_index - 1]\n",
    "            # compute the location to clear the previous EV radius\n",
    "            #print(type(last_ev), type(last_ev.location()), last_ev.location())\n",
    "            new_point = project_point_on_direction(last_ev.location(), source_point.direction(), last_ev.radius_um)\n",
    "            next_x = new_point[0]\n",
    "            next_y = new_point[1]\n",
    "        else:\n",
    "            # otherwise, we just add this starting position to the collection\n",
    "            evs_per_source_point[source_point._id] = {'evs':[]}\n",
    "            next_x = source_point.mpx\n",
    "            next_y = source_point.mpy\n",
    "            \n",
    "        # assign the location for the EV\n",
    "        ev.x_1 = next_x\n",
    "        ev.y_1 = next_y\n",
    "        \n",
    "        # displace the EV by its radius\n",
    "        ev = displace_in_direction(ev, source_point.normalx, source_point.normaly, ev.radius_um)\n",
    "\n",
    "        # keep track of the EV for this cell\n",
    "        evs_per_source_point[source_point._id]['evs'].append(ev)\n",
    "        evs.append(ev)\n",
    "        \n",
    "        # prevent evs with crazy velocities\n",
    "        if np.abs(ev.vx) >= 10 or np.abs(ev.vy) >= 10:\n",
    "            print(f'warning: {i} dirx:{sp.dirx} diry:{sp.diry}')\n",
    "        \n",
    "        if len(source_point_stack) > 0:\n",
    "            source_point = source_point_stack.pop()\n",
    "        else:\n",
    "            source_point_stack = copy.deepcopy(source_points)\n",
    "            source_point = source_point_stack.pop()\n",
    "    return evs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_evs_in_circle(n, radius, random_direction=False):\n",
    "    evs = []\n",
    "    for i in range(n):\n",
    "        a = random.random() * 2 * np.pi\n",
    "        r = radius * np.sqrt(random.random())\n",
    "        x = r * np.cos(a)\n",
    "        y = r * np.sin(a)\n",
    "        magnitude = np.sqrt(x*x + y*y)\n",
    "        \n",
    "        if random_direction:\n",
    "            a = random.random() * 2 * np.pi\n",
    "            direction = np.array([np.cos(a), np.sin(a)])\n",
    "        else:\n",
    "            direction = np.array([x/magnitude, y/magnitude])\n",
    "        \n",
    "        ev = EV(_id=i, x_1=x, y_1=y)\n",
    "        #ev.x_1 = x\n",
    "        #ev.y_1 = y\n",
    "        ev = displace_in_direction(ev, direction[0], direction[1], ev.radius_um)\n",
    "        evs.append(ev)\n",
    "    return evs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_0xml_file(evs, n, base_filename, version=0, boundaries = False, ev_ev_collisions = False, cells=None):\n",
    "\n",
    "    f0 = \"\"\"D:\\iterations\\\\v{}\\\\{}\\\\0_{}_{}.xml\"\"\".format(version,\n",
    "                                                          'bounded' if boundaries else 'unbounded', base_filename, n)\n",
    "    with open(f0, 'w') as of:\n",
    "        header = \"\"\"<states>\n",
    "    <itno>0</itno>\n",
    "    <environment>\n",
    "        <dt>{}</dt>\n",
    "        <boundaries_enabled>0</boundaries_enabled>\n",
    "        <const_ev_physics>0</const_ev_physics>\n",
    "        <const_ev_collisions>{}</const_ev_collisions>\n",
    "        <introduce_n_new_evs>0</introduce_n_new_evs>\n",
    "        <seconds_before_intoducing_new_evs>2</seconds_before_intoducing_new_evs>\n",
    "        <introducing_new_evs_threshold>0.85</introducing_new_evs_threshold>\n",
    "    </environment>\"\"\"\n",
    "\n",
    "        of.write(header.format(dt, '1' if ev_ev_collisions else '0'))\n",
    "        if boundaries:\n",
    "            # write the boundaries\n",
    "            for cell in cells:\n",
    "                of.write(cell.getXmlDescription())\n",
    "\n",
    "        # write the evs\n",
    "        for ev in evs:\n",
    "            of.write(ev.getXmlDescription())\n",
    "\n",
    "        # close the file\n",
    "        of.write(\"\"\"\n",
    "    </states>\"\"\")\n",
    "\n",
    "    print(f'file saved as: {f0}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For a temperature of 20C (293.15 Kelvin) water has a viscosity of 0.0010016\n",
      "const_Boltzmann * const_Temperature_K = 4.0473766476e-21\n",
      "6. * math.pi * const_water_viscosity = 0.018879715211013223\n"
     ]
    }
   ],
   "source": [
    "# Dynamic viscosity of water\n",
    "# source http://www.viscopedia.com/viscosity-tables/substances/water/\n",
    "water_viscosity = {20: 0.0010016, 21: 0.0009775, 22: 0.0009544, 23: 0.0009321, 24: 0.0009107, 25: 0.00089}\n",
    "\n",
    "tempC = 20\n",
    "const_Boltzmann = 1.3806504e-23\n",
    "# Temperature (Kelvins) 20C + 273.15\n",
    "const_Temperature_K = tempC + 273.15\n",
    "const_water_viscosity = water_viscosity[tempC]\n",
    "print(f\"For a temperature of {tempC}C ({const_Temperature_K} Kelvin) water has a viscosity of {const_water_viscosity}\")\n",
    "print(f\"const_Boltzmann * const_Temperature_K = {const_Boltzmann * const_Temperature_K}\")\n",
    "print(f\"6. * math.pi * const_water_viscosity = {6. * math.pi * const_water_viscosity}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dependencies\n",
    "After clicking on this cell, select `Cell>Run All above`, this will load the needed dependencies. Then proceed to the desired scenario  \n",
    "- [1D displacement test](#1D-displacement-test)\n",
    "- [2D source point test](#2D-source-point-test)\n",
    "- [Generate EVs in a pre processed cross section](#Load-previously-computed-positions-from-file)  \n",
    "\n",
    "or [Process a new cross section](#The-full-pipeline) \n",
    "\n",
    "# Simple tests\n",
    "\n",
    "## 1D displacement test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file saved as: D:\\iterations\\v14\\unbounded\\0_1D_displacement_test_50000.xml\n"
     ]
    }
   ],
   "source": [
    "# 1D scenario\n",
    "# Segments for experimentation\n",
    "segments = []\n",
    "coords = np.array([(-500,0),(-500,900)])\n",
    "avg_cell_size = 2.25\n",
    "\n",
    "segments = process_edge(segments, coords[0], coords[1], avg_cell_size, 0)\n",
    "len(segments)\n",
    "dt = 0.001983891\n",
    "\n",
    "base_filename = '1D_displacement_test'\n",
    "\n",
    "cells, secretory, ciliary = produce_cells(segments, False)\n",
    "# produce the EVs\n",
    "\n",
    "#n = 131072\n",
    "n = 50000\n",
    "version = 14\n",
    "new_evs = generate_evs_from_source_points(n, secretory)\n",
    "# apply an initial displacement to the EVs\n",
    "\n",
    "save_0xml_file(new_evs, n, base_filename, version=version, boundaries = False, ev_ev_collisions = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D source point test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_filename = '2D_source_point_test'\n",
    "dt = 0.001983891\n",
    "\n",
    "n = 131072\n",
    "version = 14\n",
    "radius = 45\n",
    "evs = generate_evs_in_circle(n, radius)\n",
    "\n",
    "save_0xml_file(evs, n, base_filename, version=13, boundaries=False, ev_ev_collisions=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D source point test with random direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_filename = '2D_source_point_randir'\n",
    "dt = 0.001983891\n",
    "\n",
    "n = 100000\n",
    "version = 14\n",
    "radius = 45\n",
    "evs = generate_evs_in_circle(n, radius, random_direction=True)\n",
    "\n",
    "save_0xml_file(evs, n, base_filename, version=13, boundaries=False, ev_ev_collisions=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A model of Extracellular Vesicles behaviour\n",
    "The state of the art about extracellular vesicle offers significant detail about the biological characterisations of their content, we identified some gaps regarding other systemic characteristics such as production rate, diffusion coefficient, mass, lifespan, and decay or decomposition rate. This complicated things a bit for creating our first model of EV behaviour. However, we expect that our findings will eventually help reduce this gap.\n",
    "\n",
    "## The environment\n",
    "\n",
    "After measuring the dimensions of the oviduct from histological images, we have identified the width of 6000nm and the height of 3000nm at the ampulla, whereas in the isthmus the width is 1000nm and height 1000nm. Therefore, the whole area of interest measures around 6 by 3 mm in the ampulla and 1 by 1 mm in the isthmus.\n",
    "\n",
    "Here we load the data from the environment from a set of previously prepared files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loaded from previously computed positions from file\n",
    "\n",
    "Previously computed distributions are stored in the directory 'distributions'. The files there are related to the experiment they belong by the first part of their name. As an example, the for the experiment 'ampulla_(3-1to3-5_L2-3_tr-amp)\\_2' would produce the files:  \n",
    "\n",
    "$$\n",
    "V2 files\n",
    "\\begin{cases}\n",
    "\\text{ampulla_(3-1to3-5_L2-3_tr-amp)_2__distribution_80-20_s-c__allCells.pickle} & \\text{w/ all cells}\\\\\n",
    "\\text{ampulla_(3-1to3-5_L2-3_tr-amp)_2__distribution_80-20_s-c__ciliary.pickle} & \\text{ciliary cells}\\\\\n",
    "\\text{ampulla_(3-1to3-5_L2-3_tr-amp)_2__distribution_80-20_s-c__secretory.pickle} & \\text{secretory cells}\n",
    "\\end{cases}\n",
    "$$\n",
    "*Note: the number 80, 20 corresponds to an 80/20 percent distribution of secretory/ciliary cells*\n",
    "\n",
    "$$\n",
    "\\text{v1 files}\n",
    "\\begin{cases}\n",
    "\\text{`ampulla_(3-1to3-5_L2-3_tr-amp)_2_0_boundaries.evbin`} & \\text{the boundaries} \\\\\n",
    "\\text{`ampulla_(3-1to3-5_L2-3_tr-amp)_2_0_starting_points.evbin`} & \\text{the starting points}  \n",
    "\\end{cases}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preset environments:\n",
    "presets = [{'f':'ampulla_(3-1to3-5_L2-3_tr-amp)_2', 's':'80', 'c':'20'}]\n",
    "p = 0\n",
    "base_filename = presets[p]['f']\n",
    "\n",
    "distributions = load_distributions(base_filename, presets[p]['s'], presets[p]['c'])\n",
    "len(distributions['cells']), len(distributions['ciliary']), len(distributions['secretory'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Vesicles\n",
    "Each vesicle is represented as a particle, assigning parameters such as mass, position, and velocity. Then, the interactions between the vesicles or with the oviduct were decided to be treated as perfectly plastic collisions observing the principles of conservation of momentum, angle, and energy.\n",
    "\n",
    "### Size\n",
    "As previously discussed in our literature review, it is well accepted that in general the range of EV sizes lies within 10 and 1500 nm. And depending on their biogenesis, this range becomes narrower. In the case of oviductosomes, experimental data has shown they typically have diametral size between 30 and 300 nm. However, more specific characterisation has narrowed down this range even further to be between 30 and 150 nm. Therefore, for the first version of our model, we decided to work with sizes in this range only (30 to 150 nm) using increments of 10 nm.\n",
    "\n",
    "### Mass\n",
    "From the published data of the molecular mass for an exosome with a diameter of 100nm, produced by blood plasma cells, we derived a mass of 1.34e-24 Kg per vesicle. Then, the corresponding mass for the range of sizes described was obtained by simple extrapolation. This procedure gave us the following masses.  \n",
    "\n",
    "| EV diameter (nm) | Mass (Kg) |  \n",
    "|------------------|-----------|\n",
    "|30|1.64397E-20|\n",
    "|40|2.19196E-20|\n",
    "|50|2.73995E-20|\n",
    "|60|3.28794E-20|\n",
    "|70|3.83593E-20|\n",
    "|80|4.38393E-20|\n",
    "|90|4.93192E-20|\n",
    "|100|5.47991E-20|\n",
    "|110|6.0279E-20|\n",
    "|120|6.57589E-20|\n",
    "|130|7.12388E-20|\n",
    "|140|7.67187E-20|\n",
    "|150|8.21986E-20|\n",
    "\n",
    "### Displacement\n",
    "Given that the vesicles are governed by physical laws, we estimate the corresponding diffusion rate for each vesicle based on their diameter. For this we use the Stokes-Einstein equation where $k_B$ is the Boltzman constant, $T$ is the absolute temperature. And the frictional coefficient, acting against the spherical particle, is given by the divisor. Where n is the viscosity of water at T temperature, and r is the radius of the vesicle.\n",
    "\n",
    "$$\n",
    "D = \\frac{k_B T}{6 \\pi \\eta}\n",
    "$$\n",
    "The Stokes-Einstein equation  \n",
    "\n",
    "We use the computed value for the displacement to estimate the initial velocity of the particle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Segments for experimentation\n",
    "#segments = []\n",
    "avg_cell_size = 2.25\n",
    "\n",
    "segments = identify_segments_in_polygons(polygons, avg_cell_size)\n",
    "#len(segments)\n",
    "dt = 0.001983891\n",
    "\n",
    "# produce the EVs\n",
    "\n",
    "n = 25000\n",
    "new_evs = generate_evs_from_source_points(n, distributions['secretory'])\n",
    "# apply an initial displacement to the EVs\n",
    "\n",
    "save_0xml_file(new_evs, n, base_filename, version=13, boundaries = False, ev_ev_collisions = False, \n",
    "               cells=distributions['cells'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The full pipeline\n",
    "\n",
    "## Identifying the cross section boundaries and internal regions\n",
    "\n",
    "After discarding the background colour, the pixels in the foreground with no transparency and which colour is close to white are selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pixels_to_um(x):\n",
    "    return x * (100 / 45)\n",
    "\n",
    "def pixels_to_mm(x):\n",
    "    return pixels_to_um(x)/1000\n",
    "\n",
    "def pixels_to_cm(x):\n",
    "    return pixels_to_mm(x)/10\n",
    "\n",
    "pixels_to_um(2500), 100/45"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_filename = 'ampulla_(3-1to3-5_L2-3_tr-amp)_2'\n",
    "\n",
    "img = imread('./img/'+base_filename+'.png')\n",
    "\n",
    "# scale the image to make 1 pixel = 1 um\n",
    "img_upscaled = rescale(img, 100/45, anti_aliasing = True)\n",
    "# get the foreground\n",
    "# foreground = img[:,:,3] != 0 & img[:,:,1] > 240\n",
    "foreground = (img_upscaled[:,:,3] > 0.2) & (img_upscaled[:,:,1] > 0.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because pixels are single and independent representations of colour we need to process them to obtain some useful geometrical information. This is done firstly by using [`skimage.measure.find-contours`](http://scikit-image.org/docs/0.14.x/api/skimage.measure.html#find-contours) to identify their contours. \n",
    "\n",
    "Then, these contours are converted into geometrical shapes using the [`shapely`](https://github.com/Toblerity/Shapely) package.  \n",
    "\n",
    "We can now proceed to check and display those having a perimeter $\\ge$ 20 pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_perimeter = 200\n",
    "contours, greater_than_N = find_contours_in_foreground(img, img_upscaled, foreground, save = False, show = False, min_perimeter=min_perimeter)\n",
    "# how many contours were found?, how many of them have more than N edges, and what are those?\n",
    "print(f'A total of {len(contours)} contours were identified. From these, {len(greater_than_N)} have a perimeter of at least {min_perimeter} pixels.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These will be further processed in the following steps to obtain the physical limits of the model. \n",
    "\n",
    "## EV potential starting point coordinates\n",
    "\n",
    "Given that the extra cellular vesicles are produced by the oviduct epitheluim, the edges of the identified contours provide the vesicle's initial coordinates and the hypothetical locations of the source points producing the vesicles of our interest.\n",
    "\n",
    "\n",
    "### Number of cells per edge\n",
    "\n",
    "The manner in which the number of cells per edge was decided changed after each iteration of the modelling activity. Here,  a brief description of the major changes is presented.  \n",
    "\n",
    "#### v1\n",
    "Took the middle point between each pair of vertices on the polygon as a source point. The source points and the oviduct boundaries were independent concepts. Both were loaded as separate binary files, one was targetting the visual representation in visualisation, the other was purely oriented to collision detection.  \n",
    "#### v2\n",
    "The second version performed fixed-lenght subdivision of the poligon edges. If edgeLenght > 0.7 subdivide each edge in 6 parts, each part became a source point. Still source points and boundaries were independent concepts. Two binary files were still required.  \n",
    "#### v3 (current)\n",
    "After consulting with Dr Paul Richmond on the performance issues observed in our previous models, the number of checks for collisions was flagged as the main need for optimisation. It was sugested to change the logical representation of the boundaries from a list of line segments and turn them into agents. Then, instead of checking every vesicle agent for possible collisions against every segment in the boundaries, only those within the possible displacement range would be checked.\n",
    "\n",
    "Such approach was estimated to require extensive modification to our model yet it would also enable exploring modelling other biological concepts observed in real cells. One of such concepts is the variations in cell distribution at different sections of the oviduct, which interestingly is also reported to change in close relation with the ovulation cycle. \n",
    "\n",
    "Therefore, after the changes were implemented, the new cell agents introduced in our model represent the phisical boundaries of the oviduct as obtained from the histology images but they also enable representing different types of cells with distinct behaviour.\n",
    "\n",
    "Firsty, for the number of cells in a given edge, each edge is subdivided by an average oviductal cell diameter of $5 \\mu m$ [x]. Then, in order to represent an average cell distributions of 50/50 ciliary/secretory [xx] at the ampulla, the fate of each cell is decided by sampling from an normally distributed random number generator (RNG). The cell will be ciliary if randomNumber &lt; 0.5 or secretory otherwise.  \n",
    "\n",
    "Allocating the hypothetical cells on the edges.\n",
    "\n",
    "### Shape simplification \n",
    "Due to the high complexity of the contour, most of the original edges have shorter length than the avg cells size. This works against our cell placement purposes. A solution for this is to simplify the shapes using less vertices. To do so, we use [`skimage.measure.approximate-polygon`](http://scikit-image.org/docs/0.14.x/api/skimage.measure.html#approximate-polygon).\n",
    "\n",
    "After testing with different values for the tolerance we decided to use 0.1 because the shape it produces the closest perimeter to the original while producing a large number of cells which would enable reproducing the distributions described by [x]\n",
    "\n",
    "| Tolerance | Vertices | Perimeter | Avg sized cells* | Smaller sized cells |\n",
    "|-----------|----------|-------|-----------------|---------------------|\n",
    "|0.0|76643.0|64701.93|0.0|76628.0|\n",
    "|0.1|15814.0|64701.93|20539.0|7695.0|\n",
    "|0.2|15318.0|64665.40|20539.0|7199.0|\n",
    "|0.25|15316.0|64665.22|20539.0|7197.0|\n",
    "|0.3|14816.0|64603.16|20766.0|6486.0|\n",
    "|0.35|13988.0|64505.53|21150.0|5410.0|\n",
    "|0.4|13368.0|64433.65|21479.0|4769.0|\n",
    "|0.5|11791.0|64201.80|21992.0|3178.0|\n",
    "|1.0|7154.0|63120.93|24337.0|209.0|\n",
    "*note:* Average-sized-cells = $\\sum$int(edge lenght/avg cell size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the cells fitted in the simplified contours\n",
    "tolerances = [0.0, 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 1.0]\n",
    "avg_cell_size = 2.25\n",
    "contours, greater_than_N, avg_cell_size\n",
    "for tolerance in tolerances:\n",
    "    simplify_contours(contours, greater_than_N, avg_cell_size, tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simplify the contours without comparison\n",
    "avg_cell_size = 2.25\n",
    "simplified_contours, polygons = simplify_contours(contours, greater_than_N, avg_cell_size, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the biological process under scrutiny, the EVs are released from the oviduct wall and their initial displacement is to the inside of the tube. To capture this in our model, we need to define an initial direction for every source point. This direction will be directly dependant on the spatial charachteristics of each shape. \n",
    "\n",
    "### Visualise the source point directions\n",
    "Our source point generation algorithm depends on the computed normal for each edge to provide the initial direction in which the EVs will displace. Because the polygons generated are described following the Right Hand Coordinate System (RHS)[http://www.ntu.edu.sg/home/ehchua/programming/opengl/CG_BasicsTheory.html] or Counter-Clockwise (CCW)[https://shapely.readthedocs.io/en/latest/manual.html#shapely.geometry.polygon.orient] direction, weu\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional\n",
    "display_edge_with_normals(polygons, avg_cell_size=2.25, i=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the cell distributions\n",
    "\n",
    "First identify the segments in each edge of the polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "segments = identify_segments_in_polygons(polygons, avg_cell_size)\n",
    "print(f'Computed {len(segments)} segments in all polygons')\n",
    "threshold_secretory = 0.8\n",
    "cells, secretory, ciliary = produce_cells(segments, True, threshold_secretory)\n",
    "print(f'Generated {len(cells)} cells. The distribution is: {len(secretory)} secretory and {len(ciliary)} ciliary')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export & visualise the computed positions\n",
    "\n",
    "The theoretical cell positions and their distributions are saved to the following files for future reuse\n",
    "\n",
    "secretory, ciliary all_cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "export_computed_positions(base_filename, threshold_secretory, cells, secretory, ciliary)\n",
    "# visualize_distributions(secretory, ciliary, threshold_secretory, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export Old file format\n",
    "The computed segments/boundaries need to be exported in a format compatible with V1 of the model. This is only for comparison purposes between both models. \n",
    "\n",
    "Then, the two files generated need to be copied to the directory where the .EXE file is located prior to lauching the application. Their names must be '0_boundaries.evbin' and '0_starting_points.evbin'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(vertices_per_visual_section, sum(vertices_per_visual_section))\n",
    "export_v1_bin_files(polygons, secretory)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<script>\n",
    "  $(document).ready(function(){\n",
    "    $('div.prompt').hide();\n",
    "    $('div.back-to-top').hide();\n",
    "    $('nav#menubar').hide();\n",
    "    $('.breadcrumb').hide();\n",
    "    $('.hidden-print').hide();\n",
    "  });\n",
    "</script>\n",
    "\n",
    "<footer id=\"attribution\" style=\"float:right; color:#999; background:#fff;\">\n",
    "Created with Jupyter, delivered by Fastly, rendered by Rackspace.\n",
    "</footer>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}