WIP: completely messed up. trying to fix surface normal directions

field.py

@@ -659,7 +659,7 @@ class Features3d:
         # Assign basic properties to class variables
         self.origin      = np.array(origin,dtype=np.float)
         self.spacing     = np.array(spacing,dtype=np.float)
-        self.dimensions  = input.shape
+        self.dimensions  = np.array(input.shape,dtype=np.int)
         self.periodicity = tuple(bool(x) for x in periodicity)
         # If regions are supposed to be inverted, i.e. the interior consists of values 
         # smaller than the threshold instead of larger, change the sign of the array.
@@ -712,8 +712,14 @@ 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)
+        contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False,compute_gradients=True)
         assert contour.is_all_triangles(), "Contouring produced non-triangle cells."
+
+        # print(contour.point_arrays['Gradients'].shape)
+        # print('gradient',time()-t)
+        # print(np.sum(gn>0),'of',len(gn))
+
+
         # Compute the connectivity of the triangulated surface: first we run an ordinary 
         # connectivity filter neglecting periodic wrapping.
         if report: print('[Features3d.triangulate] computing connectivity...')
@@ -769,6 +775,8 @@ class Features3d:
         # point of each cell determines the value for this cell.
         if report: print('[Features3d.triangulate] identifying cell based labels...')
         cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
+        cell_gradient = contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1]]
+        print(cell_gradient.shape)
         # 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.
@@ -780,8 +788,7 @@ class Features3d:
         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(
+        self._vertices = contour.points[contour.faces.reshape(contour.n_faces,4)[ind,1:]]
-                            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'.
@@ -791,8 +798,39 @@ class Features3d:
         C = self._vertices[:,2,:]
         cn = np.cross(B-A,C-A)
         cc = (A+B+C)/3
+        # print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]].shape)
+        cg = (contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:]+
+              contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,1,:]+
+              contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,2,:])/3
+        ng = (cg*cn).sum(axis=-1)
+        print('interpolated:')
+        print(cg)
+        print('component 0:')
+        print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:])
+        print('projection')
+        print(sum(ng>0),len(ng))
+
+        # self.cg = cg
+        # self.cc = cc
+        self.cg = contour.point_arrays['Gradients']
+        self.cc = contour.points
+
+        # .sum(axis=1))/3)*cn).sum(axis=-1)
         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]
+
+        # # Compute gradient of field at cell center to get correct cell normal direction
+        # t = time()
+        # cn = cn/np.sqrt(np.square(cn).sum())
+        # ccpix  = ((cc-self.origin)/self.spacing).transpose()
+        # ccnpix = (((cc+1e-4*cn)-self.origin)/self.spacing).transpose()
+        # val_cc  = ndimage.map_coordinates(self._input,ccpix,order=1,prefilter=False)
+        # val_ccn = ndimage.map_coordinates(self._input,ccnpix,order=1,prefilter=False)
+        # gn = (val_ccn-val_cc)
+        # print('gradient',time()-t)
+        # print(cg)
+        # print(np.sum(cg>0),'of',cg.shape)
+        # stop
         return
 
     @property
@@ -821,19 +859,17 @@ class Features3d:
         self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
         return
 
-    def discard_features(self,labels,clean_points=False,report=False):
+    def discard_features(self,features,report=False):
         '''Removes features from triangulated data.'''
-        from collections.abc import Iterable
+        features = self.list_of_features(features)
-        from time import time
-        if not isinstance(labels,Iterable): labels = [labels]
         # Get index ranges which are to be deleted and also create an array
         # which determines the size of the block to be deleted
         idx = []
-        gapsize = np.zeros((self.nfeatures,),dtype=np.int)
+        gapsize = np.zeros((self._nfeatures,),dtype=np.int)
-        for label in labels:
+        for feature in features:
-            rng_ = np.arange(self._offset[label],self._offset[label+1])
+            rng_ = np.arange(self._offset[feature],self._offset[feature+1])
             idx.append(rng_)
-            gapsize[label] = len(rng_)
+            gapsize[feature] = len(rng_)
         idx = np.concatenate(idx,axis=0)
         # Save former number of faces for report
         ncells = self._vertices.shape[0]
@@ -843,11 +879,12 @@ class Features3d:
         self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
         # Correct offset
         self._offset[1:] = self._offset[1:]-np.cumsum(gapsize)
-        self._offset = self._offset[:-len(labels)]
+        self._offset     = np.delete(self._offset,features,axis=0)
+        self._nfeatures  = self._offset.shape[0]-1
         if report:
             print('[Features3d.discard_features]',end=' ')
             print('discarded {:} features with {:} faces.'.format(
-                            len(labels),ncells-self._vertices.shape[0]))
+                            len(features),ncells-self._vertices.shape[0]))
         return
 
     def area(self,feature): 
@@ -902,9 +939,9 @@ class Features3d:
         self._kdradius = radius
         if report:
             print('[Features3d.build_kdtree]',end=' ')
-            print('KD-Tree built for axis {}.'.format(axis))
+            print('KD-Tree built for axis {}.'.format(kdaxis))
 
-    def inside_feature(self,coords):
+    def inside_feature(self,coords,report=True):
         '''Find nearest triangle from point R in direction ±dR and its orientation w.r.t dR. 
            The algorithm determines the closest intersection of a ray cast from the origin of the 
            point in the specified direction (with positiv and negative sign). The ray is supposed to 
@@ -934,29 +971,133 @@ class Features3d:
 
         cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
 
-        N = coords.shape[0]
+        if coords.shape[0]>12279:
+            print('Point coordinates:',coords[12279,:])
+
+        Ncoord = coords.shape[0]
         raydir = np.zeros((3,),dtype=self._vertices.dtype)
         raydir[self._kdaxis] = 1.0
-        in_feat = np.empty((N,),dtype=np.int)
+        in_feat  = np.empty((Ncoord,),dtype=np.int)
-        for ii in range(N):
+        is_hit   = np.empty((Ncoord,),dtype=np.int)
-            hit_idx,t,N = Features3d.ray_triangle_intersection(coords[ii],raydir,
+        hit_dir_ = np.empty((Ncoord,),dtype=np.int)
+        face_ = np.empty((Ncoord,),dtype=np.int)
+        for ii in range(Ncoord):
+            hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir,
                                                     self._vertices[cand[ii],0,:],
                                                     self._vertices[cand[ii],1,:],
                                                     self._vertices[cand[ii],2,:])
+            # if ii==12279: print('DEBUG',hit_idx,t,N)
             if hit_idx is None:
                 in_feat[ii] = -1
+                is_hit[ii]  = -1
+                hit_dir_[ii]=  0
+                face_[ii]   = -1
             else:
                 idx = np.argmin(np.abs(t))
-                if t[idx]*(N[idx,:]*raydir).sum(axis=-1)>0:
+                # print(cand[ii][hit_idx[idx]])
-                    in_feat[ii] = self.feature_by_cell(cand[ii][idx])
+                if ii==12279: 
+                    # print(cand[ii][idx],N[idx,:],N[idx,:]/np.sqrt(np.square(N[idx,:]).sum()),(N[idx,:]*raydir).sum(axis=-1),t[idx]*(N[idx,:]*raydir).sum(axis=-1))
+                    print(t[idx],hit_dir[idx],raydir)
+                    print(self._vertices[cand[ii][hit_idx[idx]],:,:])
+                hit_dir_[ii]= hit_dir[idx]
+                face_[ii]   = cand[ii][hit_idx[idx]]
+                if hit_dir[idx]<0:
+                    in_feat[ii] = self.feature_from_cellidx(cand[ii][hit_idx[idx]])
+                    is_hit[ii]  = 1
                 else:
                     in_feat[ii] = -1
-        return in_feat
+                    is_hit[ii]  = 0
+        if report:
+            print('[Features3d.inside_feature]',end=' ')
+            print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord))
+        return in_feat #,is_hit,hit_dir_,face_
 
-    def feature_by_cell(self,idx_cell):
+    def feature_from_cellidx(self,idx_cell):
         '''Gets feature ID for a given cell.'''
-        return np.argmax(idx_cell>self._offset)
+        # This is an optimized solution for monotonically increasing arrays: 
+        # as long as offset is smaller than the cell index, the boolean operation
+        # returns False. As soon as the first True is returned, i.e.
+        # as soon as idx_cell is equal or larger than offset, the argmax function
+        # short circuits and returns the index of the first occurence without 
+        # checking any value afterwards.
+        return np.argmax(self._offset>=idx_cell)
 
+    def cellidx_from_feature(self,features,concat=True):
+        '''Returns indices of cells which belong to given features.'''
+        from collections.abc import Iterable
+        if features is None:
+            idx = None
+        elif not isinstance(features,Iterable):
+            idx = slice(self._offset[features],self._offset[features+1])
+        elif concat:        
+            idx = []
+            for feature in features:
+                rng_ = np.arange(self._offset[feature],self._offset[feature+1])
+                idx.append(rng_)
+            idx = np.concatenate(idx,axis=0)
+        else:
+            idx = []
+            for feature in features:
+                rng_ = slice(self._offset[feature],self._offset[feature+1])
+                idx.append(rng_)
+        return idx
+
+    def ncells_from_feature(self,features):
+        '''Returns number of cells which belong to given features. Can be used
+           to construct new offsets.'''
+        features = self.list_of_features(features)
+        ncells = []
+        for feature in features:
+            ncells.append(self._offset[feature+1]-self._offset[feature])
+        return np.array(ncells)
+
+    def list_of_features(self,features):
+        '''Ensures that 'features' is a list.'''
+        from collections.abc import Iterable
+        if features is None:
+            return np.arange(0,self._nfeatures)
+        if not isinstance(features,Iterable): 
+            return [features]
+        else:
+            return list(features)
+
+    # @staticmethod
+    # def ray_triangle_intersection(r0,dr,v0,v1,v2):
+    #     '''Implements the Möller–Trumbore ray-triangle intersection algorithm. I modified the 
+    #        formulation of 
+    #        https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d
+    #        because it computes the cell normal on the way, which is needed to determine the 
+    #        direction of the hit, i.e. from the inside or outside.
+    #        Input:
+    #             R, dR: origin and direction of the ray as (3,) numpy arrays
+    #             A,B,C: vertices of N triangles as (N,3) numpy arrays
+    #        Returns:
+    #             hit_idx: index of the input triangles which returned a hit. [(Nhit,) int]
+    #             t:       parameter to determine intersection point (x = R+t*dR) [(Nhit,) float]
+    #             N:       normal vector of triangles which were hit [(Nhit,3) float]
+    #        All returned values are None if not hit occured.
+    #     '''
+    #     v0v1 = v1-v0
+    #     v0v2 = v2-v0
+    #     pvec = np.cross(dr,v0v2)
+    #     det  = (v0v1*pvec).sum(axis=-1) # det<0: from behind, det>0 from front ("culling")
+    #     norm = 1e0 # we should normalize the unit vector?
+    #     det[np.abs(det)<1e-6*norm] = np.nan # mask to avoid numpy runtime errors
+    #     invDet = 1./det 
+    #     tvec = r0-v0
+    #     qvec = np.cross(tvec,v0v1)
+    #     u = (tvec*pvec).sum(axis=-1)*invDet
+    #     v = (dr*qvec).sum(axis=-1)*invDet
+    #     is_hit = np.logical_and(
+    #                 np.logical_and(
+    #                     np.logical_and(
+    #                         np.isfinite(det),u>=-1e-6),
+    #                     v>=1e-6),
+    #                 (u+v)-1.0<=1e-6)
+    #     hit_idx = np.flatnonzero(is_hit)
+    #     if len(hit_idx)==0: return (None,None,None)
+    #     t = (v0v2[hit_idx,:]*qvec[hit_idx,:]).sum(axis=-1)*invDet[hit_idx]; 
+    #     return hit_idx,t,np.sign(det[hit_idx])
     @staticmethod
     def ray_triangle_intersection(R,dR,A,B,C):
         '''Implements the Möller–Trumbore ray-triangle intersection algorithm. I modified the 
@@ -988,46 +1129,75 @@ class Features3d:
         is_hit = np.logical_and(
                 np.logical_and(
                     np.logical_and(
-                    np.isfinite(det),u>=0.0),
+                    np.isfinite(det),u>=-1e-6),
-                    v>=0.0),
+                    v>=-1e-6),
-                (u+v)<=1.0)
+                (u+v)-1.0<=1e-6)
         hit_idx = np.flatnonzero(is_hit)
         if len(hit_idx)==0: return (None,None,None)
         # Intersection point is R+t*dR
         t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
-        return hit_idx,t,N[hit_idx]
+        return hit_idx,t,np.sign(t)*np.sign(det[hit_idx])
 
-    def faces_points_representation(self,labels,reduce_points=False):
+    def faces_points_representation(self,features,reduce_points=False):
         from collections.abc import Iterable
         # Select relevant cells
-        if labels is None:
+        idx = self.cellidx_from_feature(features)
-            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)
+            print('reduced:',faces.shape)
         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)
+            print('full:',faces.shape)
         return faces,points
 
-    def to_vtk(self,labels,reduce_points=False):
+    def cmap_features(self,nfeatures,name='tab10'):
+        from matplotlib.colors import ListedColormap
+        if name=='tab10': # seaborn/matplotlib tab10 
+            cmap = [(0.12156862745098039, 0.4666666666666667, 0.7058823529411765), 
+                    (1.0, 0.4980392156862745, 0.054901960784313725), 
+                    (0.17254901960784313, 0.6274509803921569, 0.17254901960784313), 
+                    (0.8392156862745098, 0.15294117647058825, 0.1568627450980392), 
+                    (0.5803921568627451, 0.403921568627451, 0.7411764705882353), 
+                    (0.5490196078431373, 0.33725490196078434, 0.29411764705882354), 
+                    (0.8901960784313725, 0.4666666666666667, 0.7607843137254902), 
+                    (0.4980392156862745, 0.4980392156862745, 0.4980392156862745), 
+                    (0.7372549019607844, 0.7411764705882353, 0.13333333333333333), 
+                    (0.09019607843137255, 0.7450980392156863, 0.8117647058823529)]
+        else:
+            raise NotImplementedError('Colormap not implemented.')
+        ncolor = len(cmap)
+        icolor = 0
+        output = np.empty((nfeatures,3))
+        for ii in range(nfeatures):
+            output[ii] = cmap[icolor]
+            icolor = (icolor+1)%ncolor
+        return ListedColormap(output)
+
+    def to_vtk(self,features,reduce_points=False,store_labels=False,remap_labels=False):
+        '''Returns a VTK/pyvista object for the isocontour of a single feature,
+           list of features, or the full isocontour (when None is passed).'''
         import pyvista as pv
-        from collections.abc import Iterable
+        faces,points = self.faces_points_representation(features,reduce_points=reduce_points)
-        faces,points = self.faces_points_representation(labels,reduce_points=reduce_points)
+        output = pv.PolyData(points,faces)
-        return pv.PolyData(points,faces)
+        if store_labels:
+            output.cell_arrays['label'] = np.empty(faces.shape[0],dtype=np.int)
+            offset = 0
+            for ii,feature in enumerate(self.list_of_features(features)):
+                ncells = self.ncells_from_feature(feature)
+                if remap_labels:
+                    output.cell_arrays['label'][offset:offset+ncells] = ii
+                else:
+                    output.cell_arrays['label'][offset:offset+ncells] = feature
+                offset += ncells
+        return output
 
 class BinaryFieldNd:
     def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):