diff --git a/.gitignore b/.gitignore index f6a336b..e8512f0 100644 --- a/.gitignore +++ b/.gitignore @@ -23,4 +23,4 @@ Makefile **/*.rs.bk **/*.swp Licence_notice.txt - +SPMIEEE-proof.pdf diff --git a/Cargo.toml b/Cargo.toml index fbbd36c..b5d23ab 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,9 +31,17 @@ doc = true [dependencies] katex-doc = "0.1.0" +log = "0.4" base64 = "0.21" # Added for Compressed MOC num = "0.4" # Added for MOC num-traits = "0.2" # Added for MOC +# Skymaps +byteorder = "1.4.3" +thiserror = "1.0" +colorous = "1.0" +mapproj = "0.3" +png = "0.17" +flate2 = "1.0" # Compression/decompression [dev-dependencies] criterion = "0.4" diff --git a/src/nested/map/astrometry/gal.rs b/src/nested/map/astrometry/gal.rs new file mode 100644 index 0000000..fbc7a13 --- /dev/null +++ b/src/nested/map/astrometry/gal.rs @@ -0,0 +1,172 @@ +use super::math::*; + +/// Computed using the ZYZ Euler transformation from the FK4 B1950 frame in which, based on +/// Blaauw (1958) and Liu (2011, A&A536, A16): +/// * North Galactic Pole FK4 B1950 longitude: 12h49m= +192.25 deg (1st rotation: z-axis of 192.25 deg) +/// * North Galactic Pole FK4 B1950 latitude : +27.4 deg (2nd rotation: x'-axis by 90 - 27.4 = 62.6 deg) +/// * Position angle, at the Galactic pole, from the equatorial pole, +/// of the great circle passing through the Galactic center: +123 deg. +/// Noting NEP (North Equatorial Pole), NGP (North Galactic Pole) and GC (Galactic Center), it is +/// the angle NEP-NGP-GP. The rotation from NGP so set the origin of the frame at the GC is thus +/// (3rd rotation: z''-axis by 180 - 123 = 57 deg). +/// # Note +/// * The additional precision comes from `bc/src/galactic.bc` and have been cross-checked +/// with a code using 80bits floats from F. Ochsenbein +const FK4B19502GAL: M3x3 = M3x3( + // Matrix from BC + XYZt( + -0.066_988_739_415_150_72, + -0.872_755_765_851_992_7, + -0.483_538_914_632_184_24, + ), + XYZt( + 0.492_728_466_075_323_6, + -0.450_346_958_019_961_33, + 0.744_584_633_283_031, + ), + XYZt( + -0.867_600_811_151_434_8, + -0.188_374_601_722_920_45, + 0.460_199_784_783_851_65, + ), +); + +/// Computed using the ZYZ Euler transformation from the **truncated** (ESA convention) coordinates +/// of the North Galactic Pole (NGP) plus the position angle of the NGP with respect to the +/// equatorial pole. See: +/// * Section 7.1 of Murray (1989), +/// * Eq. (1.5.9) and (1.5.10) of Hipparcos and Tycho Vol1, +/// * Eq. (2) of Liu (2011, A&A 536, 102). +/// * Section 2 from Liu (2011, A&A536, A16) +/// * Eq. (31) in Fox document +/// That is: +/// * NGP longitude in FK5 J2000: 12h51m26.2755s = +192.85948 deg (1st rotation: z-axis of 192.85948 deg) +/// * NGP latitude in FK5 J2000: +27.12825 deg (2nd rotation: x'-axis by 90 - 27.12825 = 62.87175 deg) +/// * Position angle, at the Galactic pole, from the equatorial pole, +/// of the great circle passing through the Galactic center: +122.93192 deg. +/// Noting NEP (North Equatorial Pole), NGP (North Galactic Pole) and GC (Galactic Center), it is +/// the angle NEP-NGP-GP. The rotation from NGP so set the origin of the frame at the GC is thus +/// (3rd rotation: z''-axis by 180 - 122.93192 = 57.06808 deg). +/// # Notes +/// * follows the ESA convention defined for Hipparcos, Tycho and Gaia +/// + no difference made between FK5 and ICRS +/// + based on truncated values computed from the conversion of the galactic center coordinates +/// (and position angle) from FK4 to FK5 +/// * The additional precision comes from `bc/src/galactic.bc` and have been cross-checked +/// with SOFA (`icrs2g.c` routine). +const FK5J20002GAL_ESA: M3x3 = M3x3( + XYZt( + -0.054_875_560_416_215_368, + -0.873_437_090_234_885, + -0.483_835_015_548_713_2, + ), + XYZt( + 0.494_109_427_875_583_65, + -0.444_829_629_960_011_2, + 0.746_982_244_497_218_8, + ), + XYZt( + -0.867_666_149_019_004_7, + -0.198_076_373_431_201_52, + 0.455_983_776_175_066_9, + ), +); + +// So far, no reliable and MIT high precision crate in rust: +// * rug is LGPG +// * sin of astro-float does not seems reliable enough +// Thus I decided to write a 'bc' script to compute the matrix coefficient + +/// Accounts from the difference between FK5 and ICRS in Mignard (2002). +/// * Eq. (18) from Liu (2011) A&A536, A16 +/// * Eq. (30) in Fox document +const ICRS2GAL: M3x3 = M3x3( + XYZt( + -0.054_875_657_712_619_678_1, + -0.873_437_051_955_779_129_8, + -0.483_835_073_616_418_380_3, + ), + XYZt( + 0.494_109_437_197_107_641_2, + -0.444_829_721_222_053_763_5, + 0.746_982_183_839_845_094_133, + ), + XYZt( + -0.867_666_137_557_162_561_5, + -0.198_076_337_275_070_594_6, + 0.455_983_813_691_152_347_6, + ), +); + +/// The Galactic coordinate system is barycentric, its xy-plane is the galactic plane (it is assumed +/// the the barycenter of the solar system is in the galactic plane, which is not correct) +/// and its x-axis points toward the galactic center. +/// It is defined from Equatorial frames (such a FK4, FK5 or ICRS) by the position of the +/// North Galactic Pole (NGC) in the equatorial frame, plus the position angle formed by +/// the North Equatorial Pole (NEP), the NGC and the great circle passing through the NEP and +/// the Galactic Center (GC). +pub struct Galactic { + /// Rotation matrix + m: M3x3, +} + +impl Galactic { + /// New Galactic frame for the conversion with the FK4 B19500 frame. + pub const fn new_for_fk4_b1950() -> Self { + Galactic { m: FK4B19502GAL } + } + + /// New Galactic frame for the conversion with the FK5 J2000 frame, or the ICRS frame + /// following the ESA (Hipparcos, Tycho, Gaia) conventions. + /// Its also correspond the SOFA icrs <-> gal conversion. + pub const fn new_for_fk5_j2000_and_icrs() -> Self { + Galactic { + m: FK5J20002GAL_ESA, + } + } + + /// New Galactic frame for the conversion with the ICRS frame accounting for the Mignard (2002) + /// offsets between FK5 and ICRS. + pub const fn new_for_icrs_including_fk5_icrs_offsets_from_mignard2002() -> Self { + Galactic { m: ICRS2GAL } + } + + /// Convert the given coordinate from Equatorial (FK4 B1950, FK5 J2000 or ICRS) to Galactic. + pub fn coo_eq2gal(&self, pos: &Coo) -> Coo { + Coo::from_xyz(&self.xyz_eq2gal(&pos.xyz())) + } + + /// Convert the given coordinate from Galactic to Equatorial (FK4 B1950, FK5 J2000 or ICRS). + pub fn coo_gal2eq(&self, pos: &Coo) -> Coo { + Coo::from_xyz(&self.xyz_gal2eq(&pos.xyz())) + } + + /// Convert the given Euclidean position from Equatorial (FK4 B1950, FK5 J2000 or ICRS) to + /// Galactic. + pub fn xyz_eq2gal(&self, pos: &XYZ) -> XYZ { + self.m.rotate(pos) + } + + /// Convert the given Euclidean position from Galactic to Equatorial + /// (FK4 B1950, FK5 J2000 or ICRS). + pub fn xyz_gal2eq(&self, pos: &XYZ) -> XYZ { + self.m.unrotate(pos) + } +} + +#[cfg(test)] +mod tests { + use crate::nested::map::astrometry::{gal::Galactic, math::Coo}; + + #[test] + fn test_gal2eq() { + let lon = 0.0; + let lat = 0.0; + let gal_center = Galactic::new_for_fk5_j2000_and_icrs().coo_gal2eq(&Coo::new(lon, lat)); + println!( + "Gal Center: ({:.6}, {:.6})", + gal_center.lon.to_degrees(), + gal_center.lat.to_degrees() + ) + } +} diff --git a/src/nested/map/astrometry/math.rs b/src/nested/map/astrometry/math.rs new file mode 100644 index 0000000..6e66aa6 --- /dev/null +++ b/src/nested/map/astrometry/math.rs @@ -0,0 +1,583 @@ +//! This module contains elements of linear algebra and spherical geometry needed to solve +//! astrometrical problems. +//! +use std::ops::{Add, Sub}; + +use num_traits::FloatConst; + +use crate::Customf64; + +pub trait HasLocalRotMatrix { + fn local_rot_matrix(&self) -> M3x3; +} + +#[derive(PartialEq, Debug)] +pub struct Coo { + /// longitude in `[0, 2pi[`, in `rad` + pub lon: f64, // usual units: deg + /// latitude in `[-pi/2, pi/2]`, in `rad` + pub lat: f64, // usual units: deg +} + +impl Coo { + pub fn new(lon: f64, lat: f64) -> Self { + Self { lon, lat } + } + + pub fn from_deg(lon: f64, lat: f64) -> Self { + Coo { + lon: lon.to_radians(), + lat: lat.to_radians(), + } + } + + pub fn from_xyz(xyz: &XYZ) -> Self { + Self::from(xyz.0, xyz.1, xyz.2) + } + + /// # Output + /// * `(lon, lat)`, in degrees + pub fn to_deg(&self) -> (f64, f64) { + (self.lon.to_degrees(), self.lat.to_degrees()) + } + + /// Create spherical coordinates from cartesian coordinates. + /// # Remark + /// The cartesian vector do not have to be normalized. + pub fn from(x: f64, y: f64, z: f64) -> Self { + const TWICE_PI: f64 = 2.0 * std::f64::consts::PI; + let mut lon = y.atan2(x); + if lon < 0.0_f64 { + lon += TWICE_PI; + } else if lon == TWICE_PI { + lon = 0.0; + } + let lat = z.atan2((x.pow2() + y.pow2()).sqrt()); + Coo { lon, lat } + } + + /// Provides the vector $\vec{r0} = (x, y, z)$, see Eq. 1.2.11 in + /// [The Hipparcos and Tycho Catalogues](https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf) + /// in radians. + /// Remark: the vector is normalized (unit vector) + pub fn xyz(&self) -> XYZ { + XYZ::from_coo(self.lon, self.lat) + } +} + +impl HasLocalRotMatrix for Coo { + /// Components of the 3x3 rotation matrix transforming a vector into the + /// reference frame to a vector into the local frame (i.e. the frame in + /// which the position of the projection center is (1, 0, 0). + /// + /// The 3 transposed vector composing this matrix correspond to + /// $\vec{r_0}$, $\vec{p_0}$, $\vec{q_0}$ of Eq. 1.2.16 in + /// [The Hipparcos and Tycho Catalogues](https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf) + /// + /// # Remark 1: + /// * $r_{22} = \cos(lon)$ + /// * $r_{21} = -\sin(lon)$ + /// * $r_{33} = \cos(lat)$ + /// * $r_{13} = \sin(lat)$ + /// + /// # Remark 2 + /// If you need to compute the `xyz` coordinates of the vector then you + /// may choose to compute the matrix from the Euclidean coordinates not + /// to have to call twice the trigonometric functions (perform first a + /// bench to see if the compiler is not optimizing by caching the trigonometric + /// function results the code with Euclidean coordinate is more complex since + /// division by 0 have to be tested). + fn local_rot_matrix(&self) -> M3x3 { + let (sin_lon, cos_lon) = self.lon.sin_cos(); + let (sin_lat, cos_lat) = self.lat.sin_cos(); + M3x3::new( + XYZt::new(cos_lat * cos_lon, cos_lat * sin_lon, sin_lat), // r0 + XYZt::new(-sin_lon, cos_lon, 0.0), // p0 + XYZt::new(-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat), // q0 + ) + } +} + +/// Can be called on array ref +/// ```rust +/// use cdshealpix::nested::map::astrometry::math::scalar; +/// +/// let a1: [f64; 3] = [1.0, 2.0, 3.0]; +/// let a2: [f64; 3] = [2.0, 3.0, 4.0]; +/// assert_eq!(scalar(&a1, &a2), 20.0); +/// ``` +pub fn scalar<'a, A, B>(lhs: A, rhs: B) -> f64 +where + A: IntoIterator, + B: IntoIterator, +{ + lhs + .into_iter() + .zip(rhs.into_iter()) + .map(|(a, b)| a * b) + .sum() +} + +pub trait HasXYZ { + fn x(&self) -> f64; + fn y(&self) -> f64; + fn z(&self) -> f64; + + fn to_coo(&self) -> Coo { + Coo::from(self.x(), self.y(), self.z()) + } + + fn norm(&self) -> f64 { + self.squared_norm().sqrt() + } + + fn squared_norm(&self) -> f64 { + self.x().pow2() + self.y().pow2() + self.z().pow2() + } + + fn scalar(&self, rhs: &T) -> f64 { + self.x() * rhs.x() + self.y() * rhs.y() + self.z() * rhs.z() + } + + /// Scalar product of the y and z coordinates ignoring the x coordinate). + /// Used when we know in advance that either `self.x()` or `rhs.x()` equals 0 + fn scalar_yz(&self, rhs: &T) -> f64 { + self.y() * rhs.y() + self.z() * rhs.z() + } +} + +/// 3 dimentional vector +#[derive(Debug, Clone)] +pub struct XYZ(pub f64, pub f64, pub f64); + +impl XYZ { + pub fn new(x: f64, y: f64, z: f64) -> XYZ { + XYZ(x, y, z) + } + + pub fn from_coo(lon: f64, lat: f64) -> XYZ { + let (sin_lon, cos_lon) = lon.sin_cos(); + let (sin_lat, cos_lat) = lat.sin_cos(); + XYZ::new(cos_lat * cos_lon, cos_lat * sin_lon, sin_lat) + } + + pub fn time(mut self, cte: f64) -> Self { + self.0 *= cte; + self.1 *= cte; + self.2 *= cte; + self + } +} + +impl Add for XYZ { + type Output = Self; + fn add(mut self, rhs: T) -> Self::Output { + self.0 += rhs.x(); + self.1 += rhs.y(); + self.2 += rhs.z(); + self + } +} + +impl Sub for XYZ { + type Output = Self; + fn sub(mut self, rhs: T) -> Self::Output { + self.0 -= rhs.x(); + self.1 -= rhs.y(); + self.2 -= rhs.z(); + self + } +} + +impl HasXYZ for XYZ { + fn x(&self) -> f64 { + self.0 + } + fn y(&self) -> f64 { + self.1 + } + fn z(&self) -> f64 { + self.2 + } +} +impl HasXYZ for &XYZ { + fn x(&self) -> f64 { + self.0 + } + fn y(&self) -> f64 { + self.1 + } + fn z(&self) -> f64 { + self.2 + } +} + +impl HasLocalRotMatrix for XYZ { + fn local_rot_matrix(&self) -> M3x3 { + let n_xy = self.x().hypot(self.y()); + let norm = n_xy.hypot(self.z()); + let (sin_lon, cos_lon, sin_lat, cos_lat) = if n_xy.eq0() { + (0.0, 1.0, 1.0, 0.0) + } else { + ( + self.y() / n_xy, + self.x() / n_xy, + self.z() / norm, + n_xy / norm, + ) + }; + if norm.eq0() { + Default::default() + } else if (norm - 1.0).eq0() { + M3x3( + XYZt(self.x(), self.y(), self.z()), // r0 + XYZt(-sin_lon, cos_lon, 0.0), // p0 + XYZt(-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat), // q0 + ) + } else { + M3x3( + XYZt(self.x() / norm, self.y() / norm, sin_lat), + XYZt(-sin_lon, cos_lon, 0.0), + XYZt(-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat), + ) + } + } +} + +/// Transpose of a 3 dimensional vector. +#[derive(Debug)] +pub struct XYZt(pub f64, pub f64, pub f64); + +impl XYZt { + pub fn new(x: f64, y: f64, z: f64) -> XYZt { + XYZt(x, y, z) + } + + pub fn from_coo(lon: f64, lat: f64) -> XYZt { + let (sin_lon, cos_lon) = lon.sin_cos(); + let (sin_lat, cos_lat) = lat.sin_cos(); + XYZt::new(cos_lat * cos_lon, cos_lat * sin_lon, sin_lat) + } + + pub fn time(mut self, cte: f64) -> Self { + self.0 *= cte; + self.1 *= cte; + self.2 *= cte; + self + } +} + +impl HasXYZ for XYZt { + fn x(&self) -> f64 { + self.0 + } + fn y(&self) -> f64 { + self.1 + } + fn z(&self) -> f64 { + self.2 + } +} + +#[derive(Debug)] +pub struct M3x3(pub XYZt, pub XYZt, pub XYZt); + +/// The defalt returns the identity matrix +impl Default for M3x3 { + fn default() -> Self { + Self( + XYZt(1.0, 0.0, 0.0), + XYZt(0.0, 1.0, 0.0), + XYZt(0.0, 0.0, 1.0), + ) + } +} + +impl M3x3 { + pub fn new(row_x: XYZt, row_y: XYZt, row_z: XYZt) -> M3x3 { + M3x3(row_x, row_y, row_z) + } + + /// Rotation matrix around the `x-axis` so that the new frame is obtain from a rotation + /// of angle `rotation_angle_rad` around the `x-axis`. + pub fn from_rotx(rotation_angle_rad: f64) -> Self { + let (s, c) = rotation_angle_rad.sin_cos(); + Self::new( + XYZt(1_f64, 0_f64, 0_f64), + XYZt(0_f64, c, -s), + XYZt(0_f64, s, c), + ) + } + + /// Rotation matrix around the `y-axis` so that the new frame is obtain from a rotation + /// of angle `rotation_angle_rad` around the `y-axis`. + pub fn from_roty(rotation_angle_rad: f64) -> Self { + let (s, c) = rotation_angle_rad.sin_cos(); + Self::new( + XYZt(c, 0_f64, s), + XYZt(0_f64, 1_f64, 0_f64), + XYZt(-s, 0_f64, c), + ) + } + + /// Rotation matrix around the `z-axis` so that the new frame is obtain from a rotation + /// of angle `rotation_angle_rad` around the `z-axis`. + pub fn from_rotz(rotation_angle_rad: f64) -> Self { + let (s, c) = rotation_angle_rad.sin_cos(); + Self::new( + XYZt(c, -s, 0_f64), + XYZt(s, c, 0_f64), + XYZt(0_f64, 0_f64, 1_f64), + ) + } + + /// Returns the rotation matrix corresponding to transform the old frame in the new frame + /// obtained by: + /// * a first rotation of `angle_z_rad` radians around the `z-axis`, leading to `x'y'z'` + /// * a second rotation of `angle_yp_rad` radians around the `y'-axis`, leading to `x''y''z''` + /// * a thrid rotation of `angle_zpp_rad` radians around the `z''-axis`, leading to `XYZ` + fn from_zyz_euler_intrinsic(angle_z_rad: f64, angle_yp_rad: f64, angle_zpp_rad: f64) -> Self { + // z1-y′-z2″ (intrinsic rotations) or z2-y-z1 (extrinsic rotations) + let (s1, c1) = angle_z_rad.sin_cos(); + let (s2, c2) = angle_yp_rad.sin_cos(); + let (s3, c3) = angle_zpp_rad.sin_cos(); + Self::new( + XYZt(c1 * c2 * c3 - s1 * s3, -c3 * s1 - c1 * c2 * s3, c1 * s2), + XYZt(c1 * s3 + c2 * c3 * s1, c1 * c3 - c2 * s1 * s3, s1 * s2), + XYZt(-c3 * s2, s2 * s3, c2), + ) + } + + pub fn xx(&self) -> f64 { + self.row_x().x() + } + pub fn xy(&self) -> f64 { + self.row_x().y() + } + pub fn xz(&self) -> f64 { + self.row_x().z() + } + pub fn yx(&self) -> f64 { + self.row_y().x() + } + pub fn yy(&self) -> f64 { + self.row_y().y() + } + pub fn yz(&self) -> f64 { + self.row_y().z() + } + pub fn zx(&self) -> f64 { + self.row_z().x() + } + pub fn zy(&self) -> f64 { + self.row_z().y() + } + pub fn zz(&self) -> f64 { + self.row_z().z() + } + + pub fn col_x(&self) -> XYZ { + XYZ(self.row_x().x(), self.row_y().x(), self.row_z().x()) + } + + pub fn col_y(&self) -> XYZ { + XYZ(self.row_x().y(), self.row_y().y(), self.row_z().y()) + } + + pub fn col_z(&self) -> XYZ { + XYZ(self.row_x().z(), self.row_y().z(), self.row_z().z()) + } + + pub fn row_x(&self) -> &XYZt { + &self.0 + } + pub fn row_y(&self) -> &XYZt { + &self.1 + } + pub fn row_z(&self) -> &XYZt { + &self.2 + } + + pub fn scale(mut self, cte: f64) -> M3x3 { + self.0 = self.0.time(cte); + self.1 = self.1.time(cte); + self.2 = self.2.time(cte); + self + } + + pub fn time(&self, rhs: &M3x3) -> M3x3 { + let rx = self.row_x(); + let ry = self.row_y(); + let rz = self.row_z(); + let cx = rhs.col_x(); + let cy = rhs.col_y(); + let cz = rhs.col_z(); + M3x3( + XYZt(rx.scalar(&cx), rx.scalar(&cy), rx.scalar(&cz)), + XYZt(ry.scalar(&cx), ry.scalar(&cy), ry.scalar(&cz)), + XYZt(rz.scalar(&cx), rz.scalar(&cy), rz.scalar(&cz)), + ) + } + + pub fn time_transpose_of(&self, rhs: &M3x3) -> M3x3 { + let rx = self.row_x(); + let ry = self.row_y(); + let rz = self.row_z(); + let cx = rhs.row_x(); + let cy = rhs.row_y(); + let cz = rhs.row_z(); + M3x3( + XYZt(rx.scalar(cx), rx.scalar(cy), rx.scalar(cz)), + XYZt(ry.scalar(cx), ry.scalar(cy), ry.scalar(cz)), + XYZt(rz.scalar(cx), rz.scalar(cy), rz.scalar(cz)), + ) + } + + pub fn transpose_time(&self, rhs: &M3x3) -> M3x3 { + let rx = self.col_x(); + let ry = self.col_y(); + let rz = self.col_z(); + let cx = rhs.col_x(); + let cy = rhs.col_y(); + let cz = rhs.col_z(); + M3x3( + XYZt(rx.scalar(&cx), rx.scalar(&cy), rx.scalar(&cz)), + XYZt(ry.scalar(&cx), ry.scalar(&cy), ry.scalar(&cz)), + XYZt(rz.scalar(&cx), rz.scalar(&cy), rz.scalar(&cz)), + ) + } + + /// Transform the input matrix in the global frame to a matrix in the local frame. + pub fn rotate_matrix(&self, m_global: &M3x3) -> M3x3 { + self.time(&m_global.time_transpose_of(&self)) + } + + /// Transform the input matrix in the local frame to a matrix in th global frame. + pub fn unrotate_matrix(&self, m_local: &M3x3) -> M3x3 { + self.transpose_time(&m_local.time(&self)) + } + + /// From global to local + /// # Input + /// - `v`: vector in the global frame + /// # Output + /// - vector in the local frame + pub fn rotate(&self, v: &XYZ) -> XYZ { + XYZ( + self.row_x().scalar(v), + self.row_y().scalar(v), + self.row_z().scalar(v), + ) + } + + /// From global to local + /// # Input + /// - `v`: vector in the global frame + /// # Output + /// - `(y, z)` the `y` and `z` coordiantes in local frame + pub fn rotate_proj_yz(&self, v: &XYZ) -> (f64, f64) { + (self.row_y().scalar(v), self.row_z().scalar(v)) + } + + /// From local to global + /// <=> rotate with the transpose of the matrix + /// # Input + /// - `v`: vector in the local frame + /// # Output + /// - vector in the global frame + pub fn unrotate(&self, v: &XYZ) -> XYZ { + XYZ( + self.col_x().scalar(v), + self.col_y().scalar(v), + self.col_z().scalar(v), + ) + } + + pub fn to_global(&self, local_rot_matrix: &M3x3) -> M3x3 { + local_rot_matrix.unrotate_matrix(&self) + } +} + +#[cfg(test)] +mod tests { + + use super::{M3x3, XYZt}; + // use crate::qty::{coo::Coo, HasLocalRotMatrix}; + + fn relative_diff(theo: f64, meas: f64) -> f64 { + (meas - theo).abs() / theo + } + + fn assert_eq_f64(theo: f64, meas: f64) { + assert!( + relative_diff(theo, meas) < 1.0e-13, + "{:e} != {:e}: rdiff: {:e}", + theo, + meas, + relative_diff(theo, meas) + ) + } + + #[test] + fn test_from_zyz_euler_intrinsic() { + // Angular distance between (0, 90) and (15 * (12 + 49/60), 27.4) + /*let z1 = (15.0 * (12.0_f64 + 49.0_f64 / 60.0_f64)).to_radians(); + let y2 = -27.4_f64.to_radians(); + let z3 = 123.0_f64.to_radians();*/ + + // FK4 to GAL + let z1 = 192.25_f64.to_radians(); + let y2 = 62.6_f64.to_radians(); // (90.0_f64 - 27.4_f64) because we rotate the North pole, not he vernal point, else - 27.4 deg + let z3 = 57_f64.to_radians(); // (180.0_f64 - 123.0_f64).to_radians(); + + let m33 = M3x3::from_zyz_euler_intrinsic(z1, y2, z3); + println!("{:?}", m33); + println!(""); + let mz1 = M3x3::from_rotz(z1); + let my2 = M3x3::from_roty(y2); + let mz3 = M3x3::from_rotz(z3); + let m33_v2 = mz1.time(&my2.time(&mz3)); + println!("{:?}", m33_v2); + + println!("----------------------------"); + + // echo "cos(d2r(192.25)) * cos(d2r(62.6)) * cos(d2r(57)) - sin(d2r(192.25)) * sin(d2r(57))" | bc -l + // -0.06698873941515072831 + /* + -0.0669887394151507166, 0.4927284660753235682, -0.8676008111514348650, + -0.8727557658519926782, -0.4503469580199613469, -0.1883746017229204474, + -0.4835389146321842459, 0.7445846332830310861, 0.4601997847838516236 + */ + + // echo "cos(d2r(192.85948)) * cos(d2r(62.87175)) * cos(d2r(57.06808)) - sin(d2r(192.85948)) * sin(d2r(57.06808))" | bc -l + + /* + import math + math.cos(math.radians(192.85948)) * math.cos(math.radians(62.87175)) * math.cos(math.radians(57.06808)) - math.sin(math.radians(192.85948)) * math.sin(math.radians(57.06808)) + + math.sin(math.radians(192.85948)) + -0.22256070189765312 + */ + + // FK5 to GAL + let z1 = 192.85948_f64.to_radians(); + let y2 = 62.87175_f64.to_radians(); // (90_f64 - (27.0 + 7.0 / 60.0 + 41.704/3600.0)).to_radians(); = 62.87175 + let z3 = 57.06808_f64.to_radians(); // 90.0 - 32.93192 + + println!(""); + println!("{:?}", z1.sin_cos()); + println!("{:?}", y2.sin_cos()); + println!("{:?}", z3.sin_cos()); + println!(""); + println!("{}", z1.cos() * y2.cos() * z3.cos() - z1.sin() * z3.sin()); + + let m33 = M3x3::from_zyz_euler_intrinsic(z1, y2, z3); + println!("{:?}", m33); + println!(""); + let mz1 = M3x3::from_rotz(z1); + let my2 = M3x3::from_roty(y2); + let mz3 = M3x3::from_rotz(z3); + let m33_v2 = mz1.time(&my2.time(&mz3)); + println!("{:?}", m33_v2); + } +} diff --git a/src/nested/map/astrometry/mod.rs b/src/nested/map/astrometry/mod.rs new file mode 100644 index 0000000..dd1528b --- /dev/null +++ b/src/nested/map/astrometry/mod.rs @@ -0,0 +1,7 @@ +//! This module contains codes duplicated from the not-yet-publish `astrometry` crate. +//! It will be replaced by the `astrometry` crate when the later will be public. +//! +//! We use it for frame transformations (such as Galactic to Equatorial). + +pub mod gal; +pub mod math; diff --git a/src/nested/map/fits/error.rs b/src/nested/map/fits/error.rs new file mode 100644 index 0000000..d9e3181 --- /dev/null +++ b/src/nested/map/fits/error.rs @@ -0,0 +1,48 @@ +use std::{io, num::ParseIntError}; + +use thiserror::Error; + +#[derive(Error, Debug)] +pub enum FitsError { + /// IO error + #[error("I/O error.")] + Io(#[from] io::Error), + #[error("Wrong FITS keyword. Expected: {expected:?}. Actual: {actual:?}).")] + UnexpectedKeyword { expected: String, actual: String }, + #[error("Value indicator not found in keyword record '{keyword_record:?}'.")] + ValueIndicatorNotFound { keyword_record: String }, + #[error("Wrong value for keyword '{keyword:}'. Expected: '{expected:}'. Actual: '{actual:}'.")] + UnexpectedValue { + keyword: String, + expected: String, + actual: String, + }, + #[error("Unsigned int value not found in keyword record '{keyword_record:}'.")] + UintValueNotFound { keyword_record: String }, + #[error("String value no found in keyword record '{keyword_record:}'.")] + StringValueNotFound { keyword_record: String }, + #[error("Parse {context:}. Error: {err:?}")] + WrongUintValue { context: String, err: ParseIntError }, + #[error("FITS not valid. Multiple Keyword '{keyword:}'.")] + MultipleKeyword { keyword: String }, + #[error("Missing keyword '{keyword:}'.")] + MissingKeyword { keyword: String }, + #[error("Incompatible keyword values for {keyword1:} and {keyword2:}.")] + IncompatibleKeywordContent { keyword1: String, keyword2: String }, + #[error("More data than the expected!")] + RemainingData, + #[error("Less data than expected!")] + PrematureEndOfData, + #[error("Unexpected number of data written!")] + UnexpectedWrittenSize, + #[error("unexpected depth. Max expected: {depth_max:}. Actual: {depth:}")] + UnexpectedDepth { depth: u8, depth_max: u8 }, + #[error("FITS not valid: '{msg:}'.")] + Custom { msg: String }, +} + +impl FitsError { + pub fn new_custom(msg: String) -> Self { + Self::Custom { msg } + } +} diff --git a/src/nested/map/fits/gz.rs b/src/nested/map/fits/gz.rs new file mode 100644 index 0000000..19f2b14 --- /dev/null +++ b/src/nested/map/fits/gz.rs @@ -0,0 +1,18 @@ +use std::io::{BufRead, BufReader, Read, Seek}; + +use flate2::read::GzDecoder; + +const GZ_MAGIC_NUM: [u8; 2] = [0x1F, 0x8B]; +const GZ_MAGIC_NUM_LEN: usize = GZ_MAGIC_NUM.len(); + +pub fn is_gz(reader: &mut BufReader) -> Result { + let mut gz_magic_bytes = [0u8; 2]; + reader.read_exact(&mut gz_magic_bytes)?; + reader.seek_relative(-(GZ_MAGIC_NUM_LEN as i64))?; + Ok(gz_magic_bytes == GZ_MAGIC_NUM) +} + +/// Returns an object implementing `BufRead` and decompressing on-th-fly the input `BufRead`. +pub fn uncompress(reader: R) -> BufReader> { + BufReader::new(GzDecoder::new(reader)) +} diff --git a/src/nested/map/fits/keywords.rs b/src/nested/map/fits/keywords.rs new file mode 100644 index 0000000..5d9267a --- /dev/null +++ b/src/nested/map/fits/keywords.rs @@ -0,0 +1,633 @@ +//! This module contains the MOC specific FITS keywords + +use std::{fmt, slice::ChunksMut, str}; + +use super::{ + error::FitsError, + read::{get_keyword, get_str_val_no_quote, parse_uint_val}, + write::write_keyword_record, +}; + +pub trait FitsCard: Sized { + const KEYWORD: &'static [u8; 8]; + + fn keyword_str() -> &'static str { + unsafe { str::from_utf8_unchecked(Self::KEYWORD) } + } + + fn keyword_string() -> String { + unsafe { String::from_utf8_unchecked(Self::KEYWORD.to_vec()) } + } + + fn parse_value(keyword_record: &[u8]) -> Result { + Self::specific_parse_value(keyword_record) + } + + fn specific_parse_value(keyword_record: &[u8]) -> Result; + + fn write_keyword_record(&self, keyword_record: &mut [u8]) -> Result<(), FitsError> { + write_keyword_record(keyword_record, Self::KEYWORD, &self.to_fits_value()); + Ok(()) + } + + /// Must be in quotes `'val'` is value type is string + fn to_fits_value(&self) -> String; + + /// Generate an error in case the parsed value does not match a pre-define list of possible values + /// To be called in `specific_parse_value`. + /// Essentially, it converts &str in String (because once the error is raised, the str in the + /// read buffer are out-of-scope. + fn predefine_val_err(parsed_value: &[u8], expected_values: &[&[u8]]) -> FitsError { + FitsError::UnexpectedValue { + keyword: Self::keyword_string(), + expected: format!( + "{:?}", + expected_values + .iter() + .map(|v| unsafe { String::from_utf8_unchecked(v.to_vec()) }) + .collect::>() + ), + actual: String::from_utf8_lossy(parsed_value).to_string(), + } + } +} + +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub enum Ordering { + Nested, + Ring, +} +impl FitsCard for Ordering { + const KEYWORD: &'static [u8; 8] = b"ORDERING"; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + match get_str_val_no_quote(keyword_record)? { + b"NESTED" => Ok(Ordering::Nested), + b"RING" => Ok(Ordering::Ring), + parsed_val => Err(Self::predefine_val_err(parsed_val, &[b"NESTED", b"RING"])), + } + } + + fn to_fits_value(&self) -> String { + String::from(match self { + Ordering::Nested => "'NESTED'", + Ordering::Ring => "'RING'", + }) + } +} + +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub enum CoordSys { + Cel, + Gal, +} +impl FitsCard for CoordSys { + const KEYWORD: &'static [u8; 8] = b"COORDSYS"; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + match get_str_val_no_quote(keyword_record)? { + b"CEL" => Ok(CoordSys::Cel), + b"GAL" => Ok(CoordSys::Gal), + parsed_val => Err(Self::predefine_val_err(parsed_val, &[b"CEL", b"GAL"])), + } + } + + fn to_fits_value(&self) -> String { + String::from(match self { + CoordSys::Cel => "'CEL'", + CoordSys::Gal => "'GAL'", + }) + } +} + +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub enum PixType { + Healpix, +} +impl FitsCard for PixType { + const KEYWORD: &'static [u8; 8] = b"PIXTYPE "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + match get_str_val_no_quote(keyword_record)? { + b"HEALPIX" => Ok(PixType::Healpix), + parsed_val => Err(Self::predefine_val_err(parsed_val, &[b"TCB"])), + } + } + + fn to_fits_value(&self) -> String { + String::from("'HEALPIX'") + } +} + +#[derive(Debug)] +pub struct Order(u8); +impl Order { + pub fn get(&self) -> u8 { + self.0 + } +} +impl FitsCard for Order { + const KEYWORD: &'static [u8; 8] = b"ORDER "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + parse_uint_val::(keyword_record).map(Self) + } + fn to_fits_value(&self) -> String { + format!("{}", &self.0) + } +} + +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub enum TForm1 { + B(Option), // u8 + I(Option), // u16 + J(Option), // u32 + K(Option), // u64 + E(Option), // f32 + D(Option), // f64 +} +impl TForm1 { + pub fn type_n_bytes(&self) -> usize { + match self { + TForm1::B(_) => size_of::(), + TForm1::I(_) => size_of::(), + TForm1::J(_) => size_of::(), + TForm1::K(_) => size_of::(), + TForm1::E(_) => size_of::(), + TForm1::D(_) => size_of::(), + } + } + pub fn n_bytes(&self) -> usize { + match self { + TForm1::B(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + TForm1::I(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + TForm1::J(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + TForm1::K(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + TForm1::E(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + TForm1::D(opt_n) => opt_n.unwrap_or(1) as usize * size_of::(), + } + } + pub fn n_pack(&self) -> u32 { + match self { + TForm1::B(opt_n) => opt_n.unwrap_or(1), + TForm1::I(opt_n) => opt_n.unwrap_or(1), + TForm1::J(opt_n) => opt_n.unwrap_or(1), + TForm1::K(opt_n) => opt_n.unwrap_or(1), + TForm1::E(opt_n) => opt_n.unwrap_or(1), + TForm1::D(opt_n) => opt_n.unwrap_or(1), + } + } +} +impl fmt::Display for TForm1 { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{} = {}", Self::keyword_str(), self.to_fits_value()) + } +} +impl FitsCard for TForm1 { + const KEYWORD: &'static [u8; 8] = b"TFORM1 "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + let val = get_str_val_no_quote(keyword_record)?; + if val.len() == 1 { + match val { + b"B" => Ok(TForm1::B(None)), + b"I" => Ok(TForm1::I(None)), + b"J" => Ok(TForm1::J(None)), + b"K" => Ok(TForm1::K(None)), + b"E" => Ok(TForm1::E(None)), + b"D" => Ok(TForm1::D(None)), + parsed_val => Err(Self::predefine_val_err( + parsed_val, + &[b"nB", b"nI", b"nJ", b"nK", b"nE", b"nD"], + )), + } + } else { + let (n, k) = val.split_at(1); + let n_str = unsafe { str::from_utf8_unchecked(n) }; + let n = n_str + .parse::() + .map_err(|err| FitsError::WrongUintValue { + context: n_str.to_string(), + err, + })?; + match k { + b"B" => Ok(TForm1::B(Some(n))), + b"I" => Ok(TForm1::I(Some(n))), + b"J" => Ok(TForm1::J(Some(n))), + b"K" => Ok(TForm1::K(Some(n))), + b"E" => Ok(TForm1::E(Some(n))), + b"D" => Ok(TForm1::D(Some(n))), + parsed_val => Err(Self::predefine_val_err( + parsed_val, + &[b"nB", b"nI", b"nJ", b"nK", b"nE", b"nD"], + )), + } + } + } + + fn to_fits_value(&self) -> String { + match self { + TForm1::B(opt_n) => format!( + "'{}B'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + TForm1::I(opt_n) => format!( + "'{}I'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + TForm1::J(opt_n) => format!( + "'{}J'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + TForm1::K(opt_n) => format!( + "'{}K'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + TForm1::E(opt_n) => format!( + "'{}E'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + TForm1::D(opt_n) => format!( + "'{}D'", + opt_n.map(|n| n.to_string()).unwrap_or(String::from("")) + ), + } + } +} + +#[derive(Debug)] +pub struct TType1(String); +impl TType1 { + pub fn get(&self) -> &str { + self.0.as_str() + } +} +impl FitsCard for TType1 { + const KEYWORD: &'static [u8; 8] = b"TTYPE1 "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + get_str_val_no_quote(keyword_record) + .map(|s| String::from_utf8_lossy(s).to_string()) + .map(Self) + } + fn to_fits_value(&self) -> String { + format!("'{}'", &self.0) + } +} + +#[derive(Debug)] +pub struct Nside(u32); +impl Nside { + pub fn get(&self) -> u32 { + self.0 + } +} +impl FitsCard for Nside { + const KEYWORD: &'static [u8; 8] = b"NSIDE "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + parse_uint_val::(keyword_record).map(Self) + } + fn to_fits_value(&self) -> String { + format!("{}", &self.0) + } +} + +#[derive(Debug)] +pub struct FirstPix(u64); +impl FirstPix { + pub fn get(&self) -> u64 { + self.0 + } +} +impl FitsCard for FirstPix { + const KEYWORD: &'static [u8; 8] = b"FIRSTPIX"; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + parse_uint_val::(keyword_record).map(Self) + } + fn to_fits_value(&self) -> String { + format!("{}", &self.0) + } +} + +#[derive(Debug)] +pub struct LastPix(u64); +impl LastPix { + pub fn get(&self) -> u64 { + self.0 + } +} +impl FitsCard for LastPix { + const KEYWORD: &'static [u8; 8] = b"LASTPIX "; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + parse_uint_val::(keyword_record).map(Self) + } + fn to_fits_value(&self) -> String { + format!("{}", &self.0) + } +} + +// Healpix map specific +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub enum IndexSchema { + Implicit, + Explicit, + Sparse, +} +impl FitsCard for IndexSchema { + const KEYWORD: &'static [u8; 8] = b"INDXSCHM"; + + fn specific_parse_value(keyword_record: &[u8]) -> Result { + match get_str_val_no_quote(keyword_record)? { + b"IMPLICIT" => Ok(IndexSchema::Implicit), + b"EXPLICIT" => Ok(IndexSchema::Explicit), + b"SPARSE" => Ok(IndexSchema::Sparse), + parsed_val => Err(Self::predefine_val_err( + parsed_val, + &[b"IMPLICIT", b"EXPLICIT", b"SPARSE"], + )), + } + } + + fn to_fits_value(&self) -> String { + String::from(match self { + IndexSchema::Implicit => "'IMPLICIT'", + IndexSchema::Explicit => "'EXPLICIT'", + IndexSchema::Sparse => "'SPARSE'", + }) + } +} + +// Usse the index in an array of Option(SkymapKeywords) for fast retrieving of the Card :) +pub trait SkymapCard: FitsCard { + const INDEX: u8; +} +impl SkymapCard for Ordering { + const INDEX: u8 = 0; +} +impl SkymapCard for CoordSys { + const INDEX: u8 = 1; +} +impl SkymapCard for Order { + const INDEX: u8 = 2; +} +impl SkymapCard for PixType { + const INDEX: u8 = 3; +} +impl SkymapCard for TForm1 { + const INDEX: u8 = 4; +} +impl SkymapCard for TType1 { + const INDEX: u8 = 5; +} +impl SkymapCard for Nside { + const INDEX: u8 = 6; +} +impl SkymapCard for FirstPix { + const INDEX: u8 = 7; +} +impl SkymapCard for LastPix { + const INDEX: u8 = 8; +} +impl SkymapCard for IndexSchema { + const INDEX: u8 = 9; +} + +#[derive(Debug)] +pub(super) struct SkymapKeywordsMap { + entries: [Option; 10], +} +impl SkymapKeywordsMap { + pub(super) fn new() -> SkymapKeywordsMap { + Self { + entries: [None, None, None, None, None, None, None, None, None, None], + } + } + + pub(super) fn insert(&mut self, entry: SkymapKeywords) -> Option { + self.entries[entry.index()].replace(entry) + } + + pub(super) fn get(&self) -> Option<&SkymapKeywords> { + self.entries[T::INDEX as usize].as_ref() + } + + #[allow(dead_code)] + pub(super) fn write_all(&self, keyword_records: &mut ChunksMut) -> Result<(), FitsError> { + for kw in self.entries.iter().filter_map(|v| v.as_ref()) { + kw.write_keyword_record(keyword_records.next().unwrap())?; + } + Ok(()) + } + + pub(super) fn check_pixtype(&self) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::PixType(PixType::Healpix)) => Ok(()), + None => Err(FitsError::MissingKeyword { + keyword: PixType::keyword_string(), + }), + _ => unreachable!(), // since there is only one elem in PixType enum + } + } + + pub(super) fn check_coordsys(&self, expected: CoordSys) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::CoordSys(actual)) => { + if *actual == expected { + Ok(()) + } else { + Err(FitsError::UnexpectedValue { + keyword: CoordSys::keyword_string(), + expected: CoordSys::Cel.to_fits_value(), + actual: actual.to_fits_value(), + }) + } + } + None => Err(FitsError::MissingKeyword { + keyword: CoordSys::keyword_string(), + }), + _ => unreachable!(), // since there is only one elem in CoorSys enum + } + } + + pub(super) fn check_ordering(&self, expected: Ordering) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::Ordering(actual)) => { + if *actual == expected { + Ok(()) + } else { + Err(FitsError::UnexpectedValue { + keyword: Ordering::keyword_string(), + expected: expected.to_fits_value(), + actual: actual.to_fits_value(), + }) + } + } + _ => Err(FitsError::MissingKeyword { + keyword: Ordering::keyword_string(), + }), + } + } + + pub(super) fn check_firstpix(&self, expected: u64) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::FirstPix(FirstPix(actual))) => { + if *actual == expected { + Ok(()) + } else { + Err(FitsError::UnexpectedValue { + keyword: FirstPix::keyword_string(), + expected: expected.to_string(), + actual: actual.to_string(), + }) + } + } + None => Err(FitsError::MissingKeyword { + keyword: FirstPix::keyword_string(), + }), + _ => unreachable!(), + } + } + + pub(super) fn check_lastpix(&self, expected: u64) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::LastPix(LastPix(actual))) => { + if *actual == expected { + Ok(()) + } else { + Err(FitsError::UnexpectedValue { + keyword: LastPix::keyword_string(), + expected: expected.to_string(), + actual: actual.to_string(), + }) + } + } + None => Err(FitsError::MissingKeyword { + keyword: LastPix::keyword_string(), + }), + _ => unreachable!(), + } + } + + pub(super) fn check_index_schema(&self, expected: IndexSchema) -> Result<(), FitsError> { + match self.get::() { + Some(SkymapKeywords::IndexSchema(actual)) => { + if *actual == expected { + Ok(()) + } else { + Err(FitsError::UnexpectedValue { + keyword: IndexSchema::keyword_string(), + expected: expected.to_fits_value(), + actual: actual.to_fits_value(), + }) + } + } + _ => Err(FitsError::MissingKeyword { + keyword: Ordering::keyword_string(), + }), + } + } +} + +#[derive(Debug)] +pub enum SkymapKeywords { + Ordering(Ordering), + CoordSys(CoordSys), + Order(Order), + PixType(PixType), + TForm1(TForm1), + TType1(TType1), + Nside(Nside), + FirstPix(FirstPix), + LastPix(LastPix), + IndexSchema(IndexSchema), +} +impl SkymapKeywords { + pub(super) fn is_skymap_kw(keyword_record: &[u8]) -> Result, FitsError> { + // I have not yet found how to match on the FitsCard::KEYWORD associated constant :o/ + match get_keyword(keyword_record) { + b"ORDERING" => Ordering::parse_value(keyword_record) + .map(SkymapKeywords::Ordering) + .map(Some), + b"COORDSYS" => CoordSys::parse_value(keyword_record) + .map(SkymapKeywords::CoordSys) + .map(Some), + b"ORDER " => Order::parse_value(keyword_record) + .map(SkymapKeywords::Order) + .map(Some), + b"PIXTYPE " => PixType::parse_value(keyword_record) + .map(SkymapKeywords::PixType) + .map(Some), + b"TFORM1 " => TForm1::parse_value(keyword_record) + .map(SkymapKeywords::TForm1) + .map(Some), + b"TTYPE1 " => TType1::parse_value(keyword_record) + .map(SkymapKeywords::TType1) + .map(Some), + b"NSIDE " => Nside::parse_value(keyword_record) + .map(SkymapKeywords::Nside) + .map(Some), + b"FIRSTPIX" => FirstPix::parse_value(keyword_record) + .map(SkymapKeywords::FirstPix) + .map(Some), + b"LASTPIX " => LastPix::parse_value(keyword_record) + .map(SkymapKeywords::LastPix) + .map(Some), + b"INDXSCHM" => IndexSchema::parse_value(keyword_record) + .map(SkymapKeywords::IndexSchema) + .map(Some), + _ => Ok(None), + } + } + + fn index(&self) -> usize { + (match self { + SkymapKeywords::Ordering(_) => Ordering::INDEX, + SkymapKeywords::CoordSys(_) => CoordSys::INDEX, + SkymapKeywords::Order(_) => Order::INDEX, + SkymapKeywords::PixType(_) => PixType::INDEX, + SkymapKeywords::TForm1(_) => TForm1::INDEX, + SkymapKeywords::TType1(_) => TType1::INDEX, + SkymapKeywords::Nside(_) => Nside::INDEX, + SkymapKeywords::FirstPix(_) => FirstPix::INDEX, + SkymapKeywords::LastPix(_) => LastPix::INDEX, + SkymapKeywords::IndexSchema(_) => IndexSchema::INDEX, + }) as usize + } + + pub(super) fn keyword(&self) -> &'static [u8; 8] { + match self { + SkymapKeywords::Ordering(_) => Ordering::KEYWORD, + SkymapKeywords::CoordSys(_) => CoordSys::KEYWORD, + SkymapKeywords::Order(_) => Order::KEYWORD, + SkymapKeywords::PixType(_) => PixType::KEYWORD, + SkymapKeywords::TForm1(_) => TForm1::KEYWORD, + SkymapKeywords::TType1(_) => TType1::KEYWORD, + SkymapKeywords::Nside(_) => Nside::KEYWORD, + SkymapKeywords::FirstPix(_) => FirstPix::KEYWORD, + SkymapKeywords::LastPix(_) => LastPix::KEYWORD, + SkymapKeywords::IndexSchema(_) => IndexSchema::KEYWORD, + } + } + + pub(super) fn keyword_str(&self) -> &str { + unsafe { str::from_utf8_unchecked(self.keyword()) }.trim_end() + } + + fn write_keyword_record(&self, keyword_record: &mut [u8]) -> Result<(), FitsError> { + match self { + SkymapKeywords::Ordering(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::CoordSys(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::Order(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::PixType(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::TForm1(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::TType1(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::Nside(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::FirstPix(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::LastPix(kw) => kw.write_keyword_record(keyword_record), + SkymapKeywords::IndexSchema(kw) => kw.write_keyword_record(keyword_record), + } + } +} diff --git a/src/nested/map/fits/mod.rs b/src/nested/map/fits/mod.rs new file mode 100644 index 0000000..def9ff7 --- /dev/null +++ b/src/nested/map/fits/mod.rs @@ -0,0 +1,5 @@ +pub mod error; +pub mod gz; +pub mod keywords; +pub mod read; +pub mod write; diff --git a/src/nested/map/fits/read.rs b/src/nested/map/fits/read.rs new file mode 100644 index 0000000..a96ebf0 --- /dev/null +++ b/src/nested/map/fits/read.rs @@ -0,0 +1,414 @@ +use std::{ + io, + io::{BufRead, BufReader, Read, Seek}, + num::ParseIntError, + slice::ChunksExact, + str::{self, FromStr}, +}; + +use byteorder::{BigEndian, ReadBytesExt}; +use log::{debug, warn}; + +use super::{ + super::skymap::{SkyMap, SkyMapArray}, + error::FitsError, + gz::{is_gz, uncompress}, + keywords::{ + CoordSys, FitsCard, IndexSchema, Nside, Order, Ordering, SkymapKeywords, SkymapKeywordsMap, + TForm1, TType1, + }, +}; +use crate::{depth, is_nside, n_hash}; + +/// We expect the FITS file to be a BINTABLE containing a map. +/// [Here](https://gamma-astro-data-formats.readthedocs.io/en/latest/skymaps/healpix/index.html) +/// a description of the format. +/// We so far implemented a subset of the format only: +/// * `INDXSCHM= 'IMPLICIT'` +/// * `ORDERING= 'NESTED '` +/// To be fast (in execution and development), we start by a non-flexible approach in which we +/// expect the BINTABLE extension to contains: +/// ```bash +/// XTENSION= 'BINTABLE' / binary table extension +/// BITPIX = 8 / array data type +/// NAXIS = 2 / number of array dimensions +/// NAXIS1 = ?? / length of dimension 1 +/// NAXIS2 = ?? / length of dimension 2 +/// PCOUNT = 0 / number of group parameters +/// GCOUNT = 1 / number of groups +/// TFIELDS = ?? / number of table fields +/// TTYPE1 = 'XXX' // SHOULD STARS WITH 'PROB', else WARNING +/// TFORM1 = 'XXX' // MUST CONTAINS D (f64) or E (f32) +/// TUNIT1 = 'pix-1 ' +/// TTYPE2 = ??? +/// TFORM2 = ??? +/// ... +/// MOC = T +/// PIXTYPE = 'HEALPIX ' / HEALPIX pixelisation +/// ORDERING= 'NESTED ' / Pixel ordering scheme: RING, NESTED, or NUNIQ +/// COORDSYS= 'C ' // WARNING if not found +/// NSIDE = ?? / MOC resolution (best nside) +/// or +/// ORDER = ?? / MOC resolution (best order), superseded by NSIDE +/// / (because NSIDE which are not power of 2 are possible in RING) +/// INDXSCHM= 'IMPLICIT' / Indexing: IMPLICIT or EXPLICIT +/// ... +/// END +/// ``` +/// +/// # Params +/// * `reader`: the reader over the FITS content +/// +/// # Info +/// Supports gz input stream +/// +pub fn from_fits_skymap(mut reader: BufReader) -> Result { + if is_gz(&mut reader)? { + // Probably need to build an explicit map: + // highly compressed skymaps are IMPLICIT with a lot of zeros. + from_fits_skymap_internal(uncompress(reader)) + } else { + from_fits_skymap_internal(reader) + } +} + +pub fn from_fits_skymap_internal(mut reader: R) -> Result { + let mut header_block = [b' '; 2880]; + consume_primary_hdu(&mut reader, &mut header_block)?; + // Read the extension HDU + let mut it80 = next_36_chunks_of_80_bytes(&mut reader, &mut header_block)?; + // See Table 10 and 17 in https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf + check_keyword_and_val(it80.next().unwrap(), b"XTENSION", b"'BINTABLE'")?; + check_keyword_and_val(it80.next().unwrap(), b"BITPIX ", b"8")?; + check_keyword_and_val(it80.next().unwrap(), b"NAXIS ", b"2")?; + let n_bytes_per_row = check_keyword_and_parse_uint_val::(it80.next().unwrap(), b"NAXIS1 ")?; + let n_rows = check_keyword_and_parse_uint_val::(it80.next().unwrap(), b"NAXIS2 ")?; + check_keyword_and_val(it80.next().unwrap(), b"PCOUNT ", b"0")?; + check_keyword_and_val(it80.next().unwrap(), b"GCOUNT ", b"1")?; + let n_cols = check_keyword_and_parse_uint_val::(it80.next().unwrap(), b"TFIELDS ")?; + let colname_1 = TType1::parse_value(it80.next().unwrap())?; + debug!("Skymap column name: {}", colname_1.get()); + let coltype_1 = TForm1::parse_value(it80.next().unwrap())?; + debug!("Skymap column type: {}", coltype_1.to_fits_value()); + + // nbits = |BITPIX|xGCOUNTx(PCOUNT+NAXIS1xNAXIS2x...xNAXISn) + // In our case (bitpix = Depends on TForm, GCOUNT = 1, PCOUNT = 0) => nbytes = n_cells * size_of(T) + // let data_size n_bytes as usize * n_cells as usize; // N_BYTES ok since BITPIX = 8 + // Read MOC keywords + let mut skymap_kws = SkymapKeywordsMap::new(); + 'hr: loop { + for kw_record in &mut it80 { + match SkymapKeywords::is_skymap_kw(kw_record)? { + Some(kw) => { + if let Some(previous_mkw) = skymap_kws.insert(kw) { + // A FITS keyword MUST BE uniq (I may be more relax here, taking the last one and not complaining) + warn!( + "Keyword '{}' found more than once in a same HDU! We use the first occurrence.", + previous_mkw.keyword_str() + ); + skymap_kws.insert(previous_mkw); + } + } + None => { + if &kw_record[0..4] == b"END " { + break 'hr; + } else { + debug!("Ignored FITS card: {}", unsafe { + str::from_utf8_unchecked(kw_record) + }) + } + } + } + } + // Read next 2880 bytes + it80 = next_36_chunks_of_80_bytes(&mut reader, &mut header_block)?; + } + // Check constant header params + skymap_kws.check_pixtype()?; // = HEALPIX + skymap_kws.check_ordering(Ordering::Nested)?; // So far we support only 'NESTED' + skymap_kws.check_coordsys(CoordSys::Cel)?; // So far we support only Celestial coordinates (not Galactic) + skymap_kws.check_index_schema(IndexSchema::Implicit)?; // So far we support only 'IMLPLICIT' + skymap_kws.check_firstpix(0)?; + // Check number of columns + // - we so far support only map having a single column of values + if n_cols != 1 { + return Err(FitsError::UnexpectedValue { + keyword: String::from("TFIELDS"), + expected: 1.to_string(), + actual: n_cols.to_string(), + }); + } + // Check whether lastpix + let depth = match skymap_kws.get::() { + Some(SkymapKeywords::Order(order)) => Ok(order.get()), + None => match skymap_kws.get::() { + Some(SkymapKeywords::Nside(nside)) => { + let nside = nside.get(); + if is_nside(nside) { + Ok(depth(nside)) + } else { + Err(FitsError::new_custom(format!( + "Nside is not valid (to be used in nested mode at least): {}", + nside + ))) + } + } + None => Err(FitsError::new_custom(String::from( + "Both keywords 'ORDER' and 'NSIDE' are missing!", + ))), + _ => unreachable!(), + }, + _ => unreachable!(), + }?; + let n_hash = n_hash(depth); + skymap_kws.check_lastpix(n_hash)?; + // Check whether TForm compatible with N bytes per row + let n_hash_2 = coltype_1.n_pack() as u64 * n_rows; + if n_hash != n_hash_2 { + return Err(FitsError::new_custom(format!( + "Number of elements {} do not match number of HEALPix cells {}", + n_hash_2, n_hash + ))); + } + // Check n_bytes_per_row + if n_bytes_per_row != coltype_1.n_bytes() as u64 { + return Err(FitsError::new_custom(format!( + "Number of bytes per row {} do not match TFORM1 = {}", + n_bytes_per_row, + coltype_1.to_fits_value() + ))); + } + + // Read data + match coltype_1 { + TForm1::B(_) => (0..n_hash) + .map(|_| reader.read_u8()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::U8(v.into_boxed_slice()), + }), + TForm1::I(_) => (0..n_hash) + .map(|_| reader.read_i16::()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::I16(v.into_boxed_slice()), + }), + TForm1::J(_) => (0..n_hash) + .map(|_| reader.read_i32::()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::I32(v.into_boxed_slice()), + }), + TForm1::K(_) => (0..n_hash) + .map(|_| reader.read_i64::()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::I64(v.into_boxed_slice()), + }), + TForm1::E(_) => (0..n_hash) + .map(|_| reader.read_f32::()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::F32(v.into_boxed_slice()), + }), + TForm1::D(_) => (0..n_hash) + .map(|_| reader.read_f64::()) + .collect::, io::Error>>() + .map(|v| SkyMap { + depth, + values: SkyMapArray::F64(v.into_boxed_slice()), + }), + } + .map_err(FitsError::Io) +} + +const VALUE_INDICATOR: &[u8; 2] = b"= "; + +/// # Params +/// - `header_block`: re-usable header block used to avoid multiple allocations +fn consume_primary_hdu( + reader: &mut R, + header_block: &mut [u8; 2880], +) -> Result<(), FitsError> { + let mut chunks_of_80 = next_36_chunks_of_80_bytes(reader, header_block)?; + // SIMPLE = 'T' => file compliant with the FITS standard + check_keyword_and_val(chunks_of_80.next().unwrap(), b"SIMPLE ", b"T")?; + chunks_of_80.next().unwrap(); // Do not check for BITPIX (we expect an empty header) + // NAXIS = 0 => we only support FITS files with no data in the primary HDU + check_keyword_and_val(chunks_of_80.next().unwrap(), b"NAXIS ", b"0")?; + // Ignore possible additional keywords + while !contains_end(&mut chunks_of_80) { + // Few chances to enter here (except if someone had a lot of things to say in the header) + chunks_of_80 = next_36_chunks_of_80_bytes(reader, header_block)?; + } + Ok(()) +} + +fn next_36_chunks_of_80_bytes<'a, R: BufRead>( + reader: &'a mut R, + header_block: &'a mut [u8; 2880], +) -> Result, FitsError> { + reader.read_exact(header_block)?; + Ok(header_block.chunks_exact(80)) +} + +fn contains_end<'a, I: Iterator>(chunks_of_80: &'a mut I) -> bool { + for kw_rc in chunks_of_80 { + debug_assert_eq!(kw_rc.len(), 80); + if &kw_rc[0..4] == b"END " { + return true; + } + } + false +} + +pub(super) fn check_keyword_and_val( + keyword_record: &[u8], + expected_kw: &[u8], + expected_val: &[u8], +) -> Result<(), FitsError> { + check_expected_keyword(keyword_record, expected_kw)?; + check_for_value_indicator(keyword_record)?; + check_expected_value(keyword_record, expected_val) +} + +fn check_keyword_and_parse_uint_val( + keyword_record: &[u8], + expected_kw: &[u8], +) -> Result +where + T: Into + FromStr, +{ + check_expected_keyword(keyword_record, expected_kw)?; + check_for_value_indicator(keyword_record)?; + parse_uint_val::(keyword_record) +} + +#[allow(dead_code)] +pub(super) fn check_keyword_and_get_str_val<'a>( + keyword_record: &'a [u8], + expected_kw: &[u8], +) -> Result<&'a str, FitsError> { + check_expected_keyword(keyword_record, expected_kw)?; + check_for_value_indicator(keyword_record)?; + // We go unsafe because FITS headers are not supposed to contain non-ASCII chars + get_str_val_no_quote(keyword_record).map(|bytes| unsafe { str::from_utf8_unchecked(bytes) }) +} + +pub(super) fn check_expected_keyword( + keyword_record: &[u8], + expected: &[u8], +) -> Result<(), FitsError> { + debug_assert!(keyword_record.len() == 80); // length of a FITS keyword-record + debug_assert!(expected.len() <= 8); // length of a FITS keyword + if &keyword_record[..expected.len()] == expected { + Ok(()) + } else { + // We know what we put in it, so unsafe is ok here + let expected = String::from(unsafe { str::from_utf8_unchecked(expected) }.trim_end()); + // Here, may contains binary data + let actual = String::from_utf8_lossy(&keyword_record[..expected.len()]) + .trim_end() + .to_string(); + // panic!("Ecpected: {}, Actual: {}", expected, String::from_utf8_lossy(&src[..])); + Err(FitsError::UnexpectedKeyword { expected, actual }) + } +} + +pub(super) fn check_for_value_indicator(keyword_record: &[u8]) -> Result<(), FitsError> { + debug_assert!(keyword_record.len() == 80); // length of a FITS keyword-record + if get_value_indicator(keyword_record) == VALUE_INDICATOR { + Ok(()) + } else { + let keyword_record = String::from_utf8_lossy(keyword_record) + .trim_end() + .to_string(); + Err(FitsError::ValueIndicatorNotFound { keyword_record }) + } +} + +pub(super) fn get_keyword(keyword_record: &[u8]) -> &[u8] { + &keyword_record[..8] +} +pub(super) fn get_value_indicator(keyword_record: &[u8]) -> &[u8] { + &keyword_record[8..10] +} +pub(super) fn get_value(keyword_record: &[u8]) -> &[u8] { + &keyword_record[10..] +} +pub(super) fn get_left_trimmed_value(keyword_record: &[u8]) -> &[u8] { + get_value(keyword_record).trim_ascii_start() +} + +pub(super) fn check_expected_value( + keyword_record: &[u8], + expected: &[u8], +) -> Result<(), FitsError> { + debug_assert!(keyword_record.len() == 80); // length of a FITS keyword-record + let src = get_value(keyword_record); + let lt_src = src.trim_ascii_start(); + if lt_src.len() >= expected.len() && <_src[..expected.len()] == expected { + Ok(()) + } else { + let keyword = String::from_utf8_lossy(&src[0..8]).trim_end().to_string(); + // We know what we put in it, so unsae is ok here + let expected = String::from(unsafe { str::from_utf8_unchecked(expected) }); + // Here, may contains binary data + let actual = String::from_utf8_lossy(<_src[..expected.len()]).to_string(); + Err(FitsError::UnexpectedValue { + keyword, + expected, + actual, + }) + } +} + +/// We know that the expected value does not contains a simple quote. +/// A trim_end is applied so that the result does not contain leading or trailing spaces. +pub(super) fn get_str_val_no_quote(keyword_record: &[u8]) -> Result<&[u8], FitsError> { + let mut it = get_left_trimmed_value(keyword_record).split_inclusive(|c| *c == b'\''); + if let Some([b'\'']) = it.next() { + if let Some([subslice @ .., b'\'']) = it.next() { + return Ok(subslice.trim_ascii_end()); + } + } + let keyword_record = String::from_utf8_lossy(keyword_record) + .trim_end() + .to_string(); + Err(FitsError::StringValueNotFound { keyword_record }) +} + +pub(super) fn parse_uint_val(keyword_record: &[u8]) -> Result +where + T: Into + FromStr, +{ + let src = get_left_trimmed_value(keyword_record); + let to = index_of_last_digit(src); + if to == 0 { + let keyword_record = String::from_utf8_lossy(keyword_record) + .trim_end() + .to_string(); + Err(FitsError::UintValueNotFound { keyword_record }) + } else { + // we go unsafe and unwrap because we already tested for regular digits + let str_val = unsafe { str::from_utf8_unchecked(&src[..to]) }; + str_val.parse::().map_err(|e| FitsError::WrongUintValue { + context: str_val.to_string(), + err: e, + }) + } +} + +pub(super) fn index_of_last_digit(src: &[u8]) -> usize { + for (i, c) in src.iter().enumerate() { + if !c.is_ascii_digit() { + return i; + } + } + src.len() +} diff --git a/src/nested/map/fits/write.rs b/src/nested/map/fits/write.rs new file mode 100644 index 0000000..49c4417 --- /dev/null +++ b/src/nested/map/fits/write.rs @@ -0,0 +1,214 @@ +use std::io::Write; + +use num_traits::ToBytes; + +use crate::{n_hash, nside}; + +use super::error::FitsError; + +/// Trait marking the type of the values written in the FITS map. +pub trait FitsCompatibleValue: ToBytes { + /// Size, in bytes, of a value. + fn fits_naxis1() -> u8; + fn fits_tform() -> &'static str; +} + +impl FitsCompatibleValue for u8 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "B" + } +} + +impl FitsCompatibleValue for i16 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "I" + } +} + +impl FitsCompatibleValue for i32 { + fn fits_naxis1() -> u8 { + 4 + } + + fn fits_tform() -> &'static str { + "J" + } +} + +impl FitsCompatibleValue for i64 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "K" + } +} + +impl FitsCompatibleValue for u32 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "J" + } +} + +impl FitsCompatibleValue for u64 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "K" + } +} + +impl FitsCompatibleValue for f32 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "E" + } +} + +impl FitsCompatibleValue for f64 { + fn fits_naxis1() -> u8 { + size_of::() as u8 + } + + fn fits_tform() -> &'static str { + "D" + } +} + +// Add read map from fits!! + +/// Write the values in a FITS HEALPix map. +/// +/// # Params +/// * `values`: array of values, the index of the value correspond to the HEALPix cell the value is +/// associated with. +pub fn write_implicit_skymap_fits( + mut writer: R, + values: &[T], +) -> Result<(), FitsError> { + let n_cells = values.len() as u64; + let depth = (values.len() / 12).trailing_zeros() as u8 >> 1; + if n_cells != n_hash(depth) { + return Err(FitsError::new_custom(format!( + "Number of cell {} not compatible with an HEALPix depth of {}. Expected: {}.", + n_cells, + depth, + n_hash(depth) + ))); + } + write_skymap_fits_header::<_, T>(&mut writer, depth)?; + for value in values { + writer.write_all(value.to_be_bytes().as_ref())?; + } + // Complete FITS block of 2880 bytes + let mod2880 = (n_cells as usize * T::fits_naxis1() as usize) % 2880; + if mod2880 != 0 { + writer.write_all(&vec![0_u8; 2880 - mod2880])?; + } + Ok(()) +} + +#[allow(dead_code)] +pub(crate) fn write_uint_keyword_record(dest: &mut [u8], keyword: &[u8; 8], val: u64) { + write_keyword_record(dest, keyword, &val.to_string()) +} + +#[allow(dead_code)] +pub(crate) fn write_str_keyword_record(dest: &mut [u8], keyword: &[u8; 8], val: &str) { + write_keyword_record(dest, keyword, &format!("'{}'", val)) +} + +/// # Params +/// - `dest` destination, must contains 80 bytes +/// - `value_part` is string, must be already quote: `'str_value'` +pub(crate) fn write_keyword_record(dest: &mut [u8], keyword: &[u8; 8], value_part: &str) { + const VALUE_INDICATOR: &[u8; 2] = b"= "; + debug_assert_eq!(dest.len(), 80); + dest[0..8].copy_from_slice(&keyword[..]); + dest[8..10].copy_from_slice(VALUE_INDICATOR); + let val_bytes = value_part.as_bytes(); + dest[10..10 + val_bytes.len()].copy_from_slice(val_bytes); +} + +fn write_uint_mandatory_keyword_record(dest: &mut [u8], keyword: &[u8; 8], val: u64) { + const VALUE_INDICATOR: &[u8; 2] = b"= "; + debug_assert_eq!(dest.len(), 80); + dest[0..8].copy_from_slice(&keyword[..]); + dest[8..10].copy_from_slice(VALUE_INDICATOR); + let val = val.to_string(); + let val_bytes = val.as_bytes(); + dest[30 - val_bytes.len()..30].copy_from_slice(val_bytes); +} + +fn write_str_mandatory_keyword_record(dest: &mut [u8], keyword: &[u8; 8], s: &str) { + const VALUE_INDICATOR: &[u8; 2] = b"= "; + debug_assert_eq!(dest.len(), 80); + dest[0..8].copy_from_slice(&keyword[..]); + dest[8..10].copy_from_slice(VALUE_INDICATOR); + let val = format!("'{:>8}'", s); + let val_bytes = val.as_bytes(); + dest[10..10 + val_bytes.len()].copy_from_slice(val_bytes); +} + +fn write_primary_hdu(writer: &mut R) -> Result<(), FitsError> { + let mut header_block = [b' '; 2880]; + header_block[0..30].copy_from_slice(b"SIMPLE = T"); + header_block[80..110].copy_from_slice(b"BITPIX = 8"); + header_block[160..190].copy_from_slice(b"NAXIS = 0"); + header_block[240..270].copy_from_slice(b"EXTEND = T"); + header_block[320..323].copy_from_slice(b"END"); + writer.write_all(&header_block[..]).map_err(FitsError::Io) +} + +fn write_skymap_fits_header( + mut writer: R, + depth: u8, +) -> Result<(), FitsError> { + let nside = nside(depth); + let n_cells = n_hash(depth); + + write_primary_hdu(&mut writer)?; + let mut header_block = [b' '; 2880]; + let mut it = header_block.chunks_mut(80); + // Write BINTABLE specific keywords in the buffer + it.next().unwrap()[0..20].copy_from_slice(b"XTENSION= 'BINTABLE'"); + it.next().unwrap()[0..30].copy_from_slice(b"BITPIX = 8"); + it.next().unwrap()[0..30].copy_from_slice(b"NAXIS = 2"); + write_uint_mandatory_keyword_record(it.next().unwrap(), b"NAXIS1 ", T::fits_naxis1() as u64); + write_uint_mandatory_keyword_record(it.next().unwrap(), b"NAXIS2 ", n_cells); + it.next().unwrap()[0..30].copy_from_slice(b"PCOUNT = 0"); + it.next().unwrap()[0..30].copy_from_slice(b"GCOUNT = 1"); + it.next().unwrap()[0..30].copy_from_slice(b"TFIELDS = 1"); + it.next().unwrap()[0..20].copy_from_slice(b"TTYPE1 = 'T '"); + write_str_mandatory_keyword_record(it.next().unwrap(), b"TFORM1 ", T::fits_tform()); + it.next().unwrap()[0..20].copy_from_slice(b"PIXTYPE = 'HEALPIX '"); + it.next().unwrap()[0..20].copy_from_slice(b"ORDERING= 'NESTED '"); + it.next().unwrap()[0..20].copy_from_slice(b"EXTNAME = 'xtension'"); + write_uint_mandatory_keyword_record(it.next().unwrap(), b"NSIDE ", nside as u64); + it.next().unwrap()[0..30].copy_from_slice(b"FIRSTPIX= 0"); + write_uint_mandatory_keyword_record(it.next().unwrap(), b"LASTPIX ", n_cells); + it.next().unwrap()[0..20].copy_from_slice(b"INDXSCHM= 'IMPLICIT'"); + it.next().unwrap()[0..20].copy_from_slice(b"OBJECT = 'FULLSKY '"); + it.next().unwrap()[0..26].copy_from_slice(b"CREATOR = 'CDS HEALPix Rust'"); + it.next().unwrap()[0..3].copy_from_slice(b"END"); + // Do write the header + writer.write_all(&header_block[..]).map_err(FitsError::Io) +} diff --git a/src/nested/map/img.rs b/src/nested/map/img.rs new file mode 100644 index 0000000..d109b20 --- /dev/null +++ b/src/nested/map/img.rs @@ -0,0 +1,624 @@ +use std::{ + error::Error, + f64::consts::PI, + fs::File, + io::{self, BufWriter, Write}, + ops::RangeInclusive, + path::Path, +}; + +use colorous::Gradient; +use mapproj::{ + img2celestial::Img2Celestial, img2proj::ReversedEastPngImgXY2ProjXY, pseudocyl::mol::Mol, + CanonicalProjection, CenteredProjection, ImgXY, LonLat, +}; + +use crate::nested::map::astrometry::math::Coo; +use crate::nested::{self, map::astrometry::gal::Galactic, n_hash}; + +/// Wanted image coordinate frame. +pub enum ImgCooFrame { + Equatorial, + Galactic, +} + +pub enum PosConversion { + SameMapAndImg, + EqMap2GalImg, + GalMap2EqImg, +} +impl PosConversion { + pub fn convert_img_pos_to_map_pos(&self) -> fn(f64, f64) -> (f64, f64) { + match self { + PosConversion::SameMapAndImg => |l, b| (l, b), + PosConversion::EqMap2GalImg => Self::gal2eq(), + PosConversion::GalMap2EqImg => Self::eq2gal(), + } + } + pub fn convert_map_pos_to_img_pos(&self) -> fn(f64, f64) -> (f64, f64) { + match self { + PosConversion::SameMapAndImg => |l, b| (l, b), + PosConversion::EqMap2GalImg => Self::eq2gal(), + PosConversion::GalMap2EqImg => Self::gal2eq(), + } + } + pub fn gal2eq() -> fn(f64, f64) -> (f64, f64) { + |lon, lat| { + const GAL: Galactic = Galactic::new_for_icrs_including_fk5_icrs_offsets_from_mignard2002(); + let coo_eq = Coo::new(lon, lat); + let coo_gal = GAL.coo_gal2eq(&coo_eq); + (coo_gal.lon, coo_gal.lat) + } + } + pub fn eq2gal() -> fn(f64, f64) -> (f64, f64) { + |lon, lat| { + const GAL: Galactic = Galactic::new_for_icrs_including_fk5_icrs_offsets_from_mignard2002(); + let coo_gal = Coo::new(lon, lat); + let coo_eq = GAL.coo_eq2gal(&coo_gal); + (coo_eq.lon, coo_eq.lat) + } + } +} + +pub trait Val: PartialOrd { + fn to_f64(&self) -> f64; +} + +impl Val for u8 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for u16 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for u32 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for u64 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for i8 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for i16 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for i32 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for i64 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for f32 { + fn to_f64(&self) -> f64 { + *self as f64 + } +} + +impl Val for f64 { + fn to_f64(&self) -> f64 { + *self + } +} + +pub enum ColorMapFunctionType { + Linear, + LinearLog, + LinearPow2, + LinearSqrt, + LinearAsinh, + Log, + Pow2, + Sqrt, + Asinh, +} + +/// +/// * `min`: minimum value corresponding to 0 in the ColorMap +/// * `max`: maximum value corresponding to 1 in the ColorMap +#[derive(Debug, Copy, Clone)] +pub enum ColorMapFunction { + Linear { + min: f64, + max: f64, + }, + // Use formula: 0.5*log(x+0.01)+1 returning value in [0, 1], x max must be 0.99 + LinearLog { + min: f64, + max: f64, + }, + LinearPow2 { + min: f64, + max: f64, + }, + LinearSqrt { + min: f64, + max: f64, + }, + LinearAsinh { + min: f64, + max: f64, + }, + Log { + min: f64, + max: f64, + fmin: f64, + fmax: f64, + }, + Pow2 { + min: f64, + max: f64, + fmin: f64, + fmax: f64, + }, + Sqrt { + min: f64, + max: f64, + fmin: f64, + fmax: f64, + }, + Asinh { + min: f64, + max: f64, + fmin: f64, + fmax: f64, + }, +} +impl ColorMapFunction { + pub fn new_from(cm_type: ColorMapFunctionType, min: f64, max: f64) -> Self { + match cm_type { + ColorMapFunctionType::Linear => Self::new_linear(min, max), + ColorMapFunctionType::LinearLog => Self::new_linear_log(min, max), + ColorMapFunctionType::LinearPow2 => Self::new_linear_pow2(min, max), + ColorMapFunctionType::LinearSqrt => Self::new_linear_sqrt(min, max), + ColorMapFunctionType::LinearAsinh => Self::new_linear_asinh(min, max), + ColorMapFunctionType::Log => Self::new_log(min, max), + ColorMapFunctionType::Pow2 => Self::new_pow2(min, max), + ColorMapFunctionType::Sqrt => Self::new_sqrt(min, max), + ColorMapFunctionType::Asinh => Self::new_asinh(min, max), + } + } + pub fn new_linear(min: f64, max: f64) -> Self { + ColorMapFunction::Linear { min, max } + } + pub fn new_linear_log(min: f64, max: f64) -> Self { + ColorMapFunction::LinearLog { min, max } + } + pub fn new_linear_pow2(min: f64, max: f64) -> Self { + ColorMapFunction::LinearPow2 { min, max } + } + pub fn new_linear_sqrt(min: f64, max: f64) -> Self { + ColorMapFunction::LinearSqrt { min, max } + } + pub fn new_linear_asinh(min: f64, max: f64) -> Self { + ColorMapFunction::LinearAsinh { min, max } + } + + pub fn new_log(min: f64, max: f64) -> Self { + ColorMapFunction::Log { + min, + max, + fmin: (min + 0.01).ln(), + fmax: (max + 0.01).ln(), + } + } + pub fn new_pow2(min: f64, max: f64) -> Self { + ColorMapFunction::Pow2 { + min, + max, + fmin: min * min, + fmax: max * max, + } + } + pub fn new_sqrt(min: f64, max: f64) -> Self { + ColorMapFunction::Sqrt { + min, + max, + fmin: min.sqrt(), + fmax: max.sqrt(), + } + } + pub fn new_asinh(min: f64, max: f64) -> Self { + ColorMapFunction::Asinh { + min, + max, + fmin: min.asinh(), + fmax: max.asinh(), + } + } + + /// Returns a value in `[0, 1]`. + pub fn value(&self, value: f64) -> f64 { + match self { + ColorMapFunction::Linear { min, max } => Self::linear_val(*min, *max, value, |v| v), + ColorMapFunction::LinearLog { min, max } => + // 1 + log10(0.01) / 2 = 0 since log10(10^-2) = -2 * log10(10) = -2 * ln(10)/ln(10) + // 1 + log10(1.00) / 2 = 1 since log10(1) = ln(1)/ln(10) and ln(1) = 0 + { + Self::linear_val(*min, *max, value, |v| 1.0 + 0.5 * (0.99 * v + 0.01).log10()) + } + ColorMapFunction::LinearPow2 { min, max } => Self::linear_val(*min, *max, value, |v| v * v), + ColorMapFunction::LinearSqrt { min, max } => { + Self::linear_val(*min, *max, value, |v| v.sqrt()) + } + ColorMapFunction::LinearAsinh { min, max } => { + Self::linear_val(*min, *max, value, |v| v.asinh() / 1.0_f64.asinh()) + } + ColorMapFunction::Log { + min, + max, + fmin, + fmax, + } => Self::val(*min, *max, *fmin, *fmax, value, |v| (v + 0.01).ln()), + ColorMapFunction::Pow2 { + min, + max, + fmin, + fmax, + } => Self::val(*min, *max, *fmin, *fmax, value, |v| v * v), + ColorMapFunction::Sqrt { + min, + max, + fmin, + fmax, + } => Self::val(*min, *max, *fmin, *fmax, value, |v| v.sqrt()), + ColorMapFunction::Asinh { + min, + max, + fmin, + fmax, + } => Self::val(*min, *max, *fmin, *fmax, value, |v| v.asinh()), + } + } + /// Reverse of the value, i.e. returns a value in `[1, 0]`. + pub fn value_rev(&self, value: f64) -> f64 { + 1.0 - self.value(value) + } + + // Function F takes a values in [0, 1] and must return a value in [0, 1]. + fn linear_val(min: f64, max: f64, value: f64, f: F) -> f64 + where + F: Fn(f64) -> f64, + { + if value > max { + 1.0 + } else if value < min { + 0.0 + } else { + let a = (value - min) / (max - min); + f(a) + } + } + fn val(min: f64, max: f64, fmin: f64, fmax: f64, value: f64, f: F) -> f64 + where + F: Fn(f64) -> f64, + { + if value > max { + 1.0 + } else if value < min { + 0.0 + } else { + (f(value) - fmin) / (fmax - fmin) + } + } +} + +// MOM: passer une liste de (NUNIQ, Val) ? + +/// Returns an RGBA array (each pixel is made of 4 successive u8: RGBA) using the Mollweide projection. +/// +/// # Params +/// * `smoc`: the Spatial MOC to be print; +/// * `size`: the `(X, Y)` number of pixels in the image; +/// * `proj_center`: the `(lon, lat)` coordinates of the center of the projection, in radians, +/// if different from `(0, 0)`; +/// * `proj_bounds`: the `(X, Y)` bounds of the projection, if different from the default values +/// which depends on the projection. For unbounded projections, de default value +/// is `(-PI..PI, -PI..PI)`. +pub fn to_img_default( + skymap_implicit: &[V], + img_size: (u16, u16), + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Option, + color_map_func_type: Option, +) -> Result, Box> { + let proj = Mol::new(); + let color_map = color_map.unwrap_or(colorous::TURBO); + let color_map_func_type = color_map_func_type.unwrap_or(ColorMapFunctionType::Linear); + to_img( + skymap_implicit, + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + ) +} + +/// Returns an RGBA array (each pixel is made of 4 successive u8: RGBA). +/// +/// # Params +/// * `smoc`: the Spatial MOC to be print; +/// * `size`: the `(X, Y)` number of pixels in the image; +/// * `proj`: a projection, if different from Mollweide; +/// * `proj_center`: the `(lon, lat)` coordinates of the center of the projection, in radians, +/// if different from `(0, 0)`; +/// * `proj_bounds`: the `(X, Y)` bounds of the projection, if different from the default values +/// which depends on the projection. For unbounded projections, de default value +/// is `(-PI..PI, -PI..PI)`. +/// * `pos_convert`: to handle a different coordinate system between the skymap and the image. +/// * `color_map`: +/// * `color_map_func`: the color map fonction, build it using: +pub fn to_img( + skymap_implicit: &[V], + img_size: (u16, u16), + proj: P, + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Gradient, + color_map_func_type: ColorMapFunctionType, +) -> Result, Box> { + // Make the map a trait SkyMap: Iterable { depth(), get(hash) } + + let n_cells = skymap_implicit.len() as u64; + let depth = (skymap_implicit.len() / 12).trailing_zeros() as u8 >> 1; + if n_cells != n_hash(depth) { + return Err( + format!( + "Number of cell {} not compatible with an HEALPix depth of {}. Expected: {}.", + n_cells, + depth, + n_hash(depth) + ) + .into(), + ); + } + + // Unwrap is ok here since we already tested 'n_cell' which must be > 0. + let mut iter = skymap_implicit.iter(); + let first_value = iter.next().unwrap(); + let (min, max) = iter.fold((first_value, first_value), |(min, max), val| { + ( + std::cmp::min_by(val, min, |a, b| a.partial_cmp(b).unwrap()), + std::cmp::max_by(val, max, |a, b| a.partial_cmp(b).unwrap()), + ) + }); + let (min, max) = (min.to_f64(), max.to_f64()); + /*let max = skymap_implicit + .iter() + .max_by(|a, b| ) + .map(|v| v.to_f64()) + .unwrap_or(0_f64);*/ + let color_map_func = ColorMapFunction::new_from(color_map_func_type, min, max); + + let (size_x, size_y) = img_size; + let mut v: Vec = Vec::with_capacity((size_x as usize * size_y as usize) << 2); + + let (proj_range_x, proj_range_y) = proj_bounds.unwrap_or(( + proj + .bounds() + .x_bounds() + .as_ref() + .cloned() + .unwrap_or_else(|| -PI..=PI), + proj + .bounds() + .y_bounds() + .as_ref() + .cloned() + .unwrap_or_else(|| -PI..=PI), + )); + + let img2proj = + ReversedEastPngImgXY2ProjXY::from((size_x, size_y), (&proj_range_x, &proj_range_y)); + let mut img2cel = Img2Celestial::new(img2proj, CenteredProjection::new(proj)); + if let Some((lon, lat)) = proj_center { + img2cel.set_proj_center_from_lonlat(&LonLat::new(lon, lat)); + } + + let hpx = nested::get(depth); + + let pos_convert = pos_convert.unwrap_or(PosConversion::SameMapAndImg); + let mappos2imgpos = pos_convert.convert_map_pos_to_img_pos(); + let imgpos2mappos = pos_convert.convert_img_pos_to_map_pos(); + + for y in 0..size_y { + for x in 0..size_x { + if let Some(lonlat) = img2cel.img2lonlat(&ImgXY::new(x as f64, y as f64)) { + let (lon, lat) = imgpos2mappos(lonlat.lon(), lonlat.lat()); + let idx = hpx.hash(lon, lat) as usize; + let val = &skymap_implicit[idx]; + // num_traits::cast::ToPrimitive contains to_f64 + let color = color_map.eval_continuous(color_map_func.value(val.to_f64())); + v.push(color.r); + v.push(color.g); + v.push(color.b); + v.push(255); + } else { + // Not in the proj area + v.push(255); + v.push(255); + v.push(255); + v.push(0); + } + } + } + // But, in case of sparse map with cells smaller than img pixels + for (idx, val) in skymap_implicit + .iter() + .enumerate() + .filter(|(_, val)| val.to_f64() > 0.0) + { + let (lon_rad, lat_rad) = nested::center(depth, idx as u64); + let (lon_rad, lat_rad) = mappos2imgpos(lon_rad, lat_rad); + if let Some(xy) = img2cel.lonlat2img(&LonLat::new(lon_rad, lat_rad)) { + let ix = xy.x() as u16; + let iy = xy.y() as u16; + if ix < img_size.0 && iy < img_size.1 { + let from = (xy.y() as usize * size_x as usize + ix as usize) << 2; // <<2 <=> *4 + if v[from] == 0 { + let color = color_map.eval_continuous(color_map_func.value(val.to_f64())); + v[from] = color.r; + v[from + 1] = color.g; + v[from + 2] = color.b; + v[from + 3] = 128; + } + } + } + } + Ok(v) +} + +/// # Params +/// * `smoc`: the Spatial MOC to be print; +/// * `size`: the `(X, Y)` number of pixels in the image; +/// * `proj`: a projection, if different from Mollweide; +/// * `proj_center`: the `(lon, lat)` coordinates of the center of the projection, in radians, +/// if different from `(0, 0)`; +/// * `proj_bounds`: the `(X, Y)` bounds of the projection, if different from the default values +/// which depends on the projection. For unbounded projections, de default value +/// is `(-PI..PI, -PI..PI)`. +/// * `writer`: the writer in which the image is going to be written +pub fn to_png( + skymap_implicit: &[V], + img_size: (u16, u16), + proj: Option

, + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Option, + color_map_func_type: Option, + writer: W, +) -> Result<(), Box> { + let (xsize, ysize) = img_size; + let data = if let Some(proj) = proj { + let color_map = color_map.unwrap_or(colorous::TURBO); + let color_map_func_type = color_map_func_type.unwrap_or(ColorMapFunctionType::Linear); + to_img( + skymap_implicit, + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + ) + } else { + to_img_default( + skymap_implicit, + img_size, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + ) + }?; + let mut encoder = png::Encoder::new(writer, xsize as u32, ysize as u32); // Width is 2 pixels and height is 1. + encoder.set_color(png::ColorType::Rgba); + encoder.set_depth(png::BitDepth::Eight); + let mut writer = encoder.write_header()?; + writer.write_image_data(&data).expect("Wrong encoding"); + Ok(()) +} + +/// # Params +/// * `smoc`: the Spatial MOC to be print; +/// * `size`: the `(X, Y)` number of pixels in the image; +/// * `proj`: a projection, if different from Mollweide; +/// * `proj_center`: the `(lon, lat)` coordinates of the center of the projection, in radians, +/// if different from `(0, 0)`; +/// * `proj_bounds`: the `(X, Y)` bounds of the projection, if different from the default values +/// which depends on the projection. For unbounded projections, de default value +/// is `(-PI..PI, -PI..PI)`. +/// * `path`: the path of th PNG file to be written. +/// * `view`: set to true to visualize the saved image. +#[cfg(not(target_arch = "wasm32"))] +pub fn to_png_file( + skymap_implicit: &[V], + img_size: (u16, u16), + proj: Option

, + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Option, + color_map_func_type: Option, + path: &Path, + view: bool, +) -> Result<(), Box> { + // Brackets are important to be sure the file is closed before trying to open it. + { + let file = File::create(path)?; + let mut writer = BufWriter::new(file); + to_png( + skymap_implicit, + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + &mut writer, + )?; + } + if view { + show_with_default_app(path.to_string_lossy().as_ref())?; + } + Ok(()) +} + +// Adapted from https://github.com/igiagkiozis/plotly/blob/master/plotly/src/plot.rs +#[cfg(target_os = "linux")] +pub fn show_with_default_app(path: &str) -> Result<(), io::Error> { + use std::process::Command; + Command::new("xdg-open").args([path]).output()?; + // .map_err(|e| e.into())?; + Ok(()) +} + +#[cfg(target_os = "macos")] +pub fn show_with_default_app(path: &str) -> Result<(), io::Error> { + use std::process::Command; + Command::new("open").args(&[path]).output()?; + Ok(()) +} + +#[cfg(target_os = "windows")] +pub fn show_with_default_app(path: &str) -> Result<(), io::Error> { + use std::process::Command; + Command::new("cmd") + .arg("/C") + .arg(format!(r#"start {}"#, path)) + .output()?; + Ok(()) +} diff --git a/src/nested/map/mod.rs b/src/nested/map/mod.rs new file mode 100644 index 0000000..700ff5f --- /dev/null +++ b/src/nested/map/mod.rs @@ -0,0 +1,32 @@ +pub mod astrometry; +pub mod fits; +pub mod img; +pub mod skymap; + +#[cfg(test)] +mod tests { + use super::{ + img::{to_png_file, ColorMapFunctionType, PosConversion}, + skymap::*, + }; + use mapproj::pseudocyl::mol::Mol; + + #[test] + fn test_skymap() { + let mut path = "test/resources/skymap/skymap.fits"; + let skymap = SkyMap::from_fits_file(path).unwrap(); + skymap + .to_png_file::( + (1600, 800), + None, + None, + None, + Some(PosConversion::EqMap2GalImg), + None, + Some(ColorMapFunctionType::LinearLog), //Some(ColorMapFunctionType::LinearSqrt) + "test/resources/skymap/skymap.png", + false, + ) + .unwrap(); + } +} diff --git a/src/nested/map/skymap.rs b/src/nested/map/skymap.rs new file mode 100644 index 0000000..00bb972 --- /dev/null +++ b/src/nested/map/skymap.rs @@ -0,0 +1,173 @@ +use std::{ + error::Error, + fs::File, + io::{BufReader, BufWriter, Read, Seek, Write}, + ops::{Deref, RangeInclusive}, + path::Path, +}; + +use colorous::Gradient; +use mapproj::CanonicalProjection; + +// Make a trait SkyMap: Iterable { depth(), get(hash) } +// => Make a map from a FITS file without loading data in memory (bot implicit and explicit) + +use super::{ + fits::{error::FitsError, read::from_fits_skymap, write::write_implicit_skymap_fits}, + img::{show_with_default_app, to_png, ColorMapFunctionType, PosConversion}, +}; + +// Implicit map +pub struct SkyMap { + pub depth: u8, + pub values: SkyMapArray, +} + +pub enum SkyMapArray { + U8(Box<[u8]>), + I16(Box<[i16]>), + I32(Box<[i32]>), + I64(Box<[i64]>), + F32(Box<[f32]>), + F64(Box<[f64]>), +} + +impl SkyMap { + pub fn from_fits_file>(path: P) -> Result { + File::open(path) + .map_err(FitsError::Io) + .map(BufReader::new) + .and_then(SkyMap::from_fits) + } + pub fn from_fits(reader: BufReader) -> Result { + from_fits_skymap(reader) + } + + pub fn to_fits(&self, writer: W) -> Result<(), FitsError> { + match &self.values { + SkyMapArray::U8(b) => write_implicit_skymap_fits(writer, b.deref()), + SkyMapArray::I16(b) => write_implicit_skymap_fits(writer, b.deref()), + SkyMapArray::I32(b) => write_implicit_skymap_fits(writer, b.deref()), + SkyMapArray::I64(b) => write_implicit_skymap_fits(writer, b.deref()), + SkyMapArray::F32(b) => write_implicit_skymap_fits(writer, b.deref()), + SkyMapArray::F64(b) => write_implicit_skymap_fits(writer, b.deref()), + } + } + + pub fn to_png_file>( + &self, + img_size: (u16, u16), + proj: Option

, + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Option, + color_map_func_type: Option, + path: W, + view: bool, + ) -> Result<(), Box> { + File::create(path.as_ref()) + .map_err(|e| e.into()) + .map(BufWriter::new) + .and_then(|mut writer| { + self.to_png( + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + &mut writer, + ) + }) + .and_then(|()| { + if view { + show_with_default_app(path.as_ref().to_string_lossy().as_ref()).map_err(|e| e.into()) + } else { + Ok(()) + } + }) + } + + pub fn to_png( + &self, + img_size: (u16, u16), + proj: Option

, + proj_center: Option<(f64, f64)>, + proj_bounds: Option<(RangeInclusive, RangeInclusive)>, + pos_convert: Option, + color_map: Option, + color_map_func_type: Option, + writer: W, + ) -> Result<(), Box> { + match &self.values { + SkyMapArray::U8(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + SkyMapArray::I16(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + SkyMapArray::I32(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + SkyMapArray::I64(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + SkyMapArray::F32(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + SkyMapArray::F64(b) => to_png( + b.deref(), + img_size, + proj, + proj_center, + proj_bounds, + pos_convert, + color_map, + color_map_func_type, + writer, + ), + } + } +} diff --git a/test/resources/skymap/skymap.fits b/test/resources/skymap/skymap.fits new file mode 100644 index 0000000..a49acc2 Binary files /dev/null and b/test/resources/skymap/skymap.fits differ diff --git a/test/resources/skymap/skymap.png b/test/resources/skymap/skymap.png new file mode 100644 index 0000000..94bf979 Binary files /dev/null and b/test/resources/skymap/skymap.png differ