09.24.filling-in By Joint Interpolation Of Vector fields And Grey Levels

Filling-In by Joint Interpolation of Ve tor Fields and Gray Levels

C. Ballester,  M. Bertalmio, y V: Caselles; G: Sapiro;y

and

J: Verdera;

De ember 10, 2000

Abstra t

A variational approa h for lling-in regions of missing data in digital images is introdu ed in this paper. The approa h is based on joint interpolation of the image gray-levels and gradient/isophotes dire tions, smoothly extending in an automati fashion the isophote lines into the holes of missing data. This interpolation is omputed by solving the variational problem via its gradient des ent ow, whi h leads to a set of oupled se ond order partial dierential equations, one for the gray-levels and one for the gradient orientations. The pro ess underlying this approa h an be onsidered as an interpretation of the Gestaltist's prin iple of good ontinuation. No limitations are imposed on the topology of the holes, and all regions of missing data an be simultaneously pro essed, even if they are surrounded by

ompletely dierent stru tures. Appli ations of this te hnique in lude the restoration of old photographs and removal of superimposed text like dates, subtitles, or publi ity. Examples of these appli ations are given. We on lude the paper with a number of theoreti al results on the proposed variational approa h and its orresponding gradient des ent ow. Keywords:

Interpolation, lling-in, image gradients, image gray-levels, variational approa h,

partial dierential equations, Gestalt prin iples. EDICS: 2-INTR  Dept.

de Te nologia, University of Pompeu-Fabra, Passeig de Cir umvala io, 8, 08003 Bar elona, Spain,

f oloma.ballester,vi ent. aselles,joan.verderagte n.upf.es y Department

of Ele tri al and Computer Engineering, University of Minnesota, Minneapolis, MN 55455, USA,

fmar elo,guillege e.umn.edu

1

1 Introdu tion

Filling-in missing data in digital images has a number of fundamental appli ations. They range from removing obje ts from a s ene all the way to re-tou hing damaged paintings and photographs. The basi idea is to ll-in the gap of missing data in a form that it is non-dete table by an ordinary observer. In art, this pro ess is alled inpainting [39, 16, 25, 6, 7℄. Sin e the early days of art and photography, lling-in and inpainting has been done by professional artists. Imitating their performan e with semi-automati digital te hniques is urrently an a tive area of resear h, see for example [23, 27, 28, 8℄ and the papers dis ussed below. The goal of this work is to introdu e a novel algorithm for automati ally lling-in gaps in the image. In this arti le we follow the suggestions in the on lusions se tion in [6℄ and introdu e an energy fun tional based on an interpretation of Gestaltist's good ontinuation prin iple. Suppose that we are given an image u0 : R ! [a; b℄, where R is a square of IR2 and a < b, and is an open bounded subset of R with Lips hitz ontinuous boundary. We shall all the hole or gap. We want to ll-in the hole based on the geometri and photometri information outside the hole. For that we use what we all a band around , i.e., we onsider an open region ~ of R su h that  ~ ( is the losure of the set). The band we refer to will be the set B = ~ n . To ll-in the hole we use the information of u0 ontained in B , mainly gray level and the ve tor eld of isophotes (level sets) dire tions of u0 in B .1 We attempt to ontinue the level sets of u in B inside taking into a

ount the prin iple of good ontinuation. We propose an energy fun tional whi h takes into a

ount these prin iples interpreted in a suitable way. The energy fun tional we propose has to be minimized with respe t to two variables : a ve tor eld  whi h represents the dire tions of the level lines of u, and the gray level u.  and u are onstrained in the band B by their known values there. The use of the ve tor eld of dire tions  is one of the main points of the algorithm presented in this paper, whi h permits the level sets to smoothly

ontinue inside the hole. We are then ontinuing both the geometri and photometri properties of the image inside the hole. Let us nally say that the only user intera tion required by the algorithm here introdu ed is 1

The width of the band is su h that it onveys the boundary information, mainly the gray lavel and isophotes

dire tion, and numeri aly it depends on the minimal sten il needed by the implementation.

2

to mark the regions to be lled-in. Although a number of te hniques exist for the semi-automati dete tion of image defe ts (mainly in lms), addressing this is not part of the s ope of this paper. Sin e the algorithm here presented an be used not just to restore damaged photographs but also to remove undesired obje ts and writings on the image, the regions must be marked by the user, sin e they depend on his/her subje tive sele tion. 1.1 Closely related approa hes

Before pro eeding with the detailed des ription of our algorithm, let us omment on related work. Note that image denoising is dierent to lling-in, sin e the regions of missing data are usually large. That is, regions o

upied by top to bottom s rat hes along several lm frames, long ra ks in photographs, superimposed large fonts, and so on, are of signi ant larger size than the type of noise assumed in ommon image enhan ement algorithms. In addition, in

ommon image enhan ement appli ations, the pixels ontain both information about the real data and the noise (e.g., image plus noise for additive noise), while in our appli ation there is no signi ant information in the region to be inpainted. A very a tive area related to the work here presented is the restoration of damaged lms. The basi idea here is to use information from past and future frames to restore the urrent one, e.g., [23, 28℄. Of ourse, this general approa h an not be used when dealing with still images. Another area related to the work here des ribed is texture synthesis. The basi idea here is to sele t a texture and synthesize it inside the region to be lled-in (the hole). Although outstanding texture synthesis results have been reported in the literature, e.g., [22, 14, 19, 36℄, these algorithms require the user to sele t the texture to be opied into the hole. For images where the region to be repla ed overs several dierent stru tures, the user would need to go through the tremendous work of segmenting them and sear hing orresponding repla ements throughout the pi ture. Although part of this sear h an be done automati ally, this is extremely time onsuming and requires the non-trivial sele tion of many riti al parameters, e.g., [14℄. Last, a number of fundamental works on diso

lusion and line ontinuation have been reported in the literature, and these are the losest to our approa h. A pioneering ontribution in this area is des ribed in [33℄. The authors presented a te hnique for removing o

lusions with 3

the goal of image segmentation. Sin e the region to be lled-in an be onsidered as o

luding obje ts, removing o

lusions is analogous to image inpainting. The basi idea suggested by the authors is to onne t T-jun tions at the o

luding boundaries of obje ts with elasti a minimizing urves (see later in this paper for the exa t de nition of elasti a urves). The te hnique was mainly developed for simple images obtained from a segmentation, with only a few obje ts with onstant gray-levels. Thus, they ended up by onne ting T -jun tions at the same gray level. (Other resear hers, e.g., D. Ja obs, R. Basri, S. Zu ker, et , have followed this interesting resear h area, mainly developing te hniques for smooth urve ontinuation.) Masnou and Morel [31, 32℄ re ently extended these ideas, presenting a very elegant and inspiring formal variational formulation for diso

lusion and a parti ular pra ti al algorithm implementing some of the ideas in this formulation. The algorithm lls-in by joining with geodesi urves the points of the isophotes arriving at the boundary of the region to be inpainted. The holes in their algorithm are limited to having simple topology. In addition, the angle with whi h the level lines arrive at the boundary of the holes are not (well) preserved, and the algorithm uses straight lines to join equal gray value pixels. The present work was motivated, in part, by their work. In some sense, as shown by the work of Masnou and Morel [31, 32℄ and the present paper, operators as the one proposed in [33℄ an be used to interpolate level-sets, and then re onstru t from them the interpolated gray-level image. Re ently, we have addressed the on ept of smooth ontinuation of information in the levellines dire tion in [6, 7℄. We proposed an algorithm, inspired in partial dierential equations, that propagates the image Lapla ian in this dire tion. The algorithm attempts to imitate basi approa hes used by professional restorators. The algorithm also introdu es the importan e of propagating both the gradient dire tion (geometry) and gray-values (photometry) of the image in a band surrounding the hole to be lled-in. It is part of the goal of the urrent paper to adopt some of the ideas of [6, 7℄, while deviating from the parti ular model in order to be able to de ne a formal variational approa h to the lling-in problem. The work in [6, 7℄ inspired a very elegant approa h to the lling-in problem re ently reported in [9℄ (this work was performed independently to the one reported in this paper).2 The authors 2

We thank the authors for providing us with a preliminary report of their work.

4

present a lear and very intuitive axiomati approa h to the problem. The main algorithm they propose after an interesting dis ussion of the inpainting problem is to minimize the Total Variation (TV), [35℄, of the image inside the hole (they also use, as proposed in [6, 7℄ and further studied in this paper, a band surrounding the region). They address in addition the interpolation and lling-in in the presen e of noise, a very important additional ontribution. As in the work of Masnou and Morel, their interpolation is limited to reating straight isophotes, not ne essarily smoothly ontinued from the hole boundary, and mainly is developed (as the authors learly state) for small holes. Although straight onne tions give visually pleasant results for small holes, it is important to develop a theory that permits interpolation of level lines a ross large gaps, where onne ting with straight lines will be unpleasent even for simple images. As we will argue later in this paper, in order to obtain su h a smooth interpolation and ontinuation of isophotes, it is ne essary to go into high order PDE's or systems of PDE's, as done in [6, 7℄ and here. Note that onsidering the angle of arrival of the level lines at the gap, and pursuing a smooth interpolation of it as done here, it is also supported by resear h in per eption, from the Gestalt to more re ent work, e.g., [34℄. To on lude this se tion, the interested reader is referred to the works of Nitzberg-MumfordShiota, Masnou-Morel, and Chan-Shen (as well as our previous work) to study other interesting and very related te hniques for lling-in. Let us on lude with the plan of the paper. In Se tion 2 we introdu e the problem, the fun tional spa es and the energy fun tional for image inpainting. To larify the meaning of the fun tional, we also dis uss the parti ular ase where we interpolate the gray level, knowing the ve tor eld of dire tions. Se tion 3 is devoted to numeri al experiments. Se tion 4 ontains some on lusions. Finally, Appendix 1 and 2 are devoted, respe tively, to the proof of existen e of minimizers for the energy fun tional introdu ed in this paper and to details on the numeri s for the proposed equations.

5

2 Joint interpolation of ve tor elds and gray values

Let u0 : R ! IR be an image de ned on a domain R of IR2, whi h we may suppose to be a square. Let ; ~ be two open bounded domains in IR2 with Lips hitz boundary and suppose that  ~  R. To simplify our presentation we shall assume that ~ does not tou h the boundary of the image domain R. Let B := ~ n . B will be alled the band around . Suppose that a fun tion u0 is given in B whi h, for the moment being, we shall assume to be smooth in the losure, B , of B (later we shall assume that u0 is of bounded variation, i.e., u0 2 BV (B )). Let 0 be the ve tor eld of dire tions of the gradient of u0 on B , i.e., 0 is a ve tor eld with values in IR2 satisfying 0(x)
ru0(x) = jru0 (x)j and j0(x)j  1 (ideally 1 or 0, see below). This is the information we shall use on B . We pose the image inpainting problem in the following form: Can we extend (in a reasonable way) the pair of fun tions (u0; 0 ) from the band ~ n to a pair of fun tions (u; ) de ned inside ? Of ourse, we will have to pre ise what we mean by a reasonable way. We shall dis uss and analyze a variational formulation of this lling-in problem and dis uss possible energy fun tionals, and their orresponding gradient des ent ows, whi h give a solution to it. The data are given on the band B and we should onstraint the solution (u; ) to be near the data on B . The ve tor eld  should satisfy jj  1 on and should be related to u by trying to impose that 
ru = jruj, i.e., we should impose that  is related to the ve tor eld of dire tions of the gradient of u. The ondition j(x)j  1 should be interpreted as a relaxation of this. Indeed, it may happen that (x) = 0 ( at regions) and then we annot normalize the ve tor eld to a unit ve tor. We should have in mind that the ideal ase would be that  = jrr j , u being a smooth fun tion with ru(x) 6= 0 for all x 2 . Finally, we should impose that the ve tor eld 0 in the band tries to ontinue smoothly to  inside . We shall impose this by observing that, in ase  represents the dire tions of the normals to the level lines of u, i.e., of the urves u(x1 ; x2) = ,  2 IR, then a term like div() represents its urvature. Motivated by the prin iple of smooth ontinuation, our energy fun tional should ontain terms integrating div(). Indeed, olle ting all the observations above, we propose to minimize a fun tional of the u u

6

form

Z

Z

Minimize ~ jdiv()j (a + bjrk  uj)dx +  ~ (jruj p

(1)


ru)dx;

where a, b, and  are positive onstants, k is a smoothing kernel, and ; u are submitted to

onstraints. We need to give a sense to the integrals appearing in the above expression and to make pre ise the admissible lass of fun tions where the fun tional has to be minimized. For that, we need to introdu e some fun tion spa es. This is done in the following se tion. On e this has been formally addressed, we will further dis uss the underlying on epts of the above fun tional. 2.1 Fun tion spa es

Let us rst re all the de nition of BV fun tions and total variation. Let Q be an open set. A fun tion u 2 L1(Q) whose partial derivatives in the sense of distributions are measures with nite total variation in Q is alled a fun tion of bounded variation. The lass of su h fun tions will be denoted by BV (Q). Thus u 2 BV (Q) if and only if there are Radon measures 1 ; : : : ;  de ned in Q with nite total mass in Q and N

Z

Q

uDi 'dx =

Z

Q

(2)

'di

for all ' 2 C01(Q), i = 1; : : : ; N . Thus the gradient of u is a ve tor valued measure with nite total variation Z

k ru k= supf u div ' dx : ' 2 C01(Q; IRn ); j'(x)j  1 for x 2 Qg:

(3)

Q

The spa e BV (Q) is endowed with the norm k u kBV =k u kL1 (Q) + k ru k :

(4)

We say that a measurable set E  Q has nite perimeter in Q if its indi ator fun tion  2 BV (Q). If u 2 BV (Q) almost all its level sets [u  ℄ = fx 2 Q : u(x)  g are sets of nite perimeter. For sets of nite perimeter E one an de ne the essential boundary   E , whi h is re tifable with nite H 1 measure, and ompute the normal to the level set at H 1 almost all points of   E . Thus at almost all points of almost all level sets of u 2 BV (Q) we may de ne E

N

N

7

a normal ve tor (x). This ve tor eld of normals  an be also de ned (hen e extended to all Q) as the Radon-Nikodym derivative of the measure ru with respe t to jruj, i.e., it formally satis es 
ru = jruj and, also, jj  1 a.e.. For further information on erning fun tions of bounded variation we refer to [1, 17, 40℄. Let us now introdu e the fun tion spa es for . Let Q be an open bounded subset of IR2 with a Lips hitz boundary. We de ne W 1;p(div; Q) = f 2 Lp(Q)2 : div() 2 Lp (Q)g;

and

M (div; Q) = f 2 L1 (Q)2 : div()

1  p < 1;

is a Radon measure in Qg:

The Tra e Theorem ([2℄,[10℄) guarantees that the normal omponent 
nj , is well de ned for ve tor elds  in W 1 (div; Q), or in M (div; Q). To simplify our notation we shall assume that W 1 1 (div; Q) represents the spa e M (div; Q). Next, we shall give a sense to the integrals of bounded ve tor elds with divergen e in L integrated with respe t to the gradient of a BV fun tion. For that, we shall need some results from [2℄ (see also [26℄ and [11℄). Let Q be an open bounded subset of IR with Lips hitz

ontinuous boundary. Let p  1 and q  1 be su h that 1 + 1 = 1. Following [2℄, let Q

;p

;

p

n

p

q

X (Q)p = fz 2 L1 (Q; IRn ) : div(z ) 2 Lp (Q)g:

(5)

If z 2 X (Q) and w 2 BV (Q) \ L (Q) we de ne the fun tional (z; rw) : C01(Q) ! IR by the formula < (z; rw); ' >= w ' div(z ) dx w z
r' dx: (6) q

p

Z

Z

Q

Q

Then (z; rw) is a Radon measure in Q, Z

Q

(z; rw) =

Z

Q

(7)

z
rw dx

for all w 2 W 1 1(Q) \ L (Q) and ;

p

Z B



Z

Z

(z; rw)  j(z; rw)j  kzk1 krwk B

B

for any Borel set B  Q. 8

(8)

In [2℄, a weak tra e on Q of the normal omponent of z 2 X (Q) is de ned. Con retely, it is proved that there exists a linear operator : X (Q) ! L1(Q) su h that p

p

k (z)k1  kzk1

for all x 2 Q if z 2 C 1(Q; IR );

(z )(x) = z (x)
n(x)

N

where n(x) denotess the outer unit normal at x 2 Q. We shall denote (z)(x) by [z; n℄(x). Moreover, the following Green's formula, relating the fun tion [z; n℄ and the measure (z; rw), for z 2 X (Q) and w 2 BV (Q) \ L (Q), is established: q

p

Z

Q

w div(z ) dx +

Z

Q

(z; rw) =

Z

Q

[z; n℄w dH

N

(9)

1:

If no onfusion arises, we shall denote z
rw instead of (z; ru) for z 2 X (Q) , w 2 BV (Q) \ L (Q). These results will be used in the proof of Theorem 2 given in Appendix 1. p

q

2.2 The energy derivation and interpretation

One of the key on epts above was the band around the hole. The band is of lo al hara ter but in prin iple it ould be extended to all the known part of the image. Obviously, what happens at distant parts an be independent or not from what happens at the hole, but, in our onstru tion below, we suppose that only a narrow band around the hole in uen es what happens inside the hole. Could we ll-in without the band ? To dis uss this suppose that we are given the image of Figure 1.a whi h is a gray band on a bla k ba kground partially o

luded by a square . We suppose that the sides of the square hole are orthogonal to the level lines of the original image. In these onditions, the normal omponent of the ve tor eld 0 outside is null at  . Thus if the boundary data is just 0
nj , we would have that 0
nj = 0. In parti ular, the ve tor eld  = 0 satis es this ondition. If we are not able to propagate  inside this may be ome an unpleasant situation, sin e this would mean that we do no propagate the values of u at the boundary. If we write the fun tional (1) with  = 0, it turns out to be the Total Variation [35℄. The de ision of extending the gray band or lling the hole with the bla k level would be taken as a fun tion of the perimeter of the dis ontinuities of the fun tion in the hole. Then the result of interpolating Figure 1.a using Total Variation would be that of Figure 1.b and not the one in 



9

Figure 1. , be ause the interpolating lines in Figure 1.b are shorter than the ones in Figure 1. . To over ome this situation we introdu e the band around the hole. The introdu tion of the band permits us to ee tively in orporate in the fun tional the information given by the ve tor eld  outside . In Figure 1.b we display the result of the interpolation with  = 0 on . In Figure 1. we display the result of the interpolation using the fun tional we shall ompletely des ribe below, whi h takes into a

ount the band B and omputes the ve tor eld  in ~ = [ B .

Figure 1.a.

Figure 1.b.

Figure 1. .

Thus, let B be a band around with a Lips hitz boundary ontaining the boundary of

(see Figure 2). As we made expli it above, B = ~ n . Given the band B and the fun tion u0 of bounded variation in B , we de ne the spa e BV ( ~ ; B; u0 ) = fu 2 BV ( ~ ) : u = u0

Let

! IR2 be a ve tor eld of dire tions of the gradient of u0, i.e., j0 j  1 and 0
ru0 = jru0 j as measures in B (therefore, a.e.). In pra ti e we shall onstraint the ve tor 0

:

in B g:

B

eld  to be the ve tor eld of dire tions of u0 only indire tly, through the fun tional. We ould also introdu e this as a onstraint or with a penalty term ( 0)2 (see also [9, 35℄ for penalties of this form in the gray values). R

B

Figure 2. 10

Combining the previous elements, the band, the relations between  and u, and the smoothness term on , we propose to interpolate the pair (; u) in by minimizing the fun tional: Z

Minimize ~ jdiv()j (a + bjrk  uj)dx  2 W 1 (div; ~ ) u 2 BV ( ~ ; B; u0 ) (10) jj  1; juj  ku0 k 1 ( ) jruj 
ru = 0in ~ 
nj ~ = 0
nj ~ ; where a;  > 0, b  0, k denotes a regularizing kernel of lass C 1 su h that k(x) > 0 a.e. and n(x) denotes the outer unit normal at x 2  ~ . The previous fun tional is oer ive and admits a minimum in the lass of fun tions des ribed above if p > 1. The ase p = 1 is under study. The fun tional an be interpreted as a formulation of the prin iple of good ontinuation and amodal ompletion as formulated in the Gestalt theory of vision. The following remarks

ontain heuristi arguments whi h may help to understand our hoi e. In next subse tion we shall explain in more detail the role of the term oupling  and u. p

;p

L



B



Remarks.

1. The onstant b is  0. If u is the hara teristi fun tion of the region en losed by a urve C then a term like jdiv()j jruj (11)

~ Z

p

is related to jj ds, where  is the Eu lidean urvature (of the level-sets). If p = 2, this term appears in Euler's elasti a, R

C

p

Z

C

( + 2 )ds;

; > 0:

(12)

Euler's elasti a (12) was proposed in [33℄ as a te hnique for removing o

lusions with the goal of image segmentation, sin e this riterion yields smooth, short, and not too urvy

urves. In terms of hara teristi fun tions, Euler's elasti a an be written as Z

ru jruj  + div jr uj



11



 

2

:

(13)

In [5℄, it was shown that this fun tional is not lower semi ontinuous. The fun tional proposed by Masnou and Morel [31℄, [32℄ an be interpreted as a relaxation of it, sin e it integrates fun tionals like the elasti a (plus the angle that the urve makes with the

orresponding level line arriving at the boundary) along the level lines of the fun tion u. Our fun tional an be also onsidered as a relaxed formulation of the energy of the elasti a. For that, we introdu ed  as a independent variable, and we tried to ouple it to u by imposing that 
ru = jruj (in ase ru(x) 6= 0, then (x) = jrr (( ))j ). Thus the main part of the fun tional writes jruj (a + b jdiv ()j ) : (14) u x u x

Z

p

with p  1, p = 2 orresponding to the elasti a. The hosen urve tries to minimize the length plus the angular variation of the urve. To onne t urves at distan e we should dismiss the length of the urve and hoose a to be small (but > 0, see next remark). Finally, let us say that for mathemati al reasons we have onvolved the ru term of (11) to be able to prove the existen e of a minimum for (10). From a theoreti al point of view, this may invalidate our previous omments. But, from a pra ti al point of view, it gives a weight to the urve of dis ontinuities of the image. 2. The onstant a has to be > 0. Otherwise we do not get an L bound on div(). Now, let us omment on the two terms ontaining div(). Heuristi ally, if we do not ompute  in a proper way, in a ontinuous image like in Figure 1,  ould be zero ex ept on a set of

urves. Then  = 0 a.e. on B (or on ~ ) and a term like p

Z

~

jdiv()jp dx

(15)

would produ e a null value sin e div() = 0. On the other hand, a term like (11) would integrate a power of the urvature on the level line orresponding to the dis ontinuity of the image and it would guarantee that the fun tional is not null. This argument is only heuristi and not ompletely justi ed. Indeed, we believe that in su h example as in Figure 1, a term like (15) would indu e a regularizing ee t on  and the support of  would not be a urve any more. In that ase, the integral (15) would not be null. 12

3. Related to the question dis ussed in the last omment is the possibility to ompute a regularized ve tor eld of dire tions for images whi h are onstant ex ept at jump dis ontinuities. A dire t omputation of the ve tor eld outside the hole in an image like Figure 1.a gives a null ve tor eld at all points ex ept the points on the level line separating the bla k from the white region. This may not be a good starting point to extend reasonably the ve tor eld  inside . To initialize the algorithm of steepest des ent, a regularization of  outside may be onstru ted as the ve tor eld of dire tions of the image U0 (t; x) obtained by regularizing u0(x), x 2 B , with the equation u t

ru = div jr uj

u 

=0





in Q = (0; 1)  B (16)

in S = (0; 1)  B

u(0; x) = u0 (x)

for x 2 B: As it is shown in [3℄, this equation permits a regularization of the ve tor eld of dire tions of the gradient of u, i.e., there is a ve tor eld z, jzj  1, su h that u = div(z) and z
ru = jruj. Moreover, for ea h t > 0, div(z (t)) 2 L (B ) if u0 2 L (B ) for all p  1. In this way, we initialize the steepest des ent algorithm des ribed in Se tion 3 with a regularized ve tor eld . This again raises a question, namely, if this ad-ho regularization is really needed or a regularization takes pla e with the algorithm itself, if we use an impli it numeri al s heme to solve (10). t

R

R

B

p

B

p

4. The bound juj  ku0 k 1 ( ) an be repla ed by a onstant depending on ku0 k 1 ( ) . The

onstraint that u = u0 in B ould be relaxed by adding a penalty term like (u u0 )2. Similarly, we ould add a penalty term to onstraint  to be near 0 inside B . In this

ase, we should regularize 0 in B using the equation des ribed in Remark 3. This type of approa h is addressed in the work of Chan and Shen mentioned before [9℄. L

B

L

B

R

B

5. In pra ti e, fun tional (10) is used to interpolate shapes, i.e., to interpolate level sets. The image is de omposed into upper level sets [u0  ℄, whi h are interpolated using (10) to produ e the level sets X u of a fun tion u, whi h is re onstru ted inside by using the 

13

re onstru tion formula

u(x) = supf : x 2 X ug:

To guarantee that the re onstru ted level sets orrespond to the level sets of a fun tion u, they should satisfy that X +1u  X u. In pra ti e, we for e our solution to satisfy this property. In prin iple, our fun tional (10) ould be used dire tly to interpolate fun tions. But, dis ontinuities of the image have a ontribution to the energy whi h is proportional to the jump. This gives dierent weights to dis ontinuities of dierent sizes and, as a onsequen e, they are not treated in the same manner. This is not reasonable if we want to interpolate the shapes of the image, independently of their ontrast. When taking level sets, we treat all shapes equally, and the parameters of the fun tional should only weight geometri quantities (like length, total urvature) and de ide whi h interpolation is taken as a fun tion of them. This approa h is less diusive than dire tly interpolating the gray levels. Theorem 2 in Appendix 1 proves the existen e of minimizers for our model and

an be applied to both ases, binary and gray level images. It guarantees that there are minima of our fun tional. We do not yet know the qualitative properties of those minima. Even if at the intuitive level the main quantities that appear in our fun tional are length and a measure of total urvature (like in the elasti a), we do not have a rigorous proof that this is so. The fun tional was introdu ed on a heuristi basis, but relaxations may when omputed on a o

ur as they o

ur in (16), where div( jrr j ) may represent

at region ([3, 37℄). This requires further study and we shall pursue it elsewhere. 



u u

perimeter area

6. The hoi e made in Remark 5 of de omposing the image u0 into upper level sets, interpolating them and re onstru ting the fun tion u, introdu es a la k of symmetry. Indeed, we are giving more weight to upper level sets than to lower level sets. This an be seen in Figure 3. Figure 3.a displays the image to be interpolated. It is lear that several reasonable solutions are possible and no one of them is preferable to the others. The hoi e we made gives Figure 3.b as solution, favoring that the obje t whose level is 210 goes above the obje t whose level is 0. But, in that ase, the \true" information is la king and we 14

sele ted one of the possible reasonable solution.

Figure 3.a.

Figure 3.b.

2.3 Interpolating gray values along the integral urves of a ve tor eld

Our purpose in this se tion is to further dis uss the term Z

~

(17)

jruj 
ru:

in fun tional (10) (see also [24℄ for a related, L2 and Poisson-equation based, approa h of gray value re onstru tion from image gradients). We shall see that (when  is known), when minimizing (17), we are onstru ting the fun tion u whose values on the boundary are given and whose dire tion of the gradient is given by . We shall dis uss this from a general point of view. Thus, suppose that is an open bounded domain with a Lips hitz boundary and ' 2 L1( ). Let  : ! IR2 be a ve tor eld whose smoothness will be detailed below. We ask the following question: an we interpolate the boundary data ' along the integral urves of  ? In the ase dis ussed in last Se tion,  = ?, and we propagate the boundary data ' along the integral

urves of ?. Heuristi ally,  is orthogonal to the level lines of u. Coming ba k to our general dis ussion, we want to onstru t a fun tion u : ! IR su h that uj = ' and u being onstant along the integral urves of  , i.e., the solutions of the system of ordinary dierential equations 

dX dt

This amounts to say that

=  (X ):

(18)


ru = 0;

(19)

a rst order transport equation whose hara teristi urves are the solution of the system (18). Let us dis uss the diÆ ulties posed by this formulation. First of all, existen e and uniqueness 15

of solutions of (18) is guaranteed when  is a Lips hitz ve tor eld, a very strong regularity assumption, whi h ex ludes any singularity for  . More general existen e results have been obtained in [13℄ via the study of transport equations, indeed, via formulations analog to (19). Typi ally, they are assuming that  is in some Sobolev spa e like 2 W 1 1(IR ), with some other integrability assumptions, and div( ) 2 L1(IR ). These results have been further extended in [12℄, [30℄. In parti ular, P.L. Lions in [30℄ proves a.e. existen e of solutions of (18) for ve tor elds whi h are pie ewise in W 1 1 in a pre ise sense de ned by the author. As observed in these papers, it is not known if the previous result is true for BV ve tor elds. On the other hand, even for a ve tor eld in W 1 1, for whi h we have existen e a.e. of solutions of (18), the problem of onstru ting u satisfying (19) and su h that uj = ' is not obvious. Indeed, onsider a smooth ve tor eld  de ned on a simple domain, like D := fx 2 IR2 : jxj  1g and suppose that the integral urves of  are urves that foliate D and su h that at any point of D we start a urve that ends in another point of D. Then the only possibility to extend u to so that uj = ' in a lassi al sense is that ' takes the same values at the beginning and endpoints of the integral urves of  . A possibility to over ome this diÆ ulty, would be to use the vanishing vis osity method, i.e., to solve the ellipti equation ; lo

N

N

;

;






ru + 
u = 0;

(20)

and let  ! 0. Then, we hope the sequen e u to onverge to some bounded fun tion u whi h solves the problem in a distributional sense. We do not have further information on the regularity of u. On the other hand we do not know in whi h sense the boundary onditions hold. Let us onsider the problem from the algorithmi point of view, i.e., we want to design an ee tive algorithm to solve it. Sin e the problem may be ill-posed, be ause of in ompatibility of boundary data joining two integral urves of the ve tor eld  , we propose a variational formulation of the problem. Let  =  ?. Assume that jj  1. If a solution exists, then  should point in the normal dire tion to the level lines of u. We impli itly assume that  should be

onstru ted as the ve tor eld normal to the level urves of u. Then, formally, 
ru = jruj. Thus, it seems reasonable to minimize the fun tional 

F (u) =

Z



jruj

16

Z




ru;

(exa tly the one introdu ed above) de ned in the set of fun tions of bounded variation BV ( ) whose tra e at the boundary is given by '. Let us formally integrate by parts in the se ond term of F (u) to obtain F (u) =

Z



Z

Z

jruj + div()
u






~nu;

Sin e u;  are known at the boundary, minimizing F amounts to minimize E (u) =

Z

Z



jruj + div()
u:

Let us make pre ise the lass of admissible fun tions where E is minimized. We assume that div() 2 L1 ( ) and ' 2 L1( ). It seems reasonable to impose that the solution u is a bounded fun tion with an L1 bound given by k'k1 (or a onstant related to k'k1 and the size of ). Then the se ond integral in the de nition of E (u) is well de ned. The rst integral requires the use of the spa e of bounded variation fun tions. Thus our admissible lass is A = fu 2 BV ( ) : ju(x)j  k'k1 a.e. uj = 'g. We propose 

Z

Minimize jruj +

u2A

Z



div()
u

(21)

As it is well known ([21℄, [15℄) the solution of this problem has to be understood in a weak sense as the solution of the problem Z

Minimize jruj + u 2 BV ( )

Z



div()
u +

Z





ju 'jdH 1

(22)

juj  k'k1 :

Then we have the following result. 2 L1lo ( )2, with div() 2 L1( ), ' 2 L1( ). u 2 BV ( ) su h that ju(x)j  k'k1 a.e. minimizing (22).

Theorem 1

Proof.

Let 

The result is ontained in ([21℄), Theorem 1.4.

This lari es the role of the term (17) in (10). 17

Then there is a fun tion

3 Numeri al experiments

To minimize (21) we use the steepest des ent method. For that, we formally ompute the Euler-Lagrange equation for u, namely, div



ru + div() = 0 jruj 

(23)

supplemented with Diri hlet boundary onditions for u. In pra ti e, we use the evolution equation ru div() u = div jruj 



t

with Diri hlet boundary data and initial ondition onstru ted as an ad-ho interpolation that will be orre ted by the equation. General existen e results whi h an be adapted to this ase

an be found in ([4℄). Note that the ve tor eld  is assumed to be known in this ase. This limits the usefulness of this model. But we present some experiments below to illustrate the role of this term. To minimize the fun tional (10) we use the steepest des ent method. If we denote the energy term by E (; u), the steepest des ent equations are and

t =

r E (; u) in ~

(24)

ut =

ruE (; u) in ;

(25)

supplemented with the orresponding boundary data and initial onditions. The onstraints on (; u) an be in orporated either by penalization or by brute for e after ea h time step. Given X  IR2 we denote by  the hara teristi fun tion of X , i.e.,  (x) = 1 if x 2 X , otherwise,  (x) = 0. To simplify our notation, let us write g = bjdiv()j , h = a + bjrk  uj. Then X

X

p

X

r E (; u) = pr[hjdiv()jp 2 div()℄ (ru + ru0 B ) = 0

and

rk  u ruE (; u) = div k  g jr k  uj 





 div



ru +  div() = 0; jruj 

(26) (27)

In our experiments, we take k a Gaussian kernel with small varian e, say one or two pixels. In pra ti e, one an also dismiss the kernel k. These equations have to be omplemented with the 18

orresponding boundary onditions for  and u spe i ed by the admissible lass, i.e., we spe ify the normal omponent of  in  ~ and the Diri hlet boundary ondition for u in  , sin e u = u0 in B . The initial onditions are ad-ho interpolations, for instan e, we an take u inside as the average value of u0 in B ,  inside being the dire tion of the gradient of u. One an also take a geodesi propagation inside of the values of u0 in B , with  being again the dire tion of the gradient of u. The exa t numeri al implementation of these ows is given in Appendix 2. In the experiments below, this algorithm is used to interpolate level sets, following the approa h in [31℄, [32℄. The image in B is de omposed into level sets and we get a family of binary images u0 = [  ℄,  = 0; 1; 2; :::; 255. These fun tions are interpolated inside

and we obtain a family of level sets X u. Then the fun tion u is re onstru ted using the re onstru tion formula 

u0





u(x) = supf 2 f0; 1; :::; 255g : x 2 X ug:

As observed in Remark 5 of Se tion 2.2, we for e our solution to satisfy the monotoni ity property of the level sets, i.e., that X +1 u  X u. This is imposed in the initialization of the level set X u and is maintained at ea h iteration of the algorithm by taking the supremum of the urrent solution with the hara teristi fun tion of X u. With this approa h, we diminish the diusive ee ts of the above algorithm and we better apture the shapes and dis ontinuities on the interpolated image. The onstraints on  and kuk1 an be introdu ed after ea h iteration of the above equations. We also omment that the onstraint 
ru = jruj, whi h was introdu ed as a penalization term, ould also be introdu ed by brute for e after ea h time step iteration of the algorithm. 







19

Figure 4.a.

Figure 4.b.

Figure 4. . Let us des ribe the experiments. First, in Figure 4 we display some experiment to illustrate fun tional (21). Figure 4.a displays the full image without the hole. Figure 4.b displays the image with the hole. The ve tor eld  has been omputed on Figure 4.a and we see in Figure 4. the result of interpolating the gray level knowing the ve tor eld inside . We see that the shape of the eye is re overed but not the gray level. This is not a surprise sin e the gray level inside the eye annot be re overed from the gray level on the boundary of . The algorithm is able to apture the shapes inside the eye by integrating the ve tor eld . 20

In the following experiments we show the results of the joint interpolation of gray level and the ve tor eld of dire tions using fun tional (10). The experiments have been done with p = 1 and/or p = 2. The results are quite similar. Unless expli itly stated, we display the results obtained with p = 1. Figure 5.a displays an image made of four ir les overed by a square. In Figure 5.b we display the result of the interpolation. In Figure 6.a, we display an example where the hole is not simply onne ted. The interpolation is displayed in Figure 6.b. Figure 7.a is the image of Lena with two holes, a lower one in the hat and an upper one. Figure 7.b displays the result of the interpolation. Figure 8 displays a zoom of the region around the lower hole. In Figure 9.a we display a level set of u0 orresponding to the region around the lower hole. Figures 9.b and 9. display the orresponding interpolation with p = 1 and p = 2, respe tively. Figure 10.a displays an image with text to be removed. Figure 10.b displays the orresponding re onstru tion result. Figure 11.a displays a portion of an image with text. Figure 11.b displays the orresponding re onstru tion result, obtained with p = 2.

Figure 5.a.

Figure 5.b.

Figure 6.a.

Figure 6.b.

21

Figure 7.a.

Figure 7.b.

Figure 8.a.

Figure 8.b.

Figure 9.a.

Figure 9.b.

Figure 9. .

22

Figure 10.a.

Figure 10.b.

Figure 11.a.

Figure 11.b. 23

4 Con luding remarks

In this paper we have proposed a formal variational approa h for lling-in regions of missing data in still images. The basi idea is to smoothly extend inside the hole both the ve tor eld obtained from the image gradient and the orresponding gray values. We have presented a number of examples and showed theoreti al results regarding the proposed formulation. A number of resear h dire tions are suggested by the work here presented. First of all, we need to omplement this algorithm by a te hnique apable of lling-in textured regions. Se ondly, the extension of the framework to the lling-in of other type of missing imagery data is of great interest for a number of appli ations. Last, we would like to study these ideas for interpolation in video data. These topi s will be the subje t of subsequent reports. Appendix 1: Existen e of minimizers

Re all that ~ = [ B is an open bounded set whose boundary is Lips hitz. For simpli ity, let us de ne the lass B of admissible pairs (; u) where  2 W 1 (div; ~ ), u 2 BV ( ~ ; B; u0 ), jj  1, 
ru = jruj, juj  ku0k 1 ( ) and 
nj ~ = 0
nj ~ . ;p

L

Theorem 2

B





If p > 1, there is a minimum (; u) 2 B for the problem (10).

Let us denote by E (; u) the energy de ned in (10). Let ( sequen e for E (; u). Sin e a > 0, we have that

Proof.

Z

~

n

; un )

be a minimizing

jdiv(n )jp

is bounded. Sin e j j  1, we have that  is weakly relatively ompa t in all spa es L ( ~ )2 for all 1  q < 1 and we may assume that  !  weakly in L ( ~ )2 for all 1  q < 1 and in W 1 (div; ~ ). Now, using the onstraint n

q

n

q

n

;p

Z

and integrating by parts the term Z

~

~

R

jrunj =

~ n
run ,

Z

~

n
run ;

we obtain Z

jrunj =

Z

=

Z

div(n )un + ~ [n; n℄un 

Z div(n )un + ~ [0 ; n℄u0 ~

~





24



The integration by parts is possible by results of Anzellotti ([2℄) given above. From the above identity, we obtain jru j  kdiv( )k ku k 0 + ~ ju0 j

~

sin e u = u0 in  ~ , where p0 is the exponent onjugated to p. Then, modulo a subsequen e,we may assume that u onverges to some fun tion u in L1( ~ ). Note that u 2 BV ( ~ ; B; u0). Sin e we have an L1 bound on u , we also have that u onverges to u in L ( ~ ) for all 1  q < 1. Then rk  u ! rk  u uniformly in ~ . In parti ular, we obtain Z

Z

n

n

p

n p



n

n

n

q

n

n

Z

Z

~

jdiv(n )jp (a + bjrk  unj)dx; jdiv()j (a + bjrk  uj)dx  liminf n ~ p



i.e,

E (; u)  liminf E (n ; un ): n

To prove that (; u) is a minimizer of E (; u) it suÆ es to prove that 
ru = jruj. For that, sin e div( ) weakly onverges to div() in L ( ~ ) and u ! u in L 0 ( ~ ), passing to the limit in p

n

Z

~

we get that

R

n
run =

Z

~

p

n

Z

div(n )un + ~ [0 ; n℄u0 

onverges to

~ n
run

Z

Z

~

Z

div()u + ~ [0 ; n℄u0 = ~ 
ru:



Now, Z

Z

Z


run jrunj = liminf jruj  liminf n n ~ n ~ ~



=

Z

~



Z


ru  ~ jruj;



i.e., 
ru = jruj. The pair (; u) is a minimum of E in the lass B of admissible fun tions for E. Appendix 2: Numeri al implementation

To solve equations (24) and (25), we use a impli it dis retization in time. To be pre ise, we write r E (; 0; u; v) = pr[h( + jdiv(0 )j 2)div()℄ (ru + ru0 ) = 0 (28) 

p

25

B

and 

ruE (; 0; u; v) = div k  g p

rk  u  + jrk  vj2

!

!



 div

ru +  div() = 0: (29)  + jrvj2

p

Then, we use the dis retization in time given by and

n+1 n =
tr E (n+1 ; n ; un ; un );

(30)

un+1 un =
tru E (n+1 ; n+1; un+1 ; un ):

(31)

Finally, we make the hange of variables 

+1 = n+1

n

n, vn+1 = un+1 un

and we have

 n+1 =
tr E ( n+1 + n; n; un ; un );

(32)

vn+1 =
truE (n+1 ; n+1; vn+1 + un ; un ):

(33)

Now, sin e 
nj ~ =  +1
nj ~ and u j = u +1j , then the normal omponent of  +1 and the value of v +1 are zero on the boundary, and we may use a onjugate gradient method to solve (32) and (33). The onstraint jj  1 is in orporated by brute for e after ea h time step. We an also set  = 0 and in orporate the onstraint that jruj = 
ru by brute for e after ea h time step. n



n



n



n



n

n

A knowledgments

We thank S. Betelu, A. Bertozzi, T. Chan, C. Kenney, P.L. Lions, J.M. Morel, S. Osher, E. Simon elli, and J. Shen for interesting onversations on image inpainting and lling-in. This work was partially supported by a grant from the OÆ e of Naval Resear h ONR-N00014-97-1-0509, the OÆ e of Naval Resear h Young Investigator Award, the Presidential Early Career Awards for S ientists and Engineers (PECASE), a National S ien e Foundation CAREER Award, by the National S ien e Foundation Learning and Intelligent Systems Program (LIS), and the TMR European proje t \Vis osity Solutions and their appli ations," referen e FMRX-CT98-0234.

26

Referen es

[1℄ L. Ambrosio, N. Fus o and D. Pallara, tinuity Problems, Forth oming book.

Fun tions of Bounded Variation and Free Dis on-

[2℄ G. Anzellotti, \Pairings between measures and bounded fun tions and ompensated Compa tness," Ann. di Mat. Pura ed Appl. 135, 293-318, 1993. [3℄ F. Andreu, C. Ballester, V. Caselles, J. M. Mazon, \Minimizing Total Variation ow," Dierential and Integral Equations, to appear. [4℄ F. Andreu, C Ballester, V. Caselles and J.M. Mazon, \The Diri hlet problem for the Total Variation ow," preprint. [5℄ G. Bellettini, G. Dal Maso and M. Paolini, \Semi ontinuity and relaxation properties of a

urvature depending fun tional in 2D," Ann. S uola Normale Sup. di Pisa, Cl. S i. 20(4), pp. 247-297, 1993. [6℄ M. Bertalmio, G. Sapiro, C. Ballester and V. Caselles, \Image inpainting," University of Minnesota IMA TR, De ember 1999 (available at www.ima.umn.edu). [7℄ M. Bertalmio, G. Sapiro, C. Ballester and V. Caselles, \Image inpainting," to appear in Computer Graphi s, SIGGRAPH 2000, July 2000. [8℄ C. Braverman. Photoshop retou hing handbook. IDG Books Worldwide, 1998. [9℄ T. Chan and J. Shen, \Mathemati al models for lo al deterministi inpaintings," UCLA CAM Report 00-11, Mar h 2000 (available at www.math.u la.edu). [10℄ R. Dautray and J.L. Lions, Analyse Mathematique et Cal ul Numerique pour les S ien es et les Te hniques, Tome 2, Masson, 1985.  [11℄ F. Demengel, Introdu tion aux Equations aux Derivees Partielles Elliptiques, Diderot Ed. 1999.

[12℄ B. Desjardins, \A few remarks on ordinary dierential equations," Comm. in Partial Di. Equations 21, pp. 1667-1703, 1996. 27

[13℄ R. J. DiPerna and P. L. Lions, \Ordinary dierential equations, transport theory and Sobolev spa es," Inventiones Math. 98, pp. 511-547, 1989. [14℄ A. Efros and T. Leung, \Texture synthesis by non-parametri sampling," Pro . IEEE International Conferen e Computer Vision, pp. 1033-1038, Corfu, Gree e, September 1999. [15℄ I. Ekeland and R. Temam, Amsterdam, 1976.

Convex Analysis and Variational Problems,

North Holland,

[16℄ G. Emile-Male. The Restorer's Handbook of Easel Painting. Van Nostrand Reinhold, New York, 1976. [17℄ L.C. Evans and R.F. Gariepy, Measure Theory and Fine Properties of Fun tions, Studies in Advan ed Math., CRC Press, 1992. [18℄ G. P. Galdi, An Introdu tion to the Mathemati al Theory of the Navier-Stokes Equations, Volume I, Springer Verlag, 1994. [19℄ D. Heeger and J. Bergen. \Pyramid based texture analysis/synthesis," Computer Graphi s, pp. 229-238, SIGGRAPH 95, 1995. [20℄ P. G. de Gennes and J. Prost, The Physi s of Liquid Crystals, Oxford S ien e Publi ations, 1998. [21℄ M. Giaquinta, G. Modi a and J. Sou ek, \Fun tionals with linear growth in the al ulus of variations I," Comment. Math. Univ. Carolina 20, pp. 143-156, 1979. [22℄ A. Hirani and T. Totsuka. \Combining frequen y and spatial domain information for fast intera tive image noise removal," Computer Graphi s, pp. 269-276, SIGGRAPH 96, 1996. [23℄ L. Joyeux, O. Buisson, B. Besserer, S. Boukir. \Dete tion and removal of line s rat hes in motion pi ture lms," Pro eedings of CVPR'99, IEEE Int. Conf. on Computer Vision and Pattern Re ognition, Fort Collins, Colorado, USA, June 1999. [24℄ C. Kenney and J. Langan, \A new image pro essing primitive: Re onstru ting images from modi ed ow elds," preprint, University of California, Santa Barbara. 28

[25℄ D. King. The Commissar Vanishes. Henry Holt and Company, 1997. [26℄ R. Kohn and R. Temam, \Dual spa es of stresses and strains, with appli ations to Hen ky plasti ity," Appl. Math. Optimization 10 (1983), 1-35. [27℄ A.C. Kokaram, R.D. Morris, W.J. Fitzgerald, P.J.W. Rayner. \Dete tion of missing data in image sequen es," IEEE Transa tions on Image Pro essing 11(4), pp. 1496-1508, 1995. [28℄ A.C. Kokaram, R.D. Morris, W.J. Fitzgerald, P.J.W. Rayner. \Interpolation of missing data in image sequen es," IEEE Transa tions on Image Pro essing 11(4), pp. 1509-1519, 1995. [29℄ A. Li hnewski and R. Temam, \Pseudosolutions of the time dependent minimal surfa e problem," Journal of Dierential Equations 30, pp. 340-364, 1978. [30℄ P.L. Lions, \Sur les equations dierentielles ordinaires et les equations de transport," Comptes Rendus A ad. S ien es 326, pp. 833-838, 1998. [31℄ S. Masnou, Filtrage et Deso

lusion d'Images par Methodes d'Ensembles de Niveau, These, Universite Paris-Dauphine, 1998. [32℄ S. Masnou and J.M. Morel, \Level lines based diso

lusion," 5th IEEE International Conferen e on Image Pro essing, Chi ago, Illinois, O tober 4-7, 1998. [33℄ M. Nitzberg, D. Mumford, and T. Shiota, Filtering, Segmentation, and Depth, SpringerVerlag, Berlin, 1993. [34℄ D. L. Ringa h and R. Shapley, \Spatial and temporal properties of illusory ontours and amodal ompletion", Vision Resear h 36, pp. 3037-3050, 1996. [35℄ L. Rudin, S. Osher and E. Fatemi. \Nonlinear total variation based noise removal algorithms," Physi a D 60, pp. 259-268, 1992. [36℄ E. Simon elli and J. Portilla. \Texture hara terization via joint statisti s of wavelet oeÆ ient magnitudes," 5th IEEE Int'l Conf. on Image Pro essing, Chi ago, IL. O t 4-7, 1998. 29

[37℄ D. Strong and T. Chan, \Exa t solutions to Total Variation regularization problems," UCLA CAM Report 96-41, O tober 1996. [38℄ R. Temam, \On the ontinuity of the tra e of ve tor fun tions with bounded deformation," Appl. Analysis 11, pp. 291-302, 1993. [39℄ S. Walden. The Ravished Image. St. Martin's Press, New York, 1985. [40℄ W. P. Ziemer, Weakly Dierentiable Fun tions, GTM 120, Springer Verlag, 1989.

30