1. 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[1, 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, β, which details such probabilities by the relative directions of incoming and out-going photons[3]. These inherent optical properties (IOP), 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 image f 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).

2. 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:

where parameters B=1-2(10^-D), A=1.260-0.375τ ^{0.710+0.489/ς}-(1.378+0.053ς)10^{-(0.268+0.088)τ}, E=(13.75-0.501ς)-(0.626-0.357ς)τ+(0.01258+0.00354ς)τ ², $C=\frac{1}{3}\left[1-{\left(\frac{{\theta }^{\frac{3}{2}}+E}{E}\right)}^{\frac{2}{3}}\right]$ , and D=(0.018+0.011ς+0.001725τ)τ. Note in the above equations that the optical length is defined as τ=cr (c being the total attenuation coefficient) while the single scattering albedo is given as ω ₀=b/c·ς=1/(1-ω ₀) 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].

A different approach was used by Dolin et al [12], where numerical approximations were applied in his analytical formulation of the PSF. We simplify this model by assuming little backscattering compared to the total scattering. Briefly, a modified BSF E(θ,r) can be expressed via PSF kernel G, in the order of direct beam, single and multiple scattering terms as:

Notice that a dimensionless scaling factor q is used, which can be shown to approximate the mean scattering angle [12]. P₀ and δ 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≫b_b, the backscattering coefficient, and b≫b_b, such that b±2b_b≈b, and a±2b_b≈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(θ), and B₁ is a constant. Here m is not a constant, but rather a function of IOP and τ. 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 τ (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 IOPs in search of an even simpler solution for real-time imaging needs.

3. 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 θ at range r [14, 15], therefore

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 θ. There are many possible and defendable choices of β of natural waters. Here we take a simple formulation from Wells [15], by using a phase function in the following form:

where θ₀ relates to the mean scattering angle. Substituting (15) into (14) for integration will give inverse power relationships to the scattering angle θ. 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=1/ω₀-2τθ₀. 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.

4. 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 IOPs 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 IOPs (phase function, absorption and scattering coefficients), the radius of the sphere R, and an angular resolution Δθ, 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).

Table 1. Optical parameters used to compare different PSF models as shown in Figs. 1–2.

View Table

A comparison of PSFs from different methods outlined, namely, the current study, MC simulated PSF from measured MVSM phase functions, Wells’ numerical results, Dolin’s model, and Duntley’s are shown in Figs. 1 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 (θ₀=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.

 

Fig. 1. Comparison of PSFs from different methods from a-f: “Current” (current study), “Duntley” and “Dolin” (empirical relationship in text), and “MVSM” (Monte Carlo simulated PSF based on measured volume scattering function during NATO experiment) and Wells numerical integrated PSF. Normalized to 100 mrad. See text for details.

Download Full Size | PDF

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, 1 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(-5θ^1.2) 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. ω₀=0.87, τ=5.5), or ±1 in PSF log slope. Examination of the parameter related to θ₀ 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 θ₀ (Fig. 2) demonstrated the flexibility of current model, and also suggest that the reason behind Voss’ choice of θ₀=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 τ=15, but with much less computation required.

 

Fig. 2. An example of PSF comparison under different θ₀ values. Notice that the variation of θ₀ does not affect Duntley or MVSM and all rest parameters of other models are the same as in Fig.1. Normalized to 100 mrad.

Download Full Size | PDF

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 21, 2002, Mississippi Bight, c=1.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 (τ=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 [7, 9]. These examples further demonstrate the importance and validity of the correct PSF model to underwater image restorations.

 

Fig. 3. A diver captured image and restorations with different PSFs. Image (top left, 3a) was taken during August 21, 2002 diving at Mississippi Bight with total attenuation c=1.5m^-1 at 5ft distance. 3b (top right) shows the restored image using the correct PSF model. 3c (bottom left) shows the image restored assuming a longer distance. 3d (bottom right) shows the result of a blind deconvolution with a Gaussian PSF.

Download Full Size | PDF

5. 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.