Image Anal Stereol 2002;21:201-206 Original Research Paper MATHEMATICAL MORPHOLOGY IN THE CIELAB SPACE Allan Hanbury1 and Jean Serra2 1 Pattern Recognition and Image Processing Group, Vienna University of Technology, Favoritenstraße 9/1832, A-1040 Vienna, Austria, 2Centre de Morphologie Mathematique, Ecole des Mines de Paris, 35 rue Saint-Honore, 77305 Fontainebleau Cedex, France e-mail: hanbury©prip.tuwien.ac.at, serra@cmm.ensmp.fr (Accepted September 27, 2002) ABSTRACT The use of mathematical morphology in the CIE L*a*b* colour space is discussed. It is possible to impose a total order on the colour vectors in this space by using a weighting function and lexicographical order. An order analogous to one by colour saturation is suggested by making use of a weighting function based on an electrostatic potential. This weighting function assigns a lower weight to colour vectors near the colours with maximum chroma, and higher weights to colour vectors near the lightness axis. The use of morphological operators with the colour vector order imposed by this function is demonstrated. Finally, a top-hat operator making use of the Euclidean colour distance in the L*a*b* space is introduced. Keywords: CIE L*a*b* colour space, lexicographical order, mathematical morphology. INTRODUCTION Much research has been carried out on the application of mathematical morphology to colour images, mostly based on defining orders for the vectors in the RGB colour space (Comer and Delp, 1999). Use of these formulations usually presents the disadvantage of having to arbitrarily choose one of the red, green or blue channels to play a dominant role in the ordering, although attempts have been made to overcome this limitation through the use of, for example, bit-interlacing (Chanussot and Lambert, 1998). The application of mathematical morphology in a colour space which has an angular hue component (Hanbury and Serra, 2001b) can overcome this disadvantage, allowing a non-constrained choice of the dominant hue, or permitting the implementation of rotationally invariant operators independent of the hue. In 1976, the International Commission on Illumination (CIE), introduced the L*a*b* and L*u*v* colour spaces, which were designed to be perceptually uniform and device-independent (i.e. independent of the devices which produce and display the image). These spaces also allow one to take into account the illumination characteristics of an image. Due to its device-independence, the L*a*b* space is ideal for exchanging colorimetric measurements between different observers. Nevertheless, a transformation from the RGB space to the L*a*b* space results in an irregularly shaped gamut of colours, its shape being dependent on the illumination conditions. It is possible to show that in order to avoid introducing false colours when applying morphological operators in a vector space, it is sufficient to find an injective function mapping the vectors to a complete lattice. In other words, one needs a function which associates a unique numerical weight with each vector in the colour space, allowing the vectors to be ordered by using their weights. As it is very difficult to find such a function, the requirements can be made less demanding through the use of a lexicographical order (Hanbury and Serra, 2001a). In the L*a*b* space, a measure of colour saturation does not exist, but is it possible to introduce one? In this article, we consider the use of a weighting function to simulate a colour order by saturation. MATERIALS AND METHODS THE L*a*b* SPACE The principal advantage of the L*a*b* space is its perceptual uniformity. Two colours which appear similar to a human observer lie close together in the L*a*b* space (their separation is measured by the Euclidean distance). The transformation from the RGB to the L*a*b* colour space is done by first transforming to the CIE XYZ space, and then to the L*a*b* space (Wyszecki and Stiles, 1982). When one captures an image using an RGB device, each colour in the image is specified as a triplet of values giving the amount of each of the device-specific primaries (red, green and blue) in the colour. The XYZ coordinates of these device-specific primaries can vary between devices, hence the RGB image is device-dependent. When one transforms from RGB to XYZ 201 H ANBURY A ET AL: Mathematical Morphology in the CIELAB Space coordinates, the coordinates of the device-specific primaries are taken into account, thereby making the transformed coordinates device-independent. The white point (nominally white object-colour stimulus) of the image, which is the colour obtained when the R, G and B coordinates assume their maximum values, depends on the illumination of the scene captured. The white point is taken into account by both the RGB to XYZ, and XYZ to L*a*b* transformations. If one knows the illumination conditions used when acquiring the image, then the white point can be directly specified. If the illumination conditions are unknown, a hypothesis must be made, the most common being to choose the CIE D65 daylight illuminant, which has been done for the transformation in this article. Alternatively, algorithms exist for estimating the white point of an image (Risson, 2001). In the L*a*b* colour space: L* represents the lightness; a* encodes the red-green sensation, with positive a* indicating a red colour, and negative a* a green colour; and b* encodes the yellow-blue sensation. The grey-levels or colourless points are located on the lightness axis (a* = 0. b* = 0), with black at L* = 0, and white at L* = 100. In a cylindrical coordinate representation of the space, if we take the lightness V as the "vertical" axis, then in the polar orthogonal space, one has the hue h*, measured anti-clockwise from the positive a* axis (different to the HLS or HSV hue H); and the chroma C*, the perpendicular distance to the L* axis. These cylindrical coordinates are a more convenient representation for the solution of some problems. The colour difference A£^ between two colours, each expressed in terms of L*, a* and b* is given by the Euclidean metric AEaV (AL+)2 + (Aa*)2 + (Ai>+)2 (1) In the cylindrical representation, the Euclidean distance between two colours {L\.h\.C\) and (L|;^,q),is ^ [(AT)2 + (C?)2+(C|)2 IClC^cosikl-hl)]'2 . (2) Fig. 1. (a) The a*b* histogram which results from transforming an RGB cube to the L*a*b* space, (b) The values of the extrema of the chroma C* (+) and their corresponding lightness L* (x) as a function of hue h*. We now consider the shape characteristics of the colour gamut in the L*a*b* space. To visualise the distribution of the points in the RGB space when transformed to the L*a*b* space, a two-dimensional chrominance histogram is calculated by transforming from an RGB colour cube containing points equally spaced by V256 throughout the region [0,1] x [0,1] x [0.1]. For each point [R. G. B], the resulting coordinates a* and b* are rounded to the nearest integer, and the bin [a%£*] of the histogram is incremented. A plot of the histogram is shown in Fig. la. In this image, the grey-level at each point [a*, b*] indicates the number of pixels in the RGB cube mapped to this point. In the full three-dimensional L*a*b* space, these pixels would be mapped to different lightness values. However, this two-dimensional histogram allows one to make some useful observations: Firstly, the largest number 202 Image Anal Stereol 2002;21:201-206 of points are found at the origin, the position of the lightness axis; and secondly, the distribution of colours is not circular — the maximum value assumed by C* depends on h*. The extremal points of the colour gamut in the L*a*b* space are those points which are furthest away from the lightness axis, i.e. the points with the maximum values of C*. These extrema are visible on the histogram in Fig. la. It is, however, interesting to find the lightness value corresponding to each of these extremal points. For the transformation from the RGB cube, for each integer value of h*, the point with the largest chroma C* was found. In Fig. lb, the value of C* is plotted for the extremal points corresponding to each integer value of h*, along with the lightness L* of the extremal point. These functions are henceforth denoted as Cext (h*) and Lext (h*), where the values can be read off the graph for integer values of h*, and interpolated for non-integer values. MATHEMATICAL MORPHOLOGY We aim to order the colour vectors in the L*a*b* space in a way analogous to an order by colour saturation. Hence, equal prominence should be given to the vectors closest to the extremal points of the colour gamut, and the position of the lightness axis should be taken into account. Due to the irregular shape of the colour gamut, one cannot simply define the supremum to be the colour with the largest value of C*, as the maximum possible value of C* varies with hue. We also cannot use the Euclidean distance of colours from the edges of the colour gamut, as the lightness axis is not in the geometrical centre of the colour gamut, and hence the maximum value of this distance function does not coincide with the lightness axis. Use of two-dimensional homothetic distance functions centred on the lightness axis in the planes perpendicular to this axis does not allow one to the take extremal points into consideration, and gives an equal weight to all the points on the edge of the colour gamut. We wish to define a weighting function which assigns a weight to each colour as a function of its distance to an extremal point. The weighting function w associates a weight wf with each colour vector ft w: .:fi=m.M,c:) W; where a lower weight implies a colour closer to the extremal points, and a higher weight indicates a colour closer to the lightness axis. The best solution found is to take the weighting as the value of an electrostatic Vr potential function in the L*a*b* space obtained by placing "charges" at various astutely chosen positions. Note that the electrostatic potential model was chosen for convenience. We are in no way modeling a physical situation, only making use of a concept which is well understood to simplify the problem at hand. We therefore do not make use of any units or constants from electrostatic theory. The values and positions of the charges in the L*a*b* space are, of course, heuristically chosen to produce the most useful potential function. As the numerical values of the potentials are not important, the magnitudes of the charges can be adjusted to produce values in a useful range. The potential is set up in the L*a*b* space by placing a line of positive "charge" on the (vertical) lightness axis, thereby ensuring that the surrounding greys have the highest potential; and placing negative "charges" at the extremal points of the colour gamut, imposing minima on the potential function. The potential due to a line charge can be determined analytically. Given a line charge of length /, and a point on the line at a distance x from the line centre, the potential V+ at a point with a perpendicular distance of d from point x on the tine charge is Vr+ = Aln a + Va2 + d2 (3) -jt and where A is the charge per unit distance, a = b = | +je. Note that for a numerical implementation, it is not advisable to combine the two log terms into one, as this leads to numerical instability. For the lightness axis, we aim to give a slightly lower potential to lighter greys. This is accomplished by using a positively charged line of length 200 lightness units, with the L* = 0 point of the lightness axis placed at the centre of this line. To summarise, at the point i with coordinates (Lj* .k-.Cf), the potential due to the positive line charge is calculated using Eq. 3 with X = 1, / mdd = Cf 200, x = L*t The line of negative "charge" at the extremal points is approximated by placing equidistantly spaced charges of equal magnitude at positions given by the functions Lext(A*) and Cext(h*). Because of the variations in the distance between these charges and the lightness axis, it is necessary to take into account the spatial distances between the charges, one cannot simply separate them by equal angular distances. In summary, we construct the line of extremal charge as a set of n point charges of magnitude — q, having coordinates [L^^.^.C^^)] with i € {1,...,«}. The values of the 0,- are chosen so that the charges are equidistantly spaced, that is to say that the Euclidean distance between each pair of neighbouring charges is 203 HANBURY A et AL: Mathematical Morphology in the CIELAB Space equal to a constant dext. The potential V_ at point i due to these charges is V- = -t^~, (4) where r is the Euclidean distance from an extremal point j to the point i at which the potential is being calculated. The distance r is given by Eq. 2. The positions 0, of the charges were calculated with dext = 2, and the charges were assigned unit magnitude (q = 1). The values of dext and q give the level of approximation of the line of negative charge by the set of points. Values of rfext = 1 and q = V2, for example, would give a better approximation, but with a longer calculation time. In a practical application, the calculation time is not necessarily critical — one could initially calculate the weights for all the colours in the L*a*b* gamut once, and then use a look-up table to associate them with each pixel in an image. The weight at point i is then calculated as the sum of these two potentials (Eqs. 3 and 4), i.e. wt = V+ + V_. A diagram of the equipotential lines due to the suggested charge distribution on a slice along the lightness axis in the L*a*b* space is shown in Fig. 2. Colours along each line have identical weights. Fig. 2. The equipotential lines for a vertical slice through the L*a*b* space. The hue value on the left of the lightness axis is 180°, and on the right 0°. The potential decreases from the centre outward. In order to apply a morphological operator to an image in the L*a*b* space, we first calculate the weight for each colour vector in the space. A weight image corresponding to a colour image can then be produced by replacing each colour pixel with its corresponding weight, as shown in Fig. 3. In the weight image, the grey-level represents the weight of the colour vector, and hence darker pixels indicate colours which are closer to the extremal points. Fig. 3. (a) The example image and (b) its weight image. With the potential function approach, we have created a lattice of equipotential surfaces. The colour vectors making up an equipotential surface have not yet been ordered. When defining colour morphological operators, it is advantageous that no pairs of vectors exist for which an order is not explicitly specified. To obtain such a total ordering of the vectors, we make use of the lexicographical order. The order of the equipotential surfaces is placed at the first level. For two vectors in the same equipotential surface (i.e. having equal weights), we say that the one with a higher lightness value is larger than the other. If the lightness values are also identical, then the hue must be taken into account. This lexicographical order for two arbitrary colour vectors f and /¦ is {W: < W:. Or • J ' Wi = Wj and L*>L), or w, = w, and L?=L* 1 J t J and (A,-rÄ0)<(Äy.-rA0). (5) where hQ is an arbitrarily chosen hue origin, and a1 +a2 represents the acute angle between angles a1 and a2, calculated as {|flj — a2\ if jflj — a2\ < 180° . 360° - K -a2\ if \ar-a2\ > 180° '. (6) 204 Image Anal Stereol 2002;21:201-206 In a practical application of a morphological operator to a "normal" image, the second and third levels of the lexicographical order are rarely used (Hanbury and Serra, 2001a). It is, nevertheless, possible to influence the level of usage of these lower levels by quantising the weights in the first level. If the weights are quantised so as to be represented by integer values from 0 to N, then it is clear that when N has a smaller value, the lower levels of the lexicographical order have a larger influence on the result. Once these orders have been defined, the morphological operators are defined in the standard way. The vector erosion of an image / at point x by structuring element B is %/(*) = {/(y): /(y) = inf[/(z)],z € 5X} , (7) and the corresponding dilation is obtained by replacing the inf by a sup. An opening yB is an erosion followed by a dilation, and a closing