diff --git a/Notebooks/Node SDAPs.ipynb b/Notebooks/Node SDAPs.ipynb new file mode 100644 index 0000000..a64b619 --- /dev/null +++ b/Notebooks/Node SDAPs.ipynb @@ -0,0 +1,285 @@ +{ + "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 Self-Dual Attribute Profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "l = ld2dap.LoadTIFF(layers_files[5:7])\n", + "t = ld2dap.Treshold(1e4)\n", + "a = ld2dap.SelfDualAttributeProfiles(area = [1e3, 1e6], sd = [.5])\n", + "d = ld2dap.ShowFig(stack_id='all', symb=True)\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[1]), o.data.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Algorithm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Metadata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = {'area': 2, 'sd': 3, 'moi': 2}\n", + "raster_count = 2\n", + "\n", + "att_len = list(A.values())\n", + "att_len = np.tile(att_len, raster_count)\n", + "display(att_len)\n", + "\n", + "start = np.cumsum(att_len)[:-1]\n", + "start = np.hstack(([0], start))\n", + "start\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Data\n", + "\n", + "Duplicate origin in attributes to respect Stack construction\n", + "\n", + "#### Insert a raster in a stack" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stack = o.data.copy()\n", + "raster = stack[:,:,0]\n", + "\n", + "stack.shape, raster.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "where = 3\n", + "nstack = np.insert(stack, where, raster, axis=2)\n", + "nstack.shape, (nstack[:,:,0] == nstack[:,:,where]).all()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Insert same raster in multiple places" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "raster_broadcast = np.tile(raster, (2, 1, 1))\n", + "raster_broadcast = np.rollaxis(raster_broadcast, 0, 3)\n", + "raster_broadcast.shape, (raster_broadcast[:,:,1] == raster).all()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nstack = np.insert(stack, (3), raster_broadcast, axis=2)\n", + "nstack.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check if offset is ok:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(stack[:,:,3] == nstack[:,:,4]).all()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Repeat multiple origins" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "raster_count = 2\n", + "att_len = {'area': 3, 'sd': 0, 'moi': 0}\n", + "# Ignore\n", + "where = np.array(list(att_len.values()))\n", + "where = where[where != 0] - 1\n", + "where[0] += 1\n", + "count = sum(where)\n", + "# /Ignore\n", + "origins_index = np.arange(raster_count) * count\n", + "display(origins_index)\n", + "origins = stack[:,:,origins_index]\n", + "origins.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.arange(16*9*4).reshape((16,9,4))\n", + "A[:,:,[0,1]].shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Index where origin should be placed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "raster_count = 4\n", + "\n", + "att_len = {'area': 3, 'sd': 2, 'moi': 0}\n", + "where = np.array(list(att_len.values()))\n", + "where = where[where != 0] - 1\n", + "\n", + "where[0] += 1\n", + "display(where)\n", + "count = sum(where)\n", + "display(count)\n", + "\n", + "where = np.cumsum(where[:-1])\n", + "display(where)\n", + "\n", + "offset = np.repeat(np.arange(raster_count) * count, where.size)\n", + "display(offset)\n", + "where = np.tile(where, (raster_count)) + offset\n", + "where" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = list()\n", + "A.size" + ] + } + ], + "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 +} diff --git a/Notebooks/Untitled1.ipynb b/Notebooks/Untitled1.ipynb new file mode 100644 index 0000000..2fd6442 --- /dev/null +++ b/Notebooks/Untitled1.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/ld2dap/SelfDualAttributeProfiles.py b/ld2dap/SelfDualAttributeProfiles.py index 31c5650..48309c3 100644 --- a/ld2dap/SelfDualAttributeProfiles.py +++ b/ld2dap/SelfDualAttributeProfiles.py @@ -10,13 +10,14 @@ from .core import Filter, Stack -## TODO: dep +# TODO: dep import sys import numpy as np sys.path.append('../triskele/python') import triskele -class AttributeProfiles(Filter): + +class SelfDualAttributeProfiles(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 @@ -47,70 +48,83 @@ class AttributeProfiles(Filter): def _process_len(self): att_len = dict() + att_len_cs = dict() + cs = 0 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 + al = len(values) + 1 if values is not None else 0 + att_len[att] = al + att_len_cs[att] = cs + cs += al + self.logger.debug('Attribute length: {}'.format(att_len)) + self.logger.debug('Attribute length CS: {}'.format(att_len_cs)) + return att_len, att_len_cs def _process(self, data, metadata): - t = triskele.Triskele(data, verbose=False) - att_min = t.filter(tree='tos-tree', area=self.area, + t = triskele.Triskele(data, verbose=False) + attributes = t.filter(tree='tos-tree', area=self.area, standard_deviation=self.sd, moment_of_inertia=self.moi) - ## Create new metadata - ### Pre-process descriptions + # Create new data and metadata + metadata_new = list() + + # Pre-process descriptions att_desc = self._process_desc() att_symb = self._process_symb() - ### Compute stack offsets and att length - att_len = self._process_len() + # Compute stack offsets and att length + att_len, att_len_cs = 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() + # Duplicate origin in data to respect Stack structure + # Compute insert index in where + where = np.array(list(att_len.values())) + where = where[where != 0] - 1 + where[0] += 1 + count = sum(where) + where = np.cumsum(where[:-1]) + offset = np.repeat(np.arange(len(metadata)) * count, + where.size) # Can't nest this + where = np.tile(where, (len(metadata)) + offset + + # TODO: + # - Repeat multiple origins + # - Duplicate origins to match where + # - Insert origins where + + data_new = None + for stack in metadata: if stack.end - stack.begin > 1: + self.logger.err('Nested filtering, raising error') 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: + if att_len[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]) + stack_new = Stack(raster_offset * stack.begin + att_len_cs[att], + att_len[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)) + self.logger.debug('Desc: {} + {}'.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)) + self.logger.debug('Symb: {} + {}'.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] + sd = [.1, .9] - ap = AttributeProfiles(area, sd) + ap = SelfDualAttributeProfiles(area, sd) print(ap._process_desc()) - print(ap._process_offset()) diff --git a/ld2dap/__init__.py b/ld2dap/__init__.py index 29b4120..af223f4 100644 --- a/ld2dap/__init__.py +++ b/ld2dap/__init__.py @@ -9,6 +9,7 @@ # TODO details from .AttributeProfiles import AttributeProfiles +from .SelfDualAttributeProfiles import SelfDualAttributeProfiles from .Treshold import Treshold from .LoadTIFF import LoadTIFF from .SaveFig import SaveFig diff --git a/ld2dap/core/Input.py b/ld2dap/core/Input.py index 77dec73..4499e17 100644 --- a/ld2dap/core/Input.py +++ b/ld2dap/core/Input.py @@ -10,6 +10,7 @@ from .Node import Node + class Input(Node): def __init__(self, name='__CHILD__'): super().__init__('I:{}'.format(name)) diff --git a/test.py b/test.py index 5dda6e4..3836fcf 100755 --- a/test.py +++ b/test.py @@ -11,6 +11,7 @@ #from core import Input, Output, Filter from ld2dap import LoadTIFF, SaveFig, Treshold, ShowFig, Differential from ld2dap import AttributeProfiles as APs +from ld2dap import SelfDualAttributeProfiles as SDAPs from ld2dap.core import logger import numpy as np @@ -19,7 +20,7 @@ def main(): logger.setup_logging() i = LoadTIFF(['Data/test.tiff', 'Data/test2.tiff']) t = Treshold(1e4) - ap = APs([100,1e3,1e4]) + ap = SDAPs([100,1e3], [.5,.7,.9], [.1,.3]) o = SaveFig('Res/test.png') s = ShowFig(0, True) d = Differential()