Merge branch 'master' of http://git.normalized.xyz/Florent/LD2DAPs

Notebooks/HVR Interpolations.ipynb

@@ -0,0 +1,204 @@
+{
+ "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
+}

Notebooks/Histogram Voxels Rasterisation.ipynb

@@ -389,7 +389,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Function "
+    "## Test Function "
    ]
   },
   {
@@ -462,56 +462,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# TMP"
+    "# HVR II\n",
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "dfc_c1 = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n",
-    "auto_filter(dfc_c1)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.imshow(dfc_c1)\n",
-    "plt.show()\n",
-    "plt.imsave('../Res/dfc_c1.png', dfc_c1)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def rasterize_cache(spatial, values, resolution=1., method='linear', reverse_alt=False, cache_dir='/tmp'):\n",
-    "    \"\"\"Cache layer for rasterize\"\"\"\n",
-    "    \n",
-    "    cache_dir = Path(cache_dir)\n",
-    "    name = '{}_{}_{}_{}{}'.format(data['name'], field, str(resolution).replace('.', '_'), method,\n",
-    "                                  '_reversed' if reverse_alt else '')\n",
-    "    png_file = cache_dir.joinpath(Path(name + '.png'))\n",
-    "    tif_file = cache_dir.joinpath(Path(name + '.tif'))\n",
-    "    \n",
-    "    if tif_file.exists() :\n",
-    "        print ('WARNING: Loading cached result {}'.format(tif_file))\n",
-    "        raster = triskele.read(tif_file)\n",
-    "    else:         \n",
-    "        raster = rasterizer.rasterize(spatial, values, resolution, method, reverse_alt)\n",
-    "        triskele.write(tif_file, raster)\n",
-    "        plt.imsave(png_file, raster)\n",
-    "        \n",
-    "    return raster\n",
     "\n",
-    "rasterize_cache(spatial_data, C1['z'], 10, cache_dir='../Res/TMP')"
+    "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."
    ]
   },
   {
@@ -520,9 +477,420 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "td = Path(Path('/tmp'))\n",
+    "infile = laspy.file.File('../Data/lidar/C1/272056_3289689.las')\n",
-    "tf = Path('test.png')\n",
+    "for a in infile.point_format:\n",
-    "td.joinpath(tf).j"
+    "    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()"
    ]
   },
   {

Notebooks/Raster Factory.ipynb

@@ -40,10 +40,26 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=.5, dtype=np.float32)\n",
+    "C1 = ra.bulk_load('../Data/lidar/C1', 'C1', filter_treshold=0, dtype=np.float32)\n",
-    "C2 = ra.bulk_load('../Data/lidar/C2', 'C2', filter_treshold=.5, dtype=np.float32)\n",
+    "#C2 = ra.bulk_load('../Data/lidar/C2', 'C2', filter_treshold=.5, dtype=np.float32)\n",
-    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=.5, dtype=np.float32)\n",
+    "#C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=.5, dtype=np.float32)\n",
-    "C123 = ra.bulk_load('../Data/lidar', 'C123', filter_treshold=.5, dtype=np.float32)"
+    "#C123 = ra.bulk_load('../Data/lidar', 'C123', filter_treshold=.5, dtype=np.float32)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# TMP"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "raster = ra.rasterize_cache(C1, 'intensity', 1., 'nearest', True, cache_dir='../Res/')"
    ]
   },
   {