401 lines
9.6 KiB
Plaintext
401 lines
9.6 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Setup\n",
|
|
"\n",
|
|
"## Create legacy HVR"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import sys\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"sys.path.append(\"..\")\n",
|
|
"import rasterizer\n",
|
|
"import raster_assistant as ra\n",
|
|
"\n",
|
|
"sys.path.append('../triskele/python/')\n",
|
|
"import triskele\n",
|
|
"\n",
|
|
"figsize = np.array((16, 3)) * 1.5"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"data = ra.bulk_load('../Data/lidar/C3/', filter_treshold=0)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"raster = rasterizer.vhr(data.spatial, data.z, dtype=np.float32)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.figure(figsize=figsize)\n",
|
|
"plt.imshow(raster)\n",
|
|
"plt.show()\n",
|
|
"plt.imsave('../Res/tmp.png', raster)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Noise Profile"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.hist(data.z, bins=1000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.hist(data.intensity, bins=1000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.figure(figsize=(10,10))\n",
|
|
"plt.hist2d(data.intensity, data.z, bins=1000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"treshold = .5\n",
|
|
"tresholds_i = np.percentile(data.intensity, [treshold, 100 - treshold])\n",
|
|
"tresholds_z = np.percentile(data.z, [treshold, 100 - treshold])\n",
|
|
"\n",
|
|
"intensity_filter = np.logical_or(data.intensity < tresholds_i[0], data.intensity > tresholds_i[1])\n",
|
|
"z_filter = np.logical_or(data.z < tresholds_z[0], data.z > tresholds_z[1])\n",
|
|
"all_filter = np.logical_or(intensity_filter, z_filter)\n",
|
|
"inter_filter = np.logical_and(z_filter, intensity_filter)\n",
|
|
"\n",
|
|
"data.x.size, intensity_filter.sum(), z_filter.sum(), all_filter.sum(), inter_filter.sum()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The z and intensity noise seems uncorrelated."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Noise rasters\n",
|
|
"\n",
|
|
"### Denoise altitude for voxels computation..."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"Z = data.spatial[:,2]\n",
|
|
"ra.auto_filter(Z)\n",
|
|
"data.spatial = np.hstack((data.spatial[:,:2], Z.reshape(-1, 1)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Compute noise rasters"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"all_noise = rasterizer.rasterize(data.spatial, all_filter.astype(np.float16), .5, 'nearest', dtype=np.float16)\n",
|
|
"z_noise = rasterizer.rasterize(data.spatial, z_filter.astype(np.float16), .5, 'nearest', dtype=np.float16)\n",
|
|
"i_noise = rasterizer.rasterize(data.spatial, intensity_filter.astype(np.float16), .5, 'nearest', dtype=np.float16)\n",
|
|
"inter_noise = rasterizer.rasterize(data.spatial, inter_filter.astype(np.float16), .5, 'nearest', dtype=np.float16)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Display"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, axs = plt.subplots(4, figsize=figsize * 4)\n",
|
|
"for i, fig_title in enumerate(zip((all_noise, z_noise, i_noise, inter_noise), ('All', 'Altitude', 'Intensity', 'Intersection'))):\n",
|
|
" axs[i].imshow(fig_title[0].astype(np.float32))\n",
|
|
" axs[i].set_title(fig_title[1])\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Z and Int cutoff visualisation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, axs = plt.subplots(2, 2, figsize=(200,200))\n",
|
|
"for ax, filter_title in zip(axs.flatten(), zip((all_filter, z_filter, intensity_filter, inter_filter), ('All', 'Altitude', 'Intensity', 'Intersection'))):\n",
|
|
" ax.set_title(filter_title[1])\n",
|
|
" ax.hist2d(data.intensity[np.logical_not(filter_title[0])], data.z[np.logical_not(filter_title[0])], bins=1000)\n",
|
|
"fig.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Reconstruct denoised rasters\n",
|
|
"\n",
|
|
"#### All filter"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"raster_all = rasterizer.vhr(data.spatial[np.logical_not(all_filter)], data.z[np.logical_not(all_filter)], dtype=np.float32)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.figure(figsize=figsize)\n",
|
|
"plt.imshow(raster_all)\n",
|
|
"plt.show()\n",
|
|
"plt.imsave('../Res/tmp_all.png', raster_all)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Intensity filter"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"raster_i = rasterizer.vhr(data.spatial[np.logical_not(intensity_filter)], data.z[np.logical_not(intensity_filter)], dtype=np.float32)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.figure(figsize=figsize)\n",
|
|
"plt.imshow(raster_i)\n",
|
|
"plt.show()\n",
|
|
"plt.imsave('../Res/tmp_i.png', raster_i)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Z filter"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"raster_z = rasterizer.vhr(data.spatial[np.logical_not(z_filter)], data.z[np.logical_not(z_filter)], dtype=np.float32)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.figure(figsize=figsize)\n",
|
|
"plt.imshow(raster_z)\n",
|
|
"plt.show()\n",
|
|
"plt.imsave('../Res/tmp_z.png', raster_z)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Extremum clipping"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def extremum_filter(data, treshold=.5):\n",
|
|
" tresholds_i = np.percentile(data.intensity, [treshold, 100 - treshold])\n",
|
|
" tresholds_z = np.percentile(data.z, [treshold, 100 - treshold])\n",
|
|
"\n",
|
|
" intensity_filter = np.logical_or(data.intensity < tresholds_i[0], data.intensity > tresholds_i[1])\n",
|
|
" z_filter = np.logical_or(data.z < tresholds_z[0], data.z > tresholds_z[1])\n",
|
|
" all_filter = np.logical_or(intensity_filter, z_filter)\n",
|
|
" #inter_filter = np.logical_and(z_filter, intensity_filter)\n",
|
|
"\n",
|
|
" return all_filter"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"tresholds = (0., .1, .2, .5, 1.)\n",
|
|
"\n",
|
|
"fig, axs = plt.subplots(len(tresholds), figsize=figsize * len(tresholds))\n",
|
|
"\n",
|
|
"for ax, t in zip(axs, tresholds):\n",
|
|
" f = extremum_filter(data, t)\n",
|
|
" raster = rasterizer.vhr(data.spatial[np.logical_not(f)], data.z[np.logical_not(f)], dtype=np.float32)\n",
|
|
" ax.imshow(raster)\n",
|
|
" ax.set_title(str(t))\n",
|
|
" plt.imsave('../Res/extremum_clip_z_{}.png'.format(t), raster)\n",
|
|
"plt.show()\n",
|
|
" "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"0.1 % seems optimal."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Reproductivity"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"data2 = ra.bulk_load('../Data/lidar/C1/', filter_treshold=0)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"tresholds = (0., .1, .2, .5, 1.)\n",
|
|
"\n",
|
|
"fig, axs = plt.subplots(len(tresholds), figsize=figsize * len(tresholds))\n",
|
|
"\n",
|
|
"for ax, t in zip(axs, tresholds):\n",
|
|
" f = extremum_filter(data2, t)\n",
|
|
" raster = rasterizer.vhr(data2.spatial[np.logical_not(f)], data2.z[np.logical_not(f)], dtype=np.float32)\n",
|
|
" ax.imshow(raster)\n",
|
|
" ax.set_title(str(t))\n",
|
|
" plt.imsave('../Res/extremum_clip_z_{}.png'.format(t), raster)\n",
|
|
"plt.show()\n",
|
|
" "
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|