it works now

field.py

@@ -650,16 +650,15 @@ 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',cellvol_normal_component=2,
+                        contour_method='flying_edges',
                         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.dimensions  = np.array(input.shape,dtype=np.int)
+        self.origin      = np.array(origin,dtype=np.float)-self.spacing
         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.
@@ -668,15 +667,20 @@ 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.
+        # 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.
         if has_ghost:
-            self._input = sign_invert*input
+            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
         else:
             pw = tuple((0,1) if x else (0,0) for x in periodicity)
-            self._input = np.pad(sign_invert*input,pw,mode='wrap')
+            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')
         # Triangulate
         self.triangulate(contour_method=contour_method,
-                         cellvol_normal_component=cellvol_normal_component,
                          report=report)
         # Set some state variable
         self._kdtree   = None
@@ -689,108 +693,145 @@ class Features3d:
     @classmethod
     def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False,
                             keep_input=False,contour_method='flying_edges',
-                            cellvol_normal_component=2,report=False):
+                            correct_normals=True,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',cellvol_normal_component=2,report=False):
+    def triangulate(self,contour_method='flying_edges',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._input.shape
+        datavtk.dimensions = self.dimensions
         datavtk.origin     = self.origin
         datavtk.spacing    = self.spacing
         datavtk.point_arrays['data'] = self._input.ravel('F')
-        # Compute full contour: ensure we only have triangles to efficiently extract points
+        # Compute full contour: due to the padding, the contour will result in closed object.
-        # later on.
+        # 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.
         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."
-        # Compute the connectivity of the triangulated surface: first we run an ordinary 
+        points    = contour.points
-        # connectivity filter neglecting periodic wrapping.
+        faces     = contour.faces.reshape(-1,4)
-        if report: print('[Features3d.triangulate] computing connectivity...')
+        gradients = contour.point_arrays['Gradients'][faces[:,1],:]
-        alg = vtk.vtkPolyDataConnectivityFilter() # ~twice as fast as default vtkConnectivityFilter
+        # 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()
         alg.SetInputData(contour)
         alg.SetExtractionModeToAllRegions()
         alg.SetColorRegions(True)
         alg.Update()
         contour = pv.filters._get_output(alg)
-        # Now determine the boundary points, i.e. points on the boundary which have a corresponding
+        # RegionIds are now stored as point data. To efficiently convert it to cell data, the first 
-        # 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...')
-        cell_labels  = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
+        self._points = np.array(points)
-        offset       = np.zeros((self._nfeatures+1,),dtype=np.int64)
+        self._nfeatures = nfeatures
-        offset[1:]   = np.cumsum(np.bincount(cell_labels))
+        idx   = np.flatnonzero(bd_face_tag<0)
-        ind          = np.argsort(cell_labels)
+        idxbd = np.flatnonzero(bd_face_tag>=0)
-        self._offset = offset
+        self._offset                = np.zeros((2*nfeatures+1,),dtype=np.int64)
-        self._faces  = contour.faces.reshape(contour.n_faces,4)[ind,:]
+        self._offset[1:nfeatures+1] = np.cumsum(np.bincount(regionid[idx],minlength=nfeatures))
-        self._points = contour.points
+        self._offset[nfeatures+1:]  = np.cumsum(np.bincount(regionid[idxbd],minlength=nfeatures))
-        # Compute the volume and area per cell. For the volume computation, an arbitrary component 
+        self._offset[nfeatures+1:] += self._offset[nfeatures]
-        # of the normal has to be chosen which defaults to the z-component and is set by 
+        self._faces = np.concatenate((
-        # 'cellvol_normal_component'.
+                        faces[idx,:][np.argsort(regionid[idx]),:],
-        if report: print('[Features3d.triangulate] calculating area and volume per cell...')
+                        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...')
         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 = (contour.point_arrays['Gradients'][self._faces[:,1],:]*cn).sum(axis=-1)>0
+        idx = (gradients*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])]
-        # self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])]
+        cn[idx,:] = -cn[idx,:]
-        # cn[idx,:] = -cn[idx,:]
         # Compute area and signed volume per cell
         cc = (A+B+C)/3
-        self._cell_areas   = 0.5*np.sqrt(np.square(cn).sum(axis=1))
+        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]
+        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 )
         return
 
     @property
@@ -800,23 +841,19 @@ class Features3d:
     def nfeatures(self): return self._nfeatures
 
     @property
-    def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
+    def areas(self): return self._areas
 
     @property
-    def cell_volumes(self): return np.split(self._cell_volumes,self._offset[1:-1])
+    def volumes(self): return self._volumes
-
-    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
-        self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
+        features = np.flatnonzero(np.abs(self._volumes)<threshold)
+        if len(features)>0:
+            self.discard_features(features,report=report)
         return
 
     def discard_features(self,features,clean_points=False,report=False):
@@ -836,8 +873,8 @@ class Features3d:
         nfaces = self._faces.shape[0]
         # Delete indexed elements from arrays
         self._faces   = np.delete(self._faces,idx,axis=0)
-        self._cell_areas   = np.delete(self._cell_areas,idx,axis=0)
+        self._areas   = np.delete(self._areas,features,axis=0)
-        self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
+        self._volumes = np.delete(self._volumes,features,axis=0)
         # Correct offset
         self._offset[1:] = self._offset[1:]-np.cumsum(gapsize)
         self._offset     = np.delete(self._offset,features,axis=0)
@@ -865,49 +902,42 @@ 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._cell_areas)
+            return np.sum(self.areas)
         else: 
-            return np.sum(self._cell_areas[self._offset[feature:feature+2]])
+            return self.areas[feature]
-
-    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._cell_volumes)
+            return np.sum(self.volumes)
         else: 
-            return np.sum(self._cell_volumes[self._offset[feature:feature+2]])
+            return self.volumes[feature]
-
-    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 import spatial
+        from scipy.spatial import KDTree
-        if kdaxis==0:
+        from numba import jit
-            min1 = np.amin(self._points[self._faces[:,1:],1],axis=1)
+        @jit(nopython=True,cache=True)
-            max1 = np.amax(self._points[self._faces[:,1:],1],axis=1)
+        def __get_center_and_search_radius(points,faces,kdaxis):
-            min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
+            nfaces = faces.shape[0]
-            max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
+            center = np.zeros((nfaces,2),dtype=points.dtype)
-        elif kdaxis==1:
+            radius = 0.0
-            min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
+            if   kdaxis==0: i1,i2 = 1,2
-            max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
+            elif kdaxis==1: i1,i2 = 0,2
-            min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
+            elif kdaxis==2: i1,i2 = 0,1
-            max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
+            for iface in range(nfaces):
-        elif kdaxis==2:
+                min1 = min(min(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1])
-            min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
+                max1 = max(max(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1])
-            max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
+                min2 = min(min(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2])
-            min2 = np.amin(self._points[self._faces[:,1:],1],axis=1)
+                max2 = max(max(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2])
-            max2 = np.amax(self._points[self._faces[:,1:],1],axis=1)
+                center[iface,0] = 0.5*(max1+min1)
-        else:
+                center[iface,1] = 0.5*(max2+min2)
-            raise ValueError("Invalid kdaxis.")
+                radius = max(radius,(max1-min1)**2+(max2-min2)**2)
-        center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1)
+            radius = 0.5*np.sqrt(radius)
-        radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2))
+            return center,radius
-        self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
+        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)
         self._kdaxis   = kdaxis
         self._kdradius = radius
@@ -934,6 +964,9 @@ 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."
 
@@ -943,44 +976,68 @@ class Features3d:
         elif self._kdaxis==1: query_axis = [0,2]
         elif self._kdaxis==2: query_axis = [0,1]
 
-        cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
+        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__)
 
-        Ncoord = coords.shape[0]
         raydir = np.zeros((3,),dtype=self._points.dtype)
         raydir[self._kdaxis] = 1.0
-        in_feat  = np.empty((Ncoord,),dtype=np.int)
+        #
-        is_hit   = np.empty((Ncoord,),dtype=np.int)
+        @jit(nopython=True,error_model='numpy',cache=True)
-        hit_dir_ = np.empty((Ncoord,),dtype=np.int)
+        def __ray_triangle_intersect(R,dR,A,B,C):
-        face_ = np.empty((Ncoord,),dtype=np.int)
+            '''Implements the Möller–Trumbore ray-triangle intersection algorithm. Taken from
-        print('point 12156 = ',coords[12156])
+               https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d'''
-
+            E1     = B-A
-        for ii in range(Ncoord):
+            E2     = C-A
-            hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir,
+            N      = np.cross(E1,E2)
-                                                    self._points[self._faces[cand[ii],1],:],
+            norm   = np.sqrt(np.sum(N*N))
-                                                    self._points[self._faces[cand[ii],2],:],
+            det    = -np.sum(dR*N)
-                                                    self._points[self._faces[cand[ii],3],:])
+            invdet = 1.0/det
-            # if ii==12279: print('DEBUG',hit_idx,t,N)
+            AO     = R-A
-            if ii==12156: print('DEBUG',hit_idx,t,hit_dir)
+            DAO    = np.cross(AO,dR)
-            if hit_idx is None:
+            u      =  np.sum(E2*DAO)*invdet
-                in_feat[ii] = -1
+            v      = -np.sum(E1*DAO)*invdet
-            else:
+            t      =  np.sum(AO*N)  *invdet
-                idx = np.argmin(np.abs(t))
+            return (np.abs(det)>=1e-7*norm and t>=0.0 and u>=0 and v>=0 and (u+v)-1.0<=0), t
-                if ii==12156:
+        # @jit(nopython=True,cache=True)
-                    for kk in hit_idx:
+        def __get_nearest_face(coords,raydir,cand_arr,cand_num,faces,points):
-                        print(cand[ii][kk])
+            N = coords.shape[0]
-                    # print(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:])
+            nearface = np.zeros(N,dtype=np.int64)
-                    # print(cand[ii][hit_idx[idx]])
+            isinside = np.zeros(N,dtype=np.bool_)
-                    # ifeat = self.feature_from_face(cand[ii][hit_idx[idx]])
+            jj = 0
-                    # print(self._offset[ifeat],self._offset[ifeat+1])
+            for ii in range(N):
-                    # stop
+                t = np.inf
-                # if hit_dir[idx]>0:
+                f = -1
-                #     in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]])
+                n = 0
-                # else:
+                for kk in range(cand_num[ii]):
-                #     in_feat[ii] = -1
+                    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
         if report:
             print('[Features3d.inside_feature]',end=' ')
-            print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord))
+            print('{} of {} points are located inside of features.'.format(sum(isinside),coords.shape[0]))
-        return in_feat
+        return output
 
     def feature_from_face(self,idx_cell):
         '''Gets feature ID for a given cell.'''
@@ -990,39 +1047,48 @@ 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
+        return (np.searchsorted(self._offset,idx_cell,side='right')-1)%self._nfeatures
 
-    def faces_from_feature(self,features,concat=True):
+    def faces_from_feature(self,features,incl_boundary=False):
         '''Returns indices of cells which belong to given features.'''
-        from collections.abc import Iterable
+        features = self.list_of_features(features)
-        if features is None:
+        faces = []
-            idx = slice(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 = np.arange(self._offset[feature],self._offset[feature+1])
-                idx.append(rng_)
+            faces.append(self._faces[idx,:])
-            if len(idx)==1:
+        if incl_boundary:
-                idx = np.array(idx)
-            else:
-                idx = np.concatenate(idx,axis=0)
-        else:
-            idx = []
             for feature in features:
-                rng_ = slice(self._offset[feature],self._offset[feature+1])
+                idx = np.arange(self._offset[self._nfeatures+feature],
-                idx.append(rng_)
+                                self._offset[self._nfeatures+feature+1])
-        return idx
+                faces.append(self._faces[idx,:])
+        return np.concatenate(faces,axis=0)
 
-    def nfaces_of_feature(self,features):
+    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:
+            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)
+
+    def nfaces_of_feature(self,features,incl_boundary=False):
         '''Returns number of cells which belong to given features. Can be used
            to construct new offsets.'''
         features = self.list_of_features(features)
-        ncells = []
+        nfaces = []
         for feature in features:
-            ncells.append(self._offset[feature+1]-self._offset[feature])
+            nfaces.append(self._offset[feature+1]-self._offset[feature])
-        return np.array(ncells)
+        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)
 
     def list_of_features(self,features):
         '''Ensures that 'features' is a list.'''
@@ -1034,84 +1100,6 @@ 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 
@@ -1135,10 +1123,13 @@ class Features3d:
             icolor = (icolor+1)%ncolor
         return ListedColormap(output)
 
-    def to_vtk(self,features):
+    def to_vtk(self,features,incl_regionid=False,incl_boundary=False):
         import pyvista as pv
-        idx = self.faces_from_feature(features)
+        faces  = self.faces_from_feature(features,incl_boundary=incl_boundary)
-        return pv.PolyData(self._points,self._faces[idx,:])
+        output = pv.PolyData(self._points,faces)
+        if incl_regionid:
+            output.cell_arrays['RegionId'] = self.regionid_from_feature(features,incl_boundary=incl_boundary)
+        return output
 
     def to_vtk2(self,features):
         import pyvista as pv