D Aiger
email: d.aiger@esiee.fr

H Talbot
email: h.talbot@esiee.fr

The Phase Only Transform for unsupervised surface defect detection

We present a simple, fast, and effective method for detecting defects on textured surfaces. Our method is unsupervised and contains no learning stage or information on the texture being inspected. The method is based on the Phase Only Transform (PHOT) which correspond to the Discrete Fourier Transform (DFT), normalized by the magnitude. The PHOT removes any regularities, at arbitrary scales, from the image while preserving only irregular patterns considered to represent defects. The localization is obtained by the inverse transform followed by adaptive thresholding using a simple standard statistical method. The main computational requirement is thus to apply the DFT on the input image. The method is also easy to implement in a few lines of code. Despite its simplicity, the methods is shown to be effective and generic as tested on various inputs, requiring only one parameter for sensitivity. We provide theoretical justification based on a simple model and show results on various kinds of patterns. We also discuss some limitations.

Introduction

Vision-based inspection of surfaces has many real-world applications, for instance industrial wood, steel, ceramic and silicon wafers, fruits, aircraft surfaces and many more. It is in high demand in industry in order to replace the subjective and repetitive process of manual inspection. A comprehensive survey on recent developments in vision based surface inspection using image processing techniques, particularly those that are based on texture analysis methods, was proposed by Xie. [START_REF] Xie | A review of recent advances in surface defect detection using texture analysis techniques[END_REF] According to this work, one can divide the methods for surface defect detection into four categories, namely: statistical approaches, structural approaches, filter-based methods, and model-based approaches. A significant differentiating factor in visual inspection approaches is that of supervised classification versus novelty detection. For applications where both normal and defective samples can be easily obtained, supervised classification based approaches are usually favored. However, when defects are unpredictable and defective samples are unavailable, novelty detection is more desirable.

Texture is one of the most important characteristics in identifying defects or flaws. Much effort was invested in extracting useful texture features. [START_REF] Haralick | Statistical and structural approaches to texture[END_REF][START_REF] Reed | A review of recent texture segmentation and feature extraction techniques[END_REF][START_REF] Tuceryan | Texture analysis[END_REF] Statistical texture analysis methods investigate the spatial distribution of pixel values. In structural approaches, texture is characterized by primitives or texture elements, and the spatial arrangement of these primitives. [START_REF] Nevatia | Structural analysis of natural textures[END_REF] The goals of structural approaches are to extract texture primitives, and to model the spatial arrangement. Filter based approaches share a common characteristic, which is to apply filter banks on the image and compute the energy of the filter responses. These methods can be divided into spatial domain, frequency domain, and joint spatial/spatial-frequency domain techniques. Model based methods include, among many others, fractal models, [START_REF] Freeman | The fractal geometry of nature[END_REF] autoregressive models, [START_REF] Comer | Segmentation of textured images using a multiresolution gaussian autoregressive model[END_REF] random field models, [START_REF] Li | Markov random filed modeling in image analysis[END_REF] and the texem model. [START_REF] Xie | Texem: Texture exemplars for defect detection on random textured surfaces[END_REF] In a novelty detection task, the task is to identify whether an input pattern is an expected part of the data or unknown. As for defect detection, it involves assigning a normal or abnormal label to a pattern (e.g. a surface or a pixel). In contrast to supervised classification, novelty detection only needs the normal samples for training purposes and usually uses a distance measure and a threshold for decision making. Recently, Markou and Singh [START_REF] Markou | Novelty detection: A review part 1: Statistical approaches[END_REF][START_REF] Markou | Novelty detection: A review part 2: Neural network based approaches[END_REF] gave a detailed review of novelty detection approaches, using statistical and neural network based approaches. Statistical parametric approaches are commonly used in visual inspection. [START_REF] Millan | Detection of local defects in textile webs using gabor filters[END_REF][START_REF] Kumar | Fabric defect segmentation using multichannel blob detectors[END_REF][START_REF] Bennamoun | Optimal gabor filters for textile flaw detection[END_REF][START_REF] Mirmehdi | Restructured eigenfilter matching for novelty detection in random textures[END_REF] A fundamental assumption is that the data distribution is Gaussian in nature, thus, it can be easily statistically modeled by means and covariances.

Working in the frequency domain is closely related to our contribution. Many methods apply filtering in the frequency domain, particularly when no straightforward kernel can be found in the spatial domain. The image is transformed into the Fourier domain, multiplied with the filter function and then re-transformed into the spatial domain. Coggins and Jain 16 used ring filters and orientation filters for feature extraction. D'Astous and Jernigan [START_REF] Jernigan | Texture discrimination based on detailed measures of the power spectrum[END_REF] used peak features, such as strength and area, and power dis-tribution features, such as power spectrum eigenvalues and circularity, to discriminate textures. Tsai and Heish 18 used the Fourier transform (FT) to reconstruct textile images for defect detection. Chan and Pang 19 extracted harmonic peaks from horizontal and vertical power spectrum slices. The phase of the DFT was used for matching images by the Phase Only Correlation (POC). [START_REF] Kuglin | The phase correlation image alignment method[END_REF] It was also used for image coding. [START_REF] Luo | Phase-only transform based shaping for error concealment coding[END_REF] The main focus of this chapter is to develop an unsupervised method for defect detection. We concentrated on a method that does not require any prior information or learning stage. This solves the practical problem of collecting usable examples of good patterns and defective ones. In addition, in many inspection applications, the pattern of the inspected surface is not known a-priori. As we work in the frequency domain, we have the advantages of a global view, solving the problems of selecting good kernel sizes. On the other hand, our method localizes defects without the need for any post processing. The idea is very simple: instead of trying to detect peaks in the Fourier transform of the input images, it simply removes all regularities in the image of various sizes and patterns at once by normalizing the FT of the input image by its magnitude. Through this operation, only the phase information remains while all regular patterns at all scales are removed. The localization of the defects is achieved by transforming back to the spatial domain. We show that since no analysis is being performed in the Fourier domain, removing regularities by normalizing the magnitude serves as a multiscale regularity removal, therefore, random textures are effectively removed as well.

The importance of the phase

In their important paper, Oppenheim and Lim investigated the importance of phase in signals. [START_REF] Oppenheim | The importance of phase in signals[END_REF] We review here some of their insights. In the Fourier representation of signals, spectral magnitude and phase tend to play different roles and in some situations, many of the important features of a signal are preserved even if only the phase is retained. This is not true in general for the spectral magnitude. This observation about phase has been made in a number of different contexts and applications. In general, reconstructing an object from the magnitude only is not of much value in representing the original object, whereas reconstructions from the phase only have many important features in common with the original objects. A phase-only image has Fourier transform phase equal to that of the original image and a FT magnitude of unity. of Lenna with magnitude only and phase only. It shows that although the rebuilt image contains the same frequency terms as the original input, the magnitude-only transform bears no resemblance to the original image, it does not allow recognition due to the missing phase. Conversely, a phaseonly reconstruction shows a picture similar to the shape of the original Lenna. It appears that the phase contains significant information, especially regarding the edge location. In addition, the phase only transform removes periodicity and regularity so it does more than just preserving the edges. In Section 1.3.2 we model the problem, give some theoretical justifications and explain why it works.

The new method 1.3.1. Applying the Phase Only Transform

In the context of this chapter, our goal is not to reconstruct a signal from its FT. In fact, we wish to do almost the opposite, namely, to filter out most patterns. We wish to eliminate some parts of the signal while preserving others. The above discussion provides a way to do it for our application of defect detection in images. Our purpose is to eliminate the so-called regular patterns while preserving correspondingly "rare" events in the image which can be considered to be defects. The fact that phase-only reconstruction preserves much of the correlation between signals would suggest that the location of events tends to be preserved. If we assume that in our application, a non defective region contains non-localized structures (e.g. regular patterns or homogeneous regions) and that a defect is well localized, it is reasonable to use the phase to filter all non localized patterns (see Section 1.3.2).

Our algorithm is very simple and can be implemented in a few lines of code. The first step is to apply the DFT on the input (real) image.

The discrete two-dimensional Fourier transform of an image array F (u, v) is defined in series form as:

F(u, v) = 1 N N -1 j=0 N -1 k=0 F (j, k) exp{ -2πi N (uj + vk)} where i = √ -1.
The indices (u, v) are called the spatial frequencies of the transformation. The result is a matrix of complex numbers in the frequency domain,

F(u, v) = R(u, v) + iI(u, v)
or in magnitude and phase-angle form,

F(u, v) = M(u, v) exp{iΦ(u, v)} where M(u, v) = R 2 (u, v) + I 2 (u, v) and Φ(u, v) = arctan{ I(u, v) R(u, v) }
By normalizing every complex number by dividing both the real and imaginary parts by M(u, v) we essentially remove all regular patterns at every scales at once. Note that we don't have to analyze the Fourier image. This normalization works on all frequencies at once and eliminates the regularities. The algorithm of the first stage can now be written (algorithm 1.1):

Algorithm 1.1 The Phase Only Transform Require: input image I(u, v) compute F(u, v){I} for all (u,v) do F(u, v) = F (u,v) M(u,v) end for O(u, v) = F -1 (u, v)
The resulted image is O(u, v). The spectral magnitude of images tends to fall off at high frequencies, the phase-only signal will experience a highfrequency emphasis which will accentuate narrow events without modifying their position. It is reasonable to identify 1 M(u,v) as generally emphasizing high frequencies over low frequencies.

Theoretical justification

We justify here why the Phase-Only Transform (PHOT) works. For this we express the fact that for texture, the range of value taken by the phase is small compared to that for a defect in the texture.

Phase of regular texture

We limit ourselves to 1-D, as the discussion extends readily to n-D due to the FT separability, and we carry out the discussion in the continuous domain for simplicity. We need some definitions to start: Definition 1 (projection). Let φ be a function defined over a 1-D domain D. The projection p of φ is the indicator function of the domain of φ(y) for all y.

p D [φ(x)](y) = 1 if there is an x over D such that φ(x) -y = 0 0 otherwise.

We note that if D is bounded and φ takes discrete values, i.e. if φ is zero almost everywhere, then p D [φ(x)] is also zero almost everywhere. Conversely, if φ is monotonic and non-constant, then it Lebesgue measure is non-zero.

Definition 2 (Phase integral excursion). The integral excursion E of the phase of a real signal is a measure of the range of values that the phase of its FT actually takes. More precisely, let S(x) be a real signal. Let

F[S](ω) = +∞ -∞ S(x)e -iωx dx = a S (ω)e iφS (ω)
be its FT with ω real belonging to [0, 2π]. Assuming φ S (x), the phase of S, to be either monotonic or discrete on that domain, we define

E[S] = 2π 0 p D [φ(x)](y)dy. (1.1)
Our phase integral excursion is a measure of the density of the values of φ S (ω) projected onto the y-axis. We limit ourselves to phases that are monotonic or discrete as this is sufficient for our discussion, and as in these cases, we can define p D implicitly. To continue, we need a simple definition of texture.

Definition 3 (regular texture). We defined a regular texture as a signal representable by a convergent Fourier sum:

S(x) = +N n=-N c n exp(inx), N ≥ 0. (1.2)
We know that Fourier sums can represent any square-integrable bounded periodic signal with arbitrary precision almost everywhere, which is what we require to represent regular textures. We now have then the following theorem:

Theorem 1 (Phase integral excursion of a regular texture). The phase integral excursion of any regular texture over a finite domain is zero.

Proof:

The FT of S reduces to a superposition of a finite number of Dirac peaks, therefore φ(x ∈ D) is zero almost everywhere, and so is its projection. Finally, so is its phase integral excursion, as the integral of a function with only countably non-zero values is zero.

Phase of a defect

A defect in a texture can be defined as an abrupt change in its regularity.

To study this we model it by a box function.

Definition 4 (Random box function). A random box function H(a, b)

is a function which has the following form:

H(a, b)[x] =    0 if x < a 1 if a ≤ x ≤ b 0 if x > b
We use the following well-known properties of the Fourier transform to derive the expression of the FT of H(a, b).

• Translation invariance : F[f (x -x 0 )](ω) = e -ix0ω F[f (x)](ω). • Scale invariance : F[f (ax)](ω) = 1 |a| F[f (x)]( ω a ). • The expression of the centered box function : H(-1 2 , 1 2 ) = sinc( ω 2 )
The expression of the FT of the random box function is therefore :

F[H(a, b)](ω) = e -i( a+b 2 )ω sinc( ω 2(b -a)
) .

(1.

3)

The phase of this FT is simply is monotonic and non constant, and so, even including phase wraparound, p D [φ(x)] is 1 on a measurable set. Its integral over the range of φ is therefore non-zero, and so is the integral phase excursion. We note that since a and b are random, the probability of a + b = 0 is zero.

φ H(a,b) [F[H(a, b)]](ω) = -( a + b 2 )ω. ( 1
Let us assume a regular texture on the one hand, and a regular texture with a defect in the other. Theorem 1 tells us that the former has a FT with a phase composed of only a few different values. The latter might be viewed as a superposition of a regular texture and a random box function with random values for a and b. Theorem 2 tells us that its FT features a phase composed of uncountably infinite different values.

We now show that the phase-only transform can readily distinguish between these two cases even in the discrete setting.

1D examples

In this section we show a few examples on 1D signals and give some insights about the behavior of the Phase Only Transform. We refer to the PHOT here, as the signal that is transformed back to the spatial domain, after being normalized by the magnitude. As already shown by the 2D example, most of the information on edges and sharp peaks is contained in the phase. If a signal contains a single peak or edge and a flat region, the phase part of the FFT must be significant, because the sum of many trigonometric functions is needed to construct the flat part. On the other hand, if a signal is constructed of a sum of pure sine or cosine functions of various frequencies with zero or little phase content, the PHOT will be almost zero. This is true not only for signals that are periodic within a finite support. Figure 1.2 shows such a signal. In Figure 1.3 we see a sharp peak that requires large phase content. We conclude that signals (not necessarily periodic) that have a small phase content would yield a smooth PHOT, while those with large phase content representing a peak or an edge yield a large peak in the PHOT which corresponds to the location of the peak or edge in the input signal. Assuming that a defected signal is composed of sum of sine function of various frequencies and a peak, the result of the PHOT is a collection of peaks in the spatial domain that are localized in the original defect location while the part that is corresponding to the first term is eliminated. Figure 1.4 shows a small defect (peak), composed with a sine (or cosine) wave. In Figure 1.5 we show another example on a signal that appears non-periodic due to the limited domain, yet, is composed from a sum of trigonometric functions which are all removed, while the defect remains.

Our model of an input signal is thus composed of two terms, a nondefected term, A(x) which is a sum of sine or cosine functions with relatively small phase content, and a defect term, B(x) which is assumed to be a peak or step edge, thus contains large phase content:

S(x) = A(x) + B(x)
Since the PHOT eliminates the sum of (low phase content) sines, we are left mainly with B(x), as expected from section 1.3.2. The inverse transform then yields the localization of the defect in the spatial domain.

Thresholding using Mahalanobis distance

In order to be able to use a totally unsupervised method with no learning component, we have to assume that for each input image the majority of the image pixels are intact (see Section 1.4). In this case, we can use simple statistics. We use the result of the PHOT as a probability map of a pixel being a defect. As commonly used, we assume a Gaussian distribution and use the Mahalanobis distance. We compute the mean and variance of the distribution from the image obtained by the PHOT. Since we normalize each of the FFT basis when we reconstruct the PHOT image, the global mean and standard deviation of the image are now both 1/N where N is the number of pixels. However, since the noise can be significant, we first smooth the PHOT image by a Gaussian filter and only then compute the mean and variance (we have used Gaussian of σ = 3.0). The user provides D. Aiger and H. Talbot a value in sense of Mahalanobis distance. We threshold the PHOT result such that every pixel with a distance larger than this value is considered as a defected pixel. Figure 1.6 shows an input image, the PHOT result interpreted as Mahalanobis distance from the mean and the thresholding result using a Mahalanobis distance of 4.0. Of course more sophisticated statistical methods can be used instead.

Characteristics and limitations of the Phase Only Transform

The most appealing characteristic of the PHOT is that it removes any regularities from the image without the need to identify peaks in the Fourier domain. Only spikes that do not correspond to a sum of trigonometric functions inside the image domain are left. Note that the regularities should not be presented in the entire image. Every large enough regular patterns are removed by the transform by normalizing the resulted complex number by its magnitude. In this sense, our method is different from those that work only on periodic patterns. Figure 1.7 shows an example of image that has several subpatterns that are regular but the entire image is not. The only parameter in the threshold on the Mahalanobis distance and it is exactly the same in Figures 1.6 and 1.7. The result shows that the PHOT has no difficulty in detecting defects in this image. The results look very similar to the human perception of "novel pattern". The entire image is not regular but contains patterns that in some way similarly perceived. We should note here that this can be also considered as a limitation of the method, since large defects can be viewed as regular subpattern, thus might be removed by the PHOT.

As can be expected, if we use 2D FFT on the image, every periodicity or regularity (or homogeneity) is removed by the PHOT. This contains also large defected patterns and 1D structures. For example, a defect structured as a line or scratch in the image, would not be well detected as can be seen in Figure 1.8. On the other hand, the same characteristic, can be used (to our advantage) to obtain defect detection on multiple patterns where nothing has to be known by the algorithm in advance ("blind" defect detection). In Figure 1.9 the results of our algorithm on a image that contains two totally different regularities are shown. It can be observed that the boundary regions between regularities were removed by the PHOT. This means that 1D long defected patterns may not be detected. A way to solve this problem is to apply the PHOT on lines instead on the entire 2D image. This would work however only in a highly regular patterns. We will investigate this direction in the future.

Complexity and real time performance

In many inspection system that apply defect detection algorithms for quality assurance, the time performance of the algorithm is critical as it might be used in a real manufacturing process. As can be easily concluded from our algorithm, the complexity is O(n log n) where n is the number of pixels in the input image. This, of course, comes from the DFT that we have to apply. The further processing and statistics is obviously linear with n. For very large or continuously inspected patterns, one can apply the algorithm on partial sub-windows without affecting the detection performance substantially . It is also very simple to implement the algorithm on parallel machines by decomposing the input. We successfully implemented the algorithm on a GPU (Graphics Processing Unit) using the Nvidia CUDA language. The FFT is also quite fast in practice and effective parallelization exists using Intel's SSE2 and SSE3 instructions, as well as on DSPs.

Results

We implemented the algorithm using C++ and Visual Studio. the results on a large set of images are shown in Figures 1.10. All the results were obtained using the same parameter for thresholding the Mahalanobis distance (4.0). No other parameter is needed for our algorithm. The sensitivity of the algorithm can be changed by the user by altering the Mahalanobis threshold.

Multiple sub-patterns and arbitrary patterns

As already mentioned in Section 1.4, our method does not require that the entire inspected pattern be regular. It can process many sub-patterns simultaneously. In fact, the PHOT is a detector for novel patterns. It emphasizes patterns that do not appear much in the image. It is worth noting that we do not assume anything about the size of the pattern, so it can vary. In Figure 1.11 an image containing many texture patches of different size and regularities is proceeded and the result (using Mahalanobis threshold 4.0) is shown on the right. The synthetic defect almost invisible by eye in the image is detected since it is novel. Another spike on top of the image is also detected. In Figure 1.12 a scene that contains a house with a textured roof is shown. The image contains textures as well as homogeneous and irregular regions. The synthetic defect as well as the novel pattern of the lamp on the right are well detected.

Images with no defects

We tested our simple adaptive threshold on input images which are texture patches without any defect. The purpose of this test is to verify that the method does not produce false positives. We used exactly the same parameter as in all other tests, namely, a Mahalanobis distance = 4.0. In Figure 1.13 we show two texture patches which are not quite regular (to make the test more difficult), their PHOT results and the output using threshold equal to 4.0. It can be seen that no false positive defects were produced for either inputs. It can be observed in the PHOT result (middle), how the strength of the response is related to the perception of "novelty". Although no pixel exceeds distance 4.0, some regions have larger response correlated to the measure of their regularity.

Other potential applications

The main application of the Phase Only Transform presented in this chapter is defect detection, however, as the PHOT detects novel patterns in an image it can be also used for other applications. Salient regions are generally regarded as the candidates of attention focus in human eyes, which is the key stage in object detection. The phase spectrum plays a key role for saliency detection. [START_REF] Guo | Spatio-temporal saliency detection using phase spectrum of quaternion fourier transform[END_REF] The saliency map can be calculated by the image's Phase spectrum of Fourier Transform alone. It was shown, similarly to the analysis in this chapter, that phase information specifies where each of the sinusoidal components resides within the image. The locations with less periodicity or less homogeneity in either the vertical or horizontal orientation show were the object candidates are located. In, [START_REF] Guo | Spatio-temporal saliency detection using phase spectrum of quaternion fourier transform[END_REF] each pixel of the image is represented by a quaternion that consists of color, intensity and motion feature. The Phase spectrum is then used to obtain the spatio-temporal saliency map, which considers not only salient spatial features like color, orientation and etc. in a single frame but also temporal feature between frames like motion. Two examples from 23 is shown in Figure 1. 14 Another possible application is to measure the amount of "rectification" in images containing repeated patterns (like textures) that were taken in perspective. This can subsequently allow a rectification algorithm that maximizes this measure. For example, in Figure 1.15, the left image contains more homogeneity than the unrectified image to the right, thus a measure that is based on, say, the integration of the PHOT of the image would be much larger for the image to the right. Minimizing this measure (maximizing the homogeneity) would achieve rectification. The effect of repeated pattern on the PHOT is clearly observed.

Conclusions

A novel method for defect detection on surface patches was presented. The main advantage of the new algorithm is its extreme simplicity (it consists manly of a standard forward and inverse FFT), its generality to work for various pattern without prior knowledge and the fact that it is unsupervised. We gave theoretical justification for a reasonable model. We show results on a large set of inputs and the results are very similar to the perception of defects where no prior information is given. The new algorithm has only one parameter which is the sensitivity of the algorithm. It is an advantage in real inspection systems, where ease of use is important. The algorithm is also fast in practice and can be used in real time systems. Moreover, parallelization of the algorithm can be easily obtained by simply subdividing the input.

Figure 1 .

 1 Fig. 1.1. The Phase only inverse Vs. Magnitude Only inverse -left: original, middle: magnitude only, right: phase only

. 4 )

 4 We now have the following theorem : Theorem 2 (Phase excursion of the random box function). The phase excursion of the random box function is almost surely non-zero. Proof: Ignoring phase wraparound over 2π, If a + b = 0, then φ H(a,b)

Fig. 1 .

 1 Fig. 1.2. A signal with little phase content -the PHOT is almost flat. Top: signal and its PHOT. Middle: magnitude of frequencies. Bottom: phase of frequencies.

Fig. 1 .

 1 Fig. 1.3. A small defect in a sum of sine curves. Top: signal and its PHOT. Middle: magnitude of frequencies. Bottom: phase of frequencies.

Fig. 1 .

 1 Fig. 1.4. A defect in a single sine curve. Top: signal and its PHOT. Middle: magnitude of frequencies. Bottom: phase of frequencies.

Fig. 1 . 5 .

 15 Fig. 1.5. Non periodic (but with little phase) signal and a defect (large phase content). Top: signal and its PHOT. Middle: magnitude of frequencies. Bottom: phase of frequencies.

Fig. 1 . 6 .Fig. 1 . 7 .

 1617 Fig. 1.6. Image in the process of defect detection: left -input image, middle -Mahalanobis distance from the mean (multiply by 30 for visualization), right -thresholding using distance 4.0

Fig. 1 . 8 .

 18 Fig. 1.8. Limitation of 2D transform: Scratches could not be detected as they are 1D regular.

Fig. 1 .

 1 Fig. 1.10. Results on various patterns: in each of the three columns, left -input image, right -results by thresholding using distance 4.0

Fig. 1 .Fig. 1 .

 11 Fig. 1.11. Multiple textures of various size and regularities and a synthetic defect: top -input image, bottom -thresholding using distance 4.0

Fig. 1 .Fig. 1 .

 11 Fig. 1.13. Images with no defects: left -input image, middle -PHOT result (multiplied by 30 for visualization), right -result using threshold of 4.0


Michel Couprie

Chapter 10

Distance, granulometry, skeleton

In this chapter, we present a series of concepts and operators based on the notion of distance. As often with mathematical morphology, there exists more than one way to present ideas, that are simultaneously equivalent and complementary. Here, our problem is to present methods to characterize sets of points based on metric, geometry and topology considerations. An important concept is that of the skeleton, which is of fundamental importance in pattern recognition, and has many practical applications.

Skeletons

The notion of skeleton, which we will define shortly in more than one way, is very useful in many applications. Intuitively, the skeleton of an object (of a set of points) X is a set of "median" lines, constituted of points that are located at an equal distance from distinct areas of the border of X. Historically, a definition of the general concept was proposed by Listing in 1861 [LIS 61], under the name of "cyclomatic diagram". This diagram results from the contraction of a closed curve in the Euclidean plane, under topology preservation constraints. This notion is essentially a homology concept.

A more precise and geometric definition, that of the median axis was proposed by Blum almost exactly one hundred years later [BLU 61, BLU 67]. This definition requires the notion of maximal ball.

Figure 10.1. Illustration of a maximal ball in the Euclidean case. The Ball Bmax is maximal because no other ball can contain it and simultaneously be included in the set. In contrast, the ball B2 is not maximal.

Maximal balls

Maximal balls are a simple geometric concept that will allow us to introduce the notion of skeleton of a set of points. DEFINITION 10.1 ( Ball).-The ball B r (x) of radius r and center x is the set of points y of R n such that d(x, y) < r, where d is a distance (e.g. the Euclidean distance or the city-block distance) DEFINITION 10.2 ( Maximal ball).-A maximal ball B max r of radius r for set X is a ball included in X such that ∀s > r, if B max r ⊂ B s then B s ⊂ X.

This notion is illustrated on Fig 10 .1. We remark that on this figure, where the set under study exhibit a regular border, the maximal ball B max is tangent to the border of the set in two distinct points. This is common but not the general case: sometime there are more than two points, and even a non-countable number. The point C on this figure is the center of B max . It is clearly located at the same distance of the two points of B max that are tangent to the border of the set.

Firefronts

Almost simultaneously to the definition of Blum appear those of Calabi and Hartnett [CAL 68], as well as that of Montanari [MON 68], that use the notion of firefront. DEFINITION 10.3 ( Firefront and quench function).-Suppose a fire is lit on the border of X and propagates toward the interior of X at a constant rate, in an isotropic manner. Assuming that a burnt point does not lit itself again, then the locus of points of X where several firefronts meet are called the quench points. The function that at each quench point associates the time at which the point stopped burning is equivalent to the distance to the border of the set. This function is called the quench function.

This definition is useful conceptually, because it forms the basis of many algorithmic ideas, in particular regarding skeletonization methods that use the idea of "eating up" at the border of the object. PROPERTY 10.4.-For a sufficiently regular compact set X1 in the Euclidean space, points p are centers of maximal balls of X if and only if they are quench points.

It is not immediately obvious that the two notions are equivalent, but it is not very difficult to prove. Intuitively, the meeting point C of several firefronts is located at the same distance of n distinct points A 1 , A 2 , . . . , A n of the border of X, from where the firefronts originate, due to the constant propagation rate. Because of the isotropic propagation, the segments [CA i ] all have the same length, are normal to the border of X and are radii of a ball B C,Ai , centered in C. In addition, this ball is necessarily included in X, otherwise another distinct point from the border or X would have extinguished C before the firefronts coming from the points A I . Finally, any ball that would include B C,AI and included in X must also be normal to ∂X on all points A i . It is therefore necessarily centered in C and of radius [CA i ], therefore the ball B C,Ai is maximal.

Conversely, let B C a maximal ball of X, centered in C. Since it is maximal, it intersects ∂X on at least two points A and B. Indeed, if we suppose that B C does not intersect ∂X at any point, then by definition, we can include a circle strictly in X B C , and so B C is not maximal. This is still true if B C intersects ∂X on a single point A, assuming all the following curves are tangent in A: the border of X, that of B C and that of the circle. Finally, the two firefronts issuing from A and B meet in C, and so C is a quench point.

The definition of firefronts is useful as an illustration but not very rigorous, and linked to the continuous domain. It may be useful to consider the properties of the skeleton in this domain.

Properties of the skeleton in the continuum

With standard disclaimers, the preceding definitions apply to the continuous domain, in this case arbitrary dimension manifolds, equipped with a Riemannian metric.

Here we shall only describe the R n case with the Euclidean metric. We will consider as before the set X to be skeletonized as a compact set with a border everywhere twice differentiable. To avoid some limit cases, we shall consider that the interior of X, denoted X is connected.

In this context, the two definitions of the skeleton given earlier are equivalent, and exhibit the following properties: PROPERTY 10.5.-Properties of the skeleton in the continuous domain -Non-continuity: a continuous deformation of the set can induce non-continuous changes in the associated skeleton. For instance, the euclidean skeleton of a disk is a single point : its center. However, any infinitesimal change on the border of the disk will result in at least one branch in the skeleton.

-Homotopy: the skeleton is homotopic to the original set.

-Negligibility: the Lebesgue measure of the skeleton is zero. This translates into the property that the skeleton does not contain any neighborhood.

-Invariance by isometric transform: applying any isometric transform (i.e: a translation or rotation) before or after the computation of the skeleton leads to the same result.

-Links with the distance transform:

-The function that associates any point of the skeleton with its distance to the border of the set, also known as the quench function seen earlier can be used to reconstruct the original set perfectly. This is achieved simply by dilating each point of the skeleton with the ball of radius equal to the value of the quench function at that point.

-The skeleton is centered in the object, in the sense that each point of the skeleton is located at equal distance from at least two distinct points from the border of the object.

-For each point belonging to the skeleton, the line of greatest ascent for the distance transform issuing from that point is in the direction of the skeleton. We say the upstream of every skeleton point also belongs to the skeleton.

-For every point of the skeleton, the slope of the upstream is strictly less than 1.

-The centers of maximal balls are special points of the distance transform: they form the ridge point of this transform.

These properties are explored in detail in [MAT 88a,MAT 88c]. Many are demonstrated in [RIV 87]. The non-continuity of the skeleton translates in practice into a high sensitivity to noise. Just one point added to the border of an object can induce the formation of a new branch of the skeleton. A small arbitrary change of X may involve a change of topology and radical changes in the appearance of the skeleton. In general we try to regulate these effects by filtering either on the shape of X, or on the skeleton itself [VIN 91a,ATT 95].

The homotopy of the skeleton is not guaranteed by all discrete skeletonization algorithms. Most of them however, try to ensure this important property, which is discussed in Section 10.6. Preservation of the homotopy of the skeleton is a property used in pattern recognition.

The fact that the skeleton is negligible in the continuum is sometimes reflected in need to be thin, i.e., it can be erased by an erosion with the elementary ball for the connectivity of the grid used for computing it. Some algorithms attempt to produce a skeleton of thickness of one a pixel, which is not necessarily compatible with other properties of the skeleton, such as centering.

Skeletons in discrete spaces

We can not calculate, much less use in a practical way, the continuous skeleton, but we will seek, in the discrete domain, to reproduce all or part of its characteristics.

First, note that there exist skeletonization algorithms that start from the definition of the continuum, and formulate the problem in this area. For example, we can model a forest fire through an evolution of the contour by Partial Differential Equation (PDE), and find the skeleton points by detecting points of shocks, which are points in space where only a weak solution of the PDE exists, or maintaining the Center of Maximal Balls (CMB) as anchors for the evolution of the contour [LEY 91, SID 99, TOR 03]. Ultimately, despite its theoretical elegance, these methods are still a discrete skeleton due do the fact that the PDE used for the calculation are discretized. It is difficult to prove anything about the results because the solutions are only approximate. Finally things can become even more complicated in 3D and more.

We can not reproduce in a discrete setting all the properties of continuous skeletons, but we may instead choose a particular compromise as needed.

We will try to calculate, for a discrete object X, an object S(X) -homotopic to the set X;

-thin: that is to say, a minimum number of pixel or voxel large, or an object erasable the elementary erosion;

-which contains the centers of discrete maximal balls;

-centered in the object X, that is to say that every point of S(X) should be at equal distance, for the underlying grid, of two distinct points of the boundary of X.

The latter property is usually impossible to satisfy strictly. For example, for the 8-connectivity in the plane Z 2 , for a set X that consists of a 10×10 square, the centers of maximal balls, which are within this grid squares of size e n × n with n odd, are not located at equal distance from all points of the border of X.

In the following sections we present two approaches that allows to define the skeleton of a shape in a discrete space. The first approach, described in Section 10.3, relies on the granulometry concept, familiar in mathematical morphology. It has the advantage of inducing interesting generalizations based on such family of structuring elements non-symmetrical or not connected. The second approach involves using discrete distances to define the notions of ball and maximal ball relatively to a given object (Section 10.4). This latter approach allows the use of the exact Euclidean distance, which provides for the skeleton less sensitivity to the orientation of the object in space, a full rotation invariance being impossible to guarantee.

Granulometric families and skeletons

Granulometric family

We first start by defining a way to characterize the size of a set of points. We shall see later that this characterization is rich in consequences.

DEFINITION 10.6.-We first define the concept of a granulometric family using the following steps:

1) we define a family of erosions {α i } called elementary by convex structuring elements indexed by an integer i;

2) we construct the family of structuring elements 3) these structuring elements define a granulometric family {κ i }; 4) in most cases we use a homogeneous family, that is to say when all elementary erosions are identical: ∀i,

κ i = α i . . . α 2 α 1 ; (a) (b) (c) (d)
α i = α 1 , implying κ i = (α 1 ) i .
Heterogeneous families can be used to obtain a finer control over the granulometry, in particular, they allow the use of more isotropic structuring elements. These families are illustrated in Figure 10.3. DEFINITION 10.7.-We have the following property, by construction: ∀j ≤ i, γ κi • γ κj = γ κi . We say that each κ i is open by all κ j for j ≤ i. This property is called "absorption property".

Applications of granulometries

Granulometric families allows to define incremental operations with desirable properties.

-We can trivially extend the definition to grayscale.

-For example, binary or greyscale granulometries can be used to characterize texture or to obtain information about the image without segmentation.

-We give some examples of applications of granulometries in Chapter 1 and in Chapters 13 and 19.

Ultimate eroded formula

The ultimate eroded consists of specific points of the skeleton, and is defined as follows:

DEFINITION 10.8 (Ultimate eroded).-

U (X) = i∈N U i (X), where U i (X) = ε κi (X) \ δ ∞ εκ i (X) [ε κi+1 (X)] (10.1)
Here, denotes δ ∞ A (B) the reconstruction of the set B in the mask A, as introduced in Chapter 1. It uses the basic structuring element of the grid to make the operation (see below).

The ultimate eroded has the following properties:

-Increasingness -Anti-extensivity -Idempotence (it is thus an opening) -All the U i are disjoint -If we dilate a point of U i by the structuring element κ i , the result is called a maximal ball.

The definition of a discrete ball is close to the one of the continuum: in the continuum we have a continuous family of balls.

The process of obtaining the ultimate eroded is illustrated in Figure 10.4. Each point of an ultimate eroded corresponds to a center of maximal ball, so convex as defined by the family of structuring elements α i . The connected components of the ultimate erodeds identify and possibly separate approximately convex parts of a particle [SER 82]. The ultimate eroded of each component is the last step before his disappearance by successive erosions, hence its name.

Lantuéjoul formula

From its definition, it follows that the ultimate eroded is a subset of the skeleton. We can adapt the formula of ultimate eroded to produce a different definition of the skeleton. This formula is due to Christian Lantuéjoul [LAN 78, LAN 80b]: DEFINITION 10.9 (Lantuéjoul formula).-The skeleton S(X) of a set X is the set of points defined with the following formula: Here, γ 1 is the opening by the unit ball, that is achieved with the smallest structuring element of the chosen granulometric family. In the continuum it may be infinitesimal [MAT 92], but this is not a necessity. In the discrete framework, it is done with α 1 , which is often the basic structuring element of the grid (e.g. a square 3 × 3 in 8-connectivity). This definition is inspired from the ultimate eroded, by replacing the opening by reconstruction by the unit opening. Residuals of the unit opening of each successive erosion are the elements of Lantuéjoul skeleton. This procedure is illustrated in two parts, first Figure 10.5 is the illustration of the part A\γ 1 (A) of the equation 10.2, which corresponds to residues of the opening. Then in Figure 10.6, we show the union of residues for the family erosions A = ε κi (X) with X being a rectangle. In both cases we take a disk as structuring element.

S(X) = i∈N S i (X), where S i (X) = ε κi (X) \ γ 1 [ε κi (X)] (10.2)
In a purely discrete setting, the Lantuéjoul formula still applies. In Figure 10.7 we have the example of the skeleton of Lantuéjoul of a simple discrete figure, taking as unit ball the four-connected neighborhood (the center point plus its four relative nearest neighbors). The key here is to choose a granulometric family ε κi , for example with a constant α i . We have the important following property:

PROPERTY 10.10.-Lantuéjoul skeleton is composed of the center of the maximal balls.

Without bringing a formal proof, we can see why this is so. The LANTUÉJOUL residues at iteration i of the set X are the points kept by an erosion of size i and removed through an opening of size i + 1. They are thus the center of a ball of size i contained in X, but no ball of size i + 1 or higher can be centered at this point.

Let us observe that any set of points not reconstructed by the opening by reconstruction in the formula of ultimate erodeds would also not be reconstructed with the Lantuéjoul formula, thus the following property holds.

PROPERTY 10.11.-The ultimate eroded is a subset of the Lantuéjoul skeleton.

Discrete distances

This Section provides a brief introduction to the main discrete distances used in image processing. For a more detailed treatment, the reader can refer to [THI 07]. From one of these distances, we deduce immediately the notions of ball and maximal ball for a given object X, and we call medial axis of X all the centers of the maximal balls for X.

Let us consider a set E (in the sequel we generally take E = Z n , with, in most of the examples, n = 2), a distance on E is an application d from E × E to R + that satisfies:

∀x, y ∈ E, d(x, y) = d(y, x) (symmetry) (10.3) ∀x, y ∈ E, d(x, y) = 0 ⇔ x = y (separation) (10.4) ∀x, y, z ∈ E, d(x, z) ≤ d(x, y) + d(y, z) (triangular inegality) (10.5) If we consider a subset Y of E and a point x in E, the distance from x to Y is defined by d(x, Y ) = min{d(x, y) | y ∈ Y }.
Let X be a strict subset of E, called "object", the distance map of X is the application D X from E to R + defined by:

∀x ∈ E, D X (x) = d(x, X) (10.6)
where X denotes the complement of X in E. By abuse of language, we will retain the name "distance map" even if d does not verify the triangular inequality.

In digital image processing, we have extensively used the distances known under the names of "city-block distance" (or "Manhattan distance") and "chessboard distance" (in 2D), because these distances are the easiest to calculate [ROS 68]. We denote them respectively d 4 and d 8 , referring to the number of points which are distance 1 of a given point. They are defined by:

d 4 (x, y) = n i=1 |y i -x i | (10.7) d 8 (x, y) = n sup i=1 |y i -x i | (10.8)
The major problem of these distances is their non-rotational invariance. In practice this means that if one performs a distance measurement on an object taken from a digital image, we can obtain significantly different results depending on the orientation of the object during the shooting. Note that a full rotation invariance can not be achieved when one uses discrete images, however we would like that the effects of a rotation on the measurement of distances do not exceed the size of the discretization step.

In an attempt to overcome this defect of distances d 4 and d 8 , the so-called chamfer distances have been introduced and studied [MON 68, BOR 84]. To define these distances we must assume that the set E is equipped with a module structure (see [THI 07]). We call chamfer mask a finite set of pairs consisting of a displacement v k and a weight

p k : M = {(v k , p k ) | v k ∈ E, p k ∈ R + , k ∈ {1, .
. . , m}} (10.9) satisfying the following conditions:

i) each of the v k and of the p k is not null, ii) M has a central symmetry, that is to say that if (v, p) is in M then (-v, p) is also in M , iii) M contains a base of E, that is to say that for all x in E, there exists a m-uplet (a 1 . . . a m ) of positive or null integers such that 1≤k≤m a k v k = x.

Let x and y two points of E, we set

d M (x, y) = min{ 1≤k≤m a k p k | y = x + 1≤k≤m a k v k , a k ∈ N} (10.10)
For example, the following mask allows to retrieve the d 4 distance: M 4 = {((1, 0), 1), ((0, 1), 1), ((-1, 0), 1), ((0, -1), 1)}. The mask M 5,7 = {((1, 0), 5), ((1, 1), 7), ((0, 1), 5), ((-1, 1), 7), ((-1, 0), 5), ((-1, -1), 7), ((0, -1), 5), ((1, -1), 7)} is often used in 2D. By varying the number of pairs in the mask and the weights used, it is possible to obtain more or less accurate approximations of the Euclidean distance, which we discuss below. We can compare on Figure 10.8 the forms of the balls obtained from different distances.

However, the lowest sensitivity to the effects of rotations can only be achieved by using the Euclidean distance d E :

d 2 E (x, y) = (y -x) 2 = n i=1 (y i -x i ) 2
(10.11)

d E (x, y) = d 2 E (x, y) (10.12)
The Euclidean squared distance d 2 E , which is not a distance because it does not verify the triangukar inequality is however sufficient in many applications, and easier to handle computationally than the Euclidean distance: indeed, it only involves integers when the points are integer coordinates. Moreover, if we can calculate a distance map relative to d 2 E , it is easy to deduce a map of Euclidean distance.

To calculate a distance map relative to d 2 E , the naive algorithm (direct application of the definition) has a quadratic complexity in the number of points in the image, which is very inefficient in practice. In 1980, Danielsson [DAN 80] proposes a linear time algorithm for computing an very correct approximation for this kind of distance map for two-dimensional images. However, it is only relatively recently that this problem has received a truly satisfactory solution: an algorithm giving an exact result, linear in computation time, and that can be generalized whatever the dimension n of the space. This algorithm, published in 1996 by Hirata [HIR 96], based on earlier work [SAI 94] of Saito and Toriwaki, and was found independently by Meijster et al. [MEI 00]. Another approach [MAU 03], based on the notion of Voronoï diagram, also provides an efficient and accurate algorithm. In both cases, these algorithms called "separable", compute recursively a distance map of dimension n from maps of dimension n -1 computed independently of each other (line by line, plan by plan, etc.) We describe in detail the first of these algorithms in the case n = 2.

Let us set

E = [0 . . . ℓ -1] × [0 . . . h -1].
To compute the distance map D X relative to a subset X of E, we have to compute, for any point x of E, the quantity

d X (x) = min{(y 1 -x 1 ) 2 + (y 2 -x 2 ) 2 | y ∈ X} (10.13) Let us set, for any j ∈ [0 . . . h -1] d X1 (x 1 , j) = min{(y 1 -x 1 ) 2 | (y 1 , j) ∈ X} (10.14)
We can write Eq. ( 10.13) under the form

d X (x) = min{d X1 (x 1 , j) + (j -x 2 ) 2 | 0 ≤ j < h} (10.15)
The calculation of the values of d X1 can be done independently, thanks to Eq. (10.14), for each image line, and the d X1 are stored in a table. Then, we calculate the values d X by Eq.( 10.15) independently column by column.

To interpret Eq. (10.15) geometrically, note that when x 1 is fixed, as is the case when performing the calculations for a given column, and for a given j, the expression d X1 (x 1 , j) + (jx 2 ) 2 defines a parabola whose minimum is reached for x 2 = j and is d X1 (x 1 , j). The calculations relative to Eq. (10.15) are therefore tantamount to finding the lower envelope of a family of parabolas (see Figure 10.9).

The algorithm described in [HIR 96, MEI 00] consists in calculating the points of integer abscissa of the envelope in two passes, one considering the ascending halfparabolas and the other, the descending half-parabolas. A stack structure allows to act on each parabola with only two operations: one of stacking, the other of unstacking, which ensures the linearity of the algorithm. Moreover, its implementation requires only twenty lines of code.

Following the same approach, D. Coeurjolly proposed algorithms to calculate with an optimal complexity a subset of the exact Euclidean medial axis, sufficient to reconstruct the original object, to make this reconstruction, and to calculate the function that for each point x of an object X associates the set of points of X which are minimum distance of xd [COE 03, COU 07a]. This fonction, called projection on X, plays an essential role in defining and calculating the bisector function, which is the subject of Section 10.5 of this chapter.

Bisector function

We denote by E the set R n or Z n . Let S be a non-empty subset pf E, and let x ∈ E. The projection of x on S, denoted by Π S (x), is the set of all points y in S that are closest to x; more precisely, Π S (x) = {y ∈ S, ∀z ∈ S, d(y, x) ≤ d(z, x)}.

The bisector angle of a point x in X can be defined, in the continuous framework, as the maximal unsigned angle formed by x (as the vertex) and any two points in the projection of x on X [MEY 79, TAL 92]. In particular, if # Π X (x) = 1, then the bisector angle of x is zero. The bisector function of X is the function which associates to each point x of X, its bisector angle in X.

This very definition of the bisector function was used in [ATT 96] in order to provide a filtering criterion for skeletons based on Voronoi diagrams in the continuous plane. It has been also adapted to the discrete case in [TAL 92,MAL 98,COU 07a]. We give here the definition proposed in [COU 07a].

DEFINITION 10.12.-Let X ⊂ E, and let x ∈ X. The extended projection of x on X, denoted by Π e X (x), is the union of the sets Π X (y), for all y in Γ 4 (x) such that d 2 (y, X) ≤ d 2 (x, X). The (discrete) bisector angle of x in X, denoted by θ X (x), is the maximal unsigned angle between the vectors → xy, → xz, for all y, z in Π e X (x). In particular, if #Π e X (x) = 1, then θ X (x) = 0. The (discrete) bisector function of X, denoted by θ X , is the function which associates to each point x of X, its discrete bisector angle in X.

Thanks to the algorithm introduced in [COU 07a], the Voronoi mapping Π X can be computed in optimal time. For each object point x, we must then compute Π e X (x) using the adjacency relation Γ and the distance map D 2 X . The last step to obtain the bisector angle consists in the computation of the maximum unsigned angle between all the pairs of vectors { of the extended projections for a given shape is usually quite small; thus considering all possible pairs constitutes the best choice. However, sub-quadratic algorithms exist for this task.

In Figure 10.10, we show a set X together with its medial axis (a) and the discrete bisector function θ X (b). We illustrate the use of this function to eliminate spurious points of the medial axis: in (c), we show the points of the medial axis (in black) which have a bisector angle greater than 40 degrees. A zoomed detail of both axes is shown in (d). Notice that only the bisector angles of the medial axis points need to be computed for this application. Figure 10.11 shows the bisector function of a more complex 2D shape.

The definition and computation of this discrete bisector function can be straightforwardly extended to Z 3 . To conclude this section, we present in Figure 10.12 an illustration of the bisector function of a three-dimensional object (a vertebra).

Homotopic transformations

As we saw before, the discrete medial axis is not topologically equivalent, in general, to the original object. Some algorithms used to compute the medial axis proceed by iterated thinning, that is, iterative removal of points from the object. The notion of simple point permits to garantee that such a transformation preserves topology.

Intuitively, a point of an object (a subset of Z n ) is called simple if it can be deleted from this object without altering topology. This notion, pionneered by Duda, Hart, Munson [DUD 67], Golay [GOL 69] and Rosenfeld [ROS 70], has since been the subject of an abundant literature (see e.g., [KON 89]). In particular, local characterizations of simple points have been proposed in 2D, 3D and even 4D, on which efficient implementation of thinning procedures are based [COU 09b].

In Figure 10.13, the points (or pixels) x, y, z, t are not simple: the removal of x from the set X of pixels would create a new connected component of the complement X of X ; the removal of y would merge two connected components of X ; the removal of z would split a connected component of X ; and the removal of t would delete a connected component of X. On the other hand, the pixels a, b and c are simple pixels. We see that, in 2D, the notion of connectedness (for both X and X) suffices to characterize simple pixels.

Things are more difficult in 3D. Consider the example of the set X depicted in Figure 10.14, removing the voxel x or the voxel y from X would not split, merge, create or suppress any component of X nor any component of X. However neither x nor y is simple, for the deletion of x (resp. y) causes the suppression (resp. creation) of a tunnel. The notion of tunnel can be formalized thanks to the fundamental group introduced by Poincaré. This group is a topological invariant that is preserved by any continuous deformation. However, it is still possible to characterize 3D simple points by local conditions that are only based on connectedness (see [BER 94a,BER 94b]), we give one such characterization in the sequel. The fact that a global notion, the one of simple point, can be characterized by a local test, is remarkable. It is even more surprising that connectedness alone suffices to provide such a local characterization in 3D. However, this is no longer true in 4D and higher dimensions. In [COU 09b], a definition of simple point based on the collapse operation is presented, and new local characterizations in 2D, 3D and 4D spaces are introduced. This article also analyzes the difficulties that make it impossible to extend this kind of characterizations to dimension 5 and higher.

Neighborhoods, connectedness

We recall here basic definitions of digital topology for binary images [KON 89].

A point x ∈ Z d (d = 2, 3) is defined by (x 1 , . . . , x d ) with x i ∈ Z.

We consider the neighborhood relations Γ 4 and Γ 8 defined, for any point x ∈ Z 2 , by:

Γ 4 (x) = {y ∈ Z 2 | |y 1 -x 1 | + |y 2 -x 2 | ≤ 1}, Γ 8 (x) = {y ∈ Z 2 | max(|y 1 -x 1 |, |y 2 -x 2 |) ≤ 1},
and the neighborhood relations Γ 6 , and Γ 26 and Γ 18 defined, for any point x ∈ Z 3 , by:

Γ 6 (x) = {y ∈ Z 3 | |y 1 -x 1 | + |y 2 -x 2 | + |y 3 -x 3 | ≤ 1}, Γ 26 (x) = {y ∈ Z 3 | max(|y 1 -x 1 |, |y 2 -x 2 |, |y 3 -x 3 |) ≤ 1}, Γ 18 (x) = {y ∈ Γ 26 (x) | |y 1 -x 1 | + |y 2 -x 2 | + |y 3 -x 3 | ≤ 2}.
In the sequel, we denote by n a number such that n ∈ {4, 8, 6, 18, 26}. We define

Γ * n (x) = Γ n (x) \ {x}. The point y ∈ E said to be n-adjacent to the point x ∈ E if y ∈ Γ * n (x). An n-path is a sequence of points x 0 . . . x k such that x i is n-adjacent to x i-1 for i = 1 . . . k.
Let X ⊆ E, we say that two points x, y of X are n-connected in X if there exists a n-path in X between these two points. This defines an equivalence relation. The equivalence classes of this relation are the n-connected components of X. A subset X of E is said to be n-connected if it is composed of exactly one n-connected component.

The set of all n-connected components of X is denoted by C n (X). A subset Y of E is said to be n-adjacent to a point x ∈ E if there exists a point y ∈ Y that is n-adjacent to x. The set of n-connected components of X that are n-adjacent to x is denoted by C x n (X). Remark that C n (X) and C x n (X) are sets of subsets of X, not sets of points. Furthermore, if S is a set, we denote by # S the number of its elements.

Connectivity numbers and simple points

Let us first define connectivity numbers in the 2D case. We see later that the 3D case is more complex. Let X ⊆ Z 2 and x ∈ Z 2 . The connectivity number T n (x, X) is defined by: T

n (x, X) = # C x n [Γ * 8 (x) ∩ X]
If we use the n-connectivity for X, then we must use the n-connectivity for X, for example in 2D the 4-connectivity for X is associated with the 8-connectivity for 10.6.3. Homotopic thinning It remains to decide in which order simple points will be deleted. We present, with Algorithm 1, a thinning strategy which consists of controlling this order based on a priority function. This function associates, to each point x of X, an integer or real number P (x), which represents the priority of point x. The points of X having the lowest values of P will be considered first. Some points x can be given the priority P (x) = +∞, meaning that such points cannot be deleted; in other words, these points with infinite priority are anchor points (see [DAV 81,VIN 91a,PUD 98]) for the thinnning.

Algorithm 1: Guided thinning Data: X ⊆ E, P a function from X to Z or R Result: X repeat Let x be an element of X such that x is simple for X and P is minimal ; X = X \ {x} ; until stability;

We will discuss later on (sec. 10.6.5) some specific problems that occur whenever the priority function is an Euclidean distance map.

Sequential and parallel thinning algorithms

Homotopic transforms discussed in this chapter are sequential by nature, in the sense that after each simple point deletion, this modification must be taken into account when testing the simplicity of other points. In other words, preserving topology is not guaranteed if one deletes several simple points simultaneously: for example in Figure 10.13, deleting both simple points a and a ′ would merge two connected components of the background.

Consequently, some arbitrary choices must sometimes be done with respect to the order in which simple points are treated. This may lead to different results for a same object. On the other hand, parallel thinning strategies, which are not covered by this chapter, produce skeletons that are uniquely defined.

Nevertheless, general conditions that permit to guarantee topology preservation while simultaneously deleting several simple points are much more difficult to establish than single simple point characterizations.

Skeleton based on the Euclidean distance

The skeletonization methods which are based on homotopic thinnings, in the sense of section 10.6.3, provide a formal guarantee that the skeleton and the original object have the same topology. The simplest such method consists in computing an ultimate homotopic skeleton of the object X constrained by the medial axis of X, that is, removing iteratively simple points from X which do not belong to MA(X), taking the distance map as a priority function in order to select first the points which are closest to the background. This can be done using the guided thinning procedure described previously (Algorithm 1), with P = D X and Y = MA(X).

The drawback of this method has been well analyzed in [TAL 92]. Roughly speaking the method does not guarantee that points of the homotopic skeleton outside the medial axis are "well centered" in the object. In Figure 10.15, we give an example that illustrates the kind of problem that may occur.

In Figure 10.15b, numbers correspond to squared Euclidean distances from each object point to the nearest background point. The circled point with value 1, is one of the points that belong to the constraint set Y (see Algorithm 1), that is, a point of the medial axis. Assume that all points with value less to or equal to 8 have been treated by the algorithm. At this step, points in gray are still in the object X, as well as the two circled points (the point valued 1 for it belongs to Y , and the one valued 4 for it has been detected as non-simple when it has been examined). All other points are outside of X. Clearly, the point v valued 8, adjacent to z valued 4, will be selected before its neighbor w valued 9, and since it will be simple at this stage, it will be removed from X. This behaviour will be reproduced during subsequent steps, creating a diagonal skeleton branch. This contradicts a property of skeleton in the continuous space, that asserts that such a branch should follow a steepest descent path on the distance map. To check this, let us compute the slopes of segments zv and zw in our example configuration: ( √ 8 -√ 4)/1 ≈ 0.83, and ( √ 9 -√ 4)/ √ 2 ≈ 0.71. Thus, point v should be kept in the skeleton rather than point w, following this criterion .

To solve this problem, a strategy proposed in [COU 07a] consists of defining a priority function that takes into account both the distance map and an auxiliary function defined in the neighborhood of each dynamically detected skeleton point. Let x be such a point, to any neighbor y of x that is still in X but not in Y , we associate the value p y = D X (x) + (D X (y) -D X (x))/d(x, y), with D X (x) = D 2 X (x). The new priority function, for any point y, is defined by min(p y , D X (y)). We see that (D X (y) -D X (x))/d(x, y) is the slope of xy, thus the neighbors of x will be treated in increasing order of slope, since the value p y is always less than or equal to the value D X (y) (for all x, y in Z 2 or Z 3 with x = y, we have d(x, y) ≥ 1).

For example, in the previous case, we have D X (v) = √ 8 ≈ 2.83, D X (w) = 3, p v = √ 4 + ( √ 8 -√ 4)/1 = √ 8 and p w = √ 4 + ( √ 9 -√ 4)/ √ 2 ≈ 2.71; thus the point w will be selected before v with this strategy. The algorithm is described below.

The time complexity of this algorithm depends on the data structure used to represent the sets Q and R. Specifically, this data structure must allow for efficiently choosing (p, x) at the beginning of the "while"loop, and also for fast insertion of new couples. If one uses e.g., a balanced binary tree [COR 09], the overall complexity of the algorithm is O(n log n), where n is the number of points in the image.

Conclusion

In this chapter we introduced the concepts of skeleton, medial axis, distance, granulometry and thinning, and some links between these concepts.

The definition of a series of dilations of a point by a family of structuring elements defines a granulometric family, by analogy with screening in geology. This family produces a notion of distance, and an operator: the distance function. The latter, in addition to its intrisic interest, can be used to define some characteristic points of a

Figure 10 . 2 .

 102 Figure 10.2. Illustration of the concept of firefront. Fire propagates from the border of the object toward its interior, in an isotropic manner. The point C is a quench point and is a center of maximal ball.

Figure 10 . 3 .

 103 Figure 10.3. Granulometric families that are homogeneous (a, b, c) and heterogeneous (d)

Figure 10 . 4 .

 104 Figure 10.4. Process for obtaining the ultimate eroded. In the first row (line A), we have the successive erosions A1, A2, etc., with a unit ball as structuring element. In the second row (line B), the successive binary reconstructions of Ai+1 in Ai. So B2 is the reconstruction of A3 in A2. Third row (line C) we have the residues Ci = Ai -Bi. Finally, the ultimate eroded is the union ∪iCi .

Figure 10 Figure 10 . 6 .

 10106 Figure 10.5. A residue using the Lantuéjoul formula, i.e. the result of A\γ1(A): (a) part of an object A, in light grey, (b) erosion by a unit ball, (c) result of the erosion, in light gray on a dark background, (d) dilation process by the same ball, (e) result of the opening, and its residue in black

Figure 10 . 7 .

 107 Figure 10.7. Lantuéjoul skeleton in the discrete case: in (a) the silhouette of a bear, in (b) its Lantuéjoul skeleton, with as unit ball the four-connected neighborhood, in (c) the ultimate eroded with the same unit ball

Figure 10 . 8 .

 108 Figure 10.8. Discrete balls obtained from the d4, d8 and dM 5,7 distances and from the Euclidean distance dE

Figure 10 . 9 .

 109 Figure10.9. Parabolas obtained for: j = 3, dX1(x1, j) = 4; j = 4, dX1(x1, j) = 1; j = 5, dX1(x1, j) = 2; j = 6, dX1(x1, j) = 5.

Figure 10 .Figure 10 .

 1010 Figure 10.10. (a): a set X and its medial axis [RÉM 05] (in black); (b): the bisector function θX (dark colors correspond to wide angles); (c): filtered medial axis, based on the values of θX ; (d): detail of the non-filtered and filtered medial axis.

Figure 10 .

 10 Figure 10.12. (a): a view of a subset X of Z 3 (vertebra). (b,c,d): the bisector function, illustrated in an "X-ray" manner: the gray level of a point corresponds to the average of the bisector angles on a straight line parallel to one of the three axes.

Figure 10 . 13 .

 1013 Figure 10.13. Illustration of 2D simple pixels. The set X is made of the pixels in gray, a, b, c are simple whereas x, y, z, t are not simple.

  Numerous attempts have been done during the last 40 years to solve this problem [COU 06a]. The notion of minimal non-simple set, introduced by C. Ronse [RON 88], allows for testing whether a

  Numbers indicate the distance map values. . simple points can be removed while preserving topology. The notion of P-simple point proposed in[BER 95] has the same goal and constitutes fulthermore a general algorithmic scheme for designing 3D thinning algorithms, such that topology preservation is guaranteed by construction. Recently, a general framework for the study of parallel thinning in arbitrary dimension has been developped by G.Bertrand [BER 07a]. This framework, centered on the notion of critical kernel, generalizes both the one of minimal non-simple sets and the one pf P-simple points. The interested reader can find a complete review on critical kernels on the site http://www.esiee.fr/~info/ck .

	Figure 10.15. Illustration of the geometric distorsion of a skeleton. (a): A part of an object
	(in gray), skeleton points found by Algorithm 1 with the Euclidean distance map as priority
	function (circles), steepest descent path with regard to the same map (squares). (b): One step
	of thinning (see text).

Here we admit we need to define the notion of tangent disk to the border ∂X of X. As a regularity condition, we can for instance take that ∂X be everywhere twice differentiable.

X (and vice-versa), and in 3D the 6-connectivity for X is associated with the 18-or the 26-connectivity for X (and vice-versa). This is necessary to have a correspondence between the topological characteristics of X and those of X (see [KON 89]). Furthermore, it is sometimes necessary, in 3D, to make a distinction between the 6connectivity that is associated with the 18-connectivity, and the 6-connectivity that is associated with the 26-connectivity. In order to make this distinction explicit, we use the symbol 6 + to denote the 6-connectivity that is associated with the 18connectivity (see [BER 94a]). To summarize, we have the following possibilities in 2D: (n, n) = (4, 8) or (8, 4); and in 3D: (n, n) = (6, 26), (26, 6), (6 + , 18) or (18, 6 + ).

In the 3D case, the definition of connectivity numbers is based on the notion of geodesic neighborhood. Let X ⊆ Z 3 and x ∈ Z 3 , the geodesic n-neighborhood of x in X of order k is the set Γ k n (x, X) recursively defined by: Γ

is the set composed of points y of Γ * 26 (x) ∩ X such that there exists a n-path π from x to y, the length of which is less than or equal to k, under the condition that all points of π but x must belong to Γ * 26 (x) ∩ X. The geodesic neighborhoods G n (x, X) are defined by: G

We can now give a definition of connectivity numbers in 3D [MAL 93,BER 94a,BER 94b]. DEFINITION 10.13.-Let X ⊆ Z 3 and x ∈ Z 3 . The connectivity number T n (x, X) is defined by: T

Notice that a formulation in terms of geodesic neighborhoods also permits to retrieve the definition of connectivity numbers in 2D.

If we use the n-connectivity for X and the n-connectivity for X, the connectivity numbers T n (x, X) and T n (x, X) describe the topological characteristics of point x in the object X. In particular, connectivity numbers can be used to detect whether a point is simple or not [BER 94a, BER 94b], in 2D and in 3D: PROPERTY 10.14.-Let X ⊆ E and x ∈ X. The point x is n-simple if and only if T n (x, X) = 1 and T n (x, X) = 1.

To give an intuitive interpretation of this caracterization, a point is simple if and only if there is in its neighborhood exactly one "object"component and one "background"component.

Algorithm 2: Euclidean skeleton

Data: X ⊆ E, D X the Euclidean distance map of X, Y ⊆ X Result: Z Z =X; Q ={(D X (x), x); where x is an arbitrary point of X \ Y }; R ={(p x , x); where x is an arbitrary point of X \ Y adjacent to Y , and where

binary form, the centers of maximal balls. We showed how to compute these centers in continuous and discrete spaces, by the means of a granulometric family.

To retrieve certain desirable properties of the continuous skeleton in the discrete space, in particular the conservation of topological characteristics, we have been led to introduce the concepts of simple point and homotopic thinning. These concepts apply in dimensions greater than 2, in grayscale images and also in association to discrete Euclidean metric.

With certain precautions, it is now possible do define a skeleton with good robustness properties and rotation invariance. Studying the stability of skeletons is still a subject of active research, with contributions from different disciplines, particular from computational geometry. Attali, Boissonnat and Edelsbrunner identify in [ATT 09] prominent contributions in this area. These include in particular the λ-medial axis introduced by Chazal and Lieutier in [CHA 05], having a continuity property limited to certain values of the parameter λ, which has recently been adapted to the discrete framework in [CHA 09].