BLIND DECONVOLUTION ON UNDERWATER IMAGES FOR GAS BUBBLE MEASUREMENT

Marine gas seeps, such as in the Panarea area near Sicily (McGinnis et al., 2011), emit large amounts of methane and carbon-dioxide, greenhouse gases. Better understanding their impact on the climate and the marine environment requires precise measurements of the gas flux. Camera based bubble measurement systems suffer from defocus blur caused by a combination of small depth of field, insufficient lighting and from motion blur due to rapid bubble movement. These adverse conditions are typical for open sea underwater bubble images. As a consequence so called ’bubble boxes’ have been built, which use elaborate setups, specialized cameras and high power illumination. A typical value of light power used is 1000W (Leifer et al., 2003). In this paper we propose the compensation of defocus and motion blur in underwater images by using blind deconvolution techniques. The quality of the images can be greatly improved, which will relax requirements on bubble boxes, reduce their energy consumption and widen their usability.


INTRODUCTION
In areas like the Panarea area near Sicily (McGinnis et al., 2011) or the Tommelitten (Schneider von Deimling et al., 2011) in the North Sea, large amounts of carbon dioxide and methane are emitted from underwater gas seeps.
Because of the significance of underwater gas emissions on marine life (Ishimatsu et al., 2004), (Dando and Hovland, 1992) and because of the impact on the climate, accurate measurements of the properties of the occurring gas bubbles are necessary, to be able to accurately estimate their impact (Liang et al., 2011).Different ways to measure gas bubbles have been pursued in recent years.Most prevailing are acoustic (McGinnis et al., 2006) and optical measurement.We focus on camera based measurements, because they allow more accurate observations of gas bubbles.High quality images allow flow measurement with an error of less than 6 percent, as confirmed by practical tests (Zelenka, 2014).
Underwater imaging conditions are far from ideal for this kind of observations.They typically suffer from low light, due to power restrictions.The flamboyant bubbles can reach velocities of up to 25cm per second (Leifer et al., 2003).To prevent motion blur, high shutter speed and low integration times are necessary, which results in even less light on the image sensor (for detailed experiments see (Leifer et al., 2003)).Using a fast lens, which means a high aperture, is not a solution, because this lowers the depth of field and leads to potentially defocused images.To overcome this problems, so called bubble boxes have been developed (Thomanek et al., 2010).A bubble box as shown in Figure 1 includes a camera, channeling elements to concentrate the bubbles in the cameras depth of focus, and a backlight illumination.Although front side illumination is known (Zielinski et al., 2010), for quantitative bubble measurement backlight illumination is prevalent due to higher light efficiency and bubble rim visibility.
This does not solve the problem entirely, as strong lighting is still an issue, as the 1000 Watt illuminations in (Wang and Dong,

Camera Backlight Illumination
Figure 1: Drawing of a typical bubblebox 2009) and (Leifer et al., 2003) demonstrate.Because of the imaging conditions (see above), the bubble box also requires a high speed camera, such as in (Cheng and Burkhardt, 2006) (Leifer et al., 2003).Both requirements are costly and unwieldy in expedition use.
Another factor often neglected in the literature is the presence of light scattering in the underwater ocean settings, caused by impurities in the water.This effect disturbs the light propagation and leads to blurry imaging and lowered accuracy, similar to defocus blur.
In this paper a different approach is developed, which is the application of suitable blind deconvolution techniques.The goal is not to avoid any blur or defocus, e.g. using high speed cameras and high power illuminations, but to compensate for these effects.This allows a simplified buildup of a bubble box or images of higher quality in existing builds, while allowing more flexible experiments.
Blind deconvolution techniques in general underwater settings have been used with success in (Fan et al., 2010a), (Fan et al., 2010b) and with the adapted Richardson-Lucy algorithm in (Wu et al., 2013), while the application on bubbly flow images is novel.
In this paper we focus on two different algorithms, the Richardson-Lucy standard algorithm and a more recent technique.
The Richardson-Lucy Algorithm by (Richardson, 1972) and (Lucy, 1974) with adaption to blind deconvolution by (Fish et al., 1995) is a fast standard deconvolution algorithm with a matlab implementation (TheMathworks, 2015), which includes additional dampening.In recent years significant progress in the field of blind deconvolution has been made (Levin et al., 2011).Of particular interest are the gradient sparsity methods, because the gradient sparsity prior, which was originally established for natural images, also applies very well on underwater gas bubble images.
Which ideally consist of a uniform background with sparsely distributed bubbles with high gradient rims.
We compare the heavy-tailed gradient sparsity MAP(Maximum a posteriori) blind deconvolution method by (Kotera et al., 2013) with the preceding standard technique and test whether they are suitable and applicable to underwater gas bubble images.

BLIND DECONVOLUTION TECHNIQUES
We model the observed image O with: where S is the undisturbed signal disturbed by B, the blur kernel or point spread function (PSF) and n models added noise during the image formation.The aim of deconvolution algorithms is the reconstruction of the undisturbed signal S from the observation O.If the PSF is unknown, a blind deconvolution algorithm with the added difficulty of estimating an unknown B is necessary.
We use blind deconvolution techniques, because the exact blur model of an underwater image is often unknown.Moreover motion blur is dependent on bubble speed, the very parameter to be measured.Note that this reconstruction problem of two unknowns from a single observation is inherently ill posed.A more in depth analysis of this problem can be found in (Levin et al., 2011) and (Perrone and Favaro, 2014).
The Richardson-Lucy deconvolution algorithm views this image formation model from a Bayesian perspective.For discrete event A, B with an index k and j the Bayes theorem states: From this it follows: An iterative algorithm with iteration index i is defined by: Applied to the deconvolution problem of Equation 1, this means associating the observed image O with P (B), the undisturbed signal S with P (A) and the blur kernel B with P (B|A), similar to (Fish et al., 1995).In convolution notation this can be written as: If the blur kernel is to be estimated, the analog update formula is: The Richardson-Lucy blind deconvolution algorithm (Fish et al., 1995) approaches the reconstruction problem of image S and PSF B by alternating iterations in which the image and the PSF are refined.In this paper we use the standard matlab inplementation, which includes an additional dampening, to decrease noise propagation in the course of the iterations (TheMathworks, 2015).
In the following the blind deconvolution algorithm by (Kotera et al., 2013) is presented.It relies on a gradient sparsity prior that is maximized in the MAP(Maximum a posteriori)-estimation framework.
Given the Bayesian theorem and the image model Equation 1, neglecting noise, the goal is finding the MAP(Maximum A Posteriori) of , where P (S), etc. denote the probabilities and conditional probabilities of O,S, and B, while assuming a Gaussian distribution with parameter µ of Q(S) and R(B) are regularizers of S and B respectively, to guide the optimization towards likely solutions.For the image S the regularizer Q regards the spatial gradients Dx and Dy of S: with p ∈ (0, 1) and index i over the entirety of the partial derivatives Dx and Dy.The variable p defines the norm enforced on the gradients.
The blur kernel B is regularized with where with the running index i over all Bi in B to enforce sparsity and non-negativity.
Instead of solving the MAP-problem directly, the negative logarithm of the target function is minimized, a common technique for MAP-estimation (Meintrup and Schaeffler, 2005).Constants are disregarded, because they do not influence the minimum: = − log(P (O|S, B)P (S)P (B)) In (Kotera et al., 2013) this optimization problem is solved using an augmented Lagrangian method in a multi-scale approach by alternatingly optimizing image and blur kernel.

APPLICATION AND RESULTS
In a laboratory setting a Point Grey Grasshopper 2 camera , backlight illumination and an air seep in a water tank are used in a similar build as shown in Figure 1 to obtain images of bubbles with varying acquisition parameters.The light power and camera setting such as focus, aperture and shutter speed are adjustable and provide a good foundation for experimentation with different kinds and strengths of blur.Images acquired with this system are used in the following.

Defocus blur
The two main sources of blur identified in a bubble box are defocus blur and motion blur.
First we will concentrate on defocus blur, which occurs if the bubbles are outside of the depth of field, centered at the focal plane.The depth of field determines the utility of the bubble meter, because it limits the distance range in which bubbles can be measured.A large depth of field requires a lens with a small aperture.However due to the difficult lighting conditions in underwater settings discussed in the introduction, the aperture should be as large as reasonably possible.Solving this conflict by selecting an appropriate lens is a design decision requiring a compromise.
In the following paragraph we show how the blind deconvolution techniques presented in the previous section apply to this setting and how they mitigate blur.A comparison of the results from Richardson-Lucy and gradient sparsity MAP-estimation (Kotera et al., 2013) on an image with mild defocus blur is shown in Figure 2 shows sharper edges and good deblurring.While the Richardson-Lucy algorithm is much faster, it suffers from strong ringing artifacts and increased noise.
The restoration using MAP-estimation leads to much better results, which less artifacts, less increase in noise and sharper bubble rims.Therefore it is our method of choice for the following experiments.
Using deconvolution on a noisy image as in Figure 3

Motion blur
Motion blur occurs if due to a combination of bubble motion and shutter speed, the light from a bubble is spread over several pixel.Choosing a very higher shutter speed, with a very short integration time can solve this problem, but the light intensity available for the image sensor is also determined by the integration time.Therefore this leads to dark and noisy images, unless a high power illumination is chosen.
In Figure 4 a series of images acquired with different shutter speeds is shown.In this image series automatic gain control together with strong illumination and a large aperture lens were used to allow capturing images with constant brightness, even though the integration time differs significantly.Note that such a large aperture lens also negatively influences the depth of field.
The left image was acquired with a very short integration time of only 0.0004s, set by the camera control software, shows no motion blur, while the image in the middle shows some motion blur with an estimated blur radius of approx.10 pixel and the image on the right shows very strong motion blur with a blur radius of approx.15 pixel.While applying the deconvolution algorithm on images with different intensities of motion blur, it becomes clear that there are limits to how strong the blur can be for the deconvolution to achieve a full compensation of the blur effects.We define a deconvolution successful only if the result is visually plausible and the blurred edges are fully restored without introducing artifacts, that may hinder an evaluation of the images.
Figure 5 shows the deconvolution of image with medium motion blur, the same image as in Figure 4(b).The rim of the deconvolution result is sharp, and resembles the low integration time image in Figure 4(a).
Next we measure the bubbles before and after deconvolution.For the bubble shown in Figure 6(a) which is also part of Figure 5  The initialization is shown in green, the snake is shown in red, while the fitted ellipse in blue.
deconvolution from 42 pixel to 37 pixel.This change is expected as blur around the rim of the bubble is removed, although the result is not as good as in artificial images (see Figure 9) and the restored contour could be even smaller.
Furthermore we measure the bubble size with automated algorithms (Zelenka, 2014).The different steps of the bubble segmentation are shown in Figures 6(c,d).A snake, an active contour, is initialized with the green bounding rectangle and iteratively contracts towards higher image gradients under smoothness constraints.The final state of the snake is shown in red.To accurately measure the size the blue ellipse is fitted into the snake.The automatically measured bubble size remains constant with approx.38 pixel height, which is confirmed by the segmentations in Figure 6.The weak blur on the outer rim of the bubble is disregarded and does not disturb the image gradient based algorithm, which is robust enough to cope with this blur.
For stronger motion blur shown in Figure 7 and Figure 8 the motion blur is too strong to be removed with a same parameters used in Figure 5. On the original image, the snake method fails to initialize and attempts to use multiple snakes for a single bubble (see Figure 8) and due to the blur even a manual size measurement is hard.The deconvolution with standard parameters does not fully restore the original edges but creates new edges further around the bubble.
Changing the deconvolution parameters, allowing a larger PSF and weighting the blur kernel regularization in the blur kernel estimation the 2.5 times more than the data term, compensates the stronger blur (see Figure 7).It restores the sharp rim of the bubbles, similar in sharpness to non-blur image in Figure 4(a) and allows an accurate automatic measurement as in Figure 8.Because of this advantage, we think that the strong deconvolution artifacts, apparent by the white halo around the bubbles, can be tolerated.demonstrates that a successful deconvolution with gradient sparsity MAP can provide an accurate restoration of the non-blurred image and can therefore be used for further evaluation.

CONCLUSIONS
In our experiments we have observed, how motion and defocus blur can influence the image and how deconvolution can restore them.Blur can lead to an overestimation of the bubble size or prevent automatic bubble measurements.A high measurement accuracy is necessary, because the bubble radius influences the volume with the third power and bubble shape parameters like ellipticity are important to gas dissolution models like (McGinnis et al., 2006) and (Liang et al., 2011).
The two different algorithms showed a strong difference in performance.The standard blind Richardson-Lucy deconvolution, provided a good initialization, shows acceptable results for small defocus blur, even though ringing artifacts occur.The gradient sparsity MAP deconvolution can compensate much stronger blur and even some motion blur without artifacts.For strong motion blur, strong restoration artifacts have to be accepted.
Nevertheless with deconvolution, on images which could not be evaluated beforehand, the blur is removed, sharp contours are restored and the automatic evaluation is enabled.This allows clear measurements and lowered overestimation, which are strong advantages.The complete potential of applying blind deconvolution on bubble images is shown in the results seen in the compensation of artificial motion blur.
In case of non-changing conditions, such as bubble speed and water parameters, the restored blur kernel varies little.If the imaging conditions allow it, we propose using the blind-deconvolution every few images to reconstruct the blur kernel from an initial image and reuse it for in the following images.This means for the consecutive images only non-blind deconvolution is required and the computational costs can be significantly reduced.The gradient sparsity MAP deconvolution takes 9s on the image in Figure 2, while the non-blind deconvolution Richardson-Lucy algorithm with the blur kernel set by the blind deconvolution is much faster with 1s.
Currently only spatially invariant blur is considered.For future work, a spatially variant blur kernel estimation could be beneficial, because although the motion blur of the bubble in an image is similar in direction and magnitude, due to the variations in bubble motion it is not identical.Furthermore, with a small depth of field not all bubbles are in focus and a bubble may even change between focus and defocus over time.As future work, this means that the blur model could be estimated for every bubble individually and incorporated or partially derived from the measured bubble motion.
In summary, we show in this paper, that blind deconvolution gives significant improvements in image quality and extends the usable depth of field in an underwater environment.We see that a certain amount of motion and defocus blur can be permitted with deconvolution, while maintaining high measurement accuracy .The compensation of image blur caused by defocus and motion permits more accurate flow measurements and relaxed requirements on underwater observatory equipment.This allows the design of a bubble box with higher measurement accuracy in a simplified setup, larger effective depth of field, and due to a higher tolerance to motion blur, less powerful illumination and camera.

Figure 2 :
Figure 2: Deconvolution of defocus blur.Both Richardson-Lucy and gradient sparsity MAP blind deconvolution achieve usable result, however the Richardson-Lucy method requires a good initialization and shows inferior bubble edge sharpness.
(a) reveals a negative property of the deconvolution.It increases the already present noise as shown in Figure 3(b).

Figure 3 :
Figure 3: Noisy blurred image and its gradient sparsity deconvolution result.

Figure 4 :
Figure 4: Bubble images acquired with different integration times with automatic gain control.
Figure 5: Medium motion blur image and gradient sparsity deconvolution.The dotted blue rectangle shows the cropped segment of Figure 6.

Figure 6 :
Figure 6: Measuring the bubble size with the snake algorithm.The initialization is shown in green, the snake is shown in red, while the fitted ellipse in blue.
Figure 7: Strong motion blur image and gradient sparsity deconvolution with changed weighting on the regularization allowing a more flexible blur kernel.The dotted blue rectangle shows the cropped segment of Figure 8.

Figure 8 :Figure 9 :
Figure 8: Snake based measurement of a bubble from Figure 7, colors as in Figure 6.To test whether the result of blind deconvolution is accurate, which means it does not over-or under-compensate the blur effects, tests on blurred data with known ground truth is necessary.In the following experiment a sharp image is artificially blurred, then the deconvolution algorithm is applied and the result is compared with the sharp image.The resulting images are shown in Figure 9 with the original image on the left, with the blurred image in the middle and the deconvolution on the right side.This