Skip to content

Registration

register(msims, transform_key=None, 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, n_parallel_pairwise_regs=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 n_parallel_pairwise_regs : int, optional Number of parallel pairwise registrations to run. Setting this is specifically useful for limiting memory usage. By default None (all pairwise registrations are run in parallel) 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
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
def register(
    msims: list[MultiscaleSpatialImage],
    transform_key: str = None,
    reg_channel_index: int = None,
    reg_channel: str = None,
    new_transform_key: str = None,
    registration_binning: dict[str, int] = None,
    overlap_tolerance: Union[int, dict[str, int]] = 0.0,
    pairwise_reg_func=phase_correlation_registration,
    pairwise_reg_func_kwargs: dict = None,
    groupwise_resolution_method="global_optimization",
    groupwise_resolution_kwargs: dict = None,
    pre_registration_pruning_method="alternating_pattern",
    post_registration_do_quality_filter: bool = False,
    post_registration_quality_threshold: float = 0.2,
    plot_summary: bool = False,
    pairs: list[tuple[int, int]] = None,
    scheduler=None,
    n_parallel_pairwise_regs: int = None,
    return_dict: bool = 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
    n_parallel_pairwise_regs : int, optional
        Number of parallel pairwise registrations to run. Setting this is specifically
        useful for limiting memory usage.
        By default None (all pairwise registrations are run in parallel)
    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,
        n_parallel_pairwise_regs=n_parallel_pairwise_regs,
    )

    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