vedo.tetmesh

  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3
  4try:
  5    import vedo.vtkclasses as vtk
  6except ImportError:
  7    import vtkmodules.all as vtk
  8
  9import numpy as np
 10import vedo
 11from vedo import utils
 12from vedo.base import BaseGrid
 13from vedo.mesh import Mesh
 14from vedo.file_io import download, loadUnStructuredGrid
 15
 16
 17__docformat__ = "google"
 18
 19__doc__ = """
 20Work with tetrahedral meshes.
 21
 22![](https://vedo.embl.es/images/volumetric/82767107-2631d500-9e25-11ea-967c-42558f98f721.jpg)
 23"""
 24
 25__all__ = ["TetMesh", "delaunay3d"]
 26
 27
 28##########################################################################
 29def delaunay3d(mesh, radius=0, tol=None):
 30    """
 31    Create 3D Delaunay triangulation of input points.
 32
 33    Arguments:
 34
 35        radius : (float)
 36            specify distance (or "alpha") value to control output.
 37            For a non-zero values, only tetra contained within the circumsphere
 38            will be output.
 39
 40        tol : (float)
 41            Specify a tolerance to control discarding of closely spaced points.
 42            This tolerance is specified as a fraction of the diagonal length of
 43            the bounding box of the points.
 44    """
 45    deln = vtk.vtkDelaunay3D()
 46    if utils.is_sequence(mesh):
 47        pd = vtk.vtkPolyData()
 48        vpts = vtk.vtkPoints()
 49        vpts.SetData(utils.numpy2vtk(mesh, dtype=np.float32))
 50        pd.SetPoints(vpts)
 51        deln.SetInputData(pd)
 52    else:
 53        deln.SetInputData(mesh.GetMapper().GetInput())
 54    deln.SetAlpha(radius)
 55    deln.AlphaTetsOn()
 56    deln.AlphaTrisOff()
 57    deln.AlphaLinesOff()
 58    deln.AlphaVertsOff()
 59    deln.BoundingTriangulationOff()
 60    if tol:
 61        deln.SetTolerance(tol)
 62    deln.Update()
 63    m = TetMesh(deln.GetOutput())
 64    m.pipeline = utils.OperationNode(
 65        "delaunay3d", c="#e9c46a:#edabab", parents=[mesh],
 66    )
 67    return m
 68
 69
 70def _buildtetugrid(points, cells):
 71    ug = vtk.vtkUnstructuredGrid()
 72
 73    if len(points) == 0:
 74        return ug
 75    if not utils.is_sequence(points[0]):
 76        return ug
 77
 78    if len(cells) == 0:
 79        return ug
 80
 81    if not utils.is_sequence(cells[0]):
 82        tets = []
 83        nf = cells[0] + 1
 84        for i, cl in enumerate(cells):
 85            if i in (nf, 0):
 86                k = i + 1
 87                nf = cl + k
 88                cell = [cells[j + k] for j in range(cl)]
 89                tets.append(cell)
 90        cells = tets
 91
 92    source_points = vtk.vtkPoints()
 93    varr = utils.numpy2vtk(points, dtype=np.float32)
 94    source_points.SetData(varr)
 95    ug.SetPoints(source_points)
 96
 97    source_tets = vtk.vtkCellArray()
 98    for f in cells:
 99        ele = vtk.vtkTetra()
100        pid = ele.GetPointIds()
101        for i, fi in enumerate(f):
102            pid.SetId(i, fi)
103        source_tets.InsertNextCell(ele)
104    ug.SetCells(vtk.VTK_TETRA, source_tets)
105    return ug
106
107
108##########################################################################
109class TetMesh(BaseGrid, vtk.vtkVolume):
110    """The class describing tetrahedral meshes."""
111
112    def __init__(
113        self,
114        inputobj=None,
115        c=("r", "y", "lg", "lb", "b"),  # ('b','lb','lg','y','r')
116        alpha=(0.5, 1),
117        alpha_unit=1,
118        mapper="tetra",
119    ):
120        """
121        Arguments:
122            inputobj : (vtkDataSet, list, str)
123                list of points and tet indices, or filename
124            alpha_unit : (float)
125                opacity scale
126            mapper : (str)
127                choose a visualization style in `['tetra', 'raycast', 'zsweep']`
128        """
129        vtk.vtkVolume.__init__(self)
130        BaseGrid.__init__(self)
131
132        self.useArray = 0
133
134        # inputtype = str(type(inputobj))
135        # print('TetMesh inputtype', inputtype)
136
137        ###################
138        if inputobj is None:
139            self._data = vtk.vtkUnstructuredGrid()
140
141        elif isinstance(inputobj, vtk.vtkUnstructuredGrid):
142            self._data = inputobj
143
144        elif isinstance(inputobj, vtk.vtkRectilinearGrid):
145            r2t = vtk.vtkRectilinearGridToTetrahedra()
146            r2t.SetInputData(inputobj)
147            r2t.RememberVoxelIdOn()
148            r2t.SetTetraPerCellTo6()
149            r2t.Update()
150            self._data = r2t.GetOutput()
151
152        elif isinstance(inputobj, vtk.vtkDataSet):
153            r2t = vtk.vtkDataSetTriangleFilter()
154            r2t.SetInputData(inputobj)
155            # r2t.TetrahedraOnlyOn()
156            r2t.Update()
157            self._data = r2t.GetOutput()
158
159        elif isinstance(inputobj, str):
160            if "https://" in inputobj:
161                inputobj = download(inputobj, verbose=False)
162            ug = loadUnStructuredGrid(inputobj)
163            tt = vtk.vtkDataSetTriangleFilter()
164            tt.SetInputData(ug)
165            tt.SetTetrahedraOnly(True)
166            tt.Update()
167            self._data = tt.GetOutput()
168
169        elif utils.is_sequence(inputobj):
170            # if "ndarray" not in inputtype:
171            #     inputobj = np.array(inputobj)
172            self._data = _buildtetugrid(inputobj[0], inputobj[1])
173
174        ###################
175        if "tetra" in mapper:
176            self._mapper = vtk.vtkProjectedTetrahedraMapper()
177        elif "ray" in mapper:
178            self._mapper = vtk.vtkUnstructuredGridVolumeRayCastMapper()
179        elif "zs" in mapper:
180            self._mapper = vtk.vtkUnstructuredGridVolumeZSweepMapper()
181        elif isinstance(mapper, vtk.vtkMapper):
182            self._mapper = mapper
183        else:
184            vedo.logger.error(f"Unknown mapper type {type(mapper)}")
185            raise RuntimeError()
186
187        self._mapper.SetInputData(self._data)
188        self.SetMapper(self._mapper)
189        self.color(c).alpha(alpha)
190        if alpha_unit:
191            self.GetProperty().SetScalarOpacityUnitDistance(alpha_unit)
192
193        # remember stuff:
194        self._color = c
195        self._alpha = alpha
196        self._alpha_unit = alpha_unit
197
198        self.pipeline = utils.OperationNode(
199            self, comment=f"#tets {self._data.GetNumberOfCells()}",
200            c="#9e2a2b",
201        )
202        # -----------------------------------------------------------
203
204    def _repr_html_(self):
205        """
206        HTML representation of the TetMesh object for Jupyter Notebooks.
207
208        Returns:
209            HTML text with the image and some properties.
210        """
211        import io
212        import base64
213        from PIL import Image
214
215        library_name = "vedo.tetmesh.TetMesh"
216        help_url = "https://vedo.embl.es/docs/vedo/tetmesh.html"
217
218        arr = self.thumbnail()
219        im = Image.fromarray(arr)
220        buffered = io.BytesIO()
221        im.save(buffered, format="PNG", quality=100)
222        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
223        url = "data:image/png;base64," + encoded
224        image = f"<img src='{url}'></img>"
225
226        bounds = "<br/>".join(
227            [
228                utils.precision(min_x,4) + " ... " + utils.precision(max_x,4)
229                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
230            ]
231        )
232
233        help_text = ""
234        if self.name:
235            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
236        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
237        if self.filename:
238            dots = ""
239            if len(self.filename) > 30:
240                dots = "..."
241            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
242
243        pdata = ""
244        if self._data.GetPointData().GetScalars():
245            if self._data.GetPointData().GetScalars().GetName():
246                name = self._data.GetPointData().GetScalars().GetName()
247                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
248
249        cdata = ""
250        if self._data.GetCellData().GetScalars():
251            if self._data.GetCellData().GetScalars().GetName():
252                name = self._data.GetCellData().GetScalars().GetName()
253                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
254
255        pts = self.points()
256        cm = np.mean(pts, axis=0)
257
258        allt = [
259            "<table>",
260            "<tr>",
261            "<td>", image, "</td>",
262            "<td style='text-align: center; vertical-align: center;'><br/>", help_text,
263            "<table>",
264            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
265            "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>",
266            "<tr><td><b> nr. points&nbsp/&nbsptets </b></td><td>"
267            + str(self.npoints) + "&nbsp/&nbsp" + str(self.ncells) + "</td></tr>",
268            pdata,
269            cdata,
270            "</table>",
271            "</table>",
272        ]
273        return "\n".join(allt)
274
275
276    def _update(self, data):
277        self._data = data
278        self.mapper().SetInputData(data)
279        self.mapper().Modified()
280        return self
281
282    def clone(self):
283        """Clone the `TetMesh` object to yield an exact copy."""
284        ugCopy = vtk.vtkUnstructuredGrid()
285        ugCopy.DeepCopy(self._data)
286
287        cloned = TetMesh(ugCopy)
288        pr = vtk.vtkVolumeProperty()
289        pr.DeepCopy(self.GetProperty())
290        cloned.SetProperty(pr)
291
292        # assign the same transformation to the copy
293        cloned.SetOrigin(self.GetOrigin())
294        cloned.SetScale(self.GetScale())
295        cloned.SetOrientation(self.GetOrientation())
296        cloned.SetPosition(self.GetPosition())
297
298        cloned.mapper().SetScalarMode(self.mapper().GetScalarMode())
299        cloned.name = self.name
300
301        cloned.pipeline = utils.OperationNode(
302            "clone", c="#edabab", shape="diamond", parents=[self],
303        )
304        return cloned
305
306    def compute_quality(self, metric=7):
307        """
308        Calculate functions of quality for the elements of a triangular mesh.
309        This method adds to the mesh a cell array named "Quality".
310
311        Arguments:
312            metric : (int)
313                type of estimators:
314                    - EDGE RATIO, 0
315                    - ASPECT RATIO, 1
316                    - RADIUS RATIO, 2
317                    - ASPECT FROBENIUS, 3
318                    - MIN_ANGLE, 4
319                    - COLLAPSE RATIO, 5
320                    - ASPECT GAMMA, 6
321                    - VOLUME, 7
322                    - ...
323
324        See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html)
325        for an explanation of the meaning of each metric..
326        """
327        qf = vtk.vtkMeshQuality()
328        qf.SetInputData(self._data)
329        qf.SetTetQualityMeasure(metric)
330        qf.SaveCellQualityOn()
331        qf.Update()
332        self._update(qf.GetOutput())
333        return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality"))
334
335    def compute_tets_volume(self):
336        """Add to this mesh a cell data array containing the tetrahedron volume."""
337        csf = vtk.vtkCellSizeFilter()
338        csf.SetInputData(self._data)
339        csf.SetComputeArea(False)
340        csf.SetComputeVolume(True)
341        csf.SetComputeLength(False)
342        csf.SetComputeVertexCount(False)
343        csf.SetVolumeArrayName("TetVolume")
344        csf.Update()
345        self._update(csf.GetOutput())
346        return utils.vtk2numpy(csf.GetOutput().GetCellData().GetArray("TetVolume"))
347
348    def check_validity(self, tol=0):
349        """
350        Return an array of possible problematic tets following this convention:
351        ```python
352        Valid               =  0
353        WrongNumberOfPoints = 01
354        IntersectingEdges   = 02
355        IntersectingFaces   = 04
356        NoncontiguousEdges  = 08
357        Nonconvex           = 10
358        OrientedIncorrectly = 20
359        ```
360
361        Arguments:
362            tol : (float)
363                This value is used as an epsilon for floating point
364                equality checks throughout the cell checking process.
365        """
366        vald = vtk.vtkCellValidator()
367        if tol:
368            vald.SetTolerance(tol)
369        vald.SetInputData(self._data)
370        vald.Update()
371        varr = vald.GetOutput().GetCellData().GetArray("ValidityState")
372        return utils.vtk2numpy(varr)
373
374    def threshold(self, name=None, above=None, below=None, on="cells"):
375        """
376        Threshold the tetrahedral mesh by a cell scalar value.
377        Reduce to only tets which satisfy the threshold limits.
378
379        - if `above = below` will only select tets with that specific value.
380        - if `above > below` selection range is flipped.
381
382        Set keyword "on" to either "cells" or "points".
383        """
384        th = vtk.vtkThreshold()
385        th.SetInputData(self._data)
386
387        if name is None:
388            if self.celldata.keys():
389                name = self.celldata.keys()[0]
390                th.SetInputArrayToProcess(0, 0, 0, 1, name)
391            elif self.pointdata.keys():
392                name = self.pointdata.keys()[0]
393                th.SetInputArrayToProcess(0, 0, 0, 0, name)
394            if name is None:
395                vedo.logger.warning("cannot find active array. Skip.")
396                return self
397        else:
398            if on.startswith("c"):
399                th.SetInputArrayToProcess(0, 0, 0, 1, name)
400            else:
401                th.SetInputArrayToProcess(0, 0, 0, 0, name)
402
403        if above is not None and below is not None:
404            if above > below:
405                if vedo.vtk_version[0] >= 9:
406                    th.SetInvert(True)
407                    th.ThresholdBetween(below, above)
408                else:
409                    vedo.logger.error("in vtk<9, above cannot be larger than below. Skip")
410                    return self
411            else:
412                th.ThresholdBetween(above, below)
413
414        elif above is not None:
415            th.ThresholdByUpper(above)
416
417        elif below is not None:
418            th.ThresholdByLower(below)
419
420        th.Update()
421        return self._update(th.GetOutput())
422
423    def decimate(self, scalars_name, fraction=0.5, n=0):
424        """
425        Downsample the number of tets in a TetMesh to a specified fraction.
426        Either `fraction` or `n` must be set.
427
428        Arguments:
429            fraction : (float)
430                the desired final fraction of the total.
431            n : (int)
432                the desired number of final tets
433
434        .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets.
435        """
436        decimate = vtk.vtkUnstructuredGridQuadricDecimation()
437        decimate.SetInputData(self._data)
438        decimate.SetScalarsName(scalars_name)
439
440        if n:  # n = desired number of points
441            decimate.SetNumberOfTetsOutput(n)
442        else:
443            decimate.SetTargetReduction(1 - fraction)
444        decimate.Update()
445        self._update(decimate.GetOutput())
446        self.pipeline = utils.OperationNode(
447            "decimate", comment=f"array: {scalars_name}",
448            c="#edabab", parents=[self],
449        )
450        return self
451
452    def subdvide(self):
453        """
454        Increase the number of tets of a `TetMesh`.
455        Subdivide one tetrahedron into twelve for every tetra.
456        """
457        sd = vtk.vtkSubdivideTetra()
458        sd.SetInputData(self._data)
459        sd.Update()
460        self._update(sd.GetOutput())
461        self.pipeline = utils.OperationNode(
462            "subdvide", c="#edabab", parents=[self],
463        )
464        return self
465
466    def isosurface(self, value=True):
467        """
468        Return a `vedo.Mesh` isosurface.
469
470        Set `value` to a single value or list of values to compute the isosurface(s).
471        """
472        if not self._data.GetPointData().GetScalars():
473            self.map_cells_to_points()
474        scrange = self._data.GetPointData().GetScalars().GetRange()
475        cf = vtk.vtkContourFilter()  # vtk.vtkContourGrid()
476        cf.SetInputData(self._data)
477
478        if utils.is_sequence(value):
479            cf.SetNumberOfContours(len(value))
480            for i, t in enumerate(value):
481                cf.SetValue(i, t)
482            cf.Update()
483        else:
484            if value is True:
485                value = (2 * scrange[0] + scrange[1]) / 3.0
486            cf.SetValue(0, value)
487            cf.Update()
488
489        clp = vtk.vtkCleanPolyData()
490        clp.SetInputData(cf.GetOutput())
491        clp.Update()
492        msh = Mesh(clp.GetOutput(), c=None).phong()
493        msh.mapper().SetLookupTable(utils.ctf2lut(self))
494        msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self])
495        return msh
496
497    def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)):
498        """
499        Return a 2D slice of the mesh by a plane passing through origin and assigned normal.
500        """
501        strn = str(normal)
502        if strn   ==  "x": normal = (1, 0, 0)
503        elif strn ==  "y": normal = (0, 1, 0)
504        elif strn ==  "z": normal = (0, 0, 1)
505        elif strn == "-x": normal = (-1, 0, 0)
506        elif strn == "-y": normal = (0, -1, 0)
507        elif strn == "-z": normal = (0, 0, -1)
508        plane = vtk.vtkPlane()
509        plane.SetOrigin(origin)
510        plane.SetNormal(normal)
511
512        cc = vtk.vtkCutter()
513        cc.SetInputData(self._data)
514        cc.SetCutFunction(plane)
515        cc.Update()
516        msh = Mesh(cc.GetOutput()).flat().lighting("ambient")
517        msh.mapper().SetLookupTable(utils.ctf2lut(self))
518        msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self])
519        return msh