Vector Field Restoration By The Method Of Convex ... - Semantic Scholar
Vector Field Restoration By The Method Of Convex Projections Patrice Y. Simardand Guy E. Maillouxy June 15, 2000 Abstract
In this paper, the theory of image restoration by projections onto closed convex sets is applied to the restoration of vector elds. A set of useful projection operators is presented together with a linear time numerical implementation. These projection operators can be used to restore from partial information the velocity or deformation elds computed between successive views of a scene. They also nd applications in the restoration of vector elds of physical quantities as those encountered in mechanics, hydrodynamics or electromagnetism. The method is compared with the variational approach and illustrated by restoring simulated vector elds.
P. Y. Simard was with the Institut de G enie Biomedical Ecole Polytechnique, C. P. 6079, Succ A, Montreal, P.Q., Canada, H3C 3A7. He is now with the Department of Computer Science, University of Rochester, Rochester, NY 14627. y G. E. Mailloux is with the Institut de G enie Biomedical Ecole Polytechnique, C. P. 6079, Succ A, Montreal, P.Q., Canada, H3C 3A7.
1
List of symbols caligraphic H lambda in nity nabla integral partial norm in summation delta bigcup longarrow bigcap mu
H
1 r R @
k 2P S = T)
2
1 Introduction The aim of image restoration is to make as good an estimate as possible of an original image subjected to some form of degradation. Usually, images are restored by inverse ltering or constrained deconvolution techniques. These techniques are applicable if the observed image is a linear function of the original image and if the degradation process is partially or fully known. More recent restoration techniques have generated considerable interest by using a priori knowledge of the original image or of the imaging system to place constraints on the solution of iterative restoration algorithms. The most widely known algorithm of this type is perhaps the Gerchberg-Papoulis (GP) algorithm [1, 2]. In this algorithm the a priori constraints operate on the magnitude of the image signal and on the magnitude of its Fourier transform. This method was used for the restoration of signals of known nite duration or extent and to extrapolate band-limited functions. Youla [3] and Youla and Webb [4] have shown that linear iterative algorithms of the GP type have a natural geometric interpretation and they generalized them by using an alternating orthogonal projection method in Hilbert space. They also developed general conditions for convergence and addressed the problem of noise-sensitivity of the restoration procedure. Brie y, they envisaged every known property of an original image, or signal, F as restricting it to lie in a well-de ned closed convex set in a Hilbert space H. Thus m such properties place F in the intersection C0 of the corresponding closed convex set Ci , i = 1; 2; :::;m, and the problem of restoring F from its m properties is included in that of nding at least one point belonging to C0 . To this end they derived rules, or projection operators Pi, to assign to every F 2 H its nearest neighbor Pi(F) in Ci and demonstrated that, for every F 2 H, the sequence obtained by standard recursion: Fn+1 = T(Fn ) (1) Where T = Tm Tm?1 :::T1, n = 0 to 1 and F0 = F and: Ti = 1 + i (Pi ? 1) (2) and i 's are relaxation constants in the interval 0 < i < 2 converges weakly to a point of C0. A signi cant advantage of the theory of image restoration by the method of projections onto convex sets (POCS) over GP and related algorithms is that it enables a large number of a priori constraints to be incorporated in the restoration algorithms. Youla and Webb [4] provided a set of eleven constrains of physical signi cance together with their associated projection operators. Examples of convex sets corresponding to these projection operators are: the set of all functions that vanish outside a prescribed region, the set of all functions in H whose Fourier transforms assume a prescribed value over a closed region in the Fourier plane (the projection over these two sets correspond to the GP 3
algorithm), the set of all non negative functions in H that satisfy a given energy constraint , the subset of all functions in H that are non negative, the set of all function in H whose amplitude must lie in a prescribed closed interval, etc... The method of POCS was used to restore scalar elds (images) from incomplete information [5, 6, 7] but it can also be used to restore vector elds. In this case, new useful projection operators related to vector eld properties can be developed. In a recent paper published by the authors [8], the theory and numerical implementation of a projection operator to restore divergence free vector elds from incomplete information was presented. In this paper, projection operators related to the divergence, rotational, and hyperbolic components will rst be presented together whith a linear time numerical implementation. Projection operators related to the restoration of the optical ow will then be considered and compared to the Horn and Schunck method of determining optical ow [9]. Finally, a general method to develop new projection operators will be presented.
2 Useful projection operators Let G = (gx ; gy ) be the discrete vector eld to be restored by alternating projections from the initial partial information F = (fx ; fy ) de ned on each point (i; j) of an N = n n pixel plane. Adopting the following spatial distribution of the points: (fx (i; j + 1); fy (i; j + 1)) (fx (i + 1; j + 1); fy (i + 1; j + 1)) (3) (fx (i; j); fy (i; j)) (fx (i + 1; j); fy (i + 1; j)) the discrete partial derivatives with respect to x and y can be de ned by: dfx (i; j) dfdxx (i; j) dy dfy (i; j) dfdxy (i; j) dy
= = = =
(fx (i + 1; j + 1) ? fx (i; j + 1) + (fx (i + 1; j) ? fx (i; j))=2 (fx (i + 1; j + 1) ? fx (i + 1; j) + (fx (i; j + 1) ? fx (i; j))=2 (fy (i + 1; j + 1) ? fy (i; j + 1) + (fy (i + 1; j) ? fy (i; j))=2 (4) (fy (i + 1; j + 1) ? fy (i + 1; j) + (fy (i; j + 1) ? fy (i; j))=2
Since any 2-D vector eld F can be locally approximated by a linear vector eld ~ a rotational eld R~ and F~ which is a linear combination of a divergence eld D, ~ ~ two hyperbolic elds H1 and H2, that is [8]: F~ = aD~ + bR~ + cH~ 1 + dH~2 (5) the development of projection operators for the restoration of vector elds of known divergence, rotational or hyperbolic components suggests itself and will be presented rst. 4
2.1 Projection operators related to the divergence, the rotational, and the hyperbolic components
These projection operators, denoted respectively P1 , P2 , P3 , P4 , nd interesting applications in dynamic scene analysis for restoring optical ows. For example, the divergence free projection operator P1 can be used to restore motion elds due to translational, rotational and shearing motions or any combination of these three motions. They are not restricted to optical ows but can also be used on any kind of elds encountered in electromagnetism, uid dynamics, continuum mechanics etc. where divergence free, irrotational and shearing free elds are frequent. For example, projection operator P1 can be used for restoring in 3-D the deformation eld of an incompressible body from its 2-D projections on the image plane, the projection operators P3 and P4 can be used for restoring the velocity elds of undeformable bodies, and the projection operator P2 can be used for restoring zero rotational elds or elds for which the curl is know on some points. The four corresponding sets, being vector spaces, are convex and can be de ned as: C1 = fF j div F = Dg C2 = fF j rot F = Rg C3 = fF j hyp1 F = H1g C4 = fF j hyp2 F = H2g
div F(i; j) = + dfdxx (i; j) + dfdyy (i; j) rot F(i; j) = ? dfdyx (i; j) + dfdxy (i; j) hyp1 F(i; j) = + dfdxx (i; j) ? dfdyy (i; j) hyp2 F(i; j) = + dfdyx (i; j) + dfdxy (i; j)
(6)
Where D; R; H1; H2 are scalar function of i; j and the equations for div, rot, hyp1 and hyp2 are discretizations of their corresponding continuous equations. The projection of an arbitrary F onto each of these sets is realized by:
P : The projection operator to restore a vector eld of known divergence D(i; j) 1
is de ned by the equation
P1(F) = G = (gx ; gy )
(7)
where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) + w(i ? 1; j) ? w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by divF(i; j) ? D(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (8) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n 5
P : The projection operator to restore a vector eld of known curl R(i; j) is 2
de ned by the equation
P2(F) = G = (gx ; gy )
(9)
where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) + w(i ? 1; j ? 1) ? w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) ? w(i ? 1; j) + w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by rotF(i; j) ? R(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (10) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n
P3: The projection operator to restore a vector eld of known hyperbolic component H1(i; j) is de ned by the equation P3(F) = G = (gx ; gy ) (11) where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (?w(i; j) ? w(i ? 1; j) + w(i; j ? 1) + w(i ? 1; j ? 1))=2 and w is a scalar function de ned by hyp1 F(i; j) ? H1(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (12) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n P4: The projection operator to restore a vector eld of known hyperbolic component H2(i; j) is de ned by the equation P4(F) = G = (gx ; gy ) (13) where: gx (i; j) = fx (i; j) + (+w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) ? w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) ? w(i ? 1; j) + w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by hyp2 F(i; j) ? H2(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (14) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n The demonstration for P1 is given in appendix A. The demonstrations for P2,P3 , and P4 are very similar. Note that Eqs. 8, 10,12 and 14 require solving a system of (n ? 1) (n ? 1) equations. It will be shown in section 3 how this problem can be circumvented. 6
2.2 Projection operators related to the restoration of the optical ow
Vector elds are used in dynamic image analysis when the optical ow (or the distribution of apparent velocities of movement of brightness pattern in an image) is computed from successive views of a scene [9, 10]. The Horn and Schunck method of determining the optical ow [9] is among the most popular. This method will rst be brie y presented. It will then be shown how this method (or similar methods) can be formulated in terms of POCS and what are the advantages of doing so. Following Horn and Schunck, let the motion image be represented by a continuous variation of image brightness as a function of position and time: E(x; y; t). By assuming that the brightness of a particular pont in a pattern is constant when the pattern moves, they set: Ex fx + Ey fy + Et = 0 (15) where Ex , Ey and Et denotes the partial derivatives of E(x; y; t) with respect to x; y and t respectively, and fx , fy are the optical ow components. If Eq. 15 is rewritten as rE:F = ?Et (16) where F = (fx ; fy ) is the optical ow velocity, and rE is the spatial gradient of the image, one readily sees that the component of movement in the direction of the brightness gradient can be computed by using information about the spatial and temporal brightness gradient of the images. Equation 15 constrains the velocity but does not determine the component of movement in the direction of the iso-brightness contours, at right angles to the brightness gradient. Therefore, the velocity eld F = (fx ; fy ) cannot be computed locally without the introduction of an additional constraint. Since, in the absence of occluding edges, it can be assumed that the velocity eld of the image varies smoothly almost everywhere, this additional constraint can be introduced in the form of a smoothness constraint by minimizing, for example, the square of the magnitude of the gradient of the optical ow velocity, that is by minimizing Z Z
(@fx =@x)2 + (@fy =@x)2 + (@fx =@y)2 + (@fy =@y)2 dx dy
(17)
The problem is then to minimize the sum of the errors in the equation for the rate of change of image brightness (Eq. 15) and the measure of the departure from smoothness in the optical ow velocity (Eq. 17). The resulting variational problem consists in nding the optical velocity eld F = (fx ; fy ) which minimizesZ the functional Z a(Ex fx + Ey fy + Et)+ (@fx =@x)2 + (@fy =@x)2 + (@fx =@y)2 + (@fy =@y)2 dx dy (18) 7
where a is a weighting factor. The two corresponding Euler equations r2fx = a(Ex fx + Ey fy + Et)Ex (19) r2fy = a(Ex fx + Ey fy + Et)Ey where r2 = @ 2 =@x2 + @ 2 =@y2 are then obtained using the calculus of variation and are solved for each pixel by the Gauss-Seidel method of iterative solution of a nite di erence approximation. n n fxn+1 = f nx ? Ex1+f xa+(EEx2y+f yE+y2 )Et Ex (20) n n fyn+1 = f ny ? Ex1+f xa+(EEx2y+f Ey +y2 )Et Ey
where f x and f y are local averages as de ned in [9]. It will now be shown how the problem of nding the optical ow velocity from Eq. 15 and 17 can be solved by the method of POCS. A projection operator that restores a velocity eld from its projection on the brightness gradient of an image can be developed in the following way: Consider the set C5 of elds having a known inner product KN (i; j) with a given vector eld N(nx ; ny ) C5 = fF = (fx ; fy ) j nx :fx + ny :fy + KN = 0g
(21)
C5, being de ned at each point (i; j) by the equation of a line, is convex. Note that nx, ny , and KN in Eq. 21 respectively correspond to Ex , Ey , and Et in Eq. 15. It can be shown (appendix B) that the projection operator P5 of F onto C5 is realized by:
P5: The projection operator to restore a vector eld of which the inner product KN (i; j) with another vector eld N(nx ; ny ) is known is de ned by the equation P5(F) = G = (gx; gy )
(22)
where:
gx = fx ? nx (fx :nnx +2x +fyn:n2y y +KN ) gy = fy ? ny (fx :nnx +2x +fyn:n2y y +KN ) The constraint of smoothness can be formulated in terms of POCS in the following way: Let C6 be the set of vector elds de ned by: (
x C6 = F = (fx ; fy ) j df dx
2
y + df dx
8
2
+ dfdyx
2
+ dfdyy
2
KS 2
)
(23)
Where KS is a scalar function of i; j. It can be demonstrated that C6 is convex (appendix C). C6 is the intersection of C6ij where
2 2 C6ij = F = (fx ; fy ) j dfdxx (i; j) + dfdxy (i; j) +
dfx (i; j) dy
2
+
dfy (i; j) dy
2
KS (i; j) 2
(24)
Adopting the following pixel labeling: (fx12 ; fy12) (fx22 ; fy22) (25) (fx11 ; fy11) (fx21 ; fy21) and letting: R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (26) it can be shown (appendix D) that the local projection operator P6ij of F onto C6ij is realized by:
P6ij: The projection operator to restore a vector eld subjected to the smoothness constraint of Eq. 24 is de ned by P6ij (F) = G = (gx ; gy ) (27) where G = F if R2 KS2 and gx22 = (fx22 (KS2 =R2 + 1) + fx11(1 ? KS2 =R2))=2 gx11 = (fx11 (KS2 =R2 + 1) + fx22(1 ? KS2 =R2))=2 gx21 = (fx21 (KS2 =R2 + 1) + fx12(1 ? KS2 =R2))=2 gx12 = (fx12 (KS2 =R2 + 1) + fx21(1 ? KS2 =R2))=2 (28) gy22 = (fy22(KS2 =R2 + 1) + fy11 (1 ? KS2 =R2))=2 2 2 2 2 gy11 = (fy11(KS =R + 1) + fy22 (1 ? KS =R ))=2 gy21 = (fy21(KS2 =R2 + 1) + fy12 (1 ? KS2 =R2))=2 gy12 = (fy12(KS2 =R2 + 1) + fy21 (1 ? KS2 =R2))=2 otherwise. It will be shown in section 3 that this method of local projections corresponds to projecting onto C6. The optical ow velocity eld between two images is obtained by rst computing the component of movement in the direction of the brightness gradient of the image (Eq. 15) and then by iteratively alternating projections of this component onto C5 and C6 by the corresponding projection operators P5 and P6 instead of by solving Eq. 19 and 20 iteratively for each pixel. The computation time of the two methods is about the same. the advantages of POCS will become apparent after the following development. 9
2.3 Theoretical comparison of POCS with the Horn and Schunck method
This comparison can be achieved by rst demonstrating that the Horn and Schunck algorithm (Eq. 15 to 20) corresponds to an unique projection of the zero vector eld onto C5 by using a non-Euclidian norm. Consider a Hilbert space H with elements F; G; H, etc : : :and an inner product between a pair of vector F(fx ; fy ) and G(gx ; gy ) de ned by Z Z @fx @gx + @fx @gx + @fy @gy + @fy @gy dx dy (29) F:G = @x @x @y @y @x @x @y @y and the corresponding norm by Z Z
2 2 2 2 @f @f @f @f x x y y kF k = @x + @y + @x + @y dx dy (30) Now, the rule which assign to every G 2 H its nearest neighbor P5 (F) in C5 de nes the non-linear projection operator P5 : H ?! C5 unambiguously by means of the minimality criterion 2
kF ? P5(F)k = Gmin kF ? Gk 2C5
(31)
The problem is to nd G 2 C5 which satis es this equation. Using the norm 30, Eq. 31 can be written as: Z Z
@(fx ? gx ) 2 + @(fx ? gx) 2 + @x @y 5 5 2 @(fy ? gy ) + @(fy ? gy ) 2 dx dy (32) @x @y Letting F be the null vector eld, then Eq. 32 becomes: min kF ? Gk2 = Gmin G2C 2C
Z Z
@gx 2 + @gx 2 + @gy 2 + @gy 2 dx dy (33) min k G k = min G2C5 G2C5 @x @y @x @y Note that this equation is exactly the smoothness constraint introduced by Horn and Schunck (Eq. 17). In other words, the Horn and Schunck approach corresponds to projecting a null vector eld onto C5 by using the norm of Eq. 31. The problem with this approach is that it only nds a function G in C5 which minimizes Eq. 17 the smoothness constraint. In most cases this constraint is too severe because the only eld which truly minimizes it is a translational eld. The restoration of the optical ow by the method of POCS using P5 and P6 presents the advantage that the smoothness can be restricted to be equal to or less than 2
10
KS (i; j) instead of being minimized. The value of KS (i; j) can even be \tuned" until the algorithm converges. This can be done because, if the intersection of C5 and C6 is an empty set, the alternate iterative projections onto C5 and C6 will oscillate between two values. The amplitude of this oscillation is equal to the minimum distance between the two sets. This distance can be successively reduced by modifying the KS (i; j) and, as soon as the intersection of C5 and C6 is non-empty, the algorithm will converge.
2.4 A general method to develop projection operators
Many other projection operators for special applications can be developed by the method of Lagrange multipliers, by rst demonstrating that the constraint condition de nes a convex set and then by nding the shortest distance between F and G subject to this constraint (see also appendices A, B and D). To illustrate this method consider the restoration of a vector eld whose ow is known along one axis. The corresponding convex set is de ned by: C7 = fF = (fx ; fy ) j
n X i=1
fx (i; j) = K(j)g
(34)
which, since summation is a linear operator, is convex. We want to nd the shortest distance: X (fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 (35) i;j
subject to the constraint condition: n X i=1
gx(i; j) = K(j)
(36)
According to the method of Lagrange multipliers, we consider: n X n ? X i=1 j =1
(fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 +
n X j =1
w(j)
n X i=1
gx(i; j) ? K(j)
!
(37)
Setting the partial derivative with respect to gx and gy equal to zero, we nd: ?2(fx (i; j) ? gx (i; j)) + w(j) = 0 (38) ?2(fy (i; j) ? gy (i; j)) = 0 Solving these equations we nd that: w(j) = 2 (
n X i=1
fx (i; j) ? K(j))
11
(39)
and the projection operator of F onto C7 is realized by: P7(F) = G = (gx; gy ) where:
(40)
P gx (i; j) = fx (i; j) ? ( ni=1 fx (i; j) ? K(j)) gy (i; j) = fy (i; j)
3 A linear time implementation of global projection operators Equations (7,9,11,13) of the rst four projection operators, once written in matrix form, can be solved by the method of Hockney [11]. For an N = n n pixel image, a system of (n ? 1) (n ? 1) equations has to be solved at each iteration. For a large image, this requires a large memory (2N) and long computation time (O(KN log2 N); K being the number of iteration). The smoothness constraint as stated in Eq. 17 also requires the solution of a system of (n ? 1) (n ? 1) equations. Fortunately, this problem can be circumvented. A constraint being de ned on all the domain of a vector eld F, is also valid for each subregions Ri of (F). The constraint restricted to a region Ri de nes a convex set Ci. If the union of the regions Ri is equal to (F), then the intersection of all Ci is equal to C, that is \ [ (41) Ri = (F) =) Ci = C i
i
Alternating projections of F onto each Ci by the corresponding projection operators Pi converge to a eld which belongs to C. If the Ri are m m pixels regions then the projection operators Pi requires the solution of a system of (m ? 1)2 equations. Solving (n ? 1)2 =(m ? 1)2 systems of (m ? 1)2 equations is faster and easier than solving one system of (n ? 1)2 equations. If m = 2, (m ? 1)2 independent equations have to be solved. For m = 2, the local projection operator for the divergence, the rotational and the hyperbolic components are realized by computing for each points: gx11 = fx11 + (divF ? D)=4 gx12 = fx12 + (divF ? D)=4 gx21 = fx21 ? (divF ? D)=4 gx22 = fx22 ? (divF ? D)=4 (42) gy11 = fy11 + (divF ? D)=4 gy12 = fy12 ? (divF ? D)=4 gy21 = fy21 + (divF ? D)=4 gy22 = fy22 ? (divF ? D)=4
12
gx11 = fx11 ? (rotF ? R)=4 gx12 = fx12 + (rotF ? R)=4 gx21 = fx21 ? (rotF ? R)=4 gx22 = fx22 + (rotF ? R)=4 (43) gy11 = fy11 + (rotF ? R)=4 gy12 = fy12 + (rotF ? R)=4 gy21 = fy21 ? (rotF ? R)=4 gy22 = fy22 ? (rotF ? R)=4 gx11 = fx11 + (hyp1 F ? H1)=4 gx12 = fx12 + (hyp1 F ? H1)=4 gx21 = fx21 ? (hyp1 F ? H1)=4 gx22 = fx22 ? (hyp1 F ? H1)=4 (44) gy11 = fy11 ? (hyp1 F ? H1)=4 gy12 = fy12 + (hyp1 F ? H1)=4 gy21 = fy21 ? (hyp1 F ? H1)=4 gy22 = fy22 + (hyp1 F ? H1)=4 gx11 = fx11 + (hyp2 F ? H2)=4 gx12 = fx12 ? (hyp2 F ? H2)=4 gx21 = fx21 + (hyp2 F ? H2)=4 gx22 = fx22 ? (hyp2 F ? H2)=4 (45) gy11 = fy11 + (hyp2 F ? H2)=4 gy12 = fy12 + (hyp2 F ? H2)=4 gy21 = fy21 ? (hyp2 F ? H2)=4 gy22 = fy22 ? (hyp2 F ? H2)=4 These equations were obtained by setting n = 2 in Eqs. 8, 10,12 and 14. The performance of this local projection method will be illustrated by an example: a linear hyperbolic eld de ned on 16 16 pixels has to be restored from its normal components only by using the a priori knowledge that its divergence is zero. Thus, two convex sets can be used: C1, the set of divergence free vector elds, and C5, the set of vector elds with known normal components. The hyperbolic eld was restored in three ways. The rst restoration was done by projecting alternatively onto C1 and C5. Projecting onto C1 , amounts to solving a system of 16 16 equations (this was done by the method of Hockney). The second restoration was done by considering C1, as the intersection of 9 convex sets of vectors elds that are divergence free on 6 6 pixel regions. This amounts to solving 9 systems of 5 5 equations. The third restoration was done by considering C1 as the intersection of 15 15 convex sets that are divergence free on 2 2 pixel regions. The time per iteration of this latter method is less than for the other two methods. All the restorations were done with a relaxation parameter equal to 1.9. The three methods are compared by computing the
13
percent error ek at the kth iteration [5, 8]: Fk k (46) ek = 100 k Hk ? Hk here Fk is the vector eld restored at iteration k and H is the original hyperbolic eld. From gure 1 it can be seen that the percent error curves of the three methods are quite comparable. The second restoration curve is between the two others and is not shown for clarity. The memory required by the local projection method (2 pixels) is small and the computation time is only O(Kl N) instead of O(Kg N log2 N) as for the method of Hockney, where Kl and Kg are the number of iterations of each method. As can be seen from Fig 1, Kl is approximately equal to Kg . Therefore the time complexity has been reduced from O(KN log2 N) to O(KN), where K is the number of iterations. The reason why the local projection method is so ecient is that all 2 2 pixel regions that do not have pixels in common are independent, or orthogonal, convex sets, and that projecting onto orthogonal convex sets amounts to projecting directly onto their intersection (the global constraint). These results show the advantages of the local projection method which is moreover much easier to implement than the Hockney method.
4 Examples of application This section presents two examples of vector eld restoration by the method of POCS. Other simple projection operators will also be presented to show how they can be easily derived to face practical applications. EXAMPLE I: The restoration of velocity eld (optical ow). The method of POCS is not restricted to optical ow due to motion of rigid bodies but can also be used to restore any velocity elds describing, for example, the interframe motion of uids or deformable bodies. The method can thus be illustrated on a divergence free vector eld (Fig. 2) which may represent the motion of a uid. This eld has to be restored from the partial information shown in Figure 3 which represents the components of movement in the direction of the brightness gradient. The latter was generated randomly at each point and then smoothed by a 15 15 pixel mean lter. Figure 3 thus represents the movement in the direction of a correlated brightness gradient of a random image submitted to the motion described by the velocity eld of Figure 2. The vector eld of Figure 2 can be restored from Figure 3 by projection onto three convex sets: C1, the set of divergence free vector elds; C5 , the set of vector elds whose velocity component in the direction of the brightness gradient corresponds at each point to that shown in Figure 3; and C6, the 14
set of vector elds respecting the smoothness constraint. The corresponding projection operators were de ned in section 2 as P1 , P5 and P6 respectively. When the latter projection operator is used, the optimal values of Ks (see equation 28) has to be determined as explained below. In order to measure the contribution of each projection operator, restorations have been obtained: a) by projecting onto C5 and C6 ; b) by projecting onto C1 and C5; and c) by projecting onto C1, C5 and C6. (Restoring by projection onto C1 and C6 only is, of course, meaningless). These results are compared to those obtained with the Horn and Schunck method (Figs. 4, 5, 6) and are summarized in Figure 7 which presents the curves of percent errors as a function of the number of iterations for each restoration.
4.1 Restoration by projecting onto C5 and C6
The value of Ks (x; y) (see equation 28) is generally not known for real images. So when the P6 projection operator is used its optimal value has to be determined rst. If Ks (x; y) is set to a constant equal to the maximal expected value, the restored eld can not be smooth enough. Conversely, if Ks is set too low, the restored eld can be too smooth. Fortunately, this parameter can be set to the value which will restore the smoothest eld while respecting the other constraints (C5 in this case). Geometrically, this can be seen by realizing that the size of C6 increases with Ks . If Ks is too small, the intersection with C5 will be empty, and the POCS algorithm will oscillate between the two convex sets instead of converging to an intersection point. If Ks is too large, the intersection with C5 will also be too large, the point of convergence can be too far from the desired solution and the restored eld will not be smooth enough. The optimal value of Ks can be determined by measuring, after some iterations, the RMS error between the elds restored before and after projecting onto C6. A periodic large error means that the algorithm oscillates between C6 and C5 instead of converging to a point of intersection. This means that there is no eld satisfying the two constraints. The amplitude of the RMS error is proportional to the distance between the two sets and indicates by how much Ks has to be increased. If the two convex sets intersect, the RMS error will continuously decrease with the number of iterations instead of being periodic. The smallest Ks value for which C5 and C6 intersect will restore the smoothest eld respecting the constraint of C5. The optimal Ks value can thus be determined by starting with a small Ks and by progressively increasing it until the two convex sets intersect. For restoring the eld of Figure 2 by projecting onto C5 and C6 an optimal value of Ks = 0:4 was determined in this way. The restored eld is presented in Figure 4. The percent error curves (not to be confused with the RMS error) between the original and restored elds as a function of the number of iterations with Ks = 0:4 and Ks = 1:2 are presented in Figure 7.
15
4.2 Restoration by projecting onto C1 and C5
The results are improved by projecting onto C1 and C5 instead of C5 and C6 (Fig. 7), meaning that, in this particular example, the divergence free constraint (C1) is more useful than the smoothness constraint (C6 ).
4.3 Restoration by projecting onto C1, C5 and C6
The best results are, of course, obtained by projecting onto C1 , C5 and C6 (Fig. 5 and Fig. 7), since the eld to be restored is both divergence free and smooth. In this case the optimal Ks value, determined as before, was found to be 1:2. Restoration obtained by using the same projection operators but with Ks = 0:4, yields results comparable to those obtained by projecting onto C5 and C6 with Ks = 0:4, which means that the intersection of C6 (with Ks = 0:4) with C1 and C5 is empty, or that the elds satisfying both the constraints of C1 and C5 are less smooth than the elds satisfying C5 only. It is interesting to note that the maximum value of Ks in the original eld is 1:112).
4.4 Comparison with the Horn and Schunck method
The results obtained with the Horn ans Schunck method are shown in Figures 6 and 7. They are comparable to those obtained by projecting onto C5 and C6 with Ks = 0:4 (Figs. 4 and 7). This is not surprising because both methods use the same constraints, although they are not implemented in the same way. These two methods are comparable in terms of memory and computation time and both can easily be multi-grid implemented. The only di erence of POCS, in this case, is that it constraints the smoothness to be smaller or equal to Ks (equation 28), (whereas the Horn and Schunck method minimizes this constraint), and that it allows, by adjusting the Ks parameter, to nd the smoothest eld which satis es C5 . Projecting onto C1 instead of C6 improved slightly the results (Fig. 7) meaning that the divergence free constraint is a very useful constraint that can even replace the smoothness constraint in this example. Using this constraint together with C5 and C6 (Figs. 5 and 7) gives the best results because all the a priori knowledge about the original eld (Fig. 2) are then used. EXAMPLE II: The restoration of a plane ow The eld of Figure 8 is an example of a plane ow with stagnation point (Hiemenz ow) [12]. The ow arrives from the y axis, impinges on a at wall placed at y = 0, divides into two streams on the wall and leaves in both directions leading to a stagnation point at x = y = 0. The eld of Figure 8 can be restored by the method of POCS by using only the following informations: the divergence of the eld is zero, the ow (the line integral of the divergence) is known to be zero on the external boundaries of 16
the eld, the fy components of the eld are known only on the upper boundary, the eld is known to be zero at y = 0, the eld is known to be symmetric about the y axis, and the eld is known to vary smoothly (smoothness constraint). The projection operators corresponding to the divergence free and smoothness constraints are those described in section 2 as P1 and P6 . The other four projection operators are very easy to implement: after each projection onto C1 and C6 (Ks = 0:03) the ow on the external boundaries is set to zero by subtracting from each point on the boundaries the mean vector of the points on the boundaries, the points at y = 0 are set to zero, the fy components of the points on the upper boundaries are replaced by the known one, and the eld is made symmetric by averaging the vectors of each pair of points symmetric about the y axis. The eld restored after 100 iterations is shown in Figure 9 and the percent error per iteration in Figure 10. Restorations were also made without the smoothness, symmetry and ow constraints to show the relative importance of these factors. From Figure 10, it can be seen that the symmetry constraint does not improve the restoration very much, but that the smoothness constraint is very ecient. The ow constraint, which is theoretically redundant because it is implied by the divergence free constraint, has the interesting property to reduce the percent error in the rst iterations.
5 Conclusion The main contribution of this paper was to show that the method of POCS, as originally developed by Youla and Webb [3, 4] for scalar elds (or images) can be applied to the restoration of vector elds. Projection operators were developed to restore elds of known divergence, rotational, or hyperbolic components (P1 to P4). It was shown how these projection operators can eciently be implemented by the method of local projections (section 3). With this method, computations which rst appear to be prohibitively long are reduced to a very manageable size. This method can be used whenever a convex set can itself be represented as the intersection of many local convex sets and will be very useful in developing other projection operators. Projection operators to restore vector elds from their inner product with a known eld, or to restore elds subjected to a smoothness constraint were also presented (P5 and P6 ). These two last projection operators correspond to the a priori knowledge used by Horn and Schunck to determine the optical ow. Finally, a general method to derive new projection operators was proposed. When applied to machine vision problems, such as the determination of the optical ow, which are generally solved by the calculus of variations, the method of POCS presents some interesting aspects: 1) the constraints are introduced one by one and are implemented by nding the corresponding projection operator. This is easier than nding the Euler equations of a complex functional 17
which includes all the constraints; 2) if a constraint must be smaller or equal than a value K (instead of being only equal to K) it is possible to nd the optimal K which also satis es the other constraints by checking, after some iterations, if the algorithm is not oscillating between two convex sets instead of converging to a point which belongs to their intersection (see section 4, example 1); 3) since the algorithm converges to a point of intersection of all the convex sets, the restored eld really has all the properties described by each of these convex sets; 4) some a priori knowledge, for example, that the Fourier transforms of the vector components are band-limited, or that the vector components assume prescribed values over some pixels of the images, are easier to formulate in terms of POCS than with variational approaches. Of course, the more a priori knowledge introduced in the algorithm as convex sets, the more precise the restoration. Examples of applications of this method on real noisy images (echographic textures and two-dimensional echocardiograms) are presented in reference [13].
18
Figure 1: Hyperbolic eld restoration. Error versus iteration number by solving one system of 15 15 equations, and 15 15 systems of one equation. Figure 2: A simulated divergence free vector eld Figure 3: Projection of the vector eld shown in Figure 2 along a correlated brightness gradient Figure 4: Vector eld restored by projecting onto C5 and C6(Ks = 0:4) Figure 5: Vector eld restored by projecting onto C1 , C5 , and C6 Figure 6: Vector eld restored by the Horn and Schunck method Figure 7: Error versus iteration number by projecting onto C5 and C6 (Ks = 1:2), onto C5 and C6(Ks = 0:4), onto C1 and C5, onto C1 and C5 and C6 (Ks = 1:2), and by Horn and Schunck method Figure 8: A plane ow with stagnation point Figure 9: Vector eld restored after 100 iterations Figure 10: Error versus iteration number by restoring with the six constraints, with ve constraints (without the symmetry constraint or without the ow constraint), and with four constraints (without the symmetry and smoothness constraints) 19
A Demonstration for P1 A vector eld has to be restored from the partial information F = (fx ; fy ) and the a priori knowledge that its divergence must be D(i; j). The problem is thus to nd the shortest distance between F and C1; that is, the minimum of: n X n X i=1 j =1
(fx (i; j) ? gx(i; j))2 + (fy (i; j) ? gy (i; j))2
(47)
subject to the constraint: divG(i; j) = D(i; j) According to the method of Lagrange multipliers we consider: n X n X i=1 j =1
(48)
(fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 +
n X n X i=1 j =1
w(i; j) (divG(i; j) ? D(i; j))
(49)
By setting to zero the partial derivatives with respect to gx and gy we obtain: 2(fx (i; j) ? gx (i; j)) = +w(i ? 1; j) + w(i ? 1; j ? 1) ? w(i; j) ? w(i; j ? 1) 2(fy (i; j) ? gy (i; j)) = ?w(i; j) ? w(i ? 1; j) + w(i; j ? 1) + w(i ? 1; j ? 1) (50) where w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::;n and j = 1; 2; :::;n. Adding these equations we nd: divF(i; j) ? D(i; j) = 4 w(i; j) ?(w(i ? 1; j ? 1) + w(i ? 1; j + 1) + w(i + 1; j ? 1) + w(i + 1; j + 1) (51)
B Demonstration for P5 We want to nd, on each point, the minimum of (fx ? gx)2 + (fy ? gy )2
(52)
subject to the constraint gx :nx + gy :ny + KN = 0 According to the method of Lagrange multipliers we consider (fx ? gx)2 + (fy ? gy )2 + w(gx :nx + gy :ny + KN ) 20
(53) (54)
By setting to zero the partial derivatives with respect to gx and gy we obtain: 2(?fx + gx ) + w:nx = 0 (55) 2(?fy + gy ) + w:ny = 0 Solving these equations we nd that: +fx :nx + fy :ny + KN = w(n2x + n2y )=2 (56) the solution is thus: gx = fx ? nx(fx :nx + fy :ny + KN )=(n2x + n2y ) (57) gy = fy ? ny (fx :nx + fy :ny + KN )=(n2x + n2y )
C Convexity of C6 Let g1 and g2 belong to C6 . we have to prove that 8 2 [0; 1]; g1 + (1 ? )g2 2 C6 g1, g2 2 C6 implies that the two points: ( @g@x1x ; @g@x1y ; @g@y1x ; @g@y1y ) and ( @g@x2x ; @g@x2y ; @g@y2x ; @g@y2y ) belong to the hypersphere:
(58)
x2 + y2 + z 2 + t2 KS2 : (59) Since an hypersphere is convex, the point with the following coordinates ( @f@x1x + (1 ? ) @f@x2x ; @f@y1x + (1 ? ) @f@y2x ; @f@x1y + (1 ? ) @f@x2y ; @f@y1y + (1 ? ) @f@y2y ) (60) also belongs to the hypersphere and we can write: 2 2 @f @f @f @f 1x 2x 2x 1x + + (1 ? ) + (1 ? ) @x @x @y @y 2 2 @f @f @f @f 1y 2y 2y 1y + @x + (1 ? ) @x + @y + (1 ? ) @y KS2 so C6 is convex.
(61)
D Demonstration for P6ij Using the same pixel labeling as in Eq. 25, we want to nd the shortest distance: (fx11 ? gx11)2 + (fx12 ? gx12)2 + (fx21 ? gx21)2 + (fx22 ? gx22)2 +(fy11 ? gy11)2 + (fy12 ? gy12)2 + (fy21 ? gy21)2 + (fy22 ? gy22)2 (62) 21
subject to the constraint condition: (gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 KS2 (63) This constraint can be written: (gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 = L2 (64) with 0 L KS . According to the method of Lagrange multipliers we consider: (fx11 ? gx11)2 + (fx12 ? gx12)2 + (fx21 ? gx21)2 + (fx22 ? gx22)2 +(fy11 ? gy11)2 + (fy12 ? gy12)2 + (fy21 ? gy21)2 + (fy22 ? gy22)2 +w((gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 ? L2 ) (65) By setting to zero the partial derivatives with respect to gx and gy we obtain: ?fx11 + gx11 + 2w(+gx11 ? gx22) = 0 ?fx22 + gx22 + 2w(+gx22 ? gx11) = 0 ?fx12 + gx12 + 2w(+gx12 ? gx21) = 0 ?fx21 + gx21 + 2w(+gx21 ? gx12) = 0 (66) ?fy11 + gy11 + 2w(+gy11 ? gy22) = 0 ?fy22 + gy22 + 2w(+gy22 ? gy11) = 0 ?fy12 + gy12 + 2w(+gy12 ? gy21) = 0 ?fy21 + gy21 + 2w(+gy21 ? gy12) = 0 By substitution of equations 66 into Eq. 62 we see that we have to minimize w2. From equations 66 we can write: fx22 ? fx12 + fx21 ? fx11 = (1 + 4w)(gx22 ? gx12 + gx21 ? gx11) fx22 ? fx21 + fx12 ? fx11 = (1 + 4w)(gx22 ? gx21 + gx12 ? gx11) (67) fy22 ? fy12 + fy21 ? fy11 = (1 + 4w)(gy22 ? gy12 + gy21 ? gy11) fy22 ? fy21 + fy12 ? fy11 = (1 + 4w)(gy22 ? gy21 + gy12 ? gy11) By adding these equations we obtain: R2 = (1 + 4w)2L2 (68) where R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (69) To solve Eq. 68 for 0 L KS and w2 minimum we have to consider two cases: If R2 KS2 then Eq. 68 is veri ed for w = 0, R2 = L2 and P6(F) = 22
G = (gx; gy ) = F. In the other case, w2 is minimized for L = KS and w = (R2=KS2 ? 1)=4. And the solution is: gx22 = (fx22 (KS2 =R2 + 1) + fx11(1 ? KS2 =R2))=2 gx11 = (fx11 (KS2 =R2 + 1) + fx22(1 ? KS2 =R2))=2 gx21 = (fx21 (KS2 =R2 + 1) + fx12(1 ? KS2 =R2))=2 gx12 = (fx12 (KS2 =R2 + 1) + fx21(1 ? KS2 =R2))=2 (70) gy22 = (fy22(KS2 =R2 + 1) + fy11 (1 ? KS2 =R2))=2 gy11 = (fy11(KS2 =R2 + 1) + fy22 (1 ? KS2 =R2))=2 gy21 = (fy21(KS2 =R2 + 1) + fy12 (1 ? KS2 =R2))=2 gy12 = (fy12(KS2 =R2 + 1) + fy21 (1 ? KS2 =R2))=2 R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (71) Acknowledgment: this work was supported in part by the National Sciences and Engineering Research council of Canada and the Canadian Heart Foundation through a Research Scholarship to G. E. Mailloux, and in part by the NIH research grant no. NS22407.
References [1] R.W. Gerchberg, Super-resolution through error energy reduction, Opt. Acta, vol.21, pp. 709-720, 1974. [2] A. Papoulis, A new algorithm in spectral analysis and band-limited extrapolation, IEEE Trans. Circuits Syst., vol. CAS-22, pp. 735-742, 1975. [3] D.C. Youla, Generalized image restoration by the method of alternating projections, IEEE Trans. Circuits Syst., vol. CAS-25, pp. 694-702, 1978. [4] D.C. Youla and H. Webb, Image restoration by the method of convex projections: Part 1:-Theory, IEEE Trans. Med. Imaging, vol. MI-1, pp. 81-94, 1982. [5] M.I. Sezan and H. Stark, Image restoration by the method of convex projections: Part 2-Application and numerical results, IEEE Trans. Med. Imaging, vol. MI-1, pp. 99-101, 1982. [6] M.I. Sezan and H. Stark, Tomographic image reconstruction from incomplete view data by convex projection and direct Fourier inversion, IEEE Trans. Med. Imaging, vol. MI-3, pp. 91-98, 1984. [7] R.T. Chin, C-L Yeh and W.S. Olson, Restoration of multichannel microwave radiometric images, IEEE Trans. Pattern Anal. Machine Intell., vol. PAMI-7, pp.475-484, 1985. 23
[8] P.Y. Simard and G.E. Mailloux, A projection operator for the restoration of divergence free vector elds, IEEE Trans. Pattern Anal. Machine Intell. PAMI-10, pp. 248-256, 1988. [9] B.K.P. Horn and B.G. Schunck, Determining optical ow, Artif. Intell., vol. 17, pp. 185-203, 1983. [10] E.C. Hildreth, Computation underlying the measurement of visual motion, Artif. Intell., vol. 23, pp. 309-354, 1984. [11] R. Vichnevetsky. Computer Methods for Partial Di erential Equations. vol. 1. Englewoods Cli s, Prentice-Hall, 1981. [12] H. Schlichting. Boundary Layer Theory. New-York: mcGraw-Hill, 1960. [13] G. E. Mailloux, F. Langlois, P. Y. Simard and M. Bertrand, Restoration of the velocity eld of the heart from two-dimentional echocardiograms. IEEE Trans. Med. Imaging, vol. 8, 2, pp.143-153, 1989.
24
In this paper, the theory of image restoration by projections onto closed convex sets is applied to the restoration of vector elds. A set of useful projection operators is presented together with a linear time numerical implementation. These projection operators can be used to restore from partial information the velocity or deformation elds computed between successive views of a scene. They also nd applications in the restoration of vector elds of physical quantities as those encountered in mechanics, hydrodynamics or electromagnetism. The method is compared with the variational approach and illustrated by restoring simulated vector elds.
P. Y. Simard was with the Institut de G enie Biomedical Ecole Polytechnique, C. P. 6079, Succ A, Montreal, P.Q., Canada, H3C 3A7. He is now with the Department of Computer Science, University of Rochester, Rochester, NY 14627. y G. E. Mailloux is with the Institut de G enie Biomedical Ecole Polytechnique, C. P. 6079, Succ A, Montreal, P.Q., Canada, H3C 3A7.
1
List of symbols caligraphic H lambda in nity nabla integral partial norm in summation delta bigcup longarrow bigcap mu
H
1 r R @
k 2P S = T)
2
1 Introduction The aim of image restoration is to make as good an estimate as possible of an original image subjected to some form of degradation. Usually, images are restored by inverse ltering or constrained deconvolution techniques. These techniques are applicable if the observed image is a linear function of the original image and if the degradation process is partially or fully known. More recent restoration techniques have generated considerable interest by using a priori knowledge of the original image or of the imaging system to place constraints on the solution of iterative restoration algorithms. The most widely known algorithm of this type is perhaps the Gerchberg-Papoulis (GP) algorithm [1, 2]. In this algorithm the a priori constraints operate on the magnitude of the image signal and on the magnitude of its Fourier transform. This method was used for the restoration of signals of known nite duration or extent and to extrapolate band-limited functions. Youla [3] and Youla and Webb [4] have shown that linear iterative algorithms of the GP type have a natural geometric interpretation and they generalized them by using an alternating orthogonal projection method in Hilbert space. They also developed general conditions for convergence and addressed the problem of noise-sensitivity of the restoration procedure. Brie y, they envisaged every known property of an original image, or signal, F as restricting it to lie in a well-de ned closed convex set in a Hilbert space H. Thus m such properties place F in the intersection C0 of the corresponding closed convex set Ci , i = 1; 2; :::;m, and the problem of restoring F from its m properties is included in that of nding at least one point belonging to C0 . To this end they derived rules, or projection operators Pi, to assign to every F 2 H its nearest neighbor Pi(F) in Ci and demonstrated that, for every F 2 H, the sequence obtained by standard recursion: Fn+1 = T(Fn ) (1) Where T = Tm Tm?1 :::T1, n = 0 to 1 and F0 = F and: Ti = 1 + i (Pi ? 1) (2) and i 's are relaxation constants in the interval 0 < i < 2 converges weakly to a point of C0. A signi cant advantage of the theory of image restoration by the method of projections onto convex sets (POCS) over GP and related algorithms is that it enables a large number of a priori constraints to be incorporated in the restoration algorithms. Youla and Webb [4] provided a set of eleven constrains of physical signi cance together with their associated projection operators. Examples of convex sets corresponding to these projection operators are: the set of all functions that vanish outside a prescribed region, the set of all functions in H whose Fourier transforms assume a prescribed value over a closed region in the Fourier plane (the projection over these two sets correspond to the GP 3
algorithm), the set of all non negative functions in H that satisfy a given energy constraint , the subset of all functions in H that are non negative, the set of all function in H whose amplitude must lie in a prescribed closed interval, etc... The method of POCS was used to restore scalar elds (images) from incomplete information [5, 6, 7] but it can also be used to restore vector elds. In this case, new useful projection operators related to vector eld properties can be developed. In a recent paper published by the authors [8], the theory and numerical implementation of a projection operator to restore divergence free vector elds from incomplete information was presented. In this paper, projection operators related to the divergence, rotational, and hyperbolic components will rst be presented together whith a linear time numerical implementation. Projection operators related to the restoration of the optical ow will then be considered and compared to the Horn and Schunck method of determining optical ow [9]. Finally, a general method to develop new projection operators will be presented.
2 Useful projection operators Let G = (gx ; gy ) be the discrete vector eld to be restored by alternating projections from the initial partial information F = (fx ; fy ) de ned on each point (i; j) of an N = n n pixel plane. Adopting the following spatial distribution of the points: (fx (i; j + 1); fy (i; j + 1)) (fx (i + 1; j + 1); fy (i + 1; j + 1)) (3) (fx (i; j); fy (i; j)) (fx (i + 1; j); fy (i + 1; j)) the discrete partial derivatives with respect to x and y can be de ned by: dfx (i; j) dfdxx (i; j) dy dfy (i; j) dfdxy (i; j) dy
= = = =
(fx (i + 1; j + 1) ? fx (i; j + 1) + (fx (i + 1; j) ? fx (i; j))=2 (fx (i + 1; j + 1) ? fx (i + 1; j) + (fx (i; j + 1) ? fx (i; j))=2 (fy (i + 1; j + 1) ? fy (i; j + 1) + (fy (i + 1; j) ? fy (i; j))=2 (4) (fy (i + 1; j + 1) ? fy (i + 1; j) + (fy (i; j + 1) ? fy (i; j))=2
Since any 2-D vector eld F can be locally approximated by a linear vector eld ~ a rotational eld R~ and F~ which is a linear combination of a divergence eld D, ~ ~ two hyperbolic elds H1 and H2, that is [8]: F~ = aD~ + bR~ + cH~ 1 + dH~2 (5) the development of projection operators for the restoration of vector elds of known divergence, rotational or hyperbolic components suggests itself and will be presented rst. 4
2.1 Projection operators related to the divergence, the rotational, and the hyperbolic components
These projection operators, denoted respectively P1 , P2 , P3 , P4 , nd interesting applications in dynamic scene analysis for restoring optical ows. For example, the divergence free projection operator P1 can be used to restore motion elds due to translational, rotational and shearing motions or any combination of these three motions. They are not restricted to optical ows but can also be used on any kind of elds encountered in electromagnetism, uid dynamics, continuum mechanics etc. where divergence free, irrotational and shearing free elds are frequent. For example, projection operator P1 can be used for restoring in 3-D the deformation eld of an incompressible body from its 2-D projections on the image plane, the projection operators P3 and P4 can be used for restoring the velocity elds of undeformable bodies, and the projection operator P2 can be used for restoring zero rotational elds or elds for which the curl is know on some points. The four corresponding sets, being vector spaces, are convex and can be de ned as: C1 = fF j div F = Dg C2 = fF j rot F = Rg C3 = fF j hyp1 F = H1g C4 = fF j hyp2 F = H2g
div F(i; j) = + dfdxx (i; j) + dfdyy (i; j) rot F(i; j) = ? dfdyx (i; j) + dfdxy (i; j) hyp1 F(i; j) = + dfdxx (i; j) ? dfdyy (i; j) hyp2 F(i; j) = + dfdyx (i; j) + dfdxy (i; j)
(6)
Where D; R; H1; H2 are scalar function of i; j and the equations for div, rot, hyp1 and hyp2 are discretizations of their corresponding continuous equations. The projection of an arbitrary F onto each of these sets is realized by:
P : The projection operator to restore a vector eld of known divergence D(i; j) 1
is de ned by the equation
P1(F) = G = (gx ; gy )
(7)
where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) + w(i ? 1; j) ? w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by divF(i; j) ? D(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (8) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n 5
P : The projection operator to restore a vector eld of known curl R(i; j) is 2
de ned by the equation
P2(F) = G = (gx ; gy )
(9)
where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) + w(i ? 1; j ? 1) ? w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) ? w(i ? 1; j) + w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by rotF(i; j) ? R(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (10) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n
P3: The projection operator to restore a vector eld of known hyperbolic component H1(i; j) is de ned by the equation P3(F) = G = (gx ; gy ) (11) where: gx (i; j) = fx (i; j) + (?w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) + w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (?w(i; j) ? w(i ? 1; j) + w(i; j ? 1) + w(i ? 1; j ? 1))=2 and w is a scalar function de ned by hyp1 F(i; j) ? H1(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (12) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n P4: The projection operator to restore a vector eld of known hyperbolic component H2(i; j) is de ned by the equation P4(F) = G = (gx ; gy ) (13) where: gx (i; j) = fx (i; j) + (+w(i ? 1; j) ? w(i ? 1; j ? 1) + w(i; j) ? w(i; j ? 1))=2 gy (i; j) = fy (i; j) + (+w(i; j) ? w(i ? 1; j) + w(i; j ? 1) ? w(i ? 1; j ? 1))=2 and w is a scalar function de ned by hyp2 F(i; j) ? H2(i; j) = 4 w(i; j) ? w(i ? 1; j ? 1) ? w(i ? 1; j + 1) ? w(i + 1; j ? 1) ? w(i + 1; j + 1) (14) and w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::; n j = 1; 2; :::; n The demonstration for P1 is given in appendix A. The demonstrations for P2,P3 , and P4 are very similar. Note that Eqs. 8, 10,12 and 14 require solving a system of (n ? 1) (n ? 1) equations. It will be shown in section 3 how this problem can be circumvented. 6
2.2 Projection operators related to the restoration of the optical ow
Vector elds are used in dynamic image analysis when the optical ow (or the distribution of apparent velocities of movement of brightness pattern in an image) is computed from successive views of a scene [9, 10]. The Horn and Schunck method of determining the optical ow [9] is among the most popular. This method will rst be brie y presented. It will then be shown how this method (or similar methods) can be formulated in terms of POCS and what are the advantages of doing so. Following Horn and Schunck, let the motion image be represented by a continuous variation of image brightness as a function of position and time: E(x; y; t). By assuming that the brightness of a particular pont in a pattern is constant when the pattern moves, they set: Ex fx + Ey fy + Et = 0 (15) where Ex , Ey and Et denotes the partial derivatives of E(x; y; t) with respect to x; y and t respectively, and fx , fy are the optical ow components. If Eq. 15 is rewritten as rE:F = ?Et (16) where F = (fx ; fy ) is the optical ow velocity, and rE is the spatial gradient of the image, one readily sees that the component of movement in the direction of the brightness gradient can be computed by using information about the spatial and temporal brightness gradient of the images. Equation 15 constrains the velocity but does not determine the component of movement in the direction of the iso-brightness contours, at right angles to the brightness gradient. Therefore, the velocity eld F = (fx ; fy ) cannot be computed locally without the introduction of an additional constraint. Since, in the absence of occluding edges, it can be assumed that the velocity eld of the image varies smoothly almost everywhere, this additional constraint can be introduced in the form of a smoothness constraint by minimizing, for example, the square of the magnitude of the gradient of the optical ow velocity, that is by minimizing Z Z
(@fx =@x)2 + (@fy =@x)2 + (@fx =@y)2 + (@fy =@y)2 dx dy
(17)
The problem is then to minimize the sum of the errors in the equation for the rate of change of image brightness (Eq. 15) and the measure of the departure from smoothness in the optical ow velocity (Eq. 17). The resulting variational problem consists in nding the optical velocity eld F = (fx ; fy ) which minimizesZ the functional Z a(Ex fx + Ey fy + Et)+ (@fx =@x)2 + (@fy =@x)2 + (@fx =@y)2 + (@fy =@y)2 dx dy (18) 7
where a is a weighting factor. The two corresponding Euler equations r2fx = a(Ex fx + Ey fy + Et)Ex (19) r2fy = a(Ex fx + Ey fy + Et)Ey where r2 = @ 2 =@x2 + @ 2 =@y2 are then obtained using the calculus of variation and are solved for each pixel by the Gauss-Seidel method of iterative solution of a nite di erence approximation. n n fxn+1 = f nx ? Ex1+f xa+(EEx2y+f yE+y2 )Et Ex (20) n n fyn+1 = f ny ? Ex1+f xa+(EEx2y+f Ey +y2 )Et Ey
where f x and f y are local averages as de ned in [9]. It will now be shown how the problem of nding the optical ow velocity from Eq. 15 and 17 can be solved by the method of POCS. A projection operator that restores a velocity eld from its projection on the brightness gradient of an image can be developed in the following way: Consider the set C5 of elds having a known inner product KN (i; j) with a given vector eld N(nx ; ny ) C5 = fF = (fx ; fy ) j nx :fx + ny :fy + KN = 0g
(21)
C5, being de ned at each point (i; j) by the equation of a line, is convex. Note that nx, ny , and KN in Eq. 21 respectively correspond to Ex , Ey , and Et in Eq. 15. It can be shown (appendix B) that the projection operator P5 of F onto C5 is realized by:
P5: The projection operator to restore a vector eld of which the inner product KN (i; j) with another vector eld N(nx ; ny ) is known is de ned by the equation P5(F) = G = (gx; gy )
(22)
where:
gx = fx ? nx (fx :nnx +2x +fyn:n2y y +KN ) gy = fy ? ny (fx :nnx +2x +fyn:n2y y +KN ) The constraint of smoothness can be formulated in terms of POCS in the following way: Let C6 be the set of vector elds de ned by: (
x C6 = F = (fx ; fy ) j df dx
2
y + df dx
8
2
+ dfdyx
2
+ dfdyy
2
KS 2
)
(23)
Where KS is a scalar function of i; j. It can be demonstrated that C6 is convex (appendix C). C6 is the intersection of C6ij where
2 2 C6ij = F = (fx ; fy ) j dfdxx (i; j) + dfdxy (i; j) +
dfx (i; j) dy
2
+
dfy (i; j) dy
2
KS (i; j) 2
(24)
Adopting the following pixel labeling: (fx12 ; fy12) (fx22 ; fy22) (25) (fx11 ; fy11) (fx21 ; fy21) and letting: R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (26) it can be shown (appendix D) that the local projection operator P6ij of F onto C6ij is realized by:
P6ij: The projection operator to restore a vector eld subjected to the smoothness constraint of Eq. 24 is de ned by P6ij (F) = G = (gx ; gy ) (27) where G = F if R2 KS2 and gx22 = (fx22 (KS2 =R2 + 1) + fx11(1 ? KS2 =R2))=2 gx11 = (fx11 (KS2 =R2 + 1) + fx22(1 ? KS2 =R2))=2 gx21 = (fx21 (KS2 =R2 + 1) + fx12(1 ? KS2 =R2))=2 gx12 = (fx12 (KS2 =R2 + 1) + fx21(1 ? KS2 =R2))=2 (28) gy22 = (fy22(KS2 =R2 + 1) + fy11 (1 ? KS2 =R2))=2 2 2 2 2 gy11 = (fy11(KS =R + 1) + fy22 (1 ? KS =R ))=2 gy21 = (fy21(KS2 =R2 + 1) + fy12 (1 ? KS2 =R2))=2 gy12 = (fy12(KS2 =R2 + 1) + fy21 (1 ? KS2 =R2))=2 otherwise. It will be shown in section 3 that this method of local projections corresponds to projecting onto C6. The optical ow velocity eld between two images is obtained by rst computing the component of movement in the direction of the brightness gradient of the image (Eq. 15) and then by iteratively alternating projections of this component onto C5 and C6 by the corresponding projection operators P5 and P6 instead of by solving Eq. 19 and 20 iteratively for each pixel. The computation time of the two methods is about the same. the advantages of POCS will become apparent after the following development. 9
2.3 Theoretical comparison of POCS with the Horn and Schunck method
This comparison can be achieved by rst demonstrating that the Horn and Schunck algorithm (Eq. 15 to 20) corresponds to an unique projection of the zero vector eld onto C5 by using a non-Euclidian norm. Consider a Hilbert space H with elements F; G; H, etc : : :and an inner product between a pair of vector F(fx ; fy ) and G(gx ; gy ) de ned by Z Z @fx @gx + @fx @gx + @fy @gy + @fy @gy dx dy (29) F:G = @x @x @y @y @x @x @y @y and the corresponding norm by Z Z
2 2 2 2 @f @f @f @f x x y y kF k = @x + @y + @x + @y dx dy (30) Now, the rule which assign to every G 2 H its nearest neighbor P5 (F) in C5 de nes the non-linear projection operator P5 : H ?! C5 unambiguously by means of the minimality criterion 2
kF ? P5(F)k = Gmin kF ? Gk 2C5
(31)
The problem is to nd G 2 C5 which satis es this equation. Using the norm 30, Eq. 31 can be written as: Z Z
@(fx ? gx ) 2 + @(fx ? gx) 2 + @x @y 5 5 2 @(fy ? gy ) + @(fy ? gy ) 2 dx dy (32) @x @y Letting F be the null vector eld, then Eq. 32 becomes: min kF ? Gk2 = Gmin G2C 2C
Z Z
@gx 2 + @gx 2 + @gy 2 + @gy 2 dx dy (33) min k G k = min G2C5 G2C5 @x @y @x @y Note that this equation is exactly the smoothness constraint introduced by Horn and Schunck (Eq. 17). In other words, the Horn and Schunck approach corresponds to projecting a null vector eld onto C5 by using the norm of Eq. 31. The problem with this approach is that it only nds a function G in C5 which minimizes Eq. 17 the smoothness constraint. In most cases this constraint is too severe because the only eld which truly minimizes it is a translational eld. The restoration of the optical ow by the method of POCS using P5 and P6 presents the advantage that the smoothness can be restricted to be equal to or less than 2
10
KS (i; j) instead of being minimized. The value of KS (i; j) can even be \tuned" until the algorithm converges. This can be done because, if the intersection of C5 and C6 is an empty set, the alternate iterative projections onto C5 and C6 will oscillate between two values. The amplitude of this oscillation is equal to the minimum distance between the two sets. This distance can be successively reduced by modifying the KS (i; j) and, as soon as the intersection of C5 and C6 is non-empty, the algorithm will converge.
2.4 A general method to develop projection operators
Many other projection operators for special applications can be developed by the method of Lagrange multipliers, by rst demonstrating that the constraint condition de nes a convex set and then by nding the shortest distance between F and G subject to this constraint (see also appendices A, B and D). To illustrate this method consider the restoration of a vector eld whose ow is known along one axis. The corresponding convex set is de ned by: C7 = fF = (fx ; fy ) j
n X i=1
fx (i; j) = K(j)g
(34)
which, since summation is a linear operator, is convex. We want to nd the shortest distance: X (fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 (35) i;j
subject to the constraint condition: n X i=1
gx(i; j) = K(j)
(36)
According to the method of Lagrange multipliers, we consider: n X n ? X i=1 j =1
(fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 +
n X j =1
w(j)
n X i=1
gx(i; j) ? K(j)
!
(37)
Setting the partial derivative with respect to gx and gy equal to zero, we nd: ?2(fx (i; j) ? gx (i; j)) + w(j) = 0 (38) ?2(fy (i; j) ? gy (i; j)) = 0 Solving these equations we nd that: w(j) = 2 (
n X i=1
fx (i; j) ? K(j))
11
(39)
and the projection operator of F onto C7 is realized by: P7(F) = G = (gx; gy ) where:
(40)
P gx (i; j) = fx (i; j) ? ( ni=1 fx (i; j) ? K(j)) gy (i; j) = fy (i; j)
3 A linear time implementation of global projection operators Equations (7,9,11,13) of the rst four projection operators, once written in matrix form, can be solved by the method of Hockney [11]. For an N = n n pixel image, a system of (n ? 1) (n ? 1) equations has to be solved at each iteration. For a large image, this requires a large memory (2N) and long computation time (O(KN log2 N); K being the number of iteration). The smoothness constraint as stated in Eq. 17 also requires the solution of a system of (n ? 1) (n ? 1) equations. Fortunately, this problem can be circumvented. A constraint being de ned on all the domain of a vector eld F, is also valid for each subregions Ri of (F). The constraint restricted to a region Ri de nes a convex set Ci. If the union of the regions Ri is equal to (F), then the intersection of all Ci is equal to C, that is \ [ (41) Ri = (F) =) Ci = C i
i
Alternating projections of F onto each Ci by the corresponding projection operators Pi converge to a eld which belongs to C. If the Ri are m m pixels regions then the projection operators Pi requires the solution of a system of (m ? 1)2 equations. Solving (n ? 1)2 =(m ? 1)2 systems of (m ? 1)2 equations is faster and easier than solving one system of (n ? 1)2 equations. If m = 2, (m ? 1)2 independent equations have to be solved. For m = 2, the local projection operator for the divergence, the rotational and the hyperbolic components are realized by computing for each points: gx11 = fx11 + (divF ? D)=4 gx12 = fx12 + (divF ? D)=4 gx21 = fx21 ? (divF ? D)=4 gx22 = fx22 ? (divF ? D)=4 (42) gy11 = fy11 + (divF ? D)=4 gy12 = fy12 ? (divF ? D)=4 gy21 = fy21 + (divF ? D)=4 gy22 = fy22 ? (divF ? D)=4
12
gx11 = fx11 ? (rotF ? R)=4 gx12 = fx12 + (rotF ? R)=4 gx21 = fx21 ? (rotF ? R)=4 gx22 = fx22 + (rotF ? R)=4 (43) gy11 = fy11 + (rotF ? R)=4 gy12 = fy12 + (rotF ? R)=4 gy21 = fy21 ? (rotF ? R)=4 gy22 = fy22 ? (rotF ? R)=4 gx11 = fx11 + (hyp1 F ? H1)=4 gx12 = fx12 + (hyp1 F ? H1)=4 gx21 = fx21 ? (hyp1 F ? H1)=4 gx22 = fx22 ? (hyp1 F ? H1)=4 (44) gy11 = fy11 ? (hyp1 F ? H1)=4 gy12 = fy12 + (hyp1 F ? H1)=4 gy21 = fy21 ? (hyp1 F ? H1)=4 gy22 = fy22 + (hyp1 F ? H1)=4 gx11 = fx11 + (hyp2 F ? H2)=4 gx12 = fx12 ? (hyp2 F ? H2)=4 gx21 = fx21 + (hyp2 F ? H2)=4 gx22 = fx22 ? (hyp2 F ? H2)=4 (45) gy11 = fy11 + (hyp2 F ? H2)=4 gy12 = fy12 + (hyp2 F ? H2)=4 gy21 = fy21 ? (hyp2 F ? H2)=4 gy22 = fy22 ? (hyp2 F ? H2)=4 These equations were obtained by setting n = 2 in Eqs. 8, 10,12 and 14. The performance of this local projection method will be illustrated by an example: a linear hyperbolic eld de ned on 16 16 pixels has to be restored from its normal components only by using the a priori knowledge that its divergence is zero. Thus, two convex sets can be used: C1, the set of divergence free vector elds, and C5, the set of vector elds with known normal components. The hyperbolic eld was restored in three ways. The rst restoration was done by projecting alternatively onto C1 and C5. Projecting onto C1 , amounts to solving a system of 16 16 equations (this was done by the method of Hockney). The second restoration was done by considering C1, as the intersection of 9 convex sets of vectors elds that are divergence free on 6 6 pixel regions. This amounts to solving 9 systems of 5 5 equations. The third restoration was done by considering C1 as the intersection of 15 15 convex sets that are divergence free on 2 2 pixel regions. The time per iteration of this latter method is less than for the other two methods. All the restorations were done with a relaxation parameter equal to 1.9. The three methods are compared by computing the
13
percent error ek at the kth iteration [5, 8]: Fk k (46) ek = 100 k Hk ? Hk here Fk is the vector eld restored at iteration k and H is the original hyperbolic eld. From gure 1 it can be seen that the percent error curves of the three methods are quite comparable. The second restoration curve is between the two others and is not shown for clarity. The memory required by the local projection method (2 pixels) is small and the computation time is only O(Kl N) instead of O(Kg N log2 N) as for the method of Hockney, where Kl and Kg are the number of iterations of each method. As can be seen from Fig 1, Kl is approximately equal to Kg . Therefore the time complexity has been reduced from O(KN log2 N) to O(KN), where K is the number of iterations. The reason why the local projection method is so ecient is that all 2 2 pixel regions that do not have pixels in common are independent, or orthogonal, convex sets, and that projecting onto orthogonal convex sets amounts to projecting directly onto their intersection (the global constraint). These results show the advantages of the local projection method which is moreover much easier to implement than the Hockney method.
4 Examples of application This section presents two examples of vector eld restoration by the method of POCS. Other simple projection operators will also be presented to show how they can be easily derived to face practical applications. EXAMPLE I: The restoration of velocity eld (optical ow). The method of POCS is not restricted to optical ow due to motion of rigid bodies but can also be used to restore any velocity elds describing, for example, the interframe motion of uids or deformable bodies. The method can thus be illustrated on a divergence free vector eld (Fig. 2) which may represent the motion of a uid. This eld has to be restored from the partial information shown in Figure 3 which represents the components of movement in the direction of the brightness gradient. The latter was generated randomly at each point and then smoothed by a 15 15 pixel mean lter. Figure 3 thus represents the movement in the direction of a correlated brightness gradient of a random image submitted to the motion described by the velocity eld of Figure 2. The vector eld of Figure 2 can be restored from Figure 3 by projection onto three convex sets: C1, the set of divergence free vector elds; C5 , the set of vector elds whose velocity component in the direction of the brightness gradient corresponds at each point to that shown in Figure 3; and C6, the 14
set of vector elds respecting the smoothness constraint. The corresponding projection operators were de ned in section 2 as P1 , P5 and P6 respectively. When the latter projection operator is used, the optimal values of Ks (see equation 28) has to be determined as explained below. In order to measure the contribution of each projection operator, restorations have been obtained: a) by projecting onto C5 and C6 ; b) by projecting onto C1 and C5; and c) by projecting onto C1, C5 and C6. (Restoring by projection onto C1 and C6 only is, of course, meaningless). These results are compared to those obtained with the Horn and Schunck method (Figs. 4, 5, 6) and are summarized in Figure 7 which presents the curves of percent errors as a function of the number of iterations for each restoration.
4.1 Restoration by projecting onto C5 and C6
The value of Ks (x; y) (see equation 28) is generally not known for real images. So when the P6 projection operator is used its optimal value has to be determined rst. If Ks (x; y) is set to a constant equal to the maximal expected value, the restored eld can not be smooth enough. Conversely, if Ks is set too low, the restored eld can be too smooth. Fortunately, this parameter can be set to the value which will restore the smoothest eld while respecting the other constraints (C5 in this case). Geometrically, this can be seen by realizing that the size of C6 increases with Ks . If Ks is too small, the intersection with C5 will be empty, and the POCS algorithm will oscillate between the two convex sets instead of converging to an intersection point. If Ks is too large, the intersection with C5 will also be too large, the point of convergence can be too far from the desired solution and the restored eld will not be smooth enough. The optimal value of Ks can be determined by measuring, after some iterations, the RMS error between the elds restored before and after projecting onto C6. A periodic large error means that the algorithm oscillates between C6 and C5 instead of converging to a point of intersection. This means that there is no eld satisfying the two constraints. The amplitude of the RMS error is proportional to the distance between the two sets and indicates by how much Ks has to be increased. If the two convex sets intersect, the RMS error will continuously decrease with the number of iterations instead of being periodic. The smallest Ks value for which C5 and C6 intersect will restore the smoothest eld respecting the constraint of C5. The optimal Ks value can thus be determined by starting with a small Ks and by progressively increasing it until the two convex sets intersect. For restoring the eld of Figure 2 by projecting onto C5 and C6 an optimal value of Ks = 0:4 was determined in this way. The restored eld is presented in Figure 4. The percent error curves (not to be confused with the RMS error) between the original and restored elds as a function of the number of iterations with Ks = 0:4 and Ks = 1:2 are presented in Figure 7.
15
4.2 Restoration by projecting onto C1 and C5
The results are improved by projecting onto C1 and C5 instead of C5 and C6 (Fig. 7), meaning that, in this particular example, the divergence free constraint (C1) is more useful than the smoothness constraint (C6 ).
4.3 Restoration by projecting onto C1, C5 and C6
The best results are, of course, obtained by projecting onto C1 , C5 and C6 (Fig. 5 and Fig. 7), since the eld to be restored is both divergence free and smooth. In this case the optimal Ks value, determined as before, was found to be 1:2. Restoration obtained by using the same projection operators but with Ks = 0:4, yields results comparable to those obtained by projecting onto C5 and C6 with Ks = 0:4, which means that the intersection of C6 (with Ks = 0:4) with C1 and C5 is empty, or that the elds satisfying both the constraints of C1 and C5 are less smooth than the elds satisfying C5 only. It is interesting to note that the maximum value of Ks in the original eld is 1:112).
4.4 Comparison with the Horn and Schunck method
The results obtained with the Horn ans Schunck method are shown in Figures 6 and 7. They are comparable to those obtained by projecting onto C5 and C6 with Ks = 0:4 (Figs. 4 and 7). This is not surprising because both methods use the same constraints, although they are not implemented in the same way. These two methods are comparable in terms of memory and computation time and both can easily be multi-grid implemented. The only di erence of POCS, in this case, is that it constraints the smoothness to be smaller or equal to Ks (equation 28), (whereas the Horn and Schunck method minimizes this constraint), and that it allows, by adjusting the Ks parameter, to nd the smoothest eld which satis es C5 . Projecting onto C1 instead of C6 improved slightly the results (Fig. 7) meaning that the divergence free constraint is a very useful constraint that can even replace the smoothness constraint in this example. Using this constraint together with C5 and C6 (Figs. 5 and 7) gives the best results because all the a priori knowledge about the original eld (Fig. 2) are then used. EXAMPLE II: The restoration of a plane ow The eld of Figure 8 is an example of a plane ow with stagnation point (Hiemenz ow) [12]. The ow arrives from the y axis, impinges on a at wall placed at y = 0, divides into two streams on the wall and leaves in both directions leading to a stagnation point at x = y = 0. The eld of Figure 8 can be restored by the method of POCS by using only the following informations: the divergence of the eld is zero, the ow (the line integral of the divergence) is known to be zero on the external boundaries of 16
the eld, the fy components of the eld are known only on the upper boundary, the eld is known to be zero at y = 0, the eld is known to be symmetric about the y axis, and the eld is known to vary smoothly (smoothness constraint). The projection operators corresponding to the divergence free and smoothness constraints are those described in section 2 as P1 and P6 . The other four projection operators are very easy to implement: after each projection onto C1 and C6 (Ks = 0:03) the ow on the external boundaries is set to zero by subtracting from each point on the boundaries the mean vector of the points on the boundaries, the points at y = 0 are set to zero, the fy components of the points on the upper boundaries are replaced by the known one, and the eld is made symmetric by averaging the vectors of each pair of points symmetric about the y axis. The eld restored after 100 iterations is shown in Figure 9 and the percent error per iteration in Figure 10. Restorations were also made without the smoothness, symmetry and ow constraints to show the relative importance of these factors. From Figure 10, it can be seen that the symmetry constraint does not improve the restoration very much, but that the smoothness constraint is very ecient. The ow constraint, which is theoretically redundant because it is implied by the divergence free constraint, has the interesting property to reduce the percent error in the rst iterations.
5 Conclusion The main contribution of this paper was to show that the method of POCS, as originally developed by Youla and Webb [3, 4] for scalar elds (or images) can be applied to the restoration of vector elds. Projection operators were developed to restore elds of known divergence, rotational, or hyperbolic components (P1 to P4). It was shown how these projection operators can eciently be implemented by the method of local projections (section 3). With this method, computations which rst appear to be prohibitively long are reduced to a very manageable size. This method can be used whenever a convex set can itself be represented as the intersection of many local convex sets and will be very useful in developing other projection operators. Projection operators to restore vector elds from their inner product with a known eld, or to restore elds subjected to a smoothness constraint were also presented (P5 and P6 ). These two last projection operators correspond to the a priori knowledge used by Horn and Schunck to determine the optical ow. Finally, a general method to derive new projection operators was proposed. When applied to machine vision problems, such as the determination of the optical ow, which are generally solved by the calculus of variations, the method of POCS presents some interesting aspects: 1) the constraints are introduced one by one and are implemented by nding the corresponding projection operator. This is easier than nding the Euler equations of a complex functional 17
which includes all the constraints; 2) if a constraint must be smaller or equal than a value K (instead of being only equal to K) it is possible to nd the optimal K which also satis es the other constraints by checking, after some iterations, if the algorithm is not oscillating between two convex sets instead of converging to a point which belongs to their intersection (see section 4, example 1); 3) since the algorithm converges to a point of intersection of all the convex sets, the restored eld really has all the properties described by each of these convex sets; 4) some a priori knowledge, for example, that the Fourier transforms of the vector components are band-limited, or that the vector components assume prescribed values over some pixels of the images, are easier to formulate in terms of POCS than with variational approaches. Of course, the more a priori knowledge introduced in the algorithm as convex sets, the more precise the restoration. Examples of applications of this method on real noisy images (echographic textures and two-dimensional echocardiograms) are presented in reference [13].
18
Figure 1: Hyperbolic eld restoration. Error versus iteration number by solving one system of 15 15 equations, and 15 15 systems of one equation. Figure 2: A simulated divergence free vector eld Figure 3: Projection of the vector eld shown in Figure 2 along a correlated brightness gradient Figure 4: Vector eld restored by projecting onto C5 and C6(Ks = 0:4) Figure 5: Vector eld restored by projecting onto C1 , C5 , and C6 Figure 6: Vector eld restored by the Horn and Schunck method Figure 7: Error versus iteration number by projecting onto C5 and C6 (Ks = 1:2), onto C5 and C6(Ks = 0:4), onto C1 and C5, onto C1 and C5 and C6 (Ks = 1:2), and by Horn and Schunck method Figure 8: A plane ow with stagnation point Figure 9: Vector eld restored after 100 iterations Figure 10: Error versus iteration number by restoring with the six constraints, with ve constraints (without the symmetry constraint or without the ow constraint), and with four constraints (without the symmetry and smoothness constraints) 19
A Demonstration for P1 A vector eld has to be restored from the partial information F = (fx ; fy ) and the a priori knowledge that its divergence must be D(i; j). The problem is thus to nd the shortest distance between F and C1; that is, the minimum of: n X n X i=1 j =1
(fx (i; j) ? gx(i; j))2 + (fy (i; j) ? gy (i; j))2
(47)
subject to the constraint: divG(i; j) = D(i; j) According to the method of Lagrange multipliers we consider: n X n X i=1 j =1
(48)
(fx (i; j) ? gx (i; j))2 + (fy (i; j) ? gy (i; j))2 +
n X n X i=1 j =1
w(i; j) (divG(i; j) ? D(i; j))
(49)
By setting to zero the partial derivatives with respect to gx and gy we obtain: 2(fx (i; j) ? gx (i; j)) = +w(i ? 1; j) + w(i ? 1; j ? 1) ? w(i; j) ? w(i; j ? 1) 2(fy (i; j) ? gy (i; j)) = ?w(i; j) ? w(i ? 1; j) + w(i; j ? 1) + w(i ? 1; j ? 1) (50) where w(0; j) = w(i; 0) = w(i; n) = w(n; j) = 0 for i = 1; 2; :::;n and j = 1; 2; :::;n. Adding these equations we nd: divF(i; j) ? D(i; j) = 4 w(i; j) ?(w(i ? 1; j ? 1) + w(i ? 1; j + 1) + w(i + 1; j ? 1) + w(i + 1; j + 1) (51)
B Demonstration for P5 We want to nd, on each point, the minimum of (fx ? gx)2 + (fy ? gy )2
(52)
subject to the constraint gx :nx + gy :ny + KN = 0 According to the method of Lagrange multipliers we consider (fx ? gx)2 + (fy ? gy )2 + w(gx :nx + gy :ny + KN ) 20
(53) (54)
By setting to zero the partial derivatives with respect to gx and gy we obtain: 2(?fx + gx ) + w:nx = 0 (55) 2(?fy + gy ) + w:ny = 0 Solving these equations we nd that: +fx :nx + fy :ny + KN = w(n2x + n2y )=2 (56) the solution is thus: gx = fx ? nx(fx :nx + fy :ny + KN )=(n2x + n2y ) (57) gy = fy ? ny (fx :nx + fy :ny + KN )=(n2x + n2y )
C Convexity of C6 Let g1 and g2 belong to C6 . we have to prove that 8 2 [0; 1]; g1 + (1 ? )g2 2 C6 g1, g2 2 C6 implies that the two points: ( @g@x1x ; @g@x1y ; @g@y1x ; @g@y1y ) and ( @g@x2x ; @g@x2y ; @g@y2x ; @g@y2y ) belong to the hypersphere:
(58)
x2 + y2 + z 2 + t2 KS2 : (59) Since an hypersphere is convex, the point with the following coordinates ( @f@x1x + (1 ? ) @f@x2x ; @f@y1x + (1 ? ) @f@y2x ; @f@x1y + (1 ? ) @f@x2y ; @f@y1y + (1 ? ) @f@y2y ) (60) also belongs to the hypersphere and we can write: 2 2 @f @f @f @f 1x 2x 2x 1x + + (1 ? ) + (1 ? ) @x @x @y @y 2 2 @f @f @f @f 1y 2y 2y 1y + @x + (1 ? ) @x + @y + (1 ? ) @y KS2 so C6 is convex.
(61)
D Demonstration for P6ij Using the same pixel labeling as in Eq. 25, we want to nd the shortest distance: (fx11 ? gx11)2 + (fx12 ? gx12)2 + (fx21 ? gx21)2 + (fx22 ? gx22)2 +(fy11 ? gy11)2 + (fy12 ? gy12)2 + (fy21 ? gy21)2 + (fy22 ? gy22)2 (62) 21
subject to the constraint condition: (gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 KS2 (63) This constraint can be written: (gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 = L2 (64) with 0 L KS . According to the method of Lagrange multipliers we consider: (fx11 ? gx11)2 + (fx12 ? gx12)2 + (fx21 ? gx21)2 + (fx22 ? gx22)2 +(fy11 ? gy11)2 + (fy12 ? gy12)2 + (fy21 ? gy21)2 + (fy22 ? gy22)2 +w((gx22 ? gx12 + gx21 ? gx11)2 + (gx22 ? gx21 + gx12 ? gx11)2 +(gy22 ? gy12 + gy21 ? gy11)2 + (gy22 ? gy21 + gy12 ? gy11)2 ? L2 ) (65) By setting to zero the partial derivatives with respect to gx and gy we obtain: ?fx11 + gx11 + 2w(+gx11 ? gx22) = 0 ?fx22 + gx22 + 2w(+gx22 ? gx11) = 0 ?fx12 + gx12 + 2w(+gx12 ? gx21) = 0 ?fx21 + gx21 + 2w(+gx21 ? gx12) = 0 (66) ?fy11 + gy11 + 2w(+gy11 ? gy22) = 0 ?fy22 + gy22 + 2w(+gy22 ? gy11) = 0 ?fy12 + gy12 + 2w(+gy12 ? gy21) = 0 ?fy21 + gy21 + 2w(+gy21 ? gy12) = 0 By substitution of equations 66 into Eq. 62 we see that we have to minimize w2. From equations 66 we can write: fx22 ? fx12 + fx21 ? fx11 = (1 + 4w)(gx22 ? gx12 + gx21 ? gx11) fx22 ? fx21 + fx12 ? fx11 = (1 + 4w)(gx22 ? gx21 + gx12 ? gx11) (67) fy22 ? fy12 + fy21 ? fy11 = (1 + 4w)(gy22 ? gy12 + gy21 ? gy11) fy22 ? fy21 + fy12 ? fy11 = (1 + 4w)(gy22 ? gy21 + gy12 ? gy11) By adding these equations we obtain: R2 = (1 + 4w)2L2 (68) where R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (69) To solve Eq. 68 for 0 L KS and w2 minimum we have to consider two cases: If R2 KS2 then Eq. 68 is veri ed for w = 0, R2 = L2 and P6(F) = 22
G = (gx; gy ) = F. In the other case, w2 is minimized for L = KS and w = (R2=KS2 ? 1)=4. And the solution is: gx22 = (fx22 (KS2 =R2 + 1) + fx11(1 ? KS2 =R2))=2 gx11 = (fx11 (KS2 =R2 + 1) + fx22(1 ? KS2 =R2))=2 gx21 = (fx21 (KS2 =R2 + 1) + fx12(1 ? KS2 =R2))=2 gx12 = (fx12 (KS2 =R2 + 1) + fx21(1 ? KS2 =R2))=2 (70) gy22 = (fy22(KS2 =R2 + 1) + fy11 (1 ? KS2 =R2))=2 gy11 = (fy11(KS2 =R2 + 1) + fy22 (1 ? KS2 =R2))=2 gy21 = (fy21(KS2 =R2 + 1) + fy12 (1 ? KS2 =R2))=2 gy12 = (fy12(KS2 =R2 + 1) + fy21 (1 ? KS2 =R2))=2 R2 = +(fx22 ? fx12 + fx21 ? fx11 )2 + (fx22 ? fx21 + fx12 ? fx11)2 +(fy22 ? fy12 + fy21 ? fy11 )2 + (fy22 ? fy21 + fy12 ? fy11)2 (71) Acknowledgment: this work was supported in part by the National Sciences and Engineering Research council of Canada and the Canadian Heart Foundation through a Research Scholarship to G. E. Mailloux, and in part by the NIH research grant no. NS22407.
References [1] R.W. Gerchberg, Super-resolution through error energy reduction, Opt. Acta, vol.21, pp. 709-720, 1974. [2] A. Papoulis, A new algorithm in spectral analysis and band-limited extrapolation, IEEE Trans. Circuits Syst., vol. CAS-22, pp. 735-742, 1975. [3] D.C. Youla, Generalized image restoration by the method of alternating projections, IEEE Trans. Circuits Syst., vol. CAS-25, pp. 694-702, 1978. [4] D.C. Youla and H. Webb, Image restoration by the method of convex projections: Part 1:-Theory, IEEE Trans. Med. Imaging, vol. MI-1, pp. 81-94, 1982. [5] M.I. Sezan and H. Stark, Image restoration by the method of convex projections: Part 2-Application and numerical results, IEEE Trans. Med. Imaging, vol. MI-1, pp. 99-101, 1982. [6] M.I. Sezan and H. Stark, Tomographic image reconstruction from incomplete view data by convex projection and direct Fourier inversion, IEEE Trans. Med. Imaging, vol. MI-3, pp. 91-98, 1984. [7] R.T. Chin, C-L Yeh and W.S. Olson, Restoration of multichannel microwave radiometric images, IEEE Trans. Pattern Anal. Machine Intell., vol. PAMI-7, pp.475-484, 1985. 23
[8] P.Y. Simard and G.E. Mailloux, A projection operator for the restoration of divergence free vector elds, IEEE Trans. Pattern Anal. Machine Intell. PAMI-10, pp. 248-256, 1988. [9] B.K.P. Horn and B.G. Schunck, Determining optical ow, Artif. Intell., vol. 17, pp. 185-203, 1983. [10] E.C. Hildreth, Computation underlying the measurement of visual motion, Artif. Intell., vol. 23, pp. 309-354, 1984. [11] R. Vichnevetsky. Computer Methods for Partial Di erential Equations. vol. 1. Englewoods Cli s, Prentice-Hall, 1981. [12] H. Schlichting. Boundary Layer Theory. New-York: mcGraw-Hill, 1960. [13] G. E. Mailloux, F. Langlois, P. Y. Simard and M. Bertrand, Restoration of the velocity eld of the heart from two-dimentional echocardiograms. IEEE Trans. Med. Imaging, vol. 8, 2, pp.143-153, 1989.
24
