implemented feature search. it should work now and perform well.
This commit is contained in:
parent
cee9a661c8
commit
6e3564c3c1
129
field.py
129
field.py
|
|
@ -678,6 +678,10 @@ class Features3d:
|
||||||
self.triangulate(contour_method=contour_method,
|
self.triangulate(contour_method=contour_method,
|
||||||
cellvol_normal_component=cellvol_normal_component,
|
cellvol_normal_component=cellvol_normal_component,
|
||||||
report=report)
|
report=report)
|
||||||
|
# Set some state variable
|
||||||
|
self._kdtree = None
|
||||||
|
self._kdaxis = None
|
||||||
|
self._kdradius = None
|
||||||
# Drop input if requested: this may free memory in case it was copied/temporary.
|
# Drop input if requested: this may free memory in case it was copied/temporary.
|
||||||
if not keep_input: self._input = None
|
if not keep_input: self._input = None
|
||||||
return
|
return
|
||||||
|
|
@ -870,7 +874,37 @@ class Features3d:
|
||||||
'''Returns an array with volumes of all features.'''
|
'''Returns an array with volumes of all features.'''
|
||||||
return np.add.reduceat(self._cell_volumes,self._offset[:-1])
|
return np.add.reduceat(self._cell_volumes,self._offset[:-1])
|
||||||
|
|
||||||
def find_feature(self,pts):
|
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
|
||||||
|
if kdaxis==0:
|
||||||
|
min1 = np.amin(self._vertices[:,:,1],axis=1)
|
||||||
|
max1 = np.amax(self._vertices[:,:,1],axis=1)
|
||||||
|
min2 = np.amin(self._vertices[:,:,2],axis=1)
|
||||||
|
max2 = np.amax(self._vertices[:,:,2],axis=1)
|
||||||
|
elif kdaxis==1:
|
||||||
|
min1 = np.amin(self._vertices[:,:,0],axis=1)
|
||||||
|
max1 = np.amax(self._vertices[:,:,0],axis=1)
|
||||||
|
min2 = np.amin(self._vertices[:,:,2],axis=1)
|
||||||
|
max2 = np.amax(self._vertices[:,:,2],axis=1)
|
||||||
|
elif kdaxis==2:
|
||||||
|
min1 = np.amin(self._vertices[:,:,0],axis=1)
|
||||||
|
max1 = np.amax(self._vertices[:,:,0],axis=1)
|
||||||
|
min2 = np.amin(self._vertices[:,:,1],axis=1)
|
||||||
|
max2 = np.amax(self._vertices[:,:,1],axis=1)
|
||||||
|
else:
|
||||||
|
raise ValueError("Invalid ray axis.")
|
||||||
|
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:
|
||||||
|
print('[Features3d.build_kdtree]',end=' ')
|
||||||
|
print('KD-Tree built for axis {}.'.format(axis))
|
||||||
|
|
||||||
|
def inside_feature(self,coords):
|
||||||
'''Find nearest triangle from point R in direction ±dR and its orientation w.r.t dR.
|
'''Find nearest triangle from point R in direction ±dR and its orientation w.r.t dR.
|
||||||
The algorithm determines the closest intersection of a ray cast from the origin of the
|
The algorithm determines the closest intersection of a ray cast from the origin of the
|
||||||
point in the specified direction (with positiv and negative sign). The ray is supposed to
|
point in the specified direction (with positiv and negative sign). The ray is supposed to
|
||||||
|
|
@ -889,80 +923,39 @@ 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]
|
||||||
'''
|
'''
|
||||||
# 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."
|
||||||
from time import time
|
|
||||||
from scipy import spatial
|
|
||||||
|
|
||||||
# pts = np.random.random((10000,3))
|
if not self._kdtree:
|
||||||
|
self.build_kdtree()
|
||||||
|
if self._kdaxis==0: query_axis = [1,2]
|
||||||
|
elif self._kdaxis==1: query_axis = [0,2]
|
||||||
|
elif self._kdaxis==2: query_axis = [0,1]
|
||||||
|
|
||||||
t = time()
|
cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
|
||||||
|
|
||||||
xmin = np.amin(self._vertices[:,:,0],axis=1)
|
N = coords.shape[0]
|
||||||
xmax = np.amax(self._vertices[:,:,0],axis=1)
|
raydir = np.zeros((3,),dtype=self._vertices.dtype)
|
||||||
zmin = np.amin(self._vertices[:,:,2],axis=1)
|
raydir[self._kdaxis] = 1.0
|
||||||
zmax = np.amax(self._vertices[:,:,2],axis=1)
|
in_feat = np.empty((N,),dtype=np.int)
|
||||||
print(time()-t)
|
for ii in range(N):
|
||||||
# print(xmin.shape,xmin2.shape)
|
hit_idx,t,N = Features3d.ray_triangle_intersection(coords[ii],raydir,
|
||||||
# print(np.all(np.isclose(xmin,xmin2)))
|
self._vertices[cand[ii],0,:],
|
||||||
|
self._vertices[cand[ii],1,:],
|
||||||
center = np.stack((0.5*(xmax+xmin),0.5*(zmax+zmin)),axis=1)
|
self._vertices[cand[ii],2,:])
|
||||||
radius = 0.5*np.amax(np.sqrt((xmax-xmin)**2+(zmax-zmin)**2))
|
|
||||||
print(time()-t)
|
|
||||||
del xmin,xmax,zmin,zmax
|
|
||||||
print('Preparation for KD tree in',time()-t)
|
|
||||||
|
|
||||||
t = time()
|
|
||||||
kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False)
|
|
||||||
# kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
|
||||||
# kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False)
|
|
||||||
# kd = spatial.KDTree(self._points[:,1:],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
|
||||||
print('KD tree built in',time()-t)
|
|
||||||
|
|
||||||
# t = time()
|
|
||||||
# map_ = self._points.shape[0]*[[]]
|
|
||||||
# for ii,face in enumerate(self._faces[:,1:].ravel()):
|
|
||||||
# map_[face].append(ii)
|
|
||||||
# faces_connected_to_vertex = {}
|
|
||||||
# for face_index,face in enumerate(self._faces[:,1:]):
|
|
||||||
# for vertex_index in face:
|
|
||||||
# faces_connected_to_vertex.setdefault(vertex_index,[]).append(face_index)
|
|
||||||
# print('Inverted cells-faces',time()-t)
|
|
||||||
# radius = np.sqrt(np.sum(self.spacing[1:]**2))
|
|
||||||
|
|
||||||
t = time()
|
|
||||||
bla = kd.query_ball_point(pts[:,[0,2]],radius)
|
|
||||||
print('KD tree query',time()-t)
|
|
||||||
|
|
||||||
t__ = time()
|
|
||||||
Npts = pts.shape[0]
|
|
||||||
print(pts.shape,Npts)
|
|
||||||
dR = (0,1,0)
|
|
||||||
is_inside = np.empty((Npts,),dtype=bool)
|
|
||||||
print(is_inside.shape)
|
|
||||||
for ii in range(Npts):
|
|
||||||
hit_idx,t,N = Features3d.ray_triangle_intersection(pts[ii],dR,
|
|
||||||
self._vertices[bla[ii],0,:],
|
|
||||||
self._vertices[bla[ii],1,:],
|
|
||||||
self._vertices[bla[ii],2,:])
|
|
||||||
if hit_idx is None:
|
if hit_idx is None:
|
||||||
is_inside[ii] = False
|
in_feat[ii] = -1
|
||||||
else:
|
else:
|
||||||
idx = np.argmin(np.abs(t))
|
idx = np.argmin(np.abs(t))
|
||||||
is_inside[ii] = t[idx]*(N[idx,:]*dR).sum(axis=-1)>0
|
if t[idx]*(N[idx,:]*raydir).sum(axis=-1)>0:
|
||||||
print('inside check',time()-t__)
|
in_feat[ii] = self.feature_by_cell(cand[ii][idx])
|
||||||
|
else:
|
||||||
|
in_feat[ii] = -1
|
||||||
|
return in_feat
|
||||||
|
|
||||||
# print(is_inside)
|
def feature_by_cell(self,idx_cell):
|
||||||
# The provided vertices do not need to form a closed surface! This means they can be
|
'''Gets feature ID for a given cell.'''
|
||||||
# filtered before being passed to this function.
|
return np.argmax(idx_cell>self._offset)
|
||||||
|
|
||||||
# print(radius)
|
|
||||||
# print(len(kd.query_ball_point((0.5,0.0),radius,return_sorted=True)),self._points.shape[0])
|
|
||||||
|
|
||||||
|
|
||||||
# yset = set(x.data for x in sorted(treey.at(0.0)))
|
|
||||||
# zset = set(x.data for x in sorted(treez.at(0.0)))
|
|
||||||
return is_inside
|
|
||||||
|
|
||||||
@staticmethod
|
@staticmethod
|
||||||
def ray_triangle_intersection(R,dR,A,B,C):
|
def ray_triangle_intersection(R,dR,A,B,C):
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue