Source code for starfish.core._display

import __main__
import warnings
from collections import OrderedDict
from typing import Iterable, List, Optional, Set, Tuple, Union

import numpy as np
from packaging.version import parse as parse_version

from starfish.core.imagestack.dataorder import AXES_DATA, N_AXES
from starfish.core.imagestack.imagestack import ImageStack
from starfish.core.intensity_table.intensity_table import IntensityTable
from starfish.core.morphology.binary_mask import BinaryMaskCollection
from starfish.core.types import Axes, Features

try:
    from napari import Viewer
except ImportError:
    Viewer = None


NAPARI_VERSION = "0.4.18"  # when changing this, update docs in display
INTERACTIVE = not hasattr(__main__, "__file__")


def _normalize_axes(axes: Iterable[Union[Axes, str]]) -> List[str]:
    """convert all axes to strings"""

    def _normalize(axes):
        if isinstance(axes, str):
            return axes
        elif isinstance(axes, Axes):
            return axes.value
        else:
            raise TypeError(f"passed axes must be str or Axes objects, not {type(axes)}")

    return list(_normalize(i) for i in axes)


def _max_intensity_table_maintain_dims(
    intensity_table: IntensityTable, dimensions: Set[Axes],
) -> IntensityTable:
    """
    Maximum project an IntensityTable over dimensions, retaining singletons in place of dimensions
    reduced by the max operation.

    For example, xarray.max(Axes.CH.value) would usually produce a 2-d (features, rounds) output.
    This function ensures that the output is instead (features, 1, rounds), by inserting singleton
    dimensions in place of any projected axes.

    Parameters
    ----------
    intensity_table : IntensityTable
        Input intensities
    dimensions : Set[Axes]
        Dimensions to project

    Returns
    -------
    IntensityTable :
        full-dimensionality (3-d) IntensityTable

    """
    str_dimensions = _normalize_axes(dimensions)
    initial_dimensions = OrderedDict(intensity_table.sizes)
    projected_intensities = intensity_table.max(str_dimensions)
    expanded_intensities = projected_intensities.expand_dims(str_dimensions)
    return expanded_intensities.transpose(*tuple(initial_dimensions.keys()))  # type: ignore


def _mask_low_intensity_spots(
    intensity_table: IntensityTable,
    intensity_threshold: float
) -> np.ndarray:
    """
    returns a flattened boolean mask (c order flattening) that is True where
    (feature, round, channel) combinations whose intensities are above or equal to intensity
    threshold and where values are not nan
    """
    with warnings.catch_warnings():
        # we expect to have invalid values (np.nan) in IntensityTable, hide the warnings from users
        warnings.simplefilter("ignore", RuntimeWarning)
        # reshape the values into a 1-d vector of length product(features, rounds, channels)
        return np.ravel(intensity_table.values >= intensity_threshold)


def _spots_to_markers(intensity_table: IntensityTable) -> Tuple[np.ndarray, np.ndarray]:
    """
    convert each (r, c) combination for each spot into a 5-tuple (x, y, r, c, z) for plotting
    with napari.
    """
    # get sizes which will be used to create the coordinates that correspond to the spots
    n_rounds = intensity_table.sizes[Axes.ROUND.value]
    n_channels = intensity_table.sizes[Axes.CH.value]
    n_features = intensity_table.sizes[Features.AXIS]
    code_length = n_rounds * n_channels

    # create 5-d coordinates for plotting in (x, y, round. channel, z)
    n_markers = np.product(intensity_table.shape)
    coords = np.zeros((n_markers, 5), dtype=np.uint16)

    # create the coordinates.
    # X, Y, and Z are repeated once per (r, ch) pair (code_length).
    # the cartesian product of (r, c) are created, and tiled once per feature (n_features)
    # we ensure that (ch, round) cycle in c-order to match the order of the linearized
    # array, used below for masking.
    coords[:, 0] = np.tile(np.tile(range(n_rounds), n_channels), n_features)
    coords[:, 1] = np.tile(np.repeat(range(n_channels), n_rounds), n_features)
    coords[:, 2] = np.repeat(intensity_table[Features.AXIS][Axes.ZPLANE.value].values, code_length)
    coords[:, 3] = np.repeat(intensity_table[Features.AXIS][Axes.Y.value].values, code_length)
    coords[:, 4] = np.repeat(intensity_table[Features.AXIS][Axes.X.value].values, code_length)

    sizes = np.repeat(intensity_table.radius.values, code_length)
    rc = np.zeros((sizes.shape[0], 2), dtype=int)
    z = sizes[:, np.newaxis]
    yx = np.tile(sizes[:, np.newaxis], (1, 2))
    sizes = np.concatenate((rc, yx, z), axis=1)

    return coords, sizes


[docs] def display( stack: Optional[ImageStack] = None, spots: Optional[IntensityTable] = None, masks: Optional[BinaryMaskCollection] = None, viewer: Optional[Viewer] = None, project_axes: Optional[Set[Axes]] = None, mask_intensities: float = 0., radius_multiplier: int = 1, z_multiplier: float = 1 ): """ Display an image stack, detected spots, and masks using `napari <https://github.com/napari/napari>`. Parameters ---------- stack : ImageStack ImageStack to display spots : IntensityTable IntensityTable containing spot information that was generated from the submitted stack. masks : BinaryMaskCollection Segmentation instance masks used to annotate the submitted stack. viewer : napari.Viewer Napari viewer to append the ImageStack and/or spots to. If None, creates a new viewer. Note: appending is only supported in interactive environments. project_axes : Optional[Set[Axes]] If provided, both the ImageStack and the Spots will be maximum projected along the selected axes. Useful for displaying spots across coded assays where spots may not appear in specific rounds or channels. mask_intensities : Float hide markers that correspond to intensities below this threshold value; note that any marker that is np.nan will ALSO be masked, allowing users to pre-mask the intensity table (see documentation on IntensityTable.where() for more details) (default 0, no masking) radius_multiplier : int Multiplies the radius of the displayed spots (default 1, no scaling) z_multiplier : int Multiplies the radius of the spots in z, to account for anisotropy. Important note, anisotropy of spots is DEPRECATED in napari > 0.4.17 and would be returned at a later version. Examples -------- 1. Display a stack to evaluate a filtering result. Just pass any ImageStack! >>> from starfish import display >>> display(stack) 2. Display spots of a single-molecule FISH experiment: smFISH will produce IntensityTables where most values are np.NaN. These will be masked automatically, so passing an ImageStack + IntensityTable will display spots in the rounds and channels that they are detected >>> from starfish import display >>> display(stack, intensities) 3. Diagnose spot calls within each round and channel of a coded experiment. A user might want to evaluate if spot calling is failing for a specific round/channel pair. To accomplish this, pass the intensity threshold used by the spot called to eliminate sub-threshold spots from channels. The user can additionally filter the IntensityTable to further mask additional spots (any np.NaN value will not be displayed) >>> from starfish import display >>> mask_intensities = 0.3 # this was the spot calling intensity threshold >>> display(stack, intensities, mask_intensities=mask_intensities) 4. Evaluate spot calls across rounds and channels by visualizing spots on a max projection of The rounds and channels. >>> from starfish import display, Axes >>> display(stack, intensities, project_axes={Axes.CH, Axes.ROUND}) 5. Compare the image before (raw_stack) and after (filtered_stack) filtering by displaying two stacks in the same Viewer. >>> from starfish import display >>> viewer = display(raw_stack) >>> viewer = display(stack=filtered_stack, viewer=viewer) Notes ----- - To use in ipython, use the `%gui qt` magic. - napari axes are labeled with the ImageStack axis names - Requires napari 0.4.18: use `pip install starfish[napari]` to install all necessary requirements """ if stack is None and spots is None and masks is None: raise TypeError("expected a stack, spots, or masks; got nothing") try: import napari except ImportError: raise ImportError(f"Requires at least napari {NAPARI_VERSION}. " "Run `pip install starfish[napari]` " "to install the necessary requirements.") try: version = napari.__version__ except Exception as e: raise RuntimeError("Could not identify napari version") from e if parse_version(version) < parse_version(NAPARI_VERSION): raise ValueError(f"Incorrect version {version} of napari installed." "Run `pip install starfish[napari]` " "to install the necessary requirements.") # instantiate the napari viewer if viewer is None: from qtpy.QtWidgets import QApplication app = QApplication.instance() or QApplication([]) # Get the list of axis labels axis_labels = [] for ix in range(N_AXES): for axis_name, axis_data in AXES_DATA.items(): if ix == axis_data.order: axis_labels.append(axis_name.value) break else: raise Exception(f"No AXES_DATA entry for the {ix}th axis") axis_labels.append(Axes.Y.value) axis_labels.append(Axes.X.value) viewer = Viewer(axis_labels=axis_labels) new_viewer = True elif isinstance(viewer, Viewer): new_viewer = False else: raise TypeError("viewer must be a napari.Viewer or None") if stack is not None: if project_axes is not None: stack = stack.reduce(project_axes, func="max") viewer.add_image(stack.xarray.values, rgb=False, name="stack") if spots is not None: if project_axes is not None: spots = _max_intensity_table_maintain_dims(spots, project_axes) # TODO ambrosejcarr guard rails: # 1. verify that x, y, z fit inside the imagestack (weird projections) # 2. warn if whole tiles are missing spot calls (possible miss-use) coords, sizes = _spots_to_markers(spots) # mask low-intensity values mask = _mask_low_intensity_spots(spots, mask_intensities) if not np.sum(mask): warnings.warn(f"No spots passed provided intensity threshold of {mask_intensities}") else: coords = coords[mask] sizes = sizes[mask] # adjust z-size sizes[:, 4] *= z_multiplier viewer.add_points( coords, face_color="red", edge_color="red", symbol="ring", size=np.mean(sizes[:, 2:4], axis=1) * radius_multiplier, n_dimensional=True, name="spots" ) if masks: viewer.add_labels(masks.to_label_image().label_image, name="masks") if new_viewer and not INTERACTIVE: app.exec_() # create blocking process to persist windows return viewer