3D Surface Reconstruction of Underwater Objects

In this paper, we propose a novel technique to reconstruct 3D surface of an underwater object using stereo images. Reconstructing the 3D surface of an underwater object is really a challenging task due to degraded quality of underwater images. There are various reason of quality degradation of underwater images i.e., non-uniform illumination of light on the surface of objects, scattering and absorption effects. Floating particles present in underwater produces Gaussian noise on the captured underwater images which degrades the quality of images. The degraded underwater images are preprocessed by applying homomorphic, wavelet denoising and anisotropic filtering sequentially. The uncalibrated rectification technique is applied to preprocessed images to rectify the left and right images. The rectified left and right image lies on a common plane. To find the correspondence points in a left and right images, we have applied dense stereo matching technique i.e., graph cut method. Finally, we estimate the depth of images using triangulation technique. The experimental result shows that the proposed method reconstruct 3D surface of underwater objects accurately using captured underwater stereo images.


INTRODUCTION
Optical imaging of the ocean floor offers scientists high levels of detail and ease of interpretation. However, light in underwater suffers from significant attenuation and backscatter, limiting the practical coverage of a single image to only a few square meters. For many scientific surveys, however, the area of interest is large, and can only be covered by hundreds or thousands of images acquired from a robotic vehicle or towed sled. Such surveys are required to study hydrothermal vents and spreading ridges in geology, ancient shipwrecks and settlements in archeology, forensic studies of modern shipwrecks and airplane accidents, and surveys of benthic ecosystems and species in biology.
In the field of 3D reconstruction, terrestrial applications have encouraged extensive work over the last three decades; on the other hand a limited amount of underwater applications have been explored primarily for mapping and positioning. Shape from X is a generic name for techniques that extract shape from images. Normally, in underwater environment the optical sensing techniques include Shape from Stereopsis [24], Shape from Photometric Stereo [15], Shape from Motion [7] and Active Stereo [14]. In the field of non-optic underwater sensing, acoustic cameras are employed for 3D mosaic reconstruction [2], whereas both acoustic and optic cameras are often used providing scene information that cannot be recovered from each sensor alone.
Three dimensional scene structures captured by a camera may be detected and acquired observing the apparent motion of brightness patterns from images. The primary visual motion cue useful for shape acquisition is the perceived movement of brightness patterns, known as optical flow [6] which is an approximation of the 3D world motion field. The 3D reconstruction from differential motion cues requires accurate optical flow computation. In the last fifteen years, theoretical developments in visual motion studies have established a unified framework for the treatment of the Structure from Motion (SFM) and Structure from Stereo (SFS) problems [3], also known as 3D Reconstruction from Multiple Views. 3D reconstruction from multiple views involves extracting target features from one image, matching and tracking these features across two or more images, and using triangulation to determine the position of the 3D target points relative to the camera. Visual motion methods have been well-studied, requiring densely-sampled image sequences. Instead, the trade-off in stereo vision is between the stereo correspondence problem and a more accurate and robust 3D reconstruction [19]. A large amount of works has addressed the correspondence problem, attempting to overcome the various difficulties of the large-displacement correspondence problem: occlusions, large rotations and disparities, photometric and projective distortions [12,23].
The problem of outliers has been solved by the deployment of robust estimation methods [25]. The motion and structure may be estimate by using several recursive schemes. Extended Kalman Filter is the most popular approach for jointly estimation of motion and structure with satisfactory results. However, the computation cost of this approach grows cubically with the amount of features, causing a great bottleneck for real-time performance. Stereo matching is one of the most active research areas in computer vision. Stereo matching is a hard problem due to ambiguity in un-textured and occluded areas. Only dominant features, such as points of interest, can be matched reliably. This motivates the development of progressive approaches [19]. The reduced local disparity in search range makes progressive approaches very efficient in computation and robust. However, the seed initialization remains a computational bottleneck, although, robustness can be improved by enforcing the left-right symmetry constraint.
Recently, Graph Cuts and Belief Propagation furnish combinatorial optimization frameworks in which the performance for global stereo algorithms are considerably improved according to an evaluation framework applied to a standard reference stereo data set [19]. However, the underlying brightness constancy assumption of combinatorial optimization methods severely limits the range of their applications. The earliest attempts for 3D reconstruction use methods based on volume intersection, as Shape from Silhouette [10]. Traditional methods, such as stereo, handle large visibility changes between images by solving the correspondence problem between images. The most prominent approaches in 3-D reconstruction are the Voxel Coloring [22] and the Space Carving [9]. These approaches use the color consistency to distinguish surface points from 32 other points in a scene. Cameras with an unoccluded view of a non-surface point see surfaces beyond the point, and hence inconsistent colors, in the direction of the point. Initially, the environment is represented as a discretized set of voxels, and then the algorithm is applied to color the voxels that are part of a surface in the scene. Another promising approach in 3D reconstruction is the Marching Cubes algorithm for rendering iso-surfaces from volumetric scan data [11]. The algorithm produces a triangle mesh surface representation by connecting the patches from all cubes on the iso-surface boundary.
In underwater scenario, research is directed at exploring the use of vision, potentially in conjunction of other sensors, to automatically control unmanned submersibles, including positioning and navigation by utilizing a photomosaic as a two-dimensional visual map. Recent activities combine video imagery taken from multiple views of a scene to derive size and depth measurements and 3D reconstructions. These activities support (semi-) autonomous or operator-supervised missions pertaining to automatic vision guided station keeping, location finding and navigation, survey and mapping, trajectory following and online reconstruction of a composite image, search and inspection of subsea structures. These tasks require an accurate estimation of camera position, together with fast, accurate correspondence determination, particularly for real-time registration. Common sources of error include non-planar seafloor, moving objects, illumination variations, transect superposition, positioning drift. A number of studies over the last several years have also addressed the 3D reconstruction for various applications. Khamene and Negahdaripour incorporate cues from stereo, motion and shading flow for 3-D reconstruction in underwater [8]; Majidi and Negahdaripour suggest the use of 3D reconstruction for global alignment of 3D sensor positions [13]; Nicosevici et al introduce 3D reconstruction from motion video and representation of the surface topography by piecewise planar surfaces for the construction of ortho-mosaics [16]. Hogue et al have developed a stereo vision-inertial sensing device deployed to reconstruct complex 3D structures in both the aquatic and terrestrial domains [5]. The sensor temporally combines 3D information, obtained using stereo vision algorithms with a 3 DOF inertial sensor. The resulting point cloud model is then converted to a volumetric representation and a textured polygonal mesh is extracted using the Marching Cubes algorithm [11].
In this paper, we have used stereo pair of cameras to capture the underwater images. We rectify captured underwater images using uncalibrated rectification technique. The dense stereo matching technique is applied to find correspondence points in both left and right images. We estimate the depth of an underwater image by using triangulation method. The dense depth map is converted to triangular surface meshes using the Delaunay triangulation algorithm. Thus, a texture mapping can be applied with relative easy and efficiency to the object.
The proposed method helps in the 3D surface reconstruction of underwater objects, which will be further helpful in recognition of underwater objects like mines, shipwrecks, pipelines and telecommunication cables, etc.
The organization of the paper is as follows: Section 2 presents the detailed description of steps involved in the proposed method like pre-processing of underwater images, uncalibrated stereo image rectification, stereo matching (correspondence) and depth estimation. Section 3 presents the experimental results. The Section 4 draws a conclusion.

2.3D RECONSTRUCTION
The proposed 3D surface reconstruction method for underwater objects involves several sequential steps. These steps are explained with detailed description of the algorithm.

Preprocessing
A major difficulty to process underwater images comes from light attenuation. Light attenuation limits the visibility distance, at about twenty meters in clear water and five meters or less in turbid water. The light attenuation process is caused by the absorption (which removes light energy) and scattering (which changes the direction of light path). Dealing with this difficulty, underwater imaging faces to many problems: first the rapid attenuation of light requires attaching a light source to the underwater vehicle providing the necessary lighting. Unfortunately, artificial lights tend to illuminate the scene in a non uniform fashion producing a bright spot in the center of the image and poorly illuminated area surrounding. Then, the floating particles highly variable in kind and concentration, increase absorption and scattering effects: they blur image features (forward scattering), modify colors and produce bright artifacts known as "marine snow".
We have adopted the preprocessing technique proposed by Stephene Bazielle et al. (2006) [1] to preprocess captured underwater images. This method involves sequence of preprocessing steps: (1) Removing potential moir'e effect (2) Resizing and extending symmetrically the image to get a squared image whose size is a power of two (1) Removing potential moir'e effect. A moir'e effect has the appearance of a wavy repetitive pattern on the image. It is not an underwater perturbation, and it is often considered as aliasing phenomena. Sampling moir'e mainly occurs in the analog to digital conversion process. Moir'e pattern is removed via spectral analysis by detecting peaks in the Fourier transform and deleting them assuming that they represent the moir'e effect [21]. Only few images suffer from moir'e degradation but removing it is important because the following processes enhance contrast so enhance the moir'e effect and consequently highly degrade results.
(2) Resizing and extending symmetrically the image to get a squared image whose size is a power of two. Symmetric extension prevents from potential border effects and resizing to squared image speeds up the following process by enabling to use fast Fourier transform and fast wavelet transform algorithms.
(3) Converting color space from RGB to YCbCr (Luminance Chrominance). This color space conversion allows us to work only on one channel instead of processing the three RGB channels. In YCbCr color space we process only the luminance channel (Y) corresponding to intensity (4) Homomorphic Filtering. The homomorphic filtering is used to correct non uniform illumination and to enhance contrasts in the image. It's a frequency filtering, preferred to others techniques because it corrects non uniform lighting and sharpens the edges at the same time. We consider that image is a function of the product of the illumination and the reflectance as shown below.
where ( , ) f x y is the image sensed by the camera, ( , ) i x y the illumination multiplicative factor, and ( , ) r x y the reflectance function. If we take into account this model, we can assume that the illumination factor changes slowly through the view field, therefore it represents low frequencies in the Fourier transform of the image. On the contrary reflectance is associated with high frequency components. By multiplying these components by a high-pass filter we can then suppress the low frequencies i.e. the non uniform illumination in the image. The algorithm can be decomposed as follows: -Separation of the illumination and reflectance components by taking the logarithm of the image. The logarithm converts the multiplicative effect into an additive one: ( , ) ln( ( , )) ln( ( , ) ( , )) ln( ( , )) ln( ( , )).
-Computation of the Fourier transform of the log-image ( , ) ( , ) ( , ). x -The High-pass filter is applied to the Fourier transform which decreases the contribution of low frequencies (illumination) and also amplifies the contribution of mid and high frequencies (reflectance), sharpening the edges of the object in the image: ( , ) ( , ) ( , ) ( , ) ( , ), -For each subband of each level except for the lowpass residual: Computation of the signal variance using (Eq. 7) and modification of the noisy wavelet coefficient according to the (Eq. 8).
where i y is the child and 1 i y  its parent.
-Inversion of the multiscale decomposition to reconstruct the filtered image.
(6) Anisotropic Filtering. This filter smooths the image in homogeneous area, preserves edges and enhance them. It is used to smooth textures and reduce artifacts by deleting small edges amplified by homomorphic filtering. It also removes or attenuates unwanted artifacts and remaining noise. We follow the algorithm proposed by Perona and Malik [17]. This algorithm is automatic so it uses constant parameters selected manually. The previous step of denoising is very important to obtain good results with anisotropic filtering. Anisotropic filtering is usually used as long as result is not satisfactory. In our case we loop only few times set to constant value, to preserve a short computation time. One loop of the algorithm can be decomposed as follows: For each pixel -Computation of the nearest-neighbor differences and computation of the diffusion coefficient in the four directions North, South, East, West. Many possibilities exist for this calculation; the easiest way is as follows:   (7) Adjusting image intensity. This step increases contrast by adjusting image intensity values. It suppresses eventually outliers pixels to improve contrast stretching. It then stretches contrast to use the whole range of intensity channel and if necessary it saturates some low or high values.
(8) Converting from YCbCr to RGB and reverse symmetric extension. After this step luminance channel has been preprocessed, so to regain colors we convert back the image the RGB space, and cut out the symmetric extension part of the image to recover the image with original size.
(9) Equalizing color mean. In underwater imaging color channels are rarely balanced correctly. This step enables to suppress predominant color by equalizing RGB channels means. It is rather used to produce a more pleasant image than to better segmentation. Because segmentation is in general performed on gray image and color equalization does really not change the gray image.

Uncalibrated Stereo Image Rectification
Given a pair of stereo images, rectification determines a transformation of each image plane such that pairs of conjugate epipolar lines become collinear and parallel to one of the image axes. The important advantage of rectification is that computing stereo correspondences is reduced to a 1-D search problem along the horizontal raster lines of the rectified images. In the case of uncalibrated cameras, there are more degrees of freedom in choosing the rectifying transformation and a few competing methods are present in the literature. Each aims at producing a "good" rectification by minimizing a measure of distortion, but none is clearly superior to the others, not to mention the fact that there is no agreement on what the distortion criterion should be. In this paper we adopt Quasi-Euclidean epipolar rectification method for uncalibrated images proposed by Andrea Fussiello et al. (2010) [4].
Geometrically, in the Euclidean frame, rectification is achieved by a suitable rotation of both image planes. The correspondent image transformation is the collineation induced by the plane at infinity. As a result, the plane at infinity is the locus of zero-disparity in the rectified stereo pair. This is signified by saying that Euclidean rectification is done with respect to the plane at infinity. In the uncalibrated case the reference plane is generic, as any plane can play the role of the infinity plane in the projective space. Our uncalibrated rectification can be seen as referred to a plane that approximates the plane at infinity.
We assume that intrinsic parameters are unknown and that a number of corresponding points j j r l mm  are available. The method seeks the collineations that make the original points satisfy the epipolar geometry of a rectified image pair.
The fundamental matrix of a rectified pair has a very specific form, namely it is the skew-symmetric matrix associated with the cross-product by the vector 1 (1,0,0) The left-hand side of (Eq. 12) is an algebraic error, i.e., it has no geometrical meaning, so we used instead the Sampson error that is a first order approximation of the geometric error. The matrix can be considered as the fundamental matrix between the original images, therefore, in our case, the squared Sampson error for the th j correspondence is defined as: where () i  is the th i component of the normalized vector. As this equation must hold for any j , one obtains a system of non-linear equations { 0} j E  in the unknown r H and l H . A least-squares solution can obtained with the Levenberg-Marquardt algorithm, but the way in which r H and l H are parameterized is crucial, and characterizes our approach with respect to the previous ones. We force the rectifying collineations to have the same structure as in the calibrated (Euclidean) case, i.e. to be collineations induced by the plane at infinity, namely 11 ,, r nr r or l nl l ol The old intrinsic parameters ( , ) [] u  , provided that the second and third row of nr K and nl K are the same. Hence it is not necessary to include the matrices nr K and nl K in the parameterization.
Each collineation depends in principle on five (intrinsic) plus three (rotation) unknown parameters. The rotation of one camera along its X-axis, however, can be eliminated. Consider the matrix 1 1 [ ] .
T T or r l ol Let r R and l R be the same matrices as r R and l R after premultiplying with an arbitrary (but the same for both) rotation matrix about the X-axis. It is easy to verify that 11 [ ] [ ] .
TT r l r l R u R R u R    Geometrically, this coincide with rotating a rectified pair around the baseline, which do not alter the rectification, but, in a real camera, it affects the portion of the scene that is imaged. Accordingly, we set to zero the rotation around the X-axis of the left camera.
We further reduce the number of parameters by making as educated guess on the old intrinsic parameters: no skew, principal point in the center of the image, aspect ratio equal to one. The only remaining unknowns are the focal lengths of both cameras. Assuming that they are identical and equal to  , we get: which varies in [ 1 , 1].


The minimization of the cost function is carried out using Levenburg-Marquardt, starting with all the unknown variables set to zero. When  converges outside the boundaries of the interval [ 1 , 1]  a random restart is attempted. If the problem persists the minimization is carried out with fixed 0.

 
Finally, the new intrinsic parameters ( nr K and nl K ) are set equal to the old ones: nr nl ol K K K , modulo a shift of the principal point, that might be necessary to center the rectified images in the customary image coordinate frame. Horizontal translation has no effect on the rectification, whereas vertical translation must be the same for both images.

Dense Stereo Matching
A large number of stereo matching algorithms exist and they can be classified into two main categories: local and global methods, according to the principle they are based upon. Other methods called cooperative algorithms use local and global approach at the same time. The difficulty is to choose an algorithm to perform a dense 3D reconstruction taking into account rendering, metrologic quality, computing speed and complexity of the scene. The taxonomy of [19] provides information about the overall performance of the principal algorithms (textureless regions, depth discontinuity regions, occluded regions). Finally, the graph-cut method gives excellent results, performing better in textureless areas and near discontinuities, and outperforming the other optimization methods. The major downsides are the relatively high computation time, and the need precisely tuned parameters, whose values are often image-dependent. This algorithm remains however a very interesting choice for our application, since the quality of the rendering process is a higher priority compared to execution time.
Roy S. and Cox I.J. [18] were the first ones to use this algorithm in the context of multi-camera stereovision. In order to explain the graph-cut method, we will concentrate on the case of graphs with only two terminals. Figure 1 show a simple example of a two terminal graph, which can be used to minimize an energy function on a 33  image with two labels. The two terminals are usually called source s, and sink t. They correspond to the set of labels (different depths) that can be assigned to pixels of the image.

Fig 1: Graph -Cut example
In the general case of graph-cut theory, the goal is to find a cut that has a minimum cost among all cuts, by minimizing an energy function.
Let function f be the disparity function associated to each pixel of an image. We search labeling f that minimizes the energy. To define this energy on a photo-consistence criterion (similarity between intensities of a pixel p in the first image and the pixel (p + fp) in the second image) called data term. A second term, called smoothness term, penalizes discontinuities between neighborhood pixels. Thus, the energy can be written as: (17) where term p D is the data term and ( , ) pq V f f is the smoothness term penalty between adjacent pixels.
In (kolmogorov, zabih), the energy minimization considers the input images symmetrically, handles visibility properly, and imposes spatial smoothness while preserving discontinuities.
The disparity map gives the dense correspondence map between the stereo images. Thus, the depth map is computed by triangulation with matched point pairs and camera parameters. The mathematical equation for estimation of depth map is Where, z is the depth, f is focal length of the camera, B is Baseline between two cameras and d is the disparity.
The disparity map gives us a dense correspondence map between the stereo images. Thus, the depth map is computed by triangulation with matched point is considered independently. Therefore smoothing the surface is important to obtain a spatial coherence. In this paper, a spatial coherence is achieved by filtering the depth map.

EXPERIMENTAL RESULTS
To see the effectiveness of the proposed approach on underwater images, experiment is done in a small water body. The scene includes several objects at distance [1m, 2m], near the corner of the water body. The depth at the bottom of the water body was approximately 2m. The stereo setup consists of two identical digital still cameras, which are Canon D10 water proof. The two cameras are installed on a vertical stereo mount and it is kept in water body. The Figure 2 shows the captured underwater stereo images. The underwater images are preprocessed by applying sequentially homomorphic filtering, wavelet denoising and anisotropic filtering techniques. The Figure 3 shows the preprocessed underwater images. The Figure 4 shows a result of uncalibrated rectification technique on left and right images. The Figure 5 shows the dense disparity map obtained with the implementation of the graph-cut method presented in (kolmogorov and zabih, 2002). Figure 6 shows the final reconstruction with texture mapping which can be used for 3D quantitative imaging. The dense depth map with texture mapping gives a 3D model with a good visual impression.

CONCLUSIONS
In this paper, we have presented a complete 3D surface reconstruction of small-scale underwater objects to supervised exploration of ocean floors. The first step concern with acquisition of stereo images by using stereo vision system. We preprocess the acquired stereo images sequentially by applying suitable filter. The preprocessed image is then rectified by using uncalibrated rectification technique. The rectified images are then used for stereo matching using Graph cut method. Finally, depth map is generated using triangulation technique. The experimental results shows that the proposed method reconstruct 3D surface of an underwater object accurately.