summaryrefslogtreecommitdiffstats
path: root/third_party/rust/plane-split/src/polygon.rs
diff options
context:
space:
mode:
Diffstat (limited to 'third_party/rust/plane-split/src/polygon.rs')
-rw-r--r--third_party/rust/plane-split/src/polygon.rs654
1 files changed, 654 insertions, 0 deletions
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<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<A> {
+ /// Points making the polygon.
+ pub points: [Point3D<f64>; 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<A: Copy> Clone for Polygon<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<A> Polygon<A>
+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<f64>; 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() < 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<f64>, 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<f64>,
+ transform: Transform3D<f64>,
+ anchor: A,
+ ) -> Option<Self> {
+ 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<f64>,
+ transform: &Transform3D<f64>,
+ inv_transform: &Transform3D<f64>,
+ anchor: A,
+ ) -> Option<Self> {
+ 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<f64>) -> Point2D<f64> {
+ //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<f64>) -> Option<Polygon<A>> {
+ 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<f64>) -> 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<Line> {
+ 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> {
+ 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<f64>),
+ second: (usize, Point3D<f64>),
+ ) -> (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) -> (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 >= 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<f64>,
+ ) -> (Option<Self>, Option<Self>) {
+ 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<A>; 2]>,
+ back: &mut SmallVec<[Polygon<A>; 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);
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-24 11:13:19 +0000