
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/registration/plot_stitching.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_registration_plot_stitching.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_registration_plot_stitching.py:


===========================================
Assemble images with simple image stitching
===========================================

This example demonstrates how a set of images can be assembled under
the hypothesis of rigid body motions.

.. GENERATED FROM PYTHON SOURCE LINES 10-61

.. code-block:: Python


    from matplotlib import pyplot as plt
    import numpy as np
    from skimage import data, util, transform, feature, measure, filters, metrics


    def match_locations(img0, img1, coords0, coords1, radius=5, sigma=3):
        """Match image locations using SSD minimization.

        Areas from `img0` are matched with areas from `img1`. These areas
        are defined as patches located around pixels with Gaussian
        weights.

        Parameters
        ----------
        img0, img1 : 2D array
            Input images.
        coords0 : (2, m) array_like
            Centers of the reference patches in `img0`.
        coords1 : (2, n) array_like
            Centers of the candidate patches in `img1`.
        radius : int
            Radius of the considered patches.
        sigma : float
            Standard deviation of the Gaussian kernel centered over the patches.

        Returns
        -------
        match_coords: (2, m) array
            The points in `coords1` that are the closest corresponding matches to
            those in `coords0` as determined by the (Gaussian weighted) sum of
            squared differences between patches surrounding each point.
        """
        y, x = np.mgrid[-radius : radius + 1, -radius : radius + 1]
        weights = np.exp(-0.5 * (x**2 + y**2) / sigma**2)
        weights /= 2 * np.pi * sigma * sigma

        match_list = []
        for r0, c0 in coords0:
            roi0 = img0[r0 - radius : r0 + radius + 1, c0 - radius : c0 + radius + 1]
            roi1_list = [
                img1[r1 - radius : r1 + radius + 1, c1 - radius : c1 + radius + 1]
                for r1, c1 in coords1
            ]
            # sum of squared differences
            ssd_list = [np.sum(weights * (roi0 - roi1) ** 2) for roi1 in roi1_list]
            match_list.append(coords1[np.argmin(ssd_list)])

        return np.array(match_list)









.. GENERATED FROM PYTHON SOURCE LINES 62-66

Data generation
---------------

For this example, we generate a list of slightly tilted noisy images.

.. GENERATED FROM PYTHON SOURCE LINES 66-85

.. code-block:: Python


    img = data.moon()

    angle_list = [0, 5, 6, -2, 3, -4]
    center_list = [(0, 0), (10, 10), (5, 12), (11, 21), (21, 17), (43, 15)]

    img_list = [
        transform.rotate(img, angle=a, center=c)[40:240, 50:350]
        for a, c in zip(angle_list, center_list)
    ]
    ref_img = img_list[0].copy()

    img_list = [
        util.random_noise(filters.gaussian(im, sigma=1.1), var=5e-4, rng=seed)
        for seed, im in enumerate(img_list)
    ]

    psnr_ref = metrics.peak_signal_noise_ratio(ref_img, img_list[0])








.. GENERATED FROM PYTHON SOURCE LINES 86-96

Image registration
------------------

.. note::
   This step is performed using the approach described in
   :ref:`sphx_glr_auto_examples_transform_plot_matching.py`, but any other
   method from the :ref:`sphx_glr_auto_examples_registration` section can be
   applied, depending on your problem.

Reference points are detected over all images in the list.

.. GENERATED FROM PYTHON SOURCE LINES 96-105

.. code-block:: Python


    min_dist = 5
    corner_list = [
        feature.corner_peaks(
            feature.corner_harris(img), threshold_rel=0.001, min_distance=min_dist
        )
        for img in img_list
    ]








.. GENERATED FROM PYTHON SOURCE LINES 106-109

The Harris corners detected in the first image are chosen as
references. Then the detected points on the other images are
matched to the reference points.

.. GENERATED FROM PYTHON SOURCE LINES 109-117

.. code-block:: Python


    img0 = img_list[0]
    coords0 = corner_list[0]
    matching_corners = [
        match_locations(img0, img1, coords0, coords1, min_dist)
        for img1, coords1 in zip(img_list, corner_list)
    ]








.. GENERATED FROM PYTHON SOURCE LINES 118-126

Once all the points are registered to the reference points, robust
relative affine transformations can be estimated using the RANSAC method.

.. hint::
   :func:`skimage.measure.ransac` uses random sampling to estimate a model's
   parameters. Sometimes, the number of trials may not be enough to arrive
   at a successful estimation. To make this example deterministic, we pass a
   seed ``42`` to the parameter ``rng`` that is known to be successful.

.. GENERATED FROM PYTHON SOURCE LINES 126-155

.. code-block:: Python


    src = np.array(coords0)
    trfm_list = [
        measure.ransac(
            (dst, src),
            transform.EuclideanTransform,
            min_samples=3,
            residual_threshold=2,
            max_trials=100,
            rng=42,
        )[0].params
        for dst in matching_corners
    ]

    fig, ax_list = plt.subplots(6, 2, figsize=(6, 9), sharex=True, sharey=True)
    for idx, (im, trfm, (ax0, ax1)) in enumerate(zip(img_list, trfm_list, ax_list)):
        ax0.imshow(im, cmap="gray", vmin=0, vmax=1)
        ax1.imshow(transform.warp(im, trfm), cmap="gray", vmin=0, vmax=1)

        if idx == 0:
            ax0.set_title("Tilted images")
            ax0.set_ylabel(f"Reference Image\n(PSNR={psnr_ref:.2f})")
            ax1.set_title("Registered images")

        ax0.set(xticklabels=[], yticklabels=[], xticks=[], yticks=[])
        ax1.set_axis_off()

    fig.tight_layout()




.. image-sg:: /auto_examples/registration/images/sphx_glr_plot_stitching_001.png
   :alt: Tilted images, Registered images
   :srcset: /auto_examples/registration/images/sphx_glr_plot_stitching_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 156-166

Image assembling
----------------

A composite image can be obtained using the positions of the
registered images relative to the reference one. To do so, we define
a global domain around the reference image and position the other
images in this domain.

A global transformation is defined to move the reference image in the
global domain image via a simple translation:

.. GENERATED FROM PYTHON SOURCE LINES 166-172

.. code-block:: Python

    margin = 50
    height, width = img_list[0].shape
    out_shape = height + 2 * margin, width + 2 * margin
    glob_trfm = np.eye(3)
    glob_trfm[:2, 2] = -margin, -margin








.. GENERATED FROM PYTHON SOURCE LINES 173-176

Finally, the relative position of the other images in the global
domain are obtained by composing the global transformation with the
relative transformations:

.. GENERATED FROM PYTHON SOURCE LINES 176-199

.. code-block:: Python

    global_img_list = [
        transform.warp(
            img, trfm.dot(glob_trfm), output_shape=out_shape, mode="constant", cval=np.nan
        )
        for img, trfm in zip(img_list, trfm_list)
    ]

    all_nan_mask = np.all([np.isnan(img) for img in global_img_list], axis=0)
    global_img_list[0][all_nan_mask] = 1.0

    composite_img = np.nanmean(global_img_list, 0)
    psnr_composite = metrics.peak_signal_noise_ratio(
        ref_img, composite_img[margin : margin + height, margin : margin + width]
    )

    fig, ax = plt.subplots(1, 1)

    ax.imshow(composite_img, cmap="gray", vmin=0, vmax=1)
    ax.set_axis_off()
    ax.set_title(f"Reconstructed image (PSNR={psnr_composite:.2f})")
    fig.tight_layout()

    plt.show()



.. image-sg:: /auto_examples/registration/images/sphx_glr_plot_stitching_002.png
   :alt: Reconstructed image (PSNR=36.89)
   :srcset: /auto_examples/registration/images/sphx_glr_plot_stitching_002.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.060 seconds)


.. _sphx_glr_download_auto_examples_registration_plot_stitching.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_stitching.ipynb <plot_stitching.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_stitching.py <plot_stitching.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_stitching.zip <plot_stitching.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
