beehive 0.1.1

//! A parallelogram prism aligned to X and Y axes.

use std::iter::{Extend, FromIterator, FusedIterator, Iterator};
use std::ops::{Add, AddAssign, Div, Mul, Neg, Sub, SubAssign};

use num_traits::{AsPrimitive, One, Zero};

use super::coords::{Point, PointAxial, Vector, VectorAxial};

/// A parallelogram prism aligned to X and Y axes.
///
/// `T` should generally be a signed integer type. Floating point numbers have
/// precision problems that may lead to invalid coordinates after operations.
#[derive(Copy, Clone, Eq, PartialEq, Hash, Debug, Default)]
#[cfg_attr(feature = "serde-1", derive(Serialize, Deserialize))]
pub struct QuadPrism<T> {
    pub low: PointAxial<T>,
    pub high: PointAxial<T>,
}

fn sort2<T: PartialOrd>(a: T, b: T) -> (T, T) {
    use std::cmp::Ordering as O;
    match a.partial_cmp(&b) {
        Some(O::Less) | Some(O::Equal) => (a, b),
        Some(O::Greater) | None => (b, a),
    }
}

fn sort2_ref<'a, T: PartialOrd>(a: &'a T, b: &'a T) -> (&'a T, &'a T) {
    use std::cmp::Ordering as O;
    match a.partial_cmp(b) {
        Some(O::Less) | Some(O::Equal) => (a, b),
        Some(O::Greater) | None => (b, a),
    }
}

fn bounds<T: PartialOrd>(a: PointAxial<T>, b: PointAxial<T>) -> (PointAxial<T>, PointAxial<T>) {
    let (min_x, max_x) = sort2(a.x, b.x);
    let (min_y, max_y) = sort2(a.y, b.y);
    let (min_w, max_w) = sort2(a.w, b.w);

    (
        PointAxial {
            x: min_x,
            y: min_y,
            w: min_w,
        },
        PointAxial {
            x: max_x,
            y: max_y,
            w: max_w,
        },
    )
}

fn bounds_in_place<T: PartialOrd>(a: &mut PointAxial<T>, b: &mut PointAxial<T>) {
    if a.x > b.x {
        std::mem::swap(&mut a.x, &mut b.x);
    }
    if a.y > b.y {
        std::mem::swap(&mut a.y, &mut b.y);
    }
    if a.w > b.w {
        std::mem::swap(&mut a.w, &mut b.w);
    }
}

fn clamp<T: PartialOrd>(low: T, high: T, val: T) -> T {
    if val < low {
        low
    } else if val > high {
        high
    } else {
        val
    }
}

impl<T: PartialOrd + One + Add<Output = T>> QuadPrism<T> {
    /// Creates directly from low and high coordinate points, excluding `high`.
    /// If any component of `high` is lower than `low`, the resulting prism
    /// will contain no points.
    pub fn new(low: PointAxial<T>, high: PointAxial<T>) -> Self {
        QuadPrism { low, high }
    }

    /// Creates a quad prism that contains both `a` and `b`, inclusive.
    pub fn from_points_axial(a: PointAxial<T>, b: PointAxial<T>) -> Self {
        let (low, high) = bounds(a, b);

        QuadPrism {
            low,
            high: high + VectorAxial::one(),
        }
    }

    /// Creates a quad prism that contains both `a` and `b`, inclusive.
    pub fn from_points(a: Point<T>, b: Point<T>) -> Self {
        Self::from_points_axial(a.into_axial(), b.into_axial())
    }

    /// Creates a quad prism from a base point and a size vector. Negative
    /// `size` allowed.
    pub fn from_base_size_axial(base: PointAxial<T>, size: VectorAxial<T>) -> Self
    where
        T: Copy + Zero + One + AddAssign + SubAssign + Neg<Output = T>,
    {
        use std::cmp::Ordering as O;

        let mut canon_base = base;
        let mut canon_size = size;

        macro_rules! canonicalize_comp {
            ($f:ident) => {
                if let Some(O::Less) = size.$f.partial_cmp(&T::zero()) {
                    canon_size.$f += T::one();
                    canon_base.$f += canon_size.$f;
                    canon_size.$f = -canon_size.$f + T::one();
                }
            };
        }

        canonicalize_comp!(x);
        canonicalize_comp!(y);
        canonicalize_comp!(w);

        let (low, high) = bounds(canon_base, canon_base + canon_size);

        QuadPrism { low, high }
    }

    /// Creates a quad prism from a base point and a size vector. Negative
    /// `size` allowed.
    pub fn from_base_size(base: Point<T>, size: Vector<T>) -> Self
    where
        T: Copy + Zero + One + AddAssign + SubAssign + Neg<Output = T>,
    {
        Self::from_base_size_axial(base.into_axial(), size.into_axial())
    }

    /// Creates a quad prism from a center point and a size vector. Negative
    /// `size` allowed.
    pub fn from_center_size_axial(mut center: PointAxial<T>, size: VectorAxial<T>) -> Self
    where
        T: Copy + Zero + One + AddAssign + SubAssign + Neg<Output = T> + Div<Output = T>,
    {
        let half = size / (T::one() + T::one());
        center -= half;
        Self::from_base_size_axial(center, size)
    }

    /// Creates a quad prism from a center point and a size vector. Negative
    /// `size` allowed.
    pub fn from_center_size(center: Point<T>, size: Vector<T>) -> Self
    where
        T: Copy + Zero + One + AddAssign + SubAssign + Neg<Output = T> + Div<Output = T>,
    {
        Self::from_center_size_axial(center.into_axial(), size.into_axial())
    }

    /// Expand to include the point.
    pub fn expand_into_point_axial(&mut self, mut point: PointAxial<T>) {
        bounds_in_place(&mut self.low, &mut point);
        point = point + VectorAxial::one();
        bounds_in_place(&mut point, &mut self.high);
    }

    /// Expand to include the point.
    pub fn expand_into_point(&mut self, point: Point<T>) {
        self.expand_into_point_axial(point.into_axial());
    }

    /// Limit the prism's span on X-axis to 1. Returns `None` if `x` is out of
    /// bounds.
    pub fn column(self, x: T) -> Option<Self>
    where
        T: Copy,
    {
        if x >= self.low.x && x < self.high.x {
            Some(QuadPrism {
                low: PointAxial { x, ..self.low },
                high: PointAxial {
                    x: x + T::one(),
                    ..self.high
                },
            })
        } else {
            None
        }
    }

    /// Limit the prism's span on Y-axis to 1. Returns `None` if `y` is out of
    /// bounds.
    pub fn row(self, y: T) -> Option<Self>
    where
        T: Copy,
    {
        if y >= self.low.y && y < self.high.y {
            Some(QuadPrism {
                low: PointAxial { y, ..self.low },
                high: PointAxial {
                    y: y + T::one(),
                    ..self.high
                },
            })
        } else {
            None
        }
    }

    /// Limit the prism's span on W-axis to 1. Returns `None` if `w` is out of
    /// bounds.
    pub fn layer(self, w: T) -> Option<Self>
    where
        T: Copy,
    {
        if w >= self.low.w && w < self.high.w {
            Some(QuadPrism {
                low: PointAxial { w, ..self.low },
                high: PointAxial {
                    w: w + T::one(),
                    ..self.high
                },
            })
        } else {
            None
        }
    }
}

impl<T: Copy + PartialOrd + Mul<Output = T> + Zero + Sub<Output = T>> QuadPrism<T> {
    /// Returns the size vector of the prism, clamped to non-negative X, Y,
    /// and W components.
    pub fn size_axial(&self) -> VectorAxial<T> {
        let v = self.high - self.low;
        VectorAxial {
            x: sort2(T::zero(), v.x).1,
            y: sort2(T::zero(), v.y).1,
            w: sort2(T::zero(), v.w).1,
        }
    }

    /// Returns the size vector of the prism, clamped to non-negative X, Y,
    /// and W components.
    ///
    /// Z will be non-positive as per the plane constraint.
    pub fn size(&self) -> Vector<T> {
        self.size_axial().to_cube()
    }

    /// Returns the volume of the prism.
    pub fn volume(&self) -> i64
    where
        T: AsPrimitive<i64>,
    {
        let size = self.size_axial().cast();
        size.x * size.y * size.w
    }

    /// Returns the center point of the prism.
    pub fn center_axial(&self) -> PointAxial<T>
    where
        T: One + Div<Output = T>,
    {
        let two = T::one() + T::one();
        PointAxial {
            x: (self.low.x + self.high.x - T::one()) / two,
            y: (self.low.y + self.high.y - T::one()) / two,
            w: (self.low.w + self.high.w - T::one()) / two,
        }
    }

    /// Returns the center point of the prism.
    pub fn center(&self) -> Point<T>
    where
        T: One + Div<Output = T>,
    {
        self.center_axial().to_cube()
    }
}

impl<T: PartialOrd> QuadPrism<T> {
    /// Returns `true` if the prism contains no points.
    pub fn is_empty(&self) -> bool {
        self.high.x <= self.low.x || self.high.y <= self.low.y || self.high.w <= self.low.w
    }

    /// Returns `true` if the point is contained in the prism.
    pub fn contains_axial(&self, pt: &PointAxial<T>) -> bool {
        pt.x >= self.low.x
            && pt.y >= self.low.y
            && pt.w >= self.low.w
            && pt.x < self.high.x
            && pt.y < self.high.y
            && pt.w < self.high.w
    }

    /// Returns `true` if the point is contained in the prism.
    pub fn contains(&self, pt: &Point<T>) -> bool {
        *pt.x() >= self.low.x
            && *pt.y() >= self.low.y
            && *pt.w() >= self.low.w
            && *pt.x() < self.high.x
            && *pt.y() < self.high.y
            && *pt.w() < self.high.w
    }

    /// Returns the intersection between two prisms. The result may be empty.
    ///
    /// It's possible to check for intersection by calling `is_empty` on the
    /// result. See also `intersects`.
    pub fn intersection(self, other: QuadPrism<T>) -> QuadPrism<T> {
        let (_, low) = bounds(self.low, other.low);
        let (high, _) = bounds(self.high, other.high);
        QuadPrism { low, high }
    }

    /// Returns `true` if the prisms intersect, without consuming them.
    ///
    /// Performance should be the same as using `intersection` and `is_empty`.
    /// The difference is in semantics and convenience.
    pub fn intersects(&self, other: &QuadPrism<T>) -> bool {
        sort2_ref(&self.low.x, &other.low.x).1 < sort2_ref(&self.high.x, &other.high.x).0
            && sort2_ref(&self.low.y, &other.low.y).1 < sort2_ref(&self.high.y, &other.high.y).0
            && sort2_ref(&self.low.w, &other.low.w).1 < sort2_ref(&self.high.w, &other.high.w).0
    }

    /// Returns `true` if this prism completely covers `other`, without
    /// consuming them.
    ///
    /// Performance should be the same as using `intersection` and `eq`.
    /// The difference is in semantics and convenience.
    pub fn covers(&self, other: &QuadPrism<T>) -> bool
    where
        T: PartialEq,
    {
        sort2_ref(&self.low.x, &other.low.x).1 == &other.low.x
            && sort2_ref(&self.low.y, &other.low.y).1 == &other.low.y
            && sort2_ref(&self.low.w, &other.low.w).1 == &other.low.w
            && sort2_ref(&self.high.x, &other.high.x).0 == &other.high.x
            && sort2_ref(&self.high.y, &other.high.y).0 == &other.high.y
            && sort2_ref(&self.high.w, &other.high.w).0 == &other.high.w
    }

    /// Returns a prism that covers both prisms. The result may be non-empty
    /// if both prisms are empty, but placed at different positions.
    ///
    /// It is guaranteed that:
    ///
    /// - The results are equal whichever way the function is called.
    /// - `covers` on the result returns `true` for both input prisms.
    pub fn union(self, other: QuadPrism<T>) -> QuadPrism<T> {
        let (low, _) = bounds(self.low, other.low);
        let (_, high) = bounds(self.high, other.high);
        QuadPrism { low, high }
    }

    /// Returns a prism that covers all points covered by either prism. Differs
    /// from `union` that `covers` isn't guaranteed to return `true`.
    ///
    /// It is guaranteed that:
    ///
    /// - Unions between empty prisms, wherever they're placed, are empty.
    /// - Unions between an empty prism and a non-empty one, is equal to the
    /// non-empty one.
    /// - When both prisms are non-empty, the results are equal whichever way
    /// the function is called.
    pub fn union_short_circuit(self, other: QuadPrism<T>) -> QuadPrism<T> {
        if self.is_empty() {
            other
        } else if other.is_empty() {
            self
        } else {
            self.union(other)
        }
    }

    /// Clamps a point to the bounds of this prism. Returns `None` if the prism
    /// is empty.
    pub fn clamp(self, v: Point<T>) -> Option<Point<T>>
    where
        T: Copy + Sub<Output = T> + One + Zero,
    {
        self.clamp_axial(v.into_axial()).map(PointAxial::to_cube)
    }

    /// Clamps a point to the bounds of this prism.
    pub fn clamp_axial(self, v: PointAxial<T>) -> Option<PointAxial<T>>
    where
        T: Sub<Output = T> + One,
    {
        if self.is_empty() {
            return None;
        }

        Some(PointAxial {
            x: clamp(self.low.x, self.high.x - T::one(), v.x),
            y: clamp(self.low.y, self.high.y - T::one(), v.y),
            w: clamp(self.low.w, self.high.w - T::one(), v.w),
        })
    }

    /// Returns an iterator through all points in a prism from -X to +X,
    /// -Y to +Y, bottom-up, in that order.
    pub fn points_axial(&self) -> PointsAxial<T>
    where
        T: Copy + PartialOrd + AddAssign + One,
    {
        PointsAxial {
            prism: self,
            current: self.low,
        }
    }

    /// Returns an iterator through all points in a prism from -X to +X,
    /// -Y to +Y, bottom-up, in that order.
    pub fn points(&self) -> Points<T>
    where
        T: Copy + PartialOrd + AddAssign + One,
    {
        Points(PointsAxial {
            prism: self,
            current: self.low,
        })
    }
}

/// Iterator through all points in a prism from -X to +X, -Y to +Y, bottom-up,
/// in that order.
#[derive(Clone, Eq, PartialEq, Hash, Debug)]
pub struct PointsAxial<'a, T> {
    prism: &'a QuadPrism<T>,
    current: PointAxial<T>,
}

impl<'a, T> Iterator for PointsAxial<'a, T>
where
    T: Copy + PartialOrd + AddAssign + One,
{
    type Item = PointAxial<T>;

    fn next(&mut self) -> Option<Self::Item> {
        if !self.prism.contains_axial(&self.current) {
            return None;
        }
        let result = self.current;
        self.current.x += T::one();
        if self.current.x >= self.prism.high.x {
            self.current.x = self.prism.low.x;
            self.current.y += T::one();
        }
        if self.current.y >= self.prism.high.y {
            self.current.y = self.prism.low.y;
            self.current.w += T::one();
        }
        Some(result)
    }
}

impl<'a, T> FusedIterator for PointsAxial<'a, T> where T: Copy + PartialOrd + AddAssign + One {}

/// Iterator through all points in a prism from -X to +X, -Y to +Y, bottom-up,
/// in that order.
#[derive(Clone, Eq, PartialEq, Hash, Debug)]
pub struct Points<'a, T>(PointsAxial<'a, T>);

impl<'a, T> Iterator for Points<'a, T>
where
    T: Copy + PartialOrd + Sub<Output = T> + AddAssign + Zero + One,
{
    type Item = Point<T>;

    fn next(&mut self) -> Option<Self::Item> {
        self.0.next().map(PointAxial::to_cube)
    }
}

impl<'a, T> FusedIterator for Points<'a, T> where
    T: Copy + PartialOrd + Sub<Output = T> + AddAssign + Zero + One
{
}

impl<T> Extend<PointAxial<T>> for QuadPrism<T>
where
    T: PartialOrd + One + Add<Output = T>,
{
    fn extend<I: IntoIterator<Item = PointAxial<T>>>(&mut self, it: I) {
        for pt in it {
            self.expand_into_point_axial(pt);
        }
    }
}

impl<T> Extend<Point<T>> for QuadPrism<T>
where
    T: PartialOrd + One + Add<Output = T>,
{
    fn extend<I: IntoIterator<Item = Point<T>>>(&mut self, it: I) {
        for pt in it {
            self.expand_into_point(pt);
        }
    }
}

impl<T> FromIterator<PointAxial<T>> for QuadPrism<T>
where
    T: Copy + Default + PartialOrd + One + Add<Output = T>,
{
    fn from_iter<I: IntoIterator<Item = PointAxial<T>>>(it: I) -> Self {
        let mut iter = it.into_iter();

        let mut prism = {
            if let Some(first) = iter.next() {
                QuadPrism::from_points_axial(first, first)
            } else {
                return QuadPrism::default();
            }
        };

        prism.extend(iter);
        prism
    }
}

impl<T> FromIterator<Point<T>> for QuadPrism<T>
where
    T: Copy + Default + PartialOrd + One + Add<Output = T>,
{
    fn from_iter<I: IntoIterator<Item = Point<T>>>(it: I) -> Self {
        it.into_iter().map(Point::into_axial).collect()
    }
}

#[cfg(feature = "rand-07")]
mod distribution_impl {
    use super::{Point, PointAxial, QuadPrism};

    use std::ops::Sub;

    use num_traits::Zero;
    use rand::distributions::{uniform::SampleUniform, Distribution, Uniform};
    use rand::Rng;

    /// Distribution for sampling uniformly within the quad-prism.
    impl<T> Distribution<Point<T>> for QuadPrism<T>
    where
        T: Copy + Zero + Sub<Output = T> + SampleUniform,
    {
        fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Point<T> {
            let axial: PointAxial<_> = rng.sample(self);
            axial.to_cube()
        }
    }

    /// Distribution for sampling uniformly within the quad-prism.
    impl<T> Distribution<PointAxial<T>> for QuadPrism<T>
    where
        T: Copy + SampleUniform,
    {
        fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> PointAxial<T> {
            let x = Uniform::new(self.low.x, self.high.x).sample(rng);
            let y = Uniform::new(self.low.y, self.high.y).sample(rng);
            let w = Uniform::new(self.low.w, self.high.w).sample(rng);
            PointAxial::new(x, y, w)
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    type HA = super::super::coords::PointAxial<i16>;
    type VA = super::super::coords::VectorAxial<i16>;

    #[test]
    fn can_construct_from_points() {
        fn test_from_points_sanity(a: HA, b: HA, expected_volume: i64) {
            let from_axial = QuadPrism::from_points_axial(a, b);
            let from_cube = QuadPrism::from_points(a.to_cube(), b.to_cube());
            let from_expand = {
                let mut q = QuadPrism::from_points_axial(a, a);
                q.expand_into_point_axial(b);
                q
            };
            let from_expand_cube = {
                let mut q = QuadPrism::from_points(a.to_cube(), a.to_cube());
                q.expand_into_point(b.to_cube());
                q
            };
            assert_eq!(from_axial, from_cube);
            assert_eq!(from_axial, from_expand);
            assert_eq!(from_axial, from_expand_cube);
            assert!(from_axial.contains_axial(&a));
            assert!(from_axial.contains_axial(&b));
            assert!(from_cube.contains_axial(&a));
            assert!(from_cube.contains_axial(&b));
            assert!(from_axial.contains(&a.to_cube()));
            assert!(from_axial.contains(&b.to_cube()));
            assert!(from_cube.contains(&a.to_cube()));
            assert!(from_cube.contains(&b.to_cube()));
            assert_eq!(expected_volume, from_axial.volume());
            assert_eq!(expected_volume, from_cube.volume());
        }

        test_from_points_sanity(HA::new(0, 0, 0), HA::new(0, 0, 0), 1);
        test_from_points_sanity(HA::new(1, 2, 3), HA::new(-1, -2, -3), 3 * 5 * 7);
        test_from_points_sanity(HA::new(-1, 0, 1), HA::new(2, 1, -1), 4 * 2 * 3);
        test_from_points_sanity(HA::new(1, -1, 5), HA::new(0, -1, 6), 2 * 2);
    }

    #[test]
    fn can_construct_from_base_size() {
        fn test_from_base_size_sanity(base: HA, size: VA, expected_volume: i64) {
            let from_axial = QuadPrism::from_base_size_axial(base, size);
            let from_cube = QuadPrism::from_base_size(base.to_cube(), size.to_cube());
            assert_eq!(from_axial, from_cube);
            if expected_volume > 0 {
                assert!(from_axial.contains_axial(&base));
                assert!(!from_axial.contains_axial(&(base + size)));
                assert!(from_cube.contains_axial(&base));
                assert!(!from_cube.contains_axial(&(base + size)));
                assert!(from_axial.contains(&base.to_cube()));
                assert!(!from_axial.contains(&(base + size).to_cube()));
                assert!(from_cube.contains(&base.to_cube()));
                assert!(!from_cube.contains(&(base + size).to_cube()));
            }
            assert_eq!(expected_volume, from_axial.volume());
            assert_eq!(expected_volume, from_cube.volume());
        }

        test_from_base_size_sanity(HA::new(0, 0, 0), VA::new(-1, 2, 1), 2);
        test_from_base_size_sanity(HA::new(1, 2, 3), VA::new(-7, 0, 2), 0);
        test_from_base_size_sanity(HA::new(-1, 0, 1), VA::new(3, -7, -2), 3 * 7 * 2);
        test_from_base_size_sanity(HA::new(1, -1, 5), VA::new(5, 1, -5), 5 * 5);
    }

    #[test]
    fn can_construct_from_center_size() {
        assert_eq!(
            QuadPrism::from_base_size_axial(HA::new(1, 1, 1), VA::new(3, 3, 1)),
            QuadPrism::from_center_size_axial(HA::new(2, 2, 1), VA::new(3, 3, 1)),
        );

        assert_eq!(
            QuadPrism::from_base_size(HA::new(1, 1, 1).to_cube(), VA::new(3, 3, 1).to_cube()),
            QuadPrism::from_center_size(HA::new(2, 2, 1).to_cube(), VA::new(3, 3, 1).to_cube()),
        );
    }

    #[test]
    fn can_intersect() {
        fn test_intersection(
            a: QuadPrism<i16>,
            b: QuadPrism<i16>,
            expected: Option<QuadPrism<i16>>,
        ) {
            let intersection_a = a.intersection(b);
            let intersection_b = b.intersection(a);
            assert_eq!(intersection_a, intersection_b);
            if let Some(exp) = expected {
                assert_eq!(exp, intersection_a);
                assert!(!intersection_a.is_empty());
            } else {
                assert!(intersection_a.is_empty());
            }
            assert_eq!(expected.is_some(), a.intersects(&b));
            assert_eq!(expected.is_some(), b.intersects(&a));
        }

        test_intersection(
            QuadPrism::from_points_axial(HA::new(0, 0, 0), HA::new(3, 3, 0)),
            QuadPrism::from_points_axial(HA::new(-3, -3, 0), HA::new(0, 0, 0)),
            Some(QuadPrism::from_points_axial(
                HA::new(0, 0, 0),
                HA::new(0, 0, 0),
            )),
        );

        test_intersection(
            QuadPrism::from_points_axial(HA::new(-2, -1, -3), HA::new(3, 3, 3)),
            QuadPrism::from_points_axial(HA::new(-3, -3, 2), HA::new(2, 1, -2)),
            Some(QuadPrism::from_points_axial(
                HA::new(-2, -1, -2),
                HA::new(2, 1, 2),
            )),
        );

        test_intersection(
            QuadPrism::from_points_axial(HA::new(2, 2, 0), HA::new(3, 3, 0)),
            QuadPrism::from_points_axial(HA::new(-3, -3, 0), HA::new(-2, -2, 0)),
            None,
        );
    }

    #[test]
    fn can_cover() {
        fn test_cover(
            expected: bool,
            expected_reverse: bool,
            a: QuadPrism<i16>,
            b: QuadPrism<i16>,
        ) {
            assert_eq!(expected, a.covers(&b));
            assert_eq!(expected_reverse, b.covers(&a));
            if expected {
                assert_eq!(b, a.intersection(b));
            }
            if expected_reverse {
                assert_eq!(a, b.intersection(a));
            }
        }

        test_cover(
            true,
            true,
            QuadPrism::from_points_axial(HA::new(0, 0, 0), HA::new(1, 1, 1)),
            QuadPrism::from_points_axial(HA::new(1, 1, 1), HA::new(0, 0, 0)),
        );

        test_cover(
            true,
            false,
            QuadPrism::from_points_axial(HA::new(-5, -6, -7), HA::new(5, 6, 7)),
            QuadPrism::from_points_axial(HA::new(-2, -3, -2), HA::new(1, 1, 1)),
        );

        test_cover(
            false,
            false,
            QuadPrism::from_points_axial(HA::new(-2, -5, -2), HA::new(1, 4, 1)),
            QuadPrism::from_points_axial(HA::new(-5, -1, -7), HA::new(5, 6, 3)),
        );

        test_cover(
            true,
            false,
            QuadPrism::from_points_axial(HA::new(-5, -6, -7), HA::new(5, 6, 7)),
            QuadPrism::from_points_axial(HA::new(-2, -2, -2), HA::new(-1, -1, -1)).intersection(
                QuadPrism::from_points_axial(HA::new(2, 2, 2), HA::new(1, 1, 1)),
            ),
        );
    }

    #[test]
    fn can_iter() {
        fn test_iter(qp: QuadPrism<i16>, points: &[HA]) {
            itertools::assert_equal(qp.points_axial(), points.iter().cloned());
            itertools::assert_equal(qp.points(), points.iter().cloned().map(HA::to_cube));
        }

        test_iter(
            QuadPrism::from_points_axial(HA::new(-1, -1, 0), HA::new(0, 0, 1)),
            &[
                HA::new(-1, -1, 0),
                HA::new(0, -1, 0),
                HA::new(-1, 0, 0),
                HA::new(0, 0, 0),
                HA::new(-1, -1, 1),
                HA::new(0, -1, 1),
                HA::new(-1, 0, 1),
                HA::new(0, 0, 1),
            ],
        );

        test_iter(
            QuadPrism::from_base_size_axial(HA::new(0, 0, 0), VA::new(0, 0, 0)),
            &[],
        );
        test_iter(
            QuadPrism::from_base_size_axial(HA::new(0, 0, 0), VA::new(1, 1, 1)),
            &[HA::new(0, 0, 0)],
        );
    }

    #[test]
    fn can_union() {
        fn test_union(
            a: QuadPrism<i16>,
            b: QuadPrism<i16>,
            expected: QuadPrism<i16>,
            short_circuit_expected: Option<QuadPrism<i16>>,
        ) {
            let union_a = a.union(b);
            let union_b = b.union(a);
            assert_eq!(expected, union_a);
            assert_eq!(union_a, union_b);
            assert!(union_a.covers(&a));
            assert!(union_a.covers(&b));

            let union_sc_a = a.union_short_circuit(b);
            let union_sc_b = b.union_short_circuit(a);
            if let Some(exp) = short_circuit_expected {
                assert_eq!(exp, union_sc_a);
                assert_eq!(union_sc_a, union_sc_b);
                assert!(!union_sc_a.is_empty());
            } else {
                assert!(union_sc_a.is_empty());
                assert!(union_sc_b.is_empty());
            }
        }

        // Unions between empty prisms
        test_union(
            QuadPrism::from_base_size_axial(HA::new(-3, -3, -2), VA::zero()),
            QuadPrism::from_base_size_axial(HA::new(3, 3, 2), VA::zero()),
            QuadPrism::from_points_axial(HA::new(-3, -3, -2), HA::new(2, 2, 1)),
            None,
        );

        // Unions between one empty prism and a non-empty prism
        test_union(
            QuadPrism::default(),
            QuadPrism::from_points_axial(HA::new(-3, -3, -3), HA::new(-2, -2, -2)),
            QuadPrism::from_points_axial(HA::new(-3, -3, -3), HA::new(-1, -1, -1)),
            Some(QuadPrism::from_points_axial(
                HA::new(-3, -3, -3),
                HA::new(-2, -2, -2),
            )),
        );

        // Normal union
        test_union(
            QuadPrism::from_points_axial(HA::new(-2, -1, -3), HA::new(3, 3, 3)),
            QuadPrism::from_points_axial(HA::new(-3, -3, 2), HA::new(2, 1, -2)),
            QuadPrism::from_points_axial(HA::new(-3, -3, -3), HA::new(3, 3, 3)),
            Some(QuadPrism::from_points_axial(
                HA::new(-3, -3, -3),
                HA::new(3, 3, 3),
            )),
        );
    }

    #[test]
    fn can_get_row_column_layer() {
        let prism = QuadPrism::from_points_axial(HA::new(-2, -3, -4), HA::new(1, 2, 3));

        assert_eq!(None, prism.column(-3));
        assert_eq!(None, prism.column(2));
        assert_eq!(
            Some(QuadPrism::from_points_axial(
                HA::new(-1, -3, -4),
                HA::new(-1, 2, 3)
            )),
            prism.column(-1),
        );

        assert_eq!(None, prism.row(-4));
        assert_eq!(None, prism.row(3));
        assert_eq!(
            Some(QuadPrism::from_points_axial(
                HA::new(-2, 1, -4),
                HA::new(1, 1, 3)
            )),
            prism.row(1),
        );

        assert_eq!(None, prism.layer(-5));
        assert_eq!(None, prism.layer(4));
        assert_eq!(
            Some(QuadPrism::from_points_axial(
                HA::new(-2, -3, 2),
                HA::new(1, 2, 2)
            )),
            prism.layer(2),
        );
    }

    #[test]
    fn can_clamp() {
        let prism = QuadPrism::from_points_axial(HA::new(-2, -3, -4), HA::new(1, 2, 3));
        assert_eq!(
            Some(HA::new(-2, 2, -4)),
            prism.clamp_axial(HA::new(-3, 3, -5))
        );
        assert_eq!(
            Some(HA::new(1, -3, 3)),
            prism.clamp_axial(HA::new(2, -4, 4))
        );
        assert_eq!(
            Some(HA::new(-1, 0, 1)),
            prism.clamp_axial(HA::new(-1, 0, 1))
        );

        let prism = QuadPrism::default();
        assert_eq!(None, prism.clamp_axial(HA::new(0, 0, 0)));
        assert_eq!(None, prism.clamp_axial(HA::new(1, -1, 2)));
    }
}