gitea-mirror/aladin-lite
1
0
Fork 0
mirror of https://github.com/cds-astro/aladin-lite.git synced 2025-12-12 15:49:18 -08:00
Code Issues Packages Projects Releases Wiki Activity

enhance grid wip

Browse Source
This commit is contained in:
Matthieu Baumann
2023-06-08 16:49:25 +02:00
parent 163dd7d762
commit 526cf51c4c
8 changed files with 557 additions and 450 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/core/src/math/projection/coo_space.rs
View File

@@ -9,4 +9,5 @@ pub type XYNDC = Vector2<f64>;
pub type XYClip = Vector2<f64>;
pub type XYZWorld = Vector3<f64>;
pub type XYZWWorld = Vector4<f64>;
pub type XYZWModel = Vector4<f64>;
pub type XYZWModel = Vector4<f64>;
pub type XYZModel = Vector3<f64>;

58
src/core/src/math/sph_geom/mod.rs
View File

@@ -1,2 +1,60 @@
pub mod bbox;
pub mod region;
use super::{PI, TWICE_PI};
#[inline]
pub fn is_in_lon_range(lon0: f64, lon1: f64, lon2: f64) -> bool {
// First version of the code:
// ((v2.lon() - v1.lon()).abs() > PI) != ((v2.lon() > coo.lon()) != (v1.lon() > coo.lon()))
//
// Lets note
// - lonA = v1.lon()
// - lonB = v2.lon()
// - lon0 = coo.lon()
// When (lonB - lonA).abs() <= PI
// => lonB > lon0 != lonA > lon0 like in PNPOLY
// A B lonA <= lon0 && lon0 < lonB
// --[++++[--
// B A lonB <= lon0 && lon0 < lonA
//
// But when (lonB - lonA).abs() > PI, then the test should be
// => lonA >= lon0 == lonB >= lon0
// <=> !(lonA >= lon0 != lonB >= lon0)
// A | B (lon0 < lonB) || (lonA <= lon0)
// --[++|++[--
// B | A (lon0 < lonA) || (lonB <= lon0)
//
// Instead of lonA > lon0 == lonB > lon0,
// i.e. !(lonA > lon0 != lonB > lon0).
// A | B (lon0 <= lonB) || (lonA < lon0)
// --]++|++]--
// B | A (lon0 <= lonA) || (lonB < lon0)
//
// So the previous code was bugged in this very specific case:
// - `lon0` has the same value as a vertex being part of:
// - one segment that do not cross RA=0
// - plus one segment crossing RA=0.
// - the point have an odd number of intersections with the polygon
// (since it will be counted 0 or 2 times instead of 1).
let dlon = lon2 - lon1;
if dlon < 0.0 {
(dlon >= -PI) == (lon2 <= lon0 && lon0 < lon1)
} else {
(dlon <= PI) == (lon1 <= lon0 && lon0 < lon2)
}
}
// Returns the longitude length between two longitudes
// lon1 lies in [0; 2\pi[
// lon2 lies in [0; 2\pi[
#[inline]
pub fn distance_from_two_lon(lon1: f64, lon2: f64) -> f64 {
// Cross the primary meridian
if lon1 > lon2 {
lon2 + TWICE_PI - lon1
} else {
lon2 - lon1
}
}

14
src/core/src/math/sph_geom/region.rs
View File

@@ -6,6 +6,7 @@ use cdshealpix::sph_geom::coo3d::Coo3D;
use cdshealpix::sph_geom::coo3d::Vec3;
use cdshealpix::sph_geom::ContainsSouthPoleMethod;
use crate::math::angle::ToAngle;
use al_core::{inforec, info, log};
pub enum Region {
AllSky,
@@ -30,6 +31,7 @@ pub enum PoleContained {
}
#[derive(Debug)]
pub enum Intersection {
// The segment is fully included into the region
Included,
@@ -87,13 +89,19 @@ impl Region {
}
pub fn intersects_parallel(&self, lat: f64) -> Intersection {
use crate::math::lonlat::LonLat;
match self {
// The polygon is included inside the region
Region::AllSky => Intersection::Included,
Region::Polygon { polygon, .. } => {
let vertices = polygon.intersect_parallel(lat)
.iter()
.map(|v| XYZWModel::new(v.y(), v.z(), v.x(), 1.0))
.map(|v| {
let v = XYZWModel::new(v.y(), v.z(), v.x(), 1.0);
let l = v.lonlat();
al_core::info!(l);
v
})
.collect::<Vec<_>>();
if !vertices.is_empty() {
@@ -134,8 +142,8 @@ impl Region {
}
pub fn intersects_meridian(&self, lon: f64) -> Intersection {
let n_pole_lonlat = LonLatT::new(HALF_PI.to_angle(), lon.to_angle());
let s_pole_lonlat = LonLatT::new(MINUS_HALF_PI.to_angle(), lon.to_angle());
let n_pole_lonlat = LonLatT::new(lon.to_angle(), (HALF_PI - 1e-4).to_angle());
let s_pole_lonlat = LonLatT::new(lon.to_angle(), (MINUS_HALF_PI + 1e-4).to_angle());
self.intersects_great_circle_arc(&s_pole_lonlat, &n_pole_lonlat)
}

457
src/core/src/renderable/grid/mod.rs
View File

@@ -1,9 +1,13 @@
use crate::math::MINUS_HALF_PI;
use crate::math::lonlat::LonLat;
use crate::math::projection::coo_space::XYNDC;
use crate::math::projection::coo_space::XYScreen;
use crate::math::projection::coo_space::XYZWModel;
use crate::math::sph_geom::region::Intersection;
use cdshealpix::nested::center;
use web_sys::WebGl2RenderingContext;
use cdshealpix::sph_geom::coo3d::{Coo3D};
use crate::renderable::line;
use al_core::{log, inforec, info};
@@ -52,7 +56,6 @@ use wasm_bindgen::JsValue;
use super::labels::RenderManager;
use super::TextRenderManager;
use super::line;
use al_api::resources::Resources;
impl ProjetedGrid {
@@ -184,7 +187,7 @@ impl ProjetedGrid {
fn force_update(&mut self, camera: &CameraViewPort, projection: &ProjectionType) {
self.text_renderer.begin_frame();
//let text_height = text_renderer.text_size();
let lines = lines(camera, &self.text_renderer, projection, &self.fmt);
let lines = lines(camera, projection, &self.fmt);
self.offsets.clear();
self.sizes.clear();
@@ -217,7 +220,7 @@ impl ProjetedGrid {
.flatten()
.flat_map(|v| [v.x as f32, v.y as f32])
.collect::<Vec<_>>();
//self.lines = lines;
self.num_vertices = vertices.len() >> 1;
#[cfg(feature = "webgl2")]
@@ -296,6 +299,7 @@ impl ProjetedGrid {
use crate::shader::ShaderId;
use core::num;
use std::borrow::Cow;
use std::path::is_separator;
@@ -309,34 +313,21 @@ use core::ops::Range;
#[derive(Debug)]
struct Label {
position: Vector2<f64>,
// The position
position: XYScreen,
// the string content
content: String,
// in radians
rot: f64,
}
impl Label {
fn meridian(
fov: &FieldOfView,
mut lon: f64,
m1: &Vector3<f64>,
fn from_meridian(
lonlat: &LonLatT<f64>,
camera: &CameraViewPort,
sp: Option<&Vector2<f64>>,
text_renderer: &TextRenderManager,
projection: &ProjectionType,
fmt: &angle::SerializeFmt
) -> Option<Self> {
let lat = camera.get_center().lonlat().lat();
// Do not plot meridian labels when the center of fov
// is above 80deg
if fov.is_allsky() {
// If allsky label plotting mode
// check if we are not too near of a pole
// If so, do not plot the meridian labels because
// they can overlap
if lat.abs() > ArcDeg(80.0) {
return None;
}
}
let fov = camera.get_field_of_view();
let d = if fov.contains_north_pole() {
Vector3::new(0.0, 1.0, 0.0)
} else if fov.contains_south_pole() {
@@ -345,23 +336,26 @@ impl Label {
Vector3::new(0.0, 1.0, 0.0)
};
let m2 = ((m1 + d * 1e-3).normalize()).extend(1.0);
let m1: Vector3<_> = lonlat.vector();
let m2 = (m1 + d * 1e-3).normalize();
//let s1 = projection.model_to_screen_space(&(system.to_icrs_j2000::<f64>() * m1), camera, reversed_longitude)?;
let s1 = projection.model_to_screen_space(&m1.extend(1.0), camera)?;
let s2 = projection.model_to_screen_space(&m2, camera)?;
let s2 = projection.model_to_screen_space(&m2.extend(1.0), camera)?;
//let s2 = projection.model_to_screen_space(&(system.to_icrs_j2000::<f64>() * m2), camera, reversed_longitude)?;
let ds = (s2 - s1).normalize();
let mut lon = m1.lon().to_radians();
if lon < 0.0 {
lon += TWICE_PI;
}
let content = fmt.to_string(Angle(lon));
let content = fmt.to_string(lon.to_angle());
let position = if !fov.is_allsky() {
//let dim = ctx2d.measure_text(&content).unwrap_abort();
let dim = text_renderer.get_width_pixel_size(&content);
let dim = 100.0;
let k = ds * (dim * 0.5 + 10.0);
s1 + k
@@ -369,8 +363,6 @@ impl Label {
s1
};
//position += dv * text_height * 0.5;
// rot is between -PI and +PI
let rot = if ds.y > 0.0 {
ds.x.acos()
@@ -396,44 +388,41 @@ impl Label {
})
}
fn parallel(
fov: &FieldOfView,
lat: f64,
m1: &Vector3<f64>,
fn from_parallel(
lonlat: &LonLatT<f64>,
camera: &CameraViewPort,
// in pixels
text_renderer: &TextRenderManager,
projection: &ProjectionType,
) -> Option<Self> {
let mut d = Vector3::new(-m1.z, 0.0, m1.x).normalize();
let _system = camera.get_system();
let center = camera.get_center().truncate();
//let center = (system.to_gal::<f64>() * camera.get_center()).truncate();
if center.dot(d) < 0.0 {
d = -d;
}
let m2 = (m1 + d * 1e-3).normalize();
let m1: Vector3<_> = lonlat.vector();
let s1 =
//projection.model_to_screen_space(&(system.to_icrs_j2000::<f64>() * m1.extend(1.0)), camera, reversed_longitude)?;
projection.model_to_screen_space(&m1.extend(1.0), camera)?;
let s2 =
//projection.model_to_screen_space(&(system.to_icrs_j2000::<f64>() * m2.extend(1.0)), camera, reversed_longitude)?;
projection.model_to_screen_space(&m2.extend(1.0), camera)?;
let mut t = Vector3::new(-m1.z, 0.0, m1.x).normalize();
let center = camera.get_center().truncate();
let dot_t_center = center.dot(t);
if dot_t_center.abs() < 1e-4 {
t = -t;
} else {
t = dot_t_center.signum() * t;
}
let m2 = (m1 + t * 1e-3).normalize();
let s1 = projection.model_to_screen_space(&m1.extend(1.0), camera)?;
let s2 = projection.model_to_screen_space(&m2.extend(1.0), camera)?;
let ds = (s2 - s1).normalize();
let content = angle::SerializeFmt::DMS.to_string(Angle(lat));
let content = angle::SerializeFmt::DMS.to_string(lonlat.lat());
let fov = camera.get_field_of_view();
let position = if !fov.is_allsky() && !fov.contains_pole() {
let dim = text_renderer.get_width_pixel_size(&content);
let dim = 100.0;
let k = ds * (dim * 0.5 + 10.0);
//let k = Vector2::new(0.0, 0.0);
s1 + k
} else {
s1
};
//position += dv * text_height * 0.5;
// rot is between -PI and +PI
let rot = if ds.y > 0.0 {
@@ -442,13 +431,9 @@ impl Label {
-ds.x.acos()
};
let rot = if ds.y > 0.0 {
if rot > HALF_PI {
-PI + rot
} else {
rot
}
} else if rot < -HALF_PI {
let rot = if ds.y > 0.0 && rot > HALF_PI {
-PI + rot
} else if ds.y <= 0.0 && rot < -HALF_PI {
PI + rot
} else {
rot
@@ -476,112 +461,162 @@ use crate::math::{
};
impl GridLine {
/*fn meridian(
fn meridian(
lon: f64,
lat: &Range<f64>,
sp: Option<&Vector2<f64>>,
camera: &CameraViewPort,
//text_height: f64,
text_renderer: &TextRenderManager,
projection: &ProjectionType,
fmt: &angle::SerializeFmt
) -> Option<Self> {
let fov = camera.get_field_of_view();
match fov {
FieldOfViewType::Polygon { poly, .. } => {
// Polygon case
let lon = if lon < 0.0 { lon + TWICE_PI } else { lon };
if fov.contains_both_poles() {
let camera_center = camera.get_center();
// The arc length must be < PI, so we create an arc from [(lon, -PI/2); (lon, PI/2)[
// see the cdshealpix doc:
// https://docs.rs/cdshealpix/latest/cdshealpix/sph_geom/struct.Polygon.html#method.intersect_great_circle_arc
let a = Coo3D::from_sph_coo(lon, -HALF_PI);
let b = Coo3D::from_sph_coo(lon, HALF_PI - 1e-6);
let center_lat = camera_center.lat().to_radians();
let label_lat = if center_lat > 70.0_f64.to_radians() {
70.0_f64.to_radians()
} else if center_lat < -70.0_f64.to_radians() {
(-70.0_f64).to_radians()
} else {
center_lat
};
// For those intersecting, perform the intersection
let inter_vertices = poly.intersect_great_circle_arc(&a, &b)
.iter()
.map(|v| Vector3::new(v.y(), v.z(), v.x()))
.collect::<Vec<_>>();
let lonlat = LonLatT::new(
lon.to_angle(),
label_lat.to_angle()
);
if fov.contains_both_poles() {
// Like the allsky case, we draw all the meridians
let n_pole_lonlat = LonLatT::new(Angle(lon), Angle(HALF_PI));
let s_pole_lonlat = LonLatT::new(Angle(lon), Angle(-HALF_PI));
let vertices = SmallCircle::new(&n_pole_lonlat, &s_pole_lonlat).project(camera, projection);
// Allsky case
// We do an approx saying allsky fovs intersect all meridian
// but this is not true for example for the orthographic projection
// Some meridians may not be visible
let center_lat = camera.get_center().lat();
let pos_label: Vector3<f64> = LonLatT::new(Angle(lon), center_lat).vector();
let label = Label::meridian(fov, lon, &pos_label, camera, sp, text_renderer, projection, fmt);
let label = Label::from_meridian(&lonlat, camera, projection, fmt);
// Draw the full parallel
let vertices = line::great_circle_arc::project(lon, -HALF_PI, lon, HALF_PI, camera, projection);
Some(GridLine { vertices, label })
} else {
let i = fov.intersects_meridian(lon);
match i {
Intersection::Included => {
// Longitude fov >= PI
let camera_center = camera.get_center();
let center_lat = camera_center.lat().to_radians();
let lonlat = LonLatT::new(
lon.to_angle(),
center_lat.to_angle()
);
let label = Label::from_meridian(&lonlat, camera, projection, fmt);
// Draw the full parallel
let vertices = line::great_circle_arc::project(lon, -HALF_PI, lon, HALF_PI, camera, projection);
Some(GridLine { vertices, label })
} else if fov.contains_north_pole() {
// Only one intersection
let n_pole_lonlat = LonLatT::new(Angle(lon), Angle(HALF_PI));
let i_lonlat = inter_vertices[0].lonlat();
},
Intersection::Intersect { vertices } => {
let num_intersections = vertices.len();
let (vertices, label) = match num_intersections {
1 => {
let v1 = &vertices[0];
let lonlat1 = v1.lonlat();
let lat1 = lonlat1.lat().to_radians();
let vertices = SmallCircle::new(&n_pole_lonlat, &i_lonlat).project(camera, projection);
let line_vertices = if fov.contains_north_pole() {
line::great_circle_arc::project(
lon,
lat1,
lon,
HALF_PI,
camera,
projection
)
} else {
line::great_circle_arc::project(
lon,
lat1,
lon,
MINUS_HALF_PI,
camera,
projection
)
};
let pos_label = &inter_vertices[0];
let label = Label::meridian(fov, lon, pos_label, camera, sp, text_renderer, projection, fmt);
let label = Label::from_meridian(&lonlat1, camera, projection, fmt);
(line_vertices, label)
},
2 => {
// full intersection
let v1 = &vertices[0];
let v2 = &vertices[1];
let lat1 = v1.lat().to_radians();
let lat2 = v2.lat().to_radians();
let line_vertices = line::great_circle_arc::project(
lon,
lat1,
lon,
lat2,
camera,
projection
);
let label = Label::from_meridian(&v1.lonlat(), camera, projection, fmt);
(line_vertices, label)
},
_ => {
let mut vertices = vertices.into_vec();
// One segment over two will be in the field of view
vertices.push(Vector4::new(0.0, 1.0, 0.0, 1.0));
vertices.push(Vector4::new(0.0, -1.0, 0.0, 1.0));
vertices.sort_by(|i1, i2| {
i1.y.total_cmp(&i2.y)
});
let v1 = &vertices[0];
let v2 = &vertices[1];
// meridian are part of great circles so the mean between v1 & v2 also lies on it
let vm = (v1 + v2).truncate().normalize();
let vertices = if !fov.contains_south_pole() {
&vertices[1..]
} else {
&vertices
};
let line_vertices = vertices.iter().zip(vertices.iter().skip(1))
.step_by(2)
.map(|(i1, i2)| {
line::great_circle_arc::project(
lon,
i1.lat().to_radians(),
lon,
i2.lat().to_radians(),
camera,
projection
)
})
.flatten()
.collect::<Vec<_>>();
let label = Label::from_meridian(&v1.lonlat(), camera, projection, fmt);
(line_vertices, label)
}
};
Some(GridLine { vertices, label })
} else if fov.contains_south_pole() {
// Only one intersection
let i_lonlat = inter_vertices[0].lonlat();
let s_pole_lonlat = LonLatT::new(Angle(lon), Angle(-HALF_PI));
let vertices = SmallCircle::new(&i_lonlat, &s_pole_lonlat).project(camera, projection);
let pos_label = &inter_vertices[0];
let label = Label::meridian(fov, lon, pos_label, camera, sp, text_renderer, projection, fmt);
Some(GridLine { vertices, label })
} else {
if inter_vertices.len() >= 2 {
// There must be at least 2 points intersecting the fov
let i1_lonlat = inter_vertices[0].lonlat();
let i2_lonlat = inter_vertices[1].lonlat();
let vertices = SmallCircle::new(&i1_lonlat, &i2_lonlat).project(camera, projection);
let pos_label = &inter_vertices[0];
let label = Label::meridian(fov, lon, pos_label, camera, sp, text_renderer, projection, fmt);
Some(GridLine { vertices, label })
} else {
None
}
}
},
FieldOfViewType::Allsky => {
let n_pole_lonlat = LonLatT::new(Angle(lon), Angle(HALF_PI));
let s_pole_lonlat = LonLatT::new(Angle(lon), Angle(-HALF_PI));
let vertices = SmallCircle::new(&s_pole_lonlat, &n_pole_lonlat).project(camera, projection);
// Allsky case
// We do an approx saying allsky fovs intersect all meridian
// but this is not true for example for the orthographic projection
// Some meridians may not be visible
let center_lat = camera.get_center().lat();
let pos_label: Vector3<f64> = LonLatT::new(Angle(lon), center_lat).vector();
let label = Label::meridian(fov, lon, &pos_label, camera, sp, text_renderer, projection, fmt);
Some(GridLine { vertices, label })
},
Intersection::Empty => {
None
},
}
}
}*/
}
fn parallel(
lat: f64,
camera: &CameraViewPort,
text_renderer: &TextRenderManager,
projection: &ProjectionType,
) -> Option<Self> {
let fov = camera.get_field_of_view();
@@ -589,15 +624,16 @@ impl GridLine {
// Longitude fov >= PI
let camera_center = camera.get_center();
let center_lon = camera_center.lon();
let label_vertex = LonLatT::new(
let lonlat = LonLatT::new(
center_lon,
lat.to_angle()
).vector();
);
let label = Label::parallel(fov, lat, &label_vertex, camera, text_renderer, projection);
let label = Label::from_parallel(&lonlat, camera, projection);
// Draw the full parallel
let vertices = crate::renderable::line::parallel::project(lat, center_lon.to_radians(), center_lon.to_radians() + TWICE_PI, camera, projection);
let mut vertices = line::parallel_arc::project(lat, center_lon.to_radians(), center_lon.to_radians() + PI, camera, projection);
vertices.append(&mut line::parallel_arc::project(lat, center_lon.to_radians() + PI, center_lon.to_radians() + TWICE_PI, camera, projection));
Some(GridLine { vertices, label })
} else {
// Longitude fov < PI
@@ -606,46 +642,35 @@ impl GridLine {
Intersection::Included => {
let camera_center = camera.get_center();
let center_lon = camera_center.lon();
let label_vertex = LonLatT::new(
let lonlat = LonLatT::new(
center_lon,
lat.to_angle()
).vector();
);
let label = Label::parallel(fov, lat, &label_vertex, camera, text_renderer, projection);
let label = Label::from_parallel(&lonlat, camera, projection);
// Draw the full parallel
let vertices = crate::renderable::line::parallel::project(lat, center_lon.to_radians(), center_lon.to_radians() + TWICE_PI, camera, projection);
let vertices = line::parallel_arc::project(lat, center_lon.to_radians(), center_lon.to_radians() + TWICE_PI, camera, projection);
Some(GridLine { vertices, label })
},
Intersection::Intersect { vertices: inter_vertices } => {
let is_intersecting_zero_meridian = fov.is_intersecting_zero_meridian();
Intersection::Intersect { vertices } => {
let v1 = &vertices[0];
let v2 = &vertices[1];
let lonlat1 = inter_vertices[0].lonlat();
let lonlat2 = inter_vertices[1].lonlat();
let mut lon1 = v1.lon().to_radians();
let mut lon2 = v2.lon().to_radians();
let center_lonlat = camera.get_center().lonlat();
let lon0 = center_lonlat.lon().to_radians();
let mut lon1 = lonlat1.lon().to_radians();
let mut lon2 = lonlat2.lon().to_radians();
let lon_len = if lon1 > lon2 {
lon2 + TWICE_PI - lon1
} else {
lon2 - lon1
};
let lon_len = crate::math::sph_geom::distance_from_two_lon(lon1, lon2);
// The fov should be contained into PI length
if lon_len >= PI {
std::mem::swap(&mut lon1, &mut lon2);
}
let lat_deg = lat.to_degrees();
al_core::info!(lat_deg, ":", lon0, lon1, lon2);
let vertices = crate::renderable::line::parallel::project(lat, lon1, lon2, camera, projection);
let label = Label::parallel(fov, lat, &inter_vertices[0].truncate(), camera, text_renderer, projection);
let line_vertices = line::parallel_arc::project(lat, lon1, lon2, camera, projection);
let label = Label::from_parallel(&v1.lonlat(), camera, projection);
Some(GridLine { vertices, label })
Some(GridLine { vertices: line_vertices, label })
},
Intersection::Empty => {
None
@@ -656,52 +681,51 @@ impl GridLine {
}
const GRID_STEPS: &[f64] = &[
0.0000000000048481366,
0.000000000009696273,
0.000000000024240684,
0.000000000048481368,
0.000000000096962736,
0.00000000024240682,
0.00000000048481363,
0.0000000009696273,
0.0000000024240685,
0.000000004848137,
0.000000009696274,
0.000000024240684,
0.00000004848137,
0.00000009696274,
0.00000024240686,
0.0000004848137,
0.0000009696274,
0.0000024240685,
0.000004848137,
0.000009696274,
0.000024240684,
0.000048481368,
0.000072722054,
0.00014544411,
0.00029088822,
0.00058177643,
0.0014544411,
0.0029088822,
0.004363323,
0.008726646,
0.017453292,
0.034906585,
0.08726646,
0.17453292,
0.34906584,
0.0000000000048481367,
0.000000000009696274,
0.000000000024240685,
0.000000000048481369,
0.000000000096962737,
0.00000000024240683,
0.00000000048481364,
0.0000000009696274,
0.0000000024240686,
0.000000004848138,
0.000000009696275,
0.000000024240685,
0.00000004848138,
0.00000009696275,
0.00000024240687,
0.0000004848138,
0.0000009696275,
0.0000024240686,
0.000004848138,
0.000009696275,
0.000024240685,
0.000048481369,
0.000072722055,
0.00014544412,
0.00029088823,
0.00058177644,
0.0014544412,
0.0029088823,
0.004363324,
0.008726647,
0.017453293,
0.034906586,
0.08726647,
0.17453293,
0.34906585,
std::f64::consts::FRAC_PI_4,
];
fn lines(
camera: &CameraViewPort,
//text_height: f64,
text_renderer: &TextRenderManager,
projection: &ProjectionType,
fmt: &angle::SerializeFmt,
) -> Vec<GridLine> {
// Get the screen position of the nearest pole
let _system = camera.get_system();
let fov = camera.get_field_of_view();
let sp = if fov.contains_pole() {
if fov.contains_north_pole() {
@@ -746,28 +770,23 @@ fn lines(
stop_theta -= 1e-3;
}
/*while theta < stop_theta {
while theta < stop_theta {
if let Some(line) =
GridLine::meridian(theta, &bbox.get_lat(), sp.as_ref(), camera, text_renderer, projection, fmt)
GridLine::meridian(theta, &bbox.get_lat(), sp.as_ref(), camera, projection, fmt)
{
lines.push(line);
}
theta += step_lon;
}*/
// Add parallels
//let step_lat = select_grid_step(bbox, bbox.get_lat_size() as f64, NUM_LINES);
}
let mut alpha = bbox.lat_min() - (bbox.lat_min() % step_lat);
if alpha == -HALF_PI {
alpha += step_lat;
}
let stop_alpha = bbox.lat_max();
/*if stop_alpha == HALF_PI {
stop_alpha -= 1e-3;
}*/
while alpha < stop_alpha {
if let Some(line) = GridLine::parallel(alpha, camera, text_renderer, projection) {
if let Some(line) = GridLine::parallel(alpha, camera, projection) {
lines.push(line);
}
alpha += step_lat;

317
src/core/src/renderable/line/great_circle_arc.rs
View File

@@ -1,140 +1,217 @@
use crate::{math::{lonlat::LonLatT, projection::{ProjectionType, coo_space::XYNDC}}, camera::CameraViewPort};
use crate::math::angle::Angle;
use crate::math::projection::coo_space::XYZWModel;
use cgmath::Vector2;
use cgmath::Vector3;
use crate::ProjectionType;
use crate::CameraViewPort;
use cgmath::Zero;
use cgmath::InnerSpace;
use crate::math::angle::ToAngle;
use crate::math::lonlat::LonLat;
use al_core::{log, info, inforec};
use crate::coo_space::XYZWModel;
use crate::coo_space::XYNDC;
use crate::coo_space::XYZModel;
use crate::math::{TWICE_PI, PI};
use crate::ArcDeg;
const MAX_ANGLE_BEFORE_SUBDIVISION: Angle<f64> = Angle(0.10943951023); // 12 degrees
const MAX_ITERATION: usize = 3;
use crate::LonLatT;
const MAX_ANGLE_BEFORE_SUBDIVISION: f64 = 0.1;
const MAX_ITERATION: usize = 5;
pub fn project(lonlat1: &LonLatT<f64>, lonlat2: &LonLatT<f64>, camera: &CameraViewPort, projection: &ProjectionType) -> Vec<XYNDC> {
// First longitude between [0; 2*pi[
let mut lon1 = lonlat1.lon().to_radians();
// Parallel latitude between [-0.5*pi; 0.5*pi]
let lat1 = lonlat1.lat().to_radians();
// Second longitude between [0; 2*pi[
let mut lon2 = lonlat2.lon().to_radians();
// Parallel latitude between [-0.5*pi; 0.5*pi]
let lat2 = lonlat2.lat().to_radians();
let is_intersecting_zero_meridian = lon1 > lon2;
// Requirement:
// * Latitudes between [-0.5*pi; 0.5*pi]
// * Longitudes between [0; 2\pi[
// * (lon1 - lon2).abs() < PI so that is can only either cross the preimary meridian or opposite primary meridian
// (the latest is handled because of the longitudes intervals)
pub fn project(lon1: f64, lat1: f64, lon2: f64, lat2: f64, camera: &CameraViewPort, projection: &ProjectionType) -> Vec<XYNDC> {
let mut vertices = vec![];
if is_intersecting_zero_meridian {
// Make the longitudes lie between [-PI; PI];
if lon1 > PI {
lon1 -= TWICE_PI;
}
let lonlat1 = LonLatT::new(lon1.to_angle(), lat1.to_angle());
let lonlat2 = LonLatT::new(lon2.to_angle(), lat2.to_angle());
if lon2 > PI {
lon2 -= TWICE_PI;
let v1: Vector3<_> = lonlat1.vector();
let v2: Vector3<_> = lonlat2.vector();
let p1 = projection.model_to_normalized_device_space(&v1.extend(1.0), camera);
let p2 = projection.model_to_normalized_device_space(&v2.extend(1.0), camera);
match (p1, p2) {
(Some(_), Some(_)) => {
project_line(&mut vertices, &v1, &v2, camera, projection, 0);
},
(None, Some(_)) => {
let (v1, v2) = sub_valid_domain(v2, v1, projection, camera);
project_line(&mut vertices, &v1, &v2, camera, projection, 0);
},
(Some(_), None) => {
let (v1, v2) = sub_valid_domain(v1, v2, projection, camera);
project_line(&mut vertices, &v1, &v2, camera, projection, 0);
},
(None, None) => {
}
}
let mut ndc_vertices: Vec<XYNDC> = vec![];
let start_world_vertex = LonLatT::new(lon1.to_angle(), lat1.to_angle()).vector();
let end_world_vertex = LonLatT::new(lon2.to_angle(), lat2.to_angle()).vector();
let ndc_v1 = projection.model_to_normalized_device_space(&start_world_vertex, camera);
let ndc_v2 = projection.model_to_normalized_device_space(&end_world_vertex, camera);
if let (Some(start_ndc_vertex), Some(end_ndc_vertex)) = (ndc_v1, ndc_v2) {
subdivide(
&mut ndc_vertices,
start_world_vertex,
start_ndc_vertex,
end_world_vertex,
end_ndc_vertex,
camera,
projection,
0
);
}
ndc_vertices
vertices
}
fn subdivide(
ndc_vertices: &mut Vec<XYNDC>,
start_world_vertex: XYZWModel,
start_ndc_vertex: XYNDC,
// Precondition:
// * angular distance between valid_lon and invalid_lon is < PI
// * valid_lon and invalid_lon are well defined, i.e. they can be between [-PI; PI] or [0, 2PI] depending
// whether they cross or not the zero meridian
fn sub_valid_domain(valid_v: XYZModel, invalid_v: XYZModel, projection: &ProjectionType, camera: &CameraViewPort) -> (XYZModel, XYZModel) {
let d_alpha = camera.get_aperture().to_radians() * 0.02;
end_world_vertex: XYZWModel,
end_ndc_vertex: XYNDC,
camera: &CameraViewPort,
projection: &ProjectionType,
iter: usize,
) {
if iter > MAX_ITERATION {
ndc_vertices.push(start_ndc_vertex);
ndc_vertices.push(end_ndc_vertex);
return;
}
// Project them. We are always facing the camera
let mid_world_vertex = (start_world_vertex + end_world_vertex).normalize();
if let Some(mid_ndc_vertex) = projection.model_to_normalized_device_space(&mid_world_vertex, camera) {
let ab = mid_ndc_vertex - start_ndc_vertex;
let bc = end_ndc_vertex - mid_ndc_vertex;
let ab_l = ab.magnitude2();
let bc_l = bc.magnitude2();
let ab = ab.normalize();
let bc = bc.normalize();
let theta = crate::math::vector::angle2(&ab, &bc);
// nearly colinear vectors
if theta.abs() < MAX_ANGLE_BEFORE_SUBDIVISION {
if crate::math::vector::det(&ab, &bc).abs() < 1e-2 {
ndc_vertices.push(start_ndc_vertex);
ndc_vertices.push(end_ndc_vertex);
} else {
// not colinear
ndc_vertices.push(start_ndc_vertex);
ndc_vertices.push(mid_ndc_vertex);
ndc_vertices.push(mid_ndc_vertex);
ndc_vertices.push(end_ndc_vertex);
}
} else if ab_l.min(bc_l) / ab_l.max(bc_l) < 0.1 {
if ab_l < bc_l {
ndc_vertices.push(start_ndc_vertex);
ndc_vertices.push(mid_ndc_vertex);
} else {
ndc_vertices.push(mid_ndc_vertex);
ndc_vertices.push(end_ndc_vertex);
}
let mut vv = valid_v;
let mut vi = invalid_v;
while crate::math::vector::angle3(&vv, &vi).to_radians() > d_alpha {
let vm = (vv + vi).normalize();
// check whether is it defined or not
if let Some(_) = projection.model_to_normalized_device_space(&vm.extend(1.0), camera) {
vv = vm;
} else {
// Subdivide a->b and b->c
subdivide(
ndc_vertices,
start_world_vertex,
start_ndc_vertex,
mid_world_vertex,
mid_ndc_vertex,
camera,
projection,
iter + 1
);
subdivide(
ndc_vertices,
mid_world_vertex,
mid_ndc_vertex,
end_world_vertex,
end_ndc_vertex,
camera,
projection,
iter + 1
);
vi = vm;
}
}
}
// Return the valid interval found by dichotomy
(vv, valid_v)
}
fn project_line(
vertices: &mut Vec<XYNDC>,
v1: &XYZModel,
v2: &XYZModel,
camera: &CameraViewPort,
projection: &ProjectionType,
iter: usize,
) -> bool {
let p1 = projection.model_to_normalized_device_space(&v1.extend(1.0), camera);
let p2 = projection.model_to_normalized_device_space(&v2.extend(1.0), camera);
if iter < MAX_ITERATION {
// Project them. We are always facing the camera
let vm = (v1 + v2).normalize();
let pm = projection.model_to_normalized_device_space(&vm.extend(1.0), camera);
match (p1, pm, p2) {
(Some(p1), Some(pm), Some(p2)) => {
let d12 = crate::math::vector::angle3(v1, v2).to_radians();
// Subdivide when until it is > 30 degrees
if d12 > 30.0_f64.to_radians() {
subdivide(
vertices,
v1,
v2,
&vm,
p1,
p2,
pm,
camera,
projection,
iter
);
} else {
// enough to stop the recursion
let ab = pm - p1;
let bc = p2 - pm;
let ab_u = ab.normalize();
let bc_u = bc.normalize();
let dot_abbc = crate::math::vector::dot(&ab_u, &bc_u);
let theta_abbc = dot_abbc.acos();
if theta_abbc.abs() < 5.0_f64.to_radians() {
let det_abbc = crate::math::vector::det(&ab_u, &bc_u);
if det_abbc.abs() < 1e-2 {
vertices.push(p1);
vertices.push(p2);
} else {
// not colinear but enough to stop
vertices.push(p1);
vertices.push(pm);
vertices.push(pm);
vertices.push(p2);
}
} else {
let ab_l = ab.magnitude2();
let bc_l = bc.magnitude2();
let r = (ab_l - bc_l).abs() / (ab_l + bc_l);
if r > 0.8 {
if ab_l < bc_l {
vertices.push(p1);
vertices.push(pm);
} else {
vertices.push(pm);
vertices.push(p2);
}
} else {
// Subdivide a->b and b->c
subdivide(
vertices,
v1,
v2,
&vm,
p1,
p2,
pm,
camera,
projection,
iter
);
}
}
}
true
},
_ => false
}
} else {
false
}
}
fn subdivide(
vertices: &mut Vec<XYNDC>,
v1: &XYZModel,
v2: &XYZModel,
vm: &XYZModel,
p1: XYNDC,
p2: XYNDC,
pm: XYNDC,
camera: &CameraViewPort,
projection: &ProjectionType,
iter: usize
) {
// Subdivide a->b and b->c
if !project_line(
vertices,
v1,
vm,
camera,
projection,
iter + 1
) {
vertices.push(p1);
vertices.push(pm);
}
if !project_line(
vertices,
vm,
v2,
camera,
projection,
iter + 1
) {
vertices.push(pm);
vertices.push(p2);
}
}

0
src/core/src/renderable/line/meridian.rs
View File

3
src/core/src/renderable/line/mod.rs
View File

@@ -1,7 +1,6 @@
/// This module handles the lines rendering code
pub mod meridian;
pub mod parallel;
pub mod great_circle_arc;
pub mod parallel_arc;
use al_core::WebGlContext;
use al_core::VertexArrayObject;

155
src/core/src/renderable/line/parallel.rs → src/core/src/renderable/line/parallel_arc.rs
View File

@@ -7,117 +7,52 @@ use cgmath::Zero;
use cgmath::InnerSpace;
use crate::math::angle::ToAngle;
use crate::math::lonlat::LonLat;
use al_core::{log, info, inforec};
use crate::coo_space::XYNDC;
use crate::math::{TWICE_PI, PI};
use crate::ArcDeg;
use crate::LonLatT;
const MAX_ANGLE_BEFORE_SUBDIVISION: Angle<f64> = Angle(0.174533); // 12 degrees
const MAX_ANGLE_BEFORE_SUBDIVISION: Angle<f64> = Angle(0.10);
const MAX_ITERATION: usize = 4;
// Requirement:
// * Parallel latitude between [-0.5*pi; 0.5*pi]
// * First longitude between [0; 2\pi[
// * Second lon length between [0; 2\pi[
// * Either (lon1, lat) or (lon2, lat) is defined
pub fn project(lat: f64, mut lon1: f64, lon2: f64, camera: &CameraViewPort, projection: &ProjectionType) -> Vec<XYNDC> {
let lon_len = if lon1 > lon2 {
lon2 + TWICE_PI - lon1
} else {
lon2 - lon1
};
// * (lon1 - lon2).abs() < PI
pub fn project(lat: f64, mut lon1: f64, mut lon2: f64, camera: &CameraViewPort, projection: &ProjectionType) -> Vec<XYNDC> {
let mut vertices = vec![];
if lon_len > PI + 1e-6 {
// lon mid lies in [0; 2\pi[
let lon_mid = if lon1 + PI >= TWICE_PI {
lon1 -= TWICE_PI;
lon1 + PI
} else {
lon1 + PI
};
let lon_len = crate::math::sph_geom::distance_from_two_lon(lon1, lon2);
let mut lon2 = lon1 + lon_len;
let first_half_meridian = project(lat, lon1, lon_mid, camera, projection);
let mut second_half_meridian = project(lat, lon_mid, lon2, camera, projection);
let mut vertices = first_half_meridian;
vertices.append(&mut second_half_meridian);
vertices
} else {
let mut lon2 = lon1 + lon_len;
let mut vertices = vec![];
// Can only cross the 0 meridian but not 0 and 180 ones
if lon2 > TWICE_PI {
// it crosses the zero meridian
lon2 -= TWICE_PI;
// lon1 is > PI because the lon len is <= PI
lon1 -= TWICE_PI;
}
// We know (lon1, lat) can be projected as it is a requirement of that method
let v1 = crate::math::lonlat::proj(&LonLatT::new(lon1.to_angle(), lat.to_angle()), projection, camera);
let v2 = crate::math::lonlat::proj(&LonLatT::new(lon2.to_angle(), lat.to_angle()), projection, camera);
match (v1, v2) {
(Some(v1), Some(v2)) => {
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(None, Some(v2)) => {
let (lon1, lon2) = sub_valid_domain(lat, lon2, lon1, projection, camera);
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(Some(v1), None) => {
let (lon1, lon2) = sub_valid_domain(lat, lon1, lon2, projection, camera);
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(None, None) => {}
}
vertices
// Can only cross the 0 meridian but not 0 and 180 ones
if lon2 > TWICE_PI {
// it crosses the zero meridian
lon2 -= TWICE_PI;
// lon1 is > PI because the lon len is <= PI
lon1 -= TWICE_PI;
}
}
#[inline]
pub fn is_in_lon_range(lon0: f64, lon1: f64, lon2: f64) -> bool {
// First version of the code:
// ((v2.lon() - v1.lon()).abs() > PI) != ((v2.lon() > coo.lon()) != (v1.lon() > coo.lon()))
//
// Lets note
// - lonA = v1.lon()
// - lonB = v2.lon()
// - lon0 = coo.lon()
// When (lonB - lonA).abs() <= PI
// => lonB > lon0 != lonA > lon0 like in PNPOLY
// A B lonA <= lon0 && lon0 < lonB
// --[++++[--
// B A lonB <= lon0 && lon0 < lonA
//
// But when (lonB - lonA).abs() > PI, then the test should be
// => lonA >= lon0 == lonB >= lon0
// <=> !(lonA >= lon0 != lonB >= lon0)
// A | B (lon0 < lonB) || (lonA <= lon0)
// --[++|++[--
// B | A (lon0 < lonA) || (lonB <= lon0)
//
// Instead of lonA > lon0 == lonB > lon0,
// i.e. !(lonA > lon0 != lonB > lon0).
// A | B (lon0 <= lonB) || (lonA < lon0)
// --]++|++]--
// B | A (lon0 <= lonA) || (lonB < lon0)
//
// So the previous code was bugged in this very specific case:
// - `lon0` has the same value as a vertex being part of:
// - one segment that do not cross RA=0
// - plus one segment crossing RA=0.
// - the point have an odd number of intersections with the polygon
// (since it will be counted 0 or 2 times instead of 1).
let dlon = lon2 - lon1;
if dlon < 0.0 {
(dlon >= -PI) == (lon2 <= lon0 && lon0 < lon1)
} else {
(dlon <= PI) == (lon1 <= lon0 && lon0 < lon2)
// We know (lon1, lat) can be projected as it is a requirement of that method
let v1 = crate::math::lonlat::proj(&LonLatT::new(lon1.to_angle(), lat.to_angle()), projection, camera);
let v2 = crate::math::lonlat::proj(&LonLatT::new(lon2.to_angle(), lat.to_angle()), projection, camera);
match (v1, v2) {
(Some(v1), Some(v2)) => {
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(None, Some(v2)) => {
let (lon1, lon2) = sub_valid_domain(lat, lon2, lon1, projection, camera);
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(Some(v1), None) => {
let (lon1, lon2) = sub_valid_domain(lat, lon1, lon2, projection, camera);
subdivide_multi(&mut vertices, lat, lon1, lon2, camera, projection);
},
(None, None) => {}
}
vertices
}
// Precondition:
@@ -179,11 +114,11 @@ fn subdivide(
camera: &CameraViewPort,
projection: &ProjectionType,
iter: usize,
) {
if iter < MAX_ITERATION {
let p1 = crate::math::lonlat::proj(&LonLatT::new(lon1.to_angle(), lat.to_angle()), projection, camera);
let p2 = crate::math::lonlat::proj(&LonLatT::new(lon2.to_angle(), lat.to_angle()), projection, camera);
) -> bool {
let p1 = crate::math::lonlat::proj(&LonLatT::new(lon1.to_angle(), lat.to_angle()), projection, camera);
let p2 = crate::math::lonlat::proj(&LonLatT::new(lon2.to_angle(), lat.to_angle()), projection, camera);
if iter < MAX_ITERATION {
// Project them. We are always facing the camera
let lon0 = (lon1 + lon2)*0.5;
let pm = crate::math::lonlat::proj(&LonLatT::new(lon0.to_angle(), lat.to_angle()), projection, camera);
@@ -223,7 +158,7 @@ fn subdivide(
}
} else {
// Subdivide a->b and b->c
subdivide(
if !subdivide(
vertices,
lat,
lon1,
@@ -231,9 +166,12 @@ fn subdivide(
camera,
projection,
iter + 1
);
) {
vertices.push(p1);
vertices.push(pm);
}
subdivide(
if !subdivide(
vertices,
lat,
lon0,
@@ -241,10 +179,17 @@ fn subdivide(
camera,
projection,
iter + 1
);
) {
vertices.push(pm);
vertices.push(p2);
}
}
true
},
_ => {}
_ => false
}
} else {
false
}
}