diff --git a/Cargo.toml b/Cargo.toml index 4fba42c..c84398b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ [package] name = "moc" -version = "0.11.3" +version = "0.12.0" authors = [ "F.-X. Pineau ", "Matthieu Baumann " @@ -59,6 +59,9 @@ rayon = "1.5.1" flate2 = "1.0.22" # Compression/decompression png = "0.17" # For png slab = { version = "0.4", optional = true } +stc = { git = "https://github.com/cds-astro/cds-stc-rust" } +thiserror = "1.0" # For STC-S (quick-error should be replaced by thiserror everywhere) + [dev-dependencies] rand = "0.8.3" diff --git a/README.md b/README.md index 9de0526..3643bb3 100644 --- a/README.md +++ b/README.md @@ -114,6 +114,40 @@ pub fn or( whether the index is a single index, a range lower bound or a range upper bound * [ ] Make a PostgresQL wrapper using e.g. [pgx](https://github.com/zombodb/pgx/)? + +## WARNING about the STC-S to MOC function + +STC-S parsing is ensured by the [STC crate](https://github.com/cds-astro/cds-stc-rust). + +Current discrepancies between the STC standard and this implementation: + +* The `DIFFERENCE` operation has been implemented as being a `symmetric difference` + + why? probably because: + 1. I am biased towards Boolean algebra, it as `XOR` + (exclusive `OR` or symmetric difference) but no `Difference` + 2. I read parts of the STC standard after the STC-S implementation + 3. `XOR` is already implemented in [cdshleapix](https://github.com/cds-astro/cds-healpix-rust), but `DIFFERENCE` is not. + + has stated in the STC standard: `R1 – R2 = R1 AND (NOT R2))`; + but also: `R1 - R2 = R1 AND (R1 XOR R2)`, and + `XOR = (R1 OR R2) AND (NOT (R1 AND R2))` is more complex that `DIFFERENCE` + (so is worth having implented?). +* For `Polygon`: we do not use the STC convention + + we support self-intersecting polygons + + we generally return the smallest area polygon (use `NOT` to get its complement!) + + one convention could be to use an additional (last) provided points as a control point + - note that for convex polygons, the control point could be the vertices gravity center + - in a GUI, a user could define the inner part of the polygon by a final click + + why? + 1. efficient algorithms dealing with polygons supports self-intersecting polygons + 2. to support arbitrary defined polygons by a user clicking in a viewer such as Aladin or Aladin Lite + 3. [cdshleapix](https://github.com/cds-astro/cds-healpix-rust) is based on self-intersecting polygons +* For `Box`: a position angle can be added as a last parameter, right after `bsize`. + +So far, we reject STC-S having: +* a frame different from `ICRS` +* a flavor different from `Spher2` +* units different from `degrees` + ## License Like most projects in Rust, this project is licensed under either of diff --git a/src/deser/mod.rs b/src/deser/mod.rs index c2e3ea1..6482662 100644 --- a/src/deser/mod.rs +++ b/src/deser/mod.rs @@ -13,3 +13,4 @@ pub mod json; pub mod fits; pub mod gz; pub mod img; +pub mod stcs; diff --git a/src/deser/stcs.rs b/src/deser/stcs.rs new file mode 100644 index 0000000..8b62467 --- /dev/null +++ b/src/deser/stcs.rs @@ -0,0 +1,415 @@ +//! MOC creation from an STC-S string. + +use thiserror::Error; + +use nom::{ + error::{convert_error, VerboseError}, + Err, +}; + +use stc::{ + Stc, + space::{ + common::{ + region::{BoxParams, CircleParams, ConvexParams, EllipseParams, PolygonParams}, + FillFrameRefposFlavor, Flavor, Frame, FromPosToVelocity, SpaceUnit, + }, + position::Position, + positioninterval::PositionInterval, + }, + visitor::{StcVisitResult, CompoundVisitor, SpaceVisitor, impls::donothing::VoidVisitor}, +}; + +use healpix::nested::{ + bmoc::{BMOCBuilderUnsafe, BMOC}, + box_coverage, cone_coverage_approx_custom, elliptical_cone_coverage_custom, polygon_coverage, +}; + +use crate::{moc::range::RangeMOC, qty::Hpx}; + +const HALF_PI: f64 = 0.5 * std::f64::consts::PI; +const PI: f64 = std::f64::consts::PI; +const TWICE_PI: f64 = 2.0 * std::f64::consts::PI; + +#[derive(Error, Debug)] +pub enum Stc2MocError { + #[error("Frame other than ICRS not supported (yet). Found: {found:?}")] + FrameIsNotICRS { found: Frame }, + #[error("Flavor other than Spher2 not supported (yet). Found: {found:?}")] + FlavorIsNotSpher2 { found: Flavor }, + #[error("Units ther than 'deg' not (yet?!) supported. Found: {found:?}")] + UnitsNotSupported { found: Vec }, + #[error("Convex shape not (yet?!) supported.")] + ConvexNotSupported, + #[error("Simple position not supported.")] + SimplePositionNotSupported, + #[error("Position interval not supported.")] + PositionIntervalNotSupported, + #[error("invalid header (expected {expected:?}, found {found:?})")] + WrongNumberOfParams { expected: u8, found: u8 }, + #[error("Longitude value out of bounds. Expected: [0, 360[. Actual: {value:?}")] + WrongLongitude { value: f64 }, + #[error("Latitude value out of bounds. Expected: [-90, 90[. Actual: {value:?}")] + WrongLatitude { value: f64 }, + #[error("STC-S string parsing not complete. Remaining: {rem:?}")] + ParseHasRemaining { rem: String }, + #[error("STC-S string parsing incomplete: {msg:?}")] + ParseIncomplete { msg: String }, + #[error("STC-S string parsing error: {msg:?}")] + ParseFailure { msg: String }, + #[error("STC-S string parsing failure: {msg:?}")] + ParseError { msg: String }, + #[error("No space sub-phrase found in STC-S string")] + NoSpaceFound, + #[error("Custom error: {msg:?}")] + Custom { msg: String }, +} + +#[derive(Debug, Clone)] +struct Stc2Moc { + depth: u8, + delta_depth: u8, +} +impl Stc2Moc { + fn new(depth: u8, delta_depth: Option) -> Self { + Self { depth, delta_depth: delta_depth.unwrap_or(2) } + } +} +impl CompoundVisitor for Stc2Moc { + type Value = BMOC; + type Error = Stc2MocError; + + fn visit_allsky(&mut self) -> Result { + // Ok(BMOC::new_allsky(self.depth)) + Ok(new_allsky(self.depth)) + } + + fn visit_circle(&mut self, circle: &CircleParams) -> Result { + // Get params + let lon_deg = circle.center().get(0).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty circle longitude"), + })?; + let lat_deg = circle.center().get(1).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty circle latitude"), + })?; + let radius_deg = circle.radius(); + // Convert params + let lon = lon_deg2rad(*lon_deg)?; + let lat = lat_deg2rad(*lat_deg)?; + let r = radius_deg.to_radians(); + if r <= 0.0 || PI <= r { + Err(Stc2MocError::Custom { + msg: format!("Radius out of bounds. Expected: ]0, 180[. Actual: {}.", r), + }) + } else { + Ok(cone_coverage_approx_custom( + self.depth, + self.delta_depth, + lon, + lat, + r, + )) + } + } + + fn visit_ellipse(&mut self, ellipse: &EllipseParams) -> Result { + // Get params + let lon_deg = ellipse + .center() + .get(0) + .ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty ellipse longitude"), + })?; + let lat_deg = ellipse + .center() + .get(1) + .ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty ellipse latitude"), + })?; + let a_deg = ellipse.radius_a(); + let b_deg = ellipse.radius_b(); + let pa_deg = ellipse.pos_angle(); + // Convert params + let lon = lon_deg2rad(*lon_deg)?; + let lat = lat_deg2rad(*lat_deg)?; + let a = a_deg.to_radians(); + let b = b_deg.to_radians(); + let pa = pa_deg.to_radians(); + if a <= 0.0 || HALF_PI <= a { + Err(Stc2MocError::Custom { + msg: format!( + "Semi-major axis out of bounds. Expected: ]0, 90[. Actual: {}.", + a_deg + ), + }) + } else if b <= 0.0 || a <= b { + Err(Stc2MocError::Custom { + msg: format!( + "Semi-minor axis out of bounds. Expected: ]0, {}[. Actual: {}.", + a_deg, b_deg + ), + }) + } else if pa <= 0.0 || PI <= pa { + Err(Stc2MocError::Custom { + msg: format!( + "Position angle out of bounds. Expected: [0, 180[. Actual: {}.", + pa_deg + ), + }) + } else { + Ok(elliptical_cone_coverage_custom( + self.depth, + self.delta_depth, + lon, + lat, + a, + b, + pa, + )) + } + } + + fn visit_box(&mut self, skybox: &BoxParams) -> Result { + // Get params + let lon_deg = skybox.center().get(0).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty ellipse longitude"), + })?; + let lat_deg = skybox.center().get(1).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty ellipse latitude"), + })?; + let mut a_deg = skybox.bsize().get(0).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty bsize on latitude"), + })?; + let mut b_deg = skybox.bsize().get(0).ok_or_else(|| Stc2MocError::Custom { + msg: String::from("Empty bsize on longitude"), + })?; + let mut pa_deg = skybox.bsize().get(0).copied().unwrap_or(90.0); + if a_deg < b_deg { + std::mem::swap(&mut b_deg, &mut a_deg); + pa_deg = 90.0 - pa_deg; + } + // Convert params + let lon = lon_deg2rad(*lon_deg)?; + let lat = lat_deg2rad(*lat_deg)?; + let a = a_deg.to_radians(); + let b = b_deg.to_radians(); + let pa = pa_deg.to_radians(); + if a <= 0.0 || HALF_PI <= a { + Err(Stc2MocError::Custom { + msg: format!( + "Box semi-major axis out of bounds. Expected: ]0, 90[. Actual: {}.", + a_deg + ), + }) + } else if b <= 0.0 || a <= b { + Err(Stc2MocError::Custom { + msg: format!( + "Box semi-minor axis out of bounds. Expected: ]0, {}[. Actual: {}.", + a_deg, b_deg + ), + }) + } else if !(0.0..PI).contains(&pa) { + Err(Stc2MocError::Custom { + msg: format!( + "Position angle out of bounds. Expected: [0, 180[. Actual: {}.", + pa_deg + ), + }) + } else { + Ok(box_coverage(self.depth, lon, lat, a, b, pa)) + } + } + + fn visit_polygon(&mut self, polygon: &PolygonParams) -> Result { + let vertices_deg = polygon.vertices(); + let vertices = vertices_deg + .iter() + .step_by(2) + .zip(vertices_deg.iter().skip(1).step_by(2)) + .map(|(lon_deg, lat_deg)| { + let lon = lon_deg2rad(*lon_deg)?; + let lat = lat_deg2rad(*lat_deg)?; + Ok((lon, lat)) + }) + .collect::, Stc2MocError>>()?; + Ok(polygon_coverage(self.depth, vertices.as_slice(), true)) + } + + fn visit_convex(&mut self, _convex: &ConvexParams) -> Result { + Err(Stc2MocError::ConvexNotSupported) + } + + fn visit_not(&mut self, bmoc: Self::Value) -> Result { + Ok(bmoc.not()) + } + + fn visit_union(&mut self, bmocs: Vec) -> Result { + let n = bmocs.len(); + bmocs + .into_iter() + .reduce(|acc, curr| acc.or(&curr)) + .ok_or_else(|| Stc2MocError::Custom { + msg: format!( + "Wrong number of elements in union. Expected: >=2. Actual: {} ", + n + ), + }) + } + + fn visit_intersection(&mut self, bmocs: Vec) -> Result { + let n = bmocs.len(); + bmocs + .into_iter() + .reduce(|acc, curr| acc.and(&curr)) + .ok_or_else(|| Stc2MocError::Custom { + msg: format!( + "Wrong number of elements in intersection. Expected: >=2. Actual: {} ", + n + ), + }) + } + + fn visit_difference( + &mut self, + left_bmoc: Self::Value, + right_bmoc: Self::Value, + ) -> Result { + // Warning: we interpret 'difference' as being a 'symmetrical difference', i.e. xor (not minus) + Ok(left_bmoc.xor(&right_bmoc)) + } +} + +impl SpaceVisitor for Stc2Moc { + type Value = RangeMOC>; + type Error = Stc2MocError; + type C = Self; + + fn new_compound_visitor( + &self, + fill_frame_refpos_flavor: &FillFrameRefposFlavor, + from_pos_to_velocity: &FromPosToVelocity, + ) -> Result { + // Check ICRS frame + let frame = fill_frame_refpos_flavor.frame(); + if frame != Frame::ICRS { + return Err(Stc2MocError::FrameIsNotICRS { found: frame }); + } + // Check SPHER2 flavor + let flavor = fill_frame_refpos_flavor.flavor(); + if let Some(flavor) = flavor { + if flavor != Flavor::Spher2 { + return Err(Stc2MocError::FlavorIsNotSpher2 { found: flavor }); + } + } + // Check units + let opt_units = from_pos_to_velocity.unit().cloned(); + if let Some(units) = opt_units { + for unit in units.iter().cloned() { + if unit != SpaceUnit::Deg { + return Err(Stc2MocError::UnitsNotSupported { found: units }); + } + } + } + Ok(self.clone()) + } + + fn visit_position_simple(self, _: &Position) -> Result { + Err(Stc2MocError::SimplePositionNotSupported) + } + + fn visit_position_interval(self, _: &PositionInterval) -> Result { + Err(Stc2MocError::PositionIntervalNotSupported) + } + + fn visit_allsky(self, bmoc: BMOC) -> Result { + Ok(Self::Value::from(bmoc)) + } + + fn visit_circle(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_ellipse(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_box(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_polygon(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_convex(self, _: BMOC) -> Result { + unreachable!() // because an error is raised before calling this + } + + fn visit_not(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_union(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_intersection(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } + + fn visit_difference(self, bmoc: BMOC) -> Result { + Ok(bmoc.into()) + } +} + +fn lon_deg2rad(lon_deg: f64) -> Result { + let mut lon = lon_deg.to_radians(); + if lon == TWICE_PI { + lon = 0.0; + } + if !(0.0..TWICE_PI).contains(&lon) { + Err(Stc2MocError::WrongLongitude { value: lon_deg }) + } else { + Ok(lon) + } +} + +fn lat_deg2rad(lat_deg: f64) -> Result { + let lat = lat_deg.to_radians(); + if !(-HALF_PI..=HALF_PI).contains(&lat) { + Err(Stc2MocError::WrongLatitude { value: lat_deg }) + } else { + Ok(lat) + } +} + +// TODO: remove when newt version of cdshealpix will be published +fn new_allsky(depth: u8) -> BMOC { + let mut builder = BMOCBuilderUnsafe::new(depth, 12); + builder.push_all(0, 0, 11, true); + builder.to_bmoc() +} + + +pub fn stcs2moc(depth: u8, delta_depth: Option, stcs_ascii: &str) -> Result>, Stc2MocError> { + match Stc::parse::>(stcs_ascii.trim()) { + Ok((rem, stcs)) => { + if !rem.is_empty() { + return Err(Stc2MocError::ParseHasRemaining { rem: rem.to_string() }) + } + let stc2moc_visitor = Stc2Moc::new(depth, delta_depth); + let StcVisitResult { space, .. } = stcs.accept(VoidVisitor, stc2moc_visitor, VoidVisitor, VoidVisitor); + match space { + None => Err(Stc2MocError::NoSpaceFound), + Some(space_res) => space_res, + } + } + Err(err) => { + Err(match err { + Err::Incomplete(_) => Stc2MocError::ParseIncomplete { msg: String::from("Incomplete parsing.") }, + Err::Error(e) => Stc2MocError::ParseIncomplete { msg: convert_error(stcs_ascii, e) }, + Err::Failure(e) => Stc2MocError::ParseIncomplete { msg: convert_error(stcs_ascii, e) }, + }) + } + } +}