From 0723c0b769c3a8b21afae87aaa70270d6f7ff509 Mon Sep 17 00:00:00 2001
From: Michael Krayer <michael.krayer@posteo.de>
Date: Fri, 13 Aug 2021 01:16:27 +0200
Subject: [PATCH] storing triangulation data as long list of vertices now.
 Seems to have improved performance due to better memory access pattern.

---
 field.py | 148 +++++++++++++++++++++++++------------------------------
 1 file changed, 68 insertions(+), 80 deletions(-)

diff --git a/field.py b/field.py
index 8a0da86..82af5f9 100644
--- a/field.py
+++ b/field.py
@@ -696,6 +696,7 @@ class Features3d:
         import pyvista as pv
         import vtk
         from scipy import ndimage, spatial
+        from time import time
         # Check if '_input' is available: might have been dropped after initialization
         assert self._input is not None, "'_input' not available. Initialize object with keep_input=True flag."
         # Wrap data for VTK using pyvista
@@ -707,22 +708,21 @@ class Features3d:
         # Compute full contour: ensure we only have triangles to efficiently extract points
         # later on.
         if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method))
-        contour_ = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False)
-        assert contour_.is_all_triangles(), "Contouring produced non-triangle cells."
+        contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False)
+        assert contour.is_all_triangles(), "Contouring produced non-triangle cells."
         # Compute the connectivity of the triangulated surface: first we run an ordinary 
         # connectivity filter neglecting periodic wrapping.
         if report: print('[Features3d.triangulate] computing connectivity...')
         alg = vtk.vtkPolyDataConnectivityFilter() # ~twice as fast as default vtkConnectivityFilter
-        alg.SetInputData(contour_)
+        alg.SetInputData(contour)
         alg.SetExtractionModeToAllRegions()
         alg.SetColorRegions(True)
         alg.Update()
-        contour_ = pv.filters._get_output(alg)
+        contour = pv.filters._get_output(alg)
         # Now determine the boundary points, i.e. points on the boundary which have a corresponding
         # vertex on the other side of the domain
         if report: print('[Features3d.triangulate] correcting connectivity for periodic boundaries...')
-        points = contour_.points
-        bd_points_xyz = np.zeros((0,3),dtype=points.dtype)
+        bd_points_xyz = np.zeros((0,3),dtype=contour.points.dtype)
         bd_points_idx = np.zeros((0,),dtype=np.int)
         for axis in range(3):
             if not self.periodicity[axis]: continue
@@ -732,19 +732,19 @@ class Features3d:
             period[axis]  = bound_pos-datavtk.origin[axis]
             bound_dist    = 1e-5*datavtk.spacing[axis]
             # Lower boundary
-            li = np.abs(points[:,axis]-datavtk.origin[axis])<bound_dist
+            li = np.abs(contour.points[:,axis]-datavtk.origin[axis])<bound_dist
             bd_points_idx = np.append(bd_points_idx,np.nonzero(li)[0],axis=0)
-            bd_points_xyz = np.append(bd_points_xyz,points[li,:],axis=0)
+            bd_points_xyz = np.append(bd_points_xyz,contour.points[li,:],axis=0)
             # Upper boundary
-            li = np.abs(points[:,axis]-bound_pos)<bound_dist
+            li = np.abs(contour.points[:,axis]-bound_pos)<bound_dist
             bd_points_idx = np.append(bd_points_idx,np.nonzero(li)[0],axis=0)
-            bd_points_xyz = np.append(bd_points_xyz,points[li,:]-period,axis=0)
+            bd_points_xyz = np.append(bd_points_xyz,contour.points[li,:]-period,axis=0)
         # Construct a KD Tree for efficient neighborhood search
         kd = spatial.KDTree(bd_points_xyz,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
         # Construct a map to get new labels for regions. We search for pairs of points with
         # a very small distance inbetween. Then we check if their labels differ and choose 
         # the lower label as the new joint one.
-        point_labels = contour_.point_arrays['RegionId']
+        point_labels = contour.point_arrays['RegionId']
         nfeatures    = np.max(point_labels)+1
         map_         = np.array(range(nfeatures),dtype=point_labels.dtype)
         overlap_dist = 1e-4*np.sqrt(np.square(datavtk.spacing).sum())
@@ -760,37 +760,35 @@ class Features3d:
         # 
         map_         = np.unique(map_,return_inverse=True)[1]
         point_labels = map_[point_labels]
-        nfeatures    = np.max(map_)+1 # starts with zero
+        self._nfeatures = np.max(map_)+1 # starts with zero
         # Labels are now stored as point data. To efficiently convert it to cell data, the first 
         # point of each cell determines the value for this cell.
         if report: print('[Features3d.triangulate] identifying cell based labels...')
-        ncells = contour_.n_faces
-        faces  = contour_.faces.reshape(ncells,4)[:,:]
-        cell_labels = point_labels[faces[:,1]]
+        cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
+        # We are done with VTK for now. Since we are only dealing with triangles, let us
+        # store them in a structure which is quicker to access and already sorted based 
+        # on labels.
+        # _vertices has dim(ncells,nvert,ndim), where nvert=3 (three vertices per triangle)
+        # and ndim=3 (three dimensional coordinates)
+        # _offset is the cumulative number of cells for each feature and the offset to
+        # access the _vertices array by feature.
+        if report: print('[Features3d.triangulate] storing vertices...')
+        ind = np.argsort(cell_labels)
+        self._offset = np.zeros((self._nfeatures+1,),dtype=np.int64)
+        self._offset[1:] = np.cumsum(np.bincount(cell_labels))
+        self._vertices = np.ascontiguousarray(
+                            contour.points[contour.faces.reshape(contour.n_faces,4)[ind,1:]])
         # Compute the volume and area per cell. For the volume computation, an arbitrary component 
         # 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...')
-        A = points[faces[:,1],:]
-        B = points[faces[:,2],:]
-        C = points[faces[:,3],:]
+        A = self._vertices[:,0,:]
+        B = self._vertices[:,1,:]
+        C = self._vertices[:,2,:]
         cn = np.cross(B-A,C-A)
         cc = (A+B+C)/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
-        # and group them. Therefore we sort the array holding the faces by feature label
-        # and save the offset in another array to easily extract them whenever necessary.
-        if report: print('[Features3d.triangulate] finalizing internal data...')
-        offset = np.zeros((nfeatures+1,),dtype=np.int64)
-        offset[1:] = np.cumsum(np.bincount(cell_labels))
-        ind    = np.argsort(cell_labels)
-        self._offset   = offset
-        self._faces    = faces[ind,:]
-        self._points   = points
-        self._cell_areas   = area[ind]
-        self._cell_volumes = vol[ind]
-        self._nfeatures = nfeatures
+        self._cell_areas   = 0.5*np.sqrt(np.square(cn).sum(axis=1))
+        self._cell_volumes = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
         return
 
     @property
@@ -799,12 +797,6 @@ class Features3d:
     @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 cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
 
@@ -822,7 +814,7 @@ class Features3d:
         if is_relative:
             vol_domain = np.prod(self.spacing*self.dimensions)
             threshold = threshold*vol_domain
-        self.discard_features(np.flatnonzero(self.volumes()<threshold),report=report)
+        self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
         return
 
     def discard_features(self,labels,clean_points=False,report=False):
@@ -840,9 +832,9 @@ class Features3d:
             gapsize[label] = len(rng_)
         idx = np.concatenate(idx,axis=0)
         # Save former number of faces for report
-        nfaces = self._faces.shape[0]
+        ncells = self._vertices.shape[0]
         # Delete indexed elements from arrays
-        self._faces        = np.delete(self._faces,idx,axis=0)
+        self._vertices     = np.delete(self._vertices,idx,axis=0)
         self._cell_areas   = np.delete(self._cell_areas,idx,axis=0)
         self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
         # Correct offset
@@ -851,9 +843,7 @@ class Features3d:
         if report:
             print('[Features3d.discard_features]',end=' ')
             print('discarded {:} features with {:} faces.'.format(
-                            len(labels),nfaces-self._faces.shape[0]))
-        if clean_points:
-            self.clean_points(report=report)
+                            len(labels),ncells-self._vertices.shape[0]))
         return
 
     def area(self,feature): 
@@ -862,7 +852,7 @@ class Features3d:
         if feature is None: 
             return np.sum(self._cell_areas)
         else: 
-            return np.sum(self.cell_areas[feature-1])
+            return np.sum(self._cell_areas[self._offset[feature:feature+2]])
 
     def areas(self): 
         '''Returns an array with surface areas of all features.'''
@@ -874,7 +864,7 @@ class Features3d:
         if feature is None: 
             return np.sum(self._cell_volumes)
         else: 
-            return np.sum(self.cell_volumes[feature-1])
+            return np.sum(self._cell_volumes[self._offset[feature:feature+2]])
 
     def volumes(self):
         '''Returns an array with volumes of all features.'''
@@ -907,7 +897,6 @@ class Features3d:
 
         # pts = np.random.random((10000,3))
 
-        print(self._points[self._faces[:,1:],0].shape)
         t = time()
 
         # xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
@@ -915,8 +904,8 @@ class Features3d:
         # zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
         # zmax = np.amax(self._points[self._faces[:,1:],2],axis=1)
         # print(time()-t)
-        xmin,xmax = minmax(self._points[self._faces[:,1:],0])
-        zmin,zmax = minmax(self._points[self._faces[:,1:],2])
+        xmin,xmax = minmax(self._vertices[:,:,0])
+        zmin,zmax = minmax(self._vertices[:,:,2])
         print(time()-t)
         # print(xmin.shape,xmin2.shape)
         # print(np.all(np.isclose(xmin,xmin2)))
@@ -957,9 +946,9 @@ class Features3d:
         print(is_inside.shape)
         for ii in range(Npts):
             hit_idx,t,N = Features3d.ray_triangle_intersection(pts[ii],dR,
-                                    self._points[self._faces[bla[ii],1]],
-                                    self._points[self._faces[bla[ii],2]],
-                                    self._points[self._faces[bla[ii],3]])
+                                    self._vertices[bla[ii],0,:],
+                                    self._vertices[bla[ii],1,:],
+                                    self._vertices[bla[ii],2,:])
             if hit_idx is None:
                 is_inside[ii] = False
             else:
@@ -1019,37 +1008,36 @@ class Features3d:
         t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
         return hit_idx,t,N[hit_idx]
 
-    def clean_points(self,report=False):
-        nfaces_ = self._faces.shape[0]
-        ind,inv = np.unique(self._faces[:,1:],return_inverse=True)
-        if report:
-            print('[Features3d.clean_points]',end=' ') 
-            print('removed {:} orphan points.'.format(self._points.shape[0]-len(ind)))
-        self._faces[:,1:4] = inv.reshape(nfaces_,3)
-        self._points = self._points[ind,:]
+    def faces_points_representation(self,labels,reduce_points=False):
+        from collections.abc import Iterable
+        # Select relevant cells
+        if labels is None:
+            idx = None
+        else:
+            if not isinstance(labels,Iterable): labels = [labels]
+            idx = []
+            for label in labels:
+                rng_ = np.arange(self._offset[label],self._offset[label+1])
+                idx.append(rng_)
+            idx = np.concatenate(idx,axis=0)
+        if reduce_points:
+            points,faces_ = np.unique(self._vertices[idx,:,:].reshape(-1,3),return_inverse=True,axis=0)
+            ncells = faces_.shape[0]//3
+            faces = np.empty((ncells,4),dtype=np.int)
+            faces[:,0] = 3
+            faces[:,1:] = faces_.reshape(ncells,3)
+        else:
+            points = self._vertices[idx,:,:].reshape((-1,3))
+            ncells = points.shape[0]//3
+            faces = np.empty((ncells,4),dtype=np.int)
+            faces[:,0] = 3
+            faces[:,1:] = np.arange(0,3*ncells).reshape(ncells,3)
+        return faces,points
 
     def to_vtk(self,labels,reduce_points=False):
         import pyvista as pv
         from collections.abc import Iterable
-        if labels is None:
-            return pv.PolyData(self._points,self._faces)
-        if not isinstance(labels,Iterable): 
-            labels = [labels]
-        points = self._points
-        faces  = []
-        for label in labels:
-            sl_ = slice(self._offset[label],self._offset[label+1])
-            faces.append(self._faces[sl_,:])
-        faces = np.concatenate(faces,axis=0)
-        # if reduce_points:
-        #     nfaces_ = self._faces[idx].shape[0]
-        #     ind,inv = np.unique(self._faces[:,1:4],return_inverse=True)
-        #     faces_  = np.zeros((nfaces_,4),dtype=self._faces.dtype)
-        #     faces_[:,0]   = 3
-        #     faces_[:,1:4] = inv.reshape(nfaces_,3)
-        #     points_ = self._points[ind,:]
-        #     return pv.PolyData(points_,faces_)
-        # else:
+        faces,points = self.faces_points_representation(labels,reduce_points=reduce_points)
         return pv.PolyData(points,faces)
 
 class BinaryFieldNd: