From c71f86a33bb354237cab9393ed138a4e3666c70d Mon Sep 17 00:00:00 2001
From: Michael Krayer <michael.krayer@posteo.de>
Date: Fri, 6 Aug 2021 00:34:32 +0200
Subject: [PATCH] implemented isolation/triangulation. needs to be tested

---
 field.py | 72 +++++++++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 64 insertions(+), 8 deletions(-)

diff --git a/field.py b/field.py
index 1ae057b..6755f3d 100644
--- a/field.py
+++ b/field.py
@@ -615,6 +615,7 @@ class BinaryFieldNd:
         self.labels      = None
         self.nlabels     = 0
         self.wrap        = tuple(self._ndim*[None])
+        self._feat_slice = None
         self.set_structure(False)
         self.set_periodicity(self._ndim*[False])
 
@@ -653,6 +654,7 @@ class BinaryFieldNd:
             self.labels,self.nlabels,self.wrap = self._labels_periodic()
         else:
             self.labels,self.nlabels = ndimage.label(self.data,structure=self.structure)
+        self._feat_slice = ndimage.find_objects(self.labels)
 
     def _labels_periodic(self,map_to_zero=False):
         '''Label features in an array while taking into account periodic wrapping.
@@ -735,6 +737,9 @@ class BinaryFieldNd:
             self.data = self._labels_periodic(map_to_zero=True)[0]>0
         # Invert to get the final result
         np.logical_not(self.data,self.data)
+        # If labels have been computed already, recompute them to stay consistent
+        if self.labels is not None:
+            self.label()
 
     def probe(self,idx):
         '''Returns whether or not a point at idx is True or False.'''
@@ -757,6 +762,9 @@ class BinaryFieldNd:
             self.label()
         if label is None:
             return np.bincount(self.labels.ravel())
+            # TBD: compare performance with
+            # sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))
+            # http://scipy-lectures.org/advanced/image_processing/auto_examples/plot_find_object.html
         else:
             return np.sum(self.labels==label)
     
@@ -775,27 +783,75 @@ class BinaryFieldNd:
             self.label()
         selection = self._select_feature(selection)
         # Map tagged regions to zero in order to discard them
-        map_ = np.array(range(0,self.count+1),dtype=self.label.dtype)
+        map_ = np.array(range(self.nlabels+1),dtype=self.labels.dtype)
         map_[selection] = 0
         # Remove gaps from target mapping
-        map_ = np.unique(map_,return_inverse=True)[1]
+        idx_,map_ = np.unique(map_,return_index=True,return_inverse=True)[1:3]
         # Discard regions
         self.labels  = map_[self.labels]
         self.nlabels = np.max(map_)
-        self.data    = self.labels>0
+        for axis in range(self._ndim):
+            if self.wrap[axis] is not None:
+                self.wrap[axis] = self.wrap[axis][idx_]
+        self.data = self.labels>0
+
+    def isolate_feature(self,selection):
+        from scipy import ndimage
+        if self.labels is None:
+            self.label()
+        selection = self._select_feature(selection)
+        output = []
+        for lab_ in selection:
+            # Extract feature of interest
+            if lab_==0:
+                data_ = np.logical_not(self.data)
+            else:
+                data_ = (self.labels[self._feat_slice[lab_-1]]==lab_)
+            # If feature is wrapped periodically, duplicate it and extract
+            # largest one
+            iswrapped = False
+            rep_ = self._ndim*[1]
+            for axis in range(self._ndim):
+                if self.wrap[axis] is not None and self.wrap[axis][lab_]:
+                    rep_[axis] = 2
+                    iswrapped = True
+            if iswrapped:
+                data_  = np.tile(data_,rep_)
+                l_,nl_ = ndimage.label(data_,structure=self.structure)
+                vol_   = np.bincount(l_.ravel())
+                il_    = np.argmax(vol_[1:])+1
+                sl_    = ndimage.find_objects(l_==il_)[0]
+                data_  = data_[sl_]
+            # Add to output
+            output.append(data_)
+        return output
+
+    def triangulate_feature(self,selection,origin=(0,0,0),spacing=(1,1,1)):
+        assert self._ndim==3, "Triangulation requires 3D data."
+        if self.labels is None:
+            self.label()
+        selection = self._select_feature(selection)
+        output = []
+        pw = tuple((1,1) for ii in range(self._ndim))
+        for lab_ in selection:
+            data_ = np.pad(self.isolate_feature(lab_)[0],pw,mode='constant')
+            output.append(Field3d(data_,origin,spacing).vtk_contour(0.5))
+        return output
 
     def _select_feature(self,selection):
-        dtype = self.label.dtype
-        if np.issubdtype(type(selection),np.integer):
-            return np.array(selection,dtype=dtype)
+        dtype = self.labels.dtype
+        if selection is None:
+            return np.array(range(1,self.nlabels+1))
+        elif np.issubdtype(type(selection),np.integer):
+            return np.array(selection,dtype=dtype,ndmin=1)
         elif isinstance(selection,(list,tuple,np.ndarray)):
             selection = np.array(selection)
             if selection.dtype==np.dtype('bool'):
-                assert selection.ndim==1 and selection.shape[0]==self.count+1,\
+                assert selection.ndim==1 and selection.shape[0]==self.nlabels+1,\
                         "Boolean indexing must provide count+1 values."
             else:
                 selection = selection.astype(dtype)
-                assert np.max(selection)<=self.count and np.min(selection)>=0,\
+                assert np.max(selection)<=self.nlabels and np.min(selection)>=0,\
                         "Entry in selection is out-of-bounds."
             return selection
         else: