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

.. only:: html

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

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

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

.. _sphx_glr_spam_examples_DIC_plot_multiModalRegistration.py:


Multimodal registration
==================================

A simple example to register two images acquired with different modalities

.. GENERATED FROM PYTHON SOURCE LINES 7-18

.. code-block:: default
   :lineno-start: 8


    import matplotlib as mpl
    import matplotlib.pyplot as plt
    import numpy
    import spam.datasets
    import spam.deformation

    # sphinx_gallery_thumbnail_number = 4
    import spam.DIC
    import spam.helpers








.. GENERATED FROM PYTHON SOURCE LINES 19-21

Load the two images
#####################

.. GENERATED FROM PYTHON SOURCE LINES 21-31

.. code-block:: default
   :lineno-start: 21

    xr = spam.datasets.loadConcreteXr().astype("<f4")
    ne = spam.datasets.loadConcreteNe().astype("<f4")

    #
    #
    # plt.figure()
    # plt.imshow(xr[halfSlice, :, :])
    # plt.figure()
    # plt.imshow(ne[halfSlice, :, :])








.. GENERATED FROM PYTHON SOURCE LINES 32-37

Preparation of the images
##########################################
 We start to crop the images to:
   1. select a region of interest: ``crop``
   2. keep some margin to feed the transformation: ``margin``

.. GENERATED FROM PYTHON SOURCE LINES 37-61

.. code-block:: default
   :lineno-start: 37

    cropRatio = 0.1
    crop = (
        slice(int(cropRatio * xr.shape[0]), int((1 - cropRatio) * xr.shape[0])),
        slice(int(cropRatio * xr.shape[1]), int((1 - cropRatio) * xr.shape[1])),
        slice(int(cropRatio * xr.shape[2]), int((1 - cropRatio) * xr.shape[2])),
    )
    cropPx = int(cropRatio * numpy.mean(xr.shape[0]))
    marginRatio = 0.1
    marginPx = int(marginRatio * numpy.mean(xr.shape[0]))
    cropWithMargin = (
        slice(
            int((cropRatio + marginRatio) * xr.shape[0]),
            int((1 - (cropRatio + marginRatio)) * xr.shape[0]),
        ),
        slice(
            int((cropRatio + marginRatio) * xr.shape[1]),
            int((1 - (cropRatio + marginRatio)) * xr.shape[1]),
        ),
        slice(
            int((cropRatio + marginRatio) * xr.shape[2]),
            int((1 - (cropRatio + marginRatio)) * xr.shape[2]),
        ),
    )








.. GENERATED FROM PYTHON SOURCE LINES 62-64

We rescale the two images between 0 and the number of bins int the joint histogram.
For sake of efficiency they *have* to be saved as 8bits integers (``<u1``).

.. GENERATED FROM PYTHON SOURCE LINES 64-72

.. code-block:: default
   :lineno-start: 64

    bins = 128
    xrMin = xr.min()
    xrMax = xr.max()
    neMax = ne.max()
    neMin = ne.min()
    xr = numpy.array(bins * (xr - xrMin) / (xrMax - xrMin)).astype("<u1")
    ne = numpy.array(bins * (ne - neMin) / (neMax - neMin)).astype("<u1")








.. GENERATED FROM PYTHON SOURCE LINES 73-75

To see the spatial incoherency we can plot them both on the same image with a checkerbord pattern,
where two adjacent squares represents the two images.

.. GENERATED FROM PYTHON SOURCE LINES 75-82

.. code-block:: default
   :lineno-start: 76


    halfSlice = xr.shape[0] // 2
    # checker = spam.helpers.checkerBoard(xr[halfSlice], ne[halfSlice], n=3)
    # plt.figure()
    # plt.imshow(checker, cmap="Greys")
    # plt.colorbar()








.. GENERATED FROM PYTHON SOURCE LINES 83-84

We apply and initial guess transformation.

.. GENERATED FROM PYTHON SOURCE LINES 84-91

.. code-block:: default
   :lineno-start: 84

    PhiGuess = spam.deformation.computePhi({"t": [0.0, 0.0, 0.0], "r": [15.0, 0.0, 0.0]})
    tmp = spam.deformation.decomposePhi(PhiGuess)
    neTmp = spam.DIC.applyPhi(ne.copy(), Phi=PhiGuess).astype("<u1")
    print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*tmp["t"]))
    print("Rotations   : {:.4f}, {:.4f}, {:.4f}".format(*tmp["r"]))






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Translations: 0.0000, 0.0000, 0.0000
    Rotations   : 15.0000, 0.0000, 0.0000




.. GENERATED FROM PYTHON SOURCE LINES 92-98

Registration algorithm
################################################
 We first get the parameters of the fitted gaussians and the joint histogram.
 x axis corresponds to the first image grey levels.
 y axis corresponds to the second image grey levels.
 The Gaussian parameters are parameters of the two ellipsoids that fit the two peaks.

.. GENERATED FROM PYTHON SOURCE LINES 98-112

.. code-block:: default
   :lineno-start: 98

    print("STEP 1: Get gaussian parameters")
    nPhases = 2
    gaussianParameters, jointHistogram = spam.DIC.gaussianMixtureParameters(xr[cropWithMargin], neTmp[cropWithMargin], BINS=bins, NPHASES=nPhases)
    plt.figure()
    tmp = jointHistogram.copy()
    tmp[jointHistogram <= 0] = numpy.nan
    tmp = numpy.log(tmp)
    plt.imshow(tmp.T, origin="lower", extent=[0.0, bins, 0.0, bins])
    plt.xlabel("x-ray grey levels")
    plt.ylabel("neutron grey levels")
    plt.colorbar()
    for gp in gaussianParameters:
        plt.plot(gp[1], gp[2], "b*")




.. image-sg:: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_001.png
   :alt: plot multiModalRegistration
   :srcset: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    STEP 1: Get gaussian parameters
            im1 from 29.00 to 128.00
            im2 from 23.00 to 122.00
            p normalisation: 1.00
            Find maxima
                    Min distance between maxima: 5
                    00 maxima found over the 2 needed with 1.00e-01 times of the total count
                    00 maxima found over the 2 needed with 1.00e-02 times of the total count
                    02 maxima found over the 2 needed with 1.00e-03 times of the total count
                    Maximum 1:       mu1=108.00      mu2=86.00       Phi=8.63e-03
                    Maximum 2:       mu1=109.00      mu2=67.00       Phi=3.59e-03

            Phase  1:                mu1=108.00      mu2=86.00       Phi=8.63e-03
    /usr/local/lib/python3.8/site-packages/spam/DIC/multimodal.py:675: RuntimeWarning: overflow encountered in cast
      field[nx, ny] = float(GLOBALphi) * numpy.exp(-computeLambda(a, b, c, x, GLOBALmu1, y, GLOBALmu2))
                    Fit:             a=0.18  b=0.02  c=0.04  Hessian: 0.01
                    Cov:             1,1=5.8691      1,2=-3.3674     2,2=30.3347

            Phase  2:                mu1=109.00      mu2=67.00       Phi=3.59e-03
                    Fit:             a=0.15  b=-0.00         c=0.00  Hessian: 0.00
                    Cov:             1,1=6.6323      1,2=0.8708      2,2=323.2468





.. GENERATED FROM PYTHON SOURCE LINES 113-116

Then we create the phase diagram based on the joint histogram.
Each peak corresponds to a phase (1 and 2).
The grey background (points to far away from a peak) is ignored (0).

.. GENERATED FROM PYTHON SOURCE LINES 116-133

.. code-block:: default
   :lineno-start: 116

    print("STEP 2: Create phase repartition")
    phaseDiagram, actualVoxelCoverage = spam.DIC.phaseDiagram(gaussianParameters, jointHistogram, voxelCoverage=0.99, BINS=bins)
    plt.figure()
    plt.imshow(
        phaseDiagram.T,
        origin="lower",
        extent=[0.0, bins, 0.0, bins],
        vmin=-0.5,
        vmax=nPhases + 0.5,
        cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
    )
    plt.xlabel("x-ray grey levels")
    plt.ylabel("neutron grey levels")
    plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
    for gp in gaussianParameters:
        plt.plot(gp[1], gp[2], "b*")




.. image-sg:: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_002.png
   :alt: plot multiModalRegistration
   :srcset: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    STEP 2: Create phase repartition
            9-sigma: voxel coverage = 99.15 (> 99.00%)
            8-sigma: voxel coverage = 98.99 (< 99.00%) -> Returning this phase diagram.
    /builds/ttk/spam/examples/DIC/plot_multiModalRegistration.py:125: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
      cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),




.. GENERATED FROM PYTHON SOURCE LINES 134-135

And and we use both Gaussian parameters and phase diagram as an input of the registration algorithm

.. GENERATED FROM PYTHON SOURCE LINES 135-149

.. code-block:: default
   :lineno-start: 135

    print("STEP 3: Registration")
    registration = spam.DIC.multimodalRegistration(
        xr,
        ne,
        phaseDiagram,
        gaussianParameters,
        BINS=bins,
        PhiInit=PhiGuess,
        verbose=True,
        margin=marginPx,
        maxIterations=50,
        deltaPhiMin=0.005,
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    STEP 3: Registration
            Initial state        LL = 13699.66
            Iteration Number 000 LL = 13699.66 dPhi = 0.0000 Tr =  0.000,  0.000,  0.000 Ro =  15.000,  0.000,  0.000 Zo =  1.000,  1.000,  1.000 
            Iteration Number 001 LL = 13531.89 dPhi = 0.0585 Tr = -0.041,  0.003,  0.041 Ro =  15.137,  0.020, -0.010 Zo =  1.000,  0.999,  0.998 
            Iteration Number 002 LL = 13384.98 dPhi = 0.0550 Tr = -0.077,  0.007,  0.082 Ro =  15.282,  0.039, -0.021 Zo =  0.999,  0.997,  0.996 
            Iteration Number 003 LL = 13279.62 dPhi = 0.0533 Tr = -0.115,  0.011,  0.120 Ro =  15.433,  0.059, -0.034 Zo =  0.998,  0.995,  0.994 
            Iteration Number 004 LL = 13139.97 dPhi = 0.0531 Tr = -0.154,  0.015,  0.154 Ro =  15.594,  0.082, -0.044 Zo =  0.998,  0.994,  0.993 
            Iteration Number 005 LL = 13046.46 dPhi = 0.0538 Tr = -0.196,  0.019,  0.187 Ro =  15.764,  0.106, -0.058 Zo =  0.997,  0.993,  0.991 
            Iteration Number 006 LL = 12939.69 dPhi = 0.0549 Tr = -0.242,  0.022,  0.217 Ro =  15.938,  0.129, -0.068 Zo =  0.996,  0.992,  0.990 
            Iteration Number 007 LL = 12861.05 dPhi = 0.0573 Tr = -0.292,  0.026,  0.244 Ro =  16.118,  0.151, -0.078 Zo =  0.996,  0.991,  0.989 
            Iteration Number 008 LL = 12733.84 dPhi = 0.0576 Tr = -0.342,  0.030,  0.271 Ro =  16.307,  0.170, -0.090 Zo =  0.995,  0.991,  0.988 
            Iteration Number 009 LL = 12623.39 dPhi = 0.0587 Tr = -0.394,  0.034,  0.296 Ro =  16.499,  0.191, -0.105 Zo =  0.994,  0.990,  0.988 
            Iteration Number 010 LL = 12548.71 dPhi = 0.0623 Tr = -0.450,  0.043,  0.322 Ro =  16.698,  0.211, -0.114 Zo =  0.994,  0.990,  0.987 
            Iteration Number 011 LL = 12469.26 dPhi = 0.0634 Tr = -0.507,  0.053,  0.347 Ro =  16.901,  0.231, -0.125 Zo =  0.993,  0.990,  0.987 
            Iteration Number 012 LL = 12403.03 dPhi = 0.0638 Tr = -0.565,  0.064,  0.369 Ro =  17.103,  0.254, -0.133 Zo =  0.993,  0.990,  0.987 
            Iteration Number 013 LL = 12341.03 dPhi = 0.0642 Tr = -0.623,  0.075,  0.391 Ro =  17.307,  0.279, -0.141 Zo =  0.992,  0.990,  0.987 
            Iteration Number 014 LL = 12287.28 dPhi = 0.0613 Tr = -0.679,  0.087,  0.412 Ro =  17.506,  0.302, -0.145 Zo =  0.992,  0.990,  0.987 
            Iteration Number 015 LL = 12228.45 dPhi = 0.0599 Tr = -0.733,  0.100,  0.434 Ro =  17.708,  0.325, -0.147 Zo =  0.992,  0.990,  0.987 
            Iteration Number 016 LL = 12196.08 dPhi = 0.0552 Tr = -0.782,  0.113,  0.454 Ro =  17.904,  0.345, -0.143 Zo =  0.992,  0.990,  0.987 
            Iteration Number 017 LL = 12155.41 dPhi = 0.0530 Tr = -0.828,  0.127,  0.474 Ro =  18.094,  0.364, -0.140 Zo =  0.992,  0.990,  0.987 
            Iteration Number 018 LL = 12102.12 dPhi = 0.0485 Tr = -0.871,  0.141,  0.492 Ro =  18.273,  0.383, -0.133 Zo =  0.992,  0.991,  0.987 
            Iteration Number 019 LL = 12051.18 dPhi = 0.0449 Tr = -0.909,  0.156,  0.507 Ro =  18.437,  0.402, -0.125 Zo =  0.992,  0.991,  0.988 
            Iteration Number 020 LL = 12021.00 dPhi = 0.0424 Tr = -0.946,  0.171,  0.521 Ro =  18.586,  0.422, -0.116 Zo =  0.992,  0.991,  0.988 
            Iteration Number 021 LL = 12014.26 dPhi = 0.0391 Tr = -0.980,  0.184,  0.535 Ro =  18.720,  0.441, -0.108 Zo =  0.992,  0.992,  0.988 
            Iteration Number 022 LL = 12006.01 dPhi = 0.0350 Tr = -1.009,  0.197,  0.548 Ro =  18.841,  0.461, -0.103 Zo =  0.992,  0.992,  0.988 
            Iteration Number 023 LL = 11982.48 dPhi = 0.0312 Tr = -1.035,  0.208,  0.560 Ro =  18.950,  0.477, -0.096 Zo =  0.992,  0.992,  0.988 
            Iteration Number 024 LL = 11965.45 dPhi = 0.0272 Tr = -1.058,  0.217,  0.572 Ro =  19.047,  0.489, -0.091 Zo =  0.992,  0.993,  0.989 
            Iteration Number 025 LL = 11948.00 dPhi = 0.0243 Tr = -1.079,  0.224,  0.581 Ro =  19.134,  0.500, -0.088 Zo =  0.992,  0.993,  0.989 
            Iteration Number 026 LL = 11934.43 dPhi = 0.0221 Tr = -1.097,  0.231,  0.590 Ro =  19.210,  0.510, -0.084 Zo =  0.992,  0.993,  0.989 
            Iteration Number 027 LL = 11924.29 dPhi = 0.0202 Tr = -1.115,  0.238,  0.597 Ro =  19.277,  0.519, -0.081 Zo =  0.992,  0.994,  0.989 
            Iteration Number 028 LL = 11917.08 dPhi = 0.0175 Tr = -1.130,  0.243,  0.603 Ro =  19.336,  0.527, -0.077 Zo =  0.992,  0.994,  0.989 
            Iteration Number 029 LL = 11912.87 dPhi = 0.0146 Tr = -1.144,  0.246,  0.608 Ro =  19.386,  0.534, -0.074 Zo =  0.992,  0.995,  0.989 
            Iteration Number 030 LL = 11902.09 dPhi = 0.0130 Tr = -1.155,  0.249,  0.612 Ro =  19.428,  0.540, -0.073 Zo =  0.992,  0.995,  0.989 
            Iteration Number 031 LL = 11891.75 dPhi = 0.0117 Tr = -1.166,  0.251,  0.616 Ro =  19.467,  0.546, -0.072 Zo =  0.992,  0.995,  0.989 
            Iteration Number 032 LL = 11882.82 dPhi = 0.0096 Tr = -1.175,  0.253,  0.618 Ro =  19.501,  0.551, -0.072 Zo =  0.992,  0.995,  0.990 
            Iteration Number 033 LL = 11869.45 dPhi = 0.0081 Tr = -1.183,  0.254,  0.619 Ro =  19.528,  0.555, -0.072 Zo =  0.992,  0.996,  0.990 
            Iteration Number 034 LL = 11873.17 dPhi = 0.0064 Tr = -1.189,  0.255,  0.619 Ro =  19.552,  0.558, -0.073 Zo =  0.992,  0.996,  0.990 
            Iteration Number 035 LL = 11883.93 dPhi = 0.0053 Tr = -1.194,  0.256,  0.620 Ro =  19.573,  0.560, -0.074 Zo =  0.992,  0.996,  0.990 
            Iteration Number 036 LL = 11891.06 dPhi = 0.0045 Tr = -1.199,  0.256,  0.620 Ro =  19.590,  0.563, -0.075 Zo =  0.992,  0.996,  0.990 

             -> Converged




.. GENERATED FROM PYTHON SOURCE LINES 150-153

Final transformation
#########################################################
 We can now apply the final transformation

.. GENERATED FROM PYTHON SOURCE LINES 153-157

.. code-block:: default
   :lineno-start: 153

    neReg = spam.DIC.applyPhi(ne, Phi=registration["Phi"])
    print("Translations: {:.4f}, {:.4f}, {:.4f}".format(*registration["transformation"]["t"]))
    print("Rotations   : {:.4f}, {:.4f}, {:.4f}".format(*registration["transformation"]["r"]))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Translations: -1.1988, 0.2564, 0.6200
    Rotations   : 19.5900, 0.5626, -0.0752




.. GENERATED FROM PYTHON SOURCE LINES 158-159

And check the validity of the result with a checkerboard pattern mixing the two images

.. GENERATED FROM PYTHON SOURCE LINES 159-164

.. code-block:: default
   :lineno-start: 159

    checker = spam.helpers.checkerBoard(xr[halfSlice], neReg[halfSlice], n=3)
    plt.figure()
    plt.imshow(checker, cmap="Greys")
    plt.colorbar()




.. image-sg:: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_003.png
   :alt: plot multiModalRegistration
   :srcset: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.colorbar.Colorbar object at 0x7f92aade73a0>



.. GENERATED FROM PYTHON SOURCE LINES 165-167

From the phase diagram a segemntation can also directly be obtained
We can check that phase 1 corresponds to the mortar matrix and phase 2 to the aggregates

.. GENERATED FROM PYTHON SOURCE LINES 167-177

.. code-block:: default
   :lineno-start: 167

    phaseField = registration["phaseField"]
    plt.figure()
    plt.imshow(
        phaseField[halfSlice, :, :],
        vmin=-0.5,
        vmax=nPhases + 0.5,
        cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),
    )
    plt.colorbar(ticks=numpy.arange(0, nPhases + 1))
    plt.show()



.. image-sg:: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_004.png
   :alt: plot multiModalRegistration
   :srcset: /spam_examples/DIC/images/sphx_glr_plot_multiModalRegistration_004.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    /builds/ttk/spam/examples/DIC/plot_multiModalRegistration.py:173: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
      cmap=mpl.cm.get_cmap("Set1_r", nPhases + 1),





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

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


.. _sphx_glr_download_spam_examples_DIC_plot_multiModalRegistration.py:

.. only:: html

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




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

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

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

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


.. only:: html

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

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