Stitch

Shift

view_stitch_search

Function to plot results of exhaustive search to find overlap between tile t and its neighbours. Useful for debugging the stitch section of the pipeline. White in the color plot refers to the value of score_thresh for this search.

Parameters:

Name	Type	Description	Default
nb	Notebook	Notebook containing results of the experiment. Must contain find_spots page.	required
t	int	Want to look at overlap between tile t and its north/east neighbour.	required
direction	Optional[str]	Direction of overlap interested in - either 'south'/'north' or 'west'/'east'. If None, then will look at both directions.	None

Source code in coppafish/plot/stitch/shift.py

def view_stitch_search(nb: Notebook, t: int, direction: Optional[str] = None):
    """
    Function to plot results of exhaustive search to find overlap between tile `t` and its neighbours.
    Useful for debugging the `stitch` section of the pipeline.
    White in the color plot refers to the value of `score_thresh` for this search.

    Args:
        nb: Notebook containing results of the experiment. Must contain `find_spots` page.
        t: Want to look at overlap between tile `t` and its north/east neighbour.
        direction: Direction of overlap interested in - either `'south'`/`'north'` or `'west'`/`'east'`.
            If `None`, then will look at both directions.
    """
    # NOTE that directions should actually be 'north' and 'east'
    if direction is None:
        directions = ['south', 'west']
    elif direction.lower() == 'south' or direction.lower() == 'west':
        directions = [direction.lower()]
    elif direction.lower() == 'north':
        directions = ['south']
    elif direction.lower() == 'east':
        directions = ['west']
    else:
        raise ValueError(f"direction must be either 'south' or 'west' but {direction} given.")
    direction_label = {'south': 'north', 'west': 'east'}  # label refers to actual direction

    config = nb.get_config()['stitch']
    # determine shifts to search over
    shifts = get_shifts_to_search(config, nb.basic_info)
    if not nb.basic_info.is_3d:
        config['nz_collapse'] = None
        config['shift_widen'][2] = 0  # so don't look for shifts in z direction
        config['shift_max_range'][2] = 0

    # find shifts between overlapping tiles
    c = nb.basic_info.ref_channel
    r = nb.basic_info.ref_round
    t_neighb = {'south': [], 'west': []}
    z_scale = nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy

    # align to south neighbour followed by west neighbour
    t_neighb['south'] = np.where(np.sum(nb.basic_info.tilepos_yx == nb.basic_info.tilepos_yx[t, :] + [1, 0],
                                        axis=1) == 2)[0]
    t_neighb['west'] = np.where(np.sum(nb.basic_info.tilepos_yx == nb.basic_info.tilepos_yx[t, :] + [0, 1],
                                       axis=1) == 2)[0]
    fig = []
    for j in directions:
        if t_neighb[j] in nb.basic_info.use_tiles:
            print(f'Finding shift between tiles {t} and {t_neighb[j][0]} ({direction_label[j]} overlap)')
            shift, score, score_thresh, debug_info = \
                compute_shift(spot_yxz(nb.find_spots.spot_details, t, r, c),
                              spot_yxz(nb.find_spots.spot_details, t_neighb[j][0], r, c),
                              config['shift_score_thresh'],
                              config['shift_score_thresh_multiplier'],
                              config['shift_score_thresh_min_dist'],
                              config['shift_score_thresh_max_dist'],
                              config['neighb_dist_thresh'], shifts[j]['y'],
                              shifts[j]['x'], shifts[j]['z'],
                              config['shift_widen'], config['shift_max_range'],
                              z_scale, config['nz_collapse'],
                              config['shift_step'][2])
            title = f'Overlap between t={t} and neighbor in {direction_label[j]} (t={t_neighb[j][0]}). ' \
                    f'YXZ Shift = {shift}.'
            fig = fig + [view_shifts(debug_info['shifts_2d'], debug_info['scores_2d'], debug_info['shifts_3d'],
                                     debug_info['scores_3d'], shift, debug_info['min_score_2d'],
                                     debug_info['shift_2d_initial'],
                                     score_thresh, debug_info['shift_thresh'], config['shift_score_thresh_min_dist'],
                                     config['shift_score_thresh_max_dist'], title, False)]
    if len(fig) > 0:
        plt.show()
    else:
        warnings.warn(f"Tile {t} has no overlapping tiles in nb.basic_info.use_tiles.")

view_shifts

Function to plot scores indicating number of neighbours between 2 point clouds corresponding to particular shifts applied to one of them. I.e. you can use this to view the output from coppafish/stitch/shift/compute_shift function.

Parameters:

Name	Type	Description	Default
shifts_2d	np.ndarray	int [n_shifts_2d x 2]. shifts_2d[i] is the yx shift which achieved scores_2d[i] when considering just yx shift between point clouds. I.e. first step of finding optimal shift is collapsing 3D point cloud to just a few planes and then applying a yx shift to these planes.	required
scores_2d	np.ndarray	float [n_shifts_2d]. scores_2d[i] is the score corresponding to shifts_2d[i]. It is approximately the number of neighbours between the two point clouds after the shift was applied.	required
shifts_3d	Optional[np.ndarray]	int [n_shifts_3d x 3]. shifts_3d[i] is the yxz shift which achieved scores_3d[i] when considering the yxz shift between point clouds. YX shift is in units of YX pixels. Z shift is in units of z-pixels. If None, only 2D image plotted.	None
scores_3d	Optional[np.ndarray]	float [n_shifts_3d]. scores_3d[i] is the score corresponding to shifts_3d[i]. It is approximately the number of neighbours between the two point clouds after the shift was applied.	None
best_shift	Optional[np.ndarray]	int [y_shift, x_shift, z_shift]. Best shift found by algorithm. YX shift is in units of YX pixels. Z shift is in units of z-pixels. Will be plotted as black cross on image if provided.	None
score_thresh_2d	Optional[float]	Threshold returned by compute_shift function for 2d calculation, if score is above this, it indicates an accepted 2D shift. If given, a red-white-blue colorbar will be used with white corresponding to score_thresh_2d in the 2D plot	None
best_shift_initial	Optional[np.ndarray]	int [y_shift, x_shift]. Best yx shift found by in first search of algorithm. I.e. score_thresh computation based on this. Will show as green x if given.	None
score_thresh_3d	Optional[float]	Threshold returned by compute_shift function for 3d calculation. If given, a red-white-blue colorbar will be used with white corresponding to score_thresh_3d in the 3D plots.	None
thresh_shift	Optional[np.ndarray]	int [y_shift, x_shift, z_shift]. yx shift corresponding to score_thresh. Will show as green + in both 2D and 3D plots if given.	None
thresh_min_dist	Optional[int]	shift_thresh is the shift with the max score in an annulus a distance between thresh_min_dist and thresh_max_dist away from best_shift_initial. Annulus will be shown in green if given.	None
thresh_max_dist	Optional[int]	shift_thresh is the shift with the max score in an annulus a distance between thresh_min_dist and thresh_max_dist away from best_shift_initial. Annulus will be shown in green if given.	None
title	Optional[str]	Title to show.	None
show	bool	If True, will call plt.show(), else will return fig.	True

Source code in coppafish/plot/stitch/shift.py

def view_shifts(shifts_2d: np.ndarray, scores_2d: np.ndarray, shifts_3d: Optional[np.ndarray] = None,
                scores_3d: Optional[np.ndarray] = None, best_shift: Optional[np.ndarray] = None,
                score_thresh_2d: Optional[float] = None, best_shift_initial: Optional[np.ndarray] = None,
                score_thresh_3d: Optional[float] = None, thresh_shift: Optional[np.ndarray] = None,
                thresh_min_dist: Optional[int] = None,
                thresh_max_dist: Optional[int] = None, title: Optional[str] = None, show: bool = True):
    """
    Function to plot scores indicating number of neighbours between 2 point clouds corresponding to particular shifts
    applied to one of them. I.e. you can use this to view the output from
    `coppafish/stitch/shift/compute_shift` function.

    Args:
        shifts_2d: `int [n_shifts_2d x 2]`.
            `shifts_2d[i]` is the yx shift which achieved `scores_2d[i]` when considering just yx shift between
            point clouds.
            I.e. first step of finding optimal shift is collapsing 3D point cloud to just a few planes and then
            applying a yx shift to these planes.
        scores_2d: `float [n_shifts_2d]`.
            `scores_2d[i]` is the score corresponding to `shifts_2d[i]`. It is approximately the number of neighbours
            between the two point clouds after the shift was applied.
        shifts_3d: `int [n_shifts_3d x 3]`.
            `shifts_3d[i]` is the yxz shift which achieved `scores_3d[i]` when considering the yxz shift between
            point clouds. YX shift is in units of YX pixels. Z shift is in units of z-pixels.
            If None, only 2D image plotted.
        scores_3d: `float [n_shifts_3d]`.
            `scores_3d[i]` is the score corresponding to `shifts_3d[i]`. It is approximately the number of neighbours
            between the two point clouds after the shift was applied.
        best_shift: `int [y_shift, x_shift, z_shift]`.
            Best shift found by algorithm. YX shift is in units of YX pixels. Z shift is in units of z-pixels.
            Will be plotted as black cross on image if provided.
        score_thresh_2d: Threshold returned by `compute_shift` function for 2d calculation, if `score` is above this,
            it indicates an accepted 2D shift. If given, a red-white-blue colorbar will be used with white corresponding
            to `score_thresh_2d` in the 2D plot
        best_shift_initial: `int [y_shift, x_shift]`.
            Best yx shift found by in first search of algorithm. I.e. `score_thresh` computation based on this.
            Will show as green x if given.
        score_thresh_3d: Threshold returned by `compute_shift` function for 3d calculation. If given, a red-white-blue
            colorbar will be used with white corresponding to `score_thresh_3d` in the 3D plots.
        thresh_shift: `int [y_shift, x_shift, z_shift]`.
            yx shift corresponding to `score_thresh`. Will show as green + in both 2D and 3D plots if given.
        thresh_min_dist: `shift_thresh` is the shift with the max score in an annulus a distance between
            `thresh_min_dist` and `thresh_max_dist` away from `best_shift_initial`.
            Annulus will be shown in green if given.
        thresh_max_dist: `shift_thresh` is the shift with the max score in an annulus a distance between
            `thresh_min_dist` and `thresh_max_dist` away from `best_shift_initial`.
            Annulus will be shown in green if given.
        title: Title to show.
        show: If `True`, will call `plt.show()`, else will `return fig`.
    """
    image_2d, extent_2d = get_plot_images_from_shifts(np.rint(shifts_2d).astype(int), scores_2d)
    image_2d = interpolate_array(image_2d, 0)  # replace 0 with nearest neighbor value
    fig = plt.figure(figsize=(12, 8))
    if score_thresh_2d is None:
        score_thresh_2d = (image_2d.min() + image_2d.max()) / 2
        cmap_2d = 'virids'
    else:
        cmap_2d = 'bwr'
    v_max = np.max([image_2d.max(), 1.2 * score_thresh_2d])
    v_min = image_2d.min()
    if cmap_2d == 'bwr':
        cmap_extent = np.max([v_max - score_thresh_2d, score_thresh_2d - v_min])
        # Have equal range above and below score_thresh to not skew colormap
        v_min = score_thresh_2d - cmap_extent
        v_max = score_thresh_2d + cmap_extent
    cmap_norm = matplotlib.colors.TwoSlopeNorm(vmin=v_min, vcenter=score_thresh_2d, vmax=v_max)
    if shifts_3d is not None:
        images_3d, extent_3d = get_plot_images_from_shifts(np.rint(shifts_3d).astype(int), scores_3d)
        images_3d = interpolate_array(images_3d, 0)  # replace 0 with nearest neighbor value
        if score_thresh_3d is None:
            score_thresh_3d = (images_3d.min() + images_3d.max()) / 2
            cmap_3d = 'virids'
        else:
            cmap_3d = 'bwr'
        v_max_3d = np.max([images_3d.max(), 1.2 * score_thresh_3d])
        v_min_3d = images_3d.min()
        if cmap_3d == 'bwr':
            cmap_extent = np.max([v_max_3d - score_thresh_3d, score_thresh_3d - v_min_3d])
            # Have equal range above and below score_thresh to not skew colormap
            v_min_3d = score_thresh_3d - cmap_extent
            v_max_3d = score_thresh_3d + cmap_extent
        cmap_norm_3d = matplotlib.colors.TwoSlopeNorm(vmin=v_min_3d, vcenter=score_thresh_3d, vmax=v_max_3d)
        n_cols = images_3d.shape[2]
        if n_cols > 13:
            # If loads of z-planes, just show the 13 with the largest score
            n_cols = 13
            max_score_z = images_3d.max(axis=(0, 1))
            use_z = np.sort(np.argsort(max_score_z)[::-1][:n_cols])
        else:
            use_z = np.arange(n_cols)

        plot_3d_height = int(np.ceil(n_cols / 4))
        plot_2d_height = n_cols - plot_3d_height
        ax_2d = plt.subplot2grid(shape=(n_cols, n_cols), loc=(0, 0), colspan=n_cols, rowspan=plot_2d_height)
        ax_3d = [plt.subplot2grid(shape=(n_cols, n_cols), loc=(plot_2d_height + 1, i), rowspan=plot_3d_height) for i in
                 range(n_cols)]
        for i in range(n_cols):
            # share axes for 3D plots
            ax_3d[i].get_shared_y_axes().join(ax_3d[i], *ax_3d)
            ax_3d[i].get_shared_x_axes().join(ax_3d[i], *ax_3d)
            im_3d = ax_3d[i].imshow(images_3d[:, :, use_z[i]], extent=extent_3d[:4], aspect='auto', cmap=cmap_3d,
                                    norm=cmap_norm_3d)
            z_plane = int(np.rint(extent_3d[4] + use_z[i] + 0.5))
            ax_3d[i].set_title(f'Z = {z_plane}')
            if thresh_shift is not None and z_plane == thresh_shift[2]:
                # Indicate threshold shift on correct 3d plot
                ax_3d[i].plot(thresh_shift[1], thresh_shift[0], '+', color='lime', label='Thresh shift')
            if i > 0:
                ax_3d[i].tick_params(labelbottom=False, labelleft=False)
            if best_shift is not None:
                if z_plane == best_shift[2]:
                    ax_3d[i].plot(best_shift[1], best_shift[0], 'kx')

        fig.supxlabel('X')
        fig.supylabel('Y')
        ax_3d[0].invert_yaxis()
        cbar_gap = 0.05
        cbar_3d_height = plot_3d_height / (plot_3d_height + plot_2d_height)
        cbar_2d_height = plot_2d_height / (plot_3d_height + plot_2d_height)
        cbar_ax = fig.add_axes([0.9, 0.07, 0.03, cbar_3d_height - 2*cbar_gap])  # left, bottom, width, height
        fig.colorbar(im_3d, cax=cbar_ax)
        cbar_ax.set_ylim(np.clip(v_min_3d, 0, v_max_3d), v_max_3d)
    else:
        n_cols = 1
        ax_2d = plt.subplot2grid(shape=(n_cols, n_cols), loc=(0, 0))
        ax_2d.set_xlabel('X')
        ax_2d.set_ylabel('Y')

    im_2d = ax_2d.imshow(image_2d, extent=extent_2d, aspect='auto', cmap=cmap_2d, norm=cmap_norm)
    if best_shift is not None:
        ax_2d.plot(best_shift[1], best_shift[0], 'kx', label='Best shift')
    if best_shift_initial is not None and best_shift_initial != best_shift:
        ax_2d.plot(best_shift_initial[1], best_shift_initial[0], 'x', color='lime', label='Best shift initial')
    if thresh_shift is not None:
        ax_2d.plot(thresh_shift[1], thresh_shift[0], '+', color='lime', label='Thresh shift')
    if thresh_min_dist is not None and best_shift_initial is not None:
        ax_2d.add_patch(plt.Circle((best_shift_initial[1], best_shift_initial[0]), thresh_min_dist, color='lime',
                                   fill=False))
    if thresh_max_dist is not None and best_shift_initial is not None:
        ax_2d.add_patch(plt.Circle((best_shift_initial[1], best_shift_initial[0]), thresh_max_dist, color='lime',
                                   fill=False))
    ax_2d.invert_yaxis()
    if title is None:
        title = 'Approx number of neighbours found for all shifts'
    ax_2d.set_title(title)
    ax_2d.legend(facecolor='b')
    fig.subplots_adjust(left=0.07, right=0.85, bottom=0.07, top=0.95)
    if shifts_3d is None:
        cbar_ax = fig.add_axes([0.9, 0.07, 0.03, 0.9])  # left, bottom, width, height
    else:
        cbar_ax = fig.add_axes([0.9, 0.07 + cbar_3d_height, 0.03, cbar_2d_height - 2.5 * cbar_gap])
    fig.colorbar(im_2d, cax=cbar_ax)
    cbar_ax.set_ylim(np.clip(v_min, 0, v_max), v_max)
    if show:
        plt.show()
    else:
        return fig

Point Clouds

view_stitch_overlap

This plots point clouds of neighbouring tiles with:

• No overlap
• Initial guess at shift using config['stitch']['expected_overlap']
• Overlap determined in stitch stage of the pipeline (using nb.stitch.south_shifts or nb.stitch.west_shifts)
• Their final global coordinate system positions (using nb.stitch.tile_origin)

Parameters:

Name	Type	Description	Default
nb	Notebook	Notebook containing at least stitch page.	required
t	int	Want to look at overlap between tile t and its north or east neighbour.	required
direction	str	Direction of overlap interested in - either 'south'/'north' or 'west'/'east'.	'south'

Source code in coppafish/plot/stitch/point_clouds.py

def view_stitch_overlap(nb: Notebook, t: int, direction: str = 'south'):
    """
    This plots point clouds of neighbouring tiles with:

    * No overlap
    * Initial guess at shift using `config['stitch']['expected_overlap']`
    * Overlap determined in stitch stage of the pipeline (using `nb.stitch.south_shifts` or `nb.stitch.west_shifts`)
    * Their final global coordinate system positions (using `nb.stitch.tile_origin`)

    Args:
        nb: *Notebook* containing at least `stitch` page.
        t: Want to look at overlap between tile `t` and its north or east neighbour.
        direction: Direction of overlap interested in - either `'south'`/`'north'` or `'west'`/`'east'`.
    """
    # NOTE that directions should actually be 'north' and 'east'
    if direction.lower() == 'south' or direction.lower() == 'west':
        direction = direction.lower()
    elif direction.lower() == 'north':
        direction = 'south'
    elif direction.lower() == 'east':
        direction = 'west'
    else:
        raise ValueError(f"direction must be either 'south' or 'west' but {direction} given.")

    direction_label = {'south': 'north', 'west': 'east'}  # label refers to actual direction
    if direction == 'south':
        t_neighb = np.where(np.sum(nb.basic_info.tilepos_yx == nb.basic_info.tilepos_yx[t, :] + [1, 0],
                                   axis=1) == 2)[0]
        if t_neighb not in nb.basic_info.use_tiles:
            warnings.warn(f"Tile {t} has no overlapping tiles in the south direction so changing to west.")
            direction = 'west'
        else:
            no_overlap_shift = np.array([-nb.basic_info.tile_sz, 0, 0])  # assuming no overlap between tiles
            found_shift = nb.stitch.south_shifts[np.where(nb.stitch.south_pairs[:,0] == t)[0]][0]
    if direction == 'west':
        t_neighb = np.where(np.sum(nb.basic_info.tilepos_yx == nb.basic_info.tilepos_yx[t, :] + [0, 1],
                                   axis=1) == 2)[0]
        if t_neighb not in nb.basic_info.use_tiles:
            raise ValueError(f"Tile {t} has no overlapping tiles in the west direction.")
        no_overlap_shift = np.array([0, -nb.basic_info.tile_sz, 0])  # assuming no overlap between tiles
        found_shift = nb.stitch.west_shifts[np.where(nb.stitch.west_pairs[:, 0] == t)[0]][0]

    config = nb.get_config()['stitch']
    t_neighb = t_neighb[0]
    r = nb.basic_info.ref_round
    c = nb.basic_info.ref_channel
    point_clouds = []
    # add global coordinates of neighbour tile as point cloud that is always present.
    point_clouds = point_clouds + [spot_yxz(nb.find_spots.spot_details, t_neighb, r, c) +
                                   nb.stitch.tile_origin[t_neighb]]

    local_yxz_t = spot_yxz(nb.find_spots.spot_details, t, r, c)
    # Add point cloud for tile t assuming no overlap
    point_clouds = point_clouds + [local_yxz_t + nb.stitch.tile_origin[t_neighb] + no_overlap_shift]
    # Add point cloud assuming expected overlap
    initial_shift = (1-config['expected_overlap']) * no_overlap_shift
    point_clouds = point_clouds + [local_yxz_t + nb.stitch.tile_origin[t_neighb] + initial_shift]
    # Add point cloud for tile t with found shift
    point_clouds = point_clouds + [local_yxz_t + nb.stitch.tile_origin[t_neighb] + found_shift]
    # Add point cloud for tile t in global coordinate system
    point_clouds = point_clouds + [local_yxz_t + nb.stitch.tile_origin[t]]

    neighb_dist_thresh = config['neighb_dist_thresh']
    z_scale = nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy
    pc_labels = [f'Tile {t_neighb}', f'Tile {t} - No overlap',
                 f"Tile {t} - {int(config['expected_overlap']*100)}% overlap", f'Tile {t} - Shift', f'Tile {t} - Final']
    view_point_clouds(point_clouds, pc_labels, neighb_dist_thresh, z_scale,
                      f'Overlap between tile {t} and tile {t_neighb} in the {direction_label[direction]}')
    plt.show()

view_stitch

This plots all the reference spots found (ref_round/ref_channel) in the global coordinate system created in the stitch stage of the pipeline.

It also indicates which of these spots are duplicates (detected on a tile which is not the tile whose centre they are closest to). These will be removed in the get_reference_spots step of the pipeline so we don't double count the same spot.

Parameters:

Name	Type	Description	Default
nb	Notebook	Notebook containing at least stitch page.	required

Source code in coppafish/plot/stitch/point_clouds.py

def view_stitch(nb: Notebook):
    """
    This plots all the reference spots found (`ref_round`/`ref_channel`) in the global coordinate system created
    in the `stitch` stage of the pipeline.

    It also indicates which of these spots are duplicates (detected on a tile which is not the tile whose centre
    they are closest to).
    These will be removed in the `get_reference_spots` step of the pipeline so we don't double count the same spot.

    Args:
        nb: *Notebook* containing at least `stitch` page.
    """
    is_ref = np.all((nb.find_spots.spot_details[:, 1] == nb.basic_info.ref_round,
                     nb.find_spots.spot_details[:, 2] == nb.basic_info.ref_channel), axis=0)
    local_yxz = nb.find_spots.spot_details[is_ref, -3:]
    tile = nb.find_spots.spot_details[is_ref, 0]
    local_yxz = local_yxz[np.isin(tile, nb.basic_info.use_tiles)]
    tile = tile[np.isin(tile, nb.basic_info.use_tiles)]

    # find duplicate spots as those detected on a tile which is not tile centre they are closest to
    tile_origin = nb.stitch.tile_origin
    not_duplicate = get_non_duplicate(tile_origin, nb.basic_info.use_tiles, nb.basic_info.tile_centre,
                                      local_yxz, tile)

    global_yxz = local_yxz + nb.stitch.tile_origin[tile]
    global_yxz[:, 2] = np.rint(global_yxz[:, 2])  # make z coordinate an integer
    config = nb.get_config()['stitch']
    neighb_dist_thresh = config['neighb_dist_thresh']
    z_scale = nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy

    point_clouds = [global_yxz[not_duplicate], global_yxz[np.invert(not_duplicate)]]
    pc_labels = ['Not Duplicate', 'Duplicate']
    if nb.has_page('ref_spots'):
        # Add point cloud indicating those spots that were not saved in ref_spots page
        # because they were shifted outside the tile bounds on at least one round/channel
        # and thus spot_color could not be found.
        local_yxz = local_yxz[not_duplicate]  # Want to find those which were not duplicates but still removed
        tile = tile[not_duplicate]
        local_yxz_saved = nb.ref_spots.local_yxz
        tile_saved = nb.ref_spots.tile
        local_yxz_saved = local_yxz_saved[np.isin(tile_saved, nb.basic_info.use_tiles)]
        tile_saved = tile_saved[np.isin(tile_saved, nb.basic_info.use_tiles)]
        global_yxz_ns = np.zeros((0, 3)) # not saved in ref_spots spots
        for t in nb.basic_info.use_tiles:
            missing_ind = -100
            local_yxz_t = local_yxz[tile == t]
            removed_ind = np.where(numpy_indexed.indices(local_yxz_saved[tile_saved == t], local_yxz_t,
                                                         missing=missing_ind) == missing_ind)[0]
            global_yxz_ns = np.append(global_yxz_ns, local_yxz_t[removed_ind] + nb.stitch.tile_origin[t], axis=0)
        global_yxz_ns[:, 2] = np.rint(global_yxz_ns[:, 2])
        point_clouds += [global_yxz_ns]
        pc_labels += ["No Spot Color"]

    # Sometimes can be empty point cloud, so remove these
    use_pc = [len(pc) > 0 for pc in point_clouds]
    pc_labels = [pc_labels[i] for i in range(len(use_pc)) if use_pc[i]]
    point_clouds= [point_clouds[i] for i in range(len(use_pc)) if use_pc[i]]
    vpc = view_point_clouds(point_clouds, pc_labels, neighb_dist_thresh, z_scale,
                            "Reference Spots in the Global Coordinate System")

    tile_sz = nb.basic_info.tile_sz
    for t in nb.basic_info.use_tiles:
        rect = matplotlib.patches.Rectangle((tile_origin[t, 1], tile_origin[t, 0]), tile_sz, tile_sz,
                                            linewidth=1, edgecolor='w', facecolor='none', linestyle=':')
        vpc.ax.add_patch(rect)
        vpc.ax.text(tile_origin[t, 1] + 20, tile_origin[t, 0] + 20, f"Tile {t}",
                    size=6, color='w', ha='left', weight='light')
    plt.show()

view_point_clouds

Plots two point clouds. point_clouds[0] always plotted but you can change the second point cloud using radio buttons.

Parameters:

Name	Type	Description	Default
point_clouds	List	List of point clouds, each of which is yxz coordinates float [n_points x 3]. point_clouds[0] is always plotted. YX coordinates are in units of yx pixels. Z coordinates are in units of z pixels. Radio buttons used to select other point_cloud plotted.	required
pc_labels	List	List of labels to appear in legend/radio-buttons for each point cloud. Must provide one for each point_cloud.	required
neighb_dist_thresh	float	If distance between neighbours is less than this, a white line will connect them.	5
z_scale	float	pixel_size_z / pixel_size_y i.e. used to convert z coordinates from z-pixels to yx pixels.	1
super_title	Optional[str]	Optional title for plot	None

Source code in coppafish/plot/stitch/point_clouds.py

def __init__(self, point_clouds: List, pc_labels: List, neighb_dist_thresh: float = 5, z_scale: float = 1,
             super_title: Optional[str] = None):
    """
    Plots two point clouds. `point_clouds[0]` always plotted but you can change the second point cloud using
    radio buttons.

    Args:
        point_clouds: List of point clouds, each of which is yxz coordinates `float [n_points x 3]`.
            `point_clouds[0]` is always plotted.
            YX coordinates are in units of yx pixels. Z coordinates are in units of z pixels.
            Radio buttons used to select other point_cloud plotted.
        pc_labels: List of labels to appear in legend/radio-buttons for each point cloud.
            Must provide one for each `point_cloud`.
        neighb_dist_thresh: If distance between neighbours is less than this, a white line will connect them.
        z_scale: `pixel_size_z / pixel_size_y` i.e. used to convert z coordinates from z-pixels to yx pixels.
        super_title: Optional title for plot
    """
    n_point_clouds = len(point_clouds)
    if len(point_clouds) != len(pc_labels):
        raise ValueError(f'There are {n_point_clouds} point clouds but {len(pc_labels)} labels')
    self.fig, self.ax = plt.subplots(1, 1, figsize=(15, 8))
    subplots_adjust = [0.07, 0.775, 0.095, 0.89]
    self.fig.subplots_adjust(left=subplots_adjust[0], right=subplots_adjust[1], bottom=subplots_adjust[2],
                             top=subplots_adjust[3])

    # Find neighbours between all point clouds and the first.
    self.neighb = [None]
    n_matches_str = 'Number Of Matches'
    for i in range(1, n_point_clouds):
        tree = KDTree(point_clouds[i] * [1, 1, z_scale])
        dist, neighb = tree.query(point_clouds[0] * [1, 1, z_scale],
                                  distance_upper_bound=neighb_dist_thresh * 3)
        neighb[dist > neighb_dist_thresh] = -1  # set too large distance neighb_ind to -1.
        n_matches_str = n_matches_str + f'\n- {pc_labels[i]}: {np.sum(neighb >= 0)}'
        self.neighb = self.neighb + [neighb]

    # Add text box to indicate number of matches between first point cloud and each of the others.
    n_matches_ax = self.fig.add_axes([0.8, subplots_adjust[3] - 0.75, 0.15, 0.2])
    plt.axis('off')
    n_matches_ax.text(0.05, 0.95, n_matches_str)

    #  round z to closest z-plane for plotting
    for i in range(n_point_clouds):
        point_clouds[i][:, 2] = np.rint(point_clouds[i][:, 2])

    # get yxz axis limits taking account all point clouds
    pc_min_lims = np.zeros((n_point_clouds, 3))
    pc_max_lims = np.zeros_like(pc_min_lims)
    for i in range(n_point_clouds):
        pc_min_lims[i] = np.min(point_clouds[i], axis=0)
        pc_max_lims[i] = np.max(point_clouds[i], axis=0)

    self.z_planes = np.arange(int(np.min(pc_min_lims[:, 2])), int(np.max(pc_max_lims[:, 2]) + 1))
    self.nz = len(self.z_planes)
    self.z_ind = 0
    self.z = self.z_planes[self.z_ind]
    self.z_thick = 0
    self.active_pc = [0, 1]
    self.in_z = [np.array([val[:, 2] >= self.z - self.z_thick, val[:, 2] <= self.z + self.z_thick]).all(axis=0)
                 for val in point_clouds]
    self.point_clouds = point_clouds
    self.pc_labels = np.array(pc_labels)
    self.pc_shapes = ['rx', 'bo']
    alpha = [1, 0.7]
    self.pc_plots = [self.ax.plot(point_clouds[self.active_pc[i]][self.in_z[i], 1],
                                  point_clouds[self.active_pc[i]][self.in_z[i], 0],
                                  self.pc_shapes[i], label=self.pc_labels[i], alpha=alpha[i])[0]
                     for i in range(2)]

    self.neighb_yx = None
    self.neighb_plot = None
    self.update_neighb_lines()

    self.ax.legend(loc='upper right')
    self.ax.set_ylabel('Y')
    self.ax.set_xlabel('X')
    self.ax.set_ylim(np.min(pc_min_lims[:, 0]), np.max(pc_max_lims[:, 0]))
    self.ax.set_xlim(np.min(pc_min_lims[:, 1]), np.max(pc_max_lims[:, 1]))

    if self.nz > 1:
        # If 3D, add text box to change number of z-planes collapsed onto single plane
        # and add scrolling to change z-plane
        self.ax.set_title(f'Z = {int(self.z)}', size=10)
        self.fig.canvas.mpl_connect('scroll_event', self.z_scroll)
        text_ax = self.fig.add_axes([0.8, 0.095, 0.15, 0.04])
    else:
        # For some reason in 2D, still need the text box otherwise buttons don't do work
        # But shift it off-screen and make small
        text_ax = self.fig.add_axes([40, 40, 0.00001, 0.00001])
    self.text_box = TextBox(text_ax, 'Z-Thick', self.z_thick, color='k', hovercolor=[0.2, 0.2, 0.2])
    self.text_box.cursor.set_color('r')
    # change text box title to be above not to the left of box
    label = text_ax.get_children()[0]  # label is a child of the TextBox axis
    if self.nz == 1:
        label.set_position([40, 40])  # shift label off-screen in 2D
    else:
        label.set_position([0.5, 2])  # [x,y] - change here to set the position
    # centering the text
    label.set_verticalalignment('top')
    label.set_horizontalalignment('center')
    self.text_box.on_submit(self.text_update)

    if n_point_clouds >= 3:
        # If 3 or more point clouds, add radio button to change the second point cloud shown.
        buttons_ax = self.fig.add_axes([0.8, subplots_adjust[3] - 0.2, 0.15, 0.2])
        plt.axis('off')
        self.buttons = RadioButtons(buttons_ax, self.pc_labels[1:], 0, activecolor='w')
        for i in range(n_point_clouds-1):
            self.buttons.circles[i].set_color('w')
            self.buttons.circles[i].set_color('w')
        self.buttons.set_active(0)
        self.buttons.on_clicked(self.button_update)
    if super_title is not None:
        plt.suptitle(super_title, x=(0.07 + 0.775) / 2)

Diagnostics

view_stitch_shift_info

For all north/south and east/west shifts computed in the stitch section of the pipeline, this plots the values of the shifts found and the score compared to the score_thresh.

For each direction, there will be 3 plots:

• y shift vs x shift for all pairs of neighbouring tiles
• z shift vs x shift for all pairs of neighbouring tiles
• score vs score_thresh for all pairs of neighbouring tiles (a green score = score_thresh line is plotted in this).

In each case, the markers in the plots are numbers. These numbers indicate the tile, the shift was applied to, to take it to its north or east neighbour i.e. nb.stitch.south_pairs[:, 0] or nb.stitch.west_pairs[:, 0]. The number will be blue if score > score_thresh and red otherwise.

Parameters:

Name	Type	Description	Default
nb	Notebook	Notebook containing at least the stitch page.	required
outlier	bool	If True, will plot nb.stitch.south_shift_outlier instead of nb.stitch.south_shift. In this case, only tiles for which the two are different are plotted for each round.	False

Source code in coppafish/plot/stitch/diagnostics.py

def view_stitch_shift_info(nb: Notebook, outlier: bool = False):
    """
    For all north/south and east/west shifts computed in the `stitch` section
    of the pipeline, this plots the values of the shifts found and the `score` compared to
    the `score_thresh`.

    For each direction, there will be 3 plots:

    * y shift vs x shift for all pairs of neighbouring tiles
    * z shift vs x shift for all pairs of neighbouring tiles
    * `score` vs `score_thresh` for all pairs of neighbouring tiles
    (a green score = score_thresh line is plotted in this).

    In each case, the markers in the plots are numbers.
    These numbers indicate the tile, the shift was applied to,
    to take it to its north or east neighbour i.e. `nb.stitch.south_pairs[:, 0]`
    or `nb.stitch.west_pairs[:, 0]`.
    The number will be blue if `score > score_thresh` and red otherwise.

    Args:
        nb: Notebook containing at least the `stitch` page.
        outlier: If `True`, will plot `nb.stitch.south_shift_outlier` instead of
            `nb.stitch.south_shift`. In this case, only tiles for which
            the two are different are plotted for each round.
    """
    if nb.basic_info.is_3d:
        ndim = 3
    else:
        ndim = 2
    shift_info = {}
    if len(nb.stitch.south_shifts) > 0:
        shift_info['South'] = {}
        shift_info['South']['tile'] = nb.stitch.south_pairs[:, 0]
        shift_info['South']['score_thresh'] = nb.stitch.south_score_thresh
        if outlier:
            shift_info['South']['shift'] = nb.stitch.south_outlier_shifts[:, :ndim]
            shift_info['South']['score'] = nb.stitch.south_outlier_score
        else:
            shift_info['South']['shift'] = nb.stitch.south_shifts[:, :ndim]
            shift_info['South']['score'] = nb.stitch.south_score
    if len(nb.stitch.west_shifts) > 0:
        shift_info['West'] = {}
        shift_info['West']['tile'] = nb.stitch.west_pairs[:, 0]
        shift_info['West']['score_thresh'] = nb.stitch.west_score_thresh
        if outlier:
            shift_info['West']['shift'] = nb.stitch.west_outlier_shifts[:, :ndim]
            shift_info['West']['score'] = nb.stitch.west_outlier_score
        else:
            shift_info['West']['shift'] = nb.stitch.west_shifts[:, :ndim]
            shift_info['West']['score'] = nb.stitch.west_score
    if outlier:
        title_start = "Outlier "
    else:
        title_start = ""
    shift_info_plot(shift_info, f"{title_start}Shifts found in stitch part of pipeline between each tile and the "
                                f"neighbouring tile in the direction specified")

shift_info_plot

If shift_info contains \(n\) keys, this will produce an \(n\) column x 3 row grid of subplots. For each key in shift_info dictionary, there are 3 plots:

• y shift vs x shift
• z shift vs x shift
• score vs score_thresh or n_matches vs error

In each case, the markers in the plots are numbers. These numbers are given by shift_info[key][tile]. The number will be blue if score > score_thresh and red otherwise.

Parameters:

Name	Type	Description	Default
shift_info	dict	Dictionary containing \(n\) dictionaries. Each of these dictionaries contains (either (score and score_thresh) or (n_matches, n_matches_thresh and error)): shift - float [n_tiles x 3]. \(yxz\) shifts for each tile. tile - int [n_tiles]. Indicates tile each shift was found for. score - float [n_tiles]. Indicates score found for each shift (approx number of matches between point clouds). score_thresh - float [n_tiles]. If score<score_thresh, it will be shown in red. n_matches - int [n_tiles]. Indicates score found for each shift (approx number of matches between point clouds). n_matches_thresh - int [n_tiles]. If n_matches<n_matches_thresh, it will be shown in red. error - float [n_tiles]. Average distance between neighbours. x_lim - float [n_plots x 2]. Can optionally specify number axis limits for each plot. y_lim - float [n_plots x 2]. Can optionally specify number axis limits for each plot.	required
title	Optional[str]	Overall title for the plot.	None
score_plot_thresh	int	Only shifts with score (or n_matches) > score_plot_thresh are shown.	0
fig	Optional[plt.Figure]	Can provide previous figure to plot on.	None
ax	Optional[np.ndarray]	Can provide array of plt.Axes to plot on.	None
return_ax	bool	If True, ax will be returned and plt.show() will not be run.	False

Source code in coppafish/plot/stitch/diagnostics.py

def shift_info_plot(shift_info: dict, title: Optional[str] = None, score_plot_thresh: int = 0,
                    fig: Optional[plt.Figure] = None, ax: Optional[np.ndarray] = None, return_ax: bool = False):
    """
    If `shift_info` contains $n$ keys, this will produce an $n$ column x 3 row grid of subplots.
    For each key in `shift_info` dictionary, there are 3 plots:

    * y shift vs x shift
    * z shift vs x shift
    * `score` vs `score_thresh` or `n_matches` vs `error`

    In each case, the markers in the plots are numbers.
    These numbers are given by `shift_info[key][tile]`.
    The number will be blue if `score > score_thresh` and red otherwise.

    Args:
        shift_info: Dictionary containing $n$ dictionaries.
            Each of these dictionaries contains (either (`score` and `score_thresh`) or (`n_matches`,
            `n_matches_thresh` and `error`)):

            * shift - `float [n_tiles x 3]`. $yxz$ shifts for each tile.
            * tile - `int [n_tiles]`. Indicates tile each shift was found for.
            * score - `float [n_tiles]`. Indicates `score` found for each shift (approx number of matches between
                point clouds).
            * score_thresh - `float [n_tiles]`. If `score<score_thresh`, it will be shown in red.
            * n_matches - `int [n_tiles]`. Indicates `score` found for each shift (approx number of matches between
                point clouds).
            * n_matches_thresh - `int [n_tiles]`. If `n_matches<n_matches_thresh`, it will be shown in red.
            * error - `float [n_tiles]`. Average distance between neighbours.
            * x_lim - `float [n_plots x 2]`. Can optionally specify number axis limits for each plot.
            * y_lim - `float [n_plots x 2]`. Can optionally specify number axis limits for each plot.
        title: Overall title for the plot.
        score_plot_thresh: Only shifts with `score` (or `n_matches`) > `score_plot_thresh` are shown.
        fig: Can provide previous figure to plot on.
        ax: Can provide array of plt.Axes to plot on.
        return_ax: If `True`, ax will be returned and `plt.show()` will not be run.
    """
    n_cols = len(shift_info)
    col_titles = list(shift_info.keys())
    n_rows = len(shift_info[col_titles[0]]['shift'][0])  # 2 if 2D shift, 3 if 3D.
    if fig is None:
        fig, ax = plt.subplots(n_rows, n_cols, figsize=(15, 7))
        fig.subplots_adjust(hspace=0.4, bottom=0.08, left=0.06, right=0.97, top=0.9)
    if n_cols == 1 and len(ax.shape) == 1:
        ax = ax[:, np.newaxis]
    for i in range(n_cols):
        shift_info_i = shift_info[col_titles[i]]
        # good tiles are blue, bad tiles are red
        n_tiles = len(shift_info[col_titles[i]]['tile'])
        tile_color = np.full(n_tiles, 'b')
        if 'score' in shift_info_i:
            tile_color[(shift_info_i['score'] < shift_info_i['score_thresh']).flatten()] = 'r'
            skip_tile = shift_info_i['score'] <= score_plot_thresh  # don't plot if score = 0
        elif 'n_matches' in shift_info_i:
            tile_color[(shift_info_i['n_matches'] < shift_info_i['n_matches_thresh']).flatten()] = 'r'
            skip_tile = shift_info_i['n_matches'] <= score_plot_thresh  # don't plot if n_matches = 0
        else:
            raise ValueError(f"shift_info must contain either a 'score' or 'n_matches' key")
        for t in range(n_tiles):
            if skip_tile[t]:
                continue
            ax[0, i].text(shift_info_i['shift'][t, 1], shift_info_i['shift'][t, 0],
                          str(shift_info_i['tile'][t]), color=tile_color[t], fontsize=12,
                          ha='center', va='center')
        if 'x_lim' in shift_info_i:
            ax[0, i].set_xlim(shift_info_i['x_lim'][0])
        else:
            ax[0, i].set_xlim([np.min(shift_info_i['shift'][:, 1]) - 3, np.max(shift_info_i['shift'][:, 1]) + 3])
        if 'y_lim' in shift_info_i:
            ax[0, i].set_ylim(shift_info_i['y_lim'][0])
        else:
            ax[0, i].set_ylim([np.min(shift_info_i['shift'][:, 0]) - 3, np.max(shift_info_i['shift'][:, 0]) + 3])
        if i == int(np.ceil(n_cols/2)-1):
            ax[0, i].set_xlabel('X Shift')
        if i == 0:
            ax[0, i].set_ylabel('Y Shift')
        ax[0, i].set_title(col_titles[i])

        row_ind = 1
        if n_rows == 3:
            for t in range(n_tiles):
                if skip_tile[t]:
                    continue
                ax[row_ind, i].text(shift_info_i['shift'][t, 1], shift_info_i['shift'][t, 2],
                                    str(shift_info_i['tile'][t]), color=tile_color[t], fontsize=12,
                                    ha='center', va='center')
            if 'x_lim' in shift_info_i:
                ax[row_ind, i].set_xlim(shift_info_i['x_lim'][row_ind])
            else:
                ax[row_ind, i].set_xlim([np.min(shift_info_i['shift'][:, 1]) - 3,
                                         np.max(shift_info_i['shift'][:, 1]) + 3])
            if 'y_lim' in shift_info_i:
                ax[row_ind, i].set_ylim(shift_info_i['y_lim'][row_ind])
            else:
                ax[row_ind, i].set_ylim([np.min(shift_info_i['shift'][:, 2]) - 1,
                                         np.max(shift_info_i['shift'][:, 2]) + 1])
            if i == int(np.ceil(n_cols/2)-1):
                ax[row_ind, i].set_xlabel('X Shift')
            if i == 0:
                ax[row_ind, i].set_ylabel('Z Shift')
            row_ind += 1

        if 'score' in shift_info_i:
            # Plot line so that all with score > score_thresh are above it
            if 'x_lim' in shift_info_i and 'y_lim' in shift_info_i:
                score_min = np.min(shift_info_i['y_lim'][row_ind, 0], shift_info_i['x_lim'][row_ind, 0]) - 5
                score_max = np.max(shift_info_i['y_lim'][row_ind, 1], shift_info_i['x_lim'][row_ind, 1]) + 5
            else:
                score_min = np.min(np.vstack([shift_info_i['score_thresh'], shift_info_i['score']])) - 5
                score_max = np.max(np.vstack([shift_info_i['score_thresh'], shift_info_i['score']])) + 5
            ax[row_ind, i].set_xlim([score_min, score_max])
            ax[row_ind, i].set_ylim([score_min, score_max])
            ax[row_ind, i].plot([score_min, score_max], [score_min, score_max], 'lime', linestyle=':', linewidth=2,
                                alpha=0.5)
            for t in range(n_tiles):
                if skip_tile[t]:
                    continue
                ax[row_ind, i].text(shift_info_i['score_thresh'][t], shift_info_i['score'][t],
                                    str(shift_info_i['tile'][t]), color=tile_color[t], fontsize=12,
                                    ha='center', va='center')
            if i == int(np.ceil(n_cols/2)-1):
                ax[row_ind, i].set_xlabel('Score Threshold')
            if i == 0:
                ax[row_ind, i].set_ylabel('Score')
        elif 'n_matches' in shift_info_i:
            if 'x_lim' in shift_info_i:
                ax[row_ind, i].set_xlim(shift_info_i['x_lim'][row_ind])
            else:
                ax[row_ind, i].set_xlim([np.min(shift_info_i['error']) - 0.1, np.max(shift_info_i['error']) + 0.1])
            if 'y_lim' in shift_info_i:
                ax[row_ind, i].set_ylim(shift_info_i['y_lim'][row_ind])
            else:
                ax[row_ind, i].set_ylim([np.min(shift_info_i['n_matches']) - 100,
                                         np.max(shift_info_i['n_matches']) + 100])
            for t in range(n_tiles):
                if skip_tile[t]:
                    continue
                ax[row_ind, i].text(shift_info_i['error'][t], shift_info_i['n_matches'][t],
                                    str(shift_info_i['tile'][t]), color=tile_color[t], fontsize=12,
                                    ha='center', va='center')
            if i == int(np.ceil(n_cols/2)-1):
                ax[row_ind, i].set_xlabel('Error')
            if i == 0:
                ax[row_ind, i].set_ylabel(r'$n_{matches}$')
    if title is not None:
        plt.suptitle(title)
    if return_ax:
        return ax
    else:
        plt.show()