Formally Verifying a Vertical Cell Decomposition Algorithm

For the proofs, we assume that we work with a type R of numbers with the structure of a real field, as understood by the Mathematical Components library. Among other operations, this mathematical structure provides a boolean test for equality between numbers and boolean test for comparison. The library then provides a lot of tools with their properties, like a boolean test for list membership, sorting, etc. In practice, this means that the field of rational numbers is suitable for all computations performed in this development.

Points are then described as pairs of numbers in R, and edges are described as pairs of points. These data types inherit the equality tests and the polymorphic functions for list membership of the Mathematical Components library can again be used for these data types. One peculiarity of our development is that we chose to describe the edges as pairs of points with an extra condition in our proofs: the first point has to have a first coordinate that is strictly smaller than the first coordinate of the second point. This is an example of design where some data invariant is hard-coded in the data structure. As a result, the algorithm cannot process vertical obstacles, because the data structure representing the obstacles simply does not accept such cases. In the last part of the paper, we discuss ways to circumvent this limitation concerning vertical edges.

The definition we use to introduce edges in our development is as follows

Record edge := Bedge {left_pt : pt; right_pt : pt;
    _ : left_pt.x < right_pt.x}.

This definition does not give a name to the condition that the left point has a smaller first coordinate. Such a name is given later in our development.

The most important basic block used in the algorithm is the computation of whether a point p is above or below an edge, or whether it is aligned with the two extremities of the edge. To compute this condition, it is well known in computational geometry [6] that the following determinant can be used, where l and r are the names for the left and right extremities of the edge:

This is a determinant computation, and it is well known that this determinant computes twice the signed area of the triangle lrp. The function area computes this determinant. Because it is a signed area, the value is positive if p is above and negative if p is below. We use the notations $p$ <<= $g$ and $p$ <<< $g$ to express that a point $p$ is in the half-plane below an edge $g$ or strictly below an edge, respectively.

Because the algorithm uses a vertical sweeping line, our proofs require that we check when a point is in the vertical cylinder containing the edge. This is simply done by comparing the point’s first coordinate with the first coordinates of the edge’s extremities. We write valid_edge $g$ $p$ to express that $p$ is in that cylinder. We use the function point_edge, with the notation p === g, to express that point p is included in the obstacle g.

When a point is valid for an edge, the vertical projection of that point on that edge exists. It is computed by a function called vertical_projection. Computing the vertical projection of a point on an oblique edge requires that one computes a division. This is the one place in the whole algorithm where division is used. The function vertical_projection takes a point and an edge as arguments and returns an option type: the projection is only computed if the point is in the vertical cylinder given by the edge. We prove a few properties of this projection function: the point it returns is part of the edge on which the projection happens, two points with the same second coordinate are projected to the same point, and the input point and its projection are vertically aligned.

We define a relation between edges, denoted <| and called edge_below which expresses that the first edge is below the second one. Edge $g_1$ is below Edge $g_2$ if both extremities of $g_1$ are in the half plane below $g_2$ or if both extremities of $g_2$ are in the half plane above $g_1$. During the operation of the algorithm, we use this relation to sort edges going out of a given event. For edges that all have the same left extremity, the <| relation is transitive, but in general this relation is not transitive.

Each cell has a high edge, a low edge, a left side, and a right side. For the lateral sides, we record the ordered sequence of points that correspond to events in contact with this cell. These points are unsafe. These two sequences of points are called the left points and right points of the cell (record projections named left_pts and right_pts in the code). These sequences are ordered and contain no duplication, so that the open interval between two points in a sequence is a door, a safe passage to another cell.

In the finished closed cells, the side point sequences are not empty, as they always contain the projections of some event on the low and high edge, which may be the same point.

As an intermediate data structure, we use cells where the right-hand side is empty, called open cells, but in the code it is never ambiguous whether some cell is open or closed.

The left points sequence of a cell is intended to contain points that are vertically aligned, the first point is on the high edge and the last point is on the low edge. The same properties hold for the right points. We defined a predicate named open_cell_side_limit_ok to describe the properties of the left points, and a predicate named closed_cell_side_limit_ok to describe the conjunction of the properties for the left and the right points.

We define the left limit of a cell to be the first coordinate of the first left point, and similarly for the right limit.

The safe part of a cell is defined as the union of its strict interior (points whose first coordinate is strictly between the left limit and the right limit, strictly above the low edge, and strictly below the high edge) and its doors (points on each of the sides that are strictly above the low edge and strictly below the high edge and different from the side points). In what follows, the safe part of a cell is described by a predicate named safe_part.

When the next event to be processed is not vertically aligned with the previous one, we need to scan the current sequence of open cells to know which of these cells will be affected by the current event. These cells simply are the cells for which the current event is below the high edge and above the low edge. We call these cells the contact cells. These cells are a sub-sequence of the current sequence of open cells. The function that computes this sequence actually returns several meaningful pieces of data:

1. 1.

   the prefix of unaffected cells,

2. 2.

   the sequence of contact cells without the last one,

3. 3.

   the last contact cell,

4. 4.

   the suffix of unaffected cells

5. 5.

   the low edge of the first contact cell,

6. 6.

   the high edge of the last contact cell

The function that computes this decomposition is called open_cells_decomposition. It takes as input a point and a sequence of cells. We proved that under some assumptions, the sequence of cells given as input is the concatenation of the prefix, the contact cells, and the suffix, and that the two edges given as output are indeed boundary edges of the first and last contact cells.

All contact cells are transformed into closed cells, by adding a right-side. When there is more than one contact cell, the first one has as its low edge an edge that is not in contact with the current event and as its high edge an edge whose right point is the current event. Similarly, the last contact cell has as its high edge an edge that is not in contact with the current event as its low edge an edge whose right point is the current event.

When there is only one contact cell, this is because the processed event is below that cell’s high edge and above that cell’s low edge. The sequence of the right points for the new cell only contains three points, the event and its projections to the cell’s two edges.

All contact cells are removed from the list of open cells, but new open cells, whose left side contains the currently processed event are added in their place. The new open cells are obtained by processing recursively the sequence of outgoing edges, sorted with respect to the <| relation (also known as edge_below).

In the second case, we know that the last open cell is not among the contact cells. We call the function open_cell_decomposition, but we only need to scan the suffix of the sequence of open cells (the cells that are above the last open cell). The rest of processing is the same as in the first case, paying attention to the fact that the prefix and the previous last open cell have to be included in the resulting sequence of open cells.

In the third case, the event is below the high edge of the last open cell, and this event cannot be the right extremity of an active edge, because there is no edge between the low and the high edge of the last open cell. For this reason, it is not necessary to re-run the open_cells_decomposition function, as it would simply return the last open cell as the only contact cell. Closing this contact cell to produce a new closed cell is not suitable, because it would create an empty closed cell (one where the lateral sides would be on the same vertical line).

No new closed cell is created, but the last closed cell needs to be updated to the current event as an extra unsafe point on the right side. This requirement explains why the scan state has a specific field for the last closed cell.

New open cells need to be created. If the current event has no outgoing edges, then it is enough to update the last open cell by adding the current event among its unsafe points on the left side. If the current event has outgoing edges, new open cells are created as in the first two cases, but the sequence of unsafe points on the left of the last open cell needs to be included in the sequence of left points of the first newly created open cell.

In the fourth case, the processed event is on the high edge of the last open cell. We know that the last open cell will be the first of contact cells, and we need to compute the other contact cells, by calling open_cells_decomposition, but this time with a sequence of cells that starts with the last open cell.

When closing cells, the first contact cell does not need to be closed and can simply discarded, because its left side is on same vertical line as its right side, and all the unsafe points were already taken into account by existing closed cells.

When creating new open cells, special attention is required for the case where the event has no outgoing edges. In this case, the last open cell needs to be replaced with an open cell that has as low edge the same low edge as the the last open cell, as high edge the high edge of the last contact cell (as computed by open_cells_decomposition) and as sequence of last points the sequence of left points from the last open cell where the projection on the high edge has been added.

The safety property will be expressed in the following manner: any point in the safe part of a closed cell is distinct from the obstacles and the events given as input. We already defined the safe part of closed cells in the section about cells.

In a similar manner, we can describe unsafe locations as the union of points lying on an obstacle, as described by the function point_on_edge and the set of orphaned events. In practice, we can define a predicate unsafe that takes as input a sequence of events and a point to describe this space.

We prove that the intersection between unsafe locations and the union of all safe parts of the resulting closed cells is empty.

However, we have to be careful because this proof is done under a collection of assumptions.

• $\mathrm{■}$

  The type of numbers provided to describe the point coordinates, the edges, alignment of points, and so on, with its operation and comparison predicates forms an ordered field structure. It means that the algorithm only uses addition, multiplication, and their inverses.

• $\mathrm{■}$

  The only allowed intersections between obstacles are at their extremities. Removing this constraint is planned as future work.

• $\mathrm{■}$

  The sequence of obstacles does not contain any duplications.

• $\mathrm{■}$

  There is a bounding box, given by a bottom and a top edge and two lateral sides, which contains all events. In particular, this bounding box has distinct lateral sides and bottom and top edges (the bottom edge being below the top edge).

When propagating this specification to the two phases, it turns out we need the following specification for the sorting phase.

1. 1.

   The union of all outgoing edges of all events in the sorted sequence of events contains the initial obstacles.

2. 2.

   For every event, all outgoing edges have the same left extremity, located at this event.

3. 3.

   For every event, the sequence of outgoing edges has no duplicates.

4. 4.

   The sequence of events is strictly sorted lexicographically.

5. 5.

   No event in the sequence appears in the middle of one of the obstacles.

6. 6.

   The right extremity of every outgoing edge is also located at an event existing in the sequence of events.

7. 7.

   There are no intersections between pairs of edges except at their extremities.

8. 8.

   All events are inside the bounding box.

The second part of the program, which processes the sequence of events does not make the assumption that events are always edge extremities. When specifying the safety property for the results of the algorithm, we simply need to make sure that these orphaned events are excluded from the safe subset of the configuration space.

A property that does not appear immediately as a safety property is that the middle of every cell is strictly included in the cell. This property is useful for users of this algorithm who wish to use a cell as a maneuvering space to move from one door of the cell to another door of the same cell, when both doors are on the same side. This property is the main reason for having a specific treatment for degenerate cases. This specific treatment actually guarantees that every closed cell has a left-hand side that is further to the left than the right-hand side.

To make sure that every cell contains at least a safe point, we actually need the property that every event has a list of outgoing edges without duplication. Our code to generate the sequence of events does not include steps to guarantee this, but we guarantee that it operates correctly if the input list of obstacles has no duplicates.

Finally, the statement we prove for the second phase is the following one, as written in the input language of the Rocq prover (by taste, we prefer to avoid special characters). In this statement, the function complete_process represents the algorithm that processes the sequence of events: it combines together the operation of constructing the first cells from the first event, and processing the last open cell to create the rightmost closed cell. The function events_to_edges collects all the edges from a sequence of events.

Lemma second_phase_safety (bottom top : edge) (evs : seq event) :
  well_formed_bounding_box bottom top ->
  {in bottom :: top ::
       events_to_edges evs &, forall e1 e2, inter_at_ext e1 e2} ->
  all (inside_box bottom top) [seq point e | e <- evs] ->
  sorted (@lexPtEv _) evs ->
  {in evs, forall ev, out_left_event ev} ->
  close_edges_from_events evs ->
  {in events_to_edges evs & evs, forall g e, non_inner g (point e)} ->
  {in evs, forall e, uniq (outgoing e)} ->
  {in complete_process bottom top evs, forall c,
    strict_inside_closed (cell_center c) c /\
    closed_cell_side_limit_ok c /\
    forall p, cell_safe_part c p -> ~ unsafe evs p}.

This statement has 8 prerequisites, mostly expressed as properties of the sequence of events. Among notations, {in F, forall x, …} represents a universal quantification over all elements satisfying F, {in F & G, forall x y, …} represents two quantifications over elements in F (for x) and G (for y), {in F &, forall x y, …} represents two quantifications both ranging over elements in F.

A more compact statement can be used if we consider the combination of the two phases, because the first phase guarantees the prerequisites on the sequence of events. The combination of the two phases is called edges_to_cells in this statement. While the concision of this statement is pleasant, it misses some of the power of the second phase.

Lemma two_phase_safety bottom top (edges : seq edge) :
  well_formed_bounding_box bottom top -> uniq edges ->
  {in bottom :: top :: edges &, forall e1 e2, inter_at_ext e1 e2} ->
  {in edges, forall g p, p === g -> inside_box bottom top p} ->
  {in edges_to_cells bottom top edges, forall c,
    strict_inside_closed (cell_center c) c /\
    closed_cell_side_limit_ok c /\
    forall p, cell_safe_part c p -> {in edges, forall g, ~ p === g}}.

This lemma only has 4 prerequisites: the bounding box given by the bottom and top edges must define a workable area of the plane, the sequence of input edges must have no duplication, considered edges may only meet at endpoints, and all edges must be inside the bounding box. With these four prerequisites, the algorithms of the first phase are able to produce a sequence of events satisfying the eight prerequisites of the second phase.

As a result, the algorithm produces a sequence of cells that satisfy three properties: each contains at least one point, they are well formed (the sequence of points on the side are vertical and ordered, the extremities of these sequences are on the low and high edges of the cell) and the safe part has an empty intersection with the union of the obstacles.

The infrastructure of definitions provided by the Mathematical Components library is not fully compatible with the extraction process provided by the Rocq prover [8]. To obtain code that is amenable for extraction, we write a first file called generic_trajectories, that does not rely on any of the mathematical structures of the library, but is parameterized by a given type supposed to represent numbers, with 4 operations (addition, subtraction, multiplication, division) and two comparisons. Also, this algorithm relies on a degraded type of edges, where the condition that the first point is on the left of the second point is absent from the data-structure.

The generic algorithm is used in all the proof files. In each of these proof files, we actually instantiate the type of numbers with a mathematical structure from the Mathematical Components library. Following the library’s idiom, we do not choose a specific field structure, we assume we are using one that exists. So, the proofs we make are assuming that the user of the code also respects that part of the specification and actually instantiates the software with a number structure that respects the properties specified by the real field structure.

We do provide a running implementation by instantiating the type of numbers with rational numbers, implemented as a pair of a numerator (an integer) and a denominator (a positive integer).

Every function imported from generic_trajectories is a function with at least 7 arguments. To make the whole setup comfortable, a header of notation definition is provided to instantiate the function to the type and operations provided by the field structure. This header of notations needs to be repeated in each working file, because the field structure named R in each file is conceptually a different structure for each file. However, theorems proved in one file can be reused in another simply because they are generalized at the end of each section and can be re-instantiated at each usage.

The first part of the algorithm consists in preparing the list of events so that processing these events in the order given by the list simulates the intuitive behavior of moving a vertical line from the left-hand side of the box to the right-hand side. Proving that this sorting algorithm guarantees the properties required by Lemma second_phase_safety is comparatively easy. However, this lemma only expresses safety with respect to the obstacles and points found in the list of events. Therefore, an extra property that is important is that the edges present in the resulting list of events are the input edges.

We will now concentrate on the second phase. We established six levels of invariants:

1. 1.

   In the first level, we state properties that concern the sequence of open cells and the relations between their edges and the remaining events (in the formal development this invariant is named open_cells_invariant).

2. 2.

   In the second level, we add properties that fall in three categories:

   1. (a)

      The cache consistency properties of the scan state: the scan state contains a last open cell, a last high edge, and a last x coordinate. It is assumed that the last high edge is the high edge of the last open cell, and that the last x coordinate is the left limit of the last open cell.

   2. (b)

      The expected properties of the remaining sequence of events: these events are ordered lexicographically, all the right extremities of their edges are also covered by events in the sequence, the field outgoing of each event is a list that only contains edges whose left-hand extremity is on this event, and this list has no duplicates.

   3. (c)

      There is no intersection between the high edges of all open cells and the outgoing edges of events in the sequence of events.

   4. (d)

      The bottom left corner of the last open cell is lexicographically smaller than any event in the sequence.

   5. (e)

      If the sequence of events is not empty, the next event is lexicographically larger than the left extremity of all edges of open cells, and lexicographically smaller than the right extremity.

   This is common_invariant.

3. 3.

   In the third level, we add properties that are specific to the last created open cell. This is last_open_cell_invariant.

4. 4.

   In the fourth level, we add mainly the properties that concern the closed cells that were built so far and their interactions with the current open cells. The main purpose of the properties available in this invariant is to contribute to the proof that all closed cells are disjoint and non-empty. This is disjoint_invariant.

5. 5.

   In the fifth level, we add mainly the properties that concern edges that were outgoing from events that have already been processed. In particular, we concentrate on properties that are useful to show that all these edges are covered by the top side of existing closed cells and existing open cells. This is edge_cover_invariant.

6. 6.

   In the sixth level, we add mainly the properties that concern the points recorded in cell sides. The aim is to show that cell sides enumerate the points that are unsafe on the lateral sides of cells. This is cell_sides_invariant.

In the end, the safety property for the space inside a given cell is guaranteed by the main argument that if edges are covered by the high edges of all closed cells, and all closed cells are disjoint, then edges cannot appear in the interior of a given cell.

When addressing a complex algorithm like the one presented here, it is possible to develop a simpler proof of the main result, obtained under stronger assumptions concerning the inputs. We performed such a simpler proof, which expresses how the algorithm works in cases where no events are vertically aligned. This proof relies on simpler invariants and leaves some branches of the code unverified. Still it was useful to us a a means to explore the structure of the more complete proof. In our formal development, the simpler proof is found in files general_position.v and safe_cells.v. In this paper, we concentrate on the second, more complete proof. The main lemma for that simple proof is named start_yields_safe_cells in file safe_cells.v.

The signed area function is instrumental to detect when a point is above or below an edge. More precisely, it makes it possible to detect when a point is above the straight line that carries the obstacle.

This area function is given by a determinant formula that we learned from Knuth [6]. It is quite important that this is only computed using a determinant, as it only needs ring operations (no division), without ever mentioning the projection of the given point on the given straight line.

For proofs, however, it is useful to have an alternative point of view, where one expresses the property of being above using a comparison of the second coordinates of the given point with the projected point on the given line.

The edge comparison relation, called edge_below and noted <| is directly related with the area function. To detect whether an edge is below another, we simply verify whether the two extremities of one are in the same half-plane with respect to the other.

It turns out that this relation is reflexive, not transitive, and not antisymmetric. We spent quite some time trying to show that it is transitive on some well-chosen subsets, but this turned out to be unproductive. For this reason, it is not very useful (although it is true) to show that the sequence of high edges of all open cells is sorted for the edge_below relation. We finally chose to express that this sequence of edges satisfies a stronger invariant: an edge $g_1$ that appears before an edge $g_2$ in this sequence is below.

All open cells intersect with the sweeping vertical line. As a consequence, it is useful to have a collection of lemmas establishing correspondences between edges below one another and the second coordinate of the projection of a point on these edges. This is quite natural to express but one must be careful to only consider points that are in the intersection of vertical cylinders above the two edges. The lemmas have names like order_edges_viz_point or under_pvert_y, sometimes with a strict qualifier added.

As a way to assess usability, we integrated this program into a larger application that produces trajectories without intersections with the obstacles between points inside the bounding box, when possible. To work properly, this larger application relies on a few more properties that we have not included in our specification yet. Here are two examples of such properties:

• $\mathrm{■}$

  The set of created cells should cover the whole interior of the bounding box

• $\mathrm{■}$

  When there exists a safe passage between two adjacent cells, these cells actually have a door in common.

According to our tests, these properties hold for our algorithm, but they have not been proved formally. If these properties were not satisfied, the program would not compute trajectories between some points that can actually be connected by such a trajectory.

The literature of formal verification has few examples of verification for algorithms coming from computational geometry. Pichardie and Bertot [11] published formally verified proofs of convex hull algorithms in two dimensions. They describe the determinant that can be used to know whether a point is on the left or the right of an oriented segment. We re-use that determinant in our work, to detect when a point is above or below an obstacle, because our obstacles can naturally be viewed as oriented segments (from left to right). Convex hulls, have been studied formally again in the work of Brun, Dufourd, and Magaud [4], but this time the aim was to illustrate the use of a general data structure to describe sub-divisions of the plane, called hypermaps. Our work makes no use of hypermaps. Another formally verified work using convex hulls is the work of Immler on zonotopes [5].

Geometry and convex hulls also play a significant role in two recent formalization efforts focused on mathematical riddles. F. Maric developed a formal proof of a conjecture of Erdös, that every set of points $2^{m-2}+1$, with $m\ge 3$ contains a convex polygon with $m$ sides, but just for the value $m=6$. [10] and B. Subercazeaux et al. developed a similar proof for the solution of a similar problem, where the convex polygon has to be empty, again for the value $m=6$ [13]. Both developments are characteristic in the use human expertise to reduce the problem to a discretized problem that is then solved using an automatic SAT solver, thus relying on proofs that cannot be check by the human eye. The reduction is formally verified using a interactive theorem prover, in the former case Isabelle/HOL and in the latter case Lean.

Motion planning is an envisioned application domain for the work presented here. Motion planning has also been studied formally, by Rizaldi and their co-authors [12]. This experiment is quite remote from ours, as they rely on an interaction between the proof assistant Isabelle and another symbolic engine to obtain trajectories paired with a proof. The work described in this article is less complete since we do not show how to compute trajectories, but it also places itself at a different level of interaction: we provide a correctness proof for an algorithm that can run independently from the proof assistant that was used to prove its correctness. Anand and Knepper [1] also concentrate on robot trajectories, but they do so for a single maneuver, checking that all numeric computations combine to guarantee that the robot will reach the prescribed destination, while we are interested in the preliminary step of decomposing the configuration space into cells where movement will be safe).

As future work, we wish to work in three directions. The first direction is to refine the algorithm to make it more efficient: for instance, each cell should contain a reference to its neighbors, so that computing the neighbors can be done in constant time, instead of the current situation where all existing cells must be visited to find the ones that have a door in common. The second direction is to make the algorithm stronger, by accepting vertical obstacles and by accepting obstacles that cross. The third direction is to work on trajectory computation, either by proving formally the correctness of a prototype that we already designed using Bezier curves, or by designing a variety of planning algorithms relying on different kinds of curves.