Mathematical Morphology II: Filtering

1. GEONET VI June 10-17 2007 Mathematical Morphology II:Filtering Adjunction opening and closing Opening by connected reconstruction Morphological filtering Leveling Case of the numerical functions Jean Serra, ESIEE

2. GEONET VI June 10-17 2007 Mathematical Morphology I:Structures Image analysis Complete lattices Remarquable elements Types of lattices Case of numerical functions Lattices of opérateurs Jean Serra, ESIEE

3. Binary opening and closing When B is a disc, the opening amends the caps, removes the small islands and opens isthmuses. When B is a disc, the closing closes the channels, fills completely the small lakes, and partly the gulfs.

4. The three basic properties of adjunction openings de and closings ed are also the axiomsfor the algebraic notion of an opening and a closing. Definition :In algebra, any transformation which is: increasing, anti-extensive and idempotent is called an (algebraic) opening, increasing, extensive and idempotent is called a (algebraic) closing. Particular cases : Here are two very easy ways for creating algebraic openings and closing, they are often used in practice: 1) Compute various opening (closing) and take their supremum (or the infimum in case of closings). 2) Use a reconstruction process (see V- 24 and VI - 14) Algebraic Opening and Closing

5. Let B be the image of lattice L under the algebraic opening g, i.e. B = g (L). Since g is idempotent, set B generates the family of invariant sets of g : b Î BÛ g (b) = B . 1/ Classe B is closed under sup. For any family {bj, j Î J } 5B, we have g (ƒbj, j Î J) 8ƒ{g (bj) , j Î J } = ƒ( bj, j Î J ) by increasingness, and the inverse inequality by anti-extensivity de g. Moreover , 0 Î B. Note that g does not commute under supremum. 2/ Therefore, g is the smallest extension to L of the identity on B, i.e. g (x) = ƒ { b : b Î B , b 7 x } , x Î L (1) . [ The right member is an invariant set of g smaller than x, but also that contains g (x). ] Invariant Elements

6. Conversely, let B be the class closed under sup. that an arbitrary class B0 generates. With each b Î B~0 associate dilation db(a) = b when a 8 b ; db(a) = 0 when not . By adjunction, it yields opening b when b 7 x , gb(x) = ƒ {dB(a) , db(a) 7 x } = (2) 0 when not . [ If b 7 x, it suffices to take a 8 b for obtaining a term = b in the ƒ ; if not, no db(a) ¹ 0 is 7 x ] . Then relation (1) becomes g (x) = ƒ{ db(a) , b Î B }. Theorem (G. Matheron , J. Serra ) : 1/ With any class B05 L corresponds a unique opening on L that has for invariant sets the b Î B0 . 2/ Every algebraic opening g on lattice L may be represented as the union of the adjunction openings linked to the invariant sets of g by relation (2). Representation of openings

7. Application : For creating openings with specific selection properties, one can use structuring elements with various shapes and take their supremum. Suprema of Openings • Theorem : • Any supremum of openings is still an opening. • Any infimum of closings is still a closing. Opening by sup. Opening by Original Opening by

8. Definition The set geodesic distance dX: Rn-Rn® R, with respect to reference set X is defined by dX(x,y) = Inf. of the lengths of the paths going from x to and included X; dX(x,y) = + T, if no such path exits. Properties 1) dX is a generalized distance, since dX(x,y) = dX(y,x) dX(x,y) = 0 Û x = y dX(x,z) £dX(x,y) + dX(y,z) 2) The geodesic distance is always larger than the Euclidean one; 3) A geodesic segment may not be unique Geodesic Distance Examples of geodesics in R2 : x ¥ z y' y X N.B. the portions of geodesics included in the interior of X are line segments

9. The notion of geodesic path is seldom used. By contrast, the notion of geodesic discs appears very often: BX,l(z) = {y, dX(z,y) £l} When the radius r increases, the discs progress as a wave front emitted from z inside the medium X. For a given radius l, the discs BX,l can be viewed as a set of structuring elements which vary from place to place. Geodesic Discs z BX,l(z) X

10. In the digital metrics on Zn , and when d(x) stands for the unit ball centered at point x, then the unit geodesic dilation is defined by the relation : dX(Y) = d(Y) Ç X The dilation of size n is then obtained by iteration : dX,n(Y) = d(n)X(Y) , with d(n)X(Y) = d( ... d(d(Y) Ç X) Ç X ... ) Ç X Note that the geodesic dilations are not translation invariant. (Binary) Digital Geodesic Dilation X Y dX(Y) d(2)X(Y)

11. Reconstruction Opening Y • Given X, the infinite dilation of Y dX,¥ (Y) = È{dX,l(Y) , l>0 } is a closing ; but if we consider it as an operation on the (now variable) reference setX, for a given marker Y, then dX,¥ (Y) turns out to be the reconstruction opening grec(X ; Y) = ~{ dX,l(Y) , l>0 }= ~{ gy,(X),y ÎY} of those connected components of set X that contain at least one point of set Y. X N.B.: As the grid spacing becomes finer and finer, the digital reconstruction opening tends towards the Euclidean one iff X is locally finite union of disjoint compact sets.

12. The notion of a connection allows to express, and to generalize reconstruction openings as follows 1) Call increasing binary criterionany mapping c: s(E) ®{0,1} such that: AÍB Þ c(A) £ c(B) 2) With each criterion c associate the trivial openinggT: s(E) ® s(E)gT (A) = A if c(A) = 1 gT (A) = Æ if c(A) = 0 3) By generalizing the geodesic case, we will say that grecis a reconstruction openingaccording to criterion c when : grec = ƒ{gTgx, xÎE } grec acts independently on the various components of the set under study, by keeping or removing them according as they satisfy the criterion, or not (e.g. area, Ferret diameter, volume..). Connection and Reconstruction Opening

13. The closing by reconstructionjrec=¥grec¥ is defined by duality. For example, in R2, if we take the criterion - « to have an area810 », then jrec(A) is the union of A and of the pores of A with an area 710; - or the criterion «to hit a given markerM », then jrec(A) is the union of A and of the pores of A included in Mc . Associated Lattices: We now considera family {greci} of openings by reconstruction, of criteria { ci}.Their inf €greci is still an opening by reconstruction, where each grain of A which is left must fulfill all criteria  ci , and where the sup ~greci is the opening where one criterion at least must be satisfied (dual results for the closings). Hence we may state: Proposition:Openings and closing by reconstruction constitute two complete lattices for the usual sup and inf . Closing by Reconstruction ; Lattices

14. Application: Filtering by Erosion-Recontruction • Firstly, the erosion X,Blsuppresses the connected components of X that cannot contain a disc of radius l; • then the opening grec(X ; Y) of marker Y = X,Bl «re-builts» all the others. c) Reconstruction of b) inside a) a) Initial image b) Eroded of a) by a disc

15. Application: Removal of the edge grains • Let Z be the set of the edges , and X be the grains under study ; • Set Y is the reconstruction of Z€X inside set X ; • the set difference between X and Y extracts the internal particles. X Z a) Initial image b) Particles that hit the edges c) Residue a) - b)

16. Application: Holes Filling Comment :efficient algorithm, except for the particles that hit the edges of the field. initial image X reconstruction of A inside XC A = part of the edge that hits XC

17. Individual Analysis of Particles X Algorithm (J.C. Klein) While set X is not empty do { ~ p := first point of the video scan; ~ Y := connected component of X reconstructed from p; ~ Processing of Y (and various measurments) ; ~ X := X \ Y } Y new X Individual extraction of particles

18. The adjunction opening and closing create a simpler function than the original. They smooth in a nonlinear way. The opening (closing) removes positive (negative) peaks that are thinner than the structuring element. The opening (closing) remains below (above) the original function. Effects on Functions

19. Definition : The top-hat oper tor is the residual between the initial function and its opening. Therefore it is ³ 0 Functions The top-hat is used to extract contrasted components with respect to the background. The basic top-hat extracts positive components and the dual top hat the negative ones. Typically, top-hats remove the slow trends, and thus performs a contrast enhancement. Top-hat operator

20. An Example of Top-hat Comment :The goal is the extraction of the aneurisms ( the small white spots). Top hat c), better than b) is far from being perfect. Here opening by reconstruction yields a correct solution (VI-16). Top hat by an hexagon opening of size 10. Negative image of the retina. Top hat by the sup of three segments openings of size 10.

21. In signal processing, the word "filters" is somewhat fuzzy and its meaning is context dependent. It may have the connotation of a convolution, or may denote any system which transforms one image into another. In mathematical morphology, the word has the following precise meaning: Definition (G. Matheron, J. Serra): A morphological filteris an increasing and idempotent transformation of a complete lattice into itself. Increasingness This requirement is the most fundamental. It ensures the preservation of the order relation, after filtering. This property implies some information loss during the filtering process. Definition of a Morphological Filter (I)

22. Idempotence By definition, an idempotent operator transforms any original signal, or image, into an invariant one for the said operator (Ch. I). This property appears often, but in an implicit way, in descriptions. An optical filter will be qualified of «red», an amplifier of «band-pass», for example. In Morphological filtering, the property becomes explicit, and compulsory. Idempotence may be reached after a single pass, or as as limit under iteration. It may concern a unique operation, or a whole sequence. Finally, note that when linear filters are idempotent, then they do not admit inverses : they loose information, which make them a bit more «morphological». Definition of a Morphological Filter (II)

23. Let z and y be two ordered filters (withz³y). By combining them, one can generate four filters which satisfy the following properties: Theorem (G.Matheron): 1) zy, yz, zyz, yzyare morphological filters (i.e. increasing and idempotent mappings); 2) 3) zyzis the smaller filter larger than zy Ú yz yzy is the larger filter smaller than zy Ù yz 4) We have the following equivalences zyz = yz Û yzy = zy Û yz³zy N.B.: 1-The composition of more than three operators does not give another filter. Since, as a filter, zy is idempotent, we havezyzy = zy. zy £ yz Alternating Filters y £ yzy £ zyz £ z

24. The more often, the two basic primitives used in the composition are an opening g and a closing j . Then the previous theorem states that, in general, there is no order relation between gj and jg: There is a "classical" case which illustrates point 4) of the previous theorem, namely that of the filters by reconstruction (Ch VI-21). Then the products gj and jg are ordered, with gj³ jg Comparison between gj and jg

25. Flat and Increasing Connected Operators • From now on, we focus exclusively i ) on the criterion s of flat zones ; ii ) and on those operators y : TEDTEthat are flatand increasing. Basic Properties : • Every binary connected (resp. and increasing) operator induces on TE, via the cross sections, a unique connected (resp. and increasing) operator ( H. Heijmans ) ; • In particular, the geodesic implementations extend to the numerical case ; • The properties of the set case, to be strong filters, to constitute semi-groups, etc.. are transmitted to the connected operators induced on TE. Note that a mapping may be anti-extensive on TE, but extensive on the latticeg of the partitions (e.g.reconstruction openings).

26. Numerical Geodesic Dilations (I) numerical geodesic dilation of g with respect to f • Let f and g be two numerical functions from Åd into T, with g £ f. The binary geodesic dilation of size l of each cross section of g inside that of f at the same level induces on g a dilation df,l(g) (S.Beucher) . • Equivalently, (L.Vincent) the sub-graph of df,l(g) is the set of those points of the sub-graph of f which are linked to that of g by - a non descending path - of length £ l. f df,l(g) g

27. Numerical Reconstruction Numerical Reconstruction of g inside f : • The reconstruction opening of f from g is the supremum of the geodesic dilations of g inside f, this sup being considered as a function of f: grec(f ; g) = ƒ{ df,l(g) , l>0 } The dual closing for the negative isjrec(f ; g) = m - grec(m- f ; m- g) • The three major applications are : -swamping, or reconstruction of a function by imposing markers for the maxima; - reconstruction from an erosion - contrast opening ,which extracts and filters the maxima. f grec(f ; g) g

28. Application to Retina Examination Comment:The aim is to extract and to localise aneurisms. Reconstruction operators ensure that one can remove exclusively the small and isolated peaks ( case study due to F. Zana and J.C.Klein ). b) closing by dilatation-reconstruction followed by opening by érosion- reconstruction a) Initial image c) difference a) minus b) followed by a threshold

29. An Example of Top-hat Comment :The goal is the extraction of the aneurisms ( the small white spots). Top hat c), better than b) is far from being perfect. Here opening by reconstruction yields a correct solution (VI-16). Top hat by an hexagon opening of size 10. Negative image of the retina. Top hat by the sup of three segments openings of size 10.

30. Goal To remove the useless maxima (or minima) of a function. Algorithm The "marker" is a bi-valued (0,m) function identifying the peaks of interest. The reconstruction process result is the largest function £ f and admitting maxima at the marked points only. It is called the swamping of f by opening (S.Beucher, F.Meyer), Reconstruction of a Function from Markers Swamping of f by markers m ( by opening )

31. Goal Both morphological and reconstruction openings reduce the functions according to size criteria which work on their cross sections. In opening by dynamics, the criterion holds on grey tones contrast (M.Grimaud). An Example of Swamping : Contrast Opening • Algorithm • Shift down the initial function f by constant c; • Rebuilt f from function f - c, i.e. • grec(f ; f-c) = ƒ{ df(n)(f-c) , n > 0 } • The associated top-hat extracts all peaks of dynamics ³ c

32. Application to Maxima Detection • The maxima of a numerical function on a space E are the connected components of E where f is constant and surrounded by lower values. • Therefore they are given by the residues of contrast opening of shift c = 1 • More generally, the residuals associated with a shift c extract the maxima surrounded by a descending zone deeper than c. They are called Extended Maxima (S. Beucher).

33. An Example of a Pyramid of Connected A.S.F.'s Flat zones connectivity, (i.e. j = 0 ). Each contour is preserved or suppressed, but never deformed : the initial partition increases under the successive filterings, which are strong and form a semi-group. ASF of size 8 ASF of size 4 ( hexagonal structuring elements) ASF of size 1 Initial Image

34. Given marker M, consider AÎs(E) . Let ~ gM(A) be the union of the grains of A that hit M or that are adjacent to it ; ~ jM (A) be the union of A and of its pores that are included in M and non adjacent to Mc Definition (F. Meyer) : leveling l is defined as theactivity supremum l = gMÿjM i.e. l(A) € A = gM€ A, and l(A) € Ac = jM€ Ac. Leveling lacts as opening gM inside A , and as closing jM inside Ac . Self-duality: The mapping (A,M)Dl(A,M) from s(E)´s(E)Ds(E) is self-dual. If M itself depends on A, i.e. if M = m(A), then l, as a function of A only, is self-dual iff m is already self-dual. The extension to functions (via their cross-sections) will be denoted by (f, g) DL (f, g) Levelling

35. Here are a few nice properties of levelling : Proposition (F.Meyer): The levelling (A,M)Dl(A,M) is an increasing mapping from s(E) ´ s(E)Ds(E); it admits the equivalent expression: l = gM~ ( ¥€jM) Proposition (G.Matheron):The two mappings ADlM(A) , given M, and MDlA(M) , given A, are idempotent (hence are connected filters on s(E) ). Proposition (J. Serra): Levelling ADlM(A) is a strong filter, and is equal to the commutative product of its two primitives l = gM 1jM = jM 1gM iff connection f is adjacency preventing.Then, l satisfies the stability relation : gx( I~ l ) = gx ~ gx l , i.e. preserves the sense of variation at the grains/pores junctions . Properties of Levelling

36. Initial image : « Joueur de fifre », by E. MANET Markers : hexagonal alternated filters, (non self-dual) An Example Marker j 1 g 2 flat zones : 53.813 Marker g1 j 2 flat zones : 53.858 Initial image, 83.776 pp flat zones : 34.835

37. If 0 and m stand for the two extreme bounds of the gray axis T, then the set complement operation is replaced by its function analogue f® m - f and we have for levelling L m - L ( m - f , m - g ) = L ( f, g ) (1) which means that f, g ® L(f, g) is always a self dual mapping. In addition, if g derives from f by a self-dual operation, i.e. g = g(f) with m -g( m - f ) = g( f ) (2) (e.g.convolution, median element), then levelling f ® L(f, g(f)) is self-dual. Observe that rel. (2) is distinct from that of invariance under complement g( m - f ) = g( f ) which is satisfied by the module of the gradient, or by the extended extrema, for example, and which does not imply self-duality for f ® L(f). Duality for Functions

38. Marker: extrema with a dynamics 8 h ( invariance under complement). An Example of Duality Initial image flat zones : 34.835 h = 80 flat zones : 57.445

39. We now fix set A and study the mapping MDlA(M) as marker M varies. Set A generates on s(E) the A-activity ordering 9A by the relations M19A M2 i.e. if M1 meets A or is adjacent to A, then M2 meets A or is adjacent to A and if M2 meets Ac or is adjacent to Ac, then M1 meets Ac or is adjacent to Ac. Proposition (J. Serra): If M19A M2 , then we have llA(M1)(M2) = llA(M2)(M1) = lA(M2) In practice, this granulometric pyramid allows to grade markers activities . Levelling as function of the marker M1 M2 A

40. Marker: Initialimage, where the h-extrema are given value zero (self-dual marker) An Example of Pyramid Initial image flat zones : 34.835 Levelling for h = 50 flat zones : 58.158 Levelling for h = 80 flat zones : 59.178

41. Marker: Gaussian convolution of size 5 of the noisy image An Example of Noise Reduction a) Initial image, plus 10.000 noise points b) Gaussian convolutionof a) c) Levelling of a) by b) flat zones : 46.900

42. Granulometry is the study of the size characteristics of sets and of functions. In physics, granulometries are generally based on sievesyl of increasing meshes l > 0. Now, by applying sieve l to set X, we obtain the over-sieve yl(X) Í X; if Y is another set containing X, the Y-over-sieve, for every l, is larger than the X-over-sieve, i.e.X Í Y Þ yl(X) Í yl(Y) ; if we compare two different meshes l and m such that l ³ m, the m-over-sieve is larger than the l-over-sieve, i.e. l ³ m Þ yl(X) Í ym(X) finally, by applying the largest mesh l to the m-over-sieve, we obtain again the l-over-sieve itself, i.e. yl ym(X) = ym yl (X) = yl (X) Such a description of the physical sieving suggests to resort to openings for an adequate formalism. The sizes of the structuring elements will play the role of the the sieves meshes. Granulometry: an Intuitive Approach

43. Matheron Axiomaticsdefines a granulometryas a family {gl} i) of openings depending on a positive parameter l, ii) and which decrease as l increases: l ³ m > 0 Þ gl £ gm. This second axiom is equivalent to the semi-groupwhere the composition of two operations is equal to the stronger one, namely gl gm= gm gl = gsup(l,m) (1) [{ l ³ m > 0 Þ gl £ gm } Þ gl = gl gl £ (gl gm Ú gm gl) £ gl ; conversely, gl = gm gl etgl £ I Þ gl £ gm hence semi-group (1).] If BletBmstand for the invariant elements of gl and of gm respectively, then we easily see that (1) Û BlÍ Bm Finally, when gl's are adjunction openings, i.e.gl(X)=XolB (withsimilar structuring elements), then the granulometry axioms are fulfilled if and only if the structuring element B is compact and convex. Granulometry: a Formal Approach

44. An Example of Granulometry By duality, the families of closings {jl , l>0} increasing in l generate anti-granulometries of law jl jm= jm jl = jsup(l,m) . Here, from left to right, closings by increasing discs.

45. On Adjunction Opening, and on Granulometry : The notion of an adjunction opening (or closing) and that of a granulometry were introduced by G.Matheron in 1967 {MAT67}, for the binary and translation invariant case. In 1975, {MAT75} linked them with algebraic closings in E.H.Moore’s sense (begining of the XXth century) and generalized the notion of a granulometry, but still in the set case.The extension to lattices, with the representation theorem, is due to J. Serra {SER88. For details and examples about granulometric spectrum, see {HUN75},{SER82,ch.10},{MAR87a}.For algebraic openings, see{RON93}. On Gradient and Top hat mappings : Morphological gradient appears for the first time with S.Beucher and F.Meyer in 1977 {BEU77} and top hat operator with F.Meyer in 1977 {MEY77}, followed by a reformulation in terms of function mappings in {SER82,ch.12}. More sophisticated versions of the gradient can be found in {BEU90}, and top hat mappings with non flat structuring elements are proposed in {STE83}, (see also {SER88, ch9}). References