A MULTI SCALE RETINEX WITH COLOR RESTORATION(MSRCR) FRAMEWORK FOR HAZE REMOVAL AND TONE MAPPING A THESIS
Submitted by SUDHARSAN P
In partial fulfillment for the award of the Degree of MASTER OF TECHNOLOGY IN ELECTRONICS AND COMMUNICATION ENGINEERING
(Signal Processing) Under the guidance of Dr Praveen Sankaran
DEPARTMENT OF ELECTRONICS AND COMMUNICATION ENGINEERING
NATIONAL INSTITUTE OF TECHNOLOGY CALICUT NIT CAMPUS PO, CALICUT KERALA, INDIA 673601. MAY 2012
CERTIFICATE
This is to certify that the thesis entitled: A Multi Scale Retinex with Color Restoration(MSRCR) framework for haze removal and tone mapping submitted by Mr. Sudharsan P to the National Institute of Technology Calicut towards partial fulfillment of the requirements for the award of the Degree of Master of Technology in Electronics and Communication Engineering (Signal Processing) is a bona fide record of the work carried out by him under my supervision and guidance.
Signed by Thesis Supervisors with names and date
Place: Date : Signature of Head of the Department
(Office seal)
ii
DECLARATION
I hereby declare that this submission is my own work and that to the best of my knowledge and belief, it contains no material previously published or written by another person nor material which has been accepted for the award of any other degree or diploma of the university or other institute of higher learning, except where due acknowledgement has been made in the text.
Place: Date:
Signature Name:
Reg.No:
iii
ABSTRACT A MULTI SCALE RETINEX WITH COLOR RESTORATION(MSRCR) FRAMEWORK FOR HAZE REMOVAL AND TONE MAPPING Sudharsan P National Institute of Technology Calicut 2012 Thesis Supervisor: Dr Praveen Sankaran Dynamic range of a scene is defined as the ratio of maximum to minimum luminance in a scene. The dynamic range of natural scenes are high and the illuminant conditions of the scene depends on the time of the day. Human Visual System(HVS) can perceive the correct colour of the objects in any illuminant conditions. Similarly HVS can also adapt itself to different dynamic ranges and perceive all the details of the scene clearly. But digital cameras are limited by the dynamic range and also the colour of the images captured depends on the illuminant conditions. Our aim is to enhance the quality of the recorded image as to how a human being would have perceived the scene. This property that we aim to achieve is called ‘color constancy’. Single Scale Retinex(SSR) was developed to achieve color constancy [1]. Multi Scale Retinex(MSR) algorithm [2] is a modification of the Single Scale Retinex algorithm. Multi Scale Retinex algorithm is an addition of Single Scale Retinex algorithm outputs obtained for different surround constants. A further modification of the Multi Scale Retinex algorithm is Multi Scale Retinex with Color Restoration(MSRCR) [2] which gives correct color rendition of the image. We propose a couple of modifications for the MSRCR algorithm such as fusion based MSRCR, automated MSRCR and applications of the MSRCR algorithm for haze removal and tone mapping. The modifications are proposed in the MSR stage. In fusion based MSR we fuse the three outputs of SSR algorithm instead of averaging them. Though this algorithm is computationally complex our aim is to find if it provides better results than the original algorithm. MSRCR is an image dependent algorithm. We propose an image independent MSRCR that finds parameters from the image under consideration and thus automate the algorithm. One of the interesting problems in image processing is restoration of degraded images. Outdoor images are often degraded by atmospheric effects such as haze, fog and rain. In this thesis we shall consider the problem of haze removal. One of the iv
popular haze removal algorithms proposed by Kaiming et al. [3] uses a dark channel prior based approach for haze removal. Though this approach gives very good results, this method is computationally complex. We extend the idea of a dark channel with a Retinex based approach to obtain a haze removal technique that gives good results and is also computationally simpler. Tone mapping is an algorithm that maps scenes captured by High Dynamic Range(HDR) cameras to Low Dynamic Range(LDR) devices without compromising on the details and color rendition. Tone mapping is achieved by both Global and Local tone mapping operators. In this thesis we have combined gamma mapping(Global operator) and Retinex(Local operator) to achieve good tone mapping results.
v
ACKNOWLEDGEMENT
First I would like to thank my advisor Dr. Praveen Sankaran for his invaluable project guidance. He has allowed me to work on my project independently with lots of freedom. It is his belief in me that has helped me to sail through the ups and downs of research in the past one year. Whenever I have faced any obstacle during the course of project he has been there to exactly pinpoint the flaw in my work and correct it. I would also like to thank the other panel members Dr. Deepthi. P. P and Dr. Sreelekha. G. for their invaluable suggestions which helped in the progress of my work. I am thankful to Dr. P. S. Sathidevi, The Head of the Department, ECE for encouraging our project work by providing us the necessary laboratory facilities. I would also like to thank my fellow members in the SPCOM lab with whom I have had numerous discussions during the course of project work. I would like to thank my family who have been a constant source of support in all my endeavours. Lastly but most importantly I would like to thank the Lord Almighty for providing me enough strength during my project work.
vi
c
Copyright, 2012, by Sudharsan P, All Rights Reserved
vii
TABLE OF CONTENTS Page List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xi
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv CHAPTERS 1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1
Image Enhancement . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.1.1
Spatial Domain Image Enhancement Techniques . . . . . . . .
1
1.1.2
Frequency domain enhancement techniques . . . . . . . . . . .
6
Colour image enhancement techniques . . . . . . . . . . . . . . . . .
7
1.2.1
Intensity Transformations . . . . . . . . . . . . . . . . . . . .
7
1.2.2
Histogram equalisation . . . . . . . . . . . . . . . . . . . . . .
7
1.3
Applications of Image Enhancement . . . . . . . . . . . . . . . . . . .
8
1.4
Thesis Flow: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
1.4.1
Fusion based Multi Scale Retinex with Color Restoration . . .
10
1.4.2
Automated Multi Scale Retinex with Color Restoration . . . .
11
1.4.3
Retinex based Haze Removal method . . . . . . . . . . . . . .
11
1.4.4
Tone Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
Literature Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
1.2
2
1
2.1
2.2
Retinex algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.1.1
Path version . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.1.2
Surround-Centre version . . . . . . . . . . . . . . . . . . . . .
14
2.1.3
Luma based MSR . . . . . . . . . . . . . . . . . . . . . . . . .
17
Image Fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
2.2.1
Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
2.2.2
Wavelet based fusion . . . . . . . . . . . . . . . . . . . . . . .
27
viii
3
2.3
Haze Removal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
2.4
Tone Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
2.4.1
HDR History . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
2.4.2
Dynamic Range of a scene . . . . . . . . . . . . . . . . . . . .
33
2.4.3
Tone mapping operators . . . . . . . . . . . . . . . . . . . . .
34
A fusion based approach for MSRCR . . . . . . . . . . . . . . . . . . . .
36
3.1
Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
3.1.1
Match Measure: . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.1.2
Activity Measure . . . . . . . . . . . . . . . . . . . . . . . . .
38
3.1.3
Combination: . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
3.1.4
Decision Map: . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
Dual Tree Complex Wavelet Transform(DTCWT) . . . . . . . . . . .
40
A novel automated approach for MSRCR . . . . . . . . . . . . . . . . . .
43
3.2 4
4.1 5
Proposed Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
A novel Retinex based approach for dehazing . . . . . . . . . . . . . . .
47
5.1
Koschmieder Model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
5.2
Dark Channel Prior Method . . . . . . . . . . . . . . . . . . . . . . .
48
5.2.1
Dark Channel Prior Equations . . . . . . . . . . . . . . . . . .
48
5.2.2
Transmission Map Refinement . . . . . . . . . . . . . . . . . .
49
Retinex Model for Haze Removal . . . . . . . . . . . . . . . . . . . .
51
5.3.1
Estimation of Atmospheric light . . . . . . . . . . . . . . . . .
52
5.3.2
A minor refinement on transmission map . . . . . . . . . . . .
52
Retinex based tone mapping . . . . . . . . . . . . . . . . . . . . . . . . .
54
5.3
6
6.1 7
Proposed Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
7.1
Automated Multi Scale Retinex . . . . . . . . . . . . . . . . . . . . .
55
7.2
Fusion based Multi Scale Retinex . . . . . . . . . . . . . . . . . . . .
55
7.3
Retinex based haze removal method . . . . . . . . . . . . . . . . . . .
58
7.3.1
Computational complexity . . . . . . . . . . . . . . . . . . . .
68
Retinex based tone mapping . . . . . . . . . . . . . . . . . . . . . . .
68
7.4
ix
8
Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . .
72
8.1
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
8.2
Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
LIST OF PUBLICATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . .
84
x
LIST OF ABBREVIATIONS
HVS
Human Visual System
SSR
Single Scale Retinex
MSR
Multi Scale Retinex
MSRCR
Multi Scale Retinex with Color Restoration
HDR
High Dynamic Range
HE
Histogram Equalisation
AHE
Adaptive Histogram Equalisation
CLAHE
Contrast Limited Adaptive Histogram Equalisation
RGB
Red Green Blue
CMY
Cyan Magenta Yellow
HSI
Hue Saturation Intensity
CCTV
Closed-circuit Television
MRA
Multi Resolution Analysis
RoLP
Ratio of Low Pass Pyramid
DWT
Discrete Wavelet Transform
SIDWT
Shift Invariant Discrete Wavelet Transform
MS
Maximum Selection
WA
Weighted Average
WBV
Window Based Verification
VLSI
Very Large Scale Integrated Circuits
CCD
Charge-Coupled Device
CRT
Cathode Ray Tube
LCD
Liquid Crystal Display
NASA
National Aeronautics and Space Administration
xi
LIST OF FIGURES Page 1.1
Gamma map
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2
a)Input image b)Gamma mapped image with γ =.5 . . . . . . . . . .
3
1.3
Contrast stretching [4] . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.4
a)Input image b) Contrast stretched image . . . . . . . . . . . . . . .
3
1.5
Intensity slicing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.6
a)Input Image b)Output image with some highlighted intensities . . .
5
1.7
a)Input Image b)Histogram equalised image . . . . . . . . . . . . . .
5
1.8
Unsharp masking . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.9
Homomorphic filtering . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.10 Medical image enhancement . . . . . . . . . . . . . . . . . . . . . . .
9
1.11 Image enhancement of old images [5] . . . . . . . . . . . . . . . . . .
9
1.12 Whiteboard [6] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.13 Forensic based image enhancement . . . . . . . . . . . . . . . . . . .
10
1.14 Haze formation model [3] . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.15 a) Hurlbert’s surround b) Land’s surround c) Moore’s surround . . .
15
2.16 Single scale retinex outputs a) c=15 b) c=80 c) c=250 . . . . . . . .
16
2.17 a) Input image b) Original MSRCR output c) Luma based MSRCR output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
2.18 Feature based fusion . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
2.19 Pixel based fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
2.20 Pyramid based fusion . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
2.21 Wavelet based fusion . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
2.22 Analysis 1-d filter bank . . . . . . . . . . . . . . . . . . . . . . . . . .
26
2.23 Synthesis 1-d filter bank . . . . . . . . . . . . . . . . . . . . . . . . .
26
2.24 Analysis 2-d filter bank . . . . . . . . . . . . . . . . . . . . . . . . . .
27
2.25 Synthesis 2-d filter bank . . . . . . . . . . . . . . . . . . . . . . . . .
28
2.26 Wavelet decomposition at first level . . . . . . . . . . . . . . . . . . .
28
2.27 Wavelet decomposition at second level . . . . . . . . . . . . . . . . .
29
3.1
DTCWT based fusion . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.2
DT-CWT analysis filter bank [7] . . . . . . . . . . . . . . . . . . . . .
41
3.3
DT-CWT synthesis filter bank [7] . . . . . . . . . . . . . . . . . . . .
42
xii
4.4
Histogram of a single scale retinex enhanced image [1] . . . . . . . . .
44
4.5
Histogram explaining choosing of clipping points using the variance .
45
4.6
Histogram explaining choosing of clipping points using frequency of 45
5.7
occurrence of pixels . . . . . . . . . . . . . . . . . . . . . . . . . . . . a)Input b)Transmission map t˜ c)Haze free image obtained using t˜ . .
5.8
a)Input b)Transmission map t c)Haze free image obtained using t . .
50
5.9
a)Input b)Transmission map t c)Haze free image obtained using t . .
53
5.10 a)Input b)Transmission map t c)Haze free image obtained using t . .
53
6.11 Proposed Tone Mapping method
. . . . . . . . . . . . . . . . . . . .
54
7.12 a) Input image b) Our output c) Software output . . . . . . . . . .
56
7.13 a) Input image b) Our output c) Software output . . . . . . . . . .
56
7.14 a) Input image b) Our output c) Software output . . . . . . . . . .
57
7.15 a) Input image b) Our output c) Software output . . . . . . . . . .
57
7.16 a) Input image b) Output due to addition c) Output due to fusion
59
7.17 a) Input image b)Output due to addition c) Output due to fusion .
60
49
7.18 Results of ‘Mountain’ Image obtained using dark channel prior method 61 7.19 Results of ‘Mountain’ Image obtained using Retinex based haze removal method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
7.20 Results of ‘Tokyo’ Image obtained using dark channel prior method .
63
7.21 Results of ‘Tokyo’ Image obtained using Retinex based haze removal method
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
7.22 Results of ‘Hongkong’ Image obtained using dark channel prior method 65 7.23 Results of ‘Hongkong’ Image obtained using Retinex based haze removal method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
7.24 Results of ‘House’ Image obtained using dark channel prior method .
67
7.25 Results of ‘House’ Image obtained using Retinex based haze removal method
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
7.26 a)Input HDR image b)Retinex output c)Gamma mapped HDR image d)Final Retinex output
. . . . . . . . . . . . . . . . . . . . . . . . .
70
7.27 a)Input HDR image b)Retinex output c)Gamma mapped HDR image d)Final Retinex output . . . . . . . . . . . . . . . . . . . . . . . . . .
xiii
71
LIST OF TABLES Page 7.1
Computational complexity of dark channel prior method . . . . . . .
68
7.2
Computational complexity of the proposed method . . . . . . . . . .
68
xiv
CHAPTER 1 INTRODUCTION An image is a two dimensional function g(x) where x represents the planar co ordinates (x1 , x2 ) and the amplitude g is called the intensity of the image at that point. A digital image has a finite number of points each with its own intensity. These elements are called picture elements or pels or pixels. Any image can be represented as a M × N matrix wherein the elements of the matrix represent the intensity values of the pixels. The intensity values are divided into L distinct levels, where L can be represented as a power of 2 as L = 2k . So each pixel effectively can be represented by k bits and the total number of bits required to represent the image is M × N × k bits. We have considered till now only gray scale image. A color image is a three dimensional function g(x1 , x2 , z) where the additional dimension is used to represent the color bands. In the third dimension only 3 color bands Red(R), Green(G) and Blue(B) are used. Any color can be represented as a weighted combination of these three colors. 1.1
IMAGE ENHANCEMENT
Image enhancement is the process of modifying an image so that it looks better. Image enhancement can be done either in spatial or frequency domain. Let us see some standard image enhancement techniques used. 1.1.1
Spatial Domain Image Enhancement Techniques
Spatial domain image enhancement techniques can act either on individual pixels separately or can act by taking the neighbourhood into account. 1.1.1.1
Point processing techniques(Gray Scale)
Point processing techniques act on individual pixels without taking the neighbourhood into account. This work is written as per IEEE Transactions format.
1
FIG. 1.1: Gamma map
Intensity transformations: A pixel of intensity r can be transferred to an intensity s by using the power law transformation s = c ∗ rγ . If γ = 1, it is a linear transformation. If γ < 1 the dark regions of the image are mapped to a wide range of output values and the bright regions of the image are compressed to a narrow range of output values. Exactly the opposite happens when γ is greater than 1. This is show in Fig. 1.1. Fig. 1.2 shows the input and gamma mapped image. Contrast stretching: Low contrast images occur due to poor illumination or lack of dynamic range in the acquisition system. In Contrast stretching the range of intensity levels in an image is expanded so that the image spans the full intensity range. Fig. 1.3 shows a transformation used for contrast stretching. Fig. 1.4 proves the advantage of using contrast stretching. Intensity slicing: In some applications there may be a requirement to highlight some specific intensity ranges as in medical imaging(to highlight the tumours), oceanography(to highlight the water flow). This can be done by intensity slicing. In this method , we scale the range of intensities that interests us to a high value and scale the remaining to a low value. In Fig. 1.5 the range of intensities that are to be highlighted are 100 to 200 and these are scaled to 220. The remaining intensity values which should not be highlighted are scaled to a low value 30. Fig. 1.5 represents the intensity slicing. Fig. 1.6 shows the usefulness of intensity slicing in highlighting
2
FIG. 1.2: a)Input image b)Gamma mapped image with γ =.5
FIG. 1.3: Contrast stretching [4]
FIG. 1.4: a)Input image b) Contrast stretched image 3
FIG. 1.5: Intensity slicing
specific regions of the image. Histogram equalisation: The histogram of an image with intensity levels in the range [0, L − 1] is a discrete function h(rk ) = nk where rk is the kth intensity value and nk is the number of pixels in the image with intensity rk [4]. Histogram can be converted into a probabilistic measure by dividing the number of pixels at each intensity level by the total number of pixels i.e for a M*N image, number of pixels at intensity rk can be nk /M N . Histogram of a dark image is skewed in the lower region of the histogram and similarly histogram of a bright image is skewed in the higher region of the histogram. These images are having low contrast. The histogram of a high contrast image covers a wide range of intensity values and follows a uniform probability distribution. So to increase the contrast of any image, one has to convert its histogram into a uniform histogram and this process is called histogram equalisation. Gonzalez [4] gives a good mathematical treatment of histogram equalization in his book. 1.1.1.2
Neighbourhood Processing Techniques
Till now we have seen enhancement techniques that act on individual points. Let us consider now some techniques that take into account the neighbouring pixels too.
4
FIG. 1.6: a)Input Image b)Output image with some highlighted intensities
FIG. 1.7: a)Input Image b)Histogram equalised image
5
FIG. 1.8: Unsharp masking
Histogram equalisation(HE) works well on images that have unimodal histograms(dark or bright images). But HE does not work on scenes with bimodal histograms(scenes that have both dark and bright regions). To overcome this disadvantage Adaptive Histogram Equalisation(AHE) [8] was introduced. AHE gives local contrast enhancement by applying HE in windows of adjustable size. Noise enhancement is the disadvantage of AHE. This disadvantage is overcome by a new algorithm called Contrast Limited Adaptive Histogram Equalisation(CLAHE) [9]. Unsharp masking and highboost filtering are two other important neighbourhood pixel dependent enhancement techniques . In unsharp masking, the image is blurred by convolving the image with a low pass gaussian filter. This blurred image is subtracted from the original image and a mask image is obtained. This mask image is added to the original image as g(x1 , x2 ) = f (x1 , x2 ) + k ∗ gmask (x1 , x2 ) where f is the original image, gmask is the masking image and g is the final sharpened image. When k = 1, this enhancement process is called as unsharp masking. When k > 1, the enhancement process is called highboost filtering. Fig. 1.8 explains the unsharp masking algorithm. 1.1.2
Frequency domain enhancement techniques
Till now we have seen spatial domain enhancement techniques. Let us see now homomorphic filtering, an important frequency domain enhancement method. Every pixel in an image can be represented as a product of illuminance and reflectance. This mode can be used for simultaneous dynamic range compression and contrast enhancement. It is not possible for us to separate the illuminance and reflectance component from each other in each pixel. Though we cannot separate the components, the advantage with this method is that the filter H can act separately on both the illuminant and reflectance components. The illuminance is generally slow varying and hence the low frequencies of the image are associated with it. The reflectance is 6
FIG. 1.9: Homomorphic filtering
fast varying and the high frequencies are associated with it. Filter H is designed in such a way that it can simultaneously attenuate the low frequencies and amplify the high frequencies. Fig. 1.9 explains the homomorphic filtering algorithm. 1.2
COLOUR IMAGE ENHANCEMENT TECHNIQUES
These techniques discussed above are primarily for gray scale images.Let us discuss how these techniques are applied on color images. 1.2.1
Intensity Transformations
Intensity transformations can be applied on Red-Green-Blue(RGB) or CyanMagenta-Yellow(CMY) or Hue-Saturation-Intensity(HSI) color space. The formula of transformation is color space dependent. More about color models are dealt in Chapter 2. 1.2.2
Histogram equalisation
If histogram equalisation is applied on each of the three color bands then the enhanced image has unnatural colors. Since histogram equalisation deals with intensity values, an image must be converted to HSI space and the intensity channel alone be histogram equalized. This procedure preserves the original colours. Many other gray scale image enhancement techniques such as homomorphic filtering, unsharp masking can be applied on the intensity channel of the color images alone instead of applying on the RGB color space. All these techniques are basically used for contrast enhancement. But they cannot simultaneously achieve dynamic range compression, contrast enhancement and true color rendition. Retinex can help us achieve all these three properties simultaneously. The principle is same as homomorphic filtering. But Retinex is able to achieve dynamic range compression because an exponential function is not used 7
in the final stage of enhancement. Instead of exponential function a non linear color restoration function is used. This function is able to mimic the exponential function and at the same time does not compromise on dynamic range compression. 1.3
APPLICATIONS OF IMAGE ENHANCEMENT
Image enhancement is being applied in various fields that are listed below. i. Medical Imaging uses Image enhancement for removing noise and sharpening the details [10], [11], [12], [13]. In medical imaging, minute details play a critical role in diagnosis of diseases. So this makes enhancement an essential pre processing step in medical image processing as shown in Fig. 1.10. ii. Image enhancement is also used in restoration of old paintings and historical photographs [14], [15]. Fig. 1.11 is a good example of this enhancement. iii. Image enhancement also helps partially visually impaired people to read fine print in books, computers and television. In this field a good amount of research has been carried out [16], [17], [18]. iv. Image enhancement is also used in forensic science. Images obtained from Closed Circuit Television(CCTV) cameras, crime scenes are enhanced to help security agencies in nabbing the culprits [5], [19]. Fig. 1.13 is a very good example of finger print image enhancement that will be useful in investigation of crimes [5]. v. Another interesting and important application is remote sensing. On days of atmospheric turbulence the images of earth taken from satellites look hazy and completely blurred. Various algorithms have been proposed to dehaze these images and restore the objects in the image [20]. vi. In oceanography, images gives us information about the flora and fauna, ocean floor, water flow, sediment concentration [21], [22]. vii. Research is also being done now in the application of image enhancement in e-learning. Image enhancement here is mainly used to enhance the readability of content that is written by the teacher on the black board [23], [6]. Fig. 1.12 is a very good example of enhancement that can be used in e-learning [6].
8
FIG. 1.10: Medical image enhancement
FIG. 1.11: Image enhancement of old images [5]
FIG. 1.12: Whiteboard [6]
9
FIG. 1.13: Forensic based image enhancement
1.4
THESIS FLOW:
Chapter 2 provides a detailed literature review of existing image fusion, haze removal, tone mapping and Retinex based image enhancement methods. In Chapter 3 we discuss a fusion based approach for MSRCR. In Chapter 4 we propose an automated approach that obtains parameters from the image under consideration. In Chapter 5 we propose a Retinex based haze removal method . In Chapter 6 we propose a Retinex based tone mapping approach. In Chapter 7 results of all the proposed algorithms are discussed in detail. 1.4.1
Fusion based Multi Scale Retinex with Color Restoration
Image fusion is the technique of fusing two images after aligning them. These two images could have been taken in two different lighting conditions, using two different sensors, at two different times or even at two different angles. Our objective is to fuse the two images in such a way that the resulting image is better than these two images. We can fuse the Single Scale Retinex outputs obtained for different surround constants instead of adding them. In this thesis we discuss Li Tao’s pixel based method [24] of fusing two images using the wavelet transform coefficients. The wavelet used is Dual Tree Complex Wavelet Transform(DT-CWT).
10
1.4.2
Automated Multi Scale Retinex with Color Restoration
In the Single Scale Retinex method the final output has a washed out appearance. Zia et al. [1] themselves have proposed an approach using the histogram of the Retinex enhanced image. Their approach is to clip of the extreme values of the histogram and scale the rest to 0-255. But the problem with this approach is the authors have not specified how to choose the clipping points. We propose in this thesis an image independent approach for choosing the clipping point. We use the frequency of occurrence of pixels in the image as a control measure to find this clipping point. 1.4.3
Retinex based Haze Removal method
Images of outdoor scenes are degraded by the atmospheric particles such as haze, fog and rain. It is important to develop restoration algorithms especially in applications like automatic car navigation or automatic landing systems because these systems are based on the assumptions that there are no atmospheric disturbances. In this thesis we shall consider the problem of haze removal. A widely used model for haze formation is I(x) = R(x)t(x) + a(1 − t(x)) where x is a pixel location, I is the observed image, R the scene radiance, a is the atmospheric light or airlight, and t is the transmission coefficient. The image observed is a combination of an attenuated version of underlying scene radiance with an additive haze layer where atmospheric light represents the color of the haze as shown in Fig. 1.14 which was obtained from Maltin’s Thesis [25]. The goal of any haze removal algorithm is to find R for which one needs the knowledge of a and t. This is an under-constrained problem. In a grey scale image, each pixel has only 1 constraint but 3 unknowns. In a color image, each pixel has 3 constraints and 7 under knowns(assuming t is same for each color channel). One of the latest haze removal algorithms proposed by Kaiming et al. [3] uses a dark channel prior based approach for haze removal. Though this approach gives very good results, this method is computationally complex. In this thesis we propose a RETINEX based approach that gives good results and is also computationally simpler. 1.4.4
Tone Mapping
The advancement in image acquisition technology has lead to HDR imaging which helps us to capture the entire dynamic range of a scene. The color rendition of HDR
11
FIG. 1.14: Haze formation model [3]
imaging is also very good. But HDR display devices have not developed in the same way and still remain costlier. So the challenge is to display these HDR images on low dynamic range display devices without compromising on the overall quality of the image. This operation is called tone mapping. K. Kim et al. [26] have proposed in their paper a novel Retinex based tone mapping technique. According to the author when Multi Scale Retinex is used for tone mapping, bright regions are over highlighted at the expense of darker regions. He proposes some minor modifications to the MSR algorithm and achieves better tone mapping. In this thesis we combine the original Retinex algorithm and gamma mapping and achieve good tone mapping.
12
CHAPTER 2 LITERATURE SURVEY The literature review mainly focusses on previous research carried out in the fields of image fusion, haze removal and tone mapping. A historical perspective of Retinex algorithm is also presented here. 2.1
RETINEX ALGORITHMS
The Retinex algorithm can be implemented either by the Path version [27] or the Surround-Centre based version [1]. 2.1.1
Path version
The algorithm calculates subsequent additions of pixel differences along a set of onedimensional random paths contained in the image [28]. A threshold operator T is used to remove the effect of illuminance. The image computed by this algorithm along a path τ at a particular pixel x is given by Rcτ (x) =
X
T [log[Ic (k+1)] − log[Ic (k)]]
(2.1)
kτ,k
where Ic (k) is the intensity at pixel k for color channel c and T is a threshold operator. The operation k < x is true if pixel k is located before pixel x along the path τ . Threshold operator : If x < threshold then, T [x] = 0
(2.2)
T [x] = 1
(2.3)
If x > threshold then,
Finally, the Retinex image is averaged across different paths M 1 X τ Rc (x) = R (x) M τ =1 c
(2.4)
where Rc (x) is the reflectance estimate at pixel x in color channel c and M is the number of paths. Rizzi et al. [29] had developed a path based algorithm of RETINEX using Brownian motion. The idea was first proposed in 1993 by Rizzi et al. [30] and 13
is inspired by the results on neuro physiological research. The drawback of the path based method is it does not converge always to a stable solution because convergence depends on the choice of paths. On the contrary surround-centre based method is a very stable algorithm. 2.1.2
Surround-Centre version
Surround-Centre Version of Retinex was proposed as Single Scale Retinex by Zia et al. [1]. Then later Zia et al. proposed a modified version of Single Scale Retinex known as Multi Scale Retinex with Color Restoration [2]. The Surround-Centre version of Retinex is computationally simpler than Path method and hence popularly used today. A color constancy algorithm must be able to simultaneously achieve the 3 properties given below [1], i dynamic range compression ii color independence from the spectral distribution of the scene illuminant and iii color and lightness rendition [1] The first property can be achieved by applying logarithmic transformations on the image [4]. The second property can be achieved by eliminating the illuminance component in the image. Every pixel in an image can be represented as a product of illuminance and reflectance i.e. S(x) = R(x) ∗ L(x)
(2.5)
where L represents illuminance, R represents reflectance and S represents the image pixel. Our aim is to eliminate L(x, y). Illumination varies slowly across the image unlike reflectance. So illuminance of an image can be obtained by low pass filtering the image. Instead of obtaining R = S/L, we use logarithmic approach to achieve the same, since applying logarithm on an image gives us dynamic range compression. Let s = log(S), r1 = log(R), l = log(L). So now equation 2.5 can be represented as, r1(x) = s(x) − l(x)
(2.6)
L can be obtained by convolving a low pass filter F with image S. Initially Land p proposed 1/r2 as the low pass filter F , where r = (x21 + x22 ), x1 and x2 are the 14
FIG. 2.15: a) Hurlbert’s surround b) Land’s surround c) Moore’s surround
pixel locations. By using this function, the first two properties of color constancy were achieved but not the third one. Moore et al. proposed an exponential “absolute value” F (x) = e−|r|/c [31]. Hurlbert proposed a gaussian function F (x) = e−r
2 /c2
[32]. These 3 models are shown in Fig. 2.15 from which we can infer that the inverse square rolls rapidly but at distant pixels it has higher response than the other two forms. Similarly at distant pixels exponential function has a higher response than gaussian function. Gaussian function acts more locally, exponential function is less local and the inverse square function is global in nature [1]. Mathematically, Retinex equation can be represented as, Ri (x) = log[Ii (x)] − log[F (x) ∗ Ii (x)]
(2.7)
where Ii (x) is the image distribution in the ith color spectral band, ‘∗’ denotes the convolution operation, F (x) is the surround function and Ri (x) is the associated retinex output. The retinex operation is performed on each spectral band. For the gaussian surround function, when the surround constant is varied from a small value to a large value, dynamic range compression is sacrificed for improved rendition as shown in Fig. 2.16. The middle of the range(50 < c < 100 pixels) represents a reasonable compromise, where shadows are fairly compensated and rendition achieves acceptable levels of image quality [1]. To compensate for this disadvantage of single scale retinex i.e. it can either achieve good dynamic range compression or good color rendition, we go for multi scale retinex method which helps us to achieve both simultaneously.
15
FIG. 2.16: Single scale retinex outputs a) c=15 b) c=80 c) c=250
2.1.2.1
Multi Scale Retinex(MSR)
In multi scale retinex we find the retinex output according to equation(1) for various values of surround constants and add all the retinex outputs by giving them equal weight according to the equation RM SRi =
N X
wn Rni
(2.8)
n=1
where N is the number of scales,Rni is the ith component of the nth scale,RM SRi is the ith spectral component of the MSR output,and wn is the weight associated with the nth scale [2].The only difference between R(x1 , x2 ) and Rn (x1 , x2 ) is that the surround function is given as Fn (x) = ke−r
2 /c2 n
(2.9)
Experimentation has proved that it is enough to use just 3 scales viz.one small scale(cn < 20) one large scale(cn > 200) and third one as an intermediate scale value.The three different retinex outputs that we get because of using the gaussian surround function wit three different scaling functions,is equal weighted and added to get the final retinex output. 16
2.1.2.2
Multi Scale Retinex with color restoration(MSRCR)
Retinex algorithm was originally developed for images that do not violate the ’greyworld’ assumptions. If the reflectances of the image in all the three color bands are same on an average then the image is said to obey grey-world assumption.The general effect of retinex processing on images with regional or global grey-world violations is ”greying out” of the image, either globally or in specific regions. This desaturation of color in some cases can be severe. So we have to find a good color restoration method that provides good color rendition for images that contain grey-world violations [2]. But one must not compromise on color constancy in our pursuit of color rendition, as color constancy is one of the prime objectives of retinex. The algorithm for color restoration is given below, RM SRCRi (x) = G[Ci (x)RM SRi (x) + b]
(2.10)
where 0
Ci (x) = f [Ii (x)]
(2.11)
is the ith band of the color restoration function(CRF) and RM SRCRi is the ith spectral band of the multiscale retinex with color restoration. The function that provides us the best overall color restoration is, 0
Ci (x) = βlog[αIi (x)] = βlog[αIi (x)] − βlog[
S X
Ii (x)]
i=1
where β is a gain constant, α controls the strength of the non linearity, G and b are final gain and offset values. The values specified for these constants by Zia et al.[1] are β = 46, α = 125, b = −30, G = 192. 2.1.3
Luma based MSR
Before discussing luma based MSR let us briefly discuss about color models. The colors perceived by humans are determined by the nature of light reflected from the object. Visible light belongs to a narrow band of frequencies in the electromagnetic spectrum (400-750 nm wavelength). An object appears white to the observer if it reflects light that is balanced in all visible wavelengths. Three quantities are used 17
to describe the quality of chromatic source are radiance, luminance and brightness. Radiance is the total amount of energy that flows from the light source. Luminance gives a measure of the amount of energy an observer perceives from a light source. Eg. An infrared source could have significant energy(radiance) but we cannot perceive it i.e luminance is zero. Brightness is a subjective measure and cannot be measured [4]. Around 6-7 million cones are present in the human eye. Approximately 65 percent of cones are sensitive to red light, 33 percent to green light and 2 percent are sensitive to blue. But these 2 percent cones are highly sensitive to blue colour. Due to the nature of the sensitivity of the cones, all colors are seen as a linear combination of primary colors Red, Green and Blue. The primary color can be added to produce the secondary colors of light such as magenta(red+blue), cyan(green+blue), yellow(red+green). Mixing the three primaries or a secondary with its opposite primary color in the right proportion to produce white light. The characteristics used to differentiate between colors are brightness, hue and saturation. Brightness characterises the achromatic notion of the intensity. Hue represents the dominant color as perceived by an observer. Saturation refers to the amount of white light mixed with a hue. Pink(red and white) and lavender(violet and white) colours are less saturated i.e degree of saturation is inversely proportional to the amount of white light added. Hue and saturation together form the chromatic component. So a color can be characterised by its brightness and chromaticity [4]. The amounts of red, green and blue needed to form any particular colour are called the tristiumlus values X, Y and Z. A color is then specified by its trichromatic coefficients defined as x=
X X +Y +Z
Y X +Y +Z Z z= X +Y +Z x+y+z =1
y=
2.1.3.1
(2.12) (2.13) (2.14) (2.15)
Color Models
The purpose of a color model is to facilitate the specification of colors in a standard way that can be followed across the world. A color model is a specification of a coordinate system where each color is represented by a unique point. The popular
18
color models are RGB(primary color model), CMY(secondary color model) and HSI model which corresponds closely with human perception of colours. The advantage of HSI is that it decouples the chromatic and achromatic space and we can apply all gray scale algorithms on the achromatic space and get good results. There are many color models apart from HSI which gives us a good separation of luma-chroma space. Some examples are YCbCr, YUV, YIQ, YPbPr, HSV. These models have been developed and successfully deployed for different applications. Let us see some transformation equations below. CMY Model:
C
1
R
M = 1 − G Y 1 B
(2.16)
YUV Model: The YUV format is used for analog National Television System Committee(NTSC) video but not for digital format. R .289 .587 .114 Y U = −.147 −.289 .436 × G B .615 −.515 −.100 V
(2.17)
RGB to YUV transformation for analog TV Y 1 0 1.14 R G = 1 −.395 −.581 × U V 1 2.032 0 B
(2.18)
YCbCr Model: For digital component video the color format YCbCr is used. For standard definition TV(SDTV). This standard has been specified by International TeleCommunication Union(ITU)
Y
0
.299
0.587
.114
R
× G Cb = 128 + −.169 −.331 .5 B 128 .5 −.419 −.081 Cr R 1 0 1.4 Y G = 1 −.343 −.711 × Cb − 128 B 1 1.765 0 Cr − 128 19
(2.19)
(2.20)
YCbCr model is used for High Definition Television(HDTV) too. Y 16 .183 0.614 .062 R Cb = 128 + −.101 −.339 .439 × G Cr 128 .439 −.399 −.04 B R 1.164 0 1.793 Y − 16 G = 1.164 −.213 −.533 × Cb − 128 B 1.164 2.112 0 Cr − 128
(2.21)
(2.22)
HSV model: Since Retinex tries to model Human Visual System(HVS), in Luma based MSR we use the Hue-Saturation-Value(HSV) Space as HSV too exactly models HVS. Below is given a RGB to HSV conversion. The Saturation component is given by, S =1−(
3 )(min(R, G, B)) R+G+B
(2.23)
The Intensity/Value component is given by, 1 I = (R + G + B) 3
(2.24)
The Hue Component is given by, H=
−1
θ = cos (
[(R
θ if B ≤ G
(2.25)
360 − θ if B > G
(2.26)
1 [(R − G) + (R − B)] 2 ) − G)2 + (R − B)(G − B)]1/2
(2.27)
The HSV to RGB conversion is given below. RG Sector:(0◦ ≤ H < 120◦ ) B = I(1 − S) R = I[1 +
Scos(H) ] cos(60◦ − H)
G = 3I − (R + B)
(2.28) (2.29) (2.30)
GB Sector:(120◦ ≤ H < 240◦ ) H = H − 120◦ 20
(2.31)
R = I(1 − S) G = I[1 +
Scos(H) ] cos(60◦ − H)
B = 3I − (R + G)
(2.32) (2.33) (2.34)
BR Sector:(240◦ ≤ H ≤ 360◦ ) H = H − 240◦
(2.35)
G = I(1 − S)
(2.36)
B = I[1 +
Scos(H) ] cos(60◦ − H)
R = 3I − (B + G) 2.1.3.2
(2.37) (2.38)
Algorithm
Though MSRCR performs well in most of the cases, it fails to render exact colour when there are large constant areas of constant colour. This is because, the “greyworld” assumption of MSRCR will be violated when there are large areas of constant colour. Applying MSRCR algorithm on the luma component of the image turned out to be a good solution for this problem. The PCA is the best transform from RGB space to luma and chroma space since the luma and chroma components are completely independent(orthogonal) of each other. But PCA is not widely used since the transform is image dependent and for bigger images, it will be computationally more intensive than other color space transforms. Meylan et al., have compared in their work [33] the difference in the outputs obtained when MSRCR is applied only on the luma components that are obtained from different colour transforms. According to this paper, the outputs obtained after applying MSRCR on the luma components of YUV and Lab transforms is as close as the outputs obtained after applying MSRCR on the luma component of PCA transform. Figure 2.17 shows the disadvantage of using the original MSRCR and how luma based MSRCR overcomes it. Even color rendition is better in luma based MSRCR for this particular example. Another advantage of luma based MSRCR is the simulation time is reduced approximately to one-thirds of the time it took for the original MSRCR since here the RETINEX algorithm is applied only on the luma component.
21
FIG. 2.17: a) Input image b) Original MSRCR output c) Luma based MSRCR output
2.2
IMAGE FUSION
Image fusion is the process of combining images to make a better enhanced image. Image fusion can be done at pixel or feature level [34]. Pixel level fusion is the fusion of a window of pixels from various images [35]-[36]. Feature based fusion needs extraction(by segmentation) of features contained in the images [37]-[38] as shown in Fig. 2.18. The features are decided by properties such as size, shape and texture. Before fusing two images that are taken from different angles or sensors it is important to align these two images properly and this is called as image registration. The advantage of feature based fusion is it does not require strictly registered images. Feature based fusion algorithms also avoid drawbacks of pixel based fusion like noise enhancement, blurring. But the biggest disadvantage of feature based fusion algorithms are they are computationally complex and hence we use a pixel based fusion algorithm in this thesis. Arithmetic fusion can be given by the expression F (x1 , x2 ) = ka A(x1 , x2 ) +
22
kb B(x1 , x2 ) where A and B are the source images and F is the fused image. Image averaging is the simplest fusion method where ka and kb =1/2 [39], [40] as shown in Fig. 2.19. But this fusion method results in loss of contrast. This method works well when the image regions are similar and worse when image regions are negative of each other. Yamamoto et al. successfully employed this fusion method for the detection of navigation obstacles [40]. Pixel based fusion in spatial domain often leads to reduction in contrast. But with introduction of pyramid transforms better results are obtained. The idea is to decompose the images to be fused into many layers of different resolutions and fuse each layer of an image with its counterpart layer from other images and finally obtain a fused image by reconstructing from these fused layers. Burt and Adelson were the first to propose a Multi Resolution Analysis(MRA) representation of images [41]. They introduced MRA as an image compression algorithm. But MRA has grown beyond its initial image compression role because of Mallat’s idea of representing wavelet transform in the form of a MRA representation. Let us discuss here about Burt’s Laplacian-Gaussian pyramid. A MRA pyramid decomposes an image into multiple resolutions at different scales. The lowest scale has the highest resolution and is of the same size of the image. The original image or the lowest scale h0 is low pass filtered to get the next scale h1 and then to get h2 and so on. The number of levels is application dependent and can be decided by the user. Since the low pass filter used is Gaussian , this pyramid is called Gaussian pyramid. Burt proposed in 1984 [42] that the pyramid can also be used for image fusion. This is the first paper that proposed image fusion using the MRA coefficients of an image. Toet presented a similar fusion algorithm using Ratio of Low pass(ROLP) pyramid [43]. ROLP is a pyramid which at every level of decomposition is the ratio of two successive gaussian pyramids. Toet proposed that contrast of the scene plays a very important role in human visual system and therefore a fusion technique that selects the highest local luminance contrast is likely to provide better details. Burt and Kolcynzski [44] proposed to fuse the coefficients of a gradient pyramid using a combination algorithm that uses both activity and match measure. We shall see these two terms in the later part of the discussion. Though pixel based fusion algorithms using the pyramid MRA have performed better than ordinary fusion algorithms, wavelet MRA based fusion algorithms provide us more benefits. The wavelet based algorithms are better than pyramid based ones for these reasons. i Gaussian based pyramids are 4/3 of the original image size and hence redundant. 23
FIG. 2.18: Feature based fusion
FIG. 2.19: Pixel based fusion
24
FIG. 2.20: Pyramid based fusion
Wavelets because of the orthogonal basis gives non redundant representation. ii The computational complexity of the Gaussian pyramid is far more than wavelet transform which can be easily implemented using filters. iii The reconstruction in wavelet transform is also simple and can be easily executed using dual filters to the one used in analysis side. Fig. 2.20 and Fig. 2.21 gives an algorithm of how pyramid and wavelet based algorithms work. 2.2.1
Wavelet
Let us see here the basics of wavelet before discussing the wavelet based fusion systems. The wavelet transform uses orthogonal basis to achieve multi resolution decomposition of images. It is a non redundant representation of the image. The basis functions are generated from the mother wavelet by translations and dilations. A wavelet can be represented as MRA by the Fig. 2.22. Fig. 2.22 and 2.23 represent the analysis and synthesis filters(reconstruction filter) for a one dimensional signal. In Fig. 2.22 H0 represents a high pass filter and H1 represents a low pass filter. The outputs of these two filters represents the low pass and high pass components of 25
FIG. 2.21: Wavelet based fusion
FIG. 2.22: Analysis 1-d filter bank
FIG. 2.23: Synthesis 1-d filter bank 26
FIG. 2.24: Analysis 2-d filter bank
the original signal at the first level of decomposition. After these two filters a downsampling by 2 is performed so that at all levels we have same number of elements as in the original signal. Similarly many more levels of decomposition can also be done. But in wavelet transform each time only the low pass filtered components are further decomposed. There is a another transform called wavelet packet transform in which both high pass and low pass filtered components are decomposed at each stage. In Fig. 2.23 F0 represents a high pass filter and F1 represents a low pass filter. First we have to up-sample the signals since there is a down-sampler in the analysis side. And then pass these up-sampled signals through the low pass and highpass filters respectively and finally add them to get the output signal. For perfect reconstruction to be possible, the conditions on filters in z domain are F0 (z) = H1 (−z) and F1 (z) = H0 (−z). For images(2-d signals), the wavelet is also a 2-d transform. We can apply the wavelet on rows and columns separately as shown in Fig. 2.24. After first level of decomposition image is divided into 4 sub bands as shown in Fig. 2.26 . After two levels of decomposition we get 8 sub bands as shown in Fig. 2.27 . 2.2.2
Wavelet based fusion
Wavelet based fusion outputs are better perceived by humans [45], [46]. In fusion there are certain requirements like
27
FIG. 2.25: Synthesis 2-d filter bank
FIG. 2.26: Wavelet decomposition at first level
28
FIG. 2.27: Wavelet decomposition at second level
i The fused image should preserve all the relevant information contained in the source images. ii The fusion process should not introduce any artefacts. iii Irrelevant features and noise must be attenuated. These three properties are very well satisfied by Wavelet based fusion algorithms. Wavelet based fusion can be expressed as I(x, y)=W −1 (φ(W (I1 (x, y)), W (I2 (x, y)))) where φ represents the fusion rule, W represents the wavelet transform and W −1 represents the inverse wavelet transform. All the wavelet based fusion algorithms differ only in the fusion principle φ. Li et al. [46] proposed that at each level and position of wavelet transform domain, maximum selection rule is used to determine which of the inputs contains more information. Burt et al. [41] combined coefficients in a weighted manner and these weights were found by using the local activity measure in each sub band. Wilson et al. [45] combined the wavelet coefficients in a weighted manner. The weights are decided based on the contrast measure.
29
Conventionally researchers have used discrete wavelet transform (DWT) for fusing images. But discrete wavelet transform inherently has the disadvantage that the transform is shift variant i.e. if the input samples shift by a few samples, the wavelet coefficients do not simply shift but their magnitudes alter. This property of shift variance is a great disadvantage in fusion which can be overcome by shift invariant wavelet transforms. Real transforms are shift variant because of the sub-sampling present in the transform. Removing this sub-sampling did give better fusion results than DWT. This transform is called as shift invariant discrete wavelet transform(SIDWT) [47]. But once dual tree complex wavelet transform (DT-CWT) was found [7], it has been widely used since then, as it gave better results than both DWT and SIDWT. The reason for former is shift invariance and the reason for latter is that DT-CWT has better directional selectivity than SIDWT. Various image fusion algorithms using wavelets are compared by Paul et al. [48]. In the case of wavelet transform fusion, the corresponding wavelet coefficients from the input images are combined using the fusion rule. The 3 popular fusion rules are discussed below.
i Maximum Selection(MS) Fusion: At each level and each point of the transform, the fused coefficient is the maximum one from all the source images. ii Weighted Average(WA) Scheme: The scheme developed by Burt et al. [41] used a normalised correlation between the two images’ subbands over a small local window. The resultant coefficient is found by using a weighted average of two images’ coefficients. iii Window based Verification(WBV) Scheme: Creates a binary decision map to choose between each pair of coefficients [46]. By using a quantitative measure to compare fusion schemes as proposed by Li et al. [46], DT-CWT when used along with WBV is found to be the best fusion algorithm by Paul et al. [48]. We have used here one such fusion method called 30
DT-CWT based Image Fusion which has been proposed by Li Tao [24] in her Ph.D dissertation. 2.3
HAZE REMOVAL
Now let us look into some models that are used for haze removal. Among current haze removal research, haze estimation methods can be broadly divided into two categories i Methods that rely on multiple images/additional information. ii Single image based method. Here are some methods that rely on additional information: i Multiple images of the scene taken using different degrees of polarisation. ii Multiple images taken under different weather conditions. iii Methods that require user supplied depth information or 3D model. Polarising filter based algorithm uses two or more images for haze removal. The algorithm works on the principle that ‘airlight’ of the image is polarized [49]. We will see the term ‘airlight’ in Chapter 5. Another algorithm by Narasimhan et al. [50] for haze removal uses multiple images taken under different weather conditions. Haze removal algorithms using apriori information needs either user supplied depth information [51] or a 3D model [52]. Recently single image based haze removal models have been proposed. These models use the haze formation physical model that has been developed by Koschmieder [3]. Narasimhan et al. [53] describe in detail the physics behind the atmosphere optics and the haze formation model. In 2008, Tan proposed a single image based haze removal algorithm [54]. He observed that a haze free image has higher local contrast than a hazy image. So he was able to obtain haze free image by maximizing the local contrast in the hazy image. But this approach is not based on haze formation physical model and hence the haze free images appear unnatural due to over saturation of colours. Tarel proposed a fast implementation of Tan’s algorithm [55]. Recently Kaiming et al. [3] proposed a dark channel prior based approach. This is a single image based haze removal method based on the assumption that for a haze free image, the minimum of the three colour bands R,G 31
and B will be zero. This method holds good for a wide range of hazy images. This approach gives us good results but it is computationally complex. In this report we propose a computationally simpler Retinex based haze removal method. 2.4
TONE MAPPING
The illuminant sources in a natural environment can be numerous, varying from a night scene to a early morning sunrise and a noon bright sunlight and so on [56]. Human Visual System can adapt itself to different illuminant conditions and we will be able to see all the details in the scene clearly [57], [58]. Our cameras cannot do this and the image that we capture are heavily dependent on illuminant conditions. When the dynamic range of the scene is high when compared to the dynamic range of the capturing device, then the scene is called a high dynamic range(HDR) scene. The dynamic range of several natural scenes were measured [59] and 1:160 was found as the average dynamic range of the scene. A HDR scene’s dynamic range will be of the order 1:1000. The dynamic range of traditional digital cameras is less and hence cannot capture all the details of the HDR scene. But today we have HDR cameras that have good dynamic range and can capture the scene as perceived by human eye. The cost of HDR cameras is reducing with the advancement of VLSI technology. But the advancement in display device technology has not been as fast as in image acquisition technology. Even display devices that can cater to HDR images are costly and hence less accessible. Still majority of our display devices remain low dynamic range(LDR). So when the images captured by HDR cameras are displayed on LDR devices, not all details of the scene are visible. The algorithm that helps us display HDR images on LDR devices without loss of details is called tone mapping. Tone mapping algorithms can be either local or global. Global tone mapping algorithms are spatially invariant i.e same function is applied on all the pixels. The usual global mapping operators are logarithmic, gamma and sigmoid mapping functions. Local tone mapping algorithms act differently at different spatial locations and are neighbourhood dependent. Global tone mapping methods are not good for HDR images because they do a linear compression of the dynamic range and hence there is a loss of lots of details in the tone mapped image. So a local tone mapping algorithm is indeed essential to improve the local contrast.
32
2.4.1
HDR History
The HDR history has begun right from the days of renaissance literature. Leonardo da Vinci introduced high contrast painting in the 15th century [60]. Photographic HDR techniques have been used from the mid 19th century. Multiple exposure techniques were used as differently exposed films were cut and composited in the development stage to create a HDR print of the scene. Ansel Adams developed dodging and burning techniques in the development of the print and hence captured the HDR scenes with details visible in both light and dark regions. Digital HDR imaging research started in 1980’s. In 1985 Sony patented a multiexposure CCD camera with one lens, beam splitter and 2 CCD image regions [61]. In 1987, Polaroid patented a camera that combines two different exposures to increase the dynamic range [62]. The multi exposure blending techniques were first developed by Mann et al. in 1995 [63] and later refined by Debevec et al. in 1997 [64].
2.4.2
Dynamic Range of a scene
Human Visual System can adapt to any dynamic range. HVS adapts to different luminance values by using a variable aperture i.e the pupil in the eye. The pupil can expand from 2mm to 7mm in diameter. The different cells in the retina determine how HVS adapts to different luminance levels. The rods are used for high luminance levels and cones for low luminance levels. The HVS state when only rods are used are called scotopic vision. When cones are used the vision is called photopic vision. When both rods and cones are used HVS state is called mesopic. Only when cone cells are active we can see the colour. So people for whom only rod cells work, then color vision will be blind for them [65]. HVS has a nonlinear logarithmic response to light intensity. High dynamic range imaging has varied applications like medical imaging, computer vision, space imaging, remote sensing. The basic idea of multi exposure image capture is to generate a HDR image from differently exposed LDR images if the dynamic range of the scene exceeds that of the capturing device. By varying the exposure and taking multiple images of the scene, different parts of the scene can be combined to a cohesive HDR image from differently exposed single images [66]. Mann et al. [63] were the first to propose that
33
differently exposed LDR images could be used to get a HDR image. Debevec et al. [64] and Mitsunanga et al. [67] refined the initial method proposed by Mann further. Many solutions have been proposed for the problem of misalignments caused due to camera shakes. Since the images have different image exposures the methods proposed to overcome the misalignments must be insensitive to intensity of the image. Ward et al. [68] uses median threshold and multi scale pyramid for determining the misalignments. Mantuik et al. [69] have proposed a Scale Invariant Feature Transform based method to find the misalignments. Image alignment methods correct the camera movements but often multi exposure images also include objects moving between images called ghosting. Kang et al. [70], Khan et al. [71], Pedone et al. [72], Gallo et al. [73], Sidibe et al. [74] have proposed various methods to overcome this ghosting problem. 2.4.3
Tone mapping operators
The aim of tone mapping operators is that HDR scene luminances must be satisfactorily mapped to LDR display devices like CRT and LCD. 2.4.3.1
Global operators
Schlichk’s [75] uniform rational quantization is a simple operator with two user controlled parameters. It includes a user controllable black level parameter of display δL0 and another user controllable parameter k. He also had experimented with local processing and segmentation based variations of the operator but found that global version performs better. The operator is designed for gray scale and YUV images. Ferschin et al. [76] developed an exponential compression operator. In this method an exponential function along with the average luminance of the image is used for tone mapping. Logarithmic mapping algorithm by Drago et al. [77] uses the assumption that the HVS has a logarithmic response to light intensities which helps preserving contrast and detail. The algorithm uses two controllable user parameters. One parameter is maximum display luminance Ld,max and another is constant parameter p.
34
2.4.3.2
Local operators
Ashikminn’s [78] operator follows the functionality of a HVS. It is a local operator with three stage processing. In the first stage, local adaptation luminance is estimated. The values are compressed into a display range in second stage. Some of the details lost are reintroduced in the third stage. Gradient domain processing is used in tone mapping. This algorithm is based on how HVS operates at local contrast level. Horn et al. [79] used the gradient image in determining the lightness and manipulating the gradients in producing a tone mapped image. Fattal et al. [80] uses gradient domain compression with user controllable parameters to select the level of amount of compression and amount of detail in the image. Gradient domain operators with compression produces images with exaggerated small details. Kim et al. [26] have used a modified Retinex based approach for tone mapping. In this thesis we explore a Retinex based approach of tone mapping.
35
CHAPTER 3 A FUSION BASED APPROACH FOR MSRCR Image fusion is a technique of combining two or more images so that the combined image is better enhanced than the combining images. In this section let us explore how fusion can be used in the MSRCR algorithm. In order to properly fuse the two images we need to study and understand the three properties that will exist between the two images. They are, i) Complimentary Information: Some scene information might appear in one image, but not in other. ii) Common information: The scene information that appears in both the source images. But these scene information need not look same in both the images. iii) Properties: The two images may vary in the properties i.e the sensors used to capture the two images may have different sensing capability, dynamic range, resolution and noise [24]. 3.1
ALGORITHM
In the Multi Scale Retinex(MSR) method as proposed by Zia et al. [1], the MSR output is given as a sum of three SSR images as according to the equation 2.8. We propose to fuse these three SSR images instead of adding them to get better enhanced image than the original method both qualitatively and quantitatively. Let us see the two terms that are used in this fusion algorithm. Match Measure: Match measure means the similarity between the two images in the spatial domain. Activity measure: Activity measure computes the pixel information from the MR(Multi Resolution) decomposition of the two sources. Steps: i)Obtain match measure from the spatial domain info of the two images. 36
FIG. 3.1: DTCWT based fusion
ii)Take DT-CWT of both the images and then compute the activity measure from the wavelet domain coefficients. iii) The match measure and activity measure are combined using a linear combination and a decision map to create the MR decomposition of the fused image. iv)The final step is to take inverse DT-CWT of the output of the previous step to get the fused image. In this thesis, we use a decomposition of 5 levels in DT-CWT. 3.1.1
Match Measure:
Match measure is an important property when you fuse images from two different sensors,because the image from one sensor will provide different information at the same pixel with respect to the image taken from another sensor. In order to properly compare two corresponding pixels, neighbourhoods surrounding the pixels should also be considered.Match measure is defined as a normalized correlation averaged over a neighbourhood of samples as in: P 2 x+m,y+nw IjA (x + m, y + n)IjB (x + m, y + n) AB mj (x, y) = P A B 2 2 x+m,y+nw ((Ij (x + m, y + n) ) + (Ij (x + m, y + n) ))
(3.1)
where IjA and IjB are successively sub sampled images at level j and w is the 37
5*5 neighbourhood.The sub sampled images should have the same matrix size as the wavelet detail coefficients. The value of mAB is an estimation of the similarity between the two images at each pixel.eg: mAB =1 indicates identical patterns,mAB lesser than 1 indicates less similarity between features, mAB =0 indicates that the gray scale value of all pixels in the neighbourhood are zero. 3.1.2
Activity Measure
Activity measure computes saliency of each pixel in the transform domain.The meaning of saliency depends on the properties of the images and the application of the image fusion process that we are doing.In our fusion process,we use the magnitude of the detail coefficients to calculate the activity measure as in: EjA
=
6 X
|Akj |
(3.2)
|Bjk |
(3.3)
k=1
EjB
=
6 X k=1
where Akj and Bjk are the complex detail coefficients of two source images at level j and orientation k. EjA and EjB are activity measures at level j for both source images, which are summations of magnitudes of detail coefficients of all six sub bands at each decomposition level. 3.1.3
Combination:
We do a linear combination of the MR decompositions of the two images to get a composite MR decomposition.Inverse transform of this gives us the fused image. The coefficients that are used in the linear combination are obtained from the activity and match measure. Fjk = CjA Akj + CjB Bjk
(3.4)
where Fjk is the k th detailed coefficient of the fused image at the j th level. CjA and CjB are the coefficients at level j that are dependent on the match measure and activity measure but independent of the orientation k. The way they depend on the match measure and activity measure is decided by the decision map.
38
The approximate coefficient is also calculated using the above equation with constants taken as 1/2 i.e a simple averaging is done between the MR decomposition coefficients of the source images A and B. 3.1.4
Decision Map:
Decision Map gives the value of the coefficients used for combining. 1st Case: mAB > .9 j CjA
EjA = A Ej + EjB
(3.5)
CjB
EjB = A Ej + EjB
(3.6)
2nd Case: .7 < mAB < .9 j
where T =
mAB j −.7 .9−.7
CjA =
EjA T + (1 − T )WA EjA + EjB
(3.7)
CjB =
EjB T + (1 − T )WB EjA + EjB
(3.8)
and
WA = 1 if EjA ≥ EjB otherwise 0 WB = 1 if EjA < EjB otherwise 0 3rd Case: mAB ≤ .7 j CjA = 1
(3.9)
CjB = 0
(3.10)
where EjA ≥ EjB
39
CjA = 0
(3.11)
CjB = 1
(3.12)
where EjA < EjB In the first case, the match measure is high which represents highly similar patterns and the fused coefficient is the linear combination of the coefficients of both input images with the weight factors determined by the relative relation between the activity measures [24]. In case 3, the match measure is very low which represents very dissimilar patterns and thus only the large activity measure one is included in the fused image. In case two, if the match measure is close to 0.9, the coefficients are close to the one discussed in case 1. If coefficients are close to 0.7, coefficients are chosen as close to s the one in case 3. Finally the fused image is obtained by performing the Inverse DT-CWT on the MR coefficients that we have obtained by using the equations above. 3.2
DUAL TREE COMPLEX WAVELET TRANSFORM(DTCWT)
In this section let us discuss about the DTCWT. Real wavelet transforms are shift variant i.e. a slight shift in the signal will alter the magnitude of the wavelet coefficients unpredictably. In Fourier transform a shift in the signal results only in the phase component of the transform changing but the magnitude component remains intact i.e it is a shift invariant transform. The basis function used for computing the fourier transform of a signal is a complex exponential whose real part(cos) and imaginary part(sin) form a Hilbert transform pair. Extending this idea, a complex wavelet transform has been designed so that the real part of the wavelet basis function forms a hilbert transform pair with its complex part. Complex wavelets can be realized using two trees of real wavelet transform filter bank structure and hence the name Dual Tree Complex Wavelet Transform. In this one tree caters to the real part of the transform and another caters to the complex part of the transform [7]. The disadvantage of using a single tree complex wavelet transform would be using filters having complex filter coefficients and no down samplers in the structure. The dual tree complex wavelet transform gives us the advantage that filter coefficients are real and the two trees obey all the properties of a real wavelet transform. The 40
FIG. 3.2: DT-CWT analysis filter bank [7]
only constraint will be in choosing the filter coefficients in such a way that the two real transform trees together must be equivalent to a complex wavelet transform, especially satisfying the shift invariance property. Detailed explanations and reasons for arriving at conditions to satisfy the above criteria is given by Nick Kingsbury [7]. Generally ten tap filters are used for implementing DT CWT. The filter coefficients can be obtained from Brooklyn’s DT-CWT tutorial web page [81]. The DT-CWT analysis and synthesis filter banks are shown in Fig. 3.2 and Fig. 3.3.
41
FIG. 3.3: DT-CWT synthesis filter bank [7]
42
CHAPTER 4 A NOVEL AUTOMATED APPROACH FOR MSRCR In MSRCR after performing color restoration it was found that the restored images were still “greyed-out”. It was found that the “greying -out” of the image happens at the single scale retinex stage itself. So any treatment for this problem had to be applied at that stage. Zia et al’s paper on MSRCR [1] does not deal with this problem. But their earlier work on SSR [2] deals with this problem and recommends a histogram based approach. In Land’s paper [82], this problem was not dealt with at all. Subsequently, Hulbert dealt with this problem in his Ph.D dissertation [83]. He proposed an automatic gain/offset approach, where the retinex values are adjusted by the absolute maximum and minimum found across all the values in all the color bands. But since this method could not give us a very good color rendition, we have used Zia et al.’s [2] canonical gain/offset method. The method shown in Fig. 4.4 actually results in clipping of some of the highest and lowest signal values but it is not a problem since those regions carry very less information. Problem with this method is that the authors have not clearly specified how to choose the lower clipping and upper clipping point. 4.1
PROPOSED METHOD
Though many modified versions of MSRCR are available in the literature [84], [85][86], none of these propose an automated approach for Retinex. We propose an automated (image independent) method to choose the upper and lower clipping points. We can choose the upper and lower clipping points using two methods, i. by using variance of the histogram as a control measure or ii. by using the frequency of occurrence of pixels as shown in the histogram as a control measure. Our initial approach was to use variance as a control measure. A particular test image was taken. After performing single scale retinex on that image, the histogram of the enhanced (single scale retinexed) image was plotted and variance was found from that. The clipping point was chosen as ‘x’ times the variance where ‘x’ can take 43
FIG. 4.4: Histogram of a single scale retinex enhanced image [1]
any value from 1 to 5 as shown in Fig. 4.5. The output image was obtained after clipping the histogram and rescaling the clipped region to 0 to 255 as shown in Fig. 4.4. But after testing this method across various images, we came to a conclusion that a unique ‘x’ value would not work for all images. So the procedure of finding clipping points cannot be automated if variance is chosen as a control measure. In the second method, control measure depends on the frequency of occurrence of pixels. For most of the images, the histogram of the enhanced image is similar to a gaussian (but not exact). The frequency of occurrence of pixel value ‘0’ in the enhanced image is found. Let this value be ‘max’ as shown in the Fig. 4.6. Clipping point or clipping pixel is chosen as the pixel that has ‘y’ times ‘max’ frequency of occurrence in the image a shown in the Fig. 4.6. After testing across many images, y = 0.02 was found to be an optimum value that can be used for many types of images. This approach has removed the image dependency that the previous method has and this is a really great advantage in real time applications [87], [88] where the user would not have time to choose the optimum clipping points for a particular image. 44
FIG. 4.5: Histogram explaining choosing of clipping points using the variance
FIG. 4.6: Histogram explaining choosing of clipping points using frequency of occurrence of pixels
45
The reason we are not able to automate our method using variance based control approach is because the histogram is not perfectly gaussian in nature. In this type of histograms, pixels which are at equidistant values from either side of the mean point but do not have the same percentage of occurrence of pixels. The values ‘a1’ and ‘a2’ are not equal in the Fig. 4.5. We mentioned that this histogram based approach is applied at SSR stage i.e. after SSR output is obtained for each surround constant values c = 15, 80, 250. But while testing we found that better outputs are obtained if we apply this technique after MSR instead of applying at SSR stage itself. And also the algorithm is computationally faster if we apply this histogram based correction once after MSR, instead of applying thrice in SSR stage.
46
CHAPTER 5 A NOVEL RETINEX BASED APPROACH FOR DEHAZING Images photographed in hazy or foggy conditions capture the scene correctly as the observer sees it. But there are applications like the landing system of aeroplanes or automatic car driving systems where it is important to see the objects that are obscured in the haze or fog. Haze removal methods can be divided into three broad categories viz. algorithms that use i) multiple images, ii) apriori information and iii) single image. We have made a detailed study literature on Haze removal in Chapter I. In 2011 Kaiming et al. [3] proposed a dark channel prior based approach. This is a single image based haze removal method based on the assumption that for a haze free image, the minimum of the three colour bands R,G and B will be zero. This method holds good for a wide range of hazy images. But the problem with this method is it is computationally complex. In this paper we propose a retinex based approach for haze removal. This method uses the Koschmieder haze formation method and also reduces the computational time when compared to the dark channel prior method. In Section 5.1 we will see the Koschmieder haze formation model. In Section 5.2 we shall see Kaiming et al.’s dark channel prior method. In Section 5.3 we shall discuss the proposed retinex based haze removal approach. 5.1
KOSCHMIEDER MODEL
The optical model used to describe the formation of haze proposed by Koschmieder is given as I c (x) = J c (x)t(x) + (1 − t(x))Ac
(5.13)
where I c (x) represents a hazy image. c represents Red, Green and Blue colour bands. J c (x) represents a haze free image in the three colour bands. t(x) represents the transmission map. We assume the transmission remains constant across all three color bands. Ac represents the global atmospheric light in three colour bands. Ac is a constant that affects all pixels in the image in the same manner. The challenge has been to obtain t(x) and Ac from the hazy image I c (x) and hence find the haze free image J c (x). The term J c (x)t(x) is called direct attenuation and it is a multiplicative distortion of the scene radiance J c (x). The term (1−t(x))Ac is called airlight and it is an additive 47
distortion of the scene radiance. In a homogeneous atmosphere, transmission t(x) can be expressed as t(x) = exp (−βd(x))
(5.14)
where d(x) represents the depth of the scene and β represents the scattering coefficient of the atmosphere. 5.2
DARK CHANNEL PRIOR METHOD
Dark channel prior method [3] is based on the principle that in a haze free image the intensity of at least one of the three channels is very low(dark) and close to zero. The dark channel can be described as J dark (x) = min ( min J c (y)) yΩ(x) cr,g,b
(5.15)
where J is any arbitrary image, x and y represents the pixels of the image. J c is a color band of J and Ω(x) is a local patch area centered at pixel x of the image. A filter of size Ω(x)=15 is typically used. We obtain the minimum of each pixel from the 3 color bands and create a ‘minimum’ image and then obtain the dark channel image by applying a minimum filter on this ‘minimum’ image. For a haze free image , J dark tends to 0 at all pixel locations since in a haze free image it is believed that atleast one of the color channels will be dark. Kaimming et al. [3] have used haze free images from the Flickr website and have come out with three reasons for the dark channel to exist i) shadows i.e shadows of cars, buildings ii) colourful objects i.e a green plant or a blue sky will have low intensity values in the other channels iii) dark objects. Since any haze free image will have a mixture of these objects, dark channel of a haze free image is expected to be of low intensity. 5.2.1
Dark Channel Prior Equations
Lets for now assume the atmospheric light Ac is given. We shall discuss in the Section about the estimation of atmospheric light. Equation 5.13 is normalised I c (x) J c (x) = t(x) + (1 − t(x)) Ac Ac Let us apply the minimum operator on both sides of the equation
48
(5.16)
FIG. 5.7: a)Input b)Transmission map t˜ c)Haze free image obtained using t˜
c I c (x) ˜ min(min J (x) ) + (1 − t(x)) ˜ ) = t(x) (5.17) c c yΩx yΩx Ac Ac J is a haze free image. So the dark channel of J is close to zero due to dark
min(min
channel prior. min(min yΩx
c
J c (x) )=0 Ac
(5.18)
By using equation 5.18 and 5.17 c ˜ = 1 − min(min I (x) ) t(x) c yΩx Ac
minyΩx (minc
I c (x) ) Ac
(5.19)
is the dark channel of the hazy image which gives us an estimate
of transmission map t. If t˜(x) is used for estimating the haze free image using equation 5.13 directly then we get a haze free image with blocking artifacts and halos as shown in Fig. 5.7 c) 5.2.2
Transmission Map Refinement
To remove these artifacts the transmission map has to be refined. So Kaiming et al. [3] used an image matting algorithm to overcome this problem. In this report we explain the overall refinement algorithm without going much deeper into details. For a detailed understanding of this algorithm one can refer to Levin et al. ’s [89] paper. Interactive image matting is a problem where the foreground and background objects are separated using some limited user input. User input is a coarse image mask called trimap where white pixels belong to foreground, dark pixels belong to background 49
FIG. 5.8: a)Input b)Transmission map t c)Haze free image obtained using t
and grey pixels belong to either. The image matting algorithm refines the trimap i.e resolves where the grey pixels belong to . The Image Matting equation can be expressed as I(x) = F(x)α + B(x)(1 − α)
(5.20)
where x is a pixel location, I is the input image F is the foreground image and B is the background image and α is the foreground opacity. From equation 5.13 and equation 5.20 we infer that transmission map t is same as α. The estimated coarse transmission map t˜ derived in equation 5.19 can be treated a user input trimap in interactive image matting. Levin et al., [89] proposed a closed form solution to natural image matting. The cost function in equation 5.21 can be minimized to obtain the transmission map t(x) ˜ can be written as vectors t and t˜ and used in equation 5.21. t(x) and t(x) E(t) = tT Lt + λ(t − ˜t)T (t − ˜ t)
(5.21)
The matrix L is called Laplacian Matrix. For an image of size M*N, the size of Laplacian Matrix is MN*MN. The (i, j) element of the matrix be defined as L(i, j) =
X k/i,jwk
(δij −
X 1 (1 + (Ii − µk )T ( + U3 )−1 (Ij − µk ))) |wk | |wk | k
(5.22)
Minimising equation 5.21 we get (L + λU )t = λ˜ t 50
(5.23)
where U is an identity matrix of same size as L and λ is a small value(10−4 ). From equation 5.23 we can find the refined transmission map t and hence the haze free image can be obtained using equation 5.13. Fig. 5.8 shows the refined transmission map and the haze free image obtained. From fig. 5.8 and fig. 5.7 it is clear that the refined transmission map does not have any blocking artifacts and hence the haze free image obtained also is free of any artifacts. This method has two disadvantages. i. Equation 5.22 is computationally very complex ii. The Laplacian matrix of size MN*MN is very big. So inversion of such a matrix as required by Equation 5.23 is not only computationally complex but consumes a lot of memory too. We infact faced ”Out of Memory” problems while implementing the above haze removal algorithm for an image of size 600*400 dimension using MATLAB. So we have come up with a Retinex based haze removal method that not only gives good results but is also computationally very fast and faces no memory issues. Section 5.3 deals about how basic Retinex model can be applied for haze removal. 5.3
RETINEX MODEL FOR HAZE REMOVAL
Fig. 5.7 b) shows that the transmission map is more or less a constant image. From Retinex theory it is known that the illumination of an image does not vary much across the image [1]. So the idea is to obtain an estimate of the transmission map using the Retinex algorithm. Transmission map is assumed to be constant across all three color bands. So we first convert the Image into HSV space and then use the luminance component alone to find the transmission map. V is the luminance component of the image. The transmission map can be obtained by using the below equation ˜ = F (x) ∗ V (x) t(x)
(5.24)
where F is the gaussian surround filter given by equation 5.25 and ∗ represents convolution. There is no need for any refinement on this transmission map. Using this transmission map, we can directly obtain the haze free image from equation 5.13.
51
F (x) = ke−r
2 /c2
(5.25)
But there is a problem in this method. The surround function F is given by equation 5.25. But the value c cannot be chosen arbitrarily. For some hazy images, a value of c = 15 gives a good haze free image. But for some other hazy images c = 80 or c = 250 gives good haze free images. With the limited hazy images database we have, we have not been able to find a reason as to why for some hazy images c = 15 works and for some images c = 80 or c = 250 works. 5.3.1
Estimation of Atmospheric light
In equation 5.13 atmospheric light Ac is estimated as the brightest pixel in the hazy image [55]. Though this atmospheric light estimation gives reasonably good results, Kaiming et al. [3] have proposed a better method of atmospheric light estimation. Their method uses the transmission map itself for estimation. Since the method has been well discussed in the original paper we are not detailing it here. We can use the transmission map that we have estimated using equation 5.24 and hence estimate the atmospheric light as proposed by Kaiming et al.
5.3.2
A minor refinement on transmission map
When we use transmission map t obtained using equation 5.24, the haze free image obtained looks dark as shown in Fig. 5.9 . So to avoid this dark image, we do a minor modification on the transmission map t obtained using equation 5.24. This new transmission map obtained using equation 5.26 is used for obtaining the haze free image. ˜ t(x) = 1 − t(x)
(5.26)
The transmission map and hence the haze free image obtained using equation 5.26 are shown in Fig. 5.10. But it is important to note that for estimating atmospheric light we use the transmission map obtained using equation 5.24 and finally for obtaining haze free image we use transmission map obtained by using equation 5.26 .
52
FIG. 5.9: a)Input b)Transmission map t c)Haze free image obtained using t
FIG. 5.10: a)Input b)Transmission map t c)Haze free image obtained using t
53
CHAPTER 6 RETINEX BASED TONE MAPPING The advancement in image acquisition technology has lead to HDR imaging which helps us to capture the entire dynamic range of a scene. The color rendition of HDR imaging is also very good. But HDR display devices have not developed in the same way and still remain costlier. So the challenge is to display these HDR images on low dynamic range display devices without compromising on the overall quality of the image. This operation is called tone mapping. K. Kim et al. proposed in their paper a new Retinex based tone mapping technique. According to the author when Multi Scale Retinex is used for tone mapping, bright regions are over highlighted at the expense of darker regions. He proposes some minor modifications to the MSR algorithm and achieves better tone mapping. 6.1
PROPOSED METHOD
We achieve better results by combining the local and global tone mapping. We first read the HDR image by using the ‘hdrread’ command in Matlab and then apply a gamma operator on the image, following the equation s = rγ where s and r are the new and old pixel values respectively. We use a γ value of .5. Kim et al. [26] have used only MSR for tone mapping. But we apply our automated MSRCR algorithm on this new image to get the final tone mapped image. The results are better because i) we have combined global and local tone mapping ii)Kim has used only MSR for tone mapping whereas we have brought in color restoration too by applying MSRCR. We have used here the luma based MSRCR approach.
FIG. 6.11: Proposed Tone Mapping method 54
CHAPTER 7 RESULTS AND DISCUSSIONS 7.1
AUTOMATED MULTI SCALE RETINEX
In this section we intend to compare the outputs that we have obtained using our automated method with that of the original method. We use the standard Retinex test images that are given in the website of NASA for comparison. The enhanced outputs provided for each test image in the website have been obtained by adjusting the parameters depending on the image. This Retinex algorithm is also commercially provided by True View Imaging Company as a part of Photoflair software. The demo version of this software is available for free. This software also uses an image independent (automated) approach for Retinex. Their method gives very good results for most of the images but very bad results for certain images (especially dark or shadowed images). But our automated approach gives very good enhancement for all types of images. We present here four images for comparing the software output with our automated approach. The image in Fig. 7.12 a) is obtained from the NASA Retinex website. Fig. 7.12 shows clearly that our method is better than software output. The red colour uniform and the subjects are better visible in our enhanced image than in the software enhanced image. The image in Fig. 7.13 a) was taken by us under poor lighting conditions. For this image too, our enhanced image is much better than the software’s. The image in Fig. 7.14 a) is a popular image in image enhancement applications. For this image too, our enhanced image is much better than the software’s. The image in Fig. 7.15 a) is also taken from NASA’s website. For this image, both the software and our proposed method perform equally well. In Fig. 7.12 and Fig. 7.13 artifacts are observed in the enhanced image. This problem can be dealt by using the tone mapping approach i.e applying a gamma map on the input image and then apply the MSRCR on the gamma mapped image. 7.2
FUSION BASED MULTI SCALE RETINEX
In this section we present the results obtained using Fusion based Multi Scale Retinex. From Fig. 7.16 and 7.17 it is clear that the outputs obtained using the original Multi
55
FIG. 7.12: a) Input image b) Our output c) Software output
FIG. 7.13: a) Input image b) Our output c) Software output 56
FIG. 7.14: a) Input image b) Our output c) Software output
FIG. 7.15: a) Input image b) Our output c) Software output 57
Scale Retinex and Fusion based Multi Scale Retinex are perceptually same. 7.3
RETINEX BASED HAZE REMOVAL METHOD
Fig. 7.18, Fig. 7.20, Fig. 7.22 and Fig. 7.24 gives the input image and the associated transmission maps and haze free images obtained before and after matting(Kaiming et al.’s dark channel prior method). We observe that the haze free image obtained before matting has blocking artifacts and the haze free image obtained after matting is free of these artifacts. In Fig. 7.19, Fig. 7.21, Fig. 7.23 and Fig. 7.25 the transmission maps and associated haze free images obtained for surround constant values of c = 15, 80 and 250 are shown. i) ‘Mountain’ Image: From Fig. 7.18 and Fig. 7.19 we can infer that in retinex based method a good haze free image is obtained for c=250 but the surprising fact is that the transmission map obtained for c=250 and the map obtained using dark channel prior do not appear same. In this example for c= 15 or 80 we observe that the colour of the soil is over saturated. ii) ‘Tokyo’ Image: From Fig. 7.20 and Fig. 7.21 we can infer that in retinex based method a good haze free image is obtained for c=80 or 250. For c=15, the colour of the building looks unnatural. iii) ‘Hongkong’ Image: From Fig. 7.22 and Fig. 7.23 we can infer that in retinex based method it is difficult to judge which is the best output. In Fig. 7.23 for c=15, the ‘black’ color of the buildings(in the background) in the haze free image looks same as in the haze free image obtained using dark channel prior method but the colour of the foreground buildings looks a bit saturated. For c=250, color of the foreground building is well preserved but the ’black’ colour of the buildings looks saturated. So a reasonable compromise is the haze free image obtained for c=80. iv) ‘House’ Image: From Fig. 7.24 and Fig. 7.25 we can infer that in retinex based method a good haze free image is obtained for c=15. From the above observations it is clear that for different images under consideration, different ’c’ values give good haze free images. So the future challenge is to automate the algorithm by estimating that ’c’ value which will give good haze free image from some image statistics.
58
FIG. 7.16: a) Input image b) Output due to addition c) Output due to fusion
59
FIG. 7.17: a) Input image b)Output due to addition c) Output due to fusion
60
FIG. 7.18: Results of ‘Mountain’ Image obtained using dark channel prior method
61
FIG. 7.19: Results of ‘Mountain’ Image obtained using Retinex based haze removal method
62
FIG. 7.20: Results of ‘Tokyo’ Image obtained using dark channel prior method
63
FIG. 7.21: Results of ‘Tokyo’ Image obtained using Retinex based haze removal method
64
FIG. 7.22: Results of ‘Hongkong’ Image obtained using dark channel prior method
65
FIG. 7.23: Results of ‘Hongkong’ Image obtained using Retinex based haze removal method
66
FIG. 7.24: Results of ‘House’ Image obtained using dark channel prior method
FIG. 7.25: Results of ‘House’ Image obtained using Retinex based haze removal method
67
7.3.1
Computational complexity
In this section we shall compare the computational complexity of Kaiming’s dark channel prior method and our RETINEX based approach. We have used MATLAB simulation for the comparison. After testing the dark channel prior method on images of different sizes we have observed the simulation time and listed it in table 7.1.
TABLE 7.1: Computational complexity of dark channel prior method Size of the image Time taken(seconds) 64*64 10.18 128*128 151.36 256*256 2379.75 From the table 7.1 it is clear that Kaiming et al. ’s method takes a lot of time even for a small 256*256 image. And from this table we can also infer that if a N*N image takes t seconds for simulation then a 2N*2N image takes nearly 15 ∗ t seconds for simulation. It is also not possible to simulate using this method even in our modern computers an image of size 400 ∗ 400 because of ‘Memory’ issues. The reason is for an image of size N ∗ N , we need a Laplacian Matrix L of size N 2 ∗ N 2 . In contrast our Retinex based method can be used to remove haze of image of any size and also the computational time is very less as shown in the table 7.2.
TABLE 7.2: Computational complexity of the proposed method Size of the image Time taken(seconds) 128*128 2.09 256*256 2.466 512*512 6.848 1024*1024 32.737
7.4
RETINEX BASED TONE MAPPING
In Fig. 7.26 and 7.27 the first image corresponds to the original HDR image. The second image b) is the output image obtained after applying Retinex on the image a). The third image c) is obtained after applying the gamma operator on the input 68
HDR image. The fourth image d) is obtained after applying Retinex on image c). The enhanced output obtained by apply Retinex without gamma mapping or gamma mapping without Retinex does not give good outputs. But when gamma mapping and Retinex are combined appropriately, tone mapping result is good.
69
FIG. 7.26: a)Input HDR image b)Retinex output c)Gamma mapped HDR image d)Final Retinex output
70
FIG. 7.27: a)Input HDR image b)Retinex output c)Gamma mapped HDR image d)Final Retinex output
71
CHAPTER 8 CONCLUSION AND FUTURE WORK In this chapter we shall briefly describe about our contributions and give a possible future direction for research in this field. 8.1
CONCLUSION
We have developed a fusion based approach for the Multi Scale Retinex with Color Restoration(MSRCR) algorithm. Multi Scale Retinex(MSR) is obtained by averaging the Single Scale Retinex(SSR) outputs. We propose to fuse these SSR outputs. We have used a wavelet and pixel based approach for fusion. First we obtain the wavelet coefficients of the images for 5 levels of decomposition. The fusion of the coefficients is done by using a decision rule and a combination map. To implement this two measures are used i.e i) match measure and ii) activity measure. Match measure computes the correlation between the images in the spatial domain. Activity measure gives the energy of the image at each level of decomposition. Once the fused wavelet coefficients are obtained, then inverse wavelet transform is applied to get the fused image. We used Dual Tree Complex Wavelet Transform. We have proposed an automated approach for MSRCR. In MSRCR the enhanced outputs have a “washed out” look. To avoid this problem, we use the histogram of the enhanced SSR image. For most of the images the histogram has a gaussian like shape. The approach was to clip some intensities on either side of the histogram and scale the remaining pixels to 0-255. This approach gives very good enhanced images. In this thesis we have developed an approach to determine the clipping points. We use the frequency of occurrence of pixels as a parameter from the histogram to obtain the clipping points. After testing across a good number of images we have found that it is enough to clip on either side of the histogram intensities that account for 2 percent of the total number of intensity values and scale the rest to 0-255. NASA’ s MSRCR algorithm has been implemented as a software package called ‘Photo Flair’ by ‘TrueView’ company, U. S. A. The demo version of the software can be downloaded for free from the company’s website and used. The software also seems to use an automated approach since they do not ask users to manually input any of the parameters. The software works well for most of the images except for 72
dark ones i.e for scenes whose dynamic range is high and the camera is unable to capture all the details. Our automated approach works well on all the images we have tested including the dark ones. We have proposed a Retinex based haze removal algorithm. Dark channel prior based haze removal is the most popular dehazing algorithm. This method is based on the assumption that every haze free image has a dark channel. This method gives very good results but is computationally complex too. In our method we obtain the dark channel of the hazy image by using the Retinex algorithm. The dark channel appears as a constant image. From Retinex algorithm we know that the illuminance of an image appears as a constant i.e illuminance is slow varying. So we obtain the dark channel of the hazy image from the illuminance component of the image. This approach gives same results and is also computationally simpler than the dark channel prior method. We have proposed a novel tone mapping approach. In our approach we combine the global(Gamma) and local(Retinex) tone mapping operator. We first gamma map the image and then apply MSRCR. Our approach is simple and gives comparable results to the one obtained by Kim et al. [26]. 8.2
FUTURE WORK
In the fusion based approach for the MSRCR algorithm, the enhanced outputs obtained for various images are perceptually same as the original algorithm. One of the future areas of work in this field will be to try other fusion algorithms or wavelets instead of the ones we have used and apply them on images in which original MSRCR’s performance is bad. In the automated approach for the MSRCR algorithm we have found the optimal value for automation by testing across various images. This method holds good across a good number of images we tested. But to make sure the algorithm works across all images we have to arrive at this optimal value mathematically. Retinex based haze removal algorithm gives good results and is also computationally simpler. The challenge we are working on now is automation i.e. by using some image statistics our algorithm must be able to automatically choose a surround constant ’c’ value that will remove the haze in an image and simultaneously give good colour rendition. In the Retinex based tone mapping, the first step is to apply a gamma mapping 73
on the HDR image and then apply MSRCR. We have chosen a constant gamma value of 0.5 for all images. This can also be automated i.e. an algorithm can be developed to find an appropriate gamma value for each image using its own statistics.
74
BIBLIOGRAPHY [1] D. J. Jobson, Z. Rahman, and G. A. Woodell, “Properties and performance of a centre/surround retinex,” IEEE Transactions on Image Processing, vol. 6, no. 3, 1997. [2] D. J. Jobson, Z. Rahman, and G. A. Woodell, “A multiscale retinex for bridging the gap between color images and the human observation of scenes,” IEEE Transactions on Image Processing, vol. 6, no. 7, pp. 965–976, 1997. [3] K. He, J. Sun, and X. Tang, “Single image haze removal using dark channel prior,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 12, pp. 2341–2353, 2011. [4] R. E. Woods R. C. Gonzalez, Digital Image Processing, Third Edition, Pearson Publications, 2005. [5] S. Greenberg, M. Aladjem, D. Kogan, and I. Dimitrov, “Fingerprint image enhancement using filtering techniques,” in Proceedings of 15th IEEE International Conference on Pattern Recognition, vol. 3, pp. 322–325, 2000. [6] Z. Zhang and L. W. He, “Whiteboard scanning and image enhancement,” Digital Signal Processing, vol. 17, no. 2, pp. 414–432, 2007. [7] I. W. Selesnick, R. G. Baraniuk, and N. C. Kingsbury, “The dual-tree complex wavelet transform,” IEEE Signal Processing Magazine, vol. 22, no. 6, pp. 123– 151, 2005. [8] S. M. Pizer, J. B. Zibberman, and E. Stab, “Adaptive gray level assignment in ct scan display,” Journal of Computer Assistant Tomography. [9] S.M. Pizer, E.P. Amburn, J.D. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. ter Haar Romeny, J.B. Zimmerman, and K. Zuiderveld, “Adaptive histogram equalization and its variations,” Computer Vision, Graphics, and Image Processing, vol. 39, no. 3, pp. 355–368, 1987. [10] Z. Rahman, GA Woodell, and DJ Jobson, “Retinex image enhancement: Application to medical images,” in Proceedings of the NASA Workshop on New Partnerships in Medical Diagnostic Imaging, 2001. 75
[11] A.M. Trifas and J.M. Director-Tyler, Medical image enhancement, Louisiana State University and Agricultural & Mechanical College, 2005. [12] S. Agaian and S.A. McClendon,
“Novel medical image enhancement algo-
rithms,” in Proceedings of SPIE, vol. 7532, pp. 75320W–75320W–12, 2010. [13] I.N. Bankman, Handbook of medical imaging: processing and analysis, Academic Press, 2000. [14] Z. Shi and V. Govindaraju, “Historical document image enhancement using background light intensity normalization,” in Proceedings of the 17th IEEE International Conference on Pattern Recognition, vol. 1, pp. 473–476, 2004. [15] S. C. Pei, Y. C. Zeng, and C. H. Chang, “Virtual restoration of ancient chinese paintings using color contrast enhancement and lacuna texture synthesis,” IEEE Transactions on Image Processing, vol. 13, no. 3, pp. 416–429, 2004. [16] J. Tang, J. Kim, and E. Peli, “Image enhancement in the jpeg domain for people with vision impairment,” IEEE Transactions on Biomedical Engineering, vol. 51, no. 11, pp. 2013–2023, 2004. [17] E. Peli, R. B. Goldstein, G. M. Young, C. L. Trempe, and S. M. Buzney, “Image enhancement for the visually impaired. simulations and experimental results,” Journal on Investigative Ophthalmology and Visual Science, vol. 32, no. 8, pp. 2337–2350, 1991. [18] E. Peli, E. Lee, C. L. Trempe, and S. Buzney, “Image enhancement for the visually impaired: the effects of enhancement on face recognition,” Journal of Optical Society of America, vol. 11, no. 7, pp. 1929–1939, 1994. [19] T. Ko, “Fingerprint enhancement by spectral analysis techniques,” in 31st Proceedings of Applied Imagery Pattern Recognition Workshop, pp. 133–139, 2002. [20] G. Woodell, D.J. Jobson, Z. Rahman, and G. Hines, “Advanced image processing of aerial imagery,” in Proceedings of SPIE, vol. 6246, pp. 62460E, 2006.
76
[21] H. Y. Yang, P. Y. Chen, C. C. Huang, Y. Z. Zhuang, and Y. H. Shiau, “Low complexity underwater image enhancement based on dark channel prior,” in Proceedings of the Second International Conference on Innovations in Bio-inspired Computing and Applications (IBICA), pp. 17–20, 2011. [22] J. V. Candy and D. H. Chambers, “Model-based enhancement of internal wave images,” IEEE Journal of Oceanic Engineering, vol. 22, no. 1, pp. 1–8, 1997. [23] W. J. Heng and Q. Tan, “Content enhancement for e-learning lecture video using foreground/background separation,” in Proceedings of the IEEE Workshop on Multimedia Signal Processing, pp. 436–439, 2002. [24] L. Tao, Multi-modal enhancement techniques for visibility improvement of digital images, Ph.D. thesis, Old Dominion University, 2006. [25] Erik Maltin, “Single image haze and noise removal,” M.S. thesis, University of California Santa Cruz, 2011. [26] K. Kim, J. Bae, and J. Kim, “Natural hdr image tone mapping based on retinex,” IEEE Transactions on Consumer Electronics, vol. 57, no. 4, pp. 1807– 1814, 2011. [27] Land. E. H, “Lightness and retinex theory,” Journal of Optical Society of America, 1964. [28] Laurence Meylan, Tone mapping for high dynamic range images, Ph.D. thesis, EPFL, Switzerland, 2006. [29] A. Rizzi, D. Marini, L. Rovati, and F. Docchio, “Unsupervised corrections of unknown chromatic dominants using a brownian-path-based retinex algorithm,” Journal of Electronic Imaging, vol. 12, pp. 431, 2003. [30] D. Marini and A. Rizzi, “Color constancy and optical illusions: a computer simulation with retinex theory,” in Proceedings of the Seventh International Conference on Image Analysis and Processing, World Scientific, New Jersey, pp. 657–660, 1993. [31] A. Moore, J. Allman, and R.M. Goodman, “A real-time neural system for color constancy,” IEEE Transactions on Neural Networks, vol. 2, no. 2, pp. 237–247, 1991. 77
[32] A. C. Hurlbert, The Computation of Color, Ph.D. thesis, MIT, Cambridge, September 1989. [33] L. Meylan and S. S¨ usstrunk, “The influence of luminance on local tone mapping,” in Proceedings of the 13th Color Imaging Conference, pp. 276–282, 2005. [34] C. Pohl and J. L Van Genderen, “Review article multisensor image fusion in remote sensing: concepts, methods and applications,” International Journal of remote sensing, vol. 19, no. 5, pp. 823–854, 1998. [35] P.J. Burt and R.J. Kolczynski, “Enhanced image capture through fusion,” in Proceedings of the Fourth IEEE International Conference on Computer Vision, pp. 173–182, 1993. [36] M.M. Daniel and A.S. Willsky, “A multiresolution methodology for signal-level fusion and data assimilation with applications to remote sensing,” in Proceedings of the IEEE, vol. 85, no. 1, pp. 164–180, 1997. [37] B.V. Dasarathy, “Fuzzy evidential reasoning approach to target identity and state fusion in multisensor environments,” Journal of Optical engineering, vol. 36, pp. 669, 1997. [38] M. P. Dubuisson and A. K. Jain, “Contour extraction of moving objects in complex outdoor scenes,” International Journal of Computer Vision, vol. 14, no. 1, pp. 83–105, 1995. [39] M.E. Ulug and C.L. McCullough, “Feature and data-level fusion of infrared and visual images,” in Proceedings of SPIE, vol. 3719, pp. 312, 1999. [40] K. Yamamoto and K. Yamada, “Image processing and fusion to detect navigation obstacles,” in Proceedings of SPIE, vol. 3374, pp. 337, 1998. [41] P. Burt and E. Adelson, “The laplacian pyramid as a compact image code,” IEEE Transactions on Communications, vol. 31, no. 4, pp. 532–540, 1983. [42] P. J. Burt, The pyramid as a structure for efficient computation, Rensselaer Polytechnic Institute, 1983. [43] A. Toet, “Image fusion by a ration of low-pass pyramid,” Pattern Recognition Letters, vol. 9, no. 4, pp. 245–253, 1989. 78
[44] P. J. Burt and R. J. Kolczynski, “Enhanced image capture through fusion,” in Proceedings of Fourth International Conference on Computer Vision, pp. 173– 182, 1993. [45] T. A. Wilson, S. K. Rogers, and M. Kabrisky, “Perceptual-based image fusion for hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 35, no. 4, pp. 1007–1017, 1997. [46] H. Li, B. S Manjunath, and S .K. Mitra, “Multisensor image fusion using the wavelet transform,” Journal of Graphical models and image processing, vol. 57, no. 3, pp. 235–245, 1995. [47] O. Rockinger, “Image sequence fusion using a shift-invariant wavelet transform,” in Proceedings of IEEE International Conference on Image Processing, vol. 3, pp. 288–291, 1997. [48] P. Hill, N. Canagarajah, and D. Bull, “Image fusion using complex wavelets,” in Proeedings of the 13th British Machine Vision Conference, pp. 487–496, 2002. [49] Y. Y. Schechner, S. G. Narasimhan, and S. K. Nayar, “Instant dehazing of images using polarization,” vol. 1, pp. 325–332, 2001. [50] S. K. Nayar and S. G. Narasimhan, “Vision in bad weather,” in Proceedings of the Seventh IEEE International Conference on Computer Vision, vol. 2, pp. 820–827, 1999. [51] S.G. Narasimhan and S.K. Nayar, “Interactive (de) weathering of an image using physical models,” in Proceedings of the IEEE Workshop on Color and Photometric Methods in Computer Vision, pp. 1–8, 2003. [52] J. Kopf, B. Neubert, B. Chen, M. Cohen, D. Cohen-Or, O. Deussen, M. Uyttendaele, and D. Lischinski, “Deep photo: Model-based photograph enhancement and viewing,” ACM Transactions on Graphics (TOG), vol. 27, no. 5, pp. 116, 2008. [53] S. G Narasimhan and S. K Nayar, “Vision and the atmosphere,” International Journal on Computer Vision, vol. 48, no. 3, pp. 233–254, 2002. [54] R. T. Tan, “Visibility in bad weather from a single image,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–8, 2008. 79
[55] J. P. Tarel and N. Hautiere, “Fast visibility restoration from a single color or gray level image,” in Proceedings of the IEEE 12th International Conference on Computer Vision, pp. 2201–2208, 2009. [56] Jonson. G. M., “Spatial appearance modelling and hdr rendering,” in Proceedings of IST/SPIE Electronic Imaging, San Jose, pp. 148–156, 2005. [57] S. N. Pattanaik, J. Tumblin, H. Yee, and D. P. Greenberg, “Time-dependent visual adaptation for fast realistic image display,” in Proceedings of the 27th Annual Conference on Computer graphics and Interactive techniques, pp. 47– 54, 2000. [58] M. A. Webster, “Human colour perception and its adaptation,” Network: Computation in Neural Systems, vol. 7, no. 4, pp. 587–634, 1996. [59] R. W. G. Hunt and W. D Wright, The reproduction of colour, John Wiley & Sons, 2004. [60] J. J. McCann, “Art, science, and appearance in hdr,” Journal of the Society for Information Display, vol. 15, pp. 709, 2007. [61] S. Ochi and S. Yamanaka, “Solid state image pickup device,” Sept. 10 1985, US Patent 4,541,016. [62] L. E. Alston, D. S. Levinstone, and W. T. Plummer, “Exposure control system for an electronic imaging camera having increased dynamic range,” Mar. 3 1987, US Patent 4,647,975. [63] S. Mann and R. W Picard, “On being undigital with digital cameras extending dynamic range by combining differently exposed pictures,” in Proceedings of the 48th Annual Conference on IS and T, pp. 422–428, 1995. [64] P. E. Debevec and J. Malik, “Recovering high dynamic range radiance maps from photographs,” in Proceedings of the 24th annual conference on Computer graphics and interactive techniques, 1997. [65] J. A. Ferwerda, S. N. Pattanaik, P. Shirley, and D. P. Greenberg, “A model of visual adaptation for realistic image synthesis,” in Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp. 249– 258, 1996. 80
[66] Vladimir. S. Petrovic, Multi Sensor Pixel Level Image Fusion, Ph.D. thesis, University of Manchester, England, 2001. [67] T. Mitsunaga and S. K. Nayar, “Radiometric self calibration,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 1, 1999. [68] G. Ward, “Fast, robust image registration for compositing high dynamic range photographs from hand-held exposures,” Journal of Graphics Tools, p. 329, 2005. [69] A. Tomaszewska and R. Mantiuk, “Image registration for multi-exposure high dynamic range image acquisition,” in Proceedings of International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision (WSCG), pp. 49–56, 2007. [70] S.B. Kang, M. Uyttendaele, S. Winder, and R. Szeliski, “High dynamic range video,” ACM Transactions on Graphics, vol. 22, no. 3, pp. 319–325, 2003. [71] E.A. Khan, A. O Akyuz, and E. Reinhard, “Ghost removal in high dynamic range images,” in Proceedings of the IEEE International Conference on Image Processing, Atlanta, pp. 2005–2008, 2006. [72] M. Pedone and J. Heikkil¨a, “Constrain propagation for ghost removal in high dynamic range images,” in Proceedings of International Conference on Computer Vision Theory and Applications, 2008. [73] O. Gallo, N. Gelfandz, W. C. Chen, M. Tico, and K. Pulli, “Artifact-free high dynamic range imaging,” in Proceddings of the IEEE International Conference on Computational Photography (ICCP), pp. 1–7, 2009. [74] D. Sidibe, W. Puech, and O. Strauss, “Ghost detection and removal in high dynamic range images,” in Proceedings of the European Signal Processing Conference, Glasgow, 2009. [75] C. Schlick, “Quantization techniques for visualization of high dynamic range pictures,” in Proceedings of the Fifth European Graphics Workshop on Photorealistic Rendering Techniques, vol. 94, pp. 7–20, 1994.
81
[76] P. Ferschin, I. Tastl, and W. Purgathofer, “A comparison of techniques for the transformation of radiosity values to monitor colors,” in Proceedings of the IEEE International Conference on Image Processing, vol. 3, pp. 992–996, 1994. [77] F. Drago, K. Myszkowski, T. Annen, and N. Chiba, “Adaptive logarithmic mapping for displaying high contrast scenes,” Presented in the Eurographics Computer Graphics Forum, vol. 22, no. 3, pp. 419–426, 2003. [78] M. Ashikhmin, “A tone mapping algorithm for high contrast images,” in Proceedings of the 13th Eurographics Workshop on Rendering, pp. 145–156, 2002. [79] B. K. P. Horn, “Determining lightness from an image,” Journal on Computer Graphics and Image Processing, vol. 3, no. 4, pp. 277–299, 1974. [80] R. Fattal, D. Lischinski, and M. Werman, “Gradient domain high dynamic range compression,” ACM Transactions on Graphics, vol. 21, no. 3, pp. 249–256, 2002. [81] “http://taco.poly.edu/waveletsoftware/dt2d.html,” . [82] E. H. Land, “An alternative technique for the computation of the designator in the retinex theory of color vision,” Proceedings of the National Academy of Sciences, vol. 83, no. 10, pp. 3078, 1986. [83] A. Hurlbert, “Formal connections between lightness algorithms,” Journal of Optical Society of America, vol. 3, no. 10, pp. 1684–1693, 1986. [84] Xian. Y Yong. L, “Multi dimensional multi scale image enhancement algorithm,” in Proceedings of the Second International Conference on Signal Processing Systems(ICSPS), 2010. [85] L. Tao, R. Tompkins, and V.K. Asari, “An illuminance-reflectance model for nonlinear enhancement of color images,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition-Workshops, pp. 159–159, 2005. [86] H. Hu and G. Ni, “Color image enhancement based on the improved retinex,” in Proceedings of International Conference on Multimedia Technology (ICMT), pp. 1–4, 2010.
82
[87] W.J. Kyung, D.C. Kim, and Y.H. Ha, “Real-time multi-scale retinex to enhance night scene of vehicular camera,” in Proceedings of the 17th IEEE Korea-Japan Joint Workshop on Frontiers of Computer Vision (FCV), pp. 1–4, 2011. [88] G.D. Hines, Z. Rahman, D.J. Jobson, G.A. Woodell, and S.D. Harrah, “Realtime enhanced vision system,” in Proceedings of the Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol. 5802, pp. 127–134, 2005. [89] A. Levin, D. Lischinski, and Y. Weiss, “A closed-form solution to natural image matting,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 2, pp. 228–242, 2008.
83
LIST OF PUBLICATIONS 1. Sudharsan Parthasarathy, Praveen Sankaran, “An Automated Multi Scale Retinex with Color Restoration for Image Enhancement,” National Conference on Communications, 2012, IIT Kharagpur, India. 2.
Sudharsan Parthasarathy, Praveen Sankaran, “Fusion Based Multi Scale
Retinex with Color Restoration for Image Enhancement,” , International Conference on Computer Communication and Informatics, 2012, Coimbatore, India. 3. Sudharsan Parthasarathy, Praveen Sankaran,“A RETINEX based Haze removal method,” ,International Conference on Industrial and Information Systems, August 6-9 2012, IIT Madras(COMMUNICATED).
84
