Packing Curved Objects

1 Introduction

The packing problem consists in placing n items inside a given space, such that no two items overlap. A large number of real-world applications like warehousing, logistics or parallel computing gives a great relevance to this problem. Packing has been well studied in either academic or industrial contexts, but each time considering specific shapes like bins [Lodi et al., 2002], circles [Stephenson, 2005] or polytopes

Gilles Chabert

Mines de Nantes - LINA (UMR 6241), France gilles.chabert@mines-nantes.fr

[Stoyan et al., 2005; Egeblad et al., 2009]. A lot of problems are actually studied for polytopes, some being closely related to packing, like the collision detection [Moore and Wilhelms, 1988; Brochu et al., 2012; Mainzer and Zachmann, 2015] or the calculation of the penetration depth [Cameron and Culley, 1986; Dobkin et al., 1993; Kim et al., 2002; 2004].

Our goal is to deal with the more general situation where different shapes can be mixed, including non-convex and curved shapes.

In this context, one has to develop a generic approach for solving the packing problem. One such approach has been recently proposed in [Martinez et al., 2013]. The approach splits the problem in three parts. First, they build a violation function fij for each pair (i, j) of objects. This function takes as argument the position ⃗o = (x, y) and the orientation α for both objects and gives a measure of "how much" the two objects overlap. We shall call such function an overlapping function.

Second, all these violations are summed up to form a global violation cost f:

Third, the function f is minimized using a generic blackbox optimization software called CMA-ES [Hansen and Ostermeier, 2001].

The experimental results they obtain are very encouraging. However, there is a serious bottleneck in their approach: the overlapping functions are ad-hoc formulas for every possible combination of object shapes. For instance, in a problem where objects are either circles, squares or triangles, one must build a formula for the 6 combinations: circle-circle, circle-square, circle-triangle, square-square, square-triangle and triangle-triangle. Some formulas are easy to obtain, like for the circle-circle case. Indeed, given two circles of radii ri and rj, the overlapping can be defined as the distance between the two centers and 0 if this distance exceeds the sum of the radii, that is:

This simple formula is rather an exception. Mixing shapes considerably complicate the task - even the circle-rectangle case is not so simple-. One should be convinced that, in general, these ad-hoc formulas are difficult to write. And their number grows in O(n2 ).

We show in this paper how to replace these formulas by an algorithm. This algorithm has also been implemented and this paper reports the preliminary experimental results we have obtained on original packing problems.

The paper is organized as follows. We first formalize the problem in §2, that is, the way shapes are represented and which overlapping function we actually consider. Then, we detail in §3 and §4 an algorithm for this overlapping function. Experimental results are reported in §5, followed by concluding remarks.

2 Formalism

2.1 Object definition

We consider objects described by nonlinear inequalities. This representation choice encompasses basic geometric shapes like half-planes or ellipses. We also accept the logical and/or operators in the shape description, which means that these geometric primitives can be combined to define more complex shapes like a polygon, a half-circle or a horseshoe (see Figure 1). Note that non-linear inequalities give also freedom to represent a lot of fancy shapes but real-world objects are rather obtained by combining basic primitives.

With that representation, the shape of Object i can be seen as a regular constraint ci(⃗ p). This constraint means that the point ⃗ p belongs to the object when the latter is placed at the origin of the frame and with a null rotation angle (see Figure 1).

We will identify in the sequel a shape to a constraint and talk about "an object of shape c".

2.2 The non-overlapping constraint

An object can be translated or rotated. If Object i is placed at a position ⃗ oi and rotated by an angle αi, we will denote by ⃗qi and call the parameters of Object i the vector:

The size of ⃗q is the number of degrees of freedom. If only translation is considered, the angle is dropped and ⃗q has only two components. Otherwise, it belongs to R3 .

Given a vector of parameters ⃗qi, the points ⃗ p that belong to the translated-rotated Object i are the ones satisfying

where R is the usual rotation matrix. Therefore, Objects i and j overlap iff

The non-overlapping constraint between two objects is just the negation of the latter relation.

2.3 The packing problem

The packing problem is a set of pairwise non-overlapping constraints between n objects and an additional constraint that all objects must be packed inside some container (like the enclosing bin in the bin packing).

One appealing aspect of our representation is that we don't actually need a specific treatment for this additional constraint. This is based on a very simple property. An object is inside the container iff this object does not overlap the complementary of the container. And the complementary of a shape is simply obtained with the negation of the underlying constraint. So, in our system, the packing problem consisting in placing n objects of shapes c1, . . . , cn inside a container c is reformulated as a non-overlapping constraint between n+1 shapes c1, . . . , cn, ¬c.

Note that, for obvious symmetry breaking reasons, one has to fix the origin and the orientation of one object. A natural choice is then to fix that of the container.

2.4 The overlapping function

The overlapping function fij takes as input two vectors of parameters ⃗qi and ⃗qj and gives a measure of "how much" Objects i and j overlap. The output must fulfill two properties. First, it has to be 0 iff the two objects are disjoint. Second, it has to decrease when the two objects gets more distant.

We define the overlapping function as the distance between the current vector of parameters for Object j, and the closest vector of parameters that satisfies the non-overlapping constraint. This way, the function fij can be seen as a "distance to satisfaction" for the Object j, the parameters of Object i being considered as fixed. In fact, this function exactly matches the concept of penetration depth in the realm of computational geometry [Cameron and Culley, 1986; Dobkin et al., 1993; Kim et al., 2002; 2004].

Definition 1 (Overlapping Function).

Note that the distance ‖.‖ can be either in 2 or 3 dimensions, depending if rotation is considered or not.

As already said, the previous definition introduces a distinction between the two shapes. Object i is called the reference object and Object j the moving object. Without rotation, the roles are symmetric in the sense that:

However, the previous equality does not hold with rotation. Henceforth, we will always assume Object i (resp. j) to be the reference (resp. moving) one.

It is important to notice that our definition of the overlapping function does not coincide with the overlapping surface, as Figure 2 shows. Furthermore, taking the overlapping surface as a measure of overlapping would not be appropriate since the surface can reach a local minimum in a situation where the two objects still overlap (see Figure 2). In contrast, our overlapping function only reaches a minimum at 0, that is, when the objects are disjoint. This is an important difference as this function is eventually called by a local search that can be trapped in local minima.

2.5 The overlapping region

The implementation we will propose for the overlapping function is based on the concept of overlapping region, a concept that was already introduced in [Salas et al., 2014]. We recall this concept and then give the theoretical result that makes the connection with our definition of overlapping function.

The overlapping region simply corresponds to the set of parameters for Object j that make the objects i and j overlap, the parameters of Object i being all set to 0:

Definition 2 (Overlapping Region). Given Objects i and j (resp. the reference and moving ones), the overlapping region is

This concept is depicted in Figure 3.

We show now the relation that connects the overlapping function to the overlapping region.

Proposition 1. Given two vectors of parameters ⃗ qi and ⃗ qj , let

Then

Proof. By definition,

Let us introduce the variable

We have, equivalently,

In the one hand:

In the other hand it can be proven (see Proposition 1 in [Salas et al., 2014]) that

that is,

Hence, the overlapping function maps the current position of Objects i and j to the minimal distance between ⃗r(⃗ qi, ⃗qj ) and the boundary of the overlapping region. This is also depicted in Figure 3.

3 Main Algorithm

The algorithm we propose for the overlapping function is a direct consequence of Proposition 1. We first calculate (offline) the overlapping region and then obtain (on-line) the value of f(⃗qi, ⃗qj) from it. More precisely, the two steps are:

1. We entirely calculate the overlapping region Sij . Calculating this region means that an explicit representation of Sij is computed, that is, under the form of numerical data. The data structure we use for Sij is detailed below.
2. We calculate the minimal distance between ⃗r(⃗ qi, ⃗ qj) and Sij using the previous data. The point is that this minimization problem is far more easy to solve than the original minimization problem involved in Definition 1.

3.1 Paving of the overlapping region

Let us start with the first step. The data structure we use for Sij is called a paving (see Figure 4).

Definition 3 (Paving). A paving of a set S ⊂ Rd is a triplet (I, B, O) where I (for "inside"), O (for "outside") and B (for "boundary") are three sets of boxes verifying

To calculate this paving, we apply the strategy proposed in [Salas et al., 2014]. It ensures that the total surface of the boundary (or, in other words, the precision of the paving) is less than a given value ε. Let us briefly outline this algorithm again, since we will describe an adaptation in the next section.

The algorithm is called a paver [Chabert and Jaulin, 2009], a recursive algorithm that starts with an arbitrarily large box [qj] and alternates three steps:

1. (outer rejection test). If unsatisfiability of (2) is proven for all ⃗ qj ∈ [qj] (⃗ qi being fixed), then O := O ∪ {[qj ]}. The rejection test resorts to an embedded solver that we shall not describe further, since this part has not changed from the paper cited above.
2. (inner inflation). Pick randomly a point ˜ qj ∈ [qj] and check if it belongs to Sij. If it does, then "inflate" ˜ qj to a subbox [˜ qj] proven to entirely lie inside Sij. In contrast, this part has been significantly improved and will be the topic of Section 4. Then, I := I ∪ {[˜ qj ]}. Break the remaining part [qj]\[˜ qj ] into boxes and perform a recursive call with each box.
3. if [qj] is small enough, B := B∪{[qj]}. Otherwise bisect [qj ] into two sub-boxes and perform a recursive call with each boxes.

3.2 Distance to the boundary set

Once the paving (I, B,O) is calculated, we first set ⃗r := ⃗r(⃗ qi, ⃗ qj) using (3) and then find the closest point to ⃗r which is not inside the overlapping region. However, because the boundary of the region is uncertain, we can actually only obtain an enclosure of this minimal distance:

Clearly, the accuracy of the enclosure is related to the precision ε the paving has been generated with.

In practice, to calculate a bound, say the upper one, we consider each box [q′′ ] of the set O generated by the paver, and minimize the distance over that box. That is, we calculate:

Then, we have:

The same holds for the lower bound.

Formula (5) is nothing different than the distance between a box and a point. This distance can be easily obtained in constant time (whatever is the size of the box), either using interval arithmetics or by a direct geometric argument.

Formula (6) can however be very time-consuming if the boxes in O are stored in a flat unsorted list. One has to systematically scan the whole list to get the minimal distance. This linear complexity can be easily transformed to a logarithmic complexity using a tree-structure representation, as depicted in Figure 4. The structure is then explored using standard global optimization techniques. In particular, pending nodes are stored in a heap, the criterion associated to a node being precisely its distance to ⃗r. Nodes with a distance exceeding the current upper bound for the minimum can then be discarded (just like the bounding step in a branch & bound).

4 A new inflator for the overlapping constraint

The inflation plays a key role in the paving algorithm and thus has a strong impact on the efficiency of our packing approach.We describe in this section an original inflation technique that substantially improves packing performances.

Recall that the inflator builds, from a single vector of parameters ˜ qj, a box [qj] that is inside the overlapping region Sij .

In [Salas et al., 2014], the inflator proposed inflates in all the dimensions of ⃗q simultaneously and turns out to be not efficient when the dimension is 3 (i.e., with rotation handled). In fact, this inflator is very sensitive to the number of times a variable appears in Equation (1), and the angle α has precisely many occurrences yield by the rotation matrix R-α. It was reported that the time required to calculate 3D pavings with this inflator was prohibitive.

Our new inflator does not suffer from the multiple symbol occurrences introduced by rotation matrices. It splits the inflation in two steps: it first inflates with respect to ⃗o and then with respect to α. We describe both steps and then show how the two inflations can be combined.

4.1 Translation

The inflation w.r.t. ⃗o is based on the following principle (see Figure 5). 1 Given initial parameters ( ˜ oj, ˜ αj) we first look for a vector ˜ p that belongs to both objects, the parameters of Objects i and j being respectively ⃗0 and ˜ qj, that is, ˜ p satisfies

Then, we inflate ˜ p inside Object i using generic inflation technique for inequalities [Chabert and Beldiceanu, 2010; Araya et al., 2014]. Let u := ˜ p - ˜ oj. The resulting box [⃗ p] can be translated to ˜ oj in the sense that if oj moves inside [ ⃗ oj] := [⃗ p ] - ⃗u then ⃗ oj + ⃗u is both inside [⃗ p ] and inside the Object j with parameters (⃗ oj, ˜ αj ).

1 This part is similar to the old inflator applied in the 2D case.

4.2 Rotation

The inflation for the angle works as follows (see Figure 6). Given initial parameters ( ˜ oj, ˜ αj) we look for the two angles α and α that force vector ˜ p (the same as before) to meet the boundary of Object j (this step is detailed in the next subsection). Then, it is clear that for every angle αj ∈ [α, α] there is an overlapping. However, some caution is required. First, there are generally more than 2 candidate angles (see Figure 6.d), especially if Object j is non-convex. The angle values that form the narrowest interval around ˜ αj are kept. Second, there may be no angle at all. This means that Object j can make a complete turn while containing ˜ p (or, in other words, ˜ p is in the inscribed circle of the object). Third, we have to take into account the 2π-periodicity. Indeed, the search space for angles is a bounded domain like [-π, π] or [0, 2π]. This has the unpleasant consequence to make the inflated domain disconnected, e.g., if ˜ αj = π2 , α = 0 and α = 3π 2 then the inflated interval [0, 3π 2 ] becomes [-π,- π2 ] ∪ [0, π]. So, in this case, we only keep the subinterval containing ˜ αj, that is, [-π, - π2 ].

4.3 Shape Boundary

Let us explain now how the angles α and α are calculated. First, we have to generate a constraint c′j for the boundary of Object j. In the simplest case of a single inequality,

the boundary of the object is described by the equality f(⃗ p) = 0.

In the case of a conjunction of inequalities,

The boundary is described by a disjunction of n systems:

In the most general case of a disjunction of conjunctions of inequalities,

The boundary is described by a disjunction of (n1+. . .+nm) conjunctions, each block of conjunctions being obtained with the previous pattern (see Figure 7).

The boundary constraint c′j can therefore be automatically generated by a simple symbolic processing.

The angle α and α are then two solutions of the system c′j (R-α( ˜ p - ˜ oj)) which is a disjunction of 1-dimensional equations (of variable α), i.e., something very easy to solve.

Note that if the shape is a disjunction with two subparts that plainly intersect (contrary to the figure), the boundary is only a subset of the generated system, in which case the angle inflation may not be optimal.

4.4 Global Inflation

We finally show that the two inflations can be combined, that is, the cartesian product of [⃗ oj] and [αj] = [α, α] is a valid inflation of ˜ qj .

The inflation w.r.t. translation produces a box [⃗ oj] that satisfies

Together, they give

5 Experimental Results

The goal of these preliminary experiments is to validate the theory and to have a first impression of how the algorithm behaves in practice. The ambition is neither to solve real-world problems nor to make an extensive performance analysis.

5.1 Setup

We have created 5 cases of increasing difficulty. In each case, n objects have to be placed with n = 10, 20 and 30. Each time, the size of the container has been adjusted manually so that the density of the resulting packing looks similar.

Case 1. This is a standard circle packing problem. The objectives are, first, to check that our algorithm also succeeds in solving this problem and, second, to measure its overhead as compared to a dedicated approach. In this circle packing instance, the radius of the ith circle is equal to √i, 1 ≤ i ≤ n. The radius of the enclosing circle is set to 2.1 (n=10), 2.35 (n=20) and 2.48 (n=30).

Case 2. The purpose of this case is to introduce a new shape, that is, for which (to our knowledge) no dedicated solver exists. The problem is otherwise kept as simple as possible: (1) there is only one shape in addition to the enclosing shape, (2) this shape is described by a single inequality and (3) only translation is permitted (no rotation). The shape we have chosen is an ellipse of radii 1 and 0.5. The ellipses have to be packed inside a circle of radius 2.9 (n=10), 4 (n=20) and 5 (n=30).

Case 3. We increase the complexity by allowing rotation in the same ellipse packing instance. With this new degree of freedom, the enclosing circle can be chosen a little bit smaller. The radius of the circle is set to 2.7 (n=10), 3.9 (n=20) and 4.9 (n=30).

Case 4. We increase again the complexity by replacing the ellipse with a non-convex shape which description involves several inequalities. The shape is the horseshoe depicted in Figure 1. The enclosing shape is an ellipse of radii 2/4 (n= 10), 3.9/1.8 (n=20) and 5/2.9 (n=30).

Case 5. The last case is a mixed packing problem where n/2 ellipses and n/2 horseshoes have to be packed inside a circle of radius 2.9 (n = 10), 4 (n = 20) and 6 (n = 30). Again, both translation and rotation are considered.

5.2 Results

Remind that the packing process is in two steps. The paving step, performed off-line, and the packing step. Since some pavings are used in several cases, paving and packing times are reported separately. We also only give the average time for the 13×12 circle-circle pavings calculated in Case 1 and for the other pavings involving the enclosing shape (since this shape changes with n). Tables 1 and 2 summarize the paving and packing times. Figure 8 shows pictures of the results. Every paving has been calculated with a precision ε set to 0.1.

First, we can see that the circle packing problem has been solved in a few seconds only. However, the enclosing circle is not optimal. In [Martinez et al., 2013], the circles are packed in a circle of radius 1.95 (n=10), 2.15 (n=20) and 2.27 (n=30). This gap is due to the precision ε. Using this precision for paving amounts to consider that objects to pack are "enlarged by" ε. The radius we have set for the enclosing circle is actually the smallest one for which packing is successful, given the precision ε of 0.1.

Packing for the other cases has been done in the order of the minute. Above all, we see in Figure 8 that the holes introduced by non-convex objects like the horseshoe are not lost. The paving time is in the order of the hour.

6 Conclusion

We have proposed a generic algorithm for packing non convex curved objects. The approach is not limited to polyhedral shapes and does not make convex approximation. Cavities introduced by non-convexity are used to place objects.

The approach is split in two steps. The paving step calculates off-line the set of relative positions and orientations of one object with respect to a second one that make them overlap. The packing step minimizes a fitness function which is a sum of distances between points and boundaries of the previously calculated sets. These distances are fast to calculate thanks to a tree-structured representation of the sets.

The other contribution is a new inflator for the overlapping constraint, i.e., the paving step. This new inflator uses separate techniques for inflating with respect to the position and the angle, which leads to a much larger inflation (and therefore a better packing time) than with previously existing techniques.

Preliminary experimental results are encouraging. One possible future experiment could be to see how both paving Case 5

Figure 8: Packing results.

and packing time increases as ε decreases and to mix this numerical approach with the available ad-hoc formulas.