dcel 0.8.12

// SPDX-FileCopyrightText: 2026 dcel contributors
//
// SPDX-License-Identifier: MIT OR Apache-2.0

use std::collections::BTreeSet;
use std::marker::PhantomData;

use maplike::{Clear, Get, Insert, KeyedCollection, Push, Remove, StableRemove};
use rstar::{
    AABB, Envelope, Point,
    primitives::{GeomWithData, Rectangle},
};
use rstared::AsRefRTree;
use stable_vec::StableVec;

#[cfg(feature = "undoredo")]
use undoredo::{ApplyDelta, Delta, FlushDelta, Recorder};

use crate::{
    Dcel, EdgeId, Face, FaceId, HalfEdge, HalfEdgeId, Vertex, VertexId,
    track::{HalfEdgesCounter, VertexesCounter},
};

#[derive(Clone, Debug)]
pub struct RTreedDcel<
    P: Point,
    VW = P,
    HEW = (),
    FW = (),
    VC = Vec<Vertex<VW>>,
    HEC = Vec<HalfEdge<HEW>>,
    FC = Vec<Face<FW>>,
    ER = AsRefRTree<GeomWithData<Rectangle<P>, EdgeId>>,
    FR = AsRefRTree<GeomWithData<Rectangle<P>, FaceId>>,
> {
    dcel: Dcel<VW, HEW, FW, VC, HEC, FC>,
    edges_rtree: ER,
    faces_rtree: FR,
    point_marker: PhantomData<P>,
}

pub type RTreedStableDcel<
    P,
    VW = P,
    HEW = (),
    FW = (),
    ER = AsRefRTree<GeomWithData<Rectangle<P>, EdgeId>>,
    FR = AsRefRTree<GeomWithData<Rectangle<P>, FaceId>>,
> = RTreedDcel<
    P,
    VW,
    HEW,
    FW,
    StableVec<Vertex<VW>>,
    StableVec<HalfEdge<HEW>>,
    StableVec<Face<FW>>,
    ER,
    FR,
>;

#[cfg(feature = "undoredo")]
pub type RecordingRTreedStableDcel<
    P,
    VW = P,
    HEW = (),
    FW = (),
    ER = AsRefRTree<GeomWithData<Rectangle<P>, EdgeId>>,
    FR = AsRefRTree<GeomWithData<Rectangle<P>, FaceId>>,
> = RTreedDcel<
    P,
    VW,
    HEW,
    FW,
    Recorder<StableVec<Vertex<VW>>>,
    Recorder<StableVec<HalfEdge<HEW>>>,
    Recorder<StableVec<Face<FW>>>,
    Recorder<ER>,
    Recorder<FR>,
>;

impl<
    P: Point,
    VW,
    HEW,
    FW: Default,
    VC: Default,
    HEC: Default,
    FC: Default + Push<usize, Value = Face<FW>>,
    ER: Default,
    FR: Default,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    #[inline]
    pub fn new() -> Self {
        Self {
            dcel: Dcel::new(),
            edges_rtree: Default::default(),
            faces_rtree: Default::default(),
            point_marker: PhantomData,
        }
    }
}

impl<P: Point, VW, HEW, FW, VC, HEC, FC, ER, FR> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR> {
    #[inline]
    pub fn dcel(&self) -> &Dcel<VW, HEW, FW, VC, HEC, FC> {
        &self.dcel
    }

    #[inline]
    pub fn edges_rtree(&self) -> &ER {
        &self.edges_rtree
    }

    #[inline]
    pub fn faces_rtree(&self) -> &FR {
        &self.faces_rtree
    }
}

impl<
    P: Point,
    VW: Clone + Into<P> + Eq + Ord,
    HEW: Clone + Default,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn insert_mesh(
        &mut self,
        face_polygons: impl IntoIterator<Item = impl IntoIterator<Item = VW>>,
    ) -> Vec<FaceId> {
        let faces = self.dcel.insert_mesh(face_polygons);
        let mut deduplicating_edge_set = BTreeSet::new();

        for face in &faces {
            self.add_face_to_rtree(*face);

            for edge in self.dcel.face_edges(*face).collect::<Vec<EdgeId>>() {
                if deduplicating_edge_set.insert(edge) {
                    self.add_edge_to_rtree(edge);
                }
            }
        }

        faces
    }

    pub fn insert_mesh_in_face(
        &mut self,
        face_polygons: impl IntoIterator<Item = impl IntoIterator<Item = VW>>,
    ) -> Vec<FaceId> {
        let faces = self.dcel.insert_mesh_in_face(face_polygons);
        let mut deduplicating_edge_set = BTreeSet::new();

        for face in &faces {
            self.add_face_to_rtree(*face);

            for edge in self.dcel.face_edges(*face).collect::<Vec<EdgeId>>() {
                if deduplicating_edge_set.insert(edge) {
                    self.add_edge_to_rtree(edge);
                }
            }
        }

        faces
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone + Default,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn insert_polygon(&mut self, vertex_weights: impl IntoIterator<Item = VW>) -> FaceId {
        let vertex_weights: Vec<VW> = vertex_weights.into_iter().collect();
        let face = self.dcel.insert_polygon(vertex_weights.clone());
        self.add_face_with_edges_to_rtrees_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(vertex_weights),
        );

        face
    }

    pub fn insert_polygon_in_face(
        &mut self,
        outer_face: FaceId,
        vertex_weights: impl IntoIterator<Item = VW>,
    ) -> FaceId {
        let vertex_weights: Vec<VW> = vertex_weights.into_iter().collect();

        let face = self
            .dcel
            .insert_polygon_in_face(outer_face, vertex_weights.clone());
        self.add_face_with_edges_to_rtrees_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(vertex_weights),
        );

        face
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone + Default,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn insert_edge(
        &mut self,
        from: VertexId,
        to: VertexId,
    ) -> ((HalfEdgeId, HalfEdgeId), FaceId) {
        self.split_face_by_edge(
            from,
            to,
            self.dcel.find_vertices_common_face(from, to).unwrap(),
        )
    }

    pub fn insert_edge_chain(
        &mut self,
        from: VertexId,
        to: VertexId,
        vertex_weights: impl IntoIterator<Item = VW>,
    ) -> (Vec<(HalfEdgeId, HalfEdgeId)>, FaceId) {
        self.split_face_by_edge_chain(
            from,
            to,
            vertex_weights,
            self.dcel.find_vertices_common_face(from, to).unwrap(),
        )
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    fn insert_edge_chain_with_edge_weights(
        &mut self,
        from: VertexId,
        to: VertexId,
        vertex_weights: impl IntoIterator<Item = VW>,
        edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
    ) -> (Vec<(HalfEdgeId, HalfEdgeId)>, FaceId) {
        self.split_face_by_edge_chain_with_all_weights(
            from,
            to,
            vertex_weights,
            edge_weights,
            self.dcel.find_vertices_common_face(from, to).unwrap(),
            FW::default(),
        )
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn insert_polygon_with_all_weights(
        &mut self,
        vertex_weights: impl IntoIterator<Item = VW>,
        edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
        face_weight: FW,
    ) -> FaceId {
        let vertex_weights: Vec<VW> = vertex_weights.into_iter().collect();

        let face = self.dcel.insert_polygon_with_all_weights(
            vertex_weights.clone(),
            edge_weights,
            face_weight,
        );
        self.add_face_with_edges_to_rtrees_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(vertex_weights),
        );

        face
    }

    pub fn insert_polygon_in_face_with_all_weights(
        &mut self,
        outer_face: FaceId,
        vertex_weights: impl IntoIterator<Item = VW>,
        edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
        face_weight: FW,
    ) -> FaceId {
        let vertex_weights: Vec<VW> = vertex_weights.into_iter().collect();

        let face = self.dcel.insert_polygon_in_face_with_all_weights(
            outer_face,
            vertex_weights.clone(),
            edge_weights,
            face_weight,
        );
        self.add_face_with_edges_to_rtrees_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(vertex_weights),
        );

        face
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + StableRemove<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + StableRemove<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + StableRemove<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn remove_face(&mut self, face: FaceId) -> (Vec<EdgeId>, Vec<FaceId>) {
        let mut edges: Vec<EdgeId> = self.dcel.face_edges(face).collect();
        edges.extend(self.dcel.face_spokes(face));
        let interspokes: Vec<FaceId> = self.dcel.face_interspokes(face).collect();

        // Remove the absorbed edges.
        for edge in edges {
            let endpoints = self.dcel.edge_endpoints(edge);

            self.edges_rtree.remove(&GeomWithData::new(
                Rectangle::from_corners(
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.0).clone()),
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.1).clone()),
                ),
                edge,
            ));
        }

        // Remove the absorbing face from the faces R-tree before it changes,
        // which would otherwise invalidate its bbox and make them impossible to
        // access anymore.
        self.faces_rtree.remove(&GeomWithData::new(
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(face)
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
            face,
        ));

        // Remove the absorbed faces from the faces R-tree.
        for interspoke in interspokes {
            self.faces_rtree.remove(&GeomWithData::new(
                Self::rectangle_from_vertex_weights(
                    self.dcel
                        .face_vertices(interspoke)
                        .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
                ),
                interspoke,
            ));
        }

        let (edges, faces) = self.dcel.remove_face(face);

        // Insert the absorbing face back in the R-tree now that its bbox has
        // done changing.
        self.add_face_to_rtree(face);

        (edges, faces)
    }

    pub fn remove_edge(&mut self, edge: EdgeId) -> FaceId {
        let endpoints = self.dcel.edge_endpoints(edge);

        // Remove the absorbed edge.
        self.edges_rtree.remove(&GeomWithData::new(
            Rectangle::from_corners(
                Into::<P>::into(self.dcel.vertex_weight(endpoints.0).clone()),
                Into::<P>::into(self.dcel.vertex_weight(endpoints.1).clone()),
            ),
            edge,
        ));

        // Remove the absorbing face from the faces R-tree before it changes,
        // which would otherwise invalidate its bbox and make them impossible to
        // access anymore.
        self.faces_rtree.remove(&GeomWithData::new(
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(self.dcel.incident_face(edge.lesser()))
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
            self.dcel.incident_face(edge.lesser()),
        ));

        // Remove the absorbed faces from the faces R-tree.
        self.faces_rtree.remove(&GeomWithData::new(
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(self.dcel.opposite_face(edge.lesser()))
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
            self.dcel.opposite_face(edge.lesser()),
        ));

        let absorbing_face = self.dcel.remove_edge(edge);

        // Insert the absorbing face back in the R-tree now that its bbox has
        // done changing.
        self.add_face_to_rtree(absorbing_face);

        absorbing_face
    }
}

impl<
    P: Point,
    VW,
    HEW,
    FW: Default,
    VC: Clear,
    HEC: Clear,
    FC: Clear + Push<usize, Value = Face<FW>>,
    ER: Clear,
    FR: Clear,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn clear(&mut self) {
        self.dcel.clear();
        self.edges_rtree.clear();
        self.faces_rtree.clear();
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + StableRemove<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + StableRemove<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + StableRemove<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn merge_faces_around_vertex(&mut self, inner_vertex: VertexId) {
        let absorbing_face = self
            .dcel
            .incident_face(self.dcel.vertex_representative(inner_vertex));
        self.absorb_faces_around_vertex(absorbing_face, inner_vertex);
    }

    pub fn absorb_faces_around_vertex(&mut self, absorbing_face: FaceId, inner_vertex: VertexId) {
        let initial_half_edge = self.dcel.vertex_representative(inner_vertex);

        let inner_edges: Vec<EdgeId> = self.dcel.spokes(initial_half_edge).collect();
        let perimeter_half_edges: Vec<HalfEdgeId> = self
            .dcel
            .vertex_rim_half_edges(inner_vertex)
            .collect::<Vec<HalfEdgeId>>();

        self.absorb_faces_over_edges_and_vertices_in_perimeter(
            absorbing_face,
            self.dcel
                .interspokes(initial_half_edge)
                .filter(|face| face.id() != absorbing_face.id())
                .collect::<Vec<FaceId>>(),
            inner_edges,
            [inner_vertex],
            &perimeter_half_edges,
        );
    }

    pub fn merge_faces_over_edges_and_vertices(
        &mut self,
        faces: impl IntoIterator<Item = FaceId>,
        edges: impl IntoIterator<Item = EdgeId>,
        vertices: impl IntoIterator<Item = VertexId>,
    ) {
        let mut faces = faces.into_iter();
        let absorbing_face = faces.next().unwrap();

        self.absorb_faces_over_edges_and_vertices(
            absorbing_face,
            faces.filter(|face| face.id() != absorbing_face.id()),
            edges,
            vertices,
        );
    }

    pub fn absorb_faces_over_edges_and_vertices(
        &mut self,
        absorbing_face: FaceId,
        faces: impl IntoIterator<Item = FaceId>,
        edges: impl IntoIterator<Item = EdgeId>,
        vertices: impl IntoIterator<Item = VertexId>,
    ) {
        let edges: Vec<EdgeId> = edges.into_iter().collect();
        let excluded_half_edges: Vec<HalfEdgeId> = edges
            .iter()
            .flat_map(|edge| [edge.lesser(), edge.greater()])
            .collect();

        // Find an initial half-edge that is not in the excluded list. Otherwise,
        // the circulator could end up starting from an excluded half-edge, which would
        // result in an infinite loop, as the termination condition depends on
        // returning to the initial edge.
        let initial_half_edge = self
            .dcel
            .face_half_edges(absorbing_face)
            .find(|half_edge| !excluded_half_edges.contains(half_edge))
            .unwrap();

        let perimeter_half_edges: Vec<HalfEdgeId> = self
            .dcel
            .circulate_half_edges_with_excludes(initial_half_edge, excluded_half_edges.clone())
            .collect();

        self.absorb_faces_over_edges_and_vertices_in_perimeter(
            absorbing_face,
            faces,
            edges,
            vertices,
            &perimeter_half_edges,
        );
    }

    pub fn merge_faces(&mut self, faces: impl IntoIterator<Item = FaceId>) {
        let mut faces = faces.into_iter();
        let absorbing_face = faces.next().unwrap();

        self.absorb_faces(
            absorbing_face,
            faces.filter(|face| face.id() != absorbing_face.id()),
        );
    }

    pub fn absorb_faces(
        &mut self,
        absorbing_face: FaceId,
        faces: impl IntoIterator<Item = FaceId>,
    ) {
        let mut half_edges_counter = HalfEdgesCounter::new();
        let mut vertex_weights_counter = VertexesCounter::new();
        let faces: Vec<FaceId> = faces.into_iter().collect();

        // To detect the merged edges and vertices correctly, the absorbing face
        // has to be visited in addition to the absorbed faces.
        half_edges_counter.visit_face_half_edges(&self.dcel, absorbing_face);
        vertex_weights_counter.visit_face_vertices(&self.dcel, absorbing_face);

        for &face in &faces {
            half_edges_counter.visit_face_half_edges(&self.dcel, face);
            vertex_weights_counter.visit_face_vertices(&self.dcel, face);
        }

        self.absorb_faces_over_edges_and_vertices_in_perimeter(
            absorbing_face,
            faces,
            half_edges_counter
                .inner_edges(&self.dcel)
                .collect::<Vec<EdgeId>>(),
            vertex_weights_counter
                .visited_vertices()
                .filter(|&vertex| {
                    self.dcel
                        .spokes_reverse(self.dcel.vertex_representative(vertex))
                        .all(|edge| half_edges_counter.is_inner_edge(edge))
                })
                // PERF: Needless collect?
                .collect::<Vec<VertexId>>(),
            &half_edges_counter
                .outer_edges(&self.dcel)
                .map(|(half_edge, _)| half_edge)
                .collect::<Vec<HalfEdgeId>>(),
        );
    }

    fn absorb_faces_over_edges_and_vertices_in_perimeter(
        &mut self,
        absorbing_face: FaceId,
        faces_to_absorb: impl IntoIterator<Item = FaceId>,
        edges_to_remove: impl IntoIterator<Item = EdgeId>,
        vertices_to_remove: impl IntoIterator<Item = VertexId>,
        perimeter_half_edges: &[HalfEdgeId],
    ) {
        let faces_to_absorb: Vec<FaceId> = faces_to_absorb.into_iter().collect();
        let edges_to_remove: Vec<EdgeId> = edges_to_remove.into_iter().collect();

        for &face_to_absorb in &faces_to_absorb {
            self.faces_rtree.remove(&GeomWithData::new(
                Self::rectangle_from_vertex_weights(
                    self.dcel
                        .face_vertices(face_to_absorb)
                        .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
                ),
                face_to_absorb,
            ));
        }

        for &edge_to_remove in &edges_to_remove {
            let endpoints = self.dcel.edge_endpoints(edge_to_remove);

            self.edges_rtree.remove(&GeomWithData::new(
                Rectangle::from_corners(
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.0).clone()),
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.1).clone()),
                ),
                edge_to_remove,
            ));
        }

        // Remove the absorbing face from the faces R-tree before its shape
        // changes, which would otherwise invalidate its bbox and make it
        // impossible to access anymore.
        self.faces_rtree.remove(&GeomWithData::new(
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(absorbing_face)
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
            absorbing_face,
        ));

        self.dcel.absorb_faces_over_edges_and_vertices_in_perimeter(
            absorbing_face,
            faces_to_absorb,
            edges_to_remove,
            vertices_to_remove,
            perimeter_half_edges,
        );

        // Insert the absorbing face back in the R-tree now that its bbox is
        // done changing.
        self.add_face_to_rtree(absorbing_face);
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone + Default,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn split_edge_by_vertex(
        &mut self,
        edge_to_split: EdgeId,
        vertex: VW,
    ) -> (VertexId, EdgeId) {
        let original_endpoints = self.dcel.edge_endpoints(edge_to_split);

        // Remove the edge to split from the edges R-tree before its shape
        // changes, which would otherwise invalidate its bbox and make it
        // impossible to access anymore.
        self.edges_rtree.remove(&GeomWithData::new(
            Rectangle::from_corners(
                Into::<P>::into(self.dcel.vertex_weight(original_endpoints.0).clone()),
                Into::<P>::into(self.dcel.vertex_weight(original_endpoints.1).clone()),
            ),
            edge_to_split,
        ));

        let (new_vertex, new_edge) = self.dcel.split_edge_by_vertex(edge_to_split, vertex);
        self.add_edges_to_rtree([new_edge]);

        // Insert the split edge back in the edges R-tree now that its bbox is
        // done changing.
        self.add_edges_to_rtree([edge_to_split]);

        (new_vertex, new_edge)
    }

    pub fn split_face_by_edge(
        &mut self,
        from: VertexId,
        to: VertexId,
        face_to_split: FaceId,
    ) -> ((HalfEdgeId, HalfEdgeId), FaceId) {
        let mut result = ((HalfEdgeId::new(0), HalfEdgeId::new(0)), FaceId::new(0));

        self.update_split_face_rtree(face_to_split, |dcel| {
            let (new_edge, new_face) = dcel.split_face_by_edge(from, to, face_to_split);
            let edge_id = EdgeId::new(new_edge.0, new_edge.1);
            result = (new_edge, new_face);
            (vec![edge_id], vec![new_face])
        });

        result
    }

    pub fn split_face_by_edge_chain(
        &mut self,
        from: VertexId,
        to: VertexId,
        vertex_weights: impl IntoIterator<Item = VW>,
        face_to_split: FaceId,
    ) -> (Vec<(HalfEdgeId, HalfEdgeId)>, FaceId) {
        let mut result = (vec![], FaceId::new(0));

        self.update_split_face_rtree(face_to_split, |dcel| {
            let (new_edges, new_face) =
                dcel.split_face_by_edge_chain(from, to, vertex_weights, face_to_split);
            let edge_ids = new_edges
                .iter()
                .map(|(forward, backward)| EdgeId::new(*forward, *backward))
                .collect::<Vec<_>>();
            result = (new_edges, new_face);
            (edge_ids, vec![new_face])
        });

        result
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn split_face_by_edge_chain_with_all_weights(
        &mut self,
        from: VertexId,
        to: VertexId,
        vertex_weights: impl IntoIterator<Item = VW>,
        edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
        face_to_split: FaceId,
        new_face_weight: FW,
    ) -> (Vec<(HalfEdgeId, HalfEdgeId)>, FaceId) {
        let mut result = (vec![], FaceId::new(0));

        self.update_split_face_rtree(face_to_split, |dcel| {
            let (new_edges, new_face) = dcel.split_face_by_edge_chain_with_all_weights(
                from,
                to,
                vertex_weights,
                edge_weights,
                face_to_split,
                new_face_weight,
            );
            let edge_ids = new_edges
                .iter()
                .map(|(forward, backward)| EdgeId::new(*forward, *backward))
                .collect::<Vec<_>>();
            result = (new_edges, new_face);
            (edge_ids, vec![new_face])
        });

        result
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone + Default,
    FW: Clone + Default,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    /// Partition a face into triangles by inserting a vertex inside and then
    /// adding edges between it and the original face's vertices.
    ///
    /// The original face is reused for the first triangle. New faces are
    /// created for all the other triangles.
    ///
    /// Returns the new vertex id together with the ids of all the newly created
    /// edges and faces (the reused already existing edges and face are not
    /// included).
    pub fn triangulate_face_around_point(
        &mut self,
        perimeter_face: FaceId,
        inner_vertex_weight: VW,
    ) -> (VertexId, Vec<EdgeId>, Vec<FaceId>) {
        let mut result = (VertexId::new(0), vec![], vec![]);

        self.update_split_face_rtree(perimeter_face, |dcel| {
            result = dcel.triangulate_face_around_point(perimeter_face, inner_vertex_weight);

            (result.1.clone(), result.2.clone())
        });

        result
    }

    pub fn fan_triangulate(
        &mut self,
        perimeter_face: FaceId,
        apex: VertexId,
    ) -> (Vec<EdgeId>, Vec<FaceId>) {
        self.update_split_face_rtree(perimeter_face, |dcel| {
            dcel.fan_triangulate(perimeter_face, apex)
        })
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW: Clone,
    FW: Clone,
    VC: Get<usize, Value = Vertex<VW>> + Insert<usize> + Push<usize>,
    HEC: Get<usize, Value = HalfEdge<HEW>> + Insert<usize> + Push<usize>,
    FC: Get<usize, Value = Face<FW>> + Insert<usize> + Push<usize>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    pub fn triangulate_face_around_point_with_all_weights(
        &mut self,
        perimeter_face: FaceId,
        inner_vertex_weight: VW,
        inner_edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
        triangle_face_weights: impl IntoIterator<Item = FW>,
    ) -> (VertexId, Vec<EdgeId>, Vec<FaceId>) {
        let mut result = (VertexId::new(0), vec![], vec![]);

        self.update_split_face_rtree(perimeter_face, |dcel| {
            result = dcel.triangulate_face_around_point_with_all_weights(
                perimeter_face,
                inner_vertex_weight,
                inner_edge_weights,
                triangle_face_weights,
            );
            (result.1.clone(), result.2.clone())
        });

        result
    }

    pub fn fan_triangulate_with_all_weights(
        &mut self,
        perimeter_face: FaceId,
        apex: VertexId,
        inner_edge_weights: impl IntoIterator<Item = (HEW, HEW)>,
        triangle_face_weights: impl IntoIterator<Item = FW>,
    ) -> (Vec<EdgeId>, Vec<FaceId>) {
        self.update_split_face_rtree(perimeter_face, |dcel| {
            dcel.fan_triangulate_with_all_weights(
                perimeter_face,
                apex,
                inner_edge_weights,
                triangle_face_weights,
            )
        })
    }
}

impl<
    P: Point,
    VW: Clone + Into<P>,
    HEW,
    FW,
    VC: Get<usize, Value = Vertex<VW>>,
    HEC: Get<usize, Value = HalfEdge<HEW>>,
    FC: Get<usize, Value = Face<FW>>,
    ER: Insert<GeomWithData<Rectangle<P>, EdgeId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, EdgeId>>,
    FR: Insert<GeomWithData<Rectangle<P>, FaceId>, Value = ()>
        + Remove<GeomWithData<Rectangle<P>, FaceId>>,
> RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    fn update_split_face_rtree<F>(
        &mut self,
        face_to_update: FaceId,
        mutate_fn: F,
    ) -> (Vec<EdgeId>, Vec<FaceId>)
    where
        F: FnOnce(&mut Dcel<VW, HEW, FW, VC, HEC, FC>) -> (Vec<EdgeId>, Vec<FaceId>),
    {
        // The updated face id is reused as one of the new faces. So, we need
        // to remove it from the faces R-tree before its shape changes, as it
        // would otherwise invalidate its bbox and make it impossible to access
        // anymore.
        self.faces_rtree.remove(&GeomWithData::new(
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(face_to_update)
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
            face_to_update,
        ));

        let (new_edges, new_faces) = mutate_fn(&mut self.dcel);

        // Insert the updated face back in the faces R-tree now that its bbox
        // is done changing.
        self.add_face_to_rtree(face_to_update);

        self.add_faces_to_rtree(new_faces.clone());
        self.add_edges_to_rtree(new_edges.clone());

        (new_edges, new_faces)
    }

    fn add_faces_to_rtree(&mut self, faces: impl IntoIterator<Item = FaceId>) {
        for face in faces {
            self.add_face_to_rtree(face);
        }
    }

    fn add_face_to_rtree(&mut self, face: FaceId) {
        self.add_face_to_rtree_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(face)
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
        )
    }

    fn add_face_to_rtree_with_bbox(&mut self, face: FaceId, face_rectangle: Rectangle<P>) {
        self.faces_rtree
            .insert(GeomWithData::new(face_rectangle, face), ());
    }

    fn add_face_with_edges_to_rtrees(&mut self, face: FaceId) {
        self.add_face_with_edges_to_rtrees_with_bbox(
            face,
            Self::rectangle_from_vertex_weights(
                self.dcel
                    .face_vertices(face)
                    .map(|vertex| self.dcel.vertex_weight(vertex).clone()),
            ),
        )
    }

    fn add_face_with_edges_to_rtrees_with_bbox(
        &mut self,
        face: FaceId,
        face_rectangle: Rectangle<P>,
    ) {
        self.add_face_to_rtree_with_bbox(face, face_rectangle);
        self.add_edges_to_rtree(self.dcel.face_edges(face).collect::<Vec<EdgeId>>());
    }

    fn add_edges_to_rtree(&mut self, edges: impl IntoIterator<Item = EdgeId>) {
        for edge in edges {
            self.add_edge_to_rtree(edge);
        }
    }

    fn add_edge_to_rtree(&mut self, edge: EdgeId) {
        let endpoints = self.dcel.edge_endpoints(edge);
        self.edges_rtree.insert(
            GeomWithData::new(
                Rectangle::from_corners(
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.0).clone()),
                    Into::<P>::into(self.dcel.vertex_weight(endpoints.1).clone()),
                ),
                edge,
            ),
            (),
        );
    }
}

impl<P: Point, VW: Into<P>, HEW, FW, VC, HEC, FC, ER, FR>
    RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    fn rectangle_from_vertex_weights(weights: impl IntoIterator<Item = VW>) -> Rectangle<P> {
        Rectangle::from_aabb(weights.into_iter().fold(AABB::new_empty(), |aabb, weight| {
            aabb.merged(&AABB::from_point(Into::<P>::into(weight)))
        }))
    }
}

#[cfg(feature = "undoredo")]
impl<
    P: Point,
    VW: Clone,
    HEW: Clone,
    FW: Clone,
    VCD: Clone + KeyedCollection,
    VC: Clone + KeyedCollection + ApplyDelta<VCD>,
    HECD: Clone + KeyedCollection,
    HEC: Clone + KeyedCollection + ApplyDelta<HECD>,
    FCD: Clone + KeyedCollection,
    FC: Clone + KeyedCollection + ApplyDelta<FCD>,
    FRD: Clone + KeyedCollection,
    FR: Clone + KeyedCollection + ApplyDelta<FRD>,
    ERD: Clone + KeyedCollection,
    ER: Clone + KeyedCollection + ApplyDelta<ERD>,
> ApplyDelta<RTreedDcel<P, VW, HEW, FW, VCD, HECD, FCD, ERD, FRD>>
    for RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    fn apply_delta(&mut self, delta: &Delta<RTreedDcel<P, VW, HEW, FW, VCD, HECD, FCD, ERD, FRD>>) {
        let (removed, inserted) = delta.clone().dissolve();

        let dcel_delta = Delta::with_removed_inserted(removed.dcel, inserted.dcel);
        self.dcel.apply_delta(&dcel_delta);

        let edges_rtree_delta =
            Delta::with_removed_inserted(removed.edges_rtree, inserted.edges_rtree);
        self.edges_rtree.apply_delta(&edges_rtree_delta);

        let faces_rtree_delta =
            Delta::with_removed_inserted(removed.faces_rtree, inserted.faces_rtree);
        self.faces_rtree.apply_delta(&faces_rtree_delta);
    }
}

#[cfg(feature = "undoredo")]
impl<
    P: Point,
    VW: Clone,
    HEW: Clone,
    FW: Clone,
    VCD: Clone + KeyedCollection,
    VC: Clone + KeyedCollection + FlushDelta<VCD>,
    HECD: Clone + KeyedCollection,
    HEC: Clone + KeyedCollection + FlushDelta<HECD>,
    FCD: Clone + KeyedCollection,
    FC: Clone + KeyedCollection + FlushDelta<FCD>,
    FRD: Clone + KeyedCollection,
    FR: Clone + KeyedCollection + FlushDelta<FRD>,
    ERD: Clone + KeyedCollection,
    ER: Clone + KeyedCollection + FlushDelta<ERD>,
> FlushDelta<RTreedDcel<P, VW, HEW, FW, VCD, HECD, FCD, ERD, FRD>>
    for RTreedDcel<P, VW, HEW, FW, VC, HEC, FC, ER, FR>
{
    fn flush_delta(&mut self) -> Delta<RTreedDcel<P, VW, HEW, FW, VCD, HECD, FCD, ERD, FRD>> {
        let (removed_dcel, inserted_dcel) = self.dcel.flush_delta().dissolve();
        let (removed_edges_rtree, inserted_edges_rtree) = self.edges_rtree.flush_delta().dissolve();
        let (removed_faces_rtree, inserted_faces_rtree) = self.faces_rtree.flush_delta().dissolve();

        Delta::with_removed_inserted(
            RTreedDcel {
                dcel: removed_dcel,
                edges_rtree: removed_edges_rtree,
                faces_rtree: removed_faces_rtree,
                point_marker: PhantomData,
            },
            RTreedDcel {
                dcel: inserted_dcel,
                edges_rtree: inserted_edges_rtree,
                faces_rtree: inserted_faces_rtree,
                point_marker: PhantomData,
            },
        )
    }
}

#[cfg(all(test, feature = "rstar", feature = "stable-vec"))]
mod test {
    use rstar::{
        RTreeObject,
        primitives::{GeomWithData, Rectangle},
    };

    use crate::{
        EdgeId, FaceId, HalfEdgeId, RTreedStableDcel, VertexId, assert_face_boundary,
        init_dcel_with_3x3_hex_mesh,
    };

    #[test]
    fn test_insert_edge() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        let face_to_split = FaceId::new(5);
        let face_to_split_vertices: Vec<VertexId> =
            rtreed_dcel.dcel.face_vertices(face_to_split).collect();

        // Split face 5 in two with a single edge.
        let (_new_edge, _new_face) =
            rtreed_dcel.insert_edge(face_to_split_vertices[0], face_to_split_vertices[3]);

        // There are now eleven faces in total: one unbounded and ten bounded.
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 78);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 11);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 39);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 10);

        // The original hexagon is now split into two quads.
        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        // First quad face.
        assert_face_boundary!(rtreed_dcel.dcel, 5, 4);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);
        // Second quad face.
        assert_face_boundary!(&rtreed_dcel.dcel, 10, 4);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
        assert_face_bbox_validity(&rtreed_dcel, 10);
    }

    #[test]
    fn test_insert_polygon() {
        let mut rtreed_dcel = RTreedStableDcel::<(i32, i32)>::new();
        let face = rtreed_dcel.insert_polygon([(0, 0), (10, 0), (10, 10), (0, 10)]);

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 4);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 8);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 2);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 4);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 1);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, face.id(), 4);

        assert_face_bbox_validity(&rtreed_dcel, face.id());
    }

    // TODO: Test remove face.

    // TODO: Test remove edge.

    #[test]
    fn test_merge_faces_around_vertex() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.merge_faces_around_vertex(VertexId::new(8));

        // There are now eight faces in total: one unbounded and seven bounded.
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 29);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 70);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 8);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 35);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 7);

        // Among the remaining faces, one is now a dodecagon, and the remaining
        // six are hexagons.
        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        // Face 2 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        // Face 4 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 5, 12);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        // Face 2 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 3);
        // Face 4 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn test_absorb_faces_around_vertex() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.absorb_faces_around_vertex(FaceId::new(2), VertexId::new(8));

        // There are now eight faces in total: one unbounded and seven bounded.
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 29);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 70);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 8);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 35);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 7);

        // Among the remaining faces, one is now a dodecagon, and the remaining
        // six are hexagons.
        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 12);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        // Face 4 does not exist.
        // Face 5 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        // Face 4 does not exist.
        // Face 5 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn test_merge_faces_over_edge() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.merge_faces_over_edges_and_vertices(
            [FaceId::new(5), FaceId::new(6)],
            [EdgeId::new(HalfEdgeId::new(44), HalfEdgeId::new(45))],
            [],
        );

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 74);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 9);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 37);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 8);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 5, 10);
        // Face 6 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        // Face 6 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn test_absorb_face_into_face_over_edge() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.absorb_faces_over_edges_and_vertices(
            FaceId::new(6),
            [FaceId::new(5)],
            [EdgeId::new(HalfEdgeId::new(44), HalfEdgeId::new(45))],
            [],
        );

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 74);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 9);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 37);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 8);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        // Face 5 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 6, 10);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        // Face 5 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn merge_two_faces() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.merge_faces([FaceId::new(7), FaceId::new(8)]);

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 74);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 9);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 37);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 8);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 5, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 10);
        // Face 8 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        // Face 8 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn absorb_face_into_face() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        rtreed_dcel.absorb_faces(FaceId::new(8), [FaceId::new(7)]);

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 74);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 9);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 37);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 8);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 5, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        // Face 7 does not exist.
        assert_face_boundary!(rtreed_dcel.dcel, 8, 10);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        assert_face_bbox_validity(&rtreed_dcel, 6);
        // Face 7 does not exist.
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
    }

    #[test]
    fn test_split_face_by_edge() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        let face_to_split = FaceId::new(5);
        let face_to_split_vertices: Vec<VertexId> =
            rtreed_dcel.dcel.face_vertices(face_to_split).collect();

        // Split face 5 in two with a single edge.
        let (_new_edge, _new_face) = rtreed_dcel.split_face_by_edge(
            face_to_split_vertices[0],
            face_to_split_vertices[3],
            face_to_split,
        );

        // There are now eleven faces in total: one unbounded and ten bounded.
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 30);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 78);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 11);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 39);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 10);

        // The original hexagon is now split into two quads.
        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        // Second quad face.
        assert_face_boundary!(rtreed_dcel.dcel, 5, 4);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);
        // Second quad face.
        assert_face_boundary!(&rtreed_dcel.dcel, 10, 4);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
        assert_face_bbox_validity(&rtreed_dcel, 10);
    }

    #[test]
    fn test_split_edge_by_vertex() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        let edge = EdgeId::new(HalfEdgeId::new(38), HalfEdgeId::new(39));
        let original_endpoints = rtreed_dcel.dcel.edge_endpoints(edge);
        let original_edge_count = rtreed_dcel.edges_rtree.as_ref().size();

        let (_new_vertex, new_edge) = rtreed_dcel.split_edge_by_vertex(edge, (259, 150));

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 31);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 78);
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 31);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 78);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 10);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 39);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 9);

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 7);
        assert_face_boundary!(rtreed_dcel.dcel, 5, 7);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);

        assert_edge_bbox_validity(&rtreed_dcel, edge);
        assert_edge_bbox_validity(&rtreed_dcel, new_edge);
    }

    #[test]
    fn test_split_face_by_chain_of_two_edges() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        let face_to_split = FaceId::new(5);
        let face_to_split_vertices: Vec<VertexId> =
            rtreed_dcel.dcel.face_vertices(face_to_split).collect();

        // Split face 5 in two with a chain of two edges, with their common
        // point around the face's center.
        let (_new_edges, _new_face) = rtreed_dcel.split_face_by_edge_chain(
            face_to_split_vertices[0],
            face_to_split_vertices[3],
            [(259, 150)],
            face_to_split,
        );

        // There are now eleven faces in total: one unbounded and ten bounded.
        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 31);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 80);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 11);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 40);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 10);
        // TODO: Test number of elements here and further below.

        // The original hexagon is now split into two pentagons.

        assert_face_boundary!(rtreed_dcel.dcel, 0, 0);
        assert_face_boundary!(rtreed_dcel.dcel, 1, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 2, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 3, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 4, 6);
        // First pentagon face.
        assert_face_boundary!(rtreed_dcel.dcel, 5, 5);
        assert_face_boundary!(rtreed_dcel.dcel, 6, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 7, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 8, 6);
        assert_face_boundary!(rtreed_dcel.dcel, 9, 6);
        // Second pentagon face.
        assert_face_boundary!(rtreed_dcel.dcel, 10, 5);

        assert_face_bbox_validity(&rtreed_dcel, 1);
        assert_face_bbox_validity(&rtreed_dcel, 2);
        assert_face_bbox_validity(&rtreed_dcel, 3);
        assert_face_bbox_validity(&rtreed_dcel, 4);
        assert_face_bbox_validity(&rtreed_dcel, 5);
        assert_face_bbox_validity(&rtreed_dcel, 6);
        assert_face_bbox_validity(&rtreed_dcel, 7);
        assert_face_bbox_validity(&rtreed_dcel, 8);
        assert_face_bbox_validity(&rtreed_dcel, 9);
        assert_face_bbox_validity(&rtreed_dcel, 10);
    }

    #[test]
    fn test_triangulate_face_around_point() {
        let mut rtreed_dcel = init_dcel_with_3x3_hex_mesh!(RTreedStableDcel<(i32, i32)>);
        let (_inner_vertex, new_edges, new_faces) =
            rtreed_dcel.triangulate_face_around_point(FaceId::new(5), (260, 125));

        assert_eq!(new_faces.len(), 5);
        assert_eq!(new_edges.len(), 6);

        assert_eq!(rtreed_dcel.dcel.vertices().num_elements(), 31);
        assert_eq!(rtreed_dcel.dcel.half_edges().num_elements(), 88);
        assert_eq!(rtreed_dcel.dcel.faces().num_elements(), 15);
        assert_eq!(rtreed_dcel.edges_rtree.as_ref().size(), 44);
        assert_eq!(rtreed_dcel.faces_rtree.as_ref().size(), 14);

        assert_face_boundary!(rtreed_dcel.dcel, 5, 3);

        // All new faces are triangles.
        for face in &new_faces {
            let face_vertices: Vec<VertexId> = rtreed_dcel.dcel.face_vertices(*face).collect();
            assert_eq!(face_vertices.len(), 3);
        }

        // The new edges are in the edges rtree.
        assert_eq!(
            rtreed_dcel
                .edges_rtree
                .as_ref()
                .iter()
                .filter(|e| new_edges.contains(&e.data))
                .count(),
            6
        );

        assert_face_bbox_validity(&rtreed_dcel, 5);
        for face in &new_faces {
            assert!(
                rtreed_dcel
                    .faces_rtree
                    .as_ref()
                    .locate_in_envelope(
                        &RTreedStableDcel::<(i32, i32)>::rectangle_from_vertex_weights(
                            rtreed_dcel
                                .dcel
                                .face_vertices(*face)
                                .map(|vertex| rtreed_dcel.dcel.vertex_weight(vertex).clone()),
                        )
                        .envelope(),
                    )
                    .any(|element| element.data == *face)
            );
        }
    }

    // TODO: Test fan triangulation.

    fn assert_face_bbox_validity(rtreed_dcel: &RTreedStableDcel<(i32, i32)>, face: usize) {
        let rectangle = face_rectangle(rtreed_dcel, FaceId::new(face));
        assert!(
            rtreed_dcel
                .faces_rtree
                .as_ref()
                .locate_in_envelope(&rectangle.envelope())
                .any(|&element| element == GeomWithData::new(rectangle, FaceId::new(face),))
        );
    }

    fn assert_edge_bbox_validity(rtreed_dcel: &RTreedStableDcel<(i32, i32)>, edge: EdgeId) {
        let endpoints = rtreed_dcel.dcel.edge_endpoints(edge);
        let rectangle = Rectangle::from_corners(
            rtreed_dcel.dcel.vertex_weight(endpoints.0).clone(),
            rtreed_dcel.dcel.vertex_weight(endpoints.1).clone(),
        );

        assert!(
            rtreed_dcel
                .edges_rtree
                .as_ref()
                .locate_in_envelope(&rectangle.envelope())
                .any(|&element| element == GeomWithData::new(rectangle, edge))
        );
    }

    fn face_rectangle(
        rtreed_dcel: &RTreedStableDcel<(i32, i32)>,
        face: FaceId,
    ) -> Rectangle<(i32, i32)> {
        RTreedStableDcel::<(i32, i32)>::rectangle_from_vertex_weights(
            rtreed_dcel
                .dcel
                .face_vertices(face)
                .map(|vertex| rtreed_dcel.dcel.vertex_weight(vertex).clone()),
        )
    }
}