diff --git a/src/nested/gpu.rs b/src/nested/gpu.rs index 209fe87..0759996 100644 --- a/src/nested/gpu.rs +++ b/src/nested/gpu.rs @@ -54,6 +54,102 @@ pub fn proj_gpu(x: f32, y: f32, z: f32) -> (f32, f32) { } } +/// Returns the position in on the unit sphere `(x, y, z)` of the give position in the HEALPix +/// 2D Euclidean projection plane. +/// # Inputs +/// - `X`: coordinate along the X-axis in the projection plane, in `[-4, 4]` +/// - `Y`: coordinate along the Y-axis in the projection plane, in `[-2, 2]` +/// # Output: +/// - `x`: in `[-1.0, 1.0]` +/// - `y`: in `[-1.0, 1.0]` +/// - `z`: in `[-1.0, 1.0]` +/// # Remark +/// From the HPX projection as defined in Calabreta, use: +/// - `X /= PI / 4` +/// - `Y /= PI / 4` +pub fn unproj_gpu(mut x: f32, y: f32) -> (f32, f32, f32) { + use std::f32::consts::{FRAC_PI_2, FRAC_PI_4}; + const ONE_OVER_SQRT6: f32 = 0.408_248_290_463_863_f32; + assert!((-2f32..=2f32).contains(&y)); // check y + if y > 1.0 { + // North Polar Cap + assert!((-4f32..=4f32).contains(&x)); + let x = if x < 0.0 { 8.0 + x } else { x }; + debug_assert!(x > 0.0); + let OffsetAndPM1 { offset, pm1 } = pm1_offset_decompose(x); + let sqrt_of_three_time_one_minus_sin_of = 2.0 - y; + let x = if sqrt_of_three_time_one_minus_sin_of > 1e-6 { + deal_with_numerical_approx_in_edges(pm1 / sqrt_of_three_time_one_minus_sin_of) + } else { + pm1 + } + offset as f32; + // It would be faster, but less accurate, to use: + // let z = 1.0 - sqrt_of_three_time_one_minus_sin_of.pow2() / 3.0; + // let cos_lat = sqrt(1 - z^2); + let lat = 2.0 * f32::acos(sqrt_of_three_time_one_minus_sin_of * ONE_OVER_SQRT6) - FRAC_PI_2; + let (sin_lat, cos_lat) = lat.sin_cos(); + let (sin_lon, cos_lon) = (x * FRAC_PI_4).sin_cos(); + (cos_lon * cos_lat, sin_lon * cos_lat, sin_lat) + } else if y < -1.0 { + // South polar cap + let x = if x < 0.0 { 8.0 + x } else { x }; + debug_assert!(x > 0.0); + let OffsetAndPM1 { offset, pm1 } = pm1_offset_decompose(x); + let sqrt_of_three_time_one_minus_sin_of = 2.0 + y; + let x = if sqrt_of_three_time_one_minus_sin_of > 1e-6 { + deal_with_numerical_approx_in_edges(pm1 / sqrt_of_three_time_one_minus_sin_of) + } else { + pm1 + } + offset as f32; + // It would be faster, but less accurate, to use: + // let z = -1.0 + sqrt_of_three_time_one_minus_sin_of.pow2() / 3.0; + // let cos_lat = sqrt(1 - z^2); + let lat = FRAC_PI_2 - 2.0 * f32::acos(sqrt_of_three_time_one_minus_sin_of * ONE_OVER_SQRT6); + let (sin_lat, cos_lat) = lat.sin_cos(); + let (sin_lon, cos_lon) = (x * FRAC_PI_4).sin_cos(); + (cos_lon * cos_lat, sin_lon * cos_lat, sin_lat) + } else { + // Equatorial region + let z = y * TRANSITION_Z; // = sin(lat) + let cos_lat = if z < 1e-2 { + // sqrt(1 - x²) = 1 - x²/2 - x⁴/8 - x⁶/16 + let tmp = 0.5 * z * z; + 1.0 - tmp - 0.5 * tmp * tmp + } else { + (1.0 - z * z).sqrt() + }; + let (sin_lon, cos_lon) = (x * FRAC_PI_4).sin_cos(); + (cos_lon * cos_lat, sin_lon * cos_lat, z) + } +} + +// Decompose the given positive real value in +// --* an integer offset in [1, 3, 5, 7] (*PI/4) and +// --* a real value in [-1.0, 1.0] (*PI/4) +pub(crate) struct OffsetAndPM1 { + offset: u8, // = 1, 3, 5 or 7 + pm1: f32, // in [-1.0, 1.0] +} +#[inline] +pub(crate) fn pm1_offset_decompose(x: f32) -> OffsetAndPM1 { + let floor: u8 = x as u8; + let odd_floor: u8 = floor | 1u8; + OffsetAndPM1 { + offset: odd_floor & 7u8, // value modulo 8 = 1/3/5/7 + pm1: x - (odd_floor as f32), + } +} + +fn deal_with_numerical_approx_in_edges(x: f32) -> f32 { + if x > 1.0 { + 1.0 + } else if x < -1.0 { + -1.0 + } else { + x + } +} + /// Returns the cell number (hash value) associated with the given position on the unit sphere, /// together with the offset `(dx, dy)` on the Euclidean plane of the projected position with /// respect to the origin of the cell (South vertex). @@ -241,6 +337,24 @@ fn ij2z(mut i: u32, mut j: u32) -> u32 { i } +// Returns the absolute value of the given double together with its bit of sign +struct AbsAndSign { + abs: f32, + sign: u32, +} +#[inline] +pub(crate) fn abs_sign_decompose(x: f32) -> AbsAndSign { + /// Mask to keep only the f32 sign + pub const F32_SIGN_BIT_MASK: u32 = 0x80000000; + /// Equals !F32_SIGN_BIT_MASK (the inverse of the f32 sign mask) + pub const F32_BUT_SIGN_BIT_MASK: u32 = 0x7FFFFFFF; + let bits = f32::to_bits(x); + AbsAndSign { + abs: f32::from_bits(bits & F32_BUT_SIGN_BIT_MASK), + sign: bits & F32_SIGN_BIT_MASK, + } +} + #[cfg(test)] mod tests { @@ -591,10 +705,10 @@ mod tests { } } - /* No a real test, used to check the projection code visually + // No a real test, used to check the projection code visually #[test] fn testok_proj_gpu() { - println!("ra,dec,x,y,color"); + println!("ra,dec,x,y,color,ra_unproj,dec_unproj"); for ra in 0..=360 { for dec in 0..=180 { let ra = (ra as f32).to_radians(); @@ -603,8 +717,25 @@ mod tests { let y = ra.sin() * dec.cos(); let z = dec.sin(); let (x_proj, y_proj) = proj_gpu(x, y, z); - println!("{},{},{},{},{}", ra, dec, x_proj, y_proj, ra + (dec + 2.0) * 2.0 * PI); + let (x, y, z) = unproj_gpu(x_proj, y_proj); + let r = (x * x + y * y).sqrt(); + let mut ra_unproj = y.atan2(x).to_degrees(); + let dec_unproj = z.atan2(r).to_degrees(); + if ra_unproj < 0.0 { + ra_unproj += 360.0; + } + + println!( + "{},{},{},{},{},{},{}", + ra.to_degrees(), + dec.to_degrees(), + x_proj, + y_proj, + ra + (dec + 2.0) * 2.0 * PI, + ra_unproj, + dec_unproj + ); } } - }*/ + } } diff --git a/src/nested/mod.rs b/src/nested/mod.rs index e9400de..e5671e3 100644 --- a/src/nested/mod.rs +++ b/src/nested/mod.rs @@ -1446,7 +1446,7 @@ impl Layer { } else { starting_vertex.counter_clockwise_cycle() }; - // - make the five sides + // - make the four sides self.path_along_cell_side_internal( proj_center, &v1,