diff --git a/CHANGELOG.md b/CHANGELOG.md index f02990c6..cc64d0bb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* [impr] Improve `WCS` view export with 3rd euler rotation encoding: . Still some cases are to be handled like: crval on the equator or cylindrical with a galactic frame rotation. * [fixed] Change `RADECSYS` to `RADESYS` for `Aladin#getViewWCS` to follow fits standard deprecation * [feat] Add new method `Aladin#getViewImageBuffer` to get the current view as a PNG buffer * [feat] New line rasterizer using GL instancing. diff --git a/examples/al-init-custom-options.html b/examples/al-init-custom-options.html index 8e093010..5b580961 100644 --- a/examples/al-init-custom-options.html +++ b/examples/al-init-custom-options.html @@ -23,20 +23,12 @@ reticleSize: 64, // change reticle size showContextMenu: true, showCooGrid: true, + showFrame: true, } ); }); \ No newline at end of file diff --git a/src/core/src/app.rs b/src/core/src/app.rs index a3b19d43..f8ebae5a 100644 --- a/src/core/src/app.rs +++ b/src/core/src/app.rs @@ -52,7 +52,7 @@ pub struct App { //ui: GuiRef, shaders: ShaderManager, - camera: CameraViewPort, + pub camera: CameraViewPort, downloader: Downloader, tile_fetcher: TileFetcherQueue, @@ -94,7 +94,7 @@ pub struct App { colormaps: Colormaps, - projection: ProjectionType, + pub projection: ProjectionType, // Async data receivers fits_send: async_channel::Sender, @@ -849,16 +849,6 @@ impl App { Ok(has_camera_moved) } - pub(crate) fn reset_north_orientation(&mut self) { - // Reset the rotation around the center if there is one - self.camera - .set_rotation_around_center(Angle(0.0), &self.projection); - // Reset the camera position to its current position - // this will keep the current position but reset the orientation - // so that the north pole is at the top of the center. - self.set_center(&self.get_center()); - } - pub(crate) fn read_pixel(&self, pos: &Vector2, layer: &str) -> Result { if let Some(lonlat) = self.screen_to_world(pos) { if let Some(survey) = self.layers.get_hips_from_layer(layer) { @@ -1405,10 +1395,10 @@ impl App { pub(crate) fn world_to_screen(&self, ra: f64, dec: f64) -> Option> { let lonlat = LonLatT::new(ArcDeg(ra).into(), ArcDeg(dec).into()); - let model_pos_xyz = lonlat.vector(); + let icrs_pos = lonlat.vector(); self.projection - .view_to_screen_space(&model_pos_xyz, &self.camera) + .icrs_celestial_to_screen_space(&icrs_pos, &self.camera) } pub(crate) fn screen_to_world(&self, pos: &Vector2) -> Option> { @@ -1439,11 +1429,11 @@ impl App { LonLatT::new(ra, dec) } + /// lonlat must be given in icrs frame pub(crate) fn set_center(&mut self, lonlat: &LonLatT) { self.prev_cam_position = self.camera.get_center().truncate(); - self.camera - .set_center(lonlat, CooSystem::ICRS, &self.projection); + self.camera.set_center(lonlat, &self.projection); self.request_for_new_tiles = true; // And stop the current inertia as well if there is one @@ -1534,17 +1524,17 @@ impl App { self.inertia = Some(Inertia::new(ampl.to_radians(), axis)) } - pub(crate) fn rotate_around_center(&mut self, theta: ArcDeg) { + pub(crate) fn set_view_center_pos_angle(&mut self, theta: ArcDeg) { self.camera - .set_rotation_around_center(theta.into(), &self.projection); + .set_view_center_pos_angle(theta.into(), &self.projection); // New tiles can be needed and some tiles can be removed self.request_for_new_tiles = true; self.request_redraw = true; } - pub(crate) fn get_rotation_around_center(&self) -> &Angle { - self.camera.get_rotation_around_center() + pub(crate) fn get_north_shift_angle(&self) -> Angle { + self.camera.get_north_shift_angle() } pub(crate) fn set_fov(&mut self, fov: Angle) { diff --git a/src/core/src/camera/viewport.rs b/src/core/src/camera/viewport.rs index 27d0a7e9..e1356ecc 100644 --- a/src/core/src/camera/viewport.rs +++ b/src/core/src/camera/viewport.rs @@ -17,14 +17,16 @@ use crate::healpix::coverage::HEALPixCoverage; use crate::math::angle::ToAngle; use crate::math::{projection::coo_space::XYZWModel, projection::domain::sdf::ProjDef}; +use al_core::log::console_log; use cgmath::{Matrix4, Vector2}; pub struct CameraViewPort { // The field of view angle aperture: Angle, - center: Vector4, // The rotation of the camera rotation_center_angle: Angle, + center: Vector4, w2m_rot: Rotation, + center_rot: Angle, w2m: Matrix4, m2w: Matrix4, @@ -74,6 +76,7 @@ pub struct CameraViewPort { } use al_api::coo_system::CooSystem; use al_core::WebGlContext; +use js_sys::Math::atan2; use crate::{ coosys, @@ -101,8 +104,8 @@ impl CameraViewPort { let w2m = Matrix4::identity(); let m2w = w2m; - let center = Vector4::new(0.0, 0.0, 1.0, 1.0); - + let center_rot = Angle(0.0); + let center = Vector4::new(0.0, 0.0, 0.0, 1.0); let moved = false; let zoomed = false; @@ -122,9 +125,6 @@ impl CameraViewPort { let width = width * dpi; let height = height * dpi; - //let dpi = 1.0; - //gl.scissor(0, 0, width as i32, height as i32); - let aspect = height / width; let ndc_to_clip = Vector2::new(1.0, (height as f64) / (width as f64)); let clip_zoom_factor = 1.0; @@ -143,6 +143,7 @@ impl CameraViewPort { CameraViewPort { // The field of view angle aperture, + center_rot, center, // The rotation of the cameraq w2m_rot, @@ -477,10 +478,11 @@ impl CameraViewPort { self.update_rot_matrices(proj); } - pub fn set_center(&mut self, lonlat: &LonLatT, coo_sys: CooSystem, proj: &ProjectionType) { + /// lonlat must be given in icrs frame + pub fn set_center(&mut self, lonlat: &LonLatT, proj: &ProjectionType) { let icrs_pos: Vector4<_> = lonlat.vector(); - let view_pos = coosys::apply_coo_system(coo_sys, self.get_coo_system(), &icrs_pos); + let view_pos = CooSystem::ICRS.to(self.get_coo_system()) * icrs_pos; let rot = Rotation::from_sky_position(&view_pos); // Apply the rotation to the camera to go @@ -527,13 +529,8 @@ impl CameraViewPort { if self.reversed_longitude != reversed_longitude { self.reversed_longitude = reversed_longitude; - self.rotation_center_angle = -self.rotation_center_angle; self.update_rot_matrices(proj); } - - // The camera is reversed => it has moved - self.moved = true; - self.time_last_move = Time::now(); } pub fn get_longitude_reversed(&self) -> bool { @@ -599,14 +596,17 @@ impl CameraViewPort { self.zoomed = false; } + #[inline] pub fn get_aperture(&self) -> Angle { self.aperture } + #[inline] pub fn get_center(&self) -> &Vector4 { &self.center } + #[inline] pub fn is_allsky(&self) -> bool { self.is_allsky } @@ -619,13 +619,14 @@ impl CameraViewPort { self.coo_sys } - pub fn set_rotation_around_center(&mut self, theta: Angle, proj: &ProjectionType) { - self.rotation_center_angle = theta; + pub fn set_view_center_pos_angle(&mut self, phi: Angle, proj: &ProjectionType) { + self.center_rot = phi; + self.update_rot_matrices(proj); } - pub fn get_rotation_around_center(&self) -> &Angle { - &self.rotation_center_angle + pub fn get_north_shift_angle(&self) -> Angle { + (self.w2m.x.y).atan2(self.w2m.y.y).to_angle() } } use crate::ProjectionType; @@ -657,16 +658,14 @@ impl CameraViewPort { } fn update_center(&mut self) { - // The center position is on the 3rd column of the w2m matrix self.center = self.w2m.z; - - let axis = &self.center.truncate(); - let center_rot = Rotation::from_axis_angle(axis, self.rotation_center_angle); - + // The center position is on the 3rd column of the w2m matrix + let center_axis = &self.center.truncate(); // Re-update the model matrix to take into account the rotation // by theta around the center axis - let final_rot = center_rot * self.w2m_rot; - self.w2m = (&final_rot).into(); + let r = Rotation::from_axis_angle(center_axis, self.center_rot) * self.w2m_rot; + + self.w2m = (&r).into(); if self.reversed_longitude { self.w2m = self.w2m * ID_R; } diff --git a/src/core/src/lib.rs b/src/core/src/lib.rs index 4c0a3ef1..959d168e 100644 --- a/src/core/src/lib.rs +++ b/src/core/src/lib.rs @@ -494,20 +494,30 @@ impl WebClient { /// # Arguments /// /// * `theta` - The rotation angle in degrees - #[wasm_bindgen(js_name = setRotationAroundCenter)] - pub fn rotate_around_center(&mut self, theta: f64) -> Result<(), JsValue> { + #[wasm_bindgen(js_name = setViewCenterPosAngle)] + pub fn set_view_center_pos_angle(&mut self, theta: f64) -> Result<(), JsValue> { let theta = ArcDeg(theta); - self.app.rotate_around_center(theta); + self.app.set_view_center_pos_angle(theta); Ok(()) } /// Get the absolute orientation angle of the view - #[wasm_bindgen(js_name = getRotationAroundCenter)] - pub fn get_rotation_around_center(&mut self) -> Result { - let theta = self.app.get_rotation_around_center(); + #[wasm_bindgen(js_name = getViewCenterFromNorthPoleAngle)] + pub fn get_north_shift_angle(&mut self) -> Result { + let phi = self.app.get_north_shift_angle(); + Ok(phi.to_degrees()) + } + + #[wasm_bindgen(js_name = getNorthPoleCelestialPosition)] + pub fn get_north_pole_celestial_position(&mut self) -> Result, JsValue> { + let np = self + .app + .projection + .north_pole_celestial_space(&self.app.camera); - Ok(theta.0 * 360.0 / (2.0 * std::f64::consts::PI)) + let (lon, lat) = (np.lon().to_degrees(), np.lat().to_degrees()); + Ok(Box::new([lon, lat])) } /// Get if the longitude axis is reversed @@ -571,12 +581,6 @@ impl WebClient { Ok(Box::new([lon_deg.0, lat_deg.0])) } - /// Rest the north pole orientation to the top of the screen - #[wasm_bindgen(js_name = resetNorthOrientation)] - pub fn reset_north_orientation(&mut self) { - self.app.reset_north_orientation(); - } - /// Go from a location to another one /// /// # Arguments diff --git a/src/core/src/math/projection/mod.rs b/src/core/src/math/projection/mod.rs index 559fd93a..a64eae80 100644 --- a/src/core/src/math/projection/mod.rs +++ b/src/core/src/math/projection/mod.rs @@ -150,7 +150,16 @@ pub enum ProjectionType { //Hpx(mapproj::hybrid::hpx::Hpx), } +use crate::math::lonlat::LonLat; impl ProjectionType { + pub fn north_pole_celestial_space(&self, camera: &CameraViewPort) -> LonLatT { + // This is always defined + let np_world = self.north_pole_world_space(); + + let np_celestial = camera.get_w2m() * np_world; + np_celestial.lonlat() + } + /// Screen to model space deprojection /// Perform a screen to the world space deprojection @@ -209,41 +218,28 @@ impl ProjectionType { self.world_to_screen_space(&pos_world_space, camera) } - pub fn view_to_screen_space( + pub fn icrs_celestial_to_screen_space( &self, - pos_model_space: &XYZWModel, + icrs_celestial_pos: &XYZWModel, camera: &CameraViewPort, ) -> Option> { - self.view_to_normalized_device_space(pos_model_space, camera) + self.icrs_celestial_to_normalized_device_space(icrs_celestial_pos, camera) .map(|ndc_pos| crate::ndc_to_screen_space(&ndc_pos, camera)) } - pub fn view_to_normalized_device_space( + pub fn icrs_celestial_to_normalized_device_space( &self, - pos_model_space: &XYZWModel, + icrs_celestial_pos: &XYZWModel, camera: &CameraViewPort, ) -> Option> { let view_coosys = camera.get_coo_system(); let c = CooSystem::ICRS.to::(view_coosys); let m2w = camera.get_m2w(); - let pos_world_space = m2w * c * pos_model_space; + let pos_world_space = m2w * c * icrs_celestial_pos; self.world_to_normalized_device_space(&pos_world_space, camera) } - /*pub fn view_to_normalized_device_space_unchecked( - &self, - pos_view_space: &Vector4, - camera: &CameraViewPort, - ) -> Vector2 { - let view_coosys = camera.get_coo_system(); - let c = CooSystem::ICRS.to::(view_coosys); - - let m2w = camera.get_m2w(); - let pos_world_space = m2w * c * pos_view_space; - self.world_to_normalized_device_space_unchecked(&pos_world_space, camera) - }*/ - pub fn model_to_normalized_device_space( &self, pos_model_space: &XYZWModel, @@ -673,6 +669,21 @@ pub trait Projection { /// /// * ``pos_world_space`` - The position in the world space fn world_to_clip_space(&self, pos_world_space: &XYZWWorld) -> Option>; + + /// (`alpha_p`, `delta_p`) in the WCS II paper from Mark Calabretta. + #[inline] + fn north_pole_world_space(&self) -> XYZWWorld { + // This is always defined + self.clip_to_world_space(&XYClip::new(0.0, 1.0 - 1e-5)) + .unwrap() + } + + #[inline] + fn south_pole_world_space(&self) -> XYZWWorld { + // This is always defined + self.clip_to_world_space(&XYClip::new(0.0, -1.0 + 1e-5)) + .unwrap() + } } use mapproj::ProjXY; @@ -680,6 +691,8 @@ use mapproj::ProjXY; use self::coo_space::XYScreen; use self::coo_space::XYNDC; +use super::lonlat::LonLatT; + impl<'a, P> Projection for &'a P where P: CanonicalProjection, diff --git a/src/core/src/math/rotation.rs b/src/core/src/math/rotation.rs index d20e4690..137d4747 100644 --- a/src/core/src/math/rotation.rs +++ b/src/core/src/math/rotation.rs @@ -1,7 +1,8 @@ use crate::math; -use cgmath::Quaternion; use cgmath::{BaseFloat, InnerSpace}; +use cgmath::{Euler, Quaternion}; use cgmath::{Vector3, Vector4}; +use core::f64::consts::PI; #[derive(Clone, Copy, Debug)] // Internal structure of a rotation, a quaternion @@ -136,6 +137,31 @@ where m2w * pos_model_space } + + pub fn euler(&self) -> Euler> { + self.0.into() + } + + /// Extract the 3 euler angles from the quaternion + /// Aladin Lite rotation basis is formed by Z, X, Y axis: + /// * Z axis is pointing towards us + /// * Y is pointing upward + /// * X is defined from the right-hand rule to form a basis + /// + /// The first euler angle describes the longitude (rotation around the Y axis) <=> pitch + /// The second euler angle describes the latitude (rotation around the X' modified axis) <=> yaw + /// The third euler angle describes a rotation deviation from the north pole (rotation around the Z'' modified axis) <=> roll + /// + /// Equations come from this paper (Appendix 6): + /// https://ntrs.nasa.gov/api/citations/19770024290/downloads/19770024290.pdf + pub fn euler_YXZ(&self) -> (Angle, Angle, Angle) { + let m: Matrix4 = self.0.into(); + + let a = m.x.z.atan2(m.z.z); + let b = (-m.z.y).atan2((S::one() - m.z.y * m.z.y).sqrt()); + let c = m.x.y.atan2(m.y.y); + (Angle(a), Angle(b), Angle(c)) + } } use std::ops::Mul; diff --git a/src/core/src/renderable/hips/mod.rs b/src/core/src/renderable/hips/mod.rs index 8f042e95..db2d1797 100644 --- a/src/core/src/renderable/hips/mod.rs +++ b/src/core/src/renderable/hips/mod.rs @@ -58,7 +58,7 @@ fn is_too_large(cell: &HEALPixCell, camera: &CameraViewPort, projection: &Projec .iter() .filter_map(|(lon, lat)| { let vertex = crate::math::lonlat::radec_to_xyzw(Angle(*lon), Angle(*lat)); - projection.view_to_screen_space(&vertex, camera) + projection.icrs_celestial_to_screen_space(&vertex, camera) }) .collect::>(); diff --git a/src/js/Aladin.js b/src/js/Aladin.js index 6edb986e..16cd6312 100644 --- a/src/js/Aladin.js +++ b/src/js/Aladin.js @@ -1867,8 +1867,26 @@ export let Aladin = (function () { this.view.decreaseZoom(0.01); }; - Aladin.prototype.setRotation = function (rotation) { - this.view.setRotation(rotation); + /** + * Set the view center rotation in degrees + * + * @memberof Aladin + * @param {number} rotation - The center rotation in degrees. Positive angles rotates the + * view in the counter clockwise order (or towards the east) + */ + Aladin.prototype.setViewCenterPosAngle = function (rotation) { + this.view.setViewCenterPosAngle(rotation); + }; + + /** + * Get the view center to north pole angle in degrees. This is equivalent to getting the 3rd Euler angle + * + * @memberof Aladin + * + * @returns {number} - Angle between the position center and the north pole + */ + Aladin.prototype.getCenter2NorthPoleAngle = function () { + return this.view.wasm.getViewCenterFromNorthPoleAngle(); }; // @api @@ -2264,9 +2282,8 @@ aladin.on("positionChanged", ({ra, dec}) => { * Return the current view WCS as a key-value dictionary * Can be useful in coordination with getViewDataURL * - * NOTE + TODO : Rotations are not implemented yet - * - * @API + * @memberof Aladin + * @returns {Object} - A JS object describing the WCS of the view. */ Aladin.prototype.getViewWCS = function () { // get general view properties @@ -2276,11 +2293,11 @@ aladin.on("positionChanged", ({ra, dec}) => { const height = this.view.height; // get values common for all - let cdelt1 = fov[0] / width; + let cdelt1 = -fov[0] / width; const cdelt2 = fov[1] / height; - const projectionName = this.getProjectionName(); + const projName = this.getProjectionName(); - if (projectionName == "FEYE") + if (projName == "FEYE") return "Fish eye projection is not supported by WCS standards."; // reversed longitude case @@ -2339,8 +2356,8 @@ aladin.on("positionChanged", ({ra, dec}) => { CRPIX2: height / 2 + 0.5, CRVAL1: center[0], CRVAL2: center[1], - CTYPE1: cooType1 + projectionName, - CTYPE2: cooType2 + projectionName, + CTYPE1: cooType1 + projName, + CTYPE2: cooType2 + projName, CUNIT1: "deg ", CUNIT2: "deg ", CDELT1: cdelt1, @@ -2351,6 +2368,48 @@ aladin.on("positionChanged", ({ra, dec}) => { // the radecsys keyword if (radesys == "ICRS ") WCS.RADESYS = radesys; + const isProjZenithal = ['TAN', 'SIN', 'STG', 'ZEA'].some((p) => p === projName) + if (isProjZenithal) { + // zenithal projections + // express the 3rd euler angle for zenithal projection + let thirdEulerAngle = this.getCenter2NorthPoleAngle(); + WCS.LONPOLE = 180 - thirdEulerAngle + } else { + // cylindrical or pseudo-cylindrical projections + if (!radesys) { + console.warn('TODO: Galactic frame 3rd euler rotation not handled for cylindrical projections') + } else if (WCS.CRVAL2 === 0) { + // ref point on the equator not handled (yet) + console.warn('TODO: 3rd euler rotation is not handled for ref point located at delta_0 = 0') + } else { + // ref point not on the equator + const npLonlat = this.view.wasm.getNorthPoleCelestialPosition(); + let dLon = WCS.CRVAL1 - npLonlat[0]; + + // dlon angle must lie between -PI and PI + // For dlon angle between -PI;-PI/2 or PI/2;PI one must invert LATPOLE + if (dLon < -90 || dLon > 90) { + // so that the south pole becomes upward to the ref point + WCS.LATPOLE = -90 + } + + const toRad = Math.PI / 180 + const toDeg = 1.0 / toRad; + + // Reverse the Eq 9 from the WCS II paper from Mark Calabretta to obtain LONPOLE + // function of CRVAL2 and native coordinates of the fiducial ref point, i.e. (phi_0, theta_0) = (0, 0) + // for cylindrical projections + WCS.LONPOLE = Math.asin(Math.sin(dLon * toRad) * Math.cos(WCS.CRVAL2 * toRad)) * toDeg; + + if (WCS.CRVAL2 < 0) { + // ref point is located in the south hemisphere + WCS.LONPOLE = -180 - WCS.LONPOLE; + } + + console.log("euler rot ?:", 180 - WCS.LONPOLE) + } + } + return WCS; }; diff --git a/src/js/View.js b/src/js/View.js index 09b5a411..7139180a 100644 --- a/src/js/View.js +++ b/src/js/View.js @@ -546,6 +546,8 @@ export let View = (function () { try { const [lon, lat] = view.aladin.pix2world(xymouse.x, xymouse.y, 'icrs'); view.pointTo(lon, lat); + // reset the rotation around center view + view.setViewCenterPosAngle(0.0); } catch (err) { return; @@ -709,7 +711,7 @@ export let View = (function () { view.pinchZoomParameters.initialFov = fov; view.pinchZoomParameters.initialDistance = Math.sqrt(Math.pow(e.targetTouches[0].clientX - e.targetTouches[1].clientX, 2) + Math.pow(e.targetTouches[0].clientY - e.targetTouches[1].clientY, 2)); - view.fingersRotationParameters.initialViewAngleFromCenter = view.wasm.getRotationAroundCenter(); + view.fingersRotationParameters.initialViewAngleFromCenter = view.wasm.getNorthShiftAngle(); view.fingersRotationParameters.initialFingerAngle = Math.atan2(e.targetTouches[1].clientY - e.targetTouches[0].clientY, e.targetTouches[1].clientX - e.targetTouches[0].clientX) * 180.0 / Math.PI; return; @@ -973,7 +975,7 @@ export let View = (function () { // planetary survey case rotation -= fingerAngleDiff; } - view.wasm.setRotationAroundCenter(rotation); + view.setViewCenterPosAngle(rotation); } // zoom @@ -1577,8 +1579,8 @@ export let View = (function () { }); } - View.prototype.setRotation = function(rotation) { - this.wasm.setRotationAroundCenter(rotation); + View.prototype.setViewCenterPosAngle = function(rotation) { + this.wasm.setViewCenterPosAngle(rotation); } View.prototype.setGridOptions = function (options) {