From a923744d6d754eea314a5405c996cfd702a0733b Mon Sep 17 00:00:00 2001 From: xavArtley Date: Sat, 25 Nov 2017 10:09:51 +0100 Subject: [PATCH] Greatly improve mesh grid generation : more than 100x, look at mesh speed comparison notebook --- .gitignore | 7 + ipyvolume/pylab.py | 115 +++++++++++---- notebooks/Mesh speed comparison.ipynb | 200 ++++++++++++++++++++++++++ 3 files changed, 293 insertions(+), 29 deletions(-) create mode 100644 notebooks/Mesh speed comparison.ipynb diff --git a/.gitignore b/.gitignore index 32e06e5f..6006d82e 100644 --- a/.gitignore +++ b/.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 diff --git a/ipyvolume/pylab.py b/ipyvolume/pylab.py index b4b03c64..5a00faf7 100644 --- a/ipyvolume/pylab.py +++ b/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 + diff --git a/notebooks/Mesh speed comparison.ipynb b/notebooks/Mesh speed comparison.ipynb new file mode 100644 index 00000000..2554b1e3 --- /dev/null +++ b/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 +}