mesh3D.py 39.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#!/usr/bin/env
# encoding: utf-8
"""
Author:     Daniel Boeckenhoff
Mail:       daniel.boeckenhoff@ipp.mpg.de

part of tfields library
"""
import numpy as np
import sympy
11
import rna
12
import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
13
14

# obj imports
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
15
from tfields.lib.decorators import cached_property
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
16
17
import logging
import os
18
19


20
21
22
23
24
25
def _dist_from_plane(point, plane):
    return plane['normal'].dot(point) + plane['d']


def _segment_plane_intersection(p0, p1, plane):
    """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
26
27
28
29
30
    Get the intersection between the ray spanned by p0 and p1 with the plane.
    Args:
        p0: array of length 3
        p1: array of length 3
        plane: 
31
32
    Returns:
        points, direction
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
33
            points (list of arrays of length 3): 3d points
34
35
36
    """
    distance0 = _dist_from_plane(p0, plane)
    distance1 = _dist_from_plane(p1, plane)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
37
38
    p0_on_plane = abs(distance0) < np.finfo(float).eps
    p1_on_plane = abs(distance1) < np.finfo(float).eps
39
    points = []
40
    direction = 0
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
41
    if p0_on_plane:
42
        points.append(p0)
43

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
44
    if p1_on_plane:
45
        points.append(p1)
46
    # remove duplicate points
47
48
    if len(points) > 1:
        points = np.unique(points, axis=0)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
49
    if p0_on_plane and p1_on_plane:
50
        return points, direction
51
52

    if distance0 * distance1 > np.finfo(float).eps:
53
        return points, direction
54
55
56

    direction = np.sign(distance0)
    if abs(distance0) < np.finfo(float).eps:
57
        return points, direction
58
    elif abs(distance1) < np.finfo(float).eps:
59
        return points, direction
60
61
62
    if abs(distance0 - distance1) > np.finfo(float).eps:
        t = distance0 / (distance0 - distance1)
    else:
63
        return points, direction
64

65
    points.append(p0 + t * (p1 - p0))
66
    # remove duplicate points
67
68
69
    if len(points) > 1:
        points = np.unique(points, axis=0)
    return points, direction
70
71
72
73
74
75
76


def _intersect(triangle, plane, vertices_rejected):
    """
    Intersect a triangle with a plane. Give the info, which side of the
    triangle is rejected by passing the mask vertices_rejected
    Returns:
77
78
79
80
81
82
83
84
85
        list of list. The inner list is of length 3 and refers to the points of
        new triangles. The reference is done with varying types:
            int: reference to triangle index
            complex: reference to duplicate point. This only happens in case
                two triangles are returned. Then only in the second triangle
            iterable: new vertex

    TODO:
        align norm vectors with previous face
86
    """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
87
88
    n_true = vertices_rejected.count(True)
    lonely_bool = True if n_true == 1 else False
89
90
91
92
93
    index = vertices_rejected.index(lonely_bool)
    s0, d0 = _segment_plane_intersection(triangle[0], triangle[1], plane)
    s1, d1 = _segment_plane_intersection(triangle[1], triangle[2], plane)
    s2, d2 = _segment_plane_intersection(triangle[2], triangle[0], plane)

94
95
96
    single_index = index
    couple_indices = [j for j in range(3)
                      if not vertices_rejected[j] == lonely_bool]
97
98
99
100
101
102

    # TODO handle special cases. For now triangles with at least two points on plane are excluded
    new_points = None

    if len(s0) == 2:
        # both points on plane
103
        return new_points
104
105
    if len(s1) == 2:
        # both points on plane
106
        return new_points
107
108
    if len(s2) == 2:
        # both points on plane
109
        return new_points
110
    if lonely_bool:
111
        # two new triangles
112
        if len(s0) == 1 and len(s1) == 1:
113
114
            new_points = [[couple_indices[0], s0[0], couple_indices[1]],
                          [couple_indices[1], complex(1), s1[0]]]
115
        elif len(s1) == 1 and len(s2) == 1:
116
117
            new_points = [[couple_indices[0], couple_indices[1], s1[0]],
                          [couple_indices[0], complex(2), s2[0]]]
118
        elif len(s0) == 1 and len(s2) == 1:
119
120
            new_points = [[couple_indices[0], couple_indices[1], s0[0]],
                          [couple_indices[1], s2[0], complex(2)]]
121
    else:
122
        # one new triangle
123
        if len(s0) == 1 and len(s1) == 1:
124
            new_points = [[single_index, s1[0], s0[0]]]
125
        elif len(s1) == 1 and len(s2) == 1:
126
            new_points = [[single_index, s2[0], s1[0]]]
127
        elif len(s0) == 1 and len(s2) == 1:
128
129
            new_points = [[single_index, s0[0], s2[0]]]
    return new_points
130
131


132
133
134
135
136
137
def scalars_to_fields(scalars):
    scalars = np.array(scalars)
    if len(scalars.shape) == 1:
        return [tfields.Tensors(scalars)]
    return [tfields.Tensors(fs) for fs in scalars]

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
138

139
140
141
def fields_to_scalars(fields):
    return np.array(fields)

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
142

143
def faces_to_maps(faces, *fields):
144
    return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
145

146

147
def maps_to_faces(maps):
148
149
150
151
152
    if len(maps) == 0:
        return np.array([])
    elif len(maps) > 1:
        raise NotImplementedError("Multiple maps")
    return np.array(maps[0])
153
154


155
156
157
158
159
class Mesh3D(tfields.TensorMaps):
    # pylint: disable=R0904
    """
    Points3D child used as vertices combined with faces to build a geometrical mesh of triangles
    Examples:
160
161
        >>> import tfields
        >>> import numpy as np
162
        >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], faces=[[0, 1, 2], [1, 2, 3]])
163
164
165
166
167
        >>> m.equal([[1, 2, 3],
        ...          [3, 3, 3],
        ...          [0, 0, 0],
        ...          [5, 6, 7]])
        True
168
        >>> np.array_equal(m.faces, [[0, 1, 2], [1, 2, 3]])
169
        True
170
171

        conversion to points only
172
173
174
175
176
        >>> tfields.Points3D(m).equal([[1, 2, 3],
        ...                            [3, 3, 3],
        ...                            [0, 0, 0],
        ...                            [5, 6, 7]])
        True
177
178
179
180
181
182

        Empty instances
        >>> m = tfields.Mesh3D([]);

        going from Mesh3D to Triangles3D instance is easy and will be cached.
        >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]);
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
183
        >>> assert m.triangles().equal(tfields.Triangles3D([[ 1.,  0.,  0.],
184
185
        ...                                               [ 0.,  1.,  0.],
        ...                                               [ 0.,  0.,  0.]]))
186
187

        a list of scalars is assigned to each face
188
189
190
        >>> mScalar = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]], faceScalars=[.5]);
        >>> np.array_equal(mScalar.faceScalars, [[ 0.5]])
        True
191
192

        adding together two meshes:
193
194
195
196
197
198
199
200
201
202
        >>> m2 = tfields.Mesh3D([[1,0,0],[2,0,0],[0,3,0]],
        ...                     faces=[[0,1,2]], faceScalars=[.7])
        >>> msum = tfields.Mesh3D.merged(mScalar, m2)
        >>> msum.equal([[ 1.,  0.,  0.],
        ...             [ 0.,  1.,  0.],
        ...             [ 0.,  0.,  0.],
        ...             [ 1.,  0.,  0.],
        ...             [ 2.,  0.,  0.],
        ...             [ 0.,  3.,  0.]])
        True
203
        >>> assert np.array_equal(msum.faces, [[0, 1, 2], [3, 4, 5]])
204
205
206
207

        Saving and reading
        >>> from tempfile import NamedTemporaryFile
        >>> outFile = NamedTemporaryFile(suffix='.npz')
208
        >>> m.save(outFile.name)
209
        >>> _ = outFile.seek(0)
210
        >>> m1 = tfields.Mesh3D.load(outFile.name)
211
212
        >>> bool(np.all(m == m1))
        True
213
        >>> assert np.array_equal(m1.faces, np.array([[0, 1, 2]]))
214
215

    """
216
    def __new__(cls, tensors, *fields, **kwargs):
217
218
        if not issubclass(type(tensors), Mesh3D):
            kwargs['dim'] = 3
219
220
221
        faces = kwargs.pop('faces', None)
        faceScalars = kwargs.pop('faceScalars', [])
        maps = kwargs.pop('maps', None)
222
        if maps is not None and faces is not None:
223
224
225
            raise ValueError("Conflicting options maps and faces")
        if maps is not None:
            kwargs['maps'] = maps
226
        if len(faceScalars) > 0:
227
228
229
230
231
232
233
234
235
236
237
238
239
            map_fields = scalars_to_fields(faceScalars)
        else:
            map_fields = []
        if faces is not None:
            kwargs['maps'] = faces_to_maps(faces,
                                           *map_fields)
        obj = super(Mesh3D, cls).__new__(cls, tensors, *fields, **kwargs)
        if len(obj.maps) > 1:
            raise ValueError("Mesh3D only allows one map")
        if obj.maps and obj.maps[0].dim != 3:
            raise ValueError("Face dimension should be 3")
        return obj

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
    def _save_obj(self, path, **kwargs):
        """
        Save obj as wavefront/.obj file
        """
        obj = kwargs.pop('object', None)
        group = kwargs.pop('group', None)

        cmap = kwargs.pop('cmap', 'viridis')
        map_index = kwargs.pop('map_index', None)

        path = path.replace('.obj', '')
        directory, name = os.path.split(path)

        if not (self.faceScalars.size == 0 or map_index is None):
            scalars = self.maps[0].fields[map_index]
            min_scalar = scalars[~np.isnan(scalars)].min()
            max_scalar = scalars[~np.isnan(scalars)].max()
            vmin = kwargs.pop('vmin', min_scalar)
            vmax = kwargs.pop('vmax', max_scalar)
            if vmin == vmax:
                if vmin == 0.:
                    vmax = 1.
                else:
                    vmin = 0.
Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
264
265
            import matplotlib.colors as colors
            import matplotlib.pyplot as plt
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
266
267
268
269
270
271
            norm = colors.Normalize(vmin, vmax)
            color_map = plt.get_cmap(cmap)
        else:
            # switch for not coloring the triangles and thus not producing the materials
            norm = None

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
272
273
274
        if len(kwargs) != 0:
            raise ValueError("Unused arguments.")

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
        if norm is not None:
            mat_name = name + '_frame_{0}.mat'.format(map_index)
            scalars[np.isnan(scalars)] = min_scalar - 1
            sorted_scalars = scalars[scalars.argsort()]
            sorted_scalars[sorted_scalars == min_scalar - 1] = np.nan
            sorted_faces = self.faces[scalars.argsort()]
            scalar_set = np.unique(sorted_scalars)
            scalar_set[scalar_set == min_scalar - 1] = np.nan
            mat_path = os.path.join(directory, mat_name)
            with open(mat_path, 'w') as mf:
                for s in scalar_set:
                    if np.isnan(s):
                        mf.write("newmtl nan")
                        mf.write("Kd 0 0 0\n\n")
                    else:
                        mf.write("newmtl mtl_{0}\n".format(s))
                        mf.write("Kd {c[0]} {c[1]} {c[2]}\n\n".format(c=color_map(norm(s))))
        else:
            sorted_faces = self.faces

        # writing of the obj file
        with open(path + '.obj', 'w') as f:
            f.write("# File saved with tfields Mesh3D._save_obj method\n\n")
            if norm is not None:
                f.write("mtllib ./{0}\n\n".format(mat_name))
            if obj is not None:
                f.write("o {0}\n".format(obj))
            if group is not None:
                f.write("g {0}\n".format(group))
            for vertex in self:
                f.write("v {v[0]} {v[1]} {v[2]}\n".format(v=vertex))

            last_scalar = None
            for i, face in enumerate(sorted_faces + 1):
                if norm is not None:
                    if not last_scalar == sorted_scalars[i]:
                        last_scalar = sorted_scalars[i]
                        f.write("usemtl mtl_{0}\n".format(last_scalar))
                f.write("f {f[0]} {f[1]} {f[2]}\n".format(f=face))

    @classmethod
    def _load_obj(cls, path, *group_names):
        """
        Factory method
        Given a path to a obj/wavefront file, construct the object
        """
        import csv
        log = logging.getLogger()

        with open(path, mode='r') as f:
            reader = csv.reader(f, delimiter=' ')
            groups = []
            group = None
            vertex_no = 1
            for line in reader:
                if not line:
                    continue
                if line[0] == '#':
                    continue
                if line[0] == 'g':
                    if group:
                        groups.append(group)
                    group = dict(name=line[1], vertices={}, faces=[])
                elif line[0] == 'v':
                    if not group:
                        log.warning("No group specified. I invent one myself.")
                        group = dict(name='Group', vertices={}, faces=[])
                    vertex = list(map(float, line[1:4]))
                    group['vertices'][vertex_no] = vertex
                    vertex_no += 1
                elif line[0] == 'f':
                    face = []
                    for v in line[1:]:
                        w = v.split('/')
                        face.append(int(w[0]))
                    group['faces'].append(face)

        vertices = []
        for g in groups[:]:
            vertices.extend(g['vertices'].values())

        if len(group_names) != 0:
            groups = [g for g in groups if g['name'] in group_names]

        faces = []
        for g in groups:
            faces.extend(g['faces'])
        faces = np.add(np.array(faces), -1).tolist()

        """
        Building the class from retrieved vertices and faces
        """
        if len(vertices) == 0:
            return cls([])
        faceLenghts = [len(face) for face in faces]
        for i in reversed(range(len(faceLenghts))):
            length = faceLenghts[i]
            if length == 3:
                continue
            if length == 4:
                log.warning("Given a Rectangle. I will split it but "
                            "sometimes the order is different.")
                faces.insert(i + 1, faces[i][2:] + faces[i][:1])
                faces[i] = faces[i][:3]
            else:
                raise NotImplementedError()
        mesh = cls(vertices, faces=faces)
        if group_names:
            mesh = mesh.cleaned()
        return mesh

386
387
    @classmethod
    def plane(cls, *base_vectors, **kwargs):
388
389
390
391
392
393
394
        """
        Alternative constructor for creating a plane from
        Args:
            *base_vectors: see grid constructors in core. One base_vector has to
                be one-dimensional
            **kwargs: forwarded to __new__
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
395
        vertices = tfields.Tensors.grid(*base_vectors, **kwargs)
396
397
398

        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)
399
400
401
402
403
404
405
        fix_coord = None
        for coord in range(3):
            if len(base_vectors[coord]) > 1:
                continue
            if len(base_vectors[coord]) == 0:
                continue
            fix_coord = coord
406
407
        if fix_coord is None:
            raise ValueError("Describe a plane with one variable fiexed")
408
409
410
411
412
413
414
415
416
417
418
419
420
421

        var_coords = list(range(3))
        var_coords.pop(var_coords.index(fix_coord))

        faces = []
        base0, base1 = base_vectors[var_coords[0]], base_vectors[var_coords[1]]
        for i1 in range(len(base1) - 1):
            for i0 in range(len(base0) - 1):
                idx_top_left = len(base1) * (i0 + 0) + (i1 + 0)
                idx_top_right = len(base1) * (i0 + 0) + (i1 + 1)
                idx_bot_left = len(base1) * (i0 + 1) + (i1 + 0)
                idx_bot_right = len(base1) * (i0 + 1) + (i1 + 1)
                faces.append([idx_top_left, idx_top_right, idx_bot_left])
                faces.append([idx_top_right, idx_bot_left, idx_bot_right])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
422
        inst = cls.__new__(cls, vertices, faces=faces)
423
424
425
426
        return inst

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
427
428
        """
        Construct 'cuboid' along base_vectors
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
429
430
431
432
433
434
435
436
437
438
439
440
        Examples:
            Building symmetric geometries were never as easy:

            Approximated sphere with radius 1, translated in y by 2 units
            >>> sphere = tfields.Mesh3D.grid((1, 1, 1),
            ...                              (-np.pi, np.pi, 12),
            ...                              (-np.pi / 2, np.pi / 2, 12),
            ...                              coord_sys='spherical')
            >>> sphere.transform('cartesian')
            >>> sphere[:, 1] += 2

            Oktaeder
Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
441
442
443
444
            >>> oktaeder = tfields.Mesh3D.grid((1, 1, 1),
            ...                                (-np.pi, np.pi, 5),
            ...                                (-np.pi / 2, np.pi / 2, 3),
            ...                                coord_sys='spherical')
445
            >>> oktaeder.transform('cartesian')
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
446
447
448
449
450
451
452
453
454
455
456

            Cube with edge length of 2 units
            >>> cube = tfields.Mesh3D.grid((-1, 1, 2),
            ...                            (-1, 1, 2),
            ...                            (-5, -3, 2))

            Cylinder 
            >>> cylinder = tfields.Mesh3D.grid((1, 1, 1),
            ...                                (-np.pi, np.pi, 12),
            ...                                (-5, 3, 12),
            ...                                coord_sys='cylinder')
457
            >>> cylinder.transform('cartesian')
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
458

459
        """
460
461
462
        if not len(base_vectors) == 3:
            raise AttributeError("3 base_vectors vectors required")

463
464
465
        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)

466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
        indices = [0, -1]
        coords = range(3)
        baseLengthsAbove1 = [len(b) > 1 for b in base_vectors]
        # if one plane is given: rearrange indices and coords
        if not all(baseLengthsAbove1):
            indices = [0]
            for i, b in enumerate(baseLengthsAbove1):
                if not b:
                    coords = [i]
                    break

        base_vectors = list(base_vectors)
        planes = []
        for ind in indices:
            for coord in coords:
                basePart = base_vectors[:]
                basePart[coord] = np.array([base_vectors[coord][ind]],
                                           dtype=float)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
484
                planes.append(cls.plane(*basePart, **kwargs))
485
486
487
        inst = cls.merged(*planes, **kwargs)
        return inst

488
489
490
491
492
493
494
495
496
497
498
499
500
    @property
    def faces(self):
        return maps_to_faces(self.maps)

    @faces.setter
    def faces(self, faces):
        self.maps = faces_to_maps(faces)

    @property
    def faceScalars(self):
        return fields_to_scalars(self.maps[0].fields)

    @faceScalars.setter
501
502
    def faceScalars(self, scalars):
        self.maps[0].fields = scalars_to_fields(scalars)
503

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
504
505
    @cached_property()
    def _triangles(self):
506
507
508
509
510
511
        """
        with the decorator, this should be handled like an attribute though it is a function

        """
        if self.faces.size == 0:
            return tfields.Triangles3D([])
512
513
514
515
        tris = tfields.Tensors.merged(*[self[mp.flatten()] for mp in self.maps])
        map_fields = [mp.fields for mp in self.maps]
        fields = [tfields.Tensors.merged(*fields) for fields in zip(*map_fields)]
        return tfields.Triangles3D(tris, *fields)
516

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
517
    def triangles(self):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
518
519
520
521
522
523
524
525
        """
        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()

        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
526
527
528
529
530
531
532
        return self._triangles

    def centroids(self):
        return self.triangles().centroids()

    @cached_property()
    def _planes(self):
533
534
        if self.faces.size == 0:
            return tfields.Planes3D([])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
535
536
537
538
        return tfields.Planes3D(self.centroids(), self.triangles().norms())

    def planes(self):
        return self._planes
539

540
    def nfaces(self):
541
542
        return self.faces.shape[0]

543
    def in_faces(self, points, delta, assign_multiple=False):
544
545
        """
        Check whether points lie within triangles with Barycentric Technique
546
        see Triangles3D.in_triangles
547
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
548
        masks = self.triangles().in_triangles(points, delta,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
549
                                              assign_multiple=assign_multiple)
550
551
        return masks

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
552
    def removeFaces(self, face_delete_mask):
553
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
554
        Remove faces where face_delete_mask is True
555
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
556
557
558
        face_delete_mask = np.array(face_delete_mask, dtype=bool)
        self.faces = self.faces[~face_delete_mask]
        self.faceScalars = self.faceScalars[~face_delete_mask]
559

560
    def template(self, sub_mesh):
561
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
562
        'Manual' way to build a template that can be used with self.cut
563
        Returns:
564
565
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
566
        Examples:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
567
568
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
569
            >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
570
            ...                     maps=[mp])
571
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
572
573
574
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
575

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
576
577
        TODO:
            fields template not yet implemented
578
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
579
580
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
581
        mask = self.in_faces(cents, delta=None)
582
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
583
584
585
586
587
588
589
590
591
592
593
        if inst.maps:
            scalars = []
            for face_mask in mask:
                scalars.append(face_indices[face_mask][0])
            inst.maps[0].fields = [tfields.Tensors(scalars, dim=1)]
        else:
            inst.maps = [tfields.TensorFields([],
                                              tfields.Tensors([], dim=1),
                                              dim=3,
                                              dtype=int)
                        ]
594
595
        return inst

596
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
597
        """
598
        Partition the mesh with the cuts given and return the template
599
        """
600
601
602
603
604
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

605
606
607
608
609
610
611
612
613
614
615
        inst = self.copy()

        '''
        add the indices of the vertices and maps to the fields. They will be
        removed afterwards
        '''
        if not _in_recursion:
            inst.fields.append(tfields.Tensors(np.arange(len(inst))))
            for mp in inst.maps:
                mp.fields.append(tfields.Tensors(np.arange(len(mp))))

616
        # mask for points that do not fulfill the cut expression
617
        mask = inst.evalf(expression)
618
        # remove the points
619

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
620
        if not any(~mask):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
621
622
623
624
            # no vertex is valid
            inst = inst[mask]
        elif all(~mask):
            # all vertices are valid
625
            inst = inst[mask]
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
        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()
653
654
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
655
            add vertices and faces that are at the border of the cuts
656
            '''
657
            expression_parts = tfields.lib.symbolics.split_expression(expression)
658
            if len(expression_parts) > 1:
659
                new_mesh = inst.copy()
660
661
662
663
664
665
666
667
                if at_intersection == 'splitRough':
                    """
                    the following is, to speed up the process. Problem is, that
                    triangles can exist, where all points lie outside the cut,
                    but part of the area
                    still overlaps with the cut.
                    These are at the intersection line between two cuts.
                    """
668
669
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
670
671
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
672
                        if face_on_edge:
673
                            faceIntersMask[i] = True
674
                    new_mesh.removeFaces(-faceIntersMask)
675

676
                for exprPart in expression_parts:
677
678
679
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
680
            elif len(expression_parts) == 1:
681
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
682
683
684
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
685
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
686
687
688
                norm_sympy = np.array(plane_sympy.normal_vector).astype(float)
                d = -norm_sympy.dot(np.array(plane_sympy.p1).astype(float))
                plane = {'normal': norm_sympy, 'd': d}
689

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
690
                norm_vectors = inst.triangles().norms()
691
692
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
693
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
694
                                              coord_sys=field.coord_sys)
695
696
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
697
                new_norm_vectors = []
698
                newScalarMap = []
699
                n_new = 0
700

701
702
703
704
705
706
707
                vertices = np.array(inst)
                faces = np.array(inst.maps[0])
                fields = [np.array(field) for field in inst.fields]
                faces_fields = [np.array(field) for field in inst.maps[0].fields]

                face_delete_indices = set([])
                for i, face in enumerate(inst.maps[0]):
708
                    """
709
                    vertices_rejected is a mask for each face that is True, where
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
710
                    a point is on the rejected side of the plane
711
                    """
712
                    vertices_rejected = [~mask[f] for f in face]
713
714
715
716
717
                    if any(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                    if any(vertices_rejected) and not all(vertices_rejected):
                        # face on edge
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
718
719
                        n_true = vertices_rejected.count(True)
                        lonely_bool = True if n_true == 1 else False
720

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
721
                        triangle_points = [vertices[f] for f in face]
722
                        """
723
                        Add the intersection points and faces
724
                        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
725
                        intersection = _intersect(triangle_points, plane, vertices_rejected)
726
                        last_idx = len(vertices) - 1
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
727
                        for tri_list in intersection:
728
                            new_face = []
729
730
731
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
732
                                    new_face.append(face[item])
733
734
735
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
736
                                    new_face.append(last_idx + int(item.imag))
737
738
                                else:
                                    # new vertex
739
                                    new_face.append(len(vertices))
740
                                    vertices = np.append(vertices,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
741
                                                         [[float(x) for x in item]],
742
743
744
745
746
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
747
                            faces = np.append(faces, [new_face], axis=0)
748
                            faces_fields = [np.append(field,
749
                                                      [field[i]],
750
751
752
753
754
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
755
                                                dtype=int,
756
                                                coord_sys=inst.maps[0].coord_sys)
757
758
759
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
760
                                      coord_sys=inst.coord_sys)
761
762
763
764
765
                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]
766
            else:
767
                raise ValueError("Sympy expression is not splitable.")
768
            inst = inst.cleaned()
769
        elif at_intersection == 'remove':
770
            inst = inst[mask]
771
        else:
772
773
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
774
775
776
777
778
779
780
781
782
783
784
785
786

        if _in_recursion:
            template = None
        else:
            template_field = inst.fields.pop(-1)
            template_maps = []
            for mp in inst.maps:
                t_mp = tfields.TensorFields(tfields.Tensors(mp),
                                            mp.fields.pop(-1))
                template_maps.append(t_mp)
            template = tfields.Mesh3D(tfields.Tensors(inst),
                                      template_field,
                                      maps=template_maps)
787
        return inst, template
788
789

    def _cut_template(self, template):
790
791
792
793
794
795
796
797
        """
        Args:
            template (tfields.Mesh3D)

        Examples:
            >>> import tfields
            >>> import numpy as np

798
            Build mesh
799
800
801
            >>> mmap = tfields.TensorFields([[0, 1, 2], [0, 3, 4]],
            ...                             [[42, 21], [-42, -21]])
            >>> m = tfields.Mesh3D([[0]*3, [1]*3, [2]*3, [3]*3, [4]*3],
802
803
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
804
805
            ...                    maps=[mmap])

806
            Build template
807
808
809
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
810
            ...                    [1, 0, 3, 2, 4],
811
812
            ...                    maps=[tmap])

813
            Use template as instruction to make a fast cut
814
815
816
817
818
819
820
821
822
            >>> res = m._cut_template(t)
            >>> assert np.array_equal(res.fields,
            ...                       [[0.1, 0.0, 0.3, 0.2, 0.4],
            ...                        [-0.1, 0.0, -0.3, -0.2, -0.4]])

            >>> assert np.array_equal(res.maps[0].fields[0],
            ...                       [[-42, -21], [42, 21]])
                                   
        """
823
824
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

825
        # Redirect fields
826
        fields = []
827
        if template.fields:
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
            template_field = np.array(template.fields[0])
            if len(self) > 0:
                '''
                if new vertices have been created in the template, it is
                in principle unclear what fields we have to refer to.
                Thus in creating the template, we gave np.nan.
                To make it fast, we replace nan with 0 as a dummy and correct
                the field entries afterwards with np.nan.
                '''
                nan_mask = np.isnan(template_field)
                template_field[nan_mask] = 0  # dummy reference to index 0.
                template_field = template_field.astype(int)
                for field in self.fields:
                    projected_field = field[template_field]
                    projected_field[nan_mask] = np.nan  # correction for nan
                    fields.append(projected_field)
844

845
        # Redirect maps and their fields
846
847
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
848
            mp_fields = []
849
850
851
852
            for field in mp.fields:
                if len(template_mp) == 0 and len(template_mp.fields) == 0:
                    mp_fields.append(field[0:0])  # np.empty
                else:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
853
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
854
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
855
856
857
                                          *mp_fields)
            maps.append(new_mp)

858
859
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
860
                              maps=maps)
861
862
        return inst

863
    def cut(self, expression, coord_sys=None, at_intersection=None,
864
            return_template=False):
865
866
867
        """
        cut method for Mesh3D.
        Args:
868
869
870
871
872
873
874
875
876
877
            expression (sympy logical expression | Mesh3D):
                sympy locical expression: Sympy expression that defines planes
                    in 3D
                Mesh3D: A mesh3D will be interpreted as a template, i.e. a
                    fast instruction of how to cut the triangles.
                    It is the second part of the tuple, returned by a previous
                    cut with a sympy locial expression with 'return_template=True'.
                    We use the vertices and maps of the Mesh as the sceleton of
                    the returned mesh. The fields are mapped according to
                    indices in the template.maps[i].fields.
878
            coord_sys (coordinate system to cut in):
879
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
880
881
882
                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.
883
884
            return_template (bool): If True: return the template
                            to redo the same cut fast
885
886
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
887
            >>> import numpy as np
888
            >>> import tfields
889
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
890
            >>> cut_expr = x > 1.5
891

892
893
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
894
            ...                         (0, 0, 1))
895
896
897
898
899
900
            >>> m.fields.append(tfields.Tensors(np.linspace(0, len(m) - 1,
            ...                                             len(m))))
            >>> m.maps[0].fields.append(
            ...     tfields.Tensors(np.linspace(0,
            ...                                 len(m.maps[0]) - 1,
            ...                                 len(m.maps[0]))))
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
901
            >>> mNew = m.cut(cut_expr)
902
            >>> len(mNew)
903
            8
904
            >>> mNew.nfaces()
905
906
907
908
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
909
910
911
912
913
914
915
916
917
918
            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')
919
            >>> float(m_split[:, 0].min())
920
            1.5
921
            >>> len(m_split)
922
            15
923
            >>> m_split.nfaces()
924
925
            15

926
927
            Cut with 'return_template=True' will return the exact same mesh but
            additionally an instruction to conduct the exact same cut fast (template)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
928
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
929
            ...                                    return_template=True)
930
931
932
933
934
935
936
937
938
939
            >>> m_split_template = m.cut(template)
            >>> assert m_split.equal(m_split_2, equal_nan=True)
            >>> assert m_split.equal(m_split_template, equal_nan=True)
            >>> assert len(template.fields) == 1
            >>> assert len(m_split.fields) == 1
            >>> assert len(m_split_template.fields) == 1
            >>> assert m_split.fields[0].equal(
            ...     list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
            >>> assert m_split_template.fields[0].equal(
            ...     list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
940

941
            This seems irrelevant at first but consider, the map field or the
942
943
944
945
946
947
948
949
            tensor field changes:
            >>> m_altered_fields = m.copy()
            >>> m_altered_fields[0] += 42
            >>> assert not m_split.equal(m_altered_fields.cut(template))
            >>> assert tfields.Tensors(m_split).equal(m_altered_fields.cut(template))
            >>> assert tfields.Tensors(m_split.maps[0]).equal(m_altered_fields.cut(template).maps[0])

            The cut expression may be a sympy.BooleanFunction:
950
951
            >>> cut_expr_bool_fun = (x > 1.5) & (y < 1.5) & (y >0.2) & (z > -0.5)
            >>> m_split_bool = m.cut(cut_expr_bool_fun, at_intersection='split')
952
953
954

        Returns:
            copy of cut mesh
955
            * optional: template
956
957

        """
958
        with self.tmp_transform(coord_sys or self.coord_sys):
959
            if isinstance(expression, Mesh3D):
960
961
                template = expression
                obj = self._cut_template(template)
962
963
964
965
966
967
            else:
                at_intersection = at_intersection or "remove"
                obj, template = self._cut_sympy(expression, at_intersection=at_intersection)
        if return_template:
            return obj, template
        return obj
968

969
970
971
972
973
974
975
976
977
978
979
980
981
    def disjoint_parts(self, return_template=False):
        mp_description = self.disjoint_map(0)
        parts = self.parts(mp_description)
        if not return_template:
            return parts
        else:
            templates = []
            for i, part in enumerate(parts):
                template = part.copy()
                template.maps[0].fields[0] = tfields.Tensors(mp_description[1][i])
                templates.append(template)
            return parts, templates

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
982
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
983
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
984
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
985
986
987
988
989
990
991
992
        """
        scalars_demanded = any([v in kwargs for v in ['vmin', 'vmax', 'cmap']])
        map_index = kwargs.pop('map_index', None if not scalars_demanded else 0)
        if map_index is not None:
            if not len(self.maps[0]) == 0:
                kwargs['color'] = self.maps[0].fields[map_index]

        dim_defined = False
993
        if 'axes' in kwargs:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
994
            dim_defined = True
995
996
        if 'z_index' in kwargs:
            if kwargs['z_index'] is not None:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
997
998
999
1000
1001
1002
1003
1004
1005
1006
                kwargs['dim'] = 3
            else:
                kwargs['dim'] = 2
            dim_defined = True
        if 'dim' in kwargs:
            dim_defined = True

        if not dim_defined:
            kwargs['dim'] = 2

1007
        return rna.plotting.plot_mesh(self, self.faces, **kwargs)
1008

1009

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1010
if __name__ == '__main__':  # pragma: no cover
1011
1012
    import doctest

1013
1014
    doctest.run_docstring_examples(Mesh3D.cut, globals())
    # doctest.testmod()