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

update: inside/outside detection does not work properly

This commit is contained in:
Michael Krayer 2021-08-16 12:58:23 +02:00
parent d95453ad45
commit c5908e00f5
1 changed files with 105 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

248
field.py
View File

@ -714,12 +714,6 @@ 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."
# print(contour.point_arrays['Gradients'].shape)
# print('gradient',time()-t)
# print(np.sum(gn>0),'of',len(gn))
# Compute the connectivity of the triangulated surface: first we run an ordinary # Compute the connectivity of the triangulated surface: first we run an ordinary
# connectivity filter neglecting periodic wrapping. # connectivity filter neglecting periodic wrapping.
if report: print('[Features3d.triangulate] computing connectivity...') if report: print('[Features3d.triangulate] computing connectivity...')
@ -767,70 +761,38 @@ class Features3d:
while target_ != map_[target_]: # map it recursively while target_ != map_[target_]: # map it recursively
target_ = map_[target_] target_ = map_[target_]
map_[source_] = target_ map_[source_] = target_
# map_ = np.unique(map_,return_inverse=True)[1]
map_ = np.unique(map_,return_inverse=True)[1]
point_labels = map_[point_labels] point_labels = map_[point_labels]
self._nfeatures = np.max(map_)+1 # starts with zero self._nfeatures = np.max(map_)+1 # starts with zero
# Labels are now stored as point data. To efficiently convert it to cell data, the first # 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. # point of each cell determines the value for this cell.
if report: print('[Features3d.triangulate] identifying cell based labels...') if report: print('[Features3d.triangulate] sorting faces by labels...')
cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]] cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
cell_gradient = contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1]] offset = np.zeros((self._nfeatures+1,),dtype=np.int64)
print(cell_gradient.shape) offset[1:] = np.cumsum(np.bincount(cell_labels))
# We are done with VTK for now. Since we are only dealing with triangles, let us ind = np.argsort(cell_labels)
# store them in a structure which is quicker to access and already sorted based self._offset = offset
# on labels. self._faces = contour.faces.reshape(contour.n_faces,4)[ind,:]
# _vertices has dim(ncells,nvert,ndim), where nvert=3 (three vertices per triangle) self._points = contour.points
# and ndim=3 (three dimensional coordinates) # self._offset = offset
# _offset is the cumulative number of cells for each feature and the offset to # self._faces = contour.faces.reshape(contour.n_faces,4)
# access the _vertices array by feature. # self._points = contour.points
if report: print('[Features3d.triangulate] storing vertices...')
ind = np.argsort(cell_labels)
self._offset = np.zeros((self._nfeatures+1,),dtype=np.int64)
self._offset[1:] = np.cumsum(np.bincount(cell_labels))
self._vertices = contour.points[contour.faces.reshape(contour.n_faces,4)[ind,1:]]
# Compute the volume and area per cell. For the volume computation, an arbitrary component # 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 # of the normal has to be chosen which defaults to the z-component and is set by
# 'cellvol_normal_component'. # 'cellvol_normal_component'.
if report: print('[Features3d.triangulate] calculating area and volume per cell...') if report: print('[Features3d.triangulate] calculating area and volume per cell...')
A = self._vertices[:,0,:] A = self._points[self._faces[:,1],:]
B = self._vertices[:,1,:] B = self._points[self._faces[:,2],:]
C = self._vertices[:,2,:] C = self._points[self._faces[:,3],:]
cn = np.cross(B-A,C-A) 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
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 cc = (A+B+C)/3
# print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]].shape)
cg = (contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:]+
contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,1,:]+
contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,2,:])/3
ng = (cg*cn).sum(axis=-1)
print('interpolated:')
print(cg)
print('component 0:')
print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:])
print('projection')
print(sum(ng>0),len(ng))
# self.cg = cg
# self.cc = cc
self.cg = contour.point_arrays['Gradients']
self.cc = contour.points
# .sum(axis=1))/3)*cn).sum(axis=-1)
self._cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1)) self._cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1))
self._cell_volumes = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component] self._cell_volumes = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
# # Compute gradient of field at cell center to get correct cell normal direction
# t = time()
# cn = cn/np.sqrt(np.square(cn).sum())
# ccpix = ((cc-self.origin)/self.spacing).transpose()
# ccnpix = (((cc+1e-4*cn)-self.origin)/self.spacing).transpose()
# val_cc = ndimage.map_coordinates(self._input,ccpix,order=1,prefilter=False)
# val_ccn = ndimage.map_coordinates(self._input,ccnpix,order=1,prefilter=False)
# gn = (val_ccn-val_cc)
# print('gradient',time()-t)
# print(cg)
# print(np.sum(cg>0),'of',cg.shape)
# stop
return return
@property @property
@ -859,7 +821,7 @@ class Features3d:
self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report) self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
return return
def discard_features(self,features,report=False): def discard_features(self,features,clean_points=False,report=False):
'''Removes features from triangulated data.''' '''Removes features from triangulated data.'''
features = self.list_of_features(features) features = self.list_of_features(features)
# Get index ranges which are to be deleted and also create an array # Get index ranges which are to be deleted and also create an array
@ -872,9 +834,9 @@ class Features3d:
gapsize[feature] = len(rng_) gapsize[feature] = len(rng_)
idx = np.concatenate(idx,axis=0) idx = np.concatenate(idx,axis=0)
# Save former number of faces for report # Save former number of faces for report
ncells = self._vertices.shape[0] nfaces = self._faces.shape[0]
# Delete indexed elements from arrays # Delete indexed elements from arrays
self._vertices = np.delete(self._vertices,idx,axis=0) self._faces = np.delete(self._faces,idx,axis=0)
self._cell_areas = np.delete(self._cell_areas,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) self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
# Correct offset # Correct offset
@ -884,7 +846,20 @@ class Features3d:
if report: if report:
print('[Features3d.discard_features]',end=' ') print('[Features3d.discard_features]',end=' ')
print('discarded {:} features with {:} faces.'.format( print('discarded {:} features with {:} faces.'.format(
len(features),ncells-self._vertices.shape[0])) len(features),nfaces-self._faces.shape[0]))
if clean_points:
self.clean_points(report=report)
return
def clean_points(self,report=False):
'''Removes points which are not referenced by any face.'''
nfaces = self._faces.shape[0]
ind,inv = np.unique(self._faces[:,1:],return_inverse=True)
if report:
print('[Features3d.clean_points]',end=' ')
print('removed {:} orphan points.'.format(self._points.shape[0]-len(ind)))
self._faces[:,1:4] = inv.reshape(nfaces,3)
self._points = self._points[ind,:]
return return
def area(self,feature): def area(self,feature):
@ -915,26 +890,26 @@ class Features3d:
'''Builds a KD-tree for feature search.''' '''Builds a KD-tree for feature search.'''
from scipy import spatial from scipy import spatial
if kdaxis==0: if kdaxis==0:
min1 = np.amin(self._vertices[:,:,1],axis=1) min1 = np.amin(self._points[self._faces[:,1:],1],axis=1)
max1 = np.amax(self._vertices[:,:,1],axis=1) max1 = np.amax(self._points[self._faces[:,1:],1],axis=1)
min2 = np.amin(self._vertices[:,:,2],axis=1) min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
max2 = np.amax(self._vertices[:,:,2],axis=1) max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
elif kdaxis==1: elif kdaxis==1:
min1 = np.amin(self._vertices[:,:,0],axis=1) min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
max1 = np.amax(self._vertices[:,:,0],axis=1) max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
min2 = np.amin(self._vertices[:,:,2],axis=1) min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
max2 = np.amax(self._vertices[:,:,2],axis=1) max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
elif kdaxis==2: elif kdaxis==2:
min1 = np.amin(self._vertices[:,:,0],axis=1) min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
max1 = np.amax(self._vertices[:,:,0],axis=1) max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
min2 = np.amin(self._vertices[:,:,1],axis=1) min2 = np.amin(self._points[self._faces[:,1:],1],axis=1)
max2 = np.amax(self._vertices[:,:,1],axis=1) max2 = np.amax(self._points[self._faces[:,1:],1],axis=1)
else: else:
raise ValueError("Invalid ray axis.") raise ValueError("Invalid kdaxis.")
center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1) 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)) radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2))
self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes, self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
copy_data=False,balanced_tree=balanced_tree) copy_data=False,balanced_tree=balanced_tree)
self._kdaxis = kdaxis self._kdaxis = kdaxis
self._kdradius = radius self._kdradius = radius
if report: if report:
@ -965,54 +940,50 @@ class Features3d:
if not self._kdtree: if not self._kdtree:
self.build_kdtree() self.build_kdtree()
if self._kdaxis==0: query_axis = [1,2] if self._kdaxis==0: query_axis = [1,2]
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) cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
if coords.shape[0]>12279:
print('Point coordinates:',coords[12279,:])
Ncoord = coords.shape[0] Ncoord = coords.shape[0]
raydir = np.zeros((3,),dtype=self._vertices.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) in_feat = np.empty((Ncoord,),dtype=np.int)
is_hit = np.empty((Ncoord,),dtype=np.int) is_hit = np.empty((Ncoord,),dtype=np.int)
hit_dir_ = np.empty((Ncoord,),dtype=np.int) hit_dir_ = np.empty((Ncoord,),dtype=np.int)
face_ = np.empty((Ncoord,),dtype=np.int) face_ = np.empty((Ncoord,),dtype=np.int)
print('point 12156 = ',coords[12156])
for ii in range(Ncoord): for ii in range(Ncoord):
hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir, hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir,
self._vertices[cand[ii],0,:], self._points[self._faces[cand[ii],1],:],
self._vertices[cand[ii],1,:], self._points[self._faces[cand[ii],2],:],
self._vertices[cand[ii],2,:]) self._points[self._faces[cand[ii],3],:])
# if ii==12279: print('DEBUG',hit_idx,t,N) # if ii==12279: print('DEBUG',hit_idx,t,N)
if ii==12156: print('DEBUG',hit_idx,t,hit_dir)
if hit_idx is None: if hit_idx is None:
in_feat[ii] = -1 in_feat[ii] = -1
is_hit[ii] = -1
hit_dir_[ii]= 0
face_[ii] = -1
else: else:
idx = np.argmin(np.abs(t)) idx = np.argmin(np.abs(t))
# print(cand[ii][hit_idx[idx]]) if ii==12156:
if ii==12279: for kk in hit_idx:
# print(cand[ii][idx],N[idx,:],N[idx,:]/np.sqrt(np.square(N[idx,:]).sum()),(N[idx,:]*raydir).sum(axis=-1),t[idx]*(N[idx,:]*raydir).sum(axis=-1)) print(cand[ii][kk])
print(t[idx],hit_dir[idx],raydir) # print(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:])
print(self._vertices[cand[ii][hit_idx[idx]],:,:]) # print(cand[ii][hit_idx[idx]])
hit_dir_[ii]= hit_dir[idx] # ifeat = self.feature_from_face(cand[ii][hit_idx[idx]])
face_[ii] = cand[ii][hit_idx[idx]] # print(self._offset[ifeat],self._offset[ifeat+1])
if hit_dir[idx]<0: # stop
in_feat[ii] = self.feature_from_cellidx(cand[ii][hit_idx[idx]]) # if hit_dir[idx]>0:
is_hit[ii] = 1 # in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]])
else: # else:
in_feat[ii] = -1 # in_feat[ii] = -1
is_hit[ii] = 0
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(in_feat>=0),Ncoord))
return in_feat #,is_hit,hit_dir_,face_ return in_feat
def feature_from_cellidx(self,idx_cell): def feature_from_face(self,idx_cell):
'''Gets feature ID for a given cell.''' '''Gets feature ID for a given cell.'''
# This is an optimized solution for monotonically increasing arrays: # This is an optimized solution for monotonically increasing arrays:
# as long as offset is smaller than the cell index, the boolean operation # as long as offset is smaller than the cell index, the boolean operation
@ -1020,13 +991,13 @@ 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.argmax(self._offset>=idx_cell) return np.searchsorted(self._offset,idx_cell,side='right')-1
def cellidx_from_feature(self,features,concat=True): def faces_from_feature(self,features,concat=True):
'''Returns indices of cells which belong to given features.''' '''Returns indices of cells which belong to given features.'''
from collections.abc import Iterable from collections.abc import Iterable
if features is None: if features is None:
idx = None idx = slice(None)
elif not isinstance(features,Iterable): elif not isinstance(features,Iterable):
idx = slice(self._offset[features],self._offset[features+1]) idx = slice(self._offset[features],self._offset[features+1])
elif concat: elif concat:
@ -1034,7 +1005,10 @@ class Features3d:
for feature in features: for feature in features:
rng_ = np.arange(self._offset[feature],self._offset[feature+1]) rng_ = np.arange(self._offset[feature],self._offset[feature+1])
idx.append(rng_) idx.append(rng_)
idx = np.concatenate(idx,axis=0) if len(idx)==1:
idx = np.array(idx)
else:
idx = np.concatenate(idx,axis=0)
else: else:
idx = [] idx = []
for feature in features: for feature in features:
@ -1042,7 +1016,7 @@ class Features3d:
idx.append(rng_) idx.append(rng_)
return idx return idx
def ncells_from_feature(self,features): def nfaces_of_feature(self,features):
'''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)
@ -1136,27 +1110,8 @@ class Features3d:
if len(hit_idx)==0: return (None,None,None) if len(hit_idx)==0: return (None,None,None)
# Intersection point is R+t*dR # Intersection point is R+t*dR
t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx] 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,np.sign(t)*np.sign(det[hit_idx])
return hit_idx,t,N[hit_idx]
def faces_points_representation(self,features,reduce_points=False):
from collections.abc import Iterable
# Select relevant cells
idx = self.cellidx_from_feature(features)
if reduce_points:
points,faces_ = np.unique(self._vertices[idx,:,:].reshape(-1,3),return_inverse=True,axis=0)
ncells = faces_.shape[0]//3
faces = np.empty((ncells,4),dtype=np.int)
faces[:,0] = 3
faces[:,1:] = faces_.reshape(ncells,3)
print('reduced:',faces.shape)
else:
points = self._vertices[idx,:,:].reshape((-1,3))
ncells = points.shape[0]//3
faces = np.empty((ncells,4),dtype=np.int)
faces[:,0] = 3
faces[:,1:] = np.arange(0,3*ncells).reshape(ncells,3)
print('full:',faces.shape)
return faces,points
def cmap_features(self,nfeatures,name='tab10'): def cmap_features(self,nfeatures,name='tab10'):
from matplotlib.colors import ListedColormap from matplotlib.colors import ListedColormap
@ -1181,24 +1136,31 @@ class Features3d:
icolor = (icolor+1)%ncolor icolor = (icolor+1)%ncolor
return ListedColormap(output) return ListedColormap(output)
def to_vtk(self,features,reduce_points=False,store_labels=False,remap_labels=False): def to_vtk(self,features):
'''Returns a VTK/pyvista object for the isocontour of a single feature,
list of features, or the full isocontour (when None is passed).'''
import pyvista as pv import pyvista as pv
faces,points = self.faces_points_representation(features,reduce_points=reduce_points) idx = self.faces_from_feature(features)
output = pv.PolyData(points,faces) return pv.PolyData(self._points,self._faces[idx,:])
if store_labels:
output.cell_arrays['label'] = np.empty(faces.shape[0],dtype=np.int) def to_vtk2(self,features):
offset = 0 import pyvista as pv
for ii,feature in enumerate(self.list_of_features(features)): idx = self.faces_from_feature(features)
ncells = self.ncells_from_feature(feature) faces_ = self._faces[idx,:]
if remap_labels: A = self._points[faces_[:,1],:]
output.cell_arrays['label'][offset:offset+ncells] = ii B = self._points[faces_[:,2],:]
else: C = self._points[faces_[:,3],:]
output.cell_arrays['label'][offset:offset+ncells] = feature cn = np.cross(B-A,C-A)
offset += ncells cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True))
cc = (A+B+C)/3
# print(A[0,:],B[0,:],C[0,:],cn[0,:])
# print(A.shape,B.shape,C.shape,cn.shape,faces_.shape,self._faces.shape)
output = pv.PolyData(cc)
output['cellnormal'] = cn
return output return output
def vtk_faces(self,face_idx):
import pyvista as pv
return pv.PolyData(self._points,self._faces[face_idx,:])
class BinaryFieldNd: class BinaryFieldNd:
def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False): def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):
assert isinstance(input,np.ndarray) and input.dtype==bool,\ assert isinstance(input,np.ndarray) and input.dtype==bool,\