{
 "cells": [
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "# Tresholds for Custom Raster from DFC LiDAR data\n",
    "\n",
    "Compare our results with the DFC rasters and set the tresholds for the raster factory.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load DFC raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfc_raster = triskele.read('../Data/phase1_rasters/Intensity_C3/UH17_GI3F051_TR.tif')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(dfc_raster)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The raster from DFC dataset are noised with high value noise. We need to filter high values. We empirically set the treshold to 1e4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.clip(dfc_raster, dfc_raster.min(), 1e4, out=dfc_raster)\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(dfc_raster)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set filtering and clipping treshold to process rasters from LiDAR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data without filtering or clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=0, clip_treshold=0, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_c0 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_c0)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have high value noise, but far better than the DFC noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data without filtering and minimal clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=0, clip_treshold=0.01, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_c0_01 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_c0_01)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clipping does not remove unwanted high value noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data with minimal filtering and no clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C1', 'C3', filter_treshold=0.01, clip_treshold=0, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filtering remove high value noise, but the tone mapping is bad (too dark)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data with filtering and no clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=0.1, clip_treshold=0, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_1_c0 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_1_c0)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tone mapping is correct, but interpolation artifacts appears where too much points are removed from filtering (e.g. in the stadium)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data without filtering and with clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=0, clip_treshold=0.1, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_c0_1 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_c0_1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tone map is correct, no interpolation artifact but high noise sparkle the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data with minimal filtering and minimal clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C1', 'C3', filter_treshold=0.01, clip_treshold=0.01, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_01 = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_01)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tone map is not correct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data with minimal filtering and normal clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C3 = ra.bulk_load('../Data/lidar/C3', 'C3', filter_treshold=0.2, clip_treshold=0.1, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_1 = ra.rasterize_cache('z', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tone map is correct, no interpolation artifact and low high noise in the result. We will now on choose "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load C123 data with minimal filtering and normal clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = (0.01 + 0.01 + 0.2) / 3\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C123 = ra.bulk_load('../Data/lidar/', 'C123', filter_treshold=0.08, clip_treshold=0.1, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we process the raster with the same resolution and a nearest interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_C123 = ra.rasterize_cache('z', C123, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tone map is correct, no interpolation artifact and low high noise in the result. We will now on choose "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare interpolation method\n",
    "\n",
    "### Nearest neighbour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_1_nearest = ra.rasterize_cache('intensity', C3, .5, 'nearest', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1_nearest)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_1 = ra.rasterize_cache('intensity', C3, .5, 'linear', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cubic interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_1 = ra.rasterize_cache('intensity', C3, .5, 'cubic', False, cache_dir='../Res/enrichment_rasters')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cubic interpolation seems to create negative values, maybe at the same spots of the DFC high noise ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((raster_f0_01_c0_1 < 0) * 1.)\n",
    "plt.colorbar()\n",
    "plt.title('Cubic low noise')\n",
    "plt.show()\n",
    "\n",
    "dfc_raster_raw = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((dfc_raster_raw > 1e4) * 1.)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(np.logical_and((dfc_raster_raw > 1e4), (raster_f0_01_c0_1 < 0)) * 1)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise and Cubic low noise')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((dfc_raster_raw > 1e4) * 1 - (raster_f0_01_c0_1 < 0) * 1)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise minus Cubic low noise')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerous common noise pixel between DFC noise and our cubic interpolation.\n",
    "\n",
    "Let's try with our high noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((raster_f0_01_c0_1 > raster_f0_01_c0_1_nearest.max()) * 1.)\n",
    "plt.colorbar()\n",
    "plt.title('Cubic high noise')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((dfc_raster_raw > 1e4) * 1.)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(np.logical_and((dfc_raster_raw > 1e4), (raster_f0_01_c0_1 > raster_f0_01_c0_1_nearest.max())) * 1)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise and Cubic low noise')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow((dfc_raster_raw > 1e4) * 1 - (raster_f0_01_c0_1 > raster_f0_01_c0_1_nearest.max()) * 1)\n",
    "plt.colorbar()\n",
    "plt.title('DFC high noise minus Cubic low noise')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Very low correlation between our raster and the DFC high noise.\n",
    "\n",
    "### Filter low and high interpolated values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_f0_01_c0_1_postprocess = np.clip(raster_f0_01_c0_1, C3.intensity.min(), C3.intensity.max())\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster_f0_01_c0_1_postprocess)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = ra.rasterize_cache('intensity', C3, .5, 'cubic-clip', False, cache_dir='../Res/enrichment_rasters')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C12 = ra.bulk_load(['../Data/lidar/C1', '../Data/lidar/C2'], 'C12', filter_treshold=1., clip_treshold=0.1, dtype=np.float32)"
   ]
  },
  {
   "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
}