diff --git a/docs/api/index.rst b/docs/api/index.rst index 4d28051e..de77ef5b 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -14,6 +14,7 @@ PhasorPy library version |version|. phasorpy phasor + plot io color datasets diff --git a/docs/api/plot.rst b/docs/api/plot.rst new file mode 100644 index 00000000..7e77c5f6 --- /dev/null +++ b/docs/api/plot.rst @@ -0,0 +1,5 @@ +phasorpy.plot +------------- + +.. automodule:: phasorpy.plot + :members: diff --git a/docs/conf.py b/docs/conf.py index ebdea424..1996d030 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -35,8 +35,10 @@ 'sphinx.ext.doctest', 'sphinx.ext.viewcode', 'sphinx.ext.todo', + # don't enable intersphinx since tutorials are getting littered with links # 'sphinx.ext.intersphinx', # 'numpydoc', + 'sphinx_inline_tabs', 'sphinx_copybutton', 'sphinx_click', 'sphinx_issues', @@ -59,7 +61,6 @@ # html_favicon = '' pygments_style = 'sphinx' -intersphinx_mapping = {'https://docs.python.org/': None} # extension configurations @@ -93,9 +94,22 @@ 'filename_pattern': 'phasorpy_', 'examples_dirs': '../tutorials', 'gallery_dirs': 'tutorials', + 'reference_url': {'phasorpy': None}, } copybutton_prompt_text = ( r'>>> |\.\.\. |\$ |In \[\d*\]: | {2,5}\.\.\.: | {5,8}: ' ) copybutton_prompt_is_regexp = True + +intersphinx_mapping = { + 'python': ('https://docs.python.org/3/', None), + 'numpy': ('https://numpy.org/doc/stable/', None), + 'matplotlib': ('https://matplotlib.org/stable/', None), + 'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None), + 'pandas': ('https://pandas.pydata.org/docs/', None), + 'sklearn': ('https://scikit-learn.org/stable/', None), + 'skimage': ('https://scikit-image.org/docs/stable/', None), +} + +intersphinx_disabled_reftypes = ['*'] diff --git a/requirements_dev.txt b/requirements_dev.txt index e3fe6741..ee3f98fb 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -24,6 +24,7 @@ sphinx-issues sphinx_gallery sphinx-copybutton sphinx_click +sphinx-inline-tabs numpydoc pydata-sphinx-theme # @@ -44,8 +45,8 @@ codespell # pre-commit # # optional requirements -# scipy -# scikit-image +scipy +scikit-image # scikit-learn # pandas # zarr diff --git a/src/phasorpy/_utils.py b/src/phasorpy/_utils.py new file mode 100644 index 00000000..44db6d44 --- /dev/null +++ b/src/phasorpy/_utils.py @@ -0,0 +1,262 @@ +"""Private utility functions. + +The ``phasorpy._utils`` module provides private auxiliary and convenience +functions. + +""" + +from __future__ import annotations + +__all__: list[str] = [ + 'parse_kwargs', + 'update_kwargs', + 'kwargs_notnone', + 'scale_matrix', + 'sort_coordinates', + 'phasor_to_polar_scalar', + 'phasor_from_polar_scalar', + 'circle_line_intersection', + 'circle_circle_intersection', +] + +import math +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from ._typing import Any, Sequence, NDArray + +import numpy + + +def parse_kwargs( + kwargs: dict[str, Any], + /, + *keys: str, + _del: bool = True, + **keyvalues: Any, +) -> dict[str, Any]: + """Return dict with keys from keys|keyvals and values from kwargs|keyvals. + + If `_del` is true (default), existing keys are deleted from `kwargs`. + + >>> kwargs = {'one': 1, 'two': 2, 'four': 4} + >>> kwargs2 = parse_kwargs(kwargs, 'two', 'three', four=None, five=5) + >>> kwargs == {'one': 1} + True + >>> kwargs2 == {'two': 2, 'four': 4, 'five': 5} + True + + """ + result = {} + for key in keys: + if key in kwargs: + result[key] = kwargs[key] + if _del: + del kwargs[key] + for key, value in keyvalues.items(): + if key in kwargs: + result[key] = kwargs[key] + if _del: + del kwargs[key] + else: + result[key] = value + return result + + +def update_kwargs(kwargs: dict[str, Any], /, **keyvalues: Any) -> None: + """Update dict with keys and values if keys do not already exist. + + >>> kwargs = {'one': 1, } + >>> update_kwargs(kwargs, one=None, two=2) + >>> kwargs == {'one': 1, 'two': 2} + True + + """ + for key, value in keyvalues.items(): + if key not in kwargs: + kwargs[key] = value + + +def kwargs_notnone(**kwargs: Any) -> dict[str, Any]: + """Return dict of kwargs which values are not None. + + >>> kwargs_notnone(one=1, none=None) + {'one': 1} + + """ + return dict(item for item in kwargs.items() if item[1] is not None) + + +def scale_matrix(factor: float, origin: Sequence[float]) -> NDArray[Any]: + """Return matrix to scale homogeneous coordinates by factor around origin. + + Parameters + ---------- + factor: float + Scale factor. + origin: (float, float) + Coordinates of point around which to scale. + + Returns + ------- + matrix: ndarray + A 3x3 homogeneous transformation matrix. + + Examples + -------- + >>> scale_matrix(1.1, (0.0, 0.5)) + array([[1.1, 0, -0], + [0, 1.1, -0.05], + [0, 0, 1]]) + + """ + mat = numpy.diag((factor, factor, 1.0)) + mat[:2, 2] = origin[:2] + mat[:2, 2] *= 1.0 - factor + return mat + + +def sort_coordinates( + real: Sequence[float], + imag: Sequence[float], + /, + origin: tuple[float, float] | None = None, +) -> tuple[NDArray[Any], NDArray[Any]]: + """Return cartesian coordinates sorted counterclockwise around origin. + + Parameters + ---------- + real, imag : array_like + Coordinates to be sorted. + origin : (float, float) + Coordinates around which to sort by angle. + + Returns + ------- + real, imag : ndarray + Coordinates sorted by angle. + + Examples + -------- + >>> sort_coordinates([0, 1, 2, 3], [0, 1, -1, 0]) + (array([2, 3, 1, 0]), array([-1, 0, 1, 0])) + + """ + x = numpy.asanyarray(real) + y = numpy.asanyarray(imag) + if x.ndim != 1 or x.shape != y.shape: + raise ValueError(f'invalid {x.shape=} or {y.shape=}') + if x.size < 4: + return x, y + if origin is None: + origin = x.mean(), y.mean() + indices = numpy.argsort(numpy.arctan2(y - origin[1], x - origin[0])) + return x[indices], y[indices] + + +def phasor_to_polar_scalar( + real: float, + imag: float, + /, + *, + degree: bool = False, + percent: bool = False, +) -> tuple[float, float]: + """Return polar from scalar phasor coordinates. + + >>> phasor_to_polar_scalar(1.0, 0.0, degree=True, percent=True) + (0.0, 100.0) + + """ + phi = math.atan2(imag, real) + mod = math.hypot(imag, real) + if degree: + phi = math.degrees(phi) + if percent: + mod *= 100.0 + return phi, mod + + +def phasor_from_polar_scalar( + phase: float, + modulation: float, + /, + *, + degree: bool = False, + percent: bool = False, +) -> tuple[float, float]: + """Return phasor from scalar polar coordinates. + + >>> phasor_from_polar_scalar(0.0, 100.0, degree=True, percent=True) + (1.0, 0.0) + + """ + if degree: + phase = math.radians(phase) + if percent: + modulation /= 100.0 + real = modulation * math.cos(phase) + imag = modulation * math.sin(phase) + return real, imag + + +def circle_circle_intersection( + x0: float, y0: float, r0: float, x1: float, y1: float, r1: float, / +) -> tuple[tuple[float, float], ...]: + """Return coordinates of intersection points of two circles. + + >>> circle_circle_intersection( + ... 0.0, 0.0, math.hypot(0.6, 0.4), 0.6, 0.4, 0.2 + ... ) # doctest: +NUMBER + ((0.6868, 0.2198), (0.4670, 0.5494)) + >>> circle_circle_intersection(0.0, 0.0, 1.0, 0.6, 0.4, 0.2) + () + + """ + dx = x1 - x0 + dy = y1 - y0 + dr = math.hypot(dx, dy) + ll = (r0 * r0 - r1 * r1 + dr * dr) / (dr + dr) + dd = r0 * r0 - ll * ll + if dd < 0.0 or dr < 1e-16: + return tuple() # no solution + hd = math.sqrt(dd) / dr + ld = ll / dr + return ( + (ld * dx + hd * dy + x0, ld * dy - hd * dx + y0), + (ld * dx - hd * dy + x0, ld * dy + hd * dx + y0), + ) + + +def circle_line_intersection( + x: float, y: float, r: float, x0: float, y0: float, x1: float, y1: float, / +) -> tuple[tuple[float, float], ...]: + """Return coordinates of intersection points of circle and line. + + >>> circle_line_intersection( + ... 0.6, 0.4, 0.2, 0.0, 0.0, 0.6, 0.4 + ... ) # doctest: +NUMBER + ((0.7664, 0.5109), (0.4335, 0.2890)) + >>> circle_line_intersection(0.6, 0.4, 0.2, 0.0, 0.0, 0.6, 0.1) + () + + """ + dx = x1 - x0 + dy = y1 - y0 + dr = dx * dx + dy * dy + dd = (x0 - x) * (y1 - y) - (x1 - x) * (y0 - y) + rdd = r * r * dr - dd * dd # discriminant + if rdd < 0 or dr < 1e-16: + return tuple() # no intersection + rdd = math.sqrt(rdd) + sgn = math.copysign + return ( + ( + x + (dd * dy + sgn(1, dy) * dx * rdd) / dr, + y + (-dd * dx + abs(dy) * rdd) / dr, + ), + ( + x + (dd * dy - sgn(1, dy) * dx * rdd) / dr, + y + (-dd * dx - abs(dy) * rdd) / dr, + ), + ) diff --git a/src/phasorpy/plot.py b/src/phasorpy/plot.py new file mode 100644 index 00000000..c0eabcd0 --- /dev/null +++ b/src/phasorpy/plot.py @@ -0,0 +1,1205 @@ +"""Plot phasor coordinates and related data. + +The ``phasorpy.plot`` module provides functions and classes to visualize +phasor coordinates and related data using the matplotlib library. + +""" + +from __future__ import annotations + +__all__ = [ + 'PhasorPlot', + 'plot_phasor', + 'plot_phasor_image', + 'plot_signal_image', + 'plot_polar_frequency', +] + +import math +import os +from collections.abc import Sequence +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from ._typing import Any, ArrayLike, NDArray, Literal, BinaryIO + + from matplotlib.axes import Axes + from matplotlib.image import AxesImage + from matplotlib.figure import Figure + +import numpy +from matplotlib import pyplot +from matplotlib.font_manager import FontProperties +from matplotlib.gridspec import GridSpec +from matplotlib.lines import Line2D +from matplotlib.patches import Arc, Circle, Polygon +from matplotlib.path import Path +from matplotlib.patheffects import AbstractPathEffect + +from ._utils import ( + circle_circle_intersection, + circle_line_intersection, + parse_kwargs, + phasor_from_polar_scalar, + phasor_to_polar_scalar, + sort_coordinates, + update_kwargs, +) +from .phasor import phasor_calibrate, phasor_from_lifetime + +GRID_COLOR = '0.5' +GRID_LINESTYLE = ':' +GRID_LINESTYLE_MAJOR = '-' +GRID_LINEWIDH = 1.0 +GRID_LINEWIDH_MINOR = 0.5 +GRID_FILL = False + + +class PhasorPlot: + """Phasor plot. + + Create publication quality visualizations of phasor coordinates. + + Parameters + ---------- + allquadrants : bool, optional + Show all quandrants of phasor space. + By default, only the first quadrant with universal semicricle is shown. + ax : matplotlib axes, optional + Matplotlib axes used for plotting. + By default, a new subplot axes is created. + frequency : float, optional + Laser pulse or modulation frequency in MHz. + grid : bool, optional, default: True + Display polar grid or semicircle. + **kwargs + Additional properties to set on `ax`. + + """ + + _ax: Axes + """Matplotlib axes.""" + + _limits: tuple[tuple[float, float], tuple[float, float]] + """Axes limits (xmin, xmax), (ymin, ymax).""" + + _full: bool + """Show all quadrants of phasor space.""" + + def __init__( + self, + /, + allquadrants: bool | None = None, + ax: Axes | None = None, + *, + frequency: float | None = None, + grid: bool = True, + **kwargs, + ) -> None: + # initialize empty phasor plot + self._ax = pyplot.subplots()[1] if ax is None else ax + + self._full = bool(allquadrants) + if self._full: + xlim = (-1.05, 1.05) + ylim = (-1.05, 1.05) + xticks: tuple[float, ...] = (-1.0, -0.5, 0.0, 0.5, 1.0) + yticks: tuple[float, ...] = (-1.0, -0.5, 0.0, 0.5, 1.0) + if grid: + self.polar_grid() + else: + xlim = (-0.05, 1.05) + ylim = (-0.05, 0.7) + xticks = (0.0, 0.2, 0.4, 0.6, 0.8, 1.0) + yticks = (0.0, 0.2, 0.4, 0.6) + if grid: + self.semicircle(frequency=frequency) + + title = 'Phasor plot' + if frequency is not None: + title += f' ({frequency:g} MHz)' + + update_kwargs( + kwargs, + title=title, + xlabel='G, real', + ylabel='S, imag', + aspect='equal', + xlim=xlim, + ylim=ylim, + xticks=xticks, + yticks=yticks, + ) + self._limits = (kwargs['xlim'], kwargs['ylim']) + self._ax.set(**kwargs) + + @property + def ax(self) -> Axes: + """Matplotlib :py:class:`matplotlib.axes.Axes`.""" + return self._ax + + @property + def fig(self) -> Figure | None: + """Matplotlib :py:class:`matplotlib.figure.Figure`.""" + return self._ax.get_figure() + + def show(self) -> None: + """Display all open figures. Call :py:func:`matplotlib.pyplot.show`.""" + # self.fig.show() + pyplot.show() + + def save( + self, + file: str | os.PathLike[Any] | BinaryIO | None, + /, + **kwargs: Any, + ) -> None: + """Save current figure to file. + + Parameters + ---------- + file : str, path-like, or binary file-like + Path or Python file-like object to write the current figure to. + **kwargs + Additional keyword arguments passed to + :py:func:`matplotlib:pyplot.savefig`. + + """ + pyplot.savefig(file, **kwargs) + + def plot( + self, + real: ArrayLike, + imag: ArrayLike, + /, + fmt='o', + *, + label: str | Sequence[str] | None = None, + **kwargs: Any, + ) -> None: + """Plot imag versus real coordinates as markers and/or lines. + + Parameters + ---------- + real : array_like + Real component of phasor coordinates. + Must be one or two dimensional. + imag : array_like + Imaginary component of phasor coordinates. + Must be of same shape as `real`. + fmt : str, optional, default: 'o' + Matplotlib style format string. + label : str or sequence of str, optional + Plot label. + May be a sequence if phasor coordinates are two dimensional arrays. + **kwargs + Additional parameters passed to + :py:meth:`matplotlib.axes.Axes.plot`. + + """ + ax = self._ax + if label is not None and ( + isinstance(label, str) or not isinstance(label, Sequence) + ): + label = (label,) + for ( + i, + (re, im), + ) in enumerate( + zip(numpy.array(real, ndmin=2), numpy.array(imag, ndmin=2)) + ): + lbl = None + if label is not None: + try: + lbl = label[i] + except IndexError: + pass + ax.plot(re, im, fmt, label=lbl, **kwargs) + if label is not None: + ax.legend() + + def _histogram2d( + self, + real: ArrayLike, + imag: ArrayLike, + /, + **kwargs: Any, + ) -> tuple[NDArray[Any], NDArray[Any], NDArray[Any]]: + """Return 2D histogram of imag versus real coordinates.""" + update_kwargs(kwargs, range=self._limits) + (xmin, xmax), (ymin, ymax) = kwargs['range'] + assert xmax > xmin and ymax > ymin + bins = kwargs.get('bins', 128) + if isinstance(bins, int): + assert bins > 0 + aspect = (xmax - xmin) / (ymax - ymin) + if aspect > 1: + bins = (bins, max(int(bins / aspect), 1)) + else: + bins = (max(int(bins * aspect), 1), bins) + kwargs['bins'] = bins + return numpy.histogram2d( + numpy.asanyarray(real).reshape(-1), + numpy.asanyarray(imag).reshape(-1), + **kwargs, + ) + + def hist2d( + self, + real: ArrayLike, + imag: ArrayLike, + /, + **kwargs: Any, + ) -> None: + """Plot 2D histogram of imag versus real coordinates. + + Parameters + ---------- + real : array_like + Real component of phasor coordinates. + imag : array_like + Imaginary component of phasor coordinates. + Must be of same shape as `real`. + **kwargs + Additional parameters passed to :py:meth:`numpy.histogram2d` + and :py:meth:`matplotlib.axes.Axes.pcolormesh`. + + """ + kwargs_hist2d = parse_kwargs( + kwargs, 'bins', 'range', 'density', 'weights' + ) + h, xedges, yedges = self._histogram2d(real, imag, **kwargs_hist2d) + + update_kwargs(kwargs, cmap='Blues', norm='log') + cmin = kwargs.pop('cmin', 1) + cmax = kwargs.pop('cmax', None) + if cmin is not None: + h[h < cmin] = None + if cmax is not None: + h[h > cmax] = None + self._ax.pcolormesh(xedges, yedges, h.T, **kwargs) + + def contour( + self, + real: ArrayLike, + imag: ArrayLike, + /, + **kwargs: Any, + ) -> None: + """Plot contours of imag versus real coordinates (not implemented). + + Parameters + ---------- + real : array_like + Real component of phasor coordinates. + imag : array_like + Imaginary component of phasor coordinates. + Must be of same shape as `real`. + **kwargs + Additional parameters passed to :py:func:`numpy.histogram2d` + and :py:meth:`matplotlib.axes.Axes.contour`. + + """ + update_kwargs(kwargs, cmap='Blues', norm='log') + kwargs_hist2d = parse_kwargs( + kwargs, 'bins', 'range', 'density', 'weights' + ) + h, xedges, yedges = self._histogram2d(real, imag, **kwargs_hist2d) + xedges = xedges[:-1] + (xedges[1] - xedges[0]) + yedges = yedges[:-1] + (yedges[1] - yedges[0]) + self._ax.contour(xedges, yedges, h.T, **kwargs) + + def imshow( + self, + image: ArrayLike, + /, + **kwargs: Any, + ) -> None: + """Plot an image, for example, a 2D histogram (not implemented). + + Parameters + ---------- + image : array_like + Image to display. + **kwargs + Additional parameters passed to + :py:meth:`matplotlib.axes.Axes.imshow`. + + """ + raise NotImplementedError + + def components( + self, + real: Sequence[float], + imag: Sequence[float], + /, + fraction: Sequence[float] | None = None, + **kwargs: Any, + ) -> None: + """Plot linear combinations of phasor coordinates or ranges thereof. + + Parameters + ---------- + real : sequence of float + Real component of phasor coordinates. + imag : sequence of float + Imaginary component of phasor coordinates. + fraction: sequence of float, optional + Weight associated with each component. + If None (default), outline the polygon area of possible linear + combinations of components. + Else, draw lines from the component coordinates to the weighted + average. + **kwargs + Additional parameters passed to + :py:class:`matplotlib.patches.Polygon` or + :py:class:`matplotlib.lines.Line2D`. + + """ + if fraction is None: + update_kwargs( + kwargs, + edgecolor=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH, + fill=GRID_FILL, + ) + self._ax.add_patch( + Polygon(numpy.vstack(sort_coordinates(real, imag)).T, **kwargs) + ) + return + + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH, + ) + center_re, center_im = numpy.average( + numpy.vstack((real, imag)), axis=-1, weights=fraction + ) + for re, im in zip(real, imag): + self._ax.add_line( + Line2D([center_re, re], [center_im, im], **kwargs) + ) + # TODO: add fraction labels? + + def line( + self, + real: ArrayLike, + imag: ArrayLike, + /, + **kwargs: Any, + ) -> None: + """Draw grid line. + + Parameters + ---------- + real : array_like, shape (n, ) + Real components of line start and end coordinates. + imag : array_like, shape (n, ) + Imaginary components of line start and end coordinates. + **kwargs + Additional parameters passed to + :py:class:`matplotlib.lines.Line2D`. + + """ + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH, + ) + self._ax.add_line(Line2D(real, imag, **kwargs)) + + def circle( + self, + real: float, + imag: float, + /, + radius: float, + **kwargs: Any, + ) -> None: + """Draw grid circle of radius around center. + + Parameters + ---------- + real : float + Real component of circle center coordinate. + imag : float + Imaginary component of circle center coordinate. + radius : float + Circle radius. + **kwargs + Additional parameters passed to + :py:class:`matplotlib.patches.Circle`. + + """ + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH, + fill=GRID_FILL, + ) + self._ax.add_patch(Circle((real, imag), radius, **kwargs)) + + def polar_cursor( + self, + phase: float | None = None, + modulation: float | None = None, + phase_limit: float | None = None, + modulation_limit: float | None = None, + radius: float | None = None, + **kwargs: Any, + ) -> None: + """Plot phase and modulation grid lines and arcs. + + Parameters + ---------- + phase : float, optional + Angular component of polar coordinate in radians. + modulation : float, optional + Radial component of polar coordinate. + phase_limit : float, optional + Angular component of limiting polar coordinate (in radians). + Modulation grid arcs are drawn between `phase` and `phase_limit`. + modulation_limit : float, optional + Radial component of limiting polar coordinate. + Phase grid lines are drawn from `modulation` to `modulation_limit`. + radius : float, optional + Radius of circle limiting phase and modulation grid lines and arcs. + **kwargs + Additional parameters passed to + :py:class:`matplotlib.lines.Line2D`, + :py:class:`matplotlib.patches.Circle`, and + :py:class:`matplotlib.patches.Arc`. + + """ + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH, + fill=GRID_FILL, + ) + ax = self._ax + if radius is not None and phase is not None and modulation is not None: + x = modulation * math.cos(phase) + y = modulation * math.sin(phase) + ax.add_patch(Circle((x, y), radius, **kwargs)) + del kwargs['fill'] + p0, p1 = circle_line_intersection(x, y, radius, 0, 0, x, y) + ax.add_line(Line2D((p0[0], p1[0]), (p0[1], p1[1]), **kwargs)) + p0, p1 = circle_circle_intersection(0, 0, modulation, x, y, radius) + ax.add_patch( + Arc( + (0, 0), + modulation * 2, + modulation * 2, + theta1=math.degrees(math.atan(p0[1] / p0[0])), + theta2=math.degrees(math.atan(p1[1] / p1[0])), + fill=False, + **kwargs, + ) + ) + return + + del kwargs['fill'] + for phi in (phase, phase_limit): + if phi is not None: + if modulation is not None and modulation_limit is not None: + x0 = modulation * math.cos(phi) + y0 = modulation * math.sin(phi) + x1 = modulation_limit * math.cos(phi) + y1 = modulation_limit * math.sin(phi) + else: + x0 = 0 + y0 = 0 + x1 = math.cos(phi) + y1 = math.sin(phi) + ax.add_line(Line2D((x0, x1), (y0, y1), **kwargs)) + for mod in (modulation, modulation_limit): + if mod is not None: + if phase is not None and phase_limit is not None: + theta1 = math.degrees(min(phase, phase_limit)) + theta2 = math.degrees(max(phase, phase_limit)) + else: + theta1 = 0.0 + theta2 = 360.0 if self._full else 90.0 + ax.add_patch( + Arc( + (0, 0), + mod * 2, + mod * 2, + theta1=theta1, + theta2=theta2, + fill=False, # filling arc objects is not supported + **kwargs, + ) + ) + + def polar_grid(self, **kwargs) -> None: + """Draw polar coordinate system. + + Parameters + ---------- + **kwargs + Parameters passed to + :py:class:`matplotlib.patches.Circle` and + :py:class:`matplotlib.lines.Line2D`. + + """ + ax = self._ax + # major gridlines + kwargs_copy = kwargs.copy() + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE_MAJOR, + linewidth=GRID_LINEWIDH, + # fill=GRID_FILL, + ) + ax.add_line(Line2D([-1, 1], [0, 0], **kwargs)) + ax.add_line(Line2D([0, 0], [-1, 1], **kwargs)) + ax.add_patch(Circle((0, 0), 1, fill=False, **kwargs)) + # minor gridlines + kwargs = kwargs_copy + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE, + linewidth=GRID_LINEWIDH_MINOR, + ) + for r in (1 / 3, 2 / 3): + ax.add_patch(Circle((0, 0), r, fill=False, **kwargs)) + for a in (3, 6): + x = math.cos(math.pi / a) + y = math.sin(math.pi / a) + ax.add_line(Line2D([-x, x], [-y, y], **kwargs)) + ax.add_line(Line2D([-x, x], [y, -y], **kwargs)) + + def semicircle( + self, + frequency: float | None = None, + *, + polar_reference: tuple[float, float] | None = None, + phasor_reference: tuple[float, float] | None = None, + lifetime: Sequence[float] | None = None, + labels: Sequence[str] | None = None, + show_circle: bool = True, + **kwargs, + ) -> None: + """Draw universal semicircle. + + Parameters + ---------- + frequency : float, optional + Laser pulse or modulation frequency in MHz. + polar_reference : (float, float), optional, default: (0, 1) + Polar coordinates of zero lifetime. + phasor_reference : (float, float), optional, default: (1, 0) + Phasor coordinates of zero lifetime. + Alternative to `polar_reference`. + lifetime : sequence of float, optional + Apparent single lifetimes at which to draw ticks and labels. + Only applies when `frequency` is specified. + labels : sequence of str, optional + Tick labels. By default, the values of `lifetime`. + Only applies when `frequency` and `lifetime` are specified. + show_circle : bool, optional, default: True + Draw universal semicircle. + **kwargs + Additional parameters passed to + :py:class:`matplotlib.patches.Arc` and + :py:meth:`matplotlib.axes.Axes.plot`. + + """ + update_kwargs( + kwargs, + color=GRID_COLOR, + linestyle=GRID_LINESTYLE_MAJOR, + linewidth=GRID_LINEWIDH, + ) + if phasor_reference is not None: + polar_reference = phasor_to_polar_scalar(*phasor_reference) + if polar_reference is None: + polar_reference = (0.0, 1.0) + if phasor_reference is None: + phasor_reference = phasor_from_polar_scalar(*polar_reference) + ax = self._ax + + if show_circle: + ax.add_patch( + Arc( + (phasor_reference[0] / 2, phasor_reference[1] / 2), + polar_reference[1], + polar_reference[1], + theta1=math.degrees(polar_reference[0]), + theta2=math.degrees(polar_reference[0]) + 180.0, + fill=False, + **kwargs, + ) + ) + + if frequency is not None and polar_reference == (0.0, 1.0): + # draw ticks and labels + if lifetime is None: + lifetime = [0] + [ + 2**t + for t in range(-8, 32) + if phasor_from_lifetime(frequency, 2**t)[1] >= 0.18 + ] + unit = 'ns' + else: + unit = '' + if labels is None: + labels = [f'{tau:g}' for tau in lifetime] + try: + labels[2] = f'{labels[2]} {unit}' + except IndexError: + pass + ax.plot( + *phasor_calibrate( + *phasor_from_lifetime(frequency, lifetime), + *polar_reference, + ), + path_effects=[SemicircleTicks(labels=labels)], + **kwargs, + ) + + +class SemicircleTicks(AbstractPathEffect): + """Draw ticks on universal semicircle. + + Parameters + ---------- + size : float, optional + Length of tick in dots. + The default is ``rcParams['xtick.major.size']``. + labels : sequence of str, optional + Tick labels for each vertex in path. + **kwargs + Extra keywords passed to matplotlib's + :py:meth:`matplotlib.patheffects.AbstractPathEffect._update_gc`. + + """ + + _size: float # tick length + _labels: tuple[str, ...] # tick labels + _gc: dict[str, Any] # keywords passed to _update_gc + + def __init__( + self, + size: float | None = None, + labels: Sequence[str] | None = None, + **kwargs, + ) -> None: + super().__init__((0.0, 0.0)) + + if size is None: + self._size = pyplot.rcParams['xtick.major.size'] + else: + self._size = size + if labels is None or not labels: + self._labels = () + else: + self._labels = tuple(labels) + self._gc = kwargs + + def draw_path(self, renderer, gc, tpath, affine, rgbFace=None) -> None: + """Draw path with updated gc.""" + gc0 = renderer.new_gc() + gc0.copy_properties(gc) + + # TODO: this uses private methods of the base class + gc0 = self._update_gc(gc0, self._gc) # type: ignore + trans = affine + self._offset_transform(renderer) # type: ignore + + font = FontProperties() + # approximate half size of 'x' + fontsize = renderer.points_to_pixels(font.get_size_in_points()) / 4 + size = renderer.points_to_pixels(self._size) + origin = affine.transform([[0.5, 0.0]]) + + transpath = affine.transform_path(tpath) + polys = transpath.to_polygons(closed_only=False) + + for p in polys: + # coordinates of tick ends + t = p - origin + t /= numpy.hypot(t[:, 0], t[:, 1])[:, numpy.newaxis] + d = t.copy() + t *= size + t += p + + xyt = numpy.empty((2 * p.shape[0], 2)) + xyt[0::2] = p + xyt[1::2] = t + + renderer.draw_path( + gc0, + Path(xyt, numpy.tile([Path.MOVETO, Path.LINETO], p.shape[0])), + affine.inverted() + trans, + rgbFace, + ) + + if not self._labels: + continue + # coordinates of labels + t = d * size * 2.5 + t += p + + if renderer.flipy(): + h = renderer.get_canvas_width_height()[1] + else: + h = 0.0 + + for s, (x, y), (dx, _) in zip(self._labels, t, d): + # TODO: get rendered text size from matplotlib.text.Text? + # this did not work: + # Text(d[i,0], h - d[i,1], label, ha='center', va='center') + x = x + fontsize * len(s.split()[0]) * (dx - 1.0) + y = h - y + fontsize + renderer.draw_text(gc0, x, y, s, font, 0.0) + + gc0.restore() + + +def plot_phasor( + real: ArrayLike, + imag: ArrayLike, + /, + *, + style: Literal['plot', 'hist2d', 'contour'] | None = None, + allquadrants: bool | None = None, + frequency: float | None = None, + show: bool = True, + **kwargs: Any, +) -> None: + """Plot phasor coordinates. + + A simplified interface to the :py:class:`PhasorPlot` class. + + Parameters + ---------- + real : array_like + Real component of phasor coordinates. + imag : array_like + Imaginary component of phasor coordinates. + Must be of same shape as `real`. + style : {'plot', 'hist2d', 'contour'}, optional + Method used to plot phasor coordinates. + By default, if the number of coordinates are less than 65536 + and the arrays are less than three-dimensional, `'plot'` style is used, + else `'hist2d'`. + allquadrants : bool, optional + Show all quadrants of phasor space. + By default, only the first quadrant is shown. + frequency: float, optional + Frequency of phasor plot. + If provided, the universal circle is labeled with reference lifetimes. + show : bool, optional, default: True + Display figure. + **kwargs + Additional parguments passed to :py:class:`PhasorPlot`, + :py:meth:`PhasorPlot.plot`, :py:meth:`PhasorPlot.hist2d`, or + :py:meth:`PhasorPlot.contour` depending on `style`. + + """ + init_kwargs = parse_kwargs( + kwargs, + 'ax', + 'title', + 'xlabel', + 'ylabel', + 'xlim', + 'ylim', + 'xticks', + 'yticks', + 'grid', + ) + + real = numpy.asanyarray(real) + imag = numpy.asanyarray(imag) + plot = PhasorPlot( + frequency=frequency, allquadrants=allquadrants, **init_kwargs + ) + if style is None: + style = 'plot' if real.size < 65536 and real.ndim < 3 else 'hist2d' + if style == 'plot': + plot.plot(real, imag, **kwargs) + elif style == 'hist2d': + plot.hist2d(real, imag, **kwargs) + elif style == 'contour': + plot.contour(real, imag, **kwargs) + else: + raise ValueError(f'invalid {style=}') + if show: + plot.show() + + +def plot_phasor_image( + mean: ArrayLike | None, + real: ArrayLike, + imag: ArrayLike, + *, + harmonics: int | None = None, + percentile: float | None = None, + title: str | None = None, + show: bool = True, + **kwargs: Any, +) -> None: + """Plot phasor coordinates as images. + + Preview phasor coordinates from time-resolved or hyperspectral + image stacks as returned by :py:func:`phasorpy.phasor.phasor_from_signal`. + + The last two axes are assumed to be the image axes. + Harmonics, if any, are in the first axes of `real` and `imag`. + Other axes are averaged for display. + + Parameters + ---------- + mean : array_like + Image average. Must be two or more dimensional, or None. + real : array_like + Image of real component of phasor coordinates. + The last dimensions must match shape of `mean`. + imag : array_like + Image of imaginary component of phasor coordinates. + Must be same shape as `real`. + harmonics : int, optional + Number of harmonics to display. + If `mean` is None, a nonzero value indicates the presence of harmonics + in the first axes of `mean` and `real`. Else, the presence of harmonics + is determined from the shapes of `mean` and `real`. + By default, up to 4 harmonics are displayed. + percentile : float, optional + The (q, 100-q) percentiles of image data are covered by colormaps. + By default, the complete value range of `mean` is covered, + for `real` and `imag` the range [-1..1]. + title : str, optional + Figure title. + show : bool, optional, default: True + Display figure. + **kwargs + Additional arguments passed to :func:`matplotlib.pyplot.imshow`. + + Raises + ------ + ValueError + The shapes of `mean`, `real`, and `image` do not match. + Percentile is out of range. + + """ + update_kwargs(kwargs, interpolation='nearest') + cmap = kwargs.pop('cmap', None) + shape = None + + if mean is not None: + mean = numpy.asarray(mean) + if mean.ndim < 2: + raise ValueError(f'not an image {mean.ndim=} < 2') + shape = mean.shape + mean = numpy.mean(mean.reshape(-1, *mean.shape[-2:]), axis=0) + + real = numpy.asarray(real) + imag = numpy.asarray(imag) + if real.shape != imag.shape: + raise ValueError(f'{real.shape=} != {imag.shape=}') + if real.ndim < 2: + raise ValueError(f'not an image {real.ndim=} < 2') + + if (shape is not None and real.shape[1:] == shape) or ( + shape is None and harmonics + ): + # first image dimension contains harmonics + if real.ndim < 3: + raise ValueError(f'not a multi-harmonic image {real.shape=}') + nh = real.shape[0] # number harmonics + elif shape is None or shape == real.shape: + # single harmonic + nh = 1 + else: + raise ValueError(f'shape mismatch {real.shape[1:]=} != {shape}') + + # average extra image dimensions, but not harmonics + real = numpy.mean(real.reshape(nh, -1, *real.shape[-2:]), axis=1) + imag = numpy.mean(imag.reshape(nh, -1, *imag.shape[-2:]), axis=1) + + # for MyPy + assert isinstance(mean, numpy.ndarray) or mean is None + assert isinstance(real, numpy.ndarray) + assert isinstance(imag, numpy.ndarray) + + # limit number of displayed harmonics + nh = min(4 if harmonics is None else harmonics, nh) + + # create figure with size depending on image aspect and number of harmonics + fig = pyplot.figure(layout='constrained') + w, h = fig.get_size_inches() + aspect = min(1.0, max(0.5, real.shape[-2] / real.shape[-1])) + fig.set_size_inches(w, h * 0.4 * aspect * nh + h * 0.25 * aspect) + gs = GridSpec(nh, 2 if mean is None else 3, figure=fig) + if title: + fig.suptitle(title) + + if mean is not None: + _imshow( + fig.add_subplot(gs[0, 0]), + mean, + percentile=percentile, + vmin=None, + vmax=None, + cmap=cmap, + axis=True, + title='mean', + **kwargs, + ) + + if percentile is None: + vmin = -1.0 + vmax = 1.0 + if cmap is None: + cmap = 'coolwarm_r' + else: + vmin = None + vmax = None + + for h in range(nh): + axs = [] + ax = fig.add_subplot(gs[h, -2]) + axs.append(ax) + _imshow( + ax, + real[h], + percentile=percentile, + vmin=vmin, + vmax=vmax, + cmap=cmap, + axis=mean is None and h == 0, + colorbar=percentile is not None, + title=None if h else 'G, real', + **kwargs, + ) + + ax = fig.add_subplot(gs[h, -1]) + axs.append(ax) + pos = _imshow( + ax, + imag[h], + percentile=percentile, + vmin=vmin, + vmax=vmax, + cmap=cmap, + axis=False, + colorbar=percentile is not None, + title=None if h else 'S, imag', + **kwargs, + ) + if percentile is None and h == 0: + fig.colorbar(pos, ax=axs, shrink=0.4, location='bottom') + + if show: + pyplot.show() + + +def plot_signal_image( + signal: ArrayLike, + /, + *, + axis: int | None = None, + percentile: float | Sequence[float] | None = None, + title: str | None = None, + show: bool = True, + **kwargs: Any, +) -> None: + """Plot average image and signal along axis. + + Preview time-resolved or hyperspectral image stacks to be anayzed with + :py:func:`phasorpy.phasor.phasor_from_signal`. + + The last two axes, excluding `axis`, are assumed to be the image axes. + Other axes are averaged for image display. + + Parameters + ---------- + signal : array_like + Image stack. Must be three or more dimensional. + axis : int, optional, default: -1 + Axis over which phasor coordinates would be computed. + The default is the last axis (-1). + percentile : float or [float, float], optional + The [q, 100-q] percentiles of image data are covered by colormaps. + By default, the complete value range of `mean` is covered, + for `real` and `imag` the range [-1..1]. + title : str, optional + Figure title. + show : bool, optional, default: True + Display figure. + **kwargs + Additional arguments passed to :func:`matplotlib.pyplot.imshow`. + + Raises + ------ + ValueError + Signal is not an image stack. + Percentile is out of range. + + """ + # TODO: add option to separate channels? + # TODO: add option to plot non-images? + update_kwargs(kwargs, interpolation='nearest') + signal = numpy.asarray(signal) + if signal.ndim < 3: + raise ValueError(f'not an image stack {signal.ndim=} < 3') + + if axis is None: + axis = -1 + axis %= signal.ndim + + # for MyPy + assert isinstance(signal, numpy.ndarray) + + fig = pyplot.figure(layout='constrained') + if title: + fig.suptitle(title) + w, h = fig.get_size_inches() + fig.set_size_inches(w, h * 0.7) + gs = GridSpec(1, 2, figure=fig, width_ratios=(1, 1)) + + # histogram + axes = list(range(signal.ndim)) + del axes[axis] + ax = fig.add_subplot(gs[0, 1]) + ax.set_title(f'mean, axis {axis}') + ax.plot(signal.mean(axis=tuple(axes))) + + # image + axes = list(sorted(axes[:-2] + [axis])) + ax = fig.add_subplot(gs[0, 0]) + _imshow( + ax, + signal.mean(axis=tuple(axes)), + percentile=percentile, + shrink=0.5, + title='mean', + ) + + if show: + pyplot.show() + + +def plot_polar_frequency( + frequency: ArrayLike, + phase: ArrayLike, + modulation: ArrayLike, + *, + ax: Axes | None = None, + title: str | None = None, + show: bool = True, + **kwargs, +) -> None: + """Plot phase and modulation verus frequency. + + Parameters + ---------- + frequency : array_like, shape (n, ) + Laser pulse or modulation frequency in MHz. + phase : array_like + Angular component of polar coordinates in radians. + modulation : array_like + Radial component of polar coordinates. + ax : matplotlib axes, optional + Matplotlib axes used for plotting. + By default, a new subplot axes is created. + title : str, optional + Figure title. The default is "Multi-frequency plot". + show : bool, optional, default: True + Display figure. + **kwargs + Additional arguments passed to :py:func:`matplotlib.pyplot.plot`. + + """ + # TODO: make this customizable: labels, colors, ... + if ax is None: + ax = pyplot.subplots()[1] + if title is None: + title = 'Multi-frequency plot' + if title: + ax.set_title(title) + ax.set_xscale('log', base=10) + ax.set_xlabel('frequency (MHz)') + + phase = numpy.asarray(phase) + if phase.ndim < 2: + phase = phase.reshape(-1, 1) + modulation = numpy.asarray(modulation) + if modulation.ndim < 2: + modulation = modulation.reshape(-1, 1) + + ax.set_ylabel('phase (°)', color='tab:blue') + ax.set_yticks([0.0, 30.0, 60.0, 90.0]) + for phi in phase.T: + ax.plot(frequency, numpy.rad2deg(phi), color='tab:blue', **kwargs) + ax = ax.twinx() # type: ignore + + ax.set_ylabel('modulation (%)', color='tab:red') + ax.set_yticks([0.0, 25.0, 50.0, 75.0, 100.0]) + for mod in modulation.T: + ax.plot(frequency, mod * 100, color='tab:red', **kwargs) + if show: + pyplot.show() + + +def _imshow( + ax: Axes, + image: NDArray[Any], + /, + *, + percentile: float | Sequence[float] | None = None, + vmin: float | None = None, + vmax: float | None = None, + colorbar: bool = True, + shrink: float | None = None, + axis: bool = True, + title: str | None = None, + **kwargs, +) -> AxesImage: + """Plot image array. + + Convenience wrapper around :py:func:`matplotlib.pyplot.imshow`. + + """ + update_kwargs(kwargs, interpolation='none') + if percentile is not None: + if isinstance(percentile, Sequence): + percentile = percentile[0], percentile[1] + else: + # percentile = max(0.0, min(50, percentile)) + percentile = percentile, 100.0 - percentile + if ( + percentile[0] >= percentile[1] + or percentile[0] < 0 + or percentile[1] > 100 + ): + raise ValueError(f'{percentile=} out of range') + vmin, vmax = numpy.percentile(image, percentile) + pos = ax.imshow(image, vmin=vmin, vmax=vmax, **kwargs) + if colorbar: + if percentile is not None and vmin is not None and vmax is not None: + ticks = vmin, vmax + else: + ticks = None + fig = ax.get_figure() + if fig is not None: + if shrink is None: + shrink = 0.8 + fig.colorbar(pos, shrink=shrink, location='bottom', ticks=ticks) + if title: + ax.set_title(title) + if not axis: + ax.set_axis_off() + # ax.set_anchor('C') + return pos diff --git a/tests/test__utils.py b/tests/test__utils.py new file mode 100644 index 00000000..0287eccf --- /dev/null +++ b/tests/test__utils.py @@ -0,0 +1,134 @@ +"""Test the phasorpy._utils module.""" + +import math + +import pytest +from numpy.testing import assert_allclose + +from phasorpy._utils import ( + circle_circle_intersection, + circle_line_intersection, + kwargs_notnone, + parse_kwargs, + phasor_from_polar_scalar, + phasor_to_polar_scalar, + scale_matrix, + sort_coordinates, + update_kwargs, +) + + +def test_phasor_to_polar_scalar(): + """Test phasor_to_polar_scalar function.""" + assert phasor_to_polar_scalar(0.0, 0.0) == (0.0, 0.0) + assert_allclose( + phasor_to_polar_scalar(0.8, 0.4, degree=True, percent=True), + (26.565051, 89.442719), + atol=1e-6, + ) + assert_allclose( + phasor_to_polar_scalar(1.0, 0.0, degree=True, percent=True), + (0.0, 100.0), + atol=1e-6, + ) + + +def test_phasor_from_polar_scalar(): + """Test phasor_from_polar_scalar function.""" + assert phasor_from_polar_scalar(0.0, 0.0) == (0.0, 0.0) + assert_allclose( + phasor_from_polar_scalar( + 26.565051, 89.442719, degree=True, percent=True + ), + (0.8, 0.4), + atol=1e-6, + ) + assert_allclose( + phasor_from_polar_scalar(0.0, 100.0, degree=True, percent=True), + (1.0, 0.0), + atol=1e-6, + ) + # roundtrip + assert_allclose( + phasor_from_polar_scalar( + *phasor_to_polar_scalar(-0.4, -0.2, degree=True, percent=True), + degree=True, + percent=True, + ), + (-0.4, -0.2), + atol=1e-6, + ) + + +def test_parse_kwargs(): + """Test parse_kwargs function.""" + kwargs = {'one': 1, 'two': 2, 'four': 4} + kwargs2 = parse_kwargs(kwargs, 'two', 'three', four=None, five=5) + assert kwargs == {'one': 1} + assert kwargs2 == {'two': 2, 'four': 4, 'five': 5} + + kwargs = {'one': 1, 'two': 2, 'four': 4} + kwargs2 = parse_kwargs( + kwargs, 'two', 'three', four=None, five=5, _del=False + ) + assert kwargs == {'one': 1, 'two': 2, 'four': 4} + assert kwargs2 == {'two': 2, 'four': 4, 'five': 5} + + +def test_update_kwargs(): + """Test update_kwargs function.""" + kwargs = { + 'one': 1, + } + update_kwargs(kwargs, one=None, two=2) + assert kwargs == {'one': 1, 'two': 2} + + +def test_kwargs_notnone(): + """Test kwargs_notnone function.""" + assert kwargs_notnone(one=1, none=None) == {'one': 1} + + +def test_scale_matrix(): + """Test scale_matrix function.""" + assert_allclose( + scale_matrix(1.1, (0.0, 0.5)), + [[1.1, 0, -0], [0, 1.1, -0.05], [0, 0, 1]], + 1e-6, + ) + + +def test_sort_coordinates(): + """Test sort_coordinates function.""" + x, y = sort_coordinates([0, 1, 2, 3], [0, 1, -1, 0]) + assert_allclose(x, [2, 3, 1, 0]) + assert_allclose(y, [-1, 0, 1, 0]) + + x, y = sort_coordinates([0, 1, 2], [0, 1, -1]) + assert_allclose(x, [0, 1, 2]) + assert_allclose(y, [0, 1, -1]) + + with pytest.raises(ValueError): + sort_coordinates([0, 1, 2, 3], [0, 1, -1]) + + +def test_circle_line_intersection(): + """Test circle_line_intersection function.""" + assert_allclose( + circle_circle_intersection( + 0.0, 0.0, math.hypot(0.6, 0.4), 0.6, 0.4, 0.2 + ), + ((0.6868, 0.2198), (0.4670, 0.5494)), + 1e-3, + ) + assert not circle_circle_intersection(0.0, 0.0, 1.0, 0.6, 0.4, 0.2) + + +def test_circle_circle_intersection(): + """Test circle_circle_intersection function.""" + assert_allclose( + circle_line_intersection(0.6, 0.4, 0.2, 0.0, 0.0, 0.6, 0.4), + ((0.7664, 0.5109), (0.4335, 0.2890)), + 1e-3, + ) + assert not circle_line_intersection(0.6, 0.4, 0.2, 0.0, 0.0, 0.6, 0.1) diff --git a/tests/test_plot.py b/tests/test_plot.py new file mode 100644 index 00000000..928e6f44 --- /dev/null +++ b/tests/test_plot.py @@ -0,0 +1,373 @@ +"""Test the phasorpy.plot module.""" + +import io +import math + +import numpy +import pytest +from matplotlib import pyplot + +from phasorpy.plot import ( + PhasorPlot, + plot_phasor, + plot_phasor_image, + plot_polar_frequency, + plot_signal_image, +) + +INTERACTIVE = False # enable for interactive plotting + + +class TestPhasorPlot: + """Test PhasorPlot class.""" + + def show(self, plot): + """Show plot.""" + if INTERACTIVE: + plot.show() + pyplot.close() + + def test_init(self): + """Test __init__ and attributes.""" + plot = PhasorPlot(title='default') + self.show(plot) + + plot = PhasorPlot(frequency=80.0, title='frequency') + self.show(plot) + + plot = PhasorPlot(grid=False, title='no grid') + self.show(plot) + + plot = PhasorPlot(allquadrants=True, title='allquadrants') + self.show(plot) + + plot = PhasorPlot(title='kwargs', xlim=(-0.1, 1.1), ylim=(-0.1, 0.9)) + self.show(plot) + + fig, ax = pyplot.subplots() + plot = PhasorPlot(ax=ax, title='axes') + assert plot.ax == ax + assert plot.fig == fig + self.show(plot) + + def test_save(self): + """Test save method.""" + fh = io.BytesIO() + plot = PhasorPlot(title='save') + plot.save(fh, format='png') + assert fh.getvalue()[:6] == b'\x89PNG\r\n' + pyplot.close() + + def test_plot(self): + """Test plot method.""" + plot = PhasorPlot(title='plot') + plot.plot(0.6, 0.4, label='1') + plot.plot([0.2, 0.9], [0.4, 0.3], '.-', label='2') + plot.plot( + [[0.29, 0.3, 0.31], [0.41, 0.29, 0.3]], + [[0.31, 0.29, 0.2], [0.49, 0.5, 0.51]], + 'x', + label='3', + ) + self.show(plot) + + def test_hist2d(self): + """Test hist2d method.""" + real, imag = numpy.random.multivariate_normal( + (0.6, 0.4), [[3e-3, -1e-3], [-1e-3, 1e-3]], (256, 256) + ).T + plot = PhasorPlot(title='hist2d') + plot.hist2d(real, imag) + self.show(plot) + + plot = PhasorPlot(title='hist2d parameters', allquadrants=True) + plot.hist2d( + real, imag, bins=100, cmax=500, cmap='viridis', norm='linear' + ) + self.show(plot) + + def test_contour(self): + """Test contour method.""" + real, imag = numpy.random.multivariate_normal( + (0.6, 0.4), [[3e-3, -1e-3], [-1e-3, 1e-3]], (256, 256) + ).T + plot = PhasorPlot(title='contour') + plot.contour(real, imag) + self.show(plot) + + plot = PhasorPlot(title='contour parameters', allquadrants=True) + plot.contour(real, imag, bins=200, cmap='viridis', norm='linear') + self.show(plot) + + def test_imshow(self): + """Test imshow method.""" + plot = PhasorPlot(title='imshow') + with pytest.raises(NotImplementedError): + plot.imshow([[0]]) + self.show(plot) + + @pytest.mark.parametrize('allquadrants', (True, False)) + def test_components(self, allquadrants): + """Test components method.""" + real = [0.1, 0.2, 0.5, 0.9] + imag = [0.3, 0.4, 0.5, 0.3] + weights = [2, 1, 2, 1] + plot = PhasorPlot(title='components', allquadrants=allquadrants) + with pytest.raises(ValueError): + plot.components(0.5, 0.5) + plot.components(real, imag, fill=True, facecolor='lightyellow') + plot.components(real, imag, weights, linestyle='-', color='tab:blue') + self.show(plot) + + def test_line(self): + """Test line method.""" + plot = PhasorPlot(title='line') + plot.line([0.8, 0.4], [0.2, 0.3], color='tab:red', linestyle='--') + self.show(plot) + + def test_circle(self): + """Test circle method.""" + plot = PhasorPlot(title='circle') + plot.circle(0.5, 0.2, 0.1, color='tab:red', linestyle='-') + self.show(plot) + + @pytest.mark.parametrize('allquadrants', (True, False)) + def test_polar_cursor(self, allquadrants): + """Test polar_cursor method.""" + plot = PhasorPlot(title='polar_cursor', allquadrants=allquadrants) + plot.polar_cursor() + plot.polar_cursor(0.6435, 0.5, color='tab:blue', linestyle='-') + plot.polar_cursor(0.5236, 0.6, 0.1963, 0.8, color='tab:orange') + plot.polar_cursor(0.3233, 0.9482, radius=0.05, color='tab:green') + plot.polar_cursor(0.3, color='tab:red', linestyle='--') + self.show(plot) + + def test_polar_grid(self): + """Test polar_grid method.""" + plot = PhasorPlot(grid=False, allquadrants=True, title='polar_grid') + plot.polar_grid(color='tab:red', linestyle='-') + self.show(plot) + + def test_semicircle(self): + """Test semicircle method.""" + plot = PhasorPlot(grid=False, title='empty') + plot.semicircle() + self.show(plot) + + plot = PhasorPlot(grid=False, title='frequency') + plot.semicircle(frequency=80) + self.show(plot) + + plot = PhasorPlot(grid=False, title='no labels') + plot.semicircle(frequency=80, labels=()) + self.show(plot) + + plot = PhasorPlot(grid=False, title='no circle') + plot.semicircle(frequency=80, show_circle=False) + self.show(plot) + + plot = PhasorPlot(grid=False, title='red') + plot.semicircle(frequency=80, color='tab:red', linestyle=':') + self.show(plot) + + plot = PhasorPlot(grid=False, title='lifetime') + plot.semicircle(frequency=80, lifetime=[1, 2]) + self.show(plot) + + plot = PhasorPlot(grid=False, title='labels') + plot.semicircle( + frequency=80, lifetime=[1, 2], labels=['label 1', 'label 2'] + ) + self.show(plot) + + plot = PhasorPlot(title='polar_reference', xlim=(-0.2, 1.05)) + plot.semicircle(polar_reference=(0.9852, 0.5526)) + self.show(plot) + + plot = PhasorPlot( + frequency=80.0, title='phasor_reference', xlim=(-0.2, 1.05) + ) + plot.semicircle(frequency=80.0, phasor_reference=(0.2, 0.4)) + self.show(plot) + + +def test_plot_phasor(): + """Test plot_phasor function.""" + real, imag = numpy.random.multivariate_normal( + (0.6, 0.4), [[3e-3, -1e-3], [-1e-3, 1e-3]], 32 + ).T + plot_phasor( + real, + imag, + style='plot', + title='plot', + color='tab:red', + frequency=80.0, + show=INTERACTIVE, + ) + pyplot.close() + + _, ax = pyplot.subplots() + real, imag = numpy.random.multivariate_normal( + (0.6, 0.4), [[3e-3, -1e-3], [-1e-3, 1e-3]], (256, 256) + ).T + plot_phasor( + real, + imag, + ax=ax, + title='hist2d', + cmap='Blues', + allquadrants=True, + grid=False, + show=INTERACTIVE, + ) + pyplot.close() + + plot_phasor( + real, + imag, + style='contour', + title='contour', + cmap='viridis', + levels=4, + show=INTERACTIVE, + ) + pyplot.close() + + with pytest.raises(ValueError): + plot_phasor(0, 0, style='invalid') + + +def test_plot_polar_frequency(): + """Test plot_polar_frequency function.""" + plot_polar_frequency( + [1, 10, 100], + [0, 0.5, 1], + [1, 0.5, 0], + title='plot_polar_frequency', + show=INTERACTIVE, + ) + pyplot.close() + + _, ax = pyplot.subplots() + plot_polar_frequency( + [1, 10, 100], + [[0, 0.1], [0.5, 0.55], [1, 1]], + [[[1, 0.9], [0.5, 0.45], [0, 0]]], + ax=ax, + show=INTERACTIVE, + ) + pyplot.close() + + +def test_plot_signal_image(): + """Test plot_signal_image function.""" + shape = (7, 31, 33, 11) + data = numpy.arange(math.prod(shape)).reshape(shape) + data %= math.prod(shape[-2:]) + data = data / math.prod(shape[-2:]) + + plot_signal_image(data, title='default', show=INTERACTIVE) + pyplot.close() + plot_signal_image(data, axis=0, title='axis 0', show=INTERACTIVE) + pyplot.close() + plot_signal_image(data, axis=2, title='axis 2', show=INTERACTIVE) + pyplot.close() + plot_signal_image( + data, + cmap='hot', + percentile=[5, 95], + title='percentile', + show=INTERACTIVE, + ) + pyplot.close() + + with pytest.raises(ValueError): + # not an image + plot_signal_image(data[0, 0], show=False) + pyplot.close() + with pytest.raises(ValueError): + # percentile out of range + plot_signal_image(data, percentile=(-1, 101), show=False) + pyplot.close() + + +def test_plot_phasor_image(): + """Test plot_phasor_image function.""" + shape = (7, 11, 31, 33) + data = numpy.arange(math.prod(shape)).reshape(shape) + data %= math.prod(shape[-2:]) + data = data / math.prod(shape[-2:]) + + # 2D data + d = data[0, 0] + plot_phasor_image(d, d, d, title='mean, real, imag', show=INTERACTIVE) + pyplot.close() + plot_phasor_image(None, d, d, title='real, imag', show=INTERACTIVE) + pyplot.close() + # 4D data + d = data + plot_phasor_image(d, d, d, title='4D images', show=INTERACTIVE) + pyplot.close() + # 7 harmonics + plot_phasor_image(d[0], d, d, title='harmonics up to 4', show=INTERACTIVE) + pyplot.close() + plot_phasor_image( + None, + d, + d, + harmonics=2, + title='real and imag harmonics up to 2', + show=INTERACTIVE, + ) + pyplot.close() + + d = data[0, 0] + plot_phasor_image( + d, + d, + d, + percentile=5.0, + cmap='hot', + title='5th percentile with colormap', + show=INTERACTIVE, + ) + pyplot.close() + + d = data[0, 0, 0] + with pytest.raises(ValueError): + # not an image + plot_phasor_image(d, d, d, show=False) + pyplot.close() + with pytest.raises(ValueError): + # not an image + plot_phasor_image(None, d, d, show=False) + pyplot.close() + + d = data[0, 0] + with pytest.raises(ValueError): + # not an image + plot_phasor_image(None, d, d, harmonics=2, show=False) + pyplot.close() + + d = data[0] + with pytest.raises(ValueError): + # shape mismatch + plot_phasor_image(d, d[0], d, show=False) + pyplot.close() + with pytest.raises(ValueError): + # shape mismatch + plot_phasor_image(d, d, d[0], show=False) + pyplot.close() + with pytest.raises(ValueError): + # shape mismatch + plot_phasor_image(d, d[0, :-1], d[0, :-1], show=False) + pyplot.close() + with pytest.raises(ValueError): + # percentile out of range + plot_phasor_image(d, d, d, percentile=-1, show=False) + pyplot.close() + with pytest.raises(ValueError): + # percentile out of range + plot_phasor_image(d, d, d, percentile=50, show=False) + pyplot.close() diff --git a/tutorials/phasorpy_introduction.py b/tutorials/phasorpy_introduction.py index c512e93d..ab5b17f7 100644 --- a/tutorials/phasorpy_introduction.py +++ b/tutorials/phasorpy_introduction.py @@ -1,8 +1,11 @@ """ -Introduction -============ +Introduction to PhasorPy +======================== -This tutorial provides an introduction to using the PhasorPy library. +An introduction to using the PhasorPy library. + +PhasorPy is an open-source Python library for the analysis of fluorescence +lifetime and hyperspectral images using the :doc:`/phasor_approach`. """ @@ -22,53 +25,189 @@ # Install PhasorPy # ---------------- # -# On a command prompt, shell or terminal, run ``pip`` to download and install -# the PhasorPy library and all its dependencies from the Python Package Index:: +# To download and install the PhasorPy library and all its dependencies from +# the `Python Package Index `_ (PyPI), +# run the following command on a command prompt, shell or terminal:: +# +# python -m pip install -U "phasorpy[all]" # -# $ python -m pip install -U phasorpy +# .. note:: +# The PhasorPy library is in its early stages of development +# and has not yet been released to PyPI. +# The development version of PhasorPy can be `installed manually +# `_, +# for example, using the binary wheels from `GitHub Actions +# `_, +# or the source code on GitHub (requires a C compiler):: +# +# python -m pip install git+https://github.com/phasorpy/phasorpy.git # %% # Import phasorpy # --------------- # -# Open the Python interpreter, import the ``phasorpy`` package, and check its -# version: +# Start the Python interpreter, import the ``phasorpy`` package, +# and print its version: import phasorpy print(phasorpy.__version__) # %% -# Read data file -# -------------- +# Read signal from file +# --------------------- # -# PhasorPy provides many functions to read image and metadata from file formats -# used in microscopy. -# For example, ...: +# The :py:mod:`phasorpy.datasets` module provides access to various sample +# files, for example, a TIFF file containing a time-correlated +# single photon counting (TCSPC) histogram obtained at 80 MHz. +# +# The :py:mod:`phasorpy.io` module provides many functions to read +# time-resolved and hyperspectral image and metadata from file formats used +# in microscopy. However, here the +# `tifffile `_ library is used directly: + +# TODO: use phasorpy.io function to read histogram and metadata from PTU file + +import tifffile + +from phasorpy.datasets import fetch -... +signal = tifffile.imread(fetch('Embryo.tif')) +frequency = 80.0 + +print(signal.shape, signal.dtype) + +# %% +# Plot the spatial and histogram averages: + +from phasorpy.plot import plot_signal_image + +plot_signal_image(signal, axis=0) # %% # Calculate phasor coordinates # ---------------------------- # -# Calculate the phasor coordinates for the image stack using known reference -# coordinates: +# The :py:mod:`phasorpy.phasor` module provides functions to calculate, +# convert, and correct phasor coordinates. +# +# Phasor coordinate are calculated from the signal, a TCSPC histogram in +# this case. The histogram samples are in the first dimension (`axis=0`): + +from phasorpy.phasor import phasor_from_signal + +mean, real, imag = phasor_from_signal(signal, axis=0) + +# %% +# Plot the calculated phasor coordinates: + +from phasorpy.plot import plot_phasor_image + +plot_phasor_image(mean, real, imag) + +# %% +# Calibrate phasor coordinates +# ---------------------------- +# +# Phasor coordinates from time-resolved measurements must be calibrated +# with the polar coordinates (phase and modulation) obtained from a +# reference standard of known lifetime, acquired with the same instrument +# settings. +# +# Read the signal of the reference measurement from a file: + +reference_signal = tifffile.imread(fetch('Fluorescein_Embryo.tif')) + +# %% +# Calculate phasor coordinates from the measured reference signal: + +_, measured_real, measured_imag = phasor_from_signal(reference_signal, axis=0) + +# %% +# Average the phasor coordinates, assuming there are no spatial aberations: + +from phasorpy.phasor import phasor_center + +measured_real, measured_imag = phasor_center(measured_real, measured_imag) + +# %% +# Calculate absolute phasor coordinates from the known lifetime of the +# reference (Fluorescein, 4.2 ns): + +from phasorpy.phasor import phasor_from_lifetime -... +known_real, known_imag = phasor_from_lifetime(frequency, 4.2) # %% -# Display phasor plot -# ------------------- +# Calculate polar coordinates (phase shift and relative modulation) to +# correct phasor coordinates: + +from phasorpy.phasor import polar_from_reference_phasor + +phase0, modulation0 = polar_from_reference_phasor( + measured_real, measured_imag, known_real, known_imag +) + +# %% +# Finally, calibrate the raw phasor coordinates obtained from the signal: + +from phasorpy.phasor import phasor_calibrate + +real, imag = phasor_calibrate(real, imag, phase0, modulation0) + +# %% +# Filter phasor coordinates +# ------------------------- +# +# Applying median filter to the calibrated phasor coordinates, +# often multiple times, improves contrast: -from matplotlib import pyplot +# TODO: replace this with a ``phasor_filter`` function? +from skimage.filters import median -pyplot.imshow([[0, 1], [2, 3]]) +for _ in range(2): + real = median(real) + imag = median(imag) # %% -# Segment intensity image +# Pixels with low intensities are commonly excluded from analysis and +# visualization of phasor coordinates: + +# TODO: replace this with a ``phasor_mask`` function? +mask = mean > 1 +real = real[mask] +imag = imag[mask] + +# %% +# Plot phasor coordinates # ----------------------- -# Show the intensity image, colored according to a selection of phasor -# coordinates: +# +# The :py:mod:`phasorpy.plot` module provides functions and classes for +# plotting phasor and polar coordinates. +# +# Large number of phasor coordinates, such as obtained from imaging, +# are commonly visualized as 2D histograms: + +from phasorpy.plot import PhasorPlot -... +plot = phasorpy.plot.PhasorPlot( + frequency=frequency, title='Calibrated, filtered phasor coordinates' +) +plot.hist2d(real, imag) +plot.show() + +# %% +# For comparison, the uncalibrated, unfiltered phasor coordinates: + +plot = PhasorPlot(allquadrants=True, title='Raw phasor coordinates') +plot.hist2d(*phasor_from_signal(signal, axis=0)[1:]) +plot.semicircle() +plot.show() + +# %% +# To be continued +# --------------- +# + +# %% +# sphinx_gallery_thumbnail_number = 3 diff --git a/tutorials/phasorpy_phasor_from_lifetime.py b/tutorials/phasorpy_phasor_from_lifetime.py index adeb2643..6639faf2 100644 --- a/tutorials/phasorpy_phasor_from_lifetime.py +++ b/tutorials/phasorpy_phasor_from_lifetime.py @@ -2,6 +2,8 @@ Phasor coordinates from lifetimes ================================= +An introduction to the `phasor_from_lifetime` function. + The :py:func:`phasorpy.phasor.phasor_from_lifetime` function is used to calculate phasor coordinates as a function of frequency, single or multiple lifetime components, and the pre-exponential amplitudes @@ -10,141 +12,12 @@ """ # %% -# Import required modules and functions, and define helper functions for -# plotting phasor or polar coordinates: - -import math +# Import required modules and functions: import numpy -from phasorpy.phasor import ( - phasor_from_lifetime, - phasor_semicircle, - phasor_to_polar, -) - - -def phasor_plot( - real, - imag, - fmt: str = 'o', - *, - ax=None, - mode: str = 'lifetime', - style: str | None = None, - title: str | None = None, - xlabel: str | None = None, - ylabel: str | None = None, - bins: int | None = None, - cmap=None, - show: bool = True, - return_ax: bool = False, -): - """Plot phasor coordinates using matplotlib.""" - # TODO: move this function to phasorpy.plot - from matplotlib import pyplot - from matplotlib.lines import Line2D - - if mode == 'lifetime': - xlim = [-0.05, 1.05] - ylim = [-0.05, 0.65] - ranges = [[0, 1], [0, 0.625]] - bins = 256 if bins is None else bins - bins_list = [bins, int(bins * 0.625)] - elif mode == 'spectral': - xlim = [-1.05, 1.05] - ylim = [-1.05, 1.05] - ranges = [[0, 1], [0, 1]] - bins = 256 if bins is None else bins - bins_list = [bins, bins] - else: - raise ValueError(f'unknown {mode=!r}') - if ax is None: - ax = pyplot.subplots()[1] - if style is None: - style = 'scatter' if real.size < 1024 else 'histogram' - if real is None or imag is None: - pass - elif style == 'histogram': - ax.hist2d( - real, - imag, - range=ranges, - bins=bins_list, - cmap='Blues' if cmap is None else cmap, - norm='log', - ) - elif style == 'scatter': - for re, im in zip( - numpy.array(real, ndmin=2), numpy.array(imag, ndmin=2) - ): - ax.plot(re, im, fmt) - if mode == 'lifetime': - ax.plot(*phasor_semicircle(100), color='k', lw=0.5) - ax.set_xticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0]) - ax.set_yticks([0.0, 0.2, 0.4, 0.6]) - elif mode == 'spectral': - ax.add_patch( - pyplot.Circle((0, 0), 1, color='k', lw=0.5, ls='--', fill=False) - ) - ax.add_patch( - pyplot.Circle( - (0, 0), 2 / 3, color='0.5', lw=0.25, ls='--', fill=False - ) - ) - ax.add_patch( - pyplot.Circle( - (0, 0), 1 / 3, color='0.5', lw=0.25, ls='--', fill=False - ) - ) - ax.add_line(Line2D([-1, 1], [0, 0], color='k', lw=0.5, ls='--')) - ax.add_line(Line2D([0, 0], [-1, 1], color='k', lw=0.5, ls='--')) - for a in (3, 6): - x = math.cos(math.pi / a) - y = math.sin(math.pi / a) - ax.add_line( - Line2D([-x, x], [-y, y], color='0.5', lw=0.25, ls='--') - ) - ax.add_line( - Line2D([-x, x], [y, -y], color='0.5', lw=0.25, ls='--') - ) - ax.set_xticks([-1.0, -0.5, 0.0, 0.5, 1.0]) - ax.set_yticks([-1.0, -0.5, 0.0, 0.5, 1.0]) - - ax.set( - title='Phasor plot' if title is None else title, - xlabel='G, real' if xlabel is None else xlabel, - ylabel='S, imag' if ylabel is None else ylabel, - aspect='equal', - xlim=xlim, - ylim=ylim, - ) - if show: - pyplot.show() - if return_ax: - return ax - - -def multi_frequency_plot(frequency, phase, modulation, title=None): - """Plot phase and modulation vs frequency.""" - # TODO: move this function to phasorpy.plot - from matplotlib import pyplot - - ax = pyplot.subplots()[1] - ax.set_title('Multi-frequency plot' if title is None else title) - ax.set_xscale('log', base=10) - ax.set_xlabel('frequency (MHz)') - ax.set_ylabel('phase (°)', color='tab:blue') - ax.set_yticks([0.0, 30.0, 60.0, 90.0]) - for phi in numpy.array(phase, ndmin=2).swapaxes(0, 1): - ax.plot(frequency, numpy.rad2deg(phi), color='tab:blue') - ax = ax.twinx() - ax.set_ylabel('modulation (%)', color='tab:red') - ax.set_yticks([0.0, 25.0, 50.0, 75.0, 100.0]) - for mod in numpy.array(modulation, ndmin=2).swapaxes(0, 1): - ax.plot(frequency, mod * 100, color='tab:red') - pyplot.show() - +from phasorpy.phasor import phasor_from_lifetime, phasor_to_polar +from phasorpy.plot import PhasorPlot, plot_phasor, plot_polar_frequency # %% # Single-component lifetimes @@ -156,7 +29,7 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): lifetime = numpy.array([3.9788735, 0.9947183]) -phasor_plot(*phasor_from_lifetime(80.0, lifetime)) +plot_phasor(*phasor_from_lifetime(80.0, lifetime), frequency=80.0) # %% # Multi-component lifetimes @@ -170,7 +43,9 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): [[1, 0], [0.25, 0.75], [0.5, 0.5], [0.75, 0.25], [0, 1]] ) -phasor_plot(*phasor_from_lifetime(80.0, lifetime, fraction), fmt='o-') +plot_phasor( + *phasor_from_lifetime(80.0, lifetime, fraction), fmt='o-', frequency=80.0 +) # %% # Pre-exponential amplitudes @@ -179,9 +54,10 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): # The phasor coordinates of two lifetime components with varying # pre-exponential amplitudes are also located on a line: -phasor_plot( +plot_phasor( *phasor_from_lifetime(80.0, lifetime, fraction, preexponential=True), fmt='o-', + frequency=80.0, ) # %% @@ -189,12 +65,14 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): # ---------------------------------------------- # # Phasor coordinates can be calculated at once for many frequencies, -# lifetime components, and their fractions. -# As an example, lifetimes are passed in units of s and frequencies in Hz, -# requiring to specify a unit_conversion factor: +# lifetime components, and their fractions. As an example, random distrinutions +# of lifetimes and their fractions are plotted at three frequencies. +# Lifetimes are passed in units of s and frequencies in Hz, requiring to +# specify a `unit_conversion` factor: -samples = 100 rng = numpy.random.default_rng() + +samples = 100 lifetime_distribution = ( numpy.column_stack( ( @@ -209,7 +87,7 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): (rng.random(samples), rng.random(samples), rng.random(samples)) ) -phasor_plot( +plot_phasor( *phasor_from_lifetime( frequency=[40e6, 80e6, 160e6], lifetime=lifetime_distribution, @@ -217,6 +95,7 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): unit_conversion=1.0, ), fmt='.', + label=('40 MHz', '80 MHz', '160 MHz'), ) # %% @@ -224,24 +103,24 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): # --------------- # # The phasor coordinates of a fluorescence energy transfer donor -# with a lifetime of 4.2 ns as a function of FRET efficiency +# with a single lifetime component of 4.2 ns as a function of FRET efficiency # at a frequency of 80 MHz, with some background signal and about 90 % -# of the donors participating in energy transfer, are on a curved trajectory: +# of the donors participating in energy transfer, are on a curved trajectory. +# For comparison, when 100% donors participate in FRET and there is no +# background signal, the phasor coordinates lie on the universal semicircle: samples = 25 efficiency = numpy.linspace(0.0, 1.0, samples) -# for reference, just donor with FRET -ax = phasor_plot( +plot = PhasorPlot(frequency=80.0) +plot.plot( *phasor_from_lifetime(80.0, 4.2 * (1.0 - efficiency)), + label='100% Donor in FRET', fmt='k.', - show=False, - return_ax=True, ) - -phasor_plot( +plot.plot( *phasor_from_lifetime( - frequency=80.0, + 80.0, lifetime=numpy.column_stack( ( numpy.full(samples, 4.2), # donor-only lifetime @@ -252,9 +131,10 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): fraction=[0.1, 0.9, 0.1 / 1e9], preexponential=True, ), + label='90% Donor in FRET', fmt='o-', - ax=ax, ) +plot.show() # %% # Multi-frequency plot @@ -266,9 +146,12 @@ def multi_frequency_plot(frequency, phase, modulation, title=None): frequency = numpy.logspace(-1, 4, 32) fraction = numpy.array([[1, 0], [0.5, 0.5], [0, 1]]) -multi_frequency_plot( +plot_polar_frequency( frequency, *phasor_to_polar( *phasor_from_lifetime(frequency, [3.9788735, 0.9947183], fraction) ), ) + +# %% +# sphinx_gallery_thumbnail_number = -2 diff --git a/tutorials/phasorpy_phasorplot.py b/tutorials/phasorpy_phasorplot.py new file mode 100644 index 00000000..b6376148 --- /dev/null +++ b/tutorials/phasorpy_phasorplot.py @@ -0,0 +1,177 @@ +""" +Phasor plot +=========== + +An introduction to the `PhasorPlot` class. + +The :py:class:`phasorpy.plot.PhasorPlot` class is used to ... + +""" + +# %% +# Import required modules, functions, and classes: + +import math + +import numpy +from matplotlib import pyplot + +from phasorpy.plot import PhasorPlot, plot_phasor + +# %% +# Empty phasor plot +# ----------------- +# +# Create an empty phasor plot, showing the first quadrant and the +# universal semicricle: + +plot = PhasorPlot() +plot.show() + +# %% +# Universal semicircle +# -------------------- +# +# Create a phasor plot at a frequency of 80 MHz and custom axis limits. +# Add a second, transformed universal semicircle: + +plot = PhasorPlot(frequency=80.0, xlim=(-0.2, 1.05)) +plot.semicircle(polar_reference=(0.9852, 0.5526)) +plot.show() + +# %% +# Scatter and line plots +# ---------------------- +# +# Plot phasor coordinates as scatter and/or lines: + +plot = PhasorPlot(frequency=80.0, title='Scatter and line plots') +plot.plot(0.6, 0.4, label='1') +plot.plot([0.2, 0.9], [0.4, 0.3], '.-', label='2') +plot.plot([0.39, 0.4, 0.41], [0.21, 0.19, 0.2], 'x', label='3') +plot.show() + +# %% +# Polar cursors +# ------------- +# +# Point out certain polar coordinates, and ranges thereof: + +plot = PhasorPlot(frequency=80.0, title='Polar cursors') +plot.polar_cursor(0.6435, 0.5) +plot.polar_cursor(0.5236, 0.6, 0.1963, 0.8) +plot.polar_cursor(0.3233, 0.9482, radius=0.05) +plot.show() + +# %% +# Component mixtures +# ------------------ +# +# Show linear combinations of phasor coordinates or ranges thereof: + +real, imag, weights = [0.1, 0.2, 0.5, 0.9], [0.3, 0.4, 0.5, 0.3], [2, 1, 2, 1] + +plot = PhasorPlot(frequency=80.0, title='Component mixtures') +plot.components(real, imag, linestyle='', fill=True, facecolor='lightyellow') +plot.components(real, imag, weights) +plot.show() + +# %% +# 2D Histogram +# ------------ +# +# Plot large number of phasor coordinates as a 2D histogram: + +real, imag = numpy.random.multivariate_normal( + (0.6, 0.4), [[3e-3, -1e-3], [-1e-3, 1e-3]], (256, 256) +).T +plot = PhasorPlot(frequency=80.0, title='2D Histogram') +plot.hist2d(real, imag) +plot.show() + +# %% +# Contours +# -------- +# +# Plot the contours of the density of phasor coordinates: + +plot = PhasorPlot(frequency=80.0, title='Contours') +plot.contour(real, imag) +plot.show() + + +# %% +# Image +# ----- +# +# Plot the image of a custom-colored 2D histogram: + +plot = PhasorPlot(frequency=80.0, title='Image (not implemented yet)') +# plot.imshow(image) +plot.show() + +# %% +# Combined plots +# -------------- +# +# Multiple plots can be combined: + +real2, imag2 = numpy.random.multivariate_normal( + (0.9, 0.2), [[2e-4, -1e-4], [-1e-4, 2e-4]], 4096 +).T + +plot = PhasorPlot( + title='Combined plots', xlim=(0.35, 1.03), ylim=(0.1, 0.59), grid=False +) +plot.hist2d(real, imag, bins=64, cmap='Blues') +plot.contour(real, imag, bins=48, levels=3, cmap='summer_r', norm='log') +plot.hist2d(real2, imag2, bins=64, cmap='Oranges') +plot.plot(0.6, 0.4, '.', color='tab:blue') +plot.plot(0.9, 0.2, '.', color='tab:orange') +plot.polar_cursor(math.atan(0.4 / 0.6), math.hypot(0.6, 0.4), color='tab:blue') +plot.polar_cursor( + math.atan(0.2 / 0.9), math.hypot(0.9, 0.2), color='tab:orange' +) +plot.semicircle(frequency=80.0, color='tab:purple') +plot.show() + +# %% +# All quadrants +# ------------- +# +# Create an empty phasor plot showing all four quadrants: + +plot = PhasorPlot(allquadrants=True, title='All quadrants') +plot.show() + +# %% +# Matplotlib axes +# --------------- +# +# The PhasorPlot class can use an existing matlotlib axis. +# The `PhasorPlot.ax` attribute provides access to the underlying +# matplotlib axis, for example, to add annotations: + +ax = pyplot.subplot(1, 1, 1) +plot = PhasorPlot(ax=ax, allquadrants=True, title='Matplotlib axes') +plot.hist2d(real, imag, cmap='Blues') +plot.ax.annotate( + '0.6, 0.4', + xy=(0.6, 0.4), + xytext=(0.2, 0.2), + arrowprops=dict(arrowstyle='->'), +) +pyplot.show() + + +# %% +# plot_phasor function +# -------------------- +# +# The :py:func:`phasorpy.plot.plot_phasor` function provides a simpler +# alternative to plot phasor coordinates in a single statement: + +plot_phasor(real[0, :32], imag[0, :32], fmt='.', frequency=80.0) + +# %% +# sphinx_gallery_thumbnail_number = 9