Comparison and validation of point spread models for imaging in natural waters

It is known that scattering by particulates within natural waters is the main cause of the blur in underwater images. Underwater images can be better restored or enhanced with knowledge of the point spread function (PSF) of the water. This will extend the performance range as well as the information retrieval from underwater electro-optical systems, which is critical in many civilian and military applications, including target and especially mine detection, search and rescue, and diver visibility. A better understanding of the physical process involved also helps to predict system performance and simulate it accurately on demand. The presented effort first reviews several PSF models, including the introduction of a semianalytical PSF given optical properties of the medium, including scattering albedo, mean scattering angles and the optical range. The models under comparison include the empirical model of Duntley, a modified PSF model by Dolin el al, as well as the numerical integration of analytical forms from Wells, as a benchmark of theoretical results. For experimental results, in addition to that of Duntley, we validate the above models with measured point spread functions by applying field measured scattering properties with Monte Carlo simulations. Results from these comparisons suggest it is sufficient but necessary to have the three parameters listed above to model PSFs. The simplified approach introduced also provides adequate accuracy and flexibility for imaging applications, as shown by examples of restored underwater images. 02008 Optical Society of America OCIS codes: (010.4450) Ocean optics; (110.4850) OTF; References and links I. G. D. Gilbert and R. C. Honey, "Optical turbulence in the sea," in Underwater photo-optical instrumentation applications (SPIE, 1972), pp. 49-55. 2. D. J. Bogucki, J. A. Domaradzki, R. E. Ecke, and C. R. Truman, "Light scattering on oceanic turbulence," Appl. Opt. 43, 5662-5668 (2004). 3. C. D. Mobley. Light and Water: radiative transfer in natural waters (Academic Press, New York, 1994). 4. D. Gray, "Order-of-Scattering Point Spread and Modulation Transfer Functions for Natural Waters," to bc submitted to OE (2008). 5. 1. L. Katsev, E. P. Zege, A. S. Prikhach, and I. N. Polonsky, "Efficient technique to determine backscattered light power for various atmospheric and oceanic sounding and imaging systems," J. Opt. Soc. Am. 14, 13381346 (1997). 6. W. Hou, D. Gray, A. Weidemann, G. R. Fournier, and J. L Forand, "Automated underwater image restoration and retrieval of related optical properties," in IEEE International Geoscience and Remote Sensing Symposium, IGARSS (2007), pp. 1889-1892. 7. R. C. Gonzalez, and R. E. Woods, Digital Image Processing (Prentice Hall, Upper Saddle River, N.J., 2002). 8. L. E. Mertens, and J. F. S. Replogle, "Use of point spread and beam spread functions for analysis of imaging systems in water," J. Opt. Soc. Am. 67, 1105-1117 (1977). 9. J. Zhang, Q. Zhang, and G. He, "Blind deconvolution: multiplicative iterative algorithm," Opt. Lett. 33, 2527(2008). 10. S. Q. Duntley, "Underwater lighting by submerged lasers and incandescent sources," (Scripts Instituition of Oceanography, University of California, San Diego, 1971). II. K. Voss, "Simple empirical model of the oceanic point spread function," Appl. Opt. 30. 2647-2651 (1991). #94847 $15.00 USD Received 8 Apr 2008; revised 30 May 2008; accepted 12 Jun 2008; published 20 Jun 2008 (C) 2008 OSA 23 June 2008 / Vol. 16, No. 13 / OPTICS EXPRESS 9958


Introduction
Image formation is usually affected by both the optical system design including lenses and recording devices, as well as transfer processes associated within the formation medium. Applications to remove the effects of the medium have largely been done in atmospheric situations, where astronomical and reconnaissance needs demand a better understanding and solution. This typically involves the effects due to atmospheric turbulence. In natural waters, although the turbulence does impose its limitations under certain conditions, especially at longer range[l, 2], the largest contribution to degraded images comes from suspended particulates. Their size and index of refraction changes cause blurring in the resulting imagery, thus limiting visibility range.
In general, scattering properties of the medium determine the outcome of the image transmission. In ocean and lake environment optics, the scattering properties are conveniently described and measured by the scattering coefficient (b), which determines the probability that a photon will be scattered away from its original traveling direction per unit length by the medium molecules, constituents within the medium (ie particles [3]), and turbulence [2]. The scattering parameter (b) is an integration of the volume scattering or phase function, i, which details such probabilities by the relative directions of incoming and out-going photons [3]. These inherent optical properties (1OP), although measured frequently due to their important applications in ocean optics, especially in remote sensing, have not been applied to underwater imaging issues directly, in part because they inherently only reflect the effect of the single scattering.
A more intuitive and direct measure in imaging is the point spread function (PSF), which gives the system response to a point source, and thus includes the effect of multiple scattering [3, 4]. It is the ideal parameter to study image transmission, optical sounding [5], and retrieval of optical properties [6]. Simply put, an image of an object is the combination of the original signal, f(x,y), convolved with the entire system point spread function h(x,y), adding noise n(x,y), The system response includes those from both the imaging system itself, as well as the effects of the medium (water in our case). With known characteristics of the imaging system and correct modeling of the medium, theoretically it is possible to fully recover the original signal by reversion or deconvolution [7]. Mathematically PSF is equivalent to the beam spread function (BSF) [8] which can be modeled and measured more easily. From (1), one can see that with the knowledge of g and h, the original imagef can be restored, especially when the noise can be reduced with improved system setup or optimization[9). Additionally, with the knowledge of the PSF, one can easily simulate the system outcome. This is valuable to system design, performance prediction, as well as underwater scene simulations, as it is based on the accuracy of the physical model of the system (medium included).

Selected previous PSF models
Although established almost 40 years ago, the most thorough empirical point spread model to date is still that of Duntley [10]. Duntley reported extensive lab measurements of BSFs of simulated ocean waters with remarkably different optical lengths (0.5 to 21), and summarized his findings in a single, albeit complex empirical relationship: E D =(0.018-+ 0.011'+0.001725r)r. Note in the above equations that the optical length is defined as r = cr (c being the total attenuation coefficient) while the single scattering albedo is given as q = b/c. " = 1/( -o) is used to simplify the form. H is the spectral irradiance and is normalized by the laser power P, which gives the BSF. To add to the complexity, one should realize that the above regression results are in degrees, rather than radians. It is of little surprise that the original publication had to resolve to errata to get the exact form correct, while another citation also contained misinformation likely due to printing errors [ 11].
1 -0.072r +6.3.10 -3 . r, +9.4.10 -. r ;, Notice that a dimensionless scaling factor q is used, which can be shown to approximate the mean scattering angle [12]. Po and g relate to the direct beam initial power and can be omitted when normalizations are applied to scattered photons. The above simplified result is under the assumption that the absorption coefficient a>>bb, the backscattering coefficient, and b>>bb, such that b±2bb-b, and a±2bb=a. This holds true in most conditions applied and introduces little variation when tested. This is not surprising as our focus here is on small forward scattering angles. Once again, the complexity of this model is evident and even exceeds that of Duntley's. However, it does offer a better approximation, as we shall examine later.
Lastly, using field collected data, Voss put forth a simple-looking empirical PSF expression [11], that fits 3 different types of oceanic waters (TOTO, Pacific and Sargasso Sea). The PSF covered the range between 4 and 100 mrad and was given by: where -m is the slope of log(PSF)log(6), and BI is a constant. Here m is not a constant, but rather a function of IOP and r The above formula fits all data to less than 15% error. For imaging restoration applications, especially considering optimization in convolution[9], this is an acceptable start. The m values ranged from 0.4 to 2.0 as a function of r (0 to 10). Other than graphical results, there is no definitive relationship provided to allow comparison of this method to other approaches. Nonetheless, this is a very encouraging result, since such simple relationships can be rather beneficial to imaging needs, especially when high frame rate per pixel calculations are needed. The aim of this paper is to explore this form further by examining the PSF relationship to commonly measured lOPs in search of an even simpler solution for real-time imaging needs.

A simplified PSF model
A simple form of the PSF is highly desirable, especially when routinely measured optical properties can be used as input. We introduced such a model previously [ 13], based on the analytical form by Gordon [ 14]. For brevity only key steps are outlined here. We first note that the PSF can be approximated by examining the amount of non-disturbed light variation within a cone of half angle 0 at range r [ 14,15], therefore 0 Direct integration of (14) can be carried out for a known analytical volume scattering function. The parameter t is used for the convenience of integration and can be interpreted as the portion of the scattering angle up to 6. There are many possible and defendable choices of /8 of natural waters. Here we take a simple formulation from Wells [ 15], by using a phase function in the following form: where 0o relates to the mean scattering angle. Substituting (15) into (14) for integration will give inverse power relationships to the scattering angle 6. With the help of numerical integrations, we arrive at the following empirical form: where K is a constant and does not affect imaging needs since relative units of PSF are used, and m=l/o)-20%. The later term is necessary to include higher scattering orders. This is inline with numerical calculations over different parameter ranges. Comparisons to empirical and measured formations can now be made.

Comparison and validation of PSF models
We proceed by comparing the following PSF models: 1) our simplified model; 2) Duntley's empirical relationship from simulated seawaters; 3) a modified Dolin's PSF model from analytical results; 4) numerically calculated PSF using Wells analytical method; and 5) a PSF derived from field scattering measurements and applying Monte Carlo (MC) methods. We see that the second and third methods are straight-forward. The fourth approach, found in other publications [16,17], offers a means to use ground-truth, to accurately simulate PSFs from optical scattering properties. The details of this approach can be found in a related paper [4], but a brief description is offered below. This method uses field measured lOPs including volume scattering functions to numerically calculate BSF, thus PSF, using MC. As a proven technique, MC is essentially a numerical experiment of the physical process under consideration, thus providing accurate results although it does require extensive computation as the results are statistical in nature and each set of parameters must be run separately. In fact, the MC method has been used previously to test the range dependence of the Wells theory [18]. Given the lOPs (phase function, absorption and scattering coefficients), the radius of the sphere R, and an angular resolution A, the code produces the irradiance for each angular bin over the entire sphere.
Data were obtained during an April-May 2006 NATO trial experiment in Panama City, Florida. The amount of scattering and absorption were controlled by introducing Maalox and absorption dye respectively into a tap-water filled pool. In-water optical properties during the experiment were measured and included the absorption and attenuation coefficients (Wetlabs ac-9), particle size distributions (Sequoia Scientific LISST-100), and volume scattering functions (multi-spectral volume scattering meter or MVSM).  Figs. I and 2, and parameters listed in Table 1. These are normalized to 100mrad values. The current approach agrees well with measurements as well as empirical models under the different testing conditions shown. There are only small differences except at larger scattering angles with the MVSM, Wells and Dolin. Different optical range and scattering conditions are compared in Fig. 1, under similar scattering characteristics (Oo=0.03), which would correspond to q=7.6 in Dolin's model. We can see that Dolin and Wells models match those of measured MVSM very well under all conditions, which is not surprising as they applied few approximations during derivation. At larger angles, there are visible differences between the current model results and those of Wells and Dolin, however, it still approximates Duntley's empirical results nicely. These differences only have small effects on image restoration[6], because of the relatively small effects (log scale) in grayscale pixel values (ie, I or 2 out of 255, occurring several pixels away from the center of PSF). Initial restoration testing with images obtained during the NATO experiment show little differences using any of these PSFs under typical conditions tested [ 17]. If simplicity or computation time is of less concern, the differences with the current model at larger angles can be effectively reduced by applying a "damping factor", such as exp(-50 1.) to the right side of Eq. (16), or apply Wells' or Dolin's model directly. One may notice from the figure that PSFs from MVSM are almost always higher than Duntley's. This is likely due to the differences in the scattering agents involved, as the particle sizes will be different because of the different preparation process (recall that Duntley used high purity filtered water as their experiment base in the tank, before adding similar absorption and scattering agents). In fact, the differences in PSFs can be as high as an order of magnitude under certain conditions at small angles (eg. a=0.87, r-5.5), or ±1 in PSF log slope. Examination of the parameter related to o suggests that a higher value (0.06) is more appropriate in order to match Duntley's PSF, corresponding to less forward scattering contribution in Duntley's case (Fig. 2). This is not surprising since despite the complexity of Duntley's model, it depends on only two parameters, namely the optical range and scattering albedo. It lacks the ability to alter the shape of the volume scattering function. Variations by 0o (Fig. 2) demonstrated the flexibility of current model, and also suggest that the reason behind Voss' choice of 66=0.06 [11] may likely be the result of different particle types thus the shape of scattering functions involved. Also note that the PSF from the Wells model is more responsive than that of Dolin. Also, it is worth pointing out, that the current approach matches those from Duntley's well to at least z-15, but with much less computation required. The above derived PSF model (Eq. 16) was further validated when applied to restore blurred underwater images. A Secchi disk image taken by Navy diver's Nikon Coolpix digital camera at 5 ft away during an exercise to test underwater diver visibility models (August 2 1, 2002, Mississippi Bight, c=l.5m 1 ) was used. When the correct PSF is applied, the original image showed significant improvements, both in details and contrast [ Fig. 3(b)]. When the same formula is used for the PSF but assuming a longer optical length (r=10 instead of 2.3), we notice actual degradation in details due to over-correction, although contrast was still improved. Both restorations were done using an iterative Lucy-Richardson method [7]. When a blind deconvolution approach [7] was tested with an intentionally different PSF (Gaussian in this case), the result was apparently worse, even though traditionally (ie in astronomy applications) such algorithms are not sensitive to the shape of PSFs

Conclusion
Several types of PSF models of natural waters including empirical, analytical and measurement results under different conditions are compared, along with the introduction of a simplified semi-analytical form which shows a close approximation to the above methods. Measured PSFs were obtained using MVSM-measured phase functions and MC simulations. The results show the PSFs from a modified Dolin and Wells are very accurate even at larger angles and longer optical ranges. It is shown that three parameters are needed to adequately model the system response function of the natural waters. This explains in part the inflexibility of Duntley's empirical relationship. Although small differences exist between the simplified model and other approaches, they are not significant as far as image restoration is concerned. The simpler form of the model can be beneficial to various imaging needs in natural environments, especially to real-time processing and scene generation applications. The fact that different scattering agents can result in an order of magnitude difference in PSFs between measured MC-MVSM and Duntley's at small angles, which affects the slope of the PSF significantly, suggests additional studies are needed to better determine the effects of particle types on the shape of scattering functions.