field.py

@@ -650,15 +650,16 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0):
 class Features3d:
     def __init__(self,input,threshold,origin,spacing,periodicity,
                         invert=False,has_ghost=False,keep_input=False,
-                        contour_method='flying_edges',
+                        contour_method='flying_edges',cellvol_normal_component=2,
                         report=False):
         assert len(origin)==3,  "'origin' must be of length 3"
         assert len(spacing)==3, "'spacing' must be of length 3"
         assert len(periodicity)==3, "'periodicity' must be of length 3"
         assert isinstance(input,np.ndarray), "'input' must be numpy.ndarray."
         # Assign basic properties to class variables
+        self.origin      = np.array(origin,dtype=np.float)
         self.spacing     = np.array(spacing,dtype=np.float)
-        self.origin      = np.array(origin,dtype=np.float)-self.spacing
+        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.
@@ -667,20 +668,15 @@ class Features3d:
         # '_input' is the array which is to be triangulated. Here one trailing ghost cell is 
         # required to get closed surfaces. The data will always be copied since it probably 
         # needs to be copied for vtk anyway (needs FORTRAN contiguous array) and this way
-        # there will never be side effects on the input data. The array is padded with
-        # a very high value (cannot use NaN or Inf for vtk) in order to close the contour at
-        # the boundaries.
+        # there will never be side effects on the input data.
         if has_ghost:
-            self.dimensions = input.shape+np.array((2,2,2))
-            self._input = np.full(self.dimensions,-1e30,dtype=input.dtype,order='F')
-            self._input[1:-1,1:-1,1:-1] = sign_invert*input
+            self._input = sign_invert*input
         else:
             pw = tuple((0,1) if x else (0,0) for x in periodicity)
-            self.dimensions = input.shape+np.array((2,2,2))+np.array([sum(x) for x in pw])
-            self._input = np.full(self.dimensions,-1e30,dtype=input.dtype,order='F')
-            self._input[1:-1,1:-1,1:-1] = np.pad(sign_invert*input,pw,mode='wrap')
+            self._input = np.pad(sign_invert*input,pw,mode='wrap')
         # Triangulate
         self.triangulate(contour_method=contour_method,
+                         cellvol_normal_component=cellvol_normal_component,
                          report=report)
         # Set some state variable
         self._kdtree   = None
@@ -693,145 +689,108 @@ class Features3d:
     @classmethod
     def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False,
                             keep_input=False,contour_method='flying_edges',
-                            correct_normals=True,report=False):
+                            cellvol_normal_component=2,report=False):
         return cls(fld.data,threshold,fld.origin,fld.spacing,periodicity,
                         invert=invert,has_ghost=has_ghost,keep_input=keep_input,
                         contour_method=contour_method,
+                        cellvol_normal_component=cellvol_normal_component,
                         report=report)
 
-    def triangulate(self,contour_method='flying_edges',report=False):
+    def triangulate(self,contour_method='flying_edges',cellvol_normal_component=2,report=False):
         import pyvista as pv
         import vtk
         from scipy import ndimage, spatial
-        from scipy.spatial import KDTree
-        from numba import jit
         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
         datavtk = pv.UniformGrid()
-        datavtk.dimensions = self.dimensions
+        datavtk.dimensions = self._input.shape
         datavtk.origin     = self.origin
         datavtk.spacing    = self.spacing
         datavtk.point_arrays['data'] = self._input.ravel('F')
-        # Compute full contour: due to the padding, the contour will result in closed object.
-        # This is necessary for proper inside/outside checks and volume computation. However,
-        # for surface area computation, the boundary faces have to be removed. Therefore, we 
-        # will isolate them and store them in an extra array.
+        # 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,compute_gradients=True)
         assert contour.is_all_triangles(), "Contouring produced non-triangle cells."
-        points    = contour.points
-        faces     = contour.faces.reshape(-1,4)
-        gradients = contour.point_arrays['Gradients'][faces[:,1],:]
-        # Python for loops are horribly slow, but for loops are the easiest way to perform
-        # some of the necessary check. Therefore let's define some functions to be jit compiled
-        # by numba to speed things up.
-        @jit(nopython=True)
-        def __get_boundary_faces(faces,points,bounds,tolerance):
-            assert faces.shape[1]==4 
-            assert points.shape[1]==3
-            assert len(bounds)==6
-            assert len(tolerance)==3
-            nfaces = faces.shape[0]
-            output = np.full((nfaces,),-1,dtype=np.int8)
-            for ii in range(nfaces):
-                for jj in range(3):
-                    if   (points[faces[ii,1],jj]<bounds[2*jj]+tolerance[jj] and 
-                          points[faces[ii,2],jj]<bounds[2*jj]+tolerance[jj] and
-                          points[faces[ii,3],jj]<bounds[2*jj]+tolerance[jj]):
-                        output[ii] = 2*jj
-                    elif (points[faces[ii,1],jj]>bounds[2*jj+1]-tolerance[jj] and 
-                          points[faces[ii,2],jj]>bounds[2*jj+1]-tolerance[jj] and
-                          points[faces[ii,3],jj]>bounds[2*jj+1]-tolerance[jj]):
-                        output[ii] = 2*jj+1
-            return output
-        @jit(nopython=True)
-        def __connect_faces_periodic(face_uni_lo,face_uni_hi,points,dist,idx,search_dist):
-            N = len(face_uni_hi)
-            output = np.zeros((N,4),dtype=face_uni_lo.dtype)
-            jj = 0
-            for ii in range(N):
-                if dist[ii]<search_dist:
-                    output[jj,:] = (3,face_uni_hi[ii],face_uni_lo[idx[ii]],face_uni_hi[ii])
-                    jj += 1
-            return output[:jj,:]
-        # In order to treat periodicity when labeling the various regions, we connect the 
-        # overlapping regions with degenerate triangles. Then vtkPolyDataConnectivityFilter
-        # returns the correct result right away, and we will just ignore the trailing 
-        # virtual faces afterwards.
-        tol_overlap = 1e-5*self.spacing
-        bd_face_tag = __get_boundary_faces(faces,points,tuple(contour.bounds),tol_overlap)
-        faces_connect = [faces]
-        for axis in range(3):
-            if not self.periodicity[axis]: continue
-            sl_pln = [0,1,2]
-            del sl_pln[axis]
-            face_uni_lo = np.unique(faces[bd_face_tag==2*axis,1:].ravel())
-            face_uni_hi = np.unique(faces[bd_face_tag==2*axis+1,1:].ravel())
-            kdt = KDTree(points[np.ix_(face_uni_lo,sl_pln)],leafsize=10,compact_nodes=False,
-                                                        copy_data=False,balanced_tree=False)
-            dist,idx = kdt.query(points[np.ix_(face_uni_hi,sl_pln)],k=1,
-                                                    distance_upper_bound=tol_overlap[axis])
-            faces_connect += [__connect_faces_periodic(face_uni_lo,face_uni_hi,points,
-                                                       dist,idx,tol_overlap[axis])]
-        # Save original number of faces and add the virtual connecting faces to the polydata.
-        nfaces = faces.shape[0]
-        contour.faces = np.concatenate(faces_connect,axis=0)
-        # Compute the connectivity using pure VTK (pyvista only supports vtkConnectivityFilter
-        # which takes around twice as long).
-        alg = vtk.vtkPolyDataConnectivityFilter()
+        # 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.SetExtractionModeToAllRegions()
         alg.SetColorRegions(True)
         alg.Update()
         contour = pv.filters._get_output(alg)
-        # RegionIds are now stored as point data. To efficiently convert it to cell data, the first 
+        # 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...')
+        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
+            # Compute position of boundary, period of domain and set a threshold for "on boundary" condition
+            bound_pos     = datavtk.origin[axis]+datavtk.spacing[axis]*(datavtk.dimensions[axis]-1)
+            period        = np.zeros((3,))
+            period[axis]  = bound_pos-datavtk.origin[axis]
+            bound_dist    = 1e-5*datavtk.spacing[axis]
+            # Lower boundary
+            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,contour.points[li,:],axis=0)
+            # Upper boundary
+            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,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']
+        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())
+        for (ii,jj) in kd.query_pairs(r=overlap_dist):
+            label_ii = point_labels[bd_points_idx[ii]]
+            label_jj = point_labels[bd_points_idx[jj]]
+            if label_ii!=label_jj:
+                source_ = np.maximum(label_ii,label_jj)
+                target_ = np.minimum(label_ii,label_jj)
+                while target_ != map_[target_]: # map it recursively
+                    target_ = map_[target_]
+                map_[source_] = target_
+        map_ = np.unique(map_,return_inverse=True)[1]
+        point_labels = map_[point_labels]
+        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.
-        regionid = contour.point_arrays['RegionId'][contour.faces.reshape(-1,4)[:nfaces,1]].ravel()
-        nfeatures = np.amax(regionid)+1
-        # Store the face/vertex data in class arrays. Sort the faces by RegionId for quick access
-        # later on.
         if report: print('[Features3d.triangulate] sorting faces by labels...')
-        self._points = np.array(points)
-        self._nfeatures = nfeatures
-        idx   = np.flatnonzero(bd_face_tag<0)
-        idxbd = np.flatnonzero(bd_face_tag>=0)
-        self._offset                = np.zeros((2*nfeatures+1,),dtype=np.int64)
-        self._offset[1:nfeatures+1] = np.cumsum(np.bincount(regionid[idx],minlength=nfeatures))
-        self._offset[nfeatures+1:]  = np.cumsum(np.bincount(regionid[idxbd],minlength=nfeatures))
-        self._offset[nfeatures+1:] += self._offset[nfeatures]
-        self._faces = np.concatenate((
-                        faces[idx,:][np.argsort(regionid[idx]),:],
-                        faces[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0)
-        gradients = np.concatenate((
-                        gradients[idx,:][np.argsort(regionid[idx]),:],
-                        gradients[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0)
-        # Compute the volume and area per cell.
-        if report: print('[Features3d.triangulate] calculating surface area and volume...')
+        cell_labels  = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
+        offset       = np.zeros((self._nfeatures+1,),dtype=np.int64)
+        offset[1:]   = np.cumsum(np.bincount(cell_labels))
+        ind          = np.argsort(cell_labels)
+        self._offset = offset
+        self._faces  = contour.faces.reshape(contour.n_faces,4)[ind,:]
+        self._points = contour.points
+        # 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 = self._points[self._faces[:,1],:]
         B = self._points[self._faces[:,2],:]
         C = self._points[self._faces[:,3],:]
         cn = np.cross(B-A,C-A)
         # Check if cell normal points in direction of gradient. If not, switch vertex order.
-        idx = (gradients*cn).sum(axis=-1)>0
-        self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])]
-        cn[idx,:] = -cn[idx,:]
+        idx = (contour.point_arrays['Gradients'][self._faces[:,1],:]*cn).sum(axis=-1)>0
+        # print(idx.shape,np.sum(idx),self._faces.shape,self._faces[idx,2:].shape,self._faces[idx,3:1:-1].shape)
+        # self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])]
+        # cn[idx,:] = -cn[idx,:]
         # Compute area and signed volume per cell
         cc = (A+B+C)/3
-        cell_areas    = 0.5*np.sqrt(np.square(cn).sum(axis=1))
-        cell_volumes  = 0.5*np.mean(cn*cc,axis=1)
-        @jit(nopython=True)
-        def __sum_up_cell_data(cell_data,offset,incl_boundary):
-            nfeat = (len(offset)-1)//2
-            niter = 2*nfeat if incl_boundary else nfeat
-            output = np.zeros(nfeat,dtype=cell_data.dtype)
-            for ifeat in range(niter):
-                for ii in range(offset[ifeat],offset[ifeat+1]):
-                    output[ifeat%nfeat] += cell_data[ii]
-            return output
-        self._areas   = __sum_up_cell_data(cell_areas,  self._offset,False)
-        self._volumes = __sum_up_cell_data(cell_volumes,self._offset,True )
+        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
@@ -841,19 +800,23 @@ class Features3d:
     def nfeatures(self): return self._nfeatures
 
     @property
-    def areas(self): return self._areas
+    def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
 
     @property
-    def volumes(self): return self._volumes
+    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 
+           will have a negative volume due to the direction of their normal vector.'''
+        self.discard_features(np.flatnonzero(self.volumes()<0),report=report)
+        return
 
     def reduce_noise(self,threshold=1e-5,is_relative=True,report=False):
         '''Discards all objects with smaller volume than a threshold.'''
         if is_relative:
             vol_domain = np.prod(self.spacing*self.dimensions)
             threshold = threshold*vol_domain
-        features = np.flatnonzero(np.abs(self._volumes)<threshold)
-        if len(features)>0:
-            self.discard_features(features,report=report)
+        self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
         return
 
     def discard_features(self,features,clean_points=False,report=False):
@@ -872,9 +835,9 @@ class Features3d:
         # Save former number of faces for report
         nfaces = self._faces.shape[0]
         # Delete indexed elements from arrays
-        self._faces   = np.delete(self._faces,idx,axis=0)
-        self._areas   = np.delete(self._areas,features,axis=0)
-        self._volumes = np.delete(self._volumes,features,axis=0)
+        self._faces        = np.delete(self._faces,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
         self._offset[1:] = self._offset[1:]-np.cumsum(gapsize)
         self._offset     = np.delete(self._offset,features,axis=0)
@@ -902,43 +865,50 @@ class Features3d:
         '''Returns the surface area of feature. If feature is None, total surface 
            area of all features is returned.'''
         if feature is None: 
-            return np.sum(self.areas)
+            return np.sum(self._cell_areas)
         else: 
-            return self.areas[feature]
+            return np.sum(self._cell_areas[self._offset[feature:feature+2]])
+
+    def areas(self): 
+        '''Returns an array with surface areas of all features.'''
+        return np.add.reduceat(self._cell_areas,self._offset[:-1])
 
     def volume(self,feature):
         '''Returns volume enclosed by feature. If feature isNone, total volume of 
            all features is returned.'''
         if feature is None: 
-            return np.sum(self.volumes)
+            return np.sum(self._cell_volumes)
         else: 
-            return self.volumes[feature]
+            return np.sum(self._cell_volumes[self._offset[feature:feature+2]])
+
+    def volumes(self):
+        '''Returns an array with volumes of all features.'''
+        return np.add.reduceat(self._cell_volumes,self._offset[:-1])
 
     def build_kdtree(self,kdaxis=0,leafsize=100,compact_nodes=False,balanced_tree=False,report=False):
         '''Builds a KD-tree for feature search.'''
-        from scipy.spatial import KDTree
-        from numba import jit
-        @jit(nopython=True,cache=True)
-        def __get_center_and_search_radius(points,faces,kdaxis):
-            nfaces = faces.shape[0]
-            center = np.zeros((nfaces,2),dtype=points.dtype)
-            radius = 0.0
-            if   kdaxis==0: i1,i2 = 1,2
-            elif kdaxis==1: i1,i2 = 0,2
-            elif kdaxis==2: i1,i2 = 0,1
-            for iface in range(nfaces):
-                min1 = min(min(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1])
-                max1 = max(max(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1])
-                min2 = min(min(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2])
-                max2 = max(max(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2])
-                center[iface,0] = 0.5*(max1+min1)
-                center[iface,1] = 0.5*(max2+min2)
-                radius = max(radius,(max1-min1)**2+(max2-min2)**2)
-            radius = 0.5*np.sqrt(radius)
-            return center,radius
-        center,radius = __get_center_and_search_radius(self._points,self._faces,kdaxis)
-        self._kdtree = KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
-                                     copy_data=False,balanced_tree=balanced_tree)
+        from scipy import spatial
+        if kdaxis==0:
+            min1 = np.amin(self._points[self._faces[:,1:],1],axis=1)
+            max1 = np.amax(self._points[self._faces[:,1:],1],axis=1)
+            min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
+            max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
+        elif kdaxis==1:
+            min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
+            max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
+            min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
+            max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
+        elif kdaxis==2:
+            min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
+            max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
+            min2 = np.amin(self._points[self._faces[:,1:],1],axis=1)
+            max2 = np.amax(self._points[self._faces[:,1:],1],axis=1)
+        else:
+            raise ValueError("Invalid kdaxis.")
+        center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1)
+        radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2))
+        self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
+                                             copy_data=False,balanced_tree=balanced_tree)
         self._kdaxis   = kdaxis
         self._kdradius = radius
         if report:
@@ -964,9 +934,6 @@ class Features3d:
                 t:         parameter to determine intersection point (x = R+t*dR) [float]
                 hit_dir:   direction from which the triangle was hit, from inward/outward = +1,-1 [int]
         '''
-        from numba import jit
-        from time import time
-
         coords = np.array(coords)
         assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array."
 
@@ -976,68 +943,44 @@ class Features3d:
         elif self._kdaxis==1: query_axis = [0,2]
         elif self._kdaxis==2: query_axis = [0,1]
 
-        t__ = time()
-        candidates = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius).tolist()
-        print('query',time()-t__)
-
-        t__ = time()
-        cand_num = np.asarray(tuple(map(len,candidates)))
-        cand_arr = np.concatenate(candidates).astype(np.int)
-        print('remap',time()-t__)
+        cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
 
+        Ncoord = coords.shape[0]
         raydir = np.zeros((3,),dtype=self._points.dtype)
         raydir[self._kdaxis] = 1.0
-        #
-        @jit(nopython=True,error_model='numpy',cache=True)
-        def __ray_triangle_intersect(R,dR,A,B,C):
-            '''Implements the Möller–Trumbore ray-triangle intersection algorithm. Taken from
-               https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d'''
-            E1     = B-A
-            E2     = C-A
-            N      = np.cross(E1,E2)
-            norm   = np.sqrt(np.sum(N*N))
-            det    = -np.sum(dR*N)
-            invdet = 1.0/det
-            AO     = R-A
-            DAO    = np.cross(AO,dR)
-            u      =  np.sum(E2*DAO)*invdet
-            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)
-        def __get_nearest_face(coords,raydir,cand_arr,cand_num,faces,points):
-            N = coords.shape[0]
-            nearface = np.zeros(N,dtype=np.int64)
-            isinside = np.zeros(N,dtype=np.bool_)
-            jj = 0
-            for ii in range(N):
-                t = np.inf
-                f = -1
-                n = 0
-                for kk in range(cand_num[ii]):
-                    idx = cand_arr[jj]
-                    hit_,t_ = __ray_triangle_intersect(coords[ii],raydir,
-                                                        points[faces[idx,1],:],
-                                                        points[faces[idx,2],:],
-                                                        points[faces[idx,3],:])
-                    if hit_: 
-                        n += 1
-                        if t_<t:
-                            t = t_
-                            f = idx
-                    jj += 1
-                nearface[ii] = f
-                isinside[ii] = (n%2)!=0
-            return nearface,isinside
-        t__ = time()
-        nf,isinside = __get_nearest_face(coords,raydir,cand_arr,cand_num,self._faces,self._points)
-        print('time:',time()-t__)
-        output = self.feature_from_face(nf)
-        output[np.logical_not(isinside)] = -1
+        in_feat  = np.empty((Ncoord,),dtype=np.int)
+        is_hit   = np.empty((Ncoord,),dtype=np.int)
+        hit_dir_ = np.empty((Ncoord,),dtype=np.int)
+        face_ = np.empty((Ncoord,),dtype=np.int)
+        print('point 12156 = ',coords[12156])
+
+        for ii in range(Ncoord):
+            hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir,
+                                                    self._points[self._faces[cand[ii],1],:],
+                                                    self._points[self._faces[cand[ii],2],:],
+                                                    self._points[self._faces[cand[ii],3],:])
+            # if ii==12279: print('DEBUG',hit_idx,t,N)
+            if ii==12156: print('DEBUG',hit_idx,t,hit_dir)
+            if hit_idx is None:
+                in_feat[ii] = -1
+            else:
+                idx = np.argmin(np.abs(t))
+                if ii==12156:
+                    for kk in hit_idx:
+                        print(cand[ii][kk])
+                    # print(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:])
+                    # print(cand[ii][hit_idx[idx]])
+                    # ifeat = self.feature_from_face(cand[ii][hit_idx[idx]])
+                    # print(self._offset[ifeat],self._offset[ifeat+1])
+                    # stop
+                # if hit_dir[idx]>0:
+                #     in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]])
+                # else:
+                #     in_feat[ii] = -1
         if report:
             print('[Features3d.inside_feature]',end=' ')
-            print('{} of {} points are located inside of features.'.format(sum(isinside),coords.shape[0]))
-        return output
+            print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord))
+        return in_feat
 
     def feature_from_face(self,idx_cell):
         '''Gets feature ID for a given cell.'''
@@ -1047,48 +990,39 @@ class Features3d:
         # 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.searchsorted(self._offset,idx_cell,side='right')-1)%self._nfeatures
+        return np.searchsorted(self._offset,idx_cell,side='right')-1
 
-    def faces_from_feature(self,features,incl_boundary=False):
+    def faces_from_feature(self,features,concat=True):
         '''Returns indices of cells which belong to given features.'''
-        features = self.list_of_features(features)
-        faces = []
-        for feature in features:
-            idx = np.arange(self._offset[feature],self._offset[feature+1])
-            faces.append(self._faces[idx,:])
-        if incl_boundary:
+        from collections.abc import Iterable
+        if features is None:
+            idx = slice(None)
+        elif not isinstance(features,Iterable):
+            idx = slice(self._offset[features],self._offset[features+1])
+        elif concat:        
+            idx = []
             for feature in features:
-                idx = np.arange(self._offset[self._nfeatures+feature],
-                                self._offset[self._nfeatures+feature+1])
-                faces.append(self._faces[idx,:])
-        return np.concatenate(faces,axis=0)
-
-    def regionid_from_feature(self,features,incl_boundary=False):
-        '''Returns regionid cell array to be added to PolyData.'''
-        features = self.list_of_features(features)
-        regionid = []
-        for feature in features:
-            nfaces = self._offset[feature+1]-self._offset[feature]
-            regionid.append(np.full(nfaces,feature,dtype=np.int64))
-        if incl_boundary:
+                rng_ = np.arange(self._offset[feature],self._offset[feature+1])
+                idx.append(rng_)
+            if len(idx)==1:
+                idx = np.array(idx)
+            else:
+                idx = np.concatenate(idx,axis=0)
+        else:
+            idx = []
             for feature in features:
-                nfaces = (self._offset[self._nfeatures+feature+1]-
-                          self._offset[self._nfeatures+feature])
-                regionid.append(np.full(nfaces,feature,dtype=np.int64))
-        return np.concatenate(regionid,axis=0)
+                rng_ = slice(self._offset[feature],self._offset[feature+1])
+                idx.append(rng_)
+        return idx
 
-    def nfaces_of_feature(self,features,incl_boundary=False):
+    def nfaces_of_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)
-        nfaces = []
+        ncells = []
         for feature in features:
-            nfaces.append(self._offset[feature+1]-self._offset[feature])
-        if incl_boundary:
-            for feature in features:
-                nfaces.append(self._offset[self._nfeatures+feature+1]-
-                              self._offset[self._nfeatures+feature])
-        return np.concatenate(nfaces,axis=0)
+            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.'''
@@ -1100,6 +1034,84 @@ class Features3d:
         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 
+           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.
+        '''
+        E1     = B-A
+        E2     = C-A
+        N      = np.cross(E1,E2)
+        det    = -(dR*N).sum(axis=-1) # dot product
+        det[np.abs(det)<1e-6] = np.nan # mask to avoid numpy runtime errors
+        invdet = 1.0/det
+        AO     = R-A
+        DAO    = np.cross(AO,dR)
+        # u,v,1-u-v are the barycentric coordinates of intersection
+        u      =  (E2*DAO).sum(axis=-1)*invdet
+        v      = -(E1*DAO).sum(axis=-1)*invdet
+        # Boolean array indicating hits
+        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)
+        # Intersection point is R+t*dR
+        t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
+        # return hit_idx,t,np.sign(t)*np.sign(det[hit_idx])
+        return hit_idx,t,N[hit_idx]
+
     def cmap_features(self,nfeatures,name='tab10'):
         from matplotlib.colors import ListedColormap
         if name=='tab10': # seaborn/matplotlib tab10 
@@ -1123,13 +1135,10 @@ 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):
         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)
-        return output
+        idx = self.faces_from_feature(features)
+        return pv.PolyData(self._points,self._faces[idx,:])
 
     def to_vtk2(self,features):
         import pyvista as pv