florent/ld2daps
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
9e9dbb440e
ld2daps/Notebooks/Histogram Voxels Rasterisation.ipynb
Karamaz0V1 10aed66682 HVR Denoising
2018-07-05 12:21:05 +02:00

557 lines
14 KiB
Plaintext
Raw 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": [
"## 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": [
"# 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",
"\n",
"rasterize_cache(spatial_data, C1['z'], 10, cache_dir='../Res/TMP')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"td = Path(Path('/tmp'))\n",
"tf = Path('test.png')\n",
"td.joinpath(tf).j"
]
},
{
"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