vedo.shapes

Submodule to generate simple and complex geometric shapes

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import os
   4from functools import lru_cache
   5
   6import numpy as np
   7
   8try:
   9    import vedo.vtkclasses as vtk
  10except ImportError:
  11    import vtkmodules.all as vtk
  12
  13import vedo
  14from vedo import settings
  15from vedo.colors import cmaps_names, color_map, get_color, printc
  16from vedo import utils
  17from vedo.pointcloud import Points, merge
  18from vedo.mesh import Mesh
  19from vedo.picture import Picture
  20
  21
  22__docformat__ = "google"
  23
  24__doc__ = """
  25Submodule to generate simple and complex geometric shapes
  26
  27![](https://vedo.embl.es/images/basic/extrude.png)
  28"""
  29
  30__all__ = [
  31    "Marker",
  32    "Line",
  33    "DashedLine",
  34    "RoundedLine",
  35    "Tube",
  36    "ThickTube",
  37    "Lines",
  38    "Spline",
  39    "KSpline",
  40    "CSpline",
  41    "Bezier",
  42    "Brace",
  43    "NormalLines",
  44    "StreamLines",
  45    "Ribbon",
  46    "Arrow",
  47    "Arrows",
  48    "Arrow2D",
  49    "Arrows2D",
  50    "FlatArrow",
  51    "Polygon",
  52    "Triangle",
  53    "Rectangle",
  54    "Disc",
  55    "Circle",
  56    "GeoCircle",
  57    "Arc",
  58    "Star",
  59    "Star3D",
  60    "Cross3D",
  61    "IcoSphere",
  62    "Sphere",
  63    "Spheres",
  64    "Earth",
  65    "Ellipsoid",
  66    "Grid",
  67    "TessellatedBox",
  68    "Plane",
  69    "Box",
  70    "Cube",
  71    "Spring",
  72    "Cylinder",
  73    "Cone",
  74    "Pyramid",
  75    "Torus",
  76    "Paraboloid",
  77    "Hyperboloid",
  78    "TextBase",
  79    "Text3D",
  80    "Text2D",
  81    "CornerAnnotation",
  82    "Latex",
  83    "Glyph",
  84    "Tensors",
  85    "ParametricShape",
  86    "ConvexHull",
  87    "VedoLogo",
  88]
  89
  90##############################################
  91_reps = (
  92    (":nabla", "∇"),
  93    (":inf", "∞"),
  94    (":rightarrow", "→"),
  95    (":lefttarrow", "←"),
  96    (":partial", "∂"),
  97    (":sqrt", "√"),
  98    (":approx", "≈"),
  99    (":neq", "≠"),
 100    (":leq", "≤"),
 101    (":geq", "≥"),
 102    (":foreach", "∀"),
 103    (":permille", "‰"),
 104    (":euro", "€"),
 105    (":dot", "·"),
 106    (":int", "∫"),
 107    (":pm", "±"),
 108    (":times", "×"),
 109    (":Gamma", "Γ"),
 110    (":Delta", "Δ"),
 111    (":Theta", "Θ"),
 112    (":Lambda", "Λ"),
 113    (":Pi", "Π"),
 114    (":Sigma", "Σ"),
 115    (":Phi", "Φ"),
 116    (":Chi", "X"),
 117    (":Xi", "Ξ"),
 118    (":Psi", "Ψ"),
 119    (":Omega", "Ω"),
 120    (":alpha", "α"),
 121    (":beta", "β"),
 122    (":gamma", "γ"),
 123    (":delta", "δ"),
 124    (":epsilon", "ε"),
 125    (":zeta", "ζ"),
 126    (":eta", "η"),
 127    (":theta", "θ"),
 128    (":kappa", "κ"),
 129    (":lambda", "λ"),
 130    (":mu", "μ"),
 131    (":lowerxi", "ξ"),
 132    (":nu", "ν"),
 133    (":pi", "π"),
 134    (":rho", "ρ"),
 135    (":sigma", "σ"),
 136    (":tau", "τ"),
 137    (":varphi", "φ"),
 138    (":phi", "φ"),
 139    (":chi", "χ"),
 140    (":psi", "ψ"),
 141    (":omega", "ω"),
 142    (":circ", "°"),
 143    (":onehalf", "½"),
 144    (":onefourth", "¼"),
 145    (":threefourths", "¾"),
 146    (":^1", "¹"),
 147    (":^2", "²"),
 148    (":^3", "³"),
 149    (":,", "~"),
 150)
 151
 152
 153########################################################################
 154class Glyph(Mesh):
 155    """
 156    At each vertex of a mesh, another mesh, i.e. a "glyph", is shown with
 157    various orientation options and coloring.
 158
 159    The input can also be a simple list of 2D or 3D coordinates.
 160    Color can be specified as a colormap which maps the size of the orientation
 161    vectors in `orientation_array`.
 162    """
 163    def __init__(
 164            self,
 165            mesh,
 166            glyph,
 167            orientation_array=None,
 168            scale_by_scalar=False,
 169            scale_by_vector_size=False,
 170            scale_by_vector_components=False,
 171            color_by_scalar=False,
 172            color_by_vector_size=False,
 173            tol=0,
 174            c='k8',
 175            alpha=1,
 176            **opts,
 177        ):
 178        """
 179        Arguments:
 180            orientation_array: (list, str, vtkArray)
 181                list of vectors, `vtkArray` or name of an already existing pointdata array
 182            scale_by_scalar : (bool)
 183                glyph mesh is scaled by the active scalars
 184            scale_by_vector_size : (bool)
 185                glyph mesh is scaled by the size of the vectors
 186            scale_by_vector_components : (bool)
 187                glyph mesh is scaled by the 3 vectors components
 188            color_by_scalar : (bool)
 189                glyph mesh is colored based on the scalar value
 190            color_by_vector_size : (bool)
 191                glyph mesh is colored based on the vector size
 192            tol : (float)
 193                set a minimum separation between two close glyphs
 194                (not compatible with `orientation_array` being a list).
 195
 196        Examples:
 197            - [glyphs.py](]https://github.com/marcomusy/vedo/tree/master/examples/basic/glyphs.py)
 198            - [glyphs_arrows.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/glyphs_arrows.py)
 199
 200            ![](https://vedo.embl.es/images/basic/glyphs.png)
 201        """
 202        if len(opts): # Deprecations
 203            printc(" Warning! In Glyph() unrecognized keywords:", opts, c='y')
 204            orientation_array = opts.pop("orientationArray", orientation_array)
 205            scale_by_scalar = opts.pop("scaleByScalar", scale_by_scalar)
 206            scale_by_vector_size = opts.pop("scaleByVectorSize", scale_by_vector_size)
 207            scale_by_vector_components = opts.pop("scaleByVectorComponents", scale_by_vector_components)
 208            color_by_scalar = opts.pop("colorByScalar", color_by_scalar)
 209            color_by_vector_size = opts.pop("colorByVectorSize", color_by_vector_size)
 210            printc("          Please use 'snake_case' instead of 'camelCase' keywords", c='y')
 211
 212        lighting = None
 213        if utils.is_sequence(mesh):
 214            # create a cloud of points
 215            poly = Points(mesh).polydata()
 216        elif isinstance(mesh, vtk.vtkPolyData):
 217            poly = mesh
 218        else:
 219            poly = mesh.polydata()
 220
 221        if tol:
 222            cpd = vtk.vtkCleanPolyData()
 223            cpd.SetInputData(poly)
 224            cpd.SetTolerance(tol)
 225            cpd.Update()
 226            poly = cpd.GetOutput()
 227
 228        if isinstance(glyph, Points):
 229            lighting = glyph.property.GetLighting()
 230            glyph = glyph.polydata()
 231
 232        cmap = ""
 233        if isinstance(c, str) and c in cmaps_names:
 234            cmap = c
 235            c = None
 236        elif utils.is_sequence(c):  # user passing an array of point colors
 237            ucols = vtk.vtkUnsignedCharArray()
 238            ucols.SetNumberOfComponents(3)
 239            ucols.SetName("glyph_RGB")
 240            for col in c:
 241                cl = get_color(col)
 242                ucols.InsertNextTuple3(cl[0] * 255, cl[1] * 255, cl[2] * 255)
 243            poly.GetPointData().AddArray(ucols)
 244            poly.GetPointData().SetActiveScalars("glyph_RGB")
 245            c = None
 246
 247        gly = vtk.vtkGlyph3D()
 248        gly.SetInputData(poly)
 249        gly.SetSourceData(glyph)
 250
 251        if scale_by_scalar:
 252            gly.SetScaleModeToScaleByScalar()
 253        elif scale_by_vector_size:
 254            gly.SetScaleModeToScaleByVector()
 255        elif scale_by_vector_components:
 256            gly.SetScaleModeToScaleByVectorComponents()
 257        else:
 258            gly.SetScaleModeToDataScalingOff()
 259
 260        if color_by_vector_size:
 261            gly.SetVectorModeToUseVector()
 262            gly.SetColorModeToColorByVector()
 263        elif color_by_scalar:
 264            gly.SetColorModeToColorByScalar()
 265        else:
 266            gly.SetColorModeToColorByScale()
 267
 268        if orientation_array is not None:
 269            gly.OrientOn()
 270            if isinstance(orientation_array, str):
 271                if orientation_array.lower() == "normals":
 272                    gly.SetVectorModeToUseNormal()
 273                else:  # passing a name
 274                    poly.GetPointData().SetActiveVectors(orientation_array)
 275                    gly.SetInputArrayToProcess(0, 0, 0, 0, orientation_array)
 276                    gly.SetVectorModeToUseVector()
 277            elif utils.is_sequence(orientation_array) and not tol:  # passing a list
 278                varr = vtk.vtkFloatArray()
 279                varr.SetNumberOfComponents(3)
 280                varr.SetName("glyph_vectors")
 281                for v in orientation_array:
 282                    varr.InsertNextTuple(v)
 283                poly.GetPointData().AddArray(varr)
 284                poly.GetPointData().SetActiveVectors("glyph_vectors")
 285                gly.SetInputArrayToProcess(0, 0, 0, 0, "glyph_vectors")
 286                gly.SetVectorModeToUseVector()
 287
 288        gly.Update()
 289
 290        Mesh.__init__(self, gly.GetOutput(), c, alpha)
 291        self.flat()
 292        if lighting is not None:
 293            self.property.SetLighting(lighting)
 294
 295        if cmap:
 296            lut = vtk.vtkLookupTable()
 297            lut.SetNumberOfTableValues(512)
 298            lut.Build()
 299            for i in range(512):
 300                r, g, b = color_map(i, cmap, 0, 512)
 301                lut.SetTableValue(i, r, g, b, 1)
 302            self.mapper().SetLookupTable(lut)
 303            self.mapper().ScalarVisibilityOn()
 304            self.mapper().SetScalarModeToUsePointData()
 305            if gly.GetOutput().GetPointData().GetScalars():
 306                rng = gly.GetOutput().GetPointData().GetScalars().GetRange()
 307                self.mapper().SetScalarRange(rng[0], rng[1])
 308
 309        self.name = "Glyph"
 310
 311
 312class Tensors(Mesh):
 313    """
 314    Geometric representation of tensors defined on a domain or set of points.
 315    Tensors can be scaled and/or rotated according to the source at eache input point.
 316    Scaling and rotation is controlled by the eigenvalues/eigenvectors of the symmetrical part
 317    of the tensor as follows:
 318
 319    For each tensor, the eigenvalues (and associated eigenvectors) are sorted
 320    to determine the major, medium, and minor eigenvalues/eigenvectors.
 321    The eigenvalue decomposition only makes sense for symmetric tensors,
 322    hence the need to only consider the symmetric part of the tensor,
 323    which is `1/2*(T+T.transposed())`.
 324    """
 325    def __init__(
 326        self,
 327        domain,
 328        source="ellipsoid",
 329        use_eigenvalues=True,
 330        is_symmetric=True,
 331        three_axes=False,
 332        scale=1,
 333        max_scale=None,
 334        length=None,
 335        c=None,
 336        alpha=1,
 337    ):
 338        """
 339        Arguments:
 340            source : (str, Mesh)
 341                preset type of source shape
 342                `['ellipsoid', 'cylinder', 'cube' or any specified `Mesh`]`
 343            use_eigenvalues : (bool)
 344                color source glyph using the eigenvalues or by scalars
 345            three_axes : (bool)
 346                if `False` scale the source in the x-direction,
 347                the medium in the y-direction, and the minor in the z-direction.
 348                Then, the source is rotated so that the glyph's local x-axis lies
 349                along the major eigenvector, y-axis along the medium eigenvector,
 350                and z-axis along the minor.
 351
 352                If `True` three sources are produced, each of them oriented along an eigenvector
 353                and scaled according to the corresponding eigenvector.
 354            is_symmetric : (bool)
 355                If `True` each source glyph is mirrored (2 or 6 glyphs will be produced).
 356                The x-axis of the source glyph will correspond to the eigenvector on output.
 357            length : (float)
 358                distance from the origin to the tip of the source glyph along the x-axis
 359            scale : (float)
 360                scaling factor of the source glyph.
 361            max_scale : (float)
 362                clamp scaling at this factor.
 363
 364        Examples:
 365            - [tensors.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/tensors.py)
 366            - [tensor_grid.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/tensor_grid.py)
 367
 368            ![](https://vedo.embl.es/images/volumetric/tensor_grid.png)
 369        """
 370        if isinstance(source, Points):
 371            src = source.normalize().polydata(False)
 372        else:
 373            if "ellip" in source:
 374                src = vtk.vtkSphereSource()
 375                src.SetPhiResolution(24)
 376                src.SetThetaResolution(12)
 377            elif "cyl" in source:
 378                src = vtk.vtkCylinderSource()
 379                src.SetResolution(48)
 380                src.CappingOn()
 381            elif source == "cube":
 382                src = vtk.vtkCubeSource()
 383            src.Update()
 384
 385        tg = vtk.vtkTensorGlyph()
 386        if isinstance(domain, vtk.vtkPolyData):
 387            tg.SetInputData(domain)
 388        else:
 389            tg.SetInputData(domain.GetMapper().GetInput())
 390        tg.SetSourceData(src.GetOutput())
 391
 392        if c is None:
 393            tg.ColorGlyphsOn()
 394        else:
 395            tg.ColorGlyphsOff()
 396
 397        tg.SetSymmetric(int(is_symmetric))
 398
 399        if length is not None:
 400            tg.SetLength(length)
 401        if use_eigenvalues:
 402            tg.ExtractEigenvaluesOn()
 403            tg.SetColorModeToEigenvalues()
 404        else:
 405            tg.SetColorModeToScalars()
 406        tg.SetThreeGlyphs(three_axes)
 407        tg.ScalingOn()
 408        tg.SetScaleFactor(scale)
 409        if max_scale is None:
 410            tg.ClampScalingOn()
 411            max_scale = scale * 10
 412        tg.SetMaxScaleFactor(max_scale)
 413        tg.Update()
 414        tgn = vtk.vtkPolyDataNormals()
 415        tgn.SetInputData(tg.GetOutput())
 416        tgn.Update()
 417        Mesh.__init__(self, tgn.GetOutput(), c, alpha)
 418        self.name = "Tensors"
 419
 420
 421class Line(Mesh):
 422    """
 423    Build the line segment between points `p0` and `p1`.
 424
 425    If `p0` is already a list of points, return the line connecting them.
 426
 427    A 2D set of coords can also be passed as `p0=[x..], p1=[y..]`.
 428    """
 429    def __init__(self, p0, p1=None, closed=False, res=2, lw=1, c="k1", alpha=1):
 430        """
 431        Arguments:
 432            closed : (bool)
 433                join last to first point
 434            res : (int)
 435                resolution, number of points along the line
 436                (only relevant if only 2 points are specified)
 437            lw : (int)
 438                line width in pixel units
 439            c : (color), int, str, list
 440                color name, number, or list of [R,G,B] colors
 441            alpha : (float)
 442                opacity in range [0,1]
 443        """
 444        self.slope = []  # populated by analysis.fitLine
 445        self.center = []
 446        self.variances = []
 447
 448        self.coefficients = []  # populated by pyplot.fit()
 449        self.covariance_matrix = []
 450        self.coefficients = []
 451        self.coefficient_errors = []
 452        self.monte_carlo_coefficients = []
 453        self.reduced_chi2 = -1
 454        self.ndof = 0
 455        self.data_sigma = 0
 456        self.error_lines = []
 457        self.error_band = None
 458        self.res = res
 459
 460        if isinstance(p1, Points):
 461            p1 = p1.GetPosition()
 462            if isinstance(p0, Points):
 463                p0 = p0.GetPosition()
 464        if isinstance(p0, Points):
 465            p0 = p0.points()
 466
 467        # detect if user is passing a 2D list of points as p0=xlist, p1=ylist:
 468        if len(p0) > 3:
 469            if not utils.is_sequence(p0[0]) and not utils.is_sequence(p1[0]) and len(p0)==len(p1):
 470                # assume input is 2D xlist, ylist
 471                p0 = np.stack((p0, p1), axis=1)
 472                p1 = None
 473            p0 = utils.make3d(p0)
 474
 475        # detect if user is passing a list of points:
 476        if utils.is_sequence(p0[0]):
 477            p0 = utils.make3d(p0)
 478
 479            ppoints = vtk.vtkPoints()  # Generate the polyline
 480            ppoints.SetData(utils.numpy2vtk(np.asarray(p0), dtype=np.float32))
 481            lines = vtk.vtkCellArray()
 482            npt = len(p0)
 483            if closed:
 484                lines.InsertNextCell(npt + 1)
 485            else:
 486                lines.InsertNextCell(npt)
 487            for i in range(npt):
 488                lines.InsertCellPoint(i)
 489            if closed:
 490                lines.InsertCellPoint(0)
 491            poly = vtk.vtkPolyData()
 492            poly.SetPoints(ppoints)
 493            poly.SetLines(lines)
 494            top = p0[-1]
 495            base = p0[0]
 496            self.res = 2
 497
 498        else:  # or just 2 points to link
 499
 500            line_source = vtk.vtkLineSource()
 501            p0 = utils.make3d(p0)
 502            p1 = utils.make3d(p1)
 503            line_source.SetPoint1(p0)
 504            line_source.SetPoint2(p1)
 505            line_source.SetResolution(res - 1)
 506            line_source.Update()
 507            poly = line_source.GetOutput()
 508            top = np.asarray(p1, dtype=float)
 509            base = np.asarray(p0, dtype=float)
 510
 511        Mesh.__init__(self, poly, c, alpha)
 512        self.lw(lw)
 513        self.property.LightingOff()
 514        self.PickableOff()
 515        self.DragableOff()
 516        self.base = base
 517        self.top = top
 518        self.name = "Line"
 519
 520    def linecolor(self, lc=None):
 521        """Assign a color to the line"""
 522        # overrides mesh.linecolor which would have no effect here
 523        return self.color(lc)
 524
 525    def eval(self, x):
 526        """
 527        Calculate the position of an intermediate point
 528        as a fraction of the length of the line,
 529        being x=0 the first point and x=1 the last point.
 530        This corresponds to an imaginary point that travels along the line
 531        at constant speed.
 532
 533        Can be used in conjunction with `lin_interpolate()`
 534        to map any range to the [0,1] range.
 535        """
 536        distance1 = 0.0
 537        length = self.length()
 538        pts = self.points()
 539        for i in range(1, len(pts)):
 540            p0 = pts[i - 1]
 541            p1 = pts[i]
 542            seg = p1 - p0
 543            distance0 = distance1
 544            distance1 += np.linalg.norm(seg)
 545            w1 = distance1 / length
 546            if w1 >= x:
 547                break
 548        w0 = distance0 / length
 549        v = p0 + seg * (x - w0) / (w1 - w0)
 550        return v
 551
 552    def pattern(self, stipple, repeats=10):
 553        """
 554        Define a stipple pattern for dashing the line.
 555        Pass the stipple pattern as a string like `'- - -'`.
 556        Repeats controls the number of times the pattern repeats in a single segment.
 557
 558        Examples are: `'- -', '--  -  --'`, etc.
 559
 560        The resolution of the line (nr of points) can affect how pattern will show up.
 561
 562        Example:
 563            ```python
 564            from vedo import Line
 565            pts = [[1, 0, 0], [5, 2, 0], [3, 3, 1]]
 566            ln = Line(pts, c='r', lw=5).pattern('- -', repeats=10)
 567            ln.show(axes=1).close()
 568            ```
 569            ![](https://vedo.embl.es/images/feats/line_pattern.png)
 570        """
 571        stipple = str(stipple) * int(2 * repeats)
 572        dimension = len(stipple)
 573
 574        image = vtk.vtkImageData()
 575        image.SetDimensions(dimension, 1, 1)
 576        image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 4)
 577        image.SetExtent(0, dimension - 1, 0, 0, 0, 0)
 578        i_dim = 0
 579        while i_dim < dimension:
 580            for i in range(dimension):
 581                image.SetScalarComponentFromFloat(i_dim, 0, 0, 0, 255)
 582                image.SetScalarComponentFromFloat(i_dim, 0, 0, 1, 255)
 583                image.SetScalarComponentFromFloat(i_dim, 0, 0, 2, 255)
 584                if stipple[i] == " ":
 585                    image.SetScalarComponentFromFloat(i_dim, 0, 0, 3, 0)
 586                else:
 587                    image.SetScalarComponentFromFloat(i_dim, 0, 0, 3, 255)
 588                i_dim += 1
 589
 590        polyData = self.polydata(False)
 591
 592        # Create texture coordinates
 593        tcoords = vtk.vtkDoubleArray()
 594        tcoords.SetName("TCoordsStippledLine")
 595        tcoords.SetNumberOfComponents(1)
 596        tcoords.SetNumberOfTuples(polyData.GetNumberOfPoints())
 597        for i in range(polyData.GetNumberOfPoints()):
 598            tcoords.SetTypedTuple(i, [i / 2])
 599        polyData.GetPointData().SetTCoords(tcoords)
 600        polyData.GetPointData().Modified()
 601        texture = vtk.vtkTexture()
 602        texture.SetInputData(image)
 603        texture.InterpolateOff()
 604        texture.RepeatOn()
 605        self.SetTexture(texture)
 606        return self
 607
 608    def length(self):
 609        """Calculate length of the line."""
 610        distance = 0.0
 611        pts = self.points()
 612        for i in range(1, len(pts)):
 613            distance += np.linalg.norm(pts[i] - pts[i - 1])
 614        return distance
 615
 616    def tangents(self):
 617        """
 618        Compute the tangents of a line in space.
 619
 620        Example:
 621            ```python
 622            from vedo import *
 623            shape = load(dataurl+"timecourse1d.npy")[58]
 624            pts = shape.rotate_x(30).points()
 625            tangents = Line(pts).tangents()
 626            arrs = Arrows(pts, pts+tangents, c='blue9')
 627            show(shape.c('red5').lw(5), arrs, bg='bb', axes=1).close()
 628            ```
 629            ![](https://vedo.embl.es/images/feats/line_tangents.png)
 630        """
 631        v = np.gradient(self.points())[0]
 632        ds_dt = np.linalg.norm(v, axis=1)
 633        tangent = np.array([1 / ds_dt] * 3).transpose() * v
 634        return tangent
 635
 636    def curvature(self):
 637        """
 638        Compute the signed curvature of a line in space.
 639        The signed is computed assuming the line is about coplanar to the xy plane.
 640
 641        Example:
 642            ```python
 643            from vedo import *
 644            from vedo.pyplot import plot
 645            shape = load(dataurl+"timecourse1d.npy")[55]
 646            curvs = Line(shape.points()).curvature()
 647            shape.cmap('coolwarm', curvs, vmin=-2,vmax=2).add_scalarbar3d(c='w')
 648            shape.render_lines_as_tubes().lw(12)
 649            pp = plot(curvs, ac='white', lc='yellow5')
 650            show(shape, pp, N=2, bg='bb', sharecam=False).close()
 651            ```
 652            ![](https://vedo.embl.es/images/feats/line_curvature.png)
 653        """
 654        v = np.gradient(self.points())[0]
 655        a = np.gradient(v)[0]
 656        av = np.cross(a, v)
 657        mav = np.linalg.norm(av, axis=1)
 658        mv = utils.mag2(v)
 659        val = mav * np.sign(av[:, 2]) / np.power(mv, 1.5)
 660        val[0] = val[1]
 661        val[-1] = val[-2]
 662        return val
 663
 664    def compute_curvature(self, method=0):
 665        """
 666        Add a pointdata array named 'Curvatures' which contains
 667        the curvature value at each point.
 668
 669        Keyword method is overridden in Mesh and has no effect here.
 670        """
 671        # overrides mesh.compute_curvature
 672        curvs = self.curvature()
 673        vmin, vmax = np.min(curvs), np.max(curvs)
 674        if vmin < 0 and vmax > 0:
 675            v = max(-vmin, vmax)
 676            self.cmap("coolwarm", curvs, vmin=-v, vmax=v, name="Curvature")
 677        else:
 678            self.cmap("coolwarm", curvs, vmin=vmin, vmax=vmax, name="Curvature")
 679        return self
 680
 681    def sweep(self, direction=(1, 0, 0), res=1):
 682        """
 683        Sweep the `Line` along the specified vector direction.
 684
 685        Returns a `Mesh` surface.
 686        Line position is updated to allow for additional sweepings.
 687
 688        Example:
 689            ```python
 690            from vedo import Line, show
 691            aline = Line([(0,0,0),(1,3,0),(2,4,0)])
 692            surf1 = aline.sweep((1,0.2,0), res=3)
 693            surf2 = aline.sweep((0.2,0,1))
 694            aline.color('r').linewidth(4)
 695            show(surf1, surf2, aline, axes=1).close()
 696            ```
 697            ![](https://vedo.embl.es/images/feats/sweepline.png)
 698        """
 699        line = self.polydata()
 700        rows = line.GetNumberOfPoints()
 701
 702        spacing = 1 / res
 703        surface = vtk.vtkPolyData()
 704
 705        res += 1
 706        npts = rows * res
 707        npolys = (rows - 1) * (res - 1)
 708        points = vtk.vtkPoints()
 709        points.Allocate(npts)
 710
 711        cnt = 0
 712        x = [0.0, 0.0, 0.0]
 713        for row in range(rows):
 714            for col in range(res):
 715                p = [0.0, 0.0, 0.0]
 716                line.GetPoint(row, p)
 717                x[0] = p[0] + direction[0] * col * spacing
 718                x[1] = p[1] + direction[1] * col * spacing
 719                x[2] = p[2] + direction[2] * col * spacing
 720                points.InsertPoint(cnt, x)
 721                cnt += 1
 722
 723        # Generate the quads
 724        polys = vtk.vtkCellArray()
 725        polys.Allocate(npolys * 4)
 726        pts = [0, 0, 0, 0]
 727        for row in range(rows - 1):
 728            for col in range(res - 1):
 729                pts[0] = col + row * res
 730                pts[1] = pts[0] + 1
 731                pts[2] = pts[0] + res + 1
 732                pts[3] = pts[0] + res
 733                polys.InsertNextCell(4, pts)
 734        surface.SetPoints(points)
 735        surface.SetPolys(polys)
 736        asurface = vedo.Mesh(surface)
 737        prop = vtk.vtkProperty()
 738        prop.DeepCopy(self.GetProperty())
 739        asurface.SetProperty(prop)
 740        asurface.property = prop
 741        asurface.lighting("default")
 742        self.points(self.points() + direction)
 743        return asurface
 744
 745    def reverse(self):
 746        """Reverse the points sequence order."""
 747        pts = np.flip(self.points(), axis=0)
 748        self.points(pts)
 749        return self
 750
 751
 752class DashedLine(Mesh):
 753    """
 754    Consider using `Line.pattern()` instead.
 755
 756    Build a dashed line segment between points `p0` and `p1`.
 757    If `p0` is a list of points returns the line connecting them.
 758    A 2D set of coords can also be passed as `p0=[x..], p1=[y..]`.
 759    """
 760    def __init__(self, p0, p1=None, spacing=0.1, closed=False, lw=2, c="k5", alpha=1):
 761        """
 762        Arguments:
 763            closed : (bool)
 764                join last to first point
 765            spacing : (float)
 766                relative size of the dash
 767            lw : (int)
 768                line width in pixels
 769        """
 770        if isinstance(p1, vtk.vtkActor):
 771            p1 = p1.GetPosition()
 772            if isinstance(p0, vtk.vtkActor):
 773                p0 = p0.GetPosition()
 774        if isinstance(p0, Points):
 775            p0 = p0.points()
 776
 777        # detect if user is passing a 2D list of points as p0=xlist, p1=ylist:
 778        if len(p0) > 3:
 779            if not utils.is_sequence(p0[0]) and not utils.is_sequence(p1[0]) and len(p0)==len(p1):
 780                # assume input is 2D xlist, ylist
 781                p0 = np.stack((p0, p1), axis=1)
 782                p1 = None
 783            p0 = utils.make3d(p0)
 784            if closed:
 785                p0 = np.append(p0, [p0[0]], axis=0)
 786
 787        if p1 is not None:  # assume passing p0=[x,y]
 788            if len(p0) == 2 and not utils.is_sequence(p0[0]):
 789                p0 = (p0[0], p0[1], 0)
 790            if len(p1) == 2 and not utils.is_sequence(p1[0]):
 791                p1 = (p1[0], p1[1], 0)
 792
 793        # detect if user is passing a list of points:
 794        if utils.is_sequence(p0[0]):
 795            listp = p0
 796        else:  # or just 2 points to link
 797            listp = [p0, p1]
 798
 799        listp = np.array(listp)
 800        if listp.shape[1] == 2:
 801            listp = np.c_[listp, np.zeros(listp.shape[0])]
 802
 803        xmn = np.min(listp, axis=0)
 804        xmx = np.max(listp, axis=0)
 805        dlen = np.linalg.norm(xmx - xmn) * np.clip(spacing, 0.01, 1.0) / 10
 806        if not dlen:
 807            Mesh.__init__(self, vtk.vtkPolyData(), c, alpha)
 808            self.name = "DashedLine (void)"
 809            return
 810
 811        qs = []
 812        for ipt in range(len(listp) - 1):
 813            p0 = listp[ipt]
 814            p1 = listp[ipt + 1]
 815            v = p1 - p0
 816            vdist = np.linalg.norm(v)
 817            n1 = int(vdist / dlen)
 818            if not n1:
 819                continue
 820
 821            res = 0
 822            for i in range(n1 + 2):
 823                ist = (i - 0.5) / n1
 824                ist = max(ist, 0)
 825                qi = p0 + v * (ist - res / vdist)
 826                if ist > 1:
 827                    qi = p1
 828                    res = np.linalg.norm(qi - p1)
 829                    qs.append(qi)
 830                    break
 831                qs.append(qi)
 832
 833        polylns = vtk.vtkAppendPolyData()
 834        for i, q1 in enumerate(qs):
 835            if not i % 2:
 836                continue
 837            q0 = qs[i - 1]
 838            line_source = vtk.vtkLineSource()
 839            line_source.SetPoint1(q0)
 840            line_source.SetPoint2(q1)
 841            line_source.Update()
 842            polylns.AddInputData(line_source.GetOutput())
 843        polylns.Update()
 844
 845        Mesh.__init__(self, polylns.GetOutput(), c, alpha)
 846        self.lw(lw).lighting("off")
 847        self.base = listp[0]
 848        if closed:
 849            self.top = listp[-2]
 850        else:
 851            self.top = listp[-1]
 852        self.name = "DashedLine"
 853
 854
 855class RoundedLine(Mesh):
 856    """
 857    Create a 2D line of specified thickness (in absolute units) passing through
 858    a list of input points. Borders of the line are rounded.
 859    """
 860    def __init__(self, pts, lw, res=10, c="gray4", alpha=1):
 861        """
 862        Arguments:
 863            pts : (list)
 864                a list of points in 2D or 3D (z will be ignored).
 865            lw : (float)
 866                thickness of the line.
 867            res : (int)
 868                resolution of the rounded regions
 869
 870        Example:
 871            ```python
 872            from vedo import *
 873            pts = [(-4,-3),(1,1),(2,4),(4,1),(3,-1),(2,-5),(9,-3)]
 874            ln = Line(pts, c='r', lw=2).z(0.01)
 875            rl = RoundedLine(pts, 0.6)
 876            show(Points(pts), ln, rl, axes=1).close()
 877            ```
 878            ![](https://vedo.embl.es/images/feats/rounded_line.png)
 879        """
 880        pts = utils.make3d(pts)
 881
 882        def _getpts(pts, revd=False):
 883
 884            if revd:
 885                pts = list(reversed(pts))
 886
 887            if len(pts) == 2:
 888                p0, p1 = pts
 889                v = p1 - p0
 890                dv = np.linalg.norm(v)
 891                nv = np.cross(v, (0, 0, -1))
 892                nv = nv / np.linalg.norm(nv) * lw
 893                return [p0 + nv, p1 + nv]
 894
 895            ptsnew = []
 896            for k in range(len(pts) - 2):
 897                p0 = pts[k]
 898                p1 = pts[k + 1]
 899                p2 = pts[k + 2]
 900                v = p1 - p0
 901                u = p2 - p1
 902                du = np.linalg.norm(u)
 903                dv = np.linalg.norm(v)
 904                nv = np.cross(v, (0, 0, -1))
 905                nv = nv / np.linalg.norm(nv) * lw
 906                nu = np.cross(u, (0, 0, -1))
 907                nu = nu / np.linalg.norm(nu) * lw
 908                uv = np.cross(u, v)
 909                if k == 0:
 910                    ptsnew.append(p0 + nv)
 911                if uv[2] <= 0:
 912                    alpha = np.arccos(np.dot(u, v) / du / dv)
 913                    db = lw * np.tan(alpha / 2)
 914                    p1new = p1 + nv - v / dv * db
 915                    ptsnew.append(p1new)
 916                else:
 917                    p1a = p1 + nv
 918                    p1b = p1 + nu
 919                    for i in range(0, res + 1):
 920                        pab = p1a * (res - i) / res + p1b * i / res
 921                        vpab = pab - p1
 922                        vpab = vpab / np.linalg.norm(vpab) * lw
 923                        ptsnew.append(p1 + vpab)
 924                if k == len(pts) - 3:
 925                    ptsnew.append(p2 + nu)
 926                    if revd:
 927                        ptsnew.append(p2 - nu)
 928            return ptsnew
 929
 930        ptsnew = _getpts(pts) + _getpts(pts, revd=True)
 931
 932        ppoints = vtk.vtkPoints()  # Generate the polyline
 933        ppoints.SetData(utils.numpy2vtk(np.asarray(ptsnew), dtype=np.float32))
 934        lines = vtk.vtkCellArray()
 935        npt = len(ptsnew)
 936        lines.InsertNextCell(npt)
 937        for i in range(npt):
 938            lines.InsertCellPoint(i)
 939        poly = vtk.vtkPolyData()
 940        poly.SetPoints(ppoints)
 941        poly.SetLines(lines)
 942        vct = vtk.vtkContourTriangulator()
 943        vct.SetInputData(poly)
 944        vct.Update()
 945        Mesh.__init__(self, vct.GetOutput(), c, alpha)
 946        self.flat()
 947        self.property.LightingOff()
 948        self.name = "RoundedLine"
 949        self.base = ptsnew[0]
 950        self.top = ptsnew[-1]
 951
 952
 953class Lines(Mesh):
 954    """
 955    Build the line segments between two lists of points `start_pts` and `end_pts`.
 956    `start_pts` can be also passed in the form `[[point1, point2], ...]`.
 957    """
 958    def __init__(
 959        self,
 960        start_pts,
 961        end_pts=None,
 962        dotted=False,
 963        res=1,
 964        scale=1,
 965        lw=1,
 966        c="k4",
 967        alpha=1,
 968    ):
 969        """
 970        Arguments:
 971            scale : (float)
 972                apply a rescaling factor to the lengths.
 973            c : (color, int, str, list)
 974                color name, number, or list of [R,G,B] colors
 975            alpha : (float)
 976                opacity in range [0,1]
 977            lw : (int)
 978                line width in pixel units
 979            res : (int)
 980                resolution, number of points along the line
 981                (only relevant if only 2 points are specified)
 982
 983        Examples:
 984            - [fitspheres2.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/fitspheres2.py)
 985
 986            ![](https://user-images.githubusercontent.com/32848391/52503049-ac9cb600-2be4-11e9-86af-72a538af14ef.png)
 987        """
 988        if isinstance(start_pts, Points):
 989            start_pts = start_pts.points()
 990        if isinstance(end_pts, Points):
 991            end_pts = end_pts.points()
 992
 993        if end_pts is not None:
 994            start_pts = np.stack((start_pts, end_pts), axis=1)
 995
 996        polylns = vtk.vtkAppendPolyData()
 997
 998        if not utils.is_ragged(start_pts):
 999
1000            for twopts in start_pts:
1001                line_source = vtk.vtkLineSource()
1002                line_source.SetResolution(res)
1003                if len(twopts[0]) == 2:
1004                    line_source.SetPoint1(twopts[0][0], twopts[0][1], 0.0)
1005                else:
1006                    line_source.SetPoint1(twopts[0])
1007
1008                if scale == 1:
1009                    pt2 = twopts[1]
1010                else:
1011                    vers = (np.array(twopts[1]) - twopts[0]) * scale
1012                    pt2 = np.array(twopts[0]) + vers
1013
1014                if len(pt2) == 2:
1015                    line_source.SetPoint2(pt2[0], pt2[1], 0.0)
1016                else:
1017                    line_source.SetPoint2(pt2)
1018                polylns.AddInputConnection(line_source.GetOutputPort())
1019
1020        else:
1021
1022            polylns = vtk.vtkAppendPolyData()
1023            for t in start_pts:
1024                t = utils.make3d(t)
1025                ppoints = vtk.vtkPoints()  # Generate the polyline
1026                ppoints.SetData(utils.numpy2vtk(t, dtype=np.float32))
1027                lines = vtk.vtkCellArray()
1028                npt = len(t)
1029                lines.InsertNextCell(npt)
1030                for i in range(npt):
1031                    lines.InsertCellPoint(i)
1032                poly = vtk.vtkPolyData()
1033                poly.SetPoints(ppoints)
1034                poly.SetLines(lines)
1035                polylns.AddInputData(poly)
1036
1037        polylns.Update()
1038
1039        Mesh.__init__(self, polylns.GetOutput(), c, alpha)
1040        self.lw(lw).lighting("off")
1041        if dotted:
1042            self.GetProperty().SetLineStipplePattern(0xF0F0)
1043            self.GetProperty().SetLineStippleRepeatFactor(1)
1044
1045        self.name = "Lines"
1046
1047
1048class Spline(Line):
1049    """
1050    Find the B-Spline curve through a set of points. This curve does not necessarily
1051    pass exactly through all the input points. Needs to import `scipy`.
1052    """
1053    def __init__(self, points, smooth=0, degree=2, closed=False, res=None, easing=""):
1054        """
1055        Arguments:
1056            smooth : (float)
1057                smoothing factor.
1058                - 0 = interpolate points exactly [default].
1059                - 1 = average point positions.
1060            degree : (int)
1061                degree of the spline (1<degree<5)
1062            easing : (str)
1063                control sensity of points along the spline.
1064                Available options are
1065                `[InSine, OutSine, Sine, InQuad, OutQuad, InCubic, OutCubic, InQuart, OutQuart, InCirc, OutCirc].`
1066                Can be used to create animations (move objects at varying speed).
1067                See e.g.: https://easings.net
1068            res : (int)
1069                number of points on the spline
1070
1071        See also: `CSpline` and `KSpline`.
1072
1073        Examples:
1074            - [spline_ease.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/spline_ease.py)
1075
1076                ![](https://vedo.embl.es/images/simulations/spline_ease.gif)
1077        """
1078        from scipy.interpolate import splprep, splev
1079
1080        if isinstance(points, Points):
1081            points = points.points()
1082
1083        points = utils.make3d(points)
1084
1085        per = 0
1086        if closed:
1087            points = np.append(points, [points[0]], axis=0)
1088            per = 1
1089
1090        if res is None:
1091            res = len(points) * 10
1092
1093        points = np.array(points, dtype=float)
1094
1095        minx, miny, minz = np.min(points, axis=0)
1096        maxx, maxy, maxz = np.max(points, axis=0)
1097        maxb = max(maxx - minx, maxy - miny, maxz - minz)
1098        smooth *= maxb / 2  # must be in absolute units
1099
1100        x = np.linspace(0, 1, res)
1101        if easing:
1102            if easing == "InSine":
1103                x = 1 - np.cos((x * np.pi) / 2)
1104            elif easing == "OutSine":
1105                x = np.sin((x * np.pi) / 2)
1106            elif easing == "Sine":
1107                x = -(np.cos(np.pi * x) - 1) / 2
1108            elif easing == "InQuad":
1109                x = x * x
1110            elif easing == "OutQuad":
1111                x = 1 - (1 - x) * (1 - x)
1112            elif easing == "InCubic":
1113                x = x * x
1114            elif easing == "OutCubic":
1115                x = 1 - np.power(1 - x, 3)
1116            elif easing == "InQuart":
1117                x = x * x * x * x
1118            elif easing == "OutQuart":
1119                x = 1 - np.power(1 - x, 4)
1120            elif easing == "InCirc":
1121                x = 1 - np.sqrt(1 - np.power(x, 2))
1122            elif easing == "OutCirc":
1123                x = np.sqrt(1 - np.power(x - 1, 2))
1124            else:
1125                vedo.logger.error(f"unknown ease mode {easing}")
1126
1127        # find the knots
1128        tckp, _ = splprep(points.T, task=0, s=smooth, k=degree, per=per)
1129        # evaluate spLine, including interpolated points:
1130        xnew, ynew, znew = splev(x, tckp)
1131
1132        Line.__init__(self, np.c_[xnew, ynew, znew], lw=2)
1133        self.name = "Spline"
1134
1135
1136class KSpline(Line):
1137    """
1138    Return a [Kochanek spline](https://en.wikipedia.org/wiki/Kochanek%E2%80%93Bartels_spline)
1139    which runs exactly through all the input points.
1140    """
1141    def __init__(self, points, continuity=0, tension=0, bias=0, closed=False, res=None):
1142        """
1143        Arguments:
1144            continuity : (float)
1145                changes the sharpness in change between tangents
1146            tension : (float)
1147                changes the length of the tangent vector
1148            bias : (float)
1149                changes the direction of the tangent vector
1150            closed : (bool)
1151                join last to first point to produce a closed curve
1152            res : (int)
1153                approximate resolution of the output line.
1154                Default is 20 times the number of input points.
1155
1156        ![](https://user-images.githubusercontent.com/32848391/65975805-73fd6580-e46f-11e9-8957-75eddb28fa72.png)
1157
1158        See also: `Spline` and `CSpline`.
1159        """
1160        if isinstance(points, Points):
1161            points = points.points()
1162
1163        if not res:
1164            res = len(points) * 20
1165
1166        points = utils.make3d(points).astype(float)
1167
1168        xspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkKochanekSpline()
1169        yspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkKochanekSpline()
1170        zspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkKochanekSpline()
1171        for s in [xspline, yspline, zspline]:
1172            if bias:
1173                s.SetDefaultBias(bias)
1174            if tension:
1175                s.SetDefaultTension(tension)
1176            if continuity:
1177                s.SetDefaultContinuity(continuity)
1178            s.SetClosed(closed)
1179
1180        lenp = len(points[0]) > 2
1181
1182        for i, p in enumerate(points):
1183            xspline.AddPoint(i, p[0])
1184            yspline.AddPoint(i, p[1])
1185            if lenp:
1186                zspline.AddPoint(i, p[2])
1187
1188        ln = []
1189        for pos in np.linspace(0, len(points), res):
1190            x = xspline.Evaluate(pos)
1191            y = yspline.Evaluate(pos)
1192            z = 0
1193            if lenp:
1194                z = zspline.Evaluate(pos)
1195            ln.append((x, y, z))
1196
1197        Line.__init__(self, ln, lw=2)
1198        self.clean()
1199        self.lighting("off")
1200        self.name = "KSpline"
1201        self.base = np.array(points[0], dtype=float)
1202        self.top = np.array(points[-1], dtype=float)
1203
1204
1205class CSpline(Line):
1206    """
1207    Return a Cardinal spline which runs exactly through all the input points.
1208    """
1209    def __init__(self, points, closed=False, res=None):
1210        """
1211        Arguments:
1212            closed : (bool)
1213                join last to first point to produce a closed curve
1214            res : (int)
1215                approximate resolution of the output line.
1216                Default is 20 times the number of input points.
1217
1218        See also: `Spline` and `KSpline`.
1219        """
1220
1221        if isinstance(points, Points):
1222            points = points.points()
1223
1224        if not res:
1225            res = len(points) * 20
1226
1227        points = utils.make3d(points).astype(float)
1228
1229        xspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkCardinalSpline()
1230        yspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkCardinalSpline()
1231        zspline = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkCardinalSpline()
1232        for s in [xspline, yspline, zspline]:
1233            s.SetClosed(closed)
1234
1235        lenp = len(points[0]) > 2
1236
1237        for i, p in enumerate(points):
1238            xspline.AddPoint(i, p[0])
1239            yspline.AddPoint(i, p[1])
1240            if lenp:
1241                zspline.AddPoint(i, p[2])
1242
1243        ln = []
1244        for pos in np.linspace(0, len(points), res):
1245            x = xspline.Evaluate(pos)
1246            y = yspline.Evaluate(pos)
1247            z = 0
1248            if lenp:
1249                z = zspline.Evaluate(pos)
1250            ln.append((x, y, z))
1251
1252        Line.__init__(self, ln, lw=2)
1253        self.clean()
1254        self.lighting("off")
1255        self.name = "CSpline"
1256        self.base = points[0]
1257        self.top = points[-1]
1258
1259
1260class Bezier(Line):
1261    """
1262    Generate the Bezier line that links the first to the last point.
1263    """
1264    def __init__(self, points, res=None):
1265        """
1266        Example:
1267            ```python
1268            from vedo import *
1269            import numpy as np
1270            pts = np.random.randn(25,3)
1271            for i,p in enumerate(pts):
1272                p += [5*i, 15*sin(i/2), i*i*i/200]
1273            show(Points(pts), Bezier(pts), axes=1).close()
1274            ```
1275            ![](https://user-images.githubusercontent.com/32848391/90437534-dafd2a80-e0d2-11ea-9b93-9ecb3f48a3ff.png)
1276        """
1277        N = len(points)
1278        if res is None:
1279            res = 10 * N
1280        t = np.linspace(0, 1, num=res)
1281        bcurve = np.zeros((res, len(points[0])))
1282
1283        def binom(n, k):
1284            b = 1
1285            for t in range(1, min(k, n - k) + 1):
1286                b *= n / t
1287                n -= 1
1288            return b
1289
1290        def bernstein(n, k):
1291            coeff = binom(n, k)
1292
1293            def _bpoly(x):
1294                return coeff * x ** k * (1 - x) ** (n - k)
1295
1296            return _bpoly
1297
1298        for ii in range(N):
1299            b = bernstein(N - 1, ii)(t)
1300            bcurve += np.outer(b, points[ii])
1301        Line.__init__(self, bcurve, lw=2)
1302        self.name = "BezierLine"
1303
1304
1305class NormalLines(Mesh):
1306    """
1307    Build an `Glyph` made of the normals at cells shown as lines.
1308    """
1309    def __init__(self, mesh, ratio=1, on="cells", scale=1):
1310        """
1311        If `on="points"` normals are shown at mesh vertices.
1312        """
1313        poly = mesh.clone().compute_normals().polydata()
1314
1315        if "cell" in on:
1316            centers = vtk.vtkCellCenters()
1317            centers.SetInputData(poly)
1318            centers.Update()
1319            poly = centers.GetOutput()
1320
1321        mask_pts = vtk.vtkMaskPoints()
1322        mask_pts.SetInputData(poly)
1323        mask_pts.SetOnRatio(ratio)
1324        mask_pts.RandomModeOff()
1325        mask_pts.Update()
1326
1327        ln = vtk.vtkLineSource()
1328        ln.SetPoint1(0, 0, 0)
1329        ln.SetPoint2(1, 0, 0)
1330        ln.Update()
1331        glyph = vtk.vtkGlyph3D()
1332        glyph.SetSourceData(ln.GetOutput())
1333        glyph.SetInputData(mask_pts.GetOutput())
1334        glyph.SetVectorModeToUseNormal()
1335
1336        b = poly.GetBounds()
1337        sc = max([b[1] - b[0], b[3] - b[2], b[5] - b[4]]) / 50 * scale
1338        glyph.SetScaleFactor(sc)
1339        glyph.OrientOn()
1340        glyph.Update()
1341        Mesh.__init__(self, glyph.GetOutput())
1342        self.PickableOff()
1343        self.SetProperty(mesh.GetProperty())
1344        self.property = mesh.GetProperty()
1345        self.property.LightingOff()
1346        self.name = "NormalLines"
1347
1348
1349def _interpolate2vol(mesh, kernel=None, radius=None, bounds=None, null_value=None, dims=None):
1350    # Generate a volumetric dataset by interpolating a scalar
1351    # or vector field which is only known on a scattered set of points or mesh.
1352    # Available interpolation kernels are: shepard, gaussian, voronoi, linear.
1353    #
1354    # kernel : (str) interpolation kernel type [shepard]
1355    # radius : (float) radius of the local search
1356    # bounds : (list) bounding box of the output object
1357    # dims : (list) dimensions of the output object
1358    # null_value : (float) value to be assigned to invalid points
1359    if dims is None:
1360        dims = (25, 25, 25)
1361
1362    if bounds is None:
1363        bounds = mesh.bounds()
1364    elif isinstance(bounds, vedo.base.Base3DProp):
1365        bounds = bounds.bounds()
1366
1367    # Create a domain volume
1368    domain = vtk.vtkImageData()
1369    domain.SetDimensions(dims)
1370    domain.SetOrigin(bounds[0], bounds[2], bounds[4])
1371    deltaZ = (bounds[5] - bounds[4]) / (dims[2] - 1)
1372    deltaY = (bounds[3] - bounds[2]) / (dims[1] - 1)
1373    deltaX = (bounds[1] - bounds[0]) / (dims[0] - 1)
1374    domain.SetSpacing(deltaX, deltaY, deltaZ)
1375
1376    if radius is None:
1377        radius = 2.5*np.sqrt(deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ)
1378
1379    locator = vtk.vtkStaticPointLocator()
1380    locator.SetDataSet(mesh)
1381    locator.BuildLocator()
1382
1383    if kernel == "gaussian":
1384        kern = vtk.vtkGaussianKernel()
1385        kern.SetRadius(radius)
1386    elif kernel == "voronoi":
1387        kern = vtk.vtkVoronoiKernel()
1388    elif kernel == "linear":
1389        kern = vtk.vtkLinearKernel()
1390        kern.SetRadius(radius)
1391    else:
1392        kern = vtk.vtkShepardKernel()
1393        kern.SetPowerParameter(2)
1394        kern.SetRadius(radius)
1395
1396    interpolator = vtk.vtkPointInterpolator()
1397    interpolator.SetInputData(domain)
1398    interpolator.SetSourceData(mesh)
1399    interpolator.SetKernel(kern)
1400    interpolator.SetLocator(locator)
1401    if null_value is not None:
1402        interpolator.SetNullValue(null_value)
1403    else:
1404        interpolator.SetNullPointsStrategyToClosestPoint()
1405    interpolator.Update()
1406    return interpolator.GetOutput()
1407
1408
1409def StreamLines(
1410    domain,
1411    probe,
1412    active_vectors="",
1413    integrator="rk4",
1414    direction="forward",
1415    initial_step_size=None,
1416    max_propagation=None,
1417    max_steps=10000,
1418    step_length=None,
1419    extrapolate_to_box=(),
1420    surface_constrained=False,
1421    compute_vorticity=False,
1422    ribbons=None,
1423    tubes=(),
1424    scalar_range=None,
1425    lw=None,
1426    **opts,
1427):
1428    """
1429    Integrate a vector field on a domain (a Points/Mesh or other vtk datasets types)
1430    to generate streamlines.
1431
1432    The integration is performed using a specified integrator (Runge-Kutta).
1433    The length of a streamline is governed by specifying a maximum value either
1434    in physical arc length or in (local) cell length.
1435    Otherwise, the integration terminates upon exiting the field domain.
1436
1437    Arguments:
1438        domain : (Points, Volume, vtkDataSet)
1439            the object that contains the vector field
1440        probe : (Mesh, list)
1441            the Mesh that probes the domain. Its coordinates will
1442            be the seeds for the streamlines, can also be an array of positions.
1443        active_vectors : (str)
1444            name of the vector array to be used
1445        integrator : (str)
1446            Runge-Kutta integrator, either 'rk2', 'rk4' of 'rk45'
1447        initial_step_size : (float)
1448            initial step size of integration
1449        max_propagation : (float)
1450            maximum physical length of the streamline
1451        max_steps : (int)
1452            maximum nr of steps allowed
1453        step_length : (float)
1454            length of step integration.
1455        extrapolate_to_box : (dict)
1456            Vectors that are defined on a descrete set of points
1457            are extrapolated to a 3D domain defined by its bounding box:
1458                - bounds (list), bounding box of the domain
1459                - kernel (str), interpolation kernel `["shepard","gaussian","voronoi","linear"]`
1460                - radius (float), radius of the local search
1461                - dims (list), nr of subdivisions of the domain along x, y, and z
1462                - null_value (float), value to be assigned to invalid points
1463
1464        surface_constrained : (bool)
1465            force streamlines to be computed on a surface
1466        compute_vorticity : (bool)
1467            Turn on/off vorticity computation at streamline points
1468            (necessary for generating proper stream-ribbons)
1469        ribbons : (int)
1470            render lines as ribbons by joining them.
1471            An integer value represent the ratio of joining (e.g.: ribbons=2 groups lines 2 by 2)
1472
1473        tubes : (dict)
1474            dictionary containing the parameters for the tube representation:
1475            - ratio (int), draws tube as longitudinal stripes
1476            - res (int), tube resolution (nr. of sides, 12 by default)
1477            - max_radius_factor (float), max tube radius as a multiple of the min radius
1478            - cap (bool), capping of the tube
1479            - mode (int), radius varies based on the scalar or vector magnitude:
1480                - 0 do not vary radius
1481                - 1 vary radius by scalar
1482                - 2 vary radius by vector
1483                - 3 vary radius by absolute value of scalar
1484                - 4 vary radius by vector norm
1485
1486        scalar_range : (list)
1487            specify the scalar range for coloring
1488
1489    Examples:
1490        - [streamlines1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/streamlines1.py)
1491        - [streamlines2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/streamlines2.py)
1492
1493            ![](https://vedo.embl.es/images/volumetric/81459343-b9210d00-919f-11ea-846c-152d62cba06e.png)
1494
1495        - [streamribbons.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/streamribbons.py)
1496        - [office.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/office.py)
1497
1498            ![](https://vedo.embl.es/images/volumetric/56964003-9145a500-6b5a-11e9-9d9e-9736d90e1900.png)
1499    """
1500    if len(opts): # Deprecations
1501        printc(" Warning! In StreamLines() unrecognized keywords:", opts, c='y')
1502        initial_step_size = opts.pop("initialStepSize", initial_step_size)
1503        max_propagation = opts.pop("maxPropagation", max_propagation)
1504        max_steps = opts.pop("maxSteps", max_steps)
1505        step_length = opts.pop("stepLength", step_length)
1506        extrapolate_to_box = opts.pop("extrapolateToBox", extrapolate_to_box)
1507        surface_constrained = opts.pop("surfaceConstrained", surface_constrained)
1508        compute_vorticity = opts.pop("computeVorticity", compute_vorticity)
1509        scalar_range = opts.pop("scalarRange", scalar_range)
1510        printc("          Please use 'snake_case' instead of 'camelCase' keywords", c='y')
1511
1512    if isinstance(domain, vedo.Points):
1513        if extrapolate_to_box:
1514            grid = _interpolate2vol(domain.polydata(), **extrapolate_to_box)
1515        else:
1516            grid = domain.polydata()
1517    elif isinstance(domain, vedo.BaseVolume):
1518        grid = domain.inputdata()
1519    else:
1520        grid = domain
1521
1522    if active_vectors:
1523        grid.GetPointData().SetActiveVectors(active_vectors)
1524
1525    b = grid.GetBounds()
1526    size = (b[5] - b[4] + b[3] - b[2] + b[1] - b[0]) / 3
1527    if initial_step_size is None:
1528        initial_step_size = size / 500.0
1529    if max_propagation is None:
1530        max_propagation = size
1531
1532    if utils.is_sequence(probe):
1533        pts = utils.make3d(probe)
1534    else:
1535        pts = probe.clean().points()
1536
1537    src = vtk.vtkProgrammableSource()
1538
1539    def read_points():
1540        output = src.GetPolyDataOutput()
1541        points = vtk.vtkPoints()
1542        for x, y, z in pts:
1543            points.InsertNextPoint(x, y, z)
1544        output.SetPoints(points)
1545
1546    src.SetExecuteMethod(read_points)
1547    src.Update()
1548
1549    st = vtk.vtkStreamTracer()
1550    st.SetInputDataObject(grid)
1551    st.SetSourceConnection(src.GetOutputPort())
1552
1553    st.SetInitialIntegrationStep(initial_step_size)
1554    st.SetComputeVorticity(compute_vorticity)
1555    st.SetMaximumNumberOfSteps(max_steps)
1556    st.SetMaximumPropagation(max_propagation)
1557    st.SetSurfaceStreamlines(surface_constrained)
1558    if step_length:
1559        st.SetMaximumIntegrationStep(step_length)
1560
1561    if "f" in direction:
1562        st.SetIntegrationDirectionToForward()
1563    elif "back" in direction:
1564        st.SetIntegrationDirectionToBackward()
1565    elif "both" in direction:
1566        st.SetIntegrationDirectionToBoth()
1567
1568    if integrator == "rk2":
1569        st.SetIntegratorTypeToRungeKutta2()
1570    elif integrator == "rk4":
1571        st.SetIntegratorTypeToRungeKutta4()
1572    elif integrator == "rk45":
1573        st.SetIntegratorTypeToRungeKutta45()
1574    else:
1575        vedo.logger.error(f"in streamlines, unknown integrator {integrator}")
1576
1577    st.Update()
1578    output = st.GetOutput()
1579
1580    if ribbons:
1581        scalar_surface = vtk.vtkRuledSurfaceFilter()
1582        scalar_surface.SetInputConnection(st.GetOutputPort())
1583        scalar_surface.SetOnRatio(int(ribbons))
1584        scalar_surface.SetRuledModeToPointWalk()
1585        scalar_surface.Update()
1586        output = scalar_surface.GetOutput()
1587
1588    if tubes:
1589        radius = tubes.pop("radius", domain.GetLength()/500)
1590        res = tubes.pop("res", 24)
1591        radfact = tubes.pop("max_radius_factor", 10)
1592        ratio = tubes.pop("ratio", 1)
1593        mode = tubes.pop("mode", 0)
1594        cap = tubes.pop("mode", False)
1595        if tubes:
1596            vedo.logger.warning(f"in StreamLines unknown 'tubes' parameters: {tubes}")
1597
1598        stream_tube = vtk.vtkTubeFilter()
1599        stream_tube.SetNumberOfSides(res)
1600        stream_tube.SetRadius(radius)
1601        stream_tube.SetCapping(cap)
1602        # max tube radius as a multiple of the min radius
1603        stream_tube.SetRadiusFactor(radfact)
1604        stream_tube.SetOnRatio(int(ratio))
1605        stream_tube.SetVaryRadius(int(mode))
1606
1607        stream_tube.SetInputData(output)
1608        vname = grid.GetPointData().GetVectors().GetName()
1609        stream_tube.SetInputArrayToProcess(
1610            1, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, vname
1611        )
1612        stream_tube.Update()
1613        sta = vedo.mesh.Mesh(stream_tube.GetOutput(), c=None)
1614
1615        scals = grid.GetPointData().GetScalars()
1616        if scals:
1617            sta.mapper().SetScalarRange(scals.GetRange())
1618        if scalar_range is not None:
1619            sta.mapper().SetScalarRange(scalar_range)
1620
1621        sta.phong()
1622        sta.name = "StreamLines"
1623        #############
1624        return sta  #############
1625        #############
1626
1627    sta = vedo.mesh.Mesh(output, c=None)
1628
1629    if lw is not None and len(tubes) == 0 and not ribbons:
1630        sta.lw(lw)
1631        sta.mapper().SetResolveCoincidentTopologyToPolygonOffset()
1632        sta.lighting("off")
1633
1634    scals = grid.GetPointData().GetScalars()
1635    if scals:
1636        sta.mapper().SetScalarRange(scals.GetRange())
1637    if scalar_range is not None:
1638        sta.mapper().SetScalarRange(scalar_range)
1639
1640    sta.name = "StreamLines"
1641    return sta
1642
1643
1644class Tube(Mesh):
1645    """
1646    Build a tube along the line defined by a set of points.
1647    """
1648    def __init__(self, points, r=1, cap=True, res=12, c=None, alpha=1):
1649        """
1650        Arguments:
1651            r :  (float, list)
1652                constant radius or list of radii.
1653            res : (int)
1654                resolution, number of the sides of the tube
1655            c : (color)
1656                constant color or list of colors for each point.
1657
1658        Examples:
1659            - [ribbon.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/ribbon.py)
1660            - [tube_radii.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/tube_radii.py)
1661
1662                ![](https://vedo.embl.es/images/basic/tube.png)
1663        """
1664        if isinstance(points, vedo.Points):
1665            points = points.points()
1666
1667        base = np.asarray(points[0], dtype=float)
1668        top = np.asarray(points[-1], dtype=float)
1669
1670        vpoints = vtk.vtkPoints()
1671        idx = len(points)
1672        for p in points:
1673            vpoints.InsertNextPoint(p)
1674        line = vtk.vtkPolyLine()
1675        line.GetPointIds().SetNumberOfIds(idx)
1676        for i in range(idx):
1677            line.GetPointIds().SetId(i, i)
1678        lines = vtk.vtkCellArray()
1679        lines.InsertNextCell(line)
1680        polyln = vtk.vtkPolyData()
1681        polyln.SetPoints(vpoints)
1682        polyln.SetLines(lines)
1683
1684        tuf = vtk.vtkTubeFilter()
1685        tuf.SetCapping(cap)
1686        tuf.SetNumberOfSides(res)
1687        tuf.SetInputData(polyln)
1688        if utils.is_sequence(r):
1689            arr = utils.numpy2vtk(r, dtype=float)
1690            arr.SetName("TubeRadius")
1691            polyln.GetPointData().AddArray(arr)
1692            polyln.GetPointData().SetActiveScalars("TubeRadius")
1693            tuf.SetVaryRadiusToVaryRadiusByAbsoluteScalar()
1694        else:
1695            tuf.SetRadius(r)
1696
1697        usingColScals = False
1698        if utils.is_sequence(c):
1699            usingColScals = True
1700            cc = vtk.vtkUnsignedCharArray()
1701            cc.SetName("TubeColors")
1702            cc.SetNumberOfComponents(3)
1703            cc.SetNumberOfTuples(len(c))
1704            for i, ic in enumerate(c):
1705                r, g, b = get_color(ic)
1706                cc.InsertTuple3(i, int(255 * r), int(255 * g), int(255 * b))
1707            polyln.GetPointData().AddArray(cc)
1708            c = None
1709        tuf.Update()
1710
1711        Mesh.__init__(self, tuf.GetOutput(), c, alpha)
1712        self.phong()
1713        if usingColScals:
1714            self.mapper().SetScalarModeToUsePointFieldData()
1715            self.mapper().ScalarVisibilityOn()
1716            self.mapper().SelectColorArray("TubeColors")
1717            self.mapper().Modified()
1718
1719        self.base = base
1720        self.top = top
1721        self.name = "Tube"
1722
1723def ThickTube(pts, r1, r2, res=12, c=None, alpha=1):
1724    """
1725    Create a tube with a thickness along a line of points.
1726
1727    Example:
1728    ```python
1729    from vedo import *
1730    pts = [[sin(x), cos(x), x/3] for x in np.arange(0.1, 3, 0.3)]
1731    vline = Line(pts, lw=5, c='red5')
1732    thick_tube = ThickTube(vline, r1=0.2, r2=0.3).lw(1)
1733    show(vline, thick_tube, axes=1).close()
1734    ```
1735    ![](https://vedo.embl.es/images/feats/thick_tube.png)
1736    """
1737    def make_cap(t1, t2):
1738        newpoints = t1.points().tolist() + t2.points().tolist()
1739        newfaces = []
1740        for i in range(n-1):
1741            newfaces.append([i,   i+1, i+n])
1742            newfaces.append([i+n, i+1, i+n+1])
1743        newfaces.append([2*n-1,   0, n])
1744        newfaces.append([2*n-1, n-1, 0])
1745        capm = utils.buildPolyData(newpoints, newfaces)
1746        return capm
1747    
1748    assert r1 < r2
1749
1750    t1 = Tube(pts, r=r1, cap=False, res=res)
1751    t2 = Tube(pts, r=r2, cap=False, res=res)
1752
1753    tc1a, tc1b = t1.boundaries().split()
1754    tc2a, tc2b = t2.boundaries().split()
1755    n = tc1b.npoints
1756
1757    tc1b.join(reset=True).clean() # needed because indices are flipped
1758    tc2b.join(reset=True).clean()
1759
1760    capa = make_cap(tc1a, tc2a)
1761    capb = make_cap(tc1b, tc2b)
1762
1763    thick_tube = merge(t1, t2, capa, capb).c(c).alpha(alpha)
1764    thick_tube.base = t1.base
1765    thick_tube.top = t1.top
1766    thick_tube.name = "ThickTube"
1767    return thick_tube
1768
1769class Ribbon(Mesh):
1770    """
1771    Connect two lines to generate the surface inbetween.
1772    Set the mode by which to create the ruled surface.
1773
1774    It also works with a single line in input. In this case the ribbon
1775    is formed by following the local plane of the line in space.
1776    """
1777    def __init__(
1778        self,
1779        line1,
1780        line2=None,
1781        mode=0,
1782        closed=False,
1783        width=None,
1784        res=(200, 5),
1785        c="indigo3",
1786        alpha=1,
1787    ):
1788        """
1789        Arguments:
1790            mode : (int)
1791                If mode=0, resample evenly the input lines (based on length)
1792                and generates triangle strips.
1793
1794                If mode=1, use the existing points and walks around the
1795                polyline using existing points.
1796
1797            closed : (bool)
1798                if True, join the last point with the first to form a closed surface
1799
1800            res : (list)
1801                ribbon resolutions along the line and perpendicularly to it.
1802
1803        Examples:
1804            - [ribbon.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/ribbon.py)
1805
1806                ![](https://vedo.embl.es/images/basic/ribbon.png)
1807        """
1808
1809        if isinstance(line1, Points):
1810            line1 = line1.points()
1811
1812        if isinstance(line2, Points):
1813            line2 = line2.points()
1814
1815        elif line2 is None:
1816            #############################################
1817            ribbon_filter = vtk.vtkRibbonFilter()
1818            aline = Line(line1)
1819            ribbon_filter.SetInputData(aline.polydata())
1820            if width is None:
1821                width = aline.diagonal_size() / 20.0
1822            ribbon_filter.SetWidth(width)
1823            ribbon_filter.Update()
1824            Mesh.__init__(self, ribbon_filter.GetOutput(), c, alpha)
1825            self.name = "Ribbon"
1826            ##############################################
1827            return  ######################################
1828            ##############################################
1829
1830        line1 = np.asarray(line1)
1831        line2 = np.asarray(line2)
1832
1833        if closed:
1834            line1 = line1.tolist()
1835            line1 += [line1[0]]
1836            line2 = line2.tolist()
1837            line2 += [line2[0]]
1838            line1 = np.array(line1)
1839            line2 = np.array(line2)
1840
1841        if len(line1[0]) == 2:
1842            line1 = np.c_[line1, np.zeros(len(line1))]
1843        if len(line2[0]) == 2:
1844            line2 = np.c_[line2, np.zeros(len(line2))]
1845
1846        ppoints1 = vtk.vtkPoints()  # Generate the polyline1
1847        ppoints1.SetData(utils.numpy2vtk(line1, dtype=np.float32))
1848        lines1 = vtk.vtkCellArray()
1849        lines1.InsertNextCell(len(line1))
1850        for i in range(len(line1)):
1851            lines1.InsertCellPoint(i)
1852        poly1 = vtk.vtkPolyData()
1853        poly1.SetPoints(ppoints1)
1854        poly1.SetLines(lines1)
1855
1856        ppoints2 = vtk.vtkPoints()  # Generate the polyline2
1857        ppoints2.SetData(utils.numpy2vtk(line2, dtype=np.float32))
1858        lines2 = vtk.vtkCellArray()
1859        lines2.InsertNextCell(len(line2))
1860        for i in range(len(line2)):
1861            lines2.InsertCellPoint(i)
1862        poly2 = vtk.vtkPolyData()
1863        poly2.SetPoints(ppoints2)
1864        poly2.SetLines(lines2)
1865
1866        # build the lines
1867        lines1 = vtk.vtkCellArray()
1868        lines1.InsertNextCell(poly1.GetNumberOfPoints())
1869        for i in range(poly1.GetNumberOfPoints()):
1870            lines1.InsertCellPoint(i)
1871
1872        polygon1 = vtk.vtkPolyData()
1873        polygon1.SetPoints(ppoints1)
1874        polygon1.SetLines(lines1)
1875
1876        lines2 = vtk.vtkCellArray()
1877        lines2.InsertNextCell(poly2.GetNumberOfPoints())
1878        for i in range(poly2.GetNumberOfPoints()):
1879            lines2.InsertCellPoint(i)
1880
1881        polygon2 = vtk.vtkPolyData()
1882        polygon2.SetPoints(ppoints2)
1883        polygon2.SetLines(lines2)
1884
1885        merged_pd = vtk.vtkAppendPolyData()
1886        merged_pd.AddInputData(polygon1)
1887        merged_pd.AddInputData(polygon2)
1888        merged_pd.Update()
1889
1890        rsf = vtk.vtkRuledSurfaceFilter()
1891        rsf.CloseSurfaceOff()
1892        rsf.SetRuledMode(mode)
1893        rsf.SetResolution(res[0], res[1])
1894        rsf.SetInputData(merged_pd.GetOutput())
1895        rsf.Update()
1896
1897        Mesh.__init__(self, rsf.GetOutput(), c, alpha)
1898        self.name = "Ribbon"
1899
1900
1901class Arrow(Mesh):
1902    """
1903    Build a 3D arrow from `start_pt` to `end_pt` of section size `s`,
1904    expressed as the fraction of the window size.
1905    """
1906    def __init__(
1907        self,
1908        start_pt=(0, 0, 0),
1909        end_pt=(1, 0, 0),
1910        s=None,
1911        shaft_radius=None,
1912        head_radius=None,
1913        head_length=None,
1914        res=12,
1915        c="r4",
1916        alpha=1,
1917    ):
1918        """
1919        If `c` is a `float` less than 1, the arrow is rendered as a in a color scale
1920        from white to red.
1921
1922        .. note:: If `s=None` the arrow is scaled proportionally to its length
1923
1924        ![](https://raw.githubusercontent.com/lorensen/VTKExamples/master/src/Testing/Baseline/Cxx/GeometricObjects/TestOrientedArrow.png)
1925        """
1926        # in case user is passing meshs
1927        if isinstance(start_pt, vtk.vtkActor):
1928            start_pt = start_pt.GetPosition()
1929        if isinstance(end_pt, vtk.vtkActor):
1930            end_pt = end_pt.GetPosition()
1931
1932        axis = np.asarray(end_pt) - np.asarray(start_pt)
1933        length = np.linalg.norm(axis)
1934        if length:
1935            axis = axis / length
1936        if len(axis) < 3:  # its 2d
1937            theta = np.pi / 2
1938            start_pt = [start_pt[0], start_pt[1], 0.0]
1939            end_pt = [end_pt[0], end_pt[1], 0.0]
1940        else:
1941            theta = np.arccos(axis[2])
1942        phi = np.arctan2(axis[1], axis[0])
1943        self.source = vtk.vtkArrowSource()
1944        self.source.SetShaftResolution(res)
1945        self.source.SetTipResolution(res)
1946
1947        if s:
1948            sz = 0.02
1949            self.source.SetTipRadius(sz)
1950            self.source.SetShaftRadius(sz / 1.75)
1951            self.source.SetTipLength(sz * 15)
1952
1953        # if s:
1954        #     sz = 0.02 * s * length
1955        #     tl = sz / 20
1956        #     print(s, sz)
1957        #     self.source.SetShaftRadius(sz)
1958        #     self.source.SetTipRadius(sz*1.75)
1959        #     self.source.SetTipLength(sz*15)
1960
1961        if head_length:
1962            self.source.SetTipLength(head_length)
1963        if head_radius:
1964            self.source.SetTipRadius(head_radius)
1965        if shaft_radius:
1966            self.source.SetShaftRadius(shaft_radius)
1967
1968        self.source.Update()
1969
1970        t = vtk.vtkTransform()
1971        t.RotateZ(np.rad2deg(phi))
1972        t.RotateY(np.rad2deg(theta))
1973        t.RotateY(-90)  # put it along Z
1974        if s:
1975            sz = 800 * s
1976            t.Scale(length, sz, sz)
1977        else:
1978            t.Scale(length, length, length)
1979        tf = vtk.vtkTransformPolyDataFilter()
1980        tf.SetInputData(self.source.GetOutput())
1981        tf.SetTransform(t)
1982        tf.Update()
1983
1984        Mesh.__init__(self, tf.GetOutput(), c, alpha)
1985
1986        self.phong().lighting("plastic")
1987        self.SetPosition(start_pt)
1988        self.PickableOff()
1989        self.DragableOff()
1990        self.base = np.array(start_pt, dtype=float)
1991        self.top = np.array(end_pt, dtype=float)
1992        self.tip_index = None
1993        self.fill = True  # used by pyplot.__iadd__()
1994        self.s = s if s is not None else 1  ## used by pyplot.__iadd()
1995        self.name = "Arrow"
1996
1997    def tip_point(self, return_index=False):
1998        """Return the coordinates of the tip of the Arrow, or the point index."""
1999        if self.tip_index is None:
2000            arrpts = utils.vtk2numpy(self.source.GetOutput().GetPoints().GetData())
2001            self.tip_index = np.argmax(arrpts[:, 0])
2002        if return_index:
2003            return self.tip_index
2004        return self.points()[self.tip_index]
2005
2006
2007class Arrows(Glyph):
2008    """
2009    Build arrows between two lists of points.
2010    """
2011    def __init__(
2012        self,
2013        start_pts,
2014        end_pts=None,
2015        s=None,
2016        shaft_radius=None,
2017        head_radius=None,
2018        head_length=None,
2019        thickness=1,
2020        res=12,
2021        c=None,
2022        alpha=1,
2023    ):
2024        """
2025        Build arrows between two lists of points `start_pts` and `end_pts`.
2026         `start_pts` can be also passed in the form `[[point1, point2], ...]`.
2027
2028        Color can be specified as a colormap which maps the size of the arrows.
2029
2030        Arguments:
2031            s : (float)
2032                fix aspect-ratio of the arrow and scale its cross section
2033            c : (color)
2034                color or color map name
2035            alpha : (float)
2036                set object opacity
2037            res : (int)
2038                set arrow resolution
2039
2040        Examples:
2041            - [glyphs_arrows.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/glyphs_arrows.py)
2042
2043            ![](https://user-images.githubusercontent.com/32848391/55897850-a1a0da80-5bc1-11e9-81e0-004c8f396b43.jpg)
2044        """
2045        if isinstance(start_pts, Points):
2046            start_pts = start_pts.points()
2047        if isinstance(end_pts, Points):
2048            end_pts = end_pts.points()
2049
2050        start_pts = np.asarray(start_pts)
2051        if end_pts is None:
2052            strt = start_pts[:, 0]
2053            end_pts = start_pts[:, 1]
2054            start_pts = strt
2055        else:
2056            end_pts = np.asarray(end_pts)
2057
2058        start_pts = utils.make3d(start_pts)
2059        end_pts = utils.make3d(end_pts)
2060
2061        arr = vtk.vtkArrowSource()
2062        arr.SetShaftResolution(res)
2063        arr.SetTipResolution(res)
2064
2065        if s:
2066            sz = 0.02 * s
2067            arr.SetTipRadius(sz * 2)
2068            arr.SetShaftRadius(sz * thickness)
2069            arr.SetTipLength(sz * 10)
2070
2071        if head_radius:
2072            arr.SetTipRadius(head_radius)
2073        if shaft_radius:
2074            arr.SetShaftRadius(shaft_radius)
2075        if head_length:
2076            arr.SetTipLength(head_length)
2077
2078        arr.Update()
2079        out = arr.GetOutput()
2080
2081        orients = end_pts - start_pts
2082        Glyph.__init__(
2083            self,
2084            start_pts,
2085            out,
2086            orientation_array=orients,
2087            scale_by_vector_size=True,
2088            color_by_vector_size=True,
2089            c=c,
2090            alpha=alpha,
2091        )
2092        self.flat().lighting("plastic")
2093        self.name = "Arrows"
2094
2095
2096class Arrow2D(Mesh):
2097    """
2098    Build a 2D arrow.
2099    """
2100    def __init__(
2101        self,
2102        start_pt=(0, 0, 0),
2103        end_pt=(1, 0, 0),
2104        s=1,
2105        shaft_length=0.8,
2106        shaft_width=0.05,
2107        head_length=0.225,
2108        head_width=0.175,
2109        fill=True,
2110    ):
2111        """
2112        Build a 2D arrow from `start_pt` to `end_pt`.
2113
2114        Arguments:
2115            s : (float)
2116                a global multiplicative convenience factor controlling the arrow size
2117            shaft_length : (float)
2118                fractional shaft length
2119            shaft_width : (float)
2120                fractional shaft width
2121            head_length : (float)
2122                fractional head length
2123            head_width : (float)
2124                fractional head width
2125            fill : (bool)
2126                if False only generate the outline
2127        """
2128        self.fill = fill  ## needed by pyplot.__iadd()
2129        self.s = s  #  # needed by pyplot.__iadd()
2130
2131        if s != 1:
2132            shaft_width *= s
2133            head_width *= np.sqrt(s)
2134
2135        # in case user is passing meshs
2136        if isinstance(start_pt, vtk.vtkActor):
2137            start_pt = start_pt.GetPosition()
2138        if isinstance(end_pt, vtk.vtkActor):
2139            end_pt = end_pt.GetPosition()
2140        if len(start_pt) == 2:
2141            start_pt = [start_pt[0], start_pt[1], 0]
2142        if len(end_pt) == 2:
2143            end_pt = [end_pt[0], end_pt[1], 0]
2144
2145        headBase = 1 - head_length
2146        head_width = max(head_width, shaft_width)
2147        if head_length is None or headBase > shaft_length:
2148            headBase = shaft_length
2149
2150        verts = []
2151        verts.append([0, -shaft_width / 2, 0])
2152        verts.append([shaft_length, -shaft_width / 2, 0])
2153        verts.append([headBase, -head_width / 2, 0])
2154        verts.append([1, 0, 0])
2155        verts.append([headBase, head_width / 2, 0])
2156        verts.append([shaft_length, shaft_width / 2, 0])
2157        verts.append([0, shaft_width / 2, 0])
2158        if fill:
2159            faces = ((0, 1, 3, 5, 6), (5, 3, 4), (1, 2, 3))
2160            poly = utils.buildPolyData(verts, faces)
2161        else:
2162            lines = (0, 1, 2, 3, 4, 5, 6, 0)
2163            poly = utils.buildPolyData(verts, [], lines=lines)
2164
2165        axis = np.array(end_pt) - np.array(start_pt)
2166        length = np.linalg.norm(axis)
2167        if length:
2168            axis = axis / length
2169        theta = 0
2170        if len(axis) > 2:
2171            theta = np.arccos(axis[2])
2172        phi = np.arctan2(axis[1], axis[0])
2173        t = vtk.vtkTransform()
2174        if phi:
2175            t.RotateZ(np.rad2deg(phi))
2176        if theta:
2177            t.RotateY(np.rad2deg(theta))
2178        t.RotateY(-90)  # put it along Z
2179        t.Scale(length, length, length)
2180        tf = vtk.vtkTransformPolyDataFilter()
2181        tf.SetInputData(poly)
2182        tf.SetTransform(t)
2183        tf.Update()
2184
2185        Mesh.__init__(self, tf.GetOutput(), c='k1')
2186        self.SetPosition(start_pt)
2187        self.lighting("off")
2188        self.DragableOff()
2189        self.PickableOff()
2190        self.base = np.array(start_pt, dtype=float)
2191        self.top = np.array(end_pt, dtype=float)
2192        self.name = "Arrow2D"
2193
2194
2195class Arrows2D(Glyph):
2196    """
2197    Build 2D arrows between two lists of points.
2198    """
2199    def __init__(
2200        self,
2201        start_pts,
2202        end_pts=None,
2203        s=1,
2204        shaft_length=0.8,
2205        shaft_width=0.05,
2206        head_length=0.225,
2207        head_width=0.175,
2208        fill=True,
2209        c=None,
2210        alpha=1,
2211    ):
2212        """
2213        Build 2D arrows between two lists of points `start_pts` and `end_pts`.
2214        `start_pts` can be also passed in the form `[[point1, point2], ...]`.
2215
2216        Color can be specified as a colormap which maps the size of the arrows.
2217
2218        Arguments:
2219            shaft_length : (float)
2220                fractional shaft length
2221            shaft_width : (float)
2222                fractional shaft width
2223            head_length : (float)
2224                fractional head length
2225            head_width : (float)
2226                fractional head width
2227            fill : (bool)
2228                if False only generate the outline
2229        """
2230        if isinstance(start_pts, Points):
2231            start_pts = start_pts.points()
2232        if isinstance(end_pts, Points):
2233            end_pts = end_pts.points()
2234
2235        start_pts = np.asarray(start_pts, dtype=float)
2236        if end_pts is None:
2237            strt = start_pts[:, 0]
2238            end_pts = start_pts[:, 1]
2239            start_pts = strt
2240        else:
2241            end_pts = np.asarray(end_pts, dtype=float)
2242
2243        if head_length is None:
2244            head_length = 1 - shaft_length
2245
2246        arr = Arrow2D(
2247            (0, 0, 0),
2248            (1, 0, 0),
2249            s=s,
2250            shaft_length=shaft_length,
2251            shaft_width=shaft_width,
2252            head_length=head_length,
2253            head_width=head_width,
2254            fill=fill,
2255        )
2256
2257        orients = end_pts - start_pts
2258        orients = utils.make3d(orients)
2259
2260        pts = Points(start_pts)
2261        Glyph.__init__(
2262            self,
2263            pts,
2264            arr.polydata(False),
2265            orientation_array=orients,
2266            scale_by_vector_size=True,
2267            c=c,
2268            alpha=alpha,
2269        )
2270        self.flat().lighting("off")
2271        if c is not None:
2272            self.color(c)
2273        self.name = "Arrows2D"
2274
2275
2276class FlatArrow(Ribbon):
2277    """
2278    Build a 2D arrow in 3D space by joining two close lines.
2279    """
2280    def __init__(self, line1, line2, tip_size=1, tip_width=1):
2281        """
2282        Build a 2D arrow in 3D space by joining two close lines.
2283
2284        Examples:
2285            - [flatarrow.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/flatarrow.py)
2286
2287                ![](https://vedo.embl.es/images/basic/flatarrow.png)
2288        """
2289        if isinstance(line1, Points):
2290            line1 = line1.points()
2291        if isinstance(line2, Points):
2292            line2 = line2.points()
2293
2294        sm1, sm2 = np.array(line1[-1], dtype=float), np.array(line2[-1], dtype=float)
2295
2296        v = (sm1-sm2)/3*tip_width
2297        p1 = sm1+v
2298        p2 = sm2-v
2299        pm1 = (sm1+sm2)/2
2300        pm2 = (np.array(line1[-2])+np.array(line2[-2]))/2
2301        pm12 = pm1-pm2
2302        tip = pm12/np.linalg.norm(pm12)*np.linalg.norm(v)*3*tip_size/tip_width + pm1
2303
2304        line1.append(p1)
2305        line1.append(tip)
2306        line2.append(p2)
2307        line2.append(tip)
2308        resm = max(100, len(line1))
2309
2310        Ribbon.__init__(self, line1, line2, res=(resm, 1))
2311        self.phong()
2312        self.PickableOff()
2313        self.DragableOff()
2314        self.name = "FlatArrow"
2315
2316
2317class Triangle(Mesh):
2318    """Create a triangle from 3 points in space."""
2319    def __init__(self, p1, p2, p3, c="green7", alpha=1):
2320        """Create a triangle from 3 points in space."""
2321        Mesh.__init__(self, [[p1,p2,p3], [[0,1,2]]], c, alpha)
2322        self.GetProperty().LightingOff()
2323        self.name = "Triangle"
2324
2325
2326class Polygon(Mesh):
2327    """
2328    Build a polygon in the `xy` plane.
2329    """
2330    def __init__(self, pos=(0, 0, 0), nsides=6, r=1, c="coral", alpha=1):
2331        """
2332        Build a polygon in the `xy` plane of `nsides` of radius `r`.
2333
2334        ![](https://raw.githubusercontent.com/lorensen/VTKExamples/master/src/Testing/Baseline/Cxx/GeometricObjects/TestRegularPolygonSource.png)
2335        """
2336        t = np.linspace(np.pi / 2, 5 / 2 * np.pi, num=nsides, endpoint=False)
2337        x, y = utils.pol2cart(np.ones_like(t) * r, t)
2338        faces = [list(range(nsides))]
2339        # do not use: vtkRegularPolygonSource
2340        Mesh.__init__(self, [np.c_[x, y], faces], c, alpha)
2341        if len(pos) == 2:
2342            pos = (pos[0], pos[1], 0)
2343        self.SetPosition(pos)
2344        self.GetProperty().LightingOff()
2345        self.name = "Polygon " + str(nsides)
2346
2347
2348class Circle(Polygon):
2349    """
2350    Build a Circle of radius `r`.
2351    """
2352    def __init__(self, pos=(0, 0, 0), r=1, res=120, c="gray5", alpha=1):
2353        """
2354        Build a Circle of radius `r`.
2355        """
2356        Polygon.__init__(self, pos, nsides=res, r=r)
2357
2358        self.center = []  # filled by pointcloud.pcaEllipse
2359        self.nr_of_points = 0
2360        self.va = 0
2361        self.vb = 0
2362        self.axis1 = []
2363        self.axis2 = []
2364        self.alpha(alpha).c(c)
2365        self.name = "Circle"
2366
2367
2368class GeoCircle(Polygon):
2369    """
2370    Build a Circle of radius `r`.
2371    """
2372    def __init__(self, lat, lon, r=1, res=60, c="red4", alpha=1):
2373        """
2374        Build a Circle of radius `r` as projected on a geographic map.
2375        Circles near the poles will look very squashed.
2376
2377        See example:
2378            ```bash
2379            vedo -r earthquake
2380            ```
2381        """
2382        coords = []
2383        sinr, cosr = np.sin(r), np.cos(r)
2384        sinlat, coslat = np.sin(lat), np.cos(lat)
2385        for phi in np.linspace(0, 2 * np.pi, num=res, endpoint=False):
2386            clat = np.arcsin(sinlat * cosr + coslat * sinr * np.cos(phi))
2387            clng = lon + np.arctan2(np.sin(phi) * sinr * coslat, cosr - sinlat * np.sin(clat))
2388            coords.append([clng/np.pi + 1, clat*2/np.pi + 1, 0])
2389
2390        Polygon.__init__(self, nsides=res, c=c, alpha=alpha)
2391        self.points(coords)  # warp polygon points to match geo projection
2392        self.name = "Circle"
2393
2394
2395class Star(Mesh):
2396    """
2397    Build a 2D star shape.
2398    """
2399    def __init__(self, pos=(0, 0, 0), n=5, r1=0.7, r2=1.0, line=False, c="blue6", alpha=1):
2400        """
2401        Build a 2D star shape of `n` cusps of inner radius `r1` and outer radius `r2`.
2402
2403        If line is True then only build the outer line (no internal surface meshing).
2404
2405        Example:
2406            - [extrude.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/extrude.py)
2407
2408                ![](https://vedo.embl.es/images/basic/extrude.png)
2409        """
2410        t = np.linspace(np.pi / 2, 5 / 2 * np.pi, num=n, endpoint=False)
2411        x, y = utils.pol2cart(np.ones_like(t) * r2, t)
2412        pts = np.c_[x, y, np.zeros_like(x)]
2413
2414        apts = []
2415        for i, p in enumerate(pts):
2416            apts.append(p)
2417            if i + 1 < n:
2418                apts.append((p + pts[i + 1]) / 2 * r1 / r2)
2419        apts.append((pts[-1] + pts[0]) / 2 * r1 / r2)
2420
2421        if line:
2422            apts.append(pts[0])
2423            poly = utils.buildPolyData(apts, lines=list(range(len(apts))))
2424            Mesh.__init__(self, poly, c, alpha)
2425            self.lw(2)
2426        else:
2427            apts.append((0, 0, 0))
2428            cells = []
2429            for i in range(2 * n - 1):
2430                cell = [2 * n, i, i + 1]
2431                cells.append(cell)
2432            cells.append([2 * n, i + 1, 0])
2433            Mesh.__init__(self, [apts, cells], c, alpha)
2434
2435        if len(pos) == 2:
2436            pos = (pos[0], pos[1], 0)
2437        self.SetPosition(pos)
2438        self.property.LightingOff()
2439        self.name = "Star"
2440
2441
2442class Disc(Mesh):
2443    """
2444    Build a 2D disc.
2445    """
2446    def __init__(self, pos=(0, 0, 0), r1=0.5, r2=1, res=(1, 120), angle_range=(), c="gray4", alpha=1):
2447        """
2448        Build a 2D disc of inner radius `r1` and outer radius `r2`.
2449
2450        Set `res` as the resolution in R and Phi (can be a list).
2451
2452        Use `angle_range` to create a disc sector between the 2 specified angles.
2453
2454        ![](https://raw.githubusercontent.com/lorensen/VTKExamples/master/src/Testing/Baseline/Cxx/GeometricObjects/TestDisk.png)
2455        """
2456        if utils.is_sequence(res):
2457            res_r, res_phi = res
2458        else:
2459            res_r, res_phi = res, 12 * res
2460
2461        if len(angle_range) == 0:
2462            ps = vtk.vtkDiskSource()
2463        else:
2464            ps = vtk.vtkSectorSource()
2465            ps.SetStartAngle(angle_range[0])
2466            ps.SetEndAngle(angle_range[1])
2467
2468        ps.SetInnerRadius(r1)
2469        ps.SetOuterRadius(r2)
2470        ps.SetRadialResolution(res_r)
2471        ps.SetCircumferentialResolution(res_phi)
2472        ps.Update()
2473        Mesh.__init__(self, ps.GetOutput(), c, alpha)
2474        self.flat()
2475        self.SetPosition(utils.make3d(pos))
2476        self.name = "Disc"
2477
2478
2479class Arc(Mesh):
2480    """
2481    Build a 2D circular arc between 2 points.
2482    """
2483    def __init__(
2484        self,
2485        center,
2486        point1,
2487        point2=None,
2488        normal=None,
2489        angle=None,
2490        invert=False,
2491        res=50,
2492        c="gray4",
2493        alpha=1,
2494    ):
2495        """
2496        Build a 2D circular arc between 2 points `point1` and `point2`.
2497
2498        If `normal` is specified then `center` is ignored, and
2499        normal vector, a starting `point1` (polar vector)
2500        and an angle defining the arc length need to be assigned.
2501
2502        Arc spans the shortest angular sector point1 and point2,
2503        if `invert=True`, then the opposite happens.
2504        """
2505        if len(point1) == 2:
2506            point1 = (point1[0], point1[1], 0)
2507        if point2 is not None and len(point2) == 2:
2508            point2 = (point2[0], point2[1], 0)
2509
2510        self.base = point1
2511        self.top = point2
2512
2513        ar = vtk.vtkArcSource()
2514        if point2 is not None:
2515            self.top = point2
2516            point2 = point2 - np.asarray(point1)
2517            ar.UseNormalAndAngleOff()
2518            ar.SetPoint1([0, 0, 0])
2519            ar.SetPoint2(point2)
2520            ar.SetCenter(center)
2521        elif normal is not None and angle is not None:
2522            ar.UseNormalAndAngleOn()
2523            ar.SetAngle(angle)
2524            ar.SetPolarVector(point1)
2525            ar.SetNormal(normal)
2526        else:
2527            vedo.logger.error("incorrect input combination")
2528            return
2529        ar.SetNegative(invert)
2530        ar.SetResolution(res)
2531        ar.Update()
2532        Mesh.__init__(self, ar.GetOutput(), c, alpha)
2533        self.SetPosition(self.base)
2534        self.lw(2).lighting("off")
2535        self.name = "Arc"
2536
2537class IcoSphere(Mesh):
2538    """
2539    Create a sphere made of a uniform triangle mesh.
2540    """
2541    def __init__(self, pos=(0, 0, 0), r=1, subdivisions=3, c="r5", alpha=1):
2542        """
2543        Create a sphere made of a uniform triangle mesh
2544        (from recursive subdivision of an icosahedron).
2545
2546        Example:
2547        ```python
2548        from vedo import *
2549        icos = IcoSphere(subdivisions=3)
2550        icos.compute_quality().cmap('coolwarm')
2551        icos.show(axes=1).close()
2552        ```
2553        ![](https://vedo.embl.es/images/basic/icosphere.jpg)
2554        """
2555        subdivisions = int(min(subdivisions, 9)) # to avoid disasters
2556
2557        t = (1.0 + np.sqrt(5.0)) / 2.0
2558        points = np.array([
2559            [-1,  t,  0],
2560            [ 1,  t,  0],
2561            [-1, -t,  0],
2562            [ 1, -t,  0],
2563            [ 0, -1,  t],
2564            [ 0,  1,  t],
2565            [ 0, -1, -t],
2566            [ 0,  1, -t],
2567            [ t,  0, -1],
2568            [ t,  0,  1],
2569            [-t,  0, -1],
2570            [-t,  0,  1]
2571        ])
2572        faces = [
2573            [0, 11, 5],
2574            [0, 5, 1],
2575            [0, 1, 7],
2576            [0, 7, 10],
2577            [0, 10, 11],
2578            [1, 5, 9],
2579            [5, 11, 4],
2580            [11, 10, 2],
2581            [10, 7, 6],
2582            [7, 1, 8],
2583            [3, 9, 4],
2584            [3, 4, 2],
2585            [3, 2, 6],
2586            [3, 6, 8],
2587            [3, 8, 9],
2588            [4, 9, 5],
2589            [2, 4, 11],
2590            [6, 2, 10],
2591            [8, 6, 7],
2592            [9, 8, 1]
2593        ]
2594        Mesh.__init__(self, [points * r, faces], c=c, alpha=alpha)
2595
2596        for _ in range(subdivisions):
2597            self.subdivide(method=1)
2598            pts = utils.versor(self.points()) * r
2599            self.points(pts)
2600
2601        self.SetPosition(pos)
2602        self.name = "IcoSphere"
2603
2604
2605class Sphere(Mesh):
2606    """
2607    Build a sphere.
2608    """
2609    def __init__(self, pos=(0, 0, 0), r=1, res=24, quads=False, c="r5", alpha=1):
2610        """
2611        Build a sphere at position `pos` of radius `r`.
2612
2613        Arguments:
2614            r : (float)
2615                sphere radius
2616            res : (int, list)
2617                resolution in phi, resolution in theta is by default `2*res`
2618            quads : (bool)
2619                sphere mesh will be made of quads instead of triangles
2620
2621        [](https://user-images.githubusercontent.com/32848391/72433092-f0a31e00-3798-11ea-85f7-b2f5fcc31568.png)
2622        """
2623        if len(pos) == 2:
2624            pos = np.asarray([pos[0], pos[1], 0])
2625
2626        self.radius = r  # used by fitSphere
2627        self.center = pos
2628        self.residue = 0
2629
2630        if quads:
2631            res = max(res, 4)
2632            img = vtk.vtkImageData()
2633            img.SetDimensions(res - 1, res - 1, res - 1)
2634            rs = 1.0 / (res - 2)
2635            img.SetSpacing(rs, rs, rs)
2636            gf = vtk.vtkGeometryFilter()
2637            gf.SetInputData(img)
2638            gf.Update()
2639            Mesh.__init__(self, gf.GetOutput(), c, alpha)
2640            self.lw(0.1)
2641
2642            cgpts = self.points() - (0.5, 0.5, 0.5)
2643
2644            x, y, z = cgpts.T
2645            x = x * (1 + x * x) / 2
2646            y = y * (1 + y * y) / 2
2647            z = z * (1 + z * z) / 2
2648            _, theta, phi = utils.cart2spher(x, y, z)
2649
2650            pts = utils.spher2cart(np.ones_like(phi) * r, theta, phi)
2651            self.points(pts)
2652
2653        else:
2654            if utils.is_sequence(res):
2655                res_t, res_phi = res
2656            else:
2657                res_t, res_phi = 2 * res, res
2658
2659            ss = vtk.vtkSphereSource()
2660            ss.SetRadius(r)
2661            ss.SetThetaResolution(res_t)
2662            ss.SetPhiResolution(res_phi)
2663            ss.Update()
2664
2665            Mesh.__init__(self, ss.GetOutput(), c, alpha)
2666
2667        self.phong()
2668        self.SetPosition(pos)
2669        self.name = "Sphere"
2670
2671
2672class Spheres(Mesh):
2673    """
2674    Build a large set of spheres.
2675    """
2676    def __init__(self, centers, r=1, res=8, c="r5", alpha=1):
2677        """
2678        Build a (possibly large) set of spheres at `centers` of radius `r`.
2679
2680        Either `c` or `r` can be a list of RGB colors or radii.
2681
2682        Examples:
2683            - [manyspheres.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/manyspheres.py)
2684
2685            ![](https://vedo.embl.es/images/basic/manyspheres.jpg)
2686        """
2687
2688        if isinstance(centers, Points):
2689            centers = centers.points()
2690        centers = np.asarray(centers, dtype=float)
2691        base = centers[0]
2692
2693        cisseq = False
2694        if utils.is_sequence(c):
2695            cisseq = True
2696
2697        if cisseq:
2698            if len(centers) != len(c):
2699                vedo.logger.error(f"mismatch #centers {len(centers)} != {len(c)} #colors")
2700                raise RuntimeError()
2701
2702        risseq = False
2703        if utils.is_sequence(r):
2704            risseq = True
2705
2706        if risseq:
2707            if len(centers) != len(r):
2708                vedo.logger.error(f"mismatch #centers {len(centers)} != {len(r)} #radii")
2709                raise RuntimeError()
2710        if cisseq and risseq:
2711            vedo.logger.error("Limitation: c and r cannot be both sequences.")
2712            raise RuntimeError()
2713
2714        src = vtk.vtkSphereSource()
2715        if not risseq:
2716            src.SetRadius(r)
2717        if utils.is_sequence(res):
2718            res_t, res_phi = res
2719        else:
2720            res_t, res_phi = 2 * res, res
2721
2722        src.SetThetaResolution(res_t)
2723        src.SetPhiResolution(res_phi)
2724        src.Update()
2725
2726        psrc = vtk.vtkPointSource()
2727        psrc.SetNumberOfPoints(len(centers))
2728        psrc.Update()
2729        pd = psrc.GetOutput()
2730        vpts = pd.GetPoints()
2731
2732        glyph = vtk.vtkGlyph3D()
2733        glyph.SetSourceConnection(src.GetOutputPort())
2734
2735        if cisseq:
2736            glyph.SetColorModeToColorByScalar()
2737            ucols = vtk.vtkUnsignedCharArray()
2738            ucols.SetNumberOfComponents(3)
2739            ucols.SetName("Colors")
2740            for acol in c:
2741                cx, cy, cz = get_color(acol)
2742                ucols.InsertNextTuple3(cx * 255, cy * 255, cz * 255)
2743            pd.GetPointData().AddArray(ucols)
2744            pd.GetPointData().SetActiveScalars("Colors")
2745            glyph.ScalingOff()
2746        elif risseq:
2747            glyph.SetScaleModeToScaleByScalar()
2748            urads = utils.numpy2vtk(2 * np.ascontiguousarray(r), dtype=np.float32)
2749            urads.SetName("Radii")
2750            pd.GetPointData().AddArray(urads)
2751            pd.GetPointData().SetActiveScalars("Radii")
2752
2753        vpts.SetData(utils.numpy2vtk(centers - base, dtype=np.float32))
2754
2755        glyph.SetInputData(pd)
2756        glyph.Update()
2757
2758        Mesh.__init__(self, glyph.GetOutput(), alpha=alpha)
2759        self.SetPosition(base)
2760        self.base = base
2761        self.top = centers[-1]
2762        self.phong()
2763        if cisseq:
2764            self.mapper().ScalarVisibilityOn()
2765        else:
2766            self.mapper().ScalarVisibilityOff()
2767            self.GetProperty().SetColor(get_color(c))
2768        self.name = "Spheres"
2769
2770
2771class Earth(Mesh):
2772    """
2773    Build a textured mesh representing the Earth.
2774    """
2775    def __init__(self, style=1, r=1):
2776        """
2777        Build a textured mesh representing the Earth.
2778
2779        Example:
2780            - [geodesic.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/geodesic.py)
2781
2782                ![](https://vedo.embl.es/images/advanced/geodesic.png)
2783        """
2784        tss = vtk.vtkTexturedSphereSource()
2785        tss.SetRadius(r)
2786        tss.SetThetaResolution(72)
2787        tss.SetPhiResolution(36)
2788        Mesh.__init__(self, tss, c="w")
2789        atext = vtk.vtkTexture()
2790        pnm_reader = vtk.vtkJPEGReader()
2791        fn = vedo.io.download(vedo.dataurl + f"textures/earth{style}.jpg", verbose=False)
2792        pnm_reader.SetFileName(fn)
2793        atext.SetInputConnection(pnm_reader.GetOutputPort())
2794        atext.InterpolateOn()
2795        self.SetTexture(atext)
2796        self.name = "Earth"
2797
2798
2799class Ellipsoid(Mesh):
2800    """
2801    Build a 3D ellipsoid.
2802    """
2803    def __init__(
2804        self,
2805        pos=(0, 0, 0),
2806        axis1=(1, 0, 0),
2807        axis2=(0, 2, 0),
2808        axis3=(0, 0, 3),
2809        res=24,
2810        c="cyan4",
2811        alpha=1,
2812    ):
2813        """
2814        Build a 3D ellipsoid centered at position `pos`.
2815
2816        Arguments:
2817            axis1 : (list)
2818                First axis
2819            axis2 : (list)
2820                Second axis
2821            axis3 : (list)
2822                Third axis
2823
2824        .. note:: `axis1` and `axis2` are only used to define sizes and one azimuth angle.
2825        """
2826
2827        self.center = pos
2828        self.va_error = 0
2829        self.vb_error = 0
2830        self.vc_error = 0
2831        self.axis1 = axis1
2832        self.axis2 = axis2
2833        self.axis3 = axis3
2834        self.nr_of_points = 1  # used by pcaEllipsoid
2835
2836        if utils.is_sequence(res):
2837            res_t, res_phi = res
2838        else:
2839            res_t, res_phi = 2 * res, res
2840
2841        elli_source = vtk.vtkSphereSource()
2842        elli_source.SetThetaResolution(res_t)
2843        elli_source.SetPhiResolution(res_phi)
2844        elli_source.Update()
2845        l1 = np.linalg.norm(axis1)
2846        l2 = np.linalg.norm(axis2)
2847        l3 = np.linalg.norm(axis3)
2848        self.va = l1
2849        self.vb = l2
2850        self.vc = l3
2851        axis1 = np.array(axis1) / l1
2852        axis2 = np.array(axis2) / l2
2853        axis3 = np.array(axis3) / l3
2854        angle = np.arcsin(np.dot(axis1, axis2))
2855        theta = np.arccos(axis3[2])
2856        phi = np.arctan2(axis3[1], axis3[0])
2857
2858        t = vtk.vtkTransform()
2859        t.PostMultiply()
2860        t.Scale(l1, l2, l3)
2861        t.RotateX(np.rad2deg(angle))
2862        t.RotateY(np.rad2deg(theta))
2863        t.RotateZ(np.rad2deg(phi))
2864        tf = vtk.vtkTransformPolyDataFilter()
2865        tf.SetInputData(elli_source.GetOutput())
2866        tf.SetTransform(t)
2867        tf.Update()
2868        pd = tf.GetOutput()
2869        self.transformation = t
2870
2871        Mesh.__init__(self, pd, c, alpha)
2872        self.phong()
2873        if len(pos) == 2:
2874            pos = (pos[0], pos[1], 0)
2875        self.SetPosition(pos)
2876        self.name = "Ellipsoid"
2877
2878    def asphericity(self):
2879        """
2880        Return a measure of how different an ellipsoid is froma sphere.
2881        Values close to zero correspond to a spheric object.
2882        """
2883        a,b,c = self.va, self.vb, self.vc
2884        asp = ( ((a-b)/(a+b))**2
2885              + ((a-c)/(a+c))**2
2886              + ((b-c)/(b+c))**2 )/3. * 4.
2887        return asp
2888
2889    def asphericity_error(self):
2890        """
2891        Calculate statistical error on the asphericity value.
2892
2893        Errors on the main axes are stored in
2894        `Ellipsoid.va_error, Ellipsoid.vb_error and Ellipsoid.vc_error`.
2895        """
2896        a, b, c = self.va, self.vb, self.vc
2897        sqrtn = np.sqrt(self.nr_of_points)
2898        ea, eb, ec = a / 2 / sqrtn, b / 2 / sqrtn, b / 2 / sqrtn
2899
2900        # from sympy import *
2901        # init_printing(use_unicode=True)
2902        # a, b, c, ea, eb, ec = symbols("a b c, ea, eb,ec")
2903        # L = (
2904        #    (((a - b) / (a + b)) ** 2 + ((c - b) / (c + b)) ** 2 + ((a - c) / (a + c)) ** 2)
2905        #    / 3 * 4)
2906        # dl2 = (diff(L, a) * ea) ** 2 + (diff(L, b) * eb) ** 2 + (diff(L, c) * ec) ** 2
2907        # print(dl2)
2908        # exit()
2909        dL2 = (
2910            ea ** 2
2911            * (
2912                -8 * (a - b) ** 2 / (3 * (a + b) ** 3)
2913                - 8 * (a - c) ** 2 / (3 * (a + c) ** 3)
2914                + 4 * (2 * a - 2 * c) / (3 * (a + c) ** 2)
2915                + 4 * (2 * a - 2 * b) / (3 * (a + b) ** 2)
2916            ) ** 2
2917            + eb ** 2
2918            * (
2919                4 * (-2 * a + 2 * b) / (3 * (a + b) ** 2)
2920                - 8 * (a - b) ** 2 / (3 * (a + b) ** 3)
2921                - 8 * (-b + c) ** 2 / (3 * (b + c) ** 3)
2922                + 4 * (2 * b - 2 * c) / (3 * (b + c) ** 2)
2923            ) ** 2
2924            + ec ** 2
2925            * (
2926                4 * (-2 * a + 2 * c) / (3 * (a + c) ** 2)
2927                - 8 * (a - c) ** 2 / (3 * (a + c) ** 3)
2928                + 4 * (-2 * b + 2 * c) / (3 * (b + c) ** 2)
2929                - 8 * (-b + c) ** 2 / (3 * (b + c) ** 3)
2930            ) ** 2
2931        )
2932
2933        err = np.sqrt(dL2)
2934
2935        self.va_error = ea
2936        self.vb_error = eb
2937        self.vc_error = ec
2938        return err
2939
2940
2941class Grid(Mesh):
2942    """
2943    An even or uneven 2D grid.
2944    """
2945    def __init__(
2946        self,
2947        pos=(0, 0, 0),
2948        s=(1,1),
2949        res=(10,10),
2950        lw=1,
2951        c="k3",
2952        alpha=1,
2953    ):
2954        """
2955        Create an even or uneven 2D grid.
2956
2957        Arguments:
2958            s : (float, list)
2959                if a float is provided it is interpreted as the total size along x and y,
2960                if a list of coords is provided they are interpreted as the vertices of the grid along x and y.
2961                In this case keyword `res` is ignored (see example below).
2962            res : (list)
2963                resolutions along x and y, e.i. the number of subdivisions
2964            lw : (int)
2965                line width
2966
2967        Example:
2968            ```python
2969            from vedo import *
2970            import numpy as np
2971            xcoords = np.arange(0, 2, 0.2)
2972            ycoords = np.arange(0, 1, 0.2)
2973            sqrtx = sqrt(xcoords)
2974            grid = Grid(s=(sqrtx, ycoords)).lw(2)
2975            grid.show(axes=8)
2976
2977            # can also create a grid from np.mgrid:
2978            X, Y = np.mgrid[-12:12:1000*1j, 0:15:1000*1j]
2979            vgrid = Grid(s=(X[:,0], Y[0]))
2980            vgrid.show(axes=1).close()
2981            ```
2982            ![](https://vedo.embl.es/images/feats/uneven_grid.png)
2983        """
2984        resx, resy = res
2985        sx, sy = s
2986
2987        if len(pos) == 2:
2988            pos = (pos[0], pos[1], 0)
2989
2990        if utils.is_sequence(sx) and utils.is_sequence(sy):
2991            verts = []
2992            for y in sy:
2993                for x in sx:
2994                    verts.append([x, y, 0])
2995            faces = []
2996            n = len(sx)
2997            m = len(sy)
2998            for j in range(m - 1):
2999                j1n = (j + 1) * n
3000                for i in range(n - 1):
3001                    faces.append([i + j * n, i + 1 + j * n, i + 1 + j1n, i + j1n])
3002
3003            verts = np.array(verts)
3004            Mesh.__init__(self, [verts, faces], c, alpha)
3005
3006        else:
3007            ps = vtk.vtkPlaneSource()
3008            ps.SetResolution(resx, resy)
3009            ps.Update()
3010            poly0 = ps.GetOutput()
3011            t0 = vtk.vtkTransform()
3012            t0.Scale(sx, sy, 1)
3013            tf0 = vtk.vtkTransformPolyDataFilter()
3014            tf0.SetInputData(poly0)
3015            tf0.SetTransform(t0)
3016            tf0.Update()
3017            poly = tf0.GetOutput()
3018            Mesh.__init__(self, poly, c, alpha)
3019            self.SetPosition(pos)
3020
3021        self.wireframe().lw(lw)
3022        self.GetProperty().LightingOff()
3023        self.name = "Grid"
3024
3025
3026class Plane(Mesh):
3027    """
3028    Create a plane in space.
3029    """
3030    def __init__(self, pos=(0, 0, 0), normal=(0, 0, 1), s=(1, 1), res=(1, 1), c="gray5", alpha=1):
3031        """
3032        Create a plane of size `s=(xsize, ysize)` oriented perpendicular to vector `normal`
3033        and so that it passes through point `pos`.
3034
3035        Arguments:
3036            normal : (list)
3037                normal vector to the plane
3038        """
3039        pos = utils.make3d(pos)
3040        sx, sy = s
3041
3042        self.normal = np.asarray(normal, dtype=float)
3043        self.center = np.asarray(pos, dtype=float)
3044        self.variance = 0
3045
3046        ps = vtk.vtkPlaneSource()
3047        ps.SetResolution(res[0], res[1])
3048        tri = vtk.vtkTriangleFilter()
3049        tri.SetInputConnection(ps.GetOutputPort())
3050        tri.Update()
3051        poly = tri.GetOutput()
3052        axis = self.normal / np.linalg.norm(normal)
3053        theta = np.arccos(axis[2])
3054        phi = np.arctan2(axis[1], axis[0])
3055        t = vtk.vtkTransform()
3056        t.PostMultiply()
3057        t.Scale(sx, sy, 1)
3058        t.RotateY(np.rad2deg(theta))
3059        t.RotateZ(np.rad2deg(phi))
3060        tf = vtk.vtkTransformPolyDataFilter()
3061        tf.SetInputData(poly)
3062        tf.SetTransform(t)
3063        tf.Update()
3064        Mesh.__init__(self, tf.GetOutput(), c, alpha)
3065        self.lighting("off")
3066        self.SetPosition(pos)
3067        self.name = "Plane"
3068        self.top = self.normal
3069        self.bottom = np.array([0.0, 0.0, 0.0])
3070
3071    def contains(self, points):
3072        """
3073        Check if each of the provided point lies on this plane.
3074        `points` is an array of shape (n, 3).
3075        """
3076        points = np.array(points, dtype=float)
3077        bounds = self.points()
3078
3079        mask = np.isclose(np.dot(points - self.center, self.normal), 0)
3080
3081        for i in [1, 3]:
3082            AB = bounds[i] - bounds[0]
3083            AP = points - bounds[0]
3084            mask_l = np.less_equal(np.dot(AP, AB), np.linalg.norm(AB))
3085            mask_g = np.greater_equal(np.dot(AP, AB), 0)
3086            mask = np.logical_and(mask, mask_l)
3087            mask = np.logical_and(mask, mask_g)
3088        return mask
3089
3090
3091class Rectangle(Mesh):
3092    """
3093    Build a rectangle in the xy plane.
3094    """
3095    def __init__(self, p1=(0, 0), p2=(1, 1), radius=None, res=12, c="gray5", alpha=1):
3096        """
3097        Build a rectangle in the xy plane identified by any two corner points.
3098
3099        Arguments:
3100            p1 : (list)
3101                bottom-left position of the corner
3102            p2 : (list)
3103                top-right position of the corner
3104            radius : (float, list)
3105                smoothing radius of the corner in world units.
3106                A list can be passed with 4 individual values.
3107        """
3108        if len(p1) == 2:
3109            p1 = np.array([p1[0], p1[1], 0.0])
3110        else:
3111            p1 = np.array(p1, dtype=float)
3112        if len(p2) == 2:
3113            p2 = np.array([p2[0], p2[1], 0.0])
3114        else:
3115            p2 = np.array(p2, dtype=float)
3116
3117        self.corner1 = p1
3118        self.corner2 = p2
3119
3120        color = c
3121        smoothr = False
3122        risseq = False
3123        if utils.is_sequence(radius):
3124            risseq = True
3125            smoothr = True
3126            if max(radius) == 0:
3127                smoothr = False
3128        elif radius:
3129            smoothr = True
3130
3131        if not smoothr:
3132            radius = None
3133        self.radius = radius
3134
3135        if smoothr:
3136            r = radius
3137            if not risseq:
3138                r = [r, r, r, r]
3139            rd, ra, rb, rc = r
3140
3141            if p1[0] > p2[0]:  # flip p1 - p2
3142                p1, p2 = p2, p1
3143            if p1[1] > p2[1]:  # flip p1y - p2y
3144                p1[1], p2[1] = p2[1], p1[1]
3145
3146            px, py, _ = p2 - p1
3147            k = min(px / 2, py / 2)
3148            ra = min(abs(ra), k)
3149            rb = min(abs(rb), k)
3150            rc = min(abs(rc), k)
3151            rd = min(abs(rd), k)
3152            beta = np.linspace(0, 2 * np.pi, num=res * 4, endpoint=False)
3153            betas = np.split(beta, 4)
3154            rrx = np.cos(betas)
3155            rry = np.sin(betas)
3156
3157            q1 = (rd, 0)
3158            # q2 = (px-ra, 0)
3159            q3 = (px, ra)
3160            # q4 = (px, py-rb)
3161            q5 = (px - rb, py)
3162            # q6 = (rc, py)
3163            q7 = (0, py - rc)
3164            # q8 = (0, rd)
3165            a = np.c_[rrx[3], rry[3]]*ra + [px-ra, ra]    if ra else np.array([])
3166            b = np.c_[rrx[0], rry[0]]*rb + [px-rb, py-rb] if rb else np.array([])
3167            c = np.c_[rrx[1], rry[1]]*rc + [rc, py-rc]    if rc else np.array([])
3168            d = np.c_[rrx[2], rry[2]]*rd + [rd, rd]       if rd else np.array([])
3169
3170            pts = [q1, *a.tolist(), q3, *b.tolist(), q5, *c.tolist(), q7, *d.tolist()]
3171            faces = [list(range(len(pts)))]
3172        else:
3173            p1r = np.array([p2[0], p1[1], 0.0])
3174            p2l = np.array([p1[0], p2[1], 0.0])
3175            pts = ([0.0, 0.0, 0.0], p1r - p1, p2 - p1, p2l - p1)
3176            faces = [(0, 1, 2, 3)]
3177
3178        Mesh.__init__(self, [pts, faces], color, alpha)
3179        self.SetPosition(p1)
3180        self.property.LightingOff()
3181        self.name = "Rectangle"
3182
3183
3184class Box(Mesh):
3185    """
3186    Build a box of specified dimensions.
3187    """
3188    def __init__(self, pos=(0, 0, 0), length=1, width=2, height=3, size=(), c="g4", alpha=1):
3189        """
3190        Build a box of dimensions `x=length, y=width and z=height`.
3191        Alternatively dimensions can be defined by setting `size` keyword with a tuple.
3192
3193        If `size` is a list of 6 numbers, this will be interpreted as the bounding box:
3194        `[xmin,xmax, ymin,ymax, zmin,zmax]`
3195
3196        Examples:
3197            - [aspring.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/aspring.py)
3198
3199                ![](https://vedo.embl.es/images/simulations/50738955-7e891800-11d9-11e9-85cd-02bd4f3f13ea.gif)
3200        """
3201        if len(size) == 6:
3202            bounds = size
3203            length = bounds[1] - bounds[0]
3204            width = bounds[3] - bounds[2]
3205            height = bounds[5] - bounds[4]
3206            xp = (bounds[1] + bounds[0]) / 2
3207            yp = (bounds[3] + bounds[2]) / 2
3208            zp = (bounds[5] + bounds[4]) / 2
3209            pos = (xp, yp, zp)
3210        elif len(size) == 3:
3211            length, width, height = size
3212
3213        src = vtk.vtkCubeSource()
3214        src.SetXLength(length)
3215        src.SetYLength(width)
3216        src.SetZLength(height)
3217        src.Update()
3218        pd = src.GetOutput()
3219
3220        tc = [
3221            [0.0, 0.0],
3222            [1.0, 0.0],
3223            [0.0, 1.0],
3224            [1.0, 1.0],
3225            [1.0, 0.0],
3226            [0.0, 0.0],
3227            [1.0, 1.0],
3228            [0.0, 1.0],
3229            [1.0, 1.0],
3230            [1.0, 0.0],
3231            [0.0, 1.0],
3232            [0.0, 0.0],
3233            [0.0, 1.0],
3234            [0.0, 0.0],
3235            [1.0, 1.0],
3236            [1.0, 0.0],
3237            [1.0, 0.0],
3238            [0.0, 0.0],
3239            [1.0, 1.0],
3240            [0.0, 1.0],
3241            [0.0, 0.0],
3242            [1.0, 0.0],
3243            [0.0, 1.0],
3244            [1.0, 1.0],
3245        ]
3246        vtc = utils.numpy2vtk(tc)
3247        pd.GetPointData().SetTCoords(vtc)
3248        Mesh.__init__(self, pd, c, alpha)
3249        if len(pos) == 2:
3250            pos = (pos[0], pos[1], 0)
3251        self.SetPosition(pos)
3252        self.name = "Box"
3253
3254
3255class Cube(Box):
3256    """Build a cube."""
3257    def __init__(self, pos=(0, 0, 0), side=1, c="g4", alpha=1):
3258       """Build a cube of size `side`."""
3259       Box.__init__(self, pos, side, side, side, (), c, alpha)
3260       self.name = "Cube"
3261
3262
3263class TessellatedBox(Mesh):
3264    """
3265    Build a cubic `Mesh` made of quads.
3266    """
3267    def __init__(self, pos=(0, 0, 0), n=10, spacing=(1, 1, 1), bounds=(),
3268                 c="k5", alpha=0.5,
3269    ):
3270        """
3271        Build a cubic `Mesh` made of `n` small quads in the 3 axis directions.
3272
3273        Arguments:
3274            pos : (list)
3275                position of the left bottom corner
3276            n : (int, list)
3277                number of subdivisions along each side
3278            spacing : (float)
3279                size of the side of the single quad in the 3 directions
3280        """
3281        if utils.is_sequence(n):  # slow
3282            img = vtk.vtkImageData()
3283            img.SetDimensions(n[0] + 1, n[1] + 1, n[2] + 1)
3284            img.SetSpacing(spacing)
3285            gf = vtk.vtkGeometryFilter()
3286            gf.SetInputData(img)
3287            gf.Update()
3288            poly = gf.GetOutput()
3289        else:  # fast
3290            n -= 1
3291            tbs = vtk.vtkTessellatedBoxSource()
3292            tbs.SetLevel(n)
3293            if len(bounds):
3294                tbs.SetBounds(bounds)
3295            else:
3296                tbs.SetBounds(0, n * spacing[0], 0, n * spacing[1], 0, n * spacing[2])
3297            tbs.QuadsOn()
3298            tbs.SetOutputPointsPrecision(vtk.vtkAlgorithm.SINGLE_PRECISION)
3299            tbs.Update()
3300            poly = tbs.GetOutput()
3301        Mesh.__init__(self, poly, c=c, alpha=alpha)
3302        self.SetPosition(pos)
3303        self.lw(1).lighting("off")
3304        self.base = np.array([0.5, 0.5, 0.0])
3305        self.top = np.array([0.5, 0.5, 1.0])
3306        self.name = "TessellatedBox"
3307
3308
3309class Spring(Mesh):
3310    """
3311    Build a spring model.
3312    """
3313    def __init__(
3314        self,
3315        start_pt=(0, 0, 0),
3316        end_pt=(1, 0, 0),
3317        coils=20,
3318        r1=0.1,
3319        r2=None,
3320        thickness=None,
3321        c="gray5",
3322        alpha=1,
3323    ):
3324        """
3325        Build a spring of specified nr of `coils` between `start_pt` and `end_pt`.
3326
3327        Arguments:
3328            coils : (int)
3329                number of coils
3330            r1 : (float)
3331                radius at start point
3332            r2 : (float)
3333                radius at end point
3334            thickness : (float)
3335                thickness of the coil section
3336        """
3337        diff = end_pt - np.array(start_pt, dtype=float)
3338        length = np.linalg.norm(diff)
3339        if not length:
3340            return
3341        if not r1:
3342            r1 = length / 20
3343        trange = np.linspace(0, length, num=50 * coils)
3344        om = 6.283 * (coils - 0.5) / length
3345        if not r2:
3346            r2 = r1
3347        pts = []
3348        for t in trange:
3349            f = (length - t) / length
3350            rd = r1 * f + r2 * (1 - f)
3351            pts.append([rd * np.cos(om * t), rd * np.sin(om * t), t])
3352
3353        pts = [[0, 0, 0]] + pts + [[0, 0, length]]
3354        diff = diff / length
3355        theta = np.arccos(diff[2])
3356        phi = np.arctan2(diff[1], diff[0])
3357        sp = Line(pts).polydata(False)
3358        t = vtk.vtkTransform()
3359        t.RotateZ(np.rad2deg(phi))
3360        t.RotateY(np.rad2deg(theta))
3361        tf = vtk.vtkTransformPolyDataFilter()
3362        tf.SetInputData(sp)
3363        tf.SetTransform(t)
3364        tf.Update()
3365        tuf = vtk.vtkTubeFilter()
3366        tuf.SetNumberOfSides(12)
3367        tuf.CappingOn()
3368        tuf.SetInputData(tf.GetOutput())
3369        if not thickness:
3370            thickness = r1 / 10
3371        tuf.SetRadius(thickness)
3372        tuf.Update()
3373        Mesh.__init__(self, tuf.GetOutput(), c, alpha)
3374        self.phong()
3375        self.SetPosition(start_pt)
3376        self.base = np.array(start_pt, dtype=float)
3377        self.top = np.array(end_pt, dtype=float)
3378        self.name = "Spring"
3379
3380
3381class Cylinder(Mesh):
3382    """
3383    Build a cylinder of specified height and radius.
3384    """
3385    def __init__(
3386        self,
3387        pos=(0, 0, 0),
3388        r=1,
3389        height=2,
3390        axis=(0, 0, 1),
3391        cap=True,
3392        res=24,
3393        c="teal3",
3394        alpha=1,
3395    ):
3396        """
3397        Build a cylinder of specified height and radius `r`, centered at `pos`.
3398
3399        If `pos` is a list of 2 points, e.g. `pos=[v1,v2]`, build a cylinder with base
3400        centered at `v1` and top at `v2`.
3401
3402        ![](https://raw.githubusercontent.com/lorensen/VTKExamples/master/src/Testing/Baseline/Cxx/GeometricObjects/TestCylinder.png)
3403        """
3404        if utils.is_sequence(pos[0]):  # assume user is passing pos=[base, top]
3405            base = np.array(pos[0], dtype=float)
3406            top = np.array(pos[1], dtype=float)
3407            pos = (base + top) / 2
3408            height = np.linalg.norm(top - base)
3409            axis = top - base
3410            axis = utils.versor(axis)
3411        else:
3412            axis = utils.versor(axis)
3413            base = pos - axis * height / 2
3414            top = pos + axis * height / 2
3415
3416        cyl = vtk.vtkCylinderSource()
3417        cyl.SetResolution(res)
3418        cyl.SetRadius(r)
3419        cyl.SetHeight(height)
3420        cyl.SetCapping(cap)
3421        cyl.Update()
3422
3423        theta = np.arccos(axis[2])
3424        phi = np.arctan2(axis[1], axis[0])
3425        t = vtk.vtkTransform()
3426        t.PostMultiply()
3427        t.RotateX(90)  # put it along Z
3428        t.RotateY(np.rad2deg(theta))
3429        t.RotateZ(np.rad2deg(phi))
3430        tf = vtk.vtkTransformPolyDataFilter()
3431        tf.SetInputData(cyl.GetOutput())
3432        tf.SetTransform(t)
3433        tf.Update()
3434        pd = tf.GetOutput()
3435
3436        Mesh.__init__(self, pd, c, alpha)
3437        self.phong()
3438        self.SetPosition(pos)
3439        self.base = base + pos
3440        self.top = top + pos
3441        self.name = "Cylinder"
3442
3443
3444class Cone(Mesh):
3445    """Build a cone of specified radius and height."""
3446    def __init__(
3447        self, pos=(0, 0, 0), r=1, height=3, axis=(0, 0, 1), res=48, c="green3", alpha=1
3448    ):
3449        """Build a cone of specified radius `r` and `height`, centered at `pos`."""
3450        con = vtk.vtkConeSource()
3451        con.SetResolution(res)
3452        con.SetRadius(r)
3453        con.SetHeight(height)
3454        con.SetDirection(axis)
3455        con.Update()
3456        Mesh.__init__(self, con.GetOutput(), c, alpha)
3457        self.phong()
3458        if len(pos) == 2:
3459            pos = (pos[0], pos[1], 0)
3460        self.SetPosition(pos)
3461        v = utils.versor(axis) * height / 2
3462        self.base = pos - v
3463        self.top = pos + v
3464        self.name = "Cone"
3465
3466
3467class Pyramid(Cone):
3468    """Build a pyramidal shape."""
3469    def __init__(self, pos=(0, 0, 0), s=1, height=1, axis=(0, 0, 1), c="green3", alpha=1):
3470        """Build a pyramid of specified base size `s` and `height`, centered at `pos`."""
3471        Cone.__init__(self, pos, s, height, axis, 4, c, alpha)
3472        self.name = "Pyramid"
3473
3474
3475class Torus(Mesh):
3476    """
3477    Build a toroidal shape.
3478    """
3479    def __init__(self, pos=(0, 0, 0), r1=1, r2=0.2, res=36, quads=False, c="yellow3", alpha=1):
3480        """
3481        Build a torus of specified outer radius `r1` internal radius `r2`, centered at `pos`.
3482        If `quad=True` a quad-mesh is generated.
3483        """
3484        if utils.is_sequence(res):
3485            res_u, res_v = res
3486        else:
3487            res_u, res_v = 3*res, res
3488
3489        if quads:
3490            #https://github.com/marcomusy/vedo/issues/710
3491
3492            n = res_v
3493            m = res_u
3494
3495            theta = np.linspace(0, 2.*np.pi, n)
3496            phi   = np.linspace(0, 2.*np.pi, m)
3497            theta, phi = np.meshgrid(theta, phi)
3498            t = r1 + r2 * np.cos(theta)
3499            x = t * np.cos(phi)
3500            y = t * np.sin(phi)
3501            z = r2 * np.sin(theta)
3502            pts = np.column_stack((x.ravel(), y.ravel(), z.ravel()))
3503
3504            faces = []
3505            for j in range(m-1):
3506                j1n = (j+1) * n
3507                for i in range(n-1):
3508                    faces.append([i + j*n, i+1 + j*n, i+1 + j1n, i + j1n])
3509
3510            Mesh.__init__(self, [pts, faces], c, alpha)
3511
3512        else:
3513
3514            rs = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricTorus()
3515            rs.SetRingRadius(r1)
3516            rs.SetCrossSectionRadius(r2)
3517            pfs = vtk.vtkParametricFunctionSource()
3518            pfs.SetParametricFunction(rs)
3519            pfs.SetUResolution(res_u)
3520            pfs.SetVResolution(res_v)
3521            pfs.Update()
3522            Mesh.__init__(self, pfs.GetOutput(), c, alpha)
3523
3524        self.phong()
3525        if len(pos) == 2:
3526            pos = (pos[0], pos[1], 0)
3527        self.SetPosition(pos)
3528        self.name = "Torus"
3529
3530
3531class Paraboloid(Mesh):
3532    """
3533    Build a paraboloid.
3534    """
3535    def __init__(self, pos=(0, 0, 0), height=1, res=50, c="cyan5", alpha=1):
3536        """
3537        Build a paraboloid of specified height and radius `r`, centered at `pos`.
3538
3539        Full volumetric expression is:
3540            `F(x,y,z)=a_0x^2+a_1y^2+a_2z^2+a_3xy+a_4yz+a_5xz+ a_6x+a_7y+a_8z+a_9`
3541
3542        ![](https://user-images.githubusercontent.com/32848391/51211547-260ef480-1916-11e9-95f6-4a677e37e355.png)
3543        """
3544        quadric = vtk.vtkQuadric()
3545        quadric.SetCoefficients(1, 1, 0, 0, 0, 0, 0, 0, height / 4, 0)
3546        # F(x,y,z) = a0*x^2 + a1*y^2 + a2*z^2
3547        #         + a3*x*y + a4*y*z + a5*x*z
3548        #         + a6*x   + a7*y   + a8*z  +a9
3549        sample = vtk.vtkSampleFunction()
3550        sample.SetSampleDimensions(res, res, res)
3551        sample.SetImplicitFunction(quadric)
3552
3553        contours = vtk.vtkContourFilter()
3554        contours.SetInputConnection(sample.GetOutputPort())
3555        contours.GenerateValues(1, 0.01, 0.01)
3556        contours.Update()
3557
3558        Mesh.__init__(self, contours.GetOutput(), c, alpha)
3559        self.compute_normals().phong()
3560        self.mapper().ScalarVisibilityOff()
3561        self.SetPosition(pos)
3562        self.name = "Paraboloid"
3563
3564
3565class Hyperboloid(Mesh):
3566    """
3567    Build a hyperboloid.
3568    """
3569    def __init__(self, pos=(0,0,0), a2=1, value=0.5, res=100, c="pink4", alpha=1):
3570        """
3571        Build a hyperboloid of specified aperture `a2` and `height`, centered at `pos`.
3572
3573        Full volumetric expression is:
3574            `F(x,y,z)=a_0x^2+a_1y^2+a_2z^2+a_3xy+a_4yz+a_5xz+ a_6x+a_7y+a_8z+a_9`
3575        """
3576        q = vtk.vtkQuadric()
3577        q.SetCoefficients(2, 2, -1 / a2, 0, 0, 0, 0, 0, 0, 0)
3578        # F(x,y,z) = a0*x^2 + a1*y^2 + a2*z^2
3579        #         + a3*x*y + a4*y*z + a5*x*z
3580        #         + a6*x   + a7*y   + a8*z  +a9
3581        sample = vtk.vtkSampleFunction()
3582        sample.SetSampleDimensions(res, res, res)
3583        sample.SetImplicitFunction(q)
3584
3585        contours = vtk.vtkContourFilter()
3586        contours.SetInputConnection(sample.GetOutputPort())
3587        contours.GenerateValues(1, value, value)
3588        contours.Update()
3589
3590        Mesh.__init__(self, contours.GetOutput(), c, alpha)
3591        self.compute_normals().phong()
3592        self.mapper().ScalarVisibilityOff()
3593        self.SetPosition(pos)
3594        self.name = "Hyperboloid"
3595
3596
3597def Marker(symbol, pos=(0, 0, 0), c="k", alpha=1, s=0.1, filled=True):
3598    """
3599    Generate a marker shape. Typically used in association with `Glyph`.
3600    """
3601    if isinstance(symbol, Mesh):
3602        return symbol.c(c).alpha(alpha).lighting("off")
3603
3604    if isinstance(symbol, int):
3605        symbs = ['.','o','O', '0', 'p','*','h','D','d','v','^','>','<','s', 'x', 'a']
3606        symbol = symbol % len(symbs)
3607        symbol = symbs[symbol]
3608
3609    if symbol == ".":
3610        mesh = Polygon(nsides=24, r=s * 0.6)
3611    elif symbol == "o":
3612        mesh = Polygon(nsides=24, r=s * 0.75)
3613    elif symbol == "O":
3614        mesh = Disc(r1=s * 0.6, r2=s * 0.75, res=(1, 24))
3615    elif symbol == "0":
3616        m1 = Disc(r1=s * 0.6, r2=s * 0.75, res=(1, 24))
3617        m2 = Circle(r=s * 0.36).reverse()
3618        mesh = merge(m1, m2)
3619    elif symbol == "p":
3620        mesh = Polygon(nsides=5, r=s)
3621    elif symbol == "*":
3622        mesh = Star(r1=0.65 * s * 1.1, r2=s * 1.1, line=not filled)
3623    elif symbol == "h":
3624        mesh = Polygon(nsides=6, r=s)
3625    elif symbol == "D":
3626        mesh = Polygon(nsides=4, r=s)
3627    elif symbol == "d":
3628        mesh = Polygon(nsides=4, r=s * 1.1).scale([0.5, 1, 1])
3629    elif symbol == "v":
3630        mesh = Polygon(nsides=3, r=s).rotate_z(180)
3631    elif symbol == "^":
3632        mesh = Polygon(nsides=3, r=s)
3633    elif symbol == ">":
3634        mesh = Polygon(nsides=3, r=s).rotate_z(-90)
3635    elif symbol == "<":
3636        mesh = Polygon(nsides=3, r=s).rotate_z(90)
3637    elif symbol == "s":
3638        mesh = Mesh(
3639            [[[-1, -1, 0], [1, -1, 0], [1, 1, 0], [-1, 1, 0]], [[0, 1, 2, 3]]]
3640        ).scale(s / 1.4)
3641    elif symbol == "x":
3642        mesh = Text3D("+", pos=(0, 0, 0), s=s * 2.6, justify="center", depth=0)
3643        # mesh.rotate_z(45)
3644    elif symbol == "a":
3645        mesh = Text3D("*", pos=(0, 0, 0), s=s * 2.6, justify="center", depth=0)
3646    else:
3647        mesh = Text3D(symbol, pos=(0, 0, 0), s=s * 2, justify="center", depth=0)
3648    mesh.flat().lighting("off").wireframe(not filled).c(c).alpha(alpha)
3649    if len(pos) == 2:
3650        pos = (pos[0], pos[1], 0)
3651    mesh.SetPosition(pos)
3652    mesh.name = "Marker"
3653    return mesh
3654
3655
3656class Brace(Mesh):
3657    """
3658    Create a brace (bracket) shape.
3659    """
3660    def __init__(
3661        self,
3662        q1,
3663        q2,
3664        style="}",
3665        padding1=0,
3666        font="Theemim",
3667        comment="",
3668        justify=None,
3669        angle=0,
3670        padding2=0.2,
3671        s=1,
3672        italic=0,
3673        c="k1",
3674        alpha=1,
3675    ):
3676        """
3677        Create a brace (bracket) shape which spans from point q1 to point q2.
3678
3679        Arguments:
3680            q1 : (list)
3681                point 1.
3682            q2 : (list)
3683                point 2.
3684            style : (str)
3685                style of the bracket, eg. `{}, [], (), <>`.
3686            padding1 : (float)
3687                padding space in percent form the input points.
3688            font : (str)
3689                font type
3690            comment : (str)
3691                additional text to appear next to the brace symbol.
3692            justify : (str)
3693                specify the anchor point to justify text comment, e.g. "top-left".
3694            italic : float
3695                italicness of the text comment (can be a positive or negative number)
3696            angle : (float)
3697                rotation angle of text. Use `None` to keep it horizontal.
3698            padding2 : (float)
3699                padding space in percent form brace to text comment.
3700            s : (float)
3701                scale factor for the comment
3702
3703        Examples:
3704            - [scatter3.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/scatter3.py)
3705
3706                ![](https://vedo.embl.es/images/pyplot/scatter3.png)
3707        """
3708        if isinstance(q1, vtk.vtkActor):
3709            q1 = q1.GetPosition()
3710        if isinstance(q2, vtk.vtkActor):
3711            q2 = q2.GetPosition()
3712        if len(q1) == 2:
3713            q1 = [q1[0], q1[1], 0.0]
3714        if len(q2) == 2:
3715            q2 = [q2[0], q2[1], 0.0]
3716        q1 = np.array(q1, dtype=float)
3717        q2 = np.array(q2, dtype=float)
3718        mq = (q1 + q2) / 2
3719        q1 = q1 - mq
3720        q2 = q2 - mq
3721        d = np.linalg.norm(q2 - q1)
3722        q2[2] = q1[2]
3723
3724        if style not in "{}[]()<>|I":
3725            vedo.logger.error(f"unknown style {style}." + "Use {}[]()<>|I")
3726            style = "}"
3727
3728        flip = False
3729        if style in ["{", "[", "(", "<"]:
3730            flip = True
3731            i = ["{", "[", "(", "<"].index(style)
3732            style = ["}", "]", ")", ">"][i]
3733
3734        br = Text3D(style, font="Theemim", justify="center-left")
3735        br.scale([0.4, 1, 1])
3736
3737        angler = np.arctan2(q2[1], q2[0]) * 180 / np.pi - 90
3738        if flip:
3739            angler += 180
3740
3741        _, x1, y0, y1, _, _ = br.bounds()
3742        if comment:
3743            just = "center-top"
3744            if angle is None:
3745                angle = -angler + 90
3746                if not flip:
3747                    angle += 180
3748
3749            if flip:
3750                angle += 180
3751                just = "center-bottom"
3752            if justify is not None:
3753                just = justify
3754            cmt = Text3D(comment, font=font, justify=just, italic=italic)
3755            cx0, cx1 = cmt.xbounds()
3756            cmt.rotate_z(90 + angle)
3757            cmt.scale(1 / (cx1 - cx0) * s * len(comment) / 5)
3758            cmt.shift(x1 * (1 + padding2), 0, 0)
3759            poly = merge(br, cmt).polydata()
3760
3761        else:
3762            poly = br.polydata()
3763
3764        tr = vtk.vtkTransform()
3765        tr.RotateZ(angler)
3766        tr.Translate(padding1 * d, 0, 0)
3767        pscale = 1
3768        tr.Scale(pscale / (y1 - y0) * d, pscale / (y1 - y0) * d, 1)
3769        tf = vtk.vtkTransformPolyDataFilter()
3770        tf.SetInputData(poly)
3771        tf.SetTransform(tr)
3772        tf.Update()
3773        poly = tf.GetOutput()
3774
3775        Mesh.__init__(self, poly, c, alpha)
3776        self.SetPosition(mq)
3777        self.name = "Brace"
3778        self.base = q1
3779        self.top = q2
3780
3781
3782class Star3D(Mesh):
3783    """
3784    Build a 3D starred shape.
3785    """
3786    def __init__(self, pos=(0,0,0), r=1.0, thickness=0.1, c="blue4", alpha=1):
3787        """
3788        Build a 3D star shape of 5 cusps, mainly useful as a 3D marker.
3789        """
3790        pts = ((1.34, 0., -0.37), (5.75e-3, -0.588, thickness/10), (0.377, 0.,-0.38),
3791               (0.0116, 0., -1.35), (-0.366, 0., -0.384), (-1.33, 0., -0.385),
3792               (-0.600, 0., 0.321), (-0.829, 0., 1.19), (-1.17e-3, 0., 0.761),
3793               (0.824, 0., 1.20), (0.602, 0., 0.328), (6.07e-3, 0.588, thickness/10))
3794        fcs = [[0, 1, 2], [0, 11,10], [2, 1, 3], [2, 11, 0], [3, 1, 4], [3, 11, 2],
3795               [4, 1, 5], [4, 11, 3], [5, 1, 6], [5, 11, 4], [6, 1, 7], [6, 11, 5],
3796               [7, 1, 8], [7, 11, 6], [8, 1, 9], [8, 11, 7], [9, 1,10], [9, 11, 8],
3797               [10,1, 0],[10,11, 9]]
3798
3799        Mesh.__init__(self, [pts, fcs], c, alpha)
3800        self.RotateX(90)
3801        self.scale(r).lighting("shiny")
3802
3803        if len(pos) == 2:
3804            pos = (pos[0], pos[1], 0)
3805        self.SetPosition(pos)
3806        self.name = "Star3D"
3807
3808
3809class Cross3D(Mesh):
3810    """
3811    Build a 3D cross shape.
3812    """
3813    def __init__(self, pos=(0, 0, 0), s=1.0, thickness=0.3, c="b", alpha=1):
3814        """
3815        Build a 3D cross shape, mainly useful as a 3D marker.
3816        """
3817        c1 = Cylinder(r=thickness * s, height=2 * s)
3818        c2 = Cylinder(r=thickness * s, height=2 * s).rotate_x(90)
3819        c3 = Cylinder(r=thickness * s, height=2 * s).rotate_y(90)
3820        poly = merge(c1, c2, c3).color(c).alpha(alpha).polydata(False)
3821        Mesh.__init__(self, poly, c, alpha)
3822
3823        if len(pos) == 2:
3824            pos = (pos[0], pos[1], 0)
3825        self.SetPosition(pos)
3826        self.name = "Cross3D"
3827
3828
3829class ParametricShape(Mesh):
3830    """
3831    A set of built-in shapes mainly for illustration purposes.
3832    """
3833    def __init__(self, name, res=51, n=25, seed=1):
3834        """
3835        A set of built-in shapes mainly for illustration purposes.
3836
3837        Name can be an integer or a string in this list:
3838            `['Boy', 'ConicSpiral', 'CrossCap', 'Dini', 'Enneper',
3839            'Figure8Klein', 'Klein', 'Mobius', 'RandomHills', 'Roman',
3840            'SuperEllipsoid', 'BohemianDome', 'Bour', 'CatalanMinimal',
3841            'Henneberg', 'Kuen', 'PluckerConoid', 'Pseudosphere']`.
3842
3843        Example:
3844            ```python
3845            from vedo import *
3846            settings.immediate_rendering = False
3847            plt = Plotter(N=18)
3848            for i in range(18):
3849                ps = ParametricShape(i).color(i)
3850                plt.at(i).show(ps, ps.name)
3851            plt.interactive().close()
3852            ```
3853            <img src="https://user-images.githubusercontent.com/32848391/69181075-bb6aae80-0b0e-11ea-92f7-d0cd3b9087bf.png" width="700">
3854        """
3855
3856        shapes = ['Boy', 'ConicSpiral', 'CrossCap', 'Enneper',
3857                  'Figure8Klein', 'Klein', 'Dini', 'Mobius', 'RandomHills', 'Roman',
3858                  'SuperEllipsoid', 'BohemianDome', 'Bour', 'CatalanMinimal',
3859                  'Henneberg', 'Kuen', 'PluckerConoid', 'Pseudosphere']
3860
3861        if isinstance(name, int):
3862            name = name % len(shapes)
3863            name = shapes[name]
3864
3865        if name == "Boy":
3866            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricBoy()
3867        elif name == "ConicSpiral":
3868            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricConicSpiral()
3869        elif name == "CrossCap":
3870            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricCrossCap()
3871        elif name == "Dini":
3872            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricDini()
3873        elif name == "Enneper":
3874            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricEnneper()
3875        elif name == "Figure8Klein":
3876            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricFigure8Klein()
3877        elif name == "Klein":
3878            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricKlein()
3879        elif name == "Mobius":
3880            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricMobius()
3881            ps.SetRadius(2.0)
3882            ps.SetMinimumV(-0.5)
3883            ps.SetMaximumV(0.5)
3884        elif name == "RandomHills":
3885            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricRandomHills()
3886            ps.AllowRandomGenerationOn()
3887            ps.SetRandomSeed(seed)
3888            ps.SetNumberOfHills(n)
3889        elif name == "Roman":
3890            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricRoman()
3891        elif name == "SuperEllipsoid":
3892            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricSuperEllipsoid()
3893            ps.SetN1(0.5)
3894            ps.SetN2(0.4)
3895        elif name == "BohemianDome":
3896            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricBohemianDome()
3897            ps.SetA(5.0)
3898            ps.SetB(1.0)
3899            ps.SetC(2.0)
3900        elif name == "Bour":
3901            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricBour()
3902        elif name == "CatalanMinimal":
3903            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricCatalanMinimal()
3904        elif name == "Henneberg":
3905            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricHenneberg()
3906        elif name == "Kuen":
3907            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricKuen()
3908            ps.SetDeltaV0(0.001)
3909        elif name == "PluckerConoid":
3910            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricPluckerConoid()
3911        elif name == "Pseudosphere":
3912            ps = vtk.vtkmodules.vtkCommonComputationalGeometry.vtkParametricPseudosphere()
3913        else:
3914            vedo.logger.error(f"unknown ParametricShape {name}")
3915            return
3916
3917        pfs = vtk.vtkParametricFunctionSource()
3918        pfs.SetParametricFunction(ps)
3919        pfs.SetUResolution(res)
3920        pfs.SetVResolution(res)
3921        pfs.SetWResolution(res)
3922        pfs.SetScalarModeToZ()
3923        pfs.Update()
3924
3925        Mesh.__init__(self, pfs.GetOutput())
3926
3927        if name != 'Kuen': self.normalize()
3928        if name == 'Dini': self.scale(0.4)
3929        if name == 'Enneper': self.scale(0.4)
3930        if name == 'ConicSpiral': self.bc('tomato')
3931        self.name = name
3932
3933
3934@lru_cache(None)
3935def _load_font(font):
3936    # print('_load_font()', font)
3937
3938    if utils.is_number(font):
3939        font = list(settings.font_parameters.keys())[int(font)]
3940
3941    if font.endswith(".npz"):  # user passed font as a local path
3942        fontfile = font
3943        font = os.path.basename(font).split(".")[0]
3944
3945    elif font.startswith("https"):  # user passed URL link, make it a path
3946        try:
3947            fontfile = vedo.io.download(font, verbose=False, force=False)
3948            font = os.path.basename(font).split(".")[0]
3949        except:
3950            vedo.logger.warning(f"font {font} not found")
3951            font = settings.default_font
3952            fontfile = os.path.join(vedo.fonts_path, font + ".npz")
3953
3954    else:  # user passed font by its standard name
3955        font = font[:1].upper() + font[1:]  # capitalize first letter only
3956        fontfile = os.path.join(vedo.fonts_path, font + ".npz")
3957
3958        if font not in settings.font_parameters.keys():
3959            font = "Normografo"
3960            vedo.logger.warning(
3961                f"Unknown font: {font}\n"
3962                f"Available 3D fonts are: "
3963                f"{list(settings.font_parameters.keys())}\n"
3964                f"Using font {font} instead."
3965            )
3966            fontfile = os.path.join(vedo.fonts_path, font + ".npz")
3967
3968        if not settings.font_parameters[font]["islocal"]:
3969            font = "https://vedo.embl.es/fonts/" + font + ".npz"
3970            try:
3971                fontfile = vedo.io.download(font, verbose=False, force=False)
3972                font = os.path.basename(font).split(".")[0]
3973            except:
3974                vedo.logger.warning(f"font {font} not found")
3975                font = settings.default_font
3976                fontfile = os.path.join(vedo.fonts_path, font + ".npz")
3977
3978    #####
3979    try:
3980        font_meshes = np.load(fontfile, allow_pickle=True)["font"][0]
3981    except:
3982        vedo.logger.warning(f"font name {font} not found.")
3983        raise RuntimeError
3984    return font_meshes
3985
3986
3987@lru_cache(None)
3988def _get_font_letter(font, letter):
3989    # print("_get_font_letter", font, letter)
3990    font_meshes = _load_font(font)
3991    try:
3992        pts, faces = font_meshes[letter]
3993        return utils.buildPolyData(pts, faces)
3994    except KeyError:
3995        return None
3996
3997
3998class Text3D(Mesh):
3999    """
4000    Generate a 3D polygonal Mesh to represent a text string.
4001    """
4002    def __init__(
4003        self,
4004        txt,
4005        pos=(0, 0, 0),
4006        s=1,
4007        font="",
4008        hspacing=1.15,
4009        vspacing=2.15,
4010        depth=0,
4011        italic=False,
4012        justify="bottom-left",
4013        literal=False,
4014        c=None,
4015        alpha=1,
4016    ):
4017        """
4018        Generate a 3D polygonal `Mesh` representing a text string.
4019
4020        Can render strings like `3.7 10^9` or `H_2 O` with subscripts and superscripts.
4021        Most Latex symbols are also supported.
4022
4023        Symbols `~ ^ _` are reserved modifiers:
4024        - use ~ to add a short space, 1/4 of the default empty space,
4025        - use ^ and _ to start up/sub scripting, a space terminates their effect.
4026
4027        Monospaced fonts are: `Calco, ComicMono, Glasgo, SmartCouric, VictorMono, Justino`.
4028
4029        More fonts at: https://vedo.embl.es/fonts/
4030
4031        Arguments:
4032            pos : (list)
4033                position coordinates in 3D space
4034            s : (float)
4035                size of the text
4036            depth : (float)
4037                text thickness (along z)
4038            italic : (bool), float
4039                italic font type (can be a signed float too)
4040            justify : (str)
4041                text justification as centering of the bounding box
4042                (bottom-left, bottom-right, top-left, top-right, centered)
4043            font : (str, int)
4044                some of the available 3D-polygonized fonts are:
4045                Bongas, Calco, Comae, ComicMono, Kanopus, Glasgo, Ubuntu,
4046                LogoType, Normografo, Quikhand, SmartCouric, Theemim, VictorMono, VTK,
4047                Capsmall, Cartoons123, Vega, Justino, Spears, Meson.
4048
4049                Check for more at https://vedo.embl.es/fonts/
4050
4051                Or type in your terminal `vedo --run fonts`.
4052
4053                Default is Normografo, which can be changed using `settings.default_font`.
4054
4055            hspacing : (float)
4056                horizontal spacing of the font
4057            vspacing : (float)
4058                vertical spacing of the font for multiple lines text
4059            literal : (bool)
4060                if set to True will ignore modifiers like _ or ^
4061
4062        Examples:
4063            - [markpoint.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/markpoint.py)
4064            - [fonts.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/fonts.py)
4065            - [caption.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/caption.py)
4066
4067            ![](https://vedo.embl.es/images/pyplot/fonts3d.png)
4068
4069        .. note:: Type `vedo -r fonts` for a demo.
4070        """
4071        if len(pos) == 2:
4072            pos = (pos[0], pos[1], 0)
4073
4074        if c is None:  # automatic black or white
4075            pli = vedo.plotter_instance
4076            if pli and pli.renderer:
4077                c = (0.9, 0.9, 0.9)
4078                if pli.renderer.GetGradientBackground():
4079                    bgcol = pli.renderer.GetBackground2()
4080                else:
4081                    bgcol = pli.renderer.GetBackground()
4082                if np.sum(bgcol) > 1.5:
4083                    c = (0.1, 0.1, 0.1)
4084            else:
4085                c = (0.6, 0.6, 0.6)
4086
4087        tpoly = self._get_text3d_poly(
4088            txt,
4089            s,
4090            font,
4091            hspacing,
4092            vspacing,
4093            depth,
4094            italic,
4095            justify,
4096            literal)
4097
4098        Mesh.__init__(self, tpoly, c, alpha)
4099        self.lighting("off")
4100        self.SetPosition(pos)
4101        self.PickableOff()
4102        self.DragableOff()
4103        self.name = "Text3D"
4104        self.txt = txt
4105
4106    def text(
4107        self,
4108        txt=None,
4109        s=1,
4110        font="",
4111        hspacing=1.15,
4112        vspacing=2.15,
4113        depth=0,
4114        italic=False,
4115        justify="bottom-left",
4116        literal=False,
4117    ):
4118        """
4119        Update the font style of the text.
4120        Check [available fonts here](https://vedo.embl.es/fonts).
4121        """
4122        if txt is None:
4123            return self.txt
4124
4125        tpoly = self._get_text3d_poly(
4126            txt,
4127            s,
4128            font,
4129            hspacing,
4130            vspacing,
4131            depth,
4132            italic,
4133            justify,
4134            literal,
4135        )
4136        self._update(tpoly)
4137        self.txt = txt
4138        return self
4139
4140    def _get_text3d_poly(
4141        self,
4142        txt,
4143        s=1,
4144        font="",
4145        hspacing=1.15,
4146        vspacing=2.15,
4147        depth=0,
4148        italic=False,
4149        justify="bottom-left",
4150        literal=False,
4151    ):
4152        if not font:
4153            font = settings.default_font
4154
4155        txt = str(txt)
4156
4157        if font == "VTK":  #######################################
4158            vtt = vtk.vtkVectorText()
4159            vtt.SetText(txt)
4160            vtt.Update()
4161            tpoly = vtt.GetOutput()
4162
4163        else:  ###################################################
4164
4165            stxt = set(txt)  # check here if null or only spaces
4166            if not txt or (len(stxt) == 1 and " " in stxt):
4167                return vtk.vtkPolyData()
4168
4169            if italic is True:
4170                italic = 1
4171
4172            if isinstance(font, int):
4173                lfonts = list(settings.font_parameters.keys())
4174                font = font % len(lfonts)
4175                font = lfonts[font]
4176
4177            if font not in settings.font_parameters.keys():
4178                fpars = settings.font_parameters["Normografo"]
4179            else:
4180                fpars = settings.font_parameters[font]
4181
4182            # ad hoc adjustments
4183            mono = fpars["mono"]
4184            lspacing = fpars["lspacing"]
4185            hspacing *= fpars["hspacing"]
4186            fscale = fpars["fscale"]
4187            dotsep = fpars["dotsep"]
4188
4189            # replacements
4190            if ":" in txt:
4191                for r in _reps:
4192                    txt = txt.replace(r[0], r[1])
4193
4194            if not literal:
4195                reps2 = [
4196                    ("\_", "┭"),  # trick to protect ~ _ and ^ chars
4197                    ("\^", "┮"),  #
4198                    ("\~", "┯"),  #
4199                    ("**", "^"),  # order matters
4200                    ("e+0", dotsep + "10^"),
4201                    ("e-0", dotsep + "10^-"),
4202                    ("E+0", dotsep + "10^"),
4203                    ("E-0", dotsep + "10^-"),
4204                    ("e+", dotsep + "10^"),
4205                    ("e-", dotsep + "10^-"),
4206                    ("E+", dotsep + "10^"),
4207                    ("E-", dotsep + "10^-"),
4208                ]
4209                for r in reps2:
4210                    txt = txt.replace(r[0], r[1])
4211
4212            xmax, ymax, yshift, scale = 0, 0, 0, 1
4213            save_xmax = 0
4214
4215            notfounds = set()
4216            polyletters = []
4217            ntxt = len(txt)
4218            for i, t in enumerate(txt):
4219                ##########
4220                if t == "┭":
4221                    t = "_"
4222                elif t == "┮":
4223                    t = "^"
4224                elif t == "┯":
4225                    t = "~"
4226                elif t == "^" and not literal:
4227                    if yshift < 0:
4228                        xmax = save_xmax
4229                    yshift = 0.9 * fscale
4230                    scale = 0.5
4231                    continue
4232                elif t == "_" and not literal:
4233                    if yshift > 0:
4234                        xmax = save_xmax
4235                    yshift = -0.3 * fscale
4236                    scale = 0.5
4237                    continue
4238                elif (t in (' ', '\\n')) and yshift:
4239                    yshift = 0
4240                    scale = 1
4241                    save_xmax = xmax
4242                    if t == " ":
4243                        continue
4244                elif t == "~":
4245                    if i < ntxt - 1 and txt[i + 1] == "_":
4246                        continue
4247                    xmax += hspacing * scale * fscale / 4
4248                    continue
4249
4250                ############
4251                if t == " ":
4252                    xmax += hspacing * scale * fscale
4253
4254                elif t == "\n":
4255                    xmax = 0
4256                    save_xmax = 0
4257                    ymax -= vspacing
4258
4259                else:
4260                    poly = _get_font