{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TRYSKELE Wrapper\n",
    "\n",
    "Dependence \n",
    "--\n",
    "\n",
    "`tryskele` project folder (compiled) in jupyter's root directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import libtiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "triskele_bin = Path('../triskele/build/out/apGenerator')\n",
    "triskele_bin.exists()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process = [triskele_bin.absolute(),\n",
    "          '']\n",
    "\n",
    "tryskele = subprocess.Popen(process, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "print('STDOUT:\\n' + tryskele.stdout.read().decode() + '\\nSTDERR:\\n' + tryskele.stderr.read().decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Files\n",
    "#infile  = Path('../Data/phase1_rasters/DEM_C123_TLI/UH17_GEG05_TR.tif')\n",
    "#infile  = Path('../Res/dem_u16.tif')\n",
    "infile  = Path('../Res/dem_u8.tif')\n",
    "#infile  = Path('/home/florent/Public/WIPDir/triskele/data/10m.tif')\n",
    "outfile = Path('../Res/out.tiff')\n",
    "\n",
    "area_conf      = Path('../Res/area.txt')\n",
    "deviation_conf = Path('../Res/deviation.txt')\n",
    "inertia_conf   = Path('../Res/inertia.txt')\n",
    "\n",
    "# Attributes definition\n",
    "area      = np.array([10,100,1e3,1e4,1e5,5e5])#np.arange(20) * 100000 #* 10000#np.array([10,100,1000])\n",
    "deviation = np.array([0.9,0.99,0.999,0.9999])#,1e4,1e5,5e5])\n",
    "inertia   = 6.02911 / (np.arange(3)[::-1] + 5)#np.array([10,100,1000])\n",
    "np.savetxt(area_conf, area, fmt='%d')\n",
    "np.savetxt(deviation_conf, deviation, fmt='%f')\n",
    "np.savetxt(inertia_conf, inertia, fmt='%f')\n",
    "\n",
    "process = [triskele_bin.absolute(),\n",
    "           '-i', infile.absolute(),\n",
    "           '--tos-tree',\n",
    "           '--area', area_conf.absolute(),\n",
    "           '--standard-deviation', deviation_conf.absolute(),\n",
    "           #'--moment-of-inertia', inertia_conf.absolute(),\n",
    "           '-o', outfile.absolute()]\n",
    "\n",
    "display(' '.join(str(v) for v in process))\n",
    "\n",
    "tryskele = subprocess.Popen(process, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "print('Return code: ' + str(tryskele.returncode) + '\\nSTDOUT:\\n' + tryskele.stdout.read().decode() + '\\nSTDERR:\\n' + tryskele.stderr.read().decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process = libtiff.TIFF.open(outfile)\n",
    "raster  = process.read_image()\n",
    "\n",
    "for chan in range(raster.shape[2]):\n",
    "    plt.figure(figsize=(16,3))\n",
    "    plt.imshow(raster[:,:,chan])\n",
    "    plt.colorbar()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process = libtiff.TIFF.open(outfile)\n",
    "raster  = process.read_image()\n",
    "\n",
    "rasterf = raster.astype(np.float)\n",
    "\n",
    "for chan in range(raster.shape[2] - 1):\n",
    "    plt.figure(figsize=(16*3,3*3))\n",
    "    plt.imshow(rasterf[:,:,chan+1] - rasterf[:,:,chan])\n",
    "    plt.colorbar()\n",
    "    plt.show()\n",
    "    plt.imsave('../Res/img/area_' + str(chan) + '.png', rasterf[:,:,chan+1] - rasterf[:,:,chan])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_in = libtiff.TIFF.open(infile).read_image()\n",
    "raster_in.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(raster[:,:,0] == raster_in).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "500k pixels of area for the biggest structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasterf[:,:,0] - rasterf[:,:,9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_df = pd.DataFrame(out.reshape(-1, out.shape[2]))\n",
    "out_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(out_df, return_counts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(out_df[0] * out_df[0] , 256, log=True)\n",
    "plt.show()"
   ]
  }
 ],
 "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
}