From 26a029d407be480d791972afb5975cf62c9360a6 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Fri, 19 Apr 2024 02:47:55 +0200 Subject: Adding upstream version 124.0.1. Signed-off-by: Daniel Baumann --- third_party/rust/plane-split/src/polygon.rs | 654 ++++++++++++++++++++++++++++ 1 file changed, 654 insertions(+) create mode 100644 third_party/rust/plane-split/src/polygon.rs (limited to 'third_party/rust/plane-split/src/polygon.rs') diff --git a/third_party/rust/plane-split/src/polygon.rs b/third_party/rust/plane-split/src/polygon.rs new file mode 100644 index 0000000000..5929263346 --- /dev/null +++ b/third_party/rust/plane-split/src/polygon.rs @@ -0,0 +1,654 @@ +use crate::{is_zero, Line, Plane}; + +use euclid::{ + approxeq::ApproxEq, + default::{Point2D, Point3D, Rect, Transform3D, Vector3D}, +}; +use smallvec::SmallVec; + +use std::{iter, mem}; + +/// The projection of a `Polygon` on a line. +pub struct LineProjection { + /// Projected value of each point in the polygon. + pub markers: [f64; 4], +} + +impl LineProjection { + /// Get the min/max of the line projection markers. + pub fn get_bounds(&self) -> (f64, f64) { + 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 { + /// 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 Intersection { + /// 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 { + /// Points making the polygon. + pub points: [Point3D; 4], + /// A plane describing polygon orientation. + pub plane: Plane, + /// A simple anchoring index to allow association of the + /// produced split polygons with the original one. + pub anchor: A, +} + +impl Clone for Polygon { + 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 Polygon +where + 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; 4], anchor: A) -> Option { + 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() < f64::EPSILON || edge4.square_length() < f64::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, anchor: A) -> Self { + let min = rect.min(); + let max = rect.max(); + Polygon { + points: [ + min.to_3d(), + Point3D::new(max.x, min.y, 0.0), + max.to_3d(), + Point3D::new(min.x, max.y, 0.0), + ], + plane: Plane { + normal: Vector3D::new(0.0, 0.0, 1.0), + offset: 0.0, + }, + anchor, + } + } + + /// Construct a polygon from a rectangle with 3D transform. + pub fn from_transformed_rect( + rect: Rect, + transform: Transform3D, + anchor: A, + ) -> Option { + let min = rect.min(); + let max = rect.max(); + let points = [ + transform.transform_point3d(min.to_3d())?, + transform.transform_point3d(Point3D::new(max.x, min.y, 0.0))?, + transform.transform_point3d(max.to_3d())?, + transform.transform_point3d(Point3D::new(min.x, max.y, 0.0))?, + ]; + Self::from_points(points, anchor) + } + + /// Construct a polygon from a rectangle with an invertible 3D transform. + pub fn from_transformed_rect_with_inverse( + rect: Rect, + transform: &Transform3D, + inv_transform: &Transform3D, + anchor: A, + ) -> Option { + let min = rect.min(); + let max = rect.max(); + let points = [ + transform.transform_point3d(min.to_3d())?, + transform.transform_point3d(Point3D::new(max.x, min.y, 0.0))?, + transform.transform_point3d(max.to_3d())?, + transform.transform_point3d(Point3D::new(min.x, max.y, 0.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(&0.0) || transform.m44.approx_eq(&0.0) { + 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) -> Point2D { + //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(&self, transform: &Transform3D) -> Option> { + 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(f64::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(0.0, 0.0, T::one(), 0.0)); + //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) >= 0.0); + 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() < f64::EPSILON + || (self.points[1] - self.points[3]).square_length() < f64::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) -> LineProjection { + 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) -> Intersection { + 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 { + 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), + second: (usize, Point3D), + ) -> (Option, Option) { + //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) -> (Option, Option) { + 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 >= 0.0 && t < 1.0 { + *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, + normal: &Vector3D, + ) -> (Option, Option) { + log::debug!("\tSplitting with normal"); + // figure out which side of the split does each point belong to + let mut sides = [0.0; 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 < 0.0 && side1 >= 0.0 { + &mut cut_positive + } else if side0 > 0.0 && side1 <= 0.0 { + &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) + } + } + + /// Cut a polygon with another one. + /// + /// Write the resulting polygons in `front` and `back` if the polygon needs to be split. + pub fn cut( + &self, + poly: &Self, + front: &mut SmallVec<[Polygon; 2]>, + back: &mut SmallVec<[Polygon; 2]>, + ) -> PlaneCut { + //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); + 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), 0.0) + } + }; + + 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) => PlaneCut::Sibling, + Intersection::Coplanar | Intersection::Outside => { + if dist > 0.0 { + front.push(poly.clone()); + } else { + back.push(poly.clone()); + } + + PlaneCut::Cut + } + Intersection::Inside(line) => { + let mut poly = poly.clone(); + let (res_add1, res_add2) = poly.split_with_normal(&line, &self.plane.normal); + + 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 > 0.0 { + front.push(sub) + } else { + back.push(sub) + } + } + + PlaneCut::Cut + } + } + } + + /// Returns whether both polygon's planes are parallel. + pub fn is_aligned(&self, other: &Self) -> bool { + self.plane.normal.dot(other.plane.normal) > 0.0 + } +} + +/// The result of a polygon being cut by a plane. +/// The "cut" here is an attempt to classify a plane as being +/// in front or in the back of another one. +#[derive(Debug, PartialEq)] +pub enum PlaneCut { + /// The planes are one the same geometrical plane. + Sibling, + /// Planes are different, thus we can either determine that + /// our plane is completely in front/back of another one, + /// or split it into these sub-groups. + Cut, +} + +#[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); +} -- cgit v1.2.3