diff options
Diffstat (limited to 'third_party/rust/plane-split/src')
| -rwxr-xr-x | third_party/rust/plane-split/src/bsp.rs | 152 | ||||
| -rwxr-xr-x | third_party/rust/plane-split/src/clip.rs | 161 | ||||
| -rwxr-xr-x | third_party/rust/plane-split/src/lib.rs | 301 | ||||
| -rwxr-xr-x | third_party/rust/plane-split/src/polygon.rs | 597 |
4 files changed, 1211 insertions, 0 deletions
diff --git a/third_party/rust/plane-split/src/bsp.rs b/third_party/rust/plane-split/src/bsp.rs new file mode 100755 index 0000000000..b9bb3109d1 --- /dev/null +++ b/third_party/rust/plane-split/src/bsp.rs @@ -0,0 +1,152 @@ +use crate::{is_zero, Intersection, Plane, Polygon, Splitter}; + +use binary_space_partition::{BspNode, Plane as BspPlane, PlaneCut}; +use euclid::{approxeq::ApproxEq, Point3D, Vector3D}; +use num_traits::{Float, One, Zero}; + +use std::{fmt, iter, ops}; + +impl<T, U, A> BspPlane for Polygon<T, U, A> +where + T: Copy + + fmt::Debug + + ApproxEq<T> + + ops::Sub<T, Output = T> + + ops::Add<T, Output = T> + + ops::Mul<T, Output = T> + + ops::Div<T, Output = T> + + Zero + + Float, + U: fmt::Debug, + A: Copy + fmt::Debug, +{ + fn cut(&self, mut poly: Self) -> PlaneCut<Self> { + log::debug!("\tCutting anchor {:?} by {:?}", poly.anchor, self.anchor); + log::trace!("\t\tbase {:?}", self.plane); + + //Note: we treat `self` as a plane, and `poly` as a concrete polygon here + let (intersection, dist) = match self.plane.intersect(&poly.plane) { + None => { + let ndot = self.plane.normal.dot(poly.plane.normal); + log::debug!("\t\tNormals are aligned with {:?}", ndot); + let dist = self.plane.offset - ndot * poly.plane.offset; + (Intersection::Coplanar, dist) + } + Some(_) if self.plane.are_outside(&poly.points) => { + //Note: we can't start with `are_outside` because it's subject to FP precision + let dist = self.plane.signed_distance_sum_to(&poly); + (Intersection::Outside, dist) + } + Some(line) => { + //Note: distance isn't relevant here + (Intersection::Inside(line), T::zero()) + } + }; + + match intersection { + //Note: we deliberately make the comparison wider than just with T::epsilon(). + // This is done to avoid mistakenly ordering items that should be on the same + // plane but end up slightly different due to the floating point precision. + Intersection::Coplanar if is_zero(dist) => { + log::debug!("\t\tCoplanar at {:?}", dist); + PlaneCut::Sibling(poly) + } + Intersection::Coplanar | Intersection::Outside => { + log::debug!("\t\tOutside at {:?}", dist); + if dist > T::zero() { + PlaneCut::Cut { + front: vec![poly], + back: vec![], + } + } else { + PlaneCut::Cut { + front: vec![], + back: vec![poly], + } + } + } + Intersection::Inside(line) => { + log::debug!("\t\tCut across {:?}", line); + let (res_add1, res_add2) = poly.split_with_normal(&line, &self.plane.normal); + let mut front = Vec::new(); + let mut back = Vec::new(); + + for sub in iter::once(poly) + .chain(res_add1) + .chain(res_add2) + .filter(|p| !p.is_empty()) + { + let dist = self.plane.signed_distance_sum_to(&sub); + if dist > T::zero() { + log::trace!("\t\t\tdist {:?} -> front: {:?}", dist, sub); + front.push(sub) + } else { + log::trace!("\t\t\tdist {:?} -> back: {:?}", dist, sub); + back.push(sub) + } + } + + PlaneCut::Cut { front, back } + } + } + } + + fn is_aligned(&self, other: &Self) -> bool { + self.plane.normal.dot(other.plane.normal) > T::zero() + } +} + +/// Binary Space Partitioning splitter, uses a BSP tree. +pub struct BspSplitter<T, U, A> { + tree: BspNode<Polygon<T, U, A>>, + result: Vec<Polygon<T, U, A>>, +} + +impl<T, U, A> BspSplitter<T, U, A> { + /// Create a new BSP splitter. + pub fn new() -> Self { + BspSplitter { + tree: BspNode::new(), + result: Vec::new(), + } + } +} + +impl<T, U, A> Splitter<T, U, A> for BspSplitter<T, U, A> +where + T: Copy + + fmt::Debug + + ApproxEq<T> + + ops::Sub<T, Output = T> + + ops::Add<T, Output = T> + + ops::Mul<T, Output = T> + + ops::Div<T, Output = T> + + Zero + + One + + Float, + U: fmt::Debug, + A: Copy + fmt::Debug + Default, +{ + fn reset(&mut self) { + self.tree = BspNode::new(); + } + + fn add(&mut self, poly: Polygon<T, U, A>) { + self.tree.insert(poly); + } + + fn sort(&mut self, view: Vector3D<T, U>) -> &[Polygon<T, U, A>] { + //debug!("\t\ttree before sorting {:?}", self.tree); + let poly = Polygon { + points: [Point3D::origin(); 4], + plane: Plane { + normal: -view, //Note: BSP `order()` is back to front + offset: T::zero(), + }, + anchor: A::default(), + }; + self.result.clear(); + self.tree.order(&poly, &mut self.result); + &self.result + } +} diff --git a/third_party/rust/plane-split/src/clip.rs b/third_party/rust/plane-split/src/clip.rs new file mode 100755 index 0000000000..ee240b2f2f --- /dev/null +++ b/third_party/rust/plane-split/src/clip.rs @@ -0,0 +1,161 @@ +use crate::{Intersection, NegativeHemisphereError, Plane, Polygon}; + +use euclid::{approxeq::ApproxEq, Rect, Scale, Transform3D, Trig, Vector3D}; +use num_traits::{Float, One, Zero}; + +use std::{fmt, iter, mem, ops}; + +/// A helper object to clip polygons by a number of planes. +#[derive(Debug)] +pub struct Clipper<T, U, A> { + clips: Vec<Plane<T, U>>, + results: Vec<Polygon<T, U, A>>, + temp: Vec<Polygon<T, U, A>>, +} + +impl< + T: Copy + + fmt::Debug + + ApproxEq<T> + + ops::Sub<T, Output = T> + + ops::Add<T, Output = T> + + ops::Mul<T, Output = T> + + ops::Div<T, Output = T> + + Zero + + One + + Float, + U: fmt::Debug, + A: Copy + fmt::Debug, + > Clipper<T, U, A> +{ + /// Create a new clipper object. + pub fn new() -> Self { + Clipper { + clips: Vec::new(), + results: Vec::new(), + temp: Vec::new(), + } + } + + /// Reset the clipper internals, but preserve the allocation. + pub fn reset(&mut self) { + self.clips.clear(); + } + + /// Extract the clipping planes that define the frustum for a given transformation. + pub fn frustum_planes<V>( + t: &Transform3D<T, U, V>, + bounds: Option<Rect<T, V>>, + ) -> Result<impl Iterator<Item = Plane<T, U>>, NegativeHemisphereError> { + let mw = Vector3D::new(t.m14, t.m24, t.m34); + let plane_positive = Plane::from_unnormalized(mw, t.m44)?; + + let bounds_iter_maybe = match bounds { + Some(bounds) => { + let mx = Vector3D::new(t.m11, t.m21, t.m31); + let left = bounds.origin.x; + let plane_left = + Plane::from_unnormalized(mx - mw * Scale::new(left), t.m41 - t.m44 * left)?; + let right = bounds.origin.x + bounds.size.width; + let plane_right = + Plane::from_unnormalized(mw * Scale::new(right) - mx, t.m44 * right - t.m41)?; + + let my = Vector3D::new(t.m12, t.m22, t.m32); + let top = bounds.origin.y; + let plane_top = + Plane::from_unnormalized(my - mw * Scale::new(top), t.m42 - t.m44 * top)?; + let bottom = bounds.origin.y + bounds.size.height; + let plane_bottom = + Plane::from_unnormalized(mw * Scale::new(bottom) - my, t.m44 * bottom - t.m42)?; + + Some( + plane_left + .into_iter() + .chain(plane_right) + .chain(plane_top) + .chain(plane_bottom), + ) + } + None => None, + }; + + Ok(bounds_iter_maybe + .into_iter() + .flat_map(|pi| pi) + .chain(plane_positive)) + } + + /// Add a clipping plane to the list. The plane will clip everything behind it, + /// where the direction is set by the plane normal. + pub fn add(&mut self, plane: Plane<T, U>) { + self.clips.push(plane); + } + + /// Clip specified polygon by the contained planes, return the fragmented polygons. + pub fn clip(&mut self, polygon: Polygon<T, U, A>) -> &[Polygon<T, U, A>] { + log::debug!("\tClipping {:?}", polygon); + self.results.clear(); + self.results.push(polygon); + + for clip in &self.clips { + self.temp.clear(); + mem::swap(&mut self.results, &mut self.temp); + + for mut poly in self.temp.drain(..) { + let dist = match poly.intersect_plane(clip) { + Intersection::Inside(line) => { + let (res1, res2) = poly.split_with_normal(&line, &clip.normal); + self.results.extend( + iter::once(poly) + .chain(res1) + .chain(res2) + .filter(|p| clip.signed_distance_sum_to(p) > T::zero()), + ); + continue; + } + Intersection::Coplanar => { + let ndot = poly.plane.normal.dot(clip.normal); + clip.offset - ndot * poly.plane.offset + } + Intersection::Outside => clip.signed_distance_sum_to(&poly), + }; + + if dist > T::zero() { + self.results.push(poly); + } + } + } + + &self.results + } + + /// Clip the primitive with the frustum of the specified transformation, + /// returning a sequence of polygons in the transformed space. + /// Returns None if the transformation can't be frustum clipped. + pub fn clip_transformed<'a, V>( + &'a mut self, + polygon: Polygon<T, U, A>, + transform: &'a Transform3D<T, U, V>, + bounds: Option<Rect<T, V>>, + ) -> Result<impl 'a + Iterator<Item = Polygon<T, V, A>>, NegativeHemisphereError> + where + T: Trig, + V: 'a + fmt::Debug, + { + let planes = Self::frustum_planes(transform, bounds)?; + + let old_count = self.clips.len(); + self.clips.extend(planes); + self.clip(polygon); + // remove the frustum planes + while self.clips.len() > old_count { + self.clips.pop(); + } + + let polys = self + .results + .drain(..) + .flat_map(move |poly| poly.transform(transform)); + Ok(polys) + } +} diff --git a/third_party/rust/plane-split/src/lib.rs b/third_party/rust/plane-split/src/lib.rs new file mode 100755 index 0000000000..991766a53d --- /dev/null +++ b/third_party/rust/plane-split/src/lib.rs @@ -0,0 +1,301 @@ +/*! +Plane splitting. + +Uses [euclid](https://crates.io/crates/euclid) for the math basis. +Introduces new geometrical primitives and associated logic. + +Automatically splits a given set of 4-point polygons into sub-polygons +that don't intersect each other. This is useful for WebRender, to sort +the resulting sub-polygons by depth and avoid transparency blending issues. +*/ +#![warn(missing_docs)] + +mod bsp; +mod clip; +mod polygon; + +use euclid::{approxeq::ApproxEq, Point3D, Scale, Vector3D}; +use num_traits::{Float, One, Zero}; + +use std::ops; + +pub use self::bsp::BspSplitter; +pub use self::clip::Clipper; +pub use self::polygon::{Intersection, LineProjection, Polygon}; + +fn is_zero<T>(value: T) -> bool +where + T: Copy + Zero + ApproxEq<T> + ops::Mul<T, Output = T>, +{ + //HACK: this is rough, but the original Epsilon is too strict + (value * value).approx_eq(&T::zero()) +} + +fn is_zero_vec<T, U>(vec: Vector3D<T, U>) -> bool +where + T: Copy + + Zero + + ApproxEq<T> + + ops::Add<T, Output = T> + + ops::Sub<T, Output = T> + + ops::Mul<T, Output = T>, +{ + vec.dot(vec).approx_eq(&T::zero()) +} + +/// A generic line. +#[derive(Debug)] +pub struct Line<T, U> { + /// Arbitrary point on the line. + pub origin: Point3D<T, U>, + /// Normalized direction of the line. + pub dir: Vector3D<T, U>, +} + +impl<T, U> Line<T, U> +where + T: Copy + + One + + Zero + + ApproxEq<T> + + ops::Add<T, Output = T> + + ops::Sub<T, Output = T> + + ops::Mul<T, Output = T>, +{ + /// Check if the line has consistent parameters. + pub fn is_valid(&self) -> bool { + is_zero(self.dir.dot(self.dir) - T::one()) + } + /// Check if two lines match each other. + pub fn matches(&self, other: &Self) -> bool { + let diff = self.origin - other.origin; + is_zero_vec(self.dir.cross(other.dir)) && is_zero_vec(self.dir.cross(diff)) + } + + /// Intersect an edge given by the end points. + /// Returns the fraction of the edge where the intersection occurs. + fn intersect_edge(&self, edge: ops::Range<Point3D<T, U>>) -> Option<T> + where + T: ops::Div<T, Output = T>, + { + let edge_vec = edge.end - edge.start; + let origin_vec = self.origin - edge.start; + // edge.start + edge_vec * t = r + k * d + // (edge.start, d) + t * (edge_vec, d) - (r, d) = k + // edge.start + t * edge_vec = r + t * (edge_vec, d) * d + (start-r, d) * d + // t * (edge_vec - (edge_vec, d)*d) = origin_vec - (origin_vec, d) * d + let pr = origin_vec - self.dir * self.dir.dot(origin_vec); + let pb = edge_vec - self.dir * self.dir.dot(edge_vec); + let denom = pb.dot(pb); + if denom.approx_eq(&T::zero()) { + None + } else { + Some(pr.dot(pb) / denom) + } + } +} + +/// An infinite plane in 3D space, defined by equation: +/// dot(v, normal) + offset = 0 +/// When used for plane splitting, it's defining a hemisphere +/// with equation "dot(v, normal) + offset > 0". +#[derive(Debug, PartialEq)] +pub struct Plane<T, U> { + /// Normalized vector perpendicular to the plane. + pub normal: Vector3D<T, U>, + /// Constant offset from the normal plane, specified in the + /// direction opposite to the normal. + pub offset: T, +} + +impl<T: Clone, U> Clone for Plane<T, U> { + fn clone(&self) -> Self { + Plane { + normal: self.normal.clone(), + offset: self.offset.clone(), + } + } +} + +/// An error returned when everything would end up projected +/// to the negative hemisphere (W <= 0.0); +#[derive(Clone, Debug, Hash, PartialEq, PartialOrd)] +pub struct NegativeHemisphereError; + +impl< + T: Copy + + Zero + + One + + Float + + ApproxEq<T> + + ops::Sub<T, Output = T> + + ops::Add<T, Output = T> + + ops::Mul<T, Output = T> + + ops::Div<T, Output = T>, + U, + > Plane<T, U> +{ + /// Construct a new plane from unnormalized equation. + pub fn from_unnormalized( + normal: Vector3D<T, U>, + offset: T, + ) -> Result<Option<Self>, NegativeHemisphereError> { + let square_len = normal.square_length(); + if square_len < T::approx_epsilon() * T::approx_epsilon() { + if offset > T::zero() { + Ok(None) + } else { + Err(NegativeHemisphereError) + } + } else { + let kf = T::one() / square_len.sqrt(); + Ok(Some(Plane { + normal: normal * Scale::new(kf), + offset: offset * kf, + })) + } + } + + /// Check if this plane contains another one. + pub fn contains(&self, other: &Self) -> bool { + //TODO: actually check for inside/outside + self.normal == other.normal && self.offset == other.offset + } + + /// Return the signed distance from this plane to a point. + /// The distance is negative if the point is on the other side of the plane + /// from the direction of the normal. + pub fn signed_distance_to(&self, point: &Point3D<T, U>) -> T { + point.to_vector().dot(self.normal) + self.offset + } + + /// Compute the distance across the line to the plane plane, + /// starting from the line origin. + pub fn distance_to_line(&self, line: &Line<T, U>) -> T + where + T: ops::Neg<Output = T>, + { + self.signed_distance_to(&line.origin) / -self.normal.dot(line.dir) + } + + /// Compute the sum of signed distances to each of the points + /// of another plane. Useful to know the relation of a plane that + /// is a product of a split, and we know it doesn't intersect `self`. + pub fn signed_distance_sum_to<A>(&self, poly: &Polygon<T, U, A>) -> T { + poly.points + .iter() + .fold(T::zero(), |u, p| u + self.signed_distance_to(p)) + } + + /// Check if a convex shape defined by a set of points is completely + /// outside of this plane. Merely touching the surface is not + /// considered an intersection. + pub fn are_outside(&self, points: &[Point3D<T, U>]) -> bool { + let d0 = self.signed_distance_to(&points[0]); + points[1..] + .iter() + .all(|p| self.signed_distance_to(p) * d0 > T::zero()) + } + + //TODO(breaking): turn this into Result<Line, DotProduct> + /// Compute the line of intersection with another plane. + pub fn intersect(&self, other: &Self) -> Option<Line<T, U>> { + // compute any point on the intersection between planes + // (n1, v) + d1 = 0 + // (n2, v) + d2 = 0 + // v = a*n1/w + b*n2/w; w = (n1, n2) + // v = (d2*w - d1) / (1 - w*w) * n1 - (d2 - d1*w) / (1 - w*w) * n2 + let w = self.normal.dot(other.normal); + let divisor = T::one() - w * w; + if divisor < T::approx_epsilon() * T::approx_epsilon() { + return None; + } + let origin = Point3D::origin() + self.normal * ((other.offset * w - self.offset) / divisor) + - other.normal * ((other.offset - self.offset * w) / divisor); + + let cross_dir = self.normal.cross(other.normal); + // note: the cross product isn't too close to zero + // due to the previous check + + Some(Line { + origin, + dir: cross_dir.normalize(), + }) + } +} + +/// Generic plane splitter interface +pub trait Splitter<T, U, A> { + /// Reset the splitter results. + fn reset(&mut self); + + /// Add a new polygon and return a slice of the subdivisions + /// that avoid collision with any of the previously added polygons. + fn add(&mut self, polygon: Polygon<T, U, A>); + + /// Sort the produced polygon set by the ascending distance across + /// the specified view vector. Return the sorted slice. + fn sort(&mut self, view: Vector3D<T, U>) -> &[Polygon<T, U, A>]; + + /// Process a set of polygons at once. + fn solve(&mut self, input: &[Polygon<T, U, A>], view: Vector3D<T, U>) -> &[Polygon<T, U, A>] + where + T: Clone, + U: Clone, + A: Copy, + { + self.reset(); + for p in input { + self.add(p.clone()); + } + self.sort(view) + } +} + +/// Helper method used for benchmarks and tests. +/// Constructs a 3D grid of polygons. +#[doc(hidden)] +pub fn make_grid(count: usize) -> Vec<Polygon<f32, (), usize>> { + let mut polys: Vec<Polygon<f32, (), usize>> = Vec::with_capacity(count * 3); + let len = count as f32; + polys.extend((0..count).map(|i| Polygon { + points: [ + Point3D::new(0.0, i as f32, 0.0), + Point3D::new(len, i as f32, 0.0), + Point3D::new(len, i as f32, len), + Point3D::new(0.0, i as f32, len), + ], + plane: Plane { + normal: Vector3D::new(0.0, 1.0, 0.0), + offset: -(i as f32), + }, + anchor: 0, + })); + polys.extend((0..count).map(|i| Polygon { + points: [ + Point3D::new(i as f32, 0.0, 0.0), + Point3D::new(i as f32, len, 0.0), + Point3D::new(i as f32, len, len), + Point3D::new(i as f32, 0.0, len), + ], + plane: Plane { + normal: Vector3D::new(1.0, 0.0, 0.0), + offset: -(i as f32), + }, + anchor: 0, + })); + polys.extend((0..count).map(|i| Polygon { + points: [ + Point3D::new(0.0, 0.0, i as f32), + Point3D::new(len, 0.0, i as f32), + Point3D::new(len, len, i as f32), + Point3D::new(0.0, len, i as f32), + ], + plane: Plane { + normal: Vector3D::new(0.0, 0.0, 1.0), + offset: -(i as f32), + }, + anchor: 0, + })); + polys +} diff --git a/third_party/rust/plane-split/src/polygon.rs b/third_party/rust/plane-split/src/polygon.rs new file mode 100755 index 0000000000..a7a0b7ade3 --- /dev/null +++ b/third_party/rust/plane-split/src/polygon.rs @@ -0,0 +1,597 @@ +use crate::{is_zero, Line, Plane}; + +use euclid::{approxeq::ApproxEq, default::Point2D, Point3D, Rect, Transform3D, Trig, Vector3D}; +use num_traits::{Float, One, Zero}; + +use std::{fmt, iter, mem, ops}; + +/// The projection of a `Polygon` on a line. +pub struct LineProjection<T> { + /// Projected value of each point in the polygon. + pub markers: [T; 4], +} + +impl<T> LineProjection<T> +where + T: Copy + PartialOrd + ops::Sub<T, Output = T> + ops::Add<T, Output = T>, +{ + /// Get the min/max of the line projection markers. + pub fn get_bounds(&self) -> (T, T) { + let (mut a, mut b, mut c, mut d) = ( + self.markers[0], + self.markers[1], + self.markers[2], + self.markers[3], + ); + // bitonic sort of 4 elements + // we could not just use `min/max` since they require `Ord` bound + //TODO: make it nicer + if a > c { + mem::swap(&mut a, &mut c); + } + if b > d { + mem::swap(&mut b, &mut d); + } + if a > b { + mem::swap(&mut a, &mut b); + } + if c > d { + mem::swap(&mut c, &mut d); + } + if b > c { + mem::swap(&mut b, &mut c); + } + debug_assert!(a <= b && b <= c && c <= d); + (a, d) + } + + /// Check intersection with another line projection. + pub fn intersect(&self, other: &Self) -> bool { + // compute the bounds of both line projections + let span = self.get_bounds(); + let other_span = other.get_bounds(); + // compute the total footprint + let left = if span.0 < other_span.0 { + span.0 + } else { + other_span.0 + }; + let right = if span.1 > other_span.1 { + span.1 + } else { + other_span.1 + }; + // they intersect if the footprint is smaller than the sum + right - left < span.1 - span.0 + other_span.1 - other_span.0 + } +} + +/// Polygon intersection results. +pub enum Intersection<T> { + /// Polygons are coplanar, including the case of being on the same plane. + Coplanar, + /// Polygon planes are intersecting, but polygons are not. + Outside, + /// Polygons are actually intersecting. + Inside(T), +} + +impl<T> Intersection<T> { + /// Return true if the intersection is completely outside. + pub fn is_outside(&self) -> bool { + match *self { + Intersection::Outside => true, + _ => false, + } + } + /// Return true if the intersection cuts the source polygon. + pub fn is_inside(&self) -> bool { + match *self { + Intersection::Inside(_) => true, + _ => false, + } + } +} + +/// A convex polygon with 4 points lying on a plane. +#[derive(Debug, PartialEq)] +pub struct Polygon<T, U, A> { + /// Points making the polygon. + pub points: [Point3D<T, U>; 4], + /// A plane describing polygon orientation. + pub plane: Plane<T, U>, + /// A simple anchoring index to allow association of the + /// produced split polygons with the original one. + pub anchor: A, +} + +impl<T: Clone, U, A: Copy> Clone for Polygon<T, U, A> { + fn clone(&self) -> Self { + Polygon { + points: [ + self.points[0].clone(), + self.points[1].clone(), + self.points[2].clone(), + self.points[3].clone(), + ], + plane: self.plane.clone(), + anchor: self.anchor, + } + } +} + +impl<T, U, A> Polygon<T, U, A> +where + T: Copy + + fmt::Debug + + ApproxEq<T> + + ops::Sub<T, Output = T> + + ops::Add<T, Output = T> + + ops::Mul<T, Output = T> + + ops::Div<T, Output = T> + + Zero + + One + + Float, + U: fmt::Debug, + A: Copy, +{ + /// Construct a polygon from points that are already transformed. + /// Return None if the polygon doesn't contain any space. + pub fn from_points(points: [Point3D<T, U>; 4], anchor: A) -> Option<Self> { + let edge1 = points[1] - points[0]; + let edge2 = points[2] - points[0]; + let edge3 = points[3] - points[0]; + let edge4 = points[3] - points[1]; + + if edge2.square_length() < T::epsilon() || edge4.square_length() < T::epsilon() { + return None; + } + + // one of them can be zero for redundant polygons produced by plane splitting + //Note: this would be nicer if we used triangles instead of quads in the first place... + // see https://github.com/servo/plane-split/issues/17 + let normal_rough1 = edge1.cross(edge2); + let normal_rough2 = edge2.cross(edge3); + let square_length1 = normal_rough1.square_length(); + let square_length2 = normal_rough2.square_length(); + let normal = if square_length1 > square_length2 { + normal_rough1 / square_length1.sqrt() + } else { + normal_rough2 / square_length2.sqrt() + }; + + let offset = -points[0].to_vector().dot(normal); + + Some(Polygon { + points, + plane: Plane { normal, offset }, + anchor, + }) + } + + /// Construct a polygon from a non-transformed rectangle. + pub fn from_rect(rect: Rect<T, U>, anchor: A) -> Self { + let min = rect.min(); + let max = rect.max(); + let _0 = T::zero(); + Polygon { + points: [ + min.to_3d(), + Point3D::new(max.x, min.y, _0), + max.to_3d(), + Point3D::new(min.x, max.y, _0), + ], + plane: Plane { + normal: Vector3D::new(T::zero(), T::zero(), T::one()), + offset: T::zero(), + }, + anchor, + } + } + + /// Construct a polygon from a rectangle with 3D transform. + pub fn from_transformed_rect<V>( + rect: Rect<T, V>, + transform: Transform3D<T, V, U>, + anchor: A, + ) -> Option<Self> + where + T: Trig + ops::Neg<Output = T>, + { + let min = rect.min(); + let max = rect.max(); + let _0 = T::zero(); + let points = [ + transform.transform_point3d(min.to_3d())?, + transform.transform_point3d(Point3D::new(max.x, min.y, _0))?, + transform.transform_point3d(max.to_3d())?, + transform.transform_point3d(Point3D::new(min.x, max.y, _0))?, + ]; + Self::from_points(points, anchor) + } + + /// Construct a polygon from a rectangle with an invertible 3D transform. + pub fn from_transformed_rect_with_inverse<V>( + rect: Rect<T, V>, + transform: &Transform3D<T, V, U>, + inv_transform: &Transform3D<T, U, V>, + anchor: A, + ) -> Option<Self> + where + T: Trig + ops::Neg<Output = T>, + { + let min = rect.min(); + let max = rect.max(); + let _0 = T::zero(); + let points = [ + transform.transform_point3d(min.to_3d())?, + transform.transform_point3d(Point3D::new(max.x, min.y, _0))?, + transform.transform_point3d(max.to_3d())?, + transform.transform_point3d(Point3D::new(min.x, max.y, _0))?, + ]; + + // Compute the normal directly from the transformation. This guarantees consistent polygons + // generated from various local rectanges on the same geometry plane. + let normal_raw = Vector3D::new(inv_transform.m13, inv_transform.m23, inv_transform.m33); + let normal_sql = normal_raw.square_length(); + if normal_sql.approx_eq(&T::zero()) || transform.m44.approx_eq(&T::zero()) { + None + } else { + let normal = normal_raw / normal_sql.sqrt(); + let offset = -Vector3D::new(transform.m41, transform.m42, transform.m43).dot(normal) + / transform.m44; + + Some(Polygon { + points, + plane: Plane { normal, offset }, + anchor, + }) + } + } + + /// Bring a point into the local coordinate space, returning + /// the 2D normalized coordinates. + pub fn untransform_point(&self, point: Point3D<T, U>) -> Point2D<T> { + //debug_assert!(self.contains(point)); + // get axises and target vector + let a = self.points[1] - self.points[0]; + let b = self.points[3] - self.points[0]; + let c = point - self.points[0]; + // get pair-wise dot products + let a2 = a.dot(a); + let ab = a.dot(b); + let b2 = b.dot(b); + let ca = c.dot(a); + let cb = c.dot(b); + // compute the final coordinates + let denom = ab * ab - a2 * b2; + let x = ab * cb - b2 * ca; + let y = ab * ca - a2 * cb; + Point2D::new(x, y) / denom + } + + /// Transform a polygon by an affine transform (preserving straight lines). + pub fn transform<V>(&self, transform: &Transform3D<T, U, V>) -> Option<Polygon<T, V, A>> + where + T: Trig, + V: fmt::Debug, + { + let mut points = [Point3D::origin(); 4]; + for (out, point) in points.iter_mut().zip(self.points.iter()) { + let mut homo = transform.transform_point3d_homogeneous(*point); + homo.w = homo.w.max(T::approx_epsilon()); + *out = homo.to_point3d()?; + } + + //Note: this code path could be more efficient if we had inverse-transpose + //let n4 = transform.transform_point4d(&Point4D::new(T::zero(), T::zero(), T::one(), T::zero())); + //let normal = Point3D::new(n4.x, n4.y, n4.z); + Polygon::from_points(points, self.anchor) + } + + /// Check if all the points are indeed placed on the plane defined by + /// the normal and offset, and the winding order is consistent. + pub fn is_valid(&self) -> bool { + let is_planar = self + .points + .iter() + .all(|p| is_zero(self.plane.signed_distance_to(p))); + let edges = [ + self.points[1] - self.points[0], + self.points[2] - self.points[1], + self.points[3] - self.points[2], + self.points[0] - self.points[3], + ]; + let anchor = edges[3].cross(edges[0]); + let is_winding = edges + .iter() + .zip(edges[1..].iter()) + .all(|(a, &b)| a.cross(b).dot(anchor) >= T::zero()); + is_planar && is_winding + } + + /// Check if the polygon doesn't contain any space. This may happen + /// after a sequence of splits, and such polygons should be discarded. + pub fn is_empty(&self) -> bool { + (self.points[0] - self.points[2]).square_length() < T::epsilon() + || (self.points[1] - self.points[3]).square_length() < T::epsilon() + } + + /// Check if this polygon contains another one. + pub fn contains(&self, other: &Self) -> bool { + //TODO: actually check for inside/outside + self.plane.contains(&other.plane) + } + + /// Project this polygon onto a 3D vector, returning a line projection. + /// Note: we can think of it as a projection to a ray placed at the origin. + pub fn project_on(&self, vector: &Vector3D<T, U>) -> LineProjection<T> { + LineProjection { + markers: [ + vector.dot(self.points[0].to_vector()), + vector.dot(self.points[1].to_vector()), + vector.dot(self.points[2].to_vector()), + vector.dot(self.points[3].to_vector()), + ], + } + } + + /// Compute the line of intersection with an infinite plane. + pub fn intersect_plane(&self, other: &Plane<T, U>) -> Intersection<Line<T, U>> { + if other.are_outside(&self.points) { + log::debug!("\t\tOutside of the plane"); + return Intersection::Outside; + } + match self.plane.intersect(&other) { + Some(line) => Intersection::Inside(line), + None => { + log::debug!("\t\tCoplanar"); + Intersection::Coplanar + } + } + } + + /// Compute the line of intersection with another polygon. + pub fn intersect(&self, other: &Self) -> Intersection<Line<T, U>> { + if self.plane.are_outside(&other.points) || other.plane.are_outside(&self.points) { + log::debug!("\t\tOne is completely outside of the other"); + return Intersection::Outside; + } + match self.plane.intersect(&other.plane) { + Some(line) => { + let self_proj = self.project_on(&line.dir); + let other_proj = other.project_on(&line.dir); + if self_proj.intersect(&other_proj) { + Intersection::Inside(line) + } else { + // projections on the line don't intersect + log::debug!("\t\tProjection is outside"); + Intersection::Outside + } + } + None => { + log::debug!("\t\tCoplanar"); + Intersection::Coplanar + } + } + } + + fn split_impl( + &mut self, + first: (usize, Point3D<T, U>), + second: (usize, Point3D<T, U>), + ) -> (Option<Self>, Option<Self>) { + //TODO: can be optimized for when the polygon has a redundant 4th vertex + //TODO: can be simplified greatly if only working with triangles + log::debug!("\t\tReached complex case [{}, {}]", first.0, second.0); + let base = first.0; + assert!(base < self.points.len()); + match second.0 - first.0 { + 1 => { + // rect between the cut at the diagonal + let other1 = Polygon { + points: [ + first.1, + second.1, + self.points[(base + 2) & 3], + self.points[base], + ], + ..self.clone() + }; + // triangle on the near side of the diagonal + let other2 = Polygon { + points: [ + self.points[(base + 2) & 3], + self.points[(base + 3) & 3], + self.points[base], + self.points[base], + ], + ..self.clone() + }; + // triangle being cut out + self.points = [first.1, self.points[(base + 1) & 3], second.1, second.1]; + (Some(other1), Some(other2)) + } + 2 => { + // rect on the far side + let other = Polygon { + points: [ + first.1, + self.points[(base + 1) & 3], + self.points[(base + 2) & 3], + second.1, + ], + ..self.clone() + }; + // rect on the near side + self.points = [ + first.1, + second.1, + self.points[(base + 3) & 3], + self.points[base], + ]; + (Some(other), None) + } + 3 => { + // rect between the cut at the diagonal + let other1 = Polygon { + points: [ + first.1, + self.points[(base + 1) & 3], + self.points[(base + 3) & 3], + second.1, + ], + ..self.clone() + }; + // triangle on the far side of the diagonal + let other2 = Polygon { + points: [ + self.points[(base + 1) & 3], + self.points[(base + 2) & 3], + self.points[(base + 3) & 3], + self.points[(base + 3) & 3], + ], + ..self.clone() + }; + // triangle being cut out + self.points = [first.1, second.1, self.points[base], self.points[base]]; + (Some(other1), Some(other2)) + } + _ => panic!("Unexpected indices {} {}", first.0, second.0), + } + } + + /// Split the polygon along the specified `Line`. + /// Will do nothing if the line doesn't belong to the polygon plane. + #[deprecated(note = "Use split_with_normal instead")] + pub fn split(&mut self, line: &Line<T, U>) -> (Option<Self>, Option<Self>) { + log::debug!("\tSplitting"); + // check if the cut is within the polygon plane first + if !is_zero(self.plane.normal.dot(line.dir)) + || !is_zero(self.plane.signed_distance_to(&line.origin)) + { + log::debug!( + "\t\tDoes not belong to the plane, normal dot={:?}, origin distance={:?}", + self.plane.normal.dot(line.dir), + self.plane.signed_distance_to(&line.origin) + ); + return (None, None); + } + // compute the intersection points for each edge + let mut cuts = [None; 4]; + for ((&b, &a), cut) in self + .points + .iter() + .cycle() + .skip(1) + .zip(self.points.iter()) + .zip(cuts.iter_mut()) + { + if let Some(t) = line.intersect_edge(a..b) { + if t >= T::zero() && t < T::one() { + *cut = Some(a + (b - a) * t); + } + } + } + + let first = match cuts.iter().position(|c| c.is_some()) { + Some(pos) => pos, + None => return (None, None), + }; + let second = match cuts[first + 1..].iter().position(|c| c.is_some()) { + Some(pos) => first + 1 + pos, + None => return (None, None), + }; + self.split_impl( + (first, cuts[first].unwrap()), + (second, cuts[second].unwrap()), + ) + } + + /// Split the polygon along the specified `Line`, with a normal to the split line provided. + /// This is useful when called by the plane splitter, since the other plane's normal + /// forms the side direction here, and figuring out the actual line of split isn't needed. + /// Will do nothing if the line doesn't belong to the polygon plane. + pub fn split_with_normal( + &mut self, + line: &Line<T, U>, + normal: &Vector3D<T, U>, + ) -> (Option<Self>, Option<Self>) { + log::debug!("\tSplitting with normal"); + // figure out which side of the split does each point belong to + let mut sides = [T::zero(); 4]; + let (mut cut_positive, mut cut_negative) = (None, None); + for (side, point) in sides.iter_mut().zip(&self.points) { + *side = normal.dot(*point - line.origin); + } + // compute the edge intersection points + for (i, ((&side1, point1), (&side0, point0))) in sides[1..] + .iter() + .chain(iter::once(&sides[0])) + .zip(self.points[1..].iter().chain(iter::once(&self.points[0]))) + .zip(sides.iter().zip(&self.points)) + .enumerate() + { + // figure out if an edge between 0 and 1 needs to be cut + let cut = if side0 < T::zero() && side1 >= T::zero() { + &mut cut_positive + } else if side0 > T::zero() && side1 <= T::zero() { + &mut cut_negative + } else { + continue; + }; + // compute the cut point by weighting the opposite distances + // + // Note: this algorithm is designed to not favor one end of the edge over the other. + // The previous approach of calling `intersect_edge` sometimes ended up with "t" ever + // slightly outside of [0, 1] range, since it was computing it relative to the first point only. + // + // Given that we are intersecting two straight lines, the triangles on both + // sides of intersection are alike, so distances along the [point0, point1] line + // are proportional to the side vector lengths we just computed: (side0, side1). + let point = + (*point0 * side1.abs() + point1.to_vector() * side0.abs()) / (side0 - side1).abs(); + if cut.is_some() { + // We don't expect that the direction changes more than once, unless + // the polygon is close to redundant, and we hit precision issues when + // computing the sides. + log::warn!("Splitting failed due to precision issues: {:?}", sides); + break; + } + *cut = Some((i, point)); + } + // form new polygons + if let (Some(first), Some(mut second)) = (cut_positive, cut_negative) { + if second.0 < first.0 { + second.0 += 4; + } + self.split_impl(first, second) + } else { + (None, None) + } + } +} + +#[test] +fn test_split_precision() { + // regression test for https://bugzilla.mozilla.org/show_bug.cgi?id=1678454 + let mut polygon = Polygon::<_, (), ()> { + points: [ + Point3D::new(300.0102, 150.00958, 0.0), + Point3D::new(606.0, 306.0, 0.0), + Point3D::new(300.21954, 150.11946, 0.0), + Point3D::new(300.08844, 150.05064, 0.0), + ], + plane: Plane { + normal: Vector3D::zero(), + offset: 0.0, + }, + anchor: (), + }; + let line = Line { + origin: Point3D::new(3.0690663, -5.8472385, 0.0), + dir: Vector3D::new(0.8854436, 0.46474677, -0.0), + }; + let normal = Vector3D::new(0.46474662, -0.8854434, -0.0006389789); + polygon.split_with_normal(&line, &normal); +} |