Pedro Machado

Manh Ães

D E Castro

AND Quentin M Érigot

Boris Thibert

INTERSECTION OF PARABOLOIDS AND APPLICATION TO MINKOWSKI-TYPE PROBLEMS

In this article, we study the intersection (or union) of the convex hull of N confocal paraboloids (or ellipsoids) of revolution. This study is motivated by a Minkowskitype problem arising in geometric optics. We show that in each of the four cases, the combinatorics is given by the intersection of a power diagram with the unit sphere. We prove the complexity is O(N ) for the intersection of paraboloids and Ω(N 2 ) for the intersection and the union of ellipsoids. We provide an algorithm to compute these intersections using the exact geometric computation paradigm. This algorithm is optimal in the case of the intersection of ellipsoids and is used to solve numerically the far-field reflector problem.

Introduction

The computation of intersection of half-spaces is a well-studied problem in computational geometry, which by duality is equivalent to the computation of a convex hull. Similarly, the computation of intersections or unions of spheres is also well studied and can be done by using power diagrams [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF]. In this article, we study the computation and the complexity of the intersection of the convex hull of confocal paraboloids of revolution, showing that it is equivalent to intersecting a certain power diagram with the unit sphere. Union of convex hull of confocal paraboloids of revolution, and intersection or union of convex hull of confocal ellipsoids of revolution can be studied using the same tools. These studies are motivated by inverse problems similar to Minkowski problem that arise in geometric optics. We show how the algorithm we developed to compute the intersection of paraboloids are used to solve large instances of one of these problems.

Minkowski-type problems. A theorem of Minkowski asserts that given a family of directions (y i ) 1 i N and a family of non-negative numbers (α i ) 1 i N , one can construct a convex polytope with exactly N facets, such that the ith facet has exterior normal y i and area α i . Aurenhammer, Hoffman and Aronov [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF] studied a variant of this problem involving power diagrams and showed its equivalence with the so-called constrained least-square matching problem. This article is motivated by yet another problem of Minkowski-type that arises in geometric optics, which is called the far-field reflector problem in the literature [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF][START_REF] Caffarelli | On the numerical solution of the problem of reflector design with given far-field scattering data[END_REF]. Recall that a paraboloid of revolution is defined by three parameters: its focal point, its focal distance λ and its direction y. We assume that all paraboloids are focused at the origin, and we denote P (y, λ) the convex hull of a paraboloid of revolution with direction y and focal distance λ. We will say in the following that P (y, λ) is a solid paraboloid. Paraboloids of revolution have the well-known optical property that any ray of light emanating from the origin is reflected by the surface ∂P (y, λ) in the direction y. Assume first that one wants to send the light emited from the origin in N prescribed directions y 1 , . . . , y N . From the property of a paraboloid of revolution, this can be done by considering a surface made of pieces of paraboloids of revolution whose directions are among the (y i ). In the far-field reflector problem, one would also like to prescribe the amount of light α i that is reflected in the direction y i . A theorem of Oliker-Caffarelli [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF] ensures the existence of a solution to this problem: there exist unique (up to a common multiplicative constant) focal distances λ 1 , . . . , λ N such that the surface ∂(∩ 1 i N P (y i , λ i )) reflects exactly the amount α i in each direction y i . Other types of inverse problems in geometric optics can be formulated as Minkowski-type problems involving the union of confocal solid paraboloids, and the union or intersection of confocal ellipsoids [START_REF] Oliker | Mathematical aspects of design of beam shaping surfaces in geometrical optics[END_REF][START_REF] Kochengin | Determination of reflector surfaces from near-field scattering data[END_REF].

Contributions. Motivated by these Minkowski-type problems, our goal is to compute the union and intersection of solid confocal paraboloids and ellipsoids of revolution. Using a radial parameterization, each of these computations is equivalent to the computation of a decomposition of the unit sphere into cells, that are not necessarily connected. Our contributions are the following:

• We show that each of the four types of cells can be computed by intersecting a certain power diagram with the unit sphere (Propositions 1, 5 and 7). The approach is similar to the computation of union and intersection of balls using power diagrams in [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF], or to the computation of multiplicatively weighted power diagrams in R d-1 using power diagrams in R d [START_REF] Boissonnat | On the combinatorial complexity of euclidean Voronoi cells and convex hulls of d-dimensional spheres[END_REF]. • We show that the complexity bounds of these four diagram types are different. In the case of intersection of solid confocal paraboloids in R 3 , the complexity of the intersection diagram is O(N ) (Theorem 2). This is in contrast with the Ω(N 2 ) complexity of the intersection of a power diagram with a paraboloid in R 3 [START_REF] Boissonnat | On the combinatorial complexity of euclidean Voronoi cells and convex hulls of d-dimensional spheres[END_REF]. In the case of the union and intersection of solid confocal ellipsoids, we recover this Ω(N 2 ) complexity (Theorem 8). Finally, the case of the union of paraboloids is very different from the case of the intersection of paraboloids. Indeed, in the latter case, the corresponding cells on the sphere are connected, while in the former case the number of connected component of a single cell can be Ω(N ) (Proposition 6). The complexity of the diagram in this case is unknown. • In Section 3, we describe an algorithm for computing the intersection of a power diagram with the unit sphere. This algorithm uses the exact geometric computation paradigm and can be applied to the four types of unions and intersections. It is optimal for the union and intersection of ellipsoids, but its optimality for the case of intersection of paraboloids is open. • This algorithm is then used for the numerical resolution of the far-field reflector problem.

Using a known optimal transport formulation [21,[START_REF] Glimm | Optical design of single reflector systems and the monge-kantorovich mass transfer problem[END_REF] and similar techniques to [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF], we cast this problem into a concave maximization problem in Theorem 12. This allows us to solve instances with up to 15k paraboloids, improving by several order of magnitudes upon existing numerical implementations [START_REF] Caffarelli | On the numerical solution of the problem of reflector design with given far-field scattering data[END_REF].

Intersection of confocal paraboloids of revolution

Because of their optical properties, finite intersections of solid paraboloids of revolutions with the same focal point play a crucial role in an inverse problem called the far-field reflector problem. This inverse problem is explained in more detail in Section 4. Here we study the computation and complexity of such an intersection when the focal point lies at the origin. We call this type of intersection a paraboloid intersection diagram.

Paraboloid intersection diagram.

A paraboloid of revolution in R d with focal point at the origin is uniquely defined by two parameters: its focal distance λ and its direction, described by a unit vector y. We denote the convex hull of such a paraboloid by P (y, λ). The boundary surface ∂P (y, λ) can be parameterized in spherical coordinates by the radial map u ∈ S d-1 → ρ y,λ (u) u, where the function ρ y,λ is defined by:

(1.1) ρ y,λ : u ∈ S d-1 → λ 1 -y|u .
Given a family Y = (y i ) 1 i N of unit vectors and a family λ = (λ i ) 1 i N of positive focal distances, the boundary of the intersection of the solid paraboloids (P (y i , λ i )) 1 i N is parameterized in spherical coordinates by the function:

(1.2) ρ Y,λ (u) := min 1 i N ρ yi,λi (u) = min y∈Y λ i 1 -y i |u . Definition 1.1.
The paraboloid intersection diagram associated to a family of solid paraboloids (P (y i , λ i )) 1 i N is a decomposition of the unit sphere into N cells defined by:

PI λ Y (y i ) := {u ∈ S d-1 ; ∀j ∈ {1, . . . , N }, ρ yi,λi (u) ρ yj ,λj (u)}.
The paraboloid intersection diagram corresponds to the decomposition of the unit sphere given by the lower envelope of the functions (ρ yi,λi ) 1 i N .

1.2. Power diagram formulation. We show in this section that each cell of the paraboloid intersection diagram is the intersection of a cell of a certain power diagram with the unit sphere. We first recall the definition of a power diagram. Let P = (p i ) 1 i N be a family of points in R d and (ω i ) 1 i N a family of weights. The power diagram is a decomposition of R d into N convex cells, called power cells, defined by

Pow ω P (p i ) := x ∈ R d , ∀j ∈ {1, . . . , N } x -p i 2 + ω i x -p j 2 + ω j . Proposition 1. Let (P (y i , λ i )) 1 i N be a family of confocal paraboloids. One has ∀i ∈ {1, . . . , N } PI λ Y (y i ) = S d-1 ∩ Pow ω P (p i ), where P = (p i ) 1 i N and (ω i ) 1 i N are defined by p i = -(λ -1 i /2)y i and ω i = -λ -1 i -λ -2 i /4. Proof.
For any point u ∈ S d-1 , we have the following equivalence :

u belongs to PI λ Y (y k ) ⇐⇒ k = arg min 1 i N λ i 1 -y i |u ⇐⇒ k = arg max 1 i N λ -1 i -u|λ -1 i y i .
An easy computation gives :

max 1 i N λ -1 i -u|λ -1 i y i = max 1 i N λ -1 i -u + 1 2 λ -1 i y i 2 + u 2 + 1 4 λ -1 i y i 2 = u 2 -min 1 i N u + 1 2 λ -1 i y i 2 -λ -1 i -1 4 λ -2 i .
This implies that a unit vector u belongs to the paraboloid intersection cell PI λ Y (y i ) if and only if it lies in the power cell Pow ω P (p i ). One can remark that the paraboloid intersection diagram does not change if all the focal distances λ i are multiplied by the same positive constant. This implies that the intersection with the sphere of the power cells defined in the above proposition does not change under a uniform scaling by λ (even though the whole power cells change).

1.3. Complexity of the paraboloid intersection diagram in R 3 . In this section, we show that in dimension three, the complexity of the paraboloids intersection diagram is linear in the number of paraboloids. Theorem 2. Let (P (y i , λ i )) 1 i N be a family of solid paraboloids of R 3 . Then the number of edges, vertices and faces of its paraboloid intersection diagram is in O(N ).

The proof of this theorem strongly relies on the following proposition, which shows that each cell PI λ Y (y i ) can be transformed into a finite intersection of discs, and is thus connected. Note that while it is stated only in dimension 3, this proposition holds in any ambient dimension. Proposition 3. For any two solid paraboloids P (y, λ) and P (z, µ), the projection of the set L := (∂P (y, λ)) ∩ P (z, µ) onto the plane {y} ⊥ orthogonal to y is a disc with center and radius c[(y, λ), (z, µ)] = 2λπ y (z)/ y -z 2 r[(y, λ), (z, µ)] = 2 λµ/ y -z , where π y denotes the orthogonal projection on {y} ⊥ . Moreover, given a solid paraboloid P (y, λ), the following map is one-to-one:

F (y,λ) : (S d-1 \ {y}) × R + → {y} ⊥ × R + (z, µ) → (c[(y, λ), (z, µ)], r[(y, λ), (z, µ)]) .
Proof. The proof of the first half of this proposition can be found in [START_REF] Caffarelli | On the regularity of reflector antennas[END_REF], but we include it here for the sake of completeness. We first show that the orthogonal projection onto the plane {y} ⊥ of the intersection L ′ := ∂P (y, λ) ∩ ∂P (z, ν) is a circle. Without loss of generality, we assume that y is the last basis vector (0, . . . , 0, 1). Recall that a paraboloid of revolution ∂P (y, λ) is defined implicitly by the relation x = x|y + λ. Hence, any point x in L ′ belongs to the hyperplane defined by x|z -y = λ -µ. If we denote by z ′ = π y (z), z d = z|y , x ′ = π y (x) and x d = x|y , one has

x d = z ′ |x ′ 1 -z d + µ -λ 1 -z d .
The surface ∂P (y, λ) can be parameterized over the plane {y} ⊥ by the equation

x d = x ′ 2 /2λ- λ/2. Combining this with the relations z ′ 2 + z 2 d = 1 and y -z 2 = 2(1 -z d ), we get x ′ - 2λ y -z 2 z ′ 2 = 4λµ y -z 2 .
We deduce that the projection of L ′ onto the plane {y} ′ is a circle of center c = 2λz ′ / y -z 2 and of radius r = 2 √ λµ/ y -z . Therefore, the projection π y (L) is either the disc enclosed by this circle or its complementary. In order to exclude the latter case, we remark that the intersection P (y, λ) ∩ P (z, µ) is a compact set, because it is convex and does not include a ray (assuming y = z). Hence, the projection π y (L) ⊆ π y (P (y, λ) ∩ P (z, µ)) is also compact, and therefore it is the disc of center c and radius r.

Let us now show that the map F (y,λ) is one-to-one. For a fixed positive µ, let c(z) = c([y, λ], [z, µ]). For every point z in S 2 \ {±y}, denote π 1 y (z) = π y (z)/ π y (z) . This point belongs to the unit circle in {y} ⊥ , which coincides with the equator E y of the sphere S 2 which is equidistant to the points {±y}. Then, given any constant-speed geodesic z(.) such that z(±1) = ±y and such that z(0) = z 0 ∈ E y , i.e., z(t) = sin(tπ/2)y + cos(tπ/2)z 0 , the following formula holds

c(z(t)) = 2λ π y (z(t)) y -z(t) 2 z 0 = λ cos(tπ/2) 1 -sin(tπ/2) z 0 .
One easily checks that the function t ∈ [-1, 1) → cos(tπ/2) 1-sin(tπ/2) is increasing and maps [-1, 1) to [0, +∞). The mapping z ∈ S 2 \ {y} → c(z) ∈ {y} ⊥ thus transforms bijectively every geodesic arc joining the points -y and y into a ray joining the origin to the infinity on the plane {y} ⊥ , and is therefore bijective. From the bijectivity of c and the formula defining the radius, one deduces that the map F (y,λ) is one-to-one.

Proof of Theorem 2. Proposition 3 implies that the projection of the set

L i = {ρ Y,λ (u)u, u ∈ PI λ Y (y i )} = (∂P (y i , λ i )) ∩ j =i P (y j , λ j )
onto the plane orthogonal to y i is a finite intersection of discs, and therefore convex. Since the surface ∂P (y i , λ i ) is a graph over the plane {y i } ⊥ , we deduce that L i is connected. This implies that its radial projection on the sphere, namely the paraboloid intersection cell PI λ Y (y i ), is also connected. We denote by V (resp. E, F ) the number of vertices (resp. edges, faces). Since, each cell PI λ Y (y i ) is connected, the number of faces F is bounded by N . Moreover, since there are at least three incident edges for each vertex, we have that 3V 2E. Then, by Euler's formula, we have that V 2F -4 and E 3F -6.

Even though the complexity of the paraboloid intersection diagram is O(N ), it can not be computed faster than Ω(N log N ), as stated in the proposition below. We first define the genericity condition used in the statement of this proposition.

Definition 1.2. A family of solid paraboloids (P (y

i , λ i )) 1 i N in R 3 is called in generic posi- tion if for any subset (i k ) 1 k 4 of {1, . . . , N }, the intersection 1 k 4 ∂P (y i k , λ i k ) is empty.
Remark that the intersection of four paraboloids (∂P (y i1 , λ i k )) 1 k 4 contains a point x if and only if the projection u = x/ x of this point on the unit sphere satisfies the equations

u -p i1 2 + ω i1 = . . . = u -p i4 2 + ω i4
, where the points (p i ) and the weights (ω i ) are defined by Proposition 1. The genericity condition is then equivalent to the condition that for any quadruple of weighted points (p i k , ω i k ) 1 k 4 , the weighted circumcenter does not lie on S 2 .

Proposition 4. The complexity of the computation of the paraboloid intersection diagram is Ω(N log(N )) under the algebraic tree model, and even under an assumption of genericity.

Proof. We take a family of N real numbers (t i ) 1 i N . For every i ∈ {1, . . . , N }, we put λ i = 1 and y i = ϕ(t i ) , where the map ϕ : R → R 3 defined by ϕ(t) = ( t 2 -1 1+t 2 , 2t 1+t 2 , 0) is a parameterization of the equator S 2 ∩ {z = 0} from which we removed the point (1, 0, 0). The family of paraboloids (P (y i , λ)) 1 i N is such that every cell PI λ Y (y i ) is delimited by two half great circles between the two poles, each of these half circles being shared by two cells. We add the points y N +1 = (1, 0, 0), y N +2 = (0, 1, 0), y N +3 = (-1, 0, 0), y N +4 = (0, 0 -1, 0), y N +5 = (0, 0, 1) and y N +6 = (0, 0, -1), so that (P (y i , λ)) 1 i N +6 is in general position. More precisely, the four points y N +1 , . . . , y N +4 are added to ensure that every points of the equator is at a distance strictly less than √ 2/2 from {y 1 , . . . , y N +4 }. The cells of the two poles y N +5 and y N +6 then do not intersect the equator and we keep the property that there exists a cycle with the N + 4 vertices of {y 1 , . . . , y N +4 } in the dual of the paraboloid intersection diagram. Finding this cycle then amounts to sorting the values (t i ) 1 i N +4 . The conclusion holds from the fact that a sorting algorithm has a complexity Ω(N log(N )) under the algebraic tree model.

Other types of union and intersections

Other quadrics, such as the ellipsoid of revolution, or one sheet of a two-sheeted hyperboloid of revolution can also be parametrized over a unit sphere by the inverse of an affine map [19]. In this section, we study the combinatorics of the intersection of solid ellipsoids, one of whose focal points lie at the origin. Note that this intersection naturally appears in the near-field reflector problem, where one wants to illuminate points in the space instead of directions (as in the farfield reflector problem) [START_REF] Kochengin | Determination of reflector surfaces from near-field scattering data[END_REF]. Furthermore, for both the ellipsoid and the paraboloid cases, we also study the union of the convex hulls. We show that in these three cases, the combinatorics is still given by the intersection of a power diagram with the unit sphere. However, the complexity Figure 1. A flower might be higher. We show that there exists configuration of N ellipsoids whose intersection and union have complexity Ω(N 2 ). An algorithm that matches this lower bound is provided in Section 3.

2.1. Union of confocal paraboloids of revolution. The union of a family of solid paraboloids (P (y i , λ i )) 1 i N is star-shaped with respect to the origin. Moreover, its boundary can be parameterized by the radial function u ∈ S d-1 → (max 1 i N ρ yi,λi (u)) • u. The paraboloid union diagram is a decomposition of the sphere into cells associated to the upper envelope of the functions (ρ yi,λi

) 1 i N : PU λ Y (y i ) := {u ∈ S d-1 ; ∀j ∈ {1, . . . , N }, ρ yi,λi (u) ρ yj,λj (u)
}. As before, these cells can be seen as the intersection of certain power cells with the unit sphere.

Proposition 5. Given a family (P (y i , λ i )) 1 i N of solid paraboloids, one has for all i,

PU λ Y (y i ) = S d-1 ∩ Pow ω P (p i )
, where the points and weights are given by p i = 1 2 λ -1 i y i and

ω i = λ -1 i -1 4 λ -2 i .
Proposition 3 implies that for every i, the projection of L i = (∂P (y i , λ i ))∩∂ 1 j N P (y j , λ j ) onto the plane orthogonal to y i is a finite intersection of complements of discs. In particular, this set does not need to be connected in general, and neither does the corresponding cell. This situation can happen in practice (see Proposition 6). Consequently, one cannot use the connectedness argument as in the proof of Theorem 2 to show that the paraboloid union diagram has complexity O(N ) in dimension three. Actually, the following proposition shows that a unique cell may have Ω(N ) distinct connected component. Proposition 6. One can construct a family of paraboloids (P (y i , λ i )) 0 i N such that the paraboloid union cell PU λ Y (y 0 ) has Ω(N ) connected components. Proof. Let y 0 be an arbitrary point on the sphere, and let λ 0 = 1. Now, consider a family of disks D i in the plane H = {y 0 } ⊥ with centers and radii (c i , r i ) 1 i N , and such that the set

U = N i=1 (H \ D i ) = H \ N i=1 D i
has Ω(N ) connected components. This is possible by setting up a flower shape (see Figure 1), i.e., D 1 is the unit ball and D 2 , . . . , D N are set up in a flower shape around D 1 . By the second part of Proposition 3, one can construct paraboloids (P (y i , λ i )) 1 i N such that F (y0,λ0) (y i , λ i ) = (c i , r i ). Then, the first part of Proposition 3 shows that the paraboloid union cell PU λ Y (y 0 ) is homeomorphic to U , and has therefore Ω(N ) connected components.

One can also underline that the complexity of each cell is O(N ). This is a direct consequence of Proposition 3 and the fact that the complexity of the union of

N planar discs is O(N ) [2, Lemma 1].
2.2. Intersection and union of confocal ellipsoids of revolution. An ellipsoid of revolution whose one focal point lies at the origin is characterized by two other parameters: its second focal point y and its eccentricity e in (0, 1). We denote the convex hull of such an ellipsoid of revolution E(y, e). The surface ∂E(y, e) of this set is parameterized in spherical coordinates by the function Note that the value d is fully determined by e and y and is introduced only to simplify the computations. Let Y = (y i ) 1 i N be a family of distinct points in R d and e = (e i ) 1 i N be a family of real numbers in the interval (0, 1). The boundary of the intersection of solid ellipsoids di -

1 i N E(y i , e i ) is
e 2 i 4d 2
i .

The theorem below shows that in dimension three the complexity of these diagrams can be quadratic in the number of ellipsoids. This is in sharp contrast with the case of the paraboloid intersection diagram, where the complexity is linear in the number of paraboloids.

Theorem 8. In R 3 , there exists a configuration of confocal ellipsoids of revolution such that the number of vertices and edges in the ellipsoid intersection diagram (resp. the ellipsoids union diagram) is Ω(N 2 ).

The proof of this theorem strongly relies on the following lemma, that transfers the problem to a complexity problem of power diagrams intersected by the unit sphere. Lemma 9. Let (p i , ω i ) 1 i N be a family of weighted points such that no point p i lie at the origin. Then, there exists:

(i) a family of ellipsoids (E(y i , e i )) 1 i N , such that for all i, EI e Y (y i ) = S d-1 ∩ Pow ω P (p i ). (ii) a family of ellipsoids (E(y i , e i )) 1 i N , such that for all i, EU e Y (y i ) = S d-1 ∩ Pow ω P (p i ).

Proof. We prove only the first assertion, the second one being similar. Let i ∈ {1, . . . , N }. By inverting the equations

p i = -ei 2di yi yi and ω i = 1 di + e 2 i 4d 2 i
, we get

y i = - 4
(ω i + p i 2 ) 2 -4 p i 2 p i and e i = -2 p i ω i + p i 2 .
The condition e i ∈ (0, 1) is equivalent to ω i < -p i 2 and ω i < -p i 2 -p i . These inequalities can always be satisfied by substracting a large constant to all the weights ω i , an operation that does not change the power cells.

We recall that the Voronoi diagram of a point cloud (p i ) 1 i N of R 3 is the decomposition of the space in N convex cells defined by

Vor P (p i ) := {x ∈ R 3 , x -p i x -p j }.
The next proof is illustrated by Figure 2.

Proof of Theorem 8. Thanks to Lemma 9, it is sufficient to build an example of Voronoi diagram whose intersection with S 2 has a quadratic number of edges and vertices. We can consider without restriction that N = 2k is even. Let ε > 0 be a small number. We let p 1 , . . . , p k be k points uniformly distributed on the circle centered at the origin in the plane {z = 0}, and with radius 2 -ε. We also consider k evenly distributed points q 1 , . . . , q k on (A, B) \ {O}, where A = (0, 0, -ε/4), B = (0, 0, ε/4) and O = (0, 0, 0). Our goal is now to show that for any i and j in between 1 and k, the intersection Vor P (p i ) ∩ Vor P (q j ) ∩ S 2 is non-empty, where P = {p i } ∪ {q j }. We denote by m i,j the unique point which is equidistant from p i and q j , and which lies in the horizontal plane containing q j and in the (vertical) plane passing through p i , q j and the origin. A simple computation shows that the distance δ i,j between m i,j and q j satisfies δ i,j 1 -ε/2 + ε 2 /(64 -32ε). Taking ε small enough, this implies that the point m i,j belongs to the unit ball. The line passing through m i,j and orthogonal to the plane passing through p i , q j and the origin cuts the unit sphere at two points r ± i,j . These two points belong to the same horizontal plane as q j , and they both belong to Vor P (q j ). Moreover, the distance between m i,j and each of these two points is less than (1 -δ 2 i,j ) 1/2 . We can take ε small enough so that this distance is less than sin(π/k). This implies that the points r ± i,j also belong to the Voronoi cell Vor(p i ) and therefore to the intersection Vor P (p i ) ∩ Vor P (q j ) ∩ S 2 . Thus, the intersection between the Voronoi diagram and the sphere has at least k 2 = 1 4 N 2 spherical edges.

Computing the intersection of a power diagram with a sphere

In this section, we describe a robust and efficient algorithm to compute the intersection between the power diagram of weighted points (P, ω) and the unit sphere. We call this intersection diagram. Given the set of weighted points (p i , ω i ), we define the ith intersection cell as ID ω P (p i ) = S 2 ∩ Pow ω P (p i ). This algorithm is implemented using the Computational Geometry Algorithms Library (CGAL) [1], and bears some similarity to [START_REF] Boissonnat | On the combinatorial complexity of euclidean Voronoi cells and convex hulls of d-dimensional spheres[END_REF][START_REF] Delage | Cgal-based first prototype implementation of moebius diagram in 2D[END_REF]. We concentrate our efforts in obtaining the exact combinatorial structure of the diagram on the sphere for weighted points whose coordinates and weight are rational. Exact constructions can be obtained, but are not our focus, since in our applications, the results of this section are used as input of a numerical optimization procedure; see Section 4.

Note that the intersection diagram has complex combinatorics. If one assumes that the weighted points are given by a paraboloid intersection diagram, Proposition 3 shows that the intersection cells are connected. But even in this case, the intersection between two intersection cells ID ω P (p i ) and ID ω P (p j ) can have multiple connected components. Consequently, one cannot hope to be able to reconstruct these cells from the adjacency graph (i.e., the graph that contains An Ω(N 2 ) construction. From a total of N sites, half of them in a segment inside the sphere, and half of them on a circle outside the sphere, we obtain a total of Ω(N 2 ) edges embedded on the sphere.

the points in P as vertices, and where two vertices are connected by an arc if the two cells intersect in a non-trivial circular arc), even in this simple case. In general, the cells of the intersection diagram can be disconnected and they can have holes, as shown in Figure 3. In the next paragraph, we propose an algorithm to compute a boundary representation of these cells. Recall the following definitions concerning power diagrams: Definition 3.1. When two power cells Pow ω P (p i ) ∩ Pow ω P (p j ) intersect in a 2-dimensional set, that is a (possibly unbounded) polygon in R 3 , they determine a facet of the power diagram. Similarly, a ridge is a 1-dimensional intersection Pow ω P (p i ) ∩ Pow ω P (p j ) ∩ Pow ω P (p k ) between three power cells. A ridge can be either a ray or a segment. The boundary of each facet can be written as a finite union of ridges.

In the remainder of this section, we explain how to compute a single intersection cell ID ω P (p i ), for a given p i , without any assumption on the points and weights. This intersection cell can be quite intricate (multiple connected components, holes), and will therefore be represented by its oriented boundary. More precisely, the boundary ID ω P (p i ) is a finite union of closed curves called cycles. These cycles are oriented so that for someone walking on the sphere following a cycle, the intersection cell ID ω P (p i ) lies to the right. The cell is uniquely determined by its oriented boundary. The following theorem is the main result of this section. From Theorem 8, the size of the output in the worst case for the intersection (or union) of confocal ellipsoids is Ω(N 2 ). This implies that the algorithm described below is optimal for computing an ellipsoid intersection (or union) diagram. The optimality of the algorithm for the paraboloid intersection diagram is an open problem, which can be phrased as follows: Moreover the resulting cells are connected. Adopting the randomized incremental construction paradigm (see [20,[START_REF]Computational geometry -an introduction through randomized algorithms[END_REF] for a comprehensive study), this means that the expected complexity of the cell of the ith inserted paraboloid is O(N )/i. Besides, thanks to the connectedness of the cells, we can detect the conflicts with a breadth-first search in linear time as well, without any additional structure. Therefore, the total expected cost is n i=1 O(N )/i, which is in O(N log N ). Notice that this analysis relies on the fact that the cells are connected, which is not necessarily true in the case of union of paraboloids and in the ellipsoids cases. In our application, we always use the generic and robust algorithm based on the intersections of power cells, using CGAL.

3.1.

Predicates. Geometric algorithms often rely on predicates, i.e., estimation of finite-valued geometric quantities such as the orientation of a quadruple of points in 3D, or the number of intersections between a straight line and a surface. These predicates need to be evaluated exactly: if this is not the case, some geometric algorithm may even not terminate [START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF]. We adopt the dynamic filtering technique for the computation of our predicates in CGAL. This means that for every arithmetic operation, we also compute an error bound for the result. This can be automatized using interval arithmetic [START_REF] Brónnimann | Interval arithmetic yields efficient dynamic filters for computational geometry[END_REF]. If the error bound is too large to determine the result of the predicate, the arithmetic operation is performed again with exact rational arithmetic, allowing us to return an exact answer for the predicate. Our algorithm requires the following three predicates:

(i) has on(Point p) returns +1, 0, or -1, if the point is outside, on, or inside the unit sphere respectively. (ii) number of intersections(Ridge r) returns the number of intersections between r and the unit sphere (0, 1 or 2), where r can be a segment or a ray. By convention, segment or ray touching the sphere tangentially has two intersection points, with same coordinates. Remark that the result of the first predicate can be obtained directly from the intermediate computations of this predicate. As we use them intertwined, this shortcut is adopted in our implementation. (iii) plane crossing sphere(Facet f) returns the number of intersections between the supporting plane of f and the sphere (0, 1 or +∞). One may notice that this predicate is just the dual of the first one.

Even though the algorithm below is presented as working directly on the cells of the power diagram, our actual implementation uses CGAL, and we work with the regular triangulation of the set of points (that is, the dual of the power diagram). This triangulation is a simplicial complex with 0-, 1-, 2-and 3-simplices. In CGAL, only the 0-and 3-simplices are actually stored in memory, but it remains possible to traverse the whole structure in linear time in the number of simplices. In order to handle the boundary simplices, CGAL takes the usual approach of adding a point at infinity to the initial set of points. In practice, this means the following. A ridge of the power digram can be either a segment or a ray, and is dual of a 2-simplex in the triangulation. In CGAL, such a ridge corresponds to two adjacent tetrahedra, one of whom may contain the point at infinity. The predicates above need to be adapted in order to handle k-simplices that are incident to the point at infinity. The general predicates have been implemented, but the details are omitted here.

3.2.

Computation of the oriented boundary. The output of our algorithm is the oriented boundary of the intersection cell ID ω P (p i ). Note that this is a purely combinatorial object, and no geometric construction are performed during its computation. It is described as a sequence of vertices on the sphere, oriented arcs, and oriented cycles.

• A vertex on the sphere is the result of an intersection between a ridge of the power diagram and the unit sphere. A vertex on the sphere is not constructed explicitly, but stored as an ordered pair of extremities of the corresponding ridge. Since a ridge [x 1 , x 2 ] can intersect the sphere twice, the order determines which intersection point to consider, i.e., the pair (x 1 , x 2 ) describes the intersection point that is closest to x 1 . This convention allows us to handle infinite ridges, and to uniquely identify vertices on the sphere. In degenerate cases, vertices may correspond to different ordered pairs, but having the same coordinates.

• An oriented arc corresponds to the (oriented) intersection between Pow ω P (p i ), another cell Pow ω P (p j ) and the unit sphere. Two situations arise. Either the arc is a complete circle, in which case it is also a cycle. Or the two extreme points are vertices as defined above, that is, they are the intersection between ridges and the unit sphere. An arc is oriented so that, when walking on the sphere from the first point to the second, the cell lies on the right.

• An oriented cycle is a connected component of the oriented boundary of ID ω P (p i ). A cycle can have two or more vertices on the sphere, in this case, they are represented by a cyclic sequence of arcs. However, a cycle can also have no vertices, and in such a case, it is represented by a single full arc. The orientation of a cycle is given by the orientation of its arcs.

Notice that the boundary of a power cell is composed of several convex facets, that can possibly be unbounded. The intersection of such facets with the unit sphere gives the arcs of the intersection diagram. Our algorithm constructs the oriented boundary of the intersection cell ID ω P (p i ) by iterating over each facet of the power cell Pow ω P (p i ) and obtaining implicitly both the set of vertices on the sphere V and the set of oriented arcs E. The oriented graph G = (V, E) is called the boundary graph. Once G is constructed, the oriented cycles in the boundary ∂ ID ω p (p i ) coincide with the (oriented) connected components of G, and can be obtained by a simple traversal. In the next paragraph, we explain the construction of G = (V, E).

3.3.

Computation of the boundary graph. We start by discarding the facets that do not contribute to the intersection diagram. We iterate over all facets f of Pow ω P (p i ), and keep f only if the predicate plane crossing sphere(f ) returns ∞. (The case where it returns 1 corresponds to a trivial intersection between the facet and the sphere and can be safely ignored.) Assuming that f is not discarded, the intersection of its supporting plane with the unit sphere is a circle in 3D. If all vertices of f are on or inside the sphere (i.e., has on 0), then f is inside or only touching the circle, and is discarded. See Figure 4(a).

Given a facet f = Pow ω P (p i ) ∩ Pow ω P (p j ) that was not discarded previously, we build a sequence of vertices on the sphere by traversing f in a clockwise sequence of ridges, oriented with respect to the vector p j -p i1 . We distinguish two different kinds of vertices on the sphere, the source vertices and the target vertices, as shown on Figure 4(b). Given an oriented ridge r = [x 1 , x 2 ], we compute the value of number of intersections(r). If the ridge has zero intersections with the sphere, we skip it. If there is only one intersection between r and the unit sphere, the corresponding vertex is a source if x 1 is inside the sphere (i.e., has on(x 1) 0) and it is a target vertex in the other case. Finally, if the ridge intersects the sphere twice, the oriented pair (x 1 , x 2 ) describes the target vertex, and the pair (x 2 , x 1 ) is the source vertex. The arcs of the boundary graph G corresponding to this facet are obtained by matching each source vertex with a target vertex using the same cyclic sequence. They are represented by an ordered pair of vertex indices.

Finally, we need to consider the case where none of the ridges of f intersect the sphere. There are two possibilities: either the facet f encloses the circle π ∩ S 2 , where π is the affine plane spanned by f , or the circle is completely outside from f . In the former case, we obtain an arc without vertex in the boundary of ID ω P (p i ), as shown in Figure 4(c). Its orientation is clockwise with respect to the vector p j -p i . To detect whether this happens, we simply need to determine whether the center of the circle lies inside or outside the convex polygon f . This is a classical routine that can be performed using signed volumes. For this test, we only need to use rational numbers because the center of the circle is the projection of the origin on a plane described with rational coordinates.

Proof of Theorem 10. First we build the regular triangulation of the set of N weighted points in O(N log N + C) time. Running the algorithm described in this section for every p i makes us visiting each 1-simplex, 2-simplex, and 3-simplex of the primal no more than two, three, and four times respectively; and for each visit, it takes O(1) time, since we compute constant-degree predicates. Then we visit each generated cycle at most twice. Therefore, the complexity of our algorithm is in O(N log N + C + D), where D is the complexity of the diagram on the sphere. The fact that D = O(C) completes the proof.

Application to Minkowski-type problems involving paraboloids

This section deals with the far-field reflector problem [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF], which is a Minkowski-type problem involving intersection of confocal paraboloids of revolution. Consider a point source light located at the origin O of R d that emits lights in all directions. The intensity of the light emitted in the direction x ∈ S d-1 is denoted by ρ(x). Now, consider a hypersurface R of R d . By Snell's law, every ray x emitted by the source light that intersects R at a point whose normal vector is n is reflected in the direction y = x -2 x|n n. This means that after reflection on R, the distribution of light at the origin given by ρ is transformed into a distribution ν R on the set of directions at infinity, a set that can also be described by the unit sphere. The far-field reflector problem consists in the following inverse problem: given a density ρ on the source sphere S 2 and a distribution µ on the sphere at infinity, the problem is to find an hypersurface R such that µ R coincides with µ. For computational purposes, we assume that the target measure ν is supported on a finite set of directions Y := (y i ) 1 i N , and it is thus natural to consider a reflector made of pieces of paraboloids with directions Y . This problem is numerically solved using an optimal transport formulation due to [21,[START_REF] Glimm | Optical design of single reflector systems and the monge-kantorovich mass transfer problem[END_REF].

4.1. Far-field reflector problem. We consider a finite family Y = (y 1 , . . . , y N ) of unit vectors that describe directions at infinity, non-negative numbers (α i ) 1 i N such that 1 i N α i = 1 and a probability density ρ on the unit sphere. The far-field reflector problem consists in finding a vector of non-negative focal distances (λ i ) 1 i N , such that (FF) ∀i ∈ {1, . . . , N }, ρ(PI λ Y (y i )) = α i , where for a subset X of the sphere, ρ(X) := X ρ(u)du is the weighted area of X. The far-field reflector problem can be transformed into the maximization of a concave functional, combining ideas from [21,[START_REF] Glimm | Optical design of single reflector systems and the monge-kantorovich mass transfer problem[END_REF][START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]. Theorem 12. A vector of focal distances (λ i ) 1 i N solves the far-field reflector problem (FF) if and only if the vector (γ i ) 1 i N defined by γ i = log(λ i ) is a global maximizer of the following C 1 concave function:

(4.5) Φ(γ) := N i=1 PI exp γ Y (yi) (c(u, y i ) + γ i )ρ(u)du - N i=1 γ i α i
where c(u, v) := -log(1 -u|v ) and with the convention log(0) = -∞. Moreover, the gradient of the function Φ is given by

(4.6) ∇Φ(γ) := (ρ(PI exp(γ) Y (y i )) -α i ) 1 i N .
There is a similar formulation for a reflector defined by the union of solid confocal paraboloids [START_REF] Glimm | Optical design of single reflector systems and the monge-kantorovich mass transfer problem[END_REF]. However, there is no known variational formulation for the reflector problem involving intersection of ellipsoids. A numerical approach has been proposed in [START_REF]A rigorous method for synthesis of offset shaped reflector antennas[END_REF]. We now turn to the proof of Theorem 12. This proof combines the results of Section 5 of [3] with the optimal transport formulation of the far-field reflector problem [21,[START_REF] Glimm | Optical design of single reflector systems and the monge-kantorovich mass transfer problem[END_REF]. Let us first recall some properties of supdifferential of concave functions. Given a function Φ and λ in R N , the supdifferential of Φ at λ, denoted

(a) Initial (λ i ) 1 i N (b) Final (λ i ) 1 i N (c) Rendering (d) Final (λ i ) 1 i N (e) Reflector (f) Rendering
∂ + Φ(λ) is the set of vectors v such that ∀κ ∈ R N , Φ(κ) Φ(λ) + κ -λ|v .
A function Φ is concave if and only if for every λ, the supdifferential ∂ + Φ(λ) is nonempty. If this is the case, Φ is differentiable almost everywhere, and at points of differentiability the supdifferential ∂ + Φ(λ) coincides with the singleton {∇Φ(λ)}. Finally, λ is a global maximum of Φ if and only if ∂ + Φ(λ) contains the zero vector.

Proof of Theorem 12. We consider a vector γ in R N . First remark that c(u, y i )+γ i = log(exp(γ i )/(1-u|y i ). Therefore,

u ∈ PI exp(γ) Y (y i ) ⇐⇒ ∀j ∈ {1, . . . , N }, exp(γ i ) 1 -u|y i exp(γ j ) 1 -u|y j ⇐⇒ ∀j ∈ {1, . . . , N }, c(u, y i ) + γ i c(u, y j ) + γ j
Therefore, the function Φ can be reformulated as follows:

Φ(γ) = S d-1 min 1 i N c(u, y i ) + γ i ρ(u)du - N i=1 γ i α i .
We now define T γ as the function that maps a point u on the unit sphere to the point y i such that u belongs to PI

exp(γ) Y (y i ). Then, (4.7) Φ(γ) = S d-1 (c(u, T γ (u)) + γ Tγ (u) )ρ(u)du - N i=1 γ i α i .
Moreover, for any κ in R d , one has min 1 i N c(u, y i ) + κ i c(u, T γ (u)) + κ Tγ (u) . Integrating this inequality, and substracting (4.7), we get

Φ(κ) -Φ(γ) S d-1 (κ Tγ (u) -γ Tγ (u) )ρ(u)du - 1 i N (κ i -γ i )α i 1 i N PI exp γ Y (yi) ρ(u)du -α i (κ i -γ i ) = DΦ(γ)|κ -γ
where DΦ(γ) := (ρ(PI

exp(γ) Y (y i )) -α i ) 1 i N . (4.8)
The above inequality shows that DΦ(γ) lies in ∂ + Φ(γ), i.e., this set is never empty and the function Φ is concave. Since the vector DΦ(γ) depends continuously on γ, we deduce that the function Φ is C 1 smooth, and that ∇Φ(γ) = DΦ(γ) everywhere. By Equation (4.8), a vector λ := exp(γ) solves the far-field reflector problem (FF) if and only if DΦ(γ) = 0, i.e., if and only if γ is a global maximizer of Φ. 4.2. Implementation details. The implementation of the maximization of the functional Φ follows closely [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF]. We rely on a quasi-Newton method, which only requires being able to evaluate the value of Φ and the value of its gradient at any point γ, as given by Equations (4.5)-(4.6). The computations of these values are performed in two steps. First, we compute the boundary of the paraboloid intersection cells PI exp(γ) Y (y i ), using the algorithm described in Section 3. These cells are then tessellated, and the integrals in Equations (4.5)-(4.6) are evaluated numerically using a simple Gaussian quadrature. In the experiments illustrated in Figure 5, we constructed the measure i α i δ yi so as to approximate a picture of Gaspard Monge (projected on a part of the half-sphere S 2 + := S 2 ∩ {z 0}). The density ρ is constant in the half-sphere S 2 -and vanishes in the other half. To the best of our knowledge, the only other numerical implementation of this formulation of the far-field reflector problem has been proposed in [START_REF] Caffarelli | On the numerical solution of the problem of reflector design with given far-field scattering data[END_REF]. The authors develop an algorithm, called Supporting paraboloids which bears resemblance to Bertsekas' auction algorithm for the assignment problem [START_REF] Bertsekas | Dual coordinate step methods for linear network flow problems[END_REF]. They use it to solve the far-field reflector problem with 19 paraboloids. Using the quasi-Newton approach presented above, and the algorithm developed in Section 3, we are able to solve this problem for 15,000 paraboloids in less than 10 minutes on a desktop computer. Note that the algorithm of Section 3 would probably also allow a faster and robust implementation of the Supporting paraboloids algorithm [START_REF] Caffarelli | On the numerical solution of the problem of reflector design with given far-field scattering data[END_REF].

Conclusion

In addition to the open problems mentioned earlier, let us mention a perspective. In a recent article [START_REF] Kitagawa | An iterative scheme for solving the optimal transportation problem[END_REF], the algorithm of supporting paraboloids was extended to optimal transport problems involving a cost function c that satisfies the so-called Ma-Trudinger-Wang regularity condition. For this algorithm to be practical, one needs to compute the generalized Voronoi cells efficiently, defined for any function ψ : Y → R by Vor ψ c (y) = {x ∈ X; ∀z ∈ Y, c(x, y) + ψ(y) c(x, z) + ψ(z)}. For general costs, and even in 2D, one cannot hope to do this in time below Ω(N 2 ). However, the MTW regularity condition ensures an analog of Proposition 3 and in particular, it implies that these generalized Voronoi cells are connected. One might wonder whether a randomized iterative construction could be used in this setting to yield a construction in expected time O(N log N ) in 2D. This would open the way to practical algorithms for the resolution of optimal transport problems that are intractable to PDE-based methods. 5.1. Acknowledgements. The authors would like to thank Dominique Attali, Olivier Devillers and Francis Lazarus for interesting discussions. Olivier Devillers suggested the approach used in the proof of the lower complexity bound for intersection of ellipsoids. The first author is supported by grant FACEPE/INRIA,APQ-0055-1.03/12. The second and third author would like to acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-11-BS01-014-01 (TOMMI) and ANR-13-BS01-0008-03 (TOPDATA) respectively.

19.

, A characterization of revolution quadrics by a system of partial differential equations, Proceedings of the American Mathematical Society 138 (2010), no. 11, 4075-4080. 20. J.-R. Sack and J. Urrutia (eds.), Handbook of computational geometry, North-Holland Publishing Co., 2000. 21. X.J. Wang, On the design of a reflector antenna ii, Calculus of Variations and Partial Differential Equations 20 (2004), no. 3, 329-341.

Centro de Informática, Universidade Federal de Pernambuco, Brazil.

Laboratoire Jean Kuntzmann, CNRS/Université de Grenoble, France Laboratoire Jean Kuntzmann, Université de Grenoble, France

  σ y,e (m) := d 1 -e m| y y where d = y (1 -e 2 ) 2e .

  parameterized in spherical coordinates by the lower envelope of the functions (σ yi,ei ) 1 i N . The ellipsoid intersection diagram of this family of ellipsoids is the decomposition of the unit sphere into cells associated to the lower envelope of the functions (σ yi,ei ) 1 i N :(2.3) EI e Y (y i ) := {u ∈ S d-1 ; ∀j ∈ {1, . . . , N }, σ yi,ei (u) σ yj ,ej (u)}.Similarly, the cells of the ellipsoid union diagram are associated to the upper envelope of the functions (σ yi,ei ) 1 i N as follows:(2.4) EU e Y (y i ) := {u ∈ S d-1 ; ∀j ∈ {1, . . . , N }, σ yi,ei (u) σ yj ,ej (u)}. As in the case of paraboloids, the computation of each diagram amounts to compute the intersection of a power diagram with the unit sphere. Proposition 7. Let (E(y i , e i )) 1 i N be a family of solid confocal ellipsoids. Then, (i) The cells of the ellipsoid intersection diagram are given by EI e Y (y i ) = S d-1 ∩ Pow ω P (p i ), where p i = -ei 2di yi yi and ω i = -1 di -The cells of the ellipsoid union diagram are given by : EU e Y (y i ) = S d-1 ∩ Pow ω P (p i ), where p i = ei 2di yi yi and ω i = 1

Figure 2 .

 2 Figure 2. An Ω(N 2 ) construction. From a total of N sites, half of them in a

Theorem 10 .

 10 There is an O(N log N + C) algorithm for obtaining the intersection diagram of a set P of N weighted points in R 3 , where C is the complexity of the power diagram of P .

Figure 3 .Open Question 11 .

 311 Figure 3. Examples of diagram on the sphere. Left: an intersection diagram containing faces with no vertices, faces with holes and faces with two vertices or more. Right: a diagram corresponding to the intersection of a cube and a sphere; there are seven faces, six faces with no vertices, and a face with no vertices and six holes.

Figure 4 .

 4 Figure 4. Obtaining the arcs of ID ω P (p i ). Three distinct cases: (a) inside the circle, (b) intersecting the circle, and (c) enclosing the circle. There is still a case not depict above, f can be outside the circle as well. The arcs with zero length, such as the one in (b) can be removed from the diagram depending on its use.

Figure 5 .

 5 Figure 5. Numerical computation. Calculations were done with N = 1000 paraboloids for the first row and N = 15000 paraboloids for the second row. (a) Paraboloid intersection diagram for an initial (λi) 1 i N . (b,d) Final intersection diagram after optimization. (e) Reflector surface defined by the intersection of paraboloids. (c,f ) Simulation of the illumination at infinity from a punctual light source lighting uniformly S 2 -, using LuxRender, a physically accurate raytracer engine.

We do not need orientation predicates at this point, since, in CGAL, 3-simplices of a triangulation are positively oriented.


Régis Clouard
email: regis.clouard@greyc.ismra.fr

Christine Porquet

Abderrahim Elmoataz

Marinette Revenu

WHY BUILDING KNOWLEDGE-BASED IMAGE SEGMENTATION SYSTEMS IS SO DIFFICULT

Building general-purpose image segmentation systems is a central concern for many researchers. The challenge is of course to propose more and more flexible and reliable systems, the behaviour of which closely depends on the task to perform and the images to process. In this paper, a knowledge-based system for the dynamic building of image segmentation applications is described. This system is based on the monitoring of a library of image segmentation operators. First, the distinctive features of the segmentation domain are studied, in order to find out which problem solving paradigm is the most relevant in this casehierarchical planning -and seven behavioral rules for an architecture dedicated to solving segmentation problems are proposed. Our implementation is based on the Blackboard model. Cooperation and knowledge acquisition issues, which are essential in our approach, are also briefly discussed.

Résumé

La construction de systèmes de segmentation d'images non dédiés est au coeur des préoccupations de nombreux travaux de recherche. L'enjeu est évidemment de proposer des systèmes plus flexibles et plus fiable, dont le comportement dépend étroitement de la tâche à réaliser et des images à traiter. Dans ce papier, nous décrivons un système à base de connaissances pour la construction dynamique d'applications en segmentation d'images, basé sur le principe de pilotage d'une bibliothèque d'opérateurs. Pour cela, nous étudions les caractéristiques du domaine de la segmentation d'images dans le but de définir quel modèle de résolution de problèmes semble le plus adapté à la résolution de problèmes de segmentation d'images et quelles doivent être les règles de comportement qu'il doit respecter. Notre implantation du système est basée sur le modèle du Blackboard. Dans la dernière partie, nous discutons brièvement le modèle de coopération et d'acquisition des connaissances que nous avons retenu.

Introduction

Building general-purpose image segmentation systems is a central concern for many researchers (Crevier (1993)).

The challenge is of course to propose more and more flexible and reliable systems, the behaviour of which closely depends on the task to perform and the images to process. This objective implies special problem solving architectures, enabling to integrate and bring together various pieces of knowledge, either declarative or procedural. Actually, a purely procedural approach cannot be considered because image segmentation is a complex problem, in the sense of systemic. More formally, most image segmentation problems, such as edge detection, region extraction, shape from texture, shape from shading, surface reconstruction, belong to the mathematical class of inverse ill-posed problems [START_REF] Poggio | Computational vision and regularization theory[END_REF]). It means that there exists no unique and stable transformation function (in the sense that it is continuously depending on the observed data), that can build a specific representation of a scene, starting from any kind of observation. A little noise on data can lead to a great variability in the results.

Moreover, image segmentation is in no way an optimization problem, because one cannot exhibit a generic and indisputable cost function that could be minimized. To tackle this kind of problem, one has to reduce the number of acceptable solutions, on the one hand by introducing a priori knowledge on the solution space, and on the other hand by considering segmentation processes as being decomposed into a sequence of sub-problems, that are either well-posed problems, or problems for which classical regularization methods exist.

If there are no general enough algorithms to solve any segmentation problem, on the contrary, many dedicated segmentation problems have been solved in various domains. They contributed to bring into light a great number of specific but efficient algorithms. It is thus possible to constitute a large library of operators and then to envision the building of a segmentation application as the generation of a program, through selection, parameter adjustment and linking of such operators.

In this paper, problem solving models that can be used to build knowledge-based systems for the monitoring of a library of segmentation operators are discussed and such a system is described. First, by studying peculiar features and difficulties of the image segmentation domain, we set forth seven behavioral rules that must be followed when considering suitable problem solving paradigms. Then, through the description of the system we have developed,

we discuss the choices we made according to these seven behavioral rules.

How to monitor a library of segmentation operators

In this section, we are going to examine the peculiar features of the image segmentation domain when its comes to build a system to generate applications by monitoring a library of operators.

2.1

Monitoring a library

2.2

What kind of problem solving paradigm ?

In the domain of image segmentation, there exists a large amount of algorithms that can be coded as operators, together with numerous techniques and strategies that can be represented as pieces of knowledge. Of course, the major challenge is to know when to use them and what their actual effects on images are. The problem of generating segmentation programs by monitoring a library of operators is a complex problem in the sense that "the whole is more than the sum of the parts" [START_REF] Simon | The sciences of the Artificial. the M.I.T press[END_REF]).

In fact, given the properties of operators and the laws of their interaction, it is not a trivial matter to infer the properties of the graph.

Image segmentation possesses some distinctive features related to the images that are processed, to the nature of operators and to their linking.

First, it cannot be likened to a combinatorial problem, because it is difficult to decide when an image can be recognized as a suitable solution. Moreover, one cannot a priori determine parameter values for most operators. In fact, it is well-known that parameter values have a paramount influence on results (e.g. a binarization threshold) and they must be determined with as much precision as possible.

Second, it cannot be likened to a state space searching problem. This kind of model would suppose that descriptions of an initial and a final states are available.

But, owing to their very nature, images contain an enormous amount of noisy and incomplete data and consequently, they cannot, by themselves, constitute the basis of reasoning. This can explain why reasoning paradigms such as Means-End Analysis or Strip-like

Search cannot be considered in our case, all the more because the knowledge on operators is also imprecise and unreliable [START_REF] Matsuyama | Expert Systems for Image Processing: Knowledge-based Composition of Image Analysis Processes[END_REF]). It is very difficult to a priori forecast the modifications induced by the execution of an operator, and it is perhaps even more difficult to define how to measure differences between two states.

One can then resort to reasoning in a plan space [START_REF] Currie | O-Plan: the open planning architecture[END_REF]), the only acceptable paradigm when there is an infinite number of potential configurations. It implies that the reasoning is made on actions and their linking and not on input or output data.

What kind of behavioral rules for solving segmentation problems ?

More formally, in this section we set forth seven behavioral rules for an architecture dedicated to solving segmentation problems.

Rule 1:

Choose a problem solving model based on actions.

Contrary to segmentation states, segmentation actions can be described accurately, and they can thus constitute the basis of reasoning. Segmentation actions represent image segmentation decisions that are applied to input images in order to produce output images. A segmentation problem can be decomposed into a sequence of tasks to perform, and segmentation techniques can then be proposed to perform each of these tasks (Clement & Thonnat (1993)).

Rule 2: Use a reasoning model based on hierarchical planning of actions.

Segmenting an image is not a linear process [START_REF] Gong | VISIPLAN : A Hierarchical Planning Framework for Composing Biomedical Image Analysis Processes[END_REF]). As a matter of fact, the processing of an image cannot generally be decomposed into a linear sequence of atomic operations, for which it is Because knowledge on operators and segmentation techniques are uncertain and unreliable, one has to execute operators so as to obtain intermediate results that can then be assessed in order to check the relevancy of the current plan. Because it is difficult to a priori forecast the precise effects of an operator on a given image, and in the same way it is difficult to find the optimal values of each parameter of this operator [START_REF] Matsuyama | Expert Systems for Image Processing: Knowledge-based Composition of Image Analysis Processes[END_REF]), all the more as the quality of results are very sensitive to the values of its parameters. One has to resort to mechanisms enabling operators themselves to compute their parameter values through optimization or simply satisfaction of some cost function.

Rule 4: Choose a library of operators mainly composed of "atomic" operators.

Atomic operators are operators, the action of which cannot or should not be decomposed further. This rule has two major advantages in the framework of library monitoring [START_REF] Clouard | Resolution of image processing problems by dynamic planning within the framework of the Blackboard model[END_REF]). First, it is a means to reduce the number of parameters to be tuned. Secondly, the relatively small-grain knowledge associated to each operator also makes their selection and use easier.

Moreover, one needs some special operators for combining results either in a logical, spatial or arithmetical way.

Rule 5: Adopt a hierarchical and delocated evaluation of results.

Controlling the execution of operators is not enough to ensure the relevance of the graph. The optimization of each operator does in no way mean that the whole graph will be optimized. Because evaluating the results of a segmentation process is a problem in itself, and because there are no generic evaluation functions, one has to resort to a hierarchical evaluation at each abstraction level of the planning process. So, each action produced at each abstraction level should be associated with its own evaluation functions. Moreover, the results of operators must be brought up towards upper levels and evaluated so as to take into account more and more general criteria.

Rule 6: Describe a segmentation problem as a set of tasks to perform, together with constraints and a symbolic and numeric description of the application context.

Intentions and context must be taken into account at each reasoning level, to direct choices during the planning process, control the execution of operators and select the relevant evaluation rules. Tasks describe the nature of the problem, whereas constraints define the quality of expected results. The description of the context should include three levels of information [START_REF] Elmoataz | Mécanismes opératoires d'un segmenteur d'images non dédié: définition d'une base d'opérateurs et implantation[END_REF]): physics about image formation (e.g. type of camera, acquisition conditions), perceptual information (e.g. edge type, region type, texture...) and knowledge about the semantics of the scene (e.g. shape and size of objects and relations between objects). It is essential to notice here that, in order to remain as general as possible, the knowledge and vocabulary used should be "weeded out" of all technical references to the domain of application. That is the reason why one should restrict to vocabulary in use in image segmentation (regions, boundaries, lines, pixels...) and to mathematical vocabulary (geometry, algebra...) for describing relations and features of objects. On the contrary, notions related to the domain of application (cells, roads...) should not be mentioned just as they are but described only by means of the authorized vocabulary (shape, texture...). This is a very strong and restrictive constraint but it is the only way to exhibit knowledge that is sufficiently general and reusable.

The major interest of using this restricted vocabulary is that it is common to all users. So, it can be the base of the dialogue between them.

Rule 7:

Integrate users in the resolution loop to cooperate or collaborate with the system in order to find an acceptable solution.

Despite all the attention that must be paid to the knowledge acquisition process, it is unrealistic to hope that a generalpurpose segmentation -segmentation system will ever be totally autonomous when trying to solve new segmentation problems in little-known or unknown domains. The only feasible solution is then to make users enter the resolution loop, by letting them cooperate (or collaborate) with the system during the resolution of a problem. The user can be asked by the system when some information about the goal or the context are missing, or when the system declares being unable to solve the given request. On the other hand, the user can interact opportunistically when he/she is not satisfied by the solution proposed by the system.

Conceptual model

We are now going to focus on the choices and principles we have followed to build the conceptual model of our system, in accordance with the seven behavioral rules we have just defined.

3.1

Model of knowledge

We have chosen to study the segmentation domain under five abstraction levels (Table 1) Four types of knowledge are taken into account by the reasoning mechanism:

1. Knowledge about the domain of image segmentation and its context, enabling to understand Request Definition of the problem to be solved, given by the user in terms of goals to be reached and constraints on these goals.

Task

All the primitive tasks to be solved can be found at this level. They are either deduced from the request or from the chosen analysis strategy.

Functionality

Functionalities describe general-purpose segmentation functions that must be implemented in order to solve some task.

Procedure

Procedures correspond to classical segmentation operators, but defined independently from any specific implementation.

Operator

Our run-time library of operators is our specific implementation of procedures. An operator is characterized by its prototype and the type and domain of each of its parameters. the problem data in order to determine a suitable analysis strategy and to raise ambiguities when translating this problem into segmentation tasks.

Tasks are our means of setting image segmentation problems.

2. Knowledge about the expertise in image segmentation, used for determining which techniques will be considered to perform segmentation tasks, and how results will be evaluated.

3. Knowledge about the operators of the library, used for selecting operators and setting the values of their parameters. This knowledge is totally independent from the implementation of the run-time library.

4. Knowledge about the control of the resolution, used to direct choices and to solve resolution conflicts. In most systems, this knowledge takes a procedural form. We, on the contrary, intend to express this knowledge declaratively and explicitly.

Incremental formation of the solution

In our approach, the plan is built step by step, so that it can match the current state of the solution. Each decomposition of a goal is made according to specific constraints and characteristics of the class of images.

The same kind of approach is used for assessing the quality of the plan. The plan evaluation is done in a hierarchical and delocated manner. To each goal are associated rules for assessing its own results. At the highest level, these rules are concerned with semantic features and at the lowest level with syntactic features. This principle enables to organize and distribute rules into the hierarchy.

Moreover, as a consequence of this incremental approach, it is practically impossible to forecast the plan that will come out as final solution.

Opportunistic problem-solving behavior

There are no predefined algorithms to solve a given segmentation problem, because of the size of data and the large number of alternatives to be taken into account. At each step of the resolution, several actions are feasible, so that the solution can make progress: should we decompose some goal into subgoals, assess the result of the execution of an operator, bring back up results of an evaluation or apply the evaluation rules of some other goal? Selecting the "best" action is very important for it influences directly the quality and rapidity of convergence of the solution. This falls within the competence of the control of resolution mechanism.

Choice is done dynamically, according to the state of progress of the current solution. It hinges upon general heuristics and specific focuses of attention. Our model enables to benefit from several resolution modes, either cooperative, or competitive. Several approaches are available: the goal-directed approach (choosing the actions that create important data), the action-directed approach (choosing the intrinsically important actions), the plandirected approach (choosing an action in accordance with the current resolution strategy), either in a forwardchaining manner or in a backward-chaining manner.

To that purpose, one has to consider the data of the solution as hypotheses associated with a coefficient of plausibility and a coefficient of importance towards the current state of the solution. These coefficients are used by the control mechanism to calculate priorities in the development of the solution.

A system that combines the above-mentioned features is general, powerful and flexible and can be directly implemented by means of a Blackboard architecture [START_REF] Nii | Introduction, Blackboard architecture and applications[END_REF]).

Implementation

In this section, we are now going to describe the distinctive features of our system which is based on the Blackboard architecture: the database, the knowledge sources and the way we solve control problems [START_REF] Clouard | Resolution of image processing problems by dynamic planning within the framework of the Blackboard model[END_REF]). The Blackboard model is at once a conceptual, high-level organization of information and knowledge and a general prescription for the dynamic control and use of knowledge for incremental, opportunistic problem solving [START_REF] Nii | Introduction, Blackboard architecture and applications[END_REF]).

Syntax of segmentation requests

The input of our system is a graphical interface enabling the user to specify segmentation tasks. This interface takes the form of a hypertext through which the user can progressively refine the terms of his/her request by selecting tasks, defining constraints and describing the application context by giving values to attributes at the physical, perceptual and semantic levels.

More formally, the request is defined by the following grammar in BNF form :

Database

Goals, input data, results, hypotheses are all stored into a global database, called the blackboard, which is organized vertically following the five levels previously detailed (Fig. 1). A plan is represented as a five-level tree of goals. Each node of the tree is connected to nodes at the next lower level by sequence (THEN), and conjunction (AND) links.

Links between goals can be seen as channels for the transmission of input and output images from one node to another.

Each level is described by a list of attributes where some are common to all five levels (Table 2). Other specific attributes are defined for each level. For instance, nodes at the operator level have a parameter attribute, containing the domain of possible values for each parameter of the operator, as well as a prototype attribute giving the usage of the operator.

Knowledge Sources

The knowledge base is divided into independent and autonomous modules (Knowledge Sources) that totally ignore one another. A knowledge source (KS) contains expertise to solve one part of the global problem in a unique way. So, there are as many KS as way to solve one given problem. A KS constitutes a link between two nodes of the blackboard, one using it as its input and the other as 1. Decomposition KS : They contain the knowledge to decompose some goal into subgoals at the next lower level.

For that purpose, the description of the context and the constraints associated to the goal are used to build a set of a value giving the confidence degree in the present decision related to the current solution.

Table 2: The attributes common to all the five levels.

subgoals and to specify their constraints and the relations existing between subgoals. For each subgoals, a decomposition-KS must define the goal, the constraints, the input image or the path to get them and also the evaluation rules used to check the accuracy of the results with regard to the goals to be reached and the associated constraints.

Decomposition KS are given in the following form : (where T1',T2',T3 are the same tasks as T1,T2,T3 but with different constraints).

Execution KS

Evaluation KS :

There is one evaluation KS per level, which is triggered on a goal of the tree when the decomposition of this goal is completed. Its action part propagates the results of the decomposition of the goal upwards and applies the associated evaluation rules. The result is judged acceptable when output images are in accordance with expectations.

Control blackboard and control Knowledge Sources

The task assigned to control is to select the next KS to be executed. The control structure consists of a control blackboard and a set of control knowledge sources. The Blackboard model is strongly inspired by BB1 [START_REF] Hayes-Roth | A blackboard architecture for control[END_REF]), the control blackboard being decomposed into six following levels (Table 3). The first four levels define, at each step of the resolution, the profile of the desired actions, with regard to the current state of the solution, whereas the two last levels define the profile of the actions that are actually feasible at this step of the resolution (the actions that can contribute to the progress of the solution).

The control problem is solved by control KS, following the same principles as those dealing with the domain of application. There are no differences between control-related decisions and domain-related ones. The choice of the next KS to be executed is the result of a compromise between desired KS and executable KS. The In the framework of planning image segmentation tasks, the general strategy we use is a breadth-first strategy.

It is implemented as a sequence of focuses of attention on successive levels of the database. Specific strategies can also be considered, in order to give first priority to some parts of the graph that are judged important at this step of the resolution (e.g. a temporary depth-first strategy used during the decomposition of procedures into operators). We also can apply general-purpose heuristics to give preference to :

 control actions versus domain actions,  intrinsically important actions, such as the execution of the evaluation KS,  actions that work on important data,  actions that work on data with a high coefficient of plausibility.

 actions that have a low execution cost.

How to cooperate

At present, three types of users can cooperate with the system: experts in the domain of origin of images, experts in image segmentation and the designer-programmer of the system itself. Furthermore, two forms of cooperation are available [START_REF] Clouard | A Software Workshop for Knowledge Acquisition in Image Processing[END_REF]), the first one is essential to the knowledge acquisition process, and the second one to the resolution of applications. Consequently, our system presents two faces, whether the user chooses the "manual"option for knowledge acquisition or the "automatic" one in the case of problem solving. To this purpose, representing an action-based reasoning through a hierarchical plan offers good means of cooperation [START_REF] Willamowski | Modélisation de tâches pour la résolution de problèmes en coopération système-utilisateur[END_REF]) .

Cooperation for knowledge acquisition

The manual option is simply a plan editor, that enables the image segmentation expert to manually build applications as they are prescribed by the domain expert.

The image segmentation expert can thus create his/her own base of plans of actions. Plans are coded in a completely static manner according to the task/method/tool paradigm [START_REF] Revenu | Atelier Logiciel d'acquisition et d'intégration des connaissances en Traitement et Interpretation d'images : Une approche méthodologique[END_REF]). In our case, tools correspond to operators and methods simply correspond to the links between tasks and tools. During this procedural integration process, segmentation techniques are specified as knowledge that is directly operational for the planning module, because the plans created "manually" share the same structure as the plans automatically built by the planning module. Knowledge acquisition thus requires an intermediate stage when one has to specify the set of criteria that lead to the decomposition rules of each node of the plan, i.e. type of decomposition (AND/THEN) and principally, application and non-application conditions for this decomposition. Then, once a first integration stage is over, the system is used as an experimentation tool when several plans related to the same request can be tried and tested, in order to evaluate their performance, to find reasons for their eventual failure or to propose improvements

Get objects from regions

5.2

Cooperation for problem-solving

The automatic option corresponds to a planning process involving cooperative work between the user and the reasoning module. In this automatic mode, interventions are on the system's initiative when data are missing (by asking them to the domain expert), and on the user's initiative when he/she wants to get some explanations or to modify previous choices. When an intervention is on the system's initiative, it is the system's opinion that prevails and in all other cases, the user's opinion is prevailing.

The system only appeals to the user in two specific cases: when an input request is considered as incomplete and when visual human assessment is needed.

Conclusion

Image segmentation is a typical domain where one has to develop systems establishing a cooperation between several high-level expertises. The resolution of a segmentation problem involves three categories of experts: experts from the domain of application (biologists, geographs), experts in image segmentation and designer-programmers. Domain experts are required to express clearly and exhaustively the purpose of their application, and also to provide several criteria that will be used during the a posteriori evaluation of results. Segmentation experts have to propose some kind of answer in terms of plans of actions adapted to the specific features of the request and the nature of images.

Designer-programmers finally have to provide a framework for an efficient implementation.

Rather than aiming at developing a general-purposed system, we are interested in setting up a workbench directed towards the cooperation between various agents involved in the system so as to provide assistance to the working out of some kind of computational theory of segmentation, analogous to Marr's vision theory. This approach is motivated by the development of real-size applications, while being concerned by generic principles. We are addressing issues such as:



The development of concrete applications with classical methods, in order to precise the nature of real problems and to establish an efficient dialogue between segmentation experts and domain experts, that can be based on concrete results.



The specification of a software architecture favoring cooperation. The objective is to get a system that can manage knowledge, i.e. to represent it and make it operational. It must be a framework for the development of applications and for the acquisition of know-hows and knowledge on techniques and computer tools. Although there is still a lot of work to do, the blackboard-based system we have described in this paper, starts to bring satisfactory answers to some of these objectives.

Our research takes place within the framework of the "Image Processing and Analysis" Center of Caen. This center covers a wide domain of applications, dealing with biomedical, as well as material or geological images. It constitutes a privileged field of investigation in order to validate and refine the ideas we are advocating and to make various experts cooperate.

  corresponding to the steps generally considered in image segmentation, starting from the problem specification, up to the selection of operators. This five levels correspond to a more detailed view of the three levels of Rule 2. The first two levels (Request and Task) specify the problem completely, through the determination of all the tasks to be solved. The Request represent the problem set by the user, whereas a Task defines a part of a some segmentation strategy (e.g Eliminate the background in order to extract regions of interest). The next two levels (Functionality and Procedure) are here to determine the various solutions advocated by segmentation experts, in order to solve each task. Functionalities describe the segmentation method that is selected (e.g a pixel classification method for the Task of elimination of the image background), and Procedures precise the tools that implement these functionalities (e.g a binarization based on the boundary contrast tool for the pixel classification method). The last level (Operator) ensures the instantiation of the plan with operators from the library in use. Building the plan of actions is achieved by successive refinements of goals at one level into subgoals at the next lower level. Each level corresponds to a more or less coarse version of the solution and is decomposed into an ordered sequence of more technical subgoals. The decomposition of a goal (Fig 1.) can either be a set of subgoals that must be executed sequentially (THEN links), or a conjunction of subgoals, in which case the execution order is not specified (AND links). Input images to subgoals come either from the decomposed goal or from one of the preceding goals in a sequence. Output images of a goal are those explicitly mentioned in the decomposition of this goal. More concretely, the Request is split up into primitive Tasks by determining an analysis strategy, reformulating and eliminating ambiguities. Tasks are transformed into Functionalities, which, in turn, are broken down into Procedures. Procedures are implemented by means of Operators, for which the values of each parameter must be calculated. In the case of Operators, for which a priori determination of parameter values is difficult, a set of potential values are calculated and the operators are then executed in a trial-and-error process with the different values so as to optimize some evaluation function, according to Rule 3. Images resulting from the execution of operators are brought back upwards within the hierarchy and the output images of the request constitute the final result of the resolution. The user is in charge of the final visual evaluation and can reformulate the request if he/she is not satisfied. This reformulation consists in selecting new tasks and new constraints and must be done jointly by the image segmentation expert and the domain expert, after assessing the first results produced by the system.

Fig. 1 :

 1 Fig.1: Decomposition of a goal into subgoals

  its output (Fig 2). KS are defined in the traditional condition-action style. The condition part is responsible for determining when the KS can contribute to the problem resolution. The action part acts on the solution by creating or modifying data. The knowledge base is composed of three kinds of KS :

  TRIGGER:  a task T / goal(T)=B CONDITION: state(T) solved ACTION: if ( an attribute ycontext / val(y)=V) or ( an attribute zconstraints(T) / z=W) then Create decomposition(T)=(T1,T2,T3) else Create decomposition(T)=(T0', T1',T2',T3'))

:

  There are five KS of this kind according to five execution modes: NORMAL execution, execution of an operator in OPTIMIZATION mode with the various combinations of possible values for each parameter, iterative execution (FOR), and repetitive execution (WHILE or UNTIL).

  profile of the next KS to be executed is determined by general-purpose heuristics and local focuses of attention, resulting from general or specific strategies. Executable KS are put into the agenda. At each cycle, the scheduler has to decide whether to change the current strategy or to continue the development of the current solution on the domain blackboard.

Fig. 2 :

 2 Fig. 2 : An example of a plan on the segmentation blackboard

Table 1 :

 1 The segmentation domain is studied under five levels of abstraction.

Table 3 :

 3 The six hierarchical levels of the control blackboard.