ld2daps/Notebooks/HVR Interpolations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ~~FAILURE~~ Difficulties"
   ]
  },
  {
   "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": [
    "## 3D Interpolations\n",
    "\n",
    "### Nearest neighbour\n",
    "\n",
    "Nearest neighbour from griddata works directly with data $\\in \\!R^3$.\n",
    "\n",
    "### Linear and cubic Interpolations\n",
    "\n",
    "We tried to interpolate (linear) 3D data in a regular 3D grid with no success so far. We tried:\n",
    "\n",
    "- `scipy.interpolate.griddata` is taking forever\n",
    "- `scipy.ndimage.interpolation.map_coordinates`\n",
    "- `scipy.interpolate.RegularGridInterpolator`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage.interpolation import map_coordinates\n",
    "from scipy.interpolate import RegularGridInterpolator\n",
    "\n",
    "\n",
    "map_coordinates(np.ma.masked_invalid(voxels[:,:,0]).data, np.nonzero(fnan[:,:,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import RegularGridInterpolator\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(1, 4, 11)\n",
    "y = np.linspace(4, 7, 22)\n",
    "z = np.linspace(7, 9, 33)\n",
    "x, y, z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.meshgrid(x, y, z, indexing='ij', sparse=True)[2].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi = RegularGridInterpolator((C1.x, C1.y, C1.z), C1.intensity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = voksel[:,:,:2]\n",
    "x = np.arange(0, raster.shape[1])\n",
    "y = np.arange(0, raster.shape[0])\n",
    "z = np.arange(0, raster.shape[2])\n",
    "\n",
    "raster = np.ma.masked_invalid(raster)\n",
    "xx, yy, zz = np.meshgrid(x, y, z)\n",
    "#get only the valid values\n",
    "x1 = xx[~raster.mask]\n",
    "y1 = yy[~raster.mask]\n",
    "z1 = zz[~raster.mask]\n",
    "\n",
    "newraster = raster[~raster.mask]\n",
    "\n",
    "gdr = griddata((x1, y1, z1), newraster.ravel(), (xx, yy, zz), method='linear')"
   ]
  },
  {
   "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": [
    "step = 1\n",
    "layers = gdr.shape[2]\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[:,:,i])\n",
    "    plt.imsave('../Res/voxel_layer_{}i.png'.format(i), voksel_i[:,:,i])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voksel_il = voksel.copy()\n",
    "voksel_il = r.interpolate(voksel_il, 'linear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step = 1\n",
    "layers = voksel_il.shape[2]\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[:,:,i])\n",
    "    plt.imsave('../Res/voxel_layer_{}il.png'.format(i), voksel_il[:,:,i])\n",
    "plt.show()"
   ]
  }
 ],
 "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
}