Rita Zrour

Yukiko Kenmochi

Hugues Talbot

Ikuko Shimizu

Akihiro Sugimoto

Combinatorial optimization for fitting digital line and plane

We present a method for fitting a digital line/plane from a given set of 2D/3D grid points. In the framework of discrete geometry, a digital line or plane is defined as a set of grid points located between two parallel lines or planes separated by a small distance. Our purpose is to identify such a pair of Euclidean lines/planes that represents a given set of points. This problem is formulated as a mixed integer/linear programming problem, with the objective to maximize the number of points between the two lines/planes.

I. INTRODUCTION

In the computer vision literature, recognition of straight lines and planes has been studied for decades. Since digital images always contain outliers, most robust recognition algorithms are based on statistical methods and solved by optimization methods using, for example, regression tools such as least squares, least absolute deviation, and least median of squares linear regression [START_REF] Press | Numerical recipes (Third Edition)[END_REF]. Most of these algorithms work in a continuous space. However since the data are digital, it seems logical to work with a digital model rather than a continuous one. Recently, purely discrete formulation of many classical and novel image processing problems has been proposed under the heading of "discrete geometry" [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF]. In this domain, work is accomplished on discrete data within a discrete space. In other words, we don't have any digitization error since we treat only grid points whose coordinates are represented as integers.

In a discrete space, a digital line/plane is defined as a set of grid points located between two parallel lines/planes separated by a small distance [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF], [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF]. In discrete geometry, several algorithms were proposed for recognizing a digital line or plane for a given grid point cloud [START_REF] Debled-Rennesson | Etude et reconnaissance des droites et plans discrets[END_REF], [START_REF] Buzer | An elementary algorithm for digital line recognition in the general case[END_REF]. They are, for example, based on linear programming or computational geometry using convex hulls or other geometrical properties. Those algorithms determine if a point cloud fits to a digital line or plane. Therefore, they are very sensitive to outliers; if a digital image contains outliers, they simply return the negative answer. Recently, new methods for blurred data were introduced in the framework of discrete geometry [START_REF] Debled-Rennesson | Linear segmentation of discrete curves into fuzzy segments[END_REF], [START_REF] Provot | Recognition of blurred pieces of discrete plane[END_REF], [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF], [START_REF] Faure | Tangential cover for thick digital curves[END_REF]; for example [START_REF] Provot | Recognition of blurred pieces of discrete plane[END_REF] consists of minimizing the distance between two parallel lines/planes containing all points of a given point cloud. Therefore, outliers are not separated from a fitted digital line/plane. In other words, there is only notion of inliers so we cannot process outliers with these methods.

Our purpose in this paper is to fit a digital line or plane to a given set of 2D or 3D grid points. This fitting is accomplished by defining inliers and outliers. The work is accomplished in a discrete space so that discretisation errors are processed naturally. The fitting is accomplished by maximizing the inliers and rejecting the outliers. Figure 1 shows an example of fitting a digital line to a given set of points by rejecting the points that are considered as outliers. In order to maximize the inliers we use the L0-norm, unlike statistical methods that use L1-norm or L2-norm. L0-norm considers only the binary distances 0 or 1; inliers are thus at distance of 0 while outliers are at distance of 1. Remark that such simple distances can be used because our model is digital. We formulate the problem into a mixed integer linear programming that aims to find the parameters of the digital line or plane by minimizing the number of outliers. This formulation allows us to find an optimal solution; this constitutes an advantage over statistical methods that work on a continuous non integer space and sometimes yields only approximate solutions.

This paper is organized as follows. Section II presents discrete geometrical definitions of digital lines and planes. Section III proposes the formulation for fitting a set of points to a digital line or plane using the discrete geometrical definitions; we show that our problem is written as mixed integer linear programming problems. Section IV presents some experiments by testing the method. Section V discusses the computational issue. Section VI presents an improvement to consume less computation time. Finally, Section VII states some conclusions and perspectives. Unlike in the continous domain, in the framework of discrete geometry, we take into account digitization errors when we represent straight lines and planes in a digital image. In this section we give the definitions of digital lines and planes found in [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF].

A. Digital lines

Let R be the set of real numbers. A line L in the 2D Euclidean space is defined by:

L = {(x, y) ∈ R 2 : αx + βy + δ = 0} (1)
where α, β, δ ∈ R.

Let Z be the set of integers. Z 2 denotes the set of grid points whose coordinates are all integers in the 2D Euclidean space R 2 . We now consider a digitization of L, which is given as a grid-point subset in Z 2 , called a digital line. There are various digitization techniques for a given L [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF]. In this paper we adopt a grid-line digitization with respect to L. Such a digitization of L is given by

D(L) = {(x, y) ∈ Z 2 : 0 ≤ αx + βy + δ < ω} (2)
where ω = max (|α|, |β|). D(L) is the set of all grid points that are closest to the intersection points of L with the grid lines of Z 2 and that are in the upper/under side of L. From the definition, it is obvious that we obtain a unique digital line D(L) from a given L [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF].

B. Digital planes

A plane P in the 3D Euclidean space is defined by:

P = {(x, y, z) ∈ R 3 : αx + βy + γz + δ = 0} (3)
where α, β, γ, δ ∈ R. Z 3 denotes the set of grid points whose coordinates are all integers in the 3D Euclidean space R 3 . Concerning a digitization of P, which is given as a grid-point subset in Z 3 , called a digital plane, we adopt a grid-line digitization similarly to digital lines. Such a digitization of P is then given by

D(P) = {(x, y, z) ∈ Z 3 : 0 ≤ αx + βy + γz + δ < ω} (4)
where ω = max (|α|, |β|, |γ|). As with discrete line, D(P) is unique for every P [START_REF] Klette | Digital Geometry: Geometric methods for digital pictures analysis[END_REF].

III. DIGITAL LINE/PLANE FITTING

As stated in Section I, given a finite subset S of Z n where n = 2, 3, our purpose is to fit a digital line (2D case) or plane (3D case) to S and to estimate its parameters. In our case, we would like to reject some points considered as outliers of S.

A. Formulation

If S contains no outlier, in 2D, we are simply interested in finding the solution set (α, β, δ) such that all points in S satisfying [START_REF] Kenmochi | Discrete plane segmentation and estimation from a point cloud using local geometric patterns[END_REF]. Similarly in 3D, we are interested in finding the solution set (α, β, γ, δ) such that all points in S satisfying (4).

If S contains outliers, however, we find no solution for any above linear inequality set. For such infeasible cases, we proposed to use integer linear programming to solve the problem. Integer linear programming allows us to separate the inliers from outliers.

The 2D formulation cas be described as follows: let N be the number of elements of S and p i be a binary decision variable for each point (x i , y i ) ∈ S, where i = 1, 2, . . . , N , such that p i = 0 if a point (x i , y i ) satisfies the linear inequalities of a digital line; otherwise, p i = 1. The integer linear programming problem can be described by:

[LP 1-j] minimize i=1,...,N p i subject to -M p i ≤ αx i + βy i + δ < M p i + ω for all i = 1, . . . , N (5) 0 ≤ p i ≤ 1 and integer ( 6
)
where M is a large number constant.

In order to simplify [START_REF] Provot | Recognition of blurred pieces of discrete plane[END_REF], it is possible to divide ( 5) by ω = max (|α|, |β|); this leads us to the following inequalities:

-M ω p i ≤ α ω x i + β ω y i + δ ω < M ω p i + 1 (7)
Looking into [START_REF] Reveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF], we see that four cases can be distinguished:

1) ω = |α|, α ≥ 0 2) ω = |α|, α ≤ 0 3) ω = |β|, β ≥ 0 4) ω = |β|, β ≤ 0
These four cases can be reduced into the following two cases and ( 7) or ( 5) is rewritten as follows:

[LP 1-1] If ω = |β| (j=1), -M p i ≤ α 1 x i + y i + δ < M p i + 1 (8)
where

M = M ω , α 1 = ± α ω , δ = δ ω or -δ ω + 1.
These substitutions yield the following constraint:

-1 ≤ α 1 ≤ 1. [LP 1-2] Otherwise (j=2), -M p i ≤ x i + α 2 y i + δ < M p i + 1 (9)
where

M = M ω , α 2 = ± β ω , δ = δ ω or -δ ω + 1.
These substitutions impose the following constraints:

-1 ≤ α 2 ≤ 1.
Regarding the 3D case, its formulation is similar to the 2D one. Three cases are distinguished depending on the cases where ω = |α|, ω = |β| and ω = |γ|. The integer linear programming is described for the case ω = |α| by:

[LP 2-1] minimize i=1,...,N p i subject to -M p i ≤ x i + α 2 y i + α 3 z i + δ < M p i + 1 for all i = 1, . . . , N (10)
-1 ≤ α 2 ≤ 1, (11) -1 ≤ α 3 ≤ 1, (12) 0 ≤ p i ≤ 1 and integer ( 13
)
where M is a large number constant. If we consider the case where ω = |β|, then we replace [START_REF] Faure | Tangential cover for thick digital curves[END_REF] and ( 11) by the following inequalities respectively:

[LP 2-2] -M p i ≤ α 1 x i + y i + α 3 z i + δ < M p i + 1 (14) -1 ≤ α 1 ≤ 1. (15)
If we consider the case where ω = |γ|, then we replace [START_REF] Faure | Tangential cover for thick digital curves[END_REF] and ( 12) by the following inequalities respectively:

[LP 2-3] -M p i ≤ α 1 x i + α 2 y i + z i + δ < M p i + 1 (16) -1 ≤ α 1 ≤ 1. (17)
It is noticed that we have two formulations in 2D and three formulations in 3D. We may generalize that we have n formulations in nD.

B. Algorithm

We have to solve two linear programming problems [LP 1-1] and [LP 1-2] for the 2D case, and three problems [LP 2-1], [LP 2-1] and [LP 2-3] for the 3D case. Let F j to be the solution of the objective function of the j-th problem [LP n-j] for j = 1, 2, . . . , n where n = 1, 2. Once all the F j s are computed, the minimum among them is chosen such that F = min j=1,2,...,n (F j ). It allows us to detect an optimal solution that yields fewer outliers and fits more points into a digital line or plane.

Algorithm 1 describes the different steps of fitting a digital line and plane to a given set of points for the 2D and 3D cases.

IV. EXPERIMENTS A. 2D example

Our input data is a set of points (x, y) ∈ Z 2 that satisfies y = -6x and another set of points (x, y) ∈ Z 2 that satisfies y = -x + 5. We used the free linear programming solver, lp_solve [11] to solve each [LP n-j] problem.

As mentioned in Section III-B, we compute two F j , j = 1, 2 in 2D and then we choose the minimum between them. Figure 2 shows an example of fitting a digital line to a given set of points; between the two solutions, one is optimal since it fits more points into a digital line. Table I provides an example of a fitting 200 points oscillating between the two discrete lines 0 ≤ x + y + 5 < 1 and the line 0 ≤ x + 0.1667y + 1 < 1. When the points move from one line to another, the numbers of inliers and outliers change and the line detected is the one that contains most inliers.

We also tested our method with respect to a real image as shown in Figure 3. The image dimensions are 520 × 693. Before applying our method, edge detection and mathematical morphological filtering are done for this image; the number of points in the image after this pre-processing is 1572 points. As shown in Figure 4, our method is then applied in order to fit a digital line to the set of points; Figure 5 shows the digital line detected (colored red).

We compared our method with the Hough transform [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] that is one of well-known efficient procedures for detecting lines in images. It consists of transforming each image point (x, y) into a curved line in the parameter space (θ, ρ) specified by: ρ = x cos θ + y sin θ It has been noticed through experiments that with the Hough transform, line detection usually requires interactive parameter adjusting. The parameters adjustement is done by changing the sampling intervals ∆ρ, ∆θ and the voting threshold t. Many lines can be detected if desired. Figures 6,7,8 show the Hough transformation applied to Figure 4 with different parameter values. Depending on the parameter setting, we see in the figures that the results are changed. If the threshold is too high, no line is detected. In contrast with the Hough transform, our method computes exact optimal solution without any parameter adjustement. It is noted that with our method, an iterative procedure can be applied in order to detect many lines in the image. This iterative procedure is accomplished by repeating the procedure for an optimal solution after taking off the points corresponding to the previous optimal line.

B. 3D example

Our method is also easily adapted to the 3D case in contrast to the Hough transform that does not scale well to 3D. The Hough transform requires memory consumption in 3D. For the experiments in 3D, we applied our method to a part of the same data used in [START_REF] Kenmochi | Discrete plane segmentation and estimation from a point cloud using local geometric patterns[END_REF]. Figure 9 shows a planar surface segmentation from a range image of blocks; the image size is 160 × 120. We have thirteen sets of segmented points in different colors except for those colored in light green that are detected as edge points. Among these thirteen sets, we tried to fit a digital plane to seven sets colored in green, brown, turquoise and violet. Table II shows the plane parameters for every set of points. The "Number of outliers" in Table II indicates the number of points that do not fit in the plane and N indicates the total number of points, whether they fit or not. Figures 10,11

V. COMPUTATION TIME ANALYSIS

The formulation by a mixed integer linear programming proves the efficiency of separating outliers from inliers while minimizing the number of outliers. However this formulation suffers from an important computation time when the number of points increases. In order to fit a digital line or plane to a given set S of N points, we need to solve all the n linear integer programming problems that consist of minimizing i=1,...,N p i . In terms of the number of variables, each of the [LP n -j] problems has N variables (p i for i = 1, 2, . . . N ) and n + 1 variables (α 1 , α 2 , . . . , α n , δ ). This increases the computation time.

In section IV, for the experiments, we used a set of points (x, y) ∈ Z 2 that satisfy y = -6x and another set of points (x, y) ∈ Z 2 that satisfies y = -x + 5. It was noticed that when the number of points N increases, the time complexity increases. Figure 14 shows the time variation when changing the number of points and fixing the number of outliers to 10. It is possible to see that the computation time increases as the number of points increases.

The computation time also depends on the number of outliers. Figure 15 shows the computation time variation of [LP 1-1] and [LP 1-2] when changing the number of outliers while fixing the number of points to 200. It should be noted that the computation time depends on the percentage of outliers with respect to the number of points; when the number of outliers is around 50 percent of the intial number of points, the computation time increases sharply. Figure 15 shows the time when the number of outliers is less than 25 percent or more than 75 percent of the number of points; when it is between these two numbers, the time increases rapidly to several hours (time not shown in the diagram since it exceeds the y-axis range). Moreover, it should be noted that when we increase the dimension from two to three, the time complexity also increases, since the number of formulated [LP] problems increases and the number of variables increases. Sometimes many hours are needed to solve a [LP] problem which makes Algorithm 1 impractical.

VI. IMPROVEMENT BY ADDING A CONSTRAINT

A. Problem

In most cases, each of the [LP] problems may give different solutions and one of the solutions is optimal since it gives the minimum i=1,...,N p i . In other cases, more than one system may We cannot know in advance which system is most suitable to the set S of points. This fact forces us to solve two [LP] systems in 2D and three in 3D and then to find among them the one that gives the minimum i=1,...,N p i . The computation time is a factor that depends not only on the number of points, but it also depends on the percentage of outliers with respect to the number of points (Section V). This increases the computation time when solving the [LP] problem that does not give the optimal fitting solution. Thus, it is likely that the method can be improved if we can find the right [LP] problem among the 2 (2D) or 3 (3D) [LP] problems.

B. Adding a constraint

We present a solution that decreases the computation time since it permits to detect rapidly the right equation that minimizes the objective function p i . It consists of adding the following constraint to each [LP] problem proposed in Section III-A: i=1,...,N

p i = K ( 18
)
where K is a constant that starts from 0, meaning no outliers, and then increases until a solution is found. When a solution is found there is no need to continue such an increment since we are sure that the solution is optimal. This solution gives us the faster computation time since each time the number of outliers is fixed, by fixing the range of searching solutions for each of the n [LP] problems (nD case). This is a good approach when there is a limited number of outliers, which implies that at least one of the [LP] problems has a limited number of outliers.

In such case even if the rest of the [LP] problems features a large number of outliers, they will have no solution with the constraint added. However when the number of outliers is relatively large, the computation time may increase and this approach is inefficient. Another approach is thus needed to solve the mixed integer programming problem. Figure 16 shows the variation of time as the constant K increases while fixing the number of points. It is possible to see that when K increases, the computation time also increases.

VII. CONCLUSION AND PERSPECTIVES

In this paper, we present a fitting method of a given set of grid points to a digital line or plane by formulating the problem into a mixed integer programming. The solution satisfies all points except for some outliers. The experimental results show that our method is useful for the estimation of the digital line and plane from a point cloud by determining the inliers and rejecting the outliers. It should be noted that there is no digitization error in the framework of discrete geometry so that outliers detected are not considered to be digitization errors. However, the algorithm proposed to solve the problem is time consuming when the numbers of points increased. Several solutions are expected to decrease the computation time. These solutions should be studied and tested to see their influences on the computation time as well as on the precision of the fitting; they can be divided into three directions:

• Reducing the number of points or outliers: reducing the number of points in the point cloud can be done by applying a multi-scale approach that increases the size of image pixel or volume and reduces the number of points; this decreases the number of points and may allow us to treat outliers more quickly. • Finding some outliers in a pre-processing step: finding outliers in a pre-processing step may reduce the computation time since it decreases the parameters space.

It is accomplished by using, for example, the Hough transform. The Hough transform will enable us to reject outliers before applying our method formulated into a mixed integer programming. • Proposing a new formulation: there already exist several well-known algorithms for fitting such as least median of squares linear regression [START_REF] Press | Numerical recipes (Third Edition)[END_REF]. However, the task of formulating the problem of finding digital lines or planes using these methods remains to be done. This formulation is currently under study trying to find a better and less time consuming algorithm.

Fig. 1 .

 1 Fig. 1. Fitting a digital line to a given set of points rejecting some points considered as outliers

Algorithm 1 : 6 F = F j ; 7 set α 1

 1671 Fitting a digital line or plane to a given sets of grid points input : A given set S of N points in (n + 1)D for n = 1 or 2 output: The parameters α 1 , . . . , α n , δ of the nD hyperplane and p i of every point in S begin 1 initialize F = N ; 2 for j = 1, ..., n do 3 solve [LP n-j] ;4set F j to be the optimal value of the objective 5 function; if F j < F then , . . . , α n , δ and p i to be the optimal

Fig. 2 .

 2 Fig. 2. Two solutions are obtained by solving [LP 1-1] and [LP 1-2] when fitting a digital line to the given set of points. The first one rejects 14 outliers while the second one rejects 6 outliers. It is thus clear that the second one gives the optimal result since it rejects less points.

Fig. 3 .

 3 Fig. 3. An original image

Fig. 5 .Fig. 6 .

 56 Fig. 5. Digital line fitting (red line)

  , 12 and 13 show the digital planes obtained after the fitting for the

Fig. 7 .Fig. 8 .

 78 Fig. 7. Line fitting (colored red) using the Hough transform with ∆θ = 1 and t=100

Fig. 9 . 9 Fig. 12 .

 9912 Fig. 9. Planar surface segmentation from a range image of blocks

Fig. 14 .

 14 Fig. 14. Time variation as the number of points in the point cloud increases. Fig. 15. Time variation as the number of outliers increases.

Fig. 16 .

 16 Fig. 16. Time variation with respect to the constant K

ACKNOWLEDGMENT

This work was supported in part by ANR grant SURF 05-BLAN-0071.


Rita Zrour
email: zrourr@esiee.fr

Yukiko Kenmochi
email: y.kenmochi@esiee.fr

Hugues Talbot
email: h.talbot@esiee.fr

Lilian Buzer
email: l.buzer@esiee.fr

Yskandar Hamam
email: y.hamam@esiee.fr

Ikuko Shimizu

Akihiro Sugimoto
email: sugimoto@nii.ac.jp

Optimal consensus set for digital plane fitting

This paper presents a method for fitting a digital plane to a given set of points in a 3D image in the presence of outliers. We present a new method that uses a digital plane model rather than the conventional continuous model. We show that such a digital model allows us to efficiently examine all possible consensus sets and to guarantee the solution optimality and exactness. Our algorithm has a time complexity O(N 3 log N ) together with a space complexity O(N ) where N is the number of points.

Introduction

Plane fitting is one of essential tasks in the field of 3D computer vision. For instance, this procedure is useful in image segmentation [START_REF] Köster | Mir: An approach to robust clustering -application to range image segmentation[END_REF][START_REF] Kenmochi | Discrete plane segmentation and estimation from a point cloud using local geometric patterns[END_REF], multi-view image registration [START_REF] Shum | Principal component analysis with missing data and its application to polyhedral object modeling[END_REF][START_REF] Chen | A fast automatic method for registration of partially-overlapping range images[END_REF], as well as parameter fitting [START_REF] Hartley | Multiple view geometry in computer vision[END_REF].

Several significant, optimal methods have been proposed in the literature, such as least-squares fitting, least-absolutevalues fitting and least median of squares (LMS) [START_REF] Boyd | Convex optimization[END_REF][START_REF] Hartley | Multiple view geometry in computer vision[END_REF][START_REF] Press | Numerical recipes: the art of scientific computing[END_REF]. In these methods, we use a continuous plane model that is defined by

P = {(x, y, z) ∈ R 2 : ax + by + c + z = 0}, (1)
where a, b, c ∈ R. The fitting is then carried out by optimizing different cost functions. For instance, least-squares minimizes the sum of geometric distances from all given points to the model. The solution can be obtained analytically, however it is very sensitive to the presence of outliers. Least-absolute-values fitting uses the L 1 norm, instead of the geometric distance, for its minimization. Some efficient iterative algorithms have been proposed for the fitting. However, they are still sensitive to outliers. In contrast, LMS minimizes the median of vertical/geometric distances of all given points to the model. Thus, the fitting is robust as long as the majority of given points are not outliers [START_REF] Rousseeuw | Least median of squares regression[END_REF].

In this paper, we present a globally optimal method that, for a given arbitrary cloud of points, locates the planes minimizing the number of outliers, or equivalently maximizing the number of inliers, also called the consensus set. The idea of using such consensus sets were proposed in the RANdom SAmple Consensus (RANSAC) method [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF], which is widely used in the field of computer vision. However RANSAC (and its variations) is inherently probabilistic in its approach, and does not deterministically guarantee the optimality of consensus sets. In contrast, our method is both deterministic and optimal.

In order to guarantee the optimality of consensus sets, we follow a digital geometry methodology [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF] instead of the continuous model of [START_REF] Boyd | Convex optimization[END_REF]. This methodology is arguably preferable if given all inputs are digital images. For instance, a digital model allows us to distinguish digitizationinduced noise from actual noise produced by an image sensor. Related work using digital plane models can be found in [START_REF] Kim | On the recognition of digital planes in three-dimensional space[END_REF][START_REF] Provot | Recognition of blurred pieces of discrete planes[END_REF] for digital plane recognition, and shape polyhedrization [START_REF] Sivignon | Decomposition of a three-dimensional discrete object surface into discrete plane pieces[END_REF] in both of the absence and presence of noise. However, to the best of our knowledge, outliers, namely points that do not fit the model, have never been dealt with in the field of digital geometry.

This work is an extension of our previous work on digital line fitting [START_REF] Zrour | Optimal consensus set for digital line fitting[END_REF], but uses a significantly different method. Here, we show that, even though our digital plane fitting problem is 3D, we can treat it as a 2D problem by considering it in the dual space of the duality transform [START_REF] Berg | Overmars[END_REF]. We show that our algorithm for digital plane fitting has a time complexity of O(N 3 log N ), instead of O(N 2 log N ) for digital line fitting, where N is the number of points, while the space complexity O(N ) is not affected by the dimension increase. We also point out that there are more degenerated cases for 3D than those for 2D, if our input data are all inte-gers or rational numbers, which is often the case for digital images. We thus present a way to handle such degenerated cases with only simple modifications.

The rest of the paper is as follows: in section 2 we expose the framework of our digital model. In section 3 we prove the optimality of our result. In section 4 we provide an algorithm for the computation of the fit. Section 5 provides a method for extracting the parameters from the fit. Section 6 is devoted to results and applications. Section 7 states some conclusions and perspectives.

The problem of digital plane fitting

In this paper, we use a digital model, instead of the continuous one, for planes in a discrete space Z 3 where Z is the set of all integers. We consider that using a digital model is natural when our input data is a digital image. In addition, we guarantee the optimal solution for our problem, as we mention in the next section.

A digital plane D(P) that is the digitization of P in ( 1) is defined by the set of discrete points satisfying two inequalities:

D(P) = {(x, y, z) ∈ Z 3 : 0 ≤ ax+by +z +c ≤ w} (2)
where w is a given constant value. Geometrically, D(P) is a set of discrete points lying between two parallel planes ax + by + z + c = 0 and ax + by + z + c = w, and w specifies the vertical distance between them. From the digital geometrical point of view [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF], w should not be less than 1 if we expect that D(P) is 18-connected. More precisely, in order to keep the connectivity with the minimum distance, we need to consider the distance w not only in the z-axis direction as shown in (2) but also in the xand yaxis directions. For digital planes D(P) with the distance w in the other directions, we simply make the permutation of x, y and z in (2), and put the constraints, -1 ≤ a ≤ 1, -1 ≤ b ≤ 1. Hereafter, discussions will be made by using [START_REF] Chen | A fast automatic method for registration of partially-overlapping range images[END_REF], since they are also valid for the other models.

Using the above digital plane model, our fitting problem is then described as follows: given a finite set of discrete points,

S = {x i ∈ Z 3 : i = 1, 2, . . . , N },
we would like to find a digital plane D(P) such that D(P) contains the maximum number of points in S. Points x ∈ S are called inliers if x ∈ S∩D(P); otherwise, they are called outliers.

Digital planes and their consensus sets

In this section, we show an important property of inlier sets, also called consensus sets, of digital planes. Our algorithm presented in the next section is based on this property.

Since the size of S is finite and each element x ∈ S has finite coordinates, we notice that the number of different consensus sets for a digital plane fitting of S is finite as well. Thus, if we can find all different consensus sets C from a given S, we just verify the size of each C to find the maximum one (ones if there are several) as the optimal solution. In the following, we will show that it is possible.

Before presenting the main proposition, we give some notions related to digital planes. The two parallel planes that are given by the equations (2) are called the support planes of a digital plane. Discrete points that are on support planes are called critical points of a digital plane.

Proposition 1. Let C be a consensus set of S for a digital plane. It is possible to find a new digital plane whose consensus set is the same as C such that it has at least three critical points.

Proof. Let D be an initial digital plane that contains all points in C as its inliers. Then, the following four cases can be considered when observing the critical points of D.

(1) Suppose that D has more than two critical points, then the proposition is correct in this case.

(2) Suppose that D has two critical points p 1 and p 2 . p 1 and p 2 may be located on one side or either side of the two parallel support planes of D as illustrated in Figures 1 (a) or 2 (a), respectively. First, we take the projections p 1 (resp. p 2 ) of p 1 (resp. p 2 ) on the other support plane where p 1 (resp. p 2 ) does not exist in the z-axis direction. We then apply a rotation to D in such a way as to maintain the distance w between the support planes until finding another point p 3 in C so that p 3 becomes a critical point. To achieve this, in the illustrative case of Figure 1 (a), to the support plane where p 1 and p 2 exist, we apply a rotation around the line going through p 1 and p 2 , as illustrated in Figures 1 (b). To the other support plane, we apply a rotation around the line going through p 1 and p 2 , as illustrated in Figures 1 (b). In the case of Figure 2 (a), similarly, to the support plane where p 1 (resp. p 2 ) exists, we apply a rotation around the line going through p 1 and p 2 (resp. p 1 and p 2 ), as illustrated in Figures 2 (b). Note that we can rotate D either clockwise or counterclockwise.

(3) Suppose that D has one critical point p 1 as illustrated in Figure 3 (a). In this case, we also consider the projection of p 1 , p 1 . We then apply a rotation to each support plane until finding another point p 2 in C so that p 2 becomes a critical point, as illustrated in Figure 3 (b). The support plane where p 1 (resp. p 1 ) exists is rotated around any line going through p 1 (resp. p 1 ) on the support plane. If just one point p 2 is found as a second critical point after the rotation, then another rotation is made, as mentioned in the previous case, in order to obtain a third critical point p 3 . (4) Suppose that D has no critical point as illustrated in Figure 4 (a). In this case, we first apply a translation to D in order to find a first critical point p 1 . Note that a translation can be made in any direction while the two support planes maintain a constant distance w between them. During such a translation, if more than two points are found as critical points, then the proof is complete. If just one point p 1 is found, as illustrated in Figure 4 (b), then a rotation is made as mentioned in the previous case, in order to obtain a second critical point p 2 . If during this rotation just one point p 2 is detected, then another rotation is accomplished in order to detect a third critical point p 3 as described in the first case. From this proposition, we see that we can find a digital plane D(P) for any consensus set C of S such that it has at least three critical points. This is intuitively understandable, because when we move a digital plane D(P) in the space, its consensus set C will be changed the moment that a critical point goes out from D(P), namely, becomes an outlier, due to the motion. Indeed, such a digital plane D(P) can be constructed from a triplet of points chosen from S such that they become critical points of D(P). Consequently, we can find all C from those D(P) constructed from all possible triplets of points in S. Note that there are at most eight different D(P) constructed from each triplet of points in S.

Digital plane fitting algorithm

A naive algorithm for 3D digital plane fitting derived from Proposition 1 is to make the consensus set C for each triplet of points in S and to verify the size of C for obtaining the maximum one among all triplets. Thus, such an algorithm has a O(N 4 ) time complexity. In this section, we will present a new algorithm that has a less complexity.

Our algorithm is based on a similar idea to the one for 2D digital line fitting, presented in [START_REF] Zrour | Optimal consensus set for digital line fitting[END_REF]. The key idea for the extension to 3D digital plane fitting is treating the 3D problem as a 2D problem. We first show how to reduce the dimension from three to two, and then obtain an algorithm providing a O(N 3 log N ) time and O(N ) space complexity.

We first describe the digital plane fitting problem in the dual space of the duality transform, because our algorithm works in the dual space of the duality transform [START_REF] Berg | Overmars[END_REF], similarly to the 2D digital line fitting. We then present an algorithm to exhibit the optimal consensus set (sets if there are several) that maximizes the number of inliers of a fitted digital plane from a given set S of discrete points in 3D, step by step. We also describe special treatments for degenerated cases; it should be noted that digital images likely present many degenerated cases that must be processed separately.

Digital plane fitting in the dual space

A point p = (x, y, z) in the primal space associates to a non-vertical plane

P p = {(a, b, c) : xa + yb + c + z = 0} (3)
in the dual space. Conversely, a non-vertical plane in the primal space associates to a point in the dual space. Because a digital plane defined by ( 2) is regarded as a set of nonvertical parallel planes whose normal vectors are (a, b, 1) and whose z-intercepts are between -c and w -c, it forms a vertical line segment of length w in the dual space as illustrated in Figure 5. Since all points of S in the primal space are represented by planes in the dual space, our problem of finding the optimal consensus set for digital plane fitting in the primal space is equivalent to searching the best position of the vertical line segment of length w such that it intersects with the maximum number of planes in the dual space.

We now need a search procedure for an optimal segment. Thanks to Proposition 1, we know that, for any consensus set, there exists a digital plane featuring at least three critical points, among which at least two are on one of the support planes. Thus, taking two different points p, q from S in the primal space, we first consider all the digital planes on which both p and q are critical points on the same support plane. In the dual space, digital planes having two critical points p, q forms two strips, which will be described in Section 4.2. In Section 4.3, we explain how digital planes with two critical points p, q appear in the strips when they have a third critical point r. In Section 4.4, we then show how to obtain the maximum number of inliers and its digital plane parameters in the strips.

Strips made from a critical point pair

Let p = (x p , y p , z p ) and q = (x q , y q , z q ). In the dual space, they represent two planes P p and P q , defined by [START_REF] Berg | Overmars[END_REF]. They intersect in a line L 0 pq if p and q are chosen such that (x p -x q ) 2 + (y p -y q ) 2 = 0; otherwise, P p and P q are parallel, and no intersection line can be found. The intersection line L 0 pq is represented by the following equation:

L 0 pq = {v = (a, b, c) : v = u + td, t ∈ R}
where d = (x p , y p , 1) ∧ (x q , y q , 1) = (y p -y q , x q -x p , x p y q -x q y p ), and u = (u a , u b , u c ); u is a chosen point on L 0 pq . For example, if x p y q = x q y p , by fixing u c = 0, u a and u b are automatically found since u is on both P p and P q .

Once L 0 pq is found, then, all the digital planes on which both p and q are critical points on the same support plane in the primal space correspond to the set of all the vertical line segments of length w having one of its endpoints on L 0 pq in the dual space, as shown in Figure 6. We see in the figure that the set of such digital planes, therefore, forms two strips in the plane Q pq that contains L 0 pq and the direction parallel to the c-axis. Taking the d-axis in Q pq as the orthogonal one to the c-axis, such Q pq is illustrated in Figures 6 and7. Each strip on Q pq illustrated in Figure 7 is bounded by two parallel lines, L 0 pq and L i pq for i = 1, 2, which are represented by:

L 1 pq = {v = (a, b, c) : v = u + e + td, t ∈ R}, L 2 pq = {v = (a, b, c) : v = u -e + td, t ∈ R},
where e = (0, 0, 1).

Digital planes with a critical point triplet

Hereafter, we focus on one of the strips in Q pq , because the following discussion is valid for both strips. Let us consider the strip bounded by L 0 pq and L 1 pq , as illustrated in Figure 7. According to Proposition 1, we choose a point r ∈ S \ {p, q} to be the third critical point of a fitted digital plane such that r is not colinear with p and q; the colinear case will be handled separately as a degenerated case in Section 4.6. Any point r in the primal space is represented by the line L r in Q pq in the dual space, which is the intersection between P r and Q pq , as shown in Figure 7. We see in this figure that L r intersects each of the strip boundaries, L 0 pq and L 1 pq , if it is not parallel to L 0 pq ; the parallel case will be also dealt with separately as a degenerated case in Section 4.6. The intersections between L r and L i pq , σ i r = (a i r , b i r , c i r ), for i = 0, 1, are calculated from L i pq and P r . Geometrically, the vertical line segment in the strip, one of whose endpoints is one of the intersections σ i r , in the dual space corresponds to a digital plane with three . Four points p, q, r and s in the primal space (left), and their interpretations in the cross-section Qpq of the dual space (right). Qpq is made as the plane that contains the intersection line L 0 pq of Pp and Pq and the parallel direction to the c-axis, as illustrated in Figure 6. In Qpq, all the digital planes having p and q as critical points are represented by the strips each of which is bounded by L 0 pq and either of its parallel lines L 1 pq and L 2 pq . The other points r and s in the primal space are represented by the two lines Lr and Ls in Qpq. critical points p, q and r in the primal space. This indicates that the digital planes corresponding to the vertical line segments between the two intersections σ 0 r and σ 1 r in the strip always contain r as an inlier.

Finding the largest consensus set in a strip

In order to know the number of inliers within the digital planes with two critical points p and q, we check the intersections σ 0 r and σ 1 r of L r for all r ∈ S\{p, q} with the strip boundaries, L 0 pq and L 1 pq . We use a function f i r , for i = 0, 1, which is set to be 1 if L r enters the strip from L i pq , and -1 if L r leaves the strip from L i pq . Once the intersections σ i r = (a i r , b i r , c i r ), and the associated function f i r for i = 0, 1 are calculated for all r ∈ S \ {p, q}, we sort all the quadruples (a i r , b i r , c i r , f i r ) in increasing order by using either a i r or b i r as keys; if Q pq is not perpendicular to the a-axis, we use a i r ; otherwise, we use b i r . As for determining the location of the maximum number of inliers, another function F (a) (resp. F (b), depending on the key selection) is used; after initially setting F (a) = 2 for every a, since we already know that p and q are inliers, then the value f i r is added to F (a) in the above sorted order. By looking for the maximum value of the function F (a), we obtain the parameter set (a, b, c) corresponding to the maximum optimal consensus set for a pair of critical points p and q. In this section, we consider that all L r enter or leave a strip at different moments. The degenerated cases such that many lines L r enter or leave a strip at the same moment will be described in Section 4.6.

Algorithm

We now present Algorithm 1. Input is a set S of discrete points and a distance value w of our digital plane model. Output is a set V of parameter values (a c , b c , c c ) corresponding to the fitted digital planes of that give the optimal consensus sets. In the algorithm, we consider another strip bounded by L 0 pq and L 2 pq in Steps 5, 11 and 25. We remark that, because c i r is not used for the sorting step and can be calculated from a i r and b i r , we do not have to store it for each intersection. Simply for a candidate of the optimal consensus set, we calculate it as shown in Steps 26 and 28. Note that, depending on the strip, we calculate different c c because of the translation difference w between the two strips.

In Steps 13 and 19, we only show the case where a k is used as keys for sorting. However, if Q pq is perpendicular to the a-axis, all a k has the same value. In such a case, as mentioned above, we use b k as keys, instead of a k .

The time complexity of the algorithm is O(N 3 log N ), because we have N points in S and each pair of p and q in S needs the complexity O(N log N ) for sorting at most 2N -4 different values a i r for r ∈ S \ {p, q} and i = 0, 1.

The space complexity is O(N ) because for each sorting we have at most

2N -4 different triplets (a k , b k , f k ).
Because all inputs can be given as integers or rational numbers, all computations in Algorithm 1 can be performed using only rational numbers. This ensures that all results obtained by Algorithm 1 contain no calculation error. However, degenerated cases may occur, which will be discussed

Degenerated cases

In this section, we explain how we can deal with degenerated cases, which are not considered in Algorithm 1. The following three cases are summarized.

First, if three points p, q and r are colinear in the primal space, their associated planes P p , P q and P r have a line intersection in the dual space. Therefore, for any digital plane having p and q as its critical points also has r as its another critical point. Thus, the function F (a) initially set to 2 for the inclusion of p and q as inliers will be automatically increased by 1 because of the inclusion of r.

Secondly, suppose that p, q and r are not colinear, but there is no intersection between L 0 pq and L r in Q pq ; L r is parallel to L 0 pq . If L r is between L 0 pq and L 1 pq (resp. L 2 pq ), then we set the initial value of the function F (a) to 3 when l = 1 (resp. l = 2) because r is an inlier for any a. Otherwise, we set it to 2, as described in Algorithm 1, because r is an outlier for any a.

Lastly, when many lines L r enter or leave a strip at the same moment a, all the positive valued f i r of that moment must be added to the function F (a) at once (Step 21 in Algorithm 1), and the value F (a) is compared with the current maximum value M ax (Step 22 in Algorithm 1). Note that all the negative valued f i r of the same moment a must be added after the comparison to the value F (a). Indeed such a point r must be considered as an inlier until that moment.

Feasible digital plane parameters

Once an optimal consensus set C for digital plane fitting to a given point set S is found, we need the corresponding digital plane parameters in most applications. A continuous plane model such as (1) can be used to estimate them, for example, by applying the least squared method [START_REF] Hartley | Multiple view geometry in computer vision[END_REF] to C. However, one must be careful because this may modify the inlier set. In such a case, a new C must be recalculated from a new estimated plane, leading to an inefficient iterative procedure.

In our case, however, since we use the digital plane model of (2) instead of (1), we do not need such an estimation method, and there is no danger or changing the inlier set. We can obtain all feasible solutions for the parameters of digital planes fitted to an obtained optimal C by simply looking for all feasible solutions (a, b, c) that satisfy the inequalities of (2) for all (x, y, z) ∈ C.

Such feasible solutions of digital planes are called preimages [START_REF] Garard | Gift-wrapping based preimage computation algorithm[END_REF]. In contrast to the preimages of digital lines in 2D, the structure of the preimages of digital planes is not simple; we even do not know the maximum number of vertices or facets of a convex polyhedron, that constitutes a preimage of a digital plane. Further study is necessary.

Experiments

For our experiments, we applied our proposed method to two example data, such as a 3D discrete point cloud and a 3D binary digital image. The first example is a 3D discrete point cloud in Figure 8, which is obtained after a planar surface segmentation of a range image of blocks [START_REF] Kenmochi | Discrete plane segmentation and estimation from a point cloud using local geometric patterns[END_REF]. The number of points in the cloud is 12859, and they are segmented into thirteen planar surfaces, which are illustrated in Figure 8 with points in different colors, except for those colored in light green that are detected as edge points. For each of these thirteen sets, we fitted a digital plane. We see the corresponding planes in Figure 9, and the number of points for each segmented surface and the size of its optimal consensus set in Table 1. In Figure 10, we also see that the fitted plane for the blue segmented surface points in Figure 8: inliers are colored blue while outliers are in pink.

We also applied our method to a 3D image extracted from a polymer foam observed in X-ray micro-tomography, on which homotopic thinning and surface decomposition were applied [START_REF] Plougonven | Quantitative analysis of the deformation of polypropylene foam under dynamic loading[END_REF]. Figure 11 shows a cross section of the original image and Figure 12 shows a 3D binary image obtained after homotopic thinning and surface decomposition; the image is cut into two parts for visualization. Among around 400 sets of points forming surfaces in the entire image, we choose a part, as illustrated in Figure 13, including 17 decomposed surfaces for digital plane fitting. We show the fitted planes in Figure 14, and the number of points and the optimal consensus set size for each segmented surface in Table 2. For both the examples, we set w = 1.

Conclusions

In this paper we have exposed a new method for plane fitting on discrete data such as bitmap images and volumes, using a digital geometry (DG) approach. The DG approach allows practitioners to separate effects due to digitization on the one hand and noise on the other. Using our approach, we have proposed an optimal fitting method from the point of view of the maximal consensus set: we are guaranteed to fit discrete planes with the least amount of outliers.

Our algorithm has a complexity that are identical to parameter-less traditional plane-fitting algorithms such as least median of squares regression [START_REF] Rousseeuw | Least median of squares regression[END_REF], but improves the detection of digital planes, in the presence of outliers. Future work will include improving algorithmic complexities and more complete applications such as optimal polyhedrization, image registration considering all feasible digital plane parameters. Table 2. The number of points and the optimal consensus set size for each decomposed surface in Figure 14.

Figure 1 .Figure 2 .

 12 Figure 1. A digital plane that has two critical points p 1 and p 2 on one of its support planes (a), and its rotated digital plane that also has a third critical point p 3 (b).

Figure 4 .

 4 Figure 4. A digital plane that has no critical point (a), and its translated digital plane that has one critical point p 1 (b).

Figure 5 .

 5 Figure 5. A digital plane in the primal space (left) corresponds to a vertical line segment of length w in the dual space (right).

Figure 6 .

 6 Figure 6. All the digital planes with two critical points p and q in the primal space (left) correspond to a set of vertical line segments of length w having one of its endpoints on the intersection line of the two planes Pp and Pq (right).

Figure 7

 7 Figure 7. Four points p, q, r and s in the primal space (left), and their interpretations in the cross-section Qpq of the dual space (right). Qpq is made as the plane that contains the intersection line L 0 pq of Pp and Pq and the parallel direction to the c-axis, as illustrated in Figure6. In Qpq, all the digital planes having p and q as critical points are represented by the strips each of which is bounded by L 0 pq and either of its parallel lines L 1 pq and L 2 pq . The other points r and s in the primal space are represented by the two lines Lr and Ls in Qpq.

Algorithm 1 : 4 for l = 1 , 2 do 5 initialize F = 2 ; 6 initialize 8 foreach r ∈ S \ {p, q} do 9 calculate a i r and b i r for i = 0, 1 ; 10 if l = 2 then 11 calculate a 2 r and b 2 r and reset 12 a 1 r = a 2 r and b 1 r = b 2 16 setF

 141526891101111216 Digital plane fitting input : A set S of N discrete points, a distance w output: A set V of parameter value triplets (a c , b c , c c ) of the best fitted digital planes begin 1 initialize M ax = 0; 2 foreach {p, q} ∈ S do 3 calculate L 0 pq ; the array T [k] for 7 k = 1, . . . , 2N -4; set j = 0; the tuple, (a i r , b i r , f i r ), for i = 0, 1, 17 in T [2j + i]; j = j + 1; 18 sort all the tuple elements (a k , b k , f k ) for 19 k = 1, . . . , 2j in T with a k as keys; for k = 1, . . . , 2j do 20 = F + f k ; 21 if F > M ax then set M ax = F , 22 V = ∅; if F = M ax then 23 set a c = a k , b c = b k ; 24 if l = 1 then 25 c c = -a k x p -b k y p -z p ; 26 else 27 c c = -a k x p -b k y p -z p + w; 28 put (a c , b c , c c ) in V;

Figure 8 .

 8 Figure 8. Planar surface segmentation of a 3D discrete point cloud: the number of points is 12859, and they are segmented into thirteen planar surfaces whose points are in different colors, except for those colored in light green that are detected as edge points.

Figure 9 .

 9 Figure 9. Fitted planes of segmented planar surface in Figure 8.

Figure 10 .

 10 Figure 10. The fitted plane with its optimal consensus set for the blue segmented surface points in Figure8: inliers are colored blue while outliers are colored pink.

Figure 11 .

 11 Figure 11. A cross section of a 3D image extracted from a polymer foam observed in X-ray micro-tomography.

Figure 12 .

 12 Figure 12. The 3D binary image obtained after homotopic thinning and surface decomposition applied on the image in Figure 11: the image is cut into two parts for visualization.

Figure 13 .

 13 Figure 13. Selected decomposed surfaces, which is a part of the 3D binary image in Figure 12, for digital plane fitting.

Figure 14 .

 14 Figure 14. Fitted digital planes for decomposed surfaces shown in Figure 13.

Table 1 .

 1 The number of points for each segmented surface in Figure8and the size of its optimal consensus set.

		Number of points Opt. consensus set size
	Blue	1770	1401
	Yellow	1578	1195
	Pink	1523	935
	Pale Blue	1191	922
	Orange	699	693
	Green	573	573
	Brown	545	544
	Turquoise	536	512
	Olive	440	405
	Purple	248	245
	Violet	232	206
	Moss green	223	223
	Cream	101	97

Acknowledgement

The authors express their thanks to Mr. John Chaussard, who supplied the 3D image of decomposed surfaces extracted from a polymer foam used in the experiments. A part of this work has been supported by the French National Research Agency projects, ANR-05-BLAN-0071 (SURF) and BLAN07-2-184378 (MICROFISS).