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.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        all = [
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> average size </b></td><td>" + str(average_size) + "</td></tr>",
267            "<tr><td><b> nr. points&nbsp/&nbsptets </b></td><td>" 
268            + str(self.npoints) + "&nbsp/&nbsp" + str(self.ncells) + "</td></tr>",
269            pdata,
270            cdata,
271            "</table>",
272            "</table>",
273        ]
274        return "\n".join(all)
275
276
277    def _update(self, data):
278        self._data = data
279        self.mapper().SetInputData(data)
280        self.mapper().Modified()
281        return self
282
283    def clone(self):
284        """Clone the `TetMesh` object to yield an exact copy."""
285        ugCopy = vtk.vtkUnstructuredGrid()
286        ugCopy.DeepCopy(self._data)
287
288        cloned = TetMesh(ugCopy)
289        pr = vtk.vtkVolumeProperty()
290        pr.DeepCopy(self.GetProperty())
291        cloned.SetProperty(pr)
292
293        # assign the same transformation to the copy
294        cloned.SetOrigin(self.GetOrigin())
295        cloned.SetScale(self.GetScale())
296        cloned.SetOrientation(self.GetOrientation())
297        cloned.SetPosition(self.GetPosition())
298
299        cloned.mapper().SetScalarMode(self.mapper().GetScalarMode())
300        cloned.name = self.name
301
302        cloned.pipeline = utils.OperationNode(
303            "clone", c="#edabab", shape="diamond", parents=[self],
304        )
305        return cloned
306
307    def compute_quality(self, metric=7):
308        """
309        Calculate functions of quality for the elements of a triangular mesh.
310        This method adds to the mesh a cell array named "Quality".
311
312        Arguments:
313            metric : (int)
314                type of estimators:
315                    - EDGE RATIO, 0
316                    - ASPECT RATIO, 1
317                    - RADIUS RATIO, 2
318                    - ASPECT FROBENIUS, 3
319                    - MIN_ANGLE, 4
320                    - COLLAPSE RATIO, 5
321                    - ASPECT GAMMA, 6
322                    - VOLUME, 7
323                    - ...
324
325        See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html)
326        for an explanation of the meaning of each metric..
327        """
328        qf = vtk.vtkMeshQuality()
329        qf.SetInputData(self._data)
330        qf.SetTetQualityMeasure(metric)
331        qf.SaveCellQualityOn()
332        qf.Update()
333        self._update(qf.GetOutput())
334        return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality"))
335
336    def compute_tets_volume(self):
337        """Add to this mesh a cell data array containing the tetrahedron volume."""
338        csf = vtk.vtkCellSizeFilter()
339        csf.SetInputData(self._data)
340        csf.SetComputeArea(False)
341        csf.SetComputeVolume(True)
342        csf.SetComputeLength(False)
343        csf.SetComputeVertexCount(False)
344        csf.SetVolumeArrayName("TetVolume")
345        csf.Update()
346        self._update(csf.GetOutput())
347        return utils.vtk2numpy(csf.GetOutput().GetCellData().GetArray("TetVolume"))
348
349    def check_validity(self, tol=0):
350        """
351        Return an array of possible problematic tets following this convention:
352        ```python
353        Valid               =  0
354        WrongNumberOfPoints = 01
355        IntersectingEdges   = 02
356        IntersectingFaces   = 04
357        NoncontiguousEdges  = 08
358        Nonconvex           = 10
359        OrientedIncorrectly = 20
360        ```
361
362        Arguments:
363            tol : (float)
364                This value is used as an epsilon for floating point
365                equality checks throughout the cell checking process.
366        """
367        vald = vtk.vtkCellValidator()
368        if tol:
369            vald.SetTolerance(tol)
370        vald.SetInputData(self._data)
371        vald.Update()
372        varr = vald.GetOutput().GetCellData().GetArray("ValidityState")
373        return utils.vtk2numpy(varr)
374
375    def threshold(self, name=None, above=None, below=None, on="cells"):
376        """
377        Threshold the tetrahedral mesh by a cell scalar value.
378        Reduce to only tets which satisfy the threshold limits.
379
380        - if `above = below` will only select tets with that specific value.
381        - if `above > below` selection range is flipped.
382
383        Set keyword "on" to either "cells" or "points".
384        """
385        th = vtk.vtkThreshold()
386        th.SetInputData(self._data)
387
388        if name is None:
389            if self.celldata.keys():
390                name = self.celldata.keys()[0]
391                th.SetInputArrayToProcess(0, 0, 0, 1, name)
392            elif self.pointdata.keys():
393                name = self.pointdata.keys()[0]
394                th.SetInputArrayToProcess(0, 0, 0, 0, name)
395            if name is None:
396                vedo.logger.warning("cannot find active array. Skip.")
397                return self
398        else:
399            if on.startswith("c"):
400                th.SetInputArrayToProcess(0, 0, 0, 1, name)
401            else:
402                th.SetInputArrayToProcess(0, 0, 0, 0, name)
403
404        if above is not None and below is not None:
405            if above > below:
406                if vedo.vtk_version[0] >= 9:
407                    th.SetInvert(True)
408                    th.ThresholdBetween(below, above)
409                else:
410                    vedo.logger.error("in vtk<9, above cannot be larger than below. Skip")
411                    return self
412            else:
413                th.ThresholdBetween(above, below)
414
415        elif above is not None:
416            th.ThresholdByUpper(above)
417
418        elif below is not None:
419            th.ThresholdByLower(below)
420
421        th.Update()
422        return self._update(th.GetOutput())
423
424    def decimate(self, scalars_name, fraction=0.5, n=0):
425        """
426        Downsample the number of tets in a TetMesh to a specified fraction.
427        Either `fraction` or `n` must be set.
428
429        Arguments:
430            fraction : (float)
431                the desired final fraction of the total.
432            n : (int)
433                the desired number of final tets
434
435        .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets.
436        """
437        decimate = vtk.vtkUnstructuredGridQuadricDecimation()
438        decimate.SetInputData(self._data)
439        decimate.SetScalarsName(scalars_name)
440
441        if n:  # n = desired number of points
442            decimate.SetNumberOfTetsOutput(n)
443        else:
444            decimate.SetTargetReduction(1 - fraction)
445        decimate.Update()
446        self._update(decimate.GetOutput())
447        self.pipeline = utils.OperationNode(
448            "decimate", comment=f"array: {scalars_name}",
449            c="#edabab", parents=[self],
450        )
451        return self
452
453    def subdvide(self):
454        """
455        Increase the number of tets of a `TetMesh`.
456        Subdivide one tetrahedron into twelve for every tetra.
457        """
458        sd = vtk.vtkSubdivideTetra()
459        sd.SetInputData(self._data)
460        sd.Update()
461        self._update(sd.GetOutput())
462        self.pipeline = utils.OperationNode(
463            "subdvide", c="#edabab", parents=[self],
464        )
465        return self
466
467    def isosurface(self, value=True):
468        """
469        Return a `vedo.Mesh` isosurface.
470
471        Set `value` to a single value or list of values to compute the isosurface(s).
472        """
473        if not self._data.GetPointData().GetScalars():
474            self.map_cells_to_points()
475        scrange = self._data.GetPointData().GetScalars().GetRange()
476        cf = vtk.vtkContourFilter()  # vtk.vtkContourGrid()
477        cf.SetInputData(self._data)
478
479        if utils.is_sequence(value):
480            cf.SetNumberOfContours(len(value))
481            for i, t in enumerate(value):
482                cf.SetValue(i, t)
483            cf.Update()
484        else:
485            if value is True:
486                value = (2 * scrange[0] + scrange[1]) / 3.0
487            cf.SetValue(0, value)
488            cf.Update()
489
490        clp = vtk.vtkCleanPolyData()
491        clp.SetInputData(cf.GetOutput())
492        clp.Update()
493        msh = Mesh(clp.GetOutput(), c=None).phong()
494        msh.mapper().SetLookupTable(utils.ctf2lut(self))
495        msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self])
496        return msh
497
498    def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)):
499        """
500        Return a 2D slice of the mesh by a plane passing through origin and assigned normal.
501        """
502        strn = str(normal)
503        if strn   ==  "x": normal = (1, 0, 0)
504        elif strn ==  "y": normal = (0, 1, 0)
505        elif strn ==  "z": normal = (0, 0, 1)
506        elif strn == "-x": normal = (-1, 0, 0)
507        elif strn == "-y": normal = (0, -1, 0)
508        elif strn == "-z": normal = (0, 0, -1)
509        plane = vtk.vtkPlane()
510        plane.SetOrigin(origin)
511        plane.SetNormal(normal)
512
513        cc = vtk.vtkCutter()
514        cc.SetInputData(self._data)
515        cc.SetCutFunction(plane)
516        cc.Update()
517        msh = Mesh(cc.GetOutput()).flat().lighting("ambient")
518        msh.mapper().SetLookupTable(utils.ctf2lut(self))
519        msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self])
520        return msh