{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "import matplotlib.pyplot as plt\n",
    "import laspy\n",
    "\n",
    "triskele_path = Path('../triskele/python/')\n",
    "sys.path.append(str(triskele_path.resolve()))\n",
    "import triskele\n",
    "\n",
    "figsize = np.array((16, 3)) * 1.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LiDAR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_las(path):\n",
    "    path = Path(path)\n",
    "    las = list()\n",
    "    \n",
    "    if path.is_dir():\n",
    "        for child in path.iterdir():\n",
    "            las.extend(find_las(child))\n",
    "    \n",
    "    if path.is_file() and path.suffix == '.las':\n",
    "        las.append(path)\n",
    "    \n",
    "    return las"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def auto_filter(data, treshold=.5):\n",
    "    tresholds = np.percentile(data, [treshold, 100 - treshold])\n",
    "\n",
    "    data[data < tresholds[0]] = tresholds[0]\n",
    "    data[data > tresholds[1]] = tresholds[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bulk_load(path, name=None):\n",
    "    data = {'file': path, 'name': name}\n",
    "\n",
    "    attributes = ['x', 'y', 'z', 'intensity', 'num_returns']#, 'scan_angle_rank', 'pt_src_id', 'gps_time']\n",
    "    for a in attributes:\n",
    "        data[a] = list()\n",
    "        \n",
    "    print('Load data...')\n",
    "    for f in find_las(path):\n",
    "        print('{}: '.format(f), end='')\n",
    "        infile = laspy.file.File(f)\n",
    "        for i, a in enumerate(attributes):\n",
    "            print('\\r {}: [{:3d}%]'.format(f, int(i/len(attributes) * 100)), end='')\n",
    "            data[a].extend(getattr(infile, a))\n",
    "        infile.close()\n",
    "        print('\\r {}: [Done]'.format(f))\n",
    "\n",
    "    \n",
    "    print('Create matrices...', end='')\n",
    "    for i, a in enumerate(attributes):\n",
    "        print('\\rCreate matrices: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
    "        data[a] = np.array(data[a])\n",
    "    print('\\rCreate matrices: [Done]')\n",
    "        \n",
    "    print('Filter data...', end='')\n",
    "    for i, a in enumerate(attributes):\n",
    "        print('\\rFilter data: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
    "        if a == 'x' or a == 'y':\n",
    "            continue\n",
    "        auto_filter(data[a])\n",
    "    print('\\rFilter data: [Done]')\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = bulk_load('../Data/lidar/C1/272056_3289689.las', 'C1')\n",
    "#C2 = bulk_load('../Data/lidar/C2/', 'C2')\n",
    "#C3 = bulk_load('../Data/lidar/C3/', 'C3')\n",
    "#C123 = bulk_load('../Data/lidar/', 'C123') # Todo: Copy C1 C2 C3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define raster method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 0.5\n",
    "\n",
    "spatial_data = np.array((C1['x'], C1['y'], C1['z'])).T\n",
    "\n",
    "bins = np.rint((spatial_data.max(axis=0) - spatial_data.min(axis=0)) / resolution).astype(int)\n",
    "display(bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxels_histo, edges = np.histogramdd(spatial_data, bins=bins)\n",
    "voxels_histo = voxels_histo.astype(np.float32)\n",
    "voxels_histo.shape, voxels_histo.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxels_histo_i, edges = np.histogramdd(spatial_data, bins=bins, weights=C1['intensity'])\n",
    "voxels_histo_i = voxels_histo_i.astype(np.float32)\n",
    "voxels_histo_i.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layers = voxels_histo.shape[2]\n",
    "\n",
    "plt.figure(figsize=figsize * np.array((1, layers)))\n",
    "\n",
    "for i in range(layers):\n",
    "    plt.subplot(layers, 1, i+1)\n",
    "    plt.title('Layer {}'.format(i))\n",
    "    plt.imshow(voxels_histo_i.T[i] / voxels_histo.T[i], origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/voxel_layer', voxels_histo.T[0], origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create voxels and flatten raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layers = voxels_histo.shape[2]\n",
    "\n",
    "raster = np.zeros_like(voxels_histo.T[0])\n",
    "raster[:] = np.nan\n",
    "\n",
    "voxels = (voxels_histo_i / voxels_histo).T\n",
    "\n",
    "for i in reversed(range(layers)):\n",
    "    nan_filter = np.isnan(raster)\n",
    "    raster[nan_filter] = voxels[i][nan_filter]\n",
    "\n",
    "raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del voxels_histo, voxels_histo_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster, origin='lower')\n",
    "plt.show()\n",
    "plt.imsave('../Res/voxel_raster_i_0_5.png', raster, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolate NoData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "griddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nan_filter = np.isnan(raster)\n",
    "val_filter = np.logical_not(nan_filter)\n",
    "coords = np.argwhere(val_filter)\n",
    "values = raster[val_filter]\n",
    "grid   = np.argwhere(nan_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpol_values_linear  = griddata(coords, values, grid, method='linear')\n",
    "interpol_values_cubic   = griddata(coords, values, grid, method='cubic')\n",
    "interpol_values_nearest = griddata(coords, values, grid, method='nearest')\n",
    "\n",
    "raster_interpol_linear  = raster.copy()\n",
    "raster_interpol_cubic   = raster.copy()\n",
    "raster_interpol_nearest = raster.copy()\n",
    "\n",
    "raster_interpol_linear[nan_filter]  = interpol_values_linear\n",
    "raster_interpol_cubic[nan_filter]   = interpol_values_cubic\n",
    "raster_interpol_nearest[nan_filter] = interpol_values_nearest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_interpol_linear, origin='lower')\n",
    "plt.show()\n",
    "plt.imsave('../Res/voxel_raster_interpol_tmp_i_0_5.png', raster_interpol_linear, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why do we have high values there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_interpol_nearest.max(), raster_interpol_linear.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_interpol_nearest.min(), raster_interpol_linear.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_interpol_linear[raster_interpol_linear > 1e6] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((raster_interpol_linear - raster_interpol_nearest))\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "plt.imsave('../Res/diff_interpol.png', raster_interpol_linear - raster_interpol_nearest, origin='lower')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(raster_interpol_linear - raster_interpol_nearest).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_interpol_cubic[np.isnan(raster_interpol_cubic)] = np.nanargmin(raster_interpol_cubic)\n",
    "raster_interpol_linear[np.isnan(raster_interpol_linear)] = np.nanargmin(raster_interpol_linear)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.hist(raster_interpol_cubic.reshape(-1), bins=100, label='cubic',log=True)\n",
    "plt.hist(raster_interpol_linear[raster_interpol_linear < 1e6].reshape(-1), 100, label='linear',log=False)\n",
    "plt.hist(raster_interpol_nearest.reshape(-1), 100, label='nearest',log=False, alpha=.8)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isnan(raster_interpol_nearest).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "sys.path.append(\"..\")\n",
    "import rasterizer\n",
    "importlib.reload(rasterizer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myraster = rasterizer.rasterize(spatial_data, C1['intensity'], 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myraster_r = rasterizer.rasterize(spatial_data, C1['intensity'], 1., reverse_alt=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/myrstr.png', myraster, origin='upper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/myrstr_r.png', myraster_r, origin='upper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mid_layer = myraster - myraster_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/myrstr_layer_1.png', mid_layer, origin='upper')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HVR II\n",
    "\n",
    "We will exploit the 3D to interpolate no data a much more precise way. What we need is to record empty cells where the laser do not encounter objects when previously we let no data value. This will allow better understanding of the scene before applying interpolation strategies in 3D.\n",
    "\n",
    "## Shot vector\n",
    "\n",
    "Find a way to determine the laser path in the scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "infile = laspy.file.File('../Data/lidar/C1/272056_3289689.las')\n",
    "for a in infile.point_format:\n",
    "    print(a.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "infile.header.version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each point we have its coordinate $x, y$ and $z$ and the angle `scan_angle_rank` (signed and rounded to the nearest integer) at which the laser was output :\n",
    "\n",
    "- 0 is nadir\n",
    "- -90 is the left of the aircraft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(infile.scan_angle_rank, 100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have:\n",
    "\n",
    "- Point coordinates\n",
    "- $\\alpha$ (roll) angle\n",
    "\n",
    "What we don't have:\n",
    "\n",
    "- $\\beta$ (pitch) angle : fixed angle of the captor + pitch of the aircraft\n",
    "- $\\gamma$ (yaw) angle : yaw of the aircraft, could be derived from LiDAR point cloud coordinates and gps time ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assert nadir vectors\n",
    "\n",
    "### Create voxels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(\"..\")\n",
    "import raster_assistant as ra\n",
    "import rasterizer as r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=0, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 1.\n",
    "bins = np.rint((C1.spatial.max(axis=0) - C1.spatial.min(axis=0)) / resolution).astype(int)\n",
    "bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxels_histo, edges = np.histogramdd(C1.spatial, bins=bins)\n",
    "voxels_histo_w, edges = np.histogramdd(C1.spatial, bins=bins, weights=C1.intensity)\n",
    "\n",
    "voxels = (voxels_histo_w / voxels_histo).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxels.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prototype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voksel = voxels.copy()\n",
    "\n",
    "for z, y, x in zip(*np.nonzero(np.logical_not(np.isnan(voxels)))):\n",
    "    voksel[z:,y,x] = 0.\n",
    "    \n",
    "f = ~np.isnan(voxels)\n",
    "voksel[f] = voxels[f]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step = 7\n",
    "layers = voksel.shape[0]\n",
    "\n",
    "plt.figure(figsize=figsize * np.array((1, layers)))\n",
    "\n",
    "for i in range(0, layers, step):\n",
    "    plt.subplot(layers, 1, i / step + 1)\n",
    "    plt.title('Layer {}'.format(i))\n",
    "    plt.imshow(voksel[i])\n",
    "    plt.imsave('../Res/voxel_layer_{}.png'.format(i), voksel[i])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voksel_i = voksel.copy()\n",
    "voksel_i = r.interpolate(voksel_i, 'nearest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voksel_i_disp = voksel_i.copy()\n",
    "voksel_i_disp[voksel_i_disp > 60e3] = 60e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step = 7\n",
    "layers = voksel_i.shape[0]\n",
    "\n",
    "plt.figure(figsize=figsize * np.array((1, layers)))\n",
    "\n",
    "for i in range(0, layers, step):\n",
    "    plt.subplot(layers, 1, i / step + 1)\n",
    "    plt.title('Layer {}'.format(i))\n",
    "    plt.imshow(voksel_i_disp[i])\n",
    "    plt.imsave('../Res/voxels_layers/voxel_layer_{}i.png'.format(layers - i), voksel_i_disp[i])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layers = voksel_i.shape[0]\n",
    "\n",
    "raster = np.zeros_like(voksel_i[0])\n",
    "\n",
    "for i in reversed(range(layers)):\n",
    "    nan_filter = raster == 0\n",
    "    raster[nan_filter] = voksel_i[i][nan_filter]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(raster)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow((raster == 32.) * 1.)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/raster-vhr2_1.png', raster[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "N1 = 10\n",
    "N2 = 10\n",
    "N3 = 10\n",
    "ma = np.random.choice([0,1.,2.], size=(N1,N2,N3), p=[0.98, 0.01, 0.01])\n",
    "ma = ma[:,:,:,np.newaxis]\n",
    "display(ma.shape)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.voxels(ma != [0], edgecolor=\"k\", facecolors=ma[ma!=0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voksel_i.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write voxels in LAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ofile = laspy.file.File('../Res/voxels.las', mode='w', header=infile.header)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z, y, x = np.indices(voksel_i.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = voksel_i != 0\n",
    "f.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[f], y[f], z[f], voksel_i[f]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TODO**: Try meshgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ofile.x = x[f] * resolution\n",
    "ofile.y = y[f] * resolution\n",
    "ofile.z = z[f] * resolution\n",
    "ofile.intensity = voksel_i[f]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ofile.points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ofile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ma == [0.]).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ma[...,0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl_vox = voksel_i != 0\n",
    "mpl_col = np.zeros(voksel_i.shape + (3,))\n",
    "mpl_col[...,0] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "xmin = 0; xmax = 10\n",
    "ymin = 0; ymax = 10\n",
    "zmin = 0; zmax = 10\n",
    "\n",
    "ix, iy, iz = np.indices(mpl_vox.shape)\n",
    "fx = (ix >= xmin) & (ix < xmax)\n",
    "fy = (iy >= ymin) & (iy < ymax)\n",
    "fz = (iz >= zmin) & (iz < zmax)\n",
    "f = fx & fy & fz\n",
    "shape = (xmax - xmin, ymax - ymin, zmax - zmin)\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.voxels(mpl_vox[f].reshape(shape), edgecolor=\"k\", facecolors=mpl_col[f].reshape(shape + (3,)))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}