From 69e508d8c94d8973033817ca86f357e466fc7c4d Mon Sep 17 00:00:00 2001 From: Joshua Yanovski Date: Tue, 7 Jul 2020 18:41:37 +0200 Subject: [PATCH] Make it easy to switch to SIMD for math. --- voxygen/Cargo.toml | 5 +- voxygen/src/scene/math.rs | 626 +++++++++++++++++++++++++++++++ voxygen/src/scene/mod.rs | 689 ++--------------------------------- voxygen/src/scene/terrain.rs | 50 +-- 4 files changed, 697 insertions(+), 673 deletions(-) create mode 100644 voxygen/src/scene/math.rs diff --git a/voxygen/Cargo.toml b/voxygen/Cargo.toml index 7559fad914..0009d8c47d 100644 --- a/voxygen/Cargo.toml +++ b/voxygen/Cargo.toml @@ -38,7 +38,10 @@ specs = "0.15.1" specs-idvs = { git = "https://gitlab.com/veloren/specs-idvs.git" } # Mathematics -vek = { version = "0.11.0", features = ["serde"] } +vek = { version = "0.11.0", features = [ + # "repr_simd", + "serde" +] } # Controller gilrs = { version = "0.7", features = ["serde"] } diff --git a/voxygen/src/scene/math.rs b/voxygen/src/scene/math.rs new file mode 100644 index 0000000000..cd47a51012 --- /dev/null +++ b/voxygen/src/scene/math.rs @@ -0,0 +1,626 @@ +use core::{iter, mem}; +use hashbrown::HashMap; +use num::traits::Float; +/* pub use vek::{ + geom::repr_simd::*, mat::repr_simd::column_major::Mat4, ops::*, vec::repr_simd::*, +}; */ +pub use vek::{geom::repr_c::*, mat::repr_c::column_major::Mat4, ops::*, vec::repr_c::*}; + +pub fn aabb_to_points(bounds: Aabb) -> [Vec3; 8] { + [ + Vec3::new(bounds.min.x, bounds.min.y, bounds.min.z), + Vec3::new(bounds.max.x, bounds.min.y, bounds.min.z), + Vec3::new(bounds.max.x, bounds.max.y, bounds.min.z), + Vec3::new(bounds.min.x, bounds.max.y, bounds.min.z), + Vec3::new(bounds.min.x, bounds.min.y, bounds.max.z), + Vec3::new(bounds.max.x, bounds.min.y, bounds.max.z), + Vec3::new(bounds.max.x, bounds.max.y, bounds.max.z), + Vec3::new(bounds.min.x, bounds.max.y, bounds.max.z), + ] +} + +/// Each Vec4 should be interpreted as reprenting plane +/// equation +/// +/// a(x - x0) + b(y - y0) + c(z - z0) = 0, i.e. +/// ax + by + cz - (a * x0 + b * y0 + c * z0) = 0, i.e. +/// ax + by + cz = (a * x0 + b * y0 + c * z0), i.e. +/// (lettiing d = a * x0 + b * y0 + c * z0) +/// ax + by + cz = d +/// +/// where d is the distance of the plane from the origin. +pub fn aabb_to_planes(bounds: Aabb) -> [(Vec3, T); 6] { + let zero = T::zero(); + let one = T::one(); + let bounds = bounds.map(|e| e.abs()); + [ + // bottom + (Vec3::new(zero, -one, zero), bounds.min.y), + // top + (Vec3::new(zero, one, zero), bounds.max.y), + // left + (Vec3::new(-one, zero, zero), bounds.min.x), + // right + (Vec3::new(one, zero, zero), bounds.max.x), + // near + (Vec3::new(zero, zero, -one), bounds.min.z), + // far + (Vec3::new(zero, zero, one), bounds.max.z), + ] +} + +pub fn mat_mul_points>( + mat: Mat4, + pts: &mut [Vec3], + mut do_p: impl FnMut(Vec4) -> Vec3, +) { + pts.into_iter().for_each(|p| { + *p = do_p(mat * Vec4::from_point(*p)); + }); +} + +/// NOTE: Expects points computed from aabb_to_points. +pub fn calc_view_frust_object(pts: &[Vec3; 8]) -> Vec>> { + vec![ + // near (CCW) + vec![pts[0], pts[1], pts[2], pts[3]], + // far (CCW) + vec![pts[7], pts[6], pts[5], pts[4]], + // left (CCW) + vec![pts[0], pts[3], pts[7], pts[4]], + // right (CCW) + vec![pts[1], pts[5], pts[6], pts[2]], + // bottom (CCW) + vec![pts[4], pts[5], pts[1], pts[0]], + // top (CCW) + vec![pts[6], pts[7], pts[3], pts[2]], + ] +} + +pub fn calc_view_frustum_world_coord>( + inv_proj_view: Mat4, +) -> [Vec3; 8] { + let mut world_pts = aabb_to_points(Aabb { + min: -Vec3::one(), + max: Vec3::one(), + }); + mat_mul_points(inv_proj_view, &mut world_pts, |p| Vec3::from(p) / p.w); + world_pts +} + +pub fn point_plane_distance(point: Vec3, (norm, dist): (Vec3, T)) -> T { + norm.dot(point) - dist +} + +pub fn point_before_plane(point: Vec3, plane: (Vec3, T)) -> bool { + point_plane_distance(point, plane) > T::zero() +} + +/// Returns true if and only if the final point in the polygon (i.e. the +/// first point added to the new polygon) is outside the clipping plane +/// (this implies that the polygon must be non-degenerate). +pub fn clip_points_by_plane + core::fmt::Debug>( + points: &mut Vec>, + plane: (Vec3, T), + intersection_points: &mut Vec>, +) -> bool { + /* enum Intersection { + /// Previous point was inside the plane. + Inside, + /// Previous line segment was completely outside the plane. + Outside, + /// Previous line segment went from inside the plane to outside it. + InsideOut, + } */ + // println!("points@clip_points_by_plane before clipping by {:?}: {:?}", plane, + // points); + if points.len() < 3 { + return false; + } + // NOTE: Guaranteed to succeed since points.len() > 3. + let mut current_point = points[points.len() - 1]; + let (norm, dist) = plane; + let intersect_plane_edge = |a, b| { + let diff = b - a; + let t = norm.dot(diff); + if t == T::zero() { + None + } else { + let t = (dist - norm.dot(a)) / t; + if t < T::zero() || T::one() < t { + None + } else { + Some(diff.mul_add(Vec3::broadcast(t), a)) + } + } + }; + let last_is_outside = point_before_plane(current_point, plane); + let mut is_outside = last_is_outside; + /* // Might not actually be total, but if it is partial and the point is inside it will be + // written regardless, and if it is partial and the point is outside, it means the + // second-to-last point is inside; thus, the second-to-last point will be written regardless, + // current_point will hold the new intersection point, and is_total will be false, when the + // loop ends; thus all we need to do to take this case into account is to push current_point + // onto the points vector if (is_total || is_outside) is false at the end of the loop. + let mut is_total = true; */ + let mut old_points = Vec::with_capacity((3 * points.len()) / 2); + mem::swap(&mut old_points, points); + old_points.into_iter().for_each(|point| { + /* let prev_point = current_point; + // Swap point i with the previous point in the polygon, so it is the one we normally save + // when we return false. + mem::swap(&mut current_point, point); */ + let prev_point = mem::replace(&mut current_point, point); + /* if point_before_plane(current_point) { + // If we are an outside point, we should only calculate an intersection if the previous + // point was inside. + if + is_outside s + // point was outside. + } else { + // If we are an inside point, then we should only calculate an intersection if the previous + // point was outside. + } */ + let before_plane = point_before_plane(current_point, plane); + let prev_is_outside = mem::replace(&mut is_outside, before_plane); + // println!("points@clip_points_by_plane clipping segment by {:?} (prev={:?} / + // outside={:?}, current={:?} / outside={:?})", plane, prev_point, + // prev_is_outside, current_point, is_outside); + if !prev_is_outside { + // Push previous point. + points.push(prev_point); + } + if prev_is_outside != is_outside { + if let Some(intersection_point) = intersect_plane_edge(prev_point, current_point) { + // Push intersection point. + intersection_points.push(intersection_point); + points.push(intersection_point); + } + } + /* let prev_is_total = mem::replace( + &mut is_total, + // Save the intersection point only if we go from outside to inside or inside to + // outside, and definitely intersect the plane edge. + prev_is_outside != is_outside && + + .map(|intersection_point| { + intersection_points.push(intersection_point); + if prev_is_outside { + // If the previous point is outside, we know + *point = intersection_point; + } else { + // i o i o + // + // i o (2) + // i i/o o/i (3) + // + // i o i (3) + // i i/o o/i i (4) + // + // i o i o (4) + // i i/o o/i i i/o o/i (6) + // + // i o i o i (5) + // i i/o o/i i i/o o/i i (7) + // + // i o i o i o (6) + // i i/o o/i i i/o o/i i i/o o/i (9) + current_point = intersection_point; + } + false + }) + .is_none(), + ); + // Save the previous point if it is either inside, or has been replaced by an intersection + // point. + !prev_is_outside || prev_is_total + /* match (prev_is_outside, is_outside) { + (true, true) => { + prev_is_total + }, + (true, false) => { + // Outside to inside, so save the previous point only if it's been replaced by an + // intersection point. + do_intersection(); + prev_is_total + }, + (false, true) => { + // Inside to outside, so always save the previous point, and save the intersection + // point only if we definitively intersect the plane edge. + false + }, + (false, false) => { + // Both points inside the plane, so always save the previous point. + false + } + } */ */ + }); + /* if !(is_total || is_outside) { + points.push(current_point); + } + /* match (before_plane, is_outside) { + (true, Previous::Outside) => { + + } + } + let cur_is_outside = { + if let Intersection::Inside = is_outside { + } else { + } + let prev_is_outside = mem::replace(&mut is_outside, { + let if let Intersection::Inside = is_outside { + true + } else { + false + } point_before_plane(current_point) { + }); + match (prev_is_outside, is_outside) { + (true, Some(is_outside)) => { + // Both points outside the plane, so save the previous point only if it's been + // replaced by an intersection point. + is_outside + }, + (true, false) => { + // Outside to inside, so calculate the intersection, and save it. + intersect_points.push(*point); + false + }, + (false, true) => { + // Inside to outside, so calculate the intersection, and save it and the current + // point. + intersect_points.push(*point); + false + }, + (false, false) => { + // Both points inside the plane, so save previous point + *point = * + false + } + } + if is_outside { + if prev_is_outside { + } else { + } + } else { + if prev_is_outside { + } + } + });*/ }*/ + last_is_outside +} + +fn append_intersection_points( + polys: &mut Vec>>, + intersection_points: Vec>, + tolerance: T, +) { + // NOTE: We use decoded versions of each polygon, with rounded entries. + // + // The line segments in intersection_points are consistently ordered as follows: + // each segment represents the intersection of the cutting plane with the + // polygon from which the segment came. The polygon can thus be split into + // two parts: the part "inside" the new surface (below the plane), and the + // part "outside" it (above the plane). Thus, when oriented + // with the polygon normal pointing into the camera, and the cutting plane as + // the x axis, with the "outside" part on top and the "inside" part on the + // bottom, there is a leftmost point (the point going from outside to + // inside, counterclockwise) and a rightmost point (the point going from + // inside to outside, counterclockwise). Our consistent ordering guarantees + // that the leftmost point comes before the rightmost point in each line + // segment. + // + // Why does this help us? To see that, consider the polygon adjacent to the + // considered polygon which also has the same right intersection point (we + // know there will be exactly one of these, because we have a solid + // structure and are only considering polygons that intersect the plane + // exactly two times; this means that we are ignoring polygons that intersect + // the plane at just one point, which means the two polygons must share a + // point, not be coplanar, and both intersect the plane; no theorem here, + // but I believe there can provably be at most one such instance given that + // we have at least three polygons with such a line segment). + // + // Now, for the adjacent polygon, repeat the above process. If the intersection + // point shared by the polygons is on the right in both cases, then we can + // see that the polygon's normal must be facing in the opposite direction of + // the original polygon despite being adjacent. But this + // should be impossible for a closed object! The same applies to the leftmost + // point. + // + // What is the practical upshot of all this? It means that we can consistently + // hash each line segment by its first point, which we can look up using the + // second point of a previous line segment. This will produce a chain of + // entries terminating in the original segment we looked up. As an added + // bonus, by going from leftmost point to leftmost point, we also ensure that + // we produce a polygon whose face is oriented counterclockwise around its + // normal; this can be seen by following the right-hand rule (TODO: provide + // more rigorous proof). + let tol = tolerance.recip(); + let make_key = move |point: Vec3| { + // We use floating points rounded to tolerance in order to make our HashMap + // lookups work. Otherwise we'd have to use a sorted structure, like a + // btree, which wouldn't be the end of the world but would have + // theoretically worse complexity. NOTE: Definitely non-ideal that we + // panic if the rounded value can't fit in an i64... TODO: If necessary, + // let the caller specify how to hash these keys, since in cases where + // we know the kind of floating point we're using we can just cast to bits or + // something. + point.map(|e| { + (e * tol) + .round() + .to_i64() + .expect("We don't currently try to handle floats that won't fit in an i64.") + }) + }; + let mut lines_iter = intersection_points.chunks_exact(2).filter_map(|uv| { + let u_key = make_key(uv[0]); + let v = uv[1]; + // NOTE: The reason we need to make sure this doesn't happen is that it's + // otherwise possible for two points to hash to the same value due to + // epsilon being too low. Because of the ordering mentioned previously, + // we know we should *eventually* find a pair of points starting with + // make_key(u) and ending with a different make_key(v) in such cases, so + // we just discard all the other pairs (treating them as points rather + // than lines). + (u_key != make_key(v)).then_some((u_key, v)) + }); + // .map(|uv| (make_key(uv[0]), uv[1])) + + if let Some((last_key, first)) = lines_iter.next() + /* [last, first, rest @ ..] = &*intersection_points = &*intersection_points */ + { + let lines = lines_iter.collect::>(); + /* if rest.len() < 4 { + // You need at least 3 sides for a polygon + return; + } + let lines = rest + .chunks_exact(2) + .filter_map(|uv| { + let u_key = make_key(uv[0]); + let v = uv[1]; + (u_key != make_key(v)).then_some((u_key, v)) + }) + // .map(|uv| (make_key(uv[0]), uv[1])) + .collect::>(); */ + if lines.len() < 2 { + // You need at least 3 sides for a polygon + return; + } + // println!("lines@append_intersection_points before merging points (last={:?}, + // cur={:?}): {:?}", last, cur, lines); + // let mut poly = Vec::with_capacity(lines.len() + 1); + // poly.push(first); + // NOTE: Guaranteed to terminate, provided we have no cycles besides the one + // that touches every point (which should be the case given how these + // points were generated). + let /*mut */poly_iter = iter::successors(Some(first), |&cur| lines.get(&make_key(cur)).copied()); + /* poly.extend(poly_iter.next()); + // TODO: If we were smart and pre-tested whether (last, first) was a dup (guaranteeing we + // started on a non-dup), we would not need the take_while part. + poly.extend(poly_iter.take_while(|&cur| make_key(cur) != make_key(first))); + /* while let Some(&v) = lines.get(&make_key(cur)) { + cur = v; + poly.push(cur); + } */ */ + let poly: Vec<_> = poly_iter.collect(); + // We have to check to make sure we really went through the whole cycle. + // TODO: Consider adaptively decreasing precision until we can make the cycle + // happen. + if poly.last().copied().map(make_key) == Some(last_key) { + // Push the new polygon onto the object. + polys.push(poly); + } + } +} + +pub fn clip_object_by_plane + core::fmt::Debug>( + polys: &mut Vec>>, + plane: (Vec3, T), + tolerance: T, +) { + let mut intersection_points = Vec::new(); + polys.drain_filter(|points| { + let len = intersection_points.len(); + let outside_first = clip_points_by_plane(points, plane, &mut intersection_points); + // println!("points@clip_object_by_plane after clipping by {:?} (outside_first={:?}, intersection_points={:?}): {:?}", plane, outside_first, intersection_points, points); + // Only remember intersections that are not coplanar with this side; i.e. those + // that have segment length 2. + if len + 2 != intersection_points.len() { + intersection_points.truncate(len); + } else if !outside_first { + // Order the two intersection points consistently, so that, when considered + // counterclockwise: + // - the first point goes from the exterior of the polygon (above the cutting + // plane) to its interior. + // - the second point goes from the interior of the polygon (below the cutting + // plane) to its exterior. + // the second is always going + // + // This allows us to uniquely map each line segment to an "owning" point (the + // one going from outside to inside), which happens to also point + // the segment in a counterclockwise direction around the new + // polygon normal composed of all the lines we clipped. + intersection_points.swap(len, len + 1); + } + // Remove polygon if it was clipped away + points.is_empty() + }); + // println!("polys@clip_object_by_plane after clipping by {:?} (before appending + // interection points {:?}): {:?}", plane, intersection_points, polys); + // Add a polygon of all intersection points with the plane to close out the + // object. + append_intersection_points(polys, intersection_points, tolerance); +} + +pub fn clip_object_by_aabb + core::fmt::Debug>( + polys: &mut Vec>>, + bounds: Aabb, + tolerance: T, +) { + let planes = aabb_to_planes(bounds); + // println!("planes@clip_object_by_aabb: {:?}", planes); + planes.iter().for_each(|&plane| { + clip_object_by_plane(polys, plane, tolerance); + // println!("polys@clip_object_by_aabb (after clipping by {:?}): + // {:?}", plane, polys); + }); +} + +/// Return value is 'Some(segment)' if line segment intersects the current +/// test plane. Otherwise 'None' is returned in which case the line +/// segment is entirely clipped. +pub fn clip_test(p: T, q: T, (u1, u2): (T, T)) -> Option<(T, T)> { + /* let res = */ + if p == T::zero() { + if q >= T::zero() { Some((u1, u2)) } else { None } + } else { + let r = q / p; + if p < T::zero() { + if r > u2 { + None + } else { + Some((if r > u1 { r } else { u1 }, u2)) + } + } else { + if r < u1 { + None + } else { + Some((u1, if r < u2 { r } else { u2 })) + } + } + } /*; + // println!("clip_test@(p={:?}, q={:?}, (u1, u2)=({:?}. {:?})): + // res={:?}", p, q, u1, u2, res); res*/ +} + +pub fn intersection_line_aabb + core::fmt::Debug>( + p: Vec3, + dir: Vec3, + bounds: Aabb, +) -> Option> { + // println!("before@intersection_line_aabb: p={:?} dir={:?} bounds={:?}", p, + // dir, bounds); + /* let res = */ + clip_test(-dir.z, p.z - bounds.min.z, (T::zero(), T::infinity())) + .and_then(|t| clip_test(dir.z, bounds.max.z - p.z, t)) + .and_then(|t| clip_test(-dir.y, p.y - bounds.min.y, t)) + .and_then(|t| clip_test(dir.y, bounds.max.y - p.y, t)) + .and_then(|t| clip_test(-dir.x, p.x - bounds.min.x, t)) + .and_then(|t| clip_test(dir.x, bounds.max.x - p.x, t)) + .and_then(|(t1, t2)| { + if T::zero() <= t2 { + Some(dir.mul_add(Vec3::broadcast(t2), p)) + } else if T::zero() <= t1 { + Some(dir.mul_add(Vec3::broadcast(t1), p)) + } else { + None + } + }) /*; + //println!("after@intersection_line_aabb (p={:?} dir={:?} bounds={:?}): + // {:?}", p, dir, bounds, res); res */ +} + +pub fn include_object_light_volume< + T: Float + MulAdd + core::fmt::Debug, + I: Iterator>, +>( + obj: I, + light_dir: Vec3, + bounds: Aabb, +) -> impl Iterator> { + /* obj.filter_map(move |pt| intersection_line_aabb(pt, -light_dir, bounds)) */ + // obj.map(move |pt| intersection_line_aabb(pt, -light_dir, + // bounds).unwrap_or(pt)) + obj.flat_map(move |pt| iter::once(pt).chain(intersection_line_aabb(pt, -light_dir, bounds))) +} + +pub fn calc_focused_light_volume_points + core::fmt::Debug>( + inv_proj_view: Mat4, + _light_dir: Vec3, + scene_bounding_box: Aabb, + tolerance: T, +) -> impl Iterator> { + let world_pts = calc_view_frustum_world_coord(inv_proj_view); + // println!("world_pts: {:?}", world_pts); + let mut world_frust_object = calc_view_frust_object(&world_pts); + // println!("world_frust_object: {:?}", world_frust_object); + clip_object_by_aabb(&mut world_frust_object, scene_bounding_box, tolerance); + // println!("world_frust_object@clip_object_by_aabb: {:?}", world_frust_object); + /* let object_points = world_frust_object.into_iter().flat_map(|e| e.into_iter()); + object_points.clone().chain(include_object_light_volume(object_points, light_dir, scene_bounding_box)) */ + world_frust_object.into_iter().flat_map(|e| e.into_iter()) + /* include_object_light_volume( + world_frust_object.into_iter().flat_map(|e| e.into_iter()), + light_dir, + scene_bounding_box, + ) */ +} + +/// NOTE: Will not yield useful results if pts is empty! +pub fn fit_psr< + T: Float + MulAdd, + I: Iterator>, + F: FnMut(Vec4) -> Vec4, +>( + mat: Mat4, + pts: I, + mut do_p: F, +) -> Aabb { + let mut min = Vec4::broadcast(T::infinity()); + let mut max = Vec4::broadcast(T::neg_infinity()); + pts.map(|p| do_p(mat * Vec4::::from_point(p))) + .for_each(|p| { + min = Vec4::partial_min(min, p); + max = Vec4::partial_max(max, p); + }); + Aabb { + min: min.xyz(), + max: max.xyz(), + } + /* let mut make_p = |x: f32, y: f32, z: f32| -> Vec3 { + do_p(mat * Vec4::new(x, y, z, 1.0)) + }; + let p1 = make_p(bounds.min.x, bounds.min.y, bounds.min.z); + let p2 = make_p(bounds.max.x, bounds.min.y, bounds.min.z); + let p3 = make_p(bounds.min.x, bounds.max.y, bounds.min.z); + let p4 = make_p(bounds.max.x, bounds.max.y, bounds.min.z); + let p5 = make_p(bounds.min.x, bounds.min.y, bounds.max.z); + let p6 = make_p(bounds.max.x, bounds.min.y, bounds.max.z); + let p7 = make_p(bounds.min.x, bounds.max.y, bounds.max.z); + let p8 = make_p(bounds.max.x, bounds.max.y, bounds.max.z); + // let p1: Vec4 = mat * Vec4::new(bounds.min.x, bounds.min.y, bounds.min.z, 1.0); + // let p2: Vec4 = mat * Vec4::new(0.0, bounds.min.y, 0.0, 1.0); + // let p3: Vec4 = mat * Vec4::new(0.0, 0.0, bounds.min.z, 1.0); + // let p4: Vec4 = mat * Vec4::new(bounds.max.x, 0.0, 0.0, 1.0); + // let p5: Vec4 = mat * Vec4::new(0.0, bounds.max.y, 0.0, 1.0); + // let p6: Vec4 = mat * Vec4::new(bounds.max.x, bounds.max.y, bounds.max.z, 1.0); + // println!("p1 p6 {:?} {:?}", p1, p6); + // let xmin = p1.x.min(p6.x); + // let xmax = p1.x.max(p6.x); + // println!("p1 p2 p3 p4 p5 p6: {:?} {:?} {:?} {:?} {:?} {:?}", p1, p2, p3, p4, p5, p6); + let xmin = p1.x.min(p2.x.min(p3.x.min(p4.x.min(p5.x.min(p6.x.min(p7.x.min(p8.x))))))); + let xmax = p1.x.max(p2.x.max(p3.x.max(p4.x.max(p5.x.max(p6.x.max(p7.x.max(p8.x))))))); + // let xmin = p1.x.min(p2.x.min(p3.x.min(p4.x.min(p5.x.min(p6.x))))); + // let xmax = p1.x.max(p2.x.max(p3.x.max(p4.x.max(p5.x.max(p6.x))))); + // println!("xmin: {:?}, xmax: {:?}", xmin, xmax); + // let ymin = p1.y.min(p6.y); + // let ymax = p1.y.max(p6.y); + let ymin = p1.y.min(p2.y.min(p3.y.min(p4.y.min(p5.y.min(p6.y.min(p7.y.min(p8.y))))))); + let ymax = p1.y.max(p2.y.max(p3.y.max(p4.y.max(p5.y.max(p6.y.max(p7.y.max(p8.y))))))); + // println!("ymin: {:?}, ymax: {:?}", ymin, ymax); + + // let p1: Vec4 = view_mat * Vec4::new(scene_bounds.min.x, scene_bounds.min.y, scene_bounds.min.z, 1.0); + // let p2: Vec4 = view_mat * Vec4::new(0.0, scene_bounds.min.y, 0.0, 1.0); + // let p3: Vec4 = view_mat * Vec4::new(0.0, 0.0, scene_bounds.min.z, 1.0); + // let p4: Vec4 = view_mat * Vec4::new(scene_bounds.max.x, scene_bounds.max.y, scene_bounds.max.z, 1.0); + // let p5: Vec4 = view_mat * Vec4::new(0.0, scene_bounds.max.y, 0.0, 1.0); + // let p6: Vec4 = view_mat * Vec4::new(0.0, 0.0, scene_bounds.max.z, 1.0); + // println!("p1 p2 p3 p4 p5 p6: {:?} {:?} {:?} {:?} {:?} {:?}", p1, p2, p3, p4, p5, p6); + // println!("p1 p4 {:?} {:?}", p1, p4); + let zmin = p1.z.min(p2.z.min(p3.z.min(p4.z.min(p5.z.min(p6.z.min(p7.z.min(p8.z))))))); + let zmax = p1.z.max(p2.z.max(p3.z.max(p4.z.max(p5.z.max(p6.z.max(p7.z.max(p8.z))))))); + Aabb { + min: Vec3::new(xmin, ymin, zmin), + max: Vec3::new(xmax, ymax, zmax), + } */ +} diff --git a/voxygen/src/scene/mod.rs b/voxygen/src/scene/mod.rs index 688f076f91..133dc76601 100644 --- a/voxygen/src/scene/mod.rs +++ b/voxygen/src/scene/mod.rs @@ -1,6 +1,7 @@ pub mod camera; pub mod figure; pub mod lod; +pub mod math; pub mod simple; pub mod terrain; @@ -27,8 +28,6 @@ use common::{ terrain::{BlockKind, TerrainChunk}, vol::ReadVol, }; -use core::{iter, mem}; -use hashbrown::HashMap; use num::traits::{Float, FloatConst}; use specs::{Entity as EcsEntity, Join, WorldExt}; use vek::*; @@ -42,8 +41,8 @@ const NUM_DIRECTED_LIGHTS: usize = 1; const LIGHT_DIST_RADIUS: f32 = 64.0; // The distance beyond which lights may not emit light from their origin const SHADOW_DIST_RADIUS: f32 = 8.0; const SHADOW_MAX_DIST: f32 = 96.0; // The distance beyond which shadows may not be visible -/// The minimum sin γ we will use before switching to uniform mapping. -const EPSILON_GAMMA: f64 = 0.25; +/* /// The minimum sin γ we will use before switching to uniform mapping. +const EPSILON_GAMMA: f64 = 0.25; */ // const NEAR_PLANE: f32 = 0.5; // const FAR_PLANE: f32 = 100000.0; @@ -110,555 +109,6 @@ impl<'a> SceneData<'a> { pub fn get_moon_dir(&self) -> Vec3 { Globals::get_moon_dir(self.state.get_time_of_day()) } } -pub fn aabb_to_points(bounds: Aabb) -> [Vec3; 8] { - [ - Vec3::new(bounds.min.x, bounds.min.y, bounds.min.z), - Vec3::new(bounds.max.x, bounds.min.y, bounds.min.z), - Vec3::new(bounds.max.x, bounds.max.y, bounds.min.z), - Vec3::new(bounds.min.x, bounds.max.y, bounds.min.z), - Vec3::new(bounds.min.x, bounds.min.y, bounds.max.z), - Vec3::new(bounds.max.x, bounds.min.y, bounds.max.z), - Vec3::new(bounds.max.x, bounds.max.y, bounds.max.z), - Vec3::new(bounds.min.x, bounds.max.y, bounds.max.z), - ] -} - -/// Each Vec4 should be interpreted as reprenting plane equation -/// -/// a(x - x0) + b(y - y0) + c(z - z0) = 0, i.e. -/// ax + by + cz - (a * x0 + b * y0 + c * z0) = 0, i.e. -/// ax + by + cz = (a * x0 + b * y0 + c * z0), i.e. -/// (lettiing d = a * x0 + b * y0 + c * z0) -/// ax + by + cz = d -/// -/// where d is the distance of the plane from the origin. -pub fn aabb_to_planes(bounds: Aabb) -> [(Vec3, T); 6] { - let zero = T::zero(); - let one = T::one(); - let bounds = bounds.map(|e| e.abs()); - [ - // bottom - (Vec3::new(zero, -one, zero), bounds.min.y), - // top - (Vec3::new(zero, one, zero), bounds.max.y), - // left - (Vec3::new(-one, zero, zero), bounds.min.x), - // right - (Vec3::new(one, zero, zero), bounds.max.x), - // near - (Vec3::new(zero, zero, -one), bounds.min.z), - // far - (Vec3::new(zero, zero, one), bounds.max.z), - ] -} - -pub fn mat_mul_points>( - mat: Mat4, - pts: &mut [Vec3], - mut do_p: impl FnMut(Vec4) -> Vec3, -) { - pts.into_iter().for_each(|p| { - *p = do_p(mat * Vec4::from_point(*p)); - }); -} - -/// NOTE: Expects points computed from aabb_to_points. -pub fn calc_view_frust_object(pts: &[Vec3; 8]) -> Vec>> { - vec![ - // near (CCW) - vec![pts[0], pts[1], pts[2], pts[3]], - // far (CCW) - vec![pts[7], pts[6], pts[5], pts[4]], - // left (CCW) - vec![pts[0], pts[3], pts[7], pts[4]], - // right (CCW) - vec![pts[1], pts[5], pts[6], pts[2]], - // bottom (CCW) - vec![pts[4], pts[5], pts[1], pts[0]], - // top (CCW) - vec![pts[6], pts[7], pts[3], pts[2]], - ] -} - -pub fn calc_view_frustum_world_coord>( - inv_proj_view: Mat4, -) -> [Vec3; 8] { - let mut world_pts = aabb_to_points(Aabb { - min: -Vec3::one(), - max: Vec3::one(), - }); - mat_mul_points(inv_proj_view, &mut world_pts, |p| Vec3::from(p) / p.w); - world_pts -} - -pub fn point_plane_distance(point: Vec3, (norm, dist): (Vec3, T)) -> T { - norm.dot(point) - dist -} - -pub fn point_before_plane(point: Vec3, plane: (Vec3, T)) -> bool { - point_plane_distance(point, plane) > T::zero() -} - -/// Returns true if and only if the final point in the polygon (i.e. the first -/// point added to the new polygon) is outside the clipping plane (this implies -/// that the polygon must be non-degenerate). -pub fn clip_points_by_plane + core::fmt::Debug>( - points: &mut Vec>, - plane: (Vec3, T), - intersection_points: &mut Vec>, -) -> bool { - /* enum Intersection { - /// Previous point was inside the plane. - Inside, - /// Previous line segment was completely outside the plane. - Outside, - /// Previous line segment went from inside the plane to outside it. - InsideOut, - } */ - // println!("points@clip_points_by_plane before clipping by {:?}: {:?}", plane, - // points); - if points.len() < 3 { - return false; - } - // NOTE: Guaranteed to succeed since points.len() > 3. - let mut current_point = points[points.len() - 1]; - let (norm, dist) = plane; - let intersect_plane_edge = |a, b| { - let diff = b - a; - let t = norm.dot(diff); - if t == T::zero() { - None - } else { - let t = (dist - norm.dot(a)) / t; - if t < T::zero() || T::one() < t { - None - } else { - Some(diff.mul_add(Vec3::broadcast(t), a)) - } - } - }; - let last_is_outside = point_before_plane(current_point, plane); - let mut is_outside = last_is_outside; - /* // Might not actually be total, but if it is partial and the point is inside it will be - // written regardless, and if it is partial and the point is outside, it means the - // second-to-last point is inside; thus, the second-to-last point will be written regardless, - // current_point will hold the new intersection point, and is_total will be false, when the - // loop ends; thus all we need to do to take this case into account is to push current_point - // onto the points vector if (is_total || is_outside) is false at the end of the loop. - let mut is_total = true; */ - let mut old_points = Vec::with_capacity((3 * points.len()) / 2); - mem::swap(&mut old_points, points); - old_points.into_iter().for_each(|point| { - /* let prev_point = current_point; - // Swap point i with the previous point in the polygon, so it is the one we normally save - // when we return false. - mem::swap(&mut current_point, point); */ - let prev_point = mem::replace(&mut current_point, point); - /* if point_before_plane(current_point) { - // If we are an outside point, we should only calculate an intersection if the previous - // point was inside. - if - is_outside s - // point was outside. - } else { - // If we are an inside point, then we should only calculate an intersection if the previous - // point was outside. - } */ - let before_plane = point_before_plane(current_point, plane); - let prev_is_outside = mem::replace(&mut is_outside, before_plane); - // println!("points@clip_points_by_plane clipping segment by {:?} (prev={:?} / - // outside={:?}, current={:?} / outside={:?})", plane, prev_point, - // prev_is_outside, current_point, is_outside); - if !prev_is_outside { - // Push previous point. - points.push(prev_point); - } - if prev_is_outside != is_outside { - if let Some(intersection_point) = intersect_plane_edge(prev_point, current_point) { - // Push intersection point. - intersection_points.push(intersection_point); - points.push(intersection_point); - } - } - /* let prev_is_total = mem::replace( - &mut is_total, - // Save the intersection point only if we go from outside to inside or inside to - // outside, and definitely intersect the plane edge. - prev_is_outside != is_outside && - - .map(|intersection_point| { - intersection_points.push(intersection_point); - if prev_is_outside { - // If the previous point is outside, we know - *point = intersection_point; - } else { - // i o i o - // - // i o (2) - // i i/o o/i (3) - // - // i o i (3) - // i i/o o/i i (4) - // - // i o i o (4) - // i i/o o/i i i/o o/i (6) - // - // i o i o i (5) - // i i/o o/i i i/o o/i i (7) - // - // i o i o i o (6) - // i i/o o/i i i/o o/i i i/o o/i (9) - current_point = intersection_point; - } - false - }) - .is_none(), - ); - // Save the previous point if it is either inside, or has been replaced by an intersection - // point. - !prev_is_outside || prev_is_total - /* match (prev_is_outside, is_outside) { - (true, true) => { - prev_is_total - }, - (true, false) => { - // Outside to inside, so save the previous point only if it's been replaced by an - // intersection point. - do_intersection(); - prev_is_total - }, - (false, true) => { - // Inside to outside, so always save the previous point, and save the intersection - // point only if we definitively intersect the plane edge. - false - }, - (false, false) => { - // Both points inside the plane, so always save the previous point. - false - } - } */ */ - }); - /* if !(is_total || is_outside) { - points.push(current_point); - } - /* match (before_plane, is_outside) { - (true, Previous::Outside) => { - - } - } - let cur_is_outside = { - if let Intersection::Inside = is_outside { - } else { - } - let prev_is_outside = mem::replace(&mut is_outside, { - let if let Intersection::Inside = is_outside { - true - } else { - false - } point_before_plane(current_point) { - }); - match (prev_is_outside, is_outside) { - (true, Some(is_outside)) => { - // Both points outside the plane, so save the previous point only if it's been - // replaced by an intersection point. - is_outside - }, - (true, false) => { - // Outside to inside, so calculate the intersection, and save it. - intersect_points.push(*point); - false - }, - (false, true) => { - // Inside to outside, so calculate the intersection, and save it and the current - // point. - intersect_points.push(*point); - false - }, - (false, false) => { - // Both points inside the plane, so save previous point - *point = * - false - } - } - if is_outside { - if prev_is_outside { - } else { - } - } else { - if prev_is_outside { - } - } - });*/ }*/ - last_is_outside -} - -fn append_intersection_points( - polys: &mut Vec>>, - intersection_points: Vec>, - tolerance: T, -) { - // NOTE: We use decoded versions of each polygon, with rounded entries. - // - // The line segments in intersection_points are consistently ordered as follows: - // each segment represents the intersection of the cutting plane with the - // polygon from which the segment came. The polygon can thus be split into - // two parts: the part "inside" the new surface (below the plane), and the - // part "outside" it (above the plane). Thus, when oriented - // with the polygon normal pointing into the camera, and the cutting plane as - // the x axis, with the "outside" part on top and the "inside" part on the - // bottom, there is a leftmost point (the point going from outside to - // inside, counterclockwise) and a rightmost point (the point going from - // inside to outside, counterclockwise). Our consistent ordering guarantees - // that the leftmost point comes before the rightmost point in each line - // segment. - // - // Why does this help us? To see that, consider the polygon adjacent to the - // considered polygon which also has the same right intersection point (we - // know there will be exactly one of these, because we have a solid - // structure and are only considering polygons that intersect the plane - // exactly two times; this means that we are ignoring polygons that intersect - // the plane at just one point, which means the two polygons must share a - // point, not be coplanar, and both intersect the plane; no theorem here, - // but I believe there can provably be at most one such instance given that - // we have at least three polygons with such a line segment). - // - // Now, for the adjacent polygon, repeat the above process. If the intersection - // point shared by the polygons is on the right in both cases, then we can - // see that the polygon's normal must be facing in the opposite direction of - // the original polygon despite being adjacent. But this - // should be impossible for a closed object! The same applies to the leftmost - // point. - // - // What is the practical upshot of all this? It means that we can consistently - // hash each line segment by its first point, which we can look up using the - // second point of a previous line segment. This will produce a chain of - // entries terminating in the original segment we looked up. As an added - // bonus, by going from leftmost point to leftmost point, we also ensure that - // we produce a polygon whose face is oriented counterclockwise around its - // normal; this can be seen by following the right-hand rule (TODO: provide - // more rigorous proof). - let tol = tolerance.recip(); - let make_key = move |point: Vec3| { - // We use floating points rounded to tolerance in order to make our HashMap - // lookups work. Otherwise we'd have to use a sorted structure, like a - // btree, which wouldn't be the end of the world but would have - // theoretically worse complexity. NOTE: Definitely non-ideal that we - // panic if the rounded value can't fit in an i64... TODO: If necessary, - // let the caller specify how to hash these keys, since in cases where - // we know the kind of floating point we're using we can just cast to bits or - // something. - point.map(|e| { - (e * tol) - .round() - .to_i64() - .expect("We don't currently try to handle floats that won't fit in an i64.") - }) - }; - let mut lines_iter = intersection_points.chunks_exact(2).filter_map(|uv| { - let u_key = make_key(uv[0]); - let v = uv[1]; - // NOTE: The reason we need to make sure this doesn't happen is that it's - // otherwise possible for two points to hash to the same value due to - // epsilon being too low. Because of the ordering mentioned previously, - // we know we should *eventually* find a pair of points starting with - // make_key(u) and ending with a different make_key(v) in such cases, so - // we just discard all the other pairs (treating them as points rather - // than lines). - (u_key != make_key(v)).then_some((u_key, v)) - }); - // .map(|uv| (make_key(uv[0]), uv[1])) - - if let Some((last_key, first)) = lines_iter.next() - /* [last, first, rest @ ..] = &*intersection_points = &*intersection_points */ - { - let lines = lines_iter.collect::>(); - /* if rest.len() < 4 { - // You need at least 3 sides for a polygon - return; - } - let lines = rest - .chunks_exact(2) - .filter_map(|uv| { - let u_key = make_key(uv[0]); - let v = uv[1]; - (u_key != make_key(v)).then_some((u_key, v)) - }) - // .map(|uv| (make_key(uv[0]), uv[1])) - .collect::>(); */ - if lines.len() < 2 { - // You need at least 3 sides for a polygon - return; - } - // println!("lines@append_intersection_points before merging points (last={:?}, - // cur={:?}): {:?}", last, cur, lines); - // let mut poly = Vec::with_capacity(lines.len() + 1); - // poly.push(first); - // NOTE: Guaranteed to terminate, provided we have no cycles besides the one - // that touches every point (which should be the case given how these - // points were generated). - let /*mut */poly_iter = iter::successors(Some(first), |&cur| lines.get(&make_key(cur)).copied()); - /* poly.extend(poly_iter.next()); - // TODO: If we were smart and pre-tested whether (last, first) was a dup (guaranteeing we - // started on a non-dup), we would not need the take_while part. - poly.extend(poly_iter.take_while(|&cur| make_key(cur) != make_key(first))); - /* while let Some(&v) = lines.get(&make_key(cur)) { - cur = v; - poly.push(cur); - } */ */ - let poly: Vec<_> = poly_iter.collect(); - // We have to check to make sure we really went through the whole cycle. - // TODO: Consider adaptively decreasing precision until we can make the cycle - // happen. - if poly.last().copied().map(make_key) == Some(last_key) { - // Push the new polygon onto the object. - polys.push(poly); - } - } -} - -pub fn clip_object_by_plane + core::fmt::Debug>( - polys: &mut Vec>>, - plane: (Vec3, T), - tolerance: T, -) { - let mut intersection_points = Vec::new(); - polys.drain_filter(|points| { - let len = intersection_points.len(); - let outside_first = clip_points_by_plane(points, plane, &mut intersection_points); - // println!("points@clip_object_by_plane after clipping by {:?} (outside_first={:?}, intersection_points={:?}): {:?}", plane, outside_first, intersection_points, points); - // Only remember intersections that are not coplanar with this side; i.e. those - // that have segment length 2. - if len + 2 != intersection_points.len() { - intersection_points.truncate(len); - } else if !outside_first { - // Order the two intersection points consistently, so that, when considered - // counterclockwise: - // - the first point goes from the exterior of the polygon (above the cutting - // plane) to its interior. - // - the second point goes from the interior of the polygon (below the cutting - // plane) to its exterior. - // the second is always going - // - // This allows us to uniquely map each line segment to an "owning" point (the - // one going from outside to inside), which happens to also point - // the segment in a counterclockwise direction around the new - // polygon normal composed of all the lines we clipped. - intersection_points.swap(len, len + 1); - } - // Remove polygon if it was clipped away - points.is_empty() - }); - // println!("polys@clip_object_by_plane after clipping by {:?} (before appending - // interection points {:?}): {:?}", plane, intersection_points, polys); - // Add a polygon of all intersection points with the plane to close out the - // object. - append_intersection_points(polys, intersection_points, tolerance); -} - -pub fn clip_object_by_aabb + core::fmt::Debug>( - polys: &mut Vec>>, - bounds: Aabb, - tolerance: T, -) { - let planes = aabb_to_planes(bounds); - // println!("planes@clip_object_by_aabb: {:?}", planes); - planes.iter().for_each(|&plane| { - clip_object_by_plane(polys, plane, tolerance); - // println!("polys@clip_object_by_aabb (after clipping by {:?}): {:?}", - // plane, polys); - }); -} - -/// Return value is 'Some(segment)' if line segment intersects the current test -/// plane. Otherwise 'None' is returned in which case the line segment -/// is entirely clipped. -pub fn clip_test(p: T, q: T, (u1, u2): (T, T)) -> Option<(T, T)> { - /* let res = */ - if p == T::zero() { - if q >= T::zero() { Some((u1, u2)) } else { None } - } else { - let r = q / p; - if p < T::zero() { - if r > u2 { - None - } else { - Some((if r > u1 { r } else { u1 }, u2)) - } - } else { - if r < u1 { - None - } else { - Some((u1, if r < u2 { r } else { u2 })) - } - } - } /*; - // println!("clip_test@(p={:?}, q={:?}, (u1, u2)=({:?}. {:?})): res={:?}", - // p, q, u1, u2, res); res*/ -} - -pub fn intersection_line_aabb + core::fmt::Debug>( - p: Vec3, - dir: Vec3, - bounds: Aabb, -) -> Option> { - // println!("before@intersection_line_aabb: p={:?} dir={:?} bounds={:?}", p, - // dir, bounds); - /* let res = */ - clip_test(-dir.z, p.z - bounds.min.z, (T::zero(), T::infinity())) - .and_then(|t| clip_test(dir.z, bounds.max.z - p.z, t)) - .and_then(|t| clip_test(-dir.y, p.y - bounds.min.y, t)) - .and_then(|t| clip_test(dir.y, bounds.max.y - p.y, t)) - .and_then(|t| clip_test(-dir.x, p.x - bounds.min.x, t)) - .and_then(|t| clip_test(dir.x, bounds.max.x - p.x, t)) - .and_then(|(t1, t2)| { - if T::zero() <= t2 { - Some(dir.mul_add(Vec3::broadcast(t2), p)) - } else if T::zero() <= t1 { - Some(dir.mul_add(Vec3::broadcast(t1), p)) - } else { - None - } - }) /*; - //println!("after@intersection_line_aabb (p={:?} dir={:?} bounds={:?}): - // {:?}", p, dir, bounds, res); res */ -} - -pub fn include_object_light_volume< - T: Float + MulAdd + core::fmt::Debug, - I: Iterator>, ->( - obj: I, - light_dir: Vec3, - bounds: Aabb, -) -> impl Iterator> { - /* obj.filter_map(move |pt| intersection_line_aabb(pt, -light_dir, bounds)) */ - // obj.map(move |pt| intersection_line_aabb(pt, -light_dir, - // bounds).unwrap_or(pt)) - obj.flat_map(move |pt| iter::once(pt).chain(intersection_line_aabb(pt, -light_dir, bounds))) -} - -pub fn calc_focused_light_volume_points + core::fmt::Debug>( - inv_proj_view: Mat4, - _light_dir: Vec3, - scene_bounding_box: Aabb, - tolerance: T, -) -> impl Iterator> { - let world_pts = calc_view_frustum_world_coord(inv_proj_view); - // println!("world_pts: {:?}", world_pts); - let mut world_frust_object = calc_view_frust_object(&world_pts); - // println!("world_frust_object: {:?}", world_frust_object); - clip_object_by_aabb(&mut world_frust_object, scene_bounding_box, tolerance); - // println!("world_frust_object@clip_object_by_aabb: {:?}", world_frust_object); - /* let object_points = world_frust_object.into_iter().flat_map(|e| e.into_iter()); - object_points.clone().chain(include_object_light_volume(object_points, light_dir, scene_bounding_box)) */ - world_frust_object.into_iter().flat_map(|e| e.into_iter()) - /* include_object_light_volume( - world_frust_object.into_iter().flat_map(|e| e.into_iter()), - light_dir, - scene_bounding_box, - ) */ -} - /// Approximte a scalar field of view angle using the parameterization from /// section 4.3 of Lloyd's thesis: /// @@ -781,72 +231,6 @@ fn compute_warping_parameter_perspective( ) } -/// NOTE: Will not yield useful results if pts is empty! -pub fn fit_psr< - T: Float + MulAdd, - I: Iterator>, - F: FnMut(Vec4) -> Vec3, ->( - mat: Mat4, - pts: I, - mut do_p: F, -) -> Aabb { - let mut min = Vec3::broadcast(T::infinity()); - let mut max = Vec3::broadcast(T::neg_infinity()); - pts.map(|p| do_p(mat * Vec4::::from_point(p))) - .for_each(|p| { - min = Vec3::partial_min(min, p); - max = Vec3::partial_max(max, p); - }); - Aabb { min, max } - /* let mut make_p = |x: f32, y: f32, z: f32| -> Vec3 { - do_p(mat * Vec4::new(x, y, z, 1.0)) - }; - let p1 = make_p(bounds.min.x, bounds.min.y, bounds.min.z); - let p2 = make_p(bounds.max.x, bounds.min.y, bounds.min.z); - let p3 = make_p(bounds.min.x, bounds.max.y, bounds.min.z); - let p4 = make_p(bounds.max.x, bounds.max.y, bounds.min.z); - let p5 = make_p(bounds.min.x, bounds.min.y, bounds.max.z); - let p6 = make_p(bounds.max.x, bounds.min.y, bounds.max.z); - let p7 = make_p(bounds.min.x, bounds.max.y, bounds.max.z); - let p8 = make_p(bounds.max.x, bounds.max.y, bounds.max.z); - // let p1: Vec4 = mat * Vec4::new(bounds.min.x, bounds.min.y, bounds.min.z, 1.0); - // let p2: Vec4 = mat * Vec4::new(0.0, bounds.min.y, 0.0, 1.0); - // let p3: Vec4 = mat * Vec4::new(0.0, 0.0, bounds.min.z, 1.0); - // let p4: Vec4 = mat * Vec4::new(bounds.max.x, 0.0, 0.0, 1.0); - // let p5: Vec4 = mat * Vec4::new(0.0, bounds.max.y, 0.0, 1.0); - // let p6: Vec4 = mat * Vec4::new(bounds.max.x, bounds.max.y, bounds.max.z, 1.0); - // println!("p1 p6 {:?} {:?}", p1, p6); - // let xmin = p1.x.min(p6.x); - // let xmax = p1.x.max(p6.x); - // println!("p1 p2 p3 p4 p5 p6: {:?} {:?} {:?} {:?} {:?} {:?}", p1, p2, p3, p4, p5, p6); - let xmin = p1.x.min(p2.x.min(p3.x.min(p4.x.min(p5.x.min(p6.x.min(p7.x.min(p8.x))))))); - let xmax = p1.x.max(p2.x.max(p3.x.max(p4.x.max(p5.x.max(p6.x.max(p7.x.max(p8.x))))))); - // let xmin = p1.x.min(p2.x.min(p3.x.min(p4.x.min(p5.x.min(p6.x))))); - // let xmax = p1.x.max(p2.x.max(p3.x.max(p4.x.max(p5.x.max(p6.x))))); - // println!("xmin: {:?}, xmax: {:?}", xmin, xmax); - // let ymin = p1.y.min(p6.y); - // let ymax = p1.y.max(p6.y); - let ymin = p1.y.min(p2.y.min(p3.y.min(p4.y.min(p5.y.min(p6.y.min(p7.y.min(p8.y))))))); - let ymax = p1.y.max(p2.y.max(p3.y.max(p4.y.max(p5.y.max(p6.y.max(p7.y.max(p8.y))))))); - // println!("ymin: {:?}, ymax: {:?}", ymin, ymax); - - // let p1: Vec4 = view_mat * Vec4::new(scene_bounds.min.x, scene_bounds.min.y, scene_bounds.min.z, 1.0); - // let p2: Vec4 = view_mat * Vec4::new(0.0, scene_bounds.min.y, 0.0, 1.0); - // let p3: Vec4 = view_mat * Vec4::new(0.0, 0.0, scene_bounds.min.z, 1.0); - // let p4: Vec4 = view_mat * Vec4::new(scene_bounds.max.x, scene_bounds.max.y, scene_bounds.max.z, 1.0); - // let p5: Vec4 = view_mat * Vec4::new(0.0, scene_bounds.max.y, 0.0, 1.0); - // let p6: Vec4 = view_mat * Vec4::new(0.0, 0.0, scene_bounds.max.z, 1.0); - // println!("p1 p2 p3 p4 p5 p6: {:?} {:?} {:?} {:?} {:?} {:?}", p1, p2, p3, p4, p5, p6); - // println!("p1 p4 {:?} {:?}", p1, p4); - let zmin = p1.z.min(p2.z.min(p3.z.min(p4.z.min(p5.z.min(p6.z.min(p7.z.min(p8.z))))))); - let zmax = p1.z.max(p2.z.max(p3.z.max(p4.z.max(p5.z.max(p6.z.max(p7.z.max(p8.z))))))); - Aabb { - min: Vec3::new(xmin, ymin, zmin), - max: Vec3::new(xmax, ymax, zmax), - } */ -} - impl Scene { /// Create a new `Scene` with default parameters. pub fn new(renderer: &mut Renderer, client: &Client, settings: &Settings) -> Self { @@ -1181,11 +565,11 @@ impl Scene { }; */ // let focus_frac = focus_pos.map(|e| e.fract()); - let visible_bounds = Aabb { - min: visible_bounds.min - focus_off, - max: visible_bounds.max - focus_off, + let visible_bounds = math::Aabb:: { + min: math::Vec3::from(visible_bounds.min - focus_off), + max: math::Vec3::from(visible_bounds.max - focus_off), }; - let visible_bounds_fine = Aabb { + let visible_bounds_fine = math::Aabb { min: visible_bounds.min.map(f64::from), max: visible_bounds.max.map(f64::from), }; @@ -1205,9 +589,11 @@ impl Scene { min: scene_bounds.min.map(f64::from), max: scene_bounds.max.map(f64::from), }; */ - let inv_proj_view = (proj_mat * view_mat/* * Mat4::translation_3d(-focus_off)*/) - .map(f64::from) - .inverted(); + let inv_proj_view = math::Mat4::from_col_arrays( + (proj_mat * view_mat/* * Mat4::translation_3d(-focus_off)*/).into_col_arrays(), + ) + .map(f64::from) + .inverted(); let fov = self.camera.get_fov(); let aspect_ratio = self.camera.get_aspect_ratio(); @@ -1221,7 +607,7 @@ impl Scene { let point_shadow_aspect = point_shadow_res.x as f32 / point_shadow_res.y as f32; // Construct matrices to transform from world space to light space for the sun // and moon. - let directed_light_dir = sun_dir; + let directed_light_dir = math::Vec3::from(sun_dir); /* let light_volume = calc_focused_light_volume_points(inv_proj_view, directed_light_dir.map(f64::from), scene_bounds_fine, 1e-3) // .map(|e| e - focus_off) // NOTE: Hopefully not out of bounds. @@ -1229,7 +615,7 @@ impl Scene { .collect::>(); // println!("light_volume: {:?}", light_volume); */ // let visible_light_volume = light_volume.clone(); - let visible_light_volume = calc_focused_light_volume_points(inv_proj_view, directed_light_dir.map(f64::from), visible_bounds_fine, 1e-6) + let visible_light_volume = math::calc_focused_light_volume_points(inv_proj_view, directed_light_dir.map(f64::from), visible_bounds_fine, 1e-6) // .map(|e| e - focus_off) // NOTE: Hopefully not out of bounds. .map(|v| v.map(|e| e as f32)) @@ -1304,12 +690,12 @@ impl Scene { // let look_at = bounds0.center();//Vec3::zero();// // scene_bounds.center();//Vec3::zero(); let look_at = // bounds0.center(); - let look_at = cam_pos; // /*Vec3::zero()*/scene_bounds.center()/*cam_pos*/;// - focus_off;// focus_off; + let look_at = math::Vec3::from(cam_pos); // /*Vec3::zero()*/scene_bounds.center()/*cam_pos*/;// - focus_off;// focus_off; let _light_scale = 1.5 * /*(directed_near + directed_far) / 2.0*/radius; // We upload view matrices as well, to assist in linearizing vertex positions. // (only for directional lights, so far). let mut directed_shadow_mats = Vec::with_capacity(6); - let new_dir = view_dir; + let new_dir = math::Vec3::from(view_dir); // let new_dir: Vec3 = light_volume/*visible_light_volume*/.iter().map(|p| // p - cam_pos).sum(); let new_dir = new_dir.normalized(); @@ -1324,18 +710,18 @@ impl Scene { Vec3::from(view_mat * Vec4::from_direction(Vec3::up())) .normalized() }; */ - let up: Vec3 = { + let up: math::Vec3 = { /* (directed_light_dir) .cross(new_dir) .cross(directed_light_dir) .normalized() */ - Vec3::up() + math::Vec3::up() }; // let up = Vec3::up(); // let up: Vec3 = Vec3::from(Mat4::::look_at_rh(look_at - sun_dir, // look_at, -Vec3::from(view_dir)) * Vec4::::forward_rh()); // println!("bounds0: {:?}, scene_bounds: {:?}", bounds0, scene_bounds); - directed_shadow_mats.push(Mat4::look_at_rh( + directed_shadow_mats.push(math::Mat4::look_at_rh( look_at, look_at + directed_light_dir, /* Vec3::up()*//*Vec3::from(view_dir)*//*up*//*Vec3::down() */ up, @@ -1346,7 +732,7 @@ impl Scene { // look_at, Vec3::up())); This leaves us with four dummy slots, // which we push as defaults. directed_shadow_mats - .extend_from_slice(&[Mat4::default(); 6 - NUM_DIRECTED_LIGHTS] as _); + .extend_from_slice(&[math::Mat4::default(); 6 - NUM_DIRECTED_LIGHTS] as _); // Now, construct the full projection matrices in the first two directed light // slots. let mut shadow_mats = Vec::with_capacity(6 * (lights.len() + 1)); @@ -1482,7 +868,7 @@ impl Scene { // let mut e_p: Vec3 = Vec3::zero(); v_p.z = 0.0; */ - let mut v_p = Vec3::from(light_view_mat * Vec4::from_direction(new_dir)); + let mut v_p = math::Vec3::from(light_view_mat * math::Vec4::from_direction(new_dir)); v_p.normalize(); // let dot_prod = f64::from(v_p.z); let dot_prod = new_dir.map(f64::from).dot(directed_light_dir.map(f64::from)); @@ -1497,10 +883,10 @@ impl Scene { v_p.z = 0.0; v_p.normalize(); - let l_r: Mat4 = if /*v_p.magnitude_squared() > 1e-3*//*sin_gamma > EPISLON_GAMMA*/factor != -1.0 { - Mat4::look_at_rh(Vec3::zero(), Vec3::forward_rh(), v_p) + let l_r: math::Mat4 = if /*v_p.magnitude_squared() > 1e-3*//*sin_gamma > EPISLON_GAMMA*/factor != -1.0 { + math::Mat4::look_at_rh(math::Vec3::zero(), math::Vec3::forward_rh(), v_p) } else { - Mat4::identity() + math::Mat4::identity() }; // let factor = -1.0; // let l_r: Mat4 = Mat4::look_at_rh(/*Vec3::from(e_p) - Vec3::from(v_p)*//*Vec3::up()*/e_p, /*Vec3::from(e_p)*//*Vec3::zero()*/e_p + Vec3::forward_rh(), Vec3::from(v_p)); @@ -1511,8 +897,8 @@ impl Scene { // let l_r: Mat4 = Mat4::look_at_rh(/*Vec3::from(e_p) - Vec3::from(v_p)*/Vec3::zero(), /*Vec3::from(e_p)*/-Vec3::forward_rh(), /*Vec3::up()*/-Vec3::from(v_p)); // let l_r: Mat4 = Mat4::look_at_rh(/*Vec3::from(e_p) - Vec3::from(v_p)*/Vec3::back_rh(), /*Vec3::from(e_p)*/Vec3::zero(), /*Vec3::up()*/Vec3::from(v_p)); // let l_r: Mat4 = Mat4::identity(); - let bounds0 = fit_psr(light_view_mat, visible_light_volume.iter().copied(), |p| Vec3::from(p) / p.w); - let directed_proj_mat = Mat4::orthographic_rh_no(FrustumPlanes { + let bounds0 = math::fit_psr(light_view_mat, visible_light_volume.iter().copied(), /*|p| math::Vec3::from(p) / p.w*/math::Vec4::homogenized); + let directed_proj_mat = math::Mat4::orthographic_rh_no(FrustumPlanes { // TODO: Consider adjusting resolution based on view distance. left: bounds0.min.x, right: bounds0.max.x, @@ -1524,10 +910,10 @@ impl Scene { let light_all_mat = l_r * directed_proj_mat * light_view_mat; // let bounds1 = fit_psr(light_all_mat/* * inverse_visible*/, light_volume.iter().copied(), |p| Vec3::from(p) / p.w); - let bounds0 = fit_psr(/*l_r*/light_all_mat/* * inverse_visible*/, visible_light_volume.iter().copied(), |p| Vec3::from(p) / p.w); + let bounds0 = math::fit_psr(/*l_r*/light_all_mat/* * inverse_visible*/, visible_light_volume.iter().copied(), /*|p| math::Vec3::from(p) / p.w*/math::Vec4::homogenized); // let bounds1 = fit_psr(light_all_mat/* * inverse_visible*/, aabb_to_points(visible_bounds).iter().copied(), |p| Vec3::from(p) / p.w); // let mut light_focus_pos: Vec3 = Vec3::from(light_all_mat * Vec4::from_point(focus_pos.map(f32::fract))); - let mut light_focus_pos: Vec3 = Vec3::zero();//bounds0.center();// l_r * directed_proj_mat * light_view_mat * Vec4::from_point(focus_pos.map(|e| e.fract())); + let mut light_focus_pos: math::Vec3 = math::Vec3::zero();//bounds0.center();// l_r * directed_proj_mat * light_view_mat * Vec4::from_point(focus_pos.map(|e| e.fract())); // let mut light_focus_pos: Vec3 = bounds0.center();// l_r * directed_proj_mat * light_view_mat * Vec4::from_point(focus_pos.map(|e| e.fract())); // println!("cam_pos: {:?}, focus_pos: {:?}, light_focus_pos: {:?}, v_p: {:?} bounds: {:?}, l_r: {:?}, light_view_mat: {:?}, light_all_mat: {:?}", cam_pos, focus_pos - focus_off, light_focus_pos, v_p, /*bounds1*/bounds0, l_r, light_view_mat, light_all_mat); // let w_v = Mat4::translation_3d(-Vec3::new(xmax + xmin, ymax + ymin, /*zmax + zmin*/0.0) / 2.0); @@ -1616,7 +1002,7 @@ impl Scene { } else { light_focus_pos.y }; - let w_v: Mat4 = Mat4::translation_3d(/*-bounds1.center()*/-Vec3::new(light_focus_pos.x, light_focus_pos.y/* + (directed_near - near_dist)*/,/* - /*(directed_near - near_dist)*/directed_near*//*bounds1.center().z*//*directed_near*//*bounds1.min.z - *//*(directed_near - near_dist)*//*focus_pos.z*//*light_focus_pos.z*//*light_focus_pos.z*//*center1.z*//*center1.z.max(0.0)*/light_focus_pos.z)); + let w_v: math::Mat4 = math::Mat4::translation_3d(/*-bounds1.center()*/-math::Vec3::new(light_focus_pos.x, light_focus_pos.y/* + (directed_near - near_dist)*/,/* - /*(directed_near - near_dist)*/directed_near*//*bounds1.center().z*//*directed_near*//*bounds1.min.z - *//*(directed_near - near_dist)*//*focus_pos.z*//*light_focus_pos.z*//*light_focus_pos.z*//*center1.z*//*center1.z.max(0.0)*/light_focus_pos.z)); // let w_v: Mat4 = Mat4::translation_3d(/*-bounds1.center()*/-Vec3::new(light_focus_pos.x, light_focus_pos.y,/* - /*(directed_near - near_dist)*/directed_near*//*bounds1.center().z*//*directed_near*//*bounds1.min.z - *//*(directed_near - near_dist)*//*focus_pos.z*//*light_focus_pos.z*//*light_focus_pos.z*/center1.z + directed_near - near_dist)); // let w_v: Mat4 = Mat4::translation_3d(/*-bounds1.center()*/-Vec3::new(0.0, 0.0,/* - /*(directed_near - near_dist)*/directed_near*//*bounds1.center().z*//*directed_near*//*bounds1.min.z - *//*(directed_near - near_dist)*//*focus_pos.z*//*light_focus_pos.z*/directed_near - near_dist)); /* let w_p: Mat4 = Mat4::orthographic_rh_no/*frustum_rh_no*/(FrustumPlanes { @@ -1634,10 +1020,10 @@ impl Scene { near: directed_near, far: directed_far,// directed_near + /*zmax - zmin*/bounds1.max.z - bounds1.min.z,//directed_far, }); */ - let shadow_view_mat: Mat4 = w_v * light_all_mat; - let _bounds0 = fit_psr(/*l_r*/shadow_view_mat/* * inverse_visible*/, visible_light_volume.iter().copied(), |p| Vec3::from(p) / p.w); + let shadow_view_mat: math::Mat4 = w_v * light_all_mat; + let _bounds0 = math::fit_psr(/*l_r*/shadow_view_mat/* * inverse_visible*/, visible_light_volume.iter().copied(), /*|p| math::Vec3::from(p) / p.w*/math::Vec4::homogenized); // let factor = -1.0; - let w_p: Mat4 = { + let w_p: math::Mat4 = { if /*sin_gamma > EPISLON_GAMMA*/factor != -1.0 { // Projection for y let n = directed_near;// - near_dist; @@ -1698,7 +1084,7 @@ impl Scene { // = f (2(x₀ (n / f) - l) / (r - l) - 1) // // x(f) = 2 (x₀ f + l y₀) / (r - l) - f - Mat4::new( + math::Mat4::new( s_x, o_x, 0.0, 0.0, 0.0, s_y, 0.0, o_y, 0.0, o_z, s_z, 0.0, @@ -1717,7 +1103,7 @@ impl Scene { 0.0, 0.0, s_y, o_y, 0.0, 0.0, 1.0, 0.0, ) */ - Mat4::identity() + math::Mat4::identity() } // Mat4::identity() /* let a = (n + f) / (n - f); @@ -1762,11 +1148,11 @@ impl Scene { near: zmin,//directed_near, far: zmax,//directed_far, }); */ - let shadow_all_mat: Mat4 = w_p * shadow_view_mat/*w_v * light_all_mat*/; + let shadow_all_mat: math::Mat4 = w_p * shadow_view_mat/*w_v * light_all_mat*/; let _w_p_arr = shadow_all_mat.cols.iter().map(|e| (e.x, e.y, e.z, e.w)).collect::>(); // println!("mat4 shadow_all_mat = mat4(vec4{:?}, vec4{:?}, vec4{:?}, vec4{:?});", w_p_arr[0], w_p_arr[1], w_p_arr[2], w_p_arr[3]); - let Aabb { min: Vec3 { x: xmin, y: ymin, z: zmin }, max: Vec3 { x: xmax, y: ymax, z: zmax } } = - fit_psr(/*light_all_mat*/shadow_all_mat/*shadow_view_mat*//* * inverse_visible*/, visible_light_volume.iter().copied(), |p| Vec3::from(p) / p.w); + let math::Aabb:: { min: math::Vec3 { x: xmin, y: ymin, z: zmin }, max: math::Vec3 { x: xmax, y: ymax, z: zmax } } = + math::fit_psr(/*light_all_mat*/shadow_all_mat/*shadow_view_mat*//* * inverse_visible*/, visible_light_volume.iter().copied(), /*|p| math::Vec3::from(p) / p.w*/math::Vec4::homogenized); // fit_psr(light_all_mat/* * inverse_visible*/, aabb_to_points(visible_bounds).iter().copied(), |p| Vec3::from(p) / p.w); /* let Aabb { min: Vec3 { z: zmin, .. }, max: Vec3 { z: zmax, .. } } = fit_psr(/*light_all_mat*/shadow_all_mat/* * inverse_visible*/, light_volume.iter().copied(), |p| Vec3::from(p) / p.w); @@ -1803,6 +1189,7 @@ impl Scene { let _w_p_arr = directed_proj_mat.cols.iter().map(|e| (e.x, e.y, e.z, e.w)).collect::>(); // println!("mat4 directed_proj_mat = mat4(vec4{:?}, vec4{:?}, vec4{:?}, vec4{:?});", w_p_arr[0], w_p_arr[1], w_p_arr[2], w_p_arr[3]); + let shadow_all_mat: Mat4 = Mat4::from_col_arrays(shadow_all_mat.into_col_arrays()); let _w_p_arr = (directed_proj_mat * shadow_all_mat).cols.iter().map(|e| (e.x, e.y, e.z, e.w)).collect::>(); // println!("mat4 final_mat = mat4(vec4{:?}, vec4{:?}, vec4{:?}, vec4{:?});", w_p_arr[0], w_p_arr[1], w_p_arr[2], w_p_arr[3]); diff --git a/voxygen/src/scene/terrain.rs b/voxygen/src/scene/terrain.rs index ea11244454..0bdd24f324 100644 --- a/voxygen/src/scene/terrain.rs +++ b/voxygen/src/scene/terrain.rs @@ -7,7 +7,7 @@ use crate::{ }, }; -use super::{LodData, SceneData}; +use super::{math, LodData, SceneData}; use common::{ assets, figure::Segment, @@ -2852,7 +2852,7 @@ impl Terrain { max: focus_pos + 0.5f32, }; */ let ray_direction = scene_data.get_sun_dir(); - let collides_with_aabr = |a: Aabr, b: Aabr| { + let collides_with_aabr = |a: math::Aabr, b: math::Aabr| { a.min.partial_cmple(&b.max).reduce_and() && a.max.partial_cmpge(&b.min).reduce_and() }; if ray_direction.z < 0.0 && renderer.render_mode().shadow == render::ShadowMode::Map { @@ -2860,14 +2860,18 @@ impl Terrain { min: visible_bounding_box.min - focus_off, max: visible_bounding_box.max - focus_off, }; - let visible_bounds_fine = Aabb { - min: visible_bounding_box.min.map(f64::from), - max: visible_bounding_box.max.map(f64::from), + let focus_off = math::Vec3::from(focus_off); + let visible_bounds_fine = math::Aabb:: { + min: math::Vec3::from(visible_bounding_box.min.map(f64::from)), + max: math::Vec3::from(visible_bounding_box.max.map(f64::from)), }; - let inv_proj_view = (proj_mat * view_mat/* * Mat4::translation_3d(-focus_off)*/) - .map(f64::from) - .inverted(); - let visible_light_volume = super::calc_focused_light_volume_points( + let inv_proj_view = math::Mat4::from_col_arrays( + (proj_mat * view_mat/* * Mat4::translation_3d(-focus_off)*/).into_col_arrays(), + ) + .map(f64::from) + .inverted(); + let ray_direction = math::Vec3::::from(ray_direction); + let visible_light_volume = math::calc_focused_light_volume_points( inv_proj_view, ray_direction.map(f64::from), visible_bounds_fine, @@ -2876,8 +2880,8 @@ impl Terrain { .map(|v| v.map(|e| e as f32)) .collect::>(); - let cam_pos = Vec3::from(view_mat.inverted() * Vec4::unit_w())/* + focus_off*/; - let view_dir = (focus_pos.map(f32::fract)) - cam_pos; + let cam_pos = math::Vec4::from(view_mat.inverted() * Vec4::unit_w()).xyz()/* + focus_off*/; + /* let view_dir = (focus_pos.map(f32::fract)) - cam_pos; // let new_dir: Vec3 = light_volume/*visible_light_volume*/.iter().map(|p| // p - cam_pos).sum(); let new_dir = view_dir; @@ -2889,27 +2893,27 @@ impl Terrain { new_dir } else { Vec3::from(view_mat * Vec4::from_direction(Vec3::up())).normalized() - }; - let up: Vec3 = { + }; */ + let up: math::Vec3 = { /* (ray_direction) .cross(new_dir) .cross(ray_direction) .normalized() */ - Vec3::up() + math::Vec3::up() }; - let ray_mat = Mat4::look_at_rh( + let ray_mat = math::Mat4::look_at_rh( cam_pos, cam_pos + ray_direction, up, // Vec3::up(), ); // println!("old: {:?} new: {:?}", visible_bounding_box, visible_light_volume); - let visible_bounds = Aabr::from(super::fit_psr( + let visible_bounds = math::Aabr::from(math::fit_psr( ray_mat, /* super::aabb_to_points(visible_bounding_box).iter().copied() */ visible_light_volume.into_iter(), - |p| Vec3::from(p), /* / p.w */ + |p| p, //math::Vec3::from(p), /* / p.w */ )); /* let visible_bounds_old = Aabr::from(super::fit_psr(ray_mat, super::aabb_to_points(visible_bounding_box).iter().copied(), |p| Vec3::from(p) / p.w)); println!("old: {:?} new: {:?}", visible_bounds_old, visible_bounds); */ @@ -2929,12 +2933,16 @@ impl Terrain { chunk_pos.y + chunk_sz, chunk.z_bounds.1, ]; - let chunk_box = Aabb { - min: Vec3::from(chunk_min) - focus_off, - max: Vec3::from(chunk_max) - focus_off, + let chunk_box = math::Aabb { + min: math::Vec3::from(chunk_min) - focus_off, + max: math::Vec3::from(chunk_max) - focus_off, }; - let chunk_from_light = Aabr::from(super::fit_psr(ray_mat, super::aabb_to_points(chunk_box).iter().copied(), |p| Vec3::from(p)/* / p.w*/)); + let chunk_from_light = math::Aabr::from(math::fit_psr(ray_mat, math::aabb_to_points(chunk_box).iter().copied(), |p| p/*math::Vec3::from(p)/* / p.w*/*/)); + /* let chunk_from_light = Aabr { + min: (ray_mat * Vec4::from_point(chunk_box.min)).xy(), + max: (ray_mat * Vec4::from_point(chunk_box.max)).xy(), + }.made_valid(); */ let can_shadow_sun = collides_with_aabr(chunk_from_light, visible_bounds); /* let can_shadow_sun_old = collides_with_aabr(chunk_from_light, visible_bounds_old); if can_shadow_sun != can_shadow_sun_old {