Method of enhancing color images

[0002] This application claims benefit to provisional application No. 60/163,417 filed Nov. 3, 1999.

[0003] The present invention relates to a method of enhancing colour images.

[0004] Digital radiographic images are acquired with a large dynamic range (usually 12 bits per pixel). Because the human visual system is not able to distinguish all the available grey level differences and because the image has to be brought back to 8 bits per pixel for printing and displaying, post-processing of the image is often necessary.

[0005] Multi-resolution image enhancement is being used routinely nowadays for processing of digital radiographs, i.e. for grey tone images. A multi-scale image enhancement method for radiographic images has for example been described in European patent 527 525 and European patent application 1 001 370.

[0006] An iterative processing method of the above-described kind has been described in European patent application 610 603.

[0007] As acquisition techniques improve for all kinds of images, the need for effective dynamic range compression still increases.

[0008] For colour photographic images, especially for those taken in outdoor scenes, one is very often confronted with a very wide dynamic range. When reproduced this often results in an image having areas that are too dark or too light.

[0009] Effective dynamic range compression without loss of important information is then required. As many different colours as possible have to be visualised.

[0010] Especially for photographic images it is furthermore unacceptable that the applied image processing wound result in an unnatural appearance of the processed colour image

[0011] In the article ‘Multiscale color image enhancement’, by Toet A. in Pattern Recognition Letters, NL, North Holland Publ. Amsterdam, Vol. 13, No. 3, (Mar. 1, 1992), a multiscale image decomposition method has been disclosed for application to color images.

[0012] In this method the luminance and saturation components of an image are first decomposed into contrast detail images of different spatial scales. The contrast is defined here as the ratio of the image or a low resolution version of the image and a successive lower resolution version of the image. Next, a new set of multi-scale luminance contrast primitives is then constructed by modulating the original luminance primitives at every location in the image and at every spatial scale by their corresponding saturation contrast primitives. Reconstruction of the color image from the resulting set of multiscale primitives provides a representation of the original image in which local luminance contrast is enhanced at all resolution levels.

[0013] It is an object of the present invention to provide a method for enhancing the image quality of colour images.

[0014] The above mentioned object is realised by an image processing method as described in claim 1. Specific features for preferred embodiments of the invention are disclosed in the dependent claims.

[0015] The method of the present invention is advantageous in that effective dynamic range compression can be obtained in combination with subtle feature enhancement.

[0016] The method of the present invention comprises the steps of

[0017] decomposing at least one of the colour component images of a colour image into a multi-resolution image representation comprising bandpass detail images at multiple resolution levels and a residual image,

[0018] modifying the multi-resolution image representation of at least one colour component image at at least one resolution level so as to obtain (a) modified multi-resolution image representation(s),

[0019] reconstructing said colour component image(s) from the multi-resolution representation(s) by applying a reconstruction algorithm to the residual image and the modified and unmodified bandpass detail images of the respective colour component image(s), said reconstruction algorithm being such that if it were applied to the residual image and unmodified bandpass detail images of a colour component image, the colour component image or an approximation thereof would be obtained.

[0020] Preferably the modification is non-linear. At least one non-linear modifying function is applied to the detail image(s)

[0021] In a specific embodiment the multi-resolution image representation(s) is/are modified by multiplying the value of a detail image at pixel i and at resolution level j with a factor which is obtained by evaluating a non-linear modifying function in an argument value which depends on pixel values of detail images of said colour components at pixel i and resolution level j.

[0022] The multi-resolution image representation(s) is/are for example modified by multiplying the value of a detail image at pixel i and at resolution level j with a factor which is obtained by evaluating a non-linear modifying function in an argument value which depends on pixel values of detail images of said colour components at resolution level j and pixel i.

[0023] In one embodiment the multi-resolution representation comprises directional bandpass detail images. The multi-resolution representation is e.g. a multi-resolution gradient representation.

[0024] In a specific embodiment the argument of the modifying function is the norm of the colour gradient in a specific pixel and specified resolution level.

[0025] The modifying function may for example have the form set out in claim 8. It may depend on a pixel value in the vaginal image.

[0026] Further details on the modifying function are described further on.

[0027] The colour component images may be defined in a perceptual colour space such as the L, u, v colour space.

[0028] A gamut may be defined and the colour component images may be re-scaled to that defined gamut after reconstruction. The gamut may be the gamut of the original image.

[0029] The present invention as well as specific and/or preferred embodiments thereof will be explained in detail with reference to the accompanying drawings.

[0030] FIG. 1 illustrates in general the sequence of image decomposition, modification and reconstruction according to the method of the present invention,

[0031] FIG. 2 illustrates the decomposition—with subsampling—of an image into a sequence of detail images at different resolution levels,

[0032] FIG. 3 illustrates the reconstruction of an image from modified detail images obtained by a decomposition according to the method which was illustrated in FIG. 2 and subsequently modified,

[0033] FIG. 4 illustrates the decomposition—without subsampling—of an image into a sequence of detail images at different resolution levels,

[0034] FIG. 5 illustrates the reconstruction of an image which was decomposed according to the method illustrated in FIG. 4 and subsequently modified,

[0035] FIG. 6 shows the function y(x),

[0036] FIG. 7 is an example of a modifying function.

[0037] The grey level of grey-level images is commonly represented by a single number. Pixels of colour images however are commonly specified by three numbers. These numbers represent coordinates of each pixel relative to coordinate axes in a specific colour space.

[0038] In the following the ensemble of the coordinates referring to a specific axis for all pixels of a colour image will be referred to as ‘colour component’. A single colour component of a colour image can be seen as a grey-valued image.

[0039] In a so-called R(red) G(green) B(blue) colour space, a colour is specified by the amounts of red (R), green (G) and blue (B) primary colour present in the composite colour of each pixel. This representation is much used because of its compatibility with additive colour reproduction systems.

[0040] In more perceptual colour spaces, such as HSV, Lab, Luv, colour is represented by a brightness component (V=value in HSV and L=lightness in Lab, Luv) and two other additional components which encode the chromaticity (e.g. hue, represented by ‘H’, and saturation represented by ‘S’ in the HSV colour space).

[0041] More detailed information on colour spaces and colour representations can be found in the following references: Pratt W. K., Digital Image Processing, Wiley, New York, (1991) and Wyseszcki G., Stiles W. S., Color Science: concepts and methods, quantitative data and formulae, Wiley, New York (1982).

[0042] For enhancement purposes an RGB image representation is preferably converted to its coordinates in a perceptual colour space. This is preferred since it is the objective to modify perceptual colour differences which is easier when a perceptual representation is used.

[0043] Lab and Luv are considered to be more perceptually uniform than HSV (i.e. numerical differences correspond better to perceptual ones), and therefore they are superior for the purpose of this invention.

[0044] The HSV colour space is advantageous because of its simplicity and the fact that in this colour space chromaticity is split into hue and saturation.

[0045] This splitting can however also be realised by using polar coordinates in the ab-plane in the Lab colour space or in the uv-plane when the Luv colour space is used.

[0046] FIG. 1 illustrates the general concept of the present invention.

[0047] A colour image to be processed is represented by a set of primary colour component images for example R, G, B component images.

[0048] In a first step, which is optional, the representations of the colour components of an image are transformed into a corresponding representation in a conceptual colour space such as the Luv colour space.

[0049] In this example it was assumed that the colour image was originally furnished as an R, G, B image. The R, G, B image representation is then converted into a Luv image.

[0050] Next, at least one of the L, u, v component images, or further on called L, u and v colour planes, is decomposed into a multi-resolution representation (also called multi-scale representation) comprising a set of bandpass detail images at different resolution levels and a residual image.

[0051] The detail images each represent the details present in the image at a certain resolution level. The details of an image at resolution level j can be interpreted as the difference of information between its approximation at a resolution level j and its approximation at a lower resolution level j+1.

[0052] The detail images contain the image details at a resolution level j. Each resolution level corresponds with a spatial frequency band and the detail images at a resolution level contain the image detail in the frequency band corresponding with that resolution level, i.e. all detail images are band-pass images.

[0053] In the described embodiment the multi-resolution image representation is a multi-scale gradient representation as will be explained furtheron. The image decomposition comprises a set of detail images each consisting of a horizontal and a vertical gradient component. The image gradient captures the local variations of the image. It is computed by convolving the image with two directional high pass filters. The simplest implementation makes use of a two-element filter g=(1,−1) for the horizontal direction and its transpose for the vertical direction.

[0054] It will be clear that other kinds of multi-resolution representations may be used such as Laplacian (pyramid) representations.

[0055] In the next step the multi-resolution representation of at least one component of the image in the selected colour space is modified. For this purpose the multi-resolution representation of this colour component image is pixel-wise subjected to a modification according to at least one modifying function.

[0056] The modifying function may be identical for all resolution levels of the multi-resolution representation. Alternatively it may differ between resolution levels of the same multi-resolution representation.

[0057] Likewise it is not required to use the same modification function or the same set of modification functions for application to the multi-resolution representations of each of the colour component images in the selected colour space.

[0058] FIG. 1 illustrates a specific example wherein a modification is applied to each of the image components L, u and v. This is however not necesary. It is also possible that only a single component, for example the L component of the image, is modified while the other colour component images remain unaltered or visa versa.

[0059] It is advantageous that at the finest levels of the multi-resolution representations only the elements of the multi-resolution representation of the L colour component image are modified, while the finest resolution levels of the u and v colour component images are left unaltered. In this way the additional generation or boosting of colour noise is avoided.

[0060] The modification function is a non-linear function which provides that details of small amplitude are enlarged at the expense of the ones with larger amplitude.

[0061] As a result of this modification a set of modified detail images is obtained for each of the colour component images to which modification has been applied.

[0062] In the fourth step the modified multi-resolution representation of each of the colour component images is subjected to a reconstruction processing so as to obtain modified colour component images L′, u′, v′.

[0063] If desired, the reconstructed perceptual colour space components may be transformed again into R, G, B components.

[0064] The following is a detailed description of the different steps of the present invention. The first step is illustrated in FIG. 2. This step consists of the decomposition of a colour component image into a corresponding multi-scale gradient representation whereby detail is represented by the discrete gradient in horizontal and vertical direction.

[0065] The following explanation is given for one of the colour components of the colour image, for example for the ‘L’-component. It is however to be understood that in this embodiment a similar procedure is applied to the other colour component images of the colour image in the selected colour space, so also for the u and v component in this described specific embodiment.

[0066] A colour component image is in the following indicated by the denomination I. I thus for example represents the L-component of the colour image representation in the Luv colour space. The colour image component is a map from a grid of pixels P⊂Z¹to R.

[0067] The decomposition of an image into detail images at successively coarser scale is illustrated in FIG. 2.

[0068] The decomposition comprises first the generation of approximation images of the original image at successively coarser scale. Looking at the image at a coarser scale is formulated mathematically by smoothing and subsampling of the image.

[0069] More specifically, a one-dimensional low-pass filter h is first selected, e.g. h={fraction (1/16)}(1,4,6,4,1) (underlined element is the zeroth element).

[0070] This one-dimensional low-pass filter is used to construct a two-dimensional low-pass filter h{circle around (X)}h

[0071] This filter acts on the image by convolving it with h in the horizontal direction and then in the vertical direction.

[0072] Next a set of approximation images S_jI of the image I at scales j,0≦j≦K are generated. An approximation image S_jI at scale j whereby 0≦j≦K, of the image I is defined recursively by the following formulae

[0073] S₀I=I

[0074] S_jI=↓(h{circle around (X)}h)*S_j−1I  (1)

[0075] Here ↓ denotes subsampling of the image by a factor of 2 in the horizontal and vertical directions, keeping only pixels (2x,2y) for (x,y) in the original image lattice.

[0076] The discrete gradient {right arrow over (G)}_jI of each approximation image is then calculated by applying the formula

[0077] {right arrow over (G)}_jI=((g{circle around (X)}δ)*S_jI,(δ{circle around (X)}g)*S_jI)  (2)

[0078] where g is a gradient filter (e.g. g=({fraction (1/2)}, −½) and δ is the Kronecker delta or unity filter.

[0079] (g{circle around (X)}δ) convolves with g in the horizontal direction and (δ{circle around (X)}g) convolves with g in the vertical direction.

[0080] The set {{right arrow over (G)}₀I,{right arrow over (G)}₁I, . . . {right arrow over (G)}_nI,S_KI}

[0081] is called the multiscale gradient representation of I.

[0082] In the described specific embodiment a two-dimensional gradient representation has been selected for describing the image detail present in an image at multiple resolution levels.

[0083] Other representations could also have been used instead of the gradient representation, such as a Laplacian representation.

[0084] A representation can only be useful for the problem of image enhancement if the decomposition process is reversible, i.e. if one is able to reconstruct the image from its decomposed representation. In other words, it is necessary that the invers transform of the decomposition transform exists.

[0085] The reconstruction is illustrated in FIG. 3.

[0086] To fulfill this requirement, the following constraints must be satisfied: two additional filters p, q must exist such that

[0087] p*h+q*g=δ  (3)

[0088] the odd elements of p are zero

[0089] This last property guarantees that there exists another filter r such that for any image I

[0090] ↓(p{circle around (X)}p)*I=(r{circle around (X)}r)*(↓I)  (4)

[0091] Relation (3) can be extended to two dimensions as $\left(5\right)$ $\begin{array}{c}\delta \otimes \delta =\left(p*h+q*g\right)\otimes \left(p*h+q*g\right)\\ =\left(p*h\right)\otimes \left(p*h\right)+\left(q*g\right)\otimes \left(p*h\right)+\left(p*h\right)\otimes \left(q*g\right)\text{+}\left(q*g\right)\otimes \left(q*g\right)\\ =\left(p*h\right)\otimes \left(p*h\right)+\left(q*g\right)\otimes \left(p*h\text{+}\frac{1}{2}q*g\right)+\left(p*h\text{+}\frac{1}{2}q*g\right)\otimes \left(q*g\right)\\ =\left(p*h\right)\otimes \left(p*h\right)+\left(q*g\right)\otimes \left(\frac{\delta +p*h}{2}\right)+\left(\frac{\delta +p*h}{2}\right)\otimes \left(q*g\right)\\ =\left(p*h\right)\otimes \left(p*h\right)+\left(q*g\right)\otimes l+l\otimes \left(q*g\right)\\ =\left(p\otimes p\right)*\left(h\otimes h\right)+\left(q\otimes l\right)*\left(q\otimes \delta \right)+\left(l\otimes q\right)*\left(\delta \otimes g\right)\end{array}$

[0092] where the filter l is defined by the relation $l=\left(\frac{\delta +p*h}{2}\right)$

[0093] and therefore (using (5) and (4)): $\begin{array}{cc}\begin{array}{c}\downarrow S_{j}I=\downarrow \left(p\otimes p\right)*\left(h\otimes h\right)*S_{j}I+\downarrow \left(q\otimes l\right)*\left(q\otimes \delta \right)*S_{j}I*\downarrow \left(l\otimes q\right)*\left(\delta \otimes g\right)*S_{j}I\\ =\left(r\otimes r\right)*\downarrow \left(h\otimes h\right)*S_{j}I+\downarrow \left(q\otimes l\right)*\left(g\otimes \delta \right)*S_{j}I+\downarrow \left(l\otimes q\right)*\left(\delta \otimes g\right)*S_{j}I\\ =\left(r\otimes r\right)*S_{j+1}I+\downarrow \left(q\otimes l\right)*G_j^{\left(h\right)}I+\downarrow \left(l\otimes q\right)*G_j^{\left(v\right)}I\end{array}& \left(6\right)\end{array}$

[0094] The superscripts (h) and (v) refer to the horizontal and vertical components of the gradient.

[0095] This last equation expresses that a subsampled version ↓S_jI of the approximation image S_jI can be recovered from the coarser approximation image S_j+1I and the gradient {right arrow over (G)}_jI of S_jI.

[0096] To find the full approximation image S_jI from its subsampled version the gradient can be used again to do an exact interpolation: for a pixel (2x+1,2y) lying in the middle of the two sample pixels (2x,2y) and (2x+2,2y) we have that

[0097] S_jI(2x+1,2y)=½[S_jI(2x,2y)−2G_j^(h)I(2x,2y)+S_jI(2x+2,2y)+2G_j^(h)I(2x+1,2y)]  (7)

[0098] For a pixel (2x,2y+1) lying in the middle of the two sample pixels (2x,2y) and (2x,2y+2) we have that

[0099] S_jI(2x,2y+1)=½[S_jI(2x,2y)−2G_j^(v)I(2x,2y)+S_jI(2x,2y+2)+2G_j^(v)I(2x,2y+1)]  (8)

[0100] For a pixel (2x+1,2y+1) lying in the middle of the four sample pixels (2x,2y), (2x+2,2y), (2x+2,2y), (2x,2y+2) and (2x+2,2y+2) we then have

[0101] S_jI(2x+1,2y+1)=¼

[0102] [S_jI(2x+1,2y)−2G_j^(v)I

[0103] (2x+1,2y)+S_jI

[0104] (2x,2y+1)−2G_j^(h)I

[0105] (2x,2y+1)+S_jI

[0106] (2x+1,2y+2)+2G_j^(v)I

[0107] (2x+1,2y+1)+S_jI

[0108] (2x+2,2y+1)+2G_j^(h)I

[0109] (2x+1,2y+1)]  (9)

[0110] The above described procedure (see equations (7), (8) and (9)) is called gradient based interpolation (GBI).

[0111] Therefore an image can be reconstructed from its multiscale gradient representation by recursively applying equations (6) through (9) for j=K−1 down to 0. (see FIG. 3)

[0112] This representation can also be implemented without subsampling if, instead of subsampling the approximation image, the filters are stretched by inserting zeros between any two elements and centered:

[0113] define $h_0=\frac{1}{16}\left(1,4,\underset{\_}{6},4,1\right)g_0=\frac{1}{2}\left(\underset{\_}{1},-1\right)$ $h_1=\frac{1}{16}\left(1,0,4,0,\underset{\_}{6},0,4,0,1\right)g_1=\frac{1}{2}\left(1,\underset{\_}{0},-1\right)$ $h_2=\frac{1}{16}\left(1,0,0,0,4,0,0,0,\underset{\_}{6},0,0,0,4,0,0,0,1\right)g_2=\frac{1}{2}\left(1,0,\underset{\_}{0},0,-1\right)$

[0114] The same definitions can be made for q_j,p_j,l_j,0≦j≧K. The implementation without subsampling is organized as follows (see FIG. 4):

[0115] S₀I=I

[0116] S_jI=(h_j−1{circle around (X)}h_j−1)*S_j−1I

[0117] {right arrow over (G)}_jI=((g_j{circle around (X)}δ)*S_jI,(δ{circle around (X)}g_j)*S_jI)

[0118] S_jI=(p_j{circle around (X)}p_j)*S_j+1I+(q_j{circle around (X)}l_j)*G_j^(h)I+(l_j{circle around (X)}q_j)*G_j^(v)I  (10)

[0119] To enhance the contrast in the image the multi-resolution image representation is modified in such a way that details of small amplitude are enlarged at the expense of the ones with larger amplitude. In this embodiment this modification is performed uniformly over all scales.

[0120] For this purpose the original transform coefficients resulting from the image decomposition are mapped in a continuous, monotonous non-linear way into modified coefficients.

[0121] While the colour components are each decomposed separately into their multi-resolution representation, it is preferable to perform the modification in a coupled way, i.e. to use information about all three colour components for the modification of an individual colour component.

[0122] For each colour component I the gradient (G_j^(h)I,G^(v)I) at modification level j has been computed.

[0123] Since the three components L, u, v make up the resulting colours, it makes sense to defined the colour gradient at resolution level j and pixel i as the sextet

[0124] (G_j^(h)L(i), G_j^(v)L(i), G_j^(h)u(i), G_j^(v)u(i), G_j^(h)v(i), G_j^(v)v(i))

[0125] In the Luv colour space distances between colours can be described by the Euclidean distance.

[0126] The distance between two colors C1 represented by the triplet (L1, u1, v1) and C2 represented by the triplet (L2, u2, v2) is given by

[0127] Δ=((L₁-L₂)²+(u₁-u₂)²+(v₁-v₂)²)^1/2

[0128] The norm of the colour gradient at a resolution level j and pixel i is therefore

[0129] Γ_j(i)=(∥{right arrow over (G)}_jL(i)∥²+∥{right arrow over (G)}_jv(i)∥²)^1/2∥{right arrow over (G)}_jI(i)∥=((G_j^(h)I(i))²+(G_j^(v)I(i))²)^1/2

[0130] If only the colour gradient norm is changed and the direction is left unaltered, unwanted colour shifts as well as spatial deformations in the reconstructed colour image can be avoided.

[0131] This results in the following method for modifying the multi-scale representation of a colour image.

[0132] {right arrow over (G)}′_Lj(i)={right arrow over (G)}_jL(i)·y(Γ_j(i))

[0133] {right arrow over (G)}′_uj(i)={right arrow over (G)}_ju(i)·y(Γ_j(i))

[0134] {right arrow over (G)}′_vj(i)={right arrow over (G)}_jv(i)·y(Γ_j(i))

[0135] The function y, being shown in FIG. 6, was used to multiply against the transform coefficients. This function is defined as follows $\begin{array}{cc}y\left(x\right)={\left(\frac{m}{c}\right)}^p\mathrm{if}\left|x\right|<cy\left(x\right)={\left(\frac{m}{\left|x\right|}\right)}^p\mathrm{if}c\le \left|x\right|\le my\left(x\right)=1\mathrm{if}\left|x\right|\ge m& \left(12\right)\end{array}$

[0136] m is an upper bound for the transform coefficients, transform coefficients larger than m are not modified. The value p, 0<p<1, determines the degree of nonlinearity.

[0137] The value c is a cross-over parameter. For values larger than c the modification is nonlinear, for values smaller than c the modification is linear. This is performed in order to prevent unnecessary overenhancement of the noise in the image.

[0138] The value of c is therefore preferably related to the noise content of the image, it is preferably equal to the estimated standard deviation of the noise in the image.

[0139] The parameter m can be chosen to be level dependent.

[0140] The function may additionally depend on the pixel values in he image so that the modification may depend on the image content.

[0141] By applying these formulae only colour differences for which differences in all three components are small, are enlarged.

[0142] This kind of enlargement always occurs in the direction of He colour difference; i.e. the direction of the colour gradient is never changed only its magnitude is changed.

[0143] An enhanced colour image (L′, u′, v′) is obtained by applying the reconstruction algorithm set out higher to the three thus adapted multi-scale gradient pyramidal representations.

[0144] By adapting only colour differences, the output range was not controlled. This might have as a consequence that a proper re-scaling of the dynamic range to a predefined range or gamut is required after reconstruction.

[0145] Good results are obtained by applying a linear resealing of the L′, u′, v′ values to the original L, u, v range respectively.

[0146] It is further possible to bring colours closer to the original colours by applying the following additional mapping procedure after resealing:

[0147] For every pixel i:

[0148] u′(i) is mapped onto u′(i)+{{overscore (u)}-{overscore (u′)}}

[0149] v′(i) is mapped onto v′(i)+{{overscore (v)}-{overscore (v′)}}

[0150] In these formulae the over-bar stands for the mean value of the image over all of its pixels. This mean value can be replaced by a weighted mean over a neighbourhood of the pixel i.

[0151] It is further possible to select which parts of the image are enhanced to a larger extent than other parts.

[0152] For this purpose for each pixel i the value L(i) is non-linearly mapped prior to the enhancement processing, for example L(i) is mapped onto {tilde over (L)}=L(i)^1−q100^qprior to the enhancement processing. After this processing the inverse transformation is applied to {tilde over (L)}′. The inverse transform maps {tilde over (L)}′(i) onto {tilde over (L)}′(i)1/1−q 100 ^−q/1−q.

[0153] If q>0 a darker part in the image is stretched while a lighter part is compressed. The effect is that the lighter part will be enhanced more than the darker part (as small differences get enlarged by the procedure of the present invention).

[0154] If q<0 the darker part will be enhanced more than the lighter part. The enhancement is in this case controlled by our parameters p, c, m and q. Very satisfactory results are obtained for 0.25≦p≦0.5.