baustelle bounding box

tfields/__init__.py

@@ -21,3 +21,4 @@ from .mesh3D import Mesh3D  # NOQA
 from .mesh3D import fields_to_scalars, scalars_to_fields
 from .triangles3D import Triangles3D  # NOQA
 from .planes3D import Planes3D  # NOQA
+from .bounding_box import Node

tfields/bounding_box.py

+import tfields
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy
+import copy
+
+
+class Node(object):
+    """
+    This class allows to increase the performance with cuts in x,y and z direction
+    An extension to arbitrary cuts might be possible in the future
+
+    Args:
+        parent: Parent node of self
+        mesh: Mesh corresponding to the node
+        cut_expr: Cut that determines the seperation in left and right node
+        cuts: List of cuts for the children nodes
+
+    Examples:
+        >>> import tfields
+        >>> mesh = tfields.Mesh3D.grid((5.6, 6.2, 3),
+        ...                            (-0.25, 0.25, 4),
+        ...                            (-1, 1, 10))
+        
+        >>> cuts = {'x': [5.7, 6.1],
+        ...         'y': [-0.2, 0, 0.2],
+        ...         'z': [-0.5, 0.5]}
+
+        >>> tree = tfields.bounding_box.Node(mesh,
+        ...                                  cuts,
+        ...                                  at_intersection='keep')
+        >>> 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]
+        >>> special_leaf = tree.find_leaf([5.65, -0.21, 0])
+
+    """
+    def __init__(self, mesh, cuts,
+                 coord_sys=None, at_intersection="split",
+                 parent=None, box=None, internal_template=None, cut_expr=None):
+        self.parent = parent
+        # initialize
+        self.mesh = copy.deepcopy(mesh)
+        if self.is_root():
+            cuts = copy.deepcopy(cuts)
+        self.remaining_cuts = cuts
+
+        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])]}
+        else:
+            self.box = box
+        self.left = None
+        self.right = None
+        self.at_intersection = at_intersection
+        self._internal_template = internal_template
+        if self.is_root():
+            self._trim_to_box()
+        self.cut_expr = cut_expr
+        self.left_template = None
+        self.right_template = None
+        self.coord_sys = coord_sys
+        # start the splitting process
+        self._build()
+
+    def _build(self):
+        if not self.is_last_cut():
+            self._choose_next_cut()
+            self._split()
+
+    def is_leaf(self):
+        if self.left is None and self.right is None:
+            return True
+        else:
+            return False
+
+    def is_root(self):
+        if self.parent is None:
+            return True
+        else:
+            return False
+
+    def is_last_cut(self):
+        for key in self.remaining_cuts:
+            if len(self.remaining_cuts[key]) != 0:
+                return False
+        return True
+
+    @classmethod
+    def sort_leaves(cls, leaves_list):
+        """
+        sorting the leaves first in x, then y, then z direction
+        """
+        sorted_leaves = sorted(leaves_list,
+                               key=lambda x: (x.box['x'][1], x.box['y'][1], x.box['z'][1]))
+        return sorted_leaves
+
+    def _trim_to_box(self):
+        # 6 cuts to remove outer part of the box
+        x, y, z = sympy.symbols('x y z')
+        eps = 0.0000000001
+        x_cut = (float(self.box['x'][0] - eps) <= x) & (x <= float(self.box['x'][1] + eps))
+        y_cut = (float(self.box['y'][0] - eps) <= y) & (y <= float(self.box['y'][1] + eps))
+        z_cut = (float(self.box['z'][0] - eps) <= z) & (z <= float(self.box['z'][1] + eps))
+        section_cut = x_cut & y_cut & z_cut
+
+        self.mesh, self._internal_template = self.mesh.cut(
+            section_cut,
+            at_intersection=self.at_intersection,
+            return_template=True)
+
+    def leaves(self):
+        """
+        Recursive function to create a list of all leaves
+        Returns:
+            list: of all leaves descending from this node
+        """
+        if self.is_leaf():
+            return [self]
+        else:
+            if self.left is not None:
+                leftLeaves = self.left.leaves()
+            else:
+                leftLeaves = []
+            if self.right is not None:
+                rightLeaves = self.right.leaves()
+            else:
+                rightLeaves = []
+            return tfields.lib.util.flatten(leftLeaves + rightLeaves)
+
+    def find_leaf(self, point):
+        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)
+
+    def _split(self):
+        """
+        Split the node, if there is no cut_expr set and remaing cuts exist.
+        """
+        if self.cut_expr is None and self.remaining_cuts is None:
+            raise RuntimeError("Cannot split the mesh without cut_expr and"
+                               "remaining_cuts")
+        else:
+            # create cut expression
+            x, y, z = sympy.symbols('x y z')
+            if 'x' in self.cut_expr:
+                left_cut_expression = x <= self.cut_expr['x']
+                right_cut_expression = x >= self.cut_expr['x']
+                key = 'x'
+            elif 'y' in self.cut_expr:
+                left_cut_expression = y <= self.cut_expr['y']
+                right_cut_expression = y >= self.cut_expr['y']
+                key = 'y'
+            elif 'z' in self.cut_expr:
+                left_cut_expression = z <= self.cut_expr['z']
+                right_cut_expression = z >= self.cut_expr['z']
+                key = 'z'
+            else:
+                raise KeyError()
+
+            # split the cuts into left / right
+            left_cuts = self.remaining_cuts.copy()
+            right_cuts = self.remaining_cuts.copy()
+            left_cuts[key] = [
+                value
+                for value in self.remaining_cuts[key]
+                if value <= self.cut_expr[key]
+            ]
+            right_cuts[key] = [
+                value
+                for value in self.remaining_cuts[key]
+                if value > self.cut_expr[key]
+            ]
+            left_box = copy.deepcopy(self.box)
+            right_box = copy.deepcopy(self.box)
+            left_box[key][1] = self.cut_expr[key]
+            right_box[key][0] = self.cut_expr[key]
+
+            # actually cut!
+            left_mesh, self.left_template = self.mesh.cut(
+                left_cut_expression,
+                at_intersection=self.at_intersection,
+                return_template=True)
+            right_mesh, self.right_template = self.mesh.cut(
+                right_cut_expression,
+                at_intersection=self.at_intersection,
+                return_template=True)
+
+            # two new Nodes
+            self.left = Node(left_mesh,
+                             left_cuts,
+                             parent=self, 
+                             internal_template=self.left_template, cut_expr=None,
+                             coord_sys=self.coord_sys,
+                             at_intersection=self.at_intersection,
+                             box=left_box)
+            self.right = Node(right_mesh,
+                              right_cuts,
+                              parent=self, 
+                              internal_template=self.right_template, cut_expr=None,
+                              coord_sys=self.coord_sys,
+                              at_intersection=self.at_intersection,
+                              box=right_box)
+
+    def _choose_next_cut(self):
+        """
+
+        """
+        largest = 0
+        for key in self.remaining_cuts:
+            if len(self.remaining_cuts[key]) > largest:
+                largest = len(self.remaining_cuts[key])
+                largest_key = key
+
+        median = sorted(self.remaining_cuts[largest_key])[int(0.5 * (len(self.remaining_cuts[largest_key]) - 1))]
+        # pop median cut from remaining cuts
+        self.remaining_cuts[largest_key] = [
+            x for x in self.remaining_cuts[largest_key] if x != median
+        ]
+        self.cut_expr = {largest_key: median}
+
+    def _convert_map_index(self, index):
+        """
+        Recursively getting the map fields index from root
+        Args:
+            index (int): map field index on leaf
+                (index with respect to parent node, not to root)
+        Returns:
+            int: map field index
+        """
+        if self.is_root():
+            return index
+        else:
+            return_value = self.parent._convert_map_index(
+                self.parent._internal_template.maps[0].fields[0][int(index)]
+            )
+            return return_value
+
+    def _convert_field_index(self, index):
+        """
+        Recursively getting the fields index from root
+        Args:
+            index (int): field index on leaf
+                (index with respect to parent node, not to root)
+        Returns:
+            int: field index
+        """
+        if self.is_root():
+            return index
+        else:
+            return_value = self.parent._convert_field_index(
+                self.parent._internal_template.fields[0][int(index)]
+            )
+            return return_value
+
+    def template(self):
+        """
+        Get the global template for a leaf. This can be applied to the root mesh
+        with the cut method to retrieve exactly this leaf mesh again.
+        Returns:
+            tfields.Mesh3D: mesh with first scalars as an instruction on how to build
+                this cut (scalars point to faceIndices on mother mesh). Can be
+                used with Mesh3D.cut
+        """
+
+        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
+
+class Searcher(Node):
+    def __init__(self, mesh, n_sections=None):
+        """
+        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
+
+        if n_sections is None or any([ns is None for ns in n_sections]):
+            # project all triangels to have one point at 0, 0, 0
+            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)
+
+            # abs
+            median_radius = np.median(triangle_radii)
+            cut_length = median_radius * 2 * 2
+
+            elongation = maxima - minima
+            n_sections_auto = np.round(elongation / cut_length).astype(int)
+
+            if n_sections is not None:
+                for i, ns in enumerate(n_sections):
+                    if ns is not None:
+                        n_sections_auto[i] = int(ns)
+
+        cut = {}
+        for i, key in enumerate(['x', 'y', 'z']):
+            n_cuts = n_sections[1] + 1
+            if n_cuts <= 2:
+                values = []
+            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')
+
+    def in_faces(self, tensors, delta, assign_multiple=False):
+        """
+        TODO:
+            * check rare case of point+-delta outside box
+            * combine delta with self._delta
+        Examples:
+            >>> import tfields
+            >>> 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],
+            ...                           [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)
+
+        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)
+                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__
+            """
+        return masks
+
+
+
+if __name__ == '__main__':
+    import doctest
+    # doctest.run_docstring_examples(Mesh3D.clean, globals())
+    doctest.testmod()

tfields/lib/__init__.py

@@ -53,9 +53,15 @@ def duplicates(ar, axis=None):
         >>> tfields.duplicates(a, axis=0)
         array([0, 0, 2])
 
+        An empty sequence will not throw errors
+        >>> tfields.duplicates([], axis=0)
+        array([])
+
     Returns:
         list of int: int is pointing to first occurence of unique value
     """
+    if len(ar) == 0:
+        return np.array([])
     if axis != 0:
         raise NotImplementedError()
     sidx = np.lexsort(ar.T)

tfields/lib/decorators.py

@@ -53,6 +53,7 @@ class cached_property(object):
             if self.ttl > 0 and now - last_update > self.ttl:
                 raise AttributeError
         except (KeyError, AttributeError, TypeError):
+            # no cache entry found
             value = self.fget(inst)
             try:
                 if inst._cache is None:

tfields/mesh3D.py

@@ -77,6 +77,7 @@ 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)
@@ -503,6 +504,14 @@ class Mesh3D(tfields.TensorMaps):
         return tfields.Triangles3D(tris, *fields)
 
     def triangles(self):
+        """
+        Cached method to retrieve the triangles, belonging to this mesh
+        Examples:
+            >>> import tfields
+            >>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
+            >>> assert mesh.triangles() is mesh.triangles()
+
+        """
         return self._triangles
 
     def centroids(self):
@@ -598,8 +607,37 @@ class Mesh3D(tfields.TensorMaps):
         mask = inst.evalf(expression)
         # remove the points
 
-        if not any(~mask) or all(~mask):
+        # if not any(~mask) or all(~mask):
+        #     inst = inst[mask]
+        if not any(~mask):
             inst = inst[mask]
+        elif at_intersection == 'keep':
+            expression_parts = tfields.lib.symbolics.split_expression(expression)
+            if len(expression_parts) > 1:
+                new_mesh = inst.copy()
+                for exprPart in expression_parts:
+                    inst, _ = inst._cut_sympy(exprPart,
+                                              at_intersection=at_intersection,
+                                              _in_recursion=True)
+            elif len(expression_parts) == 1:
+                face_delete_indices = set([])
+                for i, face in enumerate(inst.maps[0]):
+                    """
+                    vertices_rejected is a mask for each face that is True, where
+                    a Point is on the rejected side of the plane
+                    """
+                    vertices_rejected = [~mask[f] for f in face]
+                    if all(vertices_rejected):
+                        # delete face
+                        face_delete_indices.add(i)
+                mask = np.full(len(inst.maps[0]), True, dtype=bool)
+                for face_idx in range(len(inst.maps[0])):
+                    if face_idx in face_delete_indices:
+                        mask[face_idx] = False
+                inst.maps[0] = inst.maps[0][mask]
+            else:
+                raise ValueError("Sympy expression is not splitable.")
+            inst = inst.cleaned()
         elif at_intersection == 'split' or at_intersection == 'splitRough':
             '''
             add vertices and faces that are at the border of the cuts
@@ -672,9 +710,9 @@ class Mesh3D(tfields.TensorMaps):
                         """
                         Add the intersection points and faces
                         """
-                        newP = _intersect(triangle_points, plane, vertices_rejected)
+                        intersection = _intersect(triangle_points, plane, vertices_rejected)
                         last_idx = len(vertices) - 1
-                        for tri_list in newP:
+                        for tri_list in intersection:
                             new_face = []
                             for item in tri_list:
                                 if isinstance(item, int):
@@ -827,8 +865,9 @@ class Mesh3D(tfields.TensorMaps):
                     indices in the template.maps[i].fields.
             coord_sys (coordinate system to cut in):
             at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
-                Options:    "remove" (Default)
-                            "split" - Create new triangles that make up the old one.
+                Options:    'remove' (Default) - remove the faces that are on the edge
+                            'keep' - keep the faces that are on the edge
+                            'split' - Create new triangles that make up the old one.
             return_template (bool): If True: return the template
                             to redo the same cut fast
         Examples:
@@ -836,7 +875,7 @@ class Mesh3D(tfields.TensorMaps):
             >>> import numpy as np
             >>> import tfields
             >>> from sympy.abc import x,y,z
-            >>> cutExpr = x > 1.5
+            >>> cut_expr = x > 1.5
 
             >>> m = tfields.Mesh3D.grid((0, 3, 4),
             ...                         (0, 3, 4),
@@ -847,7 +886,7 @@ class Mesh3D(tfields.TensorMaps):
             ...     tfields.Tensors(np.linspace(0,
             ...                                 len(m.maps[0]) - 1,
             ...                                 len(m.maps[0]))))
-            >>> mNew = m.cut(cutExpr)
+            >>> mNew = m.cut(cut_expr)
             >>> len(mNew)
             8
             >>> mNew.nfaces()
@@ -855,8 +894,16 @@ class Mesh3D(tfields.TensorMaps):
             >>> float(mNew[:, 0].min())
             2.0
 
-            Cutting with the split option will create new triangles on the edge:
-            >>> m_split = m.cut(cutExpr, at_intersection='split')
+            Cutting with the 'keep' option will leave triangles on the edge
+            untouched:
+            >>> m_keep = m.cut(cut_expr, at_intersection='keep')
+            >>> float(m_keep[:, 0].min())
+            1.0
+            >>> m_keep.nfaces()
+            12
+
+            Cutting with the 'split' option will create new triangles on the edge:
+            >>> m_split = m.cut(cut_expr, at_intersection='split')
             >>> float(m_split[:, 0].min())
             1.5
             >>> len(m_split)
@@ -866,7 +913,7 @@ class Mesh3D(tfields.TensorMaps):
 
             Cut with 'return_template=True' will return the exact same mesh but
             additionally an instruction to conduct the exact same cut fast (template)
-            >>> m_split_2, template = m.cut(cutExpr, at_intersection='split',
+            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
             ...                                    return_template=True)
             >>> m_split_template = m.cut(template)
             >>> assert m_split.equal(m_split_2, equal_nan=True)
@@ -935,6 +982,129 @@ 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):