- Email: [email protected]
Direct shape from shading with improved rate of convergence
Pattern Recognition, Vol. 30, No. 3, pp. 353-365, 1997 © 1997 Pattern Recognition Society. Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved 0031-3203/97 $17.00+.00
Pergamon
PII: S0031-3203(96)00097-0
D I R E C T S H A P E F R O M S H A D I N G W I T H I M P R O V E D RATE OF C O N V E R G E N C E l W. P. CHEUNG,* C. K. LEE and K. C. LI Department of Electronic Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong (Received 11 September 1995; in revised form 20 June 1996; received for publication 10 July 1996) Abstract--In this paper, new direct shape from shading algorithms with improved rate of convergence for inclined and perpendicular light source are presented. The algorithms are mainly based on optimal dynamic system control with dynamic programming. Unlike the previous shading algorithms which first estimate the gradient of surface instead of depth using variational approaches, our algorithms provide a direct relative depth perception. Here the maximum uphill principle for inclined and perpendicular light source is derived instead of using the existing minimum downhill criteria in order to perceive the normal shape of the object. Also the convergence conditions for the algorithms are derived and proven under zero-initial condition. We show that the learning rate is linearly proportional to the minimum value of the input image in the discrete case. Furthermore, in order to increase the rate of convergence, the momentum factor is introduced. © 1997 Pattern Recognition Society. Published by Elsevier Science Ltd. Physics-based vision
3-D shape
Shape from shading
1. I N T R O D U C T I O N
Surface shape estimation of observed objects is one of the main goals in the field of three-dimensional (3-D) object recognition. The estimation method is intimately related to early visual processing. In the past few years, several methods of deriving 3-D (more precisely 2 ½-D) surface shape have been developed using different techniques and imaging devices. They are classified into direct and indirect methods. Direct or active methods recover depth information directly using range finder (1'2) and structure light (3>. Indirect methods deternfine the relative depth by cues extracted from gray level images of observed object. This class contains method such as shape from shading, defocusing and stereopsis. Since the direct or active techniques usually involve many complex external component setup, they may not be suitable for instance in object shape estimation. Hence many researchers have focused on the latter class of method. The shape from shading (4'5) methods make use of the change in image brightness of an object to determine the relative surface depth and can be effectively applied to smooth surface regions. On the other hand, the geometric stereo (6) method provides the surface heights at a set of sparse feature points which have high intensity variations, and the heights of other points are obtained by interpolations. Actually most of the shading techniques employ variational approaches which impose constraints such as smoothness, imaging equation and light conditions. As a result these techniques can only be applied to the object surface that satisfies the constraints. Also the computa* Author to whom correspondence should be addressed. 1 This work was supported by the Polytechnic Grant. 353
tional complexity of the resulting algorithms is high and increases as the number of these constraints increases. Useful alternatives for variational shading are local shading and linear shading suggested by Pentland (7's) and Tsai and Shah (9). These approaches do not require prior information on the observed scene, imaging geometry and regularization. The assumption is that the surface shape is locally spherical, and the inclined light source is not allowed. In order to improve the estimation accuracy other approaches are developed to combine the stereo vision and shading. Some of them (1°'11) use shading to strengthen the stereo conditions for surface interpolation and segmentation. The others (12-14) incorporate photometric and geometrical stereo constraints to the error minimization process which lead to more complex computation algorithms. Recently Oliensis and Dupuis (15) derived a direct shape reconstruction method which is based on the optimal control theory. However, the surface gradient iteration must satisfy the minimum downhill criteria which makes the reconstructed surface always locally concave if no other information such as singular point (the point where the intensity is maximum) depth is given. This is because the shading problem is basically ill-posed. Also the general light direction version is very complex and the rate of convergence is not mentioned. This paper is also based on the optimal control theory and dynamic programming. However, we develop the maximum uphill principle instead of the minimum downhill criteria for both perpendicular and inclined light sources to achieve direct (locally convex) shape estimation because in real life the higher gray level image point usually corresponds to the higher level. Here the iterative con-
W.P. CHEUNG et al.
354
ditions are much simpler and we show that under zeroinitial condition, the rate of convergence is approximately linearly proportional to the minimum value of the input image for our maximum uphill conditions. Also we include the momentum factor (16) to improve the overall rate of convergence. In fact the momentum factor is borrowed from the idea of learning rate improvement of the neural network (NN) and we treat the iterative surface reconstruction process as a NN learning process. Section 2 will present the derivation of maximum uphill direct shape from shading for both perpendicular and inclined light source. Next the addition of momentum factor is discussed. Finally both synthetic and real images are used to conduct the algorithm testing and comparison.
2. D I R E C T SHAPE F R O M SHADING USING DYNAMIC PROGRAMMING
In many cases, the light source is quite far from the object so that the orthogonal light source can be assumed. That is, the direction of the light source is equal to the viewing or camera axis and Fig. 1 shows the configuration. Let z be the relative distance of the object surface point to a reference plane. Then the gradient component of that point is defined as
Oz
or
I --
~/ x / p T + q2 + 1 '
(2)
where I is the intensity of the observed image, ~/the albedo constant, and the surface normal of the corresponding image point (x,y), n = ( - p , - q , 1) T. Rearranging equation (2) in eikonal form leads to 72 1 - ~ = -IVzl 2
V= 1
2.1. Perpendicular light source
Oz
n(x,y)n~ I = ~ [In(x,y)[[
(3)
and we define the left-hand side of equation (3) as V,, i.e.
In general we assumed that the surface is diffused rather than specula-type so that the Lambertain property can be used.
P=ox' q = ~
and the normalized light source directional vector is ns = (0, 0, 1) T. By the reflectance equation of a Lambertian surface we obtain
(1)
~72 12
(4)
and its value lies between zero and negative infinity. Therefore in order to apply this shading equation to real images, special care must be taken at the points where ~ e intensity is zero. As these points will cause V jumps to negative infinity and the surface estimation process will diverge. Actually there should not be any zero intensity points because these points correspond to the sharp change of the depth value which violates the smoothness assumption. Here we adapt the dynamic programming approach to transform equation (2) to iterative form. Consider the following control problem: a "particle" initially located at (xo,Yo) moves in the plane in response to the velocity control parameters u = dx/ds = 2, v = dy/ds = ~ with
I
3.
I
n,
---
(0,0,1)
Fig. 1. Source and viewing direction for the perpendicular light source configuration.
Direct shape from shading with improved rate of convergence initial values u(0) = k(0), v(0) = •(0). These values are then to be chosen at each time T to maximize the cost function defined as T
E(xo,Yo, T) =
u2 - v2) ds,
sup ~ / ( V (u(O),v(0))
(5)
0
where sup is the supremum operator. We assumed that E is a differentiable function of the starting point and by expanding equation (5) around this point we get
e(xo,Yo, T) =
[E(x((ST),y((ST),T- ~T)
sup (,(0)#(0)) ~T
+17
v/dsI
0
and
E(xo, Yo, T) =
1
sup (,(0)#(0))
( V - u(0) 2
[OE(x,y)7 2 [OE(x,y)] 2
7-x J +L
OE(x, y, T) v(O) ~ OE(x, y, T) 1
where E(x,y) = limT_~E(x,y, T) and equation (11) is equivalent to equation (3) with E(x,y) equal to z(x,y). Note that E _ < 0 and E = 0 (maximum value) at singular points which make the surface to be locally convex and allow the perception of the direct shape of the object. Since in real life, higher intensity points usually correspond to higher depth levels. The above analysis shows that the surface evolution governed by equation (10) can be derived by the dynamic minimization of the maximum error between the negative gradient magnitude (-[Vz] 2) and the image potential function V. In equation (10), E is the intermediate surface depth variable which converges to z as time T approaches to infinity.
OT
+ 0 ((ST) }.
(7)
E(x,y,k+l) = E ( x , y , k ) + A T [ C 2 x E ~ y ' k ) ) 2 +
As fiT ---+0 we have
OE(x, y, T) _ sup 1 [V - u(O) 2 - v(0) 2] or (~(0)#(0)) OE(x,y, T) u(O)-r OE(x, y, T) v(O)}. (8) Ox
Oy
For the maximization of equation (8), we choose
OE(x, y, T) Ox '
(1l)
For the discrete form of equation (10), Tis divided into time intervals and replaced by the running index, k and x and y become space index positions. Then equation (10) becomes
+ E(xo,Yo, T ) + e r [ O E ( o Y ' r ) u(O)
u(O) T
J =-v,
2.2. Shape reconstruction algorithm based on maximum uphill principle for perpendicular light source
v(0)2)~T
2
Oy
355
v(O) T _ OE(x, y, T) Oy
-+V,
(12)
where AT is a constant which is the rate of convergence. By substituting equation (4) to equation (12) and setting z~c = Ay = 1, we have the shape reconstruction iterafive equation
E(x,y,k + 1)
E(x,y,k) + AT[AxE(x,y,k) 2
(9)
Then we have Here simple backward or forward difference cannot be used because these difference approximations may not or 2 \ ~7 / ~: " satisfy the reflectance equation (2). Hence the x and y discrete directional slopes ( AxE(x, y, k ), AyE(x, y, k ) ) (10) are redefined in the following equations in order to meet As T goes to infinity, the time derivative of equation (10) equation (2) and convergence requirements (see Propovanishes which can be shown by replacing the integration sition 1):
OE(x,y,T) __ I [ v + ( O E ( x , y , T ) ~ 2 + ( O E (x~y,T))21
AxE(x,y,k) =
Ax+E(x,,y, k) if Ax+E(x,y,k ) if 2Xx+E(x, y, k) Ax-E(x,y,k) ma~(I Ax+E(x, y, k)l, I&-E(x, y, k)I) if Ax+E(x, y, k) 0 if 2x~+E(x,y,k) ( Ay+E(x, y, k)
AyE(x, y, kl = ~ Ay_E(x, y, kl
[ Oax(IAy+E(x,Y, k) l, lay E(x,y,k)l)
limit Twith infinity in equation (5) and removing all T variables in equations (6)-(9) which then leads t o
if if if if
Ay+E(x,y, k) Ay+E(x, y, k) Ay+E(x,y,k) /ky+E(x,y,k)
>_ 0 <0 >_ 0 <_ 0
and and and
Ax_E(x,y,k ) Ax_E(x, y, k) Ax E(x, y, k) and Ax E(x,y,k)
>_O, <_O, <_O, >_O,
(14)
>_ 0 <_ 0 >_ 0 <0
and and
Ay E(x,y, k) Ay_E(x, y, k) and Ay_E(x,y,k) and Ay_E(x,y,k)
>_O, <_O, < O, > O,
(15)
Ax+E(x,y,k), Ax E(x,y,k), Ay+E(x,y,k), and Ay E(x, y, k) are backward and forward differences in x
where
W.R CHEUNG et al.
356
E(x-l,y) e(x +l,y)
E(x, y) / ~ E(x-l,y) ~
E (x,y)
A_E(x,y) selected
~
selected
~z
E(x+1,y)
A +E(x,y)
a~_e&,y) case
case 1
2
E (x,y)
E(x+l,y)
~
x x
E(x- I,y) ~
//~ A,+E(x.y)
/
£(x- l,y) ~
selected
_
,Y)
~
E(x+l,y)
ArE = 0
E(x, y) case 3
case
4
Fig. 2. Maximum uphill principle. and y directions. Figure 2 illustrates the molecular representation of the modified gradient computation in the x-direction. The y-directional representation is the same with the interchange of x and y variable.
Proposition 1. The maximum uphill principle defined by equations (14) and (15) satisfy the brightness equation. Also the shape reconstruction equation (13) converges with the zero initial condition and AgE(x,y, k) <_0 Vk, where AkE(X,y, k) ~ e(x,y, k) - E(x,y, k - 1).
Proof. Consider
the reflectance equation (3) inside an infinitesimal area where the image intensity is constant. By differentiating it with respect to an arbitrary parameter ~5, it becomes
0 =p~+q~,
(16)
which implies that the slope dE/ds is either the steepest ascent or descent inside the area in the direction 6 = tan-l(q/p), where x = scos(0) and y = ssin(O). Hence the maximum uphill equations (14) and (15) satisfy the reflectance equation and we only consider four directions in the discrete case. Next we have to show that 0< -
min(x,y)AkE(x,y, k + 1) min(x,y)AgE(x, y, k)
< 1
Vx, y a n d k = l , 2 . . . (17)
for the convergence. From equation (13), the maximum uphill principle and AkE(x, y, k) < 0 we have
minlA~E(x,y.k + l)[ > minlAxE(x;y,k)[ (x,y) ' - (x,y)
(18)
and
minl2xyE(x,y,k+ (x,y)
1)l >
minlA,E(x,y,k)l, (x,y)
(19)
which implies that the minimum value with respect to x, y of AkE(x ,y,k) is increasing with k. Hence min(x,y)AgE(x,y, k + 1) ~ min(x,y)AkE(x, y, k) and together with AkE(x,y;k) <_O, equation (18) is satisfied which implies the convergence of equation (13). Notice that in order to have AgE(x,y, 1) < 0 under zero initial condition, ~ must be greater than/max. [] Equations (13)-(15) together with Proposition 1 provide the iterative scheme for surface reconstruction in the discrete form in which E is the intermediate surface depth variable. It will converge to z if the rate of convergence is chosen according to Proposition 2. Unlike many variational approaches, this scheme requires direct simple iterative equations and only the x-y discrete directional slopes need more complex computations. Figure 2 shows that the backward or forward x discrete directional slope value is chosen whenever it reflects as uphill direction with respect to the central element. If both are upward, the larger magnitude one will be taken. If both are downward, the final slope value becomes zero and this illustrates the maximum uphill principle.
Proposition 2.
The maximum rate of convergence of equation (13) is approximately linearly proportional to the minimum image intensity.
Proof. From Proposition 1, we must have AkE(x,y,k ) < 0, then under the zero-initial condition [E(x,y, 0) ~ 0] and from equation (13) at k = l , 2 and oo,
Direct shape from shading with improved rate of convergence we obtain
357
transform equation (22) to the form of equations (5) and (6), we let Z = Ixx + lzz.
(23)
Then we have
or
~2~/ 2
AT <
L i2),
ATm~ ~ .- ~ "
(24)
and
Hence min
OZ ~ x = lx + l z p
P
(20)
OZ Q = ~yy = lzq.
(25)
(21) Substituting equations (24) and (25) into equation (22), we get
[] In order to have the maximum rate of convergence, this value should be taken. Also when the image has zero minimum image intensity, the rate of convergence is also zero which matches our requirement that I cannot be zero as it will cause negative infinity in our iteration equation (13). As a result if the given image has a lower gray level content, proper positive offsetting should be given to improve the convergent rate. 2.3. Inclined light source Here we let the normalized light source directional vector be ns = (Ix, ly, lz) T and Fig. 3 depicts the configuration. By the reflectance equation of a Lambertian surface we obtain I = ~](lz - pIx - qly) ~/p2 + q2 4- 1
(22)
Without loss of generality, we let ly be zero, i.e. the light source is inclined along the x-axis. This can be done because once we know the direction of light source, we can rotate the image such that the new image corresponds to the light source with direction along x-axis. In order to
pZj_2 l - ~ - lxP+Q2+l
/~=0,
(26)
where J=l-(~£)
2.
(27)
Similar to the perpendicular light source, we introduce two control parameters governed by the following eqnation:
(9
±=Ju-lx
1-~-
and
y=v
and
v(0) = y ( 0 ) .
(28)
with initial values u(0) = 2 ( 0 ) + / x ( 1
-
~2/I2)
J These values are to be chosen at each time T to maximize the cost function defined as T
2f
sup 1 (,(0),v(0))
E(xo,Yo, T ) =
( v _ ju2 - vZ) ds,
0
I
sO( = qx@,tz)
-- (-p,-q, 1)
Fig. 3. Source and viewing direction for the inclined light source configuration.
(29)
W.R CHEUNG et al.
358
where x0 = x(0) and Y0 = y(0). Similar to dynamic programming analysis in Section 2.1, we have
OE(x,y,r) OT
1 v+
2
\
Ox
Ax+U(x,y,k)=Ax+E(x,y,k) -lx(1
-21x(1-~--2) (OE(xy'T)\ Ox +\
Ax U(x,y, k)
and v(0) r _ OE(x,y, T)
(31)
Oy
As Tgoes to infinity and with E equal to Z, equation (30) converges to equation (26). Furthermore from equation (23), the original shape z can be recovered from final Z values. The above analysis shows that the transformed surface (Z) governed by equation (30) can be derived by the dynamic minimization of the maximum error between the negative transformed gradient magnitude (-Ju 2 - va) and the image potential function V. In equation (30), E is the intermediate transformed surface depth variable which converges to Z as time Tapproaches to infinity. Also it is assumed that l~ and the shadow effect are small. In fact, the perpendicular light source shading analysis in Section 2.1 is the special case where Ix equals to zero.
2.4. Shape reconstruction algorithm based on maximum uphill principle for inclined light source
E(x,y,k + 1) = E(x,y,k) + ATIJAxE(x,y,k)2
'
J
Since we assume that the light source is inclined along xaxis only, equation (33) combined with equation (35) is the transformed version of equation (14) in which the maximum uphill operation is applied to discrete x-directional difference of U and they reduce to the perpendicular source case as 1~ becomes zero. The convergence requirements can be found in Propositions 3 and 4.
Proposition 3. The maximum uphill principle for the inclined light source defined by equations (33) and (34) satisfies the brightness equation (26). Also the shape reconstruction equation (32) converges for small lx with zero initial condition and AtE(x, y, k) <_O.
Proof. Similar to Proposition 1, consider the reflectance equation (26) inside an infinitesimal area where the image intensity is constant. By differentiating it with respect to an arbitrary parameter 4, it becomes
O=JP~-Ix (1_~)
~÷ Q~
or
~-
~d~,
(36)
which implies that the inclined directional slope
(
21~k,1 - ~)~72AxE(x , y, k)
¢]
+ AxE(x, y, k) 2 + 1 - ~- .
(32)
Similar to the perpendicular light source case, simple backward or forward difference cannot be used because these difference approximations may not satisfy the reflectance equation (26). Hence the x and y discrete directional slopes (AxE(x, y, k), AyE(x, y, k) ) are redefined in the following equation in order to meet the reflectance equation (26) and convergence requirements (see Proposition 3):
( Ax+E(x,, y, k) J Ax_E(x, y, k) AxE(x,y,k ) = [ 0 ax(IAx+E(x'y'k)l' lax E(x,y,k)l)
AyE(x, Y, k)
~2/12)
j
A~_E(x, y, k) - l~ (1 - ~2/~2)
For the discrete form of equation (30), it becomes
-
_
(35)
J
~y
u(O)r _ OE(x,y, T) Ox '
where Ax+E(x,y, k), Ax_E(x,y, k), Ay+E(x,y, k) and Ay_E(x, y, k) are the backward and forward differences in x and y directions, also
is either the steepest ascent or descent inside the area in the direction = tan-1
[~/p _
JQ
Hence the maximum uphill equations (33)-(35) satisfy the reflectance equation (26) and we only consider four directions in the discrete case. The convergence proof is similar tO Proposition 1 except that E is replaced by U in equation (18). [] if if if if
Ax+E(x,y, k) Ax+E(x,y, k) Ax+E(x, y, k) Ax+E(x, y, k)
Ay+E(x,y, k) if Ay+E(x,y, k) if Ay+E(x,y, k) Ay_E(x, y, k) max(lAy+E(x, Y, k) l, IAy-E(x, Y, k) l) if Ay+E(x,y, k) if Ay+E(x,y,k) 0
>_0 and <_0 and >_0 and <_0 and >0 <0 >0 <0
Ax_E(x, y, ~) >_O, Ax E(x,y,k) <_0, Ax_E(x, y, k) < O, Ax E(x,y,k) >_0,
and Ay E(x,y, k) and Ay E(x, y, k) and Ay E(x, y, k) and Ay E(x, y, k)
> 0, <_O, ~_ O, > O,
(33)
(34)
Direct shape from shading with improved rate of convergence Equations (32)-(35) together with Proposition 3 provide the iterative scheme for the inclined light source surface reconstruction in the discrete form. E is the intermediate surface depth variable and will converge to Z if the rate of convergence is chosen according to Proposition 4. The true z values can be recovered from equation (23).
Proposition
4. The m a x i m u m rate of convergence of equation (32), for small lx is approximately linearly proportional to the m i n i m u m image intensity.
Proof
The proof is similar to Proposition 2. From Proposition 3, we must have A~E(x, y, k) <_O, then for small lx and from equation (32) at k = l , 2 and oc, we obtain
359
Like the situation found in the perpendicular source, in order to have the maximum rate of convergence, this value should be taken. Also the input image should not have zero minimum image intensity which corresponds to abrupt jump of the surface or the shadow region.
3. I N T R O D U C T I O N OF THE M O M E N T U M FACTOR F O R THE INCREASE OF T H E RATE OF CONVERGENCE
In order to improve the rate of convergence, we modify the iterative equation (13) through the addition of a new term to the intermediate surface variable E. Hence equation (13) becomes
AkE(x, y, k + 1)
=
aAkE(x, y, k) + (1 -
a)Ar IAxE(x, y, k) 2
+/kEy(x,y, k) 2 +
1 --
~72/[(x,y)2],(39)
or
A T < ~/12 + (J + 1)[072/12) -- 1] -(J + 1)[(r/2/I 2) - 1]
lx
(37)
Hence
(38) []
where AkE(x,y,k+ 1) is the successive surface error and c~ ff [0, 1] is the momentum factor that determines the relative contribution of the current and past surface error term. Actually we borrow the momentum concept from the weight update of the neural network learning step (16) and we treat the iterative sequence as a learning process with the error identified as the last term of the right-hand side in equation (39). The m o m e n t u m factor
(a)
(b) 12 10
s
8 6
o< 4 2
,L
0 20
40 no, of i t e r a t i o n
60
80
(c)
Fig. 4. (a) The synthetic hemisphere image, (b) surface reconstruction and (c) the surface maximum % error curve (dash line indicates the effect of momentum).
360
W.R CHEUNG et al.
(a)
(b)
(c) 12 10 8 o
•
6 4 2 0
20
40
60
80
100
no. of iteration
(d)
Fig. 5. (a) The real David image, (b) surface reconstruction, (c) the rendered surface and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
can be considered as a means to increase the rate of convergence for the following reasons. When the consecutive surface error has the same sign, the exponential weighted sum of the surface error grows large in magnitude and the surface relative depth value is adjusted by a
large amount. Similarly, when the consecutive surface error has opposite signs, this sum becomes small in magnitude and the surface depth is adjusted by a small amount. On the other hand, the momentum factor can be considered as the second-order term to increase the roll-
Direct shape from shading with improved rate of convergence
(a)
361
(b)
(c)
6
5 ~_ 4 o~
~3
k,.. 2 1
0
0
2'0
4'0
step no. (d)
60
8'0
1 O0
Fig. 6. (a) The real pepper image, (b) surface reconstruction, (c) the rendered surface and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
off of the original surface error decay curve. For the case of inclined light source, the momentum factor is added to equation (32) which becomes
AkE(x, y, k + 1) = o~AkE(x, y, k) + (1 - ~ ) A r F Z X x E ( x , y, k) 2 - % ( ~ - , 2 / F ) A x E ( x , y, k)
+ AyE(x, y, k) 2 + 1 - q2/12].
(40)
4. EXPERIMENTALRESULTS In order to check the performance of our algorithms, both synthetic and real images are used. Although the derivation is quite tedious, direct iterative equation can be used to reconstruct the shape with zero-initial condition which enables us to have a fast implementation of surface shape estimation. Equations (13)-(15) and (39) are used to reconstruct the surface of a synthetic hemisphere, the real David and pepper images under perpendicular light source configuration with E(x, y, k) as the
362
W.R CHEUNG et al.
(b)
(a)
(c) 6 5 4 3 2
0
20
40
60 Mep
80
100
no.
Fig, 7. (a) The hemisphere image with 30 ° inclined light source, (b) surface reconstraction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations and (d) the surface maximum % error curve (dash line indicates the effect of momenttma). intermediate depth variable. They all have 128 x 128 resolution with 256 gray levels. Figure 4(b) shows the reconstructed surface of a synthetic hemisphere image [Fig. 4(a)]. Figure 5(b) shows the reconstructed surthce of the real David image [Fig. 5(a)] and Fig, 6(b) shows the surface of the real pepper image [Fig. 6(a)]. Since the momentum factor will not affect the final reconstructed shape, only the surfaces built from the momentum added iterative equation are displayed. Here all the images have singular points with a 255-intensity level, the albedo, r/is set to 255 in order to meet the convergence requirement in Proposition 1, Also the maximum rate of convergence
ATrnax for the synthetic hemisphere, the real David and pepper images are found using equation (21) in Proposition 2. They are 0,15, 0.08 and 0.2, respectively. Furthermore the momentum factor a is 0.2,(16)Figures 4(e), 5(d) and 6(d) show the decay curves of the maximum percentage error: rohn(xy)'
, , y, k) e/~
100 × min(x,y) A k E ( x ' Y, 1) for the surface reconstruction with and without momentum factor and the termination condition is 1% error. The dash line indicates the improvement of the conver-
Direct shape from shading with improved rate of convergence
iiii!i!iiiilN/I
:!i~ii1li!i~ i!i!i! ~ilili!ii! ;iiii!iiii•
~'~ii!ii!ii;iii1i ii iil i~!ilililil :~iiii~i
363
ii!iiii!iiiiii iiiiiii!ii!!',i liiilii~i ii~i '"i~":
ii
i iiiiiii!iiiiiiiiiiiii~i
i;iliiiii!iii~/ (b)
(a)
(c) 18 16 14
o; l o
12
8 6 4 2 0
20
4O step no.
(d)
6'0
0
Fig. 8. (a) The half cylinder image with 30° inclined light source, (b) source reconstruction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
gence rate to which the m o m e n t u m factor is added. Depending on the complexity of the input images, the improvement ranges approximately from 10 to 20 iterative steps. In order to test our inclined light source shape reconstruction, iterative equations (32)-(35) and (40), the synthetic hemisphere [Fig. 7(a)], the half cylinder [Fig. 8(a)] and the real China doll [Fig. 9(a)] images are used. Their light source is inclined along the x-direction. Both the hemisphere and half cylinder
show a 30 ° inclination [i.e. the source directional vector is (0.5000 0 0.8660)rl and the source directional vector of the China doll is estimated (17) to be (0.1963 0 0.9805) x. Figure 7(b) and (c), Fig. 8(b) and (c) and Fig. 9(b) and (c) compare the difference of the reconstructed surface using perpendicular and inclined source iterative equations. One may notice that the reconstructed surface will decline at the low gray level regions if the inclined source iterative equations are not
364
W.E CHEUNG et al.
(b)
(c)
(d) 40 35 30 o
25 20 15 10 5 0
20
40
60 step
80
100
no.
(e) Fig. 9. (a) The real China doll image with source directional vector (0.1963 0 0.9805) T, (b) surface reconstruction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations, (d) the rendered surface and (e) the surface maximum % error curve (dash line indicates the effect of momentum).
used. Similar to the perpendicular light source case, all the images have singular points with 255 intensity level, thus the albedo, ~7 is set of 255 in order to meet the convergence requirement in Proposition 3. The maximum rate of convergence ATmax for the synthetic hemisphere, the half cylinder and the China doll images is obtained using equation (38) in Proposition 4 as 0.15, 0.15 and 0.1, respectively. The momentum factor c~ is 0.2 Figs 7(d), 8(d) and 9(e) show the decay curves of maximum percentage error with and without momentum factor. However, it is assumed that the inclination of light source (lx) is small. For a large source inclination,
images may have shadow regions which violate the reflectance equation and our surface reconstruction scheme cannot apply. Under the perpendicular light source and without the momentum factor our error curve or convergence profile is similar to the existing method (15) because the only difference is the local convexity property of our method. For the inclined light source, the existing method is much more complex and the convergence result is not clearly indicated. However, with the momentum factor added, the error curve for the perpendicular and inclined light source of our method can be improved.
Direct shape from shading with improved rate of convergence 5. CONCLUSIONS In this paper, a n e w direct shape from shading m e t h o d with i m p r o v e d convergence rate for inclined and perpendicular light source is presented. The algorithms are mainly based on optimal dynamic s y s t e m control theory with d y n a m i c programming. Several lines o f instruction in the source code can i m p l e m e n t the surface shape reconstruction as the derived instead o f using the existing m i n i m u m downhill criteria in order to make the default surface shape to be locally convex for the perception o f the normal shape o f the object. Also the convergence conditions for the algorithms are derived and proven under zero-initial condition. We showed that the learning rate is linearly proportional to the m i n i m u m value o f the input image. Moreover, in order to increase the rate o f convergence, the m o m e n t u m factor is introduced. F r o m the test results o f synthetic and real images the improvem e n t ranges approximately from 10 to 20 iteration steps o f a total 100 to converge. Acknowledgements--The authors would like to thank the anonymous reviewer for the valuable comments.
REFERENCES
1. M. Rioux, Laser range finder based on synchronized scanners, Appl. Opt. 23, 3837-3844 (1984). 2. W. R Cheung, A laser imaging system for 3D leadframe inspection, M. Phil. Thesis, Hong Kong Polytechnic (1994).
365
3. J. A. Jalkio, R. C. Kim and S. K. Case, Three-dimensional inspection using multistripe structured light, Opt. Eng. 24, 966-974 (1985). 4. B. K. P. Horn and M. J. Brooks, Variational approach to shape from shading, CVG1P 33, 174-208 (1986). 5. B. K. P. Horn, Height and gradient from shading, Int. J. Comput. Vision 5, 584-598 (1990). 6. D. De Vleeschauwer, An intensity-based coarse-to-fine approach to reliably measure binocular disparity, CVGIP Image Understanding 57, 204-218 (1993). 7. A.P. Pentland, Local shading analysis, IEEE Trans. PAMI 6, 170-187 (1987). 8. A.P. Pentland, Linear shape from shading, Int. Z Comput. Vision 4, 153-162 (1990). 9. P. S. Tsai and M. Shah, Shape from shading using linear approximation, Image Vision Comput. 13, 377-383 (1995). 10. M. T. Chiaradia, A. Distance and E. Stella, Threedimensional surface reconstruction integrating shading and sparse stereo data, Opt. Eng. 28, 935-942 (1989). 11. W. E. L. Grimson, Binocular shading and visual surface reconstruction, CVGIP 28, 1943 (1984). 12. K. M. Lee and C. C. J. Kuo, Robust shape reconstruction from combined shading and stereo information, Proc. SPIE 1771, 171-182 (1992). 13. K. M. Lee and C. C. J. Kuo, Shape from shading with perspective projection, CVGIP 59, 202-212 (1994). 14. W. P. Cheung, C. K. Lee and K. C. Li, Surface reconstruction by integrating stereo vision with shading, Proc. IEEE 1ECON 2, 967-970 (1994). 15. J. Oliensis and P. Dupuis, Direct method for reconstructing shape and shading, Proc. SPIE 1570, 116-128 (1991). 16. R. A. Jacobs, Increased rates of convergence through learning rate adaptation, Neural Networks 1, 295-307 (1988). 17. Q. Zheng and R. Chellappa, Estimation of illumination direction albedo, and shape from shading, IEEE Trans. PAMI 13, 680-702 (1991).
About the A u t h o r - - W . E CHEUNG received his B.Eng. (Hons) and M.Phil. in ElecU'onic Engineering from
Hong Kong Polytechnic, in 1991 and 1994. He was a Teaching Company Associate of Hong Kong Polytechnic from 1991 to 1993. He is currently a Ph.D. candidate of The Hong Kong Polytechnic University in the Department of Electronic Engineering. His research interests are in 3D computer vision, image analysis and neural computation.
About the A u t h o r - - C . K. LEE is currently an Assistant Professor of the Hong Kong Polytechnic University in
the Department of Electronic Engineering. His research interests are in computer vision, motion control and neural network. He is a member of IEE and IEEE.
About the A u t h o r - - K . C. LI is currently an Associate Professor of The Hong Kong Polytechnic University in
the Department of Electronic Engineering. His research interests are in computer vision and telecommunication. He is a member of IEEE and Fellow of IEE.
Pergamon
PII: S0031-3203(96)00097-0
D I R E C T S H A P E F R O M S H A D I N G W I T H I M P R O V E D RATE OF C O N V E R G E N C E l W. P. CHEUNG,* C. K. LEE and K. C. LI Department of Electronic Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong (Received 11 September 1995; in revised form 20 June 1996; received for publication 10 July 1996) Abstract--In this paper, new direct shape from shading algorithms with improved rate of convergence for inclined and perpendicular light source are presented. The algorithms are mainly based on optimal dynamic system control with dynamic programming. Unlike the previous shading algorithms which first estimate the gradient of surface instead of depth using variational approaches, our algorithms provide a direct relative depth perception. Here the maximum uphill principle for inclined and perpendicular light source is derived instead of using the existing minimum downhill criteria in order to perceive the normal shape of the object. Also the convergence conditions for the algorithms are derived and proven under zero-initial condition. We show that the learning rate is linearly proportional to the minimum value of the input image in the discrete case. Furthermore, in order to increase the rate of convergence, the momentum factor is introduced. © 1997 Pattern Recognition Society. Published by Elsevier Science Ltd. Physics-based vision
3-D shape
Shape from shading
1. I N T R O D U C T I O N
Surface shape estimation of observed objects is one of the main goals in the field of three-dimensional (3-D) object recognition. The estimation method is intimately related to early visual processing. In the past few years, several methods of deriving 3-D (more precisely 2 ½-D) surface shape have been developed using different techniques and imaging devices. They are classified into direct and indirect methods. Direct or active methods recover depth information directly using range finder (1'2) and structure light (3>. Indirect methods deternfine the relative depth by cues extracted from gray level images of observed object. This class contains method such as shape from shading, defocusing and stereopsis. Since the direct or active techniques usually involve many complex external component setup, they may not be suitable for instance in object shape estimation. Hence many researchers have focused on the latter class of method. The shape from shading (4'5) methods make use of the change in image brightness of an object to determine the relative surface depth and can be effectively applied to smooth surface regions. On the other hand, the geometric stereo (6) method provides the surface heights at a set of sparse feature points which have high intensity variations, and the heights of other points are obtained by interpolations. Actually most of the shading techniques employ variational approaches which impose constraints such as smoothness, imaging equation and light conditions. As a result these techniques can only be applied to the object surface that satisfies the constraints. Also the computa* Author to whom correspondence should be addressed. 1 This work was supported by the Polytechnic Grant. 353
tional complexity of the resulting algorithms is high and increases as the number of these constraints increases. Useful alternatives for variational shading are local shading and linear shading suggested by Pentland (7's) and Tsai and Shah (9). These approaches do not require prior information on the observed scene, imaging geometry and regularization. The assumption is that the surface shape is locally spherical, and the inclined light source is not allowed. In order to improve the estimation accuracy other approaches are developed to combine the stereo vision and shading. Some of them (1°'11) use shading to strengthen the stereo conditions for surface interpolation and segmentation. The others (12-14) incorporate photometric and geometrical stereo constraints to the error minimization process which lead to more complex computation algorithms. Recently Oliensis and Dupuis (15) derived a direct shape reconstruction method which is based on the optimal control theory. However, the surface gradient iteration must satisfy the minimum downhill criteria which makes the reconstructed surface always locally concave if no other information such as singular point (the point where the intensity is maximum) depth is given. This is because the shading problem is basically ill-posed. Also the general light direction version is very complex and the rate of convergence is not mentioned. This paper is also based on the optimal control theory and dynamic programming. However, we develop the maximum uphill principle instead of the minimum downhill criteria for both perpendicular and inclined light sources to achieve direct (locally convex) shape estimation because in real life the higher gray level image point usually corresponds to the higher level. Here the iterative con-
W.P. CHEUNG et al.
354
ditions are much simpler and we show that under zeroinitial condition, the rate of convergence is approximately linearly proportional to the minimum value of the input image for our maximum uphill conditions. Also we include the momentum factor (16) to improve the overall rate of convergence. In fact the momentum factor is borrowed from the idea of learning rate improvement of the neural network (NN) and we treat the iterative surface reconstruction process as a NN learning process. Section 2 will present the derivation of maximum uphill direct shape from shading for both perpendicular and inclined light source. Next the addition of momentum factor is discussed. Finally both synthetic and real images are used to conduct the algorithm testing and comparison.
2. D I R E C T SHAPE F R O M SHADING USING DYNAMIC PROGRAMMING
In many cases, the light source is quite far from the object so that the orthogonal light source can be assumed. That is, the direction of the light source is equal to the viewing or camera axis and Fig. 1 shows the configuration. Let z be the relative distance of the object surface point to a reference plane. Then the gradient component of that point is defined as
Oz
or
I --
~/ x / p T + q2 + 1 '
(2)
where I is the intensity of the observed image, ~/the albedo constant, and the surface normal of the corresponding image point (x,y), n = ( - p , - q , 1) T. Rearranging equation (2) in eikonal form leads to 72 1 - ~ = -IVzl 2
V= 1
2.1. Perpendicular light source
Oz
n(x,y)n~ I = ~ [In(x,y)[[
(3)
and we define the left-hand side of equation (3) as V,, i.e.
In general we assumed that the surface is diffused rather than specula-type so that the Lambertain property can be used.
P=ox' q = ~
and the normalized light source directional vector is ns = (0, 0, 1) T. By the reflectance equation of a Lambertian surface we obtain
(1)
~72 12
(4)
and its value lies between zero and negative infinity. Therefore in order to apply this shading equation to real images, special care must be taken at the points where ~ e intensity is zero. As these points will cause V jumps to negative infinity and the surface estimation process will diverge. Actually there should not be any zero intensity points because these points correspond to the sharp change of the depth value which violates the smoothness assumption. Here we adapt the dynamic programming approach to transform equation (2) to iterative form. Consider the following control problem: a "particle" initially located at (xo,Yo) moves in the plane in response to the velocity control parameters u = dx/ds = 2, v = dy/ds = ~ with
I
3.
I
n,
---
(0,0,1)
Fig. 1. Source and viewing direction for the perpendicular light source configuration.
Direct shape from shading with improved rate of convergence initial values u(0) = k(0), v(0) = •(0). These values are then to be chosen at each time T to maximize the cost function defined as T
E(xo,Yo, T) =
u2 - v2) ds,
sup ~ / ( V (u(O),v(0))
(5)
0
where sup is the supremum operator. We assumed that E is a differentiable function of the starting point and by expanding equation (5) around this point we get
e(xo,Yo, T) =
[E(x((ST),y((ST),T- ~T)
sup (,(0)#(0)) ~T
+17
v/dsI
0
and
E(xo, Yo, T) =
1
sup (,(0)#(0))
( V - u(0) 2
[OE(x,y)7 2 [OE(x,y)] 2
7-x J +L
OE(x, y, T) v(O) ~ OE(x, y, T) 1
where E(x,y) = limT_~E(x,y, T) and equation (11) is equivalent to equation (3) with E(x,y) equal to z(x,y). Note that E _ < 0 and E = 0 (maximum value) at singular points which make the surface to be locally convex and allow the perception of the direct shape of the object. Since in real life, higher intensity points usually correspond to higher depth levels. The above analysis shows that the surface evolution governed by equation (10) can be derived by the dynamic minimization of the maximum error between the negative gradient magnitude (-[Vz] 2) and the image potential function V. In equation (10), E is the intermediate surface depth variable which converges to z as time T approaches to infinity.
OT
+ 0 ((ST) }.
(7)
E(x,y,k+l) = E ( x , y , k ) + A T [ C 2 x E ~ y ' k ) ) 2 +
As fiT ---+0 we have
OE(x, y, T) _ sup 1 [V - u(O) 2 - v(0) 2] or (~(0)#(0)) OE(x,y, T) u(O)-r OE(x, y, T) v(O)}. (8) Ox
Oy
For the maximization of equation (8), we choose
OE(x, y, T) Ox '
(1l)
For the discrete form of equation (10), Tis divided into time intervals and replaced by the running index, k and x and y become space index positions. Then equation (10) becomes
+ E(xo,Yo, T ) + e r [ O E ( o Y ' r ) u(O)
u(O) T
J =-v,
2.2. Shape reconstruction algorithm based on maximum uphill principle for perpendicular light source
v(0)2)~T
2
Oy
355
v(O) T _ OE(x, y, T) Oy
-+V,
(12)
where AT is a constant which is the rate of convergence. By substituting equation (4) to equation (12) and setting z~c = Ay = 1, we have the shape reconstruction iterafive equation
E(x,y,k + 1)
E(x,y,k) + AT[AxE(x,y,k) 2
(9)
Then we have Here simple backward or forward difference cannot be used because these difference approximations may not or 2 \ ~7 / ~: " satisfy the reflectance equation (2). Hence the x and y discrete directional slopes ( AxE(x, y, k ), AyE(x, y, k ) ) (10) are redefined in the following equations in order to meet As T goes to infinity, the time derivative of equation (10) equation (2) and convergence requirements (see Propovanishes which can be shown by replacing the integration sition 1):
OE(x,y,T) __ I [ v + ( O E ( x , y , T ) ~ 2 + ( O E (x~y,T))21
AxE(x,y,k) =
Ax+E(x,,y, k) if Ax+E(x,y,k ) if 2Xx+E(x, y, k) Ax-E(x,y,k) ma~(I Ax+E(x, y, k)l, I&-E(x, y, k)I) if Ax+E(x, y, k) 0 if 2x~+E(x,y,k) ( Ay+E(x, y, k)
AyE(x, y, kl = ~ Ay_E(x, y, kl
[ Oax(IAy+E(x,Y, k) l, lay E(x,y,k)l)
limit Twith infinity in equation (5) and removing all T variables in equations (6)-(9) which then leads t o
if if if if
Ay+E(x,y, k) Ay+E(x, y, k) Ay+E(x,y,k) /ky+E(x,y,k)
>_ 0 <0 >_ 0 <_ 0
and and and
Ax_E(x,y,k ) Ax_E(x, y, k) Ax E(x, y, k) and Ax E(x,y,k)
>_O, <_O, <_O, >_O,
(14)
>_ 0 <_ 0 >_ 0 <0
and and
Ay E(x,y, k) Ay_E(x, y, k) and Ay_E(x,y,k) and Ay_E(x,y,k)
>_O, <_O, < O, > O,
(15)
Ax+E(x,y,k), Ax E(x,y,k), Ay+E(x,y,k), and Ay E(x, y, k) are backward and forward differences in x
where
W.R CHEUNG et al.
356
E(x-l,y) e(x +l,y)
E(x, y) / ~ E(x-l,y) ~
E (x,y)
A_E(x,y) selected
~
selected
~z
E(x+1,y)
A +E(x,y)
a~_e&,y) case
case 1
2
E (x,y)
E(x+l,y)
~
x x
E(x- I,y) ~
//~ A,+E(x.y)
/
£(x- l,y) ~
selected
_
,Y)
~
E(x+l,y)
ArE = 0
E(x, y) case 3
case
4
Fig. 2. Maximum uphill principle. and y directions. Figure 2 illustrates the molecular representation of the modified gradient computation in the x-direction. The y-directional representation is the same with the interchange of x and y variable.
Proposition 1. The maximum uphill principle defined by equations (14) and (15) satisfy the brightness equation. Also the shape reconstruction equation (13) converges with the zero initial condition and AgE(x,y, k) <_0 Vk, where AkE(X,y, k) ~ e(x,y, k) - E(x,y, k - 1).
Proof. Consider
the reflectance equation (3) inside an infinitesimal area where the image intensity is constant. By differentiating it with respect to an arbitrary parameter ~5, it becomes
0 =p~+q~,
(16)
which implies that the slope dE/ds is either the steepest ascent or descent inside the area in the direction 6 = tan-l(q/p), where x = scos(0) and y = ssin(O). Hence the maximum uphill equations (14) and (15) satisfy the reflectance equation and we only consider four directions in the discrete case. Next we have to show that 0< -
min(x,y)AkE(x,y, k + 1) min(x,y)AgE(x, y, k)
< 1
Vx, y a n d k = l , 2 . . . (17)
for the convergence. From equation (13), the maximum uphill principle and AkE(x, y, k) < 0 we have
minlA~E(x,y.k + l)[ > minlAxE(x;y,k)[ (x,y) ' - (x,y)
(18)
and
minl2xyE(x,y,k+ (x,y)
1)l >
minlA,E(x,y,k)l, (x,y)
(19)
which implies that the minimum value with respect to x, y of AkE(x ,y,k) is increasing with k. Hence min(x,y)AgE(x,y, k + 1) ~ min(x,y)AkE(x, y, k) and together with AkE(x,y;k) <_O, equation (18) is satisfied which implies the convergence of equation (13). Notice that in order to have AgE(x,y, 1) < 0 under zero initial condition, ~ must be greater than/max. [] Equations (13)-(15) together with Proposition 1 provide the iterative scheme for surface reconstruction in the discrete form in which E is the intermediate surface depth variable. It will converge to z if the rate of convergence is chosen according to Proposition 2. Unlike many variational approaches, this scheme requires direct simple iterative equations and only the x-y discrete directional slopes need more complex computations. Figure 2 shows that the backward or forward x discrete directional slope value is chosen whenever it reflects as uphill direction with respect to the central element. If both are upward, the larger magnitude one will be taken. If both are downward, the final slope value becomes zero and this illustrates the maximum uphill principle.
Proposition 2.
The maximum rate of convergence of equation (13) is approximately linearly proportional to the minimum image intensity.
Proof. From Proposition 1, we must have AkE(x,y,k ) < 0, then under the zero-initial condition [E(x,y, 0) ~ 0] and from equation (13) at k = l , 2 and oo,
Direct shape from shading with improved rate of convergence we obtain
357
transform equation (22) to the form of equations (5) and (6), we let Z = Ixx + lzz.
(23)
Then we have
or
~2~/ 2
AT <
L i2),
ATm~ ~ .- ~ "
(24)
and
Hence min
OZ ~ x = lx + l z p
P
(20)
OZ Q = ~yy = lzq.
(25)
(21) Substituting equations (24) and (25) into equation (22), we get
[] In order to have the maximum rate of convergence, this value should be taken. Also when the image has zero minimum image intensity, the rate of convergence is also zero which matches our requirement that I cannot be zero as it will cause negative infinity in our iteration equation (13). As a result if the given image has a lower gray level content, proper positive offsetting should be given to improve the convergent rate. 2.3. Inclined light source Here we let the normalized light source directional vector be ns = (Ix, ly, lz) T and Fig. 3 depicts the configuration. By the reflectance equation of a Lambertian surface we obtain I = ~](lz - pIx - qly) ~/p2 + q2 4- 1
(22)
Without loss of generality, we let ly be zero, i.e. the light source is inclined along the x-axis. This can be done because once we know the direction of light source, we can rotate the image such that the new image corresponds to the light source with direction along x-axis. In order to
pZj_2 l - ~ - lxP+Q2+l
/~=0,
(26)
where J=l-(~£)
2.
(27)
Similar to the perpendicular light source, we introduce two control parameters governed by the following eqnation:
(9
±=Ju-lx
1-~-
and
y=v
and
v(0) = y ( 0 ) .
(28)
with initial values u(0) = 2 ( 0 ) + / x ( 1
-
~2/I2)
J These values are to be chosen at each time T to maximize the cost function defined as T
2f
sup 1 (,(0),v(0))
E(xo,Yo, T ) =
( v _ ju2 - vZ) ds,
0
I
sO( = qx@,tz)
-- (-p,-q, 1)
Fig. 3. Source and viewing direction for the inclined light source configuration.
(29)
W.R CHEUNG et al.
358
where x0 = x(0) and Y0 = y(0). Similar to dynamic programming analysis in Section 2.1, we have
OE(x,y,r) OT
1 v+
2
\
Ox
Ax+U(x,y,k)=Ax+E(x,y,k) -lx(1
-21x(1-~--2) (OE(xy'T)\ Ox +\
Ax U(x,y, k)
and v(0) r _ OE(x,y, T)
(31)
Oy
As Tgoes to infinity and with E equal to Z, equation (30) converges to equation (26). Furthermore from equation (23), the original shape z can be recovered from final Z values. The above analysis shows that the transformed surface (Z) governed by equation (30) can be derived by the dynamic minimization of the maximum error between the negative transformed gradient magnitude (-Ju 2 - va) and the image potential function V. In equation (30), E is the intermediate transformed surface depth variable which converges to Z as time Tapproaches to infinity. Also it is assumed that l~ and the shadow effect are small. In fact, the perpendicular light source shading analysis in Section 2.1 is the special case where Ix equals to zero.
2.4. Shape reconstruction algorithm based on maximum uphill principle for inclined light source
E(x,y,k + 1) = E(x,y,k) + ATIJAxE(x,y,k)2
'
J
Since we assume that the light source is inclined along xaxis only, equation (33) combined with equation (35) is the transformed version of equation (14) in which the maximum uphill operation is applied to discrete x-directional difference of U and they reduce to the perpendicular source case as 1~ becomes zero. The convergence requirements can be found in Propositions 3 and 4.
Proposition 3. The maximum uphill principle for the inclined light source defined by equations (33) and (34) satisfies the brightness equation (26). Also the shape reconstruction equation (32) converges for small lx with zero initial condition and AtE(x, y, k) <_O.
Proof. Similar to Proposition 1, consider the reflectance equation (26) inside an infinitesimal area where the image intensity is constant. By differentiating it with respect to an arbitrary parameter 4, it becomes
O=JP~-Ix (1_~)
~÷ Q~
or
~-
~d~,
(36)
which implies that the inclined directional slope
(
21~k,1 - ~)~72AxE(x , y, k)
¢]
+ AxE(x, y, k) 2 + 1 - ~- .
(32)
Similar to the perpendicular light source case, simple backward or forward difference cannot be used because these difference approximations may not satisfy the reflectance equation (26). Hence the x and y discrete directional slopes (AxE(x, y, k), AyE(x, y, k) ) are redefined in the following equation in order to meet the reflectance equation (26) and convergence requirements (see Proposition 3):
( Ax+E(x,, y, k) J Ax_E(x, y, k) AxE(x,y,k ) = [ 0 ax(IAx+E(x'y'k)l' lax E(x,y,k)l)
AyE(x, Y, k)
~2/12)
j
A~_E(x, y, k) - l~ (1 - ~2/~2)
For the discrete form of equation (30), it becomes
-
_
(35)
J
~y
u(O)r _ OE(x,y, T) Ox '
where Ax+E(x,y, k), Ax_E(x,y, k), Ay+E(x,y, k) and Ay_E(x, y, k) are the backward and forward differences in x and y directions, also
is either the steepest ascent or descent inside the area in the direction = tan-1
[~/p _
JQ
Hence the maximum uphill equations (33)-(35) satisfy the reflectance equation (26) and we only consider four directions in the discrete case. The convergence proof is similar tO Proposition 1 except that E is replaced by U in equation (18). [] if if if if
Ax+E(x,y, k) Ax+E(x,y, k) Ax+E(x, y, k) Ax+E(x, y, k)
Ay+E(x,y, k) if Ay+E(x,y, k) if Ay+E(x,y, k) Ay_E(x, y, k) max(lAy+E(x, Y, k) l, IAy-E(x, Y, k) l) if Ay+E(x,y, k) if Ay+E(x,y,k) 0
>_0 and <_0 and >_0 and <_0 and >0 <0 >0 <0
Ax_E(x, y, ~) >_O, Ax E(x,y,k) <_0, Ax_E(x, y, k) < O, Ax E(x,y,k) >_0,
and Ay E(x,y, k) and Ay E(x, y, k) and Ay E(x, y, k) and Ay E(x, y, k)
> 0, <_O, ~_ O, > O,
(33)
(34)
Direct shape from shading with improved rate of convergence Equations (32)-(35) together with Proposition 3 provide the iterative scheme for the inclined light source surface reconstruction in the discrete form. E is the intermediate surface depth variable and will converge to Z if the rate of convergence is chosen according to Proposition 4. The true z values can be recovered from equation (23).
Proposition
4. The m a x i m u m rate of convergence of equation (32), for small lx is approximately linearly proportional to the m i n i m u m image intensity.
Proof
The proof is similar to Proposition 2. From Proposition 3, we must have A~E(x, y, k) <_O, then for small lx and from equation (32) at k = l , 2 and oc, we obtain
359
Like the situation found in the perpendicular source, in order to have the maximum rate of convergence, this value should be taken. Also the input image should not have zero minimum image intensity which corresponds to abrupt jump of the surface or the shadow region.
3. I N T R O D U C T I O N OF THE M O M E N T U M FACTOR F O R THE INCREASE OF T H E RATE OF CONVERGENCE
In order to improve the rate of convergence, we modify the iterative equation (13) through the addition of a new term to the intermediate surface variable E. Hence equation (13) becomes
AkE(x, y, k + 1)
=
aAkE(x, y, k) + (1 -
a)Ar IAxE(x, y, k) 2
+/kEy(x,y, k) 2 +
1 --
~72/[(x,y)2],(39)
or
A T < ~/12 + (J + 1)[072/12) -- 1] -(J + 1)[(r/2/I 2) - 1]
lx
(37)
Hence
(38) []
where AkE(x,y,k+ 1) is the successive surface error and c~ ff [0, 1] is the momentum factor that determines the relative contribution of the current and past surface error term. Actually we borrow the momentum concept from the weight update of the neural network learning step (16) and we treat the iterative sequence as a learning process with the error identified as the last term of the right-hand side in equation (39). The m o m e n t u m factor
(a)
(b) 12 10
s
8 6
o< 4 2
,L
0 20
40 no, of i t e r a t i o n
60
80
(c)
Fig. 4. (a) The synthetic hemisphere image, (b) surface reconstruction and (c) the surface maximum % error curve (dash line indicates the effect of momentum).
360
W.R CHEUNG et al.
(a)
(b)
(c) 12 10 8 o
•
6 4 2 0
20
40
60
80
100
no. of iteration
(d)
Fig. 5. (a) The real David image, (b) surface reconstruction, (c) the rendered surface and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
can be considered as a means to increase the rate of convergence for the following reasons. When the consecutive surface error has the same sign, the exponential weighted sum of the surface error grows large in magnitude and the surface relative depth value is adjusted by a
large amount. Similarly, when the consecutive surface error has opposite signs, this sum becomes small in magnitude and the surface depth is adjusted by a small amount. On the other hand, the momentum factor can be considered as the second-order term to increase the roll-
Direct shape from shading with improved rate of convergence
(a)
361
(b)
(c)
6
5 ~_ 4 o~
~3
k,.. 2 1
0
0
2'0
4'0
step no. (d)
60
8'0
1 O0
Fig. 6. (a) The real pepper image, (b) surface reconstruction, (c) the rendered surface and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
off of the original surface error decay curve. For the case of inclined light source, the momentum factor is added to equation (32) which becomes
AkE(x, y, k + 1) = o~AkE(x, y, k) + (1 - ~ ) A r F Z X x E ( x , y, k) 2 - % ( ~ - , 2 / F ) A x E ( x , y, k)
+ AyE(x, y, k) 2 + 1 - q2/12].
(40)
4. EXPERIMENTALRESULTS In order to check the performance of our algorithms, both synthetic and real images are used. Although the derivation is quite tedious, direct iterative equation can be used to reconstruct the shape with zero-initial condition which enables us to have a fast implementation of surface shape estimation. Equations (13)-(15) and (39) are used to reconstruct the surface of a synthetic hemisphere, the real David and pepper images under perpendicular light source configuration with E(x, y, k) as the
362
W.R CHEUNG et al.
(b)
(a)
(c) 6 5 4 3 2
0
20
40
60 Mep
80
100
no.
Fig, 7. (a) The hemisphere image with 30 ° inclined light source, (b) surface reconstraction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations and (d) the surface maximum % error curve (dash line indicates the effect of momenttma). intermediate depth variable. They all have 128 x 128 resolution with 256 gray levels. Figure 4(b) shows the reconstructed surface of a synthetic hemisphere image [Fig. 4(a)]. Figure 5(b) shows the reconstructed surthce of the real David image [Fig. 5(a)] and Fig, 6(b) shows the surface of the real pepper image [Fig. 6(a)]. Since the momentum factor will not affect the final reconstructed shape, only the surfaces built from the momentum added iterative equation are displayed. Here all the images have singular points with a 255-intensity level, the albedo, r/is set to 255 in order to meet the convergence requirement in Proposition 1, Also the maximum rate of convergence
ATrnax for the synthetic hemisphere, the real David and pepper images are found using equation (21) in Proposition 2. They are 0,15, 0.08 and 0.2, respectively. Furthermore the momentum factor a is 0.2,(16)Figures 4(e), 5(d) and 6(d) show the decay curves of the maximum percentage error: rohn(xy)'
, , y, k) e/~
100 × min(x,y) A k E ( x ' Y, 1) for the surface reconstruction with and without momentum factor and the termination condition is 1% error. The dash line indicates the improvement of the conver-
Direct shape from shading with improved rate of convergence
iiii!i!iiiilN/I
:!i~ii1li!i~ i!i!i! ~ilili!ii! ;iiii!iiii•
~'~ii!ii!ii;iii1i ii iil i~!ilililil :~iiii~i
363
ii!iiii!iiiiii iiiiiii!ii!!',i liiilii~i ii~i '"i~":
ii
i iiiiiii!iiiiiiiiiiiii~i
i;iliiiii!iii~/ (b)
(a)
(c) 18 16 14
o; l o
12
8 6 4 2 0
20
4O step no.
(d)
6'0
0
Fig. 8. (a) The half cylinder image with 30° inclined light source, (b) source reconstruction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations and (d) the surface maximum % error curve (dash line indicates the effect of momentum).
gence rate to which the m o m e n t u m factor is added. Depending on the complexity of the input images, the improvement ranges approximately from 10 to 20 iterative steps. In order to test our inclined light source shape reconstruction, iterative equations (32)-(35) and (40), the synthetic hemisphere [Fig. 7(a)], the half cylinder [Fig. 8(a)] and the real China doll [Fig. 9(a)] images are used. Their light source is inclined along the x-direction. Both the hemisphere and half cylinder
show a 30 ° inclination [i.e. the source directional vector is (0.5000 0 0.8660)rl and the source directional vector of the China doll is estimated (17) to be (0.1963 0 0.9805) x. Figure 7(b) and (c), Fig. 8(b) and (c) and Fig. 9(b) and (c) compare the difference of the reconstructed surface using perpendicular and inclined source iterative equations. One may notice that the reconstructed surface will decline at the low gray level regions if the inclined source iterative equations are not
364
W.E CHEUNG et al.
(b)
(c)
(d) 40 35 30 o
25 20 15 10 5 0
20
40
60 step
80
100
no.
(e) Fig. 9. (a) The real China doll image with source directional vector (0.1963 0 0.9805) T, (b) surface reconstruction from perpendicular source iterative equations, (c) surface reconstruction using inclined source iterative equations, (d) the rendered surface and (e) the surface maximum % error curve (dash line indicates the effect of momentum).
used. Similar to the perpendicular light source case, all the images have singular points with 255 intensity level, thus the albedo, ~7 is set of 255 in order to meet the convergence requirement in Proposition 3. The maximum rate of convergence ATmax for the synthetic hemisphere, the half cylinder and the China doll images is obtained using equation (38) in Proposition 4 as 0.15, 0.15 and 0.1, respectively. The momentum factor c~ is 0.2 Figs 7(d), 8(d) and 9(e) show the decay curves of maximum percentage error with and without momentum factor. However, it is assumed that the inclination of light source (lx) is small. For a large source inclination,
images may have shadow regions which violate the reflectance equation and our surface reconstruction scheme cannot apply. Under the perpendicular light source and without the momentum factor our error curve or convergence profile is similar to the existing method (15) because the only difference is the local convexity property of our method. For the inclined light source, the existing method is much more complex and the convergence result is not clearly indicated. However, with the momentum factor added, the error curve for the perpendicular and inclined light source of our method can be improved.
Direct shape from shading with improved rate of convergence 5. CONCLUSIONS In this paper, a n e w direct shape from shading m e t h o d with i m p r o v e d convergence rate for inclined and perpendicular light source is presented. The algorithms are mainly based on optimal dynamic s y s t e m control theory with d y n a m i c programming. Several lines o f instruction in the source code can i m p l e m e n t the surface shape reconstruction as the derived instead o f using the existing m i n i m u m downhill criteria in order to make the default surface shape to be locally convex for the perception o f the normal shape o f the object. Also the convergence conditions for the algorithms are derived and proven under zero-initial condition. We showed that the learning rate is linearly proportional to the m i n i m u m value o f the input image. Moreover, in order to increase the rate o f convergence, the m o m e n t u m factor is introduced. F r o m the test results o f synthetic and real images the improvem e n t ranges approximately from 10 to 20 iteration steps o f a total 100 to converge. Acknowledgements--The authors would like to thank the anonymous reviewer for the valuable comments.
REFERENCES
1. M. Rioux, Laser range finder based on synchronized scanners, Appl. Opt. 23, 3837-3844 (1984). 2. W. R Cheung, A laser imaging system for 3D leadframe inspection, M. Phil. Thesis, Hong Kong Polytechnic (1994).
365
3. J. A. Jalkio, R. C. Kim and S. K. Case, Three-dimensional inspection using multistripe structured light, Opt. Eng. 24, 966-974 (1985). 4. B. K. P. Horn and M. J. Brooks, Variational approach to shape from shading, CVG1P 33, 174-208 (1986). 5. B. K. P. Horn, Height and gradient from shading, Int. J. Comput. Vision 5, 584-598 (1990). 6. D. De Vleeschauwer, An intensity-based coarse-to-fine approach to reliably measure binocular disparity, CVGIP Image Understanding 57, 204-218 (1993). 7. A.P. Pentland, Local shading analysis, IEEE Trans. PAMI 6, 170-187 (1987). 8. A.P. Pentland, Linear shape from shading, Int. Z Comput. Vision 4, 153-162 (1990). 9. P. S. Tsai and M. Shah, Shape from shading using linear approximation, Image Vision Comput. 13, 377-383 (1995). 10. M. T. Chiaradia, A. Distance and E. Stella, Threedimensional surface reconstruction integrating shading and sparse stereo data, Opt. Eng. 28, 935-942 (1989). 11. W. E. L. Grimson, Binocular shading and visual surface reconstruction, CVGIP 28, 1943 (1984). 12. K. M. Lee and C. C. J. Kuo, Robust shape reconstruction from combined shading and stereo information, Proc. SPIE 1771, 171-182 (1992). 13. K. M. Lee and C. C. J. Kuo, Shape from shading with perspective projection, CVGIP 59, 202-212 (1994). 14. W. P. Cheung, C. K. Lee and K. C. Li, Surface reconstruction by integrating stereo vision with shading, Proc. IEEE 1ECON 2, 967-970 (1994). 15. J. Oliensis and P. Dupuis, Direct method for reconstructing shape and shading, Proc. SPIE 1570, 116-128 (1991). 16. R. A. Jacobs, Increased rates of convergence through learning rate adaptation, Neural Networks 1, 295-307 (1988). 17. Q. Zheng and R. Chellappa, Estimation of illumination direction albedo, and shape from shading, IEEE Trans. PAMI 13, 680-702 (1991).
About the A u t h o r - - W . E CHEUNG received his B.Eng. (Hons) and M.Phil. in ElecU'onic Engineering from
Hong Kong Polytechnic, in 1991 and 1994. He was a Teaching Company Associate of Hong Kong Polytechnic from 1991 to 1993. He is currently a Ph.D. candidate of The Hong Kong Polytechnic University in the Department of Electronic Engineering. His research interests are in 3D computer vision, image analysis and neural computation.
About the A u t h o r - - C . K. LEE is currently an Assistant Professor of the Hong Kong Polytechnic University in
the Department of Electronic Engineering. His research interests are in computer vision, motion control and neural network. He is a member of IEE and IEEE.
About the A u t h o r - - K . C. LI is currently an Associate Professor of The Hong Kong Polytechnic University in
the Department of Electronic Engineering. His research interests are in computer vision and telecommunication. He is a member of IEEE and Fellow of IEE.