Feat add uv-smooth algorithm

.codespellrc

@@ -13,4 +13,5 @@ ignore-words-list =
     process,
     technik,
     thirdparty,
-    lamba
+    lamba,
+    sav

docs/conf.py

@@ -60,6 +60,7 @@
     "sphinx_changelog",
     "sphinx_design",
     "sphinx_gallery.gen_gallery",
+    "sphinxcontrib.bibtex",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -82,6 +83,8 @@
 autodoc_typehints = "description"
 autoclass_content = "both"
 
+bibtex_bibfiles = ["references.bib"]
+
 # -- Options for intersphinx extension ---------------------------------------
 
 # Example configuration for intersphinx: refer to the Python standard library.

docs/discussions/imaging/index.rst

@@ -0,0 +1,8 @@
+Imaging
+=======
+
+Imaging methods
+
+.. toctree::
+
+   uv-smooth

docs/discussions/index.rst

@@ -2,3 +2,5 @@ Discussions
 ============
 
 .. toctree::
+
+   imaging/

docs/reference/bibliography.rst

@@ -0,0 +1,6 @@
+.. _bibliography_xrayvision:
+
+Bibliography
+------------
+
+.. bibliography::

docs/reference/index.rst

@@ -11,5 +11,7 @@ Reference
    transform
    utils
    visibility
+   uv-smooth
+   bibliography
 
    ../whatsnew/index

docs/reference/uv-smooth.rst

@@ -0,0 +1,8 @@
+.. _uv-smooth:
+
+UV-SMOOTH ('xrayvision.uv_smooth')
+**********************************
+
+The ``uv_smooth`` submodule contains the uv_smooth image reconstruction technique.
+
+.. automodapi:: xrayvision.uv_smooth

docs/references.bib

@@ -0,0 +1,11 @@
+ @article{Massone2009_uv_smooth,
+ title={Hard X-ray Imaging of Solar Flares Using Interpolated Visibilities},
+ volume={703},
+ url={http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2009ApJ...703.2004M&link_type=EJOURNAL},
+ DOI={10.1088/0004-637x/703/2/2004},
+ abstractNote={RHESSI produces solar flare images with the finest angular and spectral resolutions ever achieved at hard X-ray energies. Because this instrument uses indirect, collimator-based imaging techniques, the “native” output of which is in the form of “visibilities” (two-dimensional spatial Fourier components of the image), the development and application of robust, accurate, visibility-based image reconstruction techniques is required. Recognizing that the density of spatial-frequency (u, v) coverage by RHESSI is much sparser than that normally encountered in radio astronomy, we therefore introduce a method for image reconstruction from a relatively sparse distribution of sampled visibilities. The method involves spline interpolation at spatial frequencies less than the largest sampled frequency and the imposition of a positivity constraint on the image to reduce the ringing effects resulting from an unconstrained Fourier transform inversion procedure. Using simulated images consisting both of assumed mathematical forms and of the type of structure typically associated with solar flares, we validate the fidelity, accuracy, and robustness with which the new procedure recovers input images. The method faithfully recovers both single and multiple sources, both compact and extended, over a dynamic range of ~10:1. The performance of the method, which we term as uv_smooth, is compared with other RHESSI image reconstruction algorithms currently in use and its advantages summarized. We also illustrate the application of the method using RHESSI observations of four solar flares.},
+ number={2},
+ journal={The Astrophysical Journal},
+ author={Massone, Anna Maria and Emslie, A Gordon and Hurford, G J and Prato, Marco and Kontar, Eduard P and Piana, Michele},
+ year={2009},
+ pages={2004–2016} }

pyproject.toml

@@ -185,4 +185,3 @@ ignore = [
   directory = "trivial"
   name = "Internal Changes"
   showcontent = true
-

xrayvision/conftest.py

@@ -1,8 +1,6 @@
 # Force MPL to use non-gui backends for testing.
-import matplotlib
-
 try:
-    pass
+    import matplotlib
 except ImportError:
     pass
 else:

xrayvision/imaging.py

@@ -332,7 +332,7 @@ def generate_header(
 @apu.quantity_input()  # type: ignore[untyped-decorator]
 def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibilities:
     r"""
-    Return a Visibilities object created from the map, sampling it at give `u`, `v` coordinates.
+    Return a Visibilities object created from the map, sampling it at given `u`, `v` coordinates.
 
     Parameters
     ----------

xrayvision/simulation/instruments.py

@@ -0,0 +1,50 @@
+import astropy.units as u
+import numpy as np
+
+__all__ = ["rhessi_like_uv_coverage", "stix_like_uv_coverage"]
+
+
+def rhessi_like_uv_coverage():
+    """
+    Creates RHESSI-like u, v coverage coverage distribution.
+    """
+    # Approximate RHESSI spatial frequencies (arcsec^-1) for subcollimators 1-9,
+    # based on typical RHESSI grid pitches of 2.26, 3.92, 6.79, 11.76, 20.36, 35.27, 61.08, 105.8, 183.2 arcsec.
+    isc = np.arange(0, 9)  # index subcollimator (isc) detectors 1-9
+    resolutions = [2.26, 3.92, 6.79, 11.76, 20.36, 35.27, 61.08, 105.8, 183.2] * u.arcsec
+    radii = 1 / resolutions
+
+    # 32 evenly-spaced rotation angles per detector (simulates spacecraft rotation)
+    n_angles = 32
+    angles = np.linspace(0, 2 * np.pi, n_angles, endpoint=False)
+
+    u_comp = np.outer(radii, np.cos(angles))
+    v_comp = np.outer(radii, np.sin(angles))
+
+    isc = isc[:, None] * np.ones(u_comp.shape)
+
+    return {"u": u_comp, "v": v_comp, "isc": isc}
+
+
+def stix_like_uv_coverage():
+    """
+    Creates STIX-like u, v coverage coverage distribution.
+    """
+    resolutions = [7.1, 10.2, 14.6, 20.9, 29.8, 42.7, 61.0, 87.3, 124.9, 178.6] * u.arcsec
+    radii = 1 / resolutions
+
+    # Base orientations for sub-collimators a, b, and c
+    base_angles = np.array([150, 90, 30])
+    steps = np.array([-20, -20, -20])
+    offsets = np.cumsum(np.repeat(steps.reshape(3, 1), 9, axis=1), axis=1)
+    zeros = np.zeros((3, 1))
+    offsets = np.hstack([zeros, offsets])
+    angles_cts = base_angles[:, None] + offsets
+    angles = np.mod(angles_cts, 180) * u.deg
+
+    labels = np.vstack([f"{i}a,{i}b,{i}c".split(",") for i in range(1, 11)]).T
+
+    u_comp = radii * np.cos(angles)
+    v_comp = radii * np.sin(angles)
+
+    return {"u": u_comp, "v": v_comp, "label": labels}

xrayvision/simulation/tests/test_instruments.py

@@ -0,0 +1,11 @@
+from xrayvision.simulation.instruments import rhessi_like_uv_coverage, stix_like_uv_coverage
+
+
+def test_rhessi_like_uv_coverage():
+    out = rhessi_like_uv_coverage()
+    assert out["u"].shape == out["v"].shape == out["isc"].shape == (9, 32)
+
+
+def test_stix_like_uv_coverage():
+    out = stix_like_uv_coverage()
+    assert out["u"].shape == out["v"].shape == out["label"].shape == (3, 10)

xrayvision/tests/test_uv_smooth.py

@@ -0,0 +1,142 @@
+import astropy.units as apu
+import hissw
+import numpy as np
+import pytest
+from astropy.io import fits
+from numpy.testing import assert_allclose
+
+from xrayvision.simulation.instruments import rhessi_like_uv_coverage
+from xrayvision.uv_smooth import uv_smooth, uv_smooth_new
+from xrayvision.visibility import Visibilities, VisMeta
+
+
+# @pytest.mark.skip("needs local files")
+def test_uv_smooth():
+    hdul = fits.open("~/Downloads/hsi_vis_20020221_2357_0054_46tx3e.fits")
+    times = np.unique(hdul[-1].data["TRANGE"], axis=0)
+
+    index = np.argwhere(
+        (np.all(hdul[3].data["TRANGE"] == times[7], axis=1)) & (np.all(hdul[3].data["ERANGE"] == [12.0, 25.0], axis=1))
+    )
+    vis_data = hdul[3].data[index.squeeze()]
+
+    ###############################################################################
+    # Now lets filter by ISC or detector to remove possibly bad data in this case
+    # need to remove ISC 0 and 1.
+    vis_data = vis_data[vis_data["isc"] > 2]
+    vis_data = vis_data[vis_data["obsvis"] != 0 + 0j]
+
+    # vis_sav = readsav("/Users/sm/hsi_hsi_20020221_0006-0007_12-25.sav")
+    # vis_data = vis_sav["vis"]
+    ###############################################################################
+    # Now we can create the visibility object from the filtered visibilities.
+    meta = VisMeta({"isc": vis_data["isc"]})
+
+    vunit = apu.Unit("photon/(cm**2 s)")
+    vis = Visibilities(
+        visibilities=vis_data["obsvis"] * vunit,
+        u=vis_data["u"] / apu.arcsec,
+        v=vis_data["v"] / apu.arcsec,
+        phase_center=vis_data["xyoffset"][0] * apu.arcsec,
+        meta=meta,
+        amplitude_uncertainty=vis_data["sigamp"] * vunit,
+    )
+
+    image_orig, vis_orig, ps_orig = uv_smooth(vis)
+    iamge_new, vis_new, ps_new = uv_smooth_new(vis, shape=128, pixel_size=1)
+    print("here")
+
+
+@pytest.fixture
+def rhessi_like_gaussian_vis():
+    """
+    Synthetic RHESSI-like visibilities of a circular Gaussian source at the origin.
+
+    Samples a circular Gaussian (flux=100 ph/cm²/s, sigma=5 arcsec) at u,v points
+    arranged on annuli at spatial frequencies corresponding to RHESSI subcollimators
+    2-6, with 12 rotation samples per annulus.
+
+    Returns
+    -------
+    vis : Visibilities
+        Synthetic visibilities ready for uv_smooth.
+    flux : float
+        True source flux (ph/cm^2/s).
+    sigma : float
+        True source sigma (arcsec).
+    """
+    uv = rhessi_like_uv_coverage()
+    # detectors 3-7 or isc 2-6
+    u, v, isc = uv["u"][2:7].flatten(), uv["v"][2:7].flatten(), uv["isc"][2:7].flatten()
+
+    flux = 100.0  # ph/cm^2/s
+    sigma = 5.0 * apu.arcsec
+
+    # V(u,v) = flux * exp(-2π²σ²(u²+v²))  [source at origin, so no phase term]
+    vis_vals = flux * np.exp(-2 * np.pi**2 * sigma**2 * (u**2 + v**2))
+
+    vis = Visibilities(
+        visibilities=vis_vals * apu.Unit("ph cm-2 s-1"),
+        u=u,
+        v=v,
+        meta=VisMeta({"isc": isc}),
+    )
+    return vis, flux, sigma
+
+
+def test_uv_smooth_peak_at_origin(rhessi_like_gaussian_vis):
+    vis, flux, sigma = rhessi_like_gaussian_vis
+
+    image, *_ = uv_smooth(vis, niter=50)
+    # For a source centered at (0, 0) the peak should lie at the image center
+    peak_idx = np.unravel_index(np.argmax(image), image.shape)
+    center = np.array(image.shape) // 2
+
+    assert abs(peak_idx[0] - center[0]) <= 3, f"Peak row {peak_idx[0]} too far from image center {center[0]}"
+    assert abs(peak_idx[1] - center[1]) <= 3, f"Peak col {peak_idx[1]} too far from image center {center[1]}"
+
+
+def test_uv_smooth_matches_gaussian(rhessi_like_gaussian_vis):
+    vis, flux, sigma = rhessi_like_gaussian_vis
+    sigma = sigma.value
+
+    image_orig, *og = uv_smooth(vis, niter=50)
+    image_new, *nw = uv_smooth_new(vis, shape=128, uv_pixel_size=0.0005, niter=50)
+    image_auto, *au = uv_smooth_new(vis, shape=128, pixel_size=1.0)
+
+    im_new = image_orig.shape[0]
+
+    # Build the reference Gaussian on the same pixel grid as the uv_smooth output
+    pixel_size = nw[1]  # _uv_smooth_pixel_scale()  # arcsec/pixel
+    coords = (np.arange(im_new) - im_new // 2) * pixel_size  # arcsec
+    xx, yy = np.meshgrid(coords, coords)
+    ref_image = (100.0 / (2 * np.pi * sigma**2)) * np.exp(-0.5 * (xx**2 + yy**2) / sigma**2)
+
+    np.testing.assert_allclose(image_orig, ref_image, atol=0.025)
+    np.testing.assert_allclose(image_new, ref_image, atol=0.025)
+    np.testing.assert_allclose(image_auto, ref_image, atol=0.024)
+
+
+def test_uv_smooth_idl(rhessi_like_gaussian_vis):
+    vis, flux, sigma = rhessi_like_gaussian_vis
+
+    ssw = hissw.Environment(ssw_packages=["hessi"])
+    script = """
+visin = {u:{{u | list }}, v:{{v | list }}, obsvis:{{obsvis | list}}, isc:{{isc|list}}, xyoffset:[0,0], trange:[0,0]}
+uv_smooth, visin, map, reconstructed_map_visibilities=visout
+    """
+    out = ssw.run(
+        script=script,
+        args={
+            "u": vis.u.value.tolist(),
+            "v": vis.v.value.tolist(),
+            "obsvis": vis.visibilities.value.tolist(),
+            "isc": vis.meta["isc"].tolist(),
+        },
+    )
+    image_orig, *info = uv_smooth(vis, niter=50)
+    image_new, *_ = uv_smooth_new(vis, shape=128, uv_pixel_size=0.0005, niter=50)
+    image_auto, *_ = uv_smooth_new(vis)
+    image_idl = out["map"]["data"][0]
+    assert_allclose(image_orig, image_idl, atol=5e-5)
+    assert_allclose(image_new, image_idl, atol=5e-5)

xrayvision/tests/test_visibility.py

@@ -1,3 +1,4 @@
+import astropy.units as apu
 import numpy as np
 import pytest
 from numpy.testing import assert_array_equal
@@ -23,7 +24,8 @@ def test_visibilities():
 
 
 def test_visibility():
-    vis = Visibility(vis=1 * apu.ct, u=1 / apu.deg, v=1 / apu.deg)
+    with pytest.warns(SunpyDeprecationWarning):
+        vis = Visibility(vis=1 * apu.ct, u=1 / apu.deg, v=1 / apu.deg)
     assert vis.vis == 1 * apu.ct
     assert vis.u == 1 / apu.deg
     assert vis.v == 1 / apu.deg

xrayvision/transform.py

@@ -221,8 +221,8 @@ def idft_map(
     """
     m, n = shape
     # python array index in row, column hence y, x
-    y = generate_xy(m, phase_center=phase_center[0], pixel_size=pixel_size[0])
-    x = generate_xy(n, phase_center=phase_center[1], pixel_size=pixel_size[1])
+    y = generate_xy(m, phase_center=phase_center[1], pixel_size=pixel_size[1])  # type: ignore
+    x = generate_xy(n, phase_center=phase_center[0], pixel_size=pixel_size[0])  # type: ignore
 
     x, y = np.meshgrid(x, y)