Greatly improve mesh grid generation : more than 100x, look at mesh speed comparison notebook

Author: xavArtley

.gitignore

@@ -95,6 +95,13 @@ ENV/
 # Rope project settings
 .ropeproject
 
+# Eclipse stuff
+.settings
+.project
+
+# Pydev stuff
+.pydevproject
+
 # Webpack watch files
 js/node_modules/
 # Javascript Build files

ipyvolume/pylab.py

@@ -321,35 +321,7 @@ def reshape_color(ar):
     if v is not None:
         v = reshape(v)
     _grow_limits(np.array(x).reshape(-1), np.array(y).reshape(-1), np.array(z).reshape(-1))
-    mx = nx if wrapx else nx - 1
-    my = ny if wrapy else ny - 1
-    triangles = np.zeros(((mx) * (my) * 2, 3), dtype=np.uint32)
-    lines = np.zeros(((mx) * (my) * 4, 2), dtype=np.uint32)
-
-    def index_from2d(i, j):
-        xi = (i % nx)
-        yi = (j % ny)
-        return ny * xi + yi
-        """
-        ^ ydir
-        |
-        2 3
-        0 1  ---> x dir
-        """
-
-    for i in range(mx):
-        for j in range(my):
-            p0 = index_from2d(i, j)
-            p1 = index_from2d(i + 1, j)
-            p2 = index_from2d(i, j + 1)
-            p3 = index_from2d(i + 1, j + 1)
-            triangle_index = (i * my) + j
-            triangles[triangle_index * 2 + 0, :] = [p0, p1, p3]
-            triangles[triangle_index * 2 + 1, :] = [p0, p3, p2]
-            lines[triangle_index * 4 + 0, :] = [p0, p1]
-            lines[triangle_index * 4 + 1, :] = [p0, p2]
-            lines[triangle_index * 4 + 2, :] = [p2, p3]
-            lines[triangle_index * 4 + 3, :] = [p1, p3]
+    triangles, lines = _make_triangles_lines(x.reshape(nx,ny),y.reshape(nx,ny),z.reshape(nx,ny),wrapx,wrapy)
     # print(i, j, p0, p1, p2, p3)
     mesh = ipv.Mesh(x=x, y=y, z=z, triangles=triangles if surface else None, color=color,
                        lines=lines if wireframe else None,
@@ -1004,3 +976,88 @@ def lasso(data, other=None, fig=fig):
                     if selected:
                         scatter.selected = selected
     fig.on_lasso(lasso)
+    
+
+
+def _make_triangles_lines(x, y, z, wrapx=False, wrapy=False):
+    """Transform rectangular regular grid into triangles
+    
+    :param x: {x2d}
+    :param y: {y2d}
+    :param z: {z2d}
+    :param bool wrapx: when True, the x direction is assumed to wrap, and polygons are drawn between the end end begin points
+    :param bool wrapy: simular for the y coordinate
+    :return: triangles and lines used to plot Mesh
+    """
+    assert len(x.shape) == 2, "Array x must be 2 dimensional."
+    assert len(y.shape) == 2, "Array y must be 2 dimensional."
+    assert len(z.shape) == 2, "Array z must be 2 dimensional."
+    assert x.shape == y.shape, "Arrays x and y must have same shape."
+    assert y.shape == z.shape, "Arrays y and z must have same shape."
+    
+    nx, ny = x.shape
+
+    mx = nx if wrapx else nx - 1
+    my = ny if wrapy else ny - 1
+    
+    """
+    create all pair of indices (i,j) of the rectangular grid
+    minus last row if wrapx = False => mx
+    minus last column if wrapy = False => my
+    |  (0,0)   ...   (0,j)    ...   (0,my-1)  |
+    |    .      .      .       .       .      |
+    |  (i,0)   ...   (i,j)    ...   (i,my-1)  |
+    |    .      .      .       .       .      |
+    |(mx-1,0)  ...  (mx-1,j)  ... (mx-1,my-1) |
+    """
+    i, j = np.mgrid[0:mx, 0:my]
+
+    """
+    collapsed i and j in one dimensional array, row-major order
+    ex :
+    array([[0,  1,  2],     =>   array([0, 1, 2, 3, *4*, 5])
+           [3, *4*, 5]])
+    if we want vertex 4 at (i=1,j=1) we must transform it in i*ny+j = 4
+    """
+    i, j = np.ravel(i), np.ravel(j)
+
+    """
+    Let's go for the triangles :
+        (i,j)    -  (i,j+1)   -> y dir
+        (i+1,j)  - (i+1,j+1)
+          |
+          v
+        x dir
+
+    in flatten coordinates:
+        i*ny+j     -  i*ny+j+1
+        (i+1)*ny+j -  (i+1)*ny+j+1
+    """
+
+    t1 = (i * ny + j,
+          (i + 1) % nx * ny + j,
+          (i + 1) % nx * ny + (j + 1) % ny)
+    t2 = (i * ny + j,
+          (i + 1) % nx * ny + (j + 1) % ny,
+          i * ny + (j + 1) % ny)
+
+    """
+        %nx and %ny are used for wrapx and wrapy :
+        if (i+1)=nx => (i+1)%nx=0 => close mesh in x direction
+        if (j+1)=ny => (j+1)%ny=0 => close mesh in y direction
+    """
+    
+    nt = len(t1[0])
+    
+    triangles = np.zeros((nt * 2, 3), dtype=np.uint32)
+    triangles[0::2, 0], triangles[0::2, 1], triangles[0::2, 2] = t1
+    triangles[1::2, 0], triangles[1::2, 1], triangles[1::2, 2] = t2
+    
+    lines = np.zeros((nt * 4, 2), dtype=np.uint32)
+    lines[::4,0], lines[::4,1] = t1[:2]
+    lines[1::4,0], lines[1::4,1] = t1[0],t2[2]
+    lines[2::4,0], lines[2::4,1] = t2[2:0:-1]
+    lines[3::4,0], lines[3::4,1] = t1[1],t2[1]   
+    
+    return triangles, lines
+

notebooks/Mesh speed comparison.ipynb

@@ -0,0 +1,200 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import ipyvolume.pylab as p3\n",
+    "from numpy.testing import assert_array_equal"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "from ipyvolume.pylab import make_triangles_lines as new_make_triangles_lines\n",
+    "def old_make_triangles_lines(x,y,z,wrapx=False,wrapy=False):\n",
+    "    def dim(x):\n",
+    "        d = 0\n",
+    "        el = x\n",
+    "        while True:\n",
+    "            try:\n",
+    "                el = el[0]\n",
+    "                d += 1\n",
+    "            except Exception:\n",
+    "                break\n",
+    "        return d\n",
+    "\n",
+    "    if dim(x) == 2:\n",
+    "        nx, ny = shape = x.shape\n",
+    "    else:\n",
+    "        nx, ny = shape = x[0].shape\n",
+    "\n",
+    "    def reshape(ar):\n",
+    "        if dim(ar) == 3:\n",
+    "            return [k.reshape(-1) for k in ar]\n",
+    "        else:\n",
+    "            return ar.reshape(-1)\n",
+    "\n",
+    "    x = reshape(x)\n",
+    "    y = reshape(y)\n",
+    "    z = reshape(z)\n",
+    "\n",
+    "    mx = nx if wrapx else nx - 1\n",
+    "    my = ny if wrapy else ny - 1\n",
+    "    triangles = np.zeros(((mx) * (my) * 2, 3), dtype=np.uint32)\n",
+    "    lines = np.zeros(((mx) * (my) * 4, 2), dtype=np.uint32)\n",
+    "    \n",
+    "    \n",
+    "      \n",
+    " \n",
+    "    def index_from2d(i, j):\n",
+    "        xi = (i % nx)\n",
+    "        yi = (j % ny)\n",
+    "        return ny * xi + yi\n",
+    "        \"\"\"\n",
+    "        ^ ydir\n",
+    "        |\n",
+    "        2 3\n",
+    "        0 1  ---> x dir\n",
+    "        \"\"\"\n",
+    " \n",
+    "    for i in range(mx):\n",
+    "        for j in range(my):\n",
+    "            p0 = index_from2d(i, j)\n",
+    "            p1 = index_from2d(i + 1, j)\n",
+    "            p2 = index_from2d(i, j + 1)\n",
+    "            p3 = index_from2d(i + 1, j + 1)\n",
+    "            triangle_index = (i * my) + j\n",
+    "            triangles[triangle_index * 2 + 0, :] = [p0, p1, p3]\n",
+    "            triangles[triangle_index * 2 + 1, :] = [p0, p3, p2]\n",
+    "            lines[triangle_index * 4 + 0, :] = [p0, p1]\n",
+    "            lines[triangle_index * 4 + 1, :] = [p0, p2]\n",
+    "            lines[triangle_index * 4 + 2, :] = [p2, p3]\n",
+    "            lines[triangle_index * 4 + 3, :] = [p1, p3]\n",
+    "    \n",
+    "    return triangles, lines\n",
+    "            "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "assert_array_equal(tri_new,tri_old)\n",
+    "assert_array_equal(line_new,line_old)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "330 µs ± 574 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
+      "44.9 ms ± 135 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Small data set.\n",
+    "nx = 100\n",
+    "ny = 50\n",
+    "X = np.linspace(-5, 5, nx)\n",
+    "Y = np.linspace(-5, 5, ny)\n",
+    "X, Y = np.meshgrid(X, Y, indexing='xy')\n",
+    "R = np.sqrt(X**2 + Y**2)\n",
+    "Z = np.sin(R)\n",
+    "\n",
+    "tri_new, line_new = new_make_triangles_lines(X,Y,Z)\n",
+    "tri_old, line_old = old_make_triangles_lines(X,Y,Z)\n",
+    "\n",
+    "assert_array_equal(tri_new,tri_old)\n",
+    "assert_array_equal(line_new,line_old)\n",
+    "\n",
+    "\n",
+    "%timeit new_make_triangles_lines(X,Y,Z)\n",
+    "%timeit old_make_triangles_lines(X,Y,Z)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "12.4 ms ± 80.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
+      "1.56 s ± 2.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Larger data set.\n",
+    "nx = 400\n",
+    "ny = 425\n",
+    "X = np.linspace(-5, 5, nx)\n",
+    "Y = np.linspace(-5, 5, ny)\n",
+    "X, Y = np.meshgrid(X, Y, indexing='xy')\n",
+    "R = np.sqrt(X**2 + Y**2)\n",
+    "Z = np.sin(R)\n",
+    "\n",
+    "tri_new, line_new = new_make_triangles_lines(X,Y,Z)\n",
+    "tri_old, line_old = old_make_triangles_lines(X,Y,Z)\n",
+    "\n",
+    "assert_array_equal(tri_new,tri_old)\n",
+    "assert_array_equal(line_new,line_old)\n",
+    "\n",
+    "%timeit tri_new, line_new = new_make_triangles_lines(X,Y,Z)\n",
+    "%timeit tri_old, line_old = old_make_triangles_lines(X,Y,Z)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}