convolve_rs/
smooth.rs

1//! Convolution plus the flux scaling a brightness unit needs.
2use ndarray::{Array2, CowArray, Ix2};
3
4use crate::beam::Beam;
5use crate::convolve_uv::{ConvolveError, FftFloat, FftPlans, convolve_with_gain};
6
7/// Brightness unit of an image, which decides the flux scaling.
8#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
9pub enum BrightnessUnit {
10    /// Jy/beam: scaled by `Ω_new/Ω_old` so values stay in Jy/beam.
11    #[default]
12    JyPerBeam,
13    /// Kelvin: surface brightness is conserved, so no scaling.
14    Kelvin,
15}
16
17impl BrightnessUnit {
18    /// Classify a FITS `BUNIT` string, or `None` if it is not recognised.
19    ///
20    /// ```
21    /// use convolve_rs::BrightnessUnit;
22    ///
23    /// assert_eq!(BrightnessUnit::parse("Jy/beam"), Some(BrightnessUnit::JyPerBeam));
24    /// assert_eq!(BrightnessUnit::parse(" K "), Some(BrightnessUnit::Kelvin));
25    /// assert_eq!(BrightnessUnit::parse("Jy/pixel"), None);
26    /// ```
27    pub fn parse(bunit: &str) -> Option<Self> {
28        let u = bunit.trim().trim_matches('\'').trim().to_ascii_uppercase();
29        match u.as_str() {
30            "K" | "KELVIN" => Some(BrightnessUnit::Kelvin),
31            "JY/BEAM" | "JY BEAM-1" | "JY/BM" | "JYBEAM" => Some(BrightnessUnit::JyPerBeam),
32            _ => None,
33        }
34    }
35
36    /// Classify a FITS `BUNIT` string, warning and assuming Jy/beam if it is
37    /// not recognised.
38    pub fn from_bunit(bunit: &str) -> Self {
39        Self::parse(bunit).unwrap_or_else(|| {
40            tracing::warn!("{}", Self::unrecognised_message(bunit));
41            BrightnessUnit::JyPerBeam
42        })
43    }
44
45    /// The warning for a `BUNIT` that is not recognised.
46    pub(crate) fn unrecognised_message(bunit: &str) -> String {
47        format!(
48            "Could not determine brightness unit from BUNIT={bunit:?}; assuming \
49             Jy/beam (flux scaling applied). Use a recognised unit (e.g. 'Jy/beam' \
50             or 'K') to silence this warning."
51        )
52    }
53}
54
55/// Smooth `image` from `old_beam` to `new_beam`, scaled for `unit`.
56///
57/// Arguments are as for [`convolve_uv`](crate::convolve_uv()). The result has
58/// the input's precision and shape.
59///
60/// # Errors
61///
62/// As for [`convolve_uv`](crate::convolve_uv()).
63///
64/// ```
65/// use convolve_rs::{Beam, BrightnessUnit, smooth};
66/// use ndarray::Array2;
67///
68/// let old = Beam::from_arcsec(10.0, 10.0, 0.0)?;
69/// let new = Beam::from_arcsec(20.0, 20.0, 0.0)?;
70/// let image = Array2::<f32>::from_elem((64, 64), 1.0);
71/// let dx = 2.5 / 3600.0;
72///
73/// // Jy/beam values scale by the beam-area ratio, 4.
74/// let jy = smooth(&image, &old, &new, dx, dx, BrightnessUnit::JyPerBeam, None)?;
75/// assert!((jy[(32, 32)] - 4.0).abs() < 1e-3);
76///
77/// // Kelvin: a flat image stays flat.
78/// let k = smooth(&image, &old, &new, dx, dx, BrightnessUnit::Kelvin, None)?;
79/// assert!((k[(32, 32)] - 1.0).abs() < 1e-3);
80/// # Ok::<(), Box<dyn std::error::Error>>(())
81/// ```
82pub fn smooth<'a, T: FftFloat>(
83    image: impl Into<CowArray<'a, T, Ix2>>,
84    old_beam: &Beam,
85    new_beam: &Beam,
86    dx_deg: f64,
87    dy_deg: f64,
88    unit: BrightnessUnit,
89    plans: Option<&FftPlans<T>>,
90) -> Result<Array2<T>, ConvolveError> {
91    // The filter's DC gain is g_ratio = √(Ω_new/Ω_old); Jy/beam needs its
92    // square and Kelvin needs 1.
93    let result = convolve_with_gain(
94        image.into(),
95        old_beam,
96        new_beam,
97        (dx_deg, dy_deg),
98        plans,
99        |ratio| match unit {
100            BrightnessUnit::JyPerBeam => ratio * ratio,
101            BrightnessUnit::Kelvin => 1.0,
102        },
103    )?;
104    Ok(result.image)
105}
106
107#[cfg(test)]
108mod tests {
109    use super::*;
110    use ndarray::Array2;
111
112    #[test]
113    fn test_parse_recognised() {
114        for (s, want) in [
115            ("K", BrightnessUnit::Kelvin),
116            (" k ", BrightnessUnit::Kelvin),
117            ("'K'", BrightnessUnit::Kelvin),
118            ("Kelvin", BrightnessUnit::Kelvin),
119            ("Jy/beam", BrightnessUnit::JyPerBeam),
120            ("JY BEAM-1", BrightnessUnit::JyPerBeam),
121        ] {
122            assert_eq!(BrightnessUnit::parse(s), Some(want), "{s}");
123        }
124    }
125
126    #[test]
127    fn test_parse_unrecognised_is_none() {
128        for s in ["", "Jy/pixel", "mJy"] {
129            assert_eq!(BrightnessUnit::parse(s), None);
130        }
131    }
132
133    #[test]
134    fn test_from_bunit_falls_back_to_jy_per_beam() {
135        assert_eq!(BrightnessUnit::from_bunit("K"), BrightnessUnit::Kelvin);
136        assert_eq!(BrightnessUnit::from_bunit("wat"), BrightnessUnit::JyPerBeam);
137    }
138
139    #[test]
140    fn test_kelvin_skips_flux_scaling() {
141        let old = Beam::from_arcsec(10.0, 10.0, 0.0).unwrap();
142        let new = Beam::from_arcsec(20.0, 20.0, 0.0).unwrap();
143        let img = Array2::from_elem((32, 32), 1.0_f32);
144        let dx = 2.5 / 3600.0;
145        let jy = smooth(&img, &old, &new, dx, dx, BrightnessUnit::JyPerBeam, None).unwrap();
146        let k = smooth(&img, &old, &new, dx, dx, BrightnessUnit::Kelvin, None).unwrap();
147        assert!(
148            (jy[(16, 16)] - 4.0).abs() < 1e-3,
149            "Jy/beam {}",
150            jy[(16, 16)]
151        );
152        assert!((k[(16, 16)] - 1.0).abs() < 1e-3, "Kelvin {}", k[(16, 16)]);
153    }
154}
convolve_rs/smooth.rs

convolve_rs/
smooth.rs
