atfits_rs/
pixel.rs

1//! Pixel precision and section I/O on the primary image.
2use std::os::raw::{c_int, c_void};
3
4use fitsio::FitsFile;
5use fitsio::hdu::HduInfo;
6use fitsio::images::ImageType;
7use fitsio::sys::{self, LONGLONG};
8
9use crate::error::{AtfitsError, Result, ffcall};
10
11/// Precision a cube is processed in.
12#[derive(Debug, Clone, Copy, PartialEq, Eq)]
13pub enum PixelType {
14    F32,
15    F64,
16}
17
18impl PixelType {
19    /// `-64` maps to f64; `-32` and every integer BITPIX map to f32.
20    pub fn from_bitpix(bitpix: i64) -> Self {
21        if bitpix == -64 {
22            PixelType::F64
23        } else {
24            PixelType::F32
25        }
26    }
27}
28
29/// Element types (`f32`, `f64`) for reading and writing pixel ranges.
30///
31/// Ranges are flat pixel offsets `start..end` into the primary image. Integer
32/// images are scaled by BSCALE/BZERO; BLANK pixels read as NaN, and NaN is
33/// written as BLANK (an error if the header has no BLANK).
34pub trait CubeElem: Copy + Send + Sync + 'static {
35    /// BITPIX of a float image of this type.
36    const BITPIX: i64;
37    /// Read the whole primary image.
38    fn read_full(fptr: &mut FitsFile) -> Result<Vec<Self>>;
39    /// Read pixels `start..end`.
40    fn read_section(fptr: &mut FitsFile, start: usize, end: usize) -> Result<Vec<Self>>;
41    /// Write `data` to pixels `start..end`; `data.len()` must equal `end - start`.
42    fn write_section(fptr: &mut FitsFile, start: usize, end: usize, data: &[Self]) -> Result<()>;
43}
44
45/// Read plane `index` of the primary image, where a plane is `plane_len` pixels.
46pub fn read_plane<T: CubeElem>(
47    fptr: &mut FitsFile,
48    plane_len: usize,
49    index: usize,
50) -> Result<Vec<T>> {
51    let (start, end) = plane_range(plane_len, index)?;
52    T::read_section(fptr, start, end)
53}
54
55/// Write `data` as plane `index` of the primary image.
56pub fn write_plane<T: CubeElem>(fptr: &mut FitsFile, index: usize, data: &[T]) -> Result<()> {
57    let (start, end) = plane_range(data.len(), index)?;
58    T::write_section(fptr, start, end, data)
59}
60
61fn plane_range(plane_len: usize, index: usize) -> Result<(usize, usize)> {
62    plane_len
63        .checked_mul(index)
64        .and_then(|start| Some((start, start.checked_add(plane_len)?)))
65        .ok_or_else(|| AtfitsError::Other(format!("plane {index} overflows usize")))
66}
67
68/// Pixel count of the primary image and whether it stores integers.
69fn primary_image(fptr: &mut FitsFile) -> Result<(usize, bool)> {
70    match fptr.primary_hdu()?.info {
71        HduInfo::ImageInfo { shape, image_type } => {
72            let npix = if shape.is_empty() {
73                0
74            } else {
75                shape.iter().product()
76            };
77            let integer = !matches!(image_type, ImageType::Float | ImageType::Double);
78            Ok((npix, integer))
79        }
80        _ => Err(AtfitsError::Other("primary HDU is not an image".into())),
81    }
82}
83
84fn check_range(start: usize, end: usize, npix: usize) -> Result<()> {
85    if start > end || end > npix {
86        return Err(AtfitsError::Other(format!(
87            "pixel range {start}..{end} is outside the image ({npix} pixels)"
88        )));
89    }
90    Ok(())
91}
92
93macro_rules! impl_cube_elem {
94    ($t:ty, $bitpix:expr, $code:expr) => {
95        impl CubeElem for $t {
96            const BITPIX: i64 = $bitpix;
97
98            fn read_full(fptr: &mut FitsFile) -> Result<Vec<Self>> {
99                let (npix, _) = primary_image(fptr)?;
100                Self::read_section(fptr, 0, npix)
101            }
102
103            fn read_section(fptr: &mut FitsFile, start: usize, end: usize) -> Result<Vec<Self>> {
104                let (npix, _) = primary_image(fptr)?;
105                check_range(start, end, npix)?;
106                let mut out: Vec<$t> = vec![0.0; end - start];
107                let mut nulval = <$t>::NAN;
108                let mut anynul: c_int = 0;
109                ffcall(|s| unsafe {
110                    sys::ffgpv(
111                        fptr.as_raw(),
112                        $code as c_int,
113                        start as LONGLONG + 1,
114                        out.len() as LONGLONG,
115                        &mut nulval as *mut $t as *mut c_void,
116                        out.as_mut_ptr() as *mut c_void,
117                        &mut anynul,
118                        s,
119                    )
120                })?;
121                Ok(out)
122            }
123
124            fn write_section(
125                fptr: &mut FitsFile,
126                start: usize,
127                end: usize,
128                data: &[Self],
129            ) -> Result<()> {
130                let (npix, integer) = primary_image(fptr)?;
131                check_range(start, end, npix)?;
132                if data.len() != end - start {
133                    return Err(AtfitsError::Other(format!(
134                        "pixel range {start}..{end} needs {} values, got {}",
135                        end - start,
136                        data.len()
137                    )));
138                }
139                if !(integer && data.iter().any(|v| v.is_nan())) {
140                    // cfitsio only reads `data`; the C API lacks const.
141                    return ffcall(|s| unsafe {
142                        sys::ffppr(
143                            fptr.as_raw(),
144                            $code as c_int,
145                            start as LONGLONG + 1,
146                            data.len() as LONGLONG,
147                            data.as_ptr() as *mut c_void,
148                            s,
149                        )
150                    });
151                }
152                // cfitsio finds nulls by equality, which never holds for NaN.
153                let mut nulval = <$t>::NEG_INFINITY;
154                let mut buf: Vec<$t> = data
155                    .iter()
156                    .map(|&v| if v.is_nan() { nulval } else { v })
157                    .collect();
158                ffcall(|s| unsafe {
159                    sys::ffppn(
160                        fptr.as_raw(),
161                        $code as c_int,
162                        start as LONGLONG + 1,
163                        buf.len() as LONGLONG,
164                        buf.as_mut_ptr() as *mut c_void,
165                        &mut nulval as *mut $t as *mut c_void,
166                        s,
167                    )
168                })
169            }
170        }
171    };
172}
173impl_cube_elem!(f32, -32, sys::TFLOAT);
174impl_cube_elem!(f64, -64, sys::TDOUBLE);
atfits_rs/pixel.rs

atfits_rs/
pixel.rs
