Modelling of beam propagation and its applications for underwater imaging

Abstract

In order to process underwater imaging to the best possible level, an imaging model based on beam propagation was established. The presented model included not only the laser beam propagation affected by absorption and scattering, but also the effects of underwater turbulence and the diffraction limit of sensors. By this model approximately quantified optical transfer functions (OTFs) were studied. Thus, under this framework, the approaches of image enhancement, restoration and super-resolution reconstruction (SRR) can be extended by incorporating underwater optical properties based on OTF or point spread function (PSF) of the imaging system. Experimental results proved that the imaging range and the image quality can be effectively enhanced, which are critical in underwater imaging or detecting.

Key words： image super-resolution reconstruction (SRR); optical transfer function (OTF); point spread function (PSF); range-gated

Introduction

Underwater imaging has not been a new concept since the year 1963 when an underwater transmission window of 470-580 nm, in which the attenuation of blue-green laser is far less than those other light-wave was found by Duntley [1]. The range and quality of underwater imaging are of vital importance to military and civilian applications. However, factors such as absorption, scattering and turbulence significantly reduce the light arrived at sensors and cause blurring and distortion in resulting images. Over past years, efforts have been taken to overcome these limitations, such as laser scanning imaging [2], polarized laser imaging [3] and range-gated imaging [4] which can effectively eliminate backscatter. Appropriately increasing the intensity of laser, enhancing detection rates of sensors and reducing error rate are also effective ways, however, these will undoubtedly increase the system cost greatly.

Therefore, in order to further improve the range and quality of underwater imaging beyond the hardware limitations, digital image process is needed. Image denoising, enhancement, restoration and super-resolution reconstruction (SRR) techniques are widely used. The main challenge work with digital image processing is the choice of accurate methods. Information based on the knowledge of system [5] can effectively enhance the performance of image process such as restoration. The system response which is a typical prior knowledge can be derived from underwater optical theory or image itself by measuring resolution chart at different distances. The modelling for underwater beam propagation [6-8] based on underwater optical theory can derive the system response such as point spread function (PSF). It has been studied since 1970’s when several research groups [6] suggested that linear system theory can be applied to underwater beam transmission. But most research works use traditional PSF models which cannot fully suit particular situations such as range-gated imaging.

Based on linear transmission theory, an underwater imaging model suitable for range-gated imaging system was established in this study including laser beam propagation affected by absorption and scattering, and the effects of underwater turbulence and the diffraction limit of sensors. The model-derived modulation transfer function (MTF) and PSF are applied for the theoretical basis of underwater image processing.

Underwater imaging model

The purpose of the underwater imaging model is to predict the image intensity at each pixel as a function of illumination, reflectance properties of objects, medium, and sensor characteristics. The main difference between image formation in air and in water is that the interaction between light and medium must be considered in underwater conditions unlike similar situations in the air.

Underwater laser beam transmission

Generally speaking, a typical underwater imaging system consists of three parts: the light source (laser), the image acquisition system (CCD or Intensified-CCD (ICCD)), and the object. Light on the image plane of sensor is formed by three components including nonscattered or direct light reflected by the object, the forward-scattered light from the object, and the back-scattered light from the medium (water). So, the total irradiance on image plane is generated by the linear plus of direct, forward-scattered, and back-scattered components:

(1)E_{\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{v}\mathrm{e}}(\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l})=E_{\mathrm{d}}(\mathrm{d}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t})+E_{\mathrm{f}s}(\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{d})+E_{\mathrm{b}\mathrm{s}}(\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{d}).

Fig.1 Light paths of three components

Full size|PPT slide

Figure 1 shows the light path of the three components. The direct light component received by the image sensor is the exponential decay of that reflected by the object with distance due to the attenuation of water:

(2)E_{\mathrm{d}}=E_{\mathrm{o}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}}\exp (-kl),

where l is imaging distance between the transmitting plane and the target plane, and k is the volume attenuation coefficient. Typical values of k in clear ocean, coastal and turbid waters are 0.05, 0.20 and 0.33 m^-1, respectively [9], and E_{\mathrm{o}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}} is the scalar irradiance upon the reflectance plane of the object. According to the theory of light intensity and luminous flux, regarding the underwater target as a Lambert emitter, the direct light (Lux) arrived at receiving plane can be calculated as

(3)E_{\mathrm{d}}(x,y,L)={\displaystyle \frac{n^{\mathrm{2}}{P}_{\mathrm{0}}{T}_{\mathrm{1}}{T}_{\mathrm{2}}\rho D^{\mathrm{2}}{\cos }^{\mathrm{3}}\phi {\cos }^{\mathrm{4}}{\phi }^{\prime }}{\mathrm{16}\mathrm{\pi }{l}^{\mathrm{2}}{f}^{\mathrm{2}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\alpha }{\mathrm{2}}}\right)\exp [kl(\sec \phi +\sec {\phi }^{\prime })]}},

where

\phi =\arccos {\displaystyle \frac{l}{\sqrt{{(x+d_{\mathrm{0}})}^{\mathrm{2}}+y^{\mathrm{2}}+l^{\mathrm{2}}}}}

{\phi }^{\prime }=\arccos {\displaystyle \frac{l}{\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+l^{\mathrm{2}}}}}

(x,y) denotes the coordinates on the image plane, l denotes the imaging distance, P_{\mathrm{0}} is the peak power of laser, {\tau }_{\mathrm{0}} is the pulse width, \alpha is the divergence angle (half angle) of laser beam after expansion, so the energy of one single pulse is E_{\mathrm{0}}=P_{\mathrm{0}}*{\tau }_{\mathrm{0}}; T_{\mathrm{1}} and T_{\mathrm{2}} are the optical transmittance of the optical transmission and receiving systems respectively, D/f is the relative aperture of receiving optical system, D denotes receiving diameter, f is the focal length; d_{\mathrm{0}} is the distance between optical transmission and receiving systems, \rho is the average reflectivity of the target, k and n denote attenuation coefficient and refractive index of water respectively.

As the forward-scattered light is scattered through the water after reflected by the object, according to the Fourier optics, the forward-scattered light component can be calculated by the convolution operation of direct light and the PSF of water. An empirical expression for the PSF of water has the form [10,11]:

(4)S(x,y)=({\mathrm{e}}^{-Gl}-{\mathrm{e}}^{-kl})\times F^{-\mathrm{1}}({\mathrm{e}}^{-Bl{f}_{\theta }})+{\mathrm{e}}^{-kl}\times \delta (x,y),

where k denotes the attenuation coefficient of water, G is a empirical constant related to k by , B is an empirical damping factor of scattering, l is the distance of light propagation in the water, f_{\theta } denotes the angular frequency and F^{-\mathrm{1}} stands for the inverse Fourier transform operator. This empirical PSF is the calculating basis of our model. As a result, the forward-scattered light arrived at sensors can be calculated by

(5)E_{\mathrm{f}\mathrm{s}}(x,y,l)=\left[{\displaystyle \frac{{\mathrm{e}}^{-Gl}-{\mathrm{e}}^{-kl}}{{\mathrm{e}}^{-kl}}}\times E_{\mathrm{d}}(x,y,l)\right]\ast g(x,y,l),

where

g(x,y,l)={\displaystyle \frac{F^{-\mathrm{1}}({\mathrm{e}}^{-Bl{f}_{\theta }})}{{\displaystyle \iint F^{-\mathrm{1}}({\mathrm{e}}^{-Bl{f}_{\theta }})\mathrm{d}x\mathrm{d}y}}}

and

f_{\theta }=l\times \sqrt{f_x^{\mathrm{2}}+f_y^{\mathrm{2}}}

denotes spatial frequency, * denotes convolution operation.

The back-scattered light is usually ignored in many studies for that it is reduced by separating the sensor from light source [12] and range-gating. However, it does exist even in the gated time and cannot be fully eliminated. In order to make the model more accurate, we still consider the effect of back-scattered light. For the calculation of back-scattered light, we can divide the water into several lamellas, so the distance of each lamellar from the image plane of sensors is

(6)\mathrm{\Delta }{R}_i=i{\displaystyle \frac{l_R}{n}},i=\mathrm{1},\mathrm{2},\cdots ,n,

where n is the number of lamellas, and l_R is the imaging distance based on receiving axis unlike the l which is based on transmitting axis. Like the calculation of the forward-scattered light, according to the Fourier optics, the back-scattered light component can also be derived by the convolution operation of the direct backlight and the PSF of water. Thus, the back-scattered light of the ith layer of water arriving at the image plane has the form:

(7)E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i)+\left\{\left[{\displaystyle \frac{{\mathrm{e}}^{-Gl}-{\mathrm{e}}^{-kl}}{{\mathrm{e}}^{-kl}}}\times E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i)\right]\ast g(x,y,l)\right\},

which is composed by the ith direct component E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i) and scattered component, where,

E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i)={\displaystyle \frac{\mathrm{\pi }{e}^{-k{R}_i}·T_{\mathrm{1}}·E_{si}(x,y,l,i)·\beta (\theta )·\mathrm{\Delta }{R}_i·{\cos }^{\mathrm{3}}\alpha }{\mathrm{4}{f}_{\theta }^{\mathrm{2}}}}·{\left[{\displaystyle \frac{l-\mathrm{\Delta }{R}_i·(i-\mathrm{0.5})}{l}}\right]}^{\mathrm{2}}

and

E_{\mathrm{s}\mathrm{i}}(x,y,l,i)={\displaystyle \frac{{\mathrm{e}}^{-GL}-{\mathrm{e}}^{-kL}}{{\mathrm{e}}^{-kL}}}\times \left\{\cos \alpha ·{\displaystyle \frac{P_{\mathrm{0}}}{\mathrm{4}\mathrm{\pi }\sin \left({\displaystyle \frac{\alpha }{\mathrm{2}}}\right)·R_i^{\mathrm{2}}}}\right\}\ast g(x,y,l)

denote the light intensity of the ith water lamellar, parameters G, k, T₁, \alpha, L are defined as same as that in the calculation of direct light. \beta (\theta ) denotes the volume scatter function, which has several kinds of forms developed by different researchers, like Duntley [13], Dolin et al. [14] and Wells [15]. In the model developed by this paper, the divergence angle (half angle) of laser beam and the relative aperture of receiving optics are below 10°, so it fits the range limitation of Wells’ theory () in which the form of \beta (\theta ) can be expressed as [8]

(8)\beta (\theta )={\displaystyle \frac{k\omega {\theta }_{\mathrm{0}}}{\mathrm{2}\mathrm{\pi }{({\theta }_{\mathrm{0}}^{\mathrm{2}}+{\theta }^{\mathrm{2}})}^{\mathrm{3}/\mathrm{2}}}},

where k and \omega are the total attenuation coefficient and scattering albedo respectively, and {\theta }_{\mathrm{0}} relates to the mean scattering angle; as a result, the back-scattered light arriving at image plane can be calculated as

(9)E_{\mathrm{b}\mathrm{s}}(x,y,l)={\displaystyle \sum _{i=\mathrm{1}}^n\left\{E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i)+\left[{\displaystyle \frac{{\mathrm{e}}^{-Gl}-{\mathrm{e}}^{-kl}}{{\mathrm{e}}^{-kl}}}\times E_{\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{d}}(x,y,l,i)\right]\ast g(x,y,l)\right\}},

which denotes the summation of the back-scattered intensity from all the water lamellas.

As for range-gated imaging, the back-scattered component is just from water around the target if the gating delay is set properly, we can set \mathrm{\Delta }{R}_i=l_R, then we have

(10)l+l_R={\displaystyle \frac{ct}{n}},

(11)l_R=l\cos {\theta }_s,

(12)\mathrm{d}l={\displaystyle \frac{c\mathrm{d}t}{n(\mathrm{1}+\cos {\theta }_s)}},

where {\theta }_s is the angle between receiving axis and transmitting axis. As a result, the back-scattered component can be calculated from the time integration by setting the time as the gating time.

Calculation of contrast transmittance and modulation transfer function (MTF)

In the three components of light received by sensor, the direct and forward-scattered components carry the information of target, while the back-scattered component is the background light along with varies noise which do not carry useful information, so, the contrast transmittance of image and MTF of water can be defined as mathematical operation of light components in the image at a location x, y:

(13)C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}={\displaystyle \frac{E_{\mathrm{d}}(x,y)+E_{\mathrm{f}\mathrm{s}}(x,y)}{E_{\mathrm{d}}(x,y)+E_{\mathrm{b}\mathrm{s}}(x,y)+E_{\mathrm{f}\mathrm{s}}(x,y)}},

(14)\mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{u}\mathrm{m}}={\displaystyle \frac{\mathrm{\pi }F(C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}})}{\mathrm{4}}},

where C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}} denotes the contrast transmittance of image, and \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{u}\mathrm{m}} denotes the MTF of medium (water). In order to simplify the calculation of C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}, some constant parameters determined by environments like the average reflectivity of the target, attenuation coefficient and refractive index of water, are simplified as K₁, K₂, etc. Then, the expression of E_{\mathrm{d}}(x,y,l) can be simplified as

(15)E_{\mathrm{d}}(x,y,l)=K_{\mathrm{1}}{\displaystyle \frac{P_{\mathrm{0}}{\left({\displaystyle \frac{D}{f}}\right)}^{\mathrm{2}}}{L^{\mathrm{2}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\alpha }{\mathrm{2}}}\right)\exp [K_{\mathrm{2}}l]}}{f}_{\mathrm{1}}(x,y),

and E_{\mathrm{b}\mathrm{s}} as

(16)E_{\mathrm{b}\mathrm{s}}(x,y,l)=K_{\mathrm{3}}{\displaystyle \frac{\mathrm{\pi }·\exp (K_{\mathrm{4}}l)·{\cos }^{\mathrm{4}}\alpha ·P_{\mathrm{0}}·(l-K_{\mathrm{6}})}{l^{\mathrm{2}}\sin \left({\displaystyle \frac{\alpha }{\mathrm{2}}}\right)\exp [K_{\mathrm{5}}l]}}·f_{\mathrm{2}}(x,y),

where f_{\mathrm{1}}(x,y) and f_{\mathrm{2}}(x,y) denote parameters as a function of position. Then, C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}} can be derived as

(17)C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}={\displaystyle \frac{K_{\mathrm{7}}·P_{\mathrm{0}}·{\left({\displaystyle \frac{D}{f}}\right)}^{\mathrm{2}}}{K_{\mathrm{7}}·P_{\mathrm{0}}·{\left({\displaystyle \frac{D}{f}}\right)}^{\mathrm{2}}+K_{\mathrm{8}}·\exp (K_{\mathrm{9}}·l)·{\cos }^{\mathrm{4}}\alpha ·\sin \left({\displaystyle \frac{\alpha }{\mathrm{2}}}\right)·{(l-K_{\mathrm{10}})}^{\mathrm{2}}}}·{\displaystyle \frac{f_{\mathrm{1}}(x,y)}{f_{\mathrm{2}}(x,y)}},

where parameters are defined as same as that in calculation of the light components.

We can see that, increasing the peak power of laser pulse (P_{\mathrm{0}}) and relative aperture of receiving optics (D/f), can enhance the value of C_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}, while the increase of distance between sensor and the target (l) has a opposite effect. As a result, based on the calculation of contrast transmittance, parameters of hardware can be adjusted for the purpose of image enhancement. For instance, increasing the laser power can improve the image quality. The caluculated contrast transmittance can also be used for underwater image evaluation. Higher contrast transmittance value means that more information is contained in the image which is of better quality.

Absorption and scattering are not the only factors hinder the underwater imaging, the turbulence in underwater can also severely limits underwater visibility. Hou et al. [7] gave the optical transfer function (OTF) of underwater optical turbulence using Kolmogorov turbulence model [16] by the form:

(18)$$\begin{gathered} \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}(f,r)=\exp \left[-\mathrm{1736}·K_{\mathrm{3}}·{\lambda }^{\frac{\mathrm{1}}{\mathrm{3}}}·f^{\frac{\mathrm{5}}{\mathrm{3}}}·r\right] \\ =\exp \left[-\mathrm{3.44}{\left({\displaystyle \frac{\lambda f}{R_{\mathrm{0}}}}\right)}^{\frac{\mathrm{5}}{\mathrm{3}}}·r\right], \end{gathered}$$

where f denotes the angular spatial frequency, r is the transmission range, \lambda =\mathrm{530}\mathrm{n}\mathrm{m} is the mean wavelength for underwater transmission, R_{\mathrm{0}} denotes the seeing parameter, K_{\mathrm{3}}=B_{\mathrm{1}}\chi {ϵ}^{-\frac{\mathrm{1}}{\mathrm{3}}} is the optical turbulence strength, B_{\mathrm{1}} is constant, \chi is the dissipation rate of temperature or sanity variances, and ϵ is the kinetic energy dissipation rate. In strong turbulent environments, typical values of former parameters are ϵ={\mathrm{10}}^{-\mathrm{4}},\chi ={\mathrm{10}}^{-\mathrm{10.5}},R_{\mathrm{0}}=\mathrm{0.003}. For the situations in non-turbulent environment, the value of MTF of underwater optical turbulence can be approximately \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}=\mathrm{1}.

The diffraction limit of the optical system of sensors is also an important factor inside MTF. When only one main lens exsits in the optical system, the diffraction limit factor can be defined as [17]

(19)

where f denotes the spatial frequency, and f_{\mathrm{c}\mathrm{o}} is the optical cut off frequency at the image plane, which has the form:

(20)f_{\mathrm{c}\mathrm{o}}={\displaystyle \frac{D}{\lambda ·f_l}},

where f_l denotes the focal length, and D represents the diameter of the lens, \lambda is the wavelength of operation. By far the CCD sensor is mainly chosen for the detection in range-gated imaging, the MTF of the CCD sensor depends on the size of its pixels, which can be defined as

(21)\mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{c}\mathrm{c}\mathrm{d}}={\displaystyle \frac{\sin (\mathrm{\pi }{d}_{\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}}\mathit{f})}{\mathrm{\pi }{d}_{\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}}f}},

where d_{\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}} denotes the size of pixels, and f represents the spatial frequency.

Therefore, the MTF of the whole system can be obtained by multiplying the MTF caused by each factors described above, which can be expressed as

(22)\mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}=\mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{u}\mathrm{m}}\times \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}\times \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}\times \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{c}\mathrm{c}\mathrm{d}}.

The curve of \mathrm{M}\mathrm{T}{\mathrm{F}}_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}} as a function of spatial frequency is shown in Fig. 2.

Fig.2 Comparison of relative MTFs of different factors: (a) MTF contribution from medium; (b) MTF contribution from turbulence; (c) MTF contribution from diffraction; (d) MTF contribution from ccd sensor; (e) MTF of the whole system

Full size|PPT slide

The optical devices of the imaging system such as lenses, apertures usually have circular symmetry, therefore, the PSF and MTF of the system can be obtained using hankel transform [18] from each other which has the following form:

(23)h(\theta ,L)=\mathrm{2}\mathrm{\pi }{\displaystyle \int J_{\mathrm{0}}(\mathrm{2}\mathrm{\pi }\theta \phi )}H(\phi ,R)\phi \mathrm{d}\phi ,

(24)H(\phi ,l)=\mathrm{2}\mathrm{\pi }{\displaystyle \int J_{\mathrm{0}}(\mathrm{2}\mathrm{\pi }\theta \phi )}h(\theta ,R)\theta d\theta ,

where H(\phi ,R) denotes the MTF of optical system, and h(\theta ,R) for the PSF, \phi is the spatial frequency, and l is imaging distance. The MTF and PSF which give the system response including the imaging system as well as the effects of medium are intuitive ways for describing the image formation. As a result, they can help to enhance the performance of image processing which will be discussed in Section 3.

Image processing

In order to improve the image quality of underwater imaging systems, the hardware can be upgraded such as increasing the output power of laser, improving the detecting ability of sensor and reducing error rate, etc. However, the upgrade of hardware will bring about the increase of substantial cost, and the balance is not easy to control, such as excessive increase of laser output power will cause serious impact of backscatter, this is proved by Ref. [19]. Therefore, improving image quality from the perspective of the image itself became necessary.

Image restoration

The relation between observed blurred image f^{\prime }(x,y) and original or uncorrupted signal f(x,y) can be described as

(25)f^{\prime }(x,y)=f(x,y)*h(x,y)+n(x,y),

where h(x,y) is the PSF of the system, * denotes convolution operation, n(x,y) denotes the noise of the system. Therefore, the original signal can be recovered by reversion or deconvolution from an accurate modeling of image system and medium [20]. The PSF of the imaging system calculated in Section 2 can be applied to various kinds of image restoration like wiener or blind restoration.

Image super-resolution reconstruction (SRR)

Like the contrast, resolution is another important factor for evaluating images. Image SRR [21] offers a possibility of improving image resolution beyond the hardware limitations. It has been widely studied and used recently. It refers to reconstructing a high-resolution (HR) image from one or multiple low-resolution (LR) images using the complementary information between image sequences. SRR can be divided into categories according to frequency domain and spatial domain, including interpolation, Papoulis-Gerchberg (PG) method, iterative back projection (IBP) method, projections onto convex sets (POCS) method, etc.

To apply the former calculated PSF to image super-resolution reconstruction, we choose the POCS method for its flexibility of incorporating prior knowledge of the imaging system. The main idea of the POCS method can be described as an iterative equation:

(26)f_{n+\mathrm{1}}=[P_{\mathrm{1}}{P}_{\mathrm{2}}\cdots P_k]f_n,

(27)f_{n+\mathrm{1}}=P\left[f_n+{\displaystyle \sum _{i=\mathrm{1}}^P\lambda P(g_i-H{f}_i)}\right],

where k denotes the number of limit sets, P is the projection operator, f^{n+\mathrm{1}} and f^n denote the SR image resulted from (n+1)th iteration and nth iteration, g represents the ith low resolution image, \lambda represents relaxed operator, and H denotes blurring operator which can be equivalent to the PSF of imaging system. Therefore, the combination of the POCS method and former calculated PSF can be exciting, research about which will be introduced in next section.

Experimental setup

Our model is applied to an underwater range-gated imaging system for image restoration and reconstruction, calculated contrast transmittance is used for image quality evaluation. Figure 3 shows the schematic diagram of the experimental system which consists of a Q-switch, frequency doubled Nd:YAG laser operated at 532-nm, an ICCD with programmable timing generator as external trigger controller, both the laser and the ICCD are put into a water tank, image data collected by the ICCD is transferred to computer and displayed by software.

Fig.3 Framework of range-gated imaging system

Full size|PPT slide

The parameters in our model should be set due to experimental facilities, a typical set of parameters are shown in Table 1.

Tab.1 Main physical properties of underwater range-gated imaging system

laser power (P0)	divergence angle (α)	optical transmittance (T1, T2)			relative aperture (D/f)
107 W	36 mrad	30%			0.25
image distance (L)	laser & CCD distance (d0)	average reflectivity (\rho)	attenuation coefficient (k)	refractive index (n)
40 m	20 cm	15%	0.25 m-1	1.35

The intensity distributions of the three light components are shown in Fig. 4.

Fig.4 Intensity of the light components: (a) non-scattered light; (b) forward-scattered light; (c) back-scattered light; (d) total light

Full size|PPT slide

The experiment was conducted in a boat pond which is \mathrm{150}\mathrm{m}\times \mathrm{6}\mathrm{m}\times \mathrm{4}\mathrm{m} with an attenuation of k=\mathrm{0.25}{\mathrm{m}}^{-\mathrm{1}} and scattering albedo of \omega =\mathrm{0.85}. The angle of field of view (FOV) is about 4°, which fits the range limitation of Wells’ theory (). An obtained image (original size 720 × 576 pixels, region of interest 256 × 256 pixels) of a 2 m × 2 m object in a distance of 35 m is shown in Fig. 5(a). Image restoration was performed with deconvolution filters and the PSF calculated in Section 2, the contrast values are used for evaluation. Restored results by different deconvolution filters are shown in Figs. 5(b)-5(d) with their contrast values in Table 2. It can be seen that these filters we used can contribute to improving image quality, but not significantly.

Fig.5 Restored images: (a) original image (size 256 × 256 pixels); (b) restored by Wiener filter; (c) restored by Lagrange filter; (d) restored by Lucky-Richardson filter

Full size|PPT slide

Tab.2 Contrast values of restored images

image	original	Wiener	Lagrange	Lucy-Richardson
contrast	24.9421	24.7385	24.5032	24.9456

Blind deconvolution filter which is the most popular and widely used filter in image restoration nowadays [22,23] can also be applied with our model. An important parameter for blind deconvolution is the number of iteration. A certain number of iteration can achieve better restoration, results neither do less or too more which will cause a waste of time. Restored results of different iteration number by PSF-based blind deconvolution are shown in Figs. 6(a)–6(c).

Fig.6 Restored images (size 256 × 256 pixels) with blind deconvolution by (a) 20 iteration; (b) 50 iteration; (c) 100 iteration

Full size|PPT slide

From the visual point of view by Fig. 6, we can see that restoration with iteration time of 50 performs better than that of 20, so more iteration can enhance the performance of blind deconvolution, however an iteration time of 100 makes the restoration result more ambiguous. For underwater imaging or detecting, image segmentation and information extraction are even more important, and the effect of which depend very much on the image quality which is based on the computer vision; as a result, using objective evaluation criteria such as the contrast value to gauge the image quality is more effective than the subjective evaluation. The contrast values of the restored images are shown in Table 3.

Tab.3 Contrast values of restored images

image	20 iteration	50 iteration	100 iteration
contrast	25.0099	29.5063	27.6188

It can be seen that the blind deconvolution performs better than former filters, but ringing artifacts exist in the restored result. This could be reduced using regularization such as edge detection. From the contrast value list, we can see that more iteration does not reflect a better effect, and it can be explained by that blind deconvolution algorithm has worse noise tolerance under large iteration number.

So, the quality of image can be enhanced by image restoration, and with which we can deduce the best filter and best iteration number for image restoration.

10 frames are extracted from the test video sequences collected by ICCD for super-resolution reconstruction including the sample image we used for image restoration. Figure 7 shows the reconstruction results of various SRR methods.

Fig.7 Reconstructed images (size 512 × 512 pixels) by (a) bilinear interpolation; (b) cubic convolution interpolation; (c) PG method; (d) Iterative back projection method; (e) POCS method; (f) PSF-based POCS method

Full size|PPT slide

From the visual point of view by Fig. 7, the differences of reconstructed results are not obvious, so we can only use the objective evaluation. The contrast values of reconstructed images are shown in Table 4.

Tab.4 Contrast values of reconstructed images

image	bilinear	cubic	PG	IBP	POCS	PSF-POCS
contrast	14.7464	19.1722	29.2447	29.2302	39.2438	51.8804

From the contrast values of the reconstructed images, we can see that the results of interpolation algorithms are not desirable, while other methods can offer a relative better result, this is due to the fact that creating non-exists pixels blindly blurs the boundaries of black and white, which degrades the images of black and white stripe resolution board. Reconstructed images have ringing artifacts owing to the steep cut-off frequency which also can be reduced by regularization such as edge detection. Figure 7(f) is the result of the POCS method based on our calculated PSF. As can be obviously seen from Table 4, this method performs better than other traditional SRR methods. As a result, we can conclude that the PSF based POCS method can substantially enhance the performance of super-resolution reconstruction, which can achieve a best result currently. More future work should be conducted such as the introduction and comparison of traditional PSF models for the POCS method or other SRR methods.

Conclusions

In this paper, an underwater imaging model based on the formation of underwater images is presented, along with the retrieval of optical properties. The model includes the beam propagation, responses of medium as well as sensor, and the effects of underwater turbulence. Issues implemented on underwater imaging such as denoising, image enhancement and restoration are addressed and discussed. The model was applied to a range-gated underwater imaging system for image restoration and super-resolution reconstruction in which varies filters and methods were chosen for comparison, the results show that the calculated MTF and PSF can be used to enhance the performance of both the restoration and reconstruction. Further works can be carried out including more complex methods for image restoration or super-resolution reconstruction applications. The model can also be applied to other underwater imaging systems which have a similar principle of imaging.

Acknowledgements

This paper was supported by the National Natural Science Foundation of China (Grant No. 61008050).

References

1 Duntley S Q. Light in the sea. Journal of the Optical Society of America, 1963, 53(2): 214-233

2 Acharekar M A. Underwater laser imaging system (ULIS). In: Proceedings of SPIE. 1997, 3079: 750

3 Chang P C Y, Walker J G, Hopcraft K I, Ablitt B, Jakeman E. Polarization discrimination for active imaging in scattering media. Optics Communications, 1999, 159(1-3): 1-6

4 McLean E A, Burris H R Jr, Strand M P. Short-pulse range-gated optical imaging in turbid water. Applied Optics, 1995, 34(21): 4343-4351

PMID

5 Masters B R, Barrett H H, Myers K J. Foundations of Image Science. In: Saleh B E A, ed. Wiley Series in Pure and Applied Optics. New York: Wiley-Interscience, 2007, 31(2): 114-115

6 Mertens L E, Replogle F S Jr. Use of point spread and beam spread functions for analysis of imaging systems in water. Journal of the Optical Society of America, 1977, 67(8): 1105-1117

7 Hou W L. A simple underwater imaging model. Optics Letters, 2009, 34(17): 2688-2690

PMID

8 Hou W L, Gray D J, Weidemann A D, Arnone R A. Comparison and validation of point spread models for imaging in natural waters. Optics Express, 2008, 16(13): 9958-9965

PMID

9 Jaffe J S. Computer modeling and the design of optimal underwater imaging systems. IEEE Journal of Oceanic Engineering, 1990, 15(2): 101-111

10 Rosenfeld A, Kak A C. Digital Picture Processing. 2nd ed. New York: Academic, 1982

11 Wells W H. Loss of resolution in water as a result of multiple small-angle scattering. Journal of the Optical Society of America, 1969, 59(6): 686-691

12 Yu Y F, Liu F C. System of remote-operated-vehicle-based underwater blurred image restoration. Optical Engineering, 2007, 46(11): 116002

13 Duntley S Q. Underwater Lighting by Submerged Lasers and Incandescent Sources. San Diego: Scripts Instituition of Oceanography, University of California, 1971

14 Dolin L S, Gilbert G D, Levin I, Luchinin A. Theory of Imaging Through Wavy Sea Surface. Nizhniy Novgorod: Russian Academy of Sciences, Institution of Applied Physics, 2006

15 Wells W H. Theory of small angle scattering. AGARD Lecture Series, 1973, 61: 3.3.1-3.3.19

16 Tatarskii V I. Wave Propagation in a Turbulent Medium. Silverman R S, translated. New York: McGraw-Hill, 1961: 285

17 Bonnier D, Stephane C, Yvves L, Bruno L, Marc L, Pierre G. Modelling of active TV system for surveillance operations. In: Proceedings of SPIE–The International Society for Optical Engineering. 1999, 3698: 217-228

18 Hou W L, Gray Deric J, Weidemann Alan D, Fournier Georges R, Forand J L. Automated underwater image restoration and retrieval of related optical properties. IEEE International Geoscience and Remote Sensing Symposium, 2007: 1889-1892

19 Han H W, Zhang X H, Ge W L. Performance evaluation of underwater range-gated viewing based on image quality metric. In: Proceedings of the International Conference on ICEME’09. 2009, 4: 441

20 Gonzalez R C, Woods R E. Digital Image Processing. 2nd ed. Upper Saddle River, NJ: Prentice Hall, 2002

21 Park S C, Park M K, Kang M G. Super-resolution image reconstruction: a technical overview. IEEE Signal Processing Magazine, 2003, 20(3): 21-36

22 Zhang J L, Zhang Q H, He G M. Blind deconvolution of a noisy degraded image. Applied Optics, 2009, 48(12): 2350-2355

PMID

23 Wang Z F, Tang Y D. Semi-blind image restoration based on Chan-Vese denoising model. Chinese Optics Letters, 2008, 6(6): 405-407