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

bounding box optimized

parent cec1b993
......@@ -3,6 +3,7 @@ import numpy as np
import matplotlib.pyplot as plt
import sympy
import copy
import logging
class Node(object):
......@@ -43,17 +44,20 @@ class Node(object):
# initialize
self.mesh = copy.deepcopy(mesh)
if self.is_root():
cuts = copy.deepcopy(cuts)
cuts = copy.deepcopy(cuts) # dicts are mutable
self.remaining_cuts = cuts
log = logging.getLogger('Node')
log.debug(cuts)
self.delta = delta
if box is None:
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]}
vertices = np.array(self.mesh)
self.box = {'x': [min(vertices[:, 0]) - delta,
max(vertices[:, 0]) + delta],
'y': [min(vertices[:, 1]) - delta,
max(vertices[:, 1]) + delta],
'z': [min(vertices[:, 2]) - delta,
max(vertices[:, 2]) + delta]}
else:
self.box = box
self.left = None
......@@ -245,9 +249,6 @@ class Node(object):
box=right_box)
def _choose_next_cut(self):
"""
"""
largest = 0
for key in self.remaining_cuts:
if len(self.remaining_cuts[key]) > largest:
......@@ -328,6 +329,7 @@ class Node(object):
self._template = template
return self._template
class Searcher(Node):
def __init__(self, mesh, n_sections=None, delta=0.):
"""
......@@ -344,16 +346,18 @@ class Searcher(Node):
triangles = mesh.triangles().copy()
triangles.fields = []
indices = range(triangles.ntriangles())
a_indices = [i * 3 for i in indices]
centroids = triangles.centroids()
# radius towards the points that goes through all triangles corners
triangle_radii = np.linalg.norm(centroids - triangles[a_indices],
axis=1)
ab, ac = triangles.edges()
bc = ac - ab
side_lengths = np.concatenate([np.linalg.norm(side, axis=1)
for side in [ab, ac, bc]])
# import mplTools as mpl
# axis= tfields.plotting.gca(2)
# mpl.plotHistogram(side_lengths, axis=axis)
# mpl.plt.show()
# quit()
characteristic_side_length = np.max(side_lengths)
# abs
median_radius = np.median(triangle_radii)
cut_length = median_radius * 2 * 2
cut_length = characteristic_side_length * 1.1
elongation = maxima - minima
n_sections_auto = np.round(elongation / cut_length).astype(int)
......@@ -362,10 +366,11 @@ class Searcher(Node):
for i, ns in enumerate(n_sections):
if ns is not None:
n_sections_auto[i] = int(ns)
n_sections = n_sections_auto
cut = {}
for i, key in enumerate(['x', 'y', 'z']):
n_cuts = n_sections[1] + 1
n_cuts = n_sections[i] + 1
if n_cuts <= 2:
values = []
else:
......@@ -382,11 +387,12 @@ class Searcher(Node):
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])
>>> mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
>>> tree = tfields.bounding_box.Searcher(mesh, n_sections=[5, 5, 5])
>>> points = tfields.Tensors([[0.5, 1, 2.1],
... [0.5, 0, 0],
... [0.5, 2, 2]])
... [0.5, 2, 2.1],
... [0.5, 1.5, 2.5]])
>>> box_res = tree.in_faces(points, delta=0.0001)
>>> usual_res = mesh.in_faces(points, delta=0.0001)
>>> assert np.array_equal(box_res, usual_res)
......@@ -406,27 +412,19 @@ class Searcher(Node):
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!
leaf = self.find_leaf(point)
if leaf is None:
continue
if leaf.template.nfaces() == 0:
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
if __name__ == '__main__':
import doctest
# doctest.run_docstring_examples(Mesh3D.clean, globals())
......
......@@ -1714,8 +1714,9 @@ class TensorMaps(TensorFields):
# redirect maps
for mp_idx in range(len(self.maps)):
for f in range(len(self.maps[mp_idx])):
if tensor_index in self.maps[mp_idx][f]:
index = tfields.index(self.maps[mp_idx][f], tensor_index)
mp = np.array(self.maps[mp_idx], dtype=int)
if tensor_index in mp[f]:
index = tfields.index(mp[f], tensor_index)
inst.maps[mp_idx][f][index] = duplicate_index
return inst.removed(stale_mask)
......
......@@ -54,8 +54,7 @@ def duplicates(ar, axis=None):
array([0, 0, 2])
An empty sequence will not throw errors
>>> tfields.duplicates([], axis=0)
array([])
>>> assert np.array_equal(tfields.duplicates([], axis=0), [])
Returns:
list of int: int is pointing to first occurence of unique value
......
......@@ -52,37 +52,42 @@ def igrid(*base_vectors, **kwargs):
Examples:
Initilaize using the mgrid notation
>>> import tfields
>>> tfields.igrid((0, 1, 2j), (3, 4, 2j), (6, 7, 2j))
array([[ 0., 3., 6.],
[ 0., 3., 7.],
[ 0., 4., 6.],
[ 0., 4., 7.],
[ 1., 3., 6.],
[ 1., 3., 7.],
[ 1., 4., 6.],
[ 1., 4., 7.]])
>>> import numpy as np
>>> tfields.igrid([3, 4], np.linspace(0, 1, 2),
... (6, 7, 2), iter_order=[1, 0, 2])
array([[ 3., 0., 6.],
[ 3., 0., 7.],
[ 4., 0., 6.],
[ 4., 0., 7.],
[ 3., 1., 6.],
[ 3., 1., 7.],
[ 4., 1., 6.],
[ 4., 1., 7.]])
>>> tfields.igrid(np.linspace(0, 1, 2), np.linspace(3, 4, 2),
... np.linspace(6, 7, 2), iter_order=[2, 0, 1])
array([[ 0., 3., 6.],
[ 0., 4., 6.],
[ 1., 3., 6.],
[ 1., 4., 6.],
[ 0., 3., 7.],
[ 0., 4., 7.],
[ 1., 3., 7.],
[ 1., 4., 7.]])
>>> assert np.array_equal(tfields.igrid((0, 1, 2j),
... (3, 4, 2j),
... (6, 7, 2j)),
... [[ 0., 3., 6.],
... [ 0., 3., 7.],
... [ 0., 4., 6.],
... [ 0., 4., 7.],
... [ 1., 3., 6.],
... [ 1., 3., 7.],
... [ 1., 4., 6.],
... [ 1., 4., 7.]])
>>> assert np.array_equal(tfields.igrid([3, 4], np.linspace(0, 1, 2),
... (6, 7, 2),
... iter_order=[1, 0, 2]),
... [[ 3., 0., 6.],
... [ 3., 0., 7.],
... [ 4., 0., 6.],
... [ 4., 0., 7.],
... [ 3., 1., 6.],
... [ 3., 1., 7.],
... [ 4., 1., 6.],
... [ 4., 1., 7.]])
>>> assert np.array_equal(tfields.igrid(np.linspace(0, 1, 2),
... np.linspace(3, 4, 2),
... np.linspace(6, 7, 2),
... iter_order=[2, 0, 1]),
... [[ 0., 3., 6.],
... [ 0., 4., 6.],
... [ 1., 3., 6.],
... [ 1., 4., 6.],
... [ 0., 3., 7.],
... [ 0., 4., 7.],
... [ 1., 3., 7.],
... [ 1., 4., 7.]])
"""
iter_order = kwargs.pop('iter_order', np.arange(len(base_vectors)))
......
......@@ -27,19 +27,19 @@ def mode(array, axis=0, bins='auto', range=None):
Float usage:
>>> np.random.seed(seed=0) # deterministic random
>>> n = np.random.normal(3.1, 2., 1000)
>>> tfields.lib.stats.mode(n)
array([ 2.30838613])
>>> tfields.lib.stats.mode(n, bins='sturges')
array([ 2.81321206])
>>> tfields.lib.stats.mode(np.array([n, n]), axis=1)
array([[ 2.30838613],
[ 2.30838613]])
>>> assert np.isclose(tfields.lib.stats.mode(n), [ 2.30838613])
>>> assert np.isclose(tfields.lib.stats.mode(n, bins='sturges'),
... [ 2.81321206])
>>> assert np.allclose(tfields.lib.stats.mode(np.array([n, n]), axis=1),
... [[ 2.30838613],
... [ 2.30838613]])
>>> tfields.lib.stats.mode(np.array([n, n]), axis=0).shape
(1000, 1)
>>> tfields.lib.stats.mode(np.array([n, n]), axis=1).shape
(2, 1)
>>> tfields.lib.stats.mode(np.array([n, n]), axis=None)
array([ 2.81321206])
>>> assert np.isclose(tfields.lib.stats.mode(np.array([n, n]),
... axis=None),
... [ 2.81321206])
"""
array = np.asarray(array)
......
......@@ -606,9 +606,11 @@ class Mesh3D(tfields.TensorMaps):
mask = inst.evalf(expression)
# remove the points
# if not any(~mask) or all(~mask):
# inst = inst[mask]
if not any(~mask):
# no vertex is valid
inst = inst[mask]
elif all(~mask):
# all vertices are valid
inst = inst[mask]
elif at_intersection == 'keep':
expression_parts = tfields.lib.symbolics.split_expression(expression)
......
......@@ -538,6 +538,8 @@ class Triangles3D(tfields.TensorFields):
with warnings.catch_warnings():
warnings.filterwarnings('ignore', message="invalid value encountered in less_equal")
barycentric_bools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0))
if all(~barycentric_bools):
return barycentric_bools
if min_dist_method:
orthogonal_acceptance = np.full(barycentric_bools.shape, False)
closest_indices = np.argmin(abs(w)[barycentric_bools])
......
Markdown is supported
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