Jérôme Revaud
email: jerome.revaud@insa-lyon.fr

Guillaume Lavoué
email: glavoue@insa-lyon.fr

Atilla Baskurt
email: abaskurt@insa-lyon.fr

Improving Zernike Moments Comparison for Optimal Similarity and Rotation Angle Retrieval

Keywords: Zernike moments, scene analysis, 3D object recognition, shape !

Zernike moments constitute a powerful shape descriptor in terms of robustness and description capability. However the classical way of comparing two Zernike descriptors only takes into account the magnitude of the moments and loses the phase information. The novelty of our approach is to take advantage of the phase information in the comparison process while still preserving the invariance to rotation. This new Zernike comparator provides a more accurate similarity measure together with the optimal rotation angle between the patterns, while keeping the same complexity as the classical approach. This angle information is particularly of interest for many applications, including 3D scene understanding through images. Experiments demonstrate that our comparator outperforms the classical one in terms of similarity measure. In particular the robustness of the retrieval against noise and geometric deformation is greatly improved. Moreover, the rotation angle estimation is also more accurate than state of the art algorithms.

INTRODUCTION

Zernike moments are widely used to capture global features of an image in pattern recognition and image analysis. Firstly introduced in computer vision by Teague [START_REF] Teague | Image analysis via the general theory of moments[END_REF], this shape descriptor has proved its superiority over other moment functions [START_REF] Belkasim | Pattern recognition with moment invariants: A comparative study and new results[END_REF], [START_REF] Liao | On image-analysis by moments[END_REF] regarding to its description capability and robustness to noise or deformations. Hence rotation invariant pattern recognition using Zernike moments has been extensively studied [START_REF] Khotanzad | Invariant image recognition by zernike moments[END_REF], [START_REF] Kan | Invariant character recognition with zernike and orthogonal fourier-mellin moments[END_REF]. Even very recently, a lot of authors have been working on these moments, particularly to improve their computation time [START_REF] Chong | A comparative analysis of algorithms for fast computation of zernike moments[END_REF], [START_REF] Wee | Efficient computation of radial moment functions using symmetrical property[END_REF], [START_REF]On the computational aspects of zernike moments[END_REF], [START_REF] Kotoulas | Real-time computation of zernike moments[END_REF] or their accuracy [START_REF] Amayeh | Accurate and efficient computation of high order zernike moments[END_REF].

Practically one Zernike moment is a complex number that contains two different values: magnitude and phase, however, the usual way (i.e. used in all existing algorithms) of comparing two Zernike descriptors only considers the moments magnitudes (as it brings invariance to rotation). In the context of 2D and 2D-3D indexing and recognition, this loss of information is not harmless when comparing two different patterns, and can induce erroneous results and impreciseness, as it will be further illustrated.

Using the phase information of Zernike moments (together with the magnitude) in the comparison process seems a natural way to improve the similarity measure in terms of robustness against geometric deformation or noise particularly. However in that case the resulting comparator is not invariant anymore to rotation, unless the in-plane rotation angle between the two patterns is known. Fortunately in this paper we show that the moment phases can also be used to retrieve this rotation angle in an optimal way. Finding both information (i.e. a robust rotation-invariant similarity measure together with the optimal angle of rotation) can be of great interest for many applications including image registration [START_REF] Reddy | An fft-based technique for translation, rotation, and scale-invariant image registration[END_REF], motion estimation in video and particularly scene understanding: indeed recognizing the objects composing the image and then extracting their in-plane orientation angles may help to compute their precise 3D pose and thus to understand accurately the corresponding 3D environment. A lot of work has been done for angle/similarity recognition using keypoint-based local descriptors like SIFT [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF], however, this kind of tools works only on textured objects and fails to describe smooth shapes or drawings (i.e. sketch) for instance. In such hard description/recognition cases, global shape descriptors like Zernike moments are particularly robust, that is the reason why they have recently been used for rotation invariant 2D/3D object recognition through sketches [START_REF] Filali Ansary | A bayesian 3D search engine using adaptive views clustering[END_REF], [START_REF] Hou | Calligraphic interfaces: Classifier combination for sketch-based 3d part retrieval[END_REF]. It appears quite relevant to compute the in-plane rotation together with the similarity distance in such 2D/3D indexing scenarios, particularly for some emerging applications like sketch-based modeling [START_REF] Kara | Sketch-based 3d-shape creation for industrial styling design[END_REF]. Apart from us, one approach have focused on the sub-problem of the rotation angle estimation using Zernike moments. The method was brought by Kim and Kim [START_REF] Kim | Robust rotation angle estimator[END_REF] and it proved to be very robust with respect to noise even for circular symmetric patterns. Nevertheless, the probabilistic model used to recover the rotation angle has no concrete interpretation and does not correspond to any geometrical reality, so it does not return any similarity distance. Besides, this method is based on the hypothesis that the two patterns are the same (that is, except some noise and the rotation), which does not always hold in practice.

In this context, we have developed a new general and rigorously founded approach for comparing two Zernike descriptors that takes use of both magnitude and phase information. Our approach keeps the same complexity as the standard technique (Euclidean distance between magnitude values) but provides a more accurate rotation invariant similarity measure while retrieving the optimal in-plane rotation angle. Thanks to the adaptability of Zernike description, our approach is suited to compare any kinds of images/patterns: Binary, gray level or sketch images (i.e. drawings). We also compared our results with two state-of-theart approaches for sketch and object recognition: the geometric hashing [START_REF] Wolfson | Geometric hashing: an overview[END_REF] and the deformation tolerant generalized Hough Transform from Anelli et al. [START_REF] Anelli | Deformation tolerant generalized hough transform for sketch-based image retrieval in complex scenes[END_REF].

The following section concisely presents Zernike moments. In section 3, we lean upon drawbacks of the conventional approach to build our method. The computational efficiency is also a constraint because the resulting algorithm will be used within a matching process; we thus detail a fast implementation in section 4. Finally, we present experimental results in section 5 and an application to a real 2D/3D indexing scenario in section 6.

ZERNIKE MOMENTS

Complex Zernike functions constitute a set of orthogonal basis functions mapped over the unit circle. Zernike moments of a pattern are constructed by projecting it onto those functions. They share three main properties:

• The orthogonality: this property ensures that the contribution of each moment is unique and independent.

• The rotation invariance: the magnitude of Zernike moments is independent of any planar rotation of the pattern around its center of mass. • The information compaction: low frequencies of a pattern are mostly coded into the low order moments. As a result, relatively small descriptors are robust to noise or deformations. Mathematically, Zernike basis functions are defined with an order p and a repetition q over D = {(p, q)|0 ≤ p ≤ ∞, |q| ≤ p, |p -q| = even}.

Z pq = p + 1 π x 2 +y 2 ≤1 V * pq (x, y)f (x, y)∂x∂y (1)
where V * pq denotes the complex conjugate of V pq , itself defined as:

V pq (ρ, θ) = R pq (ρ).e iqθ (2)
and

R pq (ρ) = p k = |q| |p -k| even (-1) p-k 2 p+k 2 ! p-k 2 ! k-q 2 ! k+q 2 ! ρ k
From eq. ( 1) and ( 2), Zernike moments of a pattern rotated by an angle α around its origin are given in polar coordinates as :

Z α pq = Z pq e iqα (3)
Eq (3) proves the invariance of the magnitude of Zernike moments to rotation since |Z pq e iqα | = |Z pq |.

Thanks to the property of orthogonality, the reconstruction of the pattern can be simply expressed as the sum of every Zernike basis functions weighted by the corresponding moments:

f (x, y) = (p,q)∈D Z pq V pq (x, y) (4)

SIMILARITY MEASURE AND ROTATION AN-

GLE RETRIEVAL

The classical approach

The usual way of comparing two patterns in the Zernike space only considers the magnitudes of the moments [START_REF] Khotanzad | Invariant image recognition by zernike moments[END_REF]. In the reminder of this paper, we will denote the usual comparator of Zernike descriptors as the classical one. Formally, this comparator is nothing else than an Euclidean distance between the magnitudes:

d 2 = (p,q)∈D (|Z pq | -|Z pq |) 2 (5)
However, the advantage of losing the phase information -this allows the invariance to rotation -also brings some drawbacks: The first consequence is that the classical comparator is unable to retrieve the rotation angle between two similar patterns as this information is encoded onto the moments phases. A corollary is that two symmetrical patterns will be classified as identical since their moments magnitudes are the same. Consequently, for an application which has to differentiate symmetric patterns (e.g. 2D-3D recognition that aims to retrieve viewpoint of potentially symmetric 3D objects like cars or people), the classical comparator is inoperative.

More generally, one can assume that this missing information has a negative influence on the retrieval efficiency, especially in hard cases like noisy, deformed or occluded patterns. Our approach is based on the assumption that using the phase information may result in a more robust description. Experiments have confirmed that this hypothesis is exact whatever the number of moments and for a wide type of distortions (see section 5). Figure 1 presents a short example of image retrieval from sketch which illustrates the superiority of the proposed approach upon the classical one.

The proposed Zernike comparator

In this subsection, we present a new way of comparing Zernike moments that takes both magnitude and phase information into account. Our new comparator provides a similarity score more robust than the classical one and retrieves for free an accurate angle of rotation between the two images. The angle is considered to be optimal since the Euclidean distance between the first image and the rotated second one is minimized. Let I and J be two different images, and (J * φ ) be the J image rotated by φ. The Euclidean distance between I and (J * φ ) can be expressed as a function of the rotation angle φ as follows:

d 2 I,J (φ) = x 2 +y 2 ≤1 |I(x, y) -(J * φ )(x, y)| 2 (6)
Thus our objective is to minimize this expression so as to determine the corresponding angle φ. Eq. ( 3) has shown that if the set of moments Z J pq |(p, q) ∈ D represents the J image, then Z J pq e iqφ |(p, q) ∈ D represents J * φ . By replacing I and (J * φ ) in equation ( 6) by their exact Zernike reconstruction (4), we obtain eq. ( 7) (see next page) where Z I pq and Z J pq represent Zernike moments of images I and J, respectively, and V pq , V * uv denotes the scalar product of two Zernike basis functions over the unit disc. Thanks to the orthogonality of the basis, this product is null except for the case where (p, q) = (u, v). In that case, it can be simplified into:

V pq , V * pq = π p + 1
At first glance eq. ( 7) is not trivial to minimize, but it can be rewritten into eq. ( 8) (see next page, with |Z pq | and [Z pq ] respectively the modulus and the argument of Z pq ). Formula (8) points out that the only parameter whose the distance depends is φ, which in addition is only present into the cosine functions. As a consequence, the search for optimal distance and angle will result in minimizing a sum of cosines.

Even if in real applications only a subset of Zernike moments is used (as it brings robustness to the description), the proposed method remains valid: simply, it can be seen as a fast way of retrieving the Euclidean distance between the two blurred patterns (that is, their reconstruction from the subset of moments) using their projection in Zernike space. The next section focuses on resolving the minimization so as to insure a low complexity and a fast computing time.

d 2 I,J (φ) = x 2 +y 2 ≤1 (p,q)∈D Z I pq .V pq (x, y) - (p,q)∈D Z J pq e iqφ .V pq (x, y) 2 = x 2 +y 2 ≤1 (p,q)∈D (Z I pq -Z J pq e iqφ )V pq (x, y) 2 = x 2 +y 2 ≤1      (p,q)∈D (Z I pq -Z J pq e iqφ )V pq (x, y)     (u,v)∈D (Z I uv -Z J uv e ivφ )V uv (x, y)   *    = x 2 +y 2 ≤1    (p,q),(u,v)∈D 2 (Z I pq -Z J pq e iqφ )V pq (x, y) (Z I * uv -Z J * uv e -ivφ )V * uv (x, y)    = (p,q),(u,v)∈D 2    (Z I pq -Z J pq e iqφ )(Z I * uv -Z J * uv e -ivφ )   x 2 +y 2 ≤1 V pq (x, y)V * uv (x, y)      = (p,q),(u,v)∈D 2 (Z I pq -Z J pq e iqφ )(Z I * uv -Z J * uv e -ivφ ). V pq , V * uv if (p, q) = (u, v), 0 else = (p,q)∈D Z I pq -Z J pq e iqφ 2 V pq , V * pq ( 7
)
d 2 I,J (φ) = (p,q)∈D π p + 1 Z I pq -Z J pq e iqφ 2 = (p,q)∈D π p + 1 Z I pq e i[Z I pq ] -Z J pq e i(qφ+[Z J pq ]) 2 = (p,q)∈D π p + 1 e i[Z I pq ] Z I pq -Z J pq e i(qφ+[Z J pq ]-[Z I pq ]) 2 = (p,q)∈D π p + 1 e i[Z I pq ] 2 . Z I pq -Z J pq e i(qφ+[Z J pq ]-[Z I pq ]) 2 = (p,q)∈D π p + 1 Z I pq -Z J pq cos(qφ + [Z J pq ] -[Z I pq ]) -i Z J pq sin(qφ + [Z J pq ] -[Z I pq ]) 2 = (p,q)∈D π p + 1 Z I pq -Z J pq cos(qφ + [Z J pq ] -[Z I pq ]) 2 + Z J pq sin(qφ + [Z J pq ] -[Z I pq ]) 2 = (p,q)∈D π p + 1 Z I pq 2 + Z J pq 2 -2 Z I pq Z J pq cos(qφ + [Z J pq ] -[Z I pq ]) (8)

EFFICIENT MINIMUM SEARCH

In this section, we describe how to efficiently extract the global minimum of eq. ( 8) by restricting the search using Nyquist-Shannon sampling theorem.

Assuming that the Zernike moments of each image are known until the N th order, the sum initially comprises O(N 2 ) cosine terms. Nonetheless, the sum can be simplified by removing the constant terms and by aggregating the cosine terms that own the same frequency:

A 1 cos(qφ + B 1 ) + A 2 cos(qφ + B 2 ) = |C| cos(qφ + [C]) ( 9
)
where C is a complex that worth A 1 e iB1 + A 2 e iB2 . Then, eq. ( 8) can be equivalently expressed as a sum of N cosines :

f N (φ) = N q=1
A q cos(qφ + B q ) (10)

with A q ∈ R + and B q ∈ [-π, π[.

Restricting the search of a global minimum

One can notice that f N (φ) is a 2π-periodic function.

Usually, the general technique for finding the minimum of a periodic function is a gradient descent. Indeed, f N (φ)'s first and second derivatives are easy to compute. However, such functions generally have many local minima whereas our approach requires to find the global one. One expensive solution is then to find every local minima and maxima with the gradient descent method by following the function from φ = 0 to 2π. However, f N (φ) owns a discrete Fourier spectrum bounded by a maximal frequency f M AX N = N/2π. Hence, f N (φ) has at most N local maxima and N local minima in [0, 2π[. Moreover, the Nyquist-Shannon sampling theorem teaches us that the function can not change substantially between two consecutive sampling points taken at the Nyquist frequency F = 1/T = N/π. The minimal distance between two consecutive minima is thus bounded by π/N . The initial starting points for finding every possible minima with a gradient descent can thus be equally scattered in 2N points. Moreover, by cutting f N into 4N intervals, we ensure that only one minimum or one maximum is present in each interval (see fig. 2).

Optimized minima retrieval

Our approach for optimizing the gradient descent takes advantage of the previously formulated properties. f N is sampled by 4N points equally spread between [0, 2π[: {x n = nπ/2N | 0 ≤ n < 4N }. We compute f N 's differential, denoted as f N , for each of those points. Section 4.1 ensures that if and only if a minimum is present between two consecutive points [x n , x n+1 ], then f N (x n ) is negative and f N (x n+1 ) is positive (see fig. 2). Moreover, those differential values enable to approximate the abscissa of the local minimum. Indeed, by approximating f N between [x n , x n+1 ] as a second degree polynomial, then the minimum abscissa can be evaluated at: For our application, the gradient descent algorithm does not need to be iterated to reach a high precision since this simple approximation (represented as crosses in fig. 2) is precise enough for our purpose (cf. experimental results in 5.1). Finally, the computational complexity of our approach for the distance minimization using the approximation is O(N 2 ). Our approach has then the same computational complexity than the classical Zernike comparator (eq. ( 5)).

x minimum = x n+1 f N (x n ) -x n f N (x n+1 ) f N (x n ) -f N (x n+1 ) = x n + π 2N f N (x n ) f N (x n ) -f N (x n+1 ) (11)

EXPERIMENTAL RESULTS

This section describes several experiments that illustrates:

• The efficiency of our minimum search (section 5.1). • A comparison of our approach with the classical one in terms of similarity accuracy, and retrieval performance (section 5.2.2). The average RMS errors corresponding to the approximation of the minimum position for various N .

• A comparison of our approach with the state of the art angle estimation algorithm (section 5.3).

Efficient minimum search

In order to demonstrate the efficiency of our minimum search algorithm (see section 4), a set of 30,000 random functions f N (φ) (see eq. ( 10)) has been created. For each of three different bandwidths N = {6, 12, 24}, there are 10,000 functions with random A p and B p . Firstly, we have computed for each function the exact solution φ and the approximation φ using respectively an exhaustive gradient descent and the proposed optimization. Secondly, we have compared both results and distinguished two types of situation:

1) The global minimum is correctly found by our approach: φ ≈ φ. 2) Another minimum is found. The second case derives from situations where the function admits more than one solution: In term of Zernike description, when there are s-1 secondary minima as small as the global one, that concretely means that the pattern is s-fold symmetric. For instance, a circular symmetric (2-fold symmetric) pattern can be rotated in an equivalent way by two different angles {α, α + 180 • }. In reality, what truly determines the rightness of the retrieved angle is the vertical error on the depth of the approximated minimum.

Table 1 details the RMS errors on the angle φ and on the depth of the minimum f N ( φ) = fN for both situations. The proportion of occurrences of the second case (the retrieved minimum is not the global one) is comprised below 0.15%. In the first case, the RMS errors of the angle and the minimum depth do not run over 0.81 • and 0.2%, respectively, even for N = 6 (when the approximation is the coarser). In the second case, the RMS error of the angle is high because a different minimum is retrieved. However, the depth of this minimum is similar to that of the global one since the maximum RMS error of fN is comprised below 1% in the worst case. That means that the minimum found may not be the global one but is really close in terms of depth.

As a conclusion, our minimum search seems precise enough (less than 1°of error for most of the cases) even for the case N = 6 for which the approximation is the coarser; this correspond to 16 moments. However the appropriate number of moment depends on the application, since it will also influence the robustness of the similarity measure.

Comparative study of similarity accuracy

We have conducted experiments to test the efficiency of the proposed comparator with respect to the state of the art, in terms of similarity accuracy and robustness. To that aim, we have gathered 502 logo images (about 25% are binary images and others are gray level images). Some of them are shown in figure 3. These images were cropped and re-sized in order to fit the size of a square of 100 × 100 pixels. Then, synthetic distorted images were generated for each of these original patterns. What is denoted as distortions include additive uniform noise, non-affine geometric deformation 1 , occlusion and translation. Examples of such distorted patterns at various levels are displayed in figure 4. These four types of distortion, applied at 19 different degrees, have been used to create four corresponding databases. In each database, the set of one original pattern and its 19 gradually distorted versions constitutes one class. Hence, one can see each test database as a collection of 10,040 different patterns divided into 502 classes of 20 elements.

Databases construction

For each test database, the distortion is applied gradually among 19 levels. The four types of distortion are:

1) An additive uniform noise added to the patterns such that the SNR varies from 30dB We used bi-linear interpolation for non-integer values of displacement. 3) A partial occlusion: we chose another random pattern and we paste it on the original one such that the occluded part varies from 0 (level 0) to 47.5 percents (level 19). The white parts of the pasted pattern are set transparent. 4) A translation. Its direction is random and its amplitude is comprised between 0 (level 0) and 7 pixels (level 19).

Comparison with the classical Zernike comparator

In order to conduct an exhaustive evaluation of our approach with respect to the classical way of comparing two Zernike descriptors, we have considered the similarity measure evaluation as a classification problem like in [START_REF] Khotanzad | Invariant image recognition by zernike moments[END_REF]. We have made two types of measurement using both comparators on the four databases:

• The recall-precision graph considering the previously defined classes. This measure accounts for the classification capability of our comparator and allows its complete evaluation for indexing issues. • The error rate on first retrieval: knowing a distorted pattern, we try to retrieve the original one. Contrary to the recall-precision measurement, this one corresponds to a recognition scenario where only the first retrieval is important. The Zernike moments of the patterns are extracted from the center of mass of the patterns in the case of noise and deformation, and from the center of the pictures for the cases of occlusion and translation since it has no sense to compute a center of mass when the pattern is occluded. All experiments were performed for three different numbers of moments : 25 (up to 8 th order), 49 (12 th order) and 81 (16 th order), which are representative of the numbers of Zernike moments used in the literature. In all cases, results were about the same so we decided to only illustrate the 12 th order case.

Recall Vs Precision:

A good way of measuring the efficiency of our approach in a global indexing framework is to evaluate the recall vs. precision of both comparators (see fig. 5). The experiment consists in comparing a pattern P of the database to every others, then to order the results by similarity, and finally to count the number A of relevant patterns (i.e. from the same class than P ) in the first N retrieved patterns. The precision is then equal to A/N and the recall to A/A max , with A max the size of the corresponding class. The results (averaged over all patterns from all the classes of the corresponding database) are illustrated in fig. 5 for each database (noise, deformation, occlusion, translation); the precision is presented for each value of recall. Practically, the higher is the curve, the better is the similarity measure. The first tier (i.e. recall value for N = A max ) and second tier (i.e. recall value for N = 2A max ) measures have also been computed for each graph (see table 2).

As shown by the curves, the proposed comparator performs always better than the classical one. The gain is particularly high for additive noise and deformation, where the proposed method performs up to respectively 6.8% and 5.6% better in terms of the first tier measure.

Recognition performances:

We have also conducted recognition experiments where the objective, starting from each of the 10,040 images of each database, is to retrieve the correct corresponding original pattern among the 502 original ones. The figure 6 displays the percentage of recognition errors (an erroneous pattern is returned in first position) for the classical and proposed comparator as a function of the distortion level. Globally, the proposed comparator is making about twice less errors than the classical one for medium distortions, except for the case of occlusion where our method makes about 15% less errors only. Fig. 6. The proportion of recognition errors when trying to retrieve the original pattern (among the 502 ones) from a distorted version using the classical comparator and the proposed comparator. The percentage of errors is highlighted for each method at the 15 th distortion level.

In conclusion, the use of both magnitude and phase during comparison much improves the efficiency of the descriptor in terms of classification and recognition. Indeed Recall-Precision measures have shown that the proposed method better classifies the set of patterns than the classical comparator in spite of huge distortions. Besides there were about twice less recognition errors, in our second experiment, than with the classical method for medium distortion strengths.

Comparison with geometric hashing and generalized Hough transform

We have also compared our approach with two state-of-the-art methods for sketch and object recognition: the geometric hashing [START_REF] Wolfson | Geometric hashing: an overview[END_REF] and the recent deformation tolerant generalized Hough Transform from Anelli et al. [START_REF] Anelli | Deformation tolerant generalized hough transform for sketch-based image retrieval in complex scenes[END_REF]. The geometric hashing is a well-known indexing technique that is used to quickly find matches between two sets of features (e.g. points). On the contrary, the generalized Hough transform is an object identification technique suitable for matching a shape contour model with unsegmented images: edges are extracted in the image and each edge point casts a vote in the space of the parameters (x and y position of the object to find). Anelli et al. have added two extra steps after the voting phase: 1) the votes are clustered in order to deal with small local deformations and then 2) the shape segmentation is further verified by back-projecting the image segments on the shape model and computing the model coverage score.

• For the geometric hashing, we need feature points, hence we have extracted Harris [START_REF] Harris | Geometry from visual motion[END_REF] and DoG keypoints [START_REF] Crowley | Representation for shape based on peaks and ridges in the difference of low-pass transform[END_REF], [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF] in each pattern.

An example of pattern with its keypoints is presented on figure 7. The descriptor for a given pattern consists of a list of its 15 strongest keypoints in terms of response value. Then, we have trained a hash table for each type of distortion. In the retrieval step, we begin by extracting the 25 strongest keypoints of the pattern to identify (we take more than 15 because the presence of noise can add new keypoints), and we successively project them on the hash tables of the learned patterns for different random basis (50 times). The best scores are stored for each learned pattern and it constitutes at the end the final matching scores. We deliberately chose small values for the number of keypoints and iterations so that a comparison is tractable in terms of time (see table 3). • In the case of the deformation tolerant GHT (DT-GHT) from Anelli et al. [START_REF] Anelli | Deformation tolerant generalized hough transform for sketch-based image retrieval in complex scenes[END_REF], each pattern was first indexed (i.e. canny edges extraction, segmentation of the edges, building of a Rtable ). Then to compare two patterns, we simply search the first one (the model) into the second one. This returns a matching score between 0 and 1 (1 means perfect match).

In order to properly compare our Zernike comparator with these methods, we have evaluated the recall vs. precision scores. However the comparison steps for geometric hashing and Hough transform are both quite slow so it was not possible to process the 10, 000 2 comparisons required for the recallprecision measures. As a consequence, we deliberately sampled each database into 5 deformations levels (3 rd , 7 th , 11 th , 15 th and 19 th ) × 100 patterns instead of 502, with each type of pattern (gray level / black & white) equally sampled. Each of the four 3 gives a comparison of the similarity calculation processing time for each method. Even if our method is slower than the classical Euclidean comparison, its complexity is the same (O(order 2 ) ) and both processing times are very low (less than 60 milliseconds for 1000 comparisons for our approach). However the DT-GHT and the geometric hashing are several orders of magnitude slower (37 and 27 seconds respectively for 1000 comparisons). With a number of comparisons of 500 2 per database for the recall-precision experiment, it still required about 10 hours to process the four databases for the DT-GHT or the geometric hashing against less than one minute for Zernike experiments. and DT-GHT approaches perform uppermost due to their invariance in position (however the geometric hashing is not so good because of the interpolation noise), but it can be noticed that the proposed Zernike method does not perform so bad (almost 80% on the first tier measure). For every other distortions, our approach performs better than DT-GHT and geometric hashing in terms of the first and second tier measures ; in the case of noise, the proposed method works largely better since the noise is localized in the high frequencies and hence has a small effect on the low order Zernike moments. The geometric hashing achieves low performance in all cases except for translation since the keypoint response values and locations are strongly disturbed when the images are distorted.

Rotation angle retrieval

We have conducted experiments to demonstrate the accuracy of our approach in the specific case of rotation angle estimation between similar but distorted patterns. We have performed the experiment on the previous set of 502 logo images. The images were rotated by random angles. Then, various distortions (the same as in section 5.2) were applied to each rotated pattern. The distortion levels, however, are twice smaller than in the previous section otherwise images would have been too much different to be put in correspondence by a rotation. Then, we have estimated the rotation angle with respect to the original pattern using the proposed algorithm and the most robust and acknowledged method in the state of the art: the Kim and Kim robust estimator [START_REF] Kim | Robust rotation angle estimator[END_REF]. The total number of test images is thus 200,800: 502 (the number of original patterns) × 10 (the number of rotations) × 10 (the number of distortion levels) × 4 (the number of types of distortion).

To compare to the estimator from Kim and Kim [START_REF] Kim | Robust rotation angle estimator[END_REF], we have used the same number of Zernike moments than them: 25 moments were computed up to order 8. For each of the four databases (noise, deformation, occlusion, translation), we have computed the root mean square (RMS) error of the retrieved angles using each algorithm; results are presented in table 5. The RMS error using our method is systematically lower than for the estimator from Kim and Kim, for all the types of distortion. In particular, in the case of additive noise and deformation, the proposed method provides results whose accuracy is about twice better than with the estimator from Kim and Kim. The occlusion and translation cases yields an important RMS error: in the first case, this comes from the fact that sometimes some crucial parts of the pattern in terms of rotation evidences are occluded ; in the second case this is caused by the displacement of the rotation center.

2D-3D OBJECT RECOGNITION FROM HAND-DRAWN SKETCH

The proposed algorithm has been tested within a real industrial application: the automatic transcription of sketched storyboards into reconstructed

ABLE 5

The average RMS error of rotation angle using the estimator from Kim and Kim [START_REF] Kim | Robust rotation angle estimator[END_REF]and the proposed method. 3D scenes. An example of such scenario is illustrated in figure 9: a 3D scene has been generated from a hand-drawn storyboard. Formally, an object database contains various 3D models (as polygonal meshes); on the other side, the cartoonist draws a 2D storyboard by digitally sketching the objects of the scene he has in mind. Then the objective is to recognize each 3D model from the corresponding piece of sketch, along with its 3D viewing angle, its scale and its rotation angle in the drawing plane, so as to automatically and correctly place it in the 3D scene. To obtain this result, an important number of views of each 3D model is indexed using Zernike moments like in [START_REF] Filali Ansary | A bayesian 3D search engine using adaptive views clustering[END_REF] (about 50 views per object). We do not compute Zernike moments on the whole view, but on a bounding circle containing at least 70% of the pattern. We limit the indexed surface for the following reasons:

1) we try to reduce the surface potentially disturbed by the background during retrieval 2) such approach may bring some kind of robustness regarding small occlusions. A complete robustness to occlusion could be reached if each model view was described by several Zernike circles of various scales and positions, however, this development is out of the scope of this paper.

Once we have described all the 3D model views, we search them in the sketched story-board with a twopass process to speed up the recognition. The first pass aims at finding a set of raw correspondences and the second pass refines this set. Without loss of generality, we will now describe the search of one given model into the sketch: in the first pass, we scan the whole sketched story-board with a circular window at different scales and positions.

For each position, we describe the window using Zernike moments and we store the model view which achieves the smallest distance with the local descriptor according to our comparator. We use a standard non-maxima suppression technique in the resulting scale-space of distances and we apply a first threshold to the obtained maxima: we obtain a set of potential matches. This first search is processed for a reduced set of positions and scales in order to reduce the processing time: 8 scales, starting from 400 × 400 pixels until 1400 × 1400 pixels, and for each scale the window slides over the image with a step width of windowSize/15 (since that corresponds to a translation whose amplitude is correctly retrieved in section 5). In our experiment, the story-board dimension is 3350 × 2260 and this greedy recognition takes about 2 minutes on a 2Ghz machine. The computation of Zernike moments is made faster by precomputing a set of 100×100 Zernike filters and by applying a fast smoothing approach along scales (like pyramids of Gaussian in [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF]) to quickly sample each 100×100 window. As we saw in section 5.2, Zernike distance is not so invariant to translation, so we refine the search during a second pass: for each potential match we search locally around its position in the scale-space (the steps widths are decreased both in position and scale). This allows to find the optimal position for each potential match and thus to make the difference between true and false positive after a second thresholding. This second pass takes about 30s per potential match (there is typically 4-6 potential matches per storyboard image).

Moreover, we assume that the 3D models are approximately vertically positioned on an horizontal flat ground, which means that their 3D vertical axis are probably also vertical in the picture. Since our method provides not only the similarity but also the the rotation angle, we are able to eliminate from the first pass further false positives which could not be eliminated by the classical Zernike comparator alone. Recognition results are presented in the figure 10. Note that all the bushes are not detected: this is caused by the fact that sometimes their sketch deviates too much from their 3D model. Even if there are slight errors in the placement estimation, the recognition results are still acceptable.

For such an application, the estimations are significantly better and more stable with the proposed comparator than with the classical comparator cou-pled with the in-plane rotation estimator from [START_REF] Kim | Robust rotation angle estimator[END_REF], which is less reliable when it has to estimate a rotation angle between two different patterns (see 5.

Fig. 1 .

 1 Fig. 1. The best four retrievals for the drawing of a rooster (a) from a database of 200 binary logos. Results are ordered by similarity measure (a distance of zero means a complete similarity) for both comparators.

Fig. 2 .

 2 Fig. 2. A random function f 3 cut into 4N = 12 intervals. Each interval contains at most one local minimum. When a minimum exists in the interval, the left derivative is negative and the right derivative is positive. For each minimum, a black cross figures the minimum position approximated with second order polynomial.

Fig. 3 .

 3 Fig. 3. Some of the 502 patterns used in the experiments.

(

  level 0) to 1.5dB (level 19). The resulting values are bounded in [0,255]. 2) A non-affine smooth deformation (i.e. a geometric deformation) whose amplitude varies locally from 0 (level 0) to 25 pixels (level 19) at most. This deformation is generated by two 100x100 displacement maps Dx and Dy that are initialized to 0, except for 8 random cells where we create vectors of the given amplitude and random directions. These 8 deformation vectors are then diffused until convergence (principle of the heat diffusion). The deformation function consists in Def ormedP attern(x, y) = P attern(x + Dx(x, y), y + Dy(x, y)).

Fig. 4 .

 4 Fig. 4. Illustration of the four types of distortion. From left to right the four groups of distortions: additive noise, deformation, occlusion, translation. In each group the 6 th , 13 th and 19 th distortion levels are presented. The images are inverted (black becomes white) for visualization purposes.

Fig. 7 .

 7 Fig. 7. 15 strongest keypoints of a pattern detected by the Harris corner detector (blue crosses) and the extrema in scale-space of Dog (red circles)

Fig. 8 .

 8 Fig. 8. Recall-Precision curves on the four reduced databases for the following methods: proposed Zernike, DT-GHT, geometric hashing. Each database contains 500 patterns divided into 100 classes of one original and 4 distorted versions.

Fig. 9 .

 9 Fig.9. Example of 2D/3D object recognition application: reconstruction of a 3D scene from a handdrawn 2D storyboard associated with a 3D model database. On the top, the recognized 3D model (on the right) is the best retrieved result between the drawing (on the left) and the whole set of indexed views, using the proposed Zernike comparator. Doing so, the 3D pose and the in-plane rotation angle are retrieved in the same time and thereafter allows the 3D reconstruction.

Fig. 10 .

 10 Fig. 10. Recognition results for the circus tent, the trailer and the bushes pasted on the sketch image after edge extraction. No human intervention was needed to obtain this result.

TABLE 1

 1

  NoiseDeformation Occlusion Translation 1 st tier 2 nd tier 1 st tier 2 nd tier 1 st tier 2 nd tier 1 st tier 2 nd tier

	Classical comparator	82.9%	86.2%	77.9%	84.8%	44.2%	50.8%	68.7%	75.9%
	Proposed comparator	89.7%	92.0%	83.5%	88.9%	47.2%	53.3%	71.6%	78.3%
	Difference	+6.8%	+5.9%	+5.6%	+4.1%	+3.0%	+2.5%	+2.9%	+2.4%

TABLE 2

 2 Comparison of the first tier and second tier measures for each database.

TABLE 3

 3

	Average processing times for 1000 comparisons on
	a 2 GHz computer for similarity calculation with
	each method
	databases thus contains 500 different elements.
	5.2.3.1 Timing Analysis: Table

  5.2.3.2 Recall Vs Precision: Figure 8 and table 4 illustrate the results. In the case of translation, the geometric hashing Noise Deformation Occlusion Translation 1 st tier 2 nd tier 1 st tier 2 nd tier 1 st tier 2 nd tier 1 st tier 2 nd tier

	Zernike proposed 93.4%	95.6%	88.6%	93.6%	52.9%	60.6%	78.0%	85.2%
	DT-GHT	75.5%	79.4%	88.4%	93.6%	43.5%	49.4%	96.5%	98.4%
	Hashing	62.0%	63.6%	57.2%	60.2%	47.8%	49.9%	86.9%	87.4%

TABLE 4

 4 Comparison of the first tier and second tier measures for each comparators and for each distortion.

SUMMARY AND DISCUSSION

We have presented an efficient comparator of Zernike descriptors whose novelty is to take advantage of the phase information in the comparison process while still preserving the invariance to rotation. The provided similarity measure is more robust to distortions (especially geometrical deformation and noise). The recognition errors are also about twice less important for medium distortions than with the classical comparator (i.e. the Euclidean distance between Zernike magnitudes). Moreover, our approach has the same O(order 2 ) complexity than the classical one. Finally, it provides for free an estimation of the rotation angle that outperforms the robust estimator from Kim and Kim [16]. To conclude with, this novel theoretical contribution to the Zernike framework can apply to any application already using Zernike moments. It is worth noting that it can also apply to the pseudo-Zernike moments [22] as both theories share almost the same background.


Hichem Barki
email: hichem.barki@liris.cnrs.fr

Florence Denis
email: florence.denis@liris.cnrs.fr

Florent Dupont
email: florent.dupont@liris.cnrs.fr

Contributing vertices-based Minkowski sum computation of convex polyhedra

Keywords: Minkowski sum, contributing vertices, slope diagram, convex hull, computer-aided design

Introduction

Minkowski sum implementation is of a particular interest and used in a variety of domains such as computer-aided design and manufacturing [START_REF] Lee | Polynomial/rational approximation of Minkowski sum boundary curves[END_REF], computer animation and morphing [START_REF] Kaul | Solid-interpolating deformations: construction and animation of PIPs[END_REF], morphological image analysis [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], similarity measures for convex polyhedra [START_REF] Tuzikov | Similarity measures for convex polyhedra based on Minkowski addition[END_REF], penetration depth computation and dynamic simulation [START_REF] Kim | Fast penetration depth estimation using rasterization hardware and hierarchical refinement[END_REF], robot motion planning [START_REF] Lozano-Pérez | Spatial planning: a configuration space approach[END_REF], and solid modeling.

The Minkowski sum or addition of two sets A and B in a vector space was defined by the German mathematician Herman Minkowski (1864-1909) as a position vector addition of elements of A and elements of B:

A ⊕ B = {a + b|a ∈ A, b ∈ B} (1)
Where + denotes the vector addition of two position vectors a and b coming from the two sets A and B. The Minkowski sum can also be obtained by the following definition:

A ⊕ B = a∈A B a (2)
Where ∪ denotes the set union operation and B a denotes the set B translated by a vector a. The second definition states that the Minkowski sum of two sets

A and B is obtained by sweeping all points of A by B and taking the union of all resulting points. The sweep aims at positioning or translating B such that its origin (the common initial point of all its position vectors) coincides with each point of A and to take the union of the resulting sets, as depicted in Fig. A convex polyhedron P in R 3 can be seen as the intersection of a finite number of closed half-spaces P = i=1,...,n H - i . Where H - i is a closed half-space bounded by the plane H i having a normal vector u i , the half-space H - i is called the negative side of the plane H i (see Fig. 2 for an example in 2D).

The boundary of a polyhedron can be partitioned into two-dimensional faces (facets), one-dimensional faces (edges), and zero-dimensional facets (vertices).

The supporting plane of a particular facet of a polyhedron A is the plane where the facet lies. As an example, for the convex polygon P depicted in Fig. 2.b, the supporting line (the supporting plane becomes a supporting line in 2D) of the one-dimensional face (edge) e 3 is H 3 . Further details on support function representation of convex polyhedra and supporting planes can be found in [START_REF] Ghosh | Support function representation of convex bodies, its application in geometric computing, and some related representations[END_REF]. The rest of this paper is based on the definitions we gave above and is organized as follows:

• In section 2, we review existent techniques for the computation of Minkowski sum of polyhedra.

• In section 3, we give an overview of our work and introduce key concepts we use throughout all the paper, such as the concept of contributing vertices, the relation between contributing vertices and the sweep process, the translated facets, the corner facets, and the edge facets.

• In section 4, we present in detail our new algorithm, based on the concept of contributing vertices, for the computation of the Minkowski sum of convex polyhedra. We call it Contributing Vertices-based Minkowski Sum (CVMS) algorithm.

• In section 5 we explain how the contributing vertices concept can help improve slope diagram-based methods and present our Enhanced and Simplified Slope Diagram-based Minkowski Sum (ESSDMS) algorithm.

• Finally, we present results, performance study, and comparisons between our two algorithms and other approaches found in the literature.

Previous work

The large number of applications needing Minkowski sum has inspired researchers to propose and develop various algorithms. The most common approach used for polyhedra is based on convex hulls. For two convex polyhedra

A and B, it performs vector addition of all points of A and B and computes the convex hull of the resulting point set giving us the Minkowski sum polyhedron A ⊕ B. Convex hulls are computed by several algorithms such as: the Quickhull algorithm [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF], the incremental algorithm [START_REF] Seidel | Output-size sensitive algorithms for constructive problems in computational geometry[END_REF][START_REF] Kallay | The complexity of incremental convex hull algorithms in R d[END_REF], the gift-wrapping algorithm [START_REF] Chand | An algorithm for convex polytopes[END_REF], the divide and conquer algorithm [START_REF] Preparata | Convex hulls of finite sets of points in two and three dimensions[END_REF], and Graham's algorithm [START_REF] Graham | An efficient algorithm for determining the convex hull of a finite planar set[END_REF]. A summary of these algorithms can be found in [START_REF] O'rourke | Computational Geometry in C[END_REF].

For two non-convex polyhedra A and B, the computation of the Minkowski sum polyhedron A ⊕ B is done by decomposing each non-convex polyhedron into convex pieces, computing the the pairwise Minkowski sums of all possible pairs of convex pieces from A and B, and performing the union of the pairwise Minkowski sums. The computation of the Minkowski sum of a pair of convex pieces (or polyhedra) is an important step in this approach.

For readers interested in the Minkowski sum of non-convex polyhedra, we give some references about convex decomposition and union steps. Convex decomposition aims at decomposing a non-convex polyhedron into convex pieces.

The smaller the number of pieces, the better the decomposition. The optimal decomposition of a non-convex polyhedron into convex pieces is known to be NP-hard. More than twenty years ago, Chazelle [START_REF] Chazelle | Convex decompositions of polyhedra[END_REF] proposed an optimal decomposition algorithm which generates O(r 2 ) convex pieces in O(nr 3 ) time, where r denotes the number of reflex edges and n denotes the number of polyhedron facets. Nevertheless, no practical or robust implementation has been found in the literature for Chazelle's optimal algorithm. An experimental study of convex decomposition strategies can be found in [START_REF] Chazelle | Strategies for polyhedral surface decomposition: an experimental study[END_REF].

The union step is the most time consuming step when computing the Minkowski sum of non-convex polyhedra. It can have O(n 6 ) time complexity for non-convex polyhedra, where n is the number of polyhedra facets [START_REF] Aronov | The union of convex polyhedra in three dimensions[END_REF]. Moreover, there is no robust implementation for the union computation of convex polyhedra that handles all degeneracies [START_REF] Halperin | Robust geometric computing in motion[END_REF].

While the combinatorial complexity of the convex hull approach is O(mn)

for convex polyhedra with m and n features, it can have O(m 3 n 3 ) worst-case combinatorial complexity for non convex polyhedra [START_REF] Halperin | Robust geometric computing in motion[END_REF][START_REF] Dobkin | Computing the intersectiondepth of polyhedra[END_REF].

Recently, Hachenberger implemented a data structure in CGAL [START_REF]The CGAL project homepage[END_REF] based on Nef polyhedra theory [START_REF] Nef | Beiträge zur theorie der polyeder[END_REF]. He used this Nef polyhedra implementation [START_REF] Hachenberger | Boolean operations on 3d selective nef complexes: data structure, algorithms, optimized implementation and experiments[END_REF][START_REF] Hachenberger | Boolean operations on 3d selective nef complexes: optimized implementation and experiments[END_REF] for the decomposition of polyhedra into convex parts and the union computation of pairwise Minkowski sums [START_REF] Hachenberger | Exact Minkowski sums of polyhedra and exact and efficient decomposition of polyhedra in convex pieces[END_REF]. Unfortunately, this approach is time consuming because it is based on the decomposition of polyhedra into convex pieces.

Due to non-robustness and expense of the union operation, researchers tried to bypass it by using other representations or computing only approximations. Varadhan and Manocha [START_REF] Varadhan | Accurate Minkowski sum approximation of polyhedral models[END_REF] decomposed the polyhedra into convex pieces, computed the pairwise Minkowski sums of convex pieces (by means of convex hull algorithms), approximated the union of the pairwise Minkowski sums using an adaptively subdivided voxel grid, computed signed distance on the grid points and used isosurface extraction from the distance field [START_REF] Varadhan | Topology preserving surface extraction using adaptive subdivision[END_REF].

They guarantee that their approximation has the same topology as the exact Guibas et al. [START_REF] Guibas | A kinetic framework for computational geometry[END_REF][START_REF] Guibas | Computing convolutions by reciprocal search[END_REF] presented a kinetic framework in two dimensions. An accomplishment of this framework was to define the operation of convolution on planar tracings in 2D. Basch et al. [START_REF] Basch | Polyhedral tracings and their convolution[END_REF] extended the convolution operation to polyhedral tracings in 3D and used it to generate a superset of the Minkowski sum. They extracted the exact boundary of the Minkowski sum by using arrangement computations.

Ghosh [START_REF] Ghosh | A unified computational framework for Minkowski operations[END_REF] proposed a unified computational framework to compute the Minkowski sum of polygons and polyhedra (in 2D and 3D domains). An interesting algorithm is that proposed and implemented by Fogel and

Halperin [START_REF] Fogel | Exact and efficient construction of Minkowski sums of convex polyhedra with applications[END_REF]. The authors represented the convex polyhedra in a dual space they called the Cubical Gaussian Map (CGM) and implemented their algorithm on the base of the arrangement package of CGAL. The cubical Gaussian map is a Gaussian map embedded on a unit cube instead of the unit sphere.

The Minkowski sum of two convex polyhedra A and B is obtained by computing CGM (A) and CGM (B), the respective cubical Gaussian maps of A and B, overlaying or merging CGM (A) and CGM (B), and building the Minkowski sum polyhedron S = A ⊕ B from the resulting overlay. The overlay process has a time complexity of O(f S log(f A + f B )), where f A , f B , and f S are the respective facets number of the polyhedra A, B, and S.

Overview of our work

In this work, we are interested in computing the Minkowski sum of convex polyhedra. We present a novel algorithm based on the concept of contributing vertices associated to every facet of the two polyhedra we are working with, hence its name : Contributing Vertices-based Minkowski Sum (CVMS) algorithm.

Introduction to the CVMS algorithm

The CVMS algorithm is inspired from the definition given in equation 2. This For two convex polyhedra A and B, the same logic applies. So, sweeping a particular facet f i,A of A by B is done by parallel-translating B such that its origin (the common initial point of all its position vectors) passes through all points of f i,A . This sweep implies that all vertices of B are also parallel- translated in the same manner, they are moving in planes parallel to the supporting plane of f i,A and generating facets within supporting planes parallel to that of f i,A . But only the facet generated by the vertex v k,B of B which is at maximal distance away from the supporting plane of f i,A is taken into account because it lies on the boundary of the sum polyhedron. The other facets generated by the displacement of all other vertices of B are discarded since they lie in the interior of the sum polyhedron. This particular vertex v k,B which generated this facet of the sum polyhedron A ⊕ B is called the "contributing vertex" associated to the facet f i,A .

B A origin of B (a) (b) (c) f 1,A A f 1,A f 2,A f 3,A f 4,A v 1,B v 2,B

Computation of the Minkowski sum polyhedron from contributing vertices

To construct the Minkowski sum polyhedron, we must perform three sweep The last step is to sweep the edges of A by B. For a particular edge e j,A of A, this operation implies the construction of several facets that are the Minkowski sum of the edge e j,A and all the edges of B. But not all these facets will lie necessarily on the boundary of S. Therefore, some of these facets are discarded. To determine which of these facets are retained, we use two criteria: a visibility criterion and normal orientation criterion as explained in section 4.3. The retained facets are denoted "edge facets". Note that only edges of A or B incident to facets having distinct contributing vertices will be considered in this step. A proof of this property is given in section 4.

These three steps of sweep enable us to construct the Minkowski sum polyhedron S. We show in section 4 how to compute the contributing vertices associated to the facets of A and B. This resumes to the computation of distances from vertices to supporting planes. We also show that the sum polyhedron S is composed of three (already mentioned) types of facets: "translated facets" obtained from the facets of A and their associated contributing vertices, "corner facets" obtained from the facets of B and their associated contributing vertices (commutativity property), and "edge facets" obtained from the Minkowski sum of two non-parallel edges one from A and the other from B, having incident facets with no common contributing vertices. An outline of the CVMS algorithm is presented in algorithm 1. Further details can be found found in section 4.

The CVMS algorithm for convex polyhedra

In the rest of this paper, we will not distinguish between a vector v and a point v. We will talk about a vertex v and simply refer, in the equations, to the vector pointing from the coordinates origin o to the point v as v. The distinction between a point and the corresponding vector is clear from the context.

First, let's give some notation and definitions required for the understanding of the following notions.

Notation

We consider two convex, closed and two-manifold polyhedra A and B. A is composed of f A facets, e A edges, and v A vertices. Similarly, B is composed of f B facets, e B edges, and v B vertices. The Minkowski sum of A and B is denoted S.

Prerequisites and definitions

First, let's define the concept of contributing vertices since it is needed in the rest of definitions.

Definition 1: The contributing vertex v k,B of a particular facet f i,A with an outer normal n i,A is the vertex of B which is at maximal distance away from the supporting plane of f i,A . Formally, the contributing vertex v k,B of a particular facet f i,A satisfies:

v k,B , n i,A = max R v l,B , n i,A ∀l = 1, 2, ..., v B (3)
Where ., . denotes the scalar product. An illustration of the contributing vertex concept in 2D is presented in Fig. 4.

The contributing vertex v k,A of a particular facet f i,B with an outer normal n i,B is defined in the same manner by putting A in place of B and B in place of A in equation 3.

For some facets of A and B, we can find many contributing vertices due to the fact that these contributing vertices are at the same (maximal) distance away from the supporting planes of the considered facets. This case is treated normally and does not cause any degenerate behavior. (1) Sweeping all the facets of A by B: when a particular facet f i,A is swept by B (see Fig. 5.a), this results in the facet of the sum polyhedron S having a supporting plane parallel to that of f i,A since it is generated by the contributing vertex associated to f i,A (i.e. the vertex v k,B which is at maximal distance away from the supporting plane of f i,A ). So, the sweeping of all f A facets of A results in the f A translated facets with supporting planes parallel to those of the facets of A (see Fig. (3) Sweeping all edges of A by B: when a particular edge e j,A is swept by B, i.e. the origin of B is translated so it passes through all points of e j,A (see (1) f A facets with supporting planes parallel to those of the facets of A, these facets are named the "translated facets" of S;

n 1,A n 2,A n 3,A n 4,A n 5,A v 1,B f 2,A f 3,A f 4,A f 5,A f 1,A v 2,B v 3,B v 4,B v 5,B v 2,B v 3,B v 4,B v 5,B
(2) f B facets with supporting planes parallel to those of the facets of B, these facets are named the "corner facets" of S;

(3) At most e A e B facets that result from the Minkowski sum of some pairs of non-parallel edges of A and B, these facets are named the "edge facets" of S.

Proof: Let's sweep all f A facets of A by B. For a particular facet f i,A , i = 1, . . . , f A we keep only the facet that is the result of translating f i,A according to its contributing vertex or vertices. So, there is a one-to-one mapping between the facets of A and the facets of S with supporting planes parallel the edge e k,A will be incident to f i,S and f j,S and will never contribute in any edge facet construction.

In contrast, if another edge e k,A is incident to facets f i,A and f j,A having distinct common contributing vertices v p,B and v q,B respectively, it is clear that the construction of the translated facets f i,S and f j,S corresponding to f i,A and f j,A involves v p,B and v q,B respectively. Therefore, f i,A is translated by a vector starting from coordinates origin o and ending at v p,B and f j,A is translated by a vector starting from o and ending at v q,B . Thus, the vertices v m,A and v n,A incident to the edge e k,A are translated to the positions v m,A + v p,B and v n,A +v p,B when computing f i,S and to the positions v m,A +v q,B and v n,A + v q,B when computing f j,S . Therefore, the translated facets f i,S and f j,S are not incident to the same edge of polyhedron S (they are disjoint)

and the edge e k,A creates one or more edge facets.

The proof of property 1 can be done in a similar manner for any edge e k,B of B.

In Fig. 7, a simple example illustrates property 1. The edge e 1,A in Fig.

v 1,A v 2,A v 1,A + v 1,B v 2,A + v 1,B f 1,S f 2,S f 3,S v 2,A + v 2,B (c)

Steps of the CVMS algorithm

The CVMS algorithm begins by determining the translated and corner facets of S by computing the contributing vertices associated to all facets of A and B respectively. Then, it determines edge facets. Property 1 is used in the translated and corner facets determination steps and property 2 is used in the edge facets determination step. Given two convex and closed two-manifold polyhedra A and B, the CVMS algorithm performs as follows:

(1) Translated facets determination: for each facet f i,A :

• Find its associated contributing vertices v k,B (see equation 3).

• If there are more than two contributing vertices, the facet f j,B incident to the contributing vertices associated to f i,A lies on a supporting plane parallel to that of f i,A . The corresponding translated facet is the result of the planar Minkowski sum of f i,A and f j,B (f i,A ⊕ f j,B ). The translated facet has the same outer normal orientation as f i,A and f j,B . The facet f j,B will be ignored in the corner facets determination step since it contributes once.

• If there are exactly two contributing vertices, the edge e j,B incident to the two contributing vertices associated to f i,A lies on a supporting line parallel to the supporting plane of f i,A . The resulting translated facet is computed by a planar Minkowski sum of f i,A and e j,B (f i,A ⊕ e j,B ) and has the same outer normal orientation as the facet f i,A .

• If there is only one contributing vertex v k,B , the corresponding translated facet results from translating the facet f i,A by a position vector pointing from coordinates origin o to v k,B and has the same outer normal orientation as the facet f i,A .

(2) Corner facets determination: for each facet f i,B that have not yet contributed in the previous step:

• Find its associated contributing vertices v k,A . Note that there is no case where a facet f i,B has more than two contributing vertices since such a facet has contributed once in the translated facets determination step.

• If there are exactly two contributing vertices, the edge e j,A incident to the Therefore, an edge facet e i,A ⊕ e k,B is added to the sum polyhedron S if and only if:

n 1,A × n i,k , e i,A < 0 and n i,k × n 2,A , e i,A < 0 (5)
At this point, we have constructed the sum polyhedron S. We note that we have not talked about the origin of the polyhedron B, but it is easy to consider an origin point c that is different from the coordinates origin o by translating the sum polyhedron S by a vector beginning from c and ending at o. After that, we describe the ESSDMS algorithm and indicate the role played by contributing vertices in the improvement of performance (by eliminating unnecessary processing) and the simplification of some steps.

Original slope diagram-based Minkowski sum algorithm

Slope diagrams were proposed for the first time by Ghosh [START_REF] Ghosh | A unified computational framework for Minkowski operations[END_REF]. The principle is to represent two convex and closed two-manifold polyhedra A and B in a dual space such that the orientation of the normals of the facets, edges, and vertices of the two operands determine the Minkowski sum polyhedron S. The dual space used by Ghosh is a unit sphere. For a convex polyhedron A, the slope diagram representation SD(A) is obtained as follows:

• Each facet f i,A with an outer unit normal n i,A is represented by a spherical point corresponding to the end point of the outer unit normal n i,A . This point is referred to as a facet point.

• Each edge e j,A is represented on the unit sphere by the arc of the great circle joining the two facet points representing the two facets of A incident to the edge e j,A . The spherical arc representing the edge e j,A is called an edge arc.

By "the arc of the great circle" we mean the shorter arc among the two ones joining the two facet points. The reason for choosing the shorter arc is that for convex polyhedra, the angle between two neighboring facets is always less than 180 degrees.

• Each vertex v k,A is represented by a region on the unit sphere (a spherical polygon) bounded by the facet points corresponding to the facets of A incident to the vertex v k,A and by the edge arcs corresponding to the edges of A incident to v k,A . We call this spherical polygon a vertex polygon or a vertex region.

The slope diagram of polyhedron B is obtained by following the same steps used for polyhedron A. Fig. 9 shows a slope diagram representation of a convex and closed polyhedron.

The Minkowski sum polyhedron A ⊕ B is obtained by merging the two slope diagrams SD(A) and SD(B). The merging process aims at deducing the sum facets following the intersections or coincidences of the two slope diagrams components (facet points, edge arcs, and vertex regions). Four intersection cases are detected: In order to merge the two slope diagrams embedded on the unit sphere, Ghosh transformed them into a 2D planar form by means of stereographic projection.

f 1 f 2 f 4 v 1 v 2 v 3 v 4
The stereographic projection is a projection of a sphere from one of the points N or N onto the plane σ tangent to the sphere in the diametrically opposite point N or N (see Fig. 10). Although Ghosh formulated slope diagram for the computation of the Minkowski sum of convex and non-convex polyhedra, the existent implementations are only restricted to convex polyhedra. Furthermore, the weakness of the slope diagrams based Minkowski sum algorithms come from the embedding on the unit sphere and from the stereographic projection used to merge two slope diagrams. These two steps are relatively complex, require a non-negligible computation time, and lead to precision problems when working with built-in number types such as float or double.

To overcome these two problems, we propose the ESSDMS algorithm which is somewhat similar to the improved slope diagram algorithm presented by Wu et al. [START_REF] Wu | Improvements to algorithms for computing the Minkowski sum of 3-polytopes[END_REF] and detailed in the next section.

The ESSDMS algorithm

The ESSDMS algorithm we propose here comes with simplifications and enhancements to the algorithm of Wu's et al.. It uses the concept of contributing vertices to gain additional performance and simplify some steps of the algorithm. Its implementation is based on exact number types instead of built-in number types to avoid using position thresholds. So the ESSDMS algorithm is similar to that of Wu et al. in the following points:

• No embedding of the slope diagram representation: the slope diagram principle is exploited without explicitly embedding objects in the unit sphere. When information about facet points, edge arcs, or vertex regions is needed, it is directly extracted from the polyhedra we are working with.

• No stereographic projection: because there is no explicit embedding or representation of slope diagrams, there is nothing to project. The intersections between slope diagram components are detected by the use of geometric operations obtained directly from polyhedra.

In contrast, the ESSDMS algorithm differs from that of Wu et al. in the following points:

• Use of the contributing vertices concept: this enhancement reduces the number of edges used in the determination of the intersections between edge arcs by considering only those having incident facets associated to common contributing vertices. Moreover, it allows to avoid the point-in-polygon queries needed when testing if a facet point lies inside a vertex region (intersection of a facet point and a vertex region). If we perform all other intersection queries (intersection of two facet points, intersection of a facet point and an edge arc, and the intersection of two edge arcs), it is clear that the remaining facets correspond to facet points lying necessarily inside vertex regions. So the corresponding facets in the Minkowski sum polyhedron are translated copies of these facets. By analogy with the CVMS algorithm, we can conclude that each of these facets has only one contributing vertex.

The creation of the sum facets corresponding to the remaining facets is done by using their contributing vertices. The computation of the contributing vertices does not introduce additional performance overhead (compared to the cost of detecting intersections) since they are deduced directly from the intersection configuration.

• No round-off errors: this is an implementation issue, by using exact number types (such those provided by GMP (GNU Multi Precision) library [START_REF]The GNU MP bignum library[END_REF],

we ensure that computations are done exactly and intersections are computed correctly. This avoids the use of position tolerance values or thresholds when testing coincidences of geometric features.

Given two convex, closed and two-manifold polyhedra A and B, the four intersection cases are handled as the following: The corresponding Minkowski sum facet f S for this case is the result of a planar Minkowski sum of f i,A and f j,B :

•
f S = f i,A ⊕ f j,B (7)
Since the facets f i,A and f j,B have parallel supporting planes, the contributing vertices associated to the facet f i,A are the vertices of f j,B and the contributing vertices associated to the facet f j,B are the vertices of f i,A .

•

f S = f i,A ⊕ e j,B (10)
In this case, the contributing vertices associated to the facet f i,A represented by the spherical point f are the two vertices incident to the edge e j,B

represented by the edge arc st.

By commutativity, if the facet point f corresponds to a facet f i,B and the edge arc on which f lies corresponds to an edge e j,A , the Minkowski sum facet f S is given by equation 11 and the contributing vertices associated to the facet f i,B are the two vertices incident to the edge e j,A .

f S = f i,B ⊕ e j,A (11)
• Intersection of two edge arcs: before trying to determine if two nonparallel edge arcs st and s t representing two edges e i,A and e j,B of A and B have a common intersection point, we must check if the contributing vertices associated to the facets incident to e i,A have distinct contributing vertices.

If this is not the case, we conclude that the corresponding edge arc st will not intersect any edge arc s t . Therefore, it will no longer be considered in the subsequent edge arcs intersection queries. In a similar manner, we must also check the contributing vertices associated to the facets incident to e j,B .

If there is at least one common contributing vertex, the edge arc s t will also be discarded from the rest of the iterations.

If edge arcs st and s t pass the previous step, we continue by checking if one of the two spherical points resulting from the intersection of the two great circles supporting st and s t lies on st and s t at the same time (see The corresponding Minkowski sum facet f S for this configuration is the result of a planar Minkowski sum of two edges e i,A and e j,B represented by the two edge arcs st and s t :

f S = e i,A ⊕ e j,B (12)
• Intersection of a facet point with a vertex region: from

Implementation and performance

In this section, we describe the implementations of the CVMS algorithm, the ESSDMS algorithm, and the Convex Hull-based Minkowski Sum (CHMS) algorithm. After that, we compare the performance of these approaches with other methods found in the literature and present examples of Minkowski sum polyhedra we computed.

Implementation

The CVMS algorithm has been implemented using C++ and CGAL. We used the same environment for the implementation of the ESSDMS and the CHMS algorithms. The convex polyhedra are handled by the CGAL Polyhedron 3 data structure.

The CHMS algorithm is implemented by using the CGAL function convex hull 3

for the computation of convex hulls in 3D. This function is an implementation of the Quickhull algorithm [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF].

For the computation of the planar Minkowski sums involved in the CVMS and the ESSDMS algorithms, we used the function minkowski sum 2 provided in the 2D Minkowski sum package of CGAL. This function implements the convolution operation [START_REF] Guibas | A kinetic framework for computational geometry[END_REF].

As stated previously, we guarantee the exactness of the obtained results by using exact number types. For our implementation, we used the exact number types provided by the GNU Multi Precision (GMP) library [START_REF]The GNU MP bignum library[END_REF] under CGAL.

By using exact number types, we penalize run-time performance. In fact, exact number types are very slow compared to the built-in floating point number types. To overcome this problem, we used the lazy kernel adapter [START_REF] Fabri | A generic lazy evaluation scheme for exact geometric computations[END_REF] which speeds up exact computations.

We also re-implemented the improved 3D-MSSD algorithm of Wu et al. using the same components of CGAL. Our implementation of the improved 3D-MSSD algorithm guarantees exact results because it uses exact number types, so it does not rely on position thresholds as done in the work of Wu at al..

Complexity study

In the case of the CVMS algorithm, finding the contributing vertices for the

f A facets of A takes O(f A v B )
(f A v B + f B v A + f A + f B + e A e B ).
(f A v B + f B v A )
. By summing all these times, we conclude that the ESSDMS algorithm has a worst-case time complexity of

O(f A f B + f A e B + f B e A + e A e B + f A v B + f B v A ).
For the Convex Hull-based Minkowski Sum (CHMS) algorithm we implemented, we performed vector addition of all points of A and all points of B in order to build a point cloud C. Then, we computed the convex hull of the point cloud C using the Quickhull algorithm. The first step takes O(v A v B ) time and the second has a worst-case complexity of O((f A + f B ) 2 ) [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF]. Thus, the CHMS algorithm has a worst-case time complexity of O((f A +f B ) 2 +v A v B ).

Performance benchmark

For our benchmark, we computed the Minkowski sum of several polyhedra using the CVMS algorithm, the ESSDMS algorithm, the CHMS algorithm, and the improved 3D-MSSD algorithm we re-implemented.

We also compared the above mentioned algorithms with three other algorithms. The first one is based on Nef polyhedra embedded on the sphere and implemented by Hachenberger [START_REF] Hachenberger | Boolean operations on 3d selective nef complexes: data structure, algorithms, optimized implementation and experiments[END_REF], the second one is a method based on linear programming implemented by Weibel [START_REF] Weibel | Minkowski sums[END_REF] (following the work of Fukuda [START_REF] Fukuda | From the zonotope construction to the minkowski addition of convex polytopes[END_REF]), and the third one is the Cubical Gaussian Map-based (CGM-based) algorithm of Fogel and Halperin [START_REF] Fogel | Exact and efficient construction of Minkowski sums of convex polyhedra with applications[END_REF].

We shall note that the method based on Nef polyhedra is more powerful than necessary for the Minkowski sum computation since it aims at overlaying two arbitrary Nef polyhedra embedded on the sphere (with lower dimensional features, unbounded or bounded boundaries, etc.). A second remark concerns the implementation of Weibel which is intended for the computation of polytopes (convex polyhedra in R n ) in an arbitrary dimension, so it is not specifically optimized for convex polyhedra in R 3 .

The experiments where done on a 2 GB RAM, 2.2 GHZ Intel Core 2 Duo personal computer. We used a cube and a sphere with varying facets number as input polyhedra. A comparison between running times is given in table 1:

term in the ESSDMS time complexity). This proves that the CVMS algorithm performs better than the ESSDMS algorithm and consolidates experimental results. Table 1 shows also that our CVMS and ESSDMS algorithms are much more efficient than the CGM-based algorithm of Fogel and Halperin; the computation of the cubical Gaussian maps in their algorithm is the most time consuming step when compared to the overlay one (especially for polyhedra with high complexity).

The algorithm implemented by Weibel is the slowest among all these algorithms since it is not optimized in R 3 . Another observation that is verified in table 1 is that the overlay process of the Nef polyhedra-based algorithm is not specifically devoted to the Minkowski sum computation of convex polyhedra.

An interesting fact of table 1 is that a good convex hull-based algorithm can perform better than other algorithms. In the present case, the CHMS algorithm performs better than the improved 3D-MSSD algorithm of Wu et al..

Examples

Fig. 13 shows some examples of Minkowski sum polyhedra computed by the CVMS algorithm. In Fig. 14, we included more results to show that the CVMS algorithm works with complex or big size polyhedra (having thousands or tens of thousands of facets; see Fig. algorithm can be used as an efficient replacement for the convex hull-based algorithms in order to speed-up the computation of the pairwise Minkowski sums for non-convex polyhedra (along with convex-decomposition and union

Fig. 1 .

 1 Fig. 1. Minkowski sum as a sweep of two sets. (a) Two polygons A and B. (b) B is positioned on each point of A. (c) The Minkowski sum A ⊕ B is the union of the resulting translations of B on all points of A.

Fig. 2 .

 2 Fig. 2. The intersection of closed half-planes. (a) A finite number of closed halfplanes H - 1 , . . . , H - 5 defined by H 1 , . . . , H 5 . (b) The boundary of P is partitioned into one-dimensional faces (edges) e 1 , . . . , e 5 and zero-dimensional faces (vertices) v 1 , . . . , v 5 .

  Minkowski sum and provide a two-sided Hausdorff distance bound on the approximation. Recently, Lien [29] used a point-based representation instead of the widely used mesh-based representation to compute the Minkowski sum of polyhedra. He uniformly sampled two point sets from the boundaries of two polyhedra, constructed a point set by adding all points from the two point sets already generated, and used three filters (a collision detection, normal, and octree filter) to discard interior points that are not on the boundary of the sum polyhedron. He demonstrated that the point-based representation provides the same functionality as the mesh-based representation by using it in several applications such as motion planing.

  The polyhedra are represented in a dual space called the slope diagram. Computing the intersections between the two merged or overlaid slope diagrams gives the sum polyhedron. Although Ghosh stated that slope diagrams can be used for both convex and non-convex polyhedra, the existent implementations of slope diagram-based Minkowski sum computation are only devoted to convex polyhedra. Some researchers proposed other variants of slope diagram representation. Bekker and Roerdink [34] used a slightly different slope diagram representation that works with edges instead of facets, so their 2D representation of the slope diagram reduces the problem's dimension. Wu et al. [1] presented some improvements to the slope diagram by merging the two diagrams without an explicit embedding on the unit sphere. They also avoided the use of stereographic projections (needed in the original slope diagram algorithm proposed by Ghosh).

  equation states that the Minkowski sum set A ⊕ B is obtained by sweeping all points of A by B and taking the union of all resulting points. Since we are working with polyhedra, it is clear that it is sufficient to only sweep the boundary of A by B and discard the interior points of A and B in order to have the Minkowski sum polyhedron S = A ⊕ B. Thus, we are only interested in the boundary of the resulting union. The boundary of A is composed of facets, vertices, and edges. So, in order to obtain the Minkowski sum polyhedron S, we must sweep all facets of A by B, sweep all vertices of A by B (or place B on all vertices of A), and sweep all edges of A by B. Observe that because the Minkowski sum is commutative, sweeping the vertices of A by B (or positioning B on vertices of A) is equivalent to sweeping the facets of B by A. Therefore, in the CVMS algorithm, we need to sweep all facets of A by B, sweep all facets of B by A, sweep all edges of A by B, and retain only facets that lie on the boundary of the union of these sweeps (facets lying on the boundary of the Minkowski sum polyhedron S).

3. 2

 2 From the sweep to the definition of contributing vertices To give an idea of what is a contributing vertex and how it relates to the sweep operation, we consider two convex polygons A and B depicted in Fig. 3.a, their Minkowski sum A ⊕ B is also shown in Fig. 3.b. We take a 2D example for simplicity purposes and without loss of generality. Let's consider the face (edge) f 1,A of A and sweep it by B. Sweeping f 1,A by B resumes to parallel-translating B such that its origin passes through all points of f 1,A . During this sweep, all vertices of B (i.e. v 1,B , v 2,B , and v 3,B ) are also moving in lines parallel to the supporting line of f 1,A . The displacement of these vertices produces faces (edges) having supporting lines parallel to that of f 1,A (faces that are drawn as dashed and dotted lines in Fig. 3.c). But only the facet generated by the vertex v 2,B (the dotted line in Fig. 3.c) which is at maximal distance away from the supporting line of f 1,A is considered because it lies on the boundary of A ⊕ B. The other produced faces (drawn as dashed lines in Fig. 3.c) are simply discarded because they do not lie on the boundary of A⊕B. This particular vertex v 2,B which generated this facet of the Minkowski sum polygon A ⊕ B is called the "contributing vertex" associated to the face f 1,A .

Fig. 3 .

 3 Fig. 3. (a) Two convex polygons A and B. (b) The Minkowski sum A ⊕ B depicted as a sweep. (c) The contributing vertex v 2,B associated to the face f 1,A . The dotted line is the face that lies on the boundary of A ⊕ B. The two other faces (dashed lines) are discarded.

  steps mentioned above. The first step is to sweep all facets of A by B. From what precedes, it is clear that the computation of the contributing vertices associated to the facets of A is a mean to obtain the facets of S = A ⊕ B that result from sweeping the facets of A by B, these facets of S are called "translated facets". The second step is to sweep all the facets of B by A. Therefore, computing the contributing vertices associated to the facets of B is a mean to obtain the facets of S = A ⊕ B that result from sweeping the facets of B by A. these facets are denoted "corner facets" of S.

Algorithm 1 3 : 10 : 11 : 12 : 3 . 4

 1310111234 An outline of the The CVMS algorithm for convex polyhedra Require: two convex polyhedra A and B Ensure: the Minkowski sum polyhedron S = A ⊕ B 1: for each facet f i,A of A do 2: compute its associated contributing vertices (w.r.t. B) deduce the translated facet corresponding to it 4: end for 5: for each facet f i,B of B do 6: compute its associated contributing vertices (w.r.t. A) 7: deduce the corner facet corresponding to it 8: end for 9: for each edge e j,A of A incident to facets having distinct contributing vertices do for each edge e k,B of B incident to facets having distinct contributing vertices do if e j,A and e k,B satisfy both visibility and normal orientation criteria then construct the edge facet e j,A ⊕ e k,B Contributing vertices applied to slope diagram-based algorithms We also show how the contributing vertices concept benefits slope diagram algorithms. So, we present the Enhanced and Simplified Slope Diagrambased Minkowski sum (ESSDMS) algorithm, which is somewhat similar to the improved slope diagram algorithm (improved 3D-MSSD algorithm) pre-sented in the work of Wu et al. [1]. This similarity is due to the fact that we do not embed input polyhedra on the unit sphere and do not use stereographic projections. The merging of slope diagrams and the handling of the intersections are made possible by simple geometric operations on the outer normal orientation information extracted from the polyhedra. The contributing vertices concept allows to further simplify the algorithm and to gain considerable performance by eliminating point in polygon queries and reducing edge arcs intersection queries. The computation of the contributing vertices does not introduce additional performance overhead since they are deduced directly from the intersection configuration. Details on our ESSDMS algorithm can be found in section 5.

  Since the boundary of polyhedra is composed of facets, edges, and vertices, sweeping the boundary of A by B is equivalent to sweeping the facets of A by B, sweeping the edges of A by B, and sweeping the vertices of A by B (or positioning B on vertices of A). Remember that it is sufficient to only sweep

Fig. 4 .

 4 Fig. 4. Illustration of the contributing vertices concept for the facets of a convex polygon A to be added to the convex polygon B.

  6.d).(2) Positioning B on all vertices of A: this means positioning B so that its origin coincides with each vertex v k,A . It is clear that only some facets of B will lie on the boundary of the sum polyhedron S when B is positioned on a particular vertex of A. The remaining facets of B are discarded because they lie inside S. To find the facets of B that are part of S, we exploit the commutativity property of the Minkowski sum operation, i.e. A ⊕ B = B ⊕ A. If we sweep all facets of B by A, we observe that the generated facets are the same as the facets we want to find by positioning B on each vertex v k,A . So, positioning B on all v A vertices results in the f B corner facets with supporting planes parallel to those of the facets of B (see Fig. 6.e).

Fig. 5 .

 5 Fig. 5.b). This produces facets of S that are the result of the Minkowski sum of the edge e j,A and all edges of B. But not all these generated facets lie on the boundary of the sum polyhedron S. So, S is composed of at most e A e B edge facets (see Fig. 6.f). The determination of edge facets is explained in section 4.3.

  translated facet corresponding to f i,A Contributing vertex associated to f i,A edge e j,A swept by B edge e l,B edge facet created by B

Fig. 5 .

 5 Fig. 5. (a) The creation of a translated facet corresponding to a particular facet f i,A and its associated contributing vertex (the corner facets are determined in the same manner). (b) Sweeping an edge e j,A by B and creation of an edge facet e j,A ⊕ e l,B .

Finally

  , it remains to sweep the edges of A by B. The result of this sweep are facets that are the Minkowski sum of the edges of A and the edges of B. But since we are working with polyhedra (boundary representation), we will not keep all these e A e B facets, only those lying on the boundary of the Minkowski sum polyhedron S are retained. Therefore, the Minkowski sum polyhedron is composed of at most e A e B edge facets.

Fig. 6 .Property 1 :

 61 Fig. 6. (a) Polyhedron A. (b) Polyhedron B. (c) Sum polyhedron S (d) Translated facets. (e) Corner facets. (f) Edge facets.

7.a is incident to the facets f 1 ,

 1 A and f 2,A . The vertices v 1,A and v 2,A are those vertices incident to e 1,A . The facets f 1,A and f 2,A have a common contributing vertex v 1,B (see Fig. 7.b). Therefore, the translated facets f 1,S and f 2,S corresponding to f 1,A and f 2,A are incident to the same edge of S delimited by two vertices with positions v 1,A + v 1,B and v 2,A + v 1,B (see Fig.

  7.c). In contrast, the edges e 2,A and e 3,A in Fig.7.a are incident to facets having distinct contributing vertices, so they contribute to the edge facets creation (Fig.7.c).Property 2: A facet of each polyhedron A or B contributes only once in the Minkowski sum polyhedron S; this has been proved in earlier works[START_REF] Kaul | Solid-interpolating deformations: construction and animation of PIPs[END_REF][START_REF] Ghosh | A unified computational framework for Minkowski operations[END_REF][START_REF] Bekker | An efficient algorithm to calculate the Minkowski sum of convex 3d polyhedra[END_REF]. This property allows to reduce progressively the number of facets considered in the subsequent steps of the algorithm. For example, if a facet contributed to the creation of a translation facet of S, it can be omitted from the corner facets determination step.The two mentioned properties enable a significant gain of performance by reducing the number of facets and edges to be considered when computing translated, corner, and edge facets.

Fig. 7 .

 7 Fig. 7. (a) Some edges e 1,A , e 2,A , and e 3,A of A. (b) Contributing vertices associated to f 1,A , f 2,A , and f 3,A . (c) Edge facets created by e 2,A , and e 3,A . The edge e 1,A will not produce edge facets since incident facets f 1,A and f 2,A have the same contributing vertex.

( 3 ) 4 )•

 34 two contributing vertices associated to f i,B lies on a supporting line parallel to the supporting plane of f i,B . The resulting corner facet is computed by a planar Minkowski sum of f i,B and e j,A (f i,B ⊕ e j,A ) and has the same outer normal orientation as the facet f i,B . • If there is only one contributing vertex v k,A , the corresponding corner facet results from translating the facet f i,B by a position vector pointing from coordinates origin o to v k,A and has the same outer normal orientation as the facet f i,B . Edge facets determination: we use two criteria for this purpose: the visibility of facets of B with respect to a visibility direction (defined by an edge e i,A ; see Fig. 8.a) allows to find horizon edges of B. Horizon edges are edges that separate invisible facets from other facets B with respect to the considered visibility direction. By other facets of B, we mean facets that are either visible or having supporting planes parallel to the visibility direction e i,A (see Fig. 8.b). All these horizon edges are candidates to construct an edge facet of S, but only those verifying the second criterion of normal orientation will do (see Fig. 8.c). Thus, for each edge e i,A incident to facets having distinct contributing vertices: • Determine the visibility of all facets of B. A facet f j,B with an outer normal n j,B is invisible with respect to the visibility direction e i,A if and only if: e i,A , n j,B > 0 (The facets of B that do not satisfy the above inequality are either visible or lie on supporting planes parallel to the visibility direction. The horizon edges are frontier edges between invisible facets and other facets of B with respect to the visibility direction e i,A (see Fig. 8.b). For each horizon edge e k,B incident to facets having distinct contributing vertices: Add an edge facet e i,A ⊕ e k,B to the sum polyhedron S if its outer normal orientation n i,k = e i,A ×e k,B lies between the two outer normal orientations n 1,A and n 2,A of the facets f 1,A and f 2,A incident to the edge e i,A (n 2,A follows n 1,A in a counter-clockwise order; see Fig. 8.c).

Fig. 8 .

 8 Fig. 8. (a) An edge e i,A as a visibility direction . (b) Visibility computation of all facets of B w.r.t. e i,A and horizon edges determination. (c) Only horizon edges e 4,B , e 5,B , and e 6,B will produce edge facets e i,A ⊕ e 4,B , e i,A ⊕ e 5,B , and e i,A ⊕ e 6,B . hanced and Simplified Slope Diagrams-based Minkowski Sum (ESSDMS) algorithm. We first describe the principles behind the original slope diagram.

Fig. 9 .

 9 Fig.9. A convex polyhedron and its slope diagram representation (revised from[START_REF] Ghosh | A unified computational framework for Minkowski operations[END_REF]).

Fig. 10 .

 10 Fig.10. Stereographic projection used to merge two slope diagrams (revised from[START_REF] Ghosh | A unified computational framework for Minkowski operations[END_REF]).

6 )

 6 Intersection of two facet points: two facet points are coincident if the two outer normals a and b of the corresponding facets f i,A and f j,B of A and B respectively have the same orientation (see Fig. 11.a); this is the case when a and b satisfy: a × b = 0 and a, b > 0 (The first part of equation 6 (the cross product) checks if a and b are on the same supporting line (since they have the same origin point, they have the same supporting line if they are parallel). The second part of equation 6 (the dot product) checks if a and b have the same orientation since they can be on the same supporting line but with opposite orientations.

9 )

 9 Intersection of a facet point with an edge arc: a facet point f corresponding to a facet of one of the two slope diagrams lies on an edge arc st of the other slope diagram (s and t are facet points corresponding to the facets incident to the edge represented by the edge arc st) if it lies on the great circle passing through facet points s and t and if f lies on the shorter arc joining s and t (see Fig. 11.b). These two conditions are equivalent to: f × n = 0 where n = s × t (8) s × f, n > 0 and f × t, n > 0 (Equation 8 checks if the facet point f lies in the great circle passing through s and t, this is the case when f is perpendicular to the normal n of the great circle. Equation 9 checks if the facet point f lies on the shorter arc of the great circle joining s and t. The corresponding Minkowski sum facet f S for this configuration is the result of a planar Minkowski sum of the facet represented by the facet point f and the edge represented by the edge arc st. For example, if the facet point f corresponds to a facet f i,A of A and the edge arc on which f lies corresponds to an edge e j,B of B, the Minkowski sum facet f S of S is given by:

Fig. 11 .

 11 Fig. 11.c). Suppose n is the normal of the circle containing st, i.e. n = s × t and n the normal vector of the great circle containing s t , i.e. n = s × t . The two spherical points x and x resulting from the intersection of the two great circles are obtained by the cross products x = n × n and x = n × n (see Fig. 11.c). To check that x (or x ) lies on the edge arcs st and s t , we can consider it as a facet point and use equations 8 and 9 derived for the intersection of a facet point and an edge arc.

Fig. 11 .Fig. 11 .

 1111 Fig. 11. Intersections in slope diagrams. (a) Intersection of two facet points. (b) Intersection of a facet point and an edge arc. (c) Intersection of two edge arcs. (d) Intersection of a facet point and a vertex region.

  time. Similarly, finding the contributing vertices for the f B facets of B takes at most O(f B v A ) time. The determination of translated facets takes O(f A ) time since it is done directly by examining the number of contributing vertices associated to each facet. In the same manner, constructing corner facets takes at most O(f B ) time. For edge facets, it takes at most f (e A e B ) time to find them. By adding times of all steps, we conclude that the CVMS algorithm has a worst-case time complexity of O

For

  the ESSDMS algorithm, the computation of the intersections of two facet points requires O(f A f B ) time, the computation of the intersections of facet points and edge arcs requires O(f A e B + f B e A ) time, and the computation of the intersection of two edge arcs requires at most O(e A e B ) time. Since the computation of the intersection of a facet point and a vertex region reduces to finding the contributing vertices corresponding to the remaining facet points (that have not yet contributed into the sum), the complexity of finding such intersections is at most O

Fig. 13 .

 13 Fig 14.c shows the result of such an extreme case correctly computed by the CVMS algorithm. More results of the CVMS algorithm can be found at: http://liris.cnrs.fr/hichem.barki/mksum/CVMS-convex

Fig. 14 .

 14 Fig. 14. More results of the CVMS algorithm for: (a) two spheres having thousands and tens of thousands of facets. (b) an ellipsoid-like polyhedron and a truncated dodecahedron. (c) for an extreme degenerate case (two identical spheres with thousands of facets).

8 Acknowledgment

 8 of our future work, we are generalizing the CVMS algorithm to nonconvex polyhedra. The main idea is that the boundary of each non-convex polyhedron can be seen as the union of disjoint convex surface patches (sets of facets containing no reflex edges). So by treating polyhedra in a local manner (working on each convex patch) and by correctly handling reflex edges, it is possible to generate a superset of the Minkowski sum facets. To extract the boundary of the Minkowski sum, we plan to use arrangements and surface envelope computations for which CGAL provides powerful packages. This work has been partly supported by the Cluster 2 of the Rhône-Alpes region within the LIMA Project. We are grateful to the help of P. Hachenberger and E. Fogel who provided us the necessary material and assistance for the use of their respective Nef polyhedra-based and cubical Gaussian map-based algorithms.

n 3,A n 4,A n 5,A v 1,B f 2,A f 3,A f 4,A f 5,A f 1,A v 2,B v 3,B v 4,B v 5,B v 2,B v 3,B v 4,B v 5,B v 1,B

Fig. 12 shows the ratio of the running time for the algorithms mentioned in Sum of the facets number (fA + fB) Runtime (sec.)

CVMS ( 1) ESSDMS ( 2) CHMS ( 3) MSSD ( 4) NGM ( 5)

Fuk. ( 6) CGM [START_REF] Kim | Fast penetration depth estimation using rasterization hardware and hierarchical refinement[END_REF] (1)

(2)