J Richefeu

A Manzanera

A MORPHOLOGICAL DOMINANT POINTS DETECTION AND ITS CELLULAR IMPLEMENTATION

Keywords: dominant point, mathematical morphology, level set, skeleton, parallel cellular machine

This paper addresses the problem of detection of dominant points also called interest points or key points in the literature. We present a new algorithm to compute dominant points on a gray level image. It is founded on mathematical morphology, and it is adapted to SIMD cellular, level sets based, Boolean computation. It is also compliant with the morphological scale-space framework. Cellular implementation and results are presented.

INTRODUCTION AND PRELIMINARIES

An interest point is a location in the image where the signal changes in several directions, thus carrying more information. Examples include corners (L-points with more or less important angle) and T-junctions, as well as locations where texture varies significantly [START_REF] Schmid | Comparing and evaluating interest points[END_REF]. Such points are useful in many applications, such as: image matching (correlation) as part of motion analysis or image registration processes and 3D reconstruction in computer vision or image synthesis.

A very large set of interest point detectors has already been proposed in the literature. References for quasi-exhaustive surveys are [START_REF] Deriche | A computional approach for corner and vertex detection[END_REF] [START_REF] Schmid | Comparing and evaluating interest points[END_REF]. We can file it into three categories.

Approaches in the first group involve first extracting edges as a chain code and then searching for points of maximal curvature [START_REF] Asada | The curvature primal sketch[END_REF] [START_REF] Deriche | A computional approach for corner and vertex detection[END_REF], or performing a polygonal approximation on the chains and then searching for the line segment intersections [START_REF] Horaud | Finding geometric and relational structures in an image[END_REF].

The second group consists of approaches working directly on the gray-level image. These techniques are based either on empirical local operators [START_REF] Moravec | Visual mapping by a robot rover[END_REF] or on the measurement of the gradients and curvatures of the surface. A great deal of effort has been devoted to this second type of approach in order to detect corners and vertices [START_REF] Harris | A combined corner and edge detector[END_REF].

The last group deals with theoretical models of the signal in order to obtain sub-pixel precision. Such methods are used only for specific types of interest points like corners [START_REF] Heitger | Simulation of neural contour mechanism: from simple to end-stopped cells[END_REF].

Few authors have considered morphological dominant points. One work [START_REF] Zhang | A morphological algorithm for detecting dominant points on digital curves[END_REF] deals with binary images to detect convex corner points on closed contours. It is based on residues of morphological opening.

Our approach consists in defining an interest function based on a measure of local curvature on the level sets of the gray-level image (section 2). We then use the skeleton to approximate the computation of the high curvature points (section 3).

This algorithm is well suited to SIMD cellular processing. We present the implementation on the programmable retina (section 4) that allows efficient parallel level sets based image processing. Results are finally discussed (section 5).

THE INTEREST FUNCTION

In this section, we present the concept of morphological level set processing used for our dominant points detection. Then we use it to define a gray-level interest function from the binary local curvature.

level set processing

Mathematical morphology is a nonlinear image processing framework developed by G.Matheron and J.Serra [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. Level sets represent a natural way to extend the morphological operator from binary to gray-level images.

Let E be any set. A binary image X is a subset of E.

Let K be an interval [0, n] of N. A grayscale image is a function I from E to K.
Such a grayscale image I can be represented by its level sets {I t } t∈ [1,n] defined as follows:

I t = {p ∈ E|I(p) ≥ t} (1)
So reciprocally,

I(p) = max{t; p ∈ I t } (2)

local curvature and definition of interest function

For binary images, one natural interest function is provided by the local curvature.

In Euclidian space, the curvature is defined as the instantaneous rate of slope change [START_REF] Asada | The curvature primal sketch[END_REF]. Assume that y(x) represents a curve function, and R p is the curvature at point p, then R p is represented by its derivatives :

R p = ∂ 2 y ∂x 2 [1 + ( ∂y ∂x ) 2 ] 3/2 (3)
The curvature is reciprocal of the radius of the maximal inscribed circle tangent to the curve at point p.

Let

• X be the contour of X. The interest function ϕ X associated with the binary image X is then defined as:

ϕ X (p) = R p if p ∈ • X 0 elsewhere (4)
According to the level sets framework, we can now define the interest function φ for the gray-level image I by the sum of the binary interest functions ϕ It over the level sets I t of the image.

φ : E → K p → n t=1 ϕ It (p) (5)
It can be observed that by replacing in equation ( 4) R p by a constant value we obtain the morphological gradient. So the function φ provides a measure of interest on gray-level images by weighting the local contrast by the curvature.

THE SKELETON

As the interest points will be computed by thresholding the interest function, it can be sufficient to compute only the high curvature points. Such points can be found by using a thinning algorithm.

The skeleton was introduced by H.Blum [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF] as the medial axis which is the union of the centers of the maximal inscribed balls. According to H.Blum, the highest curvature points are located on the maximal inscribed circles whose centers are the extremities of branches of the skeleton. Furthermore, this curvature is reciprocal of the distance between the contour and the extremal point of the skeleton.

For discrete images, one popular approach to object skeletonizing is to thin an object (i.e. remove iteratively contour points while preserving topology and geometry) until a minimally connected set of pixels is achieved.

The skeleton algorithm used in our dominant point detector is the MB1-Hybrid thinning algorithm detailed in [START_REF] Manzanera | Metrical properties of a family of 2D parallel thinning algorithm[END_REF]. The advantages of this algorithm are: (1) it has a certain rotation invariance; (2) it produces few spurious branches;

(3) its computational cost is low; (4) there exist 8-connected and 4-connected versions having the same geometrical properties.

One complete iteration of the MB1-Hybrid thinning algorithm actually corresponds to one iteration of a fully-parallel algorithm followed by four iterations of a semi-parallel directional algorithm in the four directions successively (North, South, East and West). Such a combination allows better rotation invariance. Approximately, it produces a skeleton branch for every angle inferior to π/2 any orientation. Boolean templates of the algorithm and implementation details can be found in [START_REF] Manzanera | Metrical properties of a family of 2D parallel thinning algorithm[END_REF].

In order to obtain the positive high curvature points on a binary image, we compute one complete iteration of the MB1-Hybrid skeleton then we calculate the extremal points of the skeleton (those which have one neighbor only in the image). In the same way, the negative high curvature points are computed thanks to the exoskeleton (skeleton of the complementary of the image). Figure 2 shows an example of this process. The previous procedure allows us to compute the convex and concave corner points (L-points). We can also calculate the junction points (T-points) by finding the multiple points of the skeletons. Nevertheless, by duality, the junction points are associated to corner points of the complementary, and the multiple points of the skeleton (resp. exoskeleton) will be detected by the presence of extremal points of the exoskeleton (resp. skeleton), as long as different connexity levels are used for the image and its complementary (discrete Jordan's theorem in the square grid [START_REF] Kong | Digital topology: Introduction and survey[END_REF]). For example, on the configuration shown on Figure 1, the junction is detected by either the presence of 4 extremal points of the 4-connected exoskeleton (points marked "+") or by the presence of 4 extremal points of the 4-connected skeleton (points marked "x"). This procedure will be performed on every level set of the gray-level image in order to obtain the interest function Figure 2: The original binary image, 8-connected skeleton (5 iterations), 4-connected exoskeleton (5 iterations) and union of the extremal points of the skeleton (points marked "x") and the exoskeleton (points marked "+"). by summing the results.

IMPLEMENTATION

The programmable retina [START_REF] Bernard | Des rétines artificielles intelligentes[END_REF] is a CMOS array sensor in which a Boolean elementary processor has been integrated inside every pixel, so it is a cellular SIMD machine with optical input. This architecture has been shown [START_REF] Manzanera | Vision artificielle rétinienne[END_REF] to be well suited for non-linear (Boolean) and level-by-level processing, thanks to an analog-to-digital conversion by multiple thresholding, that allows data processing during the acquisition.

In such architecture, the procedure used to detect dominant points is the following one : For every time index t:

1. acquisition of level set I t by threshold ; The overall number of elementary Boolean parallel operations performed by level t is less than 150. So if the programmable retina is controlled at a 10 Mhz frequence, the computation time for all the level sets is less than 4 ms for 256 gray-levels and less than 0.5 ms for 32 gray-levels.

copy and inversion of

RESULTS

We show in this section some results obtained with our simulation software of the retina architecture. The first test image is taken from the "MOVI" motion of the INRIA database and the second is a detail of the "Goldhill" test image. A hierarchical selection of the level sets defines a morphological scale-space [START_REF] Meyer | Morphological scalescale representation with levelings[END_REF]. As our interest function is a sum over the level sets, it is necessarily decreasing as the morphological scale increases. Our detector is thus compliant with such morphological scale-space as no new dominant point appears as the scale increases. An example is shown on Figure 5 using the alternated sequential filters by reconstruction.

CONCLUSION

In this paper, we have proposed a new dominant point detection algorithm based on mathematical morphology and in particular on skeletons. The algorithm can be executed in parallel cellular machines like the programmable retina very efficiently. The MB1-Hybrid thinning algorithm used in the implementation is a good trade-off between rotation invariance, efficiency and robustness. The number of dominant points can be tuned according to the contrast thanks to thresholding and according to the size of the level sets through the scale-space framework.

This study is part of our current research in the field of motion analysis. In fact, dominant points detection is the first step involved in correlation to match different images on a sequence. In following works, we will automatically detect and study motion on a video system based on the artificial retina.

Figure 1 :

 1 Figure 1: Duality junction/corner and the discrete Jordan's theorem.

  I t and I c t = N OT (I t ); 3. one iteration of thinning of images S(I t ) (skeleton) and S(I c t ) (exoskeleton); 4. search for extremal point ep(S(I)) and ep(S(I c t )); 5. logical "OR" between the two F t = ep(S(I t )) OR ep(S(I c t ));6. accumulation of the sumF = F + F t endForThe result of the previous procedure F is the interest function, dominant points are then obtained by thresholding F .

Figure 3 :

 3 Figure 3: The resulting dominant points detection image on "MOVI" test image.

Figure 3

 3 Figure 3 shows an example of dominant points detection computed on 256 gray-levels by thresholding the interest function at the value of 13. See Figure 4 for a detail.A hierarchical selection of the level sets defines a morphological scale-space[START_REF] Meyer | Morphological scalescale representation with levelings[END_REF]. As our interest function is a sum over the level sets, it is necessarily decreasing as the morphological scale increases. Our detector is thus compliant with such morphological scale-space as no new dominant point appears as the scale increases. An example is shown on Figure5using the alternated sequential filters by reconstruction.

Figure 4 :

 4 Figure 4: A detail of the interest function computed on 256 gray levels and the corresponding dominant points for the threshold at the value of 13.

Figure 5 :

 5 Figure 5: Dominant point detection over a morphological scale-space on a detail of the "goldhill" image


Antoine Manzanera
email: manzaner@ensta.fr

Thierry M Bernard

Metrical properties of a collection of 2D parallel thinning algorithms

Keywords: digital geometry, parallel thinning algorithm, (4, 8)-median axis

This paper is dedicated to the study of metrical properties of a collection of 2D thinning algorithms that we have proposed. Here, we characterize their underlying metrics and use it to reduce the classical metrical biases that affect thinning algorithms in the square grid. We show that some algorithms from the collection lead to skeletons based on a particular geometry, corresponding to the (4,8)-median axis, which is a new shape descriptor, featuring nice robustness and conditioning properties.

Introduction

We have recently proposed a family of 2-dimensional thinning algorithms to compute the skeleton of binary discrete images. We systematically constructed them from the discretization of the evolution equation of a monotonous propagating front, under topology preservation constraints, in the different connectivity models of the square grid, and for different parallelization schemes. Logic minimization was an important issue in their genesis. The eight algorithms, referred to as MB, are summarized on Figure 1, with a comparison of their different properties. The name of every algorithm is given according to the three binary labels: (1) -fp or -dir for fully parallel or directional, (2) -1 or -2 depending on the number of directions of propagation (4 or 8 respectively), and (3) -4 or -8 depending on the topology. The Boolean definition of every algorithm is presented on Figure 1: the principle of the algorithms is to delete iteratively all pixels matching the removing condition, provided that they do not match the non removing condition. For the fully parallel algorithms, all patterns are to be considered with their π/2 rotated versions. For the MBfpx-4 algorithms (x = 1 or 2), a special convention is used: the white pixels with a red dot (resp. black pixels with a green square) are the black pixels of the original image matching (resp. not matching) the removing condition. The derivation, proof, and details of implementation of each one of these algorithms can be found in [START_REF] Manzanera | MB: A coherent collection of 2D parallel thinning algorithms[END_REF]. Note that MBfp1-8 is equivalent to the algorithm proposed in [START_REF] Eckhardt | Invariant thinning[END_REF], and MBfp1-4 is equivalent to the algorithm proposed in [START_REF] Latecki | Well-composed sets[END_REF]. Figure 1 synthesizes some combinatorial (e.g. Boolean complexity and support) and topological (e.g. P-simpleness [START_REF] Bertrand | On P-simple points[END_REF]) properties that are not addressed in this paper, but detailed in [START_REF] Manzanera | MB: A coherent collection of 2D parallel thinning algorithms[END_REF]. The aim of this paper is to study the metrical properties of these algorithms. In Section 2, we show that the geometry of the MB skeletons can be formally characterized by the type of median axis that they each contain. In particular, the (4,8)-median axis is defined as the mixed case of a generic median axis including the classical morphological skeletons for the two canonical distances of the square grid. We show that the (4,8)-median axis is the locus of the cen-ters of the maximal elements from a collection of sets called (4,8)-fuzzy balls, which are formally defined. In Section 3, we show how the different underlying metrics of the algorithms lead to different behaviors with respect to rotation invariance and noise immunity, and discuss the issue of approximating Euclidean skeletons with thinning algorithms in the square grid. Conclusions and perspective work are presented in Section 4.

(K-P)-median axis and metrical properties of MB

As shape descriptor, it is obvious that the metrical properties of a skeleton are very important ; the initial shape must be recovered at least approximately from its weighted skeleton, slight variations on the contour should not lead to significant changes in the skeleton, and the skeleton must be fairly invariant to arbitrary rotations or scalings. Nevertheless, these issues are rather poorly addressed by the thinning approaches. There are other skeletonization methods designed in Euclidean frameworks, either continuous [START_REF] Siddiqi | Hamilton-Jacobi skeletons[END_REF] or discrete [START_REF] Falcão | Multiscale skeletons by image foresting transform and its application to neuromorphometry[END_REF], that address explicitly the metrical issues, at the price of a representation change, implying a higher computational cost and a loss of regularity. We wish to give in this paper a formal description of the metrical behavior of the proposed thinning algorithms, whose results can be seen on Figure 2. The geometry of the skeletons are based on different median axis, depending on the type of parallelism and the directions of deletion. We now recall the formalism needed to introduce our generic median axis:

The discrete plane is mapped to the square grid Z 2 , a (binary) image X is a subset of Z 2 . A pixel x is an element of Z 2 . The two canonical discrete distances of the square grid are respectively the 4-distance d 4 , and the 8-distance

d 8 . If x = (x 1 , x 2 ) and y = (y 1 , y 2 ), then d 4 (x, y) = |x 1 -y 1 | + |x 2 -y 2 | and d 8 (x, y) = max(|x 1 -y 1 |, |x 2 -y 2 |).
The (K, P )-median axis of image X (K and P are equal to 4 or 8 and K ≤ P ) is defined as: formally that, if the pixels are examined in the order induced by the distance d 4 to the border, then: 4

S P K (X) = {x ∈ X; (y ∈ X and d P (x, y) = 1) ⇒ d K (x, X c ) ≥ d K (y, X c )} So the (K, P )-median
(1) S 4 4 (X) ⊂ MBfp1-8(X) and (2) S 8 4 (X) ⊂ MBfp2-8(X).

Indeed, if a pixel p matches pattern (α 1 ) (resp. (α 1 ) or (α 2 )) shown on Figure 3 (or one of their π/2 rotated versions), then it is 4-adjacent (resp. 8-adjacent) to a 4-interior point q, such that d 4 (q, X c ) = d 4 (p, X c )+1, so p ∈ S 4 4 (X) (resp. p ∈ S 8 4 (X)). These properties are verified for most images, because except in pathological cases (some examples can be seen in [START_REF] Manzanera | nD skeletonization: a unified mathematical framework[END_REF]), the MBfp thinning respect the order induced by distance d 4 . It follows that the geometry of the each MB-skeleton is determined by the geometry of the corresponding (K, P )median axis. We are now going to characterize the geometry of the (4, 8)median axis, thanks to the notion of (4, 8)-fuzzy balls, that we define further: For K = 4 or 8, the K-ball of center x and radius n is defined as

B K (x, n) = {z ∈ Z 2 , d K (x, z) ≤ n}. A ball B K (x, n) is said to be maximal in the image X if ∀(y, n ) ∈ Z 2 × N, B K (x, n) ⊂ B K (y, n ) ⊂ X ⇒ (x, n) = (y, n ).
For K = P , it is well known that the (K, K)-median axis corresponds to the union of the centers of maximal K-balls [START_REF] Lantuéjoul | La squelettisation et son application aux mesures topologiques des mosaïques polycristallines[END_REF]. We are going to prove that the (4, 8)-median axis corresponds to the union of the centers of maximal (4, 8)fuzzy balls (see Figure 4), which are recursively defined as follows:

(1) A (4, 8)-fuzzy ball of radius 1 and center x B (4,8) (x, 1) is any set verifying: B 4 (x, 1) ⊂ B (4,8) (x, 1) ⊂ B 8 (x, 1). ( 2) A (4, 8)-fuzzy ball of radius n+1 and center x is a set such that there exists F x n , a (4, 8)-fuzzy ball of center x and radius n such that: B (4,8) (x, n + 1) = y∈F x n B Ky (y, 1), where K y is 4 or 8, depending on y. To prove the identity between the (4, 8)-median axis and the locus of the centers of maximal (4, 8)-fuzzy balls, we need to use the morphological erosion and dilation:

Let b ∈ Z 2 . The translated of X by b is the set X b = {x + b; x ∈ X}. Let B ⊂ Z 2 .
The morphological dilation of X by B is defined as:

X ⊕ B = b∈B X -b = {z ∈ Z 2 ; B z ∩ X = ∅}
The morphological erosion of X by B is defined as: and as

B 4 (y, n + 1) ⊂ (B 4 (y, n) ⊕ B 8 (0, 1)) (8)
from (3), we get:

[B (4,8) (x, n) ∪ B 4 (y, n + 1)] ⊂ G y n ⊕ B 8 (0, 1) (9)
Finally, as G y n is a (4, 8)-fuzzy ball of radius n and center y, we conclude thanks to ( 6) and ( 9) that [B (4,8) (x, n) ∪ B 4 (y, n + 1)] is a (4, 8)-fuzzy ball of radius (n + 1) and center y. (1) Right inclusion. Let x be the center of a (4, 8)-fuzzy balls B (4,8) (x, n) that is maximal in X. Now suppose that there exists y ∈ (B 8 (x, 1) ∩ X) such that d 4 (y, X c ) > d 4 (x, X c ). Then we must have B 4 (y, n + 1) ⊂ X. And so:

[B (4,8) (x, n) ∪ B 4 (y, n + 1)] ⊂ X (10) But from lemma 1, [B (4,8) (x, n) ∪ B 4 (y, n + 1)
] is a (4, 8)-fuzzy ball of radius n + 1 (see Figure 5(1)), which is in contradiction with the maximality of B (4,8) (x, n). (2) An element of S (4,8) (X) is center of a maximal (4, 8)-fuzzy ball.

y x z B4(y,p) B4(y,m) B4-8(x,n) (1) (2)
(2) Left inclusion. Let x ∈ S (4,8) (X). Let B (4,8) (x, n) be the biggest (4, 8)fuzzy ball of center x contained in X. We are going to prove that B (4,8) (x, n) is maximal in X.

Suppose that there exists a (4, 8)-fuzzy ball F y m of center y and radius m such that (y, m) = (x, n) and B (4,8) (x, n) ⊂ F y m ⊂ X. We have:

B 4 (y, m) ⊂ F y m (11) B 4 (x, n) ⊂ B (4,8) (x, n) (12)
The erosion of the ball B 4 (y, m) by B 4 (0, n) is a ball B 4 (y, p) containing x. But x ∈ B 4 (y, p) B 4 (0, 1), otherwise it would mean that B 4 (x, n+1) ⊂ F y m ⊂ X, and then B (4,8) (x, n) ∪ B 4 (x, n + 1) would be a (4, 8)-fuzzy ball of center x and radius (n+1) (see lemma 1) contained in X, which is in contradiction with the fact that B (4,8) (x, n) is the biggest (4, 8)-fuzzy ball of center x contained in X. So there must exist z ∈ B 4 (x, 1) (see Figure 5(2)) such that z ∈ B 4 (y, p) B 4 (0, 1). Then B 4 (z, n + 1) ⊂ F y m ⊂ X, and so d 4 (z, X c ) > d 4 (x, X c ). As z ∈ B 4 (x, 1), we get x ∈ S (4,8) (X), which is in contradiction with our hypothesis.

2

We have now identified the relation between the mixed median axis S (4,8) (X) and a particular class of sets, the (4, 8)-fuzzy balls. These balls are a new shape description tool, which interest lies in the robustness of morphological or connected skeletons defined in the square grid, as we shall illustrate in the following section.

Consequences on the geometrical behavior of MB

The fact that the MB-2 algorithms do not distinguish different (4, 8)-fuzzy balls (Figure 6(1)), make them more robust with respect to noise (Figure 6(2)), and rotation (Figure 6(3)). Obviously, the outcome is that only partial reconstructibility is possible, unlike MBfp1-8, that allows exact re-constructibility (Figure 6(4), re-constructibility is performed over the skeleton weighted with the distance to the border, with d 4 balls for MBfp1-8, and octagonal balls as a median choice of (4, 8)-fuzzy balls for MBfp2-8). Nevertheless, the better behavior of the MB-2 algorithms with respect to rotation is limited to the fact that fewer branches are generated (Figure 6(3)), but if the number of branches is fairly stable for MB-2 on the different rotated versions, their relative positions can vary significantly: see for example Figure 7, and the changes in the hierarchy of the branches. This bias is due to the different errors of the underlying distances (d 4 or d 8 for the fully parallel or construction, it turns out that the underlying metrics of the MB hybrid thinning algorithms is generated by the octagonal discrete balls minimizing the maximal error with respect to the Euclidean distance (see an example of such "optimal" octagonal ball on Figure 8). MBhyb1 and MBhyb2 are shown on the last row of Figure 2. This method leads to a significant improvement (compare images of Figure 7 with image at scale 2 of Figure 8), but the bias with respect to Euclidean distance keeps increasing with scale. Another interesting property of the fuzzy metrics of the MB-2 algorithms Optimal octogonal ball of radius 15 Fig. 8. The MB hybrid skeleton shown on image lizard at scales 1, 2 and 3: the radius of the biggest optimal octagonal ball is respectively 41, 21 and 11.

is their ability to be conditioned by a Euclidean or pseudo-Euclidean preprocessed distance. For example, in Figure 9, a chamfer distance transform of support 5 [START_REF] Borgefors | Distance transformations in digital images[END_REF] is computed, and then the thinning algorithms are applied by imposing that the pixels deleted at the same iteration are at the same distance to the border. The difference of behaviors between MB-1 and MB-2 shown on Figure 9, and the good invariance to rotation that shows MB2 in that case is explained by the fact that Euclidean discrete balls are balls is not stable by arbitrary rotation, and this lead to important biases in case of big non Euclidean (4, 8)-fuzzy balls (big 4-balls or 8-balls, typically), so for shapes with large perpendicular straight contours.

Conclusion

We have shown in this paper the consistency in terms of metrical properties of the family of thinning algorithms we have recently proposed. The geometrical behavior, the advantages and limitations of every algorithm have been clearly identified. We have shown in [START_REF] Manzanera | nD skeletonization: a unified mathematical framework[END_REF] that the 8-connected MBfp algorithms could be naturally expressed in the n-dimensional cubic grid, and we have proved the validity of the corresponding algorithms for n = 3. The same extension can be done for the other algorithms, but their validity remains to prove. This will be the subject of future work, in order to get hopefully a unified and cleanly justified thinning methodology for n-dimensional binary images.

Fig. 1 .

 1 Fig. 1. The MB family of parallel thinning algorithms.

Fig. 2 .

 2 Fig. 2. Metrical properties of the MB skeletons.

(α 1 ) ( α 2 )Fig. 3 .

 123 Fig. 3. Deletion patterns of the MBfp-skeletons.

Fig. 4 .

 4 Fig. 4. Some (4, 8)-fuzzy balls of radius 7. The extremal cases of (4, 8)-fuzzy balls are respectively the 4-ball (on the left), and the 8-ball (on the right).

2 Theorem 1 S

 21 [START_REF] Falcão | Multiscale skeletons by image foresting transform and its application to neuromorphometry[END_REF][START_REF] Manzanera | MB: A coherent collection of 2D parallel thinning algorithms[END_REF] (X) is the locus of the centers of maximal (4, 8)-fuzzy balls in X.

Fig. 5 . ( 1 )

 51 Fig. 5. (1) The center of a maximal (4, 8)-fuzzy ball is an element of S (4,8) (X).(2) An element of S (4,8) (X) is center of a maximal (4, 8)-fuzzy ball.

8 MBfp2- 8 Fig. 9 .

 889 Fig. 9. Conditioning the MB thinning algorithms by a distance transform.

  axis is the set of the maxima of distance d K , in the P -neighborhood. Depending on the values of K and P , this leads to three different median axes, shown on row 1 of Figure2. We are going to show that these different median axes determine the metrical properties of the different algorithms. As every iteration (resp. four successive sub-iterations) of MBfp (resp. MBdir) examines the 4-contour (resp. the 8-contour), the geometry of the resulting skeleton is based on distance d 4 (resp. d 8 ), as it can be seen on Figure2. In the case of the MBfp 8-connected algorithms, the isotropy allows even to prove

	Original	S (8,8)	S (4,4)	S (4,8)
	MBdir2-4	MBdir1-4	MBfp1-4	MBfp2-4
	MBdir2-8	MBdir1-8	MBfp1-8	MBfp2-8
	MBhyb1-4	MBhyb2-4	MBhyb1-8	MBhyb2-8

We also need to prove the following lemma: [START_REF] Falcão | Multiscale skeletons by image foresting transform and its application to neuromorphometry[END_REF][START_REF] Manzanera | MB: A coherent collection of 2D parallel thinning algorithms[END_REF]-fuzzy ball of radius n and center x, and if y ∈ B 8 (x, 1), then [B (4,8) (x, n) ∪ B 4 (y, n + 1)] is a (4, 8)-fuzzy ball of radius n + 1 and center y.

Preliminary remark: it is clear, by the definition of fuzzy balls, that if F x n is a (4, 8)-fuzzy ball of center x and radius n, then any set S verifying:

and

Let y ∈ B 8 (x, 1). By induction hypothesis,

)fuzzy ball of radius n and center y.

So, from (3), we get:

On the other hand, we have:

MBfp1-8

MBfp2-8

MBdir1-8

MBdir2-8

MBfp1-8

MBfp2-8

(1)

(2)

MBfp1-8

MBfp2-8

(3) (4) directional algorithms, respectively) with respect to the Euclidean distance, depending on the angles of the objects.

A first basic idea to get a "more Euclidean" geometry for the skeleton is to


Annie Druault-Vicard

Christian Tanguy
email: christian.tanguy@orange-ftgroup.com.

Exact Failure Frequency Calculations for Extended Systems

Keywords: network availability, failure frequency, failure rate, k-out-of-n systems, generating function

This paper shows how the steady-state availability and failure frequency can be calculated in a single pass for very large systems, when the availability is expressed as a product of matrices.

We apply the general procedure to k-out-of-n:G and linear consecutive k-out-of-n:F systems, and to a simple ladder network in which each edge and node may fail. We also give the associated generating functions when the components have identical availabilities and failure rates. For large systems, the failure rate of the whole system is asymptotically proportional to its size. This paves the way to ready-to-use formulae for various architectures, as well as proof that the differential operator approach to failure frequency calculations is very useful and straightforward.

NOTATION

p i , q i [success, failure] probability of component i (q i = 1 -p i ) p, q implies p i = p, q i = q (for edges).

ρ identical availability of nodes (when ρ = p). Steady-state system availability and failure frequency are important performance indices of a repairable system [START_REF] Shooman | Probabilistic reliability: an engineering approach[END_REF], [START_REF] Singh | System reliability modelling and evaluation[END_REF], [START_REF] Colbourn | The Combinatorics of Network Reliability[END_REF], [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], from which other key parameters such as the mean time between failures, average failure rate, Birnbaum importance, etc. may be deduced. In the steady-state regime, the frequency of system failure was first calculated by a cut-set [START_REF] Singh | A new method to determine the failure frequency of a complex system[END_REF] or a tie-set approach [START_REF] Singh | Tie set approach to determine the frequency of system failure[END_REF] in the case of statistically independent failures, which will also be considered here. These approaches are based on the inclusion-exclusion (IE) principle, where the failure or repair rates (more generally, the inverses of the mean down or up times), are adequately given for each term of the relevant expansion.

When all the terms of its IE expansion are kept, the exact availability is obtained as a function of each component availability. Several papers have provided a few simple recipes, describing how the system failure frequency and the failure rate can then be derived [START_REF] Schneeweiss | Computing failure frequency, MTBF & MTTR via mixed products of availabilities and unavailabilities[END_REF], [START_REF] Shi | General formulas for calculating the steady-state frequency of system failure[END_REF], [START_REF] Yuan | Boolean algebra method to calculate network system reliability indices in terms of a proposed FMEA[END_REF]. Recent refinements have been proposed when availability expressions are obtained from various instances (SVI, GVI) of sum-of-disjoint-products algorithms [START_REF] Amari | Generic rules to evaluate system-failure frequency[END_REF]. All these formal calculations boil down to a simple fact: the failure frequency may be derived from the availability through the application of a linear differential operator [START_REF] Schneeweiss | Addendum to: Computing failure frequency via mixed products of availabilities and unavailabilities[END_REF], [START_REF] Hayashi | System failure-frequency analysis using a differential operator[END_REF]. This requires knowledge of the exact availability, which is hard to come by except for trivially small networks, and may have hindered the use of this method.

Unsurprisingly, several algorithms have been put forward, in which availability and failure frequency are computed side by side in a common procedure: triangle-star transformation [START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF], OBDD calculations [START_REF] Chang | Computing System Failure Frequencies and Reliability Importance Measures Using OBDD[END_REF], and another instance of differential operator calculations [START_REF] Hayashi | Transformation from availability expression to failure frequency expression[END_REF].

In this paper, we want to promote the differential operator method for the calculation of the failure frequency by showing it gives the exact result for numerous, widely used configurations, with an arbitrary large number of components. We take advantage of recent results establishing that the availability of recursive networks may be expressed as a product of transfer matrices that take each edge and node availabilities exactly into account [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF], [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of the Brecht-Colbourn ladder and the generalized fan[END_REF], [START_REF] Tanguy | What is the probability of connecting two points ?[END_REF].

Our paper is organized as follows. In Section II, we show how the failure frequency of a system may generally be deduced from the steady-state availability when the latter is expressed by a product of transfer matrices. We first apply this method in Section III, which is devoted to k-out-of-n systems (either k-out-of-n:G or linear consecutive k-out-of-n:F ones) with distinct components. Section IV provides a generic example for the two-terminal failure frequency of a simple ladder network, which has been solved recently for arbitrary edge and node availabilities [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF]; the same procedure could easily be used for more complex networks and their all-terminal reliability too [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of the Brecht-Colbourn ladder and the generalized fan[END_REF], [START_REF] Tanguy | What is the probability of connecting two points ?[END_REF]. In each configuration, we pay attention to the case of identical components, for which the common availability is p (for edges) and ρ (for nodes). For large systems, we show that the asymptotic failure rate has a linear dependence with size, and is given by derivatives of the largest eigenvalue of the unique transfer matrix with respect to p and ρ. We conclude by a brief outlook.

II. GENERAL PROCEDURE

In many systems, as will be explicitly shown in the following sections, the availability A (or the unavailability U) is given by an expression of the form

A = v L M n • • • M 1 v R , (1)
where M k (1 ≤ k ≤ n) is a transfer matrix, the elements of which are multilinear polynomials of individual component availabilities, and where v L and v R are two vectors in which these availabilities do not appear. The mean failure frequency ν is obtained from [START_REF] Schneeweiss | Addendum to: Computing failure frequency via mixed products of availabilities and unavailabilities[END_REF], [START_REF] Hayashi | System failure-frequency analysis using a differential operator[END_REF]

ν = i λ i p i ∂A ∂p i = i µ i q i ∂U ∂q i . (2)
In order to avoid unnecessarily heavy notation, we call M ′ k the matrix obtained by applying the linear differential operator

i λ i p i ∂ ∂p i to M k . Therefore, ν = v L M ′ n M n-1 • • • M 1 v R +v L M n M ′ n-1 • • • M 1 v R + • • • +v L M n M n-1 • • • M ′ 1 v R . (3)
Since M k 's elements are at most linear functions of each p i , the derivation of M ′ k is straightforward. For instance, a matrix element p 1 + p 2 p 3 -p 1 p 2 p 3 in M k would give rise to

λ 1 p 1 + (λ 2 + λ 3 ) p 2 p 3 -(λ 1 + λ 2 + λ 3 ) p 1 p 2 p 3
; the recipes given in [START_REF] Schneeweiss | Computing failure frequency, MTBF & MTTR via mixed products of availabilities and unavailabilities[END_REF], [START_REF] Shi | General formulas for calculating the steady-state frequency of system failure[END_REF], [START_REF] Schneeweiss | Addendum to: Computing failure frequency via mixed products of availabilities and unavailabilities[END_REF] fully apply.

Both availability and failure frequency may be obtained in a single pass in the following way. Let us initialize the procedure by setting

A 1 = M 1 v R , (4)
V 1 = M ′ 1 v R . (5)
The recursion equations are

A k = M k A k-1 , (6)
V k = M k V k-1 + M ′ k A k-1 , (7)
from which we deduce the final results

A = v L A n , (8)
ν = v L V n . (9)
We can now turn to a few 'real-life' applications.

III. k-OUT-OF-n SYSTEMS

k-out-of-n systems are widely used, in various configurations; they have therefore contributed to a huge body of literature (see [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], [START_REF] Chao | Survey of reliability studies of consecutive-k-out-of-n:F & related systems[END_REF], [START_REF] Kuo | Opinions on consecutive-k-out-of-n:F systems[END_REF] and references therein). We start our discussion with these systems because each transfer matrix actually refers to a single equipment only.

A. k-out-of-n:G systems

We first consider the simple k-out-of-n:G system, where each component has an availability

p i (1 ≤ i ≤ n).
To operate as a whole, the system needs at least k elements to function. Its availability A k,n may be written as (see [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], p. 244)

A k,n = 1 -(1, 0, • • • , 0) k Λ n Λ n-1 Λ 1        1 1 . . . 1        k , (10)
with

Λ i =           q i p i 0 • • • 0 0 q i p i 0 0 0 0 . . . . . . 0 . . . . . . 0 q i p i 0 0 • • • 0 q i           k×k . ( 11
)
We have reduced the size of the matrix to a k ×k one, instead of the original

(k +1)×(k +1),
because of the nature of v L = (1, 0, . . . , 0) k and v R in eq. ( 10).

The 'derivative' of Λ i is

Λ ′ i =           -λ i p i λ i p i 0 • • • 0 0 -λ i p i λ i p i 0 0 0 0 . . . . . . 0 . . . . . . 0 -λ i p i λ i p i 0 0 • • • 0 -λ i p i           k×k , (12)
so that the computation of the failure frequency following the method given in section II is straightforward (care should of course be taken of the minus sign in eq. ( 10)).

Let us revisit Example 7.2 of [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF] (see p. 245) for the 5-out-of-8:G system with p i = 0.90, 0.89, ..., 0.83. Assuming a unique repair rate for all components, namely µ, the failure rates

λ i are such that λ i p i = µ (1 -p i ).
From the procedure detailed in Section II, we deduce When all components are identical (p i ≡ p and λ i ≡ λ), only one transfer matrix appears.

Admittedly, A k,n is so simple that a matrix formulation is hardly necessary. Nonetheless, we can give a compact expression for the generating function

G k (z) = ∞ n=0 A k,n z n (the derivation is given in the appendix): G k (k-out-of-n:G; z) = p k z k (1 -z) (1 -(1 -p) z) k . ( 13
)
Since the generating function is a formal power-series expansion, we can apply the linear differential operator λ p ∂ ∂p to eq. ( 13) so that

G k (z) = ∞ n=0 ν k,n z n is easily found to be G(k-out-of-n:G; z) = λ k p k z k (1 -(1 -p) z) k+1 , (14)
which is another formulation of the well-known result

ν k,n = λ k   n l   p k (1-p) n-k (eq. (7.10)
of [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], p. 234).

B. Linear consecutive k-out-of-n:F systems

These systems have been studied in many papers [START_REF] Chao | Survey of reliability studies of consecutive-k-out-of-n:F & related systems[END_REF], [START_REF] Kuo | Opinions on consecutive-k-out-of-n:F systems[END_REF] and a recent textbook [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF]. The reliability A k,n -the probability of operation of a system of n components, which fails if at least k consecutive elements fail -of such a system is given by (see also eq. ( 9.48) of [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], p. 344)

A k,n = (1, 0, . . . , 0) k Λ n Λ n-1 Λ 1        1 1 . . . 1        k , (15)
with

Λ i =           p i q i 0 • • • 0 p i 0 q i 0 . . . . . . . . . 0 . . . 0 p i 0 • • • 0 q i p i 0 0 • • • 0           k×k . (16)
Here again, we have reduced the size of the matrix and the vectors with respect to their original formulation. Consequently,

Λ ′ i =           λ i p i -λ i p i 0 • • • 0 λ i p i 0 -λ i p i 0 . . . . . . . . . 0 . . . 0
λ i p i 0 • • • 0 -λ i p i λ i p i 0 0 • • • 0           k×k , (17)
leading once again to a straightforward calculation of the failure frequency.

A numerical application may be found in the Lin/Con/4/11:F model, as in Example 9.6 of [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF], where the p i 's range from 0.7 to 0.9 by steps of 0.02. Assuming again that the repair rate for each equipment is µ, we get A 4,11 = 30105385968617 30517578125000 ≈ 0.98649329 and a failure frequency ν 4,11 = 155495836041 30517578125000 µ ≈ 0.050953 µ. The corresponding failure rate is then equal to 0.0516505 µ.

For the sake of completeness, we give the generating function for identical components is (see eq. (2.2) of [START_REF] Canfield | Asymptotic reliability of consecutive k-out-of-n systems[END_REF], [START_REF]z)) a partial fraction decomposition of eq[END_REF])

G k (Lin/Con/k/n:F; z) = 1 -(1 -p) k z k 1 -z + p (1 -p) k z k+1 . (18)
Use of eq. ( 18) to obtain the failure frequency generating function is straightforward, and will not be repeated here.

IV. SIMPLE LADDER

We consider in this section the two-terminal availability of a simple ladder network, displayed in Fig. 1, where successive nodes are labelled S i or T j , and where the larger black dots mark the source s and terminal t. This network is a simplified description of a standard architecture for long-hail communication networks: it consists in primary and backup paths, plus additional connections between transit nodes enabling the so-called "local protection" policy by bypassing faulty intermediate nodes or edges. Such an architecture of "absolutely reliable nodes and unreliable edges," with up to 25 edges, was chosen as Example 5 in [START_REF] Heidtmann | Smaller sums of disjoint products by subproduct inversion[END_REF] for a comparison of different "sum of disjoint products" minimizing algorithms, or by Rauzy [START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF] as well as Kuo and collaborators in OBDD test calculations [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF], [START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF], [START_REF] Fu-Min | Analyzing network reliability with imperfect nodes using OBDD[END_REF].

We showed [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF] that the two-terminal availability has a beautiful algebraic structure [START_REF] Shier | Network Reliability and Algebraic Structures[END_REF],

since its exact expression is given by a product of 3 × 3 transfer matrices (see eqs. (ai, bi, ci) and (Si, Ti).

(resp. T n ), we find that [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF], [START_REF] Tanguy | What is the probability of connecting two points ?[END_REF]

R Sn = (1 0 0) M n • • • M 0      1 0 0      , (19)
R Tn = (0 1 0) M n • • • M 0      1 0 0      . (20)
The transfer matrix M n is given by

M n =      a n S n b n c n S n T n a n b n c n S n T n a n b n S n T n c n T n a n b n c n S n T n -a n b n S n T n -b n c n S n T n a n (1 -2 b n ) c n S n T n      . (21)
For n = 0, we must set a 0 = 1, and may choose c 0 = 0 because it does not change the final result. It is worth noting that all five availabilities of the n th "cell" or building block of the network appear in a single transfer matrix M n , which is not sparse, contrary to the matrices of Section III.

Equations [START_REF] Chao | Survey of reliability studies of consecutive-k-out-of-n:F & related systems[END_REF][START_REF] Kuo | Opinions on consecutive-k-out-of-n:F systems[END_REF][START_REF] Canfield | Asymptotic reliability of consecutive k-out-of-n systems[END_REF] apply to the most general ladder in terms of individual availabilities. If an edge or a node is missing, its reliability p i should be set to zero, and its failure rate may be considered arbitrary, because it will not alter the final result. Similarly, if a given edge or node is perfect, its reliability should be equal to one; its failure rate λ i should then, of course, be set to zero.

The associated matrix M ′ n is

M ′ n =      m ′ 11 m ′ 12 m ′ 13 m ′ 21 m ′ 22 m ′ 23 m ′ 31 m ′ 32 m ′ 33      , (22)
with

m ′ 11 = (λ an + λ Sn ) a n S n , m ′ 12 = (λ bn + λ cn + λ Sn + λ Tn ) b n c n S n T n , m ′ 13 = (λ an + λ bn + λ cn + λ Sn + λ Tn ) a n b n c n S n T n , m ′ 21 = (λ an + λ bn + λ Sn + λ Tn ) a n b n S n T n , m ′ 22 = (λ cn + λ Tn ) c n T n , m ′ 23 = m ′ 13 , m ′ 31 = -m ′ 21 , m ′ 32 = -m ′ 12 , m ′ 33 = (λ an + λ cn + λ Sn + λ Tn ) a n c n S n T n -2 m ′ 13 .
When a i = b i = c i ≡ p and S i = T i ≡ ρ the three eigenvalues ζ 0 and ζ ± of the transfer matrix are [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF]

ζ 0 = p ρ (1 -p ρ), (23)
ζ ± = p ρ 2 1 + 2 p (1 -p) ρ ± √ B , (24)
with B = 1 + 4 p 2 ρ -8 p 3 ρ 2 + 4 p 4 ρ 2 . The two-terminal availabilities are [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF]

R Tn (p, ρ) = 1 2 p -ζ n+1 0 + p ρ (1 + p ρ) ζ n+1 + -ζ n+1 - ζ + -ζ - -(1 -2 p + p ρ) p 3 ρ 3 ζ n + -ζ n - ζ + -ζ - , (25)
R Sn (p, ρ) = 1 2 p +ζ n+1 0 + p ρ (1 + p ρ) ζ n+1 + -ζ n+1 - ζ + -ζ - -(1 -2 p + p ρ) p 3 ρ 3 ζ n + -ζ n - ζ + -ζ - . (26)
These expressions are identical except for the ± sign in front of the ζ n+1 0 term. Assuming that the common link failure rate is λ while that for the nodes is ξ, the failure frequency for the S 0 → T n connection is

ν = λ p ∂R Tn (p, ρ) ∂p + ξ ρ ∂R Tn (p, ρ) ∂ρ ; (27)
a similar expression applies to R Sn (p, ρ). When n is large, both availabilities are actually of the form α + ζ n + , because the modulus of ζ + is larger than that of the remaining eigenvalues for 0 ≤ p ≤ 1 [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of a simple ladder network[END_REF]. When nodes are perfect, we have therefore in this limit

ν ≈ λ p ∂α + ∂p ζ n + + n α + ∂ζ + ∂p ζ n-1 + , (28)
so that the failure rate is

λ ≈ λ ∂ ln α + ∂ ln p + n ∂ ln ζ + ∂ ln p , (29)
with

∂ ln ζ + ∂ ln p = -1 + 4 p -6 p 2 + 4 p 3 + (3 -4 p) 1 + 4 p 2 (1 -p) 2 2 (1 -p) 1 + 4 p 2 (1 -p) 2 , (30)
∂ ln α + ∂ ln p = 4 -5 p + 8 p 2 -20 p 3 + 16 p 4 -4 p 5 -(4 -7 p + 4 p 2 -2 p 3 ) 1 + 4 p 2 (1 -p) 2 2 (1 -p) [1 + 4 p 2 (1 -p) 2 ] . (31)
The variations with p of ∂ ln ζ + /∂ ln p and ∂ ln α + /∂ ln p are displayed in Figs. 2 and3.

Since |ζ + | < 1 for p < 1, the contribution of ∂ ln α + /∂ ln p will prevail, and λ will have a linear dependence with n in the large network limit (this is a general property when the eigenvalue of highest modulus is different from unity). For very reliable components, eqs. (29-31) simplify. It is easy to show by a series expansion in the link unavailability q that the global failure rate is given (to first order) by λ → (2 n + 4) λ q; this result could also have been obtained by visual inspection and enumeration of the minimal cuts. 31)). The maximum, reached for p ≈ 0.709902, is about 0.458825.

V. CONCLUSION AND OUTLOOK

We have shown that the linear differential method for computing the failure frequency is a very simple and useful one for k-out-of-n systems as well as the two-terminal availability for recursive networks (this should hold for the all-terminal availability, too [START_REF] Tanguy | What is the probability of connecting two points ?[END_REF]). Its application is not limited to the case of extremely reliable components. Even though we restricted our discussion to expressions dealing with availabilities, a similar treatment could be performed for expressions where unavailabilities are the input data (see eq. ( 2)). For more complex networks, the size of the transfer matrix increases (for instance, it is a 13 × 13 one for the 'street 3 × n' of [START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF]) but the calculations remain straightforward. Finally, the expressions given for steady-state availabilities can also be used for time-dependent systems provided that failures and reparations are still statistically independent events, because the expressions are formally identical (the availabilities of components must be replaced by the reliabilities).

APPENDIX I PROOF OF EQ. (13)

A k,n is given by

so that

G k (z) = p z 1 -(1 -p) z G k-1 (z) = p z 1 -(1 -p) z k G 0 (z) . ( 36
)
Since A 0,n = 1, ∀n, G 0 (z) = 1/(1 -z); eq. ( 13) follows.

  λ i , µ i [failure, repair] rate of component i λ, µ common [failure, repair] rate of components A steady-state availability of the system U steady-state unavailability of the system ν mean failure frequency of the system λ mean failure rate of the system (ν = A λ) M ′ i λ i p i ∂ ∂p i M G(z) generating function for the availability G(z) generating function for the failure frequency A k,n availability of a k-out-of-n:G system I. INTRODUCTION

A 5 , 8

 58 = 615925280183 625000000000 ≈ 0.98548045 and a failure frequency ν 5,8 = 8012914359 156250000000 µ ≈ 0.051283 µ. The failure rate λ 5,8 = ν 5,8 /A 5,8 is then equal to 0.0520382 µ.

Fig. 1 .

 1 Fig. 1. Different source-terminal connections for a simple ladder. Edges and nodes are indexed by their availabilities:

Fig. 2 .

 2 Fig. 2. Variation of ∂ ln ζ+/∂ ln p with p for a simple ladder with perfect nodes (eq. (30)). The maximum, reached for p ≈ 0.251641, is about 1.13827.

Fig. 3 .

 3 Fig.3. Variation of ∂ ln α+/∂ ln p with p for a simple ladder with perfect nodes (eq. (31)). The maximum, reached for

A k,n z n implies 1 z

 1 n+1 = (1 -p) A k,n + p A k-1,n . (34) Setting G k (z) = ∞ n=0 G k (z) = (1 -p) G k (z) + p G k-1 (z) .

The singular and plural of an acronym are always spelled the same. December

28, 2006 DRAFT

December 28, 2006 DRAFT