examples/user_guide/6_Trimesh.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Triangle Meshes\n",
    "\n",
    "Along with [points](2_Points.ipynb), [timeseries](3_Timeseries.ipynb), [trajectories](4_Trajectories.ipynb), and structured [grids](5_Grids.ipynb), Datashader can rasterize large triangular meshes, such as those often used to simulate data on an irregular grid:\n",
    "\n",
    "<img src=\"../assets/images/chesbay_detail.png\" width=\"500\" height=\"500\" style=\"border-width: 1px; border-style: solid;\">\n",
    "\n",
    "Any polygon can be represented as a set of triangles, and any shape can be approximated by a polygon, so the triangular-mesh support has many potential uses. \n",
    "\n",
    "In each case, the triangular mesh represents (part of) a *surface*, not a volume, and so the result fits directly into a 2D plane rather than requiring 3D rendering. This process of rasterizing a triangular mesh means generating values along specified regularly spaced intervals in the plane.  These examples from the [Direct3D docs](https://msdn.microsoft.com/en-us/library/windows/desktop/cc627092.aspx) show how this process works, for a variety of edge cases:\n",
    "<img width=500 src=\"https://msdn.microsoft.com/dynimg/IC520311.png\"/>\n",
    "\n",
    "This diagram uses \"pixels\" and colors (grayscale), but for datashader the generated raster is more precisely interpreted as a 2D array with bins, not pixels, because the values involved are numeric rather than colors.  (With datashader, colors are assigned only in the later \"shading\" stage, not during rasterization itself.) As shown in the diagram, a pixel (bin) is treated as belonging to a given triangle if its center falls either inside that triangle or along its top or left edge.\n",
    "\n",
    "The specific algorithm used to do so is based on the approach of [Pineda (1998)](http://people.csail.mit.edu/ericchan/bib/pdf/p17-pineda.pdf), which has the following features:\n",
    "  * Classification of pixels relies on triangle convexity\n",
    "  * Embarrassingly parallel linear calculations\n",
    "  * Inner loop can be calculated incrementally, i.e. with very \"cheap\" computations\n",
    "  \n",
    "and a few assumptions: \n",
    "  * Triangles should be non overlapping (to ensure repeatable results for different numbers of cores)\n",
    "  * Triangles should be specified consistently either in clockwise or in counterclockwise order of vertices (winding). \n",
    "  \n",
    "Trimesh rasterization is not yet GPU-accelerated, but it's fast because of [Numba](http://numba.pydata.org) compiling Python into SIMD machine code instructions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiny example\n",
    "\n",
    "To start with, let's generate a tiny set of 10 vertices at random locations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, datashader as ds, pandas as pd\n",
    "import datashader.utils as du, datashader.transfer_functions as tf\n",
    "from scipy.spatial import Delaunay\n",
    "import dask.dataframe as dd\n",
    "\n",
    "n = 10\n",
    "np.random.seed(2)\n",
    "\n",
    "x = np.random.uniform(size=n)\n",
    "y = np.random.uniform(size=n)\n",
    "z = np.random.uniform(0,1.0,x.shape)\n",
    "\n",
    "pts = np.stack((x,y,z)).T\n",
    "verts = pd.DataFrame(np.stack((x,y,z)).T, columns=['x', 'y' , 'z'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have a set of random x,y locations and associated z values.  We can see the numeric values with \"head\" and plot them (with color for z) using datashader's usual points plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvs = ds.Canvas(plot_height=400,plot_width=400)\n",
    "\n",
    "tf.Images(verts.head(15), tf.spread(tf.shade(cvs.points(verts, 'x', 'y', agg=ds.mean('z')), name='Points')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a trimesh, we need to connect these points together into a non-overlapping set of triangles.  One well-established way of doing so is [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def triangulate(vertices, x=\"x\", y=\"y\"):\n",
    "    \"\"\"\n",
    "    Generate a triangular mesh for the given x,y,z vertices, using Delaunay triangulation.\n",
    "    For large n, typically results in about double the number of triangles as vertices.\n",
    "    \"\"\"\n",
    "    triang = Delaunay(vertices[[x,y]].values)\n",
    "    print('Given', len(vertices), \"vertices, created\", len(triang.simplices), 'triangles.')\n",
    "    return pd.DataFrame(triang.simplices, columns=['v0', 'v1', 'v2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time tris = triangulate(verts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of triangulation is a set of triangles, each composed of three indexes into the vertices array.  The triangle data can then be visualized by datashader's ``trimesh()`` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.Images(tris.head(15), tf.shade(cvs.trimesh(verts, tris)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, datashader will rasterize your trimesh using z values [linearly interpolated between the z values that are specified at the vertices](https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Interpolation_on_a_triangular_unstructured_grid).  The shading will then show these z values as colors, as above.  You can enable or disable interpolation as you wish:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from colorcet import rainbow as c\n",
    "tf.Images(tf.shade(cvs.trimesh(verts, tris, interpolate='nearest'), cmap=c, name='10 Vertices'),\n",
    "          tf.shade(cvs.trimesh(verts, tris, interpolate='linear'),  cmap=c, name='10 Vertices Interpolated'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More complex example\n",
    "\n",
    "The small example above should demonstrate how triangle-mesh rasterization works, but in practice datashader is intended for much larger datasets. Let's consider a sine-based function `f` whose frequency varies with radius:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rad = 0.05,1.0\n",
    "\n",
    "def f(x,y):\n",
    "    rsq = x**2+y**2\n",
    "    return np.where(np.logical_or(rsq<rad[0],rsq>rad[1]), np.nan, np.sin(10/rsq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily visualize this function by sampling it on a raster with a regular grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 400\n",
    "\n",
    "ls  = np.linspace(-1.0, 1.0, n)\n",
    "x,y = np.meshgrid(ls, ls)\n",
    "img = f(x,y)\n",
    "\n",
    "raster = tf.shade(tf.Image(img, name=\"Raster\"))\n",
    "raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, you can see pronounced aliasing towards the center of this function, as the frequency starts to exceed the sampling density of the raster.  Instead of sampling at regularly spaced locations like this, let's try evaluating the function at random locations whose density varies towards the center:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polar_dropoff(n, r_start=0.0, r_end=1.0):\n",
    "    ls = np.linspace(0, 1.0, n)\n",
    "    ex = np.exp(2-5*ls)/np.exp(2)\n",
    "    radius = r_start+(r_end-r_start)*ex\n",
    "    theta  = np.random.uniform(0.0,1.0, n)*np.pi*2.0\n",
    "    x = radius * np.cos( theta )\n",
    "    y = radius * np.sin( theta )\n",
    "    return x,y\n",
    "\n",
    "x,y = polar_dropoff(n*n, np.sqrt(rad[0]), np.sqrt(rad[1]))\n",
    "z = f(x,y)\n",
    "\n",
    "verts = pd.DataFrame(np.stack((x,y,z)).T, columns=['x', 'y' , 'z'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the x,y points and optionally color them with the z value (the value of the function f(x,y)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvs = ds.Canvas(plot_height=400,plot_width=400)\n",
    "\n",
    "tf.Images(tf.shade(cvs.points(verts, 'x', 'y'), name='Points'),\n",
    "          tf.shade(cvs.points(verts, 'x', 'y', agg=ds.mean('z')), name='PointsZ'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points are clearly covering the area of the function that needs dense sampling, and the shape of the function can (roughly) be made out when the points are colored in the plot.  But let's go ahead and triangulate so that we can interpolate between the sampled values for display:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time tris = triangulate(verts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's pre-compute the combined mesh data structure for these vertices and triangles, which for very large meshes (much larger than this one!) would save plotting time later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time mesh = du.mesh(verts,tris)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This mesh can be used for all future plots as long as we don't change the number or ordering of vertices or triangles, which saves time for much larger grids.\n",
    "\n",
    "We can now plot the trimesh to get an approximation of the function with noisy sampling locally to disrupt the interference patterns observed in the regular-grid version above and preserve fidelity where it is needed. (Usually one wouldn't do this just for the purposes of plotting a function, since the eventual display on a screen is a raster image no matter what, but having a variable grid is crucial if running a simulation where fine detail is needed only in certain regions.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.shade(cvs.trimesh(verts, tris, mesh=mesh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fine detail in the heavily sampled regions is visible when zooming in closer (without resampling the function):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.Images(*([tf.shade(ds.Canvas(x_range=r, y_range=r).trimesh(verts, tris, mesh=mesh))\n",
    "            for r in [(0.1,0.8), (0.14,0.4), (0.15,0.2)]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the central disk is being filled in above, even though the function is not defined in the center.  That's a limitation of Delaunay triangulation, which will create convex regions covering the provided vertices.  You can use other tools for creating triangulations that have holes, align along certain regions, have specified densities, etc., such as [MeshPy](https://mathema.tician.de/software/meshpy) (Python bindings for [Triangle](http://www.cs.cmu.edu/~quake/triangle.html)).\n",
    "\n",
    "\n",
    "### Aggregation functions\n",
    "\n",
    "Like other datashader methods, the ``trimesh()`` method accepts an ``agg`` argument (defaulting to ``mean()``) for a reduction function that determines how the values from multiple triangles will contribute to the value of a given pixel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.Images(tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.mean('z')),name='mean'),\n",
    "          tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.max('z')), name='max'),\n",
    "          tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.min('z')), name='min'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three plots above should be nearly identical, except near the center disk where individual pixels start to have contributions from a large number of triangles covering different portions of the function space.  In this inner ring, ``mean`` reports the average value of the surface inside that pixel, ``max`` reports the maximum value of the surface (hence being darker values in this color scheme), and ``Min`` reports the minimum value contained in each pixel.  The ``min`` and ``max`` reductions are useful when looking at a very large mesh, revealing details not currently visible. For instance, if a mesh has a deep but very narrow trough, it will still show up in the ``min`` plot regardless of your raster's resolution, while it might be missed on the ``mean`` plot.  \n",
    "\n",
    "Other reduction functions are useful for making a mask of the meshed area (``any``), for showing how many triangles are present in a given pixel (``count``), and for reporting the diversity of values within each pixel (``std`` and ``var``):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.Images(tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.any('z')), name='any'),\n",
    "          tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.count()),  name='count'),\n",
    "          tf.shade(cvs.trimesh(verts, tris, mesh=mesh, agg=ds.std('z')), name='std')).cols(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallelizing trimesh aggregation with Dask\n",
    "The trimesh aggregation process can be parallelized by providing `du.mesh` and `Canvas.trimesh` with partitioned Dask dataframes.\n",
    "\n",
    "**Note:** While the calls to `Canvas.trimesh` will be parallelized across the partitions of the Dask dataframe, the construction of the partitioned mesh using `du.mesh` is not currently parallelized.  Furthermore, it currently requires loading the entire `verts` and `tris` dataframes into memory in order to construct the partitioned mesh.  Because of these constraints, this approach is most useful for the repeated aggregation of large meshes that fit in memory on a single multicore machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "verts_ddf = dd.from_pandas(verts, npartitions=4)\n",
    "tris_ddf = dd.from_pandas(tris, npartitions=4)\n",
    "mesh_ddf = du.mesh(verts_ddf, tris_ddf)\n",
    "mesh_ddf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.shade(cvs.trimesh(verts_ddf, tris_ddf, mesh=mesh_ddf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive plots\n",
    "\n",
    "By their nature, fully exploring irregular grids needs to be interactive, because the resolution of the screen and the visual system are fixed.  Trimesh renderings can be generated as above and then displayed interactively using the datashader support in [HoloViews](http://holoviews.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews.operation.datashader import datashade\n",
    "hv.extension(\"bokeh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "HoloViews is designed to make working with data easier, including support for large or small trimeshes. With HoloViews, you first declare a ``hv.Trimesh`` object, then you apply the ``datashade()`` (or just ``aggregate()``) operation if the data is large enough to require datashader.  Notice that HoloViews expects the triangles and vertices in the *opposite* order as datashader's ``cvs.trimesh()``, because the vertices are optional for HoloViews:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wireframe = datashade(hv.TriMesh((tris,verts), label=\"Wireframe\").edgepaths)\n",
    "trimesh   = datashade(hv.TriMesh((tris,hv.Points(verts, vdims='z')), label=\"TriMesh\"), aggregator=ds.mean('z'))\n",
    "\n",
    "(wireframe + trimesh).opts(width=400, height=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can zoom in on either of these plots, but they will only update if you have a live Python server (not a static web page.  The Wireframe plot will initially look like a collection of dots (as the triangles are all tiny), but zooming in will reveal the shape (if you are just looking at the static web page, eventually you will see individual pixels in the original datashaded rasterized plot, not the full trimesh available).  Notice how a few of the \"wires\" cross the center, because Delaunay triangulation has filled in the central region; other techniques as mentioned previously would be needed to avoid those.\n",
    "\n",
    "For examples of Datashader's trimesh in use, see the [Chesapeake and Delaware Bays](https://examples.pyviz.org/bay_trimesh/bay_trimesh.html) notebook:\n",
    "\n",
    "<img src=\"../assets/images/chesapeake_farout.png\" width=\"600\">"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}