Add LocalFeatures filter

Notebooks/HAPs.ipynb

@@ -77,7 +77,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "raster_p = Path('../Data/phase1_rasters/Intensity_C3/UH17_GI3F051_TR.tif')\n",
+    "raster_p = Path('../Data/phase1_rasters/DEM+B_C123/UH17_GEM051_TR.tif')\n",
     "raster   = triskele.read(raster_p)\n",
     "DFC_filter(raster)"
    ]
@@ -105,9 +105,11 @@
    "outputs": [],
    "source": [
     "#area = np.array([25, 100, 500, 1e3, 5e3, 10e3, 50e3, 100e3, 150e3])\n",
-    "area = [100, 1e3, 1e4, 1e5]\n",
+    "area = [1e3, 1e6]\n",
     "moi  = [.5, .7, .9]\n",
     "sd   = [.1, .5, .9]\n",
+    "moi = None\n",
+    "sd = None\n",
     "\n",
     "t = triskele.Triskele(raster, verbose=False)\n",
     "attributes = t.filter(tree='tos-tree',\n",
@@ -147,7 +149,7 @@
    "outputs": [],
    "source": [
     "a = attributes\n",
-    "i = 3\n",
+    "i = 1\n",
     "show(a[:,:,i].astype(np.float) - a[:,:,i+1])"
    ]
   },
@@ -158,6 +160,15 @@
     "## Compute patches"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "patch_size = 5"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -189,7 +200,7 @@
     "            stack.append(offset)\n",
     "    return np.stack(stack, axis=-1)\n",
     "\n",
-    "offset_stack = create_patches(attributes, 5)\n",
+    "offset_stack = create_patches(attributes, patch_size)\n",
     "offset_stack.shape"
    ]
   },
@@ -263,6 +274,36 @@
     "stack_max  = np.max(offset_stack, axis=-1)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.std.__name__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "create_patches.__name__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for func in (np.mean, np.var):\n",
+    "    tmp = func(offset_stack, axis=-1)\n",
+    "    show(tmp[:,:,0])\n",
+    "    print(tmp.shape)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -314,6 +355,159 @@
     "       ]\n",
     "mshow(figs, ttls, 2, save='mstack2.png')"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### TMP TRISLEK DBG"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "normal = triskele.read('../triskele/build/out/test-default.tif')\n",
+    "sd     = triskele.read('../triskele/build/out/test-sd.tif')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "normal.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "areas = np.loadtxt('../triskele/data/areaThresholds.txt')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figs = [normal[:,:,0]]\n",
+    "ttls = ['Origin']\n",
+    "\n",
+    "for area, i in zip(areas, range(1,normal.shape[2])):\n",
+    "    figs.append(normal[:,:,i])\n",
+    "    ttls.append('SDAP area {}'.format(area))\n",
+    "\n",
+    "mshow(figs, ttls, 2, '../Res/trkldbg_triskele_ap.png')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figs = [sd[:,:,0]]\n",
+    "ttls = ['Origin']\n",
+    "\n",
+    "for area, i in zip(areas, range(1,sd.shape[2])):\n",
+    "    figs.append(sd[:,:,i])\n",
+    "    ttls.append('Triskele LFSDAP (SD) area {}'.format(area))\n",
+    "\n",
+    "mshow(figs, ttls, 2, '../Res/trkldbg_triskele_sd.png')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(normal == sd).all()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "offset_stack = create_patches(normal, 5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stack_std = np.std(offset_stack, axis=-1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figs = [stack_std[:,:,0]]\n",
+    "ttls = ['Origin std 5x5']\n",
+    "\n",
+    "for area, i in zip(areas, range(1,sd.shape[2])):\n",
+    "    figs.append(stack_std[:,:,i])\n",
+    "    ttls.append('LD2DAPs LFSDAP std 5x5 area {}'.format(area))\n",
+    "\n",
+    "mshow(figs, ttls, 2, '../Res/trkldbg_ld2dap_sd.png')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "figs = list()\n",
+    "ttls = list()\n",
+    "for d in range(stack_std.shape[-1]):\n",
+    "    if d == 0:\n",
+    "        ttls.append('Origin')\n",
+    "    else:\n",
+    "        ttls.append('STD 5x5 from area {}'.format(area[d-1]))\n",
+    "    figs.append(stack_std[:,:,d])\n",
+    "    \n",
+    "mshow(figs, ttls, 2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert True"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert False"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

Notebooks/Node LFAPs.ipynb

@@ -0,0 +1,164 @@
+{
+ "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 collections\n",
+    "\n",
+    "ld2dap_path = Path('../')\n",
+    "sys.path.append(str(ld2dap_path.resolve()))\n",
+    "import ld2dap\n",
+    "from ld2dap.core import Filter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layers_files = [\n",
+    "    '../Data/phase1_rasters/DEM+B_C123/UH17_GEM051_TR.tif',\n",
+    "    '../Data/phase1_rasters/DEM_C123_3msr/UH17_GEG051_TR.tif',\n",
+    "    '../Data/phase1_rasters/DEM_C123_TLI/UH17_GEG05_TR.tif',\n",
+    "    '../Data/phase1_rasters/DSM_C12/UH17c_GEF051_TR.tif',\n",
+    "    '../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif',\n",
+    "    '../Data/phase1_rasters/Intensity_C2/UH17_GI2F051_TR.tif',\n",
+    "    '../Data/phase1_rasters/Intensity_C3/UH17_GI3F051_TR.tif'\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Node Streaming Local FeaturesAttribute Profiles\n",
+    "\n",
+    "We made some custom LFAPs in Python, but now Triskele comes with its owns in C++! Let's try it out then create a streaming node.\n",
+    "\n",
+    "## Triskele LFAPs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "sys.path.append('../triskele/python/')\n",
+    "import triskele"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "raster = triskele.read(layers_files[0])\n",
+    "raster.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lfaps = triskele.Triskele(raster)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results = lfaps.filter(area=[100, 1e4], feature='SD')\n",
+    "results.shape, lfaps.process"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Streamin Node"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import logging\n",
+    "logger = logging.getLogger()\n",
+    "logger.setLevel(logging.DEBUG)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "l = ld2dap.LoadTIFF(layers_files)\n",
+    "t = ld2dap.Treshold(1e4)\n",
+    "a = ld2dap.LocalFeatures([np.var, np.std, ])#, sd=[.4,.6,.8], moi=[.5,.9])\n",
+    "d = ld2dap.ShowFig(stack_id='all', symb=False)\n",
+    "o = ld2dap.RawOutput()\n",
+    "\n",
+    "o.input = a\n",
+    "d.input = a\n",
+    "a.input = t\n",
+    "t.input = l\n",
+    "\n",
+    "d.run()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(o.metadata[0])"
+   ]
+  },
+  {
+   "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
+}

ld2dap/LocalFeatures.py

@@ -0,0 +1,94 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+# \file LocalFeatures.py
+# \brief TODO
+# \author Florent Guiotte <florent.guiotte@gmail.com>
+# \version 0.1
+# \date 23 août 2018
+#
+# TODO details
+
+from ld2dap.core import Filter, Stack
+import numpy as np
+
+class LocalFeatures(Filter):
+    def __init__(self, features=None, patch_size=3):
+        super().__init__(self.__class__.__name__)
+        self.features = features
+        self.patch_size = patch_size
+
+    def _process(self, data, metadata):
+        self.logger.info('Local Features Filtering')
+        stacker = Stacker()
+
+        for stack in metadata:
+            self.logger.debug('{}'.format(stack))
+
+            patches = create_patches(data[:,:,stack.begin:stack.end], self.patch_size)
+
+            for function in self.features:
+                self.logger.debug('Features function {}'.format(function.__name__))
+                features = function(patches, axis=-1)
+
+                stacker.auto_stack(features, stack, 
+                        'Local Features {}'.format(function.__name__), 
+                        'f^{{{}}}'.format(function.__name__))
+
+        return stacker.pack()
+        
+def create_patches(array, patch_size=3):
+        amp = int((patch_size - 1 ) / 2)
+
+        stack = list()
+        for i in range(-amp, amp+1):
+            ai = i if i > 0 else None
+            bi = i if i < 0 else None
+            ci = -bi if bi is not None else None
+            di = -ai if ai is not None else None
+
+            for j in range(-amp, amp+1):
+                offset = np.zeros_like(array)
+                #offset = np.empty(array.shape)
+                aj = j if j > 0 else None
+                bj = j if j < 0 else None
+                cj = -bj if bj is not None else None
+                dj = -aj if aj is not None else None
+                offset[ai:bi, aj:bj] = array[ci:di, cj:dj]
+                stack.append(offset)
+        return np.stack(stack, axis=-1)
+
+class Stacker:
+    def __init__(self):
+        self.offset = 0
+        self.data_list = list()
+        self.metadata_list = list()
+
+    def auto_stack(self, data, old_stack, desc_str, symb_str):
+        if data.shape[-1] != (old_stack.end - old_stack.begin):
+            raise ValueError("Mismatch between data shape and stack length")
+
+        new_desc = list()
+        for d in old_stack.desc:
+            new_desc.append(d.copy())
+            new_desc[-1].append(desc_str)
+
+        new_symb = list()
+        for s in old_stack.symb:
+            new_symb.append([symb_str + '('])
+            new_symb[-1].extend(s.copy())
+            new_symb[-1].append(')')
+
+        self.stack(data, new_desc, new_symb)
+
+    def stack(self, data, desc, symb):
+        size_new  = data.shape[-1]
+        stack_new = Stack(self.offset, size_new)
+        self.offset += size_new
+        stack_new.desc = desc
+        stack_new.symb = symb 
+        self.data_list.append(data)
+        self.metadata_list.append(stack_new)
+
+    def pack(self):
+        return np.dstack(self.data_list), self.metadata_list
+

ld2dap/Treshold.py

@@ -20,6 +20,7 @@ class Treshold(Filter):
 
     def _process(self, data, metadata):
         # TODO: UPGRADE STACK DEPENDANCE
+        # TODO: Verify if the previous TODO is up to date
         # This filter each raster independently
         self.logger.info('Filtering')
 

ld2dap/__init__.py

@@ -10,6 +10,7 @@
 
 from .AttributeProfiles import AttributeProfiles
 from .SelfDualAttributeProfiles import SelfDualAttributeProfiles
+from .LocalFeatures import LocalFeatures
 from .Treshold import Treshold
 from .LoadTIFF import LoadTIFF
 from .SaveFig import SaveFig