Meijster, Arnold

Wilkinson, Michael H.F.

English

An area based counterpart of the binary structural opening spectra is developed It is shown that these area opening and closing spectra can be computed using an adaptation of Tarjan's union-find algorithm These spectra provide rotation, translation, and scale invariant pattern vectors for texture analysis.</p

Proceedings - University of Groningen

   University of GroningenFast computation of morphological area pattern spectraMeijster, Arnold; Wilkinson, Michael H.F.Published in:2001 INTERNATIONAL CONFERENCE ON IMAGE PROCESSING, VOL III, PROCEEDINGSIMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite fromit. Please check the document version below.Document VersionPublisher's PDF, also known as Version of recordPublication date:2001Link to publication in University of Groningen/UMCG research databaseCitation for published version (APA):Meijster, A., & Wilkinson, M. H. F. (2001). Fast computation of morphological area pattern spectra. In 2001INTERNATIONAL CONFERENCE ON IMAGE PROCESSING, VOL III, PROCEEDINGS (pp. 668-671).University of Groningen, Johann Bernoulli Institute for Mathematics and Computer Science.CopyrightOther than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of theauthor(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-amendment.Take-down policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons thenumber of authors shown on this cover page is limited to 10 maximum.Download date: 11-10-2022FAST COMPUTATION OF MORPHOLOGICAL AREA PATTERN SPECTRAArnold MeijsterCentre for High Performance ComputingUniversity of Groningen, P.O. Box 8009700 AV, Groningen, The NetherlandsMichael H. F. WilkinsonInstitute of Mathematics and ComputingScience, University of Groningen, P.O. Box800, 9700 AV, Groningen, The NetherlandsABSTRACTAn area based counterpart of the binary structural openingspectra is developed. It is shown that these area opening andclosing spectra can be computed using an adaptation of Tar-jan’s union-find algorithm. These spectra provide rotation,translation, and scale invariant pattern vectors for textureanalysis.1. INTRODUCTIONMultiscale feature extraction from textures or images canbe done by means of size distributions or granulometries[1, 2]. These are ordered sets of morphological openingswhich remove details below given scales. Usually the scaleis defined by the width of features, though other measuresare possible [3, 4]. One method to summarize the actionof a size distribution on a particular image in a single, 1-D array is called a pattern spectrum [5, 6]. In the binarycase, this summarizes how many foreground pixels remainas a function of the scale of the openings used. Figure 1shows an example of a binary image and two morphologicalpattern spectra. Image details are present at those scaleswhere the slope of the spectrum is steep.Fast algorithms for the binary case are available [7, 8],but computing patterns spectra in the grey scale case canbe time consuming, because many types of granulometriesrequire repeated openings and summing of grey levels tocompute the spectrum. Though computationally efficientalgorithms exist for certain openings with particular classesof structuring elements [9], granulometries using Euclideandiscs to achieve rotation invariance of the method are stillcostly. Breen and Jones [4] have shown that spectra basedon attribute openings, which are connected set filters canbe computed efficiently using a pixel-heap based algorithm.Attribute openings are extensions of area operators whichform a comparatively new class of image processing opera-tors [3]. They filter out or highlight detail based on the area,rather than the width of an object. In the case of attributeopenings, the filters can be based on a far larger class ofattributes of each connected component. Rotation invariant0 2.5 5 7.5 10 12.5 1502505007501000125015001750feature width (pixels)spectrum Sα XFig. 1. A binary image and two morphological pattern spec-tra derived from granulometries based on structural open-ing by discs of increasing sizes (solid), and openings-by-reconstruction with similar structuring elements (dashed).attribute openings are easily designed, so rotation invariantspectra can be obtained quickly. In this paper we propose anew algorithm, based on Tarjan’s union-find method. Thismethod is a variant of earlier work on area and attributeopenings [10, 11]. We will focus on the case of area spectra,though the method is easily extended to other attributes.In this paper we first describe opening and closing spec-tra, and the definitions of area openings and closings. Wethen go on to define area pattern spectra, and describe thenew algorithm for their computation. Finally we give someexperimental results of computational efficiency.2. THEORY2.1. Opening and Closing Pattern SpectraLet binary images X and Y be defined as the set of fore-ground pixels, i.e., a subset of the image domain M ⊂ Rn(usually n = 2). Grey scale images are a mapping from Mto R. A size distribution or granulometry is a set of opera-tors {αr} with r from some ordered set Λ (usually Λ ⊂ R6680-7803-6725-1/01/$10.00 ©2001 IEEEor Z), with the following three propertiesαr(X) ⊆X (1)X ⊆ Y ⇒ αr(X) ⊆ αr(Y ) (2)αr(αs(X)) =αmax(r,s)(X), (3)in the binary case, and in the grey scale caseαr(f) ≤f (4)f ≤ g ⇒ αr(f) ≤ αr(g) (5)αr(αs(f)) =αmax(r,s)(f), (6)for all r, s ∈ Λ. Since (3) and (6) imply idempotence, itcan be seen that size distributions are openings. The patternspectrum SαX obtained by applying a size distribution αr tobinary image X is defined asSαX(r) = A(αr(X)). (7)where A(X) is a function denoting the Lesbesgue measureof set X in Rn, or the number of pixels in Zn. In the grey-scale case the spectrum Sαf (r) of image f is just the integralof the grey level of αr(f) over the image domain, or in thediscrete case the sum:Sαf (r) =∑x∈M(αr(f))(x). (8)2.2. Area Openings and ClosingsThe theory of area operators is given only briefly here. For amore thorough discussion the reader is referred to [3]. Herewe will first discuss binary area openings and closings, andthen the extension to the grey scale case. Binary area open-ings are based on binary connected openings. Let the setX ⊆ M denote a binary image with domain M. The bi-nary connected opening Γx(X) of X at point x ∈ M yieldsthe connected component of X containing x if x ∈ X , and∅ otherwise. Thus Γx extracts the connected component towhich x belongs, discarding all others.The binary area opening can now be defined as:Definition 1 Let X ⊆ M and r ≥ 0. The binary areaopening of X with scale parameter r is given byΓar(X) = {x ∈ X|A(Γx(X)) ≥ r} (9)The binary area closing can be defined by duality:Φar(X) = [Γar(Xc)]c (10)The definition of an area opening of a grey scale imagef is most easily derived from binary images Th(f) obtainedby thresholding f at h, which can be defined as:Th(f) = {x ∈ M|f(x) ≥ h} (11)0 0.5 100.20.40.60.810 0.5 100.20.40.60.81Fig. 2. An example of area pattern spectra: (top): imagesof diatoms showing fine internal details; (bottom): corre-sponding pattern spectra obtained from series of area open-ings with increasing size criterion. Axes have been normal-izedDefinition 2 The area opening for a mapping f : M −→ Ris given by:(γar (f))(x) = sup{h|x ∈ Γar(Th(f))}. (12)The grey scale area closing φar is defined as:(φar(f))(x) = inf{h|x ∈ Φar((Th(f))c)}, (13)The interpretation of these definitions is that the areaopening of an image assigns each point the highest thresh-old at which it still belongs to a connected foreground com-ponent of area r or larger. The area closing assigns eachpoint the lowest threshold at which it belongs to a connectedbackground component of area r or larger.Vincent [3] has shown that area openings and closingsare size distributions and anti-size distributions respectively.Therefore, area opening and closing spectra can be definedas in (8). Fig. 2 shows an example of such spectra.3. COMPUTING AREA PATTERN SPECTRALet Ns be the number of scales we wish to use for our pat-tern spectra. A naive implementation would compute Nsopenings and compute the sum of grey levels for each open-ing. This is very costly for even modest values of Ns. Wetherefore propose a new algorithm, which processes the im-age just once, rather than O(Ns) times.Before going into the details of the algorithm, we firstdefine a flat zone Lh at level h of a grey scale image f as aconnected component of the set of pixels {p ∈ M|f(p) =h}. A regional maximum Mh at level h is a flat zone nomembers of which have neighbors larger than h. A peakcomponent Ph at level h is a connected component of Th(f).At each level h there may be several such sets, which willbe indexed as Lih, Pjh and Mkh , respectively, with i, j, and k669from some index set. It can be seen that any regional max-imum Mkh is also a peak component, but the reverse is nottrue. Finally, let Λ be the set of values the attribute (in thiscase area) of the peak components can assume for a givenimage size, and B be a mapping from Λ to [0..Ns), whichis the index range of the spectrum, stored in array spec.The idea is to process the image in grey scale order, fill-ing peak components at all grey levels starting at the re-gional maxima. Each time we move to a new grey level hwe may find some P jh′ at grey level h′ > h which must bemerged with some Lkh. In this case, spec[B(A(Pkh′))] mustbe incremented by (h′−h)A(P kh′), because an area openingof r = A(P kh′) would remove this quantity from the sum ofgrey levels at that point in the image.To make this work we need an efficient way to keeptrack of the peak components, which are of course disjointsets. Tarjan [12] presents the union-find algorithm whichprovides a general method for keeping track of disjoint sets.It allows performing set-union operations on sets which arein some way equivalent, while ensuring that the end prod-uct of such a union is disjoint from any other set. Tarjanuses tree structures to represent sets. Each non-root node ina tree points to its parent, while the root is flagged in someway. Two objects x and y are members of the same set ifand only if x and y are nodes of the same tree, which isequivalent to saying that they share the same root. Thereare three important operations.• Makeset(x): Create a new singleton set {x}.• FindRoot(x): Return the root element of the setcontaining x.• Union(x,y): Compute the union of the two setscontaining x and y.Each of the sets is represented as a tree, with each pixelcontaining a pointer to its parent pixel. To store the treesfor the entire image, we use an integer array parent of thesame size as the image I (i.e. N ), in which parent[p] isthe parent of pixel p. Pixels are stored as width*y+x,withx and y the pixel’s x and y coordinates, and width the im-age width. If a pixel is a root of a tree, i.e. it has no par-ent, we flag this by setting parent[p] = -area, witharea the number of pixels in the set. In the following dis-cussion pixel p refers to the integer representation of thecoordinates, whereas I[p] is the pixel’s grey level.Finally, we need to ensure that for a peak componentPh at level h, its root element r has a grey level I[r]= h.Therefore, we always make the last pixel processed the rootof the new tree. This is done by radix-sorting the pixels, andstoring the coordinates in an array S of length N . Pixels ofthe same grey level are processed in scan line order.As can be seen in Figure 3, The Union procedure per-forms the following task: Find the root r of neighbor n;void MakeSet ( int x ){ parent[x] = -1;}void Link ( int x, int y ){ parent[y] = parent[y] + parent[x];parent[x] = y;}int FindRoot ( int x ){ if ( parent[x] >=0 ){ parent[x] = FindRoot( parent[x] );return parent[x];}else return x;}void Union ( int n, int p ){ int r=FindRoot(n);if ( r != p ){ if ( I[p]<>I[r]{ spec[B(-parent[r])] -=(I[r]-I[p]) * parent[r];}Link( r, p );}}Fig. 3. The basic operations for area spectra. The root el-ement now stores the area of the component as a negativenumber.FindRoot(p) is not necessary because the current pixelp is always a root, because if linking has occurred it hasbecome the root of the resulting tree, and otherwise it is asingleton. If the root r=p we must Link them. Beforelinking we check whether I[r]= I[p], and if so, updatethe appropriate bin in spec.The complete algorithm is shown in Fig. 4. The specarray is first set zero, after which the pixels are processedin reverse grey level order. For each pixel p a union withthose neighbours which have already been processed is per-formed. After this stage spec[i] is the amount subtractedfrom the sum of grey levels at the ith scale, but not at lowerscales. The desired spectral value Si is given bySi = Si−1 − spec[i] if i > 0 (14)For i = 0 the correct value is the sum of grey levels minusspec[0]. This is carried out by the last loop in Fig. 4.670/* array S contains sorted pixel list */for (i=0;i<Ns;i++)spec[i]=0;greysum=0;for (p=0; p<Length(S); p++){pix = S[p];greysum += I[p];MakeSet(pix);for all neighbors nb of pix doif ((I[pix] < I[nb]) ||((I[pix] == I[nb]) && (nb<pix)))Union(nb,pix);}for (i=0;i<Ns;i++){ spec[i] = greysum - spec[i];greysum = spec[i];}Fig. 4. Pseudo-code showing how to compute an area open-ing pattern spectrum using the operations of Fig. 3.Table 1. Timing results of different algorithms (in seconds)Method min median maxnaive,Ns = 256 23.0 30.7 33.3pixel-heap 0.14 0.68 1.42union-find 0.092 0.124 0.1274. RESULTSThree algorithms were implemented in C: (i) the naive im-plementation using fast area openings by union-find withNs = 256 (ii) the direct approach using the pixel-heapmethod (for details see [4]), and (iii) the direct implementa-tion using union-find. All implementation were optimizedmanually for maximum performance. Timings were per-formed on 10 natural images of 256 × 256 pixels. Theresults are shown in Table 1. As expected, the direct ap-proaches outperform the naive by a factor approaching Ns.The union-find approach outperforms the pixel-heap methodby a factor of about 5, consistent with earlier findings onattribute openings [11]. Similar results were obtained onimages of 768 × 1024 (data not shown).5. DISCUSSIONWe have shown that the new algorithm is up to ten timesfaster than the pixel-heap method on natural images. Thecomplexity of the algorithm is identical to that of a singlearea or attribute opening using union-find, which was shownto be O(N log N), whereas that of the pixel-heap method isO(N2 log N), with N the number of pixels [10, 11]. Thismethod can readily be extended to other attribute operators,in the same way as in [11] for openings and closings.6. REFERENCES[1] G. Matheron, Random Sets and Integral Geometry,John Wiley, 1975.[2] L. Vincent, “Granulometries and opening trees,” Fun-damenta Informaticae, vol. 41, pp. 57–90, 2000.[3] L. Vincent, “Morphological area openings and clos-ings for grey-scale images,” in Shape in Picture:Mathematical Description of Shape in Grey-level Im-ages, Y.-L. O, A. Toet, D. Foster, H. J. A. M. Heij-mans, and P. Meer, Eds., pp. 197–208. NATO, 1993.[4] E. J. Breen and R. Jones, “Attribute openings, thin-nings and granulometries,” Computer Vision and Im-age Understanding, vol. 64, no. 3, pp. 377–389, 1996.[5] J. Serra, Image Analysis and Mathematical Morphol-ogy, vol. 1, Academic Press, New York, 2 edition,1982.[6] P. Maragos, “Pattern spectrum and multiscale shaperepresentation,” IEEE Transactions on Pattern Anal-ysis and Machine Intelligence, vol. 11, pp. 701–715,1989.[7] P. F. M. Nacken, “Chamfer metrics in mathematicalmorphology,” Journal Mathematical Imaging and Vi-sion, vol. 4, pp. 233–253, 1994.[8] P. F. M. Nacken, “Chamfer metrics, the medial axisand mathematical morphology,” Journal Mathemati-cal Imaging and Vision, vol. 6, pp. 235–248, 1996.[9] L. Vincent, “Fast grayscale granulometry algo-rithms,” in ISMM’94, Mathematical Morphologyand its Application to Image Processing, J. Serra andP. Soille, Eds., Fontainebeau, France, 1994, pp. 265–272, Kluwer Academic Publishers.[10] A. Meijster and M. H. F. Wilkinson, “An efficientalgorithm for morphological area operators,” Tech.Rep. 99-9-07, Institute for Mathematics and Comput-ing Science, University of Groningen, submitted.[11] M. H. F. Wilkinson and J. B. T. M. Roerdink,“Fast morphological attribute operations using Tar-jan’s union-find algorithm,” in Proceedings of theISMM2000, Palo Alto, CA, June 2000, pp. 311–320.[12] R. E. Tarjan, “Efficiency of a good but not linear setunion algorithm,” Journal of the ACM, vol. 22, pp.215–225, 1975.671

Fast computation of morphological area pattern spectra

https://pure.rug.nl/ws/files/3103295/2001ProcICIPMeijster.pdf

Fast computation of morphological area pattern spectra

Abstract

Similar works

Full text

Available Versions

Proceedings - University of Groningen