AN ANALYSIS OFTHE OPTIMIZATION FORMULATION OF ELASTIC INVERSE PROBLEMS Carlos E. Rivas 1 , Paul E. Barbone 1 , Assad A. Oberai 2 , Boston University, Boston MA 02215. Abstract Inverse problems are often formulated as optimization problems. One seeks the parameter distribution which, when used in a forward model of the problem, gives the best match possible to the measured data. We consider the problem of imaging the elastic modulus distributions of soft tissues in this context. We consider in particular formulating this inverse problem as a constrained optimization problem in which we minimize the objective functional (the data mismatch) with the plane stress incom- pressible elasticity equations as constraints. The Lagrange multiplier method is applied resulting in a two-field variational formulation that is often called a mixed formula- tion. We analyze the well-posedness of this problem. Under some relatively strong mathematical conditions, we prove that the continuous optimization formulation of this inverse problem is well-posed. In practice, however, when we solve this problem approximately by using classical discretization techniques, problems like oscillations and spurious non-physical solutions appear. Finally, we introduce stabilization in the discrete optimization problem to control the spurious solutions. In the context of in- verse problems, stabilization must be contrasted to regularization. The latter is used to render a continuous problem well-posed. Stabilization is used to make well-behaved a particular discretization of a well-posed mathematical problem. We expect this research to improve computational approaches to nearly all inverse problems solved via optimization methods. In the context of elastography, this research will improve our ability to screen, diagnose, and treat breast cancer, prostate cancer, and other tissue pathologies. Introduction Inverse problems are often formulated as optimization problems [3]: One seeks the parameter distribution which, when used in a forward model of the imaging problem, gives the best match possible to the measured data. In practice, the inverse problem is solved as follows: 1. One guesses a parameter distribution. 2. One solves the “forward problem” for this guessed parameter distribution. This forward solve requires a convergent, stable numerical method. 3. One compares the predicted field to the measured field. If they agree sat- isfactorily, the guessed parameter distribution is taken to be correct. If the agreement is unsatisfactory, the parameter distribution is updated and the cycle is repeated. Step 2 requires the solution of the forward problem on a discrete mesh. One expects that as that mesh is refined, the solution approaches that of the continuous equations. This is true of a convergent numerical method. We’ve found that even if a stable, optimally convergent numerical method is used in step 2, the discrete parameter distribution may not converge to the exact solution with mesh refinement. This turns out to be a generic feature of such constrained minimization problems, and thus applies equally well to elastic- ity imaging, diffuse optical tomography, electrical impedance tomography, or any inverse problem which is formulated as a constrained minimization problem. The measured data for the inverse elasticity problem is the “strain image”. It is usually obtained from ultrasound, but other imaging modalities may also be used. Figure 1 shows a scheme of how ultrasound images are obtained and how to used them to get displacements. Figure 2 shows an example of an input data and a solution of an inverse elasticity problem. There we can see the strain image and the corresponding elastic modulus reconstruction. The inclusion is a stiff 5 mm diameter cylinder with independently measured stiffness contrast of about 2.5:1. Figure 1: Image Registration In image registration undeformed and deformed images are compared to find displacements (a) (b) Figure 2: Inverse elasticity example (a) Strain image (b) Shear modulus reconstruction Background and Motivation • Elasticity imaging or elastography is a technique which uses imaging modalities to measure and image relative displacements within tissue vol- umes. Its main goal is to map the elastic modulus (stiffness) of soft tissues. Tissue stiffness or elastic modulus distribution is believed to be clinically sig- nificant and may be correlated to tumor and other tissue pathologies[1]. • In collaboration with MedBED, and a host of researchers nationally and internationally, we have demonstrated the ability to reconstruct elastic mod- ulus distributions from image data. We have worked with simulated images, images of phantoms, excised tissues, and in vivo (clinical) data, and we have worked with ultrasound and x-ray modalities. • Stabilization techniques as [2, 4] have been developed to overcome prob- lems that polute numerical solutions when some stability conditions are vio- lated. However, these techniques only work when the corresponding continu- ous problem is mathematically well posed. State of the Art • We continue to lead the field of elasticity imaging in terms of quantitative modulus reconstructions in 2D and 3D, mathematical analysis of the inverse problem, and extensions to hyperelasticity and multiphase phenomena. • The current work builds on our experience with 2D and 3D elastic mod- ulus reconstructions from both ultrasound and x-ray imaging data. Imaging data is available to us through MedBED, and through our collaborators. Formulation of the Inverse problem Constrained optimization problem Given ¯ u(x), find u(x), μ(x) and λ(x) that minimizes: Π(u, μ, λ)= 1 2 ‖u − ¯ u‖ 2 Q + a(λ, u; μ)+ τ 2 ‖L(u, μ)‖ 2 Ω ′ + α 2 ‖∇μ‖ 2 Here, ¯ u and u are the measured and the predicted displacement field, λ is the Lagrange multiplier, μ is the shear modulus and a(λ, u; μ) is the weak form of the constrained equation. The last two components of this functional are the stabilization and the regularization terms, respectively. Q defines the norm in which the misfit between u and ¯ u is calculated and it is critical in defining the well posedness of this formulation. Newton and Quasi-Newton Methods We have used both Newton and quasi-Newton [3] iteration methods to solve the nonlinear optimization problem. Newton’s method leads to the following weak equations for the updates δu, δλ, and δμ: (δu, v ) Q +τ (L(δu, ˆ μ),L(v, ˆ μ)) + a( ˆ λ, v, δμ)+τb(δμ, v )+a(δλ, v, ˆ μ)=(f,v ) ∀v a( ˆ λ, δu, γ )+ τb(γ,δu) + c(δμ, γ ; τ,α) + a(δλ, ˆ u, γ )=(g,γ ) , ∀γ a(w, δu, ˆ μ) + a(w, ˆ u, δμ) =(h, w ) , ∀w Here v ,w , and γ are the variations of u, λ, and μ, respectively. Brezzi’s theorem If there exists positive constants c 1 and c 2 such that the following conditions are satisfied, then the problem is well posed A(δu, δμ; δu, δμ) ≥ c 1 ‖|δu, δμ|‖ 2 ∀ (δu, δμ) ∈K sup B (w ; δu, δμ) ‖|δu, δμ|‖ ≥ c 2 ‖w ‖ ∀w ∈V , (δu, δμ) ∈V×M Here, V and M are appropriate function spaces and K is the kernel space that relates δu and δμ. The operators A and B are defined as: A(δu, δμ; δu, δμ)= ‖δu‖ 2 + τ ‖L(δu, ˆ μ)‖ 2 +2a( ˆ λ, δu, δμ)+ 2τb(δμ, δu)+ c(δμ, δμ; τ,α) B (w ; δu, δμ)= a(w, δu, ˆ μ)+ a(w, ˆ u, δμ) Results (a) (b) Figure 3: Computational experiment (a) Uniaxial compression in a domain with a stiffer inclusion. (b) Pure shear with linear variation of μ 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 Strain image: ε 22 x (mm) y(mm) -6.5 -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 x 10 -3 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 μ recostruction x (mm) y(mm) 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 (a) (b) Figure 4: For a domain with stiffer inclusion (a) Input strain image (b) Output shear modulus reconstruction μ distribution x (m) y (m) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 -0.2 0 0.2 0.4 0.6 0.8 1 μ distribution x (m) y (m) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 -5 0 5 10 15 (a) (b) Figure 5: Shear modulus reconstruction without stabilization (a) For an homogeneous μ (b) For a linear distribution of μ 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 Strain image: ε 12 x(mm) y(mm) 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -3 μ distribution x(mm) y(mm) 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 (a) (b) Figure 6: For a linear distribution of μ (a) Input strain image (b) Output shear modulus reconstruction Discussion • The elastic modulus reconstructions shown in Figures 4b and 6b show clearly the stiffer inclusion and the linear distribution of the shear modu- lus, respectively. These solutions were obtained assuming a homogeneous distribution for μ as the initial guess. The GLS term employed in these re- constructions helped to reduce the effect of any spurious discrete solutions. Figures 5a,b show the reconstruction for the pure shear case without GLS. The spurious solutions are present in these reconstructions • From the theoretical point of view, we can demonstrate that the inverse problem, as formulated in this work, is well posed without regularization but only for the plane stress incompressible case. This is true only when the displacement field is noiseless and when the misfit between the predicted dis- placement and the measured displacement is evaluated in the H 2 Sobolev norm. When we use the H 1 norm, we need regularization in order to guaranty well posedness. Future Work • We plan to prove convergence and stability analytically for the plane strain and 3D elasticity inverse problems. In doing that we have to identify the cor- rect norm to minimize and the form of the regularization terms to be used • Further analysis will be done to determine a good approximation for the stabilization parameter τ such that we guarantee optimal performance. Opportunities for technology transfer • In the context of elastography, by imaging in vivo distributions of the biomechanical properties of soft tissues, we expect this research to enable clinical applications in the diagnosis and treatment of breast, prostate cancer and other soft tissue pathologies. In the context of computational methods, we will develop an accurate, efficient and stable method for solving the inverse problems. Collaborators Jeffrey C. Bamber 3 , Timothy J. Hall 4 , Isaac Harari 5 , Ricardo Leiderman 6 , Elise F. Morgan 1 , Michael S. Richards 7 Acknowledgments Three-Level Diagram. This work was supported in part by CenSSIS (The Center for Subsurface Sensing and Imaging Systems) under the Engineering Research Centers Pro- gram of the National Science Founda- tion (Award No. EEC-9986821). References [1] Jeffrey C. Bamber, Paul E. Barbone, Nigel L. Bush, David O. Cosgrove, Marvin M. Doyley, Frank G. Fueschel, Paul M. Meaney, Naomi R. Miller, Tsuyoshi Shiina, and Francois Tranquart. Progress in freehand elastogra- phy of the breast. IEICE Transactions on Information and Systems, E85- D(1):5–15, 2002. [2]Leopoldo P. Franca and Thomas J.R. Hughes. Two classes of mixed finite element methods. Computer methods in applied mechanics and engineer- ing, 69:(89–129), 1988. [3] Assad A Oberai, Nachiket K Gokhale, and Gonzalo R Feij´ oo. Solution of inverse problems in elasticity imaging using the adjoint method. Inverse Problems, 19:297–313, 2003. [4]L.P. Franca T.J.R. Hughes and M. Balestra. A new finite element formu- lation for computational fluid dynamics: V. circumventing the babuska- brezzi condition: A stable petrov-galerkin formulation of the stokes prob- lem accomodating equal-order interpolations. Computer methods in ap- plied mechanics and engineering, 59:(85–99), 1986. 1 Department of Aerospace and Mechanical Engineering, Boston University, Boston MA 2 Department of Mechanical, Aerospace and Nuclear Engineering, RPI, Troy NY. 3 Royal Marsden Hospital and Institute of Cancer Research, Sutton, Surrey, UK 4 Biomedical Engineering, University of Wisconsin, Madison, WI. 5 Department of Solids and Structures, Tel Aviv University, Tel Aviv, Israel. 6 Mechanical Engineering Program, COPPE, Universidade Federal do Rio de Janeiro, RJ Brazil. 7 Department of Radiology, University of Michigan, Ann Arbor MI.