fits_io/image/
image_data.rs

1use crate::header::BayerPattern;
2use crate::image::Normalizer;
3use image::{ImageBuffer, Luma, Primitive, Rgb};
4use std::error::Error;
5use std::ops::Deref;
6
7/// A single image plane, plus the scaling needed to turn its raw samples into
8/// physical values.
9#[derive(Debug, Clone)]
10pub struct ImageData<T: Primitive> {
11    buffer: ImageBuffer<Luma<T>, Vec<T>>,
12    normalizer: Normalizer,
13    bayer_pattern: Option<BayerPattern>,
14    width: u32,
15    height: u32,
16}
17
18impl<T: Primitive> ImageData<T> {
19    /// Builds an image from its pixels and the normaliser for them.
20    pub fn from_data(
21        width: usize,
22        height: usize,
23        normalizer: Normalizer,
24        bayer_pattern: Option<BayerPattern>,
25        data: Vec<T>,
26    ) -> Result<Self, Box<dyn Error + Send + Sync>> {
27        let buffer = ImageBuffer::<Luma<T>, Vec<T>>::from_raw(width as u32, height as u32, data)
28            .ok_or("Failed to construct image buffer")?;
29        Ok(Self {
30            buffer,
31            normalizer,
32            bayer_pattern,
33            width: width as u32,
34            height: height as u32,
35        })
36    }
37
38    /// The camera bayer pattern or None if the camera is monochrome
39    pub fn bayer_pattern(&self) -> &Option<BayerPattern> {
40        &self.bayer_pattern
41    }
42
43    /// Image width in pixels
44    pub fn width(&self) -> u32 {
45        self.width
46    }
47
48    /// Image height in pixels
49    pub fn height(&self) -> u32 {
50        self.height
51    }
52
53    /// Returns the raw image data
54    pub fn raw(&self) -> &[T] {
55        self.buffer.as_raw()
56    }
57
58    /// The scaling that maps this image's raw samples onto `0.0..=1.0`.
59    pub fn normalizer(&self) -> Normalizer {
60        self.normalizer
61    }
62
63    /// Reads one raw sample as an `f64`, or 0.0 if `(x, y)` lies outside the image.
64    fn raw_at(&self, x: u32, y: u32) -> f64 {
65        if x >= self.width || y >= self.height {
66            return 0.0;
67        }
68
69        self.buffer
70            .get_pixel(x, y)
71            .0
72            .first()
73            .and_then(|sample| sample.to_f64())
74            .unwrap_or(0.0)
75    }
76
77    /// Reads one sample normalised to `0.0..=1.0`.
78    fn normalized_at(&self, x: u32, y: u32) -> f64 {
79        self.normalizer.normalize(self.raw_at(x, y))
80    }
81
82    /// Returns a normalised version of the image, where all values are converted
83    /// into f64 in the range of 0.0 - 1.0.
84    ///
85    /// BZERO and BSCALE are applied first, so the result is a normalised view of
86    /// the *physical* values rather than of the stored samples. See
87    /// [`Normalizer`] for how the black and white points are chosen.
88    pub fn normalized(&self) -> ImageBuffer<Luma<f64>, Vec<f64>> {
89        ImageBuffer::from_fn(self.width, self.height, |x, y| {
90            Luma([self.normalized_at(x, y)])
91        })
92    }
93
94    /// Performs a superpixel demosaic and returns a normalised version.
95    ///
96    /// The superpixel algorithm treats each 2x2 Bayer tile as one output pixel:
97    /// fast, but it halves the resolution in each direction. The two green
98    /// samples in the tile are averaged.
99    ///
100    /// # Errors
101    ///
102    /// Returns an error for a monochrome image, which has no Bayer pattern to
103    /// demosaic.
104    pub fn normalized_superpixel(
105        &self,
106    ) -> Result<ImageBuffer<Rgb<f64>, Vec<f64>>, Box<dyn Error + Send + Sync>> {
107        let bayer_pattern = self
108            .bayer_pattern
109            .ok_or("Can not perform superpixel demosaic on a non rgb image")?;
110        let offsets = bayer_pattern.superpixel_offsets();
111
112        Ok(ImageBuffer::from_fn(
113            self.width / 2,
114            self.height / 2,
115            |x, y| {
116                let (x, y) = (x * 2, y * 2);
117                let sample = |(dx, dy): (u32, u32)| self.normalized_at(x + dx, y + dy);
118
119                let green = (sample(offsets.green[0]) + sample(offsets.green[1])) / 2.0;
120
121                Rgb([sample(offsets.red), green, sample(offsets.blue)])
122            },
123        ))
124    }
125}
126
127impl ImageData<f64> {
128    /// Wraps an existing buffer of physical values, deriving the black and white
129    /// points from the data itself.
130    pub fn from_buffer(buffer: ImageBuffer<Luma<f64>, Vec<f64>>) -> Self {
131        let width = buffer.width();
132        let height = buffer.height();
133        let normalizer = Normalizer::from_samples(0.0, 1.0, buffer.as_raw().iter().copied());
134
135        Self {
136            buffer,
137            normalizer,
138            bayer_pattern: None,
139            width,
140            height,
141        }
142    }
143}
144
145impl<T: Primitive> Deref for ImageData<T> {
146    type Target = ImageBuffer<Luma<T>, Vec<T>>;
147    fn deref(&self) -> &Self::Target {
148        &self.buffer
149    }
150}
151
152#[cfg(test)]
153mod tests {
154    use super::ImageData;
155    use crate::header::BayerPattern;
156    use crate::image::Normalizer;
157
158    /// The unsigned-16-bit encoding astro cameras produce: BITPIX 16, BZERO 32768.
159    fn unsigned_16_bit(width: usize, height: usize, data: Vec<i16>) -> ImageData<i16> {
160        ImageData::from_data(
161            width,
162            height,
163            Normalizer::new(32768.0, 1.0, 0.0, 65535.0),
164            None,
165            data,
166        )
167        .expect("buffer dimensions match the data")
168    }
169
170    #[test]
171    fn normalizing_spreads_samples_over_the_full_range() {
172        let image = unsigned_16_bit(2, 2, vec![i16::MIN, -1, 0, i16::MAX]);
173        let normalized = image.normalized();
174
175        assert_eq!(normalized.get_pixel(0, 0)[0], 0.0);
176        assert_eq!(normalized.get_pixel(1, 1)[0], 1.0);
177        assert!((normalized.get_pixel(1, 0)[0] - 0.5).abs() < 1e-4);
178    }
179
180    #[test]
181    fn normalizing_is_monotonic() {
182        // The inverted `MAX / pixel` form used to make darker samples brighter.
183        let image = unsigned_16_bit(4, 1, vec![i16::MIN, -16384, 16384, i16::MAX]);
184        let normalized = image.normalized();
185
186        let values: Vec<f64> = (0..4).map(|x| normalized.get_pixel(x, 0)[0]).collect();
187
188        for pair in values.windows(2) {
189            assert!(
190                pair[0] < pair[1],
191                "a brighter sample must not normalise darker: {values:?}"
192            );
193        }
194    }
195
196    #[test]
197    fn normalized_values_stay_inside_the_unit_range() {
198        let image = ImageData::from_data(
199            3,
200            1,
201            Normalizer::new(0.0, 1.0, 0.0, 255.0),
202            None,
203            vec![0_u8, 128, 255],
204        )
205        .expect("buffer dimensions match the data");
206
207        for (_, _, pixel) in image.normalized().enumerate_pixels() {
208            assert!(
209                (0.0..=1.0).contains(&pixel[0]),
210                "value out of range: {}",
211                pixel[0]
212            );
213        }
214    }
215
216    const RED: u8 = 200;
217    const GREEN: u8 = 100;
218    const BLUE: u8 = 50;
219
220    /// A 4x4 mosaic tiling the given 2x2 pattern, read left to right and top to
221    /// bottom. The tile is spelled out rather than derived from
222    /// `superpixel_offsets`, so a wrong offset table cannot hide behind it.
223    fn mosaic(pattern: BayerPattern, tile: [u8; 4]) -> ImageData<u8> {
224        let mut data = Vec::with_capacity(16);
225        for y in 0..4 {
226            for x in 0..4 {
227                data.push(tile[(y % 2) * 2 + (x % 2)]);
228            }
229        }
230
231        ImageData::from_data(
232            4,
233            4,
234            Normalizer::new(0.0, 1.0, 0.0, 255.0),
235            Some(pattern),
236            data,
237        )
238        .expect("buffer dimensions match the data")
239    }
240
241    #[test]
242    fn every_bayer_pattern_demosaics_to_the_right_channels() {
243        let cases = [
244            (BayerPattern::RGGB, [RED, GREEN, GREEN, BLUE]),
245            (BayerPattern::BGGR, [BLUE, GREEN, GREEN, RED]),
246            (BayerPattern::GRBG, [GREEN, RED, BLUE, GREEN]),
247            (BayerPattern::GBRG, [GREEN, BLUE, RED, GREEN]),
248        ];
249
250        for (pattern, tile) in cases {
251            let demosaiced = mosaic(pattern, tile)
252                .normalized_superpixel()
253                .expect("a mosaic image can be demosaiced");
254
255            assert_eq!(demosaiced.dimensions(), (2, 2), "{pattern:?}");
256
257            for (x, y, pixel) in demosaiced.enumerate_pixels() {
258                let [red, green, blue] = pixel.0;
259
260                assert!(
261                    (red - RED as f64 / 255.0).abs() < 1e-9,
262                    "{pattern:?} red at ({x}, {y}): {red}"
263                );
264                assert!(
265                    (green - GREEN as f64 / 255.0).abs() < 1e-9,
266                    "{pattern:?} green at ({x}, {y}): {green}"
267                );
268                assert!(
269                    (blue - BLUE as f64 / 255.0).abs() < 1e-9,
270                    "{pattern:?} blue at ({x}, {y}): {blue}"
271                );
272            }
273        }
274    }
275
276    #[test]
277    fn demosaicing_a_monochrome_image_is_an_error() {
278        let image = unsigned_16_bit(2, 2, vec![0, 1, 2, 3]);
279
280        assert!(image.normalized_superpixel().is_err());
281    }
282
283    #[test]
284    fn a_large_zero_offset_does_not_overflow_a_narrow_sample_type() {
285        // The old code did this arithmetic in u8 and i32, so BZERO of 32768
286        // overflowed before it could be applied.
287        let image = ImageData::from_data(
288            2,
289            2,
290            Normalizer::new(32768.0, 1.0, 32768.0, 33023.0),
291            Some(BayerPattern::RGGB),
292            vec![0_u8, 128, 200, 255],
293        )
294        .expect("buffer dimensions match the data");
295
296        let demosaiced = image
297            .normalized_superpixel()
298            .expect("a mosaic image can be demosaiced");
299
300        assert_eq!(demosaiced.dimensions(), (1, 1));
301        for value in demosaiced.get_pixel(0, 0).0 {
302            assert!((0.0..=1.0).contains(&value), "value out of range: {value}");
303        }
304    }
305
306    #[test]
307    fn odd_dimensions_do_not_read_outside_the_image() {
308        let image = ImageData::from_data(
309            3,
310            3,
311            Normalizer::new(0.0, 1.0, 0.0, 255.0),
312            Some(BayerPattern::RGGB),
313            vec![0_u8; 9],
314        )
315        .expect("buffer dimensions match the data");
316
317        let demosaiced = image
318            .normalized_superpixel()
319            .expect("a mosaic image can be demosaiced");
320
321        assert_eq!(demosaiced.dimensions(), (1, 1));
322    }
323
324    #[test]
325    fn from_buffer_derives_its_range_from_the_data() {
326        let buffer = image::ImageBuffer::from_raw(2, 2, vec![10.0_f64, 20.0, 30.0, 40.0])
327            .expect("buffer dimensions match the data");
328        let image = ImageData::from_buffer(buffer);
329
330        let normalized = image.normalized();
331
332        // The old code set BSCALE to 0.0, which flattened every image to nothing.
333        assert_eq!(normalized.get_pixel(0, 0)[0], 0.0);
334        assert_eq!(normalized.get_pixel(1, 1)[0], 1.0);
335        assert!((normalized.get_pixel(1, 0)[0] - 1.0 / 3.0).abs() < 1e-9);
336    }
337}
fits_io/image/image_data.rs

fits_io/image/
image_data.rs
