diff --git a/Notebooks/HVR Interpolations.ipynb b/Notebooks/HVR Interpolations.ipynb new file mode 100644 index 0000000..49cd88f --- /dev/null +++ b/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 +} diff --git a/Notebooks/Histogram Voxels Rasterisation.ipynb b/Notebooks/Histogram Voxels Rasterisation.ipynb index abb6432..fd6298c 100644 --- a/Notebooks/Histogram Voxels Rasterisation.ipynb +++ b/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" - ] - }, - { - "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", + "# HVR II\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", - "tf = Path('test.png')\n", - "td.joinpath(tf).j" + "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()" ] }, { diff --git a/Notebooks/None0000000.png b/Notebooks/None0000000.png deleted file mode 100644 index e69de29..0000000 diff --git a/Notebooks/Raster Factory.ipynb b/Notebooks/Raster Factory.ipynb index e26f22e..4596539 100644 --- a/Notebooks/Raster Factory.ipynb +++ b/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", - "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", - "C123 = ra.bulk_load('../Data/lidar', 'C123', filter_treshold=.5, dtype=np.float32)" + "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", + "#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)" + ] + }, + { + "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/')" ] }, {