Retinex image enhancement technique is an attempt to model human vision system. Scope of this nonlinear enhancement algorithm is immense which makes it a powerful algorithm. The most important application of Retinex algorithm is the development of visibility improvement system for helping drivers with poor vision at night and bad weather condition^1-9. The algorithm is extensively used in Langley Research Centre, NASA, and DLR, Germany for their technology demonstrator projects of development of enhanced vision system for NASA 757 aircraft. The algorithm is also extensively used in medical image processing area in the TruView Imaging Company formed by the scientists of NASA who worked extensively for the development of this algorithm^10-11. This paper presents the work done on Retinex so far and proposes a method to reduce the halo artifacts for better performance of human vision system.

Keywords:    Colour image enhancement,   colour constancy,   dynamic range compression,   Gaussian surround function,   lightness number,   halo artifacts,   multi-scale retinex

In recent days, with the increasing use of digital camera capturing and rendering a good image is not a trivial task. Recorded images differ from direct view to the lack of dynamic range compression and colour constancy. Dynamic range is the ratio of the largest value to the smallest value of a physical quantity. Dynamic range of a scene is the ratio of the brightest and darkest part of the scene. The range of human vision system is quite large. The luminance of starlight is around 0.001 cd/m² to that of a sunlit is around 100,000 cd/m², which is hundred million times higher. The human eye can accommodate a dynamic range of approximately 10,000:1 in a single view, but the imaging system has a restricted dynamic range. To render a captured scene with very high dynamic range therefore remains a challenging task.

Colour constancy is defined as the maintenance of colour appearance despite variation in the colour of nearby objects and despite variation in the spectral power distribution. A captured image of rose under different lighting condition appears differently. While the colour of rose appears to be same to us in sunlight as well as room illuminated by light source, a camera film does not see it in the same way. If colour appearance is useful feature in identifying an object, then colour appearance must be constant when the object is viewed in different context. Colour constancy is the ability of human being under natural variation of light sources, but it is not perfect. If the illuminant is strongly saturated, then we cannot predict colour correctly.

The human vision system performs the task of dynamic range compression and colour constancy almost effortlessly. Thus to improve the quality of images we have to combine dynamic range compression, colour constancy and colour and lightness rendition. Multi scale Retinex with colour restoration (MSRCR), based on Land’s Retinex theory^12-16, is an attempt to achieve these goals. This is a powerful algorithm especially for the visibility improvement of the dark regions.

Retinex theory¹⁷ was proposed in 1983 as an attempt to develop a computation model of human perception of colours. Experiments based on colour perception^17-24, depends on neural structure of human vision system. On the basis of these experiments, Land concluded that the composition of the light from an area in an image does not specify the colour of that area. He found a quantity, termed it as ‘lightness’ exist which neither changes with change in illumination nor location of object on a scene. Thus on that basis of retinal cortical system there are three independent lightness making mechanism, one for long waves (R), one for middle waves (G) and the other for short waves (B) each served by its own spectrum independently. Each system forms a separate image of the scene.

In 1986, E. Land²³ proposed a simple alternative technique for computation of the designator in retinex theory. The designator is the computed numerical measure on one waveband of the lightness seen as part of the whole field of view. Instead of overall average for the denominator, average flux from contiguous areas was used. Unlike the previous Retinex technique¹⁷ which involved some kind of comparison between the flux coming from a point in remote as well as the use of contiguous area, new algorithm produced randomness in the average because of the arbitrary reflectivity and the smallness of the population of contiguous areas. Land proposed an inverse square surround.

Thus the Retinex theory can be implemented as follows:

For each waveband the relative reflectance for each point can be calculated as the relative reflectance of i to j:

$\begin{array}{l}{R}^{\mathrm{\Lambda }}(i,j)={\displaystyle \sum \delta \log {\mathrm{I}}_{k+1}/{\mathrm{I}}_k}\\ \delta \log \frac{{\mathrm{I}}_{k+1}}{{\mathrm{I}}_k}=\{\begin{array}{l}\log {\mathrm{I}}_{k+1}/{\mathrm{I}}_{k}if\left|\log \frac{{\mathrm{I}}_{k+1}}{{\mathrm{I}}_k}\right|>Th\\ 0if\left|\log \frac{{\mathrm{I}}_{k+1}}{{\mathrm{I}}_k}\right|<Th\end{array}\end{array}$ (1)

Th is the threshold.

Role of Th in Retinex is to decrease the effect channel intensity ratios.

Average relative reflectance of an area i is given as

$R^{\mathrm{\Lambda }}=\frac{{\displaystyle {\sum }_{j=1}^N{R}^{\mathrm{\Lambda }}}(i,j)}{N}$ (2)

3.1 Path Algorithm

The Retinex algorithm proposed was the first attempt to develop a computational model for human colour constancy. In random path algorithm, next pixel position is chosen randomly from the current pixel position. David¹⁹, et al. studied this algorithm and formulated the path using stochastic method: An accumulator is designed to calculate lightness values of each pixel which is initialised to zero. As the path precedes, the path value of accumulator A gets updated. This process is repeated for n number of paths. The final value for each pixel is obtained by normalising the accumulator value by the number of paths visited. If the path is too long, then the Retinex algorithm is equivalent to its maximum value. They found the algorithm too sensitive to changes in the colour of nearby objects, and hence an inadequate model of human colour constancy. Horn formulated calculation of lightness using Laplacian operator⁸. Andrew Blake¹⁸ introduces an improved version of Horn’s algorithm. He used gradient operator instead of Laplacian operator. This improved boundary condition along the image. A. Rizzi²⁵, et al. used Brownian motion to randomly distributed path. Main drawback with path algorithm is the complexity associated with number of paths, path lengths and their trajectories.

3.2 Iterative Algorithm

Frankle-McCann Retinex in 1983 proposed a multi-resolution iterative version of the algorithm. Funt²⁶, et al. published its Matlab Code in 2004. R. Sobol²⁷ proposed an improvement in this iterative model by introducing ratio-modification operator which preserves small contrast ratio and significantly compresses large contrast.

3.3 Centre/Surround Algorithm

Land evolved the concept from random walk computation to its last form as a centre/surround form. This last form of Retinex captured attention of Jobson¹⁵, et al. who studied the properties of center/surround function. In this method, the pixel under consideration is replaced by a value which depends on the weighted average of the surrounding pixels. Study of centre/surround method and propose method to reduce artifacts present around high contrast edges are presented in this paper.

The single-scale Retinex (SSR)²⁵ for a single spectral component can be represented as

$R_i(x,y)=\log {\mathrm{I}}_i(x,y)-\log [F(x,y)*{\mathrm{I}}_i(\text{x,y})]$ (3)

where I_i(x, y) is the image distribution in the ith colour spectral band, x and y are image coordinates. * denotes the convolution operation, F(x, y) is the surround function and R_i(x, y) is the associated Retinex output.

$F(x,y)=K\exp (-r^2/C^2)$ (4)

C is a scalar value called surround space constant and K is selected such that

$r=\sqrt{x2}+y2and{\displaystyle \int {\displaystyle \int F(x,y)dx,dy=1}}$ (5)

Results obtained with different Gaussian space constant values are shown in Fig. 3.

Small kernels results in dynamic range compression. With large kernel the output looks more like natural impression of the image. Middle value of surround space constant is good for compensating shadow and to achieve acceptable rendition for the processed image.

When the dynamic range of scene exceeds the dynamic range of the recording medium, there is a loss of information which cannot be recovered. Single-scale Retinex (SSR) can either provide dynamic range compression or tonal rendition but not both simultaneously. To combine the strength of various surround space Multi-scale Retinex (MSR) was developed¹².

Multi-scale Retinex output is the weighted sum of the several different SSR outputs. Mathematically,

$R_{MSR}={\displaystyle {\sum }_{i=1}^n{W}_n{R}_{ni}}$ (6)

where n is the number of scales, R_ni is the i^th component of the n^th scale, R_MSR is the i^th spectral component of the MSR output, and W_n is the weight associated with the n^th scale. In MSR the surround function is given by

$F_n(x,y)=K\exp (-r^2/C_n^2)$ (7)

where C_n is the Gaussian surround space constant. Multi-scale Retinex combines the dynamic range compression of the Single-scale Retinex with the tonal rendition to produce an output which encompasses both as shown in Fig. 4.

Although MSR gives better results by combining dynamic range compression and colour rendition, it suffers from ‘graying-out’ of uniform zones.

Images which violate gray-world assumptions on retinex processing suffer from ‘graying-out’ of the image, either globally or in specific regions. Thus gives a washed out appearance for the processed images. To restore colour^12-13,15, MSR was modified by adding a colour restoration function given as:

$R_{MSRCR}=C_i(x,y)R_{MSRi}(x,y)$ (8)

where,

$C_i(x,y)=f[{{\mathrm{I}}^{\prime }}_i(x,y)]$ (9)

${{\mathrm{I}}^{\prime }}_i(x,y)={\mathrm{I}}_i(x,y)/{\displaystyle {\sum }_{i=1}^S{\mathrm{I}}_i(x,y)}$ (10)

where ${{\text{I}}^{\prime }}_i(x,y)$ is the i^th band of the colour restoration function (CRF), R_MSRCR is the ith spectral band of the Multi-scale Retinex with colour restoration(MSRCR), S is the number of spectral channels. The value of S is generally 3 for RGB.

Results show that the function that provides the best overall colour restoration is given by

$C_i(x,y)=\beta \log [\alpha {{\mathrm{I}}^{\prime }}_i(x,y)]$ (11)

where α controls the strength of the non-linearity and β is a gain constant. The final gain/offset adjustments are required for transition from the logarithmic to the display domain. The final version of the MSRCR can be written as

$\begin{array}{l}{R}_{MSRCR}=G[C_i(x,y)\{log{\mathrm{I}}_i(x,y)\\ -\log [F(x,y)*{\mathrm{I}}_i(x,y)]\}+b]\end{array}$ (12)

where G and b are the final gain and offset values, respectively. The value of G and b are implementation dependent.

Choice of appropriate space constant values is one of the design issue. The algorithm¹⁵ is designed for constant values:

N	c1	c2	c3	α	β	wn
3	15	80	250	125	46	1/3

The result of Multi-scale–Retinex with colour restoration processing on the test image is shown in Fig. 5.

Multi-scale–Retinex with colour restoration algorithm effectively restores colours in the processed images but shows artifacts in the images having high contrast edges, as shown in Fig. 5. In this section we propose a method of reducing these halo artifacts. This work is inspired by the work of Moore²⁸, according to which the induction can be reduced by multiplying surround by gradient of the image. Unlike Moore’s work, the authors have used a Multi-scale approach and Laplacian of Gaussian as edge detector.

$\begin{array}{l}\text{Retinex}\text{algorithm}\text{has}\text{the}\text{form:}\\ \text{Retinex}=\text{Centre}-\text{Surround}\end{array}$ (13)

The surround function is given as:

$Surround=\log [F(x,y)*{\text{I}}_i(x,y)]$ (14)

To reduce the artifacts we modify the surround weight as shown below:

$\begin{array}{l}Surround=\log [{F^{\prime }}_n(x,y)*{\text{I}}_i(x,y)\\ \cdot log{\text{I}}_i(x,y)*{\text{F}}_n(x,y)]\end{array}$ (15)

where ${F^{\prime }}_n(x,y)$ the Laplacian of Gaussian is function and ${\text{F}}_n(x,y)$ is Gaussian function of the nth scale. Proposed algorithm is given in Appendix-A.

The surround function is multiplied by edge detector, thus in smooth region the output will be almost equal to zero. Hence surround value will be almost zero for smooth regions. This would result more colour rendition and reduction of artifacts near the high contrast edges. Image in Fig. 6 is processed by the proposed method with three scale values such as 3, 15, 80. This not only reduces artifacts, but also decreases the computational time.

To evaluate the performance of the proposed method, test has been conducted on images collected from various sources given in the web site. Proposed method is designed for three scales and surround space constant values chosen are 3, 15, 80. Some of the results are presented in the paper.

The face of the lady in Fig. 7(a) is not clear in MSRCR output Fig. 7(b) compared to that obtained by our method. Thus this algorithm cannot be suitable for face detection or recognition purposes. The algorithm works well for images having dark zones.

Quality assessment plays an important role in image processing application. Many cameras had inbuilt image processing algorithms to correct the quality of image before rendering to mimic closely to the perceived image by the observers. How pleasing the overall rendered image is and how well it conforms to the observer’s expectation is measured by number of factors called attributes of the image. The quality Q of an image is determined by a number of visual attributes that includes sharpness, tone rendition and colour. In this experiment we have calculated mean square error, structure similarity (SSIM) and mean structural similarity (MSSIM).

Mean square error is given by the equation:

$MSE={\displaystyle \sum {(f_i(x,y)-{f^{\prime }}_i(x,y))}^2}/MN$

Structural similarity is given by the equation:

$\begin{array}{l}SSIM(x,y)=(2{\mu }_x{\mu }_y+C_1)(2{\sigma }_{xy}+C_2)\\ /\{({\mu }_x^2+{\mu }_y^2+C_1)({\sigma }_x^2+{\sigma }_y^2+C_2\end{array}$

Mean structural similarity index measurement is given

by:

$MSSIM(x,y)=({\displaystyle \sum SSIM(x_j,y_j))/M}$

Tables 1 and 2 shows that proposed methods give better results compared to result obtained by Rahman, et al., especially with single scale where MSE is found to be highest. But the SSIM and MSSIM values are almost the same for all the results except for C where MSSIM result is quite high. Moreover Single scale result is takes less computational time as compared to Rahman, et al. algorithm.

The time complexity of the algorithm can be subdivided into three parts:

1. computation time required for Gaussian mask,
2. time complexity for convolution and logarithm, and
3. time required for colour restoration function.

Let the size of the image be of dimension M × N.

The computation time required for Gaussian mask is equal to 16NM+M+N+1. Hence the order of Gaussian mask is O(MN), if the dimension of the image is same say N, then the O(N²). To calculate time complexity for convolution operation we considered the worst condition where for all pixels are multiplied by 3 x 3 mask i.e., 9 multiplication and 9 addition thus the computational time is equal to 18 operation on each pixel so for N × M image the operations are 18MN(approximately). Thus the order of the image is equal to O(MN) if the image is of N x N dimension then it is of order O(N²).The logarithmic calculation is of same order. Since the MSRCR uses three scales that and the scales are 15 × 15, 80 × 80 and 250, the computation time is much higher than that of our proposed method. The colour restoration function requires (14MN+N+M+1) i.e. also order of O(NM). Hence the total complexity of the algorithm is of O(MN) if the dimensions are different otherwise O(N²). This calculation is done only one colour space out of RGB, the same time it required to other two planes. Hence total complexity is given by O(3N²).

The SSR provides a good mechanism for enhancing images. However depending on the surround space constant, it can either provide good tonal rendition or dynamic range compression. The MSR, comprised of the three scales such as, small, intermediate and large overcomes these limitations. The scene which violates the gray –world assumption, desaturation of colour was found to occur The MSRCR adds a colour restoration scheme which produced good colour rendition even for severe gray-world violations. MSRCR shows artifacts around high contrast edge. Proposed method reduces the artifacts, but at the same time suffers from desaturation of colour. It effectively pulls out the details from dark areas. Since designed with scales 3, 15, 80 the computational time is also reduced. Future work is to find the optimal range of surround space constant as well as number of scales required for proper colour rendition.

The images used for evaluating the performance of the algorithm are taken from www.dragon.larc.nasa.gov.

1. Ngo, H.; Tao, Li.; Ming, Zhang; Livingston A. & Asari, V. A visibility improvement system for low vision drivers by non linear enhancement of fused visible and infrared video. In Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPRW’05), 2005, 03.

2. Hines, G.D.; Rahman, Z.; D. Jobson, D.J. & Woodell, G.A. A real-time enhancement, registration, and fusion for a multi-sensor enhanced vision system. In Enhanced and Synthetic Vision 2005, SPIE 2006, 6226.

3. Verley J.G. Autonomous Landing guidance System Validation. In Enhanced and Synthetic Vision, SPIE, 1997, 3088, 19-25.

4. Clark, S.M. & Durant-Whyte, H. Autonomous land vehicle navigation using MMW Radar. In Proceedings of the International Conference on Robotics and Automation ‘98, May 1998, 4, 3697-702.

5. Doehler, H.U. & Bollmeyer, D. Simulation of imaging radar for obstacle avoidance and enhanced vision. In Enhanced and Synthetic Vision, edited by J. G. Verly, SPIE, 1997, 3088, 64-73.

6. Ferguson. D, Radke. J. System for Adverse Weather Landing. In Proceedings of the AIAA Aircraft Design, Systems and Operations Meeting, 1993, 11-13, 1-5.

7. Hecker, P. & Doehler, H.U. Enhanced vision systems: Results of simulation and operational tests. In Enhanced and synthetic vision 1998, edited by J. G. Verly, SPIE, 1998, 3364, 115-22.

8. Horn, W.F.; Stanley, Ekiert; Jeff, Radke; Richard, R.; Tucker, C.; Harold Hannan, J. & Allen, Zak. Synthetic vision system technology demonstration methodology and results to date. SAE International Technical Paper Series No. 921971, Aerotech “92, October 1992, 5-8, 1-9.

9. Nordwall, B.D. UV Sensor proposed as pilot landing aid. Aviation Week Space Technol., 1997, 147, 81-84.

10. Rahman, Z. Retinex image enhancement: Application to medical Images. In NASA workshop on New Partnerships in Medical Diagnostic Imaging, Greenbelt, Maryland, USA. July 2001.

11. Rahman, Z.; Jobson, D.J. & Woodell, G.A. MultiScale Retinex for color image enhancement. In Proceeding of the IEEE International Conference on Image Processing(ICIP), 1996.

12. Jobson, D.J.; Rahman, Z. & Woodell, G.A. Retinex image processing: Improved fidelity for direct visual observation. In Proceedings of the IS&T Fourth Color Imaging Conference: Color Science, Systems, and Applications, 1996, pp. 124-126.

13. Rahman, Z.; Jobson, D.J. & Woodell, G.A. Multiscale Retinex for color rendition and dynamic range compression. In Applications of Digital Image Processing XIX, SPIE, 1996, 2847.

14. Jobson, D.J.; Rahman, Z. & Woodell, G.A. Properties and performance of a center/surround retinex. IEEE Trans Image Proces., 1997, 6, 451-62.

15. Jobson, D.J.; Rahman, Z. & Woodell, G.A. A multi-scale Retinex for bridging the gap between color images and the human observation of scenes. IEEE Trans. Image Proces., 1997, 6(7), 965-76.

16. Wang, W. Fast Multi-Scale Retinex algorithm for color image enhancement. In Proceedings of the 2008 International Conference on Wavelet Analysis and Pattern Recognition, Hong Kong, Aug 2008, 80-85.

17. Hecker, P. Recent advances in Retinex theory and some implications for cortical computations: Color vision and the natural image. In Proceedings of National Academy of Sciences, USA. 1983, 80, 5163-169.

18. Andrew, Blake. Boundary conditions for lightness computation in Mondrian World. Comp. Vision, Graphics Image Process., 1985, 32, 314-27.

19. Brainard, David H. & Wandell, Brian A. Analysis of retinex theory of color vision. J. Opt. Soc. Am. A., 1986, 3(10), 1651-661.

20. Land, E.H. Recent advances in Retinex theory and some implications for cortical computations: Color vision and the natural image. In Proceedings of National Academy of Sciences, USA. 1983, 80, 5163-569.

21. Rahman, Z. Jobson, D.J.; Woodell, G.A. & Hines, G.D. Multi-sensor fusion and enhancement using the Retinex image enhancement algorithm. Visual Information Processing XI, Proc. SPIE 2002, 4736.

22. Rahman, Z.; Jobson, D.J. & Woodell, G.A. Retinex processing for automatic image enhancement. J. Electro. Imaging, 2004, 13(10), 100-110.

23. Land, E.H. An alternative technique for the computation of the designator in the Retinex theory of color vision. In Proceedings of National Academy of Sciences, USA. 1986, 83, 3078-080.

24. Land, E.H. Recent advances in Retinex theory. Vision Research, 1987, 26(1), 7-21.

25. Rizzi, A.; Gatta, C. & Marini, D. A new algorithm for unsupervised global and local color correction. Pattern Recognition Letters, 2003, 24(11), 1663-1677.

26. Funt, B. & Ciuera, C. Tuning Retinex parameter. J. Electro. Imaging, 2004, 13(1), 58-64.

27. Sobol, R. Improving the Retinex algorithm for rendering wide dynamic range photograph. J. Electro. Imaging, 2004, 13.

28. Moore, A.; Allman, J. & Goodman, M. A real-time neural system for color constancy. IEEE Trans. Neural Networks, 1991, 2, 237-47.

Appendix-‘A’

Algorithm: MSRCR with Halo Reduction

Input: Colour image I( M × N)

Output: Enhanced colour image I(M × N)

It includes the following steps:

Step 1. The input image is read and stored in an array.

Step 2. The individual colour components are separated and stored in individual arrays as d (R(x,y)), green (G(x,y)) and blue (B(x,y)).

Step 3. The Gaussian functions that are required for the convolution are calculated using the following equations.

$F(x,y)=K*\exp (-(x^2+y^2)/C^2)$

where C is the Gaussian surround space constant. K determined by:

${\displaystyle \iint K*\exp (-(x^2+y^2)/C^2)dxdy=1}$

Thus K is determined as K = 1/(2 × π × c²).

Step 4. The Gaussian functions are determined for two scale values of c₁ = 3, c₂ = 15. The convolution is done for each scale to get ${\text{I}}_i(x,y)*F(x,y)$ as

$\begin{array}{l}{\text{I}}_i(x,y)*F(x,y)=\\ {\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{n=0}^{N-1}\text{I}(\text{m,n})\text{F(x}-\text{m,}\text{y}-\text{n)}}}\end{array}$

Step 5. The same procedure of convolution is repeated for $F^{\prime }(x,y)$ , where $F^{\prime }(x,y)$ is given Laplacian of Gaussian.

Step 6. The outputs of the result obtained after convolving with two functions are added and the resultant is

Step 7. The convolution result for each scale is used in

$R_i(x,y)=\log (I_i(x,y)/{I^{\prime }}_i(x,y))$

i = 1: Red component,

i = 2: Green component,

i = 3: Blue component

Step 8. The multi-scale result is weighted sum of the output s of several different SSR outputs.

${R^{\prime }}_{MSRi}={\displaystyle {\sum }_{i=1}^{N}wn{R^{\prime }}_{ni}}$

where N is the number of scales, ${R^{\prime }}_{ni}$ is the ith component of the n^th scale, ${R^{\prime }}_{MSRi}$ is the ith spectral component of the MSR with halo reduction output and W_n is the weight associated with it. W_n =1/2 n=1, 2

Step 9. To calculate colourimetric transform:

${I^{″}}_i(x,y)=I_i(x,y)/{\displaystyle {\sum }_{i=1}^S{I}_i(x,y)}$

where S is the number of spectral channels.

Step 10. The colour restoration function (CRF) in the chromaticity space is calculated as

${\mathrm{C}}_i(x,y)=\beta \log [\alpha {{\text{I}}^{″}}_i(x,y)]$

where β is a gain constant = 46, α is linearity the strength of the non linearity = 125

Step 11. The MSRCR with halo reduction can be written as:

${R^{\prime }}_{MSRCRi}=C_i(x,y){R^{\prime }}_{MSRi}(x,y)$

Step 12. The MSRCR with gain offset can be written as:

Step 13. ${\text{R}}_{MSRCRi}(x,y)=G\{{\text{C}}_i(x,y)({\text{R}}_{MSRi}(x,y)+b)\}$

Step 14. The colour components are recombined to get the processed output image, which is then stored in ‘raw’ format.