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

.. only:: html

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

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

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

.. _sphx_glr_generated_gallery_imaging_demo.py:


============
Imaging Demo
============

How to create visibility from pixel data and make images.

The example uses ``stixpy`` to obtain STIX pixel data and convert these into visibilities and ``xrayvisim``
to make the images.

Imports

.. GENERATED FROM PYTHON SOURCE LINES 13-36

.. code-block:: Python


    import logging

    import astropy.units as u
    import matplotlib.pyplot as plt
    import numpy as np
    from astropy.coordinates import SkyCoord
    from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
    from sunpy.map import Map, make_fitswcs_header
    from sunpy.time import TimeRange
    from xrayvision.clean import vis_clean
    from xrayvision.imaging import vis_to_image, vis_to_map
    from xrayvision.mem import mem, resistant_mean

    from stixpy.calibration.visibility import calibrate_visibility, create_meta_pixels, create_visibility
    from stixpy.coordinates.frames import STIXImaging
    from stixpy.coordinates.transforms import get_hpc_info
    from stixpy.imaging.em import em
    from stixpy.map.stix import STIXMap  # noqa
    from stixpy.product import Product

    logger = logging.getLogger(__name__)








.. GENERATED FROM PYTHON SOURCE LINES 37-38

Read science file as Product

.. GENERATED FROM PYTHON SOURCE LINES 38-44

.. code-block:: Python


    cpd_sci = Product(
        "http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits"
    )
    cpd_sci





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

 .. code-block:: none


    Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   0%|          | 0.00/48.4M [00:00<?, ?B/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   0%|          | 1.02k/48.4M [00:00<6:07:17, 2.20kB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   0%|          | 32.7k/48.4M [00:00<10:58, 73.5kB/s]  

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   0%|          | 137k/48.4M [00:00<02:41, 298kB/s]  

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   1%|          | 417k/48.4M [00:00<00:53, 895kB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   2%|▏         | 929k/48.4M [00:00<00:24, 1.91MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   4%|▍         | 1.97M/48.4M [00:01<00:11, 3.97MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:   8%|▊         | 4.06M/48.4M [00:01<00:05, 8.12MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:  16%|█▌        | 7.85M/48.4M [00:01<00:02, 15.4MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:  32%|███▏      | 15.7M/48.4M [00:01<00:01, 30.9MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:  51%|█████     | 24.8M/48.4M [00:01<00:00, 45.3MB/s]

    solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits:  82%|████████▏ | 39.9M/48.4M [00:01<00:00, 71.0MB/s]

                                                                                                                                             
    Files Downloaded: 100%|██████████| 1/1 [00:02<00:00,  2.47s/file]
    Files Downloaded: 100%|██████████| 1/1 [00:02<00:00,  2.47s/file]

    CompressedPixelData   <sunpy.time.timerange.TimeRange object at 0x7949a3a43610>
        Start: 2021-09-23 15:20:15
        End:   2021-09-23 15:26:39
        Center:2021-09-23 15:23:27
        Duration:0.004439814814814813 days or
               0.10655555555555551 hours or
               6.393333333333331 minutes or
               383.59999999999985 seconds
        DetectorMasks
        [0...697]: [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]

        PixelMasks
        [0...697]: [['1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1']]

        EnergyEdgeMasks
        [0]: [_,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,_]




.. GENERATED FROM PYTHON SOURCE LINES 45-46

Read background file as Product

.. GENERATED FROM PYTHON SOURCE LINES 46-52

.. code-block:: Python


    cpd_bkg = Product(
        "http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits"
    )
    cpd_bkg





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

 .. code-block:: none


    Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

    solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits:   0%|          | 0.00/121k [00:00<?, ?B/s]

    solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits:   1%|          | 1.02k/121k [00:00<00:52, 2.30kB/s]

    solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits:   7%|▋         | 8.95k/121k [00:00<00:05, 20.6kB/s]

    solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits:  47%|████▋     | 56.9k/121k [00:00<00:00, 132kB/s] 

                                                                                                                                           
    Files Downloaded: 100%|██████████| 1/1 [00:01<00:00,  1.50s/file]
    Files Downloaded: 100%|██████████| 1/1 [00:01<00:00,  1.50s/file]

    CompressedPixelData   <sunpy.time.timerange.TimeRange object at 0x7949a3e9a250>
        Start: 2021-09-23 09:59:23
        End:   2021-09-23 11:35:23
        Center:2021-09-23 10:47:23
        Duration:0.06666666666666665 days or
               1.5999999999999996 hours or
               95.99999999999997 minutes or
               5759.999999999999 seconds
        DetectorMasks
        [0]: [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]

        PixelMasks
        [0]: [['1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1']]

        EnergyEdgeMasks
        [0]: [_,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]




.. GENERATED FROM PYTHON SOURCE LINES 53-54

Set time and energy ranges which will be considered for the science and the background file

.. GENERATED FROM PYTHON SOURCE LINES 54-62

.. code-block:: Python


    time_range_sci = ["2021-09-23T15:20:00", "2021-09-23T15:23:00"]
    time_range_bkg = [
        "2021-09-23T09:00:00",
        "2021-09-23T12:00:00",
    ]  # Set this range larger than the actual observation time
    energy_range = [25, 28] * u.keV








.. GENERATED FROM PYTHON SOURCE LINES 63-64

Create the meta pixel, A, B, C, D for the science and the background data

.. GENERATED FROM PYTHON SOURCE LINES 64-73

.. code-block:: Python


    meta_pixels_sci = create_meta_pixels(
        cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
    )

    meta_pixels_bkg = create_meta_pixels(
        cpd_bkg, time_range=time_range_bkg, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
    )








.. GENERATED FROM PYTHON SOURCE LINES 74-75

Perform background subtraction

.. GENERATED FROM PYTHON SOURCE LINES 75-84

.. code-block:: Python


    meta_pixels_bkg_subtracted = {
        **meta_pixels_sci,
        "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
        "abcd_rate_error_kev_cm": np.sqrt(
            meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
        ),
    }








.. GENERATED FROM PYTHON SOURCE LINES 85-86

Create visibilities from the meta pixels

.. GENERATED FROM PYTHON SOURCE LINES 86-89

.. code-block:: Python


    vis = create_visibility(meta_pixels_bkg_subtracted)








.. GENERATED FROM PYTHON SOURCE LINES 90-91

Obtain the necessary ephemeris data create HPC 0,0 coordinate

.. GENERATED FROM PYTHON SOURCE LINES 91-97

.. code-block:: Python


    vis_tr = TimeRange(vis.meta["time_range"])
    roll, solo_xyz, pointing = get_hpc_info(vis_tr.start, vis_tr.end)
    solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian")
    center_hpc = SkyCoord(0 * u.deg, 0 * u.deg, frame=Helioprojective(obstime=vis_tr.center, observer=solo))





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

 .. code-block:: none


    Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:   0%|          | 0.00/397k [00:00<?, ?B/s]

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:   0%|          | 1.02k/397k [00:00<00:43, 9.12kB/s]

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:   2%|▏         | 8.95k/397k [00:00<00:08, 45.9kB/s]

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:  10%|█         | 40.9k/397k [00:00<00:02, 158kB/s] 

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:  28%|██▊       | 112k/397k [00:00<00:00, 355kB/s] 

    solo_ANC_stix-asp-ephemeris_20210923_V02.fits:  52%|█████▏    | 208k/397k [00:00<00:00, 546kB/s]

                                                                                                    
    Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.09file/s]
    Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.09file/s]




.. GENERATED FROM PYTHON SOURCE LINES 98-101

Calibrate the visibilities

If not given will default to sun center flare location

.. GENERATED FROM PYTHON SOURCE LINES 101-104

.. code-block:: Python


    cal_vis = calibrate_visibility(vis, flare_location=center_hpc)





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

 .. code-block:: none


    Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]
    Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.86file/s]
    Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  2.86file/s]




.. GENERATED FROM PYTHON SOURCE LINES 105-106

Selected detectors 10 to 7

.. GENERATED FROM PYTHON SOURCE LINES 106-111

.. code-block:: Python


    # order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c ....
    isc_10_7 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28]
    idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_7)).ravel()








.. GENERATED FROM PYTHON SOURCE LINES 112-113

Slice the visibilities to detectors 10 - 7

.. GENERATED FROM PYTHON SOURCE LINES 113-116

.. code-block:: Python


    vis10_7 = cal_vis[idx]








.. GENERATED FROM PYTHON SOURCE LINES 117-118

Set up image parameters

.. GENERATED FROM PYTHON SOURCE LINES 118-122

.. code-block:: Python


    imsize = [512, 512] * u.pixel  # number of pixels of the map to reconstruct
    pixel = [10, 10] * u.arcsec / u.pixel  # pixel size in arcsec








.. GENERATED FROM PYTHON SOURCE LINES 123-124

Make a full disk back projection (inverse transform) map

.. GENERATED FROM PYTHON SOURCE LINES 124-127

.. code-block:: Python


    bp_image = vis_to_image(vis10_7, imsize, pixel_size=pixel)








.. GENERATED FROM PYTHON SOURCE LINES 128-129

Obtain the necessary ephemeris data

.. GENERATED FROM PYTHON SOURCE LINES 129-139

.. code-block:: Python


    vis_tr = TimeRange(vis.meta["time_range"])
    roll, solo_xyz, pointing = get_hpc_info(vis_tr.start, vis_tr.end)
    solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian")
    coord_stix = center_hpc.transform_to(STIXImaging(obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo))
    header = make_fitswcs_header(
        bp_image, coord_stix, telescope="STIX", observatory="Solar Orbiter", scale=[10, 10] * u.arcsec / u.pix
    )
    fd_bp_map = Map((bp_image, header))








.. GENERATED FROM PYTHON SOURCE LINES 140-142

Convert the coordinates and make a map in Helioprojective and rotate so "North" is "up"
Center of STIX pointing in HPC

.. GENERATED FROM PYTHON SOURCE LINES 142-148

.. code-block:: Python

    header_hp = make_fitswcs_header(
        bp_image, center_hpc, scale=[10, 10] * u.arcsec / u.pix, rotation_angle=90 * u.deg + roll
    )
    hp_map = Map((bp_image, header_hp))
    hp_map_rotated = hp_map.rotate()








.. GENERATED FROM PYTHON SOURCE LINES 149-150

Plot the both maps

.. GENERATED FROM PYTHON SOURCE LINES 150-164

.. code-block:: Python



    fig = plt.figure(layout="constrained", figsize=(12, 6))
    ax = fig.subplot_mosaic(
        [["stix", "hpc"]], per_subplot_kw={"stix": {"projection": fd_bp_map}, "hpc": {"projection": hp_map_rotated}}
    )
    fd_bp_map.plot(axes=ax["stix"])
    fd_bp_map.draw_limb()
    fd_bp_map.draw_grid()

    hp_map_rotated.plot(axes=ax["hpc"])
    hp_map_rotated.draw_limb()
    hp_map_rotated.draw_grid()




.. image-sg:: /generated/gallery/images/sphx_glr_imaging_demo_001.png
   :alt:  2021-09-23 15:21:38,  2021-09-23 15:21:38
   :srcset: /generated/gallery/images/sphx_glr_imaging_demo_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <CoordinatesMap with 2 world coordinates:

      index aliases    type   unit    wrap   format_unit visible
      ----- ------- --------- ---- --------- ----------- -------
          0     lon longitude  deg 180.0 deg         deg     yes
          1     lat  latitude  deg      None         deg     yes

    >



.. GENERATED FROM PYTHON SOURCE LINES 165-167

Estimate the flare location and plot on top of back projection map. Note the coordinates
are automatically converted from the STIXImaging to Helioprojective

.. GENERATED FROM PYTHON SOURCE LINES 167-176

.. code-block:: Python


    max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel
    # because WCS axes and array are reversed
    max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])

    ax["stix"].plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
    ax["hpc"].plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
    fig.tight_layout()





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

 .. code-block:: none

    /home/docs/checkouts/readthedocs.org/user_builds/stixpy/checkouts/209/examples/imaging_demo.py:174: UserWarning: The figure layout has changed to tight
      fig.tight_layout()




.. GENERATED FROM PYTHON SOURCE LINES 177-178

Use estimated flare location to create more accurate visibilities

.. GENERATED FROM PYTHON SOURCE LINES 178-194

.. code-block:: Python


    meta_pixels_sci = create_meta_pixels(
        cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=max_stix, no_shadowing=True
    )

    meta_pixels_bkg_subtracted = {
        **meta_pixels_sci,
        "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
        "abcd_rate_error_kev_cm": np.sqrt(
            meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
        ),
    }

    vis = create_visibility(meta_pixels_bkg_subtracted)
    cal_vis = calibrate_visibility(vis, flare_location=max_stix)








.. GENERATED FROM PYTHON SOURCE LINES 195-197

Selected detectors 10 to 3
order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c ....

.. GENERATED FROM PYTHON SOURCE LINES 197-200

.. code-block:: Python

    isc_10_3 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28, 15, 27, 31, 6, 30, 2, 25, 5, 23, 7, 29, 1]
    idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_3)).ravel()








.. GENERATED FROM PYTHON SOURCE LINES 201-202

Create an ``xrayvsion`` visibility object

.. GENERATED FROM PYTHON SOURCE LINES 202-206

.. code-block:: Python


    cal_vis.meta["offset"] = max_stix
    vis10_3 = cal_vis[idx]








.. GENERATED FROM PYTHON SOURCE LINES 207-208

Set up image parameters

.. GENERATED FROM PYTHON SOURCE LINES 208-212

.. code-block:: Python


    imsize = [129, 129] * u.pixel  # number of pixels of the map to reconstruct
    pixel = [2, 2] * u.arcsec / u.pixel  # pixel size in arcsec








.. GENERATED FROM PYTHON SOURCE LINES 213-214

Create a back projection image with natural weighting

.. GENERATED FROM PYTHON SOURCE LINES 214-217

.. code-block:: Python


    bp_nat = vis_to_image(vis10_3, imsize, pixel_size=pixel)








.. GENERATED FROM PYTHON SOURCE LINES 218-219

Create a back projection image with uniform weighting

.. GENERATED FROM PYTHON SOURCE LINES 219-222

.. code-block:: Python


    bp_uni = vis_to_image(vis10_3, imsize, pixel_size=pixel, scheme="uniform")








.. GENERATED FROM PYTHON SOURCE LINES 223-224

Create a `sunpy.map.Map` with back projection

.. GENERATED FROM PYTHON SOURCE LINES 224-227

.. code-block:: Python


    bp_map = vis_to_map(vis10_3, imsize, pixel_size=pixel)








.. GENERATED FROM PYTHON SOURCE LINES 228-229

Crete a clean image using the clean algorithm `vis_clean`

.. GENERATED FROM PYTHON SOURCE LINES 229-237

.. code-block:: Python


    niter = 200  # number of iterations
    gain = 0.1  # gain used in each clean iteration
    beam_width = 20.0 * u.arcsec
    clean_map, model_map, resid_map = vis_clean(
        vis10_3, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
    )





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

 .. code-block:: none

    2026-01-28T17:38:08Z INFO xrayvision.clean 124: Iter: 0, strength: 1.9311863170342085, location: (np.int64(62), np.int64(66))
    2026-01-28T17:38:08Z INFO xrayvision.clean 124: Iter: 25, strength: 0.4818291259500767, location: (np.int64(84), np.int64(48))
    2026-01-28T17:38:08Z INFO xrayvision.clean 145: Largest residual negative




.. GENERATED FROM PYTHON SOURCE LINES 238-239

Create a sunpy map for the clean image in Helioprojective

.. GENERATED FROM PYTHON SOURCE LINES 239-249

.. code-block:: Python


    header = make_fitswcs_header(
        clean_map.data,
        max_stix.transform_to(Helioprojective(obstime=vis_tr.center, observer=solo)),
        telescope="STIX",
        observatory="Solar Orbiter",
        scale=pixel,
        rotation_angle=90 * u.deg + roll,
    )








.. GENERATED FROM PYTHON SOURCE LINES 250-251

Crete a map using the MEM GE algorithm `mem`

.. GENERATED FROM PYTHON SOURCE LINES 251-257

.. code-block:: Python


    snr_value, _ = resistant_mean((np.abs(vis10_3.visibilities) / vis10_3.amplitude_uncertainty).flatten(), 3)
    percent_lambda = 2 / (snr_value**2 + 90)
    mem_map = mem(vis10_3, shape=imsize, pixel_size=pixel, percent_lambda=percent_lambda)






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

 .. code-block:: none

    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 0, Chi2: 234.3644777183099
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 25, Chi2: 23.49459412239177
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 50, Chi2: 18.576285880053675
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 75, Chi2: 16.5189804086438
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 100, Chi2: 15.303427427222758
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 125, Chi2: 14.528204741189544
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 150, Chi2: 13.992367150817094
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 175, Chi2: 13.61632070888609
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 200, Chi2: 13.336933875767922
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 225, Chi2: 13.116082735499745
    2026-01-28T17:38:10Z INFO xrayvision.mem 159: Iter: 250, Chi2: 12.930022572631668
    2026-01-28T17:38:10Z INFO xrayvision.mem 510: Iter: 0, Obj function: (241.24949821906276+0j)
    2026-01-28T17:38:12Z INFO xrayvision.mem 510: Iter: 25, Obj function: (29.468841693198318+0j)
    2026-01-28T17:38:14Z INFO xrayvision.mem 510: Iter: 50, Obj function: (28.679569370220314+0j)
    2026-01-28T17:38:16Z INFO xrayvision.mem 510: Iter: 75, Obj function: (28.606898703701226+0j)
    2026-01-28T17:38:19Z INFO xrayvision.mem 510: Iter: 100, Obj function: (28.598019598996785+0j)




.. GENERATED FROM PYTHON SOURCE LINES 258-259

Crete a map using the EM algorithm `EM`

.. GENERATED FROM PYTHON SOURCE LINES 259-277

.. code-block:: Python


    em_map = em(
        meta_pixels_bkg_subtracted["abcd_rate_kev_cm"],
        cal_vis,
        shape=imsize,
        pixel_size=pixel,
        flare_location=max_stix,
        idx=idx,
    )


    clean_map = Map((clean_map.data, header)).rotate()
    bp_map = Map((bp_nat, header)).rotate()
    mem_map = Map((mem_map.data, header)).rotate()
    em_map = Map((em_map, header)).rotate()

    vmax = max([clean_map.data.max(), mem_map.data.max(), em_map.data.max()])





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

 .. code-block:: none

    2026-01-28T17:38:21Z INFO stixpy.imaging.em 180: Iteration: 25, StdDeV: 0.13934706488123141, C-stat: 0.08514307890282385
    2026-01-28T17:38:21Z INFO stixpy.imaging.em 180: Iteration: 50, StdDeV: 0.03072493391043797, C-stat: 0.05479448704103872
    2026-01-28T17:38:21Z INFO stixpy.imaging.em 180: Iteration: 75, StdDeV: 0.0059690156984681594, C-stat: 0.045450543910084615
    2026-01-28T17:38:21Z INFO stixpy.imaging.em 180: Iteration: 100, StdDeV: 0.0017602673077311253, C-stat: 0.042474557084790567
    2026-01-28T17:38:21Z INFO stixpy.imaging.em 180: Iteration: 125, StdDeV: 0.0007918833025777353, C-stat: 0.04111451622144233




.. GENERATED FROM PYTHON SOURCE LINES 278-279

Finally compare the images from each algorithm

.. GENERATED FROM PYTHON SOURCE LINES 279-298

.. code-block:: Python


    fig = plt.figure(figsize=(12, 12))
    ax = fig.subplot_mosaic(
        [
            ["bp", "clean"],
            ["mem", "em"],
        ],
        subplot_kw={"projection": clean_map},
    )
    a = bp_map.plot(axes=ax["bp"])
    ax["bp"].set_title("Back Projection")
    b = clean_map.plot(axes=ax["clean"])
    ax["clean"].set_title("Clean")
    c = mem_map.plot(axes=ax["mem"])
    ax["mem"].set_title("MEM GE")
    d = em_map.plot(axes=ax["em"])
    ax["em"].set_title("EM")
    fig.colorbar(a, ax=ax.values())
    plt.show()



.. image-sg:: /generated/gallery/images/sphx_glr_imaging_demo_002.png
   :alt: Back Projection, Clean, MEM GE, EM
   :srcset: /generated/gallery/images/sphx_glr_imaging_demo_002.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_generated_gallery_imaging_demo.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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