vedo.pyplot

Advanced plotting functionalities.

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import warnings
   4import numpy as np
   5
   6try:
   7    import vedo.vtkclasses as vtk
   8except ImportError:
   9    import vtkmodules.all as vtk
  10
  11import vedo
  12from vedo import settings
  13from vedo import addons
  14from vedo import colors
  15from vedo import utils
  16from vedo.pointcloud import merge
  17from vedo.mesh import Mesh
  18from vedo.assembly import Assembly, Group
  19from vedo import shapes
  20
  21
  22__docformat__ = "google"
  23
  24__doc__ = """
  25Advanced plotting functionalities.
  26
  27![](https://vedo.embl.es/images/pyplot/fitPolynomial2.png)
  28"""
  29
  30__all__ = [
  31    "Figure",
  32    "Histogram1D",  # uncomment to generate docs
  33    "Histogram2D",
  34    "PlotXY",
  35    "PlotBars",
  36    "plot",
  37    "histogram",
  38    "fit",
  39    "donut",
  40    "violin",
  41    "whisker",
  42    "streamplot",
  43    "matrix",
  44    "DirectedGraph",
  45]
  46
  47##########################################################################
  48def _to2d(actor, offset, z, scale):
  49
  50    poly = actor.polydata()
  51
  52    tp = vtk.vtkTransformPolyDataFilter()
  53    transform = vtk.vtkTransform()
  54    transform.Scale(scale, scale, scale)
  55    transform.Translate(-offset[0], -offset[1], 0)
  56    tp.SetTransform(transform)
  57    tp.SetInputData(poly)
  58    tp.Update()
  59    poly = tp.GetOutput()
  60
  61    mapper2d = vtk.vtkPolyDataMapper2D()
  62    mapper2d.SetInputData(poly)
  63
  64    act2d = vtk.vtkActor2D()
  65    act2d.SetMapper(mapper2d)
  66
  67    act2d.GetProperty().SetColor(actor.color())
  68    act2d.GetProperty().SetOpacity(actor.alpha())
  69    act2d.GetProperty().SetLineWidth(actor.GetProperty().GetLineWidth())
  70    act2d.GetProperty().SetPointSize(actor.GetProperty().GetPointSize())
  71
  72    act2d.PickableOff()
  73
  74    csys = act2d.GetPositionCoordinate()
  75    csys.SetCoordinateSystem(4)
  76
  77    act2d.GetProperty().SetDisplayLocationToBackground()
  78    return act2d
  79
  80##########################################################################
  81class LabelData:
  82    """Helper internal class to hold label information."""
  83    def __init__(self):
  84        """Helper internal class to hold label information."""
  85        self.text = "dataset"
  86        self.tcolor = "black"
  87        self.marker = "s"
  88        self.mcolor = "black"
  89
  90
  91##########################################################################
  92class Figure(Assembly):
  93    """Format class for figures."""
  94    def __init__(
  95            self,
  96            xlim,
  97            ylim,
  98            aspect=4/3,
  99            padding=(0.05, 0.05, 0.05, 0.05),
 100            **kwargs,
 101        ):
 102        """
 103        Create an empty formatted figure for plotting.
 104
 105        Arguments:
 106            xlim : (list)
 107                range of the x-axis as [x0, x1]
 108            ylim : (list)
 109                range of the y-axis as [y0, y1]
 110            aspect : (float, str)
 111                the desired aspect ratio of the histogram. Default is 4/3.
 112                Use `aspect="equal"` to force the same units in x and y.
 113            padding : (float, list)
 114                keep a padding space from the axes (as a fraction of the axis size).
 115                This can be a list of four numbers.
 116            xtitle : (str)
 117                title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
 118            ytitle : (str)
 119                title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
 120            grid : (bool)
 121                show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
 122            axes : (dict)
 123                an extra dictionary of options for the `vedo.addons.Axes` object
 124        """
 125
 126        self.xlim = np.asarray(xlim)
 127        self.ylim = np.asarray(ylim)
 128        self.aspect = aspect
 129        self.padding = padding
 130        if not utils.is_sequence(self.padding):
 131            self.padding = [self.padding, self.padding, self.padding, self.padding]
 132
 133        self.force_scaling_types = (
 134            shapes.Glyph,
 135            shapes.Line,
 136            shapes.Rectangle,
 137            shapes.DashedLine,
 138            shapes.Tube,
 139            shapes.Ribbon,
 140            shapes.GeoCircle,
 141            shapes.Arc,
 142            shapes.Grid,
 143            # shapes.Arrows, # todo
 144            # shapes.Arrows2D, # todo
 145            shapes.Brace,  # todo
 146        )
 147
 148        options = dict(kwargs)
 149
 150        self.title  = options.pop("title", "")
 151        self.xtitle = options.pop("xtitle", " ")
 152        self.ytitle = options.pop("ytitle", " ")
 153        number_of_divisions = 6
 154
 155        self.legend = None
 156        self.labels = []
 157        self.label = options.pop("label", None)
 158        if self.label:
 159            self.labels = [self.label]
 160
 161        self.axopts = options.pop("axes", {})
 162        if isinstance(self.axopts, (bool, int, float)):
 163            if self.axopts:
 164                self.axopts = {}
 165        if self.axopts or isinstance(self.axopts, dict):
 166            number_of_divisions = self.axopts.pop("number_of_divisions", number_of_divisions)
 167
 168            self.axopts["xtitle"] = self.xtitle
 169            self.axopts["ytitle"] = self.ytitle
 170
 171            if "xygrid" not in self.axopts:  ## modify the default
 172                self.axopts["xygrid"] = options.pop("grid", False)
 173
 174            if "xygrid_transparent" not in self.axopts:  ## modify the default
 175                self.axopts["xygrid_transparent"] = True
 176
 177            if "xtitle_position" not in self.axopts:  ## modify the default
 178                self.axopts["xtitle_position"] = 0.5
 179                self.axopts["xtitle_justify"] = "top-center"
 180
 181            if "ytitle_position" not in self.axopts:  ## modify the default
 182                self.axopts["ytitle_position"] = 0.5
 183                self.axopts["ytitle_justify"] = "bottom-center"
 184
 185            if self.label:
 186                if "c" in self.axopts:
 187                    self.label.tcolor = self.axopts["c"]
 188
 189        x0, x1 = self.xlim
 190        y0, y1 = self.ylim
 191        dx = x1 - x0
 192        dy = y1 - y0
 193        x0lim, x1lim = (x0 - self.padding[0] * dx, x1 + self.padding[1] * dx)
 194        y0lim, y1lim = (y0 - self.padding[2] * dy, y1 + self.padding[3] * dy)
 195        dy = y1lim - y0lim
 196
 197        self.axes = None
 198        if xlim[0] >= xlim[1] or ylim[0] >= ylim[1]:
 199            vedo.logger.warning(f"Null range for Figure {self.title}... returning an empty Assembly.")
 200            Assembly.__init__(self)
 201            self.yscale = 0
 202            return
 203
 204        if aspect == "equal":
 205            self.aspect = dx / dy  # so that yscale becomes 1
 206
 207        self.yscale = dx / dy / self.aspect
 208
 209        y0lim *= self.yscale
 210        y1lim *= self.yscale
 211
 212        self.x0lim = x0lim
 213        self.x1lim = x1lim
 214        self.y0lim = y0lim
 215        self.y1lim = y1lim
 216
 217        self.ztolerance = options.pop("ztolerance", None)
 218        if self.ztolerance is None:
 219            self.ztolerance = dx / 5000
 220
 221        ############## create axes
 222        if self.axopts:
 223            axes_opts = self.axopts
 224            if self.axopts is True or self.axopts == 1:
 225                axes_opts = {}
 226
 227            tp, ts = utils.make_ticks(
 228                y0lim / self.yscale,
 229                y1lim / self.yscale,
 230                number_of_divisions,
 231            )
 232            labs = []
 233            for i in range(1, len(tp) - 1):
 234                ynew = utils.lin_interpolate(tp[i], [0, 1], [y0lim, y1lim])
 235                labs.append([ynew, ts[i]])
 236
 237            if self.title:
 238                axes_opts["htitle"] = self.title
 239            axes_opts["y_values_and_labels"] = labs
 240            axes_opts["xrange"] = (x0lim, x1lim)
 241            axes_opts["yrange"] = (y0lim, y1lim)
 242            axes_opts["zrange"] = (0, 0)
 243            axes_opts["y_use_bounds"] = True
 244
 245            if "c" not in axes_opts and "ac" in options:
 246                axes_opts["c"] = options["ac"]
 247
 248            self.axes = addons.Axes(**axes_opts)
 249
 250        Assembly.__init__(self, [self.axes])
 251        self.name = "Figure"
 252
 253        vedo.last_figure = self if settings.remember_last_figure_format else None
 254        return
 255
 256    def _repr_html_(self):
 257        """
 258        HTML representation of the Figure object for Jupyter Notebooks.
 259        
 260        Returns:
 261            HTML text with the image and some properties.
 262        """
 263        import io
 264        import base64
 265        from PIL import Image
 266
 267        library_name = "vedo.pyplot.Figure"
 268        help_url = "https://vedo.embl.es/docs/vedo/pyplot.html#Figure"
 269
 270        arr = self.thumbnail(zoom=1.1)
 271
 272        im = Image.fromarray(arr)
 273        buffered = io.BytesIO()
 274        im.save(buffered, format="PNG", quality=100)
 275        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 276        url = "data:image/png;base64," + encoded
 277        image = f"<img src='{url}'></img>"
 278
 279        # statisitics
 280        bounds = "<br/>".join(
 281            [
 282                vedo.utils.precision(min_x,4) + " ... " + vedo.utils.precision(max_x,4)
 283                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 284            ]
 285        )
 286
 287        help_text = ""
 288        if self.name:
 289            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 290        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 
 291        if self.filename:
 292            dots = ""
 293            if len(self.filename) > 30:
 294                dots = "..."
 295            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 296        
 297        all = [
 298            "<table>",
 299            "<tr>", 
 300            "<td>", image, "</td>",
 301            "<td style='text-align: center; vertical-align: center;'><br/>", help_text,
 302            "<table>",
 303            "<tr><td><b> nr. of objects </b></td><td>" + str(self.GetNumberOfPaths()) + "</td></tr>",
 304            "<tr><td><b> position </b></td><td>" + str(self.GetPosition()) + "</td></tr>",
 305            "<tr><td><b> x-limits </b></td><td>" + utils.precision(self.xlim,4) + "</td></tr>",
 306            "<tr><td><b> y-limits </b></td><td>" + utils.precision(self.ylim,4) + "</td></tr>",
 307            "<tr><td><b> world bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 308            "</table>",
 309            "</table>",
 310        ]
 311        return "\n".join(all)
 312
 313    def __add__(self, *obj):
 314        # just to avoid confusion, supersede Assembly.__add__
 315        return self.__iadd__(*obj)
 316
 317    def __iadd__(self, *obj):
 318        if len(obj) == 1 and isinstance(obj[0], Figure):
 319            return self._check_unpack_and_insert(obj[0])
 320
 321        obj = utils.flatten(obj)
 322        return self.insert(*obj)
 323
 324    def _check_unpack_and_insert(self, fig):
 325
 326        if fig.label:
 327            self.labels.append(fig.label)
 328
 329        if abs(self.yscale - fig.yscale) > 0.0001:
 330
 331            colors.printc("ERROR: adding incompatible Figure. Y-scales are different:",
 332                          c='r', invert=True)
 333            colors.printc("  first  figure:", self.yscale, c='r')
 334            colors.printc("  second figure:", fig.yscale, c='r')
 335
 336            colors.printc("One or more of these parameters can be the cause:", c="r")
 337            if list(self.xlim) != list(fig.xlim):
 338                colors.printc("xlim --------------------------------------------\n",
 339                              " first  figure:", self.xlim, "\n",
 340                              " second figure:", fig.xlim, c='r')
 341            if list(self.ylim) != list(fig.ylim):
 342                colors.printc("ylim --------------------------------------------\n",
 343                              " first  figure:", self.ylim, "\n",
 344                              " second figure:", fig.ylim, c='r')
 345            if list(self.padding) != list(fig.padding):
 346                colors.printc("padding -----------------------------------------\n",
 347                              " first  figure:", self.padding,
 348                              " second figure:", fig.padding, c='r')
 349            if self.aspect != fig.aspect:
 350                colors.printc("aspect ------------------------------------------\n",
 351                              " first  figure:", self.aspect,
 352                              " second figure:", fig.aspect, c='r')
 353
 354            colors.printc("\n*Consider using fig2 = histogram(..., like=fig1)", c="r")
 355            colors.printc(" Or fig += histogram(..., like=fig)\n", c="r")
 356            return self
 357
 358        offset = self.zbounds()[1] + self.ztolerance
 359
 360        for ele in fig.unpack():
 361            if "Axes" in ele.name:
 362                continue
 363            ele.z(offset)
 364            self.insert(ele, rescale=False)
 365
 366        return self
 367
 368    def insert(self, *objs, rescale=True, as3d=True, adjusted=True, cut=True):
 369        """
 370        Insert objects into a Figure.
 371
 372        The recommended syntax is to use "+=", which calls `insert()` under the hood.
 373        If a whole Figure is added with "+=", it is unpacked and its objects are added
 374        one by one.
 375
 376        Arguments:
 377            rescale : (bool)
 378                rescale the y axis position while inserting the object.
 379            as3d : (bool)
 380                if True keep the aspect ratio of the 3d object, otherwise stretch it in y.
 381            adjusted : (bool)
 382                adjust the scaling according to the shortest axis
 383            cut : (bool)
 384                cut off the parts of the object which go beyond the axes frame.
 385        """
 386        for a in objs:
 387
 388            if a in self.actors:
 389                # should not add twice the same object in plot
 390                continue
 391
 392            if isinstance(a, vedo.Points):  # hacky way to identify Points
 393                if a.ncells == a.npoints:
 394                    poly = a.polydata(False)
 395                    if poly.GetNumberOfPolys() == 0 and poly.GetNumberOfLines() == 0:
 396                        as3d = False
 397                        rescale = True
 398
 399            if isinstance(a, (shapes.Arrow, shapes.Arrow2D)):
 400                # discard input Arrow and substitute it with a brand new one
 401                # (because scaling would fatally distort the shape)
 402                prop = a.GetProperty()
 403                prop.LightingOff()
 404                py = a.base[1]
 405                a.top[1] = (a.top[1] - py) * self.yscale + py
 406                b = shapes.Arrow2D(a.base, a.top, s=a.s, fill=a.fill).z(a.z())
 407                b.SetProperty(prop)
 408                b.y(py * self.yscale)
 409                a = b
 410
 411            # elif isinstance(a, shapes.Rectangle) and a.radius is not None:
 412            #     # discard input Rectangle and substitute it with a brand new one
 413            #     # (because scaling would fatally distort the shape of the corners)
 414            #     py = a.corner1[1]
 415            #     rx1,ry1,rz1 = a.corner1
 416            #     rx2,ry2,rz2 = a.corner2
 417            #     ry2 = (ry2-py) * self.yscale + py
 418            #     b = shapes.Rectangle([rx1,0,rz1], [rx2,ry2,rz2], radius=a.radius).z(a.z())
 419            #     b.SetProperty(a.GetProperty())
 420            #     b.y(py / self.yscale)
 421            #     a = b
 422
 423            else:
 424
 425                if rescale:
 426
 427                    if not isinstance(a, Figure):
 428
 429                        if as3d and not isinstance(a, self.force_scaling_types):
 430                            if adjusted:
 431                                scl = np.min([1, self.yscale])
 432                            else:
 433                                scl = self.yscale
 434
 435                            a.scale(scl)
 436
 437                        else:
 438                            a.scale([1, self.yscale, 1])
 439
 440                    # shift it in y
 441                    a.y(a.y() * self.yscale)
 442
 443            if cut:
 444                try:
 445                    bx0, bx1, by0, by1, _, _ = a.bounds()
 446                    if self.y0lim > by0:
 447                        a.cut_with_plane([0, self.y0lim, 0], [ 0, 1, 0])
 448                    if self.y1lim < by1:
 449                        a.cut_with_plane([0, self.y1lim, 0], [ 0,-1, 0])
 450                    if self.x0lim > bx0:
 451                        a.cut_with_plane([self.x0lim, 0, 0], [ 1, 0, 0])
 452                    if self.x1lim < bx1:
 453                        a.cut_with_plane([self.x1lim, 0, 0], [-1, 0, 0])
 454                except:
 455                    # print("insert(): cannot cut", [a])
 456                    pass
 457
 458            self.AddPart(a)
 459            self.actors.append(a)
 460
 461        return self
 462
 463    def add_label(self, text, c=None, marker="", mc="black"):
 464        """
 465        Manually add en entry label to the legend.
 466
 467        Arguments:
 468            text : (str)
 469                text string for the label.
 470            c : (str)
 471                color of the text
 472            marker : (str), Mesh
 473                a marker char or a Mesh object to be used as marker
 474            mc : (str)
 475                color for the marker
 476        """
 477        newlabel = LabelData()
 478        newlabel.text = text.replace("\n", " ")
 479        newlabel.tcolor = c
 480        newlabel.marker = marker
 481        newlabel.mcolor = mc
 482        self.labels.append(newlabel)
 483        return self
 484
 485
 486    def add_legend(self,
 487            pos="top-right",
 488            relative=True,
 489            font=None,
 490            s=1,
 491            c=None,
 492            vspace=1.75,
 493            padding=0.1,
 494            radius=0,
 495            alpha=1,
 496            bc='k7',
 497            lw=1,
 498            lc='k4',
 499            z=0,
 500        ):
 501        """
 502        Add existing labels to form a legend box.
 503        Labels have been previously filled with eg: `plot(..., label="text")`
 504
 505        Arguments:
 506            pos : (str, list)
 507                A string or 2D coordinates. The default is "top-right".
 508            relative : (bool)
 509                control whether `pos` is absolute or relative, e.i. normalized
 510                to the x and y ranges so that x and y in `pos=[x,y]` should be
 511                both in the range [0,1].
 512                This flag is ignored if a string despcriptor is passed.
 513                Default is True.
 514            font : (str, int)
 515                font name or number.
 516                Check [available fonts here](https://vedo.embl.es/fonts).
 517            s : (float)
 518                global size of the legend
 519            c : (str)
 520                color of the text
 521            vspace : (float)
 522                vertical spacing of lines
 523            padding : (float)
 524                padding of the box as a fraction of the text size
 525            radius : (float)
 526                border radius of the box
 527            alpha : (float)
 528                opacity of the box. Values below 1 may cause poor rendering
 529                because of antialiasing.
 530                Use alpha = 0 to remove the box.
 531            bc : (str)
 532                box color
 533            lw : (int)
 534                border line width of the box in pixel units
 535            lc : (int)
 536                border line color of the box
 537            z : (float)
 538                set the zorder as z position (useful to avoid overlap)
 539        """
 540        sx = self.x1lim - self.x0lim
 541        s = s * sx / 55  # so that input can be about 1
 542
 543        ds = 0
 544        texts = []
 545        mks = []
 546        for i, t in enumerate(self.labels):
 547            label = self.labels[i]
 548            t = label.text
 549
 550            if label.tcolor is not None:
 551                c = label.tcolor
 552
 553            tx = vedo.shapes.Text3D(t, s=s, c=c, justify="center-left", font=font)
 554            y0, y1 = tx.ybounds()
 555            ds = max(y1 - y0, ds)
 556            texts.append(tx)
 557
 558            mk = label.marker
 559            if isinstance(mk, vedo.Points):
 560                mk = mk.clone(deep=False).lighting("off")
 561                cm = mk.center_of_mass()
 562                ty0, ty1 = tx.ybounds()
 563                oby0, oby1 = mk.ybounds()
 564                mk.shift(-cm)
 565                mk.SetOrigin(cm)
 566                mk.scale((ty1 - ty0) / (oby1 - oby0))
 567                mk.scale([1.1, 1.1, 0.01])
 568            elif mk == "-":
 569                mk = vedo.shapes.Marker(mk, s=s * 2)
 570                mk.color(label.mcolor)
 571            else:
 572                mk = vedo.shapes.Marker(mk, s=s)
 573                mk.color(label.mcolor)
 574            mks.append(mk)
 575
 576        for i, tx in enumerate(texts):
 577            tx.shift(0, -(i + 0) * ds * vspace)
 578
 579        for i, mk in enumerate(mks):
 580            mk.shift(-ds * 1.75, -(i + 0) * ds * vspace, 0)
 581
 582        acts = texts + mks
 583
 584        aleg = Assembly(acts)  # .show(axes=1).close()
 585        x0, x1, y0, y1, _, _ = aleg.GetBounds()
 586
 587        if alpha:
 588            dx = x1 - x0
 589            dy = y1 - y0
 590
 591            if not utils.is_sequence(padding):
 592                padding = [padding] * 4
 593            padding = min(padding)
 594            padding = min(padding * dx, padding * dy)
 595            if len(self.labels) == 1:
 596                padding *= 4
 597            x0 -= padding
 598            x1 += padding
 599            y0 -= padding
 600            y1 += padding
 601
 602            box = shapes.Rectangle(
 603                [x0, y0], [x1, y1], radius=radius, c=bc, alpha=alpha,
 604            )
 605            box.shift(0, 0, -dy / 100).pickable(False)
 606            if lc:
 607                box.lc(lc).lw(lw)
 608            aleg.AddPart(box)
 609
 610        xlim = self.xlim
 611        ylim = self.ylim
 612        if isinstance(pos, str):
 613            px, py = 0, 0
 614            rx, ry = (xlim[1] + xlim[0]) / 2, (ylim[1] + ylim[0]) / 2
 615            shx, shy = 0, 0
 616            if "top" in pos:
 617                if "cent" in pos:
 618                    px, py = rx, ylim[1]
 619                    shx, shy = (x0 + x1) / 2, y1
 620                elif "left" in pos:
 621                    px, py = xlim[0], ylim[1]
 622                    shx, shy = x0, y1
 623                else:  # "right"
 624                    px, py = xlim[1], ylim[1]
 625                    shx, shy = x1, y1
 626            elif "bot" in pos:
 627                if "left" in pos:
 628                    px, py = xlim[0], ylim[0]
 629                    shx, shy = x0, y0
 630                elif "right" in pos:
 631                    px, py = xlim[1], ylim[0]
 632                    shx, shy = x1, y0
 633                else:  # "cent"
 634                    px, py = rx, ylim[0]
 635                    shx, shy = (x0 + x1) / 2, y0
 636            elif "cent" in pos:
 637                if "left" in pos:
 638                    px, py = xlim[0], ry
 639                    shx, shy = x0, (y0 + y1) / 2
 640                elif "right" in pos:
 641                    px, py = xlim[1], ry
 642                    shx, shy = x1, (y0 + y1) / 2
 643            else:
 644                vedo.logger.error(f"in add_legend(), cannot understand {pos}")
 645                raise RuntimeError
 646
 647        else:
 648
 649            if relative:
 650                rx, ry = pos[0], pos[1]
 651                px = (xlim[1] - xlim[0]) * rx + xlim[0]
 652                py = (ylim[1] - ylim[0]) * ry + ylim[0]
 653                z *= xlim[1] - xlim[0]
 654            else:
 655                px, py = pos[0], pos[1]
 656            shx, shy = x0, y1
 657
 658        aleg.pos(px - shx, py * self.yscale - shy, self.z() + sx / 50 + z)
 659
 660        self.insert(aleg, rescale=False, cut=False)
 661        self.legend = aleg
 662        aleg.name = "Legend"
 663        return self
 664
 665    def as2d(self, pos="bottom-left", scale=1, padding=0.05):
 666        """
 667        Convert the Figure into a 2D static object (a 2D assembly).
 668
 669        Still experimental.
 670
 671        Arguments:
 672            pos : (str, list)
 673                position in 2D, as atring or list (x,y).
 674                Any combination of "center", "top", "bottom", "left" and "right" will work.
 675                The center of the renderer is [0,0] while top-right is [1,1].
 676            scale : (float)
 677                scaling factor
 678            padding : (float, list)
 679                a single value or a list (xpad, ypad)
 680
 681        Returns:
 682            `Group` object.
 683        """
 684        x0, x1 = self.xbounds()
 685        y0, y1 = self.ybounds()
 686        pp = self.pos()
 687        x0 -= pp[0]
 688        x1 -= pp[0]
 689        y0 -= pp[1]
 690        y1 -= pp[1]
 691
 692        if not utils.is_sequence(padding):
 693            padding = (padding, padding)
 694        padding = np.array(padding)
 695
 696        def _unpack(p):
 697            for o in list(objs):
 698                if isinstance(o, Assembly):
 699                    objs.extend(o.unpack())
 700        objs = [self]
 701        for _ in range(3):
 702            _unpack(self)
 703
 704        if "cent" in pos:
 705            offset = [(x0 + x1) / 2, (y0 + y1) / 2]
 706            position = [0, 0]
 707            if "right" in pos:
 708                offset[0] = x1
 709                position = [1 - padding[0], 0]
 710            if "left" in pos:
 711                offset[0] = x0
 712                position = [-1 + padding[0], 0]
 713            if "top" in pos:
 714                offset[1] = y1
 715                position = [0, 1 - padding[1]]
 716            if "bottom" in pos:
 717                offset[1] = y0
 718                position = [0, -1 + padding[1]]
 719        elif "top" in pos:
 720            if "right" in pos:
 721                offset = [x1, y1]
 722                position = [1, 1] - padding
 723            elif "left" in pos:
 724                offset = [x0, y1]
 725                position = [-1 + padding[0], 1 - padding[1]]
 726            else:
 727                raise ValueError(f"incomplete position pos='{pos}'")
 728        elif "bottom" in pos:
 729            if "right" in pos:
 730                offset = [x1, y0]
 731                position = [1 - padding[0], -1 + padding[1]]
 732            elif "left" in pos:
 733                offset = [x0, y0]
 734                position = [-1, -1] + padding
 735            else:
 736                raise ValueError(f"incomplete position pos='{pos}'")
 737        else:
 738            offset = [0,0]
 739            position = pos
 740
 741        scanned = []
 742        group = Group()
 743        for a in objs:
 744            if a in scanned:
 745                continue
 746            if not isinstance(a, vedo.Points):
 747                continue
 748            if a.npoints == 0:
 749                continue
 750            if a.GetProperty().GetRepresentation() == 1:
 751                # wireframe is not renedered in 2d
 752                continue
 753            aa = _to2d(a, offset, self.z(), scale / (x1-x0) * 550)
 754            aa.SetPosition(position)
 755            group += aa
 756        return group
 757
 758
 759#########################################################################################
 760class Histogram1D(Figure):
 761    "1D histogramming."
 762
 763    def __init__(
 764            self,
 765            data,
 766            weights=None,
 767            bins=None,
 768            errors=False,
 769            density=False,
 770            logscale=False,
 771            fill=True,
 772            radius=0.075,
 773            c="olivedrab",
 774            gap=0.02,
 775            alpha=1,
 776            outline=False,
 777            lw=2,
 778            lc="k",
 779            texture="",
 780            marker="",
 781            ms=None,
 782            mc=None,
 783            ma=None,
 784
 785            # Figure and axes options:
 786            like=None,
 787            xlim=None,
 788            ylim=(0, None),
 789            aspect=4/3,
 790            padding=(0., 0., 0., 0.05),
 791
 792            title="",
 793            xtitle=" ",
 794            ytitle=" ",
 795            ac="k",
 796            grid=False,
 797            ztolerance=None,
 798            label="",
 799            **fig_kwargs,
 800        ):
 801        """
 802        Creates a `Histogram1D(Figure)` object.
 803
 804        Arguments:
 805            weights : (list)
 806                An array of weights, of the same shape as `data`. Each value in `data`
 807                only contributes its associated weight towards the bin count (instead of 1).
 808            bins : (int)
 809                number of bins
 810            density : (bool)
 811                normalize the area to 1 by dividing by the nr of entries and bin size
 812            logscale : (bool)
 813                use logscale on y-axis
 814            fill : (bool)
 815                fill bars with solid color `c`
 816            gap : (float)
 817                leave a small space btw bars
 818            radius : (float)
 819                border radius of the top of the histogram bar. Default value is 0.1.
 820            texture : (str)
 821                url or path to an image to be used as texture for the bin
 822            outline : (bool)
 823                show outline of the bins
 824            errors : (bool)
 825                show error bars
 826            xtitle : (str)
 827                title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
 828            ytitle : (str)
 829                title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
 830            padding : (float), list
 831                keep a padding space from the axes (as a fraction of the axis size).
 832                This can be a list of four numbers.
 833            aspect : (float)
 834                the desired aspect ratio of the histogram. Default is 4/3.
 835            grid : (bool)
 836                show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
 837            ztolerance : (float)
 838                a tolerance factor to superimpose objects (along the z-axis).
 839
 840        Examples:
 841            - [histo_1d_a.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_1d_a.py)
 842            - [histo_1d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_1d_b.py)
 843            - [histo_1d_c.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_1d_c.py)
 844            - [histo_1d_d.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_1d_d.py)
 845
 846            ![](https://vedo.embl.es/images/pyplot/histo_1D.png)
 847        """
 848
 849        # purge NaN from data
 850        valid_ids = np.all(np.logical_not(np.isnan(data)))
 851        data = np.asarray(data[valid_ids]).ravel()
 852
 853        # if data.dtype is integer try to center bins by default
 854        if like is None and bins is None and np.issubdtype(data.dtype, np.integer):
 855            if xlim is None and ylim == (0, None):
 856                x1, x0 = data.max(), data.min()
 857                if 0 < x1 - x0 <= 100:
 858                    bins = x1 - x0 +1 
 859                    xlim = (x0 - 0.5, x1 + 0.5)
 860
 861        if like is None and vedo.last_figure is not None:
 862            if xlim is None and ylim == (0, None):
 863                like = vedo.last_figure
 864
 865        if like is not None:
 866            xlim = like.xlim
 867            ylim = like.ylim
 868            aspect = like.aspect
 869            padding = like.padding
 870            if bins is None:
 871                bins = like.bins
 872        if bins is None:
 873            bins = 20
 874
 875        if utils.is_sequence(xlim):
 876            # deal with user passing eg [x0, None]
 877            _x0, _x1 = xlim
 878            if _x0 is None:
 879                _x0 = data.min()
 880            if _x1 is None:
 881                _x1 = data.max()
 882            xlim = [_x0, _x1]
 883
 884        fs, edges = np.histogram(data, bins=bins, weights=weights, range=xlim)
 885        binsize = edges[1] - edges[0]
 886        ntot = data.shape[0]
 887
 888        fig_kwargs["title"] = title
 889        fig_kwargs["xtitle"] = xtitle
 890        fig_kwargs["ytitle"] = ytitle
 891        fig_kwargs["ac"] = ac
 892        fig_kwargs["ztolerance"] = ztolerance
 893        fig_kwargs["grid"] = grid
 894
 895        unscaled_errors = np.sqrt(fs)
 896        if density:
 897            scaled_errors = unscaled_errors / (ntot * binsize)
 898            fs = fs / (ntot * binsize)
 899            if ytitle == " ":
 900                ytitle = f"counts / ({ntot} x {utils.precision(binsize,3)})"
 901                fig_kwargs["ytitle"] = ytitle
 902        elif logscale:
 903            se_up = np.log10(fs + unscaled_errors / 2 + 1)
 904            se_dw = np.log10(fs - unscaled_errors / 2 + 1)
 905            scaled_errors = np.c_[se_up, se_dw]
 906            fs = np.log10(fs + 1)
 907            if ytitle == " ":
 908                ytitle = "log_10 (counts+1)"
 909                fig_kwargs["ytitle"] = ytitle
 910
 911        x0, x1 = np.min(edges), np.max(edges)
 912        y0, y1 = ylim[0], np.max(fs)
 913
 914        _errors = []
 915        if errors:
 916            if density:
 917                y1 += max(scaled_errors) / 2
 918                _errors = scaled_errors
 919            elif logscale:
 920                y1 = max(scaled_errors[:, 0])
 921                _errors = scaled_errors
 922            else:
 923                y1 += max(unscaled_errors) / 2
 924                _errors = unscaled_errors
 925
 926        if like is None:
 927            ylim = list(ylim)
 928            if xlim is None:
 929                xlim = [x0, x1]
 930            if ylim[1] is None:
 931                ylim[1] = y1
 932            if ylim[0] != 0:
 933                ylim[0] = y0
 934
 935        self.title = title
 936        self.xtitle = xtitle
 937        self.ytitle = ytitle
 938        self.entries = ntot
 939        self.frequencies = fs
 940        self.errors = _errors
 941        self.edges = edges
 942        self.centers = (edges[0:-1] + edges[1:]) / 2
 943        self.mean = data.mean()
 944        self.std = data.std()
 945        self.bins = edges  # internally used by "like"
 946
 947        ############################### stats legend as htitle
 948        addstats = False
 949        if not title:
 950            if "axes" not in fig_kwargs:
 951                addstats = True
 952                axes_opts = {}
 953                fig_kwargs["axes"] = axes_opts
 954            elif fig_kwargs["axes"] is False:
 955                pass
 956            else:
 957                axes_opts = fig_kwargs["axes"]
 958                if "htitle" not in axes_opts:
 959                    addstats = True
 960
 961        if addstats:
 962            htitle = f"Entries:~~{int(self.entries)}  "
 963            htitle += f"Mean:~~{utils.precision(self.mean, 4)}  "
 964            htitle += f"STD:~~{utils.precision(self.std, 4)}  "
 965
 966            axes_opts["htitle"] = htitle
 967            axes_opts["htitle_justify"] = "bottom-left"
 968            axes_opts["htitle_size"] = 0.016
 969            axes_opts["htitle_offset"] = [-0.49, 0.01, 0]
 970
 971        if mc is None:
 972            mc = lc
 973        if ma is None:
 974            ma = alpha
 975
 976        if label:
 977            nlab = LabelData()
 978            nlab.text = label
 979            nlab.tcolor = ac
 980            nlab.marker = marker
 981            nlab.mcolor = mc
 982            if not marker:
 983                nlab.marker = "s"
 984                nlab.mcolor = c
 985            fig_kwargs["label"] = nlab
 986
 987        ############################################### Figure init
 988        Figure.__init__(self, xlim, ylim, aspect, padding, **fig_kwargs)
 989        if not self.yscale:
 990            return
 991
 992        if utils.is_sequence(bins):
 993            myedges = np.array(bins)
 994            bins = len(bins) - 1
 995        else:
 996            myedges = edges
 997
 998        bin_centers = []
 999        for i in range(bins):
1000            x = (myedges[i] + myedges[i + 1]) / 2
1001            bin_centers.append([x, fs[i], 0])
1002
1003        rs = []
1004        maxheigth = 0
1005        if not fill and not outline and not errors and not marker:
1006            outline = True  # otherwise it's empty..
1007
1008        if fill:  #####################
1009            if outline:
1010                gap = 0
1011
1012            for i in range(bins):
1013                F = fs[i]
1014                if not F:
1015                    continue
1016                p0 = (myedges[i] + gap * binsize, 0, 0)
1017                p1 = (myedges[i + 1] - gap * binsize, F, 0)
1018
1019                if radius:
1020                    if gap:
1021                        rds = np.array([0, 0, radius, radius])
1022                    else:
1023                        rd1 = 0 if i < bins - 1 and fs[i + 1] >= F else radius / 2
1024                        rd2 = 0 if i > 0 and fs[i - 1] >= F else radius / 2
1025                        rds = np.array([0, 0, rd1, rd2])
1026                    p1_yscaled = [p1[0], p1[1]*self.yscale, 0]
1027                    r = shapes.Rectangle(p0, p1_yscaled, radius=rds*binsize, res=6)
1028                    r.scale([1, 1/self.yscale, 1])
1029                    r.radius = None  # so it doesnt get recreated and rescaled by insert()
1030                else:
1031                    r = shapes.Rectangle(p0, p1)
1032
1033                if texture:
1034                    r.texture(texture)
1035                    c = "w"
1036                # if texture: # causes Segmentation fault vtk9.0.3
1037                #     if i>0 and rs[i-1].GetTexture(): # reuse the same texture obj
1038                #         r.texture(rs[i-1].GetTexture())
1039                #     else:
1040                #         r.texture(texture)
1041                #     c = 'w'
1042
1043                r.PickableOff()
1044                maxheigth = max(maxheigth, p1[1])
1045                if c in colors.cmaps_names:
1046                    col = colors.color_map((p0[0]+p1[0])/2, c, myedges[0], myedges[-1])
1047                else:
1048                    col = c
1049                r.color(col).alpha(alpha).lighting("off")
1050                r.z(self.ztolerance)
1051                rs.append(r)
1052
1053        if outline:  #####################
1054            lns = [[myedges[0], 0, 0]]
1055            for i in range(bins):
1056                lns.append([myedges[i], fs[i], 0])
1057                lns.append([myedges[i + 1], fs[i], 0])
1058                maxheigth = max(maxheigth, fs[i])
1059            lns.append([myedges[-1], 0, 0])
1060            outl = shapes.Line(lns, c=lc, alpha=alpha, lw=lw)
1061            outl.z(self.ztolerance * 2)
1062            rs.append(outl)
1063
1064        if errors:  #####################
1065            for i in range(bins):
1066                x = self.centers[i]
1067                f = fs[i]
1068                if not f:
1069                    continue
1070                err = _errors[i]
1071                if utils.is_sequence(err):
1072                    el = shapes.Line([x, err[0], 0], [x, err[1], 0], c=lc, alpha=alpha, lw=lw)
1073                else:
1074                    el = shapes.Line([x, f-err/2, 0], [x, f+err/2, 0], c=lc, alpha=alpha, lw=lw)
1075                el.z(self.ztolerance*3)
1076                rs.append(el)
1077
1078        if marker:  #####################
1079
1080            # remove empty bins (we dont want a marker there)
1081            bin_centers = np.array(bin_centers)
1082            bin_centers = bin_centers[bin_centers[:, 1] > 0]
1083
1084            if utils.is_sequence(ms):  ### variable point size
1085                mk = shapes.Marker(marker, s=1)
1086                mk.scale([1, 1 / self.yscale, 1])
1087                msv = np.zeros_like(bin_centers)
1088                msv[:, 0] = ms
1089                marked = shapes.Glyph(
1090                    bin_centers,
1091                    mk,
1092                    c=mc,
1093                    orientation_array=msv,
1094                    scale_by_vector_size=True,
1095                )
1096            else:  ### fixed point size
1097
1098                if ms is None:
1099                    ms = (xlim[1] - xlim[0]) / 100.0
1100                else:
1101                    ms = (xlim[1] - xlim[0]) / 100.0 * ms
1102
1103                if utils.is_sequence(mc):
1104                    mk = shapes.Marker(marker, s=ms)
1105                    mk.scale([1, 1 / self.yscale, 1])
1106                    msv = np.zeros_like(bin_centers)
1107                    msv[:, 0] = 1
1108                    marked = shapes.Glyph(
1109                        bin_centers,
1110                        mk,
1111                        c=mc,
1112                        orientation_array=msv,
1113                        scale_by_vector_size=True,
1114                    )
1115                else:
1116                    mk = shapes.Marker(marker, s=ms)
1117                    mk.scale([1, 1 / self.yscale, 1])
1118                    marked = shapes.Glyph(bin_centers, mk, c=mc)
1119
1120            marked.alpha(ma)
1121            marked.z(self.ztolerance * 4)
1122            rs.append(marked)
1123
1124        self.insert(*rs, as3d=False)
1125        self.name = "Histogram1D"
1126
1127    def print(self, **kwargs):
1128        """Print infos about this histogram"""
1129        txt = (
1130            f"{self.name}  {self.title}\n"
1131            f"    xtitle  = '{self.xtitle}'\n"
1132            f"    ytitle  = '{self.ytitle}'\n"
1133            f"    entries = {self.entries}\n"
1134            f"    mean    = {self.mean}\n"
1135            f"    std     = {self.std}"
1136        )
1137        colors.printc(txt, **kwargs)
1138
1139
1140#########################################################################################
1141class Histogram2D(Figure):
1142    """2D histogramming."""
1143
1144    def __init__(
1145            self,
1146            xvalues,
1147            yvalues=None,
1148            bins=25,
1149            weights=None,
1150            cmap="cividis",
1151            alpha=1,
1152            gap=0,
1153            scalarbar=True,
1154
1155            # Figure and axes options:
1156            like=None,
1157            xlim=None,
1158            ylim=(None, None),
1159            zlim=(None, None),
1160            aspect=1,
1161
1162            title="",
1163            xtitle=" ",
1164            ytitle=" ",
1165            ztitle="",
1166            ac="k",
1167            **fig_kwargs,
1168        ):
1169        """
1170        Input data formats `[(x1,x2,..), (y1,y2,..)] or [(x1,y1), (x2,y2),..]`
1171        are both valid.
1172
1173        Use keyword `like=...` if you want to use the same format of a previously
1174        created Figure (useful when superimposing Figures) to make sure
1175        they are compatible and comparable. If they are not compatible
1176        you will receive an error message.
1177
1178        Arguments:
1179            bins : (list)
1180                binning as (nx, ny)
1181            weights : (list)
1182                array of weights to assign to each entry
1183            cmap : (str, lookuptable)
1184                color map name or look up table
1185            alpha : (float)
1186                opacity of the histogram
1187            gap : (float)
1188                separation between adjacent bins as a fraction for their size
1189            scalarbar : (bool)
1190                add a scalarbar to right of the histogram
1191            like : (Figure)
1192                grab and use the same format of the given Figure (for superimposing)
1193            xlim : (list)
1194                [x0, x1] range of interest. If left to None will automatically
1195                choose the minimum or the maximum of the data range.
1196                Data outside the range are completely ignored.
1197            ylim : (list)
1198                [y0, y1] range of interest. If left to None will automatically
1199                choose the minimum or the maximum of the data range.
1200                Data outside the range are completely ignored.
1201            aspect : (float)
1202                the desired aspect ratio of the figure.
1203            title : (str)
1204                title of the plot to appear on top.
1205                If left blank some statistics will be shown.
1206            xtitle : (str)
1207                x axis title
1208            ytitle : (str)
1209                y axis title
1210            ztitle : (str)
1211                title for the scalar bar
1212            ac : (str)
1213                axes color, additional keyword for Axes can also be added
1214                using e.g. `axes=dict(xygrid=True)`
1215
1216        Examples:
1217            - [histo_2d_a.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_2d_a.py)
1218            - [histo_2d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_2d_b.py)
1219
1220            ![](https://vedo.embl.es/images/pyplot/histo_2D.png)
1221        """
1222        xvalues = np.asarray(xvalues)
1223        if yvalues is None:
1224            # assume [(x1,y1), (x2,y2) ...] format
1225            yvalues = xvalues[:, 1]
1226            xvalues = xvalues[:, 0]
1227
1228        padding = [0, 0, 0, 0]
1229
1230        if like is None and vedo.last_figure is not None:
1231            if xlim is None and ylim == (None, None) and zlim == (None, None):
1232                like = vedo.last_figure
1233
1234        if like is not None:
1235            xlim = like.xlim
1236            ylim = like.ylim
1237            aspect = like.aspect
1238            padding = like.padding
1239            if bins is None:
1240                bins = like.bins
1241        if bins is None:
1242            bins = 20
1243
1244        if isinstance(bins, int):
1245            bins = (bins, bins)
1246
1247        if utils.is_sequence(xlim):
1248            # deal with user passing eg [x0, None]
1249            _x0, _x1 = xlim
1250            if _x0 is None:
1251                _x0 = xvalues.min()
1252            if _x1 is None:
1253                _x1 = xvalues.max()
1254            xlim = [_x0, _x1]
1255
1256        if utils.is_sequence(ylim):
1257            # deal with user passing eg [x0, None]
1258            _y0, _y1 = ylim
1259            if _y0 is None:
1260                _y0 = yvalues.min()
1261            if _y1 is None:
1262                _y1 = yvalues.max()
1263            ylim = [_y0, _y1]
1264
1265        H, xedges, yedges = np.histogram2d(
1266            xvalues, yvalues, 
1267            weights=weights, bins=bins, range=(xlim, ylim),
1268        )
1269
1270        xlim = np.min(xedges), np.max(xedges)
1271        ylim = np.min(yedges), np.max(yedges)
1272        dx, dy = xlim[1] - xlim[0], ylim[1] - ylim[0]
1273
1274        fig_kwargs["title"] = title
1275        fig_kwargs["xtitle"] = xtitle
1276        fig_kwargs["ytitle"] = ytitle
1277        fig_kwargs["ac"] = ac
1278
1279        self.entries = len(xvalues)
1280        self.frequencies = H
1281        self.edges = (xedges, yedges)
1282        self.mean = (xvalues.mean(), yvalues.mean())
1283        self.std = (xvalues.std(), yvalues.std())
1284        self.bins = bins  # internally used by "like"
1285
1286        ############################### stats legend as htitle
1287        addstats = False
1288        if not title:
1289            if "axes" not in fig_kwargs:
1290                addstats = True
1291                axes_opts = {}
1292                fig_kwargs["axes"] = axes_opts
1293            elif fig_kwargs["axes"] is False:
1294                pass
1295            else:
1296                axes_opts = fig_kwargs["axes"]
1297                if "htitle" not in fig_kwargs["axes"]:
1298                    addstats = True
1299
1300        if addstats:
1301            htitle = f"Entries:~~{int(self.entries)}  "
1302            htitle += f"Mean:~~{utils.precision(self.mean, 3)}  "
1303            htitle += f"STD:~~{utils.precision(self.std, 3)}  "
1304            axes_opts["htitle"] = htitle
1305            axes_opts["htitle_justify"] = "bottom-left"
1306            axes_opts["htitle_size"] = 0.0175
1307            axes_opts["htitle_offset"] = [-0.49, 0.01, 0]
1308
1309        ############################################### Figure init
1310        Figure.__init__(self, xlim, ylim, aspect, padding, **fig_kwargs)
1311
1312        if self.yscale:
1313            ##################### the grid
1314            acts = []
1315            g = shapes.Grid(
1316                pos=[(xlim[0] + xlim[1]) / 2, (ylim[0] + ylim[1]) / 2, 0],
1317                s=(dx, dy),
1318                res=bins[:2],
1319            )
1320            g.alpha(alpha).lw(0).wireframe(False).flat().lighting("off")
1321            g.cmap(cmap, np.ravel(H.T), on="cells", vmin=zlim[0], vmax=zlim[1])
1322            if gap:
1323                g.shrink(abs(1 - gap))
1324
1325            if scalarbar:
1326                sc = g.add_scalarbar3d(ztitle, c=ac).scalarbar
1327                sc.scale([self.yscale, 1, 1])  ## prescale trick
1328                sbnds = sc.xbounds()
1329                sc.x(self.x1lim + (sbnds[1] - sbnds[0]) * 0.75)
1330                acts.append(sc)
1331            acts.append(g)
1332
1333            self.insert(*acts, as3d=False)
1334            self.name = "Histogram2D"
1335
1336
1337#########################################################################################
1338class PlotBars(Figure):
1339    """Creates a `PlotBars(Figure)` object."""
1340
1341    def __init__(
1342            self,
1343            data,
1344            errors=False,
1345            logscale=False,
1346            fill=True,
1347            gap=0.02,
1348            radius=0.05,
1349            c="olivedrab",
1350            alpha=1,
1351            texture="",
1352            outline=False,
1353            lw=2,
1354            lc="k",
1355            # Figure and axes options:
1356            like=None,
1357            xlim=(None, None),
1358            ylim=(0, None),
1359            aspect=4/3,
1360            padding=(0.025, 0.025, 0, 0.05),
1361            #
1362            title="",
1363            xtitle=" ",
1364            ytitle=" ",
1365            ac="k",
1366            grid=False,
1367            ztolerance=None,
1368            **fig_kwargs,
1369        ):
1370        """
1371        Input must be in format `[counts, labels, colors, edges]`.
1372        Either or both `edges` and `colors` are optional and can be omitted.
1373
1374        Use keyword `like=...` if you want to use the same format of a previously
1375        created Figure (useful when superimposing Figures) to make sure
1376        they are compatible and comparable. If they are not compatible
1377        you will receive an error message.
1378
1379        Arguments:
1380            errors : (bool)
1381                show error bars
1382            logscale : (bool)
1383                use logscale on y-axis
1384            fill : (bool)
1385                fill bars with solid color `c`
1386            gap : (float)
1387                leave a small space btw bars
1388            radius : (float)
1389                border radius of the top of the histogram bar. Default value is 0.1.
1390            texture : (str)
1391                url or path to an image to be used as texture for the bin
1392            outline : (bool)
1393                show outline of the bins
1394            xtitle : (str)
1395                title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
1396            ytitle : (str)
1397                title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
1398            ac : (str)
1399                axes color
1400            padding : (float, list)
1401                keep a padding space from the axes (as a fraction of the axis size).
1402                This can be a list of four numbers.
1403            aspect : (float)
1404                the desired aspect ratio of the figure. Default is 4/3.
1405            grid : (bool)
1406                show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
1407
1408        Examples:
1409            - [plot_bars.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/plot_bars.py)
1410
1411               ![](https://vedo.embl.es/images/pyplot/plot_bars.png)
1412        """
1413        ndata = len(data)
1414        if ndata == 4:
1415            counts, xlabs, cols, edges = data
1416        elif ndata == 3:
1417            counts, xlabs, cols = data
1418            edges = np.array(range(len(counts) + 1)) + 0.5
1419        elif ndata == 2:
1420            counts, xlabs = data
1421            edges = np.array(range(len(counts) + 1)) + 0.5
1422            cols = [c] * len(counts)
1423        else:
1424            m = "barplot error: data must be given as [counts, labels, colors, edges] not\n"
1425            vedo.logger.error(m + f" {data}\n     bin edges and colors are optional.")
1426            raise RuntimeError()
1427
1428        # sanity checks
1429        assert len(counts) == len(xlabs)
1430        assert len(counts) == len(cols)
1431        assert len(counts) == len(edges) - 1
1432
1433        counts = np.asarray(counts)
1434        edges = np.asarray(edges)
1435
1436        if logscale:
1437            counts = np.log10(counts + 1)
1438            if ytitle == " ":
1439                ytitle = "log_10 (counts+1)"
1440
1441        if like is None and vedo.last_figure is not None:
1442            if xlim == (None, None) and ylim == (0, None):
1443                like = vedo.last_figure
1444
1445        if like is not None:
1446            xlim = like.xlim
1447            ylim = like.ylim
1448            aspect = like.aspect
1449            padding = like.padding
1450
1451        if utils.is_sequence(xlim):
1452            # deal with user passing eg [x0, None]
1453            _x0, _x1 = xlim
1454            if _x0 is None:
1455                _x0 = np.min(edges)
1456            if _x1 is None:
1457                _x1 = np.max(edges)
1458            xlim = [_x0, _x1]
1459
1460        x0, x1 = np.min(edges), np.max(edges)
1461        y0, y1 = ylim[0], np.max(counts)
1462
1463        if like is None:
1464            ylim = list(ylim)
1465            if xlim is None:
1466                xlim = [x0, x1]
1467            if ylim[1] is None:
1468                ylim[1] = y1
1469            if ylim[0] != 0:
1470                ylim[0] = y0
1471
1472        fig_kwargs["title"] = title
1473        fig_kwargs["xtitle"] = xtitle
1474        fig_kwargs["ytitle"] = ytitle
1475        fig_kwargs["ac"] = ac
1476        fig_kwargs["ztolerance"] = ztolerance
1477        fig_kwargs["grid"] = grid
1478
1479        centers = (edges[0:-1] + edges[1:]) / 2
1480        binsizes = (centers - edges[0:-1]) * 2
1481
1482        if "axes" not in fig_kwargs:
1483            fig_kwargs["axes"] = {}
1484
1485        _xlabs = []
1486        for center, xlb in zip(centers, xlabs):
1487            _xlabs.append([center, str(xlb)])
1488        fig_kwargs["axes"]["x_values_and_labels"] = _xlabs
1489
1490        ############################################### Figure
1491        self.statslegend = ""
1492        self.edges = edges
1493        self.centers = centers
1494        self.bins = edges  # internal used by "like"
1495        Figure.__init__(self, xlim, ylim, aspect, padding, **fig_kwargs)
1496        if not self.yscale:
1497            return
1498
1499        rs = []
1500        maxheigth = 0
1501        if fill:  #####################
1502            if outline:
1503                gap = 0
1504
1505            for i in range(len(centers)):
1506                binsize = binsizes[i]
1507                p0 = (edges[i] + gap * binsize, 0, 0)
1508                p1 = (edges[i + 1] - gap * binsize, counts[i], 0)
1509
1510                if radius:
1511                    rds = np.array([0, 0, radius, radius])
1512                    p1_yscaled = [p1[0], p1[1]*self.yscale, 0]
1513                    r = shapes.Rectangle(p0, p1_yscaled, radius=rds*binsize, res=6)
1514                    r.scale([1, 1/self.yscale, 1])
1515                    r.radius = None  # so it doesnt get recreated and rescaled by insert()
1516                else:
1517                    r = shapes.Rectangle(p0, p1)
1518
1519                if texture:
1520                    r.texture(texture)
1521                    c = "w"
1522
1523                r.PickableOff()
1524                maxheigth = max(maxheigth, p1[1])
1525                if c in colors.cmaps_names:
1526                    col = colors.color_map((p0[0] + p1[0]) / 2, c, edges[0], edges[-1])
1527                else:
1528                    col = cols[i]
1529                r.color(col).alpha(alpha).lighting("off")
1530                r.name = f"bar_{i}"
1531                r.z(self.ztolerance)
1532                rs.append(r)
1533
1534        elif outline:  #####################
1535            lns = [[edges[0], 0, 0]]
1536            for i in range(len(centers)):
1537                lns.append([edges[i], counts[i], 0])
1538                lns.append([edges[i + 1], counts[i], 0])
1539                maxheigth = max(maxheigth, counts[i])
1540            lns.append([edges[-1], 0, 0])
1541            outl = shapes.Line(lns, c=lc, alpha=alpha, lw=lw).z(self.ztolerance)
1542            outl.name = f"bar_outline_{i}"
1543            rs.append(outl)
1544
1545        if errors:  #####################
1546            for x, f in centers:
1547                err = np.sqrt(f)
1548                el = shapes.Line([x, f-err/2, 0], [x, f+err/2, 0], c=lc, alpha=alpha, lw=lw)
1549                el.z(self.ztolerance * 2)
1550                rs.append(el)
1551
1552        self.insert(*rs, as3d=False)
1553        self.name = "PlotBars"
1554
1555
1556#########################################################################################
1557class PlotXY(Figure):
1558    """Creates a `PlotXY(Figure)` object."""
1559    def __init__(
1560            self,
1561            #
1562            data,
1563            xerrors=None,
1564            yerrors=None,
1565            #
1566            lw=2,
1567            lc=None,
1568            la=1,
1569            dashed=False,
1570            splined=False,
1571            #
1572            elw=2,  # error line width
1573            ec=None,  # error line or band color
1574            error_band=False, # errors in x are ignored
1575            #
1576            marker="",
1577            ms=None,
1578            mc=None,
1579            ma=None,
1580            # Figure and axes options:
1581            like=None,
1582            xlim=None,
1583            ylim=(None, None),
1584            aspect=4/3,
1585            padding=0.05,
1586            #
1587            title="",
1588            xtitle=" ",
1589            ytitle=" ",
1590            ac="k",
1591            grid=True,
1592            ztolerance=None,
1593            label="",
1594            **fig_kwargs,
1595        ):
1596        """
1597        Arguments:
1598            xerrors : (bool)
1599                show error bars associated to each point in x
1600            yerrors : (bool)
1601                show error bars associated to each point in y
1602            lw : (int)
1603                width of the line connecting points in pixel units.
1604                Set it to 0 to remove the line.
1605            lc : (str)
1606                line color
1607            la : (float)
1608                line "alpha", opacity of the line
1609            dashed : (bool)
1610                draw a dashed line instead of a continuous line
1611            splined : (bool)
1612                spline the line joining the point as a countinous curve
1613            elw : (int)
1614                width of error bar lines in units of pixels
1615            ec : (color)
1616                color of error bar, by default the same as marker color
1617            error_band : (bool)
1618                represent errors on y as a filled error band.
1619                Use `ec` keyword to modify its color.
1620            marker : (str, int)
1621                use a marker for the data points
1622            ms : (float)
1623                marker size
1624            mc : (color)
1625                color of the marker
1626            ma : (float)
1627                opacity of the marker
1628            xlim : (list)
1629                set limits to the range for the x variable
1630            ylim : (list)
1631                set limits to the range for the y variable
1632            aspect : (float, str)
1633                Desired aspect ratio.
1634                Use `aspect="equal"` to force the same units in x and y.
1635                Scaling factor is saved in Figure.yscale.
1636            padding : (float, list)
1637                keep a padding space from the axes (as a fraction of the axis size).
1638                This can be a list of four numbers.
1639            title : (str)
1640                title to appear on the top of the frame, like a header.
1641            xtitle : (str)
1642                title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
1643            ytitle : (str)
1644                title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
1645            ac : (str)
1646                axes color
1647            grid : (bool)
1648                show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
1649            ztolerance : (float)
1650                a tolerance factor to superimpose objects (along the z-axis).
1651
1652        Example:
1653            ```python
1654            import numpy as np
1655            from vedo.pyplot import plot
1656            x = np.arange(0, np.pi, 0.1)
1657            fig = plot(x, np.sin(2*x), 'r0-', aspect='equal')
1658            fig+= plot(x, np.cos(2*x), 'blue4 o-', like=fig)
1659            fig.show().close()
1660            ```
1661            ![](https://vedo.embl.es/images/feats/plotxy.png)
1662
1663        Examples:
1664            - [plot_errbars.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/plot_errbars.py)
1665            - [plot_errband.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/plot_errband.py)
1666            - [plot_pip.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/plot_pip.py)
1667
1668                ![](https://vedo.embl.es/images/pyplot/plot_pip.png)
1669
1670            - [scatter1.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/scatter1.py)
1671            - [scatter2.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/scatter2.py)
1672
1673                ![](https://vedo.embl.es/images/pyplot/scatter2.png)
1674        """
1675        line = False
1676        if lw > 0:
1677            line = True
1678        if marker == "" and not line and not splined:
1679            marker = "o"
1680
1681        if like is None and vedo.last_figure is not None:
1682            if xlim is None and ylim == (None, None):
1683                like = vedo.last_figure
1684
1685        if like is not None:
1686            xlim = like.xlim
1687            ylim = like.ylim
1688            aspect = like.aspect
1689            padding = like.padding
1690
1691        if utils.is_sequence(xlim):
1692            # deal with user passing eg [x0, None]
1693            _x0, _x1 = xlim
1694            if _x0 is None:
1695                _x0 = data.min()
1696            if _x1 is None:
1697                _x1 = data.max()
1698            xlim = [_x0, _x1]
1699
1700        # purge NaN from data
1701        validIds = np.all(np.logical_not(np.isnan(data)))
1702        data = np.array(data[validIds])[0]
1703
1704        fig_kwargs["title"] = title
1705        fig_kwargs["xtitle"] = xtitle
1706        fig_kwargs["ytitle"] = ytitle
1707        fig_kwargs["ac"] = ac
1708        fig_kwargs["ztolerance"] = ztolerance
1709        fig_kwargs["grid"] = grid
1710
1711        x0, y0 = np.min(data, axis=0)
1712        x1, y1 = np.max(data, axis=0)
1713        if xerrors is not None and not error_band:
1714            x0 = min(data[:, 0] - xerrors)
1715            x1 = max(data[:, 0] + xerrors)
1716        if yerrors is not None:
1717            y0 = min(data[:, 1] - yerrors)
1718            y1 = max(data[:, 1] + yerrors)
1719
1720        if like is None:
1721            if xlim is None:
1722                xlim = (None, None)
1723            xlim = list(xlim)
1724            if xlim[0] is None:
1725                xlim[0] = x0
1726            if xlim[1] is None:
1727                xlim[1] = x1
1728            ylim = list(ylim)
1729            if ylim[0] is None:
1730                ylim[0] = y0
1731            if ylim[1] is None:
1732                ylim[1] = y1
1733
1734        self.entries = len(data)
1735        self.mean = data.mean()
1736        self.std = data.std()
1737
1738        ######### the PlotXY marker
1739        # fall back solutions logic for colors
1740        if "c" in fig_kwargs:
1741            if mc is None:
1742                mc = fig_kwargs["c"]
1743            if lc is None:
1744                lc = fig_kwargs["c"]
1745            if ec is None:
1746                ec = fig_kwargs["c"]
1747        if lc is None:
1748            lc = 'k'
1749        if mc is None:
1750            mc = lc
1751        if ma is None:
1752            ma = la
1753        if ec is None:
1754            if mc is None:
1755                ec = lc
1756            else:
1757                ec = mc
1758
1759        if label:
1760            nlab = LabelData()
1761            nlab.text = label
1762            nlab.tcolor = ac
1763            nlab.marker = marker
1764            if line and marker == "":
1765                nlab.marker = "-"
1766            nlab.mcolor = mc
1767            fig_kwargs["label"] = nlab
1768
1769        ############################################### Figure init
1770        Figure.__init__(self, xlim, ylim, aspect, padding, **fig_kwargs)
1771
1772        if not self.yscale:
1773            return
1774
1775        acts = []
1776
1777        ######### the PlotXY Line or Spline
1778        if dashed:
1779            l = shapes.DashedLine(data, c=lc, alpha=la, lw=lw)
1780            acts.append(l)
1781        elif splined:
1782            l = shapes.KSpline(data).lw(lw).c(lc).alpha(la)
1783            acts.append(l)
1784        elif line:
1785            l = shapes.Line(data, c=lc, alpha=la).lw(lw)
1786            acts.append(l)
1787
1788        if marker:
1789
1790            pts = np.c_[data, np.zeros(len(data))]
1791
1792            if utils.is_sequence(ms):
1793                ### variable point size
1794                mk = shapes.Marker(marker, s=1)
1795                mk.scale([1, 1 / self.yscale, 1])
1796                msv = np.zeros_like(pts)
1797                msv[:, 0] = ms
1798                marked = shapes.Glyph(
1799                    pts, mk, c=mc, orientation_array=msv, scale_by_vector_size=True
1800                )
1801            else:
1802                ### fixed point size
1803                if ms is None:
1804                    ms = (xlim[1] - xlim[0]) / 100.0
1805
1806                if utils.is_sequence(mc):
1807                    fig_kwargs["marker_color"] = None  # for labels
1808                    mk = shapes.Marker(marker, s=ms)
1809                    mk.scale([1, 1 / self.yscale, 1])
1810                    msv = np.zeros_like(pts)
1811                    msv[:, 0] = 1
1812                    marked = shapes.Glyph(
1813                        pts,
1814                        mk,
1815                        c=mc,
1816                        orientation_array=msv,
1817                        scale_by_vector_size=True,
1818                    )
1819                else:
1820                    mk = shapes.Marker(marker, s=ms)
1821                    mk.scale([1, 1 / self.yscale, 1])
1822                    marked = shapes.Glyph(pts, mk, c=mc)
1823
1824            marked.name = "Marker"
1825            marked.alpha(ma)
1826            marked.z(3 * self.ztolerance)
1827            acts.append(marked)
1828
1829        ######### the PlotXY marker errors
1830        ztol = self.ztolerance
1831
1832        if error_band:
1833            yerrors = np.abs(yerrors)
1834            du = np.array(data)
1835            dd = np.array(data)
1836            du[:, 1] += yerrors
1837            dd[:, 1] -= yerrors
1838            if splined:
1839                res = len(data) * 20
1840                band1 = shapes.KSpline(du, res=res)
1841                band2 = shapes.KSpline(dd, res=res)
1842                band = shapes.Ribbon(band1, band2, res=(res, 2))
1843            else:
1844                dd = list(reversed(dd.tolist()))
1845                band = shapes.Line(du.tolist() + dd, closed=True)
1846                band.triangulate().lw(0)
1847            if ec is None:
1848                band.c(lc)
1849            else:
1850                band.c(ec)
1851            band.lighting("off").alpha(la).z(ztol/10)
1852            acts.append(band)
1853
1854        else:
1855
1856            ## xerrors
1857            if xerrors is not None:
1858                if len(xerrors) == len(data):
1859                    errs = []
1860                    for i, val in enumerate(data):
1861                        xval, yval = val
1862                        xerr = xerrors[i] / 2
1863                        el = shapes.Line((xval-xerr, yval, ztol), (xval+xerr, yval, ztol))
1864                        el.lw(elw)
1865                        errs.append(el)
1866                    mxerrs = merge(errs).c(ec).lw(lw).alpha(ma).z(2 * ztol)
1867                    acts.append(mxerrs)
1868                else:
1869                    vedo.logger.error("in PlotXY(xerrors=...): mismatch in array length")
1870
1871            ## yerrors
1872            if yerrors is not None:
1873                if len(yerrors) == len(data):
1874                    errs = []
1875                    for i, val in enumerate(data):
1876                        xval, yval = val
1877                        yerr = yerrors[i]
1878                        el = shapes.Line((xval, yval-yerr, ztol), (xval, yval+yerr, ztol))
1879                        el.lw(elw)
1880                        errs.append(el)
1881                    myerrs = merge(errs).c(ec).lw(lw).alpha(ma).z(2 * ztol)
1882                    acts.append(myerrs)
1883                else:
1884                    vedo.logger.error("in PlotXY(yerrors=...): mismatch in array length")
1885
1886        self.insert(*acts, as3d=False)
1887        self.name = "PlotXY"
1888
1889
1890def plot(*args, **kwargs):
1891    """
1892    Draw a 2D line plot, or scatter plot, of variable x vs variable y.
1893    Input format can be either `[allx], [allx, ally] or [(x1,y1), (x2,y2), ...]`
1894
1895    Use `like=...` if you want to use the same format of a previously
1896    created Figure (useful when superimposing Figures) to make sure
1897    they are compatible and comparable. If they are not compatible
1898    you will receive an error message.
1899
1900    Arguments:
1901        xerrors : (bool)
1902            show error bars associated to each point in x
1903        yerrors : (bool)
1904            show error bars associated to each point in y
1905        lw : (int)
1906            width of the line connecting points in pixel units.
1907            Set it to 0 to remove the line.
1908        lc : (str)
1909            line color
1910        la : (float)
1911            line "alpha", opacity of the line
1912        dashed : (bool)
1913            draw a dashed line instead of a continuous line
1914        splined : (bool)
1915            spline the line joining the point as a countinous curve
1916        elw : (int)
1917            width of error bar lines in units of pixels
1918        ec : (color)
1919            color of error bar, by default the same as marker color
1920        error_band : (bool)
1921            represent errors on y as a filled error band.
1922            Use `ec` keyword to modify its color.
1923        marker : (str, int)
1924            use a marker for the data points
1925        ms : (float)
1926            marker size
1927        mc : (color)
1928            color of the marker
1929        ma : (float)
1930            opacity of the marker
1931        xlim : (list)
1932            set limits to the range for the x variable
1933        ylim : (list)
1934            set limits to the range for the y variable
1935        aspect : (float)
1936            Desired aspect ratio.
1937            If None, it is automatically calculated to get a reasonable aspect ratio.
1938            Scaling factor is saved in Figure.yscale
1939        padding : (float, list)
1940            keep a padding space from the axes (as a fraction of the axis size).
1941            This can be a list of four numbers.
1942        title : (str)
1943            title to appear on the top of the frame, like a header.
1944        xtitle : (str)
1945            title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
1946        ytitle : (str)
1947            title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
1948        ac : (str)
1949            axes color
1950        grid : (bool)
1951            show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
1952        ztolerance : (float)
1953            a tolerance factor to superimpose objects (along the z-axis).
1954
1955    Example:
1956        ```python
1957        from vedo.pyplot import plot
1958        import numpy as np
1959        x = np.linspace(0, 6.28, num=50)
1960        plot(np.sin(x), 'r').plot(np.cos(x), 'bo-').show().close()
1961        ```
1962        <img src="https://user-images.githubusercontent.com/32848391/74363882-c3638300-4dcb-11ea-8a78-eb492ad9711f.png" width="600">
1963
1964    Examples:
1965        - [plot_errbars.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_errbars.py)
1966        - [plot_errband.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_errband.py)
1967        - [plot_pip.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_pip.py)
1968        
1969            ![](https://vedo.embl.es/images/pyplot/plot_pip.png)
1970        
1971        - [scatter1.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/scatter1.py)
1972        - [scatter2.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/scatter2.py)
1973
1974
1975
1976    -------------------------------------------------------------------------
1977    .. note:: mode="bar"
1978
1979    Creates a `PlotBars(Figure)` object.
1980
1981    Input must be in format `[counts, labels, colors, edges]`.
1982    Either or both `edges` and `colors` are optional and can be omitted.
1983
1984    Arguments:
1985        errors : (bool)
1986            show error bars
1987        logscale : (bool)
1988            use logscale on y-axis
1989        fill : (bool)
1990            fill bars with solid color `c`
1991        gap : (float)
1992            leave a small space btw bars
1993        radius : (float)
1994            border radius of the top of the histogram bar. Default value is 0.1.
1995        texture : (str)
1996            url or path to an image to be used as texture for the bin
1997        outline : (bool)
1998            show outline of the bins
1999        xtitle : (str)
2000            title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
2001        ytitle : (str)
2002            title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
2003        ac : (str)
2004            axes color
2005        padding : (float, list)
2006            keep a padding space from the axes (as a fraction of the axis size).
2007            This can be a list of four numbers.
2008        aspect : (float)
2009            the desired aspect ratio of the figure. Default is 4/3.
2010        grid : (bool)
2011            show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
2012
2013    Examples:
2014        - [histo_1d_a.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_a.py)
2015        - [histo_1d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_b.py)
2016        - [histo_1d_c.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_c.py)
2017        - [histo_1d_d.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_d.py)
2018
2019        ![](https://vedo.embl.es/images/pyplot/histo_1D.png)
2020
2021
2022    ----------------------------------------------------------------------
2023    .. note:: 2D functions
2024
2025    If input is an external function or a formula, draw the surface
2026    representing the function `f(x,y)`.
2027
2028    Arguments:
2029        x : (float)
2030            x range of values
2031        y : (float)
2032            y range of values
2033        zlimits : (float)
2034            limit the z range of the independent variable
2035        zlevels : (int)
2036            will draw the specified number of z-levels contour lines
2037        show_nan : (bool)
2038            show where the function does not exist as red points
2039        bins : (list)
2040            number of bins in x and y
2041
2042    Examples:
2043        - [plot_fxy.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_fxy.py)
2044
2045            ![](https://vedo.embl.es/images/pyplot/plot_fxy.png)
2046
2047
2048    --------------------------------------------------------------------
2049    .. note:: mode="complex"
2050
2051    If `mode='complex'` draw the real value of the function and color map the imaginary part.
2052
2053    Arguments:
2054        cmap : (str)
2055            diverging color map (white means `imag(z)=0`)
2056        lw : (float)
2057            line with of the binning
2058        bins : (list)
2059            binning in x and y
2060
2061    Examples:
2062        - [plot_fxy.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_fxy.py)
2063
2064            ![](https://user-images.githubusercontent.com/32848391/73392962-1709a300-42db-11ea-9278-30c9d6e5eeaa.png)
2065
2066
2067    --------------------------------------------------------------------
2068    .. note:: mode="polar"
2069
2070    If `mode='polar'` input arrays are interpreted as a list of polar angles and radii.
2071    Build a polar (radar) plot by joining the set of points in polar coordinates.
2072
2073    Arguments:
2074        title : (str)
2075            plot title
2076        tsize : (float)
2077            title size
2078        bins : (int)
2079            number of bins in phi
2080        r1 : (float)
2081            inner radius
2082        r2 : (float)
2083            outer radius
2084        lsize : (float)
2085            label size
2086        c : (color)
2087            color of the line
2088        ac : (color)
2089            color of the frame and labels
2090        alpha : (float)
2091            opacity of the frame
2092        ps : (int)
2093            point size in pixels, if ps=0 no point is drawn
2094        lw : (int)
2095            line width in pixels, if lw=0 no line is drawn
2096        deg : (bool)
2097            input array is in degrees
2098        vmax : (float)
2099            normalize radius to this maximum value
2100        fill : (bool)
2101            fill convex area with solid color
2102        splined : (bool)
2103            interpolate the set of input points
2104        show_disc : (bool)
2105            draw the outer ring axis
2106        nrays : (int)
2107            draw this number of axis rays (continuous and dashed)
2108        show_lines : (bool)
2109            draw lines to the origin
2110        show_angles : (bool)
2111            draw angle values
2112
2113    Examples:
2114        - [histo_polar.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_polar.py)
2115
2116            ![](https://user-images.githubusercontent.com/32848391/64992590-7fc82400-d8d4-11e9-9c10-795f4756a73f.png)
2117
2118
2119    --------------------------------------------------------------------
2120    .. note:: mode="spheric"
2121
2122    If `mode='spheric'` input must be an external function rho(theta, phi).
2123    A surface is created in spherical coordinates.
2124
2125    Return an `Figure(Assembly)` of 2 objects: the unit
2126    sphere (in wireframe representation) and the surface `rho(theta, phi)`.
2127
2128    Arguments:
2129        rfunc : function
2130            handle to a user defined function `rho(theta, phi)`.
2131        normalize : (bool)
2132            scale surface to fit inside the unit sphere
2133        res : (int)
2134            grid resolution of the unit sphere
2135        scalarbar : (bool)
2136            add a 3D scalarbar to the plot for radius
2137        c : (color)
2138            color of the unit sphere
2139        alpha : (float)
2140            opacity of the unit sphere
2141        cmap : (str)
2142            color map for the surface
2143
2144    Examples:
2145        - [plot_spheric.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_spheric.py)
2146
2147            ![](https://vedo.embl.es/images/pyplot/plot_spheric.png)
2148    """
2149    mode = kwargs.pop("mode", "")
2150    if "spher" in mode:
2151        return _plot_spheric(args[0], **kwargs)
2152
2153    if "bar" in mode:
2154        return PlotBars(args[0], **kwargs)
2155
2156    if isinstance(args[0], str) or "function" in str(type(args[0])):
2157        if "complex" in mode:
2158            return _plot_fz(args[0], **kwargs)
2159        return _plot_fxy(args[0], **kwargs)
2160
2161    # grab the matplotlib-like options
2162    optidx = None
2163    for i, a in enumerate(args):
2164        if i > 0 and isinstance(a, str):
2165            optidx = i
2166            break
2167    if optidx:
2168        opts = args[optidx].replace(" ", "")
2169        if "--" in opts:
2170            opts = opts.replace("--", "")
2171            kwargs["dashed"] = True
2172        elif "-" in opts:
2173            opts = opts.replace("-", "")
2174        else:
2175            kwargs["lw"] = 0
2176
2177        symbs = ['.','o','O', '0', 'p','*','h','D','d','v','^','>','<','s', 'x', 'a']
2178
2179        allcols = list(colors.colors.keys()) + list(colors.color_nicks.keys())
2180        for cc in allcols:
2181            if cc == "o":
2182                continue
2183            if cc in opts:
2184                opts = opts.replace(cc, "")
2185                kwargs["lc"] = cc
2186                kwargs["mc"] = cc
2187                break
2188
2189        for ss in symbs:
2190            if ss in opts:
2191                opts = opts.replace(ss, "", 1)
2192                kwargs["marker"] = ss
2193                break
2194
2195        opts.replace(" ", "")
2196        if opts:
2197            vedo.logger.error(f"in plot(), could not understand option(s): {opts}")
2198
2199    if optidx == 1 or optidx is None:
2200        if utils.is_sequence(args[0][0]) and len(args[0][0])>1:
2201            # print('------------- case 1', 'plot([(x,y),..])')
2202            data = np.asarray(args[0]) #(x,y)
2203            x = np.asarray(data[:, 0])
2204            y = np.asarray(data[:, 1])
2205
2206        elif len(args) == 1 or optidx == 1:
2207            # print('------------- case 2', 'plot(x)')
2208            if "pandas" in str(type(args[0])):
2209                if "ytitle" not in kwargs:
2210                    kwargs.update({"ytitle" : args[0].name.replace("_","\_")})
2211            x = np.linspace(0, len(args[0]), num=len(args[0]))
2212            y = np.asarray(args[0]).ravel()
2213
2214        elif utils.is_sequence(args[1]):
2215            # print('------------- case 3', 'plot(allx,ally)',str(type(args[0])))
2216            if "pandas" in str(type(args[0])):
2217                if "xtitle" not in kwargs:
2218                    kwargs.update({"xtitle" : args[0].name.replace("_","\_")})
2219            if "pandas" in str(type(args[1])):
2220                if "ytitle" not in kwargs:
2221                    kwargs.update({"ytitle" : args[1].name.replace("_","\_")})
2222            x = np.asarray(args[0]).ravel()
2223            y = np.asarray(args[1]).ravel()
2224
2225        elif utils.is_sequence(args[0]) and utils.is_sequence(args[0][0]):
2226            # print('------------- case 4', 'plot([allx,ally])')
2227            x = np.asarray(args[0][0]).ravel()
2228            y = np.asarray(args[0][1]).ravel()
2229
2230    elif optidx == 2:
2231        # print('------------- case 5', 'plot(x,y)')
2232        x = np.asarray(args[0]).ravel()
2233        y = np.asarray(args[1]).ravel()
2234
2235    else:
2236        vedo.logger.error(f"plot(): Could not understand input arguments {args}")
2237        return None
2238
2239    if "polar" in mode:
2240        return _plot_polar(np.c_[x, y], **kwargs)
2241
2242    return PlotXY(np.c_[x, y], **kwargs)
2243
2244
2245def histogram(*args, **kwargs):
2246    """
2247    Histogramming for 1D and 2D data arrays.
2248
2249    This is meant as a convenience function that creates the appropriate object
2250    based on the shape of the provided input data.
2251
2252    Use keyword `like=...` if you want to use the same format of a previously
2253    created Figure (useful when superimposing Figures) to make sure
2254    they are compatible and comparable. If they are not compatible
2255    you will receive an error message.
2256
2257    -------------------------------------------------------------------------
2258    .. note:: default mode, for 1D arrays
2259
2260    Creates a `Histogram1D(Figure)` object.
2261
2262    Arguments:
2263        weights : (list)
2264            An array of weights, of the same shape as `data`. Each value in `data`
2265            only contributes its associated weight towards the bin count (instead of 1).
2266        bins : (int)
2267            number of bins
2268        vrange : (list)
2269            restrict the range of the histogram
2270        density : (bool)
2271            normalize the area to 1 by dividing by the nr of entries and bin size
2272        logscale : (bool)
2273            use logscale on y-axis
2274        fill : (bool)
2275            fill bars with solid color `c`
2276        gap : (float)
2277            leave a small space btw bars
2278        radius : (float)
2279            border radius of the top of the histogram bar. Default value is 0.1.
2280        texture : (str)
2281            url or path to an image to be used as texture for the bin
2282        outline : (bool)
2283            show outline of the bins
2284        errors : (bool)
2285            show error bars
2286        xtitle : (str)
2287            title for the x-axis, can also be set using `axes=dict(xtitle="my x axis")`
2288        ytitle : (str)
2289            title for the y-axis, can also be set using `axes=dict(ytitle="my y axis")`
2290        padding : (float, list)
2291            keep a padding space from the axes (as a fraction of the axis size).
2292            This can be a list of four numbers.
2293        aspect : (float)
2294            the desired aspect ratio of the histogram. Default is 4/3.
2295        grid : (bool)
2296            show the background grid for the axes, can also be set using `axes=dict(xygrid=True)`
2297        ztolerance : (float)
2298            a tolerance factor to superimpose objects (along the z-axis).
2299
2300    Examples:
2301        - [histo_1d_a.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_a.py)
2302        - [histo_1d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_b.py)
2303        - [histo_1d_c.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_c.py)
2304        - [histo_1d_d.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_1d_d.py)
2305
2306        ![](https://vedo.embl.es/images/pyplot/histo_1D.png)
2307
2308
2309    -------------------------------------------------------------------------
2310    .. note:: default mode, for 2D arrays
2311
2312    Input data formats `[(x1,x2,..), (y1,y2,..)] or [(x1,y1), (x2,y2),..]`
2313    are both valid.
2314
2315    Arguments:
2316        bins : (list)
2317            binning as (nx, ny)
2318        weights : (list)
2319            array of weights to assign to each entry
2320        cmap : (str, lookuptable)
2321            color map name or look up table
2322        alpha : (float)
2323            opacity of the histogram
2324        gap : (float)
2325            separation between adjacent bins as a fraction for their size.
2326            Set gap=-1 to generate a quad surface.
2327        scalarbar : (bool)
2328            add a scalarbar to right of the histogram
2329        like : (Figure)
2330            grab and use the same format of the given Figure (for superimposing)
2331        xlim : (list)
2332            [x0, x1] range of interest. If left to None will automatically
2333            choose the minimum or the maximum of the data range.
2334            Data outside the range are completely ignored.
2335        ylim : (list)
2336            [y0, y1] range of interest. If left to None will automatically
2337            choose the minimum or the maximum of the data range.
2338            Data outside the range are completely ignored.
2339        aspect : (float)
2340            the desired aspect ratio of the figure.
2341        title : (str)
2342            title of the plot to appear on top.
2343            If left blank some statistics will be shown.
2344        xtitle : (str)
2345            x axis title
2346        ytitle : (str)
2347            y axis title
2348        ztitle : (str)
2349            title for the scalar bar
2350        ac : (str)
2351            axes color, additional keyword for Axes can also be added
2352            using e.g. `axes=dict(xygrid=True)`
2353
2354    Examples:
2355        - [histo_2d_a.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_2d_a.py)
2356        - [histo_2d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_2d_b.py)
2357
2358        ![](https://vedo.embl.es/images/pyplot/histo_2D.png)
2359
2360
2361    -------------------------------------------------------------------------
2362    .. note:: mode="3d"
2363
2364    If `mode='3d'`, build a 2D histogram as 3D bars from a list of x and y values.
2365
2366    Arguments:
2367        xtitle : (str)
2368            x axis title
2369        bins : (int)
2370            nr of bins for the smaller range in x or y
2371        vrange : (list)
2372            range in x and y in format `[(xmin,xmax), (ymin,ymax)]`
2373        norm : (float)
2374            sets a scaling factor for the z axis (frequency axis)
2375        fill : (bool)
2376            draw solid hexagons
2377        cmap : (str)
2378            color map name for elevation
2379        gap : (float)
2380            keep a internal empty gap between bins [0,1]
2381        zscale : (float)
2382            rescale the (already normalized) zaxis for visual convenience
2383
2384    Examples:
2385        - [histo_2d_b.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_2d_b.py)
2386
2387
2388    -------------------------------------------------------------------------
2389    .. note:: mode="hexbin"
2390
2391    If `mode='hexbin'`, build a hexagonal histogram from a list of x and y values.
2392
2393    Arguments:
2394        xtitle : (str)
2395            x axis title
2396        bins : (int)
2397            nr of bins for the smaller range in x or y
2398        vrange : (list)
2399            range in x and y in format `[(xmin,xmax), (ymin,ymax)]`
2400        norm : (float)
2401            sets a scaling factor for the z axis (frequency axis)
2402        fill : (bool)
2403            draw solid hexagons
2404        cmap : (str)
2405            color map name for elevation
2406
2407    Examples:
2408        - [histo_hexagonal.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_hexagonal.py)
2409
2410        ![](https://vedo.embl.es/images/pyplot/histo_hexagonal.png)
2411
2412
2413    -------------------------------------------------------------------------
2414    .. note:: mode="polar"
2415
2416    If `mode='polar'` assume input is polar coordinate system (rho, theta):
2417
2418    Arguments:
2419        weights : (list)
2420            Array of weights, of the same shape as the input.
2421            Each value only contributes its associated weight towards the bin count (instead of 1).
2422        title : (str)
2423            histogram title
2424        tsize : (float)
2425            title size
2426        bins : (int)
2427            number of bins in phi
2428        r1 : (float)
2429            inner radius
2430        r2 : (float)
2431            outer radius
2432        phigap : (float)
2433            gap angle btw 2 radial bars, in degrees
2434        rgap : (float)
2435            gap factor along radius of numeric angle labels
2436        lpos : (float)
2437            label gap factor along radius
2438        lsize : (float)
2439            label size
2440        c : (color)
2441            color of the histogram bars, can be a list of length `bins`
2442        bc : (color)
2443            color of the frame and labels
2444        alpha : (float)
2445            opacity of the frame
2446        cmap : (str)
2447            color map name
2448        deg : (bool)
2449            input array is in degrees
2450        vmin : (float)
2451            minimum value of the radial axis
2452        vmax : (float)
2453            maximum value of the radial axis
2454        labels : (list)
2455            list of labels, must be of length `bins`
2456        show_disc : (bool)
2457            show the outer ring axis
2458        nrays : (int)
2459            draw this number of axis rays (continuous and dashed)
2460        show_lines : (bool)
2461            show lines to the origin
2462        show_angles : (bool)
2463            show angular values
2464        show_errors : (bool)
2465            show error bars
2466
2467    Examples:
2468        - [histo_polar.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_polar.py)
2469
2470        ![](https://vedo.embl.es/images/pyplot/histo_polar.png)
2471
2472
2473    -------------------------------------------------------------------------
2474    .. note:: mode="spheric"
2475
2476    If `mode='spheric'`, build a histogram from list of theta and phi values.
2477
2478    Arguments:
2479        rmax : (float)
2480            maximum radial elevation of bin
2481        res : (int)
2482            sphere resolution
2483        cmap : (str)
2484            color map name
2485        lw : (int)
2486            line width of the bin edges
2487
2488    Examples:
2489        - [histo_spheric.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_spheric.py)
2490
2491        ![](https://vedo.embl.es/images/pyplot/histo_spheric.png)
2492    """
2493    mode = kwargs.pop("mode", "")
2494    if len(args) == 2:  # x, y
2495
2496        if "spher" in mode:
2497            return _histogram_spheric(args[0], args[1], **kwargs)
2498
2499        if "hex" in mode:
2500            return _histogram_hex_bin(args[0], args[1], **kwargs)
2501
2502        if "3d" in mode.lower():
2503            return _histogram_quad_bin(args[0], args[1], **kwargs)
2504
2505        return Histogram2D(args[0], args[1], **kwargs)
2506
2507    elif len(args) == 1:
2508
2509        if isinstance(args[0], vedo.Volume):
2510            data = args[0].pointdata[0]
2511        elif isinstance(args[0], vedo.Points):
2512            pd0 = args[0].pointdata[0]
2513            if pd0:
2514                data = pd0.ravel()
2515            else:
2516                data = args[0].celldata[0].ravel()
2517        else:
2518            try:
2519                if "pandas" in str(type(args[0])):
2520                    if "xtitle" not in kwargs:
2521                        kwargs.update({"xtitle" : args[0].name.replace("_","\_")})
2522            except:
2523                pass
2524            data = np.asarray(args[0])
2525
2526        if "spher" in mode:
2527            return _histogram_spheric(args[0][:, 0], args[0][:, 1], **kwargs)
2528
2529        if data.ndim == 1:
2530            if "polar" in mode:
2531                return _histogram_polar(data, **kwargs)
2532            return Histogram1D(data, **kwargs)
2533
2534        if "hex" in mode:
2535            return _histogram_hex_bin(args[0][:, 0], args[0][:, 1], **kwargs)
2536
2537        if "3d" in mode.lower():
2538            return _histogram_quad_bin(args[0][:, 0], args[0][:, 1], **kwargs)
2539
2540        return Histogram2D(args[0], **kwargs)
2541
2542    vedo.logger.error(f"in histogram(): could not understand input {args[0]}")
2543    return None
2544
2545
2546def fit(
2547        points,
2548        deg=1,
2549        niter=0,
2550        nstd=3,
2551        xerrors=None,
2552        yerrors=None,
2553        vrange=None,
2554        res=250,
2555        lw=3,
2556        c='red4',
2557    ):
2558    """
2559    Polynomial fitting with parameter error and error bands calculation.
2560    Errors bars in both x and y are supported.
2561
2562    Returns a `vedo.shapes.Line` object.
2563
2564    Additional information about the fitting output can be accessed with:
2565
2566    `fitd = fit(pts)`
2567
2568    - `fitd.coefficients` will contain the coefficients of the polynomial fit
2569    - `fitd.coefficient_errors`, errors on the fitting coefficients
2570    - `fitd.monte_carlo_coefficients`, fitting coefficient set from MC generation
2571    - `fitd.covariance_matrix`, covariance matrix as a numpy array
2572    - `fitd.reduced_chi2`, reduced chi-square of the fitting
2573    - `fitd.ndof`, number of degrees of freedom
2574    - `fitd.data_sigma`, mean data dispersion from the central fit assuming `Chi2=1`
2575    - `fitd.error_lines`, a `vedo.shapes.Line` object for the upper and lower error band
2576    - `fitd.error_band`, the `vedo.mesh.Mesh` object representing the error band
2577
2578    Errors on x and y can be specified. If left to `None` an estimate is made from
2579    the statistical spread of the dataset itself. Errors are always assumed gaussian.
2580
2581    Arguments:
2582        deg : (int)
2583            degree of the polynomial to be fitted
2584        niter : (int)
2585            number of monte-carlo iterations to compute error bands.
2586            If set to 0, return the simple least-squares fit with naive error estimation
2587            on coefficients only. A reasonable non-zero value to set is about 500, in
2588            this case *error_lines*, *error_band* and the other class attributes are filled
2589        nstd : (float)
2590            nr. of standard deviation to use for error calculation
2591        xerrors : (list)
2592            array of the same length of points with the errors on x
2593        yerrors : (list)
2594            array of the same length of points with the errors on y
2595        vrange : (list)
2596            specify the domain range of the fitting line
2597            (only affects visualization, but can be used to extrapolate the fit
2598            outside the data range)
2599        res : (int)
2600            resolution of the output fitted line and error lines
2601
2602    Examples:
2603        - [fit_polynomial1.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/fit_polynomial1.py)
2604
2605        ![](https://vedo.embl.es/images/pyplot/fitPolynomial1.png)
2606    """
2607    if isinstance(points, vedo.pointcloud.Points):
2608        points = points.points()
2609    points = np.asarray(points)
2610    if len(points) == 2:  # assume user is passing [x,y]
2611        points = np.c_[points[0], points[1]]
2612    x = points[:, 0]
2613    y = points[:, 1]  # ignore z
2614
2615    n = len(x)
2616    ndof = n - deg - 1
2617    if vrange is not None:
2618        x0, x1 = vrange
2619    else:
2620        x0, x1 = np.min(x), np.max(x)
2621        if xerrors is not None:
2622            x0 -= xerrors[0] / 2
2623            x1 += xerrors[-1] / 2
2624
2625    tol = (x1 - x0) / 10000
2626    xr = np.linspace(x0, x1, res)
2627
2628    # project x errs on y
2629    if xerrors is not None:
2630        xerrors = np.asarray(xerrors)
2631        if yerrors is not None:
2632            yerrors = np.asarray(yerrors)
2633            w = 1.0 / yerrors
2634            coeffs = np.polyfit(x, y, deg, w=w, rcond=None)
2635        else:
2636            coeffs = np.polyfit(x, y, deg, rcond=None)
2637        # update yerrors, 1 bootstrap iteration is enough
2638        p1d = np.poly1d(coeffs)
2639        der = (p1d(x + tol) - p1d(x)) / tol
2640        yerrors = np.sqrt(yerrors * yerrors + np.power(der * xerrors, 2))
2641
2642    if yerrors is not None:
2643        yerrors = np.asarray(yerrors)
2644        w = 1.0 / yerrors
2645        coeffs, V = np.polyfit(x, y, deg, w=w, rcond=None, cov=True)
2646    else:
2647        w = 1
2648        coeffs, V = np.polyfit(x, y, deg, rcond=None, cov=True)
2649
2650    p1d = np.poly1d(coeffs)
2651    theor = p1d(xr)
2652    fitl = shapes.Line(xr, theor, lw=lw, c=c).z(tol * 2)
2653    fitl.coefficients = coeffs
2654    fitl.covariance_matrix = V
2655    residuals2_sum = np.sum(np.power(p1d(x) - y, 2)) / ndof
2656    sigma = np.sqrt(residuals2_sum)
2657    fitl.reduced_chi2 = np.sum(np.power((p1d(x) - y) * w, 2)) / ndof
2658    fitl.ndof = ndof
2659    fitl.data_sigma = sigma  # worked out from data using chi2=1 hypo
2660    fitl.name = "LinearPolynomialFit"
2661
2662    if not niter:
2663        fitl.coefficient_errors = np.sqrt(np.diag(V))
2664        return fitl  ################################
2665
2666    if yerrors is not None:
2667        sigma = yerrors
2668    else:
2669        w = None
2670        fitl.reduced_chi2 = 1
2671
2672    Theors, all_coeffs = [], []
2673    for i in range(niter):
2674        noise = np.random.randn(n) * sigma
2675        coeffs = np.polyfit(x, y + noise, deg, w=w, rcond=None)
2676        all_coeffs.append(coeffs)
2677        P1d = np.poly1d(coeffs)
2678        Theor = P1d(xr)
2679        Theors.append(Theor)
2680    all_coeffs = np.array(all_coeffs)
2681    fitl.monte_carlo_coefficients = all_coeffs
2682
2683    stds = np.std(Theors, axis=0)
2684    fitl.coefficient_errors = np.std(all_coeffs, axis=0)
2685
2686    # check distributions on the fly
2687    # for i in range(deg+1):
2688    #     histogram(all_coeffs[:,i],title='par'+str(i)).show(new=1)
2689    # histogram(all_coeffs[:,0], all_coeffs[:,1],
2690    #           xtitle='param0', ytitle='param1',scalarbar=1).show(new=1)
2691    # histogram(all_coeffs[:,1], all_coeffs[:,2],
2692    #           xtitle='param1', ytitle='param2').show(new=1)
2693    # histogram(all_coeffs[:,0], all_coeffs[:,2],
2694    #           xtitle='param0', ytitle='param2').show(new=1)
2695
2696    error_lines = []
2697    for i in [nstd, -nstd]:
2698        el = shapes.Line(xr, theor + stds * i, lw=1, alpha=0.2, c="k").z(tol)
2699        error_lines.append(el)
2700        el.name = "ErrorLine for sigma=" + str(i)
2701
2702    fitl.error_lines = error_lines
2703    l1 = error_lines[0].points().tolist()
2704    cband = l1 + list(reversed(error_lines[1].points().tolist())) + [l1[0]]
2705    fitl.error_band = shapes.Line(cband).triangulate().lw(0).c("k", 0.15)
2706    fitl.error_band.name = "PolynomialFitErrorBand"
2707    return fitl
2708
2709
2710def _plot_fxy(
2711        z,
2712        xlim=(0, 3),
2713        ylim=(0, 3),
2714        zlim=(None, None),
2715        show_nan=True,
2716        zlevels=10,
2717        c=None,
2718        bc="aqua",
2719        alpha=1,
2720        texture="",
2721        bins=(100, 100),
2722        axes=True,
2723    ):
2724    if c is not None:
2725        texture = None  # disable
2726
2727    ps = vtk.vtkPlaneSource()
2728    ps.SetResolution(bins[0], bins[1])
2729    ps.SetNormal([0, 0, 1])
2730    ps.Update()
2731    poly = ps.GetOutput()
2732    dx = xlim[1] - xlim[0]
2733    dy = ylim[1] - ylim[0]
2734    todel, nans = [], []
2735
2736    for i in range(poly.GetNumberOfPoints()):
2737        px, py, _ = poly.GetPoint(i)
2738        xv = (px + 0.5) * dx + xlim[0]
2739        yv = (py + 0.5) * dy + ylim[0]
2740        try:
2741            with warnings.catch_warnings():
2742                warnings.simplefilter("ignore")
2743                zv = z(xv, yv)
2744                if np.isnan(zv) or np.isinf(zv) or np.iscomplex(zv):
2745                    zv = 0
2746                    todel.append(i)
2747                    nans.append([xv, yv, 0])
2748        except:
2749            zv = 0
2750            todel.append(i)
2751            nans.append([xv, yv, 0])
2752        poly.GetPoints().SetPoint(i, [xv, yv, zv])
2753
2754    if todel:
2755        cellIds = vtk.vtkIdList()
2756        poly.BuildLinks()
2757        for i in todel:
2758            poly.GetPointCells(i, cellIds)
2759            for j in range(cellIds.GetNumberOfIds()):
2760                poly.DeleteCell(cellIds.GetId(j))  # flag cell
2761        poly.RemoveDeletedCells()
2762        cl = vtk.vtkCleanPolyData()
2763        cl.SetInputData(poly)
2764        cl.Update()
2765        poly = cl.GetOutput()
2766
2767    if not poly.GetNumberOfPoints():
2768        vedo.logger.error("function is not real in the domain")
2769        return None
2770
2771    if zlim[0]:
2772        tmpact1 = Mesh(poly)
2773        a = tmpact1.cut_with_plane((0, 0, zlim[0]), (0, 0, 1))
2774        poly = a.polydata()
2775    if zlim[1]:
2776        tmpact2 = Mesh(poly)
2777        a = tmpact2.cut_with_plane((0, 0, zlim[1]), (0, 0, -1))
2778        poly = a.polydata()
2779
2780    cmap = ""
2781    if c in colors.cmaps_names:
2782        cmap = c
2783        c = None
2784        bc = None
2785
2786    mesh = Mesh(poly, c, alpha).compute_normals().lighting("plastic")
2787
2788    if cmap:
2789        mesh.compute_elevation().cmap(cmap)
2790    if bc:
2791        mesh.bc(bc)
2792    if texture:
2793        mesh.texture(texture)
2794
2795    acts = [mesh]
2796    if zlevels:
2797        elevation = vtk.vtkElevationFilter()
2798        elevation.SetInputData(poly)
2799        bounds = poly.GetBounds()
2800        elevation.SetLowPoint(0, 0, bounds[4])
2801        elevation.SetHighPoint(0, 0, bounds[5])
2802        elevation.Update()
2803        bcf = vtk.vtkBandedPolyDataContourFilter()
2804        bcf.SetInputData(elevation.GetOutput())
2805        bcf.SetScalarModeToValue()
2806        bcf.GenerateContourEdgesOn()
2807        bcf.GenerateValues(zlevels, elevation.GetScalarRange())
2808        bcf.Update()
2809        zpoly = bcf.GetContourEdgesOutput()
2810        zbandsact = Mesh(zpoly, "k", alpha).lw(1).lighting("off")
2811        zbandsact.mapper().SetResolveCoincidentTopologyToPolygonOffset()
2812        acts.append(zbandsact)
2813
2814    if show_nan and todel:
2815        bb = mesh.GetBounds()
2816        if bb[4] <= 0 and bb[5] >= 0:
2817            zm = 0.0
2818        else:
2819            zm = (bb[4] + bb[5]) / 2
2820        nans = np.array(nans) + [0, 0, zm]
2821        nansact = shapes.Points(nans, r=2, c="red5", alpha=alpha)
2822        nansact.GetProperty().RenderPointsAsSpheresOff()
2823        acts.append(nansact)
2824
2825    if isinstance(axes, dict):
2826        axs = addons.Axes(mesh, **axes)
2827        acts.append(axs)
2828    elif axes:
2829        axs = addons.Axes(mesh)
2830        acts.append(axs)
2831
2832    assem = Assembly(acts)
2833    assem.name = "PlotFxy"
2834    return assem
2835
2836
2837def _plot_fz(
2838        z,
2839        x=(-1, 1),
2840        y=(-1, 1),
2841        zlimits=(None, None),
2842        cmap="PiYG",
2843        alpha=1,
2844        lw=0.1,
2845        bins=(75, 75),
2846        axes=True,
2847    ):
2848    ps = vtk.vtkPlaneSource()
2849    ps.SetResolution(bins[0], bins[1])
2850    ps.SetNormal([0, 0, 1])
2851    ps.Update()
2852    poly = ps.GetOutput()
2853    dx = x[1] - x[0]
2854    dy = y[1] - y[0]
2855
2856    arrImg = []
2857    for i in range(poly.GetNumberOfPoints()):
2858        px, py, _ = poly.GetPoint(i)
2859        xv = (px + 0.5) * dx + x[0]
2860        yv = (py + 0.5) * dy + y[0]
2861        try:
2862            zv = z(complex(xv), complex(yv))
2863        except:
2864            zv = 0
2865        poly.GetPoints().SetPoint(i, [xv, yv, np.real(zv)])
2866        arrImg.append(np.imag(zv))
2867
2868    mesh = Mesh(poly, alpha).lighting("plastic")
2869    v = max(abs(np.min(arrImg)), abs(np.max(arrImg)))
2870    mesh.cmap(cmap, arrImg, vmin=-v, vmax=v)
2871    mesh.compute_normals().lw(lw)
2872
2873    if zlimits[0]:
2874        mesh.cut_with_plane((0, 0, zlimits[0]), (0, 0, 1))
2875    if zlimits[1]:
2876        mesh.cut_with_plane((0, 0, zlimits[1]), (0, 0, -1))
2877
2878    acts = [mesh]
2879    if axes:
2880        axs = addons.Axes(mesh, ztitle="Real part")
2881        acts.append(axs)
2882    asse = Assembly(acts)
2883    asse.name = "PlotFz"
2884    if isinstance(z, str):
2885        asse.name += " " + z
2886    return asse
2887
2888
2889def _plot_polar(
2890        rphi,
2891        title="",
2892        tsize=0.1,
2893        lsize=0.05,
2894        r1=0,
2895        r2=1,
2896        c="blue",
2897        bc="k",
2898        alpha=1,
2899        ps=5,
2900        lw=3,
2901        deg=False,
2902        vmax=None,
2903        fill=False,
2904        splined=False,
2905        nrays=8,
2906        show_disc=True,
2907        show_lines=True,
2908        show_angles=True,
2909    ):
2910    if len(rphi) == 2:
2911        rphi = np.stack((rphi[0], rphi[1]), axis=1)
2912
2913    rphi = np.array(rphi, dtype=float)
2914    thetas = rphi[:, 0]
2915    radii = rphi[:, 1]
2916
2917    k = 180 / np.pi
2918    if deg:
2919        thetas = np.array(thetas, dtype=float) / k
2920
2921    vals = []
2922    for v in thetas:  # normalize range
2923        t = np.arctan2(np.sin(v), np.cos(v))
2924        if t < 0:
2925            t += 2 * np.pi
2926        vals.append(t)
2927    thetas = np.array(vals, dtype=float)
2928
2929    if vmax is None:
2930        vmax = np.max(radii)
2931
2932    angles = []
2933    points = []
2934    for t, r in zip(thetas, radii):
2935        r = r / vmax * r2 + r1
2936        ct, st = np.cos(t), np.sin(t)
2937        points.append([r * ct, r * st, 0])
2938    p0 = points[0]
2939    points.append(p0)
2940
2941    r2e = r1 + r2
2942    lines = None
2943    if splined:
2944        lines = shapes.KSpline(points, closed=True)
2945        lines.c(c).lw(lw).alpha(alpha)
2946    elif lw:
2947        lines = shapes.Line(points)
2948        lines.c(c).lw(lw).alpha(alpha)
2949
2950    points.pop()
2951
2952    ptsact = None
2953    if ps:
2954        ptsact = shapes.Points(points, r=ps, c=c, alpha=alpha)
2955
2956    filling = None
2957    if fill and lw:
2958        faces = []
2959        coords = [[0, 0, 0]] + lines.points().tolist()
2960        for i in range(1, lines.npoints):
2961            faces.append([0, i, i + 1])
2962        filling = Mesh([coords, faces]).c(c).alpha(alpha)
2963
2964    back = None
2965    back2 = None
2966    if show_disc:
2967        back = shapes.Disc(r1=r2e, r2=r2e * 1.01, c=bc, res=(1, 360))
2968        back.z(-0.01).lighting("off").alpha(alpha)
2969        back2 = shapes.Disc(r1=r2e / 2, r2=r2e / 2 * 1.005, c=bc, res=(1, 360))
2970        back2.z(-0.01).lighting("off").alpha(alpha)
2971
2972    ti = None
2973    if title:
2974        ti = shapes.Text3D(title, (0, 0, 0), s=tsize, depth=0, justify="top-center")
2975        ti.pos(0, -r2e * 1.15, 0.01)
2976
2977    rays = []
2978    if show_disc:
2979        rgap = 0.05
2980        for t in np.linspace(0, 2 * np.pi, num=nrays, endpoint=False):
2981            ct, st = np.cos(t), np.sin(t)
2982            if show_lines:
2983                l = shapes.Line((0, 0, -0.01), (r2e * ct * 1.03, r2e * st * 1.03, -0.01))
2984                rays.append(l)
2985                ct2, st2 = np.cos(t + np.pi / nrays), np.sin(t + np.pi / nrays)
2986                lm = shapes.DashedLine(
2987                    (0, 0, -0.01), (r2e * ct2, r2e * st2, -0.01), spacing=0.25
2988                )
2989                rays.append(lm)
2990            elif show_angles:  # just the ticks
2991                l = shapes.Line(
2992                    (r2e * ct * 0.98, r2e * st * 0.98, -0.01),
2993                    (r2e * ct * 1.03, r2e * st * 1.03, -0.01),
2994                )
2995            if show_angles:
2996                if 0 <= t < np.pi / 2:
2997                    ju = "bottom-left"
2998                elif t == np.pi / 2:
2999                    ju = "bottom-center"
3000                elif np.pi / 2 < t <= np.pi:
3001                    ju = "bottom-right"
3002                elif np.pi < t < np.pi * 3 / 2:
3003                    ju = "top-right"
3004                elif t == np.pi * 3 / 2:
3005                    ju = "top-center"
3006                else:
3007                    ju = "top-left"
3008                a = shapes.Text3D(int(t * k), pos=(0, 0, 0), s=lsize, depth=0, justify=ju)
3009                a.pos(r2e * ct * (1 + rgap), r2e * st * (1 + rgap), -0.01)
3010                angles.append(a)
3011
3012    mrg = merge(back, back2, angles, rays, ti)
3013    if mrg:
3014        mrg.color(bc).alpha(alpha).lighting("off")
3015    rh = Assembly([lines, ptsact, filling] + [mrg])
3016    rh.base = np.array([0, 0, 0], dtype=float)
3017    rh.top = np.array([0, 0, 1], dtype=float)
3018    rh.name = "PlotPolar"
3019    return rh
3020
3021
3022def _plot_spheric(rfunc, normalize=True, res=33, scalarbar=True, c="grey", alpha=0.05, cmap="jet"):
3023    sg = shapes.Sphere(res=res, quads=True)
3024    sg.alpha(alpha).c(c).wireframe()
3025
3026    cgpts = sg.points()
3027    r, theta, phi = utils.cart2spher(*cgpts.T)
3028
3029    newr, inans = [], []
3030    for i in range(len(r)):
3031        try:
3032            ri = rfunc(theta[i], phi[i])
3033            if np.isnan(ri):
3034                inans.append(i)
3035                newr.append(1)
3036            else:
3037                newr.append(ri)
3038        except:
3039            inans.append(i)
3040            newr.append(1)
3041
3042    newr = np.array(newr, dtype=float)
3043    if normalize:
3044        newr = newr / np.max(newr)
3045        newr[inans] = 1
3046
3047    nanpts = []
3048    if inans:
3049        redpts = utils.spher2cart(newr[inans], theta[inans], phi[inans])
3050        nanpts.append(shapes.Points(redpts, r=4, c="r"))
3051
3052    pts = utils.spher2cart(newr, theta, phi)
3053
3054    ssurf = sg.clone().points(pts)
3055    if inans:
3056        ssurf.delete_cells_by_point_index(inans)
3057
3058    ssurf.alpha(1).wireframe(0).lw(0.1)
3059
3060    ssurf.cmap(cmap, newr)
3061    ssurf.compute_normals()
3062
3063    if scalarbar:
3064        xm = np.max([np.max(pts[0]), 1])
3065        ym = np.max([np.abs(np.max(pts[1])), 1])
3066        ssurf.mapper().SetScalarRange(np.min(newr), np.max(newr))
3067        sb3d = ssurf.add_scalarbar3d(s=(xm * 0.07, ym), c="k").scalarbar
3068        sb3d.rotate_x(90).pos(xm * 1.1, 0, -0.5)
3069    else:
3070        sb3d = None
3071
3072    sg.pickable(False)
3073    asse = Assembly([ssurf, sg] + nanpts + [sb3d])
3074    asse.name = "PlotSpheric"
3075    return asse
3076
3077
3078def _histogram_quad_bin(x,y, **kwargs):
3079    # generate a histgram with 3D bars
3080    #
3081    histo = Histogram2D(x,y, **kwargs)
3082
3083    gap = kwargs.pop("gap", 0)
3084    zscale = kwargs.pop("zscale", 1)
3085    cmap = kwargs.pop("cmap", "Blues_r")
3086
3087    gr = histo.actors[2]
3088    d = gr.diagonal_size()
3089    tol = d / 1_000_000  # tolerance
3090    if gap >= 0:
3091        gr.shrink(1 - gap - tol)
3092    gr.map_cells_to_points()
3093
3094    faces = np.array(gr.faces())
3095    s = 1 / histo.entries * len(faces) * zscale
3096    zvals = gr.pointdata["Scalars"] * s
3097
3098    pts1 = gr.points()
3099    pts2 = np.copy(pts1)
3100    pts2[:, 2] = zvals + tol
3101    newpts = np.vstack([pts1, pts2])
3102    newzvals = np.hstack([zvals, zvals]) / s
3103
3104    n = pts1.shape[0]
3105    newfaces = []
3106    for f in faces:
3107        f0, f1, f2, f3 = f
3108        f0n, f1n, f2n, f3n = f + n
3109        newfaces.extend(
3110            [
3111                [f0, f1, f2, f3],
3112                [f0n, f1n, f2n, f3n],
3113                [f0, f1, f1n, f0n],
3114                [f1, f2, f2n, f1n],
3115                [f2, f3, f3n, f2n],
3116                [f3, f0, f0n, f3n],
3117            ]
3118        )
3119
3120    msh = Mesh([newpts, newfaces]).pickable(False)
3121    msh.cmap(cmap, newzvals, name="Frequency")
3122    msh.lw(1).lighting("ambient")
3123
3124    histo.actors[2] = msh
3125    histo.RemovePart(gr)
3126    histo.AddPart(msh)
3127    return histo
3128
3129
3130def _histogram_hex_bin(
3131        xvalues,
3132        yvalues,
3133        bins=12,
3134        norm=1,
3135        fill=True,
3136        c=None,
3137        cmap="terrain_r",
3138        alpha=1,
3139    ):
3140    xmin, xmax = np.min(xvalues), np.max(xvalues)
3141    ymin, ymax = np.min(yvalues), np.max(yvalues)
3142    dx, dy = xmax - xmin, ymax - ymin
3143
3144    if utils.is_sequence(bins):
3145        n, m = bins
3146    else:
3147        if xmax - xmin < ymax - ymin:
3148            n = bins
3149            m = np.rint(dy / dx * n / 1.2 + 0.5).astype(int)
3150        else:
3151            m = bins
3152            n = np.rint(dx / dy * m * 1.2 + 0.5).astype(int)
3153
3154    src = vtk.vtkPointSource()
3155    src.SetNumberOfPoints(len(xvalues))
3156    src.Update()
3157    poly = src.GetOutput()
3158
3159    values = np.stack((xvalues, yvalues), axis=1)
3160    zs = [[0.0]] * len(values)
3161    values = np.append(values, zs, axis=1)
3162
3163    poly.GetPoints().SetData(utils.numpy2vtk(values, dtype=np.float32))
3164    cloud = Mesh(poly)
3165
3166    col = None
3167    if c is not None:
3168        col = colors.get_color(c)
3169
3170    hexs, binmax = [], 0
3171    ki, kj = 1.33, 1.12
3172    r = 0.47 / n * 1.2 * dx
3173    for i in range(n + 3):
3174        for j in range(m + 2):
3175            cyl = vtk.vtkCylinderSource()
3176            cyl.SetResolution(6)
3177            cyl.CappingOn()
3178            cyl.SetRadius(0.5)
3179            cyl.SetHeight(0.1)
3180            cyl.Update()
3181            t = vtk.vtkTransform()
3182            if not i % 2:
3183                p = (i / ki, j / kj, 0)
3184            else:
3185                p = (i / ki, j / kj + 0.45, 0)
3186            q = (p[0] / n * 1.2 * dx + xmin, p[1] / m * dy + ymin, 0)
3187            ids = cloud.closest_point(q, radius=r, return_cell_id=True)
3188            ne = len(ids)
3189            if fill:
3190                t.Translate(p[0], p[1], ne / 2)
3191                t.Scale(1, 1, ne * 10)
3192            else:
3193                t.Translate(p[0], p[1], ne)
3194            t.RotateX(90)  # put it along Z
3195            tf = vtk.vtkTransformPolyDataFilter()
3196            tf.SetInputData(cyl.GetOutput())
3197            tf.SetTransform(t)
3198            tf.Update()
3199            if c is None:
3200                col = i
3201            h = Mesh(tf.GetOutput(), c=col, alpha=alpha).flat()
3202            h.lighting("plastic")
3203            h.PickableOff()
3204            hexs.append(h)
3205            if ne > binmax:
3206                binmax = ne
3207
3208    if cmap is not None:
3209        for h in hexs:
3210            z = h.GetBounds()[5]
3211            col = colors.color_map(z, cmap, 0, binmax)
3212            h.color(col)
3213
3214    asse = Assembly(hexs)
3215    asse.SetScale(1.2 / n * dx, 1 / m * dy, norm / binmax * (dx + dy) / 4)
3216    asse.SetPosition(xmin, ymin, 0)
3217    asse.base = np.array([0, 0, 0], dtype=float)
3218    asse.top = np.array([0, 0, 1], dtype=float)
3219    asse.name = "HistogramHexBin"
3220    return asse
3221
3222
3223def _histogram_polar(
3224        values,
3225        weights=None,
3226        title="",
3227        tsize=0.1,
3228        bins=16,
3229        r1=0.25,
3230        r2=1,
3231        phigap=0.5,
3232        rgap=0.05,
3233        lpos=1,
3234        lsize=0.04,
3235        c='grey',
3236        bc="k",
3237        alpha=1,
3238        cmap=None,
3239        deg=False,
3240        vmin=None,
3241        vmax=None,
3242        labels=(),
3243        show_disc=True,
3244        nrays=8,
3245        show_lines=True,
3246        show_angles=True,
3247        show_errors=False,
3248    ):
3249    k = 180 / np.pi
3250    if deg:
3251        values = np.array(values, dtype=float) / k
3252    else:
3253        values = np.array(values, dtype=float)
3254
3255    vals = []
3256    for v in values:  # normalize range
3257        t = np.arctan2(np.sin(v), np.cos(v))
3258        if t < 0:
3259            t += 2 * np.pi
3260        vals.append(t + 0.00001)
3261
3262    histodata, edges = np.histogram(
3263        vals, weights=weights, bins=bins, range=(0, 2 * np.pi)
3264    )
3265
3266    thetas = []
3267    for i in range(bins):
3268        thetas.append((edges[i] + edges[i + 1]) / 2)
3269
3270    if vmin is None:
3271        vmin = np.min(histodata)
3272    if vmax is None:
3273        vmax = np.max(histodata)
3274
3275    errors = np.sqrt(histodata)
3276    r2e = r1 + r2
3277    if show_errors:
3278        r2e += np.max(errors) / vmax * 1.5
3279
3280    back = None
3281    if show_disc:
3282        back = shapes.Disc(r1=r2e, r2=r2e * 1.01, c=bc, res=(1, 360))
3283        back.z(-0.01)
3284
3285    slices = []
3286    lines = []
3287    angles = []
3288    errbars = []
3289
3290    for i, t in enumerate(thetas):
3291        r = histodata[i] / vmax * r2
3292        d = shapes.Disc((0, 0, 0), r1, r1 + r, res=(1, 360))
3293        delta = np.pi / bins - np.pi / 2 - phigap / k
3294        d.cut_with_plane(normal=(np.cos(t + delta), np.sin(t + delta), 0))
3295        d.cut_with_plane(normal=(np.cos(t - delta), np.sin(t - delta), 0))
3296        if cmap is not None:
3297            cslice = colors.color_map(histodata[i], cmap, vmin, vmax)
3298            d.color(cslice)
3299        else:
3300            if c is None:
3301                d.color(i)
3302            elif utils.is_sequence(c) and len(c) == bins:
3303                d.color(c[i])
3304            else:
3305                d.color(c)
3306        d.alpha(alpha).lighting("off")
3307        slices.append(d)
3308
3309        ct, st = np.cos(t), np.sin(t)
3310
3311        if show_errors:
3312            show_lines = False
3313            err = np.sqrt(histodata[i]) / vmax * r2
3314            errl = shapes.Line(
3315                ((r1 + r - err) * ct, (r1 + r - err) * st, 0.01),
3316                ((r1 + r + err) * ct, (r1 + r + err) * st, 0.01),
3317            )
3318            errl.alpha(alpha).lw(3).color(bc)
3319            errbars.append(errl)
3320
3321    labs = []
3322    rays = []
3323    if show_disc:
3324        outerdisc = shapes.Disc(r1=r2e, r2=r2e * 1.01, c=bc, res=(1, 360))
3325        outerdisc.z(-0.01)
3326        innerdisc = shapes.Disc(r1=r2e / 2, r2=r2e / 2 * 1.005, c=bc, res=(1, 360))
3327        innerdisc.z(-0.01)
3328        rays.append(outerdisc)
3329        rays.append(innerdisc)
3330
3331        rgap = 0.05
3332        for t in np.linspace(0, 2 * np.pi, num=nrays, endpoint=False):
3333            ct, st = np.cos(t), np.sin(t)
3334            if show_lines:
3335                l = shapes.Line((0, 0, -0.01), (r2e * ct * 1.03, r2e * st * 1.03, -0.01))
3336                rays.append(l)
3337                ct2, st2 = np.cos(t + np.pi / nrays), np.sin(t + np.pi / nrays)
3338                lm = shapes.DashedLine(
3339                    (0, 0, -0.01), (r2e * ct2, r2e * st2, -0.01), spacing=0.25
3340                )
3341                rays.append(lm)
3342            elif show_angles:  # just the ticks
3343                l = shapes.Line(
3344                    (r2e * ct * 0.98, r2e * st * 0.98, -0.01),
3345                    (r2e * ct * 1.03, r2e * st * 1.03, -0.01),
3346                )
3347            if show_angles:
3348                if 0 <= t < np.pi / 2:
3349                    ju = "bottom-left"
3350                elif t == np.pi / 2:
3351                    ju = "bottom-center"
3352                elif np.pi / 2 < t <= np.pi:
3353                    ju = "bottom-right"
3354                elif np.pi < t < np.pi * 3 / 2:
3355                    ju = "top-right"
3356                elif t == np.pi * 3 / 2:
3357                    ju = "top-center"
3358                else:
3359                    ju = "top-left"
3360                a = shapes.Text3D(int(t * k), pos=(0, 0, 0), s=lsize, depth=0, justify=ju)
3361                a.pos(r2e * ct * (1 + rgap), r2e * st * (1 + rgap), -0.01)
3362                angles.append(a)
3363
3364    ti = None
3365    if title:
3366        ti = shapes.Text3D(title, (0, 0, 0), s=tsize, depth=0, justify="top-center")
3367        ti.pos(0, -r2e * 1.15, 0.01)
3368
3369    for i, t in enumerate(thetas):
3370        if i < len(labels):
3371            lab = shapes.Text3D(
3372                labels[i], (0, 0, 0), s=lsize, depth=0, justify="center"  # font="VTK",
3373            )
3374            lab.pos(
3375                r2e * np.cos(t) * (1 + rgap) * lpos / 2,
3376                r2e * np.sin(t) * (1 + rgap) * lpos / 2,
3377                0.01,
3378            )
3379            labs.append(lab)
3380
3381    mrg = merge(lines, angles, rays, ti, labs)
3382    if mrg:
3383        mrg.color(bc).lighting("off")
3384
3385    acts = slices + errbars + [mrg]
3386    asse = Assembly(acts)
3387    asse.frequencies = histodata
3388    asse.bins = edges
3389    asse.name = "HistogramPolar"
3390    return asse
3391
3392def _histogram_spheric(thetavalues, phivalues, rmax=1.2, res=8, cmap="rainbow", gap=0.1):
3393
3394    x, y, z = utils.spher2cart(np.ones_like(thetavalues) * 1.1, thetavalues, phivalues)
3395    ptsvals = np.c_[x, y, z]
3396
3397    sg = shapes.Sphere(res=res, quads=True).shrink(1-gap)
3398    sgfaces = sg.faces()
3399    sgpts = sg.points()
3400
3401    cntrs = sg.cell_centers()
3402    counts = np.zeros(len(cntrs))
3403    for p in ptsvals:
3404        cell = sg.closest_point(p, return_cell_id=True)
3405        counts[cell] += 1
3406    acounts = np.array(counts, dtype=float)
3407    counts *= (rmax - 1) / np.max(counts)
3408
3409    for cell, cn in enumerate(counts):
3410        if not cn:
3411            continue
3412        fs = sgfaces[cell]
3413        pts = sgpts[fs]
3414        _, t1, p1 = utils.cart2spher(pts[:, 0], pts[:, 1], pts[:, 2])
3415        x, y, z = utils.spher2cart(1 + cn, t1, p1)
3416        sgpts[fs] = np.c_[x, y, z]
3417
3418    sg.points(sgpts)
3419    sg.cmap(cmap, acounts, on="cells")
3420    vals = sg.celldata["Scalars"]
3421
3422    faces = sg.faces()
3423    points = sg.points().tolist() + [[0.,0.,0.]]
3424    lp = len(points)-1
3425    newfaces = []
3426    newvals = []
3427    for i, f in enumerate(faces):
3428        p0,p1,p2,p3 = f
3429        newfaces.append(f)
3430        newfaces.append([p0,lp,p1])
3431        newfaces.append([p1,lp,p2])
3432        newfaces.append([p2,lp,p3])
3433        newfaces.append([p3,lp,p0])
3434        for j in range(5):
3435            newvals.append(vals[i])
3436
3437    newsg = Mesh([points, newfaces]).cmap(cmap,newvals,on='cells')
3438    newsg.compute_normals().lw(0.1).flat()
3439    newsg.name = "HistogramSpheric"
3440    return newsg
3441
3442
3443def donut(
3444        fractions,
3445        title="",
3446        tsize=0.3,
3447        r1=1.7,
3448        r2=1,
3449        phigap=0,
3450        lpos=0.8,
3451        lsize=0.15,
3452        c=None,
3453        bc="k",
3454        alpha=1,
3455        labels=(),
3456        show_disc=False,
3457    ):
3458    """
3459    Donut plot or pie chart.
3460
3461    Arguments:
3462        title : (str)
3463            plot title
3464        tsize : (float)
3465            title size
3466        r1 : (float) inner radius
3467        r2 : (float)
3468            outer radius, starting from r1
3469        phigap : (float)
3470            gap angle btw 2 radial bars, in degrees
3471        lpos : (float)
3472            label gap factor along radius
3473        lsize : (float)
3474            label size
3475        c : (color)
3476            color of the plot slices
3477        bc : (color)
3478            color of the disc frame
3479        alpha : (float)
3480            opacity of the disc frame
3481        labels : (list)
3482            list of labels
3483        show_disc : (bool)
3484            show the outer ring axis
3485
3486    Examples:
3487        - [donut.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/donut.py)
3488
3489            ![](https://vedo.embl.es/images/pyplot/donut.png)
3490    """
3491    fractions = np.array(fractions, dtype=float)
3492    angles = np.add.accumulate(2 * np.pi * fractions)
3493    angles[-1] = 2 * np.pi
3494    if angles[-2] > 2 * np.pi:
3495        print("Error in donut(): fractions must sum to 1.")
3496        raise RuntimeError
3497
3498    cols = []
3499    for i, th in enumerate(np.linspace(0, 2 * np.pi, 360, endpoint=False)):
3500        for ia, a in enumerate(angles):
3501            if th < a:
3502                cols.append(c[ia])
3503                break
3504    labs = ()
3505    if labels:
3506        angles = np.concatenate([[0], angles])
3507        labs = [""] * 360
3508        for i in range(len(labels)):
3509            a = (angles[i+1] + angles[i]) / 2
3510            j = int(a / np.pi * 180)
3511            labs[j] = labels[i]
3512
3513    data = np.linspace(0, 2 * np.pi, 360, endpoint=False) + 0.005
3514    dn = _histogram_polar(
3515        data,
3516        title=title,
3517        bins=360,
3518        r1=r1,
3519        r2=r2,
3520        phigap=phigap,
3521        lpos=lpos,
3522        lsize=lsize,
3523        tsize=tsize,
3524        c=cols,
3525        bc=bc,
3526        alpha=alpha,
3527        vmin=0,
3528        vmax=1,
3529        labels=labs,
3530        show_disc=show_disc,
3531        show_lines=0,
3532        show_angles=0,
3533        show_errors=0,
3534    )
3535    dn.name = "Donut"
3536    return dn
3537
3538
3539def violin(
3540        values,
3541        bins=10,
3542        vlim=None,
3543        x=0,
3544        width=3,
3545        splined=True,
3546        fill=True,
3547        c="violet",
3548        alpha=1,
3549        outline=True,
3550        centerline=True,
3551        lc="darkorchid",
3552        lw=3,
3553    ):
3554    """
3555    Violin style histogram.
3556
3557    Arguments:
3558        bins : (int)
3559            number of bins
3560        vlim : (list)
3561            input value limits. Crop values outside range
3562        x : (float)
3563            x-position of the violin axis
3564        width : (float)
3565            width factor of the normalized distribution
3566        splined : (bool)
3567            spline the outline
3568        fill : (bool)
3569            fill violin with solid color
3570        outline : (bool)
3571            add the distribution outline
3572        centerline : (bool)
3573            add the vertical centerline at x
3574        lc : (color)
3575            line color
3576
3577    Examples:
3578        - [histo_violin.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/histo_violin.py)
3579
3580            ![](https://vedo.embl.es/images/pyplot/histo_violin.png)
3581    """
3582    fs, edges = np.histogram(values, bins=bins, range=vlim)
3583    mine, maxe = np.min(edges), np.max(edges)
3584    fs = fs.astype(float) / len(values) * width
3585
3586    rs = []
3587
3588    if splined:
3589        lnl, lnr = [(0, edges[0], 0)], [(0, edges[0], 0)]
3590        for i in range(bins):
3591            xc = (edges[i] + edges[i + 1]) / 2
3592            yc = fs[i]
3593            lnl.append([-yc, xc, 0])
3594            lnr.append([yc, xc, 0])
3595        lnl.append((0, edges[-1], 0))
3596        lnr.append((0, edges[-1], 0))
3597        spl = shapes.KSpline(lnl).x(x)
3598        spr = shapes.KSpline(lnr).x(x)
3599        spl.color(lc).alpha(alpha).lw(lw)
3600        spr.color(lc).alpha(alpha).lw(lw)
3601        if outline:
3602            rs.append(spl)
3603            rs.append(spr)
3604        if fill:
3605            rb = shapes.Ribbon(spl, spr, c=c, alpha=alpha).lighting("off")
3606            rs.append(rb)
3607
3608    else:
3609        lns1 = [[0, mine, 0]]
3610        for i in range(bins):
3611            lns1.append([fs[i], edges[i], 0])
3612            lns1.append([fs[i], edges[i + 1], 0])
3613        lns1.append([0, maxe, 0])
3614
3615        lns2 = [[0, mine, 0]]
3616        for i in range(bins):
3617            lns2.append([-fs[i], edges[i], 0])
3618            lns2.append([-fs[i], edges[i + 1], 0])
3619        lns2.append([0, maxe, 0])
3620
3621        if outline:
3622            rs.append(shapes.Line(lns1, c=lc, alpha=alpha, lw=lw).x(x))
3623            rs.append(shapes.Line(lns2, c=lc, alpha=alpha, lw=lw).x(x))
3624
3625        if fill:
3626            for i in range(bins):
3627                p0 = (-fs[i], edges[i], 0)
3628                p1 = (fs[i], edges[i + 1], 0)
3629                r = shapes.Rectangle(p0, p1).x(p0[0] + x)
3630                r.color(c).alpha(alpha).lighting("off")
3631                rs.append(r)
3632
3633    if centerline:
3634        cl = shapes.Line([0, mine, 0.01], [0, maxe, 0.01], c=lc, alpha=alpha, lw=2).x(x)
3635        rs.append(cl)
3636
3637    asse = Assembly(rs)
3638    asse.name = "Violin"
3639    return asse
3640
3641
3642def whisker(
3643        data,
3644        s=0.25,
3645        c='k',
3646        lw=2,
3647        bc='blue',
3648        alpha=0.25,
3649        r=5,
3650        jitter=True,
3651        horizontal=False,
3652    ):
3653    """
3654    Generate a "whisker" bar from a 1-dimensional dataset.
3655
3656    Arguments:
3657        s : (float)
3658            size of the box
3659        c : (color)
3660            color of the lines
3661        lw : (float)
3662            line width
3663        bc : (color)
3664            color of the box
3665        alpha : (float)
3666            transparency of the box
3667        r : (float)
3668            point radius in pixels (use value 0 to disable)
3669        jitter : (bool)
3670            add some randomness to points to avoid overlap
3671        horizontal : (bool)
3672            set horizontal layout
3673
3674    Examples:
3675        - [whiskers.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/whiskers.py)
3676
3677            ![](https://vedo.embl.es/images/pyplot/whiskers.png)
3678    """
3679    xvals = np.zeros_like(np.asarray(data))
3680    if jitter:
3681        xjit = np.random.randn(len(xvals)) * s / 9
3682        xjit = np.clip(xjit, -s / 2.1, s / 2.1)
3683        xvals += xjit
3684
3685    dmean = np.mean(data)
3686    dq05 = np.quantile(data, 0.05)
3687    dq25 = np.quantile(data, 0.25)
3688    dq75 = np.quantile(data, 0.75)
3689    dq95 = np.quantile(data, 0.95)
3690
3691    pts = None
3692    if r:
3693        pts = shapes.Points([xvals, data], c=c, r=r)
3694
3695    rec = shapes.Rectangle([-s / 2, dq25], [s / 2, dq75], c=bc, alpha=alpha)
3696    rec.GetProperty().LightingOff()
3697    rl = shapes.Line(
3698        [[-s / 2, dq25], [s / 2, dq25], [s / 2, dq75], [-s / 2, dq75]], closed=True
3699    )
3700    l1 = shapes.Line([0, dq05, 0], [0, dq25, 0], c=c, lw=lw)
3701    l2 = shapes.Line([0, dq75, 0], [0, dq95, 0], c=c, lw=lw)
3702    lm = shapes.Line([-s / 2, dmean], [s / 2, dmean])
3703    lns = merge(l1, l2, lm, rl)
3704    asse = Assembly([lns, rec, pts])
3705    if horizontal:
3706        asse.rotate_z(-90)
3707    asse.name = "Whisker"
3708    asse.info["mean"] = dmean
3709    asse.info["quantile_05"] = dq05
3710    asse.info["quantile_25"] = dq25
3711    asse.info["quantile_75"] = dq75
3712    asse.info["quantile_95"] = dq95
3713    return asse
3714
3715
3716def streamplot(
3717    X,
3718    Y,
3719    U,
3720    V,
3721    direction="both",
3722    max_propagation=None,
3723    mode=1,
3724    lw=0.001,
3725    c=None,
3726    probes=(),
3727):
3728    """
3729    Generate a streamline plot of a vectorial field (U,V) defined at positions (X,Y).
3730    Returns a `Mesh` object.
3731
3732    Arguments:
3733        direction : (str)
3734            either "forward", "backward" or "both"
3735        max_propagation : (float)
3736            maximum physical length of the streamline
3737        lw : (float)
3738            line width in absolute units
3739        mode : (int)
3740            mode of varying the line width:
3741            - 0 - do not vary line width
3742            - 1 - vary line width by first vector component
3743            - 2 - vary line width vector magnitude
3744            - 3 - vary line width by absolute value of first vector component
3745
3746    Examples:
3747        - [plot_stream.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/plot_stream.py)
3748
3749            ![](https://vedo.embl.es/images/pyplot/plot_stream.png)
3750    """
3751    n = len(X)
3752    m = len(Y[0])
3753    if n != m:
3754        print("Limitation in streamplot(): only square grids are allowed.", n, m)
3755        raise RuntimeError()
3756
3757    xmin, xmax = X[0][0], X[-1][-1]
3758    ymin, ymax = Y[0][0], Y[-1][-1]
3759
3760    field = np.sqrt(U * U + V * V)
3761
3762    vol = vedo.Volume(field, dims=(n, n, 1))
3763
3764    uf = np.ravel(U, order="F")
3765    vf = np.ravel(V, order="F")
3766    vects = np.c_[uf, vf, np.zeros_like(uf)]
3767    vol.pointdata["StreamPlotField"] = vects
3768
3769    if len(probes) == 0:
3770        probe = shapes.Grid(pos=((n-1)/2,(n-1)/2,0), s=(n-1, n-1), res=(n-1,n-1))
3771    else:
3772        if isinstance(probes, vedo.Points):
3773            probes = probes.points()
3774        else:
3775            probes = np.array(probes, dtype=float)
3776            if len(probes[0]) == 2:
3777                probes = np.c_[probes[:, 0], probes[:, 1], np.zeros(len(probes))]
3778        sv = [(n - 1) / (xmax - xmin), (n - 1) / (ymax - ymin), 1]
3779        probes = probes - [xmin, ymin, 0]
3780        probes = np.multiply(probes, sv)
3781        probe = vedo.Points(probes)
3782
3783    stream = vedo.shapes.StreamLines(
3784        vol,
3785        probe,
3786        tubes={"radius": lw, "mode": mode,},
3787        lw=lw,
3788        max_propagation=max_propagation,
3789        direction=direction,
3790    )
3791    if c is not None:
3792        stream.color(c)
3793    else:
3794        stream.add_scalarbar()
3795    stream.lighting("off")
3796
3797    stream.scale([1 / (n - 1) * (xmax - xmin), 1 / (n - 1) * (ymax - ymin), 1])
3798    stream.shift(xmin, ymin)
3799    return stream
3800
3801
3802def matrix(
3803        M,
3804        title='Matrix',
3805        xtitle='',
3806        ytitle='',
3807        xlabels=(),
3808        ylabels=(),
3809        xrotation=0,
3810        cmap='Reds',
3811        vmin=None,
3812        vmax=None,
3813        precision=2,
3814        font='Theemim',
3815        scale=0,
3816        scalarbar=True,
3817        lc='white',
3818        lw=0,
3819        c='black',
3820        alpha=1,
3821    ):
3822    """
3823    Generate a matrix, or a 2D color-coded plot with bin labels.
3824
3825    Returns an `Assembly` object.
3826
3827    Arguments:
3828        M : (list, numpy array)
3829            the input array to visualize
3830        title : (str)
3831            title of the plot
3832        xtitle : (str)
3833            title of the horizontal colmuns
3834        ytitle : (str)
3835            title of the vertical rows
3836        xlabels : (list)
3837            individual string labels for each column. Must be of length m
3838        ylabels : (list)
3839            individual string labels for each row. Must be of length n
3840        xrotation : (float)
3841            rotation of the horizontal labels
3842        cmap : (str)
3843            color map name
3844        vmin : (float)
3845            minimum value of the colormap range
3846        vmax : (float)
3847            maximum value of the colormap range
3848        precision : (int)
3849            number of digits for the matrix entries or bins
3850        font : (str)
3851            font name. Check [available fonts here](https://vedo.embl.es/fonts).
3852
3853        scale : (float)
3854            size of the numeric entries or bin values
3855        scalarbar : (bool)
3856            add a scalar bar to the right of the plot
3857        lc : (str)
3858            color of the line separating the bins
3859        lw : (float)
3860            Width of the line separating the bins
3861        c : (str)
3862            text color
3863        alpha : (float)
3864            plot transparency
3865
3866    Examples:
3867        - [np_matrix.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/np_matrix.py)
3868
3869            ![](https://vedo.embl.es/images/pyplot/np_matrix.png)
3870    """
3871    M = np.asarray(M)
3872    n, m = M.shape
3873    gr = shapes.Grid(res=(m, n), s=(m / (m + n) * 2, n / (m + n) * 2), c=c, alpha=alpha)
3874    gr.wireframe(False).lc(lc).lw(lw)
3875
3876    matr = np.flip(np.flip(M), axis=1).ravel(order="C")
3877    gr.cmap(cmap, matr, on="cells", vmin=vmin, vmax=vmax)
3878    sbar = None
3879    if scalarbar:
3880        gr.add_scalarbar3d(title_font=font, label_font=font)
3881        sbar = gr.scalarbar
3882    labs = None
3883    if scale != 0:
3884        labs = gr.labels(
3885            on="cells",
3886            scale=scale / max(m, n),
3887            precision=precision,
3888            font=font,
3889            justify="center",
3890            c=c,
3891        )
3892        labs.z(0.001)
3893    t = None
3894    if title:
3895        if title == "Matrix":
3896            title += " " + str(n) + "x" + str(m)
3897        t = shapes.Text3D(title, font=font, s=0.04, justify="bottom-center", c=c)
3898        t.shift(0, n / (m + n) * 1.05)
3899
3900    xlabs = None
3901    if len(xlabels) == m:
3902        xlabs = []
3903        jus = "top-center"
3904        if xrotation > 44:
3905            jus = "right-center"
3906        for i in range(m):
3907            xl = shapes.Text3D(xlabels[i], font=font, s=0.02, justify=jus, c=c).rotate_z(
3908                xrotation
3909            )
3910            xl.shift((2 * i - m + 1) / (m + n), -n / (m + n) * 1.05)
3911            xlabs.append(xl)
3912
3913    ylabs = None
3914    if len(ylabels) == n:
3915        ylabels = list(reversed(ylabels))
3916        ylabs = []
3917        for i in range(n):
3918            yl = shapes.Text3D(
3919                ylabels[i], font=font, s=0.02, justify="right-center", c=c
3920            )
3921            yl.shift(-m / (m + n) * 1.05, (2 * i - n + 1) / (m + n))
3922            ylabs.append(yl)
3923
3924    xt = None
3925    if xtitle:
3926        xt = shapes.Text3D(xtitle, font=font, s=0.035, justify="top-center", c=c)
3927        xt.shift(0, -n / (m + n) * 1.05)
3928        if xlabs is not None:
3929            y0, y1 = xlabs[0].ybounds()
3930            xt.shift(0, -(y1 - y0) - 0.55 / (m + n))
3931    yt = None
3932    if ytitle:
3933        yt = shapes.Text3D(
3934            ytitle, font=font, s=0.035, justify="bottom-center", c=c
3935        ).rotate_z(90)
3936        yt.shift(-m / (m + n) * 1.05, 0)
3937        if ylabs is not None:
3938            x0, x1 = ylabs[0].xbounds()
3939            yt.shift(-(x1 - x0) - 0.55 / (m + n), 0)
3940    asse = Assembly(gr, sbar, labs, t, xt, yt, xlabs, ylabs)
3941    asse.name = "Matrix"
3942    return asse
3943
3944
3945def CornerPlot(points, pos=1, s=0.2, title="", c="b", bg="k", lines=True, dots=True):
3946    """
3947    Return a `vtkXYPlotActor` that is a plot of `x` versus `y`,
3948    where `points` is a list of `(x,y)` points.
3949
3950    Assign position following this convention:
3951
3952        - 1, topleft,
3953        - 2, topright,
3954        - 3, bottomleft,
3955        - 4, bottomright.
3956    """
3957    if len(points) == 2:  # passing [allx, ally]
3958        points = np.stack((points[0], points[1]), axis=1)
3959
3960    c = colors.get_color(c)  # allow different codings
3961    array_x = vtk.vtkFloatArray()
3962    array_y = vtk.vtkFloatArray()
3963    array_x.SetNumberOfTuples(len(points))
3964    array_y.SetNumberOfTuples(len(points))
3965    for i, p in enumerate(points):
3966        array_x.InsertValue(i, p[0])
3967        array_y.InsertValue(i, p[1])
3968    field = vtk.vtkFieldData()
3969    field.AddArray(array_x)
3970    field.AddArray(array_y)
3971    data = vtk.vtkDataObject()
3972    data.SetFieldData(field)
3973
3974    xyplot = vtk.vtkXYPlotActor()
3975    xyplot.AddDataObjectInput(data)
3976    xyplot.SetDataObjectXComponent(0, 0)
3977    xyplot.SetDataObjectYComponent(0, 1)
3978    xyplot.SetXValuesToValue()
3979    xyplot.SetAdjustXLabels(0)
3980    xyplot.SetAdjustYLabels(0)
3981    xyplot.SetNumberOfXLabels(3)
3982
3983    xyplot.GetProperty().SetPointSize(5)
3984    xyplot.GetProperty().SetLineWidth(2)
3985    xyplot.GetProperty().SetColor(colors.get_color(bg))
3986    xyplot.SetPlotColor(0, c[0], c[1], c[2])
3987
3988    xyplot.SetXTitle(title)
3989    xyplot.SetYTitle("")
3990    xyplot.ExchangeAxesOff()
3991    xyplot.SetPlotPoints(dots)
3992
3993    if not lines:
3994        xyplot.PlotLinesOff()
3995
3996    if isinstance(pos, str):
3997        spos = 2
3998        if "top" in pos:
3999            if "left" in pos:
4000                spos = 1
4001            elif "right" in pos:
4002                spos = 2
4003        elif "bottom" in pos:
4004            if "left" in pos:
4005                spos = 3
4006            elif "right" in pos:
4007                spos = 4
4008        pos = spos
4009    if pos == 1:
4010        xyplot.GetPositionCoordinate().SetValue(0.0, 0.8, 0)
4011    elif pos == 2:
4012        xyplot.GetPositionCoordinate().SetValue(0.76, 0.8, 0)
4013    elif pos == 3:
4014        xyplot.GetPositionCoordinate().SetValue(0.0, 0.0, 0)
4015    elif pos == 4:
4016        xyplot.GetPositionCoordinate().SetValue(0.76, 0.0, 0)
4017    else:
4018        xyplot.GetPositionCoordinate().SetValue(pos[0], pos[1], 0)
4019
4020    xyplot.GetPosition2Coordinate().SetValue(s, s, 0)
4021    return xyplot
4022
4023
4024def CornerHistogram(
4025        values,
4026        bins=20,
4027        vrange=None,
4028        minbin=0,
4029        logscale=False,
4030        title="",
4031        c="g",
4032        bg="k",
4033        alpha=1,
4034        pos="bottom-left",
4035        s=0.175,
4036        lines=True,
4037        dots=False,
4038        nmax=None,
4039    ):
4040    """
4041    Build a histogram from a list of values in n bins.
4042    The resulting object is a 2D actor.
4043
4044    Use `vrange` to restrict the range of the histogram.
4045
4046    Use `nmax` to limit the sampling to this max nr of entries
4047
4048    Use `pos` to assign its position:
4049        - 1, topleft,
4050        - 2, topright,
4051        - 3, bottomleft,
4052        - 4, bottomright,
4053        - (x, y), as fraction of the rendering window
4054    """
4055    if hasattr(values, "inputdata"):
4056        values = utils.vtk2numpy(values.inputdata().GetPointData().GetScalars())
4057
4058    n = values.shape[0]
4059    if nmax and nmax < n:
4060        # subsample:
4061        idxs = np.linspace(0, n, num=int(nmax), endpoint=False).astype(int)
4062        values = values[idxs]
4063
4064    fs, edges = np.histogram(values, bins=bins, range=vrange)
4065
4066    if minbin:
4067        fs = fs[minbin:-1]
4068    if logscale:
4069        fs = np.log10(fs + 1)
4070    pts = []
4071    for i in range(len(fs)):
4072        pts.append([(edges[i] + edges[i+1]) / 2, fs[i]])
4073
4074    cplot = CornerPlot(pts, pos, s, title, c, bg, lines, dots)
4075    cplot.SetNumberOfYLabels(2)
4076    cplot.SetNumberOfXLabels(3)
4077    tprop = vtk.vtkTextProperty()
4078    tprop.SetColor(colors.get_color(bg))
4079    tprop.SetFontFamily(vtk.VTK_FONT_FILE)
4080    tprop.SetFontFile(utils.get_font_path(vedo.settings.default_font))
4081    tprop.SetOpacity(alpha)
4082    cplot.SetAxisTitleTextProperty(tprop)
4083    cplot.GetProperty().SetOpacity(alpha)
4084    cplot.GetXAxisActor2D().SetLabelTextProperty(tprop)
4085    cplot.GetXAxisActor2D().SetTitleTextProperty(tprop)
4086    cplot.GetXAxisActor2D().SetFontFactor(0.55)
4087    cplot.GetYAxisActor2D().SetLabelFactor(0.0)
4088    cplot.GetYAxisActor2D().LabelVisibilityOff()
4089    return cplot
4090
4091
4092class DirectedGraph(Assembly):
4093    """
4094    Support for Directed Graphs.
4095    """
4096    def __init__(self, **kargs):
4097        """
4098        A graph consists of a collection of nodes (without postional information)
4099        and a collection of edges connecting pairs of nodes.
4100        The task is to determine the node positions only based on their connections.
4101
4102        This class is derived from class `Assembly`, and it assembles 4 Mesh objects
4103        representing the graph, the node labels, edge labels and edge arrows.
4104
4105        Arguments:
4106            c : (color)
4107                Color of the Graph
4108            n : (int)
4109                number of the initial set of nodes
4110            layout : (int, str)
4111                layout in
4112                `['2d', 'fast2d', 'clustering2d', 'circular', 'circular3d', 'cone', 'force', 'tree']`.
4113                Each of these layouts has different available options.
4114
4115        ---------------------------------------------------------------
4116        .. note:: Options for layouts '2d', 'fast2d' and 'clustering2d'
4117
4118        Arguments:
4119            seed : (int)
4120                seed of the random number generator used to jitter point positions
4121            rest_distance : (float)
4122                manually set the resting distance
4123            nmax : (int)
4124                the maximum number of iterations to be used
4125            zrange : (list)
4126                expand 2d graph along z axis.
4127
4128        ---------------------------------------------------------------
4129        .. note:: Options for layouts 'circular', and 'circular3d':
4130
4131        Arguments:
4132            radius : (float)
4133                set the radius of the circles
4134            height : (float)
4135                set the vertical (local z) distance between the circles
4136            zrange : (float)
4137                expand 2d graph along z axis
4138
4139        ---------------------------------------------------------------
4140        .. note:: Options for layout 'cone'
4141
4142        Arguments:
4143            compactness : (float)
4144                ratio between the average width of a cone in the tree,
4145                and the height of the cone.
4146            compression : (bool)
4147                put children closer together, possibly allowing sub-trees to overlap.
4148                This is useful if the tree is actually the spanning tree of a graph.
4149            spacing : (float)
4150                space between layers of the tree
4151
4152        ---------------------------------------------------------------
4153        .. note:: Options for layout 'force'
4154
4155        Arguments:
4156            seed : (int)
4157                seed the random number generator used to jitter point positions
4158            bounds : (list)
4159                set the region in space in which to place the final graph
4160            nmax : (int)
4161                the maximum number of iterations to be used
4162            three_dimensional : (bool)
4163                allow optimization in the 3rd dimension too
4164            random_initial_points : (bool)
4165                use random positions within the graph bounds as initial points
4166
4167        Examples:
4168            - [lineage_graph.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/lineage_graph.py)
4169            
4170                ![](https://vedo.embl.es/images/pyplot/graph_lineage.png)
4171            
4172            - [graph_network.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/pyplot/graph_network.py)
4173
4174                ![](https://vedo.embl.es/images/pyplot/graph_network.png)
4175        """
4176
4177        vedo.Assembly.__init__(self)
4178
4179        self.nodes = []
4180        self.edges = []
4181
4182        self._node_labels = []  # holds strings
4183        self._edge_labels = []
4184        self.edge_orientations = []
4185        self.edge_glyph_position = 0.6
4186
4187        self.zrange = 0.0
4188
4189        self.rotX = 0
4190        self.rotY = 0
4191        self.rotZ = 0
4192
4193        self.arrow_scale = 0.15
4194        self.node_label_scale = None
4195        self.node_label_justify = "bottom-left"
4196
4197        self.edge_label_scale = None
4198
4199        self.mdg = vtk.vtkMutableDirectedGraph()
4200
4201        n = kargs.pop("n", 0)
4202        for _ in range(n):
4203            self.add_node()
4204
4205        self._c = kargs.pop("c", (0.3, 0.3, 0.3))
4206
4207        self.gl = vtk.vtkGraphLayout()
4208
4209        self.font = kargs.pop("font", "")
4210
4211        s = kargs.pop("layout", "2d")
4212        if isinstance(s, int):
4213            ss = [
4214                "2d",
4215                "fast2d",
4216                "clustering2d",
4217                "circular",
4218                "circular3d",
4219                "cone",
4220                "force",
4221                "tree",
4222            ]
4223            s = ss[s]
4224        self.layout = s
4225
4226        if "2d" in s:
4227            if "clustering" in s:
4228                self.strategy = vtk.vtkClustering2DLayoutStrategy()
4229            elif "fast" in s:
4230                self.strategy = vtk.vtkFast2DLayoutStrategy()
4231            else:
4232                self.strategy = vtk.vtkSimple2DLayoutStrategy()
4233            self.rotX = 180
4234            opt = kargs.pop("rest_distance", None)
4235            if opt is not None:
4236                self.strategy.SetRestDistance(opt)
4237            opt = kargs.pop("seed", None)
4238            if opt is not None:
4239                self.strategy.SetRandomSeed(opt)
4240            opt = kargs.pop("nmax", None)
4241            if opt is not None:
4242                self.strategy.SetMaxNumberOfIterations(opt)
4243            self.zrange = kargs.pop("zrange", 0)
4244
4245        elif "circ" in s:
4246            if "3d" in s:
4247                self.strategy = vtk.vtkSimple3DCirclesStrategy()
4248                self.strategy.SetDirection(0, 0, -1)
4249                self.strategy.SetAutoHeight(True)
4250                self.strategy.SetMethod(1)
4251                self.rotX = -90
4252                opt = kargs.pop("radius", None)  # float
4253                if opt is not None:
4254                    self.strategy.SetMethod(0)
4255                    self.strategy.SetRadius(opt)  # float
4256                opt = kargs.pop("height", None)
4257                if opt is not None:
4258                    self.strategy.SetAutoHeight(False)
4259                    self.strategy.SetHeight(opt)  # float
4260            else:
4261                self.strategy = vtk.vtkCircularLayoutStrategy()
4262                self.zrange = kargs.pop("zrange", 0)
4263
4264        elif "cone" in s:
4265            self.strategy = vtk.vtkConeLayoutStrategy()
4266            self.rotX = 180
4267            opt = kargs.pop("compactness", None)
4268            if opt is not None:
4269                self.strategy.SetCompactness(opt)
4270            opt = kargs.pop("compression", None)
4271            if opt is not None:
4272                self.strategy.SetCompression(opt)
4273            opt = kargs.pop("spacing", None)
4274            if opt is not None:
4275                self.strategy.SetSpacing(opt)
4276
4277        elif "force" in s:
4278            self.strategy = vtk.vtkForceDirectedLayoutStrategy()
4279            opt = kargs.pop("seed", None)
4280            if opt is not None:
4281                self.strategy.SetRandomSeed(opt)
4282            opt = kargs.pop("bounds", None)
4283            if opt is not None:
4284                self.strategy.SetAutomaticBoundsComputation(False)
4285                self.strategy.SetGraphBounds(opt)  # list
4286            opt = kargs.pop("nmax", None)
4287            if opt is not None:
4288                self.strategy.SetMaxNumberOfIterations(opt)  # int
4289            opt = kargs.pop("three_dimensional", True)
4290            if opt is not None:
4291                self.strategy.SetThreeDimensionalLayout(opt)  # bool
4292            opt = kargs.pop("random_initial_points", None)
4293            if opt is not None:
4294                self.strategy.SetRandomInitialPoints(opt)  # bool
4295
4296        elif "tree" in s:
4297            self.strategy = vtk.vtkSpanTreeLayoutStrategy()
4298            self.rotX = 180
4299
4300        else:
4301            vedo.logger.error(f"Cannot understand layout {s}. Available layouts:")
4302            vedo.logger.error(
4303                "[2d,fast2d,clustering2d,circular,circular3d,cone,force,tree]"
4304            )
4305            raise RuntimeError()
4306
4307        self.gl.SetLayoutStrategy(self.strategy)
4308
4309        if len(kargs):
4310            vedo.logger.error(f"Cannot understand options: {kargs}")
4311
4312
4313    def add_node(self, label="id"):
4314        """Add a new node to the `Graph`."""
4315        v = self.mdg.AddVertex()  # vtk calls it vertex..
4316        self.nodes.append(v)
4317        if label == "id":
4318            label = int(v)
4319        self._node_labels.append(str(label))
4320        return v
4321
4322    def add_edge(self, v1, v2, label=""):
4323        """Add a new edge between to nodes.
4324        An extra node is created automatically if needed."""
4325        nv = len(self.nodes)
4326        if v1 >= nv:
4327            for _ in range(nv, v1 + 1):
4328                self.add_node()
4329        nv = len(self.nodes)
4330        if v2 >= nv:
4331            for _ in range(nv, v2 + 1):
4332                self.add_node()
4333        e = self.mdg.AddEdge(v1, v2)
4334        self.edges.append(e)
4335        self._edge_labels.append(str(label))
4336        return e
4337
4338    def add_child(self, v, node_label="id", edge_label=""):
4339        """Add a new edge to a new node as its child.
4340        The extra node is created automatically if needed."""
4341        nv = len(self.nodes)
4342        if v >= nv:
4343            for _ in range(nv, v + 1):
4344                self.add_node()
4345        child = self.mdg.AddChild(v)
4346        self.edges.append((v, child))
4347        self.nodes.append(child)
4348        if node_label == "id":
4349            node_label = int(child)
4350        self._node_labels.append(str(node_label))
4351        self._edge_labels.append(str(edge_label))
4352        return child
4353
4354    def build(self):
4355        """
4356        Build the `DirectedGraph(Assembly)`.
4357        Accessory objects are also created for labels and arrows.
4358        """
4359        self.gl.SetZRange(self.zrange)
4360        self.gl.SetInputData(self.mdg)
4361        self.gl.Update()
4362
4363        gr2poly = vtk.vtkGraphToPolyData()
4364        gr2poly.EdgeGlyphOutputOn()
4365        gr2poly.SetEdgeGlyphPosition(self.edge_glyph_position)
4366        gr2poly.SetInputData(self.gl.GetOutput())
4367        gr2poly.Update()
4368
4369        dgraph = Mesh(gr2poly.GetOutput(0))
4370        # dgraph.clean() # WRONG!!! dont uncomment
4371        dgraph.flat().color(self._c).lw(2)
4372        dgraph.name = "DirectedGraph"
4373
4374        diagsz = self.diagonal_size() / 1.42
4375        if not diagsz:
4376            return None
4377
4378        dgraph.SetScale(1 / diagsz)
4379        if self.rotX:
4380            dgraph.rotate_x(self.rotX)
4381        if self.rotY:
4382            dgraph.rotate_y(self.rotY)
4383        if self.rotZ:
4384            dgraph.rotate_z(self.rotZ)
4385
4386        vecs = gr2poly.GetOutput(1).GetPointData().GetVectors()
4387        self.edge_orientations = utils.vtk2numpy(vecs)
4388
4389        # Use Glyph3D to repeat the glyph on all edges.
4390        arrows = None
4391        if self.arrow_scale:
4392            arrow_source = vtk.vtkGlyphSource2D()
4393            arrow_source.SetGlyphTypeToEdgeArrow()
4394            arrow_source.SetScale(self.arrow_scale)
4395            arrow_source.Update()
4396            arrow_glyph = vtk.vtkGlyph3D()
4397            arrow_glyph.SetInputData(0, gr2poly.GetOutput(1))
4398            arrow_glyph.SetInputData(1, arrow_source.GetOutput())
4399            arrow_glyph.Update