Registration

register(msims, transform_key, reg_channel_index=None, reg_channel=None, new_transform_key=None, registration_binning=None, overlap_tolerance=0.0, pairwise_reg_func=phase_correlation_registration, pairwise_reg_func_kwargs=None, groupwise_resolution_method='global_optimization', groupwise_resolution_kwargs=None, pre_registration_pruning_method='alternating_pattern', post_registration_do_quality_filter=False, post_registration_quality_threshold=0.2, plot_summary=False, pairs=None, scheduler=None, return_dict=False)

Register a list of views to a common extrinsic coordinate system.

This function is the main entry point for registration.

1) Build a graph of pairwise overlaps between views 2) Determine registration pairs from this graph 3) Register each pair of views. Need to add option to pass registration functions here. 4) Determine the parameters mapping each view into the new extrinsic coordinate system. Currently done by determining a reference view and concatenating for reach view the pairwise transforms along the shortest paths towards the ref view.

Parameters

msims : list of MultiscaleSpatialImage Input views reg_channel_index : int, optional Index of channel to be used for registration, by default None reg_channel : str, optional Name of channel to be used for registration, by default None Overrides reg_channel_index transform_key : str, optional Extrinsic coordinate system to use as a starting point for the registration, by default None new_transform_key : str, optional If set, the registration result will be registered as a new extrinsic coordinate system in the input views (with the given name), by default None registration_binning : dict, optional Binning applied to each dimensionn during registration, by default None overlap_tolerance : float, optional Extend overlap regions considered for pairwise registration. - if 0, the overlap region is the intersection of the bounding boxes. - if > 0, the overlap region is the intersection of the bounding boxes extended by this value in all spatial dimensions. - if None, the full images are used for registration pairwise_reg_func : Callable, optional Function used for registration. pairwise_reg_func_kwargs : dict, optional Additional keyword arguments passed to the registration function groupwise_resolution_method : str, optional Method used to determine the final transform parameters from pairwise registrations: - 'global_optimization': global optimization considering all pairwise transforms - 'shortest_paths': concatenation of pairwise transforms along shortest paths groupwise_resolution_kwargs : dict, optional Additional keyword arguments passed to the groupwise optimization function pre_registration_pruning_method : str, optional Method used to eliminate registration edges (e.g. diagonals) from the view adjacency graph before registration. Available methods: - None: No pruning, useful when no regular arrangement is present. - 'alternating_pattern': Prune to edges between squares of differering colors in checkerboard pattern. Useful for regular 2D tile arrangements (of both 2D or 3D data). - 'shortest_paths_overlap_weighted': Prune to shortest paths in overlap graph (weighted by overlap). Useful to minimize the number of pairwise registrations. - 'otsu_threshold_on_overlap': Prune to edges with overlap above Otsu threshold. This is useful for regular 2D or 3D grid arrangements, as diagonal edges will be pruned. - 'keep_axis_aligned': Keep only edges that align with tile axes. This is useful for regular grid arrangements and to explicitely prune diagonals, e.g. when other methods fail. post_registration_do_quality_filter : bool, optional post_registration_quality_threshold : float, optional Threshold used to filter edges by quality after registration, by default None (no filtering) plot_summary : bool, optional If True (and new_transform_key is set), plot graphs summarising the registration process and results: 1) Cross correlation values of pairwise registrations (stack boundaries shown as before registration) 2) Residual distances between registration edges after global parameter resolution. Grey edges have been removed during glob param res (stack boundaries shown as after registration). Stack boundary positions reflect the registration result. By default False pairs : list of tuples, optional If set, initialises the view adjacency graph using the indicates pairs of view/tile indices, by default None scheduler : str, optional Dask scheduler to use for parallel computation, by default None return_dict : bool, optional If True, return a dict containing params, registration metrics and more, by default False

Returns

list of xr.DataArray Parameters mapping each view into a new extrinsic coordinate system or dict Dictionary containing the following keys: - 'params': Parameters mapping each view into a new extrinsic coordinate system - 'pairwise_registration': Dictionary containing the following - 'summary_plot': Tuple containing the figure and axis of the summary plot - 'graph': networkx graph of pairwise registrations - 'metrics': Dictionary containing the following metrics: - 'qualities': Edge registration qualities - 'groupwise_resolution': Dictionary containing the following - 'summary_plot': Tuple containing the figure and axis of the summary plot - 'graph': networkx graph of groupwise resolution - 'metrics': Dictionary containing the following metrics: - 'residuals': Edge residuals after groupwise resolution

Source code in src/multiview_stitcher/registration.py

def register(
    msims: list[MultiscaleSpatialImage],
    transform_key,
    reg_channel_index=None,
    reg_channel=None,
    new_transform_key=None,
    registration_binning=None,
    overlap_tolerance: Union[int, dict[str, int]] = 0.0,
    pairwise_reg_func=phase_correlation_registration,
    pairwise_reg_func_kwargs=None,
    groupwise_resolution_method="global_optimization",
    groupwise_resolution_kwargs=None,
    pre_registration_pruning_method="alternating_pattern",
    post_registration_do_quality_filter=False,
    post_registration_quality_threshold=0.2,
    plot_summary=False,
    pairs=None,
    scheduler=None,
    return_dict=False,
):
    """

    Register a list of views to a common extrinsic coordinate system.

    This function is the main entry point for registration.

    1) Build a graph of pairwise overlaps between views
    2) Determine registration pairs from this graph
    3) Register each pair of views.
       Need to add option to pass registration functions here.
    4) Determine the parameters mapping each view into the new extrinsic
       coordinate system.
       Currently done by determining a reference view and concatenating for reach
       view the pairwise transforms along the shortest paths towards the ref view.

    Parameters
    ----------
    msims : list of MultiscaleSpatialImage
        Input views
    reg_channel_index : int, optional
        Index of channel to be used for registration, by default None
    reg_channel : str, optional
        Name of channel to be used for registration, by default None
        Overrides reg_channel_index
    transform_key : str, optional
        Extrinsic coordinate system to use as a starting point
        for the registration, by default None
    new_transform_key : str, optional
        If set, the registration result will be registered as a new extrinsic
        coordinate system in the input views (with the given name), by default None
    registration_binning : dict, optional
        Binning applied to each dimensionn during registration, by default None
    overlap_tolerance : float, optional
        Extend overlap regions considered for pairwise registration.
        - if 0, the overlap region is the intersection of the bounding boxes.
        - if > 0, the overlap region is the intersection of the bounding boxes
            extended by this value in all spatial dimensions.
        - if None, the full images are used for registration
    pairwise_reg_func : Callable, optional
        Function used for registration.
    pairwise_reg_func_kwargs : dict, optional
        Additional keyword arguments passed to the registration function
    groupwise_resolution_method : str, optional
        Method used to determine the final transform parameters
        from pairwise registrations:
        - 'global_optimization': global optimization considering all pairwise transforms
        - 'shortest_paths': concatenation of pairwise transforms along shortest paths
    groupwise_resolution_kwargs : dict, optional
        Additional keyword arguments passed to the groupwise optimization function
    pre_registration_pruning_method : str, optional
        Method used to eliminate registration edges (e.g. diagonals) from the view adjacency
        graph before registration. Available methods:
        - None: No pruning, useful when no regular arrangement is present.
        - 'alternating_pattern': Prune to edges between squares of differering
            colors in checkerboard pattern. Useful for regular 2D tile arrangements (of both 2D or 3D data).
        - 'shortest_paths_overlap_weighted': Prune to shortest paths in overlap graph
            (weighted by overlap). Useful to minimize the number of pairwise registrations.
        - 'otsu_threshold_on_overlap': Prune to edges with overlap above Otsu threshold.
            This is useful for regular 2D or 3D grid arrangements, as diagonal edges will be pruned.
        - 'keep_axis_aligned': Keep only edges that align with tile axes. This is useful for regular grid
            arrangements and to explicitely prune diagonals, e.g. when other methods fail.
    post_registration_do_quality_filter : bool, optional
    post_registration_quality_threshold : float, optional
        Threshold used to filter edges by quality after registration,
        by default None (no filtering)
    plot_summary : bool, optional
        If True (and `new_transform_key` is set), plot graphs summarising the registration process and results:
        1) Cross correlation values of pairwise registrations
           (stack boundaries shown as before registration)
        2) Residual distances between registration edges after global parameter resolution.
           Grey edges have been removed during glob param res (stack boundaries shown as after registration).
        Stack boundary positions reflect the registration result.
        By default False
    pairs : list of tuples, optional
        If set, initialises the view adjacency graph using the indicates
        pairs of view/tile indices, by default None
    scheduler : str, optional
        Dask scheduler to use for parallel computation, by default None
    return_dict : bool, optional
        If True, return a dict containing params, registration metrics and more, by default False

    Returns
    -------
    list of xr.DataArray
        Parameters mapping each view into a new extrinsic coordinate system
    or
    dict
        Dictionary containing the following keys:
        - 'params': Parameters mapping each view into a new extrinsic coordinate system
        - 'pairwise_registration': Dictionary containing the following
            - 'summary_plot': Tuple containing the figure and axis of the summary plot
            - 'graph': networkx graph of pairwise registrations
            - 'metrics': Dictionary containing the following metrics:
                - 'qualities': Edge registration qualities
        - 'groupwise_resolution': Dictionary containing the following
            - 'summary_plot': Tuple containing the figure and axis of the summary plot
            - 'graph': networkx graph of groupwise resolution
            - 'metrics': Dictionary containing the following metrics:
                - 'residuals': Edge residuals after groupwise resolution
    """

    if pairwise_reg_func_kwargs is None:
        pairwise_reg_func_kwargs = {}

    if groupwise_resolution_kwargs is None:
        groupwise_resolution_kwargs = {}

    sims = [msi_utils.get_sim_from_msim(msim) for msim in msims]

    if "c" in msi_utils.get_dims(msims[0]):
        if reg_channel is None:
            if reg_channel_index is None:
                for msim in msims:
                    if "c" in msi_utils.get_dims(msim):
                        raise (
                            Exception("Please choose a registration channel.")
                        )
            else:
                reg_channel = sims[0].coords["c"][reg_channel_index]

        msims_reg = [
            msi_utils.multiscale_sel_coords(msim, {"c": reg_channel})
            if "c" in msi_utils.get_dims(msim)
            else msim
            for imsim, msim in enumerate(msims)
        ]
    else:
        msims_reg = msims

    g = mv_graph.build_view_adjacency_graph_from_msims(
        msims_reg,
        transform_key=transform_key,
        pairs=pairs,
        overlap_tolerance=overlap_tolerance,
    )

    if pre_registration_pruning_method is not None:
        g_reg = prune_view_adjacency_graph(
            g,
            method=pre_registration_pruning_method,
        )
    else:
        g_reg = g

    g_reg_computed = compute_pairwise_registrations(
        msims_reg,
        g_reg,
        transform_key=transform_key,
        registration_binning=registration_binning,
        overlap_tolerance=overlap_tolerance,
        pairwise_reg_func=pairwise_reg_func,
        pairwise_reg_func_kwargs=pairwise_reg_func_kwargs,
        scheduler=scheduler,
    )

    if post_registration_do_quality_filter:
        # filter edges by quality
        g_reg_computed = mv_graph.filter_edges(
            g_reg_computed,
            threshold=post_registration_quality_threshold,
            weight_key="quality",
        )

    params, groupwise_resolution_info_dict = groupwise_resolution(
        g_reg_computed,
        method=groupwise_resolution_method,
        **groupwise_resolution_kwargs,
    )

    params = [params[iview] for iview in sorted(g_reg_computed.nodes())]

    if new_transform_key is not None:
        for imsim, msim in enumerate(msims):
            msi_utils.set_affine_transform(
                msim,
                params[imsim],
                transform_key=new_transform_key,
                base_transform_key=transform_key,
            )

        if plot_summary or return_dict:
            edges = list(g_reg_computed.edges())
            fig_pair_reg, ax_pair_reg = vis_utils.plot_positions(
                msims,
                transform_key=transform_key,
                edges=edges,
                edge_color_vals=np.array(
                    [
                        g_reg_computed.get_edge_data(*e)["quality"].mean()
                        for e in edges
                    ]
                ),
                edge_label="Pairwise view correlation",
                display_view_indices=True,
                use_positional_colors=False,
                plot_title="Pairwise registration summary",
                show_plot=plot_summary,
            )

            if groupwise_resolution_info_dict is not None:
                edge_residuals = np.array(
                    [
                        groupwise_resolution_info_dict[
                            "optimized_graph_t0"
                        ].get_edge_data(*e)["residual"]
                        if e
                        in groupwise_resolution_info_dict[
                            "optimized_graph_t0"
                        ].edges
                        else np.nan
                        for e in edges
                    ]
                )
                edge_clims = [
                    np.nanmin(edge_residuals),
                    np.nanmax(edge_residuals),
                ]
                if edge_clims[0] == edge_clims[1]:
                    edge_clims = [0, 1]
                fig_group_res, ax_group_res = vis_utils.plot_positions(
                    msims,
                    transform_key=new_transform_key,
                    edges=edges,
                    edge_color_vals=edge_residuals,
                    edge_cmap="Spectral_r",
                    edge_clims=edge_clims,
                    edge_label="Remaining edge residuals [distance units]",
                    display_view_indices=True,
                    use_positional_colors=False,
                    plot_title="Global parameter resolution summary",
                    show_plot=plot_summary,
                )
            else:
                fig_group_res, ax_group_res = None, None
    else:
        fig_pair_reg, ax_pair_reg, fig_group_res, ax_group_res = [
            None,
            None,
            None,
            None,
        ]

    if return_dict:
        # limit output graphs to first timepoint (for now)
        g_reg_computed = g_reg_computed.copy()
        for e in g_reg_computed.edges:
            for k, v in g_reg_computed.edges[e].items():
                if isinstance(v, xr.DataArray) and "t" in v.coords:
                    g_reg_computed.edges[e][k] = g_reg_computed.edges[e][
                        k
                    ].sel({"t": g_reg_computed.edges[e][k].coords["t"][0]})

        return {
            "params": params,
            "pairwise_registration": {
                "summary_plot": (fig_pair_reg, ax_pair_reg),
                "graph": g_reg_computed,
                "metrics": {
                    "qualities": nx.get_edge_attributes(
                        g_reg_computed, "quality"
                    ),
                },
            },
            "groupwise_resolution": {
                "summary_plot": (fig_group_res, ax_group_res),
                "graph": groupwise_resolution_info_dict["optimized_graph_t0"],
                "metrics": {
                    "residuals": nx.get_edge_attributes(
                        groupwise_resolution_info_dict["optimized_graph_t0"],
                        "residual",
                    ),
                },
            },
        }
    else:
        return params