From 3d6ecf2bda8c95c29c3a59483ea82f4e1930b97d Mon Sep 17 00:00:00 2001
From: Michael Krayer <michael.krayer@posteo.de>
Date: Thu, 19 Aug 2021 16:31:53 +0200
Subject: [PATCH] fixed label() for BinaryField

---
 field.py | 54 ++++++++++++++++++++++++++++++------------------------
 1 file changed, 30 insertions(+), 24 deletions(-)

diff --git a/field.py b/field.py
index 5f03cfd..c6cdf21 100644
--- a/field.py
+++ b/field.py
@@ -1004,7 +1004,7 @@ class Features3d:
             v      = -np.sum(E1*DAO)*invdet
             t      =  np.sum(AO*N)  *invdet
             return (np.abs(det)>=1e-7*norm and t>=0.0 and u>=0 and v>=0 and (u+v)-1.0<=0), t
-        # @jit(nopython=True,cache=True)
+        @jit(nopython=True,cache=True)
         def __get_nearest_face(coords,raydir,cand_arr,cand_num,faces,points):
             N = coords.shape[0]
             nearface = np.zeros(N,dtype=np.int64)
@@ -1100,7 +1100,8 @@ class Features3d:
         else:
             return list(features)
 
-    def cmap_features(self,nfeatures,name='tab10'):
+    @staticmethod
+    def cmap_features(nfeatures,name='tab10'):
         from matplotlib.colors import ListedColormap
         if name=='tab10': # seaborn/matplotlib tab10 
             cmap = [(0.12156862745098039, 0.4666666666666667, 0.7058823529411765), 
@@ -1123,33 +1124,35 @@ class Features3d:
             icolor = (icolor+1)%ncolor
         return ListedColormap(output)
 
-    def to_vtk(self,features,incl_regionid=False,incl_boundary=False):
+    def to_vtk(self,features,incl_boundary=False,incl_regionid=False,remap_regionid=True):
         import pyvista as pv
         faces  = self.faces_from_feature(features,incl_boundary=incl_boundary)
         output = pv.PolyData(self._points,faces)
         if incl_regionid:
             output.cell_arrays['RegionId'] = self.regionid_from_feature(features,incl_boundary=incl_boundary)
+            if remap_regionid:
+                output.cell_arrays['RegionId'] = np.unique(output.cell_arrays['RegionId'],return_inverse=True)[1]
         return output
 
-    def to_vtk2(self,features):
-        import pyvista as pv
-        idx = self.faces_from_feature(features)
-        faces_ = self._faces[idx,:]
-        A = self._points[faces_[:,1],:]
-        B = self._points[faces_[:,2],:]
-        C = self._points[faces_[:,3],:]
-        cn = np.cross(B-A,C-A)
-        cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True))
-        cc = (A+B+C)/3
-        # print(A[0,:],B[0,:],C[0,:],cn[0,:])
-        # print(A.shape,B.shape,C.shape,cn.shape,faces_.shape,self._faces.shape)
-        output = pv.PolyData(cc)
-        output['cellnormal'] = cn
-        return output
+    # def to_vtk2(self,features):
+    #     import pyvista as pv
+    #     idx = self.faces_from_feature(features)
+    #     faces_ = self._faces[idx,:]
+    #     A = self._points[faces_[:,1],:]
+    #     B = self._points[faces_[:,2],:]
+    #     C = self._points[faces_[:,3],:]
+    #     cn = np.cross(B-A,C-A)
+    #     cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True))
+    #     cc = (A+B+C)/3
+    #     # print(A[0,:],B[0,:],C[0,:],cn[0,:])
+    #     # print(A.shape,B.shape,C.shape,cn.shape,faces_.shape,self._faces.shape)
+    #     output = pv.PolyData(cc)
+    #     output['cellnormal'] = cn
+    #     return output
 
-    def vtk_faces(self,face_idx):
-        import pyvista as pv
-        return pv.PolyData(self._points,self._faces[face_idx,:])
+    # def vtk_faces(self,face_idx):
+    #     import pyvista as pv
+    #     return pv.PolyData(self._points,self._faces[face_idx,:])
 
 class BinaryFieldNd:
     def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):
@@ -1204,7 +1207,7 @@ class BinaryFieldNd:
         if not any(self.periodicity):
             return
         # Get a mapping of labels which differ at periodic overlap
-        map_  = np.array(range(nlabels_+1),dtype=labels_.dtype)
+        map_  = np.array(range(self.nlabels+1),dtype=self._labels.dtype)
         wrap_ = self._ndim*[None]
         for axis in range(self._ndim):
             if not self.periodicity[axis]: continue
@@ -1268,8 +1271,11 @@ class BinaryFieldNd:
                 labels_keep |= set(np.unique(labels_[sl_lo])) | set(np.unique(labels_[sl_hi]))
         labels_keep.discard(0)
         self._data = np.isin(labels_,tuple(labels_keep),invert=True)
-        # If labels have been computed already, recompute them to stay consistent
-        if self._labels is not None: self.label()
+        # If labels have been computed already, set them to None because they are now invalid
+        if self._labels is not None: 
+            self._labels = None
+            self.nlabels = None
+            self.wrap    = tuple(self._ndim*[None])
         # Compute mask if it is to be returned, otherwise we are done
         if return_mask:
             np.logical_xor(mask,self._data,out=mask)