{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "import numpy as np\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 = (16, 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load LAS files data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_las(path):\n",
    "    print(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": [
    "files = ['../Data/lidar/C1/272056_3289689.las',\n",
    "         '../Data/lidar/C1/272652_3289689.las',\n",
    "         '../Data/lidar/C1/273248_3289689.las',\n",
    "         '../Data/lidar/C1/273844_3289689.las']\n",
    "\n",
    "folder = Path('../Data/lidar/')\n",
    "files  = find_las(folder)\n",
    "\n",
    "xdata = list()\n",
    "ydata = list()\n",
    "idata = list()\n",
    "\n",
    "for file in files:\n",
    "    infile = laspy.file.File(file)\n",
    "    print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n",
    "    xdata.append(infile.x)\n",
    "    ydata.append(infile.y)\n",
    "    idata.append(infile.intensity)\n",
    "    #idata.append(infile.num_returns)\n",
    "    \n",
    "    \n",
    "data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(data[:,2], 1000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extremum = 0.5\n",
    "tresholds = np.percentile(data[:,2], [extremum, 100 - extremum])\n",
    "display(tresholds)\n",
    "\n",
    "filtered_data = data[:,2].copy()\n",
    "filtered_data[data[:,2] < tresholds[0]] = tresholds[0]\n",
    "filtered_data[data[:,2] > tresholds[1]] = tresholds[1]\n",
    "\n",
    "plt.hist(filtered_data, 1000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(data[:,2] < tresholds[0]).sum(), (data[:,2] > tresholds[1]).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tresholds[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display a 2D hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 0.5\n",
    "ground_size = np.array((data[:,0].max() - data[:,0].min(), data[:,1].max() - data[:,1].min()))\n",
    "display(ground_size)\n",
    "\n",
    "bins = np.rint(ground_size / resolution).astype(int)\n",
    "bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[:,0].max(), data[:,0].min(), data[:,1].max(), data[:,1].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "hist, xedges, yedges = np.histogram2d(data[:,0], data[:,1], bins=bins)\n",
    "plt.imshow(hist.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/hist.png', hist.T, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display intensity map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=figsize)\n",
    "hist_i, xedges, yedges = np.histogram2d(data[:,0], data[:,1], bins=bins, weights=data[:,2])\n",
    "plt.imshow(hist_i.T / hist.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that we have no data at this resolution. Should interpolate ?\n",
    "\n",
    "### Better intensity display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = hist_i.T /  hist.T\n",
    "f = np.logical_not(np.isnan(i))\n",
    "plt.hist(i[f], bins=1000)\n",
    "plt.show()\n",
    "\n",
    "i[f].min(), i[f].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imax_disp = 10000\n",
    "\n",
    "idisp = i\n",
    "idisp[i > imax_disp] = imax_disp\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(idisp, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/hist_int.png', idisp, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create rasters DFC comptabible\n",
    "\n",
    "[Official description](http://www.grss-ieee.org/community/technical-committees/data-fusion/data-fusion-contest/) announce .5 m resolution on the DSMs.\n",
    "\n",
    "## Compare with existing DSMs\n",
    "\n",
    "### Shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfc_dsm = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n",
    "dfc_dsm.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idisp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Offset\n",
    "\n",
    "#### Filter DFC raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodata_filter = dfc_dsm > 1e4\n",
    "dfc_dsm[nodata_filter] = dfc_dsm[nodata_filter == False].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(dfc_dsm.reshape(-1), 1000, label='DFC')\n",
    "plt.hist(idisp[f], 1000, label='our', alpha=.8)\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visual test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_disp = np.moveaxis(np.array((np.flip(idisp, 0), dfc_dsm, np.zeros_like(idisp))), 0, 2)\n",
    "display(offset_disp.shape)\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(offset_disp / np.nanmax(offset_disp), origin='upper', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seems good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/offset.png',offset_disp / np.nanmax(offset_disp), origin='upper')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Better raster with interpolation of missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import griddata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test Griddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridx, gridy = np.mgrid[0:10, 0:10] + 10\n",
    "gridx, gridy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampled_data = np.array([[0, 0, 5], [9, 9, 5], [0, 9, 10], [9, 0, 10]])\n",
    "sampled_data[:, :2] += 10\n",
    "display(sampled_data)\n",
    "\n",
    "inter_data = griddata(sampled_data[:,:2], sampled_data[:,2], (gridx, gridy), method='cubic') # linear, nearest, cubic\n",
    "\n",
    "plt.imshow(inter_data, origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic grid computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_resolution = .01\n",
    "\n",
    "sampled_data = np.array([[0, 0, 5], [9, 9, 5], [0, 9, 10], [20, 0, 10]])\n",
    "\n",
    "def rasterize(coords, values, resolution, method='cubic'):\n",
    "    if isinstance(coords, np.ndarray):\n",
    "        xmin, xmax = coords[:,0].min(), coords[:,0].max()\n",
    "        ymin, ymax = coords[:,1].min(), coords[:,1].max()\n",
    "    else:\n",
    "        xmin, xmax = coords[0].min(), coords[0].max()\n",
    "        ymin, ymax = coords[1].min(), coords[1].max()\n",
    "\n",
    "    gridx, gridy = np.mgrid[xmin:xmax:resolution, ymin:ymax:resolution]\n",
    "    \n",
    "    return griddata(coords, values, (gridx, gridy), method=method) # linear, nearest, cubic\n",
    "\n",
    "plt.imshow(rasterize(sampled_data[:,:2], sampled_data[:,2], test_resolution, 'nearest').T, origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(data[:,2], 1000)\n",
    "plt.show()\n",
    "\n",
    "should_have_keep_idata = data[:,2].copy()\n",
    "should_have_keep_idata[data[:,2] > 1e4] = 1e4\n",
    "\n",
    "plt.hist(should_have_keep_idata, 700)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape, filtered_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treshold = -1\n",
    "#raster = rasterize(data[:treshold,:2], should_have_keep_idata[:treshold], 0.5, 'nearest')\n",
    "raster = rasterize(data[:,:2], filtered_data, 0.5, 'cubic')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster.T, origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#raster = rasterize(data[:treshold,:2], should_have_keep_idata[:treshold], 0.5, 'nearest')\n",
    "raster = rasterize((data[:,0], data[:,1]), data[:,2], 0.5, 'nearest')\n",
    "\n",
    "plt.figure(figsize=figsize)\n",
    "plt.imshow(raster.T, origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imsave('../Res/raster.png', raster.T, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save TIFF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triskele.write('../Res/my_intensity.tiff', np.flip(idisp, 0))"
   ]
  }
 ],
 "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
}