mwtkrayer
/
suspendtools
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

it works now

This commit is contained in:
Michael Krayer 2021-08-19 09:29:43 +02:00
parent e95b7dc74a
commit 928cc9c9d2
1 changed files with 166 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

394
field.py
View File

@ -650,7 +650,7 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0):
class Features3d: class Features3d:
def __init__(self,input,threshold,origin,spacing,periodicity, def __init__(self,input,threshold,origin,spacing,periodicity,
invert=False,has_ghost=False,keep_input=False, invert=False,has_ghost=False,keep_input=False,
contour_method='flying_edges',cellvol_normal_component=2, contour_method='flying_edges',
report=False): report=False):
assert len(origin)==3, "'origin' must be of length 3" assert len(origin)==3, "'origin' must be of length 3"
assert len(spacing)==3, "'spacing' must be of length 3" assert len(spacing)==3, "'spacing' must be of length 3"
@ -681,7 +681,6 @@ class Features3d:
self._input[1:-1,1:-1,1:-1] = np.pad(sign_invert*input,pw,mode='wrap') self._input[1:-1,1:-1,1:-1] = np.pad(sign_invert*input,pw,mode='wrap')
# Triangulate # Triangulate
self.triangulate(contour_method=contour_method, self.triangulate(contour_method=contour_method,
cellvol_normal_component=cellvol_normal_component,
report=report) report=report)
# Set some state variable # Set some state variable
self._kdtree = None self._kdtree = None
@ -694,14 +693,13 @@ class Features3d:
@classmethod @classmethod
def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False, def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False,
keep_input=False,contour_method='flying_edges', 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, return cls(fld.data,threshold,fld.origin,fld.spacing,periodicity,
invert=invert,has_ghost=has_ghost,keep_input=keep_input, invert=invert,has_ghost=has_ghost,keep_input=keep_input,
contour_method=contour_method, contour_method=contour_method,
cellvol_normal_component=cellvol_normal_component,
report=report) 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 pyvista as pv
import vtk import vtk
from scipy import ndimage, spatial from scipy import ndimage, spatial
@ -723,8 +721,9 @@ class Features3d:
if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method)) 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) contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False,compute_gradients=True)
assert contour.is_all_triangles(), "Contouring produced non-triangle cells." assert contour.is_all_triangles(), "Contouring produced non-triangle cells."
points = contour.points points = contour.points
faces = contour.faces.reshape(-1,4) 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 # 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 # some of the necessary check. Therefore let's define some functions to be jit compiled
# by numba to speed things up. # by numba to speed things up.
@ -734,9 +733,9 @@ class Features3d:
assert points.shape[1]==3 assert points.shape[1]==3
assert len(bounds)==6 assert len(bounds)==6
assert len(tolerance)==3 assert len(tolerance)==3
ncells = faces.shape[0] nfaces = faces.shape[0]
output = np.full((ncells,),-1,dtype=np.int8) output = np.full((nfaces,),-1,dtype=np.int8)
for ii in range(ncells): for ii in range(nfaces):
for jj in range(3): for jj in range(3):
if (points[faces[ii,1],jj]<bounds[2*jj]+tolerance[jj] and 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,2],jj]<bounds[2*jj]+tolerance[jj] and
@ -757,10 +756,12 @@ class Features3d:
output[jj,:] = (3,face_uni_hi[ii],face_uni_lo[idx[ii]],face_uni_hi[ii]) output[jj,:] = (3,face_uni_hi[ii],face_uni_lo[idx[ii]],face_uni_hi[ii])
jj += 1 jj += 1
return output[:jj,:] 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 tol_overlap = 1e-5*self.spacing
bd_face_tag = __get_boundary_faces(faces,points,tuple(contour.bounds),tol_overlap) bd_face_tag = __get_boundary_faces(faces,points,tuple(contour.bounds),tol_overlap)
#####
faces_connect = [faces] faces_connect = [faces]
for axis in range(3): for axis in range(3):
if not self.periodicity[axis]: continue if not self.periodicity[axis]: continue
@ -774,10 +775,11 @@ class Features3d:
distance_upper_bound=tol_overlap[axis]) distance_upper_bound=tol_overlap[axis])
faces_connect += [__connect_faces_periodic(face_uni_lo,face_uni_hi,points, faces_connect += [__connect_faces_periodic(face_uni_lo,face_uni_hi,points,
dist,idx,tol_overlap[axis])] dist,idx,tol_overlap[axis])]
##### # Save original number of faces and add the virtual connecting faces to the polydata.
ncells = faces.shape[0] nfaces = faces.shape[0]
contour.faces = np.concatenate(faces_connect,axis=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 = vtk.vtkPolyDataConnectivityFilter()
alg.SetInputData(contour) alg.SetInputData(contour)
alg.SetExtractionModeToAllRegions() alg.SetExtractionModeToAllRegions()
@ -786,43 +788,50 @@ class Features3d:
contour = pv.filters._get_output(alg) contour = pv.filters._get_output(alg)
# RegionIds are now stored as point data. To efficiently convert it to cell data, the first # RegionIds 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. # point of each cell determines the value for this cell.
regionid = contour.point_arrays['RegionId'][contour.faces.reshape(-1,4)[:ncells,1]].ravel() regionid = contour.point_arrays['RegionId'][contour.faces.reshape(-1,4)[:nfaces,1]].ravel()
nregion = np.amax(regionid)+1 nfeatures = np.amax(regionid)+1
# Store the face/vertex data in class arrays. Sort the faces by RegionId for quick access # Store the face/vertex data in class arrays. Sort the faces by RegionId for quick access
# later on. # later on.
if report: print('[Features3d.triangulate] sorting faces by labels...') if report: print('[Features3d.triangulate] sorting faces by labels...')
self._points = points self._points = np.array(points)
self._nfeatures = nregion self._nfeatures = nfeatures
#
idx = np.flatnonzero(bd_face_tag<0) idx = np.flatnonzero(bd_face_tag<0)
idxbd = np.flatnonzero(bd_face_tag>=0) idxbd = np.flatnonzero(bd_face_tag>=0)
self._offset = np.zeros((nregion+1,),dtype=np.int64) self._offset = np.zeros((2*nfeatures+1,),dtype=np.int64)
self._offset[1:] = np.cumsum(np.bincount(regionid[idx],minlength=nregion)) self._offset[1:nfeatures+1] = np.cumsum(np.bincount(regionid[idx],minlength=nfeatures))
self._offsetbd = np.full((nregion+1,),self._offset[-1],dtype=np.int64) self._offset[nfeatures+1:] = np.cumsum(np.bincount(regionid[idxbd],minlength=nfeatures))
self._offsetbd[1:] += np.cumsum(np.bincount(regionid[idxbd],minlength=nregion)) self._offset[nfeatures+1:] += self._offset[nfeatures]
self._faces = np.concatenate(( self._faces = np.concatenate((
faces[idx,:][np.argsort(regionid[idx]),:], faces[idx,:][np.argsort(regionid[idx]),:],
faces[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0) faces[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0)
# # gradients = np.concatenate((
# # self._gradient = (contour.point_arrays['Gradients'][self._faces[:,1],:] gradients[idx,:][np.argsort(regionid[idx]),:],
gradients[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0)
# # Compute the volume and area per cell. For the volume computation, an arbitrary component # Compute the volume and area per cell.
# # of the normal has to be chosen which defaults to the z-component and is set by if report: print('[Features3d.triangulate] calculating surface area and volume...')
# # 'cellvol_normal_component'. A = self._points[self._faces[:,1],:]
# if report: print('[Features3d.triangulate] calculating area and volume per cell...') B = self._points[self._faces[:,2],:]
# A = self._points[self._faces[:,1],:] C = self._points[self._faces[:,3],:]
# B = self._points[self._faces[:,2],:] cn = np.cross(B-A,C-A)
# C = self._points[self._faces[:,3],:] # Check if cell normal points in direction of gradient. If not, switch vertex order.
# cn = np.cross(B-A,C-A) idx = (gradients*cn).sum(axis=-1)>0
# # Check if cell normal points in direction of gradient. If not, switch vertex order. self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])]
# idx = (contour.point_arrays['Gradients'][self._faces[:,1],:]*cn).sum(axis=-1)>0 cn[idx,:] = -cn[idx,:]
# # print(idx.shape,np.sum(idx),self._faces.shape,self._faces[idx,2:].shape,self._faces[idx,3:1:-1].shape) # Compute area and signed volume per cell
# # self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])] cc = (A+B+C)/3
# # cn[idx,:] = -cn[idx,:] cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1))
# # Compute area and signed volume per cell cell_volumes = 0.5*np.mean(cn*cc,axis=1)
# cc = (A+B+C)/3 @jit(nopython=True)
# self._cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1)) def __sum_up_cell_data(cell_data,offset,incl_boundary):
# self._cell_volumes = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component] 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 return
@property @property
@ -832,23 +841,19 @@ class Features3d:
def nfeatures(self): return self._nfeatures def nfeatures(self): return self._nfeatures
@property @property
def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1]) def areas(self): return self._areas
@property @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): def reduce_noise(self,threshold=1e-5,is_relative=True,report=False):
'''Discards all objects with smaller volume than a threshold.''' '''Discards all objects with smaller volume than a threshold.'''
if is_relative: if is_relative:
vol_domain = np.prod(self.spacing*self.dimensions) vol_domain = np.prod(self.spacing*self.dimensions)
threshold = threshold*vol_domain 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 return
def discard_features(self,features,clean_points=False,report=False): def discard_features(self,features,clean_points=False,report=False):
@ -867,9 +872,9 @@ class Features3d:
# Save former number of faces for report # Save former number of faces for report
nfaces = self._faces.shape[0] nfaces = self._faces.shape[0]
# Delete indexed elements from arrays # Delete indexed elements from arrays
self._faces = np.delete(self._faces,idx,axis=0) 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 # Correct offset
self._offset[1:] = self._offset[1:]-np.cumsum(gapsize) self._offset[1:] = self._offset[1:]-np.cumsum(gapsize)
self._offset = np.delete(self._offset,features,axis=0) self._offset = np.delete(self._offset,features,axis=0)
@ -897,52 +902,43 @@ class Features3d:
'''Returns the surface area of feature. If feature is None, total surface '''Returns the surface area of feature. If feature is None, total surface
area of all features is returned.''' area of all features is returned.'''
if feature is None: if feature is None:
return np.sum(self._cell_areas) return np.sum(self.areas)
else: 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): def volume(self,feature):
'''Returns volume enclosed by feature. If feature isNone, total volume of '''Returns volume enclosed by feature. If feature isNone, total volume of
all features is returned.''' all features is returned.'''
if feature is None: if feature is None:
return np.sum(self._cell_volumes) return np.sum(self.volumes)
else: 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): def build_kdtree(self,kdaxis=0,leafsize=100,compact_nodes=False,balanced_tree=False,report=False):
'''Builds a KD-tree for feature search.''' '''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)
copy_data=False,balanced_tree=balanced_tree) self._kdtree = KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
copy_data=False,balanced_tree=balanced_tree)
self._kdaxis = kdaxis self._kdaxis = kdaxis
self._kdradius = radius self._kdradius = radius
if report: if report:
@ -968,6 +964,9 @@ class Features3d:
t: parameter to determine intersection point (x = R+t*dR) [float] 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] 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) coords = np.array(coords)
assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array." assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array."
@ -977,44 +976,68 @@ class Features3d:
elif self._kdaxis==1: query_axis = [0,2] elif self._kdaxis==1: query_axis = [0,2]
elif self._kdaxis==2: query_axis = [0,1] 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 = np.zeros((3,),dtype=self._points.dtype)
raydir[self._kdaxis] = 1.0 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: if report:
print('[Features3d.inside_feature]',end=' ') 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): def feature_from_face(self,idx_cell):
'''Gets feature ID for a given cell.''' '''Gets feature ID for a given cell.'''
@ -1024,30 +1047,20 @@ class Features3d:
# as soon as idx_cell is equal or larger than offset, the argmax function # 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 # short circuits and returns the index of the first occurence without
# checking any value afterwards. # 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,incl_boundary=False): def faces_from_feature(self,features,incl_boundary=False):
'''Returns indices of cells which belong to given features.''' '''Returns indices of cells which belong to given features.'''
from collections.abc import Iterable features = self.list_of_features(features)
faces = [] faces = []
if features is None: for feature in features:
faces.append(self._faces) idx = np.arange(self._offset[feature],self._offset[feature+1])
if incl_boundary:
faces.append(self._faces)
elif not isinstance(features,Iterable):
idx = slice(self._offset[features],self._offset[features+1])
faces.append(self._faces[idx,:]) faces.append(self._faces[idx,:])
if incl_boundary: if incl_boundary:
idx = slice(self._offsetbd[features],self._offsetbd[features+1])
faces.append(self._faces[idx,:])
else:
for feature in features: for feature in features:
idx = np.arange(self._offset[feature],self._offset[feature+1]) idx = np.arange(self._offset[self._nfeatures+feature],
self._offset[self._nfeatures+feature+1])
faces.append(self._faces[idx,:]) faces.append(self._faces[idx,:])
if incl_boundary:
for feature in features:
idx = np.arange(self._offsetbd[feature],self._offsetbd[feature+1])
faces.append(self._faces[idx,:])
return np.concatenate(faces,axis=0) return np.concatenate(faces,axis=0)
def regionid_from_feature(self,features,incl_boundary=False): def regionid_from_feature(self,features,incl_boundary=False):
@ -1059,18 +1072,23 @@ class Features3d:
regionid.append(np.full(nfaces,feature,dtype=np.int64)) regionid.append(np.full(nfaces,feature,dtype=np.int64))
if incl_boundary: if incl_boundary:
for feature in features: for feature in features:
nfaces = self._offsetbd[feature+1]-self._offsetbd[feature] nfaces = (self._offset[self._nfeatures+feature+1]-
self._offset[self._nfeatures+feature])
regionid.append(np.full(nfaces,feature,dtype=np.int64)) regionid.append(np.full(nfaces,feature,dtype=np.int64))
return np.concatenate(regionid,axis=0) return np.concatenate(regionid,axis=0)
def nfaces_of_feature(self,features): def nfaces_of_feature(self,features,incl_boundary=False):
'''Returns number of cells which belong to given features. Can be used '''Returns number of cells which belong to given features. Can be used
to construct new offsets.''' to construct new offsets.'''
features = self.list_of_features(features) features = self.list_of_features(features)
ncells = [] nfaces = []
for feature in features: 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): def list_of_features(self,features):
'''Ensures that 'features' is a list.''' '''Ensures that 'features' is a list.'''
@ -1082,84 +1100,6 @@ class Features3d:
else: else:
return list(features) 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'): def cmap_features(self,nfeatures,name='tab10'):
from matplotlib.colors import ListedColormap from matplotlib.colors import ListedColormap
if name=='tab10': # seaborn/matplotlib tab10 if name=='tab10': # seaborn/matplotlib tab10
@ -1185,10 +1125,8 @@ class Features3d:
def to_vtk(self,features,incl_regionid=False,incl_boundary=False): def to_vtk(self,features,incl_regionid=False,incl_boundary=False):
import pyvista as pv import pyvista as pv
faces = self.faces_from_feature(features,incl_boundary=incl_boundary) faces = self.faces_from_feature(features,incl_boundary=incl_boundary)
print(faces.shape)
output = pv.PolyData(self._points,faces) output = pv.PolyData(self._points,faces)
print(output.n_faces)
if incl_regionid: if incl_regionid:
output.cell_arrays['RegionId'] = self.regionid_from_feature(features,incl_boundary=incl_boundary) output.cell_arrays['RegionId'] = self.regionid_from_feature(features,incl_boundary=incl_boundary)
return output return output