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

.. only:: html

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

        New to DeepInverse? Get started with the basics with the
        :ref:`5 minute quickstart tutorial <sphx_glr_auto_examples_basics_demo_quickstart.py>`.

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

.. _sphx_glr_auto_examples_physics_demo_remote_sensing.py:


Remote sensing with satellite images
====================================

In this example we demonstrate remote sensing inverse problems for multispectral satellite imaging.

These have important applications for image restoration in environmental monitoring, urban planning, disaster recovery etc.

We will demonstrate pan-sharpening, i.e., recovering high-resolution multispectral images from measurement pairs of
low-resolution multispectral images and high-resolution panchromatic (single-band) images with the forward
operator :class:`deepinv.physics.Pansharpen`.

We will also demonstrate other inverse problems including compressive spectral imaging and hyperspectral unmixing.

We provide a convenient satellite image dataset for pan-sharpening :class:`deepinv.datasets.NBUDataset` provided in the paper :footcite:t:`meng2020large`.
which includes data from several satellites such as WorldView satellites.

.. tip::

    For remote sensing experiments, DeepInverse provides the following classes:

    - :class:`Pan-sharpening <deepinv.physics.Pansharpen>`
    - :class:`Compressive spectral imaging <deepinv.physics.CompressiveSpectralImaging>`
    - :class:`Hyperspectral unmixing <deepinv.physics.HyperSpectralUnmixing>`
    - :class:`Super resolution <deepinv.physics.Downsampling>`
    - :class:`Satellite imagery dataset <deepinv.datasets.NBUDataset>`
    - Metrics for multispectral data: :class:`QNR <deepinv.loss.metric.QNR>`, :class:`SpectralAngleMapper <deepinv.loss.metric.SpectralAngleMapper>`, :class:`ERGAS <deepinv.loss.metric.ERGAS>`

.. GENERATED FROM PYTHON SOURCE LINES 33-37

.. code-block:: Python

    import deepinv as dinv
    import torch

    device = dinv.utils.get_freer_gpu() if torch.cuda.is_available() else "cpu"







.. GENERATED FROM PYTHON SOURCE LINES 38-55

Load raw pan-sharpening measurements
------------------------------------
The dataset includes raw pansharpening measurements
containing ``(MS, PAN)`` where ``MS`` are the low-res (4-band) multispectral and ``PAN`` are the high-res
panchromatic images. Note there are no ground truth images!

.. note::

  The pan-sharpening measurements are provided as a :class:`deepinv.utils.TensorList`, since
  the pan-sharpening physics :class:`deepinv.physics.Pansharpen` is a stacked physics combining
  :class:`deepinv.physics.Downsampling` and :class:`deepinv.physics.Decolorize`.
  See the User Guide :ref:`physics_combining` for more information.

Note, for plotting purposes we only plot the first 3 bands (RGB).

Note also that the linear adjoint must assume the unknown spectral response function (SRF).


.. GENERATED FROM PYTHON SOURCE LINES 55-82

.. code-block:: Python


    DATA_DIR = dinv.utils.get_data_home()
    dataset = dinv.datasets.NBUDataset(DATA_DIR, return_pan=True, download=True)

    y = dataset[0].unsqueeze(0).to(device)  # MS (1,4,256,256), PAN (1,1,1024,1024)

    physics = dinv.physics.Pansharpen((4, 1024, 1024), factor=4, device=device)

    # Pansharpen with classical Brovey method
    x_hat = physics.A_dagger(y)  # shape (1,4,1024,1024)

    dinv.utils.plot(
        [
            y[0][:, :3],
            y[1],  # Note this will be interpolated to match high-res image size
            x_hat[:, :3],
            physics.A_adjoint(y)[:, :3],
        ],
        titles=[
            "Input MS",
            "Input PAN",
            "Pseudo-inverse \n using \n Brovey method",
            "Linear adjoint",
        ],
        dpi=1200,
    )




.. image-sg:: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_001.png
   :alt: Input MS, Input PAN, Pseudo-inverse   using   Brovey method, Linear adjoint
   :srcset: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    Downloading datasets/nbu/gaofen-1.zip

      0%|          | 0/7914941 [00:00<?, ?it/s]
     43%|████▎     | 3.25M/7.55M [00:00<00:00, 34.0MB/s]
    100%|██████████| 7.55M/7.55M [00:00<00:00, 57.3MB/s]

    Extracting:   0%|          | 0/12 [00:00<?, ?it/s]
    Extracting: 100%|██████████| 12/12 [00:00<00:00, 915.60it/s]
    Dataset has been successfully downloaded.




.. GENERATED FROM PYTHON SOURCE LINES 83-85

Evaluate performance - note we can only use QNR as we have no GT


.. GENERATED FROM PYTHON SOURCE LINES 85-90

.. code-block:: Python


    qnr = dinv.metric.QNR()
    print(qnr(x_net=x_hat, x=None, y=y, physics=physics))






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

 .. code-block:: none

    tensor([0.4133])




.. GENERATED FROM PYTHON SOURCE LINES 91-96

Simulate remote-sensing measurements
------------------------------------
We can also simulate measurements from various remote sensing inverse problems so that we have pairs of
measurements and ground truth. Now, the dataset loads ground truth images ``x``:


.. GENERATED FROM PYTHON SOURCE LINES 96-101

.. code-block:: Python


    dataset = dinv.datasets.NBUDataset(DATA_DIR, return_pan=False)

    x = dataset[0].unsqueeze(0).to(device)  # just MS of shape 1,4,256,256








.. GENERATED FROM PYTHON SOURCE LINES 102-105

For **compressive spectral imaging**, we use the coded-aperture snapshot spectral imaging (CASSI) model,
which is a popular hyperspectral imaging method. See :class:`deepinv.physics.CompressiveSpectralImaging`


.. GENERATED FROM PYTHON SOURCE LINES 105-110

.. code-block:: Python


    physics = dinv.physics.CompressiveSpectralImaging(x.shape[1:], mode="sd", device=device)
    y = physics(x)  # 1,1,256,256
    dinv.utils.plot([x[:, :3], y], titles=["Image x", "CASSI meas. y"])




.. image-sg:: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_002.png
   :alt: Image x, CASSI meas. y
   :srcset: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 111-115

For **hyperspectral unmixing**, our images are the measurements and we seek to recover abundances
given the endmember matrix in the linear mixing model.
In this toy example, we perform unmixing with 2 endmembers: one purely yellow and one purely blue.


.. GENERATED FROM PYTHON SOURCE LINES 115-125

.. code-block:: Python


    physics = dinv.physics.HyperSpectralUnmixing(
        M=torch.tensor([[0.5, 0.5, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]), device=device
    )
    abundance = physics.A_adjoint(x)  # 1,2,256,256
    dinv.utils.plot(
        [x[:, :3], abundance[:, [0]], abundance[:, [1]]],
        titles=["Mixed image", "Yellow abudance", "Blue abundance"],
    )




.. image-sg:: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_003.png
   :alt: Mixed image, Yellow abudance, Blue abundance
   :srcset: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 126-129

For the **pansharpening** physics, we assume a flat spectral response function,
but this can also be jointly learned. We simulate Gaussian noise on the panchromatic images.


.. GENERATED FROM PYTHON SOURCE LINES 129-137

.. code-block:: Python


    physics = dinv.physics.Pansharpen((4, 256, 256), factor=4, srf="flat", device=device)

    y = physics(x)

    # Pansharpen with classical Brovey method
    x_hat = physics.A_dagger(y)








.. GENERATED FROM PYTHON SOURCE LINES 138-154

Solving pan-sharpening with neural networks
-------------------------------------------
The pan-sharpening physics is compatible with the rest of the DeepInverse library
so we can solve the inverse problem using any method provided in the library.
For example, we use here the PanNet :footcite:t:`yang2017pannet` model.

This model can be trained using losses such as supervised learning using :class:`deepinv.loss.SupLoss`
or self-supervised learning using Equivariant Imaging :class:`deepinv.loss.EILoss`, which was applied to
pan-sharpening in :footcite:t:`wang2024perspective`.

For evaluation, we use the standard full-reference metrics (ERGAS, SAM) and no-reference (QNR).

.. note::

  This is a tiny example using 5 images. We demonstrate training for 1 epoch for speed, but you can train from scratch using 50 epochs.


.. GENERATED FROM PYTHON SOURCE LINES 154-158

.. code-block:: Python


    model = dinv.models.PanNet(hrms_shape=(4, 256, 256), device=device)
    x_net = model(y, physics)








.. GENERATED FROM PYTHON SOURCE LINES 159-163

Example training loss using measurement consistency on the multispectral images
and Stein's Unbiased Risk Estimate on the panchromatic images.
For metrics, we use standard full-reference and no-reference multispectral pan-sharpening metrics,
since ground-truth is now available.

.. GENERATED FROM PYTHON SOURCE LINES 163-173

.. code-block:: Python


    loss = dinv.loss.StackedPhysicsLoss(
        [dinv.loss.MCLoss(), dinv.loss.SureGaussianLoss(0.05)]
    )

    sam = dinv.metric.distortion.SpectralAngleMapper()
    ergas = dinv.metric.distortion.ERGAS(factor=4)
    qnr = dinv.metric.QNR()
    print(sam(x_hat, x), ergas(x_hat, x), qnr(x_hat, x=None, y=y, physics=physics))





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

 .. code-block:: none

    tensor([0.0545]) tensor([4.1538]) tensor([0.4737])




.. GENERATED FROM PYTHON SOURCE LINES 174-176

For training, we first load optimizer and pretrained model,
then train using the deepinv Trainer.

.. GENERATED FROM PYTHON SOURCE LINES 176-210

.. code-block:: Python


    optimizer = torch.optim.Adam(model.parameters())

    from deepinv.models.utils import get_weights_url

    file_name = "demo_nbu_pansharpen.pth"
    url = get_weights_url(model_name="demo", file_name=file_name)
    ckpt = torch.hub.load_state_dict_from_url(
        url, map_location=lambda storage, loc: storage, file_name=file_name
    )
    model.load_state_dict(ckpt["state_dict"])
    optimizer.load_state_dict(ckpt["optimizer"])

    from torch.utils.data import DataLoader

    trainer = dinv.Trainer(
        model=model,
        physics=physics,
        optimizer=optimizer,
        losses=loss,
        metrics=[sam, ergas],
        train_dataloader=DataLoader(dataset),
        epochs=1,
        online_measurements=True,
        plot_images=False,
        compare_no_learning=True,
        no_learning_method="A_dagger",
        show_progress_bar=False,
        device=device,
    )

    trainer.train()
    trainer.test(DataLoader(dataset))





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

 .. code-block:: none

    Downloading: "https://huggingface.co/deepinv/demo/resolve/main/demo_nbu_pansharpen.pth?download=true" to /home/runner/.cache/torch/hub/checkpoints/demo_nbu_pansharpen.pth

      0%|          | 0.00/955k [00:00<?, ?B/s]
    100%|██████████| 955k/955k [00:00<00:00, 13.5MB/s]
    The model has 77124 trainable parameters
    Train epoch 0: TotalLoss=0.002, SpectralAngleMapper=0.28, ERGAS=23.228
    Eval epoch 0: SpectralAngleMapper=0.279, SpectralAngleMapper no learning=0.039, ERGAS=24.093, ERGAS no learning=4.351
    Test results:
    SpectralAngleMapper no learning: 0.039 +- 0.011
    SpectralAngleMapper: 0.279 +- 0.012
    ERGAS no learning: 4.351 +- 1.179
    ERGAS: 24.093 +- 9.487

    {'SpectralAngleMapper no learning': np.float64(0.03946961760520935), 'SpectralAngleMapper no learning_std': np.float64(0.011352569126202245), 'SpectralAngleMapper': np.float64(0.2789239645004272), 'SpectralAngleMapper_std': np.float64(0.011772674749104428), 'ERGAS no learning': np.float64(4.351437377929687), 'ERGAS no learning_std': np.float64(1.1789995963538327), 'ERGAS': np.float64(24.09258117675781), 'ERGAS_std': np.float64(9.487485230098399)}



.. GENERATED FROM PYTHON SOURCE LINES 211-212

Plot sample results:

.. GENERATED FROM PYTHON SOURCE LINES 212-223

.. code-block:: Python

    dinv.utils.plot(
        [
            x[:, :3],
            y[0][:, :3],
            y[1],
            x_hat[:, :3],
            x_net[:, :3],
        ],
        titles=["x HRMS", "y LRMS", "y PAN", "Estimate (classical)", "Estimate (PanNet)"],
    )




.. image-sg:: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_004.png
   :alt: x HRMS, y LRMS, y PAN, Estimate (classical), Estimate (PanNet)
   :srcset: /auto_examples/physics/images/sphx_glr_demo_remote_sensing_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 224-227

:References:

.. footbibliography::


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

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


.. _sphx_glr_download_auto_examples_physics_demo_remote_sensing.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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