WIP on SelfDualAttributeProfiles

Notebooks/Node Streaming Benchmark.ipynb

@@ -18,6 +18,42 @@
     "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": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "l = ld2dap.LoadTIFF(layers_files[5])\n",
+    "t = ld2dap.Treshold(1e4)\n",
+    "a = ld2dap.AttributeProfiles(area = [1e3, 1e6])\n",
+    "f = ld2dap.Differential()\n",
+    "d = ld2dap.ShowFig('all', fname='../Res/aps.png', pad_inches=0, symb=True, vmax=255)\n",
+    "\n",
+    "t.input = l\n",
+    "a.input = t\n",
+    "f.input = a\n",
+    "d.input = f\n",
+    "d.run()"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -28,16 +64,15 @@
     "trsh = ld2dap.Treshold(1e4)\n",
     "aps  = ld2dap.AttributeProfiles(area=[1e4,1e5], moi=[.1,.5,.9])\n",
     "diff = ld2dap.Differential()\n",
-    "\n",
     "disp = ld2dap.ShowFig('all', True)\n",
     "out  = ld2dap.RawOutput()\n",
-    "out.input  = aps\n",
     "out2 = ld2dap.RawOutput()\n",
     "\n",
-    "out2.input = diff\n",
+    "out.input  = aps\n",
-    "disp.input = diff\n",
+    "#diff.input = aps\n",
     "\n",
-    "diff.input = aps\n",
+    "#out2.input = diff\n",
+    "disp.input = aps\n",
     "\n",
     "aps.input  = trsh\n",
     "\n",
@@ -52,7 +87,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(out.metadata[1])"
+    "'\\rho'.format()"
    ]
   },
   {

Notebooks/SDAPs.ipynb

@@ -0,0 +1,125 @@
+{
+ "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",
+    "\n",
+    "triskele_path = Path('../triskele/python/')\n",
+    "sys.path.append(str(triskele_path.resolve()))\n",
+    "import triskele\n",
+    "\n",
+    "def DFC_filter(raster):\n",
+    "    raster[raster > 1e4] = raster[raster < 1e4].max()\n",
+    "\n",
+    "    \n",
+    "def show(im, im_size=1, save=None):\n",
+    "    plt.figure(figsize=(16*im_size,3*im_size))\n",
+    "    plt.imshow(im)\n",
+    "    plt.colorbar()\n",
+    "    \n",
+    "    if save is not None:\n",
+    "        plt.savefig(save, bbox_inches='tight', pad_inches=1)\n",
+    "    \n",
+    "    plt.show()\n",
+    "\n",
+    "def mshow(Xs, titles=None, im_size=1, save=None):\n",
+    "    s = len(Xs)\n",
+    "\n",
+    "    plt.figure(figsize=(16*im_size,3*im_size*s))\n",
+    "\n",
+    "    for i in range(s):\n",
+    "        plt.subplot(s,1,i+1)\n",
+    "        plt.imshow(Xs[i])\n",
+    "        \n",
+    "        if titles is not None:\n",
+    "            plt.title(titles[i])\n",
+    "            \n",
+    "        plt.colorbar()\n",
+    "    \n",
+    "    if save is not None:\n",
+    "        plt.savefig(save, bbox_inches='tight', pad_inches=1)\n",
+    "        \n",
+    "    plt.show()"
+   ]
+  },
+  {
+   "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']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layers = list()\n",
+    "\n",
+    "for file in layers_files:\n",
+    "    print('Loading {}'.format(file))\n",
+    "    layer = triskele.read(file)\n",
+    "    DFC_filter(layer)\n",
+    "    layers.append(layer)\n",
+    "\n",
+    "layers_stack = np.stack(layers, axis=2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t = triskele.Triskele(layers_stack, verbose=False)\n",
+    "\n",
+    "attributes = t.filter(tree='min-tree',\n",
+    "                      area=[1e3,1e5], \n",
+    "                      standard_deviation=[.3,.9],\n",
+    "                      moment_of_inertia=[.1,.4]\n",
+    "                     )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "attributes.shape"
+   ]
+  }
+ ],
+ "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/SelfDualAttributeProfiles.py

@@ -0,0 +1,116 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+# \file SelfDualAttributeProfiles.py
+# \brief TODO
+# \author Florent Guiotte <florent.guiotte@gmail.com>
+# \version 0.1
+# \date 17 avril 2018
+#
+# TODO details
+
+from .core import Filter, Stack
+
+## TODO: dep
+import sys
+import numpy as np
+sys.path.append('../triskele/python')
+import triskele
+
+class AttributeProfiles(Filter):
+    def __init__(self, area=None, sd=None, moi=None):
+        super().__init__(self.__class__.__name__)
+        self.area = np.sort(area) if area is not None else None
+        self.sd   = np.sort(sd) if sd is not None else None
+        self.moi  = np.sort(moi) if moi is not None else None
+
+    def _process_desc(self):
+        att_desc = dict()
+        for att in ['area', 'sd', 'moi']:
+            att_desc[att] = list()
+            if self.__getattribute__(att) is not None:
+                att_desc[att].append(None)
+                att_desc[att].extend(
+                  ['Self-dual {} {}'.format(att, x) for x in self.__getattribute__(att)])
+
+        return att_desc
+
+    def _process_symb(self):
+        att_symb = dict()
+        for att in ['area', 'sd', 'moi']:
+            att_symb[att] = list()
+            if self.__getattribute__(att) is not None:
+                att_symb[att].append(None)
+                att_symb[att].extend(
+                  ['\rho^{{{}}}_{{{}}}'.format(att, x) for x in self.__getattribute__(att)])
+
+        return att_symb
+
+    def _process_len(self):
+        att_len = dict()
+        for att in ['area', 'sd', 'moi']:
+            values = self.__getattribute__(att)
+            att_len[att] = len(values) if values is not None else 0
+        return att_len
+
+    def _process(self, data, metadata):
+        t       = triskele.Triskele(data, verbose=False)
+        att_min = t.filter(tree='tos-tree', area=self.area,
+                              standard_deviation=self.sd,
+                              moment_of_inertia=self.moi)
+
+        ## Create new metadata
+        ### Pre-process descriptions
+        att_desc = self._process_desc()
+        att_symb = self._process_symb()
+
+        ### Compute stack offsets and att length
+        att_len       = self._process_len()
+        raster_offset = sum(att_len.values())
+
+        ### Merge old and new descriptions
+        metadata_new = list()
+        ### Re-order to create original APs
+        raster_list  = list()
+        for stack in metadata:
+            if stack.end - stack.begin > 1:
+                raise NotImplementedError('Nested filtering not implemented yet')
+
+            do = 0
+            sb = stack.begin * (raster_offset + 1)
+            for att in ['area', 'sd', 'moi']:
+                if att_offset[att] == 0:
+                    continue
+
+                al = att_len[att]
+                raster_list.append(att_min[:,:,sb+do+al:sb+do:-1])
+                raster_list.append(att_min[:,:,sb])
+                print('DEBUG: copying layer {}'.format(sb))
+                raster_list.append(att_max[:,:,sb+do+1:sb+do+al+1])
+                do += al
+
+                stack_new = Stack(dso + stack_offset * stack.begin, att_offset[att],
+                        stack.desc[0], stack.symb[0])
+
+                for old_desc, new_desc in zip(stack_new.desc, att_desc[att]):
+                    print('DESCRIPTION: {} > {}'.format(old_desc, new_desc))
+                    old_desc.append(new_desc)
+
+                for old_symb, new_symb in zip(stack_new.symb, att_symb[att]):
+                    print('symbRIPTION: {} > {}'.format(old_symb, new_symb))
+                    old_symb.append(new_symb)
+
+                metadata_new.append(stack_new)
+
+
+        data_new = np.dstack(raster_list)
+
+        return data_new, metadata_new
+
+if __name__ == '__main__':
+    area = [10, 100, 1000]
+    sd   = [.1, .9]
+
+    ap = AttributeProfiles(area, sd)
+
+    print(ap._process_desc())
+    print(ap._process_offset())

ld2dap/ShowFig.py

@@ -13,12 +13,14 @@ from ld2dap.core import Output
 import matplotlib.pyplot as plt
 
 class ShowFig(Output):
-    def __init__(self, stack_id=0, symb=False, bbox_inches='tight', pad_inches=1):
+    def __init__(self, stack_id=0, symb=False, fname=None, bbox_inches='tight', pad_inches=1, vmax=None):
         super().__init__(self.__class__.__name__)
         self.bbox_inches = bbox_inches
         self.pad_inches  = pad_inches
         self.stack_id    = stack_id
         self.symb        = symb
+        self.fname       = fname
+        self.vmax        = vmax
 
     def _process(self, data, metadata):
         if self.stack_id == 'all':
@@ -36,7 +38,7 @@ class ShowFig(Output):
 
         for i, di in enumerate(range(stack.begin, stack.end)):
             f1 = fig.add_subplot(fig_count, 1, i + 1)
-            img = f1.imshow(data[:,:,di])
+            img = f1.imshow(data[:,:,di], vmax=self.vmax)
             plt.colorbar(img)
             if self.symb:
                 plt.rc('text', usetex=True)
@@ -44,5 +46,7 @@ class ShowFig(Output):
             else:
                 f1.set_title(' > '.join(filter(None, stack.desc[i])))
 
+        if self.fname is not None:
+            fig.savefig(self.fname, bbox_inches=self.bbox_inches, pad_inches=self.pad_inches)
         plt.show()
         plt.close(fig)

ld2dap/Treshold.py

@@ -17,6 +17,7 @@ class Treshold(Filter):
         self.max_value = max_value #if max_value is not None else treshold
 
     def _process(self, data, metadata):
+        # TODO: UPGRADE RASTER DEPENDANCE
         if self.max_value is None:
             self.max_value = data[data < self.treshold].max()
 
@@ -25,4 +26,7 @@ class Treshold(Filter):
                 d.append('treshold {}'.format(self.treshold))
                 #s.append('T_{{{}}}'.format(self.treshold))
 
+            ## TODO: TMP FIX
+            #data[:,:,stack.begin][data[:,:,stack.begin] > self.treshold] = data[:,:,stack.begin][data[:,:,stack.begin]].max()
+        
         return data * (data < self.treshold) + self.max_value * (data >= self.treshold), metadata

test.py

@@ -17,7 +17,7 @@ import numpy as np
 def main():
     i = LoadTIFF(['Data/test.tiff', 'Data/test2.tiff'])
     t = Treshold(1e4)
-    ap = APs(np.array([100,1e3,1e4]))
+    ap = APs([100,1e3,1e4])
     o = SaveFig('Res/test.png')
     s = ShowFig(0, True)
     d = Differential()