Image Anal Stereol 2015;34:111-123 Original Research Paper doi: 10.5566/ias.1291 MATHEMATICAL MORPHOLOGY BASED CHARACTERIZATION OF BINARY IMAGE Raghvendra Sharma^ and B. S. Daya Sagar Systems Science and Informatics Unit, Indian Statistical Institute-Bangalore Centre, 8th Mile, Mysore Road, RV College PO, Bangalore-560059, India. e-mails: rsharma@isibang.ac.in; bsdsagar@isibang.ac.in (Received February 19, 2015; revised May 25, 2015; accepted June 10, 2015) ABSTRACT This paper reports the results of a theoretical study on morphological characterization of foreground (X) and background (Xc) of a discrete binary image. Erosion asymmetry and dilation asymmetry, defined to elaborate smoothing of an image respectively by contraction and expansion, are generalized for multiscale smoothing, and their relationships with morphological skeleton and ridge (background skeleton) transformations are discussed. Then we develop algorithms identifying image topology in terms of critical scales corresponding to close-hulls and open-skulls, along with a few other salient characteristics, as respective smoothing by expansion and contraction proceeds. For empirical demonstration of these algorithms, essentially to unravel the hidden characteristics of topological and geometrical relevance, we considered deterministic and random binary Koch quadric fractals. A shape-size based zonal quantization technique for image and its background is introduced as analytical outcome of these algorithms. The ideas presented and demonstrated on binary fractals could be easily extended to the grayscale images and fractals. Keywords: dilation asymmetry, erosion asymmetry, close-hull, open-skull, degree of stability, hull fragments, skull fragments. INTRODUCTION The two fields mathematical morphology (Matheron, 1975; Serra, 1982) and fractal geometry (Mandelbrot, 1982) evolved independently but almost in same era near 1960s. Fractal geometry offers computation of fractal dimension, which is one of the measures to quantify the degree of roughness, and has been widely employed within the context of image processing, in texture-based image classification, image compressions, and many (Kaplan, 1999; Tolle et al., 2003; Xia et al., 2006; Ji et al., 2013). However, the measures that quantify the shape-size content in varied types of images such as binary, grayscale, color and hyper-spectral images; are from the field of mathematical morphology. Interestingly, various mathematical morphological transformations are rightly appropriate to characterize fractal objects, fractal functions, fractal surfaces etc. This paper emphasizes on characterization of spatial binary objects (e.g. fractals) not by fractal dimensions but by morphological analysis. Mathematical Morphology evolved as set-theory based image analysis approach. The central idea of mathematical morphology is to examine geometrical structure of an image object with reference to an object of simple shape and size, termed as structuring element. The interactions of image object with structuring elements result in nonlinear smoothing filters, also known as morphologic transforms (Serra, 1982; Maragos and Schafer, 1987a,b). Dilation and erosion are basic and mutually dual mathematical morphologic transforms. Dilation of an image followed by erosion is closing transform, and dual of this, i.e., erosion of an image followed by dilation is opening transform. Skeleton and ridge (henceforth referred to denote background skeleton) transforms, which are dual to each other, are geometrical representations summarizing overall shape, abstract structure, orientation of foreground and background of an object respectively. Skeleton of a continuous binary image was first introduced by Blum (1967), as medial axis. Skeleto-nization is based upon a systematic use of multiscale erosions and openings with certain logical operations (Lantuejoul, 1980), while the ridge follows multiscale dilations and closings of an object. Shape-size analysis plays prime role in object recognition, and is an area of utmost importance in digital image processing and computer vision. Shape-size analysis of image by multiscale smoothing has been used in mathematical morphology since its inception. 111 Sharma S et al: Morphological characteristics of binary image Granulometry and antigranulometry (Matheron, 1975; Serra, 1982; Dougherty et al., 1989; Heijmans, 1994; Soille, 2003), primary tools of morphology, systematically apply morphological transformations and their cascade-operations for multiscale shape representation, and define shape-size content based quantitative indexes to compute spatial complexity of foreground and background regions of a set (e.g. binary fractal) or a function (e.g. grayscale image). Maragos (1989) developed shape-size descriptor called pattern spectrum to detect critical scales in an image object and to quantify various aspect of shape-size content of the image. Several researchers employed binary morphological transformations to decompose spatial objects (e.g. binary images) into non-overlapping disks (convex polygons) of various sizes, and also into skeletons (Maragos, 1989; Pitas and Venet-sanopoulos, 1990, 1992; Goutsias and Schonfeld, 1991; Reinhardt and Higgins, 1996a,b; Xu, 1996, 2001a,b,c, 2003a,b). These studies were carried out essentially for representation, and recognition of binary patterns. The main focus of the aforementioned papers was to decompose an image object based upon certain criterion e.g. homothetics of used structuring element or collection of convex polygonal etc. The current work is motivated by this seed idea of decomposition of image, but adapts an altogether different approach driven by internal image characteristics to achieve the same. This very characteristic primarily is motivated by the work of Beucher (2005) on morphological residues, and could be identified by critical scales obtained from unique topological networks of image, namely skeleton and ridge, during multiscale analysis. The focus of this paper is not only limited to shape-size based systematic zonal quantization of image and its background, but it also provides spatial modeling of image in terms of shape-size topology, texture analysis and pattern classification. In this paper, we demonstrate few aspects related to multiscale smoothing based shape-size complexity of image with primary focus only on finite discrete binary images and finite discrete binary structuring elements containing origin. Multiscale opening and closing are employed respectively for the analyses of foreground and background of a spatial object (e.g. binary fractal). It is known that the multiscale openings are scaled up to ultimate erosion of the image (Beucher, 2005, 1994), while the multiscale closings typically are scaled such that any further closing by a bigger size structuring element yields no expansion in closed version (assuming finite image and enough big background). In our experimentation during multiscale analysis of discrete image objects, we have observed that image (and its background) may have no roughness with reference to a specific discrete size. This 'no roughness' means a minimum (i.e., zero) value of corresponding quantitative indexes or pattern spectrum as referred above. In some cases, this 'no roughness' pattern stabilizes for a few suc-cessive multiscale iterations, which could be referred as stabilization of residues (Beucher, 2005). The overall analysis and results shown provide few insights on image topology that were not reported elsewhere to best of our knowledge. This topographic characteristic provides a base for shape-size based quantification of foreground and background of image. This paper is organized as follows: We define 'dilation asymmetry' and 'erosion asymmetry' in the next section. These terms are not mentioned in classical mathematical morphology and are defined in this paper to elaborate smoothing of an image by expansion and contraction respectively. After this, we generalize these terms for multiscale smoothing and show their oneness respectively with ridge and skeleton. This is followed by a discussion on shape-size complexities of image background and foreground along with their 'degree of stability', respectively via multiple close-hulls and open-skulls, under the influence of increasing cycles of morphological closing and opening. We also propose algorithms to obtain critical scales corresponding to multiple close-hulls and open-skulls of an image along with their degree of stability. Afterwards, we demonstrate experimental results of these algorithms on deterministic and random binary Koch quadric fractals. As an outcome of analysis of these algorithms, we define systematic quantization techniques to quantify image and its background into zonal fragments with open-skulls and close-hulls respectively being the quantifiers. There exist various techniques to estimate morphology based fractal dimension (e.g. Maragos and Sun, 1993; Radhakrishnan et al., 2004). We compute a scale invariant but shape-dependent morphological quantitative index and correlate this with analytical fractal dimension of binary fractal. The paper ends with conclusion of the results of our study. Notation: Z = set of all integers; Z2 = two dimensional grid of discrete points; X, A, B, C = subsets of Z2; A c C = set A is subset of set C; A - C = set difference between A and C; U = set union; X©B = dilation of X by B; X0B = erosion of X by B; X • B = closing of X by B; X o B = opening of X by B; 0 = empty set. 112 Image Anal Stereol 2015;34:101-110 MORPHOLOGICAL SMOOTHING - ASYMMETRY AND SYMMETRY OF IMAGE Morphological dilation and erosion transforms of an image X are actually defined against symmetric set of structuring element B with respect to origin (Matheron, 1975; Serra, 1982; Maragos and Schafer, 1986). For the ease of notation, but with the same spirit of definitions, we refer simply X©B and X0B respectively as morphological dilation and erosion of X by B. The cascade of erosion-dilation (resp. the cascade of dilation-erosion), in other words X o B (resp. X • B), implemented by an arbitrary size of B, provide morphological opening and closing. Closing is an extensive transform while opening is an antiextensive transform (Serra, 1982), i.e., for finite discrete binary image X and for finite discrete binary structuring element B, X o B c X c X • B. An image smoothes in closing transform such that the entire asymmetry of image with reference to struc-turing element is removed at cost of its expansion. Similarly, image smoothes in opening and entire asymmetry of image with reference to structuring element is removed at cost of its contraction. The terms 'smoothing by expansion' and 'smoothing by contraction' are used later in the document with this essence. Further in this context, we define dilation and erosion asymmetry in coming sub-sections. DILATION ASYMMETRY AND SYMMETRY We define dilation asymmetry Dasym(X), (or simply Dasym as we refer it further) of X with reference to B as: Dasym(X) = (X • B) - X. (1a) Dilation asymmetry is a measure of asymmetry or distortion (by expansion) of X by B due to closing. Dasym equals to 0 in Eq. 1a imply X • B = X, which means X has no dilation asymmetry, or is perfectly dilation symmetric to B and cannot be expanded by B in closing cycle(s). We define dilation symmetry Dsym(X), (or simply Dsym as referred further) of X with reference to B as: Dsym(X) = (X ®B) - (X • B). (1b) Dilation symmetry is a measure of symmetry of X with reference to B during the closing transform. For case of perfect dilation symmetry, i.e., X• B = X, Dsym = (X©B) - X, which implies that the entire incremented part of X due to dilation is symmetric with reference to B and hence, completely shrinks back during erosion in closing cycle. Interestingly, Eqs. 1a,b jointly illustrate that Dasym is incremented part of X due to dilation, which does not shrink back after erosion in closing cycle; while Dsym is incremented part of X due to dilation, which shrinks back after erosion as closing cycle completes. This understanding elaborates the obvious observation from Eqs. 1a,b that Dasym and Dsym are mutually disjoint, and: Dasym(X) U Dsym(X) = (X ©B) - X. (1c) Eq. 1c indicates that set union of Dasym and Dsym is a zonal ring between X©B and X. We stated earlier that entire dilation asymmetry of X with reference to B is removed at cost of its expansion in closing, which means that closed version X • B is perfectly dilation symmetric to B, i.e., Dasym(X • B) = 0. This, with replacement of X by X • B in (1a) yields (X • B) • B = X • B; concluding that closing of an image is an idempotent transformation. This also infers that dilation of X • B by B results in X©B, i.e., (X• B)© B = X©B, or (X©B) o B = X©B; which indicates that dilated version of an image is super stable with respect to openings. These properties were classically established by Serra (1982), and our purpose to refer them here is to establish a consistency of definitions in Eqs. 1a,b with them. Fig. 1 illustrates smoothing of image in closing transform and also summarizes the above mentioned properties of dilation asymmetry and dilation symmetry. Fig. 1. Image smoothing in closing transform along with properties of dilation asymmetry and dilation symmetry. EROSION ASYMMETRY AND SYMMETRY We define erosion asymmetry Easym(X), (or simply Easym as referred further) of X with reference to B as: Easym(X) = X - (X o B). (2a) Erosion asymmetry is a measure of asymmetry or distortion (by contraction) of X by B due to opening. 113 Sharma S et al: Morphological characteristics of binary image Easym equals to 0 in Eq. 2a imply X= X o B, which means X has no erosion asymmetry, or is perfectly erosion symmetric to B and cannot be contracted by B in opening cycle(s). We define erosion symmetry Esym(X), (or simply Esym as referred further) of X with reference to B as: Esym(X) = (X o B) - (X e B). (2b) Erosion symmetry is a measure of symmetry of X with reference to B during the opening transform. For case of perfect erosion symmetry, i.e., X = XoB, Esym = X - (X0B), which implies that the entire reduced part of X due to erosion is symmetric with reference to B and hence, completely expands back during dilation in opening cycle. Eqs. 2a,b illustrate that Easym is reduced part of X due to erosion, which does not expand back after dilation in opening cycle; while Esym is reduced part of X due to erosion, which expands back after dilation as opening cycle completes. This elaborates the otherwise obvious conclusion from Eqs. 2a,b that Easym and Esym are mutually disjoint, and: Easym(X) U Esym(X) = X - (X e B). (2c) Set union of Easym and Esym is a zonal ring between X and X0B, as is clear from Eq. 2c. Dilation and erosion (as well closing and opening) are mutually dual transforms. The notion of previous sub-section could easily be correlated to procure Easym(X o B) = 0, followed by (X o B) o B = X o B; deducing that opening of an image is an idempotent transform. Similarly, we procure (X o B)0B = X0B, or (X0B) • B = X0B; which points that eroded version of an image is super stable with respect to closings. The latter procurement, and replacement of X by X0B in Eq. 1a results in Dasym(X0B) = 0. Similarly, as we know that (X©B) o B = X©B, replacing X by X©B in Eq. 2a results in Easym(X©B) = 0. In summary, eroded version of an image has no dilation asymmetry, and dilated version of an image has no erosion asymmetry. These properties were originally esta-blished by Serra (1982), and our purpose to refer them here is to establish a consistency of definitions in Eqs. 2a,b with them. Fig. 2 illustrates smoothing of image in opening transform and also summarizes the above mentioned properties of erosion asymmetry and symmetry. Fig. 2. Image smoothing in opening transform along with properties of erosion asymmetry and erosion symmetry. MULTISCALE MORPHOLOGICAL SMOOTHING DILATION ASYMMETRY AND SYMMETRY - GENERIC FRAMEWORK In previous section, we discussed that X completely smoothes in closing transform, by removal of its total asymmetry with reference to B at cost of its expansion. The scope of this paper is discrete binary space, and within this scope, further smoothing of X by expansion is feasible via its closing by bigger discrete size structuring elements belonging to family of B. By convention, base structuring element B is a discrete binary prototype pattern of unit (i.e., one) size, and a finite pattern nB, given by Eq. 3, defines a family of binary patterns belonging to B and specified by discrete size n. nB = 1 ©5©...©B. (3) where n is non negative integer. By convention, nB = {(0,0)} for n = 0. The cascade of erosion-dilation, and the cascade of dilation-erosion, implemented by nB for successive values of n, provides multiscale morphological opening and closing respectively. These multi-scale openings and closings are respectively denoted by X o nB and X • nB. A closing by nB removes entire asymmetry of and up to size nB out of X, i.e., (X• mB) c (X• nB) for any non-negative integers m, n such that n > m. Now as we perform multiscale smoothing of X by expansion, we are interested in a generalized version of Eq. 1a, i.e., dilation asymmetry Dnasym(X) for nth iteration of multiscale smoothing by expansion. The most obvious generalized version seems to be Dnasym(X) = X • nB - X, where n is a positive integer. n-times 114 Image Anal Stereol 2015;34:101-110 But with this pattern of generalization, we land up in a situation where Dpasym(X) E Dqasym(X) V q > p where p, q are positive integers. Rather than this approach of generalization, a much preferable idea to analyze the smoothing pattern of image will be to have smaller disjoint dilation asymmetry sets as iterations of multiscale smoothing proceed. To visualize disjoint dilation asymmetry sets, we consider topology of zonal ring between X and X©B, given by Eq. 1c, as elementary theme. Also, Fig. 1 depicts flip and flop of component Dsym within this zonal ring, by dilation of X • B and erosion of X©B respectively. So, if we realize multiscale smoothing of X by expansion in terms of zonal rings between X©nB and X©(n+1)B where n is a non negative integer, we shall obtain disjoint patterns of dilation asym-metry and dilation symmetry as smoothing proceeds, because a pair of unique dilation asymmetry and symmetry set shall exist within the region of each (mutually disjoint) zonal ring. Here we assume empty sets of dilation asymmetry as disjoint sets, if so happens for few iterations of multiscale smoothing. To generalize Eqs. 1a,b with this approach of zonal rings, we need to consider dilated image of current iteration of multiscale smoothing as input image for next iteration. We denote iteration count by nonnegative integer n and input image for nth iteration of smoothing by Xn. For n = 0, which is starting iteration, the input image X0 is X itself. As this approach considers dilated version of current iteration being input for next iteration, Xn = X©nB. With this approach of zonal rings, the generalized versions Dnasym(X) and Dnsym(X) of dilation asymmetry and symmetry respectively for nth iteration of multiscale smoothing by expansion are: Dnasym(X) = (X © nB) • B - (X © nB). (4a) Dnsym(X) = (X © (n+1)B) - (X © nB) • B. (4b) Note that Eqs. 4a,b are trivial to achieve with replacement of X by Xn in Eqs. 1a,b. Assuming finite image and enough big background, the upper threshold of n, denoted by K, is maximum (or last) iteration in which smoothing of X by expansion occurs; indicating that Dnasym is an empty set for each n > K. If so, using (4a) we obtain X • (n+1)B = X • nB V n > K. In other words, X • (K+1)B is maximum possible closed version of X obtained by multiscale smoothing. We will shortly discuss a theoretical limit on maximum iteration of multiscale smoothing. Analogous to previous section, Dnasym(X) (or simply Dnasym as we refer it further) is measure of asymmetry or distortion (by expansion) of Xn by B due to closing. Similarly, Dnsym(X) (or simply Dnsym as referred further) is measure of symmetry of Xn with reference to B during the closing transform. Dnasym and Dnsym are mutually disjoint as is obvious from Eqs. 4a,b; and set union of Dnasym and Dnsym is a zonal ring between Xn©B and Xn as indicated by Eq. 4c. Dnasym(X) U Dnsym(X) = (X © (n+1)B) - (X © nB). (4c) Eq. 4a in fact replicates the morphological ridge transform of discrete binary image X, which is realized in this sub-section as a asymmetry of topological zonal rings expanding out of X. Henceforth, the usage of term dilation asymmetry in essence could be pondered as morphological ridge (background skeleton). Maragos (1989) proposed extended reduced skeleton transform (ERST) which may eliminate some ridge redundancy in few scenarios and is compatible with pattern spectrum. The upper threshold value K obtained from ERST based model may slightly differ in some scenarios (e.g. a reduction of one iteration) than the one obtained from ridge based model (Eq. 4a) as described above. As our focus in this paper is more towards classifying overall shape-size topology of image, we embrace with ridge based model with awareness of this slight effect. However, we will point out this possible minor impact, wherever applicable. The generalized version of results summarized in Fig. 1 of previous section could be trivially obtained for nth zonal ring with replacement of X by Xn, i.e., (Xn • B) • B = Xn • B, and (Xn©B) o B = Xn©B. Replacing Xn by X©nB, the generalized results are ((X©nB) • B) • B = (X©nB) • B, and (X©(n+1)B) o B = X©(n+1)B. As a closing by nB removes entire asymmetry of and up to size nB out of X, a more generic flavor of results is: (X • nB) • mB = (X • nB) V m < n, (5a) (X © nB) o mB = (X © nB) V m < n, (5b) where m, n are non-negative integers. We come back to pending discussion on theoretical limit for maximum iteration of smoothing by expansion. A cut-off cap Kmax on iterations could be introduced because of the fact that multiscale smoothing of X by expansion definitely stops as soon as X becomes subset of a specific discrete size binary pattern belonging to family of B as described in Eq. 3. The outcome of this viewpoint is cut-off value given by Eq. 6. Kmax = minimum n | X (z(n+1)B, or Kmax = minimum n | { (n+1)B 0 X f 0}, (6) 115 Sharma S et al: Morphological characteristics of binary image where n is non-negative integer. Clearly, Dnasym = 0 for n > Kmax. The dealing of X against (n+1)B (and not nB) in Eq. 6 makes it consistent with model used in Eq. 4a, where iterations of smoothing start at n=0, resulting in smoothing of X by (n+1)B in nth iteration. Eq. 6 is primarily designed for both X and B to be 2D (Two-Dimensional), which is prime focus of our experimentation. However, it accommodates not so practical scenarios where X is 1-D (One-Dimensional) but B is 2-D. X and B both being 1-D is a trivial case. A scenario where X is 2-D and B is 1-D is also not very common. For such scenario, a minimum value of n such that (n+1)B supersedes longest anisotropic straight line which is subset of X, could be considered as Kmax. It is evident that K is strictly lesser than Kmax. The accurate value of K precisely depends upon geometric constitutions of X and B, and happens to be far less than Kmax in most of the practical cases. EROSION ASYMMETRY AND SYMMETRY - GENERIC FRAMEWORK Analogous to previous sub-section, the disjoint patterns of erosion asymmetry and erosion symmetry could be realized by pervasion of zonal ring topology between X and X0B (refer Eq. 2c) to Xn and Xn0B, where n is non-negative integer, and Xn (= X0nB) is generic input image for nth iteration of multiscale smoothing by contraction. The generalized versions Enasym(X) and Ensym(X) (or simply Enasym and Ensym as referred further) of erosion asymmetry and symmetry respectively derived with this approach are: Enasym(X) = (X e nB) - (X e nB) o B. (7a) Ensym(X) = (X e nB) o B - (X e (n+1)B). (7b) The iterations start from n = 0. X shrinks as iterations of multiscale smoothing by contraction proceed, and it is evident that maximum possible iteration Nmax is identified by ultimate erosion of image, i.e.: Nmax = n | {(X e nB ^ 0) * (X e (n+1)B = 0)}. (8) Observe from Eq. 7a and Eq. 8 that erosion asymmetry for n = Nmax is a non-empty set equals to (XeNmaxB). Now as we have identified Nmax, going by analogy of previous sub-section, can we think of a upper threshold on n, which is denoted by N and represent the maximum (or last) iteration in which smoothing of X by contraction occurs? The answer is that depending upon geometry of X and B, the ultimate erosion of X may occur well before its complete smoothing, and accordingly N may or may not exist. So, in some scenarios, we may obtain Nmax without obtaining N. However, N may exist in some cases, and could be conceptualized as an iteration such that input image Xn is perfectly erosion symmetric to B for ne (N, Nmax); i.e., Enasym is 0 in open interval (N, Nmax). X o (N+1)B is the smallest non-empty opened version of X obtained by multiscale smoothing, if N exists. Enasym and Ensym are mutually disjoint as is perceptible from Eqs. 7a,b; and their union is a zonal ring between XenB and Xe(n+1)B. Enasym(X) U Ensym(X) = (X e nB) - (X 9 (n+1)B). (7c) Conspicuously, Eq. 7a replicates morphological skeleton transform (Lantuejoul, 1980) of discrete binary image X, which is depicted here as asymmetry of shrinking topological zonal rings. Now onwards, the usage of term erosion asymmetry in essence could be contemplated as morphological skeleton. Maragos (1989) proposed reduced skeleton transform (RST) which may eliminate some skeleton redundancy in few cases depending on choice of X and B, and is compatible with pattern spectrum. The upper threshold value N (if it exists as explained above), obtained from RST based model may slightly differ in some scenarios than the one obtained from skeleton based model (Eq. 7a). We adopt skeleton based model in this paper on similar lines as explained in previous sub-section. Apparently, the generalized version of results summarized in Fig. 2 of previous section are ((XenB) o B) o B = (XenB) o B, and (Xe(n+1)B) • B = Xe(n+1)B. Eqs. 9a,b demonstrate more generic flavor, known as absorption property, on account of erudition that an opening by nB removes entire asymmetry of and up to size nB out of X. (X o nB) o mB = (X o nB) V m < n, (9a) (X e nB) • mB = (X e nB) V m < n, (9b) where m, n are non-negative integers. SHAPE SIZE TOPOLOGY -CLOSE HULLS AND OPEN SKULLS OF IMAGE CLOSE-HULLS IN BACKGROUND SMOOTHING Let us consider multiscale closing of X in close interval [0, Kmax], which in previous section was derived as potential integral range of iterations for multiscale smoothing of X by expansion. The closed version X • kB corresponding to iteration count k is regarded as a close-hull (or only hull as we often 116 Image Anal Stereol 2015;34:101-110 refer it further) of X if condition X • kB = X • (k+1)B is satisfied. This condition could be achieved by replacing Dnasym as 0 for n = k in Eq. 4a. Consequently, an empty dilation asymmetry set for an iteration of smoothing in Eq. 4a ensures that a close-hull exists corresponding to this iteration. Apparently, the mentioned condition for a hull could be potentially satisfied for multiple iterations in [0, Kmax], which implies that there may be multiple hulls of X. An interesting point of observation is that if a hull condition is satisfied for n>0; the corresponding input image X©nB in Eq. 4a for this iteration is stable with reference to opening and closing by B, i.e., (X©nB) oB = (X©nB) = (X©nB)• B. Now, we first define 'degree of stability' of a close-hull; and then will refine close-hull condition as per this definition. The degree of stability of a hull, which is obtained corresponding to iteration count k, is defined as number of successive iterations immediately after k for which dilation asymmetry does not appear. If degree of stability for this hull is s; then X • kB = X • (k+1+i)B, 0 < i < s. In other words, there exists null asymmetry pattern for discrete size interval [(k+1)B, (k+s+1)B]. Obviously, iterations corresponding to degree of stability do not stand up for a new hull; and to avoid any such confusion, we refine condition to achieve hull as follows. The closed version X • kB corresponding to iteration count k is regarded as a hull of X, if condition in Eq. 10 is satisfied. {X • kB = X • (k+1)Bj A {X • (k-1)B £X • kB}, (10) where k is a non negative integer in [0, Kmax] as we said earlier. There is a fair possibility that first hull of image is obtained at k = 0 itself, and condition (Eq. 10) is truncated to X = X • B for this case. Let's denote iteration count for first hull by Kj and its degree of stability by S1. In general, S1 could be much less than (Kmax - K1), and hence dilation asymmetry appears again at (K1+S1+1)th iteration. If reappeared dilation asymmetry persists continuously for T1 iterations before disappearing again, then iteration count K2 for second hull will be (K1+S1+T1+1). This pattern may repeat depending upon geometry of X and B as we iterate further till Kmax and multiple hulls of X could be obtained. Let M be a positive integer denoting maximum number of hulls of X, and Ki (0 < Ki < Kmax) be the iteration count for which ith hull (1 < i < M) is obtained as per condition in Eq. 10. So, X • KMB is maximum possible hull of X, i.e., X• nB = X• KMB V n > KM, and there will be no more smoothing of X by expansion in further iterations till Kmax (if KM < Kmax) or beyond. In other words, all iterations succeeding KM belong to degree of stability of maximum hull, and hence its degree of stability could be considered as infinite assuming image has huge background. In previous section, we denoted K as maximum (or last) iteration in which smoothing of X by expansion occurs, and hence K= KM-1. Iterations K1 to KM, corresponding to hulls, are measure of scale invariant but shape dependent characteristic of background of X. Iterations Ki, along with degree of stability of hulls, which we denote by Si for ith hull (1 < i < M), exhibit shape-size topology of background of X. Algorithm (1) below, which computes total number M of hulls along with Ki and Si (1 < i < M) values for these hulls, has a linear complexity. It requires image X, structuring element B and Kmax (obtained as per Eq. 6) as input parameters; and assumes that X has enough background to hold dilations till Kmax. As degree of stability of maximum hull theoretically is infinite, Alg. 1 refrains to compute the same. Let INPUT_IMAGE, DILATED_ IMAGE and CLOSED_IMAGE denote the sets indicating input, dilated and closed image respectively for a generic iteration during background smoothing. Let n, hull_count and degree_of_stability indicates generic iteration count, hull number and degree of stability of hull respectively. Let last_iteration_closed_input_image_same be a flag indicating if closed and input images for previous iteration were same or not. Algorithm (1): Computing Background Topology of Image Initialization: INPUT_IMAGE= X; n= 0; hull_count= 0; degree_of_stability = 0; last_iteration_closed_input_image_same= false; for n=0..Kmax do DILATED_IMAGE = INPUT_IMAGE0B; CLOSED_IMAGE = DILATED_IMAGE0B; if (CLOSED_IMAGE == INPUT_IMAGE) then if (last_iteration_closed_input_image_same == true) then // Increment degree of stability of current hull. degree_of_stability ++; else // This iteration corresponds to a new hull. hull_count ++; M = hull_count; Khull_count = n; // Set flag for upcoming iterations. last_iteration_closed_input_image_same = true; end //if(last_iteration_closed_input_image_same == true) else // Iteration indicates non-empty dilation asymmetry either // between image and first hull or successive hulls. if (last_iteration_closed_input_image_same == true) then //Store & reset degree of stability of current hull. 117 Sharma S et al: Morphological characteristics of binary image Shuii_count = degree_of_stability; degree_of_stability = 0; // Reset flag for upcoming iterations. last_iteration_closed_input_image_same = false; end //if(last_iteration_closed_input_image_same == true) end // if (CLOSED_IMAGE == INPUT_IMAGE) INPUT_IMAGE = DILATED_IMAGE; end // for Alg. 1 could be easily modified for ERST model (Maragos, 1989). This model, depending upon geometry, may have slight difference in some values of K and Si for few combinations of X and B. OPEN-SKULLS IN FOREGROUND SMOOTHING Applying cognitive resemblance to close-hull, an opened version X o nB corresponding to iteration count n of multiscale smoothing by contraction is conceptualized as an open-skull (or only skull) of X, if condition in Eq. 11 is satisfied. {X o nB = X o (n+1)B} A {X o (n-1)B £ X o nB}, (11) where n is a non negative integer in close interval [0, Nmax], which is range of iterations of multiscale smoothing of X by contraction as derived in previous section. Apparently, there may be multiple iterations in [0, Nmax] satisfying condition of Eq. 11, and hence there may be multiple skulls of X. The pruned version of Eq. 11, to monitor if (first) skull is obtained at n=0 itself, is X = Xo B. In logical correlation to degree of stability of a hull, the degree of stability of a skull obtained corresponding to iteration count n, is defined as number of successive iterations immediately after n for which erosion asymmetry does not appear. If we again denote s as degree of stability of this skull, then X o nB = X o (n+1+i)B, 0 < i < s. Condition in Eq. 11 ensures that iterations corresponding to degree of stability of existing skull are not interpreted as new skull. The model proposed in Eq. 7a together with check in Eq. 11 indicates that an iteration resulting in empty erosion asymmetry set corresponds to either a new skull, or degree of stability of existing skull. For such an iteration, a parity condition, i.e., (X0nB)• B = (X0nB) = (X0nB) o B is attained, provided n is strictly greater than 0; and the same is visible from Eq. 7a and Eq. 9b. As we have discussed skull in near resemblance to hull; it is good to point a contrasting fact that (at least one) hull exists for any finite image (assuming enough big image background), however there may be cases where an image has no skull depending upon the choice of structuring element. Having said that; we continue with generic case where X may have multiple skulls in [0, Nmax]. Let M be a positive integer denoting maximum number of skulls of X, and Ni (0 < Ni < Nmax) be the iteration count for which ith skull (1 < i < M) is obtained as per Eq. 11. So, X o NMB is ultimate (last) skull of X. In previous section, we conceptualized maximum iteration N of smoothing of X by contraction. For cases where N exists, N = Nm-1, and X o nB = X o NmB V ne [Nm, Nmax]. We know that erosion asymmetry always persists for n = Nmax ; and hence the last skull has finite degree of stability by definition. For cases where N exists, degree of stability of last skull is (Nmax-1-NM). Degree of stability of last skull is less than (Nmax-1-NM) for cases where N does not exists. Iterations Ni to NM are measure of scale invariant but shape dependent characteristic of foreground of X. These iterations along with degree of stability Si (1 < i < M) of skulls represent shape-size topology of foreground of X. Algorithm (2) computes total number M of skulls of image X, along with their Ni and Si (1 < i < M) values. Let INPUT_IMAGE, ERODED_IMAGE and OPENED_IMAGE denote the sets indicating input, eroded and opened image respectively for a generic iteration during foreground smoothing. Let n, skull_count and degree_of_stability indicates generic iteration count, skull number and degree of stability of skull respectively. Let last_iteration_open_input_ image_same be a flag indicating if opened and input images for previous iteration were same. Algorithm (2): Computing Foreground Topology of Image Initialization: INPUT_IMAGE= X; n=0; skull_count= 0; degree_of_stability = 0; last_iteration_open_input_image_same= false; while (1) do ERODED_IMAGE = INPUT_IMAGE0B; if (ERODED_IMAGE == 0 ) then Nmax = n; if (skull_count > 0) then if (degree_of_stability > 0) then // All iterations in (NM, Nmax) belong to stability of // last skull. Note that default value of stability is 0. Sskull_count = degree_of_stability; end; // (degree_of_stability > 0) end // if (skull_count > 0) break; // exit of algorithm end // if (ERODED_IMAGE == 0) OPENED_IMAGE = ERODED_IMAGE0B; if (OPENED_IMAGE == INPUT_IMAGE) then if (last_iteration_open_input_image_same == true) then // Increment degree of stability of current skull. degree_of_stability ++; else // A new skull. 118 Image Anal Stereol 2015;34:101-110 skull_count ++; M = skull_count; Nskull_count = n; // Set flag for upcoming iterations. last_iteration_open_input_image_same = true; end // if (last_iteration_open_input_image_same == true) else // Iteration indicates non-empty erosion asymmetry either // between image and first skull or successive skulls. if (last_iteration_open_input_image_same == true) then //Store & reset degree of stability of current skull. Sskull_count = degree_of_stability; degree_of_stability = 0; // Reset flag for upcoming iterations. last_iteration_open_input_image_same = false; end // if (last_iteration_open_input_image_same == true) end // if (OPENED_IMAGE == INPUT_IMAGE) INPUT_IMAGE = ERODED_IMAGE; n = n+1; end // while(1) Alg. 2 could be easily modified for RST model (Maragos, 1989), which may have slight difference in some values of N and Si in few cases. In essence, here we discussed multiple hulls and skulls of an image during multiscale analysis. Critical scales in terms of iterations K and N corresponding to hulls and skulls are entrenched, and concept of a locally stable or silent zone, where hulls and skulls stabilize for certain scale-interval is introduced. This literally means that scale-intervals [K; , K^S^ and [N N;+S;] in a typical size distribution or pattern spectrum plot reflect a straight line coinciding to axis representing size of closing and opening respectively. We also furnish algorithms to compute foreground and background topology of image in terms of these salient scales. EMPIRICAL ANALYSIS HULL AND SKULL BASED CHARACTERIZATION OF FRACTAL OBJECTS We applied Alg. 1 on deterministic and random binary Koch quadric fractals of Fig. 3. Table 1 shows the results with B as symmetrical flat structuring element of primitive size 3x3, containing 4 neighbors and rhombic in shape. For deterministic fractal: M = 6, Km = 62 and hence K (= KM-1) is 61. The number of successive iterations for which dilation asymmetry persisted e.g. between first and second hull of deterministic fractal is K2-(Kj+Si+1), i.e., 7 (iteration 12 to 18 both inclusive). We discussed earlier that depending upon geometric constitutions of X and B, the value of Kmax may be much higher than K, which exactly is reflected for deterministic and random fractals with Kmax value at 220 and 180 respectively. M-1 Observe that a considerable 29 (= £ S; + M) out of i=1 total 63(=Km+1) iterations in case of deterministic fractal corresponds to empty dilation asymmetry sets, but this number is not so high for random fractal. Fig. 3. (a) Deterministic and (b) random binary Koch quadric fractals (both 1024 x 1024 pixels). Applying Alg. 2 on deterministic and random binary Koch quadric fractals of Fig. 3 with exactly same rhombic structuring element yields the results as shown in Table 2. We obtained the Nmax value for deterministic and random fractal as 63 and 54 respectively; which implies that degree of stability of last skull is less than (Nmax-1-NM) for both these cases, and hence N does not exists. A simple example where N exists is circular image with exterior spikes on its boundary. The number of successive iterations for which erosion asymmetry persisted e.g. between first and second skull of random fractal is N2-(Ni+Si+1), i.e., 1 (iteration number 20). Again, notice M that 26 (= £ Si + M) out of total 64(=Nmax+1) i=1 iterations for case of deterministic fractal corresponds to empty erosion asymmetry sets, and this number is pretty low for random fractal. The ERST and RST model (Maragos, 1989) yields no impact on results of Table 1 and Table 2 respectively. Notice that the choice of rhombic structuring element, especially for deterministic fractal, establishes a very close equilibrium between the values obtained for skulls and their background counterparts, i.e., hulls, as is evident from Tables 1 and 2. For deterministic fractal of Fig. 3a, a flat square structuring element of primitive size 3x3 yields 10 hulls and 3 skulls, while a flat octagonal structuring element (primitive size 5x5) yields 3 hulls and 1 skull. For random fractal of Fig. 3b, square structuring element yields 2 close-hulls and 1 skull. In general, the model given by Alg. 1 and Alg. 2 could be applied to any image, resulting in one or more 119 Sharma S et al: Morphological characteristics of binary image hulls and zero or more skulls depending upon choice of X and B. Table 1. Background Shape Size Topology of (A) Deterministic and (B) Random Binary Koch Quadric Fractals of Fig. 3. JA) ÎBL Hull No. i Ki Si Hull No. i Ki Si 1 11 0 1 14 1 2 19 2 2 27 0 3 25 2 3 33 2 4 35 1 4 46 0 5 42 18 5 49 œ 6 62 œ Table 2. Foreground Shape Size Topology of (A) Deterministic and (B) Random Binary Koch Quadric Fractals of Fig. 3 (A) (B) Skull No. i Ni Si Skull No. i Ni Si 1 12 1 1 19 0 2 20 2 2 21 0 3 26 2 3 37 0 4 36 0 4 41 4 5 41 13 5 47 0 6 58 2 6 49 1 HULL AND SKULL BASED QUANTIFICATION OF IMAGE TEXTURE We know that granulomeres, and their discrete derivatives referred as size distribution (Matheron, 1975) or pattern spectra (Maragos, 1989), provide information regarding shape and size, and play important role in image texture analysis and pattern classification. Iterations corresponding to hulls and skulls reflect critical scales, and being these scales the quantifiers, we define hull and skull based discrete derivatives for zonal classification (of background and foreground respectively) of image. These derivatives are termed as Hull Fragment Pattern Spectrum (HFPS) and Skull Fragment Pattern Spectrum (SFPS) and are defined as: HFPSl(X) = J A* * Ki B) - A(X)'i = 1 A(X • Kfi) - A(X • Ki-1B),2 < i ' - • N / * * * (j) (k) (1) Fig. 6a-g skull fragments (in order from centre towards boundary of object), and (h) to (l) visible hull fragments (in order starting from object boundary and growing outwards) of deterministic fractal of Fig. 3a. Rhombus of primitive size 3x3 is used as structuring element. Max. (i.e. 6th) hull fragment is too small as is conspicuous from Fig. 4a and not shown. The hull and skull fragments based model could be utilized for image analysis, especially for images with systematic disorder. Fig. 6a-l reflects an ordered pattern in sequence of zonal fragments. This is due to systematic disorder in overall texture of deterministic fractal and will certainly not exist for all images (e.g. random fractal) under consideration. CONCLUSIONS This paper re-visits morphological ridge and skeleton transforms with a flavor of dilation and erosion asymmetry, followed by discussion on multiple hulls and skulls of an image. We provide a morphological 121 Sharma S et al: Morphological characteristics of binary image model to quantify image and its background based upon potential critical-scales, which demonstrate a sort of local stability in image texture. This local stability could be interpreted as invariance of image against certain scales during successive morphological opening and closing transforms. This phenomenon sometimes may persist consistently for a scale-interval. In summary, the proposed model comes handy in topographical analysis of image as well its background, and accordingly provides a tool for zonal-decomposition. The model exhibits scale invariant but shape dependent characteristics and is very effective for analysis of images having systematic disorder in overall texture. In particular, we can consider e.g. geospatial objects under this category where topographic characteristics stabilize locally for certain scale-interval, and proposed model could be handy to study the spatio-temporal stability of such objects. Though the paper demonstrate binary images as experimental prototype, all the concepts are applicable for gray tone images as well. There exists a scope for follow-up work in consideration of few aspects, e.g. a) possible usage of multiscale index (identified by ratio KM/NM, or preferably KM/Nmax) of an image for a given shape as a measure of background versus foreground complexity b) characterizing fractal dimension via morphological operations c) applying the work in this paper to grayscale domain. ACKNOWLEDGEMENTS Authors acknowledge the support received from Indian Statistical Institute under the project no: SSIU-030-2013-16. REFERENCES Beucher S (2005). Numerical residues. Mathematical Morphology: 40 Years On Computational Imaging and Vision 30:23-32. Beucher S (1994). Watershed, hierarchical segmentation and waterfall algorithm. Mathematical Morphology and its Applications to Image Processing, Fontainebleau, Jean Serra and Pierre Soille (Eds.), Kluwer Ac. Publ. Nld:69-76. Blum H (1967). A transformation for extracting new descriptors of shape in models for the perception of speech and visual forms. W Wathen-Dunn (Ed.), Cambridge, MA: M.I.T. Press. Dougherty ER, Kraus EJ, Pelz JB (1989). Image segmentation by local morphological granulomeres. Geoscience and Remote Sensing Symposium. Goutsias J, Schonfeld D (1991). Morphological representation of discrete and binary images. IEEE T Signal Proces 39:1369-79. Heijmans HJAM (1994). Morphological image operators. New York: Academic Press. Ji H, Yang X, Ling H, Xu Y (2013). Wavelet domain multifractal analysis for static and dynamic texture classification. IEEE T Image Process 22:286-99. Kaplan LM (1999). Extended fractal analysis for texture classification and segmentation. IEEE T Image Process 8:1572-85. Lantuejoul C (1980). Skeletonization in quantitative metallography. Issues of Digital Image Processing, Haralick RM and Simon JC (Eds.), The Netherlands: Sijthoff and Noordhoff. Mandelbrot BB (1982). Fractal geometry of nature. New York: W H Freeman. Maragos P (2005). Morphological filtering for image enhancement and feature detection. The Image and Video Processing Handbook, A C Bovik (Ed.). Elsevier Academic Press, 135-56. Maragos P, Schafer RW (1987a). Morphological filterspart I: their set-theoretic analysis and relations to linear shift-invariant filters. IEEE T Acoust Speech 35:115369. Maragos P, Schafer RW (1987b). Morphological filterspart II: their relations to median, order-statistic, and stack filters. IEEE T Acoust Speech 35:1170-84. Maragos P, Fang-Kuo Sun (1993). Measuring the fractal dimension of signals: morphological covers and iterative optimization. IEEE T Signal Proces 41:108-21. Maragos P, Schafer RW(1986). Morphological skeleton representation and coding of binary images. IEEE T Acoust Speech 34:1228-44. Maragos PA (1989). Pattern spectrum and multiscale shape representation. IEEE T Pattern Anal 11(7):701-16. Matheron G (1975). Random sets and integral geometry. New York: Wiley. Pitas I, Venetsanopoulos AN (1990). Morphological shape decomposition. IEEE T Pattern Anal 12:38-45. Pitas I, Venetsanopoulos AN (1992). Morphological shape representation. Pattern Recogn 25:555-65. Radhakrishnan P, Teo Lay Lian, Daya Sagar BS (2004). Estimation of fractal dimension through morphological decomposition. ELSEVIER Chaos Soliton Fract 21: 563-72. Reinhardt JM, Higgins WE (1996a). Efficient morphological shape representation. IEEE T Image Process 5:89101. Reinhardt JM, Higgins WE (1996b). Comparison between the morphological skeleton and morphological shape decomposition. IEEE T Pattern Anal 18(9):951-57. Serra J (1982). Image analysis and mathematical morphology. London: Academic Press. Soille P (2003). Morphological image analysis - principles and applications. Springer. 122 Image Anal Stereol 2015;34:101-110 Tolle CR, McJunkin TR, Gorsich DJ (2003). Suboptimal minimum cluster volume cover-based method for measuring fractal dimension. IEEE T Pattern Anal 25:32-41. Xia Y, Feng D, Zhao R (2006). Morphology-based multi-fractal estimation for texture segmentation. IEEE T Image Process 15:614-23. Xu J (1996). Morphological decomposition of 2-D binary shapes into conditionally maximal convex polygons. Pattern Recogn 29:1075-104. Xu J (2001a). Morphological decomposition of 2-D binary shapes into convex polygons: A heuristic algorithm. IEEE T Image Process 10:61-71. Xu J (2001b). Morphological representation of 2-D binary shapes using rectangular components. Pattern Recogn 34:277-86. Xu J (2001c). Efficient morphological shape representation with overlapping disk components. IEEE T Image Process 10:1346-56. Xu J (2003a). Efficient morphological shape representation by varying overlapping levels among representative disks. Pattern Recogn 36:429-37. Xu J (2003b). A generalized discrete morphological skeleton transform with multiple structuring elements for the extraction of structural shape components. IEEE T Image Process 12:1677-86. 123