Pierre Dubosclard

Stanislas Larnier

Hubert Konik

Ariane Herbulot

Michel Devy

Automatic method for visual grading of seed food products

This paper presents an automatic method for visual grading, designed to solve the industrial problem of evaluation of seed lots. The sample is thrown in bulk onto a tray placed in a chamber for acquiring color image. An image processing method had been developed to separate and characterize each seed. The approach adopted for the segmentation step is based on the use of marked point processes and active contour, leading to tackle the problem by a technique of energy minimization.

Introduction

In agriculture, the global grain harvest reached several billion tons each year. Seed producers exchange their crops at a price determined by the quality of their production. This assessment, called grading, is performed for each set on a representative sample. The diculty of this assessment is to fully characterize the sample. To do so, it is necessary to qualify each of its elements. Historically, this has been performed manually by an operator. This method is exposed to various problems and the results can vary from one operator to another.

Alpha MOS company [START_REF]Alpha MOS[END_REF] develops systems for quality control of food products. It proposes a visual sensory system to provide an alternative to human evaluation. The assessment should be simple to implement and at least as fast as the human evaluation. The seed samples are presented in bulk, without any arrangement, but they are spread over a tray in such way that there is no overlapping between the objects to avoid occlusion (Figure 1).

The evaluation by the instrument is composed of three steps. The rst step is the acquisition: the operator places the samples on the tray in the instrument and takes an image. Then the detection step consists in nding each object in the image, to nally classify them in dierent quality classes regarding several criteria (shape, color, spot) during the classication step.

The quality of the sample can then be deduced from the result of the classication. The detection step is the main diculty. It is necessary to develop a segmentation method to isolate each object under the following constraints:

the number of seeds is unknown ( an approximative estimation can be done); the objects have quite generical geometric and chromatic features; they are randomly placed, without arrangement and with no overlapping.

In Section 2, a state of the art around the visual grading problem is presented. Section 3 describes the data acquisition system. The notion of marked point processes is introduced in Section 4. The segmentation steps are detailed in Section 5. Numerical results are presented in Section 6.

State of the art

Several studies have been conducted on the cereal seeds grading. Augustin et al. [START_REF] Agustin | Automatic milled rice quality analysis[END_REF] focused on the quality control of grain of rice, regarding dierent criteria of shape and color. From these criteria, a classication method based on neural network is used to qualify each grain. This approach gives good results for the classication of complete, broken and colors defect rice grain. However this method is applied on images with separated grain. The segmentation issue is then simplied by an operator or a mechanical system (vibrating bowl or slot) to separate the grain in front of the camera.

Other studies have been conducted on the cereal segmentation topic, mainly on wheat and rice. Yao et al. [START_REF] Yao | An automatic segmentation algorithm for touching rice grains images[END_REF] and Faessel et al. [START_REF] Faessel | Touching grain kernels separation by gap-lling[END_REF] focused on detection and separation of rice grain. They both address the problem by working on a binary image obtained by a threshold to separate the objects from the background. Yao et al. [START_REF] Yao | An automatic segmentation algorithm for touching rice grains images[END_REF] then work on the contours and search the concaves angles to connect them two at a time in order to detect objects boundaries. Faessel et al. [START_REF] Faessel | Touching grain kernels separation by gap-lling[END_REF] used a mathematical morphology method on the binary image: a skeleton operation on the background. The open lines of the skeleton, without ending, are then combined under some constraints to obtain the objects boundaries. These two methods give good results on image of touching grain with low density of objects. The computation times are short, but these methods are not adapted for images with heaps and high density of seeds.

Acquisition system

The acquisitions are made in a cabin (Figure 2) which integrates a camera and a lighting system. This cabin oers stable and reproducible acquisition condition, independently from the external lighting. Some improvements have been made on the existing system available at Alpha MOS. The lighting system and the camera have been replaced by new material to improve the quality and the stability of the color image acquisition. The lighting source retained is composed of white LEDs. These LEDs have a continuous spectrum in the visible range and were chosen for their stability over time in term of luminous intensity. As LEDs are punctual sources, a diuser is placed downstream to ensure the lighting homogeneity in the acquisition area. The image acquisition is performed at a distance of 400 mm from the object plan by a CMOS mono sensor color camera of 5 megapixel with a 5 mm lens. The chosen camera was a Basler acA2500-14gc. It oers a resolution on the object plan around 6 pixels per millimeter, which is important for our application as the objects have a size of only few millimeters. The image acquisitions presented in this paper were obtained with this system.

Marked point processes

The notion of marked point processes has widely been used to represent stochastic phenomena such as waiting queue. More recently, this approach was used to extract objects in image processing with for example, the detection of roads [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] or to count trees [START_REF] Perrin | A marked point process model for tree crown extraction in plantations[END_REF] on satellite and aerial images.

Introduction to the marked point processes

Figure 1 presents an example of seeds. The chosen approach to modelize and extract the seeds is based on the marked point processes. Indeed, the seeds can be represented by a generical simple shape and there is no arrangement between them, they are randomly disposed. These objects are dened by their positions and their geometric attributes or marks. Let χ be the space of objects such as χ = P × M , with P the space of the position and M the space of the geometric attributes describing the object. A conguration of objects from χ, noted x, is a non-arranged list of objects:

x = {x 1 , . . . , x n }, n ∈ N, x i ∈ χ, i = 1, . . . , n.
The objects to detect can be approximated by an ellipse characterized by its marks, for example its orientation, its minor axis and its major axis (Figure 3). Fig. 4: Object x i and its crown F(x i ) (a), and their equivalent disretized (b). For each object x i of the conguration x, an energy U (x i ) composed of two terms is associated. The rst term is a data term noted U d (x i ), which represents the likelihood of the marked point process regarding the data (the image in our case). This term is dened by the data of the object itself. The second term is an a priori term, noted U p (x i ), which imposes condition on the overall conguration.

In the context of object detection inside an image, the aim is to nd the most likely object conguration. This research is based on the two energies terms that are dened in the next section.

Energy U(x)

The energy associated to an object x i , noted U (x i ), is the sum of the data term U d (x i ) and a priori term U p (x i ). The energy for the conguration x is then:

U (x) = U d (x) + γU p (x)
with γ a weight coecient which is determined empiricaly.

The term U d (x) takes into account the image data for each object of x. It is computed by using the Bhattacharyya distance, noted d B , dened in [START_REF] Descamps | Automatic amingo detection using a multiple birth and death process[END_REF]:

d B (x i , F(x i )) = (µ 1 -µ 2 ) 2 4 σ 2 1 + σ 2 2 - 1 2 log 2σ 1 σ 2 σ 2 1 + σ 2 2
with an object x i ∈ χ and F(x i ) the object crown (Figure 4), (µ 1 , σ 1 ) and (µ 2 , σ 2 ) respectively the means and the variances of the radiometric values of the object and its crown. The computation of this distance provides a criterion that highlights area with important contrast between the object and its crown. It also takes into account the homogeneity of the area. Finally, the term U d (x) is dened as follows:

U d (x) = xi∈x U d (x i ) = xi∈x Q(d B (x i , F(x i )))
with Q(d B ) ∈ [-1,1] a quality function which favorizes or penalizes the objects considering a given threshold d 0 :

Q(d B ) = (1 - d B d 0 ) if d B < d 0 , Q(d B ) = exp - d B -d 0 100 -1 if d B ≥ d 0 .
The objects having an important contrast with their crown (d b > d 0 ) are then favorized and their associated data energy is negative. The U p (x) term gives information on the a priori knowledge on the target conguration, like the interactions between the objects. In the context of seed segmentation, U p (x) is a repulsive term that penalizes the objects overlapping. For each object x i ∈ χ, U p (x i ) is the sum of repulsive strengths emitted by the objects in interaction with x i , that are in overlapping with x i . These repulsive strengths are computed by counting the number of pixels that belong to the object x i and to its neighbouring objects noted V(x i ):

U p (x) = xi∈x U p (x i ) with U p (x i ) = xj ∈V(xi) A(x i ∩ x j )
and A(x i ∩ x j ) the common area of x i and x j objects.

The presented segmentation method is inspired by the multiple Birth-and-Death algorithm described by Descombes et al. in [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF]. This approach involved the marked point processes in an optimization framework. But we adapt this approach to treat our segmentation topic by adding a detection step between the birth step and the death step. The dierence with the approach of Descombes et al. is that we do not consider the optimization on the entire conguration but on specic objects. We use this Birth-and-Death dynamic as a sampler, the optimization part is realized by an active contour method detailed later. First, the initialization of the method is presented, then the dierent steps of the method are described. Birth map The rst initialization step consists of computing an image that is named birth map. This image has the same size as the input image and it associates to every position p a probability B(p) that there is an object centered at this position.

Initialization

This image is computed in two steps. The rst step consists in a binarization of the input image to separate the objects from the background, pixels of the background are set to zero and the probability associated is null. The second step is the computation of the Euclidean distance to the contours. Orientation maps The orientation maps are a set of images that associate to every dierent probabilities on the possible orientations of an object in this point. To obtain them, mean lters which have a geometrical shape and dierent orientations are used on the binary image. The geometrical shape is the mean of the possible shapes which approximate the seeds.

For example, the ellipse had been chosen in the case of rice seeds. Figure 6 presents three orientation maps obtained from the input image in Figure 5.

Active contour

The objects created in the context of the marked point processes have only simple shapes like ellipses with a limited range of axes sizes. The computational time is the main reason. But to accurately detect every object in the image, we need to obtain the most precise boundaries.

As the objects to detect have a generical shape, we decided to use the method based on an active contour with a geometric shape prior proposed by Bresson et al. [START_REF] Bresson | A variational model for object segmentation using boundary information and shape prior driven by the mumfordshah functional[END_REF]. This method follows the well-known energy functionnal model of Chen et al. [START_REF] Chen | Using prior shapes in geometric active contours in a variational framework[END_REF] where the shape prior of Leventon et al. [START_REF] Leventon | Statistical shape inuence in geodesic active contours[END_REF] is integrated. Finally, to improve the robustness of the method, Bresson et al. add a region-based energy term based on the Mumford-Shah functionnal (Vese and Chan [START_REF] Chan | Active contours without edges[END_REF]). This method is then based on three complementary terms dedicated to shape, boundary and region inside the contour.

Birth-and-Death dynamic

The method adapted for the Birth-and-Death algorithm is composed of three steps that are done iteratively. We added a detection step to the original approach in order to detect the object boundary with more accuracy, but also to be able to segment the heaps progressively from their boundaries to their cores. Birth step The rst step consists of objects and is illustrated on Figure 7.

For each point p of the input image, if there is not already an object at this position, the birth probability δB(p) is computed from the birth map, with δ a regularization coecient that handles the approximative number of objects. If an object x i is created, its axes are randomly chosen and its orientation is obtained from the set of orientation maps. For a given number of angles, the set of orientation maps provides the probability that an object has this angle. Then the data energy U p (x i ) of the object x i is computed and the object is placed in the conguration x by sorting them regarding their data energy. Once all the image is scanned, the algorithm goes to the next step.

Detection step The aim of the detection step is to validate the objects of the conguration x. Figure 8 described this step. For each object x i of the conguration taking by their data energy classement order, their data energy is compared to a threshold. If their energy is inferior to this threshold, they may be correctly placed. Then the active contour approach detailed previously is use to validate this hypothesis on one hand, and to get an accurate boundary if the object is correct on the other hand. The result of the active contour is analyzed to determine dierent criteria like area or roundness. From these criteria the object is then validated or not. If the object is validated, the object is removed from the conguration x, the birth map is updated by aecting the probability to create an object in this area at zero. The input image is also updated by turning the pixels values of the correct object to zero, so the heaps can be progressively processed from their boundaries to their cores. If the object is not correct, we only remove it from the conguration x. Death step The death step consists of cleaning the conguration x and is illustrated on Figure 9. If some object have been validated during the detection step of the current iteration, we compute the new data term of all the objects from the updated input image. We then compute the a priori energy and the death probability of each object in the conguration:

D(x p ) = δa ϕ (x ) 1 + δa ϕ (x p )
x with a ϕ (x p ) = exp (-ϕU (x p )). The object is then remove from the conguration with the probability D(x p ). Figure 10 illustrates the behaviour of the detection step in the cases of bad starting contour (a) and a good one (b). The rst leads to a nal contour with a shape distant from an ellipse and is rejected. The second leads to a nal contour with a shape similar to an ellipse and is accepted.

Figure 11 and Figure 12 presents the nal segmentation. The green color represents the contours. For the greater part, rice seeds and oats are well detected. Some improvements could be made. Further development is under consideration especially in the algorithm parameters selection. This selection could be made thanks to an automatic learning on separated seed.

Conclusion

This paper proposes an approach to perform the visual quality control of cereal seeds samples. This operation called visual grading can be treated in three steps: acquisition, segmentation and classication. An acquisition system of color images has been created to collect the data. A new segmentation approach has been developed, based on the marked point processes. The Birth-an-Death dynamic has been modied to integrate a new detection step based on an active contour with a shape prior term. The results on rice seeds are promising.

Experiments with higher density and with other type of seeds (barley, pea, wheat) are in progress. Other tests like comparison with human operator and reproductibility on the same sample in dierent congurations are also underway.

In the future, some algorithm parameters will be automatically learnt on simple images with a representation sample of separated seeds. The shape parameters would be extracted from statistics on the binarized image. The integration of a 3D data acquisition system like stereovision with two cameras is under consideration. Despite the hardware cost, such data might be useful in particular to enrich the birth map but also to provide criteria for the classication stage.

Acknowledgement

This CIFRE thesis work was made possible thanks to the involment of Alpha MOS company.

Fig. 1 : 2 :

 12 Fig. 1: Image of wheat seeds in bulk. Fig. 2: Acquisition system.

Fig. 3 :

 3 Fig. 3: Position and marks of an ellipse.

Fig. 5 :

 5 Fig. 5: From left to right : input image and birth map corresponding.

Fig. 6 :

 6 Fig. 6: Example of three orientation maps.

Fig. 7 :

 7 Fig. 7: Birth step.

Fig. 8 :

 8 Fig. 8: Detection step.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Death step.

Fig. 11 :

 11 Fig. 11: Final segmentation on rice sample.

Fig. 12 :

 12 Fig. 12: Final segmentation on oats ple.


Guyslain Naves
email: guyslain.naves@lif.univ-mrs.fr

Arnaud Spiwack
email: arnaud@spiwack.net

Balancing lists: a proof pearl ⋆

Starting with an algorithm to turn lists into full trees which uses non-obvious invariants and partial functions, we progressively encode the invariants in the types of the data, removing most of the burden of a correctness proof. The invariants are encoded using non-uniform inductive types which parallel numerical representations in a style advertised by Okasaki, and a small amount of dependent types.

Introduction

Starting with a list lst, we want to turn it into a binary tree tr of the following form (in Ocaml): type α tree = Node of α tree * α * α tree Leaf

With the constraints that lst must be the infix traversal of tr and that tr must be full, in the sense that every level except the last are required to be completely filled. Such a function turns, in particular, sorted lists into balanced binary search trees.

There are a number of folklore algorithms to achieve this result in linear time. Here we consider one of these algorithms, presented in Section 2, which repeatedly pairs up trees of height h in a list to form a list of trees of height h + 1. Our interest in this algorithm sprouts from the fact that its correctness is not obvious; in particular the invariants are complex: the main loop operates on a list of length 2 k -1 whose elements are alternately of two distinct forms.

In Sections 3 and 4 we show refinements of the algorithm where the invariants are pushed into the types, leading to a complete and short proof of correctness in Coq.

We start by giving a first, simple, implementation of the balancing algorithm. The heart of the algorithm relies on using an alternating list of length 2 k -1, where odd-position elements are trees and even-position elements are labels, of type α (indices starting from 1). A full tree of height k can be decomposed into the first k -1 levels, containing 2 k-1 -1 internal nodes, and the kth level, which contains both nodes and leaves. Thus, the 2 k-1 -1 labels in the alternating list will be used to label the internal nodes in the k -1 first levels of the balanced tree, while the 2 k-1 trees, all of height at most one at first, will constitute the level k.

Though we could encode labels as trees of height one in the alternating list, we rather use an appropriate type for the sake of readability:

type α tree_or_elt = Elt of α Tree of α tree
We decompose the problem into two parts: computing an alternating list of length 2 k -1 from an arbitrary list of labels, and then transforming this alternating list into a balanced tree. We first show how to solve the second part: turning an alternating list into a full tree.

Given an alternating list lst, by pairing the trees in lst using only one traversal of the list, we obtain an alternating list with exactly half as many trees. Each pairing requires two trees and one label used as a root. In order to build a list that is alternated, we also need a second label, that is kept as a single element. This explains why we consider at each step the four first elements of the list.

A single traversal, encoded by pass : α tree_or_elt list → α tree_or_elt_list, reduces an alternating list of length 2 k -1 3 to an alternating list of length 2 k-1 -1. By iterating this process using loop : α tree_or_elt list → α tree, we reduce the original list to a list of length one, whose one element is a balanced tree t such that the infix traversal of t is the initial list.

let rec loop = function [] → Leaf [Tree t] → t list → loop (pass list)
Notice how the invariant that alternating lists have length 2 k -1 is maintained: this is because, for k 2, we have

2 k -1 = 4 × (2 k-2 -1) + 3, hence we obtain an alternating list of length 2 × (2 k-2 -1) + 1 = 2 k-1 -1.
It remains to show how to transform a list of labels of length n into an alternating list of trees and labels. Each of the original trees has height zero or one: they are leaves or contain only one label. Because we want a list of length precisely 2 k -1, for k = 1 + ⌊log 2 n⌋, it means we need 2 k -1n leaves. This quantity is computed as the variable missing. The function pad computes the alternating list by creating as many leaves as needed, alternating them with elements, and once enough leaves are created, promotes all the odd-position labels into trees. The balancing algorithm balance: α list → α tree is thus given by the composition of loop with complete:

let complete list = let n = List.length list in let rec pow2 i = if i <= n then pow2 (2*i) else i in let missing = ( pow2
let balance list = loop (complete list)
As for the complexity of this algorithm, notice that pass and complete are both clearly in linear-time in the length of the lists on which they work, while loop recurses on lists whose length are halved at each recursive step. Hence balance is a linear-time algorithm.

Removing partial functions

The loop function of Section 2 relies on the invariant that the list argument has length 2 k -1. Additionally, all the odd-position values must be of the form Tree t, whereas all the even-position values must be of the form Elt x. If either of these invariants is broken, we would run into the assert false of pass.

It is not immediately apparent that these properties hold. If it does not take a tremendous effort to convince oneself that the balance function of Section 2 is indeed correct, a direct mechanically checked proof would not be very practical.

Length invariants

Our goal in this section is to avoid resorting to assert false. In addition to making sure that balance indeed terminates with a value, it will make it considerably simpler to implement the balancing algorithm in Coq in Section 4. To achieve this goal, it is necessary to have more precise types. Let us focus first on the length invariants: we will need to define a type which contains exactly the non-empty lists of length 2 k -1.

A data structure which holds 2 k -1 elements brings complete binary trees to mind. Even if it is possible -though not necessary convenient -to represent complete binary trees in Ocaml, they are not the appropriate structure. First, because complete binary trees are full trees and are, hence, unlikely to serve as a useful intermediate data structure to build a full tree. Second because there is a simpler -albeit more exotic -alternative.

Indeed, lists can be seen as decorated unary numbers: there is an element at each successor. Different kinds of lists can be obtained, more or less systematically, by varying the numerical representation. This idea goes back to Guibas & al. in [START_REF] Guibas | A new representation for linear lists[END_REF] and a fairly thorough exploration can be found in Okasaki [2, Chapters 9&10]. In the simplest cases, the analogous list structure corresponds to a structurally recursive exponentiation algorithm. For regular lists, a list of size n whose elements have type a can be recursively defined with the following equations:

a 0 = 1 a n+1 = a × a n
Replacing unary numbers with binary numbers, we obtain the binary exponentiation algorithm:

   a 2 0 -1 = 1 a 2n = (a 2 ) n a 2n+1 = a × (a 2 ) n Okasaki [2,
Chapter 10] uses a non-uniform inductive type to encode the latter exponentiation algorithm into a type of lists he calls binary lists. We are only interested in lists of length 2 k -1, that is a length written only with the digit 1 in binary representation. So following Okasaki, but skipping the second equation above (which corresponds to the digit 0) we define the following non-uniform inductive type, which we call power lists:

module PowerList = struct type α t = Zero TwicePlusOne of α * (α*α) t end
This type actually appears in Okasaki [2, Chapter 10] as an introduction to non-uniform binary lists. Relatedly, Okasaki [START_REF] Okasaki | From fast exponentiation to square matrices[END_REF] leverages a tail-recursive binary exponentiation algorithm to define a type capturing precisely square matrices; on the other hand, Myers [START_REF] Myers | An applicative random-access stack[END_REF] introduced a flavour of list based on skew binary numbers which are not easily captured as exponentiation.

Although the power lists may look like some sort of trees, it is not a very accurate depiction. The easiest way to picture how power lists works is to see TwicePlusOne as a fancy (::), then the lists with, respectively, 1, 3, 7, and 15 elements are as follows:

-[1] -[1;(2,3)] -[1;(2,3);((4,5),(6,7))] -[1;
(2,3);((4,5),(6,7));(((8,9),(10,11)),((12,13),(14,15)))]

Elements appear in order, like in a regular list, but they are packed twice as tightly after each TwicePlusOne.

Just like with regular lists, there is a map function for power lists. Due to the non-uniformity it is a little trickier3 than the regular list map: in the recursive steps, the argument function f needs to process two consecutive elements at a time.

module PowerList = struct . . .

let rec map : α β. (α→β) → α t → β t = fun f → function Zero → Zero TwicePlusOne (elt,lst) → let f' (x,y) = f x , f y in TwicePlusOne (f elt , map f' lst) end 3.

Alternation

In Section 2, labels are separated from trees dynamically. The pass function verifies that trees and labels are interleaved properly, and fails if they are not. In this section, instead, we consider a variant of α PowerList.t where every odd position contains a tree, and every even position contains an element. More generally, we define a type (ω,η) AlternatingPowerList.t where odd positions have type ω, and even positions have type η. Such a list should have the following pattern:

- [ω;(η,ω);((η,ω),(η,ω))] After the first element, which must have type ω, there is no difference between even and odd positions: indeed, excluding the first element, we are actually building an (η*ω) PowerList.t. Hence the definition:

module AlternatingPowerList = struct type (ω,η) t = Zero TwicePlusOne of ω * (η*ω) PowerList.t end
For brevity, let us write PL and APL for PowerList and AlternatingPowerList respectively.

Using these alternating power lists, we can define a version of the pass function free of assert false. Indeed, consider an alternating power list of length at least 3: it is of the form APL.TwicePlusOne (a, PL.TwicePlusOne ((b,c), lst)), where lst has type ((η*ω)*(η*ω)) PowerList.t. The pass function of Section 2, as it happens, manipulates its arguments by groups of four elements: basically, pass is simply a map over lst.

We hence define the function pass which joins the trees in a list of length 2 k+2 -14 , producing a list of length 2 k+1 -1. The function loop is virtually unchanged from Section 2, except it acts on power lists:

Padding

Now that there is no more assert false in the code of loop, we need to change the complete function of Section 2 so that it returns an (α tree,α) APL.t rather than a list. The heart of this section is a function which turns an α list into an (α*α tree) PL.t. The final function, which produces an (α tree,α) APL.t is a simple wrapper around the former. We want to turn a list lst of length n + 1 into a pair of its first element, converted into a tree, plus a power list of length 2 × (2 k -1) n representing its tail tail. Each element of the power list is a pair, whose first term is an element, and its second term is a tree of height at most one. In particular, the length of the returned power list is always even, so if tail has odd length, we will need to insert at least a Leaf. This suggests that we may inspect the parity of the length of tail, and insert an extra element precisely when it is odd. This leads to a slightly different padding procedure than that of Section 2, in particular the leaves are not inserted at the same position, but it is inconsequential.

An α list of even length can be turned into an (α*α) list whose length is halved. This turns out to be interesting for our recursion, since it mimics the inductive step of power lists. Also, in the case of even length, we need to distinguish the empty case from the non-empty case: the former will be turned into the empty power list APL.Zero while the latter will be turned into a power list of the form APL.TwicePlusOne((x,y),l), where x and y correspond to the two first elements of tail. These different cases are represented in the following view: The padding function itself, of_list, is at first sight far from intuitive. Let us recall that we want to turn a list of labels of arbitrary length, into a power list of pairs. A label can be thought of as a bit of weight 2 0 , while a pair of labels would be a bit of weight 2 1 , and so on. At first, all our bits have weight 2 0 and consists in one label each. We can build bits of higher weight by pairing up two bits of the same weight. A bit made up only of labels is called pure. We can also double the weight of a bit by interlacing leaves with it (with the function pad), but this gives a bit made of pairs of labels and trees, call them impure. Lastly, we can also transmute a pure bit into an impure bit of the same weight (with the function coerce), by replacing odd-position labels by trees of height one.

Each recursive step consists in taking a list of pure bits of the same weight 2 k , and outputing exactly one impure bit of size 2 k+1 , plus a list of pure bits of weight 2 k+1 , which is converted recursively. We thus obtain, successively, one bit of each weight from 2 1 to 2 l , for some l, encoding a list of length 2 l+1 -2, as expected.

At any recursive step, suppose first that the number of bits of weight 2 k is odd. As we need to compute only bits of weight 2 k+1 , one of them impure, we are forced to use pad on one bit, and to pair up the others. Suppose now that the number of bits of weight 2 k is even. In that case, we can pair them all into bits of weight 2 k+1 , and then use coerce on one of them to make the impure bit.

The last difficulty is that pad and coerce both depend on the current weight of the bits, hence we need to update them at each recursive step. pad must add leaves between every two consecutive labels, in even positions, while coerce must upgrade every even-position label into a tree of height one. This leads to the following definition: The final function, balance:α list → α tree, implements the same algorithm as Section 2 without any partial functions. What we may have lost in this section, compared to the simple algorithm, is the simplicity of the complexity analysis of the algorithm. The subtleties of the main functions require a finer analysis. Consider first the function PL.map: clearly the number of recursive calls depends only logarithmically on the number of elements in the power list. But each recursive call uses as its first argument a function twice as complex than the previous one. This leads to the following inequation over the complexity C(n, m, f ) of map, where n is the number of elements in the power list, m is the size of the elements in the power list, and f is the complexity of the mapped function:

C(n, m, f ) ≤ f (m) + C n -1 2 , 2m, k → 2 × f (k/2) + O(1) + O(1)
From there, it is easy to prove that C(n, m, f ) = n.f (m) + O(n), so that PL.map runs indeed in linear-time, and so is loop. Similarly the complexity of PL.of_list can be described by a higher-order recursive inequation (almost the same as above, except that the complexity depends on two functions and the constant term is replaced by a linear term), whose solution gives also a linear-time complexity.

Turning to Coq

There is still a property of the algorithm that the implementation of Section 3 does not make obvious: that the algorithm actually does build full trees. In this section we shall build into the type of balance that its output is indeed full.

To that effect, we will use Coq rather than Ocaml. Even if it is possible, with some effort, to represent full trees and implement the algorithm in Ocaml -and relatively easy in Haskell -a Coq implementation also gives us termination by construction. Coq forces every recursion to be structural, which will prove to be rather entertaining.

At a superficial level, a visible difference with the Ocaml implementation is that Powerlist.t and AlternatingPowerList.t must be decorated with the k such that the length is 2 k -1: it is the structural recursion parameter of the balance_powerlist function. Because it makes the code simpler, we will use a recursive definition rather than an inductive one:

Module PowerList. Fixpoint T (A:Type) (k:nat) := match k with 0 ⇒ unit:Type S k' ⇒ A * T (A * A) k' end.
End PowerList.

We will also need a version where k can be arbitrary. For that purpose we use

End PowerList

The definition of AlternatingPowerList.T and AlternatingPowerList.U are similar.

Full trees

To code full trees, we index trees by their height, and specify that leaves can happen only at height 0 or 1:

Inductive FullTree (A:Type) : nat → Type := Leaf 0 : FullTree A 0 Leaf 1 : FullTree A 1 Node {k:nat} :

FullTree A k → A → FullTree A k → FullTree A (S k).
If we omitted the constructor Leaf 1 , we would have a definition of complete binary trees: both subtrees of a node are complete binary trees of the same height.

We allow the full trees to be incomplete by letting Leaf 1 take the place of nodes on the last level. Using the type FullTree A k in place of the type α tree, the functions pass and balance_powerlist are virtually unmodified5 with respect to Section 3. Only their types change to reflect the extra information:

Definition pass {A k p} : APL.T (FullTree A (S p)) A (S (S k)) → APL.T (FullTree A (S (S p))) A (S k). Fixpoint loop {A k p} : APL.T (FullTree A (S p)) A (S k) → FullTree A (plus k (S p)) {struct k}.

The algorithm indeed builds only full trees.

Structural initialisation

The padding conversion from lists to power lists, in Section 3, is not structural due to the use of pair_up in the recursive call. To tackle this recursion, we shall make use of another intermediate structure. What we need, essentially, is that all the calls to pair_up are pre-calculated, so the intermediate structure will be like parity except that the calls to (α*α) list are replaced by inductive calls.

As it turns out, this is another non-uniform datatype which corresponds to a numerical representation. Indeed, any natural number can be written in binary with digits 1 and 2 (but not 0). In this system, for example, 8 = 1×2 2 +1×2 1 +2×2 0 is represented as 112. Here is the definition, where tpo reads "twice plus one" and tpt "twice plus two": End BinaryList.

Module BinaryList.
To turn a non-empty list into a BinaryList.T, all we need is a function cons of type A → T A → T A to add an element in front of the list. On the numerical representation side, it corresponds to adding 1. It behaves like adding 1 in the usual binary representation, except that 1-s are turned into 2-s without a carry and 2-s into 1-s while producing a carry: Module BinaryList.

. . . End BinaryList.

Note that while cons takes, in the worst case, logarithmic time with respect to the length of the list, building a list by repeatedly using cons is still linear. Indeed, as previously mentioned, cons mimics the successor algorithm for binary numbers, whose amortized complexity is well-known to be constant. We also need a function which turns a T (A * A) into a T A. This is effectively multiplication by 2. The lack of 0 among the digits 6 makes this process recursive. A simple presentation of the doubling algorithm consists in adding a 0 at the end of the number, then eliminating the 0 using the following equalities:

   0 = • x20 = x12 x10 = x02
In terms of binary lists: Module BinaryList.

. . . End BinaryList.

We can now write a structurally recursive padding function, using binary lists as the structural argument. As we do not know in advance the length of the produced list, a PowerList.U is returned. We write BL as a shorthand for BinaryList:

Module PowerList. . . . Fixpoint of_binary_list {A X} (d:A→X) (f:A * A→X) (l:BL.T A) : U X := match l with BL.zero ⇒ zero BL.tpo a l ⇒ tpo (d a) (of_binary_list (d×d) (f×f) l) BL.tpt a b l ⇒ tpo (f (a,b)) (of_binary_list (d×d) (f×f) l) end. End PowerList.
Where g×f is the function which maps (x,y) to (g x,f y).

Conclusion

The balance function of Section 4 is, by virtue of its type alone, a total function which turns lists into full binary trees. Yet, to the cost of using intermediary data-structures, it effectively implements the algorithm of Section 2.

The missing piece is to prove that the infix traversal of balance l is indeed l. The infix traversal of a (full) tree is represented in Coq with the functions We can then state the theorem:

Theorem balance_preserves_order A (l:list A) : list_of_full_tree (balance l) = l.

The proof is short and straightforward: we define a traversal function for each intermediate structure; and state a variant of balance_preserves_order for each intermediate function. Proving the intermediate lemmas is not difficult and can be mostly automatised: we use a very simple generic automated tactic, which discharges most goals. This theorem concludes our easy formal proof of the balancing algorithm.

  let join left node right = Tree (Node (left, node, right)) let rec pass = function Tree left :: Elt root :: Tree right :: Elt e :: others → join left root right :: Elt e :: pass others [Tree left; Elt root; Tree right] → [join left root right] _ → assert false

  1) -n -1 in let rec pad missing = function head::tail when missing <> 0 → Tree Leaf :: Elt head :: pad (missing -1) tail odd::even::others → join Leaf odd Leaf :: Elt even :: pad 0 others [single] → [join Leaf single Leaf] [] → [] in pad missing list

  let pass left (root,right) apl = let join_up ((single,left),(root,right)) = single, Node (left,root,right) in APL.TwicePlusOne ( Node (left,root,right) , PL.map join_up apl) let rec loop : ǫ. (ǫ tree,ǫ) APL.t → ǫ tree = function APL.Zero → Leaf APL.TwicePlusOne (tree,PL.Zero) → tree APL.TwicePlusOne (tree,PL.TwicePlusOne (pair,apl)) → loop (pass tree pair apl)

  type α parity = Empty Odd of α * (α*α) list Even of (α*α) * (α*α) list let pair_up lst = let succ elt = function Empty → Odd (elt, []) Odd (b,pairs) → Even ((elt,b), pairs) Even (bc,pairs) → Odd (elt, bc::pairs) in List.fold_right succ lst Empty

module

  PowerList = struct . . . let rec of_list : α β. (α→β) → (α*α→β) → α list → β t = fun pad coerce bits → let pad' (x,y) = (pad x, pad y) in let coerce' (x,y) = (coerce x, coerce y) in match pair_up bits with Empty → Zero Odd (a,pures) → TwicePlusOne (pad a, of_list pad' coerce' pures) Even (ab,pures) → TwicePlusOne (coerce ab, of_list pad' coerce' pures) end With that function, we can conclude our implementation. Again writing PL and APL for PowerList and AlternatingPowerList respectively: module AlternatingPowerList = struct . . . let of_list leaf up id = function [] → Zero a::l → let pad x = id x , leaf in let coerce (x,y) = id x , up y in TwicePlusOne (up a, PowerList.of_list pad coerce l) end let singleton x = Node(Leaf,x,Leaf) let balance l = loop (APL.of_list Leaf singleton (fun e→e) l)

  Coq's type of dependent pairs { n:nat & F n}. The constructor for dependent pairs is written n , x . The implicit version comes with constructorstpo stands for "twice plus one": Module PowerList. . . . Definition U (A:Type) := { k:nat & T A k }. Definition zero {A:Type} : U A := 0 , tt . Definition tpo {A:Type} (a:A) (l:U (A * A)) : U A := let ' k,l := l in S k , (a,l) .

  Inductive T (A:Type) : Type := zero tpo (a:A) (l:T (A * A)) tpt (a b: A) (l:T (A * A)).

Fixpoint

  cons {A} (a:A) (l:T A) : T A := match l with zero ⇒ tpo a zero tpo b l ⇒ tpt a b l tpt b c l ⇒ tpo a (cons (b,c) l) end. Definition of_list {A} (l:list A) : T A := List.fold_right cons zero l.

Fixpoint

  twice {A} (l:T (A * A)) : T A := match l with zero ⇒ zero tpo (a,b) l ⇒ tpt a b (twice l) tpt (a,b) cd l ⇒ tpt a b (tpo cd l) end.

  The rest follows straightforwardly, and we can define the following functions which conclude the algorithm (BL, PL, and APL stand for BinaryList, PowerList, and AlternatingPowerList respectively): Module AlternatingPowerList. . . . Definition of_binary_list {A Odd Even} (d:Odd) (f:A→Odd) (g:A→Even) (l:BL.T A) : U Odd Even := match l with BL.zero ⇒ zero BL.tpo a l ⇒ let d' x := ( g x , d ) in tpo (f a) (PL.of_binary_list d' (g×f) (BL.twice l)) BL.tpt a b l ⇒ let d' x := ( g x , d ) in tpo (f a) (PL.of_binary_list d' (g×f) (BL.tpo b l)) end. Definition of_list {A Odd Even} (d:Odd) (f:A→Odd) (g:A→Even) (l:list A) : U Odd Even := of_binary_list d f g (BL.of_list l). End AlternatingPowerList. Definition singleton {A:Type} (x:A) : FullTree A 1 := Node Leaf 0 x Leaf 0 . Definition balance {A:Type} (l:list A) : { k:nat & FullTree A k } := let ' k,l := APL.of_list Leaf 1 singleton (fun x⇒x) l in match k with 0 ⇒ fun _ ⇒ 0 , Leaf 0 S k ⇒ fun l ⇒ plus k 1 , loop l end l.

Fixpoint

  list_of_full_tree_n {A n} (t:FullTree A n) : list A := match t with Leaf 0 ⇒ [] Leaf 1 ⇒ [] Node _ t 1 x t 2 ⇒ list_of_full_tree_n t 1 + + [x] + + list_of_full_tree_n t 2 end.Definition list_of_full_tree {A} (t:{ k:nat & FullTree A k }) : list A := list_of_full_tree_n (projT2 t).

The type annotation on PowerList.map informs Ocaml that map is a non-uniform recursive function. Without the type annotation, Ocaml simply assumes that map is uniformly recursive and fails to typecheck since f and f' have different types.

To ensure that its argument list has at least three elements, pass takes the three first elements as extra arguments. In other words pass t (x,s) l is meant to be read as pass (APL.TwicePlusOne (t , PL.TwicePlusOne ((x,s),l))).

In fact, as can be seen from its type, loop only handles non-empty alternating power lists. This is due to a small technicality: the recursive step of loop is the case S (S k), but Coq does not recognise S k as a structural subterm of S (S k), so the definition from Section 3 does not verifies Coq's structural recursion criterion. As a workaround, the empty case is moved to the balance function.

The constructor zero represents an empty list of digits.

⋆ This research has received funding from the European Research Council under the FP7 grant agreement 278673, Project MemCAD