Commit cec1b993 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

bounding box is working

parent 393877d5
......@@ -32,12 +32,12 @@ class Node(object):
>>> leaves = tree.leaves()
>>> leaves = tfields.bounding_box.Node.sort_leaves(leaves)
>>> meshes = [leaf.mesh for leaf in leaves]
>>> templates = [leaf.template() for leaf in leaves]
>>> templates = [leaf.template for leaf in leaves]
>>> special_leaf = tree.find_leaf([5.65, -0.21, 0])
"""
def __init__(self, mesh, cuts,
coord_sys=None, at_intersection="split",
coord_sys=None, at_intersection="split", delta=0.,
parent=None, box=None, internal_template=None, cut_expr=None):
self.parent = parent
# initialize
......@@ -46,16 +46,22 @@ class Node(object):
cuts = copy.deepcopy(cuts)
self.remaining_cuts = cuts
self.delta = delta
if box is None:
self.box = {'x': [min(self.mesh[:, 0]), max(self.mesh[:, 0])],
'y': [min(self.mesh[:, 1]), max(self.mesh[:, 1])],
'z': [min(self.mesh[:, 2]), max(self.mesh[:, 2])]}
self.box = {'x': [min(self.mesh[:, 0]) - delta,
max(self.mesh[:, 0]) + delta],
'y': [min(self.mesh[:, 1]) - delta,
max(self.mesh[:, 1]) + delta],
'z': [min(self.mesh[:, 2]) - delta,
max(self.mesh[:, 2]) + delta]}
else:
self.box = box
self.left = None
self.right = None
self.at_intersection = at_intersection
self._internal_template = internal_template
if self.is_leaf():
self._template = None
if self.is_root():
self._trim_to_box()
self.cut_expr = cut_expr
......@@ -88,6 +94,20 @@ class Node(object):
return False
return True
def in_box(self, point):
x, y, z = point
for key in ['x', 'y', 'z']:
value = locals()[key]
if value < self.box[key][0] or self.box[key][1] < value:
return False
return True
@property
def root(self):
if self.is_root:
return self
return self.parent.root
@classmethod
def sort_leaves(cls, leaves_list):
"""
......@@ -130,17 +150,30 @@ class Node(object):
rightLeaves = []
return tfields.lib.util.flatten(leftLeaves + rightLeaves)
def find_leaf(self, point):
def find_leaf(self, point, _in_recursion=False):
"""
Returns:
Node / None:
Node: leaf note, containinig point
None: point outside root box
"""
x, y, z = point
if self.is_root():
if not self.in_box(point):
return None
else:
if not _in_recursion:
raise RuntimeError("Only root node can search for all leaves")
if self.is_leaf():
return self
x, y, z = point
if len(self.cut_expr) > 1:
raise ValueError("cut_expr is too long")
key = self.cut_expr.keys()[0]
value = locals()[key]
if value <= self.cut_expr[key]:
return self.left.find_leaf(point)
return self.right.find_leaf(point)
return self.left.find_leaf(point, _in_recursion=True)
return self.right.find_leaf(point, _in_recursion=True)
def _split(self):
"""
......@@ -262,6 +295,7 @@ class Node(object):
)
return return_value
@property
def template(self):
"""
Get the global template for a leaf. This can be applied to the root mesh
......@@ -275,34 +309,33 @@ class Node(object):
if not self.is_leaf():
raise RuntimeError("Only leaf nodes can return a template")
template = self._internal_template.copy()
template_field = []
if template.fields:
for idx in template.fields[0]:
template_field.append(self._convert_field_index(idx))
template.fields = [tfields.Tensors(template_field, dim=1,
dtype=int)]
template_map_field = []
if len(template.maps[0]) > 0:
for idx in template.maps[0].fields[0]:
template_map_field.append(self._convert_map_index(idx))
template.maps[0].fields = [tfields.Tensors(template_map_field,
dim=1, dtype=int)]
return template
if self._template is None:
template = self._internal_template.copy()
template_field = []
if template.fields:
for idx in template.fields[0]:
template_field.append(self._convert_field_index(idx))
template.fields = [tfields.Tensors(template_field, dim=1,
dtype=int)]
template_map_field = []
if len(template.maps[0]) > 0:
for idx in template.maps[0].fields[0]:
template_map_field.append(self._convert_map_index(idx))
template.maps[0].fields = [tfields.Tensors(template_map_field,
dim=1, dtype=int)]
self._template = template
return self._template
class Searcher(Node):
def __init__(self, mesh, n_sections=None):
def __init__(self, mesh, n_sections=None, delta=0.):
"""
Args:
n_sections (None or list of len 3):
None: auto search n_sections
list: number of sections in the index dimension
"""
# mesh = mesh.copy()
# mesh.maps[0].fields = [tfields.Tensors(np.arange(len(mesh.maps[0])))]
minima = mesh.min(axis=0) - 0.0001
maxima = mesh.max(axis=0) + 0.0001
......@@ -338,34 +371,37 @@ class Searcher(Node):
else:
values = np.linspace(minima[i], maxima[i], n_cuts)[1:-1]
cut[key] = values
print cut
return super(Searcher, self).__init__(mesh, cut, at_intersection='keep')
return super(Searcher, self).__init__(mesh, cut, at_intersection='keep',
delta=delta)
def in_faces(self, tensors, delta, assign_multiple=False):
def in_faces(self, tensors, delta=-1, assign_multiple=False):
"""
TODO:
* check rare case of point+-delta outside box
* combine delta with self._delta
Examples:
>>> import tfields
>>> import numpy as np
>>> mesh = tfields.Mesh3D.grid((0, 1, 2), (2, 1, 2), (2, 3, 2))
>>> tree = tfields.bounding_box.Searcher(mesh, n_sections=[None, 2, None])
>>> points = tfields.Tensors([[0.5, 1, 2],
>>> points = tfields.Tensors([[0.5, 1, 2.1],
... [0.5, 0, 0],
... [0.5, 2, 2]])
>>> box_res = tree.in_faces(points, delta=0.0001)
>>> usual_res = mesh.in_faces(points, delta=0.0001)
>>> box_res, usual_res
>>> assert np.array_equal(box_res, usual_res)
"""
if not self.is_root():
raise ValueError("in_faces may only be called by root Node.")
if self.at_intersection != 'keep':
raise ValueError("Intersection method must be 'keep' for in_faces"
"method.")
if self.mesh.nfaces() == 0:
return np.empty((tensors.shape[0], 0), dtype=bool)
if delta == -1:
delta = self.delta
masks = np.zeros((len(tensors), self.mesh.nfaces()), dtype=bool)
with tensors.tmp_transform(self.mesh.coord_sys):
......@@ -375,30 +411,18 @@ class Searcher(Node):
# "to be in triangles (before handling assign_multiple)."
# .format(nP=masks.sum(), **locals())) # SLOW!
leaf = self.find_leaf(point)
leaf_mask = masks[i]
leaf_mask[leaf.mesh.maps[0].fields[0]] |= leaf.mesh.triangles()._in_triangles(point, delta)
masks[i] = leaf_mask
if len(masks) == 0:
# empty sequence
return masks
if not assign_multiple:
"""
index of first faces containing the point for each point. This gets the
first membership index always. ignoring multiple triangle memberships
"""
faceMembershipIndices = np.argmax(masks, axis=1)
# True if point lies in any triangle
membershipBools = (masks.sum(axis=1) != 0)
masks = np.full(masks.shape, False, dtype=bool)
for j, valid in enumerate(membershipBools):
if valid:
masks[j, faceMembershipIndices[j]] = True
"""
masks is now the matrix as described in __doc__
"""
if leaf is None:
continue
leaf_mask = leaf.template.triangles()._in_triangles(point, delta)
original_face_indices = leaf.template.maps[0].fields[0][leaf_mask]
if not assign_multiple and len(original_face_indices) > 0:
original_face_indices = original_face_indices[:1]
# axis = tfields.plotting.gca(3)
# tfields.Tensors([point]).plot(axis=axis)
# leaf.mesh.plot(axis=axis, edgecolor='k')
# tfields.plotting.show()
masks[i, original_face_indices] = True
return masks
......
......@@ -77,7 +77,6 @@ def _intersect(triangle, plane, vertices_rejected):
TODO:
align norm vectors with previous face
"""
print triangle, plane, vertices_rejected
nTrue = vertices_rejected.count(True)
lonely_bool = True if nTrue == 1 else False
index = vertices_rejected.index(lonely_bool)
......@@ -982,129 +981,6 @@ class Mesh3D(tfields.TensorMaps):
return tfields.plotting.plot_mesh(self, self.faces, **kwargs)
class BoxSearcher(object):
"""
Note:
* too small boxes
Examples:
>>> import tfields
>>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
>>> box_searcher = tfields.mesh3D.BoxSearcher(mesh)
>>> box_searcher._find_part_index([1, 2, 3])
728
>>> points = tfields.Tensors([[0.5, 1, 2],
... [0.5, 0, 0],
... [0.5, 2, 2]])
>>> box_res = box_searcher.in_faces(points,
... delta=0.0001)
>>> usual_res = mesh.in_faces(points, delta=0.0001)
>>> usual_res, (box_res == usual_res)
"""
def __init__(self, mesh, nx=10, ny=10, nz=10):
self._mesh = mesh.copy()
self._mesh.maps[0].fields = [tfields.Tensors(np.arange(len(mesh.maps[0])))]
self._minima = self._mesh.min(axis=0) - 0.0001
self._maxima = self._mesh.max(axis=0) + 0.0001
self._nx = nx
self._ny = ny
self._nz = nz
self._x_space = np.linspace(self._minima[0], self._maxima[0], self._nx)
self._y_space = np.linspace(self._minima[1], self._maxima[1], self._ny)
self._z_space = np.linspace(self._minima[2], self._maxima[2], self._nz)
if self._mesh.triangles().ntriangles() != 0:
self._build_parts()
def _build_parts(self):
from sympy.abc import x, y, z
self._parts = []
for x_min, x_max in tfields.lib.util.pairwise(self._x_space):
x_min, x_max = self.v_min_v_max(x_min, x_max, 0)
x_part = self._mesh.cut((x > x_min) & (x < x_max), at_intersection='keep')
for y_min, y_max in tfields.lib.util.pairwise(self._y_space):
y_min, y_max = self.v_min_v_max(y_min, y_max, 1)
xy_part = x_part.cut((y > y_min) & (y < y_max), at_intersection='keep')
for z_min, z_max in tfields.lib.util.pairwise(self._z_space):
z_min, z_max = self.v_min_v_max(z_min, z_max, 2)
xyz_part = xy_part.cut((z > z_min) & (z < z_max), at_intersection='keep')
# cut_expression = ((x > x_min) & (x <= x_max) &
# (y > y_min) & (y <= y_max) &
# (z > z_min) & (z <= z_max))
# xyz_part = self._mesh.cut(cut_expression, at_intersection='keep')
self._parts.append(xyz_part)
def v_min_v_max(self, v_min, v_max, v_index):
if v_min == self._minima[v_index]:
v_min = -np.finfo(float).max
if v_max == self._maxima[v_index]:
v_max = np.finfo(float).max
return v_min, v_max
def _in_between(self, value, v_min, v_max, v_index):
v_min, v_max = self.v_min_v_max(v_min, v_max, v_index)
return (v_min < value and value <= v_max)
def _find_part_index(self, point):
x, y, z = point
index = 0
for x_min, x_max in tfields.lib.util.pairwise(self._x_space):
for y_min, y_max in tfields.lib.util.pairwise(self._y_space):
for z_min, z_max in tfields.lib.util.pairwise(self._z_space):
if (self._in_between(x, x_min, x_max, 0) and
self._in_between(y, y_min, y_max, 1) and
self._in_between(z, z_min, z_max, 2)):
return index
index += 1
def in_faces(self, tensors, delta, assign_multiple=False):
"""
TODO:
* check rare case of point+-delta outside box
* combine delta with self._delta
"""
if self._mesh.nfaces() == 0:
return np.empty((tensors.shape[0], 0), dtype=bool)
masks = np.zeros((len(tensors), self._mesh.nfaces()), dtype=bool)
with tensors.tmp_transform(self._mesh.coord_sys):
for i, point in enumerate(iter(tensors)):
# print(i)
# print("Status: {i} points processed. Yet {nP} are found "
# "to be in triangles (before handling assign_multiple)."
# .format(nP=masks.sum(), **locals())) # SLOW!
part_mesh = self._parts[self._find_part_index(point)]
part_mask = masks[i]
part_mask[part_mesh.maps[0].fields[0]] |= part_mesh.triangles()._in_triangles(point, delta)
masks[i] = part_mask
# print("Status: All points processed. Yet {nP} are found "
# "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
if len(masks) == 0:
# empty sequence
return masks
if not assign_multiple:
"""
index of first faces containing the point for each point. This gets the
first membership index always. ignoring multiple triangle memberships
"""
faceMembershipIndices = np.argmax(masks, axis=1)
# True if point lies in any triangle
membershipBools = (masks.sum(axis=1) != 0)
masks = np.full(masks.shape, False, dtype=bool)
for j, valid in enumerate(membershipBools):
if valid:
masks[j, faceMembershipIndices[j]] = True
"""
masks is now the matrix as described in __doc__
"""
return masks
if __name__ == '__main__': # pragma: no cover
import doctest
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment