SIAM J.I MAGING S CIENCES c2013Society for Industrial and Applied Mathematics Vol.6,No.3,pp.1345–1366A Variational Approach for Image Stitching II:Using Image Gradients ∗Michael K.Ng †and Wei Wang ‡Abstract.In [W.Wang and M.K.Ng,SIAM J.Imaging Sci.,6(2013),pp.1318–1344],we proposed and developed an image stitching algorithm by studying a variational model for automatically computing weighting mask functions on input images and stitching them together.The main aim of this paper is to further develop an image stitching algorithm using the gradients of input images.Our idea is to study a variational method for computing a stitched image by using an energy functional containing the data-fitting term based on the difference between the gradients of the stitched image and the input images,and the Laplacian regularization term based on the smoothness of weighting mask functions.The use of image gradient information allows us to automatically adjust the stitched image to handle color inconsistency across input images.In the model,we incorporate both boundary conditions of the stitched image and the weighting mask functions.The existence of a solution of the proposed energy functional is shown.We also present an alternating minimizing algorithm for solvingthe variational model numerically,and we show the convergence of this algorithm.Experimental results show that the performance of the proposed method is better than the other testing methods proposed in the literature for input images with color inconsistency.Key words.image stitching,variational model,image gradients,algorithm,weighting mask functions AMS subject classifications.65K10,68U10,91-08DOI.10.1137/rge objects often cannot be captured in a single picture under the view of camera phones or Planetary Data System cameras.Image stitching is a process that combines two or more images and blends them into one.It is the main step in the generation of panoramic images,and it is widely used in remote sensing [1,2],superresolution [3],andtexture synthesis [4].The main aim of image stitching is to find a visually acceptable orseamless blending from the input images with overlapping regions.Stitching problems usually contains two steps:image alignment and image blending.The goal of image alignment is to find corresponding point pairs in the overlapping region of two images;see,for instance,[5,6,7,8,9].Image blending combines the two aligned images seamlessly.There are two main approaches for image stitching in the literature.Optimal seam meth-ods [4,10,11,12,13,14]search for a curve in the overlapping region on which the differencesbetween two original images are minimal.Then each image is copied to the corresponding side of the curve.When the difference between two input images is zero on the curve,there will∗Received by the editors April 2,2012;accepted for publication April 10,2013;published electronically July 11,2013./journals/siims/6-3/87214.html †Centre for Mathematical Imaging and Vision and Department of Mathematics,Hong Kong Baptist University,Kowloon Tong,Hong Kong (mng@.hk ).This author’s research was supported in part by an HKRGC grant and HKBU FRG grant.‡Corresponding author.Department of Mathematics,Tongji University,Shanghai 200442,China (weiwamng@ ).This author’s research was supported by National Natural Science Foundation of China grant 11201341and China Postdoctoral Science Foundation grants 2012M511126and 2013T60459.1345D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p1346MICHAEL K.NG AND WEI WANGnot be visual seams on the curve.However,the seam is visible when there is no such curve in practice.This case happens when there is a global intensity difference between the input images.The second approach uses a weighting mask over the overlapping area to smooth thetransition between the input images.In feathering [15],the stitching image is a weighted com-bination of the input images.Other image blending methods based on combined features in different resolutions are introduced in [16,17,18].The weighting mask functions vary spatiallyas a function of the distance from the seam.Pyramid blending [19]smooths the transition in each frequency band independently.The key step of this kind of approach is how to choose a weighting mask function properly.Similarly,Zomet et al.[20]and Levin et al.[21]proposed an energy functional based on image intensity and gradient for computing image stitching.In [22,23],Suen,Lam,and Wong considered optimization problems for image stitching based on image intensity,gradient,and curvature.Thebasic idea is to find a stitching image whose pixel intensity,gradient,and curvature values are close to those of the inputs.The methodcan stitch seamlessly when uniform photometric inconsistency exists.However,the weightingmask functions associated with these energy functionals are still required to choose properly in order to obtain high visual quality of image stitching ually,suitable weighting mask functions are unknown in advance.In [24],we have proposed and developed an image stitching algorithm by studying avariational method for automatically computing weighting mask functions on input images and stitching them together.In order to enhance visual quality of stitched images where inputimages have color inconsistency or different lighting conditions,a color correction procedure is required in their postprocessing step.The main aim of this paper is to further develop an imagestitching algorithm using the gradients of inputimages.Our idea is to study a variational model to compute a stitched image by using an energy functional containing the data-fitting term based on the difference between the gradients of the stitched image and the input images,and the Laplacian regularization term based onthe smoothness of weighting mask functions.The use of image gradient information allows us to automatically adjust the stitched imageto handle color inconsistency across input images.In the model,we incorporate boundary conditions of both the stitched image and the weighting mask functions.The existence of the solution of the proposed energy functional is shown.We also present an alternating minimizing algorithm to solve the variational model numerically,and we show the convergence of this algorithm.Experimental results show that the performance of the proposed method is betterthan that of the other testing methods proposed in the literature when input images have color inconsistency or different lighting conditions.This paper is organized as follows.In section 2,we describe the proposed model and show the existence of the solution of the proposed energy functional.In section 3,we present an alternating minimizing algorithm to solve the proposed model numerically,and we show theconvergence of this algorithm.In section 4,we present numerical examples which demonstrate the effectiveness of the proposed model.Finally,the concluding remarks are given in section 5.2.The proposed model.2.1.The motivation.Let us consider the following example to demonstrate the proposed model.In Figures 1(a)–1(b),there are two images to be stitched.Assume that the input images have been aligned,i.e.,we know exactly where the overlapping region is in each inputD o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h pVARIATIONAL IMAGE STITCHING USING GRADIENTS 13471Ω0N (a)(b)(c)Figure 1.The two input images (a)I 1and (b)I 2;(c)the overlapping region Ω.image.The overlapping region Ωis shown in Figure 1(c).In [24],we study a variational approach containing an energy functional based on image intensity to determine both weighting mask functions and a stitched image.Let I 1and I 2represent the two input images.We want to determine a weighting mask function in theoverlapping region such that it is smooth based on the goal of perfect blending and also variesregularly based on our assumption of transition direction from one image to another image.The following optimization problem is studied in [24]:(2.1)min(α,u )∈Λold E old (α,u )= Ω|∇α|2dx + Ωα(u −I 1)2dx +Ω(1−α)(u −I 2)2dx ,where Λold ={(α,u )|α∈W 1,2(Ω),α∈˜Λ,0≤α≤1,u ∈L 2(Ω)}is our admissible set containing the constraints of two variables:the weighting mask functionαand the stitched image u .Here we define the gradient operator to be ∇:= ∂x ∂y .We note that the first term in (2.1)ensures thesmoothness of α,and the next two terms are interpreted as fidelity terms.When αtends to 1,u is forced to get closer to I 1.When αtendsto 0,u is forced to get closer to I 2.On the other hand,the set ˜Λfor αcoincides with the transition direction which is based on the overlapping situation of the input images in practice.As an example of the overlapping region Ωin Figure 1(c),the weighting mask function varies smoothly from 1to 0along the left-to-right direction.Then we consider Dirichlet boundary condition α=1on the vertical line marked with blue color and α=0on the vertical line marked with red color,and the Neumann boundary condition ∂α∂ n =0on the twohorizontal lines marked with green color.In particular,the set ˜Λin this example is given as follows:˜Λ= α|∂x α≤0,α=f on D,∂α∂ n =0on N ,and D ∪N =∂Ω ,where D and N refer to the boundaries for the Dirichlet and Neumann boundary conditions in ∂Ω,respectively.The situations where we want to smooth the transition from right to left D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p1348MICHAEL K.NG AND WEI WANGFigure2.The stitched image by using(left)the old model in(2.1)and(right)the proposed method.(∂xα≥0),top to bottom(∂yα≤0),or bottom to top(∂yα≥0)can be handled similarly;see[24]for a detailed discussion.In Figure2(left),we give the stitched image by using themodel in(2.1).Although wefind that both input images can be stitched smoothly together,the visual quality of the stitched image is not good where the two input images have differentlighting conditions;see Figures1(a)–1(b).In this paper,our main aim is to further develop an image stitching algorithm using the gradients of input images so that the stitched image can be adjusted automatically to handlecolor inconsistency across input images;see the stitched image computed by the new modelin Figure2(right).2.2.The gradient stitching.LetΩ1(Ω2,respectively)be the region viewed exclusively inimage I1(I2,respectively)satisfyingΩ1∩Ω2=Ω1∩Ω=Ω2∩Ω=∅,and letΩc=Ω1∪Ω2∪Ωbe the combined domain(see Figure3for the example in Figure1).First,we define a newweighting mask function in the combined domainΩc based on the weighting maskαgiven inthe overlapping regionΩas follows:(2.2)ω(x,y):=⎧⎨⎩1,(x,y)∈Ω1,α(x,y),(x,y)∈Ω,0,(x,y)∈Ω2.Then we extend I1and I2into the combined domainΩc by setting proper values inΩ2andΩ1,respectively.We consider the following minimization problem:(2.3)min(α,u)∈ΛE(α,u)=Ω|∇α|2dx+Ωcω|∇u−∇I1|2dx+Ωc(1−ω)|∇u−∇I2|2dx, whereΛ=(α,u)|α∈W1,2(Ω),α∈˜Λ,0≤α≤1,u∈W1,2(Ωc),u=I1on∂Ω1∩∂Ωc,u=I2on∂Ω2∩∂Ωc,∂u∂ n=0on N,I min≤u≤I max(2.4)is our admissible set containing the constraints.Here I min and I max are the minimum andmaximum values,respectively,in I1and I2.We remark that the gradients of the stitchedimage and the weighting mask functions are involved in(2.3);thus their boundary conditionsmust be included in(2.4)in order to guarantee the existence of the solution of(2.3).Downloaded1/31/13to218.21.115.245.RedistributionsubjecttoSIAMlicenseorcopyright;seehttp://www.siam.org/journals/ojsa.phpVARIATIONAL IMAGE STITCHING USING GRADIENTS 1349Ω1ΩΩ2Figure 3.The combined domain Ωc for the example in Figure 1.In the previous subsection,we have discussed how to set the boundary condition of α.The boundary condition for u can be defined similarly.As an example,the new weighting mask function ωin Figure 3varies smoothly from 1to 0along the left-to-right direction for thetwo input images in Figures 1(a)–1(b).Moreover,we have the following Dirichlet boundary condition:the values of the stitched image u are equal to the values of the input image I 1on the line marked with blue color,and the values of the stitched image u are equal to the values of the input image I 2on the line marked with red color.The two horizontal lines marked with green color are set for the Neumann boundary condition ∂u ∂ n =0.According to this setting,we keep the colors of I 1and I 2in the left-and right-hand sides of the stitchedimage,respectively,and fuse their colors smoothly in the overlapping region of the stitched image.We see from Figure 2(c)that the new model can provide better visual quality of the stitched image than that of the previous model.In section 4,we will show more examples ofhow to set the boundary condition of u to stitch input images together.2.3.Mathematical analysis.In this subsection,we analyze the existence of a solution for (2.3).Let us consider the unconstrained form of the problem (2.3)as follows:(2.5)min (α,u,v )∈Λ0{E1(α,u,v )},whereE 1(α,u,v )≡ Ω|∇α|2dx+Ωc ω|∇u −∇I 1|2dx+ Ωc(1−ω)|∇u −∇I 2|2dx+μ2 Ω|∂x α+v |2dx ,Λ0=(α,u )|α∈W 1,2(Ω),0≤α≤1,α=f on D,∂α∂ n =0on N ,u ∈W 1,2(Ωc ),u =I 1on ∂Ω1∩∂Ωc ,u =I 2on ∂Ω2∩∂Ωc ,∂u∂ n =0on N ,I min ≤u ≤I max ,v ≥0.We introduce a new variable v and add a penalty term to handle the constraint ∂x α≤0.If μis large enough,v ≥0guarantees ∂x α≤0.Then we have the following theorem for problem (2.5).Theorem 2.1.Assume f is continuous on ∂Ω,and 0≤f ≤1.Assume I 1,I 2∈W 1,2(Ωc ),and I 1and I 2are continuous on ∂Ω1∩∂Ωc and ∂Ω2∩∂Ωc ,respectively.Then the problem in (2.5)has at least one solution.D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p1350MICHAEL K.NG AND WEI WANGProof .It is well known that,with the requirements that the length of D be nonzero and f be continuous,the solution of the problem min α=f on D,∂α∂ n=0on NΩ|∇α|2dx exists and can be derived by solving the following equation:⎧⎨⎩Δα=0in Ω,α=fon D,∂α∂ n =0on N .The above statement can be found in Chapter 2of [25].Then we denote the solution as ˆα,and by the maximum principle,we can easily obtain 0≤ˆα≤1.We can also get ˆu satisfying the constraints in Λ0by using the same arguments.Then we choose α=ˆα,u =ˆu ,and v =1,and our energy E 1(α,u,v )will be finite.This implies that the infimum of the energy E 1must be finite.Suppose {(αn ,˜u n ,v n )}⊂Λ0is a minimizing sequence of problem (2.5).Then there existsa constant M >0such that E 1(αn ,˜u n ,v n )≤M,which implies Ω|∇αn |2dx ≤M and Ω|∂x αn +v n |2dx ≤M.Associated with the boundedness of α,we derive that {αn }is uniformly bounded in W 1,2(Ω),and {v n }is uniformly bounded in L 2(Ω).Noting that W 1,2(Ω)is compactly embedded in L 2(Ω),we deduce that up to a subsequence αn converges to α∗∈W 1,2(Ω),and v n weakly converges to v ∗∈L 2(Ω);i.e.,(2.6)αn −−−−→L 2(Ω)α∗,and v n v ∗∈L 2(Ω).Now we can derive the corresponding ωn based on αn ,and we can easily deduce ωn −−−−→L 2(Ωc )ω∗,where ω∗=⎧⎨⎩1,x ∈Ω1,α∗,x ∈Ω,0,x ∈Ω2.Then we fix αn and v n to consider the following problem:min (αn ,u,v n )∈Λ0{E 1(αn ,u,v n )}.We denote the solution of the above problem as u n .Thus,E 1(αn ,u n ,v n )≤E 1(αn ,˜u n ,v n ).This means {(αn ,u n ,v n )}is also a minimizing sequence of problem (2.5).From Theorem 2in [21],we know that the minimum of E 1(αn ,u,v n)is obtained with the image having the closest D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h pVARIATIONAL IMAGE STITCHING USING GRADIENTS 1351derivatives to ωn ∇I 1+(1−ωn )∇I 2under the L 2norm.Noting that I 1,I 2∈W 1,2(Ωc ),we have that Ω|∇u n |2dx is uniformly bounded,and I min ≤u n ≤I max ;i.e.,{u n }is uniformly bounded in W 1,2(Ωc ).We deduce that up to a subsequence u n converges to I min ≤u ∗∈W 1,2(Ωc )≤I max ,and ∇u n weakly converges to∇u ∗∈L 2(Ωc );i.e.,(2.7)u n −−−−→L 2(Ωc )u ∗,and ∇u n ∇u ∗∈L 2(Ωc ).The last convergence result can be derived in the same way as ∇αn converges;see the following discussion.Then we can rewrite the first three terms of E 1(αn ,u n ,v n ),associated with (2.6),(2.7),and from Theorem 2.1.3in [26],the dominated convergence theorem,and the convexity of the three terms in ∇α,∇u ,we have the following lower semicontinuity:lim infn →∞ Ω|∇αn |2dx + Ωc ωn |∇u n −∇I 1|2dx + Ωc(1−ωn )|∇u n −∇I 2|2dx =lim infn →∞ Ω|∇αn |2dx + Ωc ωn |∇I 1|2dx + Ωc(1−ωn )|∇I 2|2dx +Ωc |∇u n |2dx −2 Ωc ωn ∇u n ·∇I 1dx −2 Ωc(1−ωn )∇u n ·∇I 2dx ≥Ω|∇α∗|2dx + Ωc ω∗|∇I 1|2dx + Ωc (1−ω∗)|∇I 2|2dx)+Ωc |∇u ∗|2dx −2 Ωc ω∗∇u ∗·∇I 1dx −2 Ωc(1−ω∗)∇u ∗·∇I 2dx =Ω|∇α∗|2dx + Ωcω∗|∇u ∗−∇I 1|2dx +Ωc (1−ω∗)|∇u ∗−∇I 2|2dx .By using (2.6),we know that αn →α∗,a .e .in Ω.Noting that 0≤αn ≤1,we deduce0≤α∗≤1,a .e .in Ω.Then we claim thatα∗=f on D,∂α∗∂ n =0on N .Indeed,αn is bounded in BV (Ω).Then up to a subsequence we have ∇αn ∗−→∇α∗in a measure sense;i.e.,for any ϕ∈C ∞0(Ω)2,the following convergence holds:(2.8)Ωϕ·∇αn dx →Ωϕ·∇α∗dx .Noting that ∇αn is bounded in L 2(Ω),we havethat up to a subsequence it weakly converges in L 2(Ω).By combining (2.8)and the property that a function in L 2(Ω)can be approximatedby smooth functions,we obtain the following result:∇αn ∇α∗in L 2(Ω).Then we choose any φ∈C ∞(¯Ω)2,where ¯Ωmeans the closure of Ω,and we have Ωφ·∇αn dx →Ωφ·∇α∗dx .D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p1352MICHAEL K.NG AND WEI WANGBy using the formula of integration by parts, ∂Ωαn φ· n d H 1− Ωαn div φdx →∂Ωα∗φ· n d H 1−Ωα∗div φdx ,where n is the normal vector of ∂Ω,and H 1is the one-dimensional Hausdorffmeasure.Noting (2.6),we have Ωαn div φdx → Ωα∗div φdx .Then combining the boundary condition,we deriveD fφ· n d H 1= D α∗φ· n d H 1+ Nα∗φ· n d H 1.Because this is valid for any φ,and f is continuous,we obtain that α∗=f on D and ∂α∗∂ n =0on N .We remark here that we can get u ∗|∂Ω1∩∂Ωc =I 1,u ∗|∂Ω2∩∂Ω=I 2,∂u ∗∂ n =0on N byusing similar arguments.As a consequence of the lower semicontinuity for the L 2norm,lim infn →∞ Ω|∂x αn +v n |2dx ≥ Ω|∂x α∗+v ∗|2dx .Therefore,min (α,u,v )∈Λ0{E 1(α,u,v )}=lim inf n →∞E 1(αn ,u n ,v n )≥E 1(α∗,u ∗,v ∗),and we have α∗=f on D ,∂α∗∂ n =0on N ,0≤α∗≤1;u ∗|∂Ω1∩∂Ωc =I 1,u ∗|∂Ω2∩∂Ωc =I 2,∂u ∗∂ n =0on N ,I min ≤u ∗∈W1,2(Ωc )≤I max ;and v ∗≥0.This completes the proof.3.The alternating minimization algorithm.In this section,we develop a numerical method for solving problem (2.3).We first transform problem (2.3)into the following equiv-alent problem:(3.1)min α∈W 1,2(Ω),u ∈W 1,2(Ωc ),v ∈L 2(Ω){E 2(α,u,v )}subject to ∂x α+v =0,v ≥0,0≤α≤1,α=f on D,∂α∂ n=0on N ,u |∂Ω1∩∂Ωc =I 1,u |∂Ω2∩∂Ωc =I 2,∂u∂ n =0on N ,I min ≤u ≤I max ,whereE 2(α,u,v )≡ Ω|∇α|2dx + Ωc ω|∇u −∇I 1|2dx + Ωc (1−ω)|∇u −∇I 2|2dx .There are three variables α,u ,and v in the optimization problem.A simple but effectivestrategy is to use an alternating minimization scheme for solving the optimization problem as follows.D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h pVARIATIONAL IMAGE STITCHING USING GRADIENTS 1353Algorithm 1.(i)Let α0=α0be the initial mask function.(ii)At the k th iteration:•Given αk ,generate ωk based on (2.2),and compute u k +1by solving (3.2)min uE 3(u )≡Ωcωk |∇u −∇I 1|2dx + Ωc(1−ωk )|∇u −∇I 2|2dx subject tou |∂Ω1∩∂Ωc =I 1,u |∂Ω2∩∂Ωc =I 2,∂u ∂ n =0on N ,I min ≤u ≤I max .•Given u k +1,compute αk +1,vk +1by solving min α∈W 1,2(Ω),v ∈L 2(Ω)E 4(α,v )(3.3)≡Ω|∇α|2dx + Ωα(|∇u k +1−∇I 1|2−|∇u k +1−∇I 2|2)dx ,subject to∂x α+v =0,v ≥0,0≤α≤1,α=f on D,∂α∂ n =0on N .(iii)Go back to step (ii)until ||αk +1−αk ||||αk +1||≤ αand ||u k +1−u k ||||u k +1||≤ u .3.1.Solution of u -subproblem.For the minimization subproblem in (3.2),we first in-troduce notation ι1(w ):=0,I min ≤w ≤I max ,+∞otherwise ,H 1(u )=Ωc ωk |∇u −∇I 1|2dx + Ωc (1−ωk )|∇u −∇I 2|2dx .Then we can rewrite problem (3.2)as follows:(3.4)min u,w{E 5(u,w )≡ι1(w )+H 1(u )}subject to u −w =0,u |∂Ω1∩∂Ωc =I 1,u |∂Ω2∩∂Ωc =I 2,∂u ∂ n=0on N .For this constrained optimization problem,we can employ the alternating direction methodof multipliers (ADMM),which is introduced in [27,28,29]to solve this problem.By using the Lagrangian multiplier λ1to the linear constraint,the augmented Lagrangian function of (3.4)is given by L 1(u,w,λ1)=ι1(w )+H 1(u )+ λ1,u −w +μ12||u −w ||2.D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p1354MICHAEL K.NG AND WEI WANGHere the scalar product ·,· is the inner product of L 2(Ωc ).Then the ADMM iterations are described in the following steps:(3.5)w i +1=argmin wL 1(u i ,w,λi1),(3.6)u i +1=argmin u |∂Ω1∩∂Ωc =I 1,u |∂Ω2∩∂Ωc =I 2,∂u∂ n =0on N L 1(u,w i +1,λi1),(3.7)λi +11=λi 1+μ1(u i +1−w i +1).We rewrite problem (3.5)as follows:w i +1=argmin w ι1(w )+μ12w −u i −λi 1μ1 2 .It is obvious that w i +1is equal to the following projection:w i +1=max min u i +λi 1μ1,I max ,I min .The Euler–Lagrange equation for problem (3.6)is given by the following equation:(∂T x ∂x +∂T y ∂y +μ1I )u =∂T x ωk ∂x I 1+∂T y ωk ∂y I 1+∂T x (1−ωk )∂x I 2+∂Ty (1−ωk )∂y I 2+μ1w i +1−λi 1.Here we note that ∂x and ∂y are the matrix forms of the corresponding discrete convolution operators,and I is the identity matrix.As the coefficient matrix of the above linear system isgiven by the second-order differential operator with Dirichlet boundary conditions,this matrix equation can be solved by many efficient numerical solvers and preconditioning techniques;see,for instance,the book by Saad [30].3.2.Solution of α-subproblem.For the minimization subproblem in (3.3),we make use of the variable-splitting technique to solve the constrained problem.Similarly,let’s firstintroduce some notation.For simplicity,we denote g =|∇u k +1−∇I 1|2−|∇uk +1−∇I 2|2and denote ι2,ι3as the indicator functions for the sets 0≤α≤1,v ≥0,respectively.Moreprecisely,they are defined as follows:ι2(α):=0,0≤α≤1,+∞otherwiseand ι3(v ):=0,v ≥0,+∞otherwise .D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h pOn the other hand,we denoteH 2(α)= Ω|∇α|2dx + Ωαg dx;then we can rewrite problem (3.3)as follows:(3.8)min α,β,v{E 6(α,β,v )≡ι1(β)+ι2(v )+H 2(α)},subject to ∂x α+v =0,α−β=0,α=f on D,∂α∂ n =0on N .For this constrained optimization problem,we can also employ ADMM to solve this problem.By using the Lagrangian multiplier λ2,λ3to the linear constraint,the augmented Lagrangianfunction of (3.8)is given byL 2(α,β,v,λ2,λ3)=ι2(β)+ι3(v )+H 2(α)+ λ2,α−β + λ3,∂x α+v +μ22||∂x α+v ||2+μ22||α−β||2.We penalize the last two norms equally usingthe same parameter μ2.Then the ADMM iterations are described in the following steps:(3.9)(βi +1,v i +1)=argmin β,v L 2(αi ,β,v,λi 2,λi3),(3.10)αi +1=argmin α=f on D,∂α∂ n=0on N L 2(α,βi +1,v i +1,λi 2,λi3),(3.11)λi +12=λi2+μ2(αi +1−βi +1),(3.12)λi +13=λi 3+μ2(∂x αi +1+v i +1).We rewrite problem (3.9)in detail as follows:(βi +1,v i +1)=argmin β,vι2(β)+μ22 β−αi −λi 2μ2 2+ι3(v )+μ22 v +∂x αi +λi 3μ2 2 .It is obvious that βi +1and v i +1are equal to the following two projections:βi +1=max min αi +λi 2μ2,1 ,0 ,D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h pv i +1=max −∂x αi −λi 3μ2,0 .The Euler–Lagrange equation for problem (3.10)is given by the following equation:(3.13)((2+μ2)∂T x∂x +2∂Ty ∂y+μ2I )α=−∂T x (μ2v i +1+λ3)−λ2−g +μ2βi +1.The coefficient matrix of the above linear system is given by the second-order differential operator with the Dirichlet and Neumann boundary conditions.This matrix equation can also be solved by many efficient numerical solvers and preconditioning techniques;see the book by Saad [30].In particular,we consider the overlapping region to be rectangular like the example inFigure 1.In this case,the Dirichlet boundary conditions are on the blue line and the red line,and the Neumann boundary conditions are on the green lines.In [31],it was mentionedthat the linear system in (3.13)arising from the differential operator with the Dirichlet and Neumann boundary conditions can be diagonalized by using the discrete sine transform matrix in the x -direction and the discrete cosine transform matrix in the y -direction.As we can make use of fast sine and cosine transforms to implement these computation tasks,the cost of solving such a linear system is O (N log N ),where N is the number of pixels in the overlapping region.3.3.Convergence analysis.In this subsection,we study the convergence of Algorithm 1.We first note that the convergence for ADMM [32]can be used here for the proposed method-ology for solving the problem in (3.4)and (3.8).We conclude the convergence for ADMM asthe following theorem.Theorem 3.1.For the problem in (3.4)and (3.8),let u 0,λ01and α0,λ02,λ03be arbitrary,andlet μ1,μ2>0.Then the sequences {u k ,w k ,λk 1}and {(αk ,βk ,v k ,λk 2,λk 3)}generated by (3.5)–(3.7)and (3.9)–(3.12)converge to (u ∗,w ∗,λ∗1)and (α∗,β∗,v ∗,λ∗2,λ∗3),which are saddle points of L 1and L 2.Now we can show the convergence of Algorithm 1.Theorem 3.2.Suppose I 1,I 2∈W 1,2(Ωc ).Let ˆΛ= (α,u,v )|(α,u,v )∈W 1,2(Ω)×W 1,2(Ωc )×L 2(Ω),0≤α≤1,α=f on D,∂α∂ n =0on N ,∂x α+v =0,u |∂Ω1∩∂Ωc =I 1,u |∂Ω2∩∂Ωc =I 2,∂u∂ n =0on N ,I min ≤u ≤I max ,v ≥0 .Let {(αk ,u k ,v k )}be the sequence generated by Algorithm 1.Then {(αk ,u k ,v k )}converges to(α∗,u ∗,v ∗)∈ˆΛ(up to a subsequence),and for any (α,u,v )∈ˆΛ,we have E 2(α∗,u ∗,v ∗)≤E 2(α,u ∗,v ),E 2(α∗,u ∗,v ∗)≤E 2(α∗,u,v ∗).Proof .First,we deduce the following inequality from Algorithm 1:E 2(αk +1,u k +1,v k +1)≤E 2(αk ,uk +1,v k)≤E 2(αk ,u k ,v k ).D o w n l o a d e d 10/31/13 t o 218.201.115.245. R e d i s t r i b u t i o n s u b j e c t t o S I A M l i c e n s e o r c o p y r i g h t ; s e e h t t p ://w w w .s i a m .o r g /j o u r n a l s /o j s a .p h p。