From fdbe6fb70b2a4df3ce3aa5613cae417357d89b1c Mon Sep 17 00:00:00 2001
From: Michael Krayer <michael.krayer@posteo.de>
Date: Wed, 11 Aug 2021 11:31:02 +0200
Subject: [PATCH] WIP: find_region

---
 field.py | 89 +++++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 63 insertions(+), 26 deletions(-)

diff --git a/field.py b/field.py
index 6decfe7..955643c 100644
--- a/field.py
+++ b/field.py
@@ -771,11 +771,11 @@ class Features3d:
         # of the normal has to be chosen which defaults to the z-component and is set by 
         # 'cellvol_normal_component'.
         if report: print('[Features3d.triangulate] calculating area and volume per cell...')
-        X = points[faces[:,1],:]
-        Y = points[faces[:,2],:]
-        Z = points[faces[:,3],:]
-        cn = np.cross(X-Z,Y-X)
-        cc = (X+Y+Z)/3
+        U = points[faces[:,1],:]
+        V = points[faces[:,2],:]
+        W = points[faces[:,3],:]
+        cn = np.cross(U-W,V-U)
+        cc = (U+V+W)/3
         area = 0.5*np.sqrt(np.square(cn).sum(axis=1))
         vol  = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
         # Now the label is known per cell. We only need to find all cells with the same label
@@ -784,7 +784,6 @@ class Features3d:
         if report: print('[Features3d.triangulate] finalizing internal data...')
         offset = np.zeros((nfeatures+1,),dtype=np.int64)
         offset[1:] = np.cumsum(np.bincount(cell_labels))
-        print(nfeatures,np.bincount(cell_labels),offset)
         ind    = np.argsort(cell_labels)
         self._offset   = offset
         self._faces    = faces[ind,:]
@@ -797,14 +796,20 @@ class Features3d:
     @property
     def threshold(self): return self._threshold
 
+    @property
+    def nfeatures(self): return self._nfeatures
+
+    @property
+    def faces(self): return np.split(self._faces,self._offset)
+
     @property
     def points(self): return self._points
 
     @property
-    def faces(self): return self._faces
+    def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
 
     @property
-    def nfeatures(self): return self._nfeatures
+    def cell_volumes(self): return np.split(self._cell_volumes,self._offset[1:-1])
 
     def fill_holes(self,report=False):
         '''Remove triangulation which is fully enclosed by another. These regions 
@@ -851,18 +856,6 @@ class Features3d:
             self.clean_points(report=report)
         return
 
-    @property
-    def faces(self): return np.split(self._faces,self._offset)
-
-    @property
-    def points(self): return self._points
-
-    @property
-    def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
-
-    @property
-    def cell_volumes(self): return np.split(self._cell_volumes,self._offset[1:-1])
-
     def area(self,feature): 
         '''Returns the surface area of feature. If feature is None, total surface 
            area of all features is returned.'''
@@ -887,12 +880,56 @@ class Features3d:
         '''Returns an array with volumes of all features.'''
         return np.add.reduceat(self._cell_volumes,self._offset[:-1])
 
-    def find_feature(self,coords):
-        coords = np.array(coords)
-        assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as 3xN array."
-        # if method=='binary':
-        #     idx = np.round()
-        # elif method=='triangulation':
+    def find_feature(self):
+        # coords = np.array(coords)
+        # assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array."
+        from time import time
+        from scipy import spatial
+
+        t = time()
+
+        xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
+        xmax = np.amax(self._points[self._faces[:,1:],0],axis=1)
+        zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
+        zmax = np.amax(self._points[self._faces[:,1:],2],axis=1)
+
+        center = np.stack((0.5*(xmax+xmin),0.5*(zmax+zmin)),axis=1)
+        radius = np.sqrt(2.0)*np.maximum(
+            np.amax(0.5*(xmax-xmin)),
+            np.amax(0.5*(zmax-zmin)))
+        print('Preparation for KD tree in',time()-t)
+
+        t = time()
+        kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False)
+        # kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
+        # kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False)
+        # kd = spatial.KDTree(self._points[:,1:],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
+        print('KD tree built in',time()-t)
+
+        query = np.random.random((10000,2))
+
+        # t = time()
+        # map_ = self._points.shape[0]*[[]]
+        # for ii,face in enumerate(self._faces[:,1:].ravel()):
+        #     map_[face].append(ii)
+        # faces_connected_to_vertex = {}
+        # for face_index,face in enumerate(self._faces[:,1:]):
+        #     for vertex_index in face:
+        #         faces_connected_to_vertex.setdefault(vertex_index,[]).append(face_index)
+        # print('Inverted cells-faces',time()-t)
+
+        radius = np.sqrt(np.sum(self.spacing[1:]**2))
+        
+        t = time()
+        kd.query_ball_point(query,radius)
+        print('KD tree query',time()-t)
+
+        # print(radius)
+        # print(len(kd.query_ball_point((0.5,0.0),radius,return_sorted=True)),self._points.shape[0])
+
+
+        # yset = set(x.data for x in sorted(treey.at(0.0)))
+        # zset = set(x.data for x in sorted(treez.at(0.0)))
         return
 
     def clean_points(self,report=False):