Algorithms for processing polarization rich optical imaging data

Algorithms for processing polarization-rich optical imaging data A Thesis Submitted for the Degree of Master of Science (Engineering) in the Faculty of Engineering by Umesh R S Department of Electrical Engineering Indian Institute of Science Bangalore – 560 012 May 2004 i Acknowledgments Four magnanimous people have been behind the completion of this thesis. First and foremost, I must thank my advisor Dr A G Ramakrishnan, for having confidence in me and supporting me all the way, till the thesis saw its logical end. His simplicity, openness and tranquility have influenced me deeply. In no way less, stands Dr Hema Ramachandran of the Raman Research Institute, Bangalore, who provided me with all the necessary experimental facilities and took care of me as her own student, though I vexed her by absconding time and again from that wonderful place. Without her support, this thesis would not have seen light. I owe a lot of my learning to her and RRI. Srikanth gave a brotherly backing all throughout my thesis work. Without his gyaan of estimation theory and his battery of key references, this work would have taken an inordinately long time. Divya deserves special thanks for her unflinching support on the experimental front. I am indebted to her for teaching me optics, light scattering and for sharing the ups and downs of research work. On the non-technical front, I thank Yoga for making my life richer with his arguments and thoughts ranging from philosophy to eating. Hari and Bhatta ensured that we were served a full course of jokes during lunch and dinner times. Special thanks to them and to Vatsan and VK who helped me a lot during course work. My peers in Biomedical lab and Multimedia lab deserve their share of thanks, for making my stay in IISc a cherishable one. Equally, people of the Optics lab of RRI deserve special thanks too, for making me feel one among them. Besides these, Nityotsava was always a welcome distraction towards the end of my thesis work. It brought me close to some extremely good and like minded people, of whom, I must mention Anand. I should also thank Prof K R Ramakrishnan for extending his support in the absence of Dr AGR. I am grateful to Dr Ashok Rao, for having initiated me into the wonderful subject of signal processing. There are many others who have influenced my attitudes towards many things in life during this period, for whom I am indebted, but the list is too long to be produced here. But more than anybody else, I must thank amma, anna, akka, bhaava and paapu, who selflessly supported me all the while. This thesis is dedicated to them. Contents List of symbols and abbreviations viii Abstract x 1 Introduction to optical imaging 1 1.1 Imaging science . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Optical imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Light scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Organization of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Scattering of polarized light 2.1 2.2 7 Polarized light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 States of polarization . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 Stokes vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Light scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Single scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.2 Multiple scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Mueller matrices and light scattering . . . . . . . . . . . . . . . . . 20 3 Imaging through scattering media: A review 23 3.1 A brief overview of optical imaging . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Polarization based direct imaging schemes . . . . . . . . . . . . . . . . . . 27 3.2.1 Effect of scatterers on imaging . . . . . . . . . . . . . . . . . . . . . 28 3.2.1.1 The volume and density of scatterers . . . . . . . . . . . . 30 ii iii CONTENTS 3.2.1.2 The size of the scatterers . . . . . . . . . . . . . . . . . . 31 3.2.1.3 The relative refractive index of the scatterers . . . . . . . 32 3.2.2 Effect of source parameters and polarization optics . . . . . . . . . 32 3.2.3 Effect of spatial filtering and detector characteristics . . . . . . . . 34 3.2.4 Imaging schemes and processing algorithms . . . . . . . . . . . . . 35 3.2.4.1 Polarization difference imaging schemes . . . . . . . . . . 36 3.2.4.2 Polarization modulation imaging schemes . . . . . . . . . 41 3.2.4.3 Polarization diversity active imaging schemes . . . . . . . 43 4 Processing polarization-rich data 46 4.1 Signal modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 Polarization intensity imaging . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2.1 4.3 Polarization intensity imaging - the case of white noise . . . . . . . 52 4.2.1.1 Intensity imaging using PDI estimator . . . . . . . . . . . 52 4.2.1.2 Intensity imaging using PMI estimator . . . . . . . . . . . 54 4.2.1.3 The holy grail for intensity imaging . . . . . . . . . . . . . 55 4.2.1.4 Comparison of estimators . . . . . . . . . . . . . . . . . . 60 4.2.2 Polarization intensity imaging - the case of coloured noise . . . . . . 61 4.2.3 Chunked data processing and bootstrapping . . . . . . . . . . . . . 65 4.2.4 Computational complexity of the estimators . . . . . . . . . . . . . 66 4.2.5 Comparison of polarization intensity estimators . . . . . . . . . . . 68 4.2.5.1 Comparisons based on Monte-Carlo simulations . . . . . . 68 4.2.5.2 Comparisons based on performance in real noise . . . . . . 70 4.2.5.3 Comparison of estimators within and across classes . . . . 73 4.2.5.4 A weak case for the APES estimator . . . . . . . . . . . . 79 4.2.5.5 A case for the chunked PMI estimator . . . . . . . . . . . 80 4.2.5.6 Effect of chunking parameters on PII APES estimator . . 81 4.2.5.7 Effect of bootstrap length on APES estimator . . . . . . . 83 4.2.5.8 Other issues related to polarization intensity imaging . . . 87 Imaging using Degree of polarization . . . . . . . . . . . . . . . . . . . . . 89 4.3.0.9 The PDI DOLP estimators . . . . . . . . . . . . . . . . . 91 iv CONTENTS 4.3.0.10 The PMI DOLP estimator . . . . . . . . . . . . . . . . . . 92 4.3.1 4.4 4.5 Comparison of DOLP imaging estimators . . . . . . . . . . . . . . . 92 4.3.1.1 A case against the chunked PMI DOLP estimator . . . . . 97 4.3.1.2 Effect of chunking on APES DOLP estimator . . . . . . . 97 4.3.1.3 Effect of bootstrap length on APES DOLP estimator . . . 98 Polarization orientation imaging . . . . . . . . . . . . . . . . . . . . . . . . 102 4.4.1 Polarization orientation imaging - the case of white noise . . . . . . 103 4.4.2 Polarization orientation imaging - the case of coloured noise . . . . 103 4.4.3 Comparison of POI estimators . . . . . . . . . . . . . . . . . . . . . 104 4.4.3.1 A case for the APES POI estimator . . . . . . . . . . . . 105 4.4.3.2 Effect of chunking on the APES POI estimator . . . . . . 110 A potpourri of visualization parameters . . . . . . . . . . . . . . . . . . . . 110 5 Imaging results 116 6 Conclusions and topics for further research 129 6.1 Contributions of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.2 Extension to the case of circular polarization . . . . . . . . . . . . . . . . . 130 6.3 Topics for further research . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.3.1 Correlation based processing techniques . . . . . . . . . . . . . . . . 131 6.3.2 The imaging problem vis-a-vis the image transfer problem . . . . . 131 6.3.3 The vectorial point spread function of light . . . . . . . . . . . . . . 132 6.3.4 Post processing procedures . . . . . . . . . . . . . . . . . . . . . . . 133 6.3.5 Harmonic based processing scheme . . . . . . . . . . . . . . . . . . 133 A CRLB calculations - The case of white noise 136 B CRLB calculations - The case of coloured noise 140 Bibliography 145 List of Figures 1.1 Effect of increase of number of scatterers on imaging - Simulated images. . 4 2.1 Vibration ellipse with ellipticity B/A and azimuth γ . . . . . . . . . . . . . 9 3.1 A block diagrammatic approach to polarization based imaging techniques . 27 3.2 Summary of scattering regimes 4.1 The modified PMI imaging scheme . . . . . . . . . . . . . . . . . . . . . . 47 4.2 Performance of PII estimators in white noise. . . . . . . . . . . . . . . . . 76 4.3 Performance of PII estimators in AR1 noise. . . . . . . . . . . . . . . . . . 77 4.4 Performance of PII estimators in real noise. . . . . . . . . . . . . . . . . . 78 4.5 PII PMI N and its chunked variants in white and coloured noise . . . . . . 82 4.6 Effect of chunk length on PII APES C estimator in coloured noise . . . . . 84 4.7 Effect of overlap length on PII APES estimator in white noise . . . . . . . 85 4.8 Effect of overlap length on PII APES estimator in coloured noise . . . . . . 86 4.9 Effect of bootstrap length on PII APES B estimator . . . . . . . . . . . . . 88 . . . . . . . . . . . . . . . . . . . . . . . . 30 4.10 Performance of estimators of DOLP in white noise. . . . . . . . . . . . . . 94 4.11 Performance of estimators of DOLP in AR1 noise. . . . . . . . . . . . . . . 95 4.12 Performance of estimators of DOLP in real noise. . . . . . . . . . . . . . . 96 4.13 Effect of chunk length on DOLP APES C estimator in white noise. . . . . 99 4.14 Effect of overlap length on APES C DOLP estimator in white noise . . . . 100 4.15 Effect of bootstrap length on DOLP APES estimator in white noise. . . . . 101 4.16 Performance of phase estimators in white noise. . . . . . . . . . . . . . . . 106 4.17 Performance of phase estimators in AR1 noise. . . . . . . . . . . . . . . . . 107 4.18 Performance of phase estimators in real noise. . . . . . . . . . . . . . . . . 108 v LIST OF FIGURES vi 4.19 A case for APES N phase estimator in white noise . . . . . . . . . . . . . . 109 4.20 Effect of chunk length on APES C phase estimator . . . . . . . . . . . . . 111 4.21 Result obtained by fusing the visualization parameters . . . . . . . . . . . 115 5.1 Results to illustrate the bias of PII PDI N estimator . . . . . . . . . . . . . 122 5.2 Comparison of results of PII and DOLP imaging. . . . . . . . . . . . . . . 123 5.3 Comparison of PII results obtained using PDI N and PMI N. . . . . . . . . 124 5.4 Comparison of PII results for SET 1, obtained from various estimators . . 125 5.5 Comparison of DOPI results for SET 2, obtained from various estimators . 126 5.6 Comparison of results obtained with coherent and incoherent sources. . . . 127 5.7 Segmentation of a POI result. . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.1 Results of sinusoidal detection on data of SET 2 . . . . . . . . . . . . . . . 133 6.2 Result obtained by the harmonic based processing scheme . . . . . . . . . 135 List of Tables 2.1 Stokes representation of the fundamental states of polarization . . . . . . . 13 2.2 Anisotropy values of spherical polystyrene particles of different sizes, used in our imaging experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.1 Data sets used for analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.2 Statistics related to the data sets used for analysis. . . . . . . . . . . . . . 117 vii List of symbols and abbreviations PDI Polarization Difference Imaging PMI Polarization Modulation Imaging CCD Charge Coupled Device DOP Degree of Polarization DOLP Degree of Linear Polarization RRI Relative Refractive Index g Anisotropy parameter ls Scattering mean free path L Scattering slab thickness l∗ (TMFP) Transport Mean Free Path τ Optical thickness M Mueller matrix SV Stokes Vector CW Continuous wave SNR Signal to Noise Ratio ODLP Observed Degree of Linear Polarization PDAI Polarization Diversity Active Imaging MMIP Mueller Matrix Imaging Polarimetry PII Polarization Intensity Imaging DOPI Degree of Polarization Imaging POI Polarization Orientation Imaging AR Auto Regressive iid Independent and identically distributed viii LIST OF SYMBOLS AND ABBREVIATIONS M Periodicity of the rotating polaroid f Frequency N Number of data points available for analysis DFT Discrete Fourier Transform MVU Minimum Variance Unbiased (estimator) CRLB Cramer-Rao Lower Bound pdf Probability density function MLE Maximum Likelihood Estimator LSEK Least Squared Estimator of K sinusoids APES Amplitude and Phase Estimator for Sinusoids FFT Fast Fourier Transform pfa Probability of false alarm ix Abstract Optical imaging is a vibrant research area, fueled by defense and bio-medical applications. The non-ionizing nature of near-infrared and optical waves, coupled with the possibility of performing functional imaging, has ushered in new interest in bio-medical optical imaging. The technologically advanced state of optical and infrared sources and the availability of fast and inexpensive detectors have accelerated the research in this field. The challenges encountered in optical imaging can usually be attributed to the phenomenon of scattering. Scattering blurs the details of an object being imaged, thereby reducing the achievable resolution. Many imaging schemes have been devised to circumvent this degradation induced by scattering. In one such class of techniques called as ‘Direct imaging techniques’, the main aim is to reject scattered light and use only unscattered light for imaging. The criteria for rejecting and accepting different parts of the radiation as scattered and unscattered can vary, and this has given rise to different imaging schemes. Of particular interest to us, in this thesis, are the continuous-wave, direct imaging schemes, which use polarization of the received radiation to discriminate the unscattered radiation from the scattered part. The emphasis in this thesis is on the signal processing methodologies adopted in such schemes. Throughout the thesis, we have used only linear polarization for our study. However, parallels can be drawn from this study, to the case of circular polarization for most of the situations discussed. A class of continuous-wave, polarization based direct imaging techniques use simple subtraction of sets of co-polarized and cross-polarized images to obtain an image corresponding to the unscattered component of light. These schemes are called as polarization x Abstract xi difference imaging (PDI) schemes. There are other schemes which use polarization modulation (either at the source or at the receiver), followed by sinusoidal estimation, to discriminate between scattered and unscattered components. These are called as polarization modulation imaging (PMI) schemes. In this thesis, we present a mathematical framework to analyze and compare these apparently disparate imaging methodologies. Theoretically and through Monte-Carlo simulations, we have studied the relative advantages and disadvantages of the PDI and PMI schemes under both white and coloured noise conditions, concluding that, in general, the PMI scheme gives better estimates of the unscattered component. The results of simulations and experiments corroborate our theoretical arguments. The PMI scheme is shown to give asymptotically efficient estimates of the unscattered component, whereas the PDI scheme is shown to give biased estimates of the same. Moreover, it is shown that PDI scheme is a particular case of the PMI scheme. We have also suggested that the PMI scheme, which modulates the received radiation rather than the incident radiation is a more useful variant, since it can be easily adapted for passive imaging as well. We conducted PMI experiments in the Optics lab of the Raman Research Institute, Bangalore, wherein, we used a 10 mW , 632.8 nm, linearly polarized, He-Ne laser source to image opaque objects through scattering slabs of mono-disperse polystyrene microspheres dispersed in water. Polystyrene particles with diameters 2.97 µm, 0.06 µm and 0.13 µm were used in our studies. An intensified charge coupled device (CCD) camera was used to capture the images. Results showed that imaging can be performed beyond 40 optical thicknesses, for particles of 0.13 µm diameter. For larger particles, the depth to which imaging could be performed, was much lesser. Experiments were also conducted to image an object through mist, with successful results. An experiment using an incoherent, white light source, showed that using incoherent sources can yield better imaging results than coherent sources. We attribute this to the speckle noise induced by coherent sources. We have also analyzed schemes which use the degree of polarization (DOP) of the scattered light as the visualization parameter. The expressions for DOP are available in Abstract xii the literature, only for the PDI scheme, and not for the PMI scheme. By simple theoretical analysis, we have shown that the DOP information can be obtained by PMI schemes also. We have compared the various methods of estimating the DOP information and conclude that the PMI and the PDI schemes are better than all the other schemes considered. We observe that the current PMI and PDI schemes cannot discriminate different states of linear polarization, which can occur in real data. We show that this can be achieved by using the phase information of the recorded sinusoids in the PMI scheme, at no extra computational cost, when estimating the unscattered component. We call such schemes as Polarization Orientation Imaging schemes (POI). However, the PDI scheme is incapable of estimating the phase of the sinusoids, which is its major drawback. We studied the suitability of using matched-filter based estimation techniques to estimate the various visualization parameters discussed, namely, the unscattered component, the DOP and the phase. We found that the PMI and the PDI schemes show superior performance than these schemes, except that the PDI scheme cannot estimate phase. However, there are a few conditions when the matched-filter based techniques can give better results, which we have mentioned. We have shown that the three visualization parameters can be fused to form a colour image, which gives a holistic view of the scene and mentioned the advantages of such a rendition. We also report the advantages of analyzing chunks of data and bootstrapped data under various circumstances, to estimate the various visualization parameters. We have also briefly touched upon the possible post processing that can be performed on the obtained results, and as an example, shown the segmentation of a POI result. Chapter 1 Introduction to optical imaging 1.1 Imaging science From the point of view of this thesis, we define ‘imaging’ as the creation of a visual representation of some measurable property of a person, object, or a phenomenon. Thus, ‘Imaging Science’ becomes the pursuit of scientific understanding of imaging or an imaging technique. This inter-disciplinary branch strives to view phenomena normally invisible to the human eye by translating them by some means, to a visually perceptible form. What makes imaging a very interesting, diverse and an inter-disciplinary field is, that the bearer of (or what conveys) the ‘measurable property’ can take many forms. For example, Ultrasound and X-rays are both used for imaging the interiors of the human body; though both are used for the same purpose, the information the two modalities give can be different. Similarly, visible, infrared and radio waves are used for remote sensing applications depending on the features one is interested in. The intelligence in designing new imaging techniques lies in perfecting the art of choosing a suitable bearer of the ‘measurable property’ for the application on hand. In this thesis, we deal with only those techniques, where the properties of interest of objects being imaged are carried by visible and near infrared waves. This field of imaging 1 CHAPTER 1. INTRODUCTION TO OPTICAL IMAGING 2 is called as Optical Imaging. Wherever generalizations to other regions of the electromagnetic spectrum are possible, a mention will be made. Having narrowed down to some extent towards the area of interest of this thesis, we delve deeper into optical imaging itself. 1.2 Optical imaging Optical imaging has diverse applications ranging from microscopy to astrophysics. Cameras, optical telescopes, optical microscopes are all optical imaging devices. Optical imaging is also a very important part of remote sensing. There are a few important reasons which make optical imaging desirable and sometimes inevitable. In general, the wavelength used for imaging limits the achievable resolution or details in image quality [1]. So, one would like to use the shortest wavelength at disposal to image objects if one could do so. In the case of remote sensing, one would have liked to use wavelengths below visible range for imaging, if earth were to reflect these wavelengths. Radiations with wavelength below that of the visible range are absorbed strongly by earth’s atmosphere, or else, our existence would have been in jeopardy. On the other hand, the visible, infrared and radio waves are easily allowed to reach the earth’s surface. The reflected waves of these wavelengths carry information about the composition of the earth’s surface. This makes the use of wavelengths shorter than that of the visible spectrum inviable for remote sensing. This reality has made remote sensing to exceedingly depend on visible, infrared and radio wavelengths. Similarly, since X-ray is an ionizing radiation, there has been a penchant among radiologists for using visible wavelengths for imaging the interiors of the human body. Currently, interest is growing in optical imaging and we may not be far from realizing diagnostic optical imaging tools. Besides these, defense needs have further pushed the limits of optical imaging. All these diverse applications have made this field, a vibrant one. CHAPTER 1. INTRODUCTION TO OPTICAL IMAGING 3 Although optical imaging can offer solutions to many important imaging problems, it is not easy to image using optical wavelengths. Optical radiation has been known to be one of the most difficult imaging probes. One important reason is light scattering, which manifests itself in most of the imaging conditions. Though the phenomenon of absorption (which is closely linked to scattering) also affects imaging in certain conditions, many a time, it becomes a useful property in imaging. The phenomenon of scattering and absorption affect all imaging systems to varying extents, and hence, have attracted immense interest from researchers across different disciplines. In this thesis, we neglect the effects of absorption and consider imaging affected by scattering only. More specifically, we deal with two-dimensional imaging, and with imaging techniques that deal with polarization information, apart from the intensity information. Having narrowed down to the problem further, let us see some aspects of light scattering. 1.3 Light scattering In this section, we dwell on a very simplistic picture of light scattering and look at the way it affects imaging. The problem of imaging through light scattering media is also termed as ‘Imaging through turbid media’ [2] or ‘Imaging through turbidity’. Note that we are not considering the effect of turbulence induced by wind etc., which we come across in literature related to atmospheric optics. The effect of light scattering is the only topic of interest in this thesis. To understand the reasons that make imaging under scattering conditions a formidable task, we consider a simple case, where, we try to obtain the shadowgram of an opaque object resting in a scattering medium, using active illumination (i.e., the light source is controllable). If the object is in an environment in which the scattering effects are negligible, shining light from one side would cast a sharp shadow of the object on a screen placed on the other side of the object. Fig. 1.1(a) shows such a simulated shadow of a elliptical object. 4 CHAPTER 1. INTRODUCTION TO OPTICAL IMAGING If the object is placed in a scattering medium, some light will enter the region of geometric shadow of the object, due to scattering. If the scattering increases, more light will stray into the region of geometric shadow. This decreases the contrast between the shadow and the surrounding background, introducing a blurring effect on the shadowgram. Fig. 1.1(b) shows what could be the result of mild scattering. If the concentration of the scatterers is increased further, the blur keeps increasing, till a point is reached, when the shadow becomes indiscernible with respect to the background. Fig. 1.1(c) shows the effect of increase in scattering on the shadowgram shown in Fig. 1.1(a). The images in Fig. 1.1 were generated by a simulation of scattering. Generally, to qualitatively assess the effect of light scattering, researchers resort to the photon picture of light. In this setting, light inside a scattering medium can be regarded as consisting of three different components, viz., the ballistic or the unscattered photons, the diffuse or the multiply scattered photons and the quasi-ballistic photons or weakly scattered photons. Among these components, the ballistic photons which travel straight and unscattered are the ones which are capable of forming shadowgrams of an opaque inclusion. The quasi-ballistic photons, which undergo scattering only through small angles can be imagined to meander about the forward direction like a snake and for this reason are also called as ‘snake photons’. Though the snake photons also blur the image of the object slightly, they retain their direction to some extent and are hence, useful. The (a) (b) (c) Figure 1.1: Effect of increase of number of scatterers on imaging - Simulated images. (a) Shadow of an elliptical object formed with negligible scattering (b) Effect of mild scattering. (c) Effect of intesnse scattering. CHAPTER 1. INTRODUCTION TO OPTICAL IMAGING 5 diffuse photons cause most of the blurring and are the ones to be avoided while imaging. The composition of light collected by the receiver has been analyzed in different ways, though, parallels can be drawn between the description given above and the ones given in [3, 4]. By a look at Fig. 1.1, it is evident that we could have obtained a sharp shadow of the object, if we could have prevented light from entering the geometric shadow region, i.e., if we could have retained the rays passing straight through the medium and rejected the scattered rays. Some imaging methodologies strive to achieve this single goal, i.e., to capture the ballistic and snake photons and avoid collecting the diffuse photons. These schemes are referred by a collective term, namely, ‘Direct imaging’ schemes [5]. In scattered polarized light, the ballistic and snake photons retain their initial polarization state to a greater extent, as compared to the diffuse photons. Hence, polarization of the scattered photons can be used to distinguish the less scattered ones from the rest. Such schemes are called polarization based direct imaging schemes. There are many other imaging schemes which analyze the diffuse intensity collected, and these are in general called as ‘Indirect imaging’ schemes. In these schemes, imaging is viewed as an inverse source problem and the wave equation is used for obtaining the results. In this thesis, we have considered direct, polarization based imaging schemes only. In our work, we have looked at two-dimensional imaging schemes only and especially, the transmission mode imaging schemes. Three dimensional imaging schemes have not been the subject of analysis in this thesis. But, a general review of the different schemes existing to date and in particular, a review of polarization based schemes is given in the third chapter. 1.4 Organization of the thesis The next chapter contains the essential background material related to polarization optics and light scattering, that is needed to understand the latter chapters. The chapter following the preliminaries contains a review of the field of optical imaging. The emphasis CHAPTER 1. INTRODUCTION TO OPTICAL IMAGING 6 is on continuous wave imaging using polarized light, the topic addressed in this thesis. The review mainly concentrates on the analysis methods adopted in such imaging schemes. The fourth chapter provides a mathematical framework for the analysis of various imaging schemes, using which, we assess the advantages and disadvantages of the various schemes. We have also proposed some new imaging schemes and processing techniques in that chapter. In the fifth chapter, we have illustrated the results obtained from actual experimental data. The final chapter summarizes the contributions of this thesis and provides pointers for further research on the topic. Chapter 2 Scattering of polarized light In this chapter, we give a gist of the terms from the field of polarization optics and light scattering, used frequently in this thesis, along with the bare-minimum mathematics related to them. The contents of this chapter have been compiled by borrowing heavily from [6]. 2.1 Polarized light Since the thesis deals with imaging using polarized light, we begin with the most basic notions of polarization of light, its representation and measurement. Polarization is a property which arises out of the transverse nature of the electromagnetic (EM) radiation, and is related to the orientation of the plane of vibration of its electric field. The magnetic field associated with the radiation is not taken into account to denote polarization. The compelling need to study the vectorial nature of light is encountered when one tries to understand light propagation in restricted media like optical fibers or in the presence of surfaces with discontinuous refractive indices (like scattering media), where the boundary conditions need to be taken into account. In the case of a plane surface, it can be shown that the component of the electric field perpendicular to the plane of incidence (TE wave) and the one parallel to it (TM wave) are completely unrelated, and their behavior during 7 CHAPTER 2. SCATTERING OF POLARIZED LIGHT 8 reflection and refraction phenomena are independent. So, it looks as though the EM field can be treated as a superposition of two plane-polarized waves, which can form a basis for generating any other polarization state. Since the monochromatic plane wave is the fundamental entity in the description of polarized light, we examine it in greater detail. 2.1.1 States of polarization Consider a plane monochromatic wave, with angular frequency ω and wave number k. The wave number k is defined as k = 2πn ν 2π = c λ (2.1) where, n is the refractive index of the medium and ν is the frequency of the EM radiation and λ is the wavelength in the dielectric medium and is given by λ = c/nν. Suppose the wave is propagating along the z direction in a non-absorbing medium. The electric field E of the wave can be represented mathematically as E = A cos (kz − ωt) + B sin (kz − ωt) where, the real vectors A and B are independent of position in the medium. The electric field vector at any point lies in a plane, such that, the normal to the plane will be parallel to the direction of propagation. In a particular plane, say z = 0, the tip of the electric vector traces out a curve: E(z = 0) = A cos (ωt) + B sin (ωt) The above equation describes an ellipse, or more aptly, the vibration ellipse, a sketch of which is shown in Fig 2.1. If A = 0 or B = 0, the vibration ellipse is just a straight line and the wave is said to be linearly polarized or plane polarized. The non-zero vector specifies the direction of vibration. If | A |=| B | and A · B = 0, then the vibration ellipse is a circle and the wave is said to be circularly polarized. In general, a monochromatic wave is elliptically polarized. CHAPTER 2. SCATTERING OF POLARIZED LIGHT 9 −B γ −A A B Figure 2.1: Vibration ellipse with ellipticity B/A and azimuth γ As can be observed from Fig. 2.1, a given vibration ellipse can be traced out in two opposite directions, clockwise and anti-clockwise. We adopt the convention, according to which, an elliptically polarized wave is said to be right-handed if the vibration ellipse is rotating in the clockwise sense as viewed by an observer looking towards the source. Apart from the handedness, the ellipse can be characterized by its ellipticity (the ratio of the length of the semi-minor axis to the semi-major axis) and its azimuth (the angle made by the semi-major axis to the horizontal of the chosen frame of reference). In the example vibration ellipse sketched in Fig 2.1, the ellipticity is given by B/A , the azimuth is γ and the intensity of the wave is given by I = A2 + B 2 . The four parameters, viz. handedness, ellipticity, azimuth and the intensity, called as the ellipsometric parameters, completely specify the state of polarization of a wave. Generalizing the above analysis, we can say that the correlation between any two orthogonal components of the electric field decides the polarization state of the light. If the orthogonal states are completely correlated, the light is totally polarized. If they are completely uncorrelated, the light is said to be unpolarized. If they are partially correlated, light is said to be partially polarized. Only a purely monochromatic source can emit totally polarized light. If a source exhibits spectral width, or, if a purely monochromatic wave passes through depolarizing media that introduce random phase shifts between two orthogonal states of the basis, the light becomes partially polarized. This is what is commonly encountered in practice. CHAPTER 2. SCATTERING OF POLARIZED LIGHT 10 There are different ways of mathematically representing polarized light [7]. The Stokes representation turns out to be the most convenient for our purpose, since it can easily handle the representation of both fully and partially polarized light, both of which are encountered in light scattering studies. Also, the practical measurability of the parameters of the representation and the additive nature of the Stokes parameters (which is explained shortly), coupled with the Mueller matrix transformation methods to study light scattering, qualify the Stokes representation as an ideal choice to study polarization in light scattering. Hence, we study the Stokes representation in some detail. 2.1.2 Stokes vectors We introduce the Stokes parameters along with their notations through the experiments by which they can be determined. We assume that a polarization insensitive detector is used to measure the different intensities in the experiments and that the various polarizers and wave plates are ideal and do not absorb light. To define the polarization states of the beam, we use a set of orthogonal axes eˆk and ê⊥ , which we refer to as ‘horizontal’ and ‘vertical’, respectively. Then, the electric field E can be represented as E = E0 exp (ikz − iwt); E0 = Ek eˆk + E⊥ ê⊥ Now, we make the following measurements on the given beam of light. We have omitted a scaling factor in all intensity calculations from field values, since it does not hamper the understanding of the actual idea. 1. The intensity falling on the detector with no polarizers in the path is measured, and the value is given by Ek Ek∗ + E⊥ E⊥∗ 2. Let a horizontal polarizer be introduced in the path of the beam. The intensity transmitted will be Ek Ek∗ . If a vertical polarizer is introduced instead, the intensity CHAPTER 2. SCATTERING OF POLARIZED LIGHT 11 measured will be E⊥ E⊥∗ . The difference between the two measured intensities will be Ik − I⊥ = Ek Ek∗ − E⊥ E⊥∗ 3. Let a polarizer be placed in the path of the beam, with its axis aligned at +45◦ to the horizontal. We introduce a new set of basis vectors ê+ and ê− which are obtained by rotating the eˆk vector by +45◦ and −45◦ respectively. The new basis will be 1 ê+ = √ (êk + ê⊥ ); 2 1 ê− = √ (êk − ê⊥ ) 2 In this basis, the electric field vector E0 can be written as E0 = E+ ê+ + E− ê− where 1 E+ = √ (Ek + E⊥ ); 2 1 E− = √ (Ek − E⊥ ) 2 The amplitude of the transmitted wave through the polarizer aligned at +45◦ will be E+ = √1 (Ek + E⊥ ), 2 giving an intensity I+ = (Ek Ek∗ + E⊥ E⊥∗ + Ek E⊥∗ + E⊥ Ek∗ )/2. Similarly, the intensity transmitted through the polarizer aligned at −45 ◦ will be I− = (Ek Ek∗ + E⊥ E⊥∗ − Ek E⊥∗ − E⊥ Ek∗ )/2. The difference between the two intensities will be I+ − I− = Ek E⊥∗ + E⊥ Ek∗ 4. Now, we introduce the right and the left circular polarizers. The respective basis vectors are given by 1 êR = √ (êk + iê⊥ ); 2 1 êL = √ (êk − iê⊥ ) 2 In this basis, the electric field vector E0 can be written as E0 = ER êR + EL êL , where 1 ER = √ (Ek − iE⊥ ); 2 1 EL = √ (Ek + iE⊥ ) 2 CHAPTER 2. SCATTERING OF POLARIZED LIGHT 12 The intensity transmitted through the right and the left circular polarizers will be IR = (Ek Ek∗ + E⊥ E⊥∗ + iEk E⊥∗ − iE⊥ Ek∗ )/2 and IL = (Ek Ek∗ + E⊥ E⊥∗ − iEk E⊥∗ + iE⊥ Ek∗ )/2, respectively. The difference between the two gives IR − IL = i(Ek E⊥∗ − E⊥ Ek∗ ) With these four measurements, we can get the Stokes parameters or the Stokes vector as I = Ek Ek∗ + E⊥ E⊥∗ (2.2) Q = Ek Ek∗ − E⊥ E⊥∗ = Ik − I⊥ U = Ek E⊥∗ + E⊥ Ek∗ = I+ − I− V = i(Ek E⊥∗ − E⊥ Ek∗ ) = IR − IL Observe that Q and U depend upon the choice of the horizontal and the vertical axes, but I and V do not. The sign of V signifies the handedness of the ellipse; positive stands for right-handedness and negative stands for left-handedness. The methodology of transformation of the Stokes Vector, when the frame of reference changes, has been detailed in [6]. Table 2.1 shows the Stokes vectors for the commonly encountered polarization states. The angles on the top of each of the linearly polarized vectors corresponds to the angle made by the polarization axis with the horizontal. One important observation that has to be made here is that the Stokes vectors do not form a vector space. To realize this, it is sufficient to notice that the Stokes vector of the linear polarized state at 45 deg is not the sum (ignoring scaling) of the Stokes vector perpendicular and parallel to the reference frame. The Stokes vector representation of other commonly encountered polarization states are given in [7, 8]. 13 CHAPTER 2. SCATTERING OF POLARIZED LIGHT Table 2.1: Stokes representation of the fundamental states of polarization ◦ 0 1  1   0 0 ◦   90  1   −1       0  0 ◦ +45 1  0     1  0 ◦ −45  1  0     −1  0 γ 1  cos(2γ)   sin(2γ) 0      Linearly polarized light Right   1  0     0  1  Left  1  0     0  −1 Circularly polarized light In the case of an ideal, strictly monochromatic wave, the four parameters are dependent, and it can be shown that I 2 = Q2 + U 2 + V 2 . The Stokes parameters of a quasimonochromatic beam (which, in general, is partially polarized) are obtained by taking the time averaged quantities over an interval long compared with the period, in which case, it can be shown that I 2 ≥ Q2 + U 2 + V 2 . The equality holds if the light is completely polarized. Accordingly, we can define Overall degree of polarization (DOP) = p Degree of linear polarization (DOLP) = p Degree of circular polarization = Q2 + U 2 + V 2 I Q2 + U 2 I V I Next, we review the terminology used in light scattering, relevant to this thesis. (2.3) (2.4) (2.5) CHAPTER 2. SCATTERING OF POLARIZED LIGHT 2.2 14 Light scattering The physical basis for scattering by any system can be traced to the heterogeneity of the system, either at the molecular level or on the scale of aggregation of many molecules. Since heterogeneity is a rule rather than an exception in the real world, scattering is an omnipresent phenomenon. The study of light scattering begins with the understanding of scattering by single particles within the framework of EM theory, based upon which, the more complicated phenomena are studied. When an obstacle is illuminated by an EM wave, electric charges in the obstacle are set into oscillatory motion by the electric field of the incident wave. The accelerated charges re-radiate EM energy in all directions. This secondary radiation is called the radiation scattered by the obstacle. Essentially, we can put it as, scattering = excitation + re-radiation. If the frequency of the scattered light is the same as that of the incident light, then the scattering event is called as elastic or coherent scattering. If the incident EM energy is transformed into thermal energy, then the process is called as absorption. Scattering and absorption are not mutually independent processes, though, depending on the situation, one may be more prominent than the other. In this thesis, we deal only with scattering by particles and neglect the effect of absorption. Particles in a collection are electromagnetically coupled. Each particle is not only excited by the external field, but also by the resultant field scattered by all the other particles. Since the field scattered by a particle depends on the total field to which it is exposed, a rigorous theoretical treatment of scattering by many particles is a formidable task. Considerable simplification results, if we assume single scattering, i.e., the number of particles in the collection is sufficiently small and their separation sufficiently large, so that, in the neighborhood of any particle, the resultant field due to scattering by all the other particles is small compared with the external field. With this assumption, the total CHAPTER 2. SCATTERING OF POLARIZED LIGHT 15 scattered field is just the sum of the fields scattered by the individual particles, each of which is acted on by the external field in isolation from the other particles. It is difficult to state precise general conditions under which the single scattering criterion is satisfied. Next, we study the case of single scattering in some detail. 2.2.1 Single scattering The essential problem in scattering theory is to find the scattered field as a function of direction, for the given incident field and scatterer. Let us consider the case of light of arbitrary wavelength being incident on an arbitrarily shaped particle. The total flux scattered by the particle in all directions can be considered to be the flux of the incident wave falling on a virtual area σs called as the scattering cross-section of the particle. i.e., Fscat = σs Fin where Fscat and Fin refer to the total scattered and incident flux, respectively. Similarly, the absorption cross-section σa is defined as the virtual area over which the real overall absorption of the incident wave would occur. The extinction cross-section σext , (a measure of the overall extinction suffered by the incident wave in the direction of transmission), is the sum of scattering and absorption cross-sections. i.e., σext = σs + σa For non-absorbing particles, σext = σs . These cross-sections depend upon the particle size, shape, orientation, wavelength of the incident light, the relative refractive index (RRI) between the scattering particle and the surrounding medium (m) and the polarization state of the incident light. The size parameter (x) of the scattering particles takes into account the relative dimensions of the size of the particle and the wavelength (λ) of the incident radiation. CHAPTER 2. SCATTERING OF POLARIZED LIGHT 16 e.g. for a spherical scattering particle of radius a, x is defined as x= 2πa λ (2.6) The same relation can be written as x = ka, where k is the wave number of the radiation in the medium surrounding the particle (see eqn 2.1). Though extinction depends only on the scattering amplitude in the forward direction, it is the combined effect of absorption in the particle and scattering in all directions by the particle. In our work, since we used spherical scattering particles, we restrict further discussions to scattering by spherical particles only. Mie theory helps to compute the scattering properties of spherical particles. The theory calculates the angular dependence of the scattered intensity for light polarized parallel and perpendicular to the scattering plane (the plane containing the incident and scattered rays). From this, the intensity for any polarization can be calculated. For particles much smaller than the wavelength, the scattered intensity is distributed more or less uniformly in all directions, in which case, the scattering is said to be isotropic (this is also referred to as Rayleigh scattering). On the other hand, for particles larger than the wavelength, light is scattered preferentially in the forward direction and the scattering is called as anisotropic (or Mie scattering). The reason for anisotropy can be attributed to the interference between the scattered waves. For particles much smaller than the wavelength (Rayleigh scattering), the scattered waves will be in phase in all directions and hence, we see nearly an isotropic distribution of intensity around the particle. However, for large scatterers, the forward scattered waves interfere constructively and the others interfere nearly destructively, resulting in anisotropic distribution of intensity around the scatterer. The anisotropy is specified in terms of the anisotropy parameter g (dimensionless), or the scattering indicatrix (or scattering diagram) [9]. The term phase-function is also used sometimes to denote the same quantity. We shall use the terminology ‘anisotropy CHAPTER 2. SCATTERING OF POLARIZED LIGHT 17 parameter’ in further discussions. The anisotropy parameter, g, is a measure of the amount of light retained in the forward direction after a single scattering event and is mathematically given as g = hcos(θ)i, where θ is the deflection taken by a photon after scattering. g is a function of the radius of the spherical particle, the wavelength of the incident light and the refractive index contrast between the particles and the surrounding medium. It varies from nearly 0 (isotropic scattering) to 1 (strictly forward scattering). Table 2.2 gives the anisotropy values for the spherical particles used in our imaging experiments. In all the experiments, the wavelength of the source used was 0.6328µm. The refractive index of the scattering polystyrene spheres was 1.59. The polystyrene spheres were dispersed in water (refractive index = 1.33). Hence, the RRI between the scattering spheres and the surrounding medium is 1.1955. It has been observed that scattering by a single particle or collection of identical particles does not decrease the degree of polarization of fully polarized incident light, though the nature of polarization of the scattered light may change [6]. But, upon scattering by a collection of non-identical particles, the incident polarized light can be rendered partially polarized due to depolarization. These features are independent of the specific nature of the particles. This is true, only if the assumption of single scattering is valid. It has also been observed that the scattering near the forward direction is always coherent due to the presence of unscattered light. Table 2.2: Anisotropy values of spherical polystyrene particles of different sizes, used in our imaging experiments Radius(µ) 0.03 0.065 1.485 Size parameter (x) Anisotropy(g) 0.396 0.027 0.858 0.128 19.6 0.81 18 CHAPTER 2. SCATTERING OF POLARIZED LIGHT When light is scattered by a collection of particles, one cannot simply add the scattered intensities from the individual particles to obtain the resultant intensity in a given particular direction. In such cases, the total scattered intensity is determined by the square of the absolute value of the total electric field. I(t, r) ∝| E1 (t, r) + E2 (t, r) + E3 (t, r) + . . . |2 , where Ei (t, r) are contributions to the total electric field from different scattering events. We can observe that the interference effects between different fields are also taken into the summation, apart from the individual contribution from each of the fields. When scattering becomes severe, we say the light is multiply scattered. Next, we briefly review the concepts of multiple scattering, a detailed treatment of which is given in [10, 11]. 2.2.2 Multiple scattering In a medium consisting of a large number of scatterers, the incident field undergoes recurrent random scattering before it exits the medium. In such multiple scattering media, light is assumed to propagate diffusively. The interference effects are assumed to be scrambled due to many random scattering events (except in special cases like backscattering). The position and time dependent intensity is described by the diffusion equation or some other simplification of the radiative transfer equation. Based on certain assumptions and observations, a multiple scattering medium is usually characterized by a few length scales. The characteristics of different media are usually compared on the basis of these length scales. The mean free path (MFP) is one such length scale, used to characterize the scattering process. For instance, the scattering mean free path, ls , is defined as the average distance between two successive scattering events, and is given by ls = 1 nσs (2.7) CHAPTER 2. SCATTERING OF POLARIZED LIGHT 19 where n is the number of scatterers per unit volume (number density) of the scattering medium and σs is the scattering cross-section of the individual scatterers. Similarly, the absorption and extinction mean free paths can be calculated by substituting σa and σext in the above equation instead of σs . The definition of diffusive transmission is itself based on these length scales [11]; e.g. Light transmission through a semi-infinite scattering slab of thickness L is considered to be diffusive, if the following inequality is satisfied. λ ≪ ls ≪ L ≪ l a where λ is the wavelength of the incident radiation, ls and la are the scattering and absorption mean free paths of the scattering slab, respectively. The various length scales used to describe radiative light transport and their significance have been detailed in [11]. In the diffusive regime, a wave can be assumed to have undergone scattering after traveling, on an average, a distance of ls ; but it hardly means that the direction of the wave travel is randomized. Due to the anisotropy of scattering, the wave may propagate in the near forward direction even after several scattering events. To accommodate this characteristic of scattering, a transport mean free path (TMFP) is introduced, which is defined as the average distance that the light travels before its direction of propagation is randomized. The TMFP, denoted by l ∗ is given by l∗ = ls (1 − g) (2.8) Clearly, for isotropic scattering l ∗ ∼ ls since g → 0, which seems logical, because, the wave can be scattered into any angle with almost uniform probability over 4π steradian. When the scattering medium is taken to be in the form of a semi-infinite slab of thickness L along the direction of incidence of light and infinite in extent in the transverse direction, we can obtain a convenient measure of the scattering encountered by light, as given by CHAPTER 2. SCATTERING OF POLARIZED LIGHT 20 the optical thickness of the medium, which is defined as τ= L l∗ (2.9) Consider the propagation of light through a semi-infinite slab of thickness L along the +Z axis. The incoming beam decays exponentially due to multiple scattering. The unscattered beam at a depth z in the medium is given by the Beer-Lambert’s law as [12] z I(z) = I0 exp (− ) ls (2.10) where I0 is the incident intensity of the beam. This exponential decay of the unscattered photons gives an idea of the difficulty involved in performing direct imaging in highly scattering media. The ratio of unscattered to scattered photons decreases very rapidly as the penetration depth increases and puts a limitation on the depth up to which direct imaging can be performed. 2.2.3 Mueller matrices and light scattering Most of the interactions of polarized light with optical elements can be expressed as a linear relationship; i.e., the output Stokes parameters can be obtained by a linear combination of the input parameters. The same technique has been used to study scattering too. Mueller matrices have been extensively used for this purpose. In any co-ordinate system, each optical component in a particular orientation can be mathematically represented by its appropriate, real, 4 × 4 transformation matrix called the Mueller matrix.. The elements of the Mueller matrix M depend upon the properties of the interacting element, frequency of light and in case of scattering, on the scattering angle. The Mueller matrices of some of the standard optical elements are given in [8, 7, 6]. If we denote the Stokes vector (SV) of the input light as Sin and the vector obtained after transmission through the optical element as Sout , then, Sout = MSin , where M is the representative 21 CHAPTER 2. SCATTERING OF POLARIZED LIGHT Mueller matrix of the element. In the case of multiple optical elements arranged in a cascade, the total effect of the cascade is determined by the product of the individual Mueller matrices, in the correct order. e.g. The Mueller matrix of an ideal linear polarizer whose transmission axis makes an angle θ with respect to the horizontal is given by  1 cos 2θ sin 2θ 0     2 cos 2θ cos 2θ cos 2θ sin 2θ 0 1   M(θ) =   2  2 sin 2θ cos 2θ sin 2θ sin 2θ 0    0 0 0 0 Now, if light with input SV Sin = (Ii Qi Ui Vi )T is incident on this polarizer, yielding the output SV (Io Qo Uo Vo )T , then the input and output intensities are related as 1 Io = (Ii + Qi cos 2θ + Ui sin 2θ) 2 (2.11) Thus, if a linearly polarized wave is incident on a polaroid rotating at an angular frequency ω, then, the transmitted intensity is of the form 1 Io = (Ii + Qi cos(2ωt) + Ui sin(2ωt)) 2 (2.12) where, ωt = θ is the instantaneous angle made by the polaroid pass axis with the plane of polarization of the incident light. Since, for a plane polarized light, I i = Qi and Ui = Vi = 0, we are left with 1 Io = Ii (1 + cos(2ωt)) = Ii cos2 ωt 2 (2.13) This expression is called as Malus’ law. Mueller matrices have been extensively used in the computation of Stokes vectors emerging as a result of scattering and are referred to as Scattering matrices or Phase matrices. CHAPTER 2. SCATTERING OF POLARIZED LIGHT 22 The advantage of the Mueller matrix formalism is that, it gives a simple means of determining the polarization of the scattered light, given any arbitrary polarization of the incident light. Scattering can depolarize the incident beam, mix its polarization states and change its direction. The elements of a scattering matrix depend on the size and shape of the scatterer, the refractive index contrast between the scatterer and the surrounding medium, the angle of scattering and the azimuth of the scattering plane. Details regarding the derivation of the scattering Mueller matrix can be seen in [6]. The Stokes parameters of the light scattered by a collection of randomly separated particles are the sum of the Stokes parameters of the light scattered by the individual particles. Therefore, the scattering matrix for such a collection is merely the sum of the individual particle scattering matrices. For any particle or collection of particles symmetric about the direction of incidence of the beam, the elements of the scattering Mueller matrix, should be independent of the azimuthal angle φ. If unpolarized light is incident on one or more particles, the Stokes parameters of the scattered light can be shown to be in general, partially polarized. This result shows that scattering is a mechanism for polarizing light. The degree of polarization of the scattered light depend on the scattering direction and has been found to be maximum when the scattered direction is normal to incident direction. A detailed study of Mueller matrices of spherical scatterers has been given in [6]. It has been shown that the matrix contains only four independent elements. These elements can be calculated precisely using Mie theory. The relevant mathematics and algorithms for programming have also been given in [6]. With this, we come to the conclusion of all the preliminaries needed to deal with studies on scattering of polarized light by spherical particles, which is the main work reported in this thesis. Chapter 3 Imaging through scattering media: A review The literature related to imaging techniques of interest to us can be found in diverse fields from astrophysics to microscopy. Due to the vast expanse of literature, we do not emphasize on any particular field of application, but we concentrate on some common ideas and techniques that permeate all these diverse fields. Any modifications that can be incorporated based on the knowledge of the specific domain of application can push the utility of these techniques further. However, there is no denying that biomedical applications have taken a lion’s share in this review. In this review, we concentrate on direct imaging techniques that use CW monochromatic, polarized light sources to obtain two-dimensional images of objects hidden in a semi-infinite slab of scattering media, by using polarization information apart from intensity information. We adopt a top-down approach to get to the imaging schemes of interest to us. So, first, we begin with a brief overview of the field of optical imaging. 23 CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 3.1 24 A brief overview of optical imaging In recent times, researchers have shown great interest in optical imaging schemes for biomedical applications, apart from its use in defense and atmospheric optics. As a testimony to this, one can cite journals that have earmarked some issues to highlight the progress made in the field. Some of these special issues are [13, 14, 15, 16]. There is an increasing acceptance of the possibility that optical imaging can play a complementary role to the existing biomedical imaging techniques, and in due course, it may even become a main-stream imaging technology [17]. The main advantages of optical imaging schemes as compared to other imaging modalities for biomedical applications are • The non-ionizing character of visible and near-infrared radiation. • The ability to perform very fast detection of optical signals due to the presence of fast and affordable detectors, coupled with the possibility of being able to leverage a host of contrast agents to image dynamic processes [18, 19]. Researchers have been successful in bringing optical imaging techniques nearly on par with other established imaging techniques, at least for imaging soft tissues like breast [20] and neonatal head [21]. The polarization based techniques have also found a niche application area in skin tissue polarimetry [16]. A broad overview of the various optical imaging schemes being developed for medical imaging applications has been given in [5, 22]. Coming to the applications of optical imaging in other fields, one can see optical techniques being used for imaging through rain and fog [23], through haze [4], underwater imaging [24] and many more. There are interesting results [25, 26] which claim that a strong multiple scattering wall, far from being a hindrance to imaging, could serve as a thin lens, which can produce a high resolution image of an arbitrarily shaped three-dimensional object hidden in it. Simple correlation techniques applied on the speckle patterns reflected from the scattering medium are shown to do the trick. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 25 Though the problem of imaging through scattering media looks similar to that of imaging through turbulence, there is no consensus as to whether the tools used for the latter are applicable to the problem on hand. Some researchers have mentioned that the pointspread function analysis, which is the basis of speckle interferometry is not applicable to the problem of image transfer through multiply scattering media, due to the vanishingly small angle of isoplanaticity for multiple scattering [25]. But, ample effort has gone into the analysis of point-spread function of multiply scattering media [9], and there obviously seems to be some use of it, contradicting the view of [25]. This, and many other inconsistencies which we come across as we proceed, emphasize the fact that the field of optical imaging is yet to find conclusive answers to many basic problems, and this makes it a very fertile field for research. Combining the results of different optical imaging techniques can yield better results than what could be obtained by employing any one technique in isolation [26], and at times, it may be inevitable too. In this thesis, we are concerned with imaging schemes that closely resemble the technique of polarization discrimination described in the articles [5] and [27]. We have studied continuous-wave (CW), polarization based imaging schemes only and not pulse based imaging schemes. Though in [5], the polarization discrimination scheme has been said to be able to image not more than about a centimeter of soft tissue, the logical extension of what can be done with intensity alone seems to have been underplayed. Since polarization data is a superset of intensity data, one can hope to use all the intensity based techniques and also see if more information could be obtained by polarization. Essentially, polarization data should in no way hamper what could have been done with intensity alone. If this argument is valid, polarization based schemes have much more to offer to the field of imaging. In short, the polarization based schemes will be useful in those circumstances, where the scattering is neither so less, that only intensity would suffice for imaging the hidden CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 26 object, nor so high, that no polarization information is left at the receiver to get additional information about the hidden object. Many defense applications fall into this category, e.g., target detection in fog, rain, snow or haze. In fact, the development of polarization discrimination imaging technique can be attributed to research in underwater imaging [28]. In many such applications, the source information will be totally unknown (as compared to medical imaging schemes, where the sources will be completely known). In such situations, polarization information can complement the intensity information. An important issue is the comparative performance of CW imaging schemes and pulse based imaging schemes. The latter schemes use ultra-short pulses and high speed gating techniques to capture the ballistic component of the scattered light [29, 30, 5]. The method that we are interested in, uses CW sources, and banks upon polarization of the unscattered light and spatial filtering to do away with scattered photons [28, 27, 31]. Though the ultrafast shuttering techniques are very effective and simple, they are prohibitively costly, as compared to the CW techniques. So, it is important to know if the CW methodologies can give performances comparable to that of the shuttering techniques. It has been experimentally shown that CW based schemes can achieve the same level of diffusive light rejection as pulse schemes [32, 33]. However, the article does not compare the schemes varying the pulse widths or the receiver acceptance angles, and other parameters that affect the imaging performance. Nonetheless, it is an extremely important result. One should only be worried whether the generalization of the results to all pulse widths and acceptance angles at all wavelengths is really valid, which is unlikely to be the case. However, many CW polarization based imaging schemes bank upon this result to vindicate their performance vis-a-vis the pulsed imaging schemes [27, 31]. Even otherwise, CW polarization based schemes are worth studying, because, there are many circumstances (mainly in atmospheric optics), where sunlight is the source, and the information about specific objects (usually man-made objects) in a scene can be obtained by analyzing the polarization information of the sunlight reflected from the scene. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 27 Next, we review some CW polarization based, direct imaging schemes. 3.2 Polarization based direct imaging schemes Any polarization based active, direct imaging scheme can be visualized to consist of five main sections, shown as different blocks in Fig. 3.1. We see as to what parameters in each of these sections affect the image quality. Figure 3.1: A block diagrammatic approach to polarization based imaging techniques • Source and polarization state generation optics The source can be coherent or incoherent. It can be monochromatic or polychromatic. The polarization state generators can either be integrated into the source, as in the case of polarized laser sources, or could be a separate section, consisting of polarizers and wave-plates, to generate the desired states of polarization. • Semi-infinite scattering slab with hidden inclusion The scattering slab is characterized by its thickness L or optical thickness τ , the density of the scatterers specified in terms of ls or l∗ , the size of the scatterers specified in terms of the size parameter x or the anisotropy parameter g, and the relative refractive index (RRI) m between the scatterers and the surrounding medium. The opaque inclusion is usually assumed to be totally absorbing and reflection from it are neglected. For applications like biomedical imaging, the subject of interest can be the variation of some parameter across the slab. In such cases, there is no separate inclusion to be imaged. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 28 • Collection and polarization analysis optics This usually consists of a spatial filtering section comprising lenses and apertures and a polarization analysis section contain polarizers and/or wave-plates. This part of the setup discriminates photons based on their direction of arrival and the state of polarization, so that, to the extent possible, only the unscattered photons are allowed to reach the detector. This unit is characterized by the focal lengths of the lenses being used, the sizes of the apertures, orientation of polaroids and wave plates, etc. • Detector The detector is usually either a scanning lock-in amplifier or a charge coupled device (CCD), though it can take other forms too. In case a CCD is used, its gain, signal to noise ratio (SNR) at various gains, the spatial resolution, quantization depths, sensitivity and integration time etc., affect the final results. • Data analysis and interpretation Data analysis and interpretation are based on measurements of features that discriminate regions of the scattering medium from that of opaque inclusion. The modeling of the scattering phenomenon plays a major role in deciding the analysis and interpretation algorithms. We segregate the review according to these blocks. For better flow in the presentation, we begin the review by studying the influence of the scattering slab parameters on the performance of imaging. After that, we review the literature associated with Source and polarization state generation optics. Later, we review the Collection and polarization analysis methodologies. Finally, we review different polarization based, CW imaging schemes along with their associated processing techniques. 3.2.1 Effect of scatterers on imaging In general, light scattering by a collection of particles can be categorized into three possible regimes, based on the amount of unscattered light present in the received radiation, and the polarization memory of scattered light [34, 12, 35]. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 29 • Regime 1 When scatterer concentration is low, a significant amount of unscattered light is detected, and also, a major part of the scattered light maintains its original polarization. Under such circumstances, polarization based imaging techniques are really not necessary, since, intensity based techniques followed by image restoration should suffice. • Regime 2 When scatterer concentration is such that no unscattered light is detectable, but the scattered light is still partially polarized, the polarization properties of the scattered light depend on the particle size. In this regime, the DOP of scattered photons varies with the number of scattering events undergone by the photons. Hence, DOP can act as a discriminant for segregating short and long path photons. Thus, polarization based schemes are advantageous in this regime. • Regime 3 When scatterer concentration is too high, unscattered light is not detectable, and also, all the scattered light is nearly totally depolarized. In such adverse scattering conditions, neither intensity based, nor polarization based direct imaging techniques are useful. To be able to image in this regime, one must use indirect imaging techniques like optical coherence tomography [36], diffusion optical tomography [18] or other methods [5, 37]. We summarize the above discussion in Fig. 3.2 by reproducing the sketch given in [12]. Henceforth, all the advantages and applications of polarization based imaging that we mention, will pertain to scatter concentrations where, the imaging condition can be categorized as belonging to regime 2 or 1. Though the classification of the scattering regimes as shown above looks attractive, there are no clear boundaries between these regimes; rather, there is a gradual transition CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 30 Figure 3.2: Summary of scattering regimes from one regime to another. Hence, this classification is hardly of any use, unless we define some measures that determine the scattering regime to which a collection of particles belongs to. Towards this end, a few mesoscopic length scales have been defined [35, 38, 39]. But, the major problem with these definitions is that, there is no way of theoretically calculating them from the slab parameters. These values are determined only experimentally. This restricts their utility to being usable only as a priori knowledge in designing new imaging systems. The following subsections discuss three important parameters of the scatterers that determine the scattering regime and influence the quality of imaging. 3.2.1.1 The volume and density of scatterers The scattering slab thickness L, which represents the volume and the density of the scatterers (which in turn decides the values of ls (eqn 2.7) and l∗ (eqn 2.8)), determines the scattering regime [11] and the amount of unscattered light left in the received radiation (eqn 2.10). The calculation of the DOP in scattered light also utilizes these length scales and is given in [40]. The characteristic length scales of depolarization for slab geometry have also been defined in [40]. These parameters give a first hand information about CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 31 the scattering regime that one may encounter, and also about the gross behavior of the scattered light. 3.2.1.2 The size of the scatterers Though the size of the scatterers is a very important factor, it can hardly be ascertained in most of the applications. This stresses the need to model the scatterers based on the application on hand, to obtain good results [41]. We now list a few important observations regarding the effect of size of the scatterers on imaging, assuming a fixed wavelength. In the ensuing discussion, we use the terms Rayleigh and Mie regimes to refer to situations where, the particle sizes are much smaller and larger than the incident wavelengths, respectively. 1. The spatial resolution achievable using polarization based techniques depends heavily on the scattering anisotropy of the medium [42]. 2. The size of the scatterers influences the ballistic propagation and depolarization in scattering media [35]. 3. As the diameter of the scattering particle increases, the number of scattering events required to degrade the ballistic propagation decreases [35, 43]. 4. The transport mean free path length of the scattered polarized component increases with the particle diameter [35, 43, 38]. 5. For small size parameters, the characteristic length of depolarization is independent of the RRI [44]. 6. Imaging through samples with large scatterers is restricted to lower optical thicknesses than for samples with smaller scatterers [35]. 7. For media with very high anisotropy parameter, polarization based techniques do not yield any advantage over polarization-insensitive techniques for the purpose of imaging [41, 34, 12, 35]. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 3.2.1.3 32 The relative refractive index of the scatterers We list a few important observations on the effect of RRI on imaging. 1. Different media with similar ls and g values can exhibit distinctly different polarization characteristics if the RRI are different [44]. 2. For small size parameters, the depolarization length is independent of the RRI [44]. 3. For large RRI and large size parameters, the anisotropy factor is diminished. This behavior is due to the increase of backscattering as the RRI increases [44]. 4. Irrespective of the size parameter, the circular depolarization length depends strongly on the RRI [44]. With this, we complete the review of the effect of the scattering slab on imaging. We now analyze the effect of source and polarization optics on the imaging performance. 3.2.2 Effect of source parameters and polarization optics The important variables that need to be considered here are the wavelength of the source, its polarization state, its spectral width and coherence. We neglect the effect of coherence in this review. We have already seen the influence of wavelength (by way of size parameter) on imaging performance. Hence, we restrict our review to analyze only the effect of polarization state of light on the imaging performance. The key results have been listed below. 1. Depolarization of scattered light depends on the initial polarization of the photons [42]. 2. Few scattering events are needed to randomize circular polarization in Rayleigh regime [42, 45, 46, 40] and linear polarization in Mie regime [42]. Many scattering events are required to randomize circular polarization in Mie regime [42, 45, 46, 40]. 3. In dilute suspensions of microspheres, where independent scattering can be assumed, the DOLP decreases as the scatterer concentration increases. However, for dense CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 33 suspensions, the DOP increases with the scatterer concentration, and further, preferential propagation of linear over circular polarized light is observed [45, 46]. 4. For biological tissues, the widths of the point-spread functions do not depend on whether the incident light is linearly or circularly polarized [41]. 5. Polarization discrimination of ballistic and snake photons based on circularly polarized light gives rise to a wider point-spread function than the one based on linearly polarized light [41]. 6. Circularly polarized light cannot be used to discriminate between weakly and strongly scattered photons for media containing spheres of large diameter. The converse is true in media containing small spheres [41]. 7. The DOP decays at the same rate for both incident linear and circular polarization states for small detector apertures and it is independent of the anisotropy parameter [41]. 8. In Rayleigh regime, the depolarization rates for both the incident linear and circular polarization states are nearly the same. 9. Circular depolarization length strongly depends on the RRI and systematically decreases after each scattering event, whereas, linear depolarization length decreases with the randomization of directions. This is true, independent of the size parameter and the refractive indices of the scatterers and the size of the receiver’s field of view [44]. 10. The criterion to choose the initial polarization state of the probe beam, when the optical properties of the scattering medium are known, is addressed in [39]. For Rayleigh scatterers, linear polarization has been recommended, and for Mie scatterers, circular polarization has been suggested. With this, we complete the review of the ways in which the polarization state of the input light affects imaging performance. It should however be remembered that only in active imaging schemes, we can choose the polarization state. Many applications use CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 34 passive imaging schemes, in which case, the above review will not be able to help predict the imaging performance. Now, we review the literature related to the fourth section of the schematic given in Fig.3.1. 3.2.3 Effect of spatial filtering and detector characteristics The variables in the detector section that affect the imaging performance are, the parameters of the spatial filter, like, focal lengths of the lenses and the sizes of apertures [47], the detector sensitivity, resolution (both spatial and quantization levels in the case of CCDs), noise characteristics, gain and the polarization analysis optics. The detector unit plays an important role in the case of CW imaging, since, apart from polarization retained in the unscattered or less scattered photons, which are weighted appropriately by the polarization optics, what really helps in rejecting diffuse photons is, the spatial filtering. The amount of light retained after spatial filtering is very low. Hence, detectors sensitive to very low light levels are needed. Studies reveal that polarization gating methods are superior to the pinhole gating method when signal strength is weak [48]. The polarization analysis optics will usually consist of polarizers and wave-plates (either fixed or variable), to weigh the collected radiation according to the state of polarization. Essentially, they work as gates that allow only particular kind of polarization and attenuate other polarization states heavily. The action of polarizers and wave-plates on radiations of different polarizations is well expounded in [7]. However, the effect of spatial filtering on scattered light is less well documented. We now look at some details of it. In spatial Fourier filters, adjustments of the aperture size at the Fourier plane can be made to act as a variable temporal gate for light emanating from the scattering media [47]. The above study deals with pulsed imaging schemes. However, the conclusions arrived at, that the diffusive component of the scattered light can be removed by angle gating using spatial filtering, holds even in the case of CW imaging. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 35 In the same vein, there have been experimental and theoretical studies on the detection of ballistic and rejection of diffuse light in trans-illumination confocal and heterodyne imaging systems [49]. From these studies, expressions for optimum pinhole size for ballistic light detection and diffuse light rejection for confocal imaging have been derived. For large detector apertures, even for small optical depths, the blurring effects due to multiple scattering become evident, since, light from wide angles are collected and in this case, the polarization-difference imaging schemes have been found to be more effective [41]. Detectors with high sensitivity and gain are usually needed for direct imaging. As for any other application, detectors with high spatial resolution and greater quantization levels are desirable. Next, we review some polarization based, CW direct imaging techniques and associated processing algorithms. 3.2.4 Imaging schemes and processing algorithms There are mainly three variants of CW polarization based imaging schemes, depending upon the parameter of visualization adapted in the scheme. The first scheme is based on orthogonal polarization differencing. The second scheme uses polarization modulation, finally leading to sinusoidal amplitude estimates as the visualization parameter. These two schemes are in fact interchangeable, and as we will see in the following chapter, the former is a particular case of the latter. The third group of schemes is based on the principles of ellipsometry, wherein the schemes make use of Mueller matrix estimates of the scattering medium for the purpose of visualization. In this survey, we have emphasized on the first two classes only. However, we have mentioned some of the aspects of the schemes based on ellipsometry. To begin with, we present the PDI schemes. The first few schemes have been explained in detail, and for the ones that follow later, we mention only the differences, since, the differences between the implementations are very little. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 3.2.4.1 36 Polarization difference imaging schemes One of the earliest works on PDI schemes is reported in [28, 3]. These articles elucidate the ideas behind PDI, the algorithm employed for reconstruction of the target and the SNR analysis. There is a brief mention of the possible applications of the technique in [3], which vary from polarization based microscopy to runway lighting systems. The experimental setup used in [28] is a passive PDI system. The polarization information obtained after reflection from the target is used for reconstruction of the target image. The actual imaging was performed as follows. A polarizer/analyzer pair was used to record images parallel and perpendicular to the fixed analyzer axis denoted as Ik and I⊥ , respectively. These two formed the orthogonal polarization state images. In order to reduce the noise variance at each pixel, a set of 128 such images were obtained at each polarization state, added and suitably scaled. It is shown that the polarization difference image Ik − I⊥ , when enhanced, could give the features of the aluminium target, along with the locations of the two abraded patches in it. Following the processing behind the technique, it was named as polarization-difference imaging. The reason behind the success of the system is hypothesized as follows. The image Ik is formed due to both the scattered and the unscattered light, whereas, the image I⊥ is formed predominantly by the scattered light (since the unscattered light would be removed by the polarizer-analyzer combination). So, when we subtract I⊥ from Ik , what we are left with, is the image formed due to unscattered light alone. In short, the polarization difference image is the result of common-mode rejection, i.e., the scattered light, which is common to both I⊥ and Ik , is rejected by the act of differencing. In fact, this is the argument used in all PDI schemes and is perfectly valid. The efficacy of the system has been explained on the basis of the observed degree of linear polarization (ODLP ), which is defined for a region as hODLP iregion PD hIregion (x, y)i = PS hIregion (x, y)i (3.1) CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 37 where I P D (x, y) = Ik (x, y) − I⊥ (x, y) and I P S (x, y) = Ik (x, y) + I⊥ (x, y) Vindicating the applicability of the technique to real, natural settings, it is mentioned that the PDI scheme was found to give encouraging results for experimental conditions with ODLP ≤ 0.01, which is much lower than what one can expect from object surfaces in natural environment. This scheme is very simple, passive and potentially very fast. Besides this, the same analogy can be carried forward to any wavelength of the electromagnetic spectrum. In case there is a need, further processing like histogram equalization etc. can be performed on the PDI results. The only possible bottleneck could be the low light levels that one usually encounters in scattering conditions, which may slow down the process of acquiring images. This stresses the need for having highly sensitive detectors like intensified CCD or its variants for the success of the scheme. Since the processing at each pixel is independent of the other pixels, the whole scheme can be parallelized, making the analysis very fast. For assessing skin lesions in superficial epidermal and papillary dermal layers, an imaging modality which uses a video camera, wherein the mechanism of contrast is governed by the reflectance of polarized light has been described in [50, 51]. In this scheme, the nearly unpolarized, diffusively reflected light from deeper layers of skin is rejected in favor of polarization retaining light backscattered from superficial layers, where most of the skin lesions occur. Care is taken to avoid specular reflection from the skin surface. This system has demonstrated the ability to visualize the true margins of skin cancer, that were not easily discernible by dermatologists by using a simple, incoherent white-light source. The resultant images are also supposed to subtract melanin pigmentation from the images of pigmented skin lesions, revealing the underlying structure. The resultant image is obtained as not just the difference of the images as in [28], but as DOP (eqn 2.4) at each pixel. In this case, the eqn 3.1 is explained as follows. The numerator symbolizes suppression of highly scattered light whereas the denominator is supposed to ensure the cancellation of the attenuation due to pigmentation of melanin. It is also described as to how the polarization based camera could improve the estimate of the size of a sclerotic CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 38 basal cell carcinoma, thereby proving useful in performing surgical excision of such cells. The influence of particle size on active polarization discrimination imaging in underwater applications has been studied using Monte-Carlo techniques [52]. Both circular and linear polarization states have been considered in the simulations. The simulations mimic the situation of a depolarizing target being illuminated by linearly or circularly polarized light, through a scattering medium. The depolarized light reflected from the object is the quantity of interest in this problem. In order to view the object, one must reject the light scattered from the medium, but retain the light scattered from the object. It is mentioned that the light backscattered from the medium with few scattering events will retain its original polarization to a great extent, whereas the light reflected from the object would be totally depolarized. Hence, viewing the medium from a state, orthogonal to the incident state of polarization should suppress the radiation scattered from the medium, but retain the radiation reflected by the target. This is what is leveraged in this scheme. This scheme is shown to work even for trans-illumination and many off-axis geometries. The article gives the following important results. 1. In general, the effectiveness of the orthogonal polarization technique decreases with increasing scattering particle size. 2. There is a clear maximum in the contrast enhancement at a particular depth (depending on particle size). 3. The depth at which the orthogonal polarization technique is most effective increases with particle size but with an improvement factor that reduces with particle size. This, though is generally true, will fail when the radius of particles nears the wavelength of illumination. In this condition, the scattering properties will show a nonmonotonic behavior. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 39 Another reflectance mode, active polarization-difference imaging scheme, useful for underwater applications has been discussed in [53]. Enhancement in target details is obtained by illuminating the scene with a source of known state of polarization and detecting the reflected light orthogonally polarized to the incident state. The authors refer to this as polarization discrimination imaging. It has been shown that the effectiveness of this technique depends very weakly on the particle shape and on the form of illumination geometry as mentioned in [54]. The authors also endorse the observations of [52]. They also propose a image subtraction based technique that extends the visibility depth further than a factor of 2, achievable by observing the cross-polarization image. The difference between the scheme mentioned in [53] and that of [50, 51] is that, though, in both methods the co-polarized and the cross-polarized images are obtained, in [53], it is the cross-polarized image that contains the information of the target, whereas, it was the co-polarized image that carries information in the experiments described in [50, 51]. The important thing to note here is that, in [53], the co-polarized and cross-polarized images do not contain the same amount of scattered light. Hence, subtracting co-polarized image from the cross-polarized image does not provide the details of the target. It is a fraction of the co-polarized image that needs to be subtracted from the cross-polarized image to obtain the desired result. Through simulation results, it has been reported that the optimum subtraction fraction is found to be ∼ 13%. By simulations, it has also been reported that the subtraction technique can improve the visibility depth by ∼ 1 mfp (scattering mean free path), which is a significant improvement for the application envisaged. Also, there are comparisons of this scheme with two other subtracted-image enhancement schemes, both of which are variants of unsharp masking. It is reported that, among the three enhancement schemes, the polarization based scheme gave the best results. Some novel approaches to tissue optics by laser light scattering have been studied in [37]. Three different techniques have been reviewed, and imaging based on DOP as visualization parameter is one among them. The techniques have been compared with traditional CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 40 approaches to diagnostics and imaging of macroscopically inhomogeneous, multiply scattering objects through simulations and experiments as well. Here are the highlights of this article. 1. For the intermediate scattering regime (l ∗ ∼ L), the ballistic and snake photons contribute significantly to the detected light. Hence, visualization and location of the absorbing inhomogeneity by using DOP as a parameter would yield results better than using intensity alone. 2. The efficacy of the polarization based schemes increases as anisotropy decreases, as mentioned in [52]. 3. The maximal sensitivity that can be achieved with the DOP as the visualization parameter should be expected if the modal value of the effective path length is of the order of depolarization length and this can be taken as the criterion for the applicability of the technique. Results of experimental and simulation evaluations of the possibility of using linearly polarized CW laser radiation for imaging of absorbers embedded in a multiple scattering medium is given in [55]. The transmitted light is analyzed in terms of DOLP (eqn 2.4). The contrast and effective signal to noise ratio for images reconstructed using DOLP as visualization parameter have been discussed. A parameter to define the sharpness of the resulting images has been quoted from earlier literature. Here too, the polarization based schemes can give better results than intensity based schemes in the intermediate scattering regime. With this, we conclude the review of PDI schemes. Essentially we found two ways of interpreting the result. One is, by the method of difference of the co-polarized and the cross-polarized images (or their fractions) and the other is, by calculating the DOP. We compare the two schemes in the next chapter. Next, we review another polarization based imaging scheme, which uses sinusoidal estimation for visualization. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 3.2.4.2 41 Polarization modulation imaging schemes As a parallel development to the PDI schemes, some researchers were working on the polarization modulation imaging (PMI) schemes. One such scheme is reported in [31]. In this scheme, the source polarization state, which is taken to be linear, is made to rotate (i.e., modulated) at a known rate. The scattered light is filtered through a fixed analyzer, so that, the intensity due to the ballistic and snake photons gets sinusoidally modulated at twice the rate of modulation of the source polarization (see Malus’ law, eqn 2.11, 2.13). However, the diffusely scattered light does not show such a variation, since the diffusely scattered light will either be unpolarized or randomly polarized. The contribution of the ballistic component to the overall signal is usually small compared to the contribution of diffusely scattered light, and is given by Beer-Lambert’s law (eqn 2.10). The sinusoidal component in the received signal is detected using lock-in detection. The resolution of such a scheme is limited only by the diffraction phenomena (i.e., the theoretical limit) and hence is much better than the resolution that can be obtained with photon density wave imaging [18], for which, the effective wavelength is much longer. The above method was considerably improved in [27, 2, 56]. The data acquisition was hastened by replacing the scanning lock-in detector by acquisitions with a CCD camera. The apertures were replaced by a spatial filtering scheme. Since the acquired images were processed offline, the need for a lock-in amplifier was eliminated. These simplifications rendered the PMI scheme very simple and fast and comparable to the PDI schemes mentioned in [28, 3]. Now, we take a brief look at the data acquisition and processing scheme adopted in [27]. Linearly polarized light, whose plane of polarization was controlled by a rotating polaroid was made to fall on a scattering slab containing an opaque inclusion. The rotating polaroid, which was mounted on a stepper motor was moved by a known, fixed angle at each step. At each step, an image of the scattered light filtered through a spatial filter and a fixed polarization analyzer was captured by an intensified CCD. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 42 Since the ballistic and snake photons retain their initial polarization even after scattering, their contribution to the images should vary according to Malus’ law (eqn 2.13), as the plane of polarization of the incident light is rotated. So, if we take a series of images with the plane of polarization varying by fixed angular increments at each step, then the series of values at each pixel should show a sinusoidal variation if the ballistic and snake components are high enough. The diffusive component hardly contains any polarization information, and that is random. Hence, its contribution to the series of values at each pixel is noisy. Now, the task of finding the ballistic and snake photon components in the signal boils down to finding the amplitude of the sinusoid at each pixel. The opaque inclusion essentially blocks the ballistic light from reaching the CCD camera. Hence, at the pixels lying in the geometric shadow region of the inclusion, we do not expect to find any ballistic component, i.e., the sinusoidal component is almost zero. At other pixels, a sinusoidal component will be present. So, a measure of the ballistic component (viz. the sinusoidal component) recorded at each pixel location of the series of images should be able to give the shadow of the opaque inclusion. The estimation of the sinusoidal component at each pixel is performed in [27, 2, 56] with the ubiquitous tool of Fourier transform. Since the rate of change of intensity due to ballistic and snake photons is twice that of the input state of polarization (ω), we need to find the 2ω component of the Fourier transform at each pixel. The resultant gray-scale image is constructed with the intensity at each pixel being equal to the magnitude square of the 2ω component. A suitable scaling may have to be applied to the image, to be able to perceive it. It has been reported that the method has been able to give good results for objects hidden in turbid media with turbidity as high as 30l ∗ , with a spatial resolution of about 100 µm, using a continuous, 1 mW laser source. It is important to note the subtle difference in the principle of polarization discrimination as given in [28], as compared to that given in [31, 27, 2, 56]. In [28], the polarization of CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 43 the light reflected from the target is constant and at the detector, the analyzer is rotated to acquire the co-polarized and cross-polarized images. But in [31, 27] the analyzer is held fixed and the source polarization is modulated. Hence, the later technique can be used only in the case of active imaging, where the source is controllable, whereas, the former can be used even in passive imaging schemes. Perhaps, one would get the same results by having a fixed source polarization plane and a rotating analyzer, even in the imaging setup described in [31, 27]. This is an important change that can extend the utility of the PMI procedure. The two classes of imaging schemes reviewed till now, i.e., the PDI and the PMI, have utilized the polarization of the scattered light to visualize the hidden object. However, no property of the inclusion itself can be discerned by these methods. The presence or absence of an inclusion and in case an inclusion is present, its dimensions and some properties (as in the case of the two patches described in [28]) are the only information that can be obtained by these methods. However, some imaging schemes which study the polarization properties of a scene, instead of the polarization state of light, have been developed [57]. These are termed as polarization diversity active imaging (PDAI) schemes. We next study two such schemes, which are extensions of ellipsometry to active imaging and rough surface optics. We have reviewed the PDAI schemes for the sake of completeness only, and we mainly concentrate on PDI and PMI schemes in the following chapters of the thesis. 3.2.4.3 Polarization diversity active imaging schemes In PDAI scheme defined in [57], a scene is illuminated with a sequence of polarization states and the measurements of polarization state scattered from the scene are captured as images. These images are then analyzed to determine the Mueller matrix at each pixel. The Mueller matrix of the target is obtained using the Mueller Matrix Imaging Polarimeter (MMIP), which uses the dual rotation retarder technique [58]. The authors explain as to how, depolarization of light can be used to estimate target roughness and texture. The CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 44 information obtained by this scheme has been used to measure the orientation of bodies and estimate their refractive index. The defense related applications of PDAI, like de-camouflaging and target detection have been studied [59]. The main idea in this scheme is to use the DOP of the reflected light (defined differently from eqn 2.4), to distinguish man-made objects form natural objects. In this technique, the fact that man-made objects depolarize incident polarized light to a lesser extent as compared to natural objects has been leveraged. The authors demonstrate the idea by using the Dual Rotation Retarder Technique (see [57]). In this scheme, the polarization degree at every pixel location is characterized by its Mueller matrix M. This is different from the way the DOLP was defined in the PDI schemes (see eqn 2.4). The polarization degree P d is defined by the authors as P d = 100 sP 3 i=0 P3 M2ij − M200 3M200 j=0 and intensity is defined as I = M00 where, Mij stand for the elements of the Mueller matrix M. The degree of polarization P d varies from 0% (corresponding to a totally depolarizing target) to 100% (a nondepolarizing target). The DOP as defined above is independent of the actual incident and received intensities (it depends only on the Mueller matrix elements), but the DOP as defined in [28] depends on the actual received intensities. The method proposed in fact captures all the properties of the medium through which the light scatters. Hence, this may be a better representation of the DOP, than that defined in [28]. The drawback of the scheme is the huge computational cost involved in finally obtaining polarization images, as compared to the scheme described in [28]. The results are usually represented in pseudo-colors for better visualization. CHAPTER 3. IMAGING THROUGH SCATTERING MEDIA: A REVIEW 45 The technique of Mueller Matrix Imaging Polarimetry (MMIP) has also been used for characterizing tissue properties [60], through which, the results of [50] have been verified. MMIP acquires images of samples and calculates the full Mueller matrix for each pixel of the image. Since the Mueller matrix completely characterizes the polarization properties at a pixel, all the information about the different regions of the tissue can be obtained and thus, the regions can be characterized. It has been shown that MMIP can be used as a technique for characterizing various dermatological diseases. Since the processing schemes discussed so far do not need information about the wavelengths involved, all the imaging techniques discussed so far now can be extended to other regions of the electro-magnetic spectrum too. They can also be modified to handle circularly polarized light. Our survey of the various polarization based imaging techniques concludes here. Chapter 4 Processing polarization-rich data There are mainly two imaging techniques which make use of the polarization of the scattered light to distinguish different regions in a scene, namely the Polarization difference imaging (PDI) and the polarization modulation imaging (PMI). Till now, to our knowledge, there has been no comparative study of these two schemes. In this chapter, we embark on such a task. We compare the theoretical and experimental performance of the PMI scheme described in [27] and the PDI scheme described in [3]. We also propose few other processing techniques and theoretically study and compare them with these two existing schemes. We propose a minor change to the PMI scheme proposed in [31, 27], to make it suitable for both active and passive imaging. The necessary change is to keep the plane of polarization of the incident light fixed and allow the analyzer to rotate. Also, we can observe that the PDI scheme becomes a particular case of this modified PMI scheme as explained below. In the modified PMI scheme, if images are captured with angular displacements of π/2 of the analyzer, the frequency of the resulting sinusoid due to unscattered light will be half the rate at which the images are captured (the sampling rate), i.e., only two points of one whole period of a sinusoid will be sampled. If one of the sampled points (images) is chosen to be at the maximum of the sinusoid (the co-polarized image), the next sample will naturally be that of the minimum of the sinusoid (the cross-polarized image). Thus, 46 47 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA we can get the polarization difference data by properly choosing the sampling point. In the rest of the thesis, whenever we refer to the PMI scheme, we assume such a modified PMI scheme, as given in Fig. 4.1. LASER BE FP O SF RP CCD Figure 4.1: The modified PMI imaging scheme. BE: Beam Expander, FP: Fixed Polaroid, O: Object immersed in scattering medium, SF: Spatial Filter, RP: Rotating Polaroid, CCD: Charge Coupled Device In order to compare different processing schemes, we utilize the concepts of estimation theory, with the assumption that the imaging methodologies being considered are essentially different estimators for estimating the same quantity. The characteristics of the estimators that are of concern to us are the bias and the variance. An estimator is said to be unbiased if it yields the true value of the parameter being estimated on an average, over all possible values of the parameter being estimated; else, it is termed as biased [61]. The variance of the estimator is a measure of the closeness of the estimated values to the actual value and is the criterion usually used to rank the estimators. Another important characteristic of an estimator is its efficiency, i.e., how good the estimator is, in using the available data to estimate the unknown parameters. The efficient estimator is the best, minimum variance, unbiased estimator one can hope to design. From the point of view of data processing and visualization, we identify three different parameters, which can be used to interpret the polarization information contained in the received radiation. We classify the imaging methodologies based on these visualization parameters as • Polarization intensity imaging (PII) If we can segregate the information in the received radiation as corresponding to CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 48 the polarized and the unpolarized parts (which we show to be possible), we can use the amount of polarized light as a visualization parameter. We call such schemes as belonging to the category of polarization intensity imaging. The PMI described in [27] can be taken as an example of such a scheme. • Degree of polarization imaging (DOPI) In these schemes, the quantity of interest is the degree of polarization. It is a measure of the purity of the received partially polarized radiation and may convey information about the nature of obscured objects, apart from their location. In our comparisons, we consider only the degree of linear polarization (see eqn ??). This parameter has been used in [28] to visualize the polarization information. • Polarization orientation imaging (POI) This scheme is useful only when linearly polarized light is used for imaging. When the received radiation contains multiple linearly polarized states, we may be able to identify and categorize the obscured objects based on the orientation of the plane of polarization of the received radiation. To the best of our knowledge, this has not been utilized in any of the imaging schemes, though, as we will show, this can be obtained at no extra cost, when we use the PMI scheme. In order to be able to use the concepts of estimation theory, we need to model the data captured at every pixel location, in each of these imaging schemes. We resort to the Stokes vector analysis for this purpose. As a precursor to modeling the data, we would like to stress the point that the processing involved in all the imaging schemes discussed in this thesis is essentially 1-D. A typical, polarization based, direct imaging technique yields a series of images, with the images being captured at particular states of the analyzer optics. The processing is performed at a particular pixel location, by looking at the series of data values at that location, across images. Every processing operation is applied to all the time series. After estimating the visualization parameter at every pixel location, we construct the resultant image by plotting the estimated values as an image. One important factor that CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 49 can affect the resultant image, is the variance of the estimated quantities at each pixel. If the estimated quantities show a large variance, we may obtain drastically different results when we repeat experiments with the same imaging setup. Estimation theory comes to our rescue here, since the variance of the estimators can be used to judge the performance of the processing schemes. Though the actual value of the estimated quantities are very important, in many cases, it is the estimator’s variance that is crucial. The resultant image is usually segmented into two or more regions, corresponding to the hidden objects and the background. In such cases, the true values at each pixel location will not be very crucial, unless they affect the segmentation process itself. In our analysis of the processing schemes, we assume that the data values are recorded using an ideal detector. i.e., we have ignored the effect of quantization of the actual data when the information is converted to an image. Also, the detector is assumed to have infinite dynamic range. We assume that these factors affect all the processing schemes equally, and hence, the best processing scheme should remain so, even when applied on quantized and clipped data. In all further discussions, we also assume that the incident light is linearly polarized. Extension of the analysis to circular polarization wherever suitable, is usually simple, and is mentioned later. The case of passive imaging also assumes that the information about the object of interest is parameterized by the linear polarization information present in the received radiation. With this background information, we model the intensity recorded in the PMI scheme, (shown in Fig 4.1) as follows. 4.1 Signal modeling In general, the Stokes vector (SV) recorded by a point detector (or at a pixel location of a CCD camera), can be represented as [hIs i hQs i hUs i hVs i]′ where h·i represents time averaging. Though, due to finite detector area, there is spatial averaging also, CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 50 we are more interested in time averaging, since the areas involved are small and the dominant factor affecting the recorded SV turns out to be time. The statistics of the time averaged, scattered Stokes vectors is of importance to diverse applications and is yet to be studied fully [38]. But, since the integration times are usually large compared to the coherence time (a measure of the duration of constancy of the instantaneous SV of the wave field) [38], it is usually assumed by virtue of the central-limit theorem, that the time averaged Stokes parameters recorded during different sub-intervals are statistically independent, Gaussian random variables. However, for some experimental data, the power spectral density of the time series revealed an underlying coloured noise process, which can be modeled by an Auto Regressive 1 (AR1) process [69]. Hence, we do not assume the noise to be white, in modeling the observed data. With this information, we can model the intensity recorded at an arbitrary pixel location at the nth step of acquisition of PMI scheme as (from eqn 2.12) 1 Ir (n) = (Is + Qs cos 2θn + Us sin 2θn ) + v(n) 2 n = 0, 1, 2, · · · , N − 1 (4.1) where θn is the angle made by the analyzer with respect to the horizontal at the nth step and v(n) is an AR1 process, given by v(n) = av(n − 1) + w(n) (4.2) where, w(n) is a zero mean, independent and identically distributed (iid) Gaussian noise process with unknown variance σ 2 and a is the unknown AR1 coefficient. In case a = 0, the noise process becomes white. Henceforth, we deal with the analysis of a single time series, representing data at any pixel. The analysis essentially applies to all other pixel locations as well. 51 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA Assuming that the rotating polaroid takes M steps to complete exactly one rotation, θ n increments by 2π c M at each step of the polaroid. Suppose we start grabbing images starting with an arbitrary orientation φ of the rotating polaroid, the intensity recorded at the n th step can be represented as 1 Ir (n) = 2  Is + Qs cos     4πn 4πn + 2φ + Us sin + 2φ + v(n) M M (4.3) A more useful representation of the same equation would be 1 Ir (n) = 2  Is + where, α = arctan  Qs Us p Q2s + Us2 sin  4πn + 2φ + α M  + v(n) (4.4)  . We observe from eqn 4.4 that the component of the intensity that is independent of the Is , 2 which corresponds √ 2 to 2 the diffuse part of the scattered Qs +Us light. The amplitude of the sinusoidal part, i.e., corresponds to the ballistic and 2 orientation of the analyzer is snake components and is a measure of the magnitude of the polarized light in the received radiation. We denote the former by Iscat and the latter by Ibal . Using this representation, Ir (n) = Iscat + Ibal sin  4πn +β M  + v(n) n = 0, 1, 2, · · · , N − 1 (4.5) where we have replaced the term 2φ + α by a single variable, β. The discrete frequency f of Ibal is given by f = 2 . M In all our comparisons of the imaging schemes, we assume that there are N images available for analysis. i.e., in case of PDI scheme, there would be N 2 images each, corre- sponding to the co-polarized and cross-polarized data. In PMI schemes, there would be N images constituting a time series at every pixel location, which would be analyzed for estimating the sinusoidal component. We further assume that N is an integral multiple of M . Although this condition is not very strict as N increases, for the sake of analysis, we continue with this assumption. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 52 The comparison of the analysis schemes needs the knowledge of the noise statistics at each pixel location, which is seldom known a priori. Still, we assume that the noise characteristics are known and analyze the various schemes, since we get an idea of the performance of the various estimators given a particular noise condition. We do not explicitly estimate the noise variance terms, since the quantity of interest to us is the unscattered component of light and noise is a nuisance parameter. Each imaging scheme is analyzed for two cases: (i) v(n) is white, and (ii) v(n) is coloured. When we assume v(n) is white, we replace v(n) by w(n) for clarity. We now analyze imaging schemes that exploit polarization intensity. 4.2 4.2.1 Polarization intensity imaging Polarization intensity imaging - the case of white noise We observed earlier that the PMI scheme described in [27], belongs to this class. If the processing of data in PDI schemes is restricted to taking the difference of the intensities alone, the PDI schemes can also give the polarization intensity information. The difference between the maximum and minimum intensities of the PDI data gives twice the amplitude of the sinusoid buried in noise. We shall compare the performance of such a PDI scheme, with that of the PMI scheme. 4.2.1.1 Intensity imaging using PDI estimator The analysis of the PDI scheme can itself be subdivided into two cases. The first is the general case, where the orientation of the plane of polarization of the incident radiation is not known. Most of the passive imaging schemes belong to this category. In the second case, it is known exactly. Many active imaging schemes fall into this category. However, since the rotating polaroid is at the receiving end in the modified PMI scheme, it is possible to know the orientation of the plane of polarization of the incident radiation by finding the position of the analyzer which corresponds to either the maximum or the minimum CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 53 intensity recorded, even in passive imaging schemes. Let us first analyze the general case and later see as to what happens when the plane of polarization of the incident light is exactly known. We denote the co-polarized intensity recorded in a general PDI scheme by Ik and the cross-polarized intensity by I⊥ . Since the polarization orientation of the incident light is not known, we assume that the intensity recorded as Ik is obtained with the analyzer at an angle φ with respect to the horizontal. Thus, the images need not correspond to the co-polarization and cross-polarization locations. The recorded intensities would then be (from eqn 4.5) Ik (n) = Iscat + Ibal sin (β) + w(n) (4.6) I⊥ (n) = Iscat + Ibal sin (β + π) + w ′ (n) (4.7) where w(n) and w ′ (n) are zero mean iid Gaussian random variables and β = (2φ + α). In PDI scheme, the estimate of the ballistic component in the recorded data is given by N Iˆbal,P DI 2
1 X Ik (n) − I⊥ (n) = N n=1 (4.8) In many implementations of PDI schemes, the scaling factor has been found to be different. It is an arbitrary constant in some schemes [53, 28]. However, using the differencing scheme, we can theoretically estimate the ballistic component by eqn 4.8. If any arbitrary scaling is used, the estimate would be a scaled version of the actual ballistic component. The scaling would affect the variance of the estimate also. Hence, we restrict our analysis to the theoretically correct scaling value. 54 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA Now, by substituting the expressions for Ik and I⊥ from eqns 4.6 and 4.7 to eqn 4.8, and by simplifying the resulting expression, we obtain, Iˆbal,P DI = Ibal sin β + w∗ (n) (4.9) where w∗ (n) is a zero-mean, Gaussian iid noise, with variance σ2 , N where σ 2 is the variance of w(n) and w ′ (n). We can easily see from eqn 4.9 that E{Iˆbal,P DI } = Ibal sin β var{Iˆbal,P DI } = (4.10) σ2 N (4.11) where E{·} stands for expectation and var{·} stands for variance of a random variable. Clearly, the estimate Iˆbal,P DI is biased, since the estimated value depends upon β. If the co-polarized and cross-polarized images do not correspond exactly to the maximum and minimum values, the estimated value of the unscattered component will not be the true value. This is a big disadvantage of PDI scheme. Only when β = π2 , we get the true estimate of the unscattered component. The estimator has a variance of σ2 , N irrespective of the value of β. Thus, only when the plane of polarization of the incident light is exactly known, we get true estimates of the ballistic component. Next, we analyze the performance of the PMI polarization intensity estimator. 4.2.1.2 Intensity imaging using PMI estimator Suppose N data values from the PMI scheme are available for analysis. It can be seen from eqn 4.5 that the 2N th M component of the N-point DFT (Discrete Fourier Transform) would contain the information of the sinusoid. Here, we assume that 2N M is an integer, i.e., N is an integral multiple of M . The N-point DFT of a signal x(n) is given by X(k) = N −1 X n=0 x(n)e−j2πkn/N k = 0, 1, 2, ......, N − 1 (4.12) CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 55 So, as per [27], where the PMI scheme has been used, the estimate of the ballistic component in the signal is given by Iˆbal,P M I = Ir  2N M  2 (4.13) where, Ir stands for the DFT sequence of Ir (n). Though, this is one way of representing the polarization data, for the sake of comparison of various schemes, we use a properly scaled definition of the DFT to get the actual estimates of Ibal . We can easily find that the estimate of Ibal is actually given by 2 Iˆbal,P M I = Ir N  2N M  (4.14) The analysis of the bias and variance of the estimator is quite formidable, and hence we have resorted to numerical simulations for the purpose. We compare the results of the simulation with that of the PDI intensity estimator, in a later section. Next, we analyze as to what best can be done to estimate Ibal , from the basic principles of estimation theory. 4.2.1.3 The holy grail for intensity imaging In the framework of estimation theory, the best estimator, with the optimality criterion being the minimum variance of the estimated quantity, is the minimum variance, unbiased (MVU) estimator [61]. However, existence of an MVU estimator does not ensure that the estimator is efficient. In addition to being unbiased, if an estimator attains the CramerRao lower bound (CRLB), the estimator will be efficient. This is the best MVU estimator one can hope to design. We first look at the possibility of finding such an optimal estimator for the problem on hand. In general, MVU estimators do not exist for all unknown parameters under all circumstances. As a first step in determining whether such an estimator exists, we see what the CRLB for the estimate is, and check if some estimator satisfies it. Here, the quantity that CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 56 we wish to estimate is Ibal of eqn 4.5. The derivation of the CRLB for various unknown parameters of eqn 4.5 when the noise is white, is given in Appendix A. It can easily be found that the CRLB for Ibal is 2σ 2 . N We now need to find if there exists an estimator that can achieve this bound. From a theorem related to the CRLB [61], it is known that an unbiased estimator for a parameter θ exists iff ∂ ln p(Ir ; θ) = f (θ)(g(Ir ) − θ) ∂θ (4.15) where p(Ir ; θ) is the probability density function (pdf) of variable Ir , parameterized by the unknown variable vector θ (see eqn A.2) and f (·) and g(·) are some functions. If such a condition is satisfied, then the MVU estimator of θ is, θ̂ = g(Ir ) (4.16) and the minimum variance is given by   var(θ̂i ) ≥ I −1 (θ) ii where I(θ) is the Fisher information matrix. For the problem on hand, θ = Ibal and Ir (n) is as given in eqn 4.5. It can be found that the condition demanded by eqn 4.15 cannot be satisfied in this case. Hence, we abandon our search for the MVU estimator of Ibal and look at other estimators that give variances closest to the CRLB. When an MVU estimator does not exist, or if it cannot be found, one usually resorts to the Maximum Likelihood Estimator (MLE), due to its interesting properties. MLE performs optimally when large enough data points are available for analysis. It also asymptotically achieves the CRLB, and hence is asymptotically efficient and optimal. Further, if an efficient estimator exists, it is achieved by the MLE. Due to all these 57 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA reasons, we resort to the MLE for estimating Ibal . The problem at hand is similar to estimating the amplitude of a single sinusoid, except for the constant term Iscat . We explore whether we can extend the analysis of the case of a single sinusoid, detailed in [61] to that of ours. By modifying eqn 4.5 as given below, it is clear that the data can be modeled linearly. Ir (n) = Iscat + Ibal cos β sin  4πn M  + Ibal sin β cos  4πn M  + w(n) (4.17) With this modification, we can express the above equation as  Ir (0)   1 0 1  

   1 cos 4π sin 4π  M M  

  = sin 8π cos 8π  1 Ir (2)  M M  .  .. .. ..  .  . .  . .       4π(N −1) −1) 1 sin 4π(N cos Ir (N − 1) M M | {z } {z |          Ir (1) Ir H     Iscat     Ibal cos β    Ibal sin β | {z Θ }         +      } | w(0) w(1) w(2) .. .           w(N − 1) {z } W (4.18) or, equivalently, by matrix notation as Ir = HΘ + W (4.19) If we can estimate Ibal cos β and Ibal sin β, we can estimate Ibal as Ibal q = + (Ibal cos β)2 + (Ibal sin β)2 (4.20) The reason for choosing such a linear form to model the data is intentional, since the linear model gives immense advantage in designing the estimator. It has been proved that [61] if the observed data X are described by the general linear model X = HΘ + W (4.21) 58 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA where H is a known N × p matrix with N > p and of rank p, Θ is a p × 1 parameter vector to be estimated, and W is a noise vector with pdf N (0, C), then the Maximum Likelihood Estimator of Θ is
−1 T −1 Θ̂ = HT C−1 H H C X (4.22) Θ̂ is also an efficient estimator, in that it attains the CRLB and hence is the MVU estimator. The pdf of Θ̂ is Θ̂ ∼ N (Θ, (HT C−1 H)−1 ) (4.23) Since the linear model of eqn 4.19 satisfies the above conditions, we have an efficient estimator for Ibal cos β and Ibal sin β and thus for Ibal , which is the MLE of Ibal . Although the estimator is biased if data points are few; it is asymptotically unbiased [62]. Moreover, since the estimate is obtained by a non-linear transformation of the MVU estimates, the estimator cannot be efficient [61]. For the problem on hand, the rank of the matrix H can be shown to be 3 by considering just the first 3 rows and performing row reduction on the 3 × 3 matrix obtained. The noise vector has the covariance matrix σ 2 I, where I is an N × N unit matrix, and σ 2 is the unknown noise variance of the series. To obtain the MVU estimates of Ibal cos β and Ibal sin β easily, we need a simple form for the inversion of the matrix H. It can be shown that this happens if N is an integral multiple of M and then, the MLE estimate Θ̂ for various components turns out to be N −1 1 X Ir (n) Iˆscat,M V U = N n=0   N −1 2 X 4πn ˆ Ir (n) sin Ibal cos βM V U = N n=0 M (4.24) (4.25) 59 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA Iˆbal sin βM V U   N −1 2 X 4πn = Ir (n) cos N n=0 M (4.26) where, M is the periodicity (number of steps per rotation) of the rotating polaroid. The covariance matrix of the estimate Θ̂ is,  σ2 N    0  0 0 2σ 2 N 0 0    0   (4.27) 2σ 2 N Hence, we have var{Iˆscat,M V U } = σ2 N 2σ 2 var{Iˆbal cos βM V U } = var{Iˆbal sin βM V U } = N (4.28) (4.29) substituting eqns 4.25 and 4.26 into eqn 4.20, we obtain Iˆbal,M LE v u(N −1  )2 (N )2  −1 X X 4πn 4πn 2u Ir (n) cos + Ir (n) sin = t N M M n=0 n=0 (4.30) It turns out that Iˆbal,M LE is a random variable with the distribution being that of the
N square root of a gamma random variable with density Γ 1, 4σ 2 , when finite samples are used. However, due to the properties of the MLE [61], the density asymptotically tends 2 to N (Ibal , 2σ ). N It is clear that the right hand side of the above equation is the same as that of eqn 4.14. Hence, we arrive at the important result that the MLE estimate of Ibal can be obtained by the PMI scheme. It is worth observing that if β = π2 , we can obtain the MVU estimates of Ibal . Hence, if we know the exact phase relations, we can obtain the MVU estimates of Ibal , using the PMI scheme. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 60 We end our search for better estimators of Ibal here, since the MLE estimator almost always does the best job, when the MVU estimator does not exist. In all subsequent references to PMI, we imply the usage of the MLE estimator given in eqn 4.19, followed by the transformation given in eqn 4.20. This also corresponds to the LSEK (Least Squared Estimator of K sinusoids) estimator mentioned in [62], as against the DFT, mentioned as LSE1 in [62]. The advantages of LSEK over LSE1 have been mentioned in the same article, and apply to this problem also, since we are estimating both the DC component and the sinusoidal amplitude simultaneously. 4.2.1.4 Comparison of estimators The following observations can be made from our previous discussions. • The PDI estimator is biased due to its dependence on the value of β. • Though the MLE (and hence the PMI estimator) is biased when only few data points are available for analysis, it is asymptotically unbiased and efficient. • The variance of the PDI estimator (eqn 4.11) is always less than that of the PMI estimator. In fact, it is less than CRLB given in A.8. This is at the cost of the bias of the estimator, which does not improve with the amount of data available for analysis, unlike the PMI estimator. • Unbiased estimates of Ibal can be obtained if β = π2 , using both the PMI and the PDI estimators. • If we can choose β = π , 2 the PDI estimator ensures better performance than the PMI estimator, due to lesser variance. • If β cannot be known a priori, the PMI estimator is preferable over the PDI estimator. The exact phase relations can usually be known in the case of active imaging. Hence, for active imaging, PDI schemes are more useful than the PMI schemes. For passive imaging, though the PMI scheme seems to be more suitable, there are certain circumstances where CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 61 the PDI scheme may still perform better. There are applications where the parameter of interest is not the exact value of the sinusoidal amplitude, but its relative value across the scene. Since the PDI scheme gives uniformly scaled values of the sinusoidal amplitude across the scene, it may be better to use the PDI scheme since its variance is lower than the PMI scheme. At the end of this section, we present results from numerical simulations to substantiate our analysis. With this, we conclude our analysis of the various estimators for Ibal in white noise, and proceed to study some estimators for estimating the same quantity in coloured noise. 4.2.2 Polarization intensity imaging - the case of coloured noise For the theoretical analysis of the various estimators in coloured noise, we need to be specific about the kind of noise being considered, though, the estimators themselves need not have a priori knowledge of the noise characteristics. It has been found from experimental data that the noise process can be adequately modeled by an Auto Regressive (AR) process of order 1. We reproduce here, the basic equation that governs the behaviour of the data samples. Ir (n) = Iscat + Ibal sin  4πn +β M  + v(n) n = 0, 1, 2, · · · , N − 1 (4.31) The analysis of the PDI and PMI schemes, in the case of white noise gave us closed form expressions for the mean and the variance of the estimators. However, to be able to do the same in the case of coloured noise is a formidable task. Moreover, since the theoretical analysis of other estimators discussed in this section also throw up similar challenges, we have resorted to Monte-Carlo simulations to compare the estimators. First, we make some observations about the PDI and the PMI estimators, following which, we see if an MVU estimator or an MLE exists for the problem on hand, and from there, we proceed to see some estimators apart from the ones we have discussed so far, and finally compare all the estimators discussed. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 62 We know that the assumed co-polarized and cross-polarized intensities in PDI schemes are given by Ik (n) = Iscat + Ibal sin (β) + v(n) (4.32) I⊥ (n) = Iscat + Ibal sin (β + π) + v ′ (n) (4.33) In the case of white noise, the v(n) and v ′ (n) terms are independent, and hence we could easily proceed to find the mean and the variance of the estimators. However, in the case of coloured noise, these two are dependent. Further, the v(n), n = 0, 1, 2, · · · , N − 1 terms are themselves correlated and the v ′ (n) terms also show a similar behaviour. Hence, it is clear that the PDI scheme performs poorly in coloured noise, as compared to its performance in white noise. This observation is, in general, true for all estimators that we consider. The LSEK estimator, used in the PMI scheme, is also suboptimal, since, it too considers the noise to be uncorrelated [62]. Hence, we now set out to find the best possible estimator of Ibal in coloured noise, from the point of view of estimation theory. Due to the same reasons as in the case of white noise, it has been found that we can not obtain an MVU estimate of Ibal , and we have to be satisfied with the MLE estimate given by eqn 4.22. However, what makes the coloured noise case different is that, the noise covariance matrix C in this case will not be σ 2 I, and hence, the estimates depend explicitly on the noise covariance matrix terms. On the other hand, if the noise is white, we can obtain the MLE estimates of the amplitude, without having to know the noise variance per se. In the case of coloured noise, the noise covariance matrix and needs to be obtained from the observed data. This affects the performance of the MLE . Moreover, we do not directly have the noise samples, since the data contains signal plus noise information. This makes the estimation of the noise covariance matrix difficult. Hence, numerous estimators have been designed to tackle this problem [62]. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 63 We can use an intuitive, nonetheless, a very effective technique to estimate the noise samples from the data containing both signal and noise. These noise only samples can then be used to obtain the noise covariance matrix. The technique we have used has been found to improve the performance of most of the high performance estimators given in [62]. We compare them at the end of this section. Now, we see the technique of estimating the noise only samples. In the PMI imaging setup, we know the periodicity M of the sinusoidal component a priori, since the rotating polaroid is at the receiving end of the imaging setup. We use this information to achieve our objective of obtaining the noise samples alone. In the PMI scheme, N , the number of data points available for analysis, is chosen to be an integral multiple of M . By taking a N −point DFT of the observed data, the information about the sinusoid will be localized to a single DFT coefficient. Here, we refer to both positive
th and negative frequency terms, when we say a single DFT coefficient. It is the 2N M coefficient and its symmetric counterpart, that carry this information. Similarly, the I scat information is localized to the DC component. The rest of the frequency components contain only the noise information. We make an assumption that, by reducing the two DFT coefficients containing the signal information (Iscat and Ibal ) to zero, and taking an inverse DFT, we obtain noise information alone. In doing so, though the signal information is eliminated from the reconstructed samples, some noise information is also lost in the process and changes in phase will affect the reconstructed signal. However, the idea is that, an estimate of the noise covariance matrix, obtained with this reconstructed noise-only samples, would be better than that obtained as given in [62]. Once we obtain the noise samples in this manner, we can estimate the noise covariance matrix. As we saw earlier, since we usually cannot know the noise covariance matrix a priori, there have been various techniques developed to estimate the amplitude of a sinusoid in coloured noise, some of which have been reviewed and compared in [62]. Of all the CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 64 matched-filter based amplitude estimation techniques mentioned in [62], for our application, the APES (Amplitude and Phase Estimator for Sinusoids) algorithm seems to be the best suited, since we need to estimate the amplitude of one sinusoid only. The algorithm has been found to be asymptotically statistically efficient. However, what is more interesting is that, it is reported to be unbiased even when only a finite number of data points are available for estimating the amplitude of a complex sinusoid. But, since we are interested in estimating the magnitude of the sinusoid, the efficiency will be lost. Still, the performance of the APES estimator for small data lengths is what makes it attractive for our application. The design of the APES estimator has been detailed in [62, 63], and its efficient implementation has been described in [64]. Our modification of the APES estimator, using the noise covariance matrix obtained as above, shows lower variance in estimating the amplitudes, especially at lower SNR. We call it as the APESR estimator, where the appended R stands for ‘robust’. The design of this estimator is briefly described below. Let x(n), n = 0, 1, 2, · · · , N −1 denote the observed samples containing both the signal and noise information. We obtain the noise only samples z(n) from x(n), as explained earlier. We form overlapping sub-vectors of z(n), denoted by q(l) as followsq(l) = [z(l) z(l + 1) z(l + 2) · · · z(l + R − 1)]T l = 0, 1, 2, · · · L − 1. (4.34) where L = N − R + 1. A suitable choice of R for the APES estimator has been reported to be N 2 in [63], and we retain the same for this modified estimator too. The estimate of the noise covariance matrix Q̂ is obtained as L−1 Q̂ = 1X z(l)z T (l) L l=0 (4.35) This Q̂ is used to estimate the amplitude of the sinusoid. However, the estimated amplitude will be a complex number, whose absolute value corresponds to the magnitude CHAPTER 4. PROCESSING POLARIZATION-RICH DATA of the sinusoid (Ibal ), and whose angle will be equal to π 2 65 − β. Next, we discuss a couple of techniques that may improve the performance of the various estimators discussed so far. Later, we compare the performance of all the estimators. 4.2.3 Chunked data processing and bootstrapping A general technique that can be applied for all the estimators discussed so far, for both white and coloured noise, is chunked data processing. Similarly, bootstrapping is another general technique that can be used to improve the performance of estimators. In chunked data processing, instead of using all the data points to obtain a single estimate, we break the data into smaller chunks, from which we obtain multiple estimates for the same variable and average them. However, there is no guarantee that the resultant estimate obtained would be better than what can be obtained by considering all the data points in a single go. But, as we shall see, the APES estimators perform significantly better with chunked data. The PMI, PDI and the MLE estimators do not show such improvements. Moreover, chunking can be extended to overlapping windows of data, leading to a large number of relatively high variance estimates, which, upon averaging, can yield results better than the single estimate obtained from the available data. Though it looks as if the chunked data processing technique requires huge computation, surprisingly, it need not. The processing of chunked data when pushed to its limit (i.e., when we consider chunks which have all but a single data point in common), can be viewed as processing data using a small buffer. As a new sample comes in, the oldest sample leaves the buffer, thereby leading to a new realization of the size of the data buffer. A new estimate is found for every single data point in this manner, and all the values averaged. If the buffer size chosen is small, the processing is still simple, though it needs to be repeated many times. Moreover, the memory requirement will also be less, since we need not retain all the data values to estimate the parameters. This is an advantage of chunked data processing. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 66 In bootstrapping, one generates more data samples from the existing ones. As the number of samples available for analysis increases, the noise variance comes down, and hence, the performance of the estimator improves. However, the computational costs involved are much higher in this case, as compared to chunking. A good introduction to bootstrapping is given in [65]. For our problem, we can assume the data to be emanating from M/2 different sources, where each source represents a point in half a period of the rotating polaroid. If the noise is known to be white, we can form extended data by sampling with or without replacement, the data points corresponding to these positions, and concatenating the data. However, the extension to the case of coloured noise is not as simple. An application of bootstrapping for detection of a sinusoid in coloured noise is given in [66], where different types of bootstrapping algorithms and ways of implementing them in coloured noise are given. We used the block bootstrapped algorithm. First, we estimate the noise only part of the data as mentioned earlier, instead of the time domain approach given in [66]. Then, we bootstrap the noise samples to obtain better estimates of the noise covariance matrix, to be used in the MLE and the APESR estimators. It would also be interesting to know as to whether we can further improve the estimation performance by processing chunked bootstrapped data. We carried out such studies too, and the results are reported later, when comparing the different estimators. Next, we study the computational complexity of the various estimators, since some applications demand real-time processing, whereas for others, offline processing suffices. 4.2.4 Computational complexity of the estimators A basic PDI scheme with modest memory requirements can be implemented as follows. The effect of adding co-polarized or cross-polarized images can be obtained by increasing the integration time of the CCD to a long time. Then, by subtracting such a single cross-polarized image from a co-polarized image, we can obtain PDI results. Such a scheme can be implemented with the memory requirement being as little as twice the CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 67 size of co-polarized and cross-polarized images. In case long integration times cannot be allowed, the scheme can be parallelized by using multiple sensors. However, the memory requirements increase proportionally with the number of sensors. The actual processing performed by a signal processor involves only subtraction and addition of integers, and a final scaling of the values, to be rendered as images. This can be carried out at extremely fast rates, and hence the PDI scheme can easily be implemented in real time by modern signal processors. The PMI scheme has higher computational needs than the PDI scheme. The core task in PMI processing is the computation of a single DFT coefficient (or LSEK estimate), for which there are efficient methods. The problem of computing only a few DFT coefficients has been addressed by the Goertzel algorithm, the FFT (Fast Fourier Transform) pruning technique and other techniques, which have been compared in [67]. The memory requirements for implementing these techniques are the same as in the case of PDI schemes. The implementation essentially needs successive multiply-add-accumulate steps, followed by a square-root operation. Hence, the processing is a bit more involved, than in the case of PDI schemes. The PMI scheme also lends itself to parallelization and can be implemented in real time. A direct implementation of the Goertzel algorithm needs memory at least twice the size of the co-polarized or cross-polarized images. The MLE and the two APES estimators involve more computation than the PDI and PMI schemes. The bottleneck in these algorithms is the inversion of matrices. In that aspect, the MLE scheme needs more number crunching than the APES schemes, since it has to invert a matrix at least twice the size of what the APES estimators have to. An efficient implementation of the APES estimators is given in [68], along with the relative computational costs involving FFT based algorithms. Among the two APES estimators, the APESR estimator needs less computation, and hence should be preferred wherever the performances of both the estimators are similar. The MLE and the APES estimators also have the disadvantage that they need to have the whole data to begin processing, unlike the PDI and PMI schemes, which can compute in-place, sequentially, as the data arrives. Hence, the memory requirements of the former estimators are higher than that 68 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA of the latter. However, if the noise covariance matrix can be known a priori, the APES and the MLE methods estimate in real time. As mentioned in the previous subsection, the bootstrapping algorithms need considerably higher computational efforts than the normal algorithms, but chunked data processing can be simpler than the normal processing schemes. Since there are various parameters involved to compare the computational costs, like the chunk length, the overlap length and the length of the bootstrapped data, we do not get into a detailed analysis of the computational costs involved in bootstrapping and chunking. Next, we compare the performance of the various estimators based on Monte-Carlo simulations and also based on the performance of the estimators in real noisy data. 4.2.5 Comparison of polarization intensity estimators While comparing the performance of various estimators, we first give the results of Monte-Carlo simulations. We give only those results on real data, that do not match the predictions of Monte-Carlo simulations. 4.2.5.1 Comparisons based on Monte-Carlo simulations All simulations reported in this thesis have used 200 realizations of 64 data points generated according to eqn 4.31, to test the PMI, MLE, APES, and the APESR estimators. For testing the PDI estimator, we used 32 data points each with β = π 2 and 3π , 2 corresponding to the co-polarized and the cross-polarized data, respectively. We have also tested the performance of all the estimators on chunked data. We used a chunk length of 32 data points, and an overlap of 31 data points between consecutive windows, thereby obtaining 32 estimates, which were averaged to obtain the final result. We have analyzed the effect of bootstrapping on all the estimators by constructing an extended sequence of 128 data points from the given 64 data points. In this case, the dimension of the noise covariance matrix used in MLE and APESR estimators changed from 32 × 32 to 64 × 64. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 69 All the data generated according to eqn 4.31 assume the signal frequency f to be 0.125. The values of Iscat and Ibal are chosen to be 4 and 1, respectively, so that the degree of linear polarization (DOLP) corresponds to 0.25. As we shall see, the performance of the estimators does not depend on the values of Iscat or Ibal , but on the noise characteristics. Hence, even if we choose any other value of Iscat and Ibal , the bias and variance of the estimators will remain the same. Various characteristics of the estimators are studied by varying the SNR through the variance of the iid Gaussian random variable governing the noise process. For deciding the range of SNR values to be considered, we calculated the average SNR of a 10 × 20 region of eight data sets that we used. The region chosen was outside the region of the geometric shadow of the hidden object, which ensured that we obtain the ballistic sinusoidal component apart from noise. Such an analysis showed the actual SNR to vary from around -14 dB to 9 dB. Hence, we have considered the SNR range from -25 dB to +25 dB for our analysis. For calculating the SNR, we used the average of the ratios of the signal energy to noise energy. Since we knew the frequency of the signal, we could localize it to one frequency component and consider the rest of the components to contain noise information. We used zero mean data for calculating the overall energy. For the case of white noise with variance σ 2 , the SNR is given by SNR in dB = 10 log10  2 Ibal 2σ 2  (4.36) For coloured noise, we have used the local SNR as the varying parameter, to test the performance of the estimators. The local SNR is defined as [62] SNR in dB = 10 log10 2 N Ibal φ(f ) where, φ(f ) is the power spectral density of the observation noise v(n) [69] φ(f ) = |1+ σ2 , −j2πf k |2 k=1 ak e Pp CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 70 σ 2 is the variance of the white noise process, which drives the AR process, and ak are the AR coefficients governing the noise process. Though we found that the noise process in the actual data could be adequately represented by an AR1 process, we have tested the algorithms with AR2 noise process too, so that the performance of the estimators in unknown noise conditions could be better understood. For studying the performance of the estimators under AR1 noise, we chose the AR coefficient a, to be 0.50, in all the simulations. To study the performance of the estimators in AR2 noise, we chose the AR coefficients to be a1 = 0.50 and a2 = −0.125, respectively. The pole frequency for this choice of AR coefficients corresponds to a discrete frequency of 0.125, and coincides with the frequency of the sinusoid, thus creating a relatively difficult situation to estimate the sinusoidal amplitude. In our comparison of estimators, we emphasize on simulations rather than actual data, since the actual data pertain to only two experiments. However, the noise models give a generalization to the data that can be obtained in polarization based imaging experiments. 4.2.5.2 Comparisons based on performance in real noise To compare the performance of the estimators in real noise, we used data from polarization modulation imaging experiments conducted with the setup shown in Fig. 4.1. The details of the experimental setup are given below. • Source: Linearly polarized He-Ne Laser. Wavelength= 632.8 nm. Power= 10 mW. Beam Dimension (1/e2 ): 0.65 mm. • Beam Expansion Optics: It consisted of two plano-convex lenses of focal lengths 5 cm and 30 cm, giving a beam expansion factor of 6, with which we could get a beam diameter of nearly 7 mm. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 71 • Spatial Filter: The spatial filter consisted of a convex lens with a focal length of 10 cm followed by an iris with a very narrow opening placed at the focal point of the lens. The diameter of the iris was not more than 1 mm, but we could not measure it precisely. • Charge Coupled Device (CCD): EEVTM , the intensified CCD used has 8 bit resolution and gives black and white images. The gain could be varied over a large range. We grabbed 512 raw images of size 240 × 320 in each of the experiments. For all the experiments, the frequency of the rotating polaroid was 0.125. i.e., the time series at each pixel location consisted of 512 data points with frequency of Ibal being 0.125. • Hidden Object: The object used for imaging purposes was an opaque cross of nearly 1 mm thickness. For analyzing the performance of the estimators, we used data obtained from two experiments, one with the scattering medium consisting of polystyrene microspheres of diameter 2.97µ dispersed in water, and the other with the scattering medium being mist. We call the former medium as ‘medium 1’, and the latter as ‘medium 2’. The noise data generated from data recorded with these two media will henceforth be called ‘noise 1’ and ‘noise 2’, respectively. A glance of the power spectra at a few pixel locations of each set showed that the noise process corresponding to medium 1 was nearly white, while that of medium 2 was coloured. This was endorsed by the average AR coefficients calculated from 200 pixel locations, using the Levinson-Durbin recursion algorithm. Noise 1 had the first two AR coefficients as -0.0021 and +0.0183, while noise 2 showed the coefficients as +0.1561 and +0.0853, thereby supporting our observations. The pixel locations chosen for our analysis were the same as that used for calculating the average SNR of real data. In all the polarization imaging experiments conducted, the exact orientation of the incident linearly polarized light was not known and hence the sinusoidal phase was also unknown. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 72 In order to compare the performance of the estimators in real noise, we need to have control over the values of Ibal and Iscat , which cannot be easily attained in actual experiments. None of our data showed uniform distribution of Ibal and Iscat . Thus, we had to find a way of comparing the performance of the estimators. Towards that end, we first obtained 512 noise only samples from 200 pixel locations, by the method described in section 4.2.2. Out of these 512 samples, we assumed the first 64 to represent the actual noise samples. To these noise samples, we added a sinusoid of frequency 0.25, with varying amplitudes, to simulate different SNR conditions from -25 to +25 dB. The frequency (0.25) was chosen to avoid the possible changes in noise conditions that would have taken place at a frequency of 0.125, during the process of estimating the noise samples. The amplitude of the sinusoid added was chosen using the definition of SNR for white noise given in eqn 4.36 . The value of Iscat was chosen to be four times the value of Ibal , thereby giving a DOLP of 0.25. The 200 pixel locations chosen to obtain the noise only samples were the same as those from which we obtained the average SNR and average AR coefficients, that we mentioned earlier. A few points which distinguish the analysis of the estimators in real noise from the analysis in simulated noise are mentioned below. • The samples representing noise are not exactly noise samples, but were obtained from signal+noise observations. Still, we call this noise as ‘real noise’, to distinguish it from simulated noise. • For polarization magnitude and degree of polarization estimators, the variance and bias of the estimators depend on the noise properties only, and not the SNR, as is evident from the calculations of CRLB given in Appendix A. Hence, the bias and variance curves of some of the estimators remain nearly constant for all values of SNR. This does not imply that there is no improvement in the performance of the estimators. It has to be interpreted bearing in mind, the changing amplitude of the sinusoid with SNR. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 73 • Since the noise properties are unknown, we have not plotted the CRLB in our discussion on real data. 4.2.5.3 Comparison of estimators within and across classes We compare the performance of 5 different classes of estimators, namely, the PDI, PMI, the MLE and the two APES estimators. Within each class, we have variants which process chunked data, bootstrapped data and chunked bootstrapped data (estimators which first bootstrap the data and then chunk the bootstrapped data to estimate the parameters), leading to 22 estimators in all. Due to the paucity of space, we choose the best estimators from each class and compare only these five estimators. For convenience, we use the following notations to distinguish the different estimators in each class. The normal estimators are denoted by the estimator class name followed by ‘N’. Suffix ‘B’ stands for the estimator which uses the bootstrapped data. An estimator class name followed by ‘C’ denotes an estimator which performs analysis on chunked data. Similarly, the suffix ‘CB’ stands for an estimator which analyzes bootstrapped and chunked data. In the case of MLE and the APESR estimators, we implemented two different algorithms for bootstrapping. The first assumes the noise to be white and uses sampling with replacement to extend the data. The other assumes the noise to be coloured, and uses block based bootstrapping algorithm to find better estimates of the noise covariance matrix. However, we found that the block based bootstrapping algorithms did not show reliable performance, and hence we do not report those results here. The performance of all the estimators were similar in AR1 and AR2 noise conditions. Hence, we do not report our observations on the performance of the estimators in AR2 noise. In all the comparisons to follow, when we refer to all noise conditions, it is inclusive of the real noise too. We now present our observations on the performance of the various polarization intensity estimators of each class, under varying SNR for both simulated and real data. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 74 • For simulated data, all the estimators were found to be positively biased under coloured noise. The bias and variance of all the estimators dramatically increased in coloured noise as compared to white noise, indicating poorer performance of the estimators in coloured noise. • Among the variants of the PDI estimator, the PDI N showed the least bias and variance under all noise conditions. It was found to be positively biased in white noise, up to -15 dB, beyond which, it was found to be unbiased. It was also found to be positively biased at all SNR in real noise, like other PDI variants. We attribute this to the unknown phase of the sinusoid in experimental data. • Among the PMI estimators, the PMI N showed the least bias and variance under all noise conditions. The PMI N estimator was also found to be positively biased up to around -5 dB in white noise conditions, beyond which, it was found to be unbiased. In real noise, the bias of PMI N estimator was found to be higher than PMI C beyond -10 dB. • The MLE N showed better performance up to around 10 dB, beyond which, MLE C and MLE CB showed better performance in white noise. Under coloured and real noise conditions, MLE N showed the best performance at all SNR. • Among the APES estimators, the APES CB showed the best performance under all noise conditions. However, it was found to be positively biased even under white noise. For real noise conditions, APES C showed better performance than APES CB. However, the bias of APES C was a bit higher than APES N for noise 2. • Among the APESR estimators, the APESR N showed better performance than other estimators up to -10 dB in the case of white noise and up to 15 dB in the case of coloured noise, beyond which, the APESR CB was better. Hence, we choose the APESR N to be the representative of this class, since we are interested in performance at low SNR. We could not observe any clear trend in real noise. APES CB performed better than all others in noise 1, whereas, APES C performed better in noise 2. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 75 Thus, the best estimators of polarization intensity information under varying noise conditions for simulated data were PDI N, PMI N, MLE N, APES CB and APESR N. Figs. 4.2(a) and 4.2(b) show the bias and the variance of different estimators in white noise. The CRLB of Fig. 4.2(b) was calculated using eqn A.7. Figs. 4.3(a) and 4.3(b) show the bias and the performance of the estimators in AR1 noise. For this case, CRLB was calculated using eqn B.17. The best estimators of polarization intensity information from real data were PDI N, PMI N and MLE N for the first three classes. Though there was no consistency in performance among the APES and the APESR estimators, we choose APES C and APESR CB as representatives of their classes, since they showed better performance with noise 1 and noise 2, respectively. Fig. 4.4 shows the bias and variance of these estimators in real noise. We have plotted the results corresponding to that of noise 1 only. The same hold true even for results with noise 2. From the results of numerical simulations, we observe the following • The performance of the estimators is relatively independent of the noise characteristics. The PDI and the PMI estimators consistently perform better than others. • The PDI N estimator is unbiased and has the least variance among all the estimators for any noise, when β = π2 . Hence, the PDI estimator should be the preferred choice for PII, unless accurate estimates of Ibal are needed when the phase conditions are not known. The variance of the PDI estimators was lower than that of the CRLB under all noise conditions. This should not come as a surprise, since the reduction in variance is at the cost of the bias of the estimator, as we find in the case of real data (see Fig. 4.4(a)). Only when β = π/2, we get unbiased estimates of polarization intensity. Otherwise, the bias will not decrease even with increasing amount of data available for analysis. • The PMI N estimator is the next best, and has the added advantage of being asymptotically unbiased, unlike the PDI estimator. Hence, when exact phase relations are CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of PII estimators in white noise. Only the best performer in each class of estimators, namely, PDI, PMI, MLE, APES and APESR are shown. (b) Variance of PII estimators in white noise. Figure 4.2: Performance of PII estimators in white noise. 76 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of PII estimators in AR1 noise. (b) Variance of PII estimators in AR1 noise. Figure 4.3: Performance of PII estimators in AR1 noise. 77 78 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of PII estimators in real noise. 9 8 PDI N PMI N MLE N APES C APESR CB 7 Variance 6 5 4 3 2 1 0 −25 −20 −15 −10 −5 0 5 10 15 20 SNR dB (b) Variance of PII estimators in real noise. Figure 4.4: Performance of PII estimators in real noise. 25 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 79 not known, as in the case of passive imaging, and accurate estimates of Ibal are needed, the PMI estimator should be the preferred choice. • Though MLE N and other estimators perform on par with the PMI N estimator under varying noise conditions including real noise, their high computational costs favour using the PDI N or the PMI N estimators. • The APES estimators are hardly useful for PII, when the frequency of the sinusoid is exactly known. We performed similar analysis for data of lengths 32 and 128. The chunk lengths considered for the two cases were 16 and 64, respectively. The overlap was 15 for the 32-point data, and 63 for the 128-point data sets. The bootstrap lengths considered were 64 and 256, respectively. The behaviour of the estimators remained the same under all these conditions, thereby generalizing our results to varying data lengths for simulated data. For real noise 1, it was found that the behaviour of the estimators did not vary with the number of data points available for analysis. However, in case of noise 2, with increase in the number of data points available for analysis from 64 to 128, chunking improved the performance of all the estimators. Hence, PDI C, PMI C, and MLE C showed better overall performance compared to their normal counterparts. However, with N = 32, the relative performance of the estimators remained the same. Thus, we conclude that data chunking may come to help when larger number of data points are available for analysis. 4.2.5.4 A weak case for the APES estimator As we observed, the APES estimators do not seem to be much use for PII. However, as we mentioned earlier, the behaviour of the PDI N and PMI N estimators, as reported earlier, can be observed only under certain conditions. It is important to know as to how the performance of these estimators change, when these conditions are not satisfied. The analysis of the APES estimators vis-a-vis the PMI (LSEK) estimator, which we have used here, has been given in [62]. The behaviour of PDI N when the phase relationships are CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 80 not known has already been discussed. Now, let us analyze the behaviour of PMI N, when N is not an integral multiple of M . The same analogy holds true for MLE N also, since, its performance also depends on the same criterion. The first important result is that the estimates Iˆbal cos β and Iˆbal sin β will no longer be MVU estimates [61]. This may perturb the PMI and MLE estimates by large amounts from the actual values. After obtaining the data, we may find the frequency of the signal to be different from the expected value. Such a case may happen when the frequency of the rotating polaroid is not exactly known. In such cases, we may have to estimate the frequency first and then perform amplitude estimation. Then, all the useful properties of the APES estimators listed in [62] come to fore. However, we do not get into a detailed analysis of the case, since more often than not, the frequency of the rotating polaroid can be tightly controlled. We conducted simulations to understand the behaviour of the PMI N, MLE N and the APES N estimators when the frequency of the sinusoid was chosen to be 0.115, instead of 0.125. However, with this small change in frequency, we could not find a significant difference in the bias of the estimators at low SNR. We could observe that the APES N estimator performed better only at SNR above 25 dB. 4.2.5.5 A case for the chunked PMI estimator It is known that coarser DFT will pack energy of wider frequency ranges into the DFT coefficients. This follows from Parseval’s theorem, which has to do with conservation of energy under Fourier transformation. Thus, we studied the possibility of using PMI C to improve the performance of the PMI estimator, so that we can achieve lower bias, than what is possible with the PMI N estimator, when the frequency is slightly off from the expected value. Towards that end, we performed simulations wherein the frequency of Ibal was chosen to be 0.115, but still, we assumed the frequency to be 0.125, for analysis. A situation of this sort can occur due to unknown errors in the frequency of rotation of polaroid. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 81 The simulations considered chunks of length 8, 16 and 32. The overlap in all the cases was one less than chunk length. The variance of the chunked PMI estimators were all nearly the same and were very close to that of the PMI N estimator, which are shown in Figs. 4.2(b) and 4.3(b). Here again, we report the performance of the estimators in white and AR1 noise only. Fig. 4.5 shows the bias of the PMI N and its chunked variants when f = 0.115 under white and AR1 noise conditions. It can be observed that the bias in the estimated amplitude varies widely with the chunking length. Clearly, for the case of white noise, at higher SNR, the bias decreases, as the chunk length decreases, validating our hypothesis. However, for the case of coloured noise, the bias is high for lower chunk lengths at low SNR. So, chunked PMI is useful only when the SNR is high, and more so, when the noise is white. However, we are not in a position to comment about the optimum chunk length for a given frequency. It can just be commented that chunking may improve the performance of the PMI estimator, when there are small uncertainties in frequency, around the expected value. 4.2.5.6 Effect of chunking parameters on PII APES estimator Though we found the APES estimator to be not all that useful for intensity imaging, we still report the effect of chunking parameters on their performance. We observed that processing chunked data involves two parameters, namely, the chunk length and the overlap between chunks. The way these two parameters affect the estimation procedure is important, since we observed that the APES CB estimators can perform considerably better under certain conditions. We have chosen the APES CB estimator to study the effect of these parameters. We considered chunks of length 8, 16, 32 and 56, on 200 sets of 64 data points. For each chunk length, the overlap was chosen to be one less than the chunk length. All the other parameters of the simulations were similar to that used in the comparison of estimators. We assumed AR1 noise for these studies and Fig. 4.6 shows the result. It can be observed that both bias and variance of the APES estimator are reduced, as the chunk length decreases. Thus chunking can improve the performance of the APES estimator CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of PII PMI N and variants of PMI C in white noise. (b) Bias of PII PMI N and variants of PMI C in AR1 noise. Figure 4.5: Performance of PII PMI N and its chunked variants, when N is not an integral multiple of M . For this example, f = 0.115 and N = 64. 82 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 83 dramatically. However, the effects of varying the chunk length could be different for different estimators, and have to be studied on a case by case basis. Next, we report our findings on the effect of varying overlap on chunked estimation. To understand the effect of overlap, we conducted simulations with the chunk length fixed at 32, and the overlap taking values of 0, 8, 16, 24 and 31. Figs. 4.7 and 4.8 report the results of analysis for white and AR1 noise, respectively. As can be observed, in white noise, there is a clear trade-off between bias and the variance of the estimators as the overlap varies. The overlap length that gives the minimum variance shows the highest bias and vice-versa. In coloured noise, since the bias goes from positive to negative at around 15 dB, overlapping helps only at low SNR. At high SNR, overlapping gives negatively biased results. Hence, depending on the importance of the criterion, one can choose to resort to overlapping or non-overlapping technique. However, the way chunking and overlapping work in conjunction needs to be studied in detail, to understand and predict the behaviour of the estimators. Moreover, the behaviour may change with the estimators themselves, and hence, each estimator needs to be studied independently. 4.2.5.7 Effect of bootstrap length on APES estimator We studied the performance of the APES B estimator with increasing bootstrap length. For the bootstrapping algorithm, the frequency of the signal must be precisely known, since the periodicity of the signal is taken in to consideration. We chose bootstrap lengths of 96, 128, 256 and 512, which were all generated by extending 64 observed data points, assuming the noise to be white. Hence, we report the behaviour of APES B in white noise only. All the other parameters of the simulations were the same, as in previous cases. Fig. 4.9 shows the performance of the APES B estimator and its bootstrapped variant in white noise. Increased bootstrapping reduces the variance at all SNR. Bootstrapping reduces the bias of the estimator at low SNR, but at high SNR, a keen observation of the results shows that bootstrapping induces slightly higher bias than the APES N (not CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of PII APES N and variants of APES C in AR1 noise. (b) Variance of PII APES N and variants of APES C in AR1 noise. Figure 4.6: Effect of chunk length on PII APES C estimator in AR1 noise. 84 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 85 (a) Bias of PII APES N and variants of APEC C estimators in white noise. (b) Variance of PII APES N and variants of APES C estimators in white noise. Figure 4.7: Effect of overlap length on the performance of PII APES estimator in white noise. For all the cases, the chunk length chosen was 32. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 86 (a) Bias of PII APES N and variants of APES C estimators in AR1 noise. (b) Variance of PII APES N and variants of APES C estimators in AR1 noise. Figure 4.8: Effect of overlap length on the performance of PII APES estimator in AR1 noise. For all the cases, the chunk length chosen was 32. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 87 evident in the plot). The effect of bootstrapping under coloured noise also showed similar behaviour, albeit all the estimators were positively biased. We conclude our discussion on PII by commenting on other issues involved. 4.2.5.8 Other issues related to polarization intensity imaging Obtaining noise only samples In order to study the noise characteristics, it is important to have noise only realizations. Such data can be obtained by grabbing images by keeping the rotating polaroid stationary at a fixed location. An estimate of noise can be obtained by removing the dc value from each of the time series. This can be used to estimate the noise covariance matrix, which can be directly used in MLE, thereby reducing the computations considerably. Choice of f In PII, we have the option of choosing the frequency of rotation of the polaroid. Throughout our discussion, we have assumed the signal to be at a frequency of 0.125. Here are a few points as to how the frequency affects PII: • The choice of f will not make a difference in case the noise is white. • If the noise is a low pass AR1 process, we stand to benefit by choosing a higher f , since the lower noise energy at higher frequencies leads to improved performance of the estimators. Similarly, choosing a lower f would be helpful if the noise is a high-pass AR1 process. However, choosing a very high or very low frequency will affect polarization orientation imaging, which we are yet to discuss. • In case the underlying noise process is of a higher order than AR1, it is judicious to avoid the pole frequencies while choosing f , in order to obtain better estimates. For this, we need to know the properties of noise a priori, which can be obtained as discussed earlier. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 88 (a) Bias of PII APES N and variants of APES B estimator in white noise. (b) Variance of PII APES N and variants of APES B estimator in white noise. Figure 4.9: Effect of bootstrap length on PII APES B estimator in white noise. For all the cases, N = 64. 89 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA Choice of β The derivation of CRLB for Ibal in coloured noise, given in Appendix B, shows dependence on β (see eqn B.17). However, for the noise conditions considered in our simulations, the terms of the equation that depend on β were extremely small, compared to the term independent of β, i.e.a11 a22 . Hence, we could not see any significant change in bias or variance of the estimators, when β was changed and so, we do not report those results. The same is true for the CRLB of DOLP and phase too. With this, we conclude our analysis of the various estimators for PII and proceed to analyze the estimators of DOLP information. 4.3 Imaging using Degree of polarization The DOLP of scattered light was defined in eqn 2.4. An imaging scheme that uses DOLP as the visualization parameter is reported in [28], where DOLP is defined as DOLP = Ik − I ⊥ Ik + I ⊥ (4.37) where, Ik and I⊥ refer to the co-polarized and cross-polarized intensities. It seems as though there is no relationship between the above definitions of DOLP. However, by substituting β = π 2 into eqns 4.32 and 4.33, we get the DOLP as given by eqn 4.37 to be Ibal = DOLP = Iscat p Q2s + Us2 Is (4.38) The last relationship is obtained by substituting for Ibal and Iscat from eqn 4.4. This proves eqn 4.37 to be correct, in the light of eqn 2.4. A pertinent question at this juncture is, how the DOLP based images differ from polarization intensity images. The numerator of the expression for DOLP (eqn 4.37) corresponds to the polarization intensity image. However, by scaling this information with the sum of the co-polarized and cross-polarized intensities, we normalize the polarization CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 90 information by the amount of unpolarized light present at the pixel location. Hence, the resultant images indicate the purity of the polarized light received. It can be observed that the DOLP values are both lower and upper bound (0 ≤ DOLP ≤ 1). If DOLP is 0, the received radiation does not contain any polarized light. As the DOLP increases and reaches 1, it symbolizes increasing polarized component in the received radiation, finally reaching a state where the received radiation is totally polarized. Images using DOLP are useful if the scene is non-uniformly illuminated, and more so, in passive imaging involving specular objects. DOLP can give information about the nature of the objects apart from their position in a scattering medium. Smooth surfaces like metals and glass reflect light well. It is a known fact that reflection can induce varying degrees of polarization depending on the angle of incidence. Hence, in passive imaging, one may be able to distinguish metallic and shiny surfaces by looking at the degree of polarization of the received radiation. This idea was utilized to distinguish regions abraded differently, in the experiment mentioned in [28]. Though the information of DOLP has been explicitly given in the PDI schemes, there is no such mention of a measure of DOLP in the PMI schemes given in [27, 31]. It can easily be observed that, in the PMI scheme, the DOLP information can be obtained as the ratio of the amplitude of the sinusoidal component to the DC component, since Iscat corresponds to the DC component. Thus, we see that both the PDI and the PMI schemes are capable of estimating both the polarization intensity and degree of polarization parameters. In the same manner, all the estimators discussed in the previous section can also estimate the DOLP. However, there can be a host of methods that can estimate DOLP, since the sinusoidal amplitude and DC values can be estimated by different techniques. We do not consider such a situation here, and assume that the same estimators will be used to estimate both the numerator and denominator quantities, though it may be possible to obtain improved performance by varying the estimation technique for the numerator and the denominator. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 4.3.0.9 91 The PDI DOLP estimators A general PDI DOLP estimator can be mathematically expressed as, d P DI = 1 DOLP K K X k=1   N P 2K −1 n=0 N P 2K −1 n=0
 Ik (n) − I⊥ (n) 
Ik (n) + I⊥ (n) (4.39) where, Ik (n) and I⊥ (n) are given by eqns 4.32 and 4.33. Here, we have assumed N to be even for mathematical convenience, though this is not a necessary condition. Two estimators which represent the extreme cases of the above general estimator are d P DI1 = DOLP and P N2 −1 n=0 P N2 −1 n=0
Ik (n) − I⊥ (n)
Ik (n) + I⊥ (n) −1   2 Ik (n) − I⊥ (n) 2 X = N n=0 Ik (n) + I⊥ (n) (4.40) N d P DI2 DOLP (4.41) The latter has the form of a chunked estimator of DOLP. It is clear that all the estimators derived from eqn 4.39 will be biased, since the numerator corresponds to scaled versions of the polarization intensity data, which was shown to be biased. However, if β = π 2 in eqns 4.32 and 4.33, or if the relative (rather than the d P DI can be used. Howactual) values of DOLP in a scene are of interest, then, DOLP ever, if β → 0, the noise may completely obscure the DOLP information. Henceforth, we assume β = π2 , so that the numerator is an unbiased estimate of Ibal . d P DI is that of a ratio estimator, with both the numerator and the The form of DOLP denominator being Gaussian random variables which are uncorrelated in case the noise is white, but correlated in case of coloured noise. The theoretical analysis of this estimator is by no means an easy task. We can do no better than refer to [70] for an understanding of the distribution and densities of such ratios. However, to compare the estimators, we CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 92 need to know the bias and variance of these estimators. Such a study has been reported in [71]. It has been found that d P DI2 estimator is biased and DOLP d P DI1 is asymptotically unbiased. • DOLP • The variances of the estimators are equal and diminish to zero asymptotically. Hence, both the estimators are consistent. d P DI1 . It can also be observed that the memFor these reasons, it is better to use DOLP d P DI1 is lesser than that of DOLP d P DI2 , thereby strengthening ory requirements of DOLP its case. There is also a mention of another ratio estimator in [71] that is unbiased. But, it needs prior knowledge of the noise variances and covariances. Since this information is not available in our problem, we do not study that estimator. 4.3.0.10 The PMI DOLP estimator In our analysis of the PII, we found that, for β 6= π2 , we can obtain only ML estimate of Ibal and not its MVU estimate. We also observed that, the ML estimate of Iscat also corresponds to its MVU estimate in both coloured and white noise. Hence, we can obtain MLE of DOLP using the invariance property of the MLE [61]. The transformation ˆ b DOLP = IˆIbal being non-invertible, the MLE maximizes a modified likelihood function, scat as explained in [61]. By the theorem on the asymptotic behaviour of the MLE [61], we conclude that the MLE achieves the CRLB asymptotically. The CRLB for DOLP, for white and coloured noise cases have been derived in Appendices A and B, respectively. 4.3.1 Comparison of DOLP imaging estimators As in the case of PII, we used real noise and also resorted to Monte-Carlo simulations to study the performance of the various estimators of DOLP. The parameters for the simulations were the same as we reported for PII. A common observation that could be made across all estimators is that, their performance in coloured noise was very poor for SNR from around -15 dB to 5 dB. This is the region of transition from positive to negative bias for polarization intensity estimators in coloured noise, as we saw in the CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 93 previous section. However, the performance of the estimators on real noise did not show such droop in performance. This observation is perhaps due to the noise being nearly white. Hence, estimators of DOLP can still be used at these SNR. We now mention the best estimators in each class for real and simulated data. • The PDI N and the PMI N showed minimum variance and bias under all noise conditions and are hence the best estimators in their respective classes. • Among the MLE estimators, the MLE N showed better performance overall, compared to all other estimators of its class. • APES CB showed the best performance in white noise. No clear trend could be observed in the performance of the APES estimators in coloured noise. For real noise, it was observed that APES N performed better in noise 1 and APES CB performed better in noise 2. We consider APES CB as the representative for its class, due to its superior performance in white noise and in noise 2. • Among the APESR estimates, APESR CB performed better in white noise. The estimators being unreliable in coloured noise, we consider APESR CB as the representative of this class. However, in real noise, APES N performed relatively better. Figs. 4.10, 4.11 and 4.12 shows the bias and variance of the best estimators from each class in white, AR1 and noise 1, respectively. We can observe that PDI N and PMI N estimators seem to be the best for estimating DOLP too. Though MLE N and APES estimators also perform as good as PMI estimators, their computational costs are much higher. Hence, PDI N or PMI N should be preferred over them. As mentioned in the section on polarization intensity estimation (section 4.2), the behaviour of the PDI N and PMI N, as reported here, are based on certain conditions. The conditions for PDI imaging can be usually satisfied in active imaging, while the conditions on PMI scheme can be satisfied in both active and passive imaging schemes. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of estimators of DOLP in white noise. (b) Variance of estimators of DOLP in white noise. Figure 4.10: Performance of estimators of DOLP in white noise. 94 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of estimators of DOLP in AR1 noise. (b) Variance of estimators of DOLP in AR1 noise. Figure 4.11: Performance of estimators of DOLP in AR1 noise. 95 96 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of estimators of DOLP in real noise. 1 PDI N PMI N MLE N APES CB APESR CB 0.9 0.8 0.7 Variance 0.6 0.5 0.4 0.3 0.2 0.1 0 −25 −20 −15 −10 −5 0 5 10 15 20 SNR dB (b) Variance of estimators of DOLP in real noise. Figure 4.12: Performance of estimators of DOLP in real noise. 25 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 97 We also performed simulations to see if the APES estimators provide better estimates of DOLP when the frequency of the sinusoid was slightly off from the expected value. The analysis was similar to that of section 4.2.5.4. However, here too, we found the PMI N to perform better than the APES estimators at SNR below 25 dB, when the frequency of the sinusoid was 0.115, instead of the usual 0.125. The advantages of the APES estimators can only be seen at frequencies far away from that resolved by the PMI estimator and have been detailed in [62]. 4.3.1.1 A case against the chunked PMI DOLP estimator Due to the reasons we cited in section 4.2.5.5, it intuitively looks like, it may be better to use PMI C instead of PMI N, when the frequency of the sinusoid can vary slightly around the expected value. We conducted simulations on the same lines of section 4.2.5.5, to study the usefulness of PMI C in such cases. Our study showed that chunking worsened the performance of the PMI DOLP estimators. The variance of the estimator increased with increase in chunk length. However, the bias of the chunked estimators decreased as chunk length increased, beyond SNR of -20 dB. But, the gain in terms of bias was found to be very little. Hence, we conclude that chunking does not improve the performance of PMI DOLP estimator. 4.3.1.2 Effect of chunking on APES DOLP estimator To study the effect of chunking parameters on estimation of DOLP, we chose the APES estimator, which showed positive response to chunking. We conducted Monte-Carlo simulations on the lines of section 4.2.5.6 to study the effects of these parameters. By the very nature of definition of DOLP, we expected the same results, as in section 4.2.5.6, even for DOLP imaging. It was indeed found to be so in white noise, where, processing small chunks decreased the variance up to around 20 dB. The bias was also found to decrease up to SNR of -5 dB, beyond which, larger chunks showed lower bias. At high SNR, bias was found to decreases as chunk length increased, with the least being that of APES N estimator. However, in coloured noise, the behaviour of the estimators were erratic, due to the transition of bias from positive to negative values. The APES N seemed to be the CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 98 most stable estimator in coloured noise. Fig. 4.13 shows the bias and variance of APES C estimators with varying chunk lengths in white noise. Increasing the overlap length gave lower bias and variance in white noise, and hence overlapping should be useful in estimating DOLP values. Fig. 4.14 shows the results of simulations conducted to study the effect of overlapping in white noise. Since the estimators showed erratic performance in coloured noise, we do not report those results. 4.3.1.3 Effect of bootstrap length on APES DOLP estimator We conducted simulations on the lines of section 4.2.5.7 to study the effect of bootstrap length on the APES DOLP estimator. Fig. 4.15 shows the results of simulations. Here are the important observations related to that study. • Bootstrapping decreased the bias and the variance of APES DOLP estimators at low SNR. However, at high SNR (beyond 10 dB), bootstrapping induced slight positive bias. Hence, bootstrapping can be useful for estimating DOLP at low SNR (below 15 dB) in white noise. • In coloured noise, due to the unstable behaviour of the APES DOLP estimator, no clear trend could be observed and APES N was more stable than the other estimators. With this, we conclude our analysis of the various estimators of DOLP. The choice of suitable frequency should be based on the considerations mentioned in section 4.2.5.8. Varying the number of data points for analysis did not change our observations. As reported in section 4.2.5.8, varying β did not have much impact on DOLP estimation. One thing that we hypothesize, but did not test is that, one can perhaps obtain better estimates of DOLP by estimating Ibal using chunking, and estimating Iscat without chunking and taking their ratio. Such a hybrid estimator may yield better results than bootstrapping alone, or chunking alone. Next, we study the various estimators used for Polarization orientation imaging (POI). CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of APES C DOLP estimator in white noise. (b) Variance of APES C DOLP estimator in white noise. Figure 4.13: Effect of chunk length on DOLP APES C estimator in white noise. 99 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 100 (a) Bias of APES C DOLP estimator in white noise. (b) Variance of APES C DOLP estimator in white noise. Figure 4.14: Effect of overlap length on APES C estimator in white noise. For all the cases, the chunk length was 32. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of APES B DOLP estimator in white noise. (b) Variance of APES B DOLP estimator in white noise. Figure 4.15: Effect of bootstrap length on DOLP APES estimator in white noise. 101 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 4.4 102 Polarization orientation imaging The received radiation can contain different states of linear polarization due to various factors. It could be due to extraneous polarized sources deployed for other imaging purposes. It could also be due to changes in polarization state induced by reflection. The latter case is more commonly encountered in passive imaging scenarios, where specular surfaces can change the orientation of the plane of polarization of incident linearly polarized, or even induce polarization to the incident unpolarized light. In any case, the ability to distinguish different states of linear polarization adds value to the imaging methodology. Further, this can be obtained at no extra cost, while using any of the schemes, except the PDI. We next analyze as to how we can obtain the polarization state information. Assume that with the PMI scheme of Fig. 4.1, an experiment is conducted with an arbitrary orientation of the plane of polarization of the source, yielding data that follows eqn 4.5. With all the parameters being the same, we study the effect of a change in the orientation of the plane of polarization of the source by an angle γ, on the recorded data. In eqn 4.4, we assumed that the angle φ represents the orientation of the rotating polaroid with respect to the horizontal. Now, we need to replace φ by φ + γ. This change essentially leads to a change in the value of β, which corresponds to the change in phase of the recorded sinusoids. i.e., a change in the orientation of the plane of polarization manifests itself as a change in the phase of the recorded sinusoids, and does not affect the sinusoidal and DC amplitudes. The Polarization intensity imaging and the DOLP imaging cannot capture this information, since they ignore the estimated phase of the sinusoids, even though all the estimators except the PDI estimator can give this information. Thus, we can use the same estimators to obtain all the three visualization parameters, namely the polarization intensity, the DOLP and the polarization orientation (PO), which give distinct information about the scene being imaged. 103 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 4.4.1 Polarization orientation imaging - the case of white noise The CRLB for estimation of the phase of a sinusoid in white noise has been derived in Appendix A and is given by eqn A.8. However, it has been found that there does not exist an MVU phase estimator that attains the CRLB for this case in white Gaussian noise [61]. This result is applicable to our problem too. But, the MLE of phase exists due to the invariance property, and is given approximately by [61] (from eqn 4.25 and 4.26) β̂M LE = arctan Iˆbal sin βM V U Iˆbal cos βM V U ! = arctan PN −1 n=0 Ir (n) cos PN −1 n=0 Ir (n) sin
! 4πn M
4πn M (4.42) where, M is the periodicity of the rotating polaroid, and N is the number of data points available for analysis. β̂M LE can also be obtained as the argument of the complex DFT coefficient corresponding to the frequency of the sinusoid. Since the PMI estimator obtains the DFT coefficient from which we can estimate the phase of the sinusoid, we conclude that the PMI phase estimator obtains the MLE estimate of phase. The asymptotic variance of the PMI phase estimator has been shown to reach the CRLB given by eqn A.8. The transformation function used to estimate the phase is not an invertible function, since the transformation maps two parameters, namely Iˆbal sin βM V U and Iˆbal cos βM V U into one variable β̂M LE . Hence, β̂M LE actually maximizes the modified likelihood function in the manner explained in [61]. 4.4.2 Polarization orientation imaging - the case of coloured noise The CRLB for phase of the sinusoid in coloured noise has been derived in Appendix B, and is given by eqn B.18. It can be observed that the CRLB for β depends not only on the amplitude of the sinusoid Ibal (or the SNR) as in the case of white noise, but also on β itself. The value of β, for which the variance is minimized, has also been derived CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 104 in Appendix B. However, the dependence on β was found to be very weak, essentially making the CRLB invariant to β. As in the case of white noise, the PMI estimate of phase is given by the argument of the complex DFT coefficient corresponding to the frequency of the sinusoid. The APES, APESR estimators as given in [62] also give complex amplitude estimates, leading to estimates of phase. In the case of MLE, we need to use eqn 4.42. Due to the asymptotic properties of MLE, the estimator achieves the bound given in B.18 when large number of data points are available for analysis. Once again, β̂M LE actually maximizes the modified likelihood function, as in the case of white noise. 4.4.3 Comparison of POI estimators We now compare the different estimators for phase discussed so far. Before doing so, we make a few observations that are common to all POI estimators. • All the estimators, irrespective of their class, showed high and erratic bias up to −5dB in white noise, beyond which they performed better. A similar behaviour was observed in noise 2 also. In coloured noise, such a behaviour persisted up to +5 dB, thereby decreasing the reliability of the estimators. From these observations, we conclude that reliable phase estimation is possible only at relatively high SNR. • In white noise, the variance of the estimators behaved as expected and governed by CRLB. However, in coloured noise, the variance of the estimators remained nearly constant up to around +5 dB, much below the CRLB, beyond which, they behaved as expected. • Chunking seemed to help estimators in both simulated and real noise. Next, we proceed to list the best estimator in each class for both simulated and real noise. • The PMI N showed the best performance in its class in white and coloured noise. In real noise, PMI N showed the least variance but higher bias than PMI C. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 105 • MLE CB showed the best performance on simulated data under all noise conditions. In noise 1, MLE C performed better than MLE N in terms of bias. However, in terms of variance, MLE N showed the best performance. • Among the APES estimators, the APES CB showed the best performance on simulated data. However, APES C and APES N preformed better in noise 1 and 2, respectively. We use APES CB to compare with the estimators of other classes. • Among the APESR estimators, there was no clear trend in the behaviour, neither in simulated noise nor in real noise. Since APESR CB showed better performance in white noise, we use it to compare with estimators of other classes. Fig. 4.16 shows the bias and the variance of the phase estimators in white noise. Similarly, Fig. 4.17 shows the bias and variance of the same estimators in AR1 noise. Fig. 4.18 shows the performance of the various estimators in noise 1. As we can observe, for POI, there is no clear advantage for any particular estimator. Since the PMI phase estimator is simpler than the rest and shows the best performance in white noise, we conclude that PMI estimator should be the preferred choice for POI. 4.4.3.1 A case for the APES POI estimator We also studied the performance of the estimators in white noise when the frequency of the sinusoid was 0.115, instead of 0.125. Fig. 4.19 shows the bias and variance of the estimators under such conditions. It is observed that PMI N and MLE N give biased estimates of phase, which remain constant even at high SNR. However, the APES N estimator gives low bias and nearly 0 bias after an SNR of 5 dB. This is the advantage of the APES estimator, over the PMI estimator. However, the variance of the APES estimator is higher than that of the PMI N estimator. This can be reduced by resorting to chunking and overlapping. We see a strong dependence of phase on the frequency of the signal. Hence, for phase estimation, when the frequency could be off by a small margin, the APES estimators CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of phase estimators in white noise. (b) Variance of phase estimators in white noise. Figure 4.16: Performance of phase estimators in white noise. 106 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of phase estimators in AR1 noise. (b) Variance of phase estimators in AR1 noise. Figure 4.17: Performance of phase estimators in AR1 noise. 107 108 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) Bias of phase estimators in real noise. 3 PMI N MLE B1 APES CB APESR CB 2.5 Variance 2 1.5 1 0.5 0 −25 −20 −15 −10 −5 0 5 10 15 20 SNR dB (b) Variance of phase estimators in real noise. Figure 4.18: Performance of phase estimators in real noise. 25 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 109 (a) Bias of PMI N, APES N and MLE N phase estimators in white noise. (b) Variance of PMI N, APES N and MLE N phase estimators in white noise. Figure 4.19: Performance of PMI, MLE and APES estimators in white noise, when N is not an integral multiple of M . For this example, f = 0.115. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 110 should be a good choice. However, the APESR estimator performs extremely poorly under such circumstances as can be seen in Fig. 4.16 and should be avoided. The performance of all the estimators being erratic in coloured noise, we do not report those results here. Also PMI C did not perform well in this case. The bias of the PMI estimator in fact deteriorated with chunking. 4.4.3.2 Effect of chunking on the APES POI estimator Having found the utility of the APES estimator for phase estimation, we also studied the effect of chunking parameters on the performance of the APES C estimator. Simulations were performed with the same parameters, as in the case of polarization intensity imaging. It was found that using smaller chunks reduced the variance of the APES C estimator at all SNR and brought down the bias at low SNR. However, at SNR beyond 5 dB, smaller chunk lengths induced greater bias. Fig. 4.20 show the results of the simulations. Similarly, we studied the effect of overlap length on APES C estimator. It could be observed that with increased overlapping, the bias and the variance of the estimates came down by very small amounts. We also studied the effect of bootstrapping on the APES estimator. Simulations showed that bootstrapping did not show any improvement on the performance in white noise. 4.5 A potpourri of visualization parameters We defined three different visualization parameters for polarization based imaging schemes. We also studied the various estimators for each of these visualization parameters and studied their characteristics. In all the schemes discussed, the visualization parameters are represented as gray scale images. If we want to study all the visualization parameters of a scene, we need to study three different images. In this section, we propose a new scheme of rendering the polarization information, where, the visualization parameters are intuitively mapped to various aspects of a colour image, thereby giving a holistic view of the scene. We also show the advantages of such a representation. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 111 (a) Bias of APES N and variants of APES C phase estimators in white noise. (b) Variance of APES N and variants of APES C phase estimators in white noise. Figure 4.20: Effect of chunk length on the performance of APES N and variants of APES C phase estimators in white noise. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 112 The polarization magnitude information inherently has the notion of intensity, i.e., the intensity of polarized light. The DOLP signifies the purity of the polarized radiation and parallels the idea of saturation in color images. In colour images, the greater the purity of a colour, higher will be the saturation. The polarization orientation parameter distinguishes different states of linear polarization and is akin to different hues in a scene. Thus, the three visualization parameters correspond intuitively to the parameters of a normal colour image and hence can be rendered so. We next see some simulation results to illustrate this aspect. We synthesized 32 images of size 160 × 100 pixels to study the suitability of fusing the visualization parameters. The data at each pixel location was synthesized according to eqn 4.5. We considered the noise to be white for these simulations, but the results for the case of coloured noise will be no different. The synthesized images were divided into four quadrants, as shown in Fig. 4.21(a). The value of Iscat was chosen to be 170 in quadrants 1 and 4, and 240 in quadrants 2 and 3. Within each quadrant, we chose rectangular subregions of size 80 × 50, as shown in Fig. 4.21(a). The value of Ibal was chosen to be 0 at all locations, except the subregions. It was chosen to be 0.2 in the subregions of every quadrant, thus giving a DOLP of 0.0012 in the subregions of quadrants 1 and 4 and 0.008, in the subregions of quadrants 2 and 3. Before the addition of noise, the images were blurred using a Gaussian mask of size 9 × 9, to simulate the blurring due to the optical elements. Though we have not analyzed the nature of the blur in the actual experimental setup, we use the Gaussian blur only to study the probable effect of blurring. To the blurred data, we added white Gaussian noise, at every pixel location, across images. The noise at each pixel location had a variance of 0.05, leading to a SNR of nearly -4 dB in all the subregions. The reason for choosing such small values for Ibal and DOLP is to show the robustness of the polarization imaging technique even at relatively low SNR, when the individual images in the series do not convey any visual information about polarization by themselves. The information becomes evident after processing. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 113 The phases of the sinusoids in the subregions were chosen to be 2π/9, 0, 5π/18 and 15π/9, respectively. The phases were chosen arbitrarily, to illustrate the feasibility of rendering the polarization information as colour images and do not have anything specific to the processing itself. We have used the PMI estimator to estimate the various visualization parameters, since the PDI scheme is incapable of estimating the phase information. The parameters can also be obtained using the APES, APESR or the MLE estimator. Fig. 4.21(b) shows a representative image of the set of 32 images which looked nearly alike. Fig. 4.21(c) shows the histogram equalized version of the representative image. As we can observe, we do not find any discernible difference among the subregions, but, we find the difference in gray scales of the four quadrants. The same was true in the case of other images in the series too. Fig. 4.21(d) shows the result of polarization magnitude analysis. It can be observed that the four subregions have the same polarization magnitude information. Fig. 4.21(e) shows the result of DOLP analysis, wherein, we can observe slightly higher DOLP in the subregions of quadrants 1 and 4, as compared to those of quadrants 2 and 3. The magnitude of DOLP being very small, we perhaps could not have expected a drastically different result, showing the quadrants very differently. Fig. 4.21(f) shows the result of polarization orientation analysis. The subregions are clearly distinguishable from the background. However, only the subregion of the first quadrant stands apart uniquely from the others. Thus, though it is very attractive to use the POI to differentiate regions with polarization information from others that do not, it is not easy to differentiate the orientation from these results. Fig. 4.21(g) shows the result of fusing all the visualization parameters into a colour image, as explained earlier. The colour image clearly shows the different polarization orientations of the subregions. A keen observation also shows the lower saturation levels of the colors in quadrants 2 and 3, as compared to the colours of quadrants 1 and 4. CHAPTER 4. PROCESSING POLARIZATION-RICH DATA 114 This corresponds to the difference in DOLP of these subregions. It is clear that rendering the parameters as a colour image can provide better insight into the various polarization parameters, as compared to rendering them as gray scale images. 115 CHAPTER 4. PROCESSING POLARIZATION-RICH DATA (a) (b) (c) (d) (e) (f) (a) (b) (c) (d) (e) (f) (g) Break-up of synthetic images. A representative of the 32 images in the series. Result of histogram equalizing (b). PMI N magnitude estimation result. PMI N DOLP estimation result. PMI N PO estimation result. Result of fusing (d), (e) and (f). (g) Figure 4.21: Results of fusing the visualization parameters. Chapter 5 Imaging results In this chapter, we present results of the imaging experiments conducted with various scattering media. The experiments were conducted with the setup explained in section 4.2.5.2. Table 5.1 lists the various media that we used for the experiments reported in this chapter. Table 5.2 gives the statistical information of the corresponding data sets, obtained by measuring the different quantities mentioned, over a region of 10 × 20 pixels. In all the imaging experiments mentioned in this chapter, the periodicity of the rotating polaroid was 16, and hence, the periodicity of Ibal was 8, or its frequency was 0.125. Also, except with data of SET 1, in all the experiments reported here, the object (shadow of the object in this case) was not visible in any of the individual images obtained, even after processing by various techniques. In SET 1, a very faint shadow of the object could Table 5.1: Data sets used for analysis. Name SET 1 SET 2 SET 3 SET 4a Scattering particles Polystyrene spheres Polystyrene spheres Polystyrene spheres Polystyrene spheres of of of of diameter diameter diameter diameter ls (µ) 2.97µ 334.7 2.97µ 314.0 2.97µ 2.97µ - l∗ (µ) τ 1760 5.68 1652 6.05 - The medium was the same as for SET 3, but the source was incoherent white light, instead of the laser, used in obtaining SET 3. a 116 CHAPTER 5. IMAGING RESULTS 117 Table 5.2: Statistics related to the data sets used for analysis. Name Average SNR (dB) Average AR coefficients SET 1 −1.2 +0.078 and +0.023 SET 2 −6 −0.002 and +0.018 SET 3 −15 +0.055 and +0.052 SET 4 −8 +0.109 and +0.053 be perceived in some of the images that were close to co-polarization location. However, we use data from SET 1 to illustrate the various aspects of polarization based imaging, since the visibility of the shadow was very poor. Similar results could be obtained with other data sets too. All the images illustrated in this chapter are histogram equalized versions of the actual results obtained, since the actual images had very low contrast. We have emphasized on the results of PDI and PMI estimators, since we found them to be suitable for all imaging purposes. Processing results also endorsed the results of our analysis. Fig. 5.1 shows the results of polarization intensity imaging using PDI N estimator, when 8 co-polarized and 8 cross-polarized images of SET 1 were available for analysis. The periodicity of Ibal being 8, there were 4 positions (corresponding to the different phases of Ibal ), that could be used to obtain polarization intensity information using the PDI N estimator. Fig. 5.1(a) shows the shadow of the opaque cross without any scattering. This was the object to be imaged in SET 1. Fig. 5.1(b) shows a representative of a nearly co-polarized image obtained with scattering medium being present. Since we did not know the actual orientation of the polaroid with respect to the source radiation, we call it as nearly co-polarized. We can see that the shadow of the cross is faintly visible in Fig. 5.1(b). However, no visible details could be observed in the nearly cross-polarized images. Figs. 5.1(c-f) show the results obtained by using PDI N estimator, assuming different positions of the rotating polaroid to correspond to the co-polarized location. We can observe considerable improvement in the results from Fig. 5.1(c) to 5.1(f). In Fig. 5.1(c), we observe that there is no information about the hidden object. This is due CHAPTER 5. IMAGING RESULTS 118 to the bias of the PDI N estimator. Similar results could be observed in estimation of DOLP with the same data. Fig. 5.2 shows the comparison between polarization intensity and DOLP results obtained from SET 1. The figures on the left correspond to PII and the ones on the right are results of DOLP imaging. The images from top to bottom were obtained with N being 8, 32 and 128, respectively. There is a gradual improvement in the quality of results, as N increases from 8 to 128. We observe that the DOLP images have a brighter central region than polarization intensity images. Comparing the results with that of Fig. 5.1(a), we observe that the DOLP images are closer to the actual image than the polarization intensity images. Fig. 5.3 shows the comparison of polarization intensity results obtained from PDI N and PMI N with data from SET 1. The images on the right correspond to PMI N, while the others correspond to PDI N. The results in each row were obtained with the value of N being 8, 32 and 128, respectively. We can observe that when fewer images are available for analysis, PDI N gives better results than PMI N, and hence should be preferred. Beyond N = 256, not much difference could be observed in the results. Similar behaviour was seen in DOLP images too. Fig. 5.4 shows the comparison of results obtained from different estimators with N = 16, from SET 1 data. We see that PDI N gives the best result among all the estimators. PMI N and MLE N estimators give the next best result, which look very similar. However, comparing the computational complexity, one should choose PMI N instead of MLE N, in case PDI N is not suitable for some reason. The results of the APES estimators are less encouraging, as we could guess from our analysis in the previous chapter. Similar comparisons for data from SET 2 are given in Fig. 5.5. We now compare the results obtained with the scattering medium being the same, but the source being a coherent laser beam in one case, and an incoherent white source in another. SET 3 and SET 4 correspond to such data sets. However, the medium was not CHAPTER 5. IMAGING RESULTS 119 calibrated and we could not know the values of ls and l∗ of it. Figs. 5.6(a) and 5.6(b) show the shadows of a opaque cross obtained without scattering from two different sources. The former was obtained with a monochromatic laser source and the latter with an ordinary white light source. Once the scattering medium was introduced, in both the cases, we could not observe any details of the shadow in individual images by applying standard image processing techniques. The results of processing 128 images of SET 3 and SET 4 with PDI N DOLP estimator are given in Figs. 5.6(c) and 5.6(d), respectively. Similarly, Figs. 5.6(e) and 5.6(f) show the results of processing 128 images of the same sets with PMI N DOLP estimator. Similar results were observed with PII. We see that the results obtained with the incoherent source are far superior to those obtained with the laser. The main reason for this should be the low SNR of SET 3 as compared to that of SET 4, as given in Table. 5.2. The other reason for poor results in SET 3 could be due to the speckle noise that could be present in SET 3 induced by the coherent illumination of the laser. Though we do not have a concrete evidence for this, it is well known that speckle is formed when scattering occurs due to a coherent source and that it degrades imaging performance. The other argument that hints at the same reason is the results of SET 2, where, though the SNR was lower than that observed in SET 4, the results obtained are inferior to those obtained with SET 4. Both results were obtained by analyzing the same number of data points. However, this is a open question to be investigated, and there seems to be very little literature regarding this. But, it is clear that incoherent sources of very high powers will be needed to image objects through greater optical depths, as compared to laser sources, since incoherent sources have high beam divergence rates compared to lasers. Next, we show some advantages of POI imaging using results obtained from SET 4. To illustrate this, we use the estimates obtained by the PMI N estimator. As seen in section 4.4, POI is useful only if the SNR is relatively high, or a large number of data points are available for analysis. For illustrating the advantages of POI, we use 512 images of SET 4, where the SNR is relatively higher than other data sets. CHAPTER 5. IMAGING RESULTS 120 If a single linearly polarized source is used, then, ideally, the region not blocked by the object should yield the same phase information at every pixel location. Depending on the SNR, the values of estimated phase will only fluctuate slightly around the mean value. However, wherever the object blocks the ballistic light, the phase estimation would yield random results. Thus, the correlation in the intensity values of neighbouring pixels in regions receiving the ballistic component should be much higher than that of regions blocked by the object. As we show, this information can be exploited in segmenting the image into target and background regions. This sort of segmentation and hence, POI is very useful in defense related applications. Moreover, the analogy can be extended to regions containing different polarization orientations also. For testing the validity of our hypothesis, we considered rectangular blocks of size 9 × 9 and 15 × 15, centered around every pixel location, excluding the boundary pixels. We cross-correlated the data in every block with blocks around the adjacent 8 neighbors and took the maximum of the cross correlation values as the result for that pixel. We finally plotted these results as an image. Instead of taking the maximum of the correlation values, we also took the average of the correlation values, and observed similar results. We found that the resultant images had histograms which clearly showed two modes. This information can be used to segment the images into regions corresponding to the hidden object and the background. Fig. 5.7(a) shows the image of the object without scattering. Fig. 5.7(b) shows the POI result obtained from 512 images of SET 4. Fig. 5.7(c) shows the result obtained from the cross-correlation technique described above, when the block size was 9 × 9. Fig 5.7(d) shows the result with blocks of size 21 × 21. Fig. 5.7(e) shows the histogram of the image in Fig. 5.7(d). Fig. 5.7(f) shows the result of thresholding Fig. 5.7(d) at a value lying in the valley between the two modes of the histogram. As we can observe, though the exact boundary of the object is not visible, we get an idea of the presence of a hidden object in the medium. The processing can be followed by morphological operations to get better results. This sort of information is what is most of the time needed in defense applications. This is a very useful result, as far as POI applications to defense is concerned. CHAPTER 5. IMAGING RESULTS 121 Though we can segment polarization intensity and DOLP images, the POI segmentation is robust when the variation in the estimates of phase will be smaller than the variation in the intensity of Ibal and Iscat . Though Ibal and Iscat can have larger distributions naturally, it is unlikely that the phase has such large distributions. This is the advantage of using POI results for segmentation. It was also observed that the intensities of PII and DOLP results did not show any modes as such, and hence we need to do further processing to be able to get reliable segmentation results. Many other post processing techniques can be applied to the results obtained by PII, DOPI and POI. We can fuse the three results to obtain colour images in case multiple states of linear polarizations exist in the received radiation, as explained in section 4.5. We did not have any data set which had multiple states of linear polarization, and hence we do not have such results. The segmentation scheme discussed earlier is fairly rudimentary. We can perhaps develop many more such segmentation algorithms, specifically aimed at polarization imaging. We also conducted experiments with the medium being mist, and could get positive results after processing. We observed that for the 2.97µ diameter particles, when the optical thickness τ reached 6.77, we could not retrieve the shadow, even after processing 512 images. This was the limit to which we could image with those particles. However, with particles of 0.11µ diameter, we could image up to optical thickness of nearly 40. This observation endorses the result that polarization based schemes work better in media containing small sized scatterers, than in media with large scatterers [35, 52, 37]. 122 CHAPTER 5. IMAGING RESULTS (a) (b) (c) (d) (e) (f) Figure 5.1: PII to illustrate the bias of PDI N estimator; N = 16. (a) (b) (c)-(f) Image of the object without scattering. A representative image of the series. PDI N results obtained with various co-polarization locations. 123 CHAPTER 5. IMAGING RESULTS (a) (b) (c) (d) (e) (f) Figure 5.2: Comparison of results of PII and DOLP imaging. (a) (c) (e) PII using PDI N; N = 8 PII using PDI N; N = 32 PII using PDI N; N = 128 (b) (d) (f) DOPI using PDI N; N = 8 DOPI using PDI N; N = 32 DOPI using PDI N; N = 128 124 CHAPTER 5. IMAGING RESULTS (a) (b) (c) (d) (e) (f) Figure 5.3: Comparison of PII results obtained using PDI N and PMI N. (a) (c) (e) PII using PDI N; N = 8 PII using PDI N; N = 32 PII using PDI N; N = 128 (b) (d) (f) PII using PMI N; N = 8 PII using PMI N; N = 32 PII using PMI N; N = 128 125 CHAPTER 5. IMAGING RESULTS (a) Actual Image (b) PDI N (c) PMI N (d) MLE N (e) APES N (f) APESR N Figure 5.4: Comparison of PII results obtained from 16 images of SET 1. 126 CHAPTER 5. IMAGING RESULTS (a) Actual Image (b) PDI N (c) PMI N (d) MLE N (e) APES N (f) APESR N Figure 5.5: Comparison of DOPI results obtained from 128 images of SET 2. 127 CHAPTER 5. IMAGING RESULTS (a) (b) (c) (d) (e) (f) Figure 5.6: Comparison of results obtained with coherent and incoherent sources. (a) (c) (e) Image without scattering (Laser) DOPI using PDI N; N = 128 (Laser) DOPI using PMI N; N = 128 (Laser) (b) (d) (f) Image without scattering (whitle light) DOPI using PDI N; N = 128 (white light) DOPI using PMI N; N = 128 (white light) 128 CHAPTER 5. IMAGING RESULTS (a) (b) (c) (d) (e) (f) Figure 5.7: Segmentation of a POI result. (a) (c) (e) Image without scattering 9 × 9 block processing result Histogram of (d) (b) (d) (f) Resutl of PMI N POI; N = 512 15 × 15 block processing result Result of segmenting (e) Chapter 6 Conclusions and topics for further research 6.1 Contributions of the thesis In the previous chapters, we studied various visualization parameters that can be used to render polarization information contained in the received radiation. We also studied various estimators for measuring these visualization parameters. We found that the PDI and the PMI schemes are suitable for estimating the visualization parameters discussed. Interestingly, they are also the simplest of all the estimators studied. Owing to the bias of the PDI estimator and its inability to capture polarization orientation information, PMI N is the true all purpose estimator. As we saw, PDI N is a particular case of the PMI N estimator. Hence, we conclude that PMI N should be the estimator of choice to estimate all the polarization related parameters of the received radiation. We also studied as to how chunking and bootstrapping can help improve the performance of the estimators. However, the improvement gained by such schemes were marginal. There was no significant difference in the results obtained. But, their utility has been proved using simulations. We also introduced some post estimation techniques like segmentation of the POI images to segregate target and background regions and fusing the three visualization parameters to obtain a holistic view of the scene. Due to lack of 129 CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH 130 experimental data, we could not illustrate the techniques on actual data. To summarize, here are the main contributions of the thesis. • We have made an extensive literature survey of light scattering and optical imaging, to position the problem of continuous wave, polarization based direct imaging in its right place among other optical imaging techniques. We have also studied its advantages and limitations. • To the best of our knowledge, for the first time, we have classified the imaging schemes based on the visualization parameters, and pooled in most of the techniques in one place and explained how and when the various schemes are useful. • We have built a mathematical framework to compare the various schemes mentioned in the literature. The framework begins with modeling the received scattered radiation and ends with predicting the performance of the various estimators by making use of principles of estimation theory. • Through simulations and theoretically, we have analyzed the performance of various estimators for the different visualization parameters under varying noise conditions. • We have also briefly explained the possibility of using post estimation techniques like segmentation, to derive more useful information from the estimated quantities. • The introduction of polarization orientation imaging, the idea of combining the visualization parameters in to a colour image and the correlation based segmentation of POI results are to the best of our knowledge, totally new concepts in this field. • In short, we have been able to leverage the tools of signal processing, especially estimation theory, for analyzing and improving continuous wave, polarization based, direct imaging techniques. 6.2 Extension to the case of circular polarization The emphasis in this thesis has been on imaging schemes which analyze linearly polarized light. However, it can be easily observed that the same techniques can be used CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH 131 for PII and imaging using DOP of circularly polarized light. The notion of POI imaging will have to be modified in this case, since there can be only two pure states of circular polarization. We can perhaps allocate two particular hues to denote circularly polarized states, and the rest can be used to represent various linear polarization states, in a all encompassing linear and circular polarization based imaging scheme. For realizing such a scheme, suitable changes will have to be made in the imaging setup [8], in order to capture the information corresponding to the circularly polarized states. With such a scheme, two colour images can be obtained, one corresponding to the visualization parameters of the linearly polarized states, and the other, to the circularly polarized states. 6.3 6.3.1 Topics for further research Correlation based processing techniques We studied the various algorithms for processing polarization rich data. There is another processing technique for the same purpose, based on 2-D correlation. Here, a real, paraxial image of objects hidden behind a multiple-scattering barrier can be obtained from the light diffused through the barrier [26]. The processing involves correlating the speckle pattern produced by a known, reference source with that of the unknown distribution due to the object. Most of our data sets did not have prominent speckle patterns, due to averaging over long periods of time. Moreover, we did not have data with the source alone. Hence, we did not pursue this algorithm. It has been reported that highresolution images can be obtained using this simple technique. We have not addressed these schemes. 6.3.2 The imaging problem vis-a-vis the image transfer problem An interesting issue that needs to be understood well is the relationship between the imaging problem, vis-a-vis, the image transfer problem. The latter can be imagined to be an extreme case of the imaging problem, with the object being outside the scattering medium, towards the source. If the medium were to be linear, then, irrespective of the position of the object, the result obtained would be the same. However, such is not the CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH 132 case. Thus, if we can understand the relationship between these two problems, we can use the point-spread function analysis of the image transfer problem to solve the imaging problem. However, for the problem of polarized light transfer, we need to know the vectorial point spread function of light. Next, we share our ideas about this problem. 6.3.3 The vectorial point spread function of light In this thesis, we analyzed the 1-D signal processing aspects of the problem of direct, polarization based imaging schemes. However, as we observed, even with the best estimators, useful results could be obtained only up to a very small optical thickness. The main reason for this was observed to be the spread of scattered light, which can be characterized by the point spread function. Fig. 6.1(a) and 6.1(b) show the results of using detection theory to test the presence of sinusoids at various pixel locations of SET 2 (see section 5), assuming the noise to be white. The images are binary, with dark regions denoting regions of absence of sinusoids, and bright regions denoting their presence. The first image corresponds to the probability of false alarm (pfa) being 10−14 , and the second, corresponds to pfa being 10−18 . It can be observed that despite the pfa being so low, most of the regions show the presence of a sinusoid, including regions of geometric shadow, thereby emphasizing the spread in the ballistic components. It has been reported in [41], that the least spread among the polarization based schemes discussed was that of the polarization difference. An interesting exercise would be to analyze the spread function assuming the ballistic component to vary sinusoidally. Such an analysis will establish the relative strengths of each of the imaging schemes. Detailed analysis of spread function requires the knowledge of radiative transfer theory [72]. The point spread function analysis for transfer of unpolarized light (scalar light) through scattering media has been covered exhaustively in [9]. For understanding the spread of polarized light, we need to extend the analysis of [9] to the vectorial nature of light. If such an analysis can be carried out, then the comparison of various imaging schemes in terms of their abilities to resolve objects hidden in scattering media can be CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH (a) pfa=10−14 133 (b) pfa=10−18 Figure 6.1: Results of sinusoidal detection on data of SET 2; N = 512. carried out. Such an analysis is essential to improve the performance of the polarization based schemes. With the knowledge of the vectorial point spread function, we can deconvolve the results obtained using 1-D processing schemes, to obtain improved results. 6.3.4 Post processing procedures All the results that we have reported have been obtained by plotting the results obtained from 1-D analysis of the data. However, there is enough scope to improve the rendition of the results using 2-D processing schemes, based on the parameters of interest. e.g., we can perform colour based de-noising, followed by segmentation, on the colour image obtained by processing polarization data, to show the exact regions of different polarization orientations. Similarly, median filtering can be used to remove impulsive noise from the images. Segmentation followed by classification can be used to label different images of a scene. Such post-processing schemes can aid the users to pick their parameters of interest easily. We have not delved into such an analysis of post processing schemes, since most of the time, they are application dependent. 6.3.5 Harmonic based processing scheme In all the 1-D processing schemes that we have mentioned, we considered the ballistic component to show up at a single particular frequency. It was found to be the case, in CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH 134 almost all the data sets except two. In a data set taken with a mixture of 2.97µ and 0.06µ particles, we could observe a sinusoid at the first harmonic position of the expected frequency as well. Similarly, in a data set taken with the scattering medium consisting of a different concentration of the same particle, we could observe a sinusoid at the expected frequency, its first and also the second harmonics. We could observe that in the regions where the object blocked the ballistic components, the amplitude of the second harmonic was less, than in other regions. We used the ratio of the amplitude of the first harmonic to the fundamental, as the visualization parameter, after removing the outliers. It was observed that drastically better results could be obtained, than the PMI N results. The PMI N results also showed a strange behaviour, in the sense that, the shadow regions looked brighter than the background in the final result. PDI N result also showed similar behaviour. It could be observed that the data contained lot of speckle, as compared to other data sets. The object was also faintly visible before processing. However, none of these should give results as given by PMI N. The new processing scheme based on ratio of the amplitudes of the first harmonic to the fundamental, which we discussed, showed expected results. We call this scheme as the harmonic based processing scheme. Fig. 6.2(a) shows the image of the object without scattering. Fig. 6.2(b) shows a representative image of the series of images captured. Fig. 6.2(c) shows the PMI N PII result. Fig. 6.2(d) shows the result obtained with the harmonic based processing scheme, which we have proposed. We can clearly see the superiority of the result obtained with the proposed processing scheme, as compared to the PMI N result. However, we did not have more data to analyze the reasons behind the findings. It would be interesting to know as to under what conditions, the harmonic based processing scheme can come to use. CHAPTER 6. CONCLUSIONS AND TOPICS FOR FURTHER RESEARCH (a) (b) (c) (d) 135 Figure 6.2: Comparison of results of PMI N PII and the harmonic based processing scheme. The scattering medium contained 2.97µ and 0.06µ particles. (a) (b) (c) (d) Image of the object without scattering A representative of the 512 images grabbed PMI N PII result; N = 512 Result obtained with harmonic based processing Appendix A CRLB calculations - The case of white noise The equation governing the behaviour of data in polarization imaging is (eqn 4.5) Ir (n) = Iscat + Ibal sin  4πn +β M  + w(n) n = 0, 1, 2, . . . , N − 1 (A.1) where Iscat ≥ 0, Ibal ≥ 0 In the above equation, the unknown parameters of interest θ, are Iscat , Ibal and β. The variance of w(n) is unknown, but it is not of interest to us. The regularity conditions that the pdf p(Ir ; θ) has to satisfy for the bounds to exist, is given by,   ∂ ln p(Ir ; θ) E =0 ∂θ f or all θ where, N −1 1 1 X (Ir (n) − Iscat − Ibal sin γ)2 exp − p(Ir ; θ) = 2πσ 2 2σ 2 n=0 σ 2 is the noise variance of the time series, and γ = 136 4πN M + β. ! (A.2) APPENDIX A. CRLB CALCULATIONS - THE CASE OF WHITE NOISE 137 It can be easily verified that these conditions are satisfied and hence, we proceed to find the bounds for the parameters. Since multiple parameters are unknown, we find the bounds by finding the elements of the Fisher information matrix [I(θ)]ij which are given by [61]: M −1 1 X ∂s[n; θ] ∂s[n; θ] [I(θ)]ij = 2 σ n=0 ∂θi ∂θj (A.3) In this case, we have θ = [Iscat Ibal β]. Before finding the elements of the Fisher information matrix, we make the assumption that the number of data points available for analysis is an integral multiple of M , the period of the rotating polaroid. i.e., N = kM , where k is an integer. Moreover, we also assume that the periodicity of the rotation polaroid M itself is such that, the rotating polaroid makes exactly one rotation in M steps. However, these assumptions are needed only to arrive at easy expressions for the elements of the Fisher information matrix and the bounds derived will hold even otherwise, if N is quite large. Now, we calculate the elements of the Fisher information matrix, from first principles, using eqn A.3. [I(θ)]11 N N 1 X 1= 2 = 2 σ n=0 σ [I(θ)]12   N  4πn 1 X 1 × sin +β =0 = 2 σ n=0 M [I(θ)]13 [I(θ)]22   N  4πn 1 X 1 × Ibal cos +β =0 = 2 σ n=0 M   N  1 X N 4πn 2 = 2 1 × sin +β = 2 σ n=0 M 2σ APPENDIX A. CRLB CALCULATIONS - THE CASE OF WHITE NOISE [I(θ)]23 [I(θ)]33 138     N  1 X 4πn 4πn = 2 sin + β × Ibal cos +β =0 σ n=0 M M   N  4πn I2 N 1 X 2 2 +β = bal 2 = 2 Ibal cos σ n=0 M 2σ It can easily be verified that [I(θ)]21 = [I(θ)]12 [I(θ)]31 = [I(θ)]13 [I(θ)]32 = [I(θ)]23 Thus the Fisher information matrix is given by  N σ2   I(θ) =  0  0 0 0 N 2σ 2 0 0 2 N Ibal 2σ 2      (A.4) It is known from the Cramer-Rao Lower Bound theorem for vector parameters [61], that   var(θ̂i ) = [Cθ̂ ]ii ≥ I −1 (θ) ii (A.5) Hence, we obtain the bounds on the various parameters as σ2 var{Iscat } ≥ N (A.6) 2σ 2 var{Ibal } ≥ N (A.7) var{β} ≥ 2σ 2 1 = 2 Ibal N Nη where η is the SNR. (A.8) APPENDIX A. CRLB CALCULATIONS - THE CASE OF WHITE NOISE Bounds for estimating 139 Ibal Iscat In this section, we try to find the bound for estimates of the transformation d = Ibal DOLP Iscat (A.9) Assume that it is desired to estimate α = g(θ) for g, a r-dimensional function. From CRLB for transformation of vector parameters, it has been shown in [61], that Cα̂ − ∂g(θ) −1 ∂g(θ)T I (θ) ≥0 ∂θ ∂θ (A.10) In our case, g corresponds to DOLP and it is a one-dimensional transformation of d is related to the variance of Ibal parameters Ibal and Iscat . Hence, the variance of DOLP and Iscat as follows T d } ≥ ∂g(θ) I −1 (θ) ∂g(θ) var{DOLP ∂θ ∂θ (A.11) where, ∂g(θ) h = ∂θ d ∂ DOLP ∂Iscat d ∂ DOLP ∂Ibal d ∂ DOLP ∂β i = h − II2bal scat 1 Iscat 0 i (A.12) which leads to the relation, 2 2 2 d } ≥ σ Ibal + 2σ var{DOLP 4 2 N Iscat N Iscat (A.13) Appendix B CRLB calculations - The case of coloured noise The equation governing the behaviour of data in polarization imaging is (eqn 4.5) Ir (n) = Iscat + Ibal sin  4πn +β M  + v(n) n = 0, 1, 2, . . . , N − 1 (B.1) where Iscat ≥ 0, Ibal ≥ 0 and v(n) represents coloured noise. The values that we need to estimate are Iscat , Ibal and β. For analysis, we consider v(n) to be AR1 process with pdf N (0, C), where C is the noise covariance matrix. v(n) is given by v(n) = av(n − 1) + w(n) (B.2) where, a is the unknown AR1 co-efficient and w(n) is zero-mean Gaussian iid random variable with unknown variance σ 2 . Although we do not know the values of a and σ 2 , we fix these quantities in simulations. Once a and σ 2 are known, C can be calculated as follows [69]. 140 141 APPENDIX B. CRLB CALCULATIONS - THE CASE OF COLOURED NOISE We denote the autocorrelation values of the AR1 process by rv (n), n = 0, 1, 2, · · · σ2 rv (0) = 1 − a2 (B.3) and all other autocorrelation values can be found using the relationship rv (k) = rv (0)a|k| (B.4) Once the autocorrelation values are found, the noise covariance matrix of size p × p can be obtained by imposing a toeplitz structure on the matrix, with the values of the first row being the first p autocorrelation values. Another important relationship that has been used in our simulations, is the relationship between the power spectral density Px (ejω ) of the noise process and its variance, given by Px (ejω ) = σ2 | 1 + ae−jω |2 (B.5) Next, we derive the CRLB for Iscat , Ibal and β from first principles. The elements of the Fisher information matrix for the general Gaussian case are given by [61]  ∂µ(θ) [I(θ)]ij = ∂θi T C −1    1 ∂C(θ) −1 ∂C(θ) ∂µ(θ) −1 + tr C (θ) (θ) C (θ) ∂θj 2 ∂θi ∂θj  (B.6) where    ∂µ(θ)  =  ∂θi   ∂[µ(θ)]1 ∂θi ∂[µ(θ)]2 ∂θi .. . ∂[µ(θ)]N ∂θi         (B.7) For the problem on hand, the covariance matrix terms Cij do not depend on θ and hence the second term on the right hand side of the above equation becomes zero. The APPENDIX B. CRLB CALCULATIONS - THE CASE OF COLOURED NOISE 142 covariance matrix C (N ×N ) can be obtained by using the knowledge of the noise variance and the AR1 coefficient as explained earlier. With this observation, we find the elements of the Fisher information matrix, after rearranging eqn B.1 as follows, Ir (n) = Iscat + Ibal cos β sin     4πn 4πn + Ibal sin β cos + v(n) M M n = 0, 1, 2, . . . , N (B.8) The above equation can be easily represented in linear form as  Ir (0)   1 0 1     I (1)   1 sin 4π
cos 4π
 r   M M   

8π 8π  Ir (2)  =  1 sin cos    M M   .  .. .. .. .    . . .   . 

1 sin 4πN cos 4πN Ir (N ) M M {z | | {z } Ir H or, equivalently, by matrix notation as    I   scat    I cos β   bal   Ibal sin β  {z | Θ } Ir = HΘ + V   v(0)     v(1)        + v(2)       .  ..    } v(N ) | {z } (B.9) V (B.10) We can easily observe that the quantity corresponding to µ(θ) of eqn B.6 is, HΘ Hence h ∂µ(θ) =H ∂θi ∂Θ ∂θi i (B.11) Let P= h ∂Θ ∂θ1 ∂Θ ∂θ2 ∂Θ ∂θ3 i  1 0 0   =  0 cos β −Ibal sin β  0 sin β Ibal cos β      (B.12) APPENDIX B. CRLB CALCULATIONS - THE CASE OF COLOURED NOISE 143 Hence, from eqn B.6, the Fisher information matrix is given by [I(θ)] = (HP)T C −1 HP = PT HT C −1 HP (B.13) Let us denote the 3 × 3 matrix HT C −1 H by A. It can easily be observed that A will be a symmetric matrix. Moreover, due to the nature of the matrix H, we can easily verify that A22 = A33 . Thus, the general representation of A will be  a a a  11 12 13  A =  a12 a22 a23  a13 a23 a22      (B.14) Note that the elements of A can be theoretically calculated once C is known. By substituting for HT C −1 H and P in eqn B.13, we obtain the Fisher information matrix as    [I(θ)] =   a11 a12 cos β + a13 sin β Ibal (a13 cos β − a12 sin β) a12 cos β + a13 sin β a22 + a23 sin 2β a23 Ibal cos 2β Ibal (a13 cos β − a12 sin β) a23 Ibal cos 2β 2 Ibal (a22 − a23 sin 2β)      Now, we can obtain the bounds on the variances of Iscat , Ibal and β by inverting the Fisher information matrix. We do not find the inverse of the whole matrix; instead, we find only the necessary elements of the inverse of the matrix. The determinant of [I(θ)] can be obtained as the product of the determinants of PT ,A and P, due to the relation B.13. Hence, 2 det [I(θ)] = Ibal det(A) (B.15) APPENDIX B. CRLB CALCULATIONS - THE CASE OF COLOURED NOISE 144 We can see that  2 2 2  −1  a223 cos2 2β a11 Ibal (a22 − a223 sin2 2β) − Ibal I (θ) 11 = 2 Ibal det(A) a11 (a222 − a223 ) = det(A)  −1  a11 a22 + (a13 a12 − a11 a23 ) sin 2β − a213 cos2 β − a212 sin2 β I (θ) 22 = det(A)  −1  a11 a22 + (a11 a23 − a12 a13 ) sin 2β − a212 cos2 β − a213 sin2 β I (θ) 33 = 2 det(A) Ibal (B.16) (B.17) (B.18) The bounds on the various parameters are given by   var{Iˆscat } ≥ I −1 (θ) 11 (B.19)   var{Iˆbal } ≥ I −1 (θ) 22 (B.20)   var{β̂} ≥ I −1 (θ) 33 (B.21) From the above observations, we can immediately infer the following. • The CRLB for Iscat is independent of other parameters of interest. • The CRLB for Ibal is dependent on the phase of the sinusoid, which is itself an unknown parameter. • The CRLB for β is not only dependent on β, but in addition, it also depends on the amplitude of the sinusoidal component Ibal , similar to the case of CRLB for β in the case of white noise (see eqn A.8). Since the phase of the sinusoid itself is unknown in our case, we find the global minimum variance for the estimation of Ibal , by finding that β, which yields minimum variance of Ibal . This can be found to occur at 1 β = arctan 2  2(a12 a13 − a11 a23 ) a212 − a213  (B.22) APPENDIX B. CRLB CALCULATIONS - THE CASE OF COLOURED NOISE 145 It can be verified that the same value of β yields the minimum variance of β̂ too. As we can see, the value of β that yields minimum variance of Iˆbal and β̂ depends on the noise parameters. Bounds for estimating Ibal Iscat The procedure for estimating the bounds are similar to the analysis carried out in the case of white noise. From eqn A.12, we have the required bound to be b var{DOLP } ≥ T(PT HT C−1 HP)−1 TT (B.23) where T= h − II2bal scat 1 Iscat 0 i (B.24) Bibliography [1] N. S. Kopeika, A System Engineering Approach to Imaging. Bellingham, Washington, USA: SPIE Optical Engineering Press, 1998. [2] R. Hema, “Imaging through turbid media,” CURRENT SCIENCE, vol. 76, pp. 1334– 1340, May 1999. [3] J. S. Tyo, M. P. Rowe, E. N. PughJr, and N. Engheta, “Target detection in optically scattering media by polarization-difference imaging,” Applied Optics, vol. 35, pp. 1855–1870, Apr. 1996. [4] Y. Y. Schechner, S. G. Narasimhan, and S. K. Nayar, “Polarization-based vision through haze,” Applied Optics, vol. 42, pp. 511–525, Jan. 2003. [5] J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. experimental techniques,” Phys.Med.Biol., vol. 42, pp. 825–840, 1997. [6] C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles. New York: John Wiley and Sons Inc, 1998. [7] S. Huard, Polarization of light. John Wiley & Sons, 1997. [8] W. S. Bickel and W. M. Bailey, “Stokes vectors, mueller matrices, and polarized scattered light,” American Journal of Physics, vol. 53, May 1985. [9] E. P. Zege, A. P. Ivanov, and I. L. Katsev, Image Transfer Through a Scattering Medium. Berlin Heidelberg: Springer-Verlag, 1991. 146 BIBLIOGRAPHY 147 [10] vande Hulst H C, MULTIPLE LIGHT SCATTERING Tables Formulas and Applications. New York USA: Academic Press, 1980. [11] van Rossum M C W and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy mesoscopy and diffusion,” Reviews of Modern Physics, vol. 71, Jan. 1999. [12] S. P. Morgan, M. P. Knong, and M. G. Somekh, “Effects of polarization state and scatterer concentration on optical imaging through scattering media,” Applied Optics, vol. 36, pp. 1560–1565, Mar. 1997. [13] Journal of Optical Society of America, A, vol. 14, no. 1, Jan. 1997. [14] Applied Optices, vol. 36, no. 1, Jan. 1997. [15] Applied Optices, vol. 36, no. 2, Jan. 1997. [16] Journal of Biomedical Optics, vol. 7, no. 3, July 2002. [17] “Mathematics and physics of emerging biomedical imaging.” Committee on the Mathematics and Physics of Emerging Dynamic Biomedical Imaging- National Research Council, 1996. [18] D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Processing Magazine, pp. 57–75, Nov. 2001. [19] R. L. Barbour and H. L. Graber, “Diagnostic imaging with light, and beyond,” in Proceedings-19th International conference-IEEE/EMBS, pp. 2735–2740, 1997. [20] H. Jiang, N. V. Iftimia, Y. Xu, J. A. Eggert, L. L. Fajardo, and K. L. Klove, “Nearinfrared optical imaging of the breast with model-based reconstruction,” Academic Radiology, vol. 9, no. 2, pp. 186–194, 2002. [21] J. C. Hebden, A. Gibson, R. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical BIBLIOGRAPHY 148 tomography of the premature infant brain,” Physics in Medicine and Biology, vol. 47, pp. 4155–4166, 2002. [22] J. C. Hebden and D. T. Delpy, “Diagnostic imaging with light,” The British Journal of Radiology, vol. 70, pp. S206–S214, 1997. [23] L. R. Bissonette, “Imaging through fog and rain,” Optical Engineering, vol. 31, pp. 1045–1052, May 1992. [24] F. M. Caimi, “A review of recent underwater imaging methods and advancements,” in Coastal Ocean: Prospects for the 21st Century, vol. 1, pp. 69–71, IEEE, 1996. [25] I. Freund, “Image reconstruction through multiple scattering media,” Optics Communications, vol. 86, pp. 216–227, Nov. 1991. [26] I. Freund, “Correlation imaging through multiply scattering media,” PHYSICS LETTERS A, vol. 147, pp. 502–506, July 1990. [27] R. Hema and N. Andal, “Two-dimensional imaging through turbid media using a continuous wave light source,” Optics Communications, vol. 154, pp. 255–260, Sept. 1998. [28] M. P. Rowe, E. N. PughJr, J. S. Tyo, and N. Engheta, “Polarization-difference imaging: a biologically inspired technique for observation through scattering media,” OPTICS LETTERS, vol. 20, pp. 608–610, Mar. 1995. [29] S. G. Demos and R. R. Alfano, “Temporal gating in highly scattering media by the degree of optical polarization,” OPTICS LETTERS, vol. 21, p. 161, Feb. 1996. [30] S. G. Demos and R. R. Alfano, “Optical polarization imaging,” Applied Optics, vol. 36, pp. 150–155, Jan. 1997. [31] O. Emile, F. Bretenaker, and A. L. Floch, “Rotating polarization imaging in turbid media,” OPTICS LETTERS, vol. 21, pp. 1706–1708, Oct. 1996. [32] H. Horinaka, K. Hashimoto, K. Wada, and Y. Cho, “Extraction of quasi straightforward propagating photons from diffused light transmission through a scattering BIBLIOGRAPHY 149 medium by polarization modulation,” OPTICS LETTERS, vol. 20, pp. 1501–1503, July 1995. [33] H. Horinaka, K. Hashimoto, K. Wada, T. Umeda, and Y. Cho, “Optical ct imaging in highly scattering media by extraction of photons preserving initial polarization,” in International Symposium on Polarization Analysis and Applications to Device Technology (T. Yoshizawa and H. Yokota, eds.), vol. 2873, pp. 54–57, SPIE, 1996. [34] S. P. Morgan, M. P. Khong, and M. G. Somekh, “Polarization imaging through scattering media,” in Photon Propagation in Tissues (B. Chance, D. T. Delpy, and G. J. Mueller, eds.), vol. 2626, pp. 265–272, SPIE, 1995. [35] G. Jarry, E. Steimer, V. Damaschini, M. Epifanie, M. Jurczak, and R. Kaiser, “Coherence and polarization of light propagating through scattering media and biological tissues,” Applied Optics, vol. 37, pp. 7357–7367, Nov. 1998. [36] J. M. Schmitt, “Optical coherence tomography: A review,” IEEE JOURNAL OF SELECTED TOPICS IN QUANTUM ELECTRONICS, vol. 5, pp. 1205–1215, Aug. 1999. [37] D. A. Zimnyakov, V. V. Tuchin, R. A. Zdrajevsky, and Y. P. Sinichkin, “Coherent and polarization imaging: Novel approaches in tissue diagnostics by laser light scattering,” in Laser tissue interaction XI: Photochemical, Photothermal and Photomechanical (D. D. Duncan and J. O. Hollinger, eds.), vol. 3914, pp. 407–422, SPIE, 2000. [38] C. Brosseau, Fundamentals of polarized light-A statistical optics approach. New York: John Wiley and Sons Inc, 1998. [39] D. A. Zimnyakov and Y. P. Sinichkin, “A study of polarization decay as applied to improved imaging in scattering media,” Journal of Optics A:Pure Applied Optics, vol. 2, pp. 200–208, Feb. 2000. BIBLIOGRAPHY 150 [40] D. Bicout, C. Brossear, A. S. Martinez, and J. M. Schmitt, “Depolarization of multiply scattered waves by spherical diffusers: Influence of the size parameter,” PHYSICAL REVIEW E, vol. 49, pp. 1767–1770, Feb. 1994. [41] M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissues,” J. Opt. Soc. Am. A, vol. 18, pp. 948–960, Apr. 2001. [42] J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Applied Optics, vol. 31, pp. 6535–6546, Oct. 1992. [43] K. K. Bizheva, A. M. Siegel, and D. A. Boas, “Path-length-resolved dynamic light scattering in highly scattering random media: The transition to diffusing wave spectroscopy,” PHYSICAL REVIEW E, vol. 58, pp. 7664–7667, Dec. 1998. [44] A. D. Kim and M. Moscoso, “Influence of the relative refractive index on the depolarization of multiply scattered waves,” PHYSICAL REVIEW E, vol. 64, July 2001. [45] V. Shankaran, M. J. Everett, D. J. Maitland, and J. T. W. Jr., “Comparison of polarized-light propagation in biological tissue and phantoms,” OPTICS LETTERS, vol. 24, pp. 1044–1046, Aug. 1999. [46] V. Shankaran, J. T. W. Jr., and D. J. Maitland, “Polarized light propagation through tissue phantoms containing densely packed scatterers,” OPTICS LETTERS, vol. 25, pp. 239–241, Feb. 2000. [47] Q. Z. Wang, X. Liang, L. Wang, P. P. Ho, and R. R. Alfano, “Fourier spatial filter acts as a temporal gate for light propagating through turbid media,” OPTICS LETTERS, vol. 20, pp. 1498–1500, July 1995. [48] X. Gan, S. P. Schilders, and M. Gu, “Image enhancement through turbid media under a microscope by use of polarization gating methods,” J. Opt. Soc. Am. A, vol. 16, pp. 2177–2184, Sept. 1999. BIBLIOGRAPHY 151 [49] M. Kempe, A. Z. Genack, W. Rudolph, and P. Dorn, “Ballistic and diffuse light detection in confocal and heterodyne imaging systems,” J. Opt. Soc. Am. A, vol. 14, pp. 216–223, Jan. 1997. [50] S. L. Jacques and K. Lee, “Imaging tissues with a polarized light video camera,” in 1999 International Conference on Biomedical Optics (Q. Luo, B. Chance, L. V. Wang, and S. L. Jacques, eds.), vol. 3863, pp. 68–74, SPIE, 1999. [51] S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” Journal of Biomedical Optics, vol. 7, pp. 329–340, July 2002. [52] K. Turpin, J. G. Walker, P. C. Y. Chang, K. I. Hopcraft, B. Ablitt, and E. Jakeman, “The influence of particle size in active polarization imaging in scattering media,” Optics Communications, vol. 168, pp. 325–335, Sept. 1999. [53] J. G. Walker, P. C. Y. Chang, and K. I. Hopcraft, “Visibility depth improvement in active polarization imaging in scattering media,” Applied Optics, vol. 39, pp. 4993– 4941, Sept. 2000. [54] P. Chang, J. Walker, K. Hopcraft, B. Ablitt, and E. Jakeman, “Polarization discrimination for active imaging in scattering media,” Optics Communications, vol. 159, pp. 1–6, Jan. 1999. [55] Y. P. Sinichkin, D. A. Zimnyakov, and V. V. Giterman, “Direct polarization imaging of turbid tissues with cw laser source: Potentialities and restrictions,” in Coherence Domain Optical Methods in Biomedical Science and Clinical Applications III (V. V. Tuchin and J. A. Izatt, eds.), vol. 3598, pp. 258–262, SPIE, 1999. [56] G. Venkatesh, M. Sushil, R. Hema, and A. K. Sood, “Imaging in turbid media using quasi-ballistic photons,” Optics Communications, vol. 170, pp. 331–345, Nov. 1999. [57] R. A. Chipman, “Polarization diversity active imaging,” in Image Reconstruction and Restoration II (T. J. Schulz, ed.), vol. 3170, pp. 68–73, SPIE, 1997. [58] J. L. Pezzantini and R. A. Chipman, “Mueller matrix imaging polarimetry,” Optical Engineering, vol. 34, p. 1558, June 1995. BIBLIOGRAPHY 152 [59] P. Clemenceau, S. Breugnot, and L. Collot, “Polarization diversity active imaging,” in Laser Radar Technology and Applications III (G. W. Kamerman, ed.), vol. 3380, pp. 284–290, SPIE, 1998. [60] M. H. Smith, A. Lompado, and P. Burke, “Mueller matrix imaging polarimetry in dermatology,” in Biomedical Diagnostic, Guidance and Surgical Assist-Systems II (Vo-Dhin, Grundfest, and B. D. A, eds.), vol. 3911, SPIE, 2000. [61] S. M. kay, Fundamentals of Statistical Signal Processing: Estimation Theory. New Jersey: Prentice Hall International Inc., 1993. [62] P. Stoica, H. Li, and J. Li, “Amplitude estimation of sinusoidal signals: Survey, new results, and an application,” IEEE TRANSACTION ON SIGNAL PROCESSING, vol. 48, pp. 338–352, Feb. 2000. [63] J. Li and P. Stoica, “An adaptive filtering approach to spectral estimation and sar imaging,” IEEE TRANSACTIONS ON SIGNAL PROCESSING, vol. 44, pp. 1469– 1484, June 1996. [64] P. Stoica, A. Jakobsson, and J. Li, “Matched-filter bank interpretation of some spectral estimators,” Signal Processing, vol. 66, pp. 45–59, June 1998. [65] A. M. Zoubir and B. Boashash, “The bootstrap and its application in signal processing,” IEEE Signal Processing Magazine, pp. 56–76, 1998. [66] H.-T. Ong and A. M. Zoubir, “Bootstrap-based detection of signals with unknown parameters in unspecified correlated interference,” IEEE Transaction on Signal Processing, vol. 51, pp. 135–141, Jan. 2003. [67] H. V. Sorensen, C. S. Burrus, and D. L. Jones, “A new efficient algorithm for computing few dft points,” in IEEE international symposium on circuits and systems, vol. 2, pp. 1915–1918, June 1988. [68] Z.-S. Liu, H. Li, and J. Li, “Efficient implementation of capon and apes for spectral estimation,” IEEE Transactions on Aerospace and Electronic Systems, vol. 34, pp. 1314–1319, Oct. 1998. BIBLIOGRAPHY 153 [69] M. H. Hayes, Statistical Digital Signal Processing and Modeling. New York: John Wiley and Sons Inc., 1996. [70] G. Marsaglia, “Ratios of normal variables and ratios of sums of uniform variables,” American statistical association journal, vol. 60, pp. 193–204, Mar. 1965. [71] G. M. P. van Kempen and L. J. van Vliet, “Mean and variance of ratio estimators used in fluorescence ratio imaging,” Cytometry, vol. 39, pp. 300–305, 2000. [72] S. Chandrashekar, Radiative transfer. New York: Dover Publications Inc, 1960.