vedo.pyplot

Advanced plotting functionalities.

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