mesh.py 23.8 KB
Newer Older
1
2
#!/usr/bin/env python
#
3
# mesh.py - The TriangleMesh class.
4
5
6
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
7
"""This module provides the :class:`Mesh` class, which represents a
8
3D model made of triangles.
9

10
See also the following modules:
11

12
  .. autosummary::
13

14
15
16
17
     fsl.data.vtk
     fsl.data.gifti
     fsl.data.freesurfer

18
19
A handful of standalone functions are provided in this module, for doing
various things with meshes:
20
21
22
23
24
25
26

  .. autosummary::
     :nosignatures:

     calcFaceNormals
     calcVertexNormals
     needsFixing
27
28
"""

29

30
import logging
31
32
import collections

33
34
import os.path as op
import numpy   as np
35

36
37
38
import fsl.utils.meta       as meta
import fsl.utils.notifier   as notifier
import fsl.transform.affine as affine
39

40
41
42
43

log = logging.getLogger(__name__)


44
45
class Mesh(notifier.Notifier, meta.Meta):
    """The ``Mesh`` class represents a 3D model. A mesh is defined by a
46
    collection of ``N`` vertices, and ``M`` triangles.  The triangles are
Paul McCarthy's avatar
Paul McCarthy committed
47
    defined by ``(M, 3)`` indices into the list of vertices.
48
49


50
    A ``Mesh`` instance has the following attributes:
51

52

53
    ============== ======================================================
54
    ``name``       A name, typically the file name sans-suffix.
55

56
57
    ``dataSource`` Full path to the mesh file (or ``None`` if there is
                   no file associated with this mesh).
58

59
60
    ``nvertices``  The number of vertices in the mesh.

61
62
63
64
65
    ``vertices``   A ``(n, 3)`` array containing the currently selected
                   vertices. You can assign  a vertex set key to this
                   attribute to change the selected vertex set.

    ``bounds``     The lower and upper bounds
66

67
68
    ``indices``    A ``(m, 3)`` array containing the vertex indices
                   for ``m`` triangles
69

70
71
    ``normals``    A  ``(m, 3)`` array containing face normals for the
                   triangles
72

73
74
    ``vnormals``   A ``(n, 3)`` array containing vertex normals for the
                   the current vertices.
75
76
77
78
79

    ``trimesh``    (if the `trimesh <https://github.com/mikedh/trimesh>`_
                   library is present) A ``trimesh.Trimesh`` object which
                   can be used for geometric queries on the mesh.
    ============== ======================================================
80

81

82
83
84
85
86
87
88
89
    **Vertex sets**


    A ``Mesh`` object can be associated with multiple sets of vertices, but
    only one set of triangles. Vertices can be added via the
    :meth:`addVertices` method. Each vertex set must be associated with a
    unique key - you can then select the current vertex set via the
    :meth:`vertices` property. Most ``Mesh`` methods will raise a ``KeyError``
90
91
92
93
94
95
    if you have not added any vertex sets, or selected a vertex set. The
    following methods are available for managing vertex sets:

    .. autosummary::
       :nosignatures:

96
       loadVertices
97
98
       addVertices
       selectedVertices
99
       vertexSets
100

101
102
103
104
105
106
    .. note:: Internally the ``Mesh`` class may store two versions of the
              triangles, with opposite unwinding orders. It keeps track of the
              required triangle unwinding order for each vertex set, so that
              the :meth:`indices` method will return the appropriate copy for
              the currently selected vertex set.

107
108
109
110
111
112

    **Vertex data**


    A ``Mesh`` object can store vertex-wise data. The following methods can be
    used for adding/retrieving vertex data:
113
114
115
116

    .. autosummary::
       :nosignatures:

117
       loadVertexData
118
       addVertexData
119
       getVertexData
120
       vertexDataSets
121
       clearVertexData
122

123

124
    **Notification**
125

126

127
128
129
130
131
132
    The ``Mesh`` class inherits from the :class:`Notifier` class. Whenever the
    ``Mesh`` vertex set is changed, a notification is emitted via the
    ``Notifier`` interface, with a topic of ``'vertices'``. When this occurs,
    the :meth:`vertices`, :meth:`bounds`, :meth:`normals` and :attr:`vnormals`
    properties will all change so that they return data specific to the newly
    selected vertex set.
133
134


Paul McCarthy's avatar
Paul McCarthy committed
135
    **Metadata**
136

137

138
139
    The ``Mesh`` class also inherits from the :class:`Meta` class, so
    any metadata associated with the ``Mesh`` may be added via those methods.
140

141

142
    **Geometric queries**
143

144

145
146
    If the ``trimesh`` library is present, the following methods may be used
    to perform geometric queries on a mesh:
147

148
149
    .. autosummary::
       :nosignatures:
150

151
152
153
154
       rayIntersection
       planeIntersection
       nearestVertex
    """
155
156


157
158
159
160
161
162
    def __init__(self,
                 indices,
                 name='mesh',
                 dataSource=None,
                 vertices=None,
                 fixWinding=False):
163
164
165
166
        """Create a ``Mesh`` instance.

        Before a ``Mesh`` can be used, some vertices must be added via the
        :meth:`addVertices` method.
167

168
169
170
171
172
        :arg indices:    A list of indices into the vertex data, defining the
                         mesh triangles.

        :arg name:       A name for this ``Mesh``.

173
174
175
176
        :arg dataSource: The data source for this ``Mesh``.

        :arg vertices:   Initial vertex set to add - given the key
                         ``'default'``.
177
178
179

        :arg fixWinding: Ignored if ``vertices is None``. Passed through to the
                         :meth:`addVertices` method along with ``vertices``.
180
181
        """

182
183
        self.__name       = name
        self.__dataSource = dataSource
184
        self.__nvertices  = int(indices.max()) + 1
185
186
187
188
189
190
191
        self.__selected   = None

        # We potentially store two copies of
        # the indices, with opposite unwinding
        # orders. The vindices dict stores refs
        # to one or the other for each vertex
        # set.
192
        self.__indices      = np.asarray(indices, dtype=np.int32).reshape((-1, 3))
193
194
        self.__fixedIndices = None
        self.__vindices     = collections.OrderedDict()
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

        # All of these are populated
        # in the addVertices method
        self.__vertices = collections.OrderedDict()
        self.__loBounds = collections.OrderedDict()
        self.__hiBounds = collections.OrderedDict()

        # These get populated on
        # normals/vnormals accesses
        self.__faceNormals = collections.OrderedDict()
        self.__vertNormals = collections.OrderedDict()

        # this gets populated in
        # the addVertexData method
        self.__vertexData  = collections.OrderedDict()

        # this gets populated
        # in the trimesh method
        self.__trimesh = collections.OrderedDict()

215
216
217
        # Add initial vertex
        # set if provided
        if vertices is not None:
218
            self.addVertices(vertices, fixWinding=fixWinding)
219

220

221
222
223
224
225
226
227
228
229
230
231
232
233
    def __repr__(self):
        """Returns a string representation of this ``Mesh`` instance. """
        return '{}({}, {})'.format(type(self).__name__,
                                   self.name,
                                   self.dataSource)


    def __str__(self):
        """Returns a string representation of this ``Mesh`` instance.
        """
        return self.__repr__()


234
235
236
237
238
239
    @property
    def name(self):
        """Returns the name of this ``Mesh``. """
        return self.__name


Paul McCarthy's avatar
Paul McCarthy committed
240
241
242
243
244
245
    @name.setter
    def name(self, name):
        """Set the name of this ``Mesh``. """
        self.__name = name


246
247
248
249
    @property
    def dataSource(self):
        """Returns the data source of this ``Mesh``. """
        return self.__dataSource
250
251


252
253
254
255
256
257
    @property
    def nvertices(self):
        """Returns the number of vertices in the mesh. """
        return self.__nvertices


258
259
260
    @property
    def vertices(self):
        """The ``(N, 3)`` vertices of this mesh. """
261
262
263
        return self.__vertices[self.__selected]


264
    @vertices.setter
265
266
267
    def vertices(self, key):
        """Select the current vertex set - a ``KeyError`` is raised
        if no vertex set with the specified ``key`` has been added.
268
269
270
271

        When the current vertex set is changed, a notification is emitted
        through the :class:`.Notifier` interface, with the topic
        ``'vertices'``.
272
273
274
275
276
        """

        # Force a key error if
        # the key is invalid
        self.__vertices[key]
277

278
279
280
        if self.__selected != key:
            self.__selected = key
            self.notify(topic='vertices')
281

282
283
284
285

    @property
    def indices(self):
        """The ``(M, 3)`` triangles of this mesh. """
286
        return self.__vindices[self.__selected]
287
288


289
290
291
292
293
    @property
    def normals(self):
        """A ``(M, 3)`` array containing surface normals for every
        triangle in the mesh, normalised to unit length.
        """
294

295
        selected = self.__selected
296
        indices  = self.__vindices[selected]
297
298
        vertices = self.__vertices[selected]
        fnormals = self.__faceNormals.get(selected, None)
299

300
301
302
        if fnormals is None:
            fnormals = calcFaceNormals(vertices, indices)
            self.__faceNormals[selected] = fnormals
303

304
        return fnormals
305
306
307
308


    @property
    def vnormals(self):
309
        """A ``(N, 3)`` array containing normals for every vertex
310
311
312
        in the mesh.
        """

313
        selected = self.__selected
314
        indices  = self.__vindices[selected]
315
316
        vertices = self.__vertices[selected]
        vnormals = self.__vertNormals.get(selected, None)
317

318
319
320
        if vnormals is None:
            vnormals = calcVertexNormals(vertices, indices, self.normals)
            self.__vertNormals[selected] = vnormals
321

322
        return vnormals
323
324


325
326
    @property
    def bounds(self):
327
        """Returns a tuple of values which define a minimal bounding box that
328
329
        will contain all of the currently selected vertices in this
        ``Mesh`` instance. The bounding box is arranged like so:
330
331
332

            ``((xlow, ylow, zlow), (xhigh, yhigh, zhigh))``
        """
333
334
335
        lo = self.__loBounds[self.__selected]
        hi = self.__hiBounds[self.__selected]
        return lo, hi
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
    def loadVertices(self, infile, key=None, **kwargs):
        """Loads vertex data from the given ``infile``, and adds it as a vertex
        set with the given ``key``. This implementation supports loading vertex
        data from white-space delimited text files via ``numpy.loadtxt``, but
        sub-classes may override this method to support additional file types.


        :arg infile: File to load data from.

        :arg key:    Key to pass to :meth:`addVertices`. If not provided,
                     set to ``infile`` (converted to an absolute path)

        All of the other arguments are passed through to :meth:`addVertices`.

        :returns:    The loaded vertices.
        """

        infile = op.abspath(infile)

        if key is None:
            key = infile

        vertices = np.loadtxt(infile)

362
        return self.addVertices(vertices, key, **kwargs)
363
364


365
    def addVertices(self, vertices, key=None, select=True, fixWinding=False):
366
        """Adds a set of vertices to this ``Mesh``.
367

368
369
        :arg vertices:   A `(n, 3)` array containing ``n`` vertices, compatible
                         with the indices specified in :meth:`__init__`.
370

371
372
        :arg key:        A key for this vertex set. If ``None`` defaults to
                         ``'default'``.
373
374
375
376
377
378
379

        :arg select:     If ``True`` (the default), this vertex set is
                         made the currently selected vertex set.

        :arg fixWinding: Defaults to ``False``. If ``True``, the vertex
                         winding order of every triangle is is fixed so they
                         all have outward-facing normal vectors.
380
381
382
383
384

        :returns:        The vertices, possibly reshaped

        :raises:         ``ValueError`` if the provided ``vertices`` array
                         has the wrong number of vertices.
385
386
        """

387
388
389
        if key is None:
            key = 'default'

390
391
392
        vertices = np.asarray(vertices)
        lo       = vertices.min(axis=0)
        hi       = vertices.max(axis=0)
393

394
395
396
397
398
399
400
401
402
403
404
405
406
        # Don't allow vertices of
        # different size to be added
        try:
            vertices = vertices.reshape(self.nvertices, 3)

        # reshape raised an error -
        # wrong number of vertices
        except ValueError:
            raise ValueError('{}: invalid number of vertices: '
                             '{} != ({}, 3)'.format(key,
                                                    vertices.shape,
                                                    self.nvertices))

407
        self.__vertices[key] = vertices
408
        self.__vindices[key] = self.__indices
409
410
        self.__loBounds[key] = lo
        self.__hiBounds[key] = hi
411

412
413
        if select:
            self.vertices = key
414

415
416
        if fixWinding:
            indices  = self.__indices
417
            normals  = calcFaceNormals(vertices, indices)
418
            needsFix = needsFixing(vertices, indices, normals, lo, hi)
419

420
421
422
            # See needsFixing documentation
            if needsFix:

423
424
                if self.__fixedIndices is None:
                    self.__fixedIndices = indices[:, [0, 2, 1]]
425

426
427
                self.__vindices[   key] = self.__fixedIndices
                self.__faceNormals[key] = normals * -1
428
429
            else:
                self.__faceNormals[key] = normals
430

431
432
        return vertices

433

434
    def vertexSets(self):
435
436
437
        """Returns a list containing the keys of all vertex sets. """
        return list(self.__vertices.keys())

438
439
440
441
442
443

    def selectedVertices(self):
        """Returns the key of the currently selected vertex set. """
        return self.__selected


444
445
446
447
448
449
450
451
452
453
    def loadVertexData(self, infile, key=None):
        """Loads vertex-wise data from the given ``infile``, and adds it with
        the given ``key``. This implementation supports loading data from
        whitespace-delimited text files via ``numpy.loadtxt``, but sub-classes
        may override this method to support additional file types.

        :arg infile: File to load data from.

        :arg key:    Key to pass to :meth:`addVertexData`. If not provided,
                     set to ``infile`` (converted to an absolute path)
454
455

        :returns:    The loaded vertex data.
456
457
458
459
460
461
462
463
464
        """

        infile = op.abspath(infile)

        if key is None:
            key = infile

        vertexData = np.loadtxt(infile)

465
        return self.addVertexData(key, vertexData)
466
467


468
469
470
    def addVertexData(self, key, vdata):
        """Adds a vertex-wise data set to the ``Mesh``. It can be retrieved
        by passing the specified ``key`` to the :meth:`getVertexData` method.
471
472

        :returns: The vertex data, possibly reshaped.
473
        """
474

475
        nvertices = self.nvertices
476

477
478
479
        if vdata.ndim not in (1, 2) or vdata.shape[0] != nvertices:
            raise ValueError('{}: incompatible vertex data '
                             'shape: {}'.format(key, vdata.shape))
480

481
482
483
484
        vdata                  = vdata.reshape(nvertices, -1)
        self.__vertexData[key] = vdata

        return vdata
485
486


487
488
489
490
    def getVertexData(self, key):
        """Returns the vertex data for the given ``key`` from the
        internal vertex data cache. If there is no vertex data iwth the
        given key, a ``KeyError`` is raised.
491
492
        """

493
        return self.__vertexData[key]
494
495
496
497


    def clearVertexData(self):
        """Clears the internal vertex data cache - see the
498
        :meth:`addVertexData` and :meth:`getVertexData` methods.
499
        """
500
        self.__vertexData = collections.OrderedDict()
501
502


503
504
505
506
507
    def vertexDataSets(self):
        """Returns a list of keys for all loaded vertex data sets. """
        return list(self.__vertexData.keys())


508
    @property
509
510
511
512
513
    def trimesh(self):
        """Reference to a ``trimesh.Trimesh`` object which can be used for
        geometric operations on the mesh.

        If the ``trimesh`` or ``rtree`` libraries are not available, this
514
515
        function returns ``None``, and none of the geometric query methods
        will do anything.
516
517
518
519
520
521
522
523
524
525
526
527
528
529
        """

        # trimesh is an optional dependency - rtree
        # is a depedendency of trimesh which is a
        # wrapper around libspatialindex, without
        # which trimesh can't be used for calculating
        # ray-mesh intersections.
        try:
            import trimesh
            import rtree   # noqa
        except ImportError:
            log.warning('trimesh is not available')
            return None

530
531
532
533
534
535
536
        tm = self.__trimesh.get(self.__selected, None)

        if tm is None:
            tm = trimesh.Trimesh(self.vertices,
                                 self.indices,
                                 process=False,
                                 validate=False)
537

538
            self.__trimesh[self.__selected] = tm
539

540
        return tm
541
542


543
    def rayIntersection(self, origins, directions, vertices=False):
544
545
546
547
548
        """Calculate the intersection between the mesh, and the rays defined by
        ``origins`` and ``directions``.

        :arg origins:    Sequence of ray origins
        :arg directions: Sequence of ray directions
549
550
551
552
553
554
555
556
557
        :returns:        A tuple containing:

                           - A ``(n, 3)`` array containing the coordinates
                             where the mesh was intersected by each of the
                             ``n`` rays.

                           - A ``(n,)`` array containing the indices of the
                             triangles that were intersected by each of the
                             ``n`` rays.
558
559
        """

560
        trimesh = self.trimesh
561
562

        if trimesh is None:
563
            return np.zeros((0, 3)), np.zeros((0,))
564

565
        tris, rays, locs = trimesh.ray.intersects_id(
566
567
568
            origins,
            directions,
            return_locations=True,
569
            multiple_hits=False)
570

Paul McCarthy's avatar
Paul McCarthy committed
571
        if len(tris) == 0:
572
573
            return np.zeros((0, 3)), np.zeros((0,))

574
575
        # sort by ray. I'm Not sure if this is
        # needed - does trimesh do it for us?
576
        rayIdxs = np.asarray(np.argsort(rays))
577
578
        locs    = locs[rayIdxs]
        tris    = tris[rayIdxs]
579

580
        return locs, tris
581
582


583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
    def nearestVertex(self, points):
        """Identifies the nearest vertex to each of the provided points.

        :arg points: A ``(n, 3)`` array containing the points to query.

        :returns:    A tuple containing:

                      - A ``(n, 3)`` array containing the nearest vertex for
                        for each of the ``n`` input points.

                      - A ``(n,)`` array containing the indices of each vertex.

                      - A ``(n,)`` array containing the distance from each
                        point to the nearest vertex.
        """

599
        trimesh = self.trimesh
600
601

        if trimesh is None:
602
            return np.zeros((0, 3)), np.zeros((0, )), np.zeros((0, ))
603
604
605
606
607
608
609

        dists, idxs = trimesh.nearest.vertex(points)
        verts       = self.vertices[idxs, :]

        return verts, idxs, dists


610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
    def planeIntersection(self,
                          normal,
                          origin,
                          distances=False):
        """Calculate the intersection of this ``TriangleMesh`` with
        the plane defined by ``normal`` and ``origin``.

        :arg normal:    Vector defining the plane orientation

        :arg origin:    Point defining the plane location

        :arg distances: If ``True``, barycentric coordinates for each
                        intersection line vertex are calculated and returned,
                        giving their respective distance from the intersected
                        triangle vertices.

        :returns:       A tuple containing
Paul McCarthy's avatar
Paul McCarthy committed
627

628
629
630
631
632
633
                          - A ``(m, 2, 3)`` array containing ``m`` vertices:
                            of a set of lines, defining the plane intersection

                          - A ``(m,)`` array containing the indices of the
                            ``m`` triangles that were intersected.

Paul McCarthy's avatar
Paul McCarthy committed
634
                          - (if ``distances is True``) A ``(m, 2, 3)`` array
635
636
637
638
639
                            containing the barycentric coordinates of each
                            line vertex with respect to its intersected
                            triangle.
        """

640
        trimesh = self.trimesh
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672

        if trimesh is None:
            return np.zeros((0, 3)), np.zeros((0, 3))

        import trimesh.intersections as tmint
        import trimesh.triangles     as tmtri

        lines, faces = tmint.mesh_plane(
            trimesh,
            plane_normal=normal,
            plane_origin=origin,
            return_faces=True)

        if not distances:
            return lines, faces

        # Calculate the barycentric coordinates
        # (distance from triangle vertices) for
        # each intersection line

        triangles = self.vertices[self.indices[faces]].repeat(2, axis=0)
        points    = lines.reshape((-1, 3))

        if triangles.size > 0:
            dists = tmtri.points_to_barycentric(triangles, points)
            dists = dists.reshape((-1, 2, 3))
        else:
            dists = np.zeros((0, 2, 3))

        return lines, faces, dists


673
674
675
def calcFaceNormals(vertices, indices):
    """Calculates face normals for the mesh described by ``vertices`` and
    ``indices``.
676

677
678
679
680
    :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
    :arg indices:  A ``(m, 3)`` array containing the mesh triangles.
    :returns:      A ``(m, 3)`` array containing normals for every triangle in
                   the mesh.
681
    """
682

683
684
685
    v0 = vertices[indices[:, 0]]
    v1 = vertices[indices[:, 1]]
    v2 = vertices[indices[:, 2]]
686

687
    fnormals = np.cross((v1 - v0), (v2 - v0))
688
    fnormals = affine.normalise(fnormals)
689

690
    return np.atleast_2d(fnormals)
691

692

693
694
695
def calcVertexNormals(vertices, indices, fnormals):
    """Calculates vertex normals for the mesh described by ``vertices``
    and ``indices``.
696

697
698
699
700
701
702
    :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
    :arg indices:  A ``(m, 3)`` array containing the mesh triangles.
    :arg fnormals: A ``(m, 3)`` array containing the face/triangle normals.
    :returns:      A ``(n, 3)`` array containing normals for every vertex in
                   the mesh.
    """
703

704
    vnormals = np.zeros((vertices.shape[0], 3), dtype=float)
705

706
707
708
709
710
    # TODO make fast. I can't figure
    # out how to use np.add.at to
    # accumulate the face normals for
    # each vertex.
    for i in range(indices.shape[0]):
711

712
        v0, v1, v2 = indices[i]
713

714
715
716
        vnormals[v0, :] += fnormals[i]
        vnormals[v1, :] += fnormals[i]
        vnormals[v2, :] += fnormals[i]
717

718
    # normalise to unit length
719
    return affine.normalise(vnormals)
720
721


722
723
724
def needsFixing(vertices, indices, fnormals, loBounds, hiBounds):
    """Determines whether the triangle winding order, for the mesh described by
    ``vertices`` and ``indices``, needs to be flipped.
725

726
727
728
    If this function returns ``True``, the given ``indices`` and ``fnormals``
    need to be adjusted so that all face normals are facing outwards from the
    centre of the mesh. The necessary adjustments are as follows::
729

730
731
        indices[:, [1, 2]] = indices[:, [2, 1]]
        fnormals           = fnormals * -1
732

733
734
735
736
737
    :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
    :arg indices:  A ``(m, 3)`` array containing the mesh triangles.
    :arg fnormals: A ``(m, 3)`` array containing the face/triangle normals.
    :arg loBounds: A ``(3, )`` array contaning the low vertex bounds.
    :arg hiBounds: A ``(3, )`` array contaning the high vertex bounds.
738

739
740
    :returns:      ``True`` if the ``indices`` and ``fnormals`` need to be
                   adjusted, ``False`` otherwise.
741
742
    """

743
744
745
746
747
748
749
750
751
752
    # Define a viewpoint which is
    # far away from the mesh.
    camera = loBounds - (hiBounds - loBounds)

    # Find the nearest vertex
    # to the viewpoint
    dists = np.sqrt(np.sum((vertices - camera) ** 2, axis=1))
    ivert = np.argmin(dists)
    vert  = vertices[ivert]

753
754
755
756
757
    # Get all the triangles
    # that this vertex is in
    # and their face normals
    itris = np.where(indices == ivert)[0]
    norms = fnormals[itris, :]
758

759
    # Calculate the angle between each
760
    # normal, and a vector from the
761
762
763
764
765
766
767
    # vertex to the camera. If more than
    # 50% of the angles are negative
    # (== more than 90 degrees == the
    # face is facing away from the
    # camera), assume that we need to
    # flip the triangle winding order.
    angles = np.dot(norms, affine.normalise(camera - vert))
768
    return ((angles >= 0).sum() / len(itris)) < 0.5