florent/ld2daps
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
master
ld2daps/Notebooks/Histogram Voxels Rasterisation.ipynb
Karamaz0V1 ff8fc08fe3 HVR II
2018-07-24 17:28:04 +02:00

925 lines
21 KiB
Plaintext
Raw Permalink Blame History

{
"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
}
Reference in New Issue View Git Blame Copy Permalink