thm (DPM) for 
odel. Some small 
or determining the 
rithm, five small 
ed at the centre of 
[he other four are 
pper-right, lower- 
image and height 
gned as candidate 
je training regions 
or every candidate 
ning reflectance 
gorithm (DHI) for 
t image (see the 
uct an iterative 
ximation. A fixed 
frame. After the 
> model, a set of 
rs) with respect to 
1pdated improved 
ined. Using these 
ht data, a shading 
ly shaded regions. 
every candidate 
] region and the 
e candidate model 
selected as the 
1,---,m be the i th 
c models. Let 
parameters of the 
om DRP after the 
dated gradients in 
the iterations. The 
k, je») , 
letric modeli. (1) 
°d training regions 
ns is 
ietric modeli. (2) 
model MR is thus 
3) 
n fact, it is similar 
) 
algorithms. The 
ctance properties. 
1pdated improved 
ving algorithm to 
RP. In DHI, we 
. 
  
By using the photometric model 3 determined in the training 
frame, the known light source, with the brightness image and 
the approximate height image, we first estimate the surface 
reflectance property images in DRP. The number of the 
reflectance property images depends on how many parameters 
are used in the determined photometric model. E.g., if itis the 
Torrance-Sparrow model, there are three reflectance property 
images (the diffusion parameter image, the specular parameter 
image and the surface roughness parameter image). The value 
of each pixel in the property image represents the value of 
surface reflectance property at its coordinates. If it is the 
Lambertian model there is only one image containing albedo 
values. Then these estimated reflectance properties are used in 
DHI to improve the approximate height image. The improved 
height image is then used to improve the estimation of the 
reflectance property image(s), and so on. Obviously, this is an 
iteration strategy. In each iteration, a shading algorithm is 
called to obtain the shaded image. The mean square error 
between the brightness image and the shaded image is 
calculated to determine whether the iteration should be stopped. 
If the mean square error does not reduce, the program 
terminates and outputs the determined photometric model, the 
surface reflectance property image(s) and the improved height 
image. 
3.2.1 DRP - A Region Growing Algorithm to Determine 
Surface Reflectance Properties 
So far, our problem becomes to obtain the surface reflectance 
properties based on a given photometric model, a given 
brightness image and a given height image. Though this seems 
similar with the problem in Ikeuchi and Sato's paper (1991) 
and the problem in Kay and Caelli’s paper (1994), there are 
significant differences with their cases. We estimate parameters 
at each point on the object surface. Our method differs from 
that of Ikeuchi and Sato which assumed constant regions. Also, 
our method differs from that of Kay and Caelli which used 
multiple brightness images obtained with different light 
sources. 
A region growing algorithm based on least square fitting is 
proposed to estimate the inhomogeneous reflectance properties. 
In the algorithm, a set of initial reflectance properties are 
estimated in a set of kernel surface points. Surface reflectance 
properties are obtained by minimizing 
» 2 
e, 7 X [e»-ho»]. (4) 
kr zyen 
where Q, is the set of kernel points, / is the brightness image, 
Iy is the image irradiance corresponding to the photometric 
model R. For the Lambertian model, Eq. 4 becomes 
2 
e n Y [16.3) - ,:S-NG«»)] : (5) 
* xyef, 
where c, is the Lambertian diffuse reflection coefficient. S is 
the direction of light source, N is the surface normal. For the 
Torrance-Sparrow model, Eq. 4 becomes 
2 
e, = 2 [es - «4:8: NG) 7e, npo oan! j 
(6) 
x x,yef2, 
where c, is the specular coefficient, c, is the surface roughness 
International Archives of Photogrammetry and Remote Sensing. Vol. XXXI, Part B3. Vienna 1996 
  
coefficient, « is the specular angle calculated as 
ox, y)= cos”! (H-N(x,y)), (7) 
where H=(S+V)/[s+V|| with v the viewing direction. 
The minimization of Eq. 5 to obtain c, is simple. To minimize 
vanités and c,, we follow Kay and Caelli’s 
method. From the kernel set, we grow the neighbours of the 
kernel set if the neighbours have the same properties as the 
kernel set. Whether a neighbour point is to be grown depends 
on the error after growing. Let Q, «1 denote the new set of the 
Eq. 6 to obtain c 
kernel set 2, adding a neighbour point. Let ¢; 5, ¢, q, and 
c, g, be three parameters estimated in set © . The error ea 
Fale k p +1 
is 
Ca, +1 = Y (5) 7 &4,0, :S-N(x,y)- 
x,yef2, t1 
Cr, P(g, G5) )] 2 (8) 
«ih, ©) 
€0 +1 
& 
  
  
this neighbour point is added to the kernel set. In Eq. 9, th isa 
predefined threshold. Therefore, the new kernel set has one 
point more. After the neighbours of the old kernel set Q, are 
tested and the kernel set has grown, parameters c, , c, and c, 
are re-estimated and updated. If the kernel set has not grown, 
the growing for set Q, is stopped. The growing will start with 
a new kernel set until all points on the surface are estimated. 
The initial kernel set can be a small nxn window anywhere. 
However, the deviation of [70 »-lycy]. must be small by 
using the estimated parameters in the initial set to guarantee 
homogeneity. In the end, points having not been grown are 
interpolated by the successfully estimated neighbours. 
3.2.2 DHI Algorithm 
In fact, DHI is an SFS algorithm to obtain the improved height 
image. We mainly follow Zheng and Chellappa's algorithm 
(1991). Their algorithm is a constrained optimization problem 
inimiving 
f[F(p.q.Z) dx dy, (10) 
where 
F=[R(p,q)- 10x,” + 
[Rp (p.9) Px + Ra(P9) 9x = T(x)” + 
[R,(P.q) Py * Ry (CP. 4):4y - 1G Y) + 
np - Z,Y *(- 2,» Y. (11) 
In their algorithm, the intensity gradient constraint and the 
height gradient constraint are applied. Using the calculus of 
variations, minimization of Eq. 10 is equivalent to solving the 
Euler equations: 
1031