{
 "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
}