tween 
s. But 
ion of 
lation 
often 
The 
in that 
> from 
of the 
g the 
> data, 
an be 
an be 
>-scan 
ed by 
> data 
2. SHAPE FROM SHADING USING 
SIDE-SCAN SONAR IMAGES 
2.1 Inverse Problem and Regularization 
By means of an echo sounding system, a 3D seafloor surface 
Z(x,y) can be projected onto a 2D image with its intensity 
I(x,y) representing strength of reflected acoustic signals from 
the seafloor. Generally, the intensity of an image is described 
by a reflectance map R (also called backscattering model) 
(Kober and Leberl, 1991): 
I(x,y) =R(Z, Z, B, L, 9), (1) 
where 
Z,, Z, = slopes of Z(x,y) in x and y direction, 
P 7 illumination vector, 
L = viewing vector, and 
o, albedo. 
The shape from shading technique analyzes image intensity 
and computes slopes (Z,, Z,) The surface model of the 
seafloor Z(x,y) can be derived tying the seafloor surface 
slopes to the existing known depth points. In comparison to 
the forward calculation of the reflected acoustic signal 
represented by (1), the shape from shading leads to a 
mathematic inverse problem which cannot be simply solved. 
In this inverse problem, the sea floor surface is reconstructed 
from 2D image intensity. But the existence, the uniqueness, 
and the stability of the solution cannot be guaranteed without 
additional constraints. In this sense this inverse problem is 
mathematically ill-posed (Terzopoulos, 1986a, March, 1988, 
Lee and Pavlidis, 1988, Li, 1990). A well posed problem has a 
unique solution which depends on input data. The 
regularization technique transforms an ill-posed problem into 
a well-posed problem. In the inverse problem of surface 
reconstruction additional constraints are often applied in order 
to solve an ill-posed problem. The basic theory of 
regularization can be found in (Terzopoulos, 1986a, Tikhonov 
and Arsenin, 1977, Poggio et. al, 1985). 
The following objective function (also called cost function) is 
usually used for shape from shading technique: 
£- Jf (IG y) - RZ. Zu) 4 (22,322*,*Z^.)) dxdy, (2) 
where Z,, Z, and Z are second derivatives and A is a 
regularization parameter; I, R, Z, and Z, are the same as in 
(1). The first term in (2) represents a least square requirement 
which forces the difference between the image intensity I(x, y) 
and the calculated image function to a minimum by adjusting 
surface slopes (Z,, Z,). The term [[I(x, y) - R(Z,, Z,)P dxdy 
is also called a penalty functional because it increases if 
unreasonable values for the surface slopes (Z, Z,) are 
selected; furthermore € increases. This is a penalty in the sense 
of optimization because € should be minimized. 
The second term in (2) is mean curvature constraint of the 
reconstructed surface; that is, it adds a smoothness constraint 
to the least square requirement in the first term. Therefore, 
781 
equation (2) also represents a damped least squares. Since the 
smoothness constraint stabilizes the solution in the sense of 
optimization, analog to the penalty functional, the term [| 
(Z^,42Z?,,*Z? .) dxdy is also called stabilizing functional. 
The regularization parameter A adjusts the relation between 
the least square requirement and the smoothness constraint and 
it controls the roughness of the reconstructed surface model. It 
is clear if a large value of A is selected the algorithm will 
produce a smoothed surface which is not always true in 
accordance with seafloor topography. An appropriate choice 
of A can lead to surface slope estimates which approach 
seafloor topography properly and supply a reflectance map R 
close to the image intensity I(x, y). 
The solution of equation (2) by variation theory results in an 
Euler-Lagrange equation: 
À AZ, + [I(x, y) -R(Z, Z)] à -0 
A AZ + [I(x, y) - RZ, Z,)] 3 =O. 3) 
Here AZ, and AZ, are differentials of Z, and Z , and +x and 
a are partial derivatives of R with respect to Z, and Z,. A 
» " 
discrete form of (3) is (Kober and Leberl, 1991, Frankot and 
Chellappa, 1987): 
  
Z, aij) = Z, (ij) * 3/103) (IG, j) - R,,. Z,, i, )] 22 
Za. j) 7 Z, i, j) * 3/103) [IG, j) - R(Z,, Z,, i, 3] = 
4) 
with 
Z, (x,y) = [Z,,G+1,) + Z, i, j-1) * 
Z, (i-1, j) * Z, (i, j+1)] / 4, 
Z, (x, y) 7 IZ, (i, j) * Z, i, j-1) + 
Zl j) + ZG, j*1] / 4, 
Z 
n,x? 
Z,y 7 slopes of Z(x, y) in x and y direction 
in the nth iteration. 
Equation (4) is actually a relaxation formula, where slope 
estimates in (n+1)th iteration (Z,,,,, Z,,,,) can be calculated 
by (Z,,, Z,,) from the nth iteration. Therefore, computational 
reconstruction of the seafloor surface will be realized by an 
iterative procedure where the surface slopes are iteratively 
improved during the relaxation procedure. 
2.2 Improvement of Seafloor Surface Models 
Bathymetric data are often available in gridded format. That 
means, for every point (x, yj in a grid there is a 
corresponding depth z, Depths between grid points may be 
calculated by interpolation methods. In this case the 
interpolated depths are estimates made by considering 
neighboring depth points. Thus, the density of the grid