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