visualizing_point_cloud_mayavi.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:bd1fb351343a9526a9979147a28e303466c3b5713251f934ceef22152c216cb8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Visualizing point clouds with Mayavi"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook will cover the basics of visualizing point clouds with Mayavi. Mayavi is a python codeable interface for VTK, which offers the possibility to plot scientific data as well as point clouds."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import mayavi.mlab as mlab\n",
      "from tvtk.api import tvtk\n",
      "from mayavi.modules.glyph import Glyph\n",
      "from mayavi.tools.sources import MGlyphSource\n",
      "from mayavi.modules.surface import Surface\n",
      "from mayavi.filters.mask_points import MaskPoints"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Plotting the point cloud"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting the point cloud can be done very easy by using the mlab.mplot3d() command. It is that way however not easy to change the colors and add different modules to the pipeline. So this pipeline is what we will build manually. A pipeline is, quoting wikipedias definition: \n",
      "\n",
      "In 3D computer graphics, the graphics pipeline or rendering pipeline refers to the sequence of steps used to create a 2D raster representation of a 3D scene.\n",
      "\n",
      "So this is what we will do. But first import the data and create a data source:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xyz = np.load(r'D:\\rongen\\xyz_delft_center.npy', mmap_mode = 'r')\n",
      "rgb = np.load(r'D:\\rongen\\rgb_delft_center.npy', mmap_mode = 'r')\n",
      "\n",
      "def create_data_source():\n",
      "    '''\n",
      "    Function to create a data source readible by mayavi/vtk\n",
      "    '''\n",
      "    # Create the datasource\n",
      "    data_source  = MGlyphSource()\n",
      "    # Add the x, y and z data to it\n",
      "    data_source.reset(x = xyz[:, 0], y = xyz[:, 1], z = xyz[:, 2])\n",
      "    # Add a vtk_data_source to the pipeline\n",
      "    vtk_data_source = mlab.pipeline.add_dataset(data_source.dataset)\n",
      "    # To be able to change the content of the datasource afterwards, we state the following:\n",
      "    data_source.m_data = vtk_data_source\n",
      "    return vtk_data_source\n",
      "vtk_data_source = create_data_source()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "C:\\Users\\rongen\\AppData\\Local\\Continuum\\Anaconda\\lib\\site-packages\\traits\\has_traits.py:1536: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.\n",
        "  setattr( self, name, value )\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The window has poppud up already, but we still need to do some steps before we can see the point cloud.\n",
      "\n",
      "Add the colors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def add_colors(vtk_data_source):\n",
      "    '''\n",
      "    Function to add colors to data source\n",
      "    Input : the data source which has to be colored\n",
      "    '''\n",
      "    # Create an array which vtk can read\n",
      "    col = tvtk.UnsignedCharArray()\n",
      "    # Set color data to it\n",
      "    col.from_array(rgb)\n",
      "    # Change the scalar (color) data of the dataset to the color array\n",
      "    vtk_data_source.mlab_source.dataset.point_data.scalars=col\n",
      "\n",
      "add_colors(vtk_data_source)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now it is possible to add a mask to the pipeline. The mask will prevent some points from being rendered, which can come in handy when you do not want to show all points, to save for example memory. In our case, if you have sufficient memory (> 4 GB) you can probably plot all points and skip the mask. In that case you can add the glyph directly to the vtk_data_source (next step)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create an engine to which filters can be added\n",
      "engine = mlab.get_engine()\n",
      "\n",
      "def create_mask(npoints = 100000):\n",
      "    '''\n",
      "    Function to create a mask\n",
      "    Input : the number of points which are visible\n",
      "    '''\n",
      "    # Create the mask\n",
      "    mask = MaskPoints()\n",
      "    engine.add_filter(mask, vtk_data_source)\n",
      "    mask.filter.random_mode = True\n",
      "    mask.filter.maximum_number_of_points = npoints\n",
      "    return mask\n",
      "mask = create_mask()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And add a filter which contains information on how to visualize the points."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create the glyph with which the point cloud actually gets points\n",
      "def create_glyph(add_to = 'mask'):\n",
      "    '''\n",
      "    Function to add glyphs to the data\n",
      "    '''\n",
      "    glyph = Glyph()\n",
      "    if add_to == 'mask':\n",
      "        engine.add_filter(glyph, mask)\n",
      "    else:\n",
      "        engine.add_filter(glyph, vtk_data_source)\n",
      "    glyph.glyph.glyph_source.glyph_source.glyph_type = 'vertex'\n",
      "\n",
      "glyph = create_glyph()\n",
      "# Show the figure\n",
      "mlab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Improve navigation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The point cloud looks fine, but the navigation is still a bit annoying. The camera is rolling, which we do not want in our case. Also the focal point is fixed, so we can only look at one point. Lets improve these things, starting with the rolling."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#-------------------------------------------------------------------------------\n",
      "# Add a terrain interactor style, which does not roll. For this we use the tvtk\n",
      "# module which has vtk objects as traits.\n",
      "\n",
      "# The interactor style had to be inserted before the mask or glyph is added to \n",
      "# the pipeline, so we create the animation again below (Close it first)\n",
      "\n",
      "# Set the interactor to the new style\n",
      "def change_interaction():\n",
      "    '''\n",
      "    Function to change the interaction style\n",
      "    '''\n",
      "    istyle = tvtk.InteractorStyleTerrain()\n",
      "    iactor = fig.scene.interactor\n",
      "    iactor.interactor_style = istyle\n",
      "\n",
      "# figure\n",
      "fig             = mlab.gcf()\n",
      "\n",
      "# data_source\n",
      "vtk_data_source = create_data_source()\n",
      "add_colors(vtk_data_source)\n",
      "\n",
      "# interactorstyle\n",
      "change_interaction()\n",
      "\n",
      "# engine, mask and glyph\n",
      "engine          = mlab.get_engine()\n",
      "mask            = create_mask(200000)\n",
      "glyph           = create_glyph()\n",
      "\n",
      "# show\n",
      "mlab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The camera stopped rolling, but the zoom by scroll function has been disabled for some reason.. Let's build that one ourselves, while we our builing in key press event for moving around."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We will have to add observers for key press events and mouse wheel events.\n",
      "# Change the interaction function and add observers\n",
      "\n",
      "def change_interaction():\n",
      "    '''\n",
      "    Function to change the interaction style and adds observers\n",
      "    '''\n",
      "    istyle = tvtk.InteractorStyleTerrain()\n",
      "    iactor = fig.scene.interactor\n",
      "    iactor.interactor_style = istyle\n",
      "    istyle.add_observer('KeyPressEvent',           Callback)\n",
      "    istyle.add_observer('MouseWheelBackwardEvent', Callback)\n",
      "    istyle.add_observer('MouseWheelForwardEvent',  Callback)\n",
      "\n",
      "# A callback function is needed to handle the events. The abbsorvers give two\n",
      "# arguments to the callback function: the object, with all the event traits\n",
      "# and the eventname itself, which is just a string.\n",
      "\n",
      "def Callback(obj, event):\n",
      "    '''\n",
      "    Callback function, handels the interaction events\n",
      "    '''\n",
      "    # Get the geometric properties (azimuth, elevation, distance, focal point)\n",
      "    # This van easily be done be calling mlav.view() without arguments.\n",
      "    # mlab.move() without argument gives the camera position and the focal point\n",
      "    azm, elev, dist, fp = mlab.view()\n",
      "    cam, fp = mlab.move()\n",
      "    # Set the height of the focal point to 2.5 m.\n",
      "    fp[2] = 2.5\n",
      "    # Convert to radians\n",
      "    azm = azm * 2 * np.pi / 360.\n",
      "    \n",
      "    # To move the camera we will use the numpad keys.\n",
      "    if event == 'KeyPressEvent':\n",
      "        # Get the pressed key\n",
      "        key = obj.GetInteractor().GetKeyCode()\n",
      "        \n",
      "        # Let the moved distance be determined by the hieght of the camera (cam[2])\n",
      "        step = cam[2]/50.+1\n",
      "        \n",
      "        if key == '8':\n",
      "            mlab.view(focalpoint = fp + [-step*np.cos(azm),-step*np.sin(azm), 0], distance = dist)\n",
      "        if key == '4':\n",
      "            mlab.view(focalpoint = fp + [ step*np.sin(azm),-step*np.cos(azm), 0], distance = dist)\n",
      "        if key == '5':\n",
      "            mlab.view(focalpoint = fp + [ step*np.cos(azm), step*np.sin(azm), 0], distance = dist)\n",
      "        if key == '6':\n",
      "            mlab.view(focalpoint = fp + [-step*np.sin(azm), step*np.cos(azm), 0], distance = dist)\n",
      "    \n",
      "    # Make our own zoom function\n",
      "    if event == 'MouseWheelForwardEvent':\n",
      "        dist *= 0.9\n",
      "        mlab.view(focalpoint = fp, distance = dist)\n",
      "    if event == 'MouseWheelBackwardEvent':\n",
      "        dist *= 1.111\n",
      "        mlab.view(focalpoint = fp, distance = dist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig             = mlab.gcf()\n",
      "vtk_data_source = create_data_source()\n",
      "add_colors(vtk_data_source)\n",
      "change_interaction()\n",
      "engine          = mlab.get_engine()\n",
      "mask            = create_mask(1000000)\n",
      "glyph           = create_glyph()\n",
      "mlab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That works fine! So we now have a detailed point cloud, which interacts fine."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Making a movie"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next thing we are going to do is create a sequence of png snapshots, which we can compile to a movie. We will add another key to our interaction function, which will rotate the camera and take a snapshot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import time\n",
      "\n",
      "def Callback(obj, event):\n",
      "    '''\n",
      "    Callback function, handels the interaction events\n",
      "    '''\n",
      "    # Get the geometric properties (azimuth, elevation, distance, focal point)\n",
      "    # This van easily be done be calling mlav.view() without arguments.\n",
      "    # mlab.move() without argument gives the camera position and the focal point\n",
      "    azm, elev, dist, fp = mlab.view()\n",
      "    cam, fp = mlab.move()\n",
      "    # Set the height of the focal point to 2.5 m.\n",
      "    fp[2] = 2.5\n",
      "    # Convert to radians\n",
      "    azm = azm * 2 * np.pi / 360.\n",
      "    \n",
      "    # To move the camera we will use the numpad keys.\n",
      "    if event == 'KeyPressEvent':\n",
      "        # Get the pressed key\n",
      "        key = obj.GetInteractor().GetKeyCode()\n",
      "        \n",
      "        # Let the moved distance be determined by the hieght of the camera (cam[2])\n",
      "        step = cam[2]/50.+1\n",
      "        \n",
      "        if key == '8':\n",
      "            mlab.view(focalpoint = fp + [-step*np.cos(azm),-step*np.sin(azm), 0], distance = dist)\n",
      "        if key == '4':\n",
      "            mlab.view(focalpoint = fp + [ step*np.sin(azm),-step*np.cos(azm), 0], distance = dist)\n",
      "        if key == '5':\n",
      "            mlab.view(focalpoint = fp + [ step*np.cos(azm), step*np.sin(azm), 0], distance = dist)\n",
      "        if key == '6':\n",
      "            mlab.view(focalpoint = fp + [-step*np.sin(azm), step*np.cos(azm), 0], distance = dist)\n",
      "    \n",
      "        # Make a Movie\n",
      "        if key == 'm':\n",
      "            # Enable off screen rendering, which shortens the time to create snapshots greatly\n",
      "            engine.current_scene.scene.off_screen_rendering = True\n",
      "            # Loop through different azimuths (rotations)\n",
      "            i = 0\n",
      "            for azimuth in range(0,360,30):\n",
      "                # Set the view:\n",
      "                mlab.view(azimuth = azimuth, elevation = 75, distance = 500, focalpoint = np.array([84465, 447547, 2.5]))\n",
      "                # Save the scene\n",
      "                engine.scenes[0].scene.save(r'D:\\rongen\\snapshot{:03d}.png'.format(i))\n",
      "                i += 1\n",
      "                # Pause for a second, to avoid break downs\n",
      "                time.sleep(1)\n",
      "            # Disable off screen rendering again\n",
      "            engine.current_scene.scene.off_screen_rendering = False\n",
      "    \n",
      "    # Make our own zoom function\n",
      "    if event == 'MouseWheelForwardEvent':\n",
      "        dist *= 0.9\n",
      "        mlab.view(focalpoint = fp, distance = dist)\n",
      "    if event == 'MouseWheelBackwardEvent':\n",
      "        dist *= 1.111\n",
      "        mlab.view(focalpoint = fp, distance = dist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Pressing m will create a sequence of snapshots from different angles now. These snapshots will have to be combined with a movie maker program or ffmpeg or whatever you like."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig             = mlab.gcf()\n",
      "vtk_data_source = create_data_source()\n",
      "add_colors(vtk_data_source)\n",
      "change_interaction()\n",
      "engine          = mlab.get_engine()\n",
      "mask            = create_mask(100000)\n",
      "glyph           = create_glyph()\n",
      "mlab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Adding model data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A way to add modeldata in Mayavi, is to create a tvtk unstructured grid, and add the data to it. The next section will explain how to do this.\n",
      "\n",
      "Note that a lot of data files have a different are build up differently. It is therefor hard to write a script which is applicable for all grids. The next cell will therefore not work, but should be applied to your specific dataset."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The trick to make an tvtk.Unstructuredgrid is to make an administration of the\n",
      "# cells in the grid. Per cell should be known which grid points surround the\n",
      "# cell. An index of each node per cell is needed.\n",
      "\n",
      "# Let's say we have a netcdf file with the following data:\n",
      "xnode = dataset.variables['FlowElemContour_x'][:] # (N,4)-array with the x-coordinates of the cell corners\n",
      "ynode = dataset.variables['FlowElemContour_y'][:] # (N,4)-array with the y-coordinates of the cell corners\n",
      "\n",
      "xcc = dataset.variables['FlowElem_xcc'][:] # (N,1)-array with the x-coordinates of the cell centers\n",
      "ycc = dataset.variables['FlowElem_ycc'][:] # (N,1)-array with the x-coordinates of the cell centers\n",
      "\n",
      "# First we need to create an array with points.\n",
      "points = np.c_[xnode.ravel(), ynode.ravel()]\n",
      "\n",
      "# We only need every point once, so sort out the unqie points:\n",
      "b = np.ascontiguousarray(points).view(np.dtype((np.void, points.dtype.itemsize * points.shape[1])))\n",
      "_, idx = np.unique(b, return_index=True)\n",
      "unique_points = points[idx]\n",
      "\n",
      "# For some reason creating a tvtk.UnstructuredGrid goes wrong when point 0 is\n",
      "# requested. Therefore we add a dummy point, which will not be requested\n",
      "unique_points = np.vstack((unique_points[0,:]+0.5, unique_points))\n",
      "\n",
      "# Write a function to check where a point is in the unique_points array\n",
      "def checkwhere(elm):\n",
      "    c1 = np.where(elm[0] == unique_points[:,0])[0]\n",
      "    c2 = np.where(elm[1] == unique_points[:,1])[0]\n",
      "    return np.intersect1d(c1, c2)[0]\n",
      "\n",
      "# Create the indices of the cell corners for each cell\n",
      "indices = np.array([checkwhere(i) for i in np.c_[xnode.ravel(), ynode.ravel()]])\n",
      "indices = indices.reshape((len(indices)/4,4))\n",
      "\n",
      "# Number of cells\n",
      "n_cells = int(indices.shape[0])\n",
      "\n",
      "# Array with number of nodes per cell\n",
      "elemtype = np.sum(indices > 0, 1)\n",
      "\n",
      "# Array with cell types, these depend on the number of nodes per cell.\n",
      "# if these are not four per cell, script something to get the right cell\n",
      "# type for each quantity of nodes\n",
      "cell_types = np.ones(len(elemtype))*tvtk.Quad().cell_type\n",
      "\n",
      "# Array with number of nodes per cell, and index of nodes\n",
      "cells = np.c_[elemtype, indices].ravel()\n",
      "\n",
      "# Array with offset, needed to read out cells\n",
      "offsets = np.r_[[0], np.cumsum(elemtype + 1)]\n",
      "\n",
      "# Convert cells to tvtk CellArray\n",
      "cell_array = tvtk.CellArray()\n",
      "cell_array.set_cells(n_cells, cells)\n",
      "\n",
      "# Create the grid\n",
      "grid = tvtk.UnstructuredGrid(points = unique_points)\n",
      "grid.set_cells(cell_types, offsets, cell_array)\n",
      "\n",
      "# Now scalar values can be set to the cell centers or cell npoints (nodes)\n",
      "# If you have cell data available you might need to interpolate to \n",
      "# point data, or vice versa\n",
      "grid.point_data.scalars = your_point_data\n",
      "grid.point_data.scalars.name = 'my_data'\n",
      "grid.cell_data.scalars = your_point_data\n",
      "grid.cell_data.scalars.name = 'my_data'\n",
      "\n",
      "# Save the grid to an xml file\n",
      "w = tvtk.XMLUnstructuredGridWriter(input = grid, file_name = 'my_grid.xml')\n",
      "w.write()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "When you have a grid created or saved, you can add it to Mayavi as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the grid data if it is not stored in memory\n",
      "read = tvtk.XMLUnstructuredGridReader(file_name = 'path\\to\\xml\\file.xml')\n",
      "grid = read.get_output()\n",
      "\n",
      "# Show the grid:\n",
      "fig = mlab.gcf()\n",
      "engine = mlab.get_engine()\n",
      "mlab.pipeline.add_dataset(grid)\n",
      "engine.add_module(Surface())\n",
      "mlab.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In Mayavi, a lot of properties can be changed for the appearance of your grid. This can be done by clicking the Mayavi icon, which opens the pipeline. A pro tip is to click the recording button. The code of the changes yo do after clicking the recording button will be shown in the recording window. This makes it easy to find out how to change certain properties of your grids in your script directly."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}