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