Mathematical modeling in full-field optical coherence elastography
Habib Ammari, Elie Bretin, Pierre Millien, Laurent Seppecher, Jin-Keun Seo
MMathematical modeling in full-field optical coherenceelastography ∗ Habib Ammari † Elie Bretin ‡ Pierre Millien † Laurent Seppecher † Jin-Keun Seo § October 3, 2018
Abstract
We provide a mathematical analysis of and a numerical framework for full-fieldoptical coherence elastography, which has unique features including micron-scaleresolution, real-time processing, and non-invasive imaging. We develop a novel al-gorithm for transforming volumetric optical images before and after the mechanicalsolicitation of a sample with sub-cellular resolution into quantitative shear modulusdistributions. This has the potential to improve sensitivities and specificities in thebiological and clinical applications of optical coherence tomography.
Mathematics Subject Classification (MSC2000): 35R30, 35B30.Keywords: full-field optical coherence tomography, elastography, hybrid imaging, optimal control, high-resolution shear modulus imaging, biological tissues.
Optical coherence tomography (OCT) is a non-invasive and a non-ionizing imaging tech-nique that produces high-resolution images of biological tissues. It performs optical slicingin the sample, to allow three-dimensional reconstructions of internal structures. Conven-tional optical coherence time-domain and frequency-domain tomographies require trans-verse scanning of the illumination spot in one or two directions to obtain cross-sectional oren face images, respectively. Full-field OCT allows OCT to be performed without trans-verse scanning; the tomographic images are obtained by combining interferometric imagesacquired in parallel using an image sensor. Both the transverse and the axial resolutionsare of the order of 1 µ m; see [9, 10]. We refer to [11] for the mathematical modeling ofOCT. ∗ This work was supported by the ERC Advanced Grant Project MULTIMOD–267184. † Department of Mathematics and Applications, Ecole Normale Sup´erieure, 45 Rue d’Ulm, 75005 Paris,France ([email protected], [email protected], [email protected]). ‡ Institut Camille Jordan, Universit´e de Lyon, Lyon, F-69003, France ([email protected]). § Department of Computational Science and Engineering, Yonsei University, 50 Yonsei-Ro, Seodaemun-Gu, Seoul 120-749, Korea ([email protected]). a r X i v : . [ m a t h . O C ] M a y lastography is an imaging-based technique for the estimation of the elastic propertiesof tissues. Given that the mechanical properties of tissues and cells are related to theirstructure and function, changes in those properties can reflect healthy or pathologicalstates such as weakening of vessel walls or cirrhosis of the liver. Elastography can aid theidentification of suspicious lesions, the diagnosis of various diseases and the monitoringof the effectiveness of treatments (see [15, 16]). Different imaging modalities (e.g., ultra-sound and magnetic resonance imaging) can be used to measure tissue displacements andto estimate the resulting tissue stiffness and viscosity. Magnetic resonance elastographyis relatively expensive, due to the high magnetic field environment, which requires specif-ically designed equipment. Several reconstruction approaches for elastography have beenderived [3, 4, 5, 22].In [18], elastographic contrast has been combined with full-field OCT with the aimof creating a virtual palpation map at the micrometer scale. The idea is to register avolumetric optical image before and after mechanical solicitation of the sample. Based onthe assumption that the density of the optical scatterers is advected by the deformation,the displacement map can be first estimated. Then, using a quasi-incompressible modelfor the tissue elasticity, the shear modulus distribution can be reconstructed from theestimated displacement map.The OCT elastography is able to perform displacement measurements with sub-cellularresolution. It enables a more precise characterization of tissues than that achieved usingultrasound or magnetic resonance elastography; therefore, it provides a more accurateassessment of microscale variations of elastic properties. A map of mechanical proper-ties added as a supplementary contrast mechanism to morphological images could aiddiagnosis. The technique costs less than other elastography techniques.The mapping of mechanical properties was first introduced to OCT imaging by Schmitt[21], who measured displacements as small as a few micrometers in heterogeneous gelatinphantoms containing scattering particles in addition to living skin. Various subsequentapplications have employed OCT methods in elastography; these include dynamic andfull-field optical coherence elastography (see [14, 19, 20]).In all of the aforementioned techniques, transforming the OCT images before and afterthe application of a load into quantitative maps of the shear modulus is a challengingproblem.In this paper we present a mathematical and numerical framework for the OCT-elastography experiment described in [18]. Using the set of images before and aftermechanical solicitation we design a novel method to reconstruct the shear modulus dis-tribution inside the sample.To mathematically formulate the problem, let Ω ⊂ R d , d = 2 , , and let ε be theknown piecewise smooth optical index of the medium, and µ be its shear modulus. In thispaper we consider heterogeneous (unknown) shear modulus distributions. The mediumis solicited mechanically. Since compression modulus of biological media is four order ofmagnitude larger than the shear modulus, it can be shown that the displacement map u obeys the linearized equations of incompressible fluids or the Stokes system [3, 4, 5].The model problem is then the following Stokes system in a heteregeneous medium whichreads: ∇ · (cid:0) µ ( ∇ u + ∇ u T ) (cid:1) + ∇ p = 0 in Ω , ∇ · u = 0 in Ω , u = f on ∂ Ω , (1.1)2here superposed T denotes the transpose and the real-valued vector f satisfies the com-patibility condition (cid:82) ∂ Ω f · ν = 0 with ν being the outward normal at ∂ Ω .Throughout this paper, we assume that µ ∈ C , (Ω ) and f ∈ C ( ∂ Ω ) d . From [7, 12,13], (1.1) has a unique solution u ∈ C (Ω ) d . Moreover, there exists a positive constant C depending only on µ and Ω such that || u || C (Ω ) d ≤ C || f || C ( ∂ Ω ) d . Using a second OCT scan, one has access to the optical index of the deformed medium ε u ( (cid:101) x ) , ∀ (cid:101) x ∈ Ω u , where Ω u is defined byΩ u = { x + u ( x ) , x ∈ Ω } . The new optical index is linked to the original one by ε ( x ) = ε u ( x + u ( x )) , ∀ x ∈ Ω . (1.2)The goal is to reconstruct the shear modulus map µ on Ω from the functions ε and ε u . We first prove that, in two dimensions, if the direction of ∇ ε |∇ ε | is not constant in aneighborhood of x , then the displacement field u at x can be approximately reconstructed.In three dimensions, one shall assume that the vectors ∇ ε ( y ) |∇ ε ( y ) | are not coplanar for y a neighborhood of x . Hence, the reconstructed value of u ( x ) serves as an initial guessfor the minimization of the discrepancy between computed and measured changes in theoptical index. Then, we compute an element of the subgradient [8] of the discrepancyfunctional. Finally, we implement a minimization scheme to retrieve the shear modulusmap from the reconstructed displacements.The paper is organized as follows. Section 2 is devoted to some mathematical prelim-inaries. In section 3 we consider piecewise smooth ε functions and first derive a leading-order Taylor expansion of ε u as || u || C goes to zero. Then we provide an initial guess bylinearization. Finally, we prove the Fr´echet differentiability of the discrepancy functionalbetween the measured and the computed advected images. The displacement field insidethe sample can be obtained as the minimizer of such functional. Section 4 is devoted tothe reconstruction of the shear modulus from the displacement measurements. In section5 we present some numerical results to highlight the viability and the performance of theproposed algorithm. The paper ends with a short discussion. Let Ω be a bounded smooth domain in R d , d = 2 ,
3. We start by defining a class ofpiecewise smooth functions.
Definition 2.1
For any k ∈ N , α ∈ ]0 , , for any curve S of class C ,α for some < α < such that Ω \ S is a union of connected domains Ω i , i = 1 , , · · · n , we define C k,αS (cid:0) Ω (cid:1) tobe the class of functions f : Ω −→ R satisfying f | Ω i ∈ C k,αS (cid:0) Ω i (cid:1) ∀ i = 1 , · · · n. (2.1)3 efinition 2.2 We define
BV(Ω) as the subspace of L (Ω) of all the functions f whoseweak derivative Df is a finite Radon measure. In other terms, f satisfies (cid:90) Ω f ∇ · F ≤ C sup x ∈ Ω | F | , ∀ F ∈ C (Ω) d for some positive constant C with C (Ω) being the set of compactly supported C functions.The derivative of a function f ∈ BV(Ω) can be decomposed as Df = ∇ f H d + [ f ] ν s H d − S + D c f, where H d is the Lebesgue measure on Ω, H d − S is the surface Hausdorff measure on arectifiable surface S , ν S is a normal vector defined a.e. on S , ∇ f ∈ L (Ω) is the smoothderivative of f , [ f ] ∈ L ( S, H d − S ) is the jump of f across S and D c f is a vector measuresupported on a set of Hausdorff dimension less than ( d − d − Definition 2.3
We define
SBV(Ω) as the subspace of
BV(Ω) of all the functions f sat-isfying D c f = 0 . Definition 2.4
For any ≤ p ≤ + ∞ , we define SBV p (Ω) = (cid:8) f ∈ SBV(Ω) ∩ L p (Ω) , ∇ f ∈ L p (Ω) d (cid:9) . As SBV(Ω) is a good model for piecewise- W , functions, the space SBV p (Ω) canbe seen as the space of piecewise- W ,p functions. Here, W ,p (Ω) = { f ∈ L p (Ω) , ∇ f ∈ L p (Ω) d } for p ≥ ∞ (Ω) is a nice definition of piecewise Lipschitz function. Notealso that C k,αS (cid:0) Ω (cid:1) ⊂ SBV p (Ω).From now on, we assume that the optical index in the medium ε belongs to C k,αS (cid:0) Ω (cid:1) ,which is a simple but good model for a discontinuous medium. Some of the followingpropositions are true for more general maps ε ∈ SBV(Ω). In these propositions we onlyassume that ε is in SBV(Ω). Let Ω (cid:98) (Ω ∩ Ω u ) be a smooth simply connected domain. On Ω, we have ε u = ε ◦ ( I + u ) − ε = ε u ◦ ( I + u ) , where I is the d × d identity matrix. Proposition 3.1
Let ε ∈ BV(Ω) and let u ∈ C (Ω) d be such that (cid:107) u (cid:107) C (Ω) d < . Then,for any ψ ∈ C (Ω) , we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω ( ε − ε u ) ψ − (cid:90) Ω ψ u · Dε (cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:107) u (cid:107) C (Ω) d (cid:107) u (cid:107) C (Ω) d (cid:107) ψ (cid:107) C (Ω) | ε | TV(Ω) , (3.1)4 here the constant C is independent of ψ and | | TV(Ω) denotes the total variation semi-norm. Estimate (3.1) yields that ε u − ε + u · Dε (cid:107) u (cid:107) C (Ω) d weakly converges to in C (Ω) when (cid:107) u (cid:107) C (Ω) d goes to . Proof.
For each t ∈ [0 , φ t by φ − t ( x ) = x + t u ( x ). Let η > ε ( η ) be a smooth function such that (cid:107) ε − ε ( η ) (cid:107) L (Ω) →
0, and | ε ( η ) | TV(Ω) →| ε | TV(Ω) as η →
0. Analogously, we define ε ( η ) u to be the smooth approximation of ε u givenby ε ( η ) u ( x ) = ε ( η ) ◦ φ ( x ) . From ε ( η ) u ( x ) − ε ( η ) ( x ) = (cid:0) ε ( η ) ◦ φ (cid:1) ( x ) − (cid:0) ε ( η ) ◦ φ (cid:1) ( x ) , ∀ x ∈ Ω , we have ε ( η ) u ( x ) − ε ( η ) ( x ) = (cid:90) ∇ ε ( η ) ( φ t ( x )) · ∂ t φ t ( x ) dt, ∀ x ∈ Ω . Therefore, for ψ ∈ C ∞ (Ω) with C ∞ (Ω) being the set of compactly supported C ∞ functions, (cid:90) Ω (cid:2) ε ( η ) u ( x ) − ε ( η ) ( x ) + ∇ ε ( η ) ( x ) · u ( x ) (cid:3) ψ ( x ) d x = (cid:90) Ω (cid:20)(cid:90) ∇ ε ( η ) ( φ t ( x )) · ∂ t φ t ( x ) dt (cid:21) ψ ( x ) d x + (cid:90) Ω ∇ ε ( η ) ( x ) · u ( x ) ψ ( x ) d x , ∀ x ∈ Ω . (3.2)By a change of variables in the first integral and using the fact that ∂ t φ t ( x ) = − ∂ x φ t ( x ) ∂ t φ − t ( y ) | y = φ t ( x ) , we get, for all x ∈ Ω, (cid:90) (cid:20)(cid:90) Ω ∇ ε ( η ) ( φ t ( x )) · ∂ t φ t ( x ) ψ ( x ) d x (cid:21) dt = − (cid:90) (cid:90) Ω ∇ ε ( η ) ( y ) · (cid:2) ∂ x φ t ( φ − t ( y )) ∂ t φ − t ( y ) (cid:3) | det ∂ x φ − t ( y ) | ψ (cid:0) φ − t ( y ) (cid:1) d y dt. Here, det denotes the determinant of a matrix. Since ∀ ( y , t ) ∈ Ω × [0 , , ∂ t φ − t ( y ) = u ( y ) ,∂ y φ − t ( y ) = I + t ∇ u ( y ) , and ∂ x φ t ( φ − t ( y )) ∂ y φ − t ( y ) = I , we can write (cid:90) (cid:90) Ω (cid:2) ∇ ε ( η ) ( φ t ( x )) · ∂ t φ t ( x ) ψ ( x ) d x (cid:3) dt = − (cid:90) (cid:90) Ω ∇ ε ( η ) ( y ) · (cid:2) ( I + t ∇ u ( y )) − u ( y ) (cid:3) | det I + t ∇ u ( y ) | ψ (cid:0) φ − t ( y ) (cid:1) d y dt, (cid:90) Ω (cid:2) ε ( η ) u ( x ) − ε ( η ) ( x ) + ∇ ε ( η ) ( x ) · u ( x ) (cid:3) ψ ( x ) d x = (cid:90) (cid:90) Ω ∇ ε ( η ) ( x ) · u ( x ) (cid:2) ψ ( x ) − ψ (cid:0) φ − t ( x ) (cid:1) (cid:3) d x dt + (cid:90) (cid:90) Ω ∇ ε ( η ) ( x ) · (cid:0)(cid:2) ( I + t ∇ u ( x )) − | det I + t ∇ u ( x ) | − I (cid:3) u ( x ) (cid:1) ψ (cid:0) φ − t ( x ) (cid:1) d x dt. (3.3)The first term in the right-hand side of (3.3) can be estimated as follows: (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) (cid:90) Ω ∇ ε ( η ) ( x ) · u ( x ) (cid:2) ψ ( x ) − ψ (cid:0) φ − t ( x ) (cid:1) (cid:3) d x dt (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) u (cid:107) C (Ω) d (cid:107)∇ ε ( η ) (cid:107) L (Ω) d (cid:107)∇ ψ (cid:107) C (Ω) d . Let tr denote the trace of a matrix. Using the fact that( I + t ∇ u ) − = (cid:88) i =0 ( − i ( t ∇ u ) i , which follows from || u || C (Ω) d <
1, anddet ( I + t ∇ u ) = − tr t ∇ u + det t ∇ u if d = 2 , t ∇ u − (cid:2) (tr t ∇ u ) − tr ( t ∇ u ) (cid:3) + det t ∇ u if d = 3 , we get (cid:90) (cid:90) Ω ∇ ε ( η ) ( x ) · u ( x ) (cid:2) ( I + t ∇ u ( x )) − | det I + t ∇ u ( x ) | − I (cid:3) ψ (cid:0) φ − t ( x ) (cid:1) dxdt ≤ (cid:107) u (cid:107) C (Ω) d (cid:107) u (cid:107) C (Ω) d (cid:107)∇ ε ( η ) (cid:107) L (Ω) d (cid:107) ψ (cid:107) C (Ω) , which is the desired estimate for the second term in the right-hand side of (3.3).Now, we can deduce the final result by density when η →
0. Since u ∈ C (Ω) d and ψ ∈ C (Ω), we can write (cid:90) Ω ψ u · ∇ ε ( η ) = − (cid:90) Ω ∇ · ( ψ u ) ε ( η ) . Since (cid:107) ε ( η ) − ε (cid:107) L (Ω) →
0, we have (cid:90) Ω ∇ · ( ψ u ) ε ( η ) → (cid:90) Ω ∇ · ( ψ u ) ε. As | ε ( η ) | TV(Ω) → | ε | TV(Ω) , we arrive at (3.1) and the proof of the proposition is complete. (cid:3) .2 Local recovery via linearization Assuming that ε ∈ SBV (Ω), we can write Dε = ∇ ε H d + [ ε ] S ν S H d − S , where ν S is the outward normal at the oriented surface S of discontinuity of ε .The data consists of ε and ε u on Ω. In order to reconstruct u , we can use the firstorder approximation of ε − ε u : ε − ε u ≈ u · Dε, given by Proposition 3.1. These data can be decomposed into two parts: u · Dε ( · ) = u · ∇ ε H d + [ ε ] S u · ν S H d − S = d reg H d + d sing H d − S . Let w be a mollifier supported on [ − , δ >
0, we define w δ = 1 δ d w (cid:16) · δ (cid:17) , and introduce u δ ( x ) = (cid:90) Ω u ( y ) w δ ( | x − y | ) d y . Since u is smooth, for any x ∈ Ω, u δ ( x ) is a good approximation of u on the ball withcenter x and radius δ .We want to find an approximate value for u δ from the optical measurements and use itas an initial guess in an optimization procedure. For doing so, we introduce the functional J x : R d −→ R given by u (cid:55)−→ J x ( u ) = (cid:90) Ω |∇ ε ( y ) · u − d reg ( y ) | w δ ( | x − y | ) d y + (cid:90) Ω | [ ε ] S u · ν S − d sing ( y ) | w δ ( | x − y | ) d y , and look for minimizers of J x in R d . The gradient of J x can be explicitly computed asfollows: ∇ J x ( u ) = 2 (cid:90) Ω ( ∇ ε ( y ) · u − d reg ( y )) ∇ ε ( y ) w δ ( | x − y | ) d y + 2 (cid:90) Ω ([ ε ] S ( y ) u · ν ( y ) − d sing ( y )) [ ε ] S ( y ) ν ( y ) w δ ( | x − y | ) d y . In the case where ε has no jumps, J x is a quadratic functional and we have ∇ J x ( u ) = 0 ⇔ u T (cid:18)(cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ) d y (cid:19) = (cid:90) x + δB d reg ( y ) w δ ( | x − y | ) ∇ ε ( y ) d y , (3.4)where B is the ball with center 0 and radius 1.If the matrix (cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ) is invertible, then the minimizer is given by u T = (cid:18)(cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ) d y (cid:19) − (cid:90) x + δB d reg w δ ( | x − y | ) ∇ ε ( y ) d y . (3.5)7he following proposition gives a sufficient condition for the invertibilty of the matrix (cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ). Proposition 3.2
Suppose that ε has no jumps and d = 2 . Assume x + δB ⊂ Ω . Then,if all vectors ∇ ε in { y : w δ ( | y − x | ) (cid:54) = 0 } are not collinear, then the matrix (cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ) d y is invertible. Proof.
Writing ∀ y ∈ x + δB, ∇ ε ( y ) = u ( y ) e + v ( y ) e , where { e , e } is the cannonical basis of R , it follows that ∇ ε ∇ ε T ( y ) = u ( y ) e e T + v ( y ) e e T + u ( y ) v ( y ) (cid:0) e e T + e e T (cid:1) , ∀ y ∈ x + δB. Computing the convolution with respect to w δ , we get (cid:90) Ω w δ ( | y − y | ) ∇ ε ( y ) ∇ ε T ( y ) d y = (cid:18)(cid:90) Ω u ( y ) w δ ( | y − x | ) d y (cid:19) e e T + (cid:18)(cid:90) Ω v ( y ) w δ ( | y − x | ) d y (cid:19) e e T + (cid:18)(cid:90) Ω u ( y ) v ( y ) w Tδ ( | y − x | ) d y (cid:19) (cid:0) e e T + e e T (cid:1) . This matrix is not invertible if and only if (cid:18)(cid:90) Ω u ( y ) w δ ( | y − x | ) d y (cid:19) (cid:18)(cid:90) Ω v ( y ) w δ ( | y − x | ) d y (cid:19) = (cid:18)(cid:90) Ω u ( y ) v ( y ) w δ ( | y − x | ) d y (cid:19) , which is exactly the equality case in weighted Cauchy-Schwarz inequality. So, if thereexist two points y , y ∈ { y : w δ ( | y − x | ) (cid:54) = 0 } such that ∇ ε ( y ) × ∇ ε ( y ) (cid:54) = 0, then u is not proportional to v , and the matrix is invertible. (cid:3) Remark 3.3
Assuming that ∇ ε ( y ) (cid:54) = 0 for y ∈ x + δB ⊂ Ω , Proposition 3.2 gives thatthe direction of ∇ ε |∇ ε | in not constant in x + δB ⊂ Ω if and only if (cid:90) x + δB ∇ ε ( y ) ∇ ε T ( y ) d y is invertible.Hence, under the above condition on ε in the neighborhood x + δB , the displacement field u at x can be approximately reconstructed. Remark 3.4
By exactly the same arguments as those in two dimensions, one can provethat in the three-dimensional case, if all vectors ∇ ε in { y : w δ ( | y − x | ) (cid:54) = 0 } are notcoplanar, then the matrix (cid:90) Ω w δ ( | x − y | ) ∇ ε ( y ) ∇ ε T ( y ) d y is invertible.On the other hand, in the case where ε is piecewise smooth, one can first detect thesurface of jumps of ε using for example an edge detection algorithm [6, 17] and then applythe proposed local algorithm in order to have a good approximation of u in the domainswhere ε is smooth. .3 Minimization of the discrepancy functional Let ε ∈ C k,αS (cid:0) Ω (cid:1) , where S is the surface of discontinuity. For the sake of simplicity weassume that Ω \ S is the union of two connected domains Ω ∪ Ω . Therefore, ε can bewritten as ε = ε χ Ω + ε χ Ω (3.6)with ε i ∈ C (Ω i ), for i = 1 , u ∗ the applied (true) displacement on Ω (as defined in (1.1)) and (cid:101) ε the mea-sured deformed optical index given by (cid:101) ε = ε ◦ ( I + u ∗ ) − . The following result holds.
Proposition 3.5
Let ε verify (3.6), u ∗ ∈ C (Ω) d be the solution of (1.1), and (cid:101) ε = ε ◦ ( I + u ∗ ) − . Suppose that Ω (cid:98) Ω . Then, the functional I defined by I : C (Ω) d −→ R , u (cid:55)−→ I ( u ) = (cid:90) Ω | (cid:101) ε ◦ ( I + u ) − ε | d x (3.7) has a nonempty subgradient. Let ξ in the dual of C (Ω) d be given by ξ : h (cid:55)→ (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x )] h ( x ) · D (cid:101) ε ◦ ( I + u )( x ) d x . (3.8) For || h || C (Ω) d small enough, we have I ( u + h ) − I ( u ) ≥ ( ξ , h ) , where ( , ) is the duality product between C (Ω) d and its dual, which means that ξ ∈ ∂I with ∂I being the subgradient of I . Remark 3.6
It is worth emphasizing that if ε has no jump, then I is Fr´echet differentiableand ξ is its Fr´echet derivative. Remark 3.7
Under the assumptions of Proposition 3.5, if u ∗ is small enough (in C -norm), then (cid:101) ε = ε ◦ ( I + u ∗ ) − can be written as (cid:101) ε = (cid:101) ε + (cid:101) ε χ ˜Ω , (3.9) with (cid:101) ε ∈ C (Ω) and (cid:101) ε ∈ C (Ω) . In the sequel, we shall define ˜Ω i = ( I + u ∗ ) (Ω i ) and ˜ f i = ε i ◦ ( I + u ∗ ) − . For doing so, we extend ˜ f into a function (cid:101) ε defined on the wholedomain such that (cid:101) ε ∈ C (Ω) and (cid:101) ε (cid:12)(cid:12) ˜Ω = ˜ f . Then, we set (cid:101) ε = ˜ f − (cid:101) ε on ˜Ω . Finally,we extend (cid:101) ε into a compactly supported C -function on the whole domain Ω . We first prove the following lemma. 9 emma 3.8
Let u , h ∈ C (Ω) d and let (cid:101) ε be as in (3.9). Then, for (cid:107) u − u ∗ (cid:107) C (Ω) d and (cid:107) h (cid:107) C (Ω) d small enough, we have (cid:90) Ω [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x = (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x + o ( (cid:107) h (cid:107) C (Ω) d ) , (3.10) where δ ∂ ˜Ω is the Dirac distribution on ∂ ˜Ω and (cid:101) ε is defined in Remark 3.7. Proof.
We start by decomposing (cid:101) ε as follows: (cid:90) Ω [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x = (cid:90) Ω (cid:20)(cid:0)(cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u ) (cid:1) + (cid:0)(cid:101) ε ( x + u + h ) χ ˜Ω ( x + u + h ) − (cid:101) ε ( x + u ) χ ˜Ω ( x + u ) (cid:1)(cid:21) d x . Now, by developing the square, the first term can be estimated by (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω (cid:0)(cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u ) (cid:1) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) (cid:101) ε (cid:107) C (Ω) (cid:107) h (cid:107) C (Ω) d . Next, we write (cid:101) ε ( x + u + h ) χ ˜Ω ( x + u + h ) − (cid:101) ε ( x + u ) χ ˜Ω ( x + u ) = [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] χ ˜Ω ( x + u + h )+ (cid:2) χ ˜Ω ( x + u + h ) − χ ˜Ω ( x + u ) (cid:3) (cid:101) ε ( x + u ) . Since ( (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )) (cid:101) ε ( x + u ) ∈ C (Ω), Proposition 3.1 yields (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω (cid:2)(cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u ) (cid:3) (cid:2) χ ˜Ω ( x + u + h ) − χ ˜Ω ( x + u ) (cid:3) (cid:101) ε ( x + u ) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:18)(cid:90) Ω [ h · ∇ (cid:101) ε ( x + u )] d x (cid:19) / (cid:18)(cid:20)(cid:90) Ω [ h · ν (cid:101) ε ( x + u )] δ ∂ ˜Ω ( x + u ) d x (cid:21) + o ( (cid:107) h (cid:107) C (Ω) d ) (cid:19) / ≤ C (cid:107) h (cid:107) C (Ω) d . We now need to handle the last term (cid:90) Ω (cid:0) (cid:2) χ ˜Ω ( x + u + h ) − χ ˜Ω ( x + u ) (cid:3) (cid:101) ε ( x + u ) (cid:1) d x = (cid:90) Ω (cid:12)(cid:12) χ ˜Ω ( x + u + h ) − χ ˜Ω ( x + u ) (cid:12)(cid:12) (cid:101) ε ( x + u ) d x . Using Proposition 3.1, we obtain that (cid:90) Ω (cid:0) (cid:12)(cid:12) χ ˜Ω ( x + u + h ) − χ ˜Ω ( x + u ) (cid:12)(cid:12) (cid:101) ε ( x + u ) (cid:1) d x = (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x + o ( (cid:107) h (cid:107) C (Ω) d ) , which completes the proof of the lemma. (cid:3) We are now ready to prove Proposition 3.5.10 roof. If u ∈ C (Ω) and h ∈ C (Ω) , then we have I ( u + h ) − I ( u ) = (cid:90) Ω [ (cid:101) ε ( x + u + h ) + (cid:101) ε ( x + u ) − ε ( x )] [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x , and hence, I ( u + h ) − I ( u ) = (cid:90) Ω [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x + 2 (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x )] [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x . For any η >
0, let g ( η ) be a smooth, compactly supported function such that (cid:107) g ( η ) − [ (cid:101) ε ◦ ( I + u ) − ε ] (cid:107) L (Ω) < η and (cid:12)(cid:12) | g ( η ) | TV(Ω) − | (cid:101) ε ◦ ( I + u ) − ε | TV(Ω) (cid:12)(cid:12) < η ;see [2].Now, we write (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x )] [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x = (cid:90) Ω g η ( x ) [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x + (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x ) − g η ( x )] [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x . Let τ h be the translation operator. Then, τ h satisfies, for any h ∈ C (Ω) d , (cid:107) τ h [ f ] − f (cid:107) p ≤ C ( f ) (cid:107) h (cid:107) C (Ω) d , ∀ f ∈ SBV p (Ω) . (3.11)Using Cauchy-Schwartz’ inequality, we get (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x ) − g η ( x )] [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ Cη (cid:107) h (cid:107) C (Ω) d , (3.12)where C is a constant depending on (cid:101) ε, u , and Ω.We know that for a certain function ρ such that ρ ( s ) → s → (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω g η ( x ) [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] d x − (cid:90) Ω g η ( x ) h ( x ) · D ( (cid:101) ε ◦ ( I + u )) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) h (cid:107) C (Ω) d ρ ( (cid:107) h (cid:107) C (Ω) d ) . (3.13)Now, we have the following estimate: (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω g η ( x ) h ( x ) · D ( (cid:101) ε ◦ ( I + u )) d x − (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x )] h ( x ) · D ( (cid:101) ε ◦ ( I + u )) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:48) η (cid:107) h (cid:107) C (Ω) d . (3.14)Indeed, since (cid:101) ε ∈ C k,αS (cid:0) Ω (cid:1) ⊂ SBV(Ω), (cid:101) ε ◦ ( I + u ) ∈ SBV(Ω) and we can write thefollowing decomposition of D ( (cid:101) ε ◦ ( I + u )) into a continuous part and a jump part on arectifiable surface S : D ( (cid:101) ε ◦ ( I + u )) = ∇ ( (cid:101) ε ◦ ( I + u )) H d + [ (cid:101) ε ◦ ( I + u )] ν S H d − S ,
11e have that (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω (cid:2) g η ( x ) − [ (cid:101) ε ( x + u ) − ε ( x )] (cid:3) h ( x ) · ∇ ( (cid:101) ε ◦ ( I + u )) ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ C η (cid:107) h (cid:107) C (Ω) d . For the jump part, since S is a rectifiable surface and the function f η = g η − [ (cid:101) ε ◦ ( I + u ) − ε ]is piecewise continuous, it is possible to define a trace f η | S on the surface S satisfying (cid:107) f η | S (cid:107) L ( S ) ≤ C (cid:107) f η (cid:107) L (Ω) for some positive constant C depending only on S and Ω. Then we get (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) S f η h ( x ) · [ (cid:101) ε ◦ ( I + u )] ν S H d − S (cid:12)(cid:12)(cid:12)(cid:12) ≤ C η (cid:107) h (cid:107) C (Ω) d for some positive constant C independent of η and h .Now, the last term (cid:82) Ω [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] can be handled using Lemma 3.8.Doing so, we obtain (cid:90) Ω [ (cid:101) ε ( x + u + h ) − (cid:101) ε ( x + u )] = (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) + o ( (cid:107) h (cid:107) C (Ω) d ) . (3.15)Combining (3.12), (3.13), (3.14), and (3.15), we get that for every η > (cid:12)(cid:12)(cid:12)(cid:12) I ( u + h ) − I ( u ) − (cid:90) Ω [ (cid:101) ε ( x + u ) − ε ( x )] h ( x ) · D (cid:101) ε ◦ ( I + u )( x ) d x − (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x (cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:107) h (cid:107) C (Ω) d (cid:18) ρ ( (cid:107) h (cid:107) C (Ω) d ) + η (cid:19) for some positive constant C independent of h and η .Finally, it follows that I ( u + h ) − I ( u ) = ( ξ , h ) + (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x + o ( (cid:107) h (cid:107) C (Ω) d ) , where ξ is defined by (3.8). Hence, either (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x is of order of (cid:107) h (cid:107) C (Ω) d and we get I ( u + h ) − I ( u ) ≥ ( ξ , h )for (cid:107) h (cid:107) C (Ω) d small enough or (cid:90) Ω (cid:101) ε ( x + u ) | h · ν | δ ∂ ˜Ω ( x + u ) d x = o ( (cid:107) h (cid:107) C (Ω) d ) and in thiscase, I is Fr´echet differentiable and ξ is its Fr´echet derivative. The proof of Proposition3.5 is then complete. (cid:3) Remark 3.9
The minimization of the functional I gives a reconstruction of u ∗ on asubdomain Ω ⊂ Ω . In practical conditions, since u ∗ is small Ω is almost the wholedomain Ω . The values of u ∗ on the boundary are known and, since u ∗ is of class C , itis possible to deduce the values of u ∗ on Ω \ Ω by interpolation. Reconstruction of the shear modulus
The problem is now to recover the function µ the reconstructed internal data u . For doingso, we use the method described in [4]. We introduce the operator F u = F [ µ ] = ∇ · (cid:0) µ ( ∇ u + ∇ u T ) (cid:1) + ∇ p = 0 in Ω , ∇ · u = 0 in Ω , u = f on ∂ Ω , and minimize the function K given by C , (Ω ) −→ R µ (cid:55)−→ K [ µ ] = (cid:90) Ω |F [ µ ] − u | d x . According to [4], K is Fr´echet differentiable and its gradient can be explicitly computed.Let v be the solution of ∇ · (cid:0) µ ( ∇ v + ∇ v T ) (cid:1) + ∇ q = ( K [ µ ] − u ) in Ω , ∇ · v = 0 in Ω , v = 0 on ∂ Ω . Then, ∇K ( µ )[ h ] = (cid:90) Ω h ( ∇ v + ∇ v T ) : ( ∇ u + ∇ u T ) d x . A gradient descent method can be applied in order to reconstruct µ from u . We refer to[4] for more details. We take Ω = [0 , and discretize it on a 300 ×
300 grid, and generate a random Gaussianprocess to model the optical index ε of the medium as shown in Figure 5.1. Given a shearmodulus µ map on Ω; see Figure 5 (left), we solve (1.1) on Ω via a finite element methodcompute the displacement field u . We then compute the displaced optical index ε u byusing a spline interpolation approach and proceed to recover the shear modulus from thedata ε and ε u on the grid by the method described in the paper.Using (3.5), we first compute the initial guess u δ for the displacement field as the least-square solution to minimization of J x . Figure 5.2 shows the kernel w δ used to compute u δ . As one can see δ needs to be large enough so the matrix w δ (cid:63) (cid:0) ∇ ε ∇ ε T (cid:1) is invertible ateach point x , which is basically saying that δ must be bigger than the correlation lengthof ε . Figure 5.3 shows the conditioning of the matrix w δ (cid:63) (cid:0) ∇ ε ∇ ε T (cid:1) . Figure 5.4 shows thetrue displacement u ∗ , the result of the first order approximation (i.e., the initial guess) u δ and then the result of the optimization process using a gradient descent method tominimize the discrepancy functional I .Once the displacement inside the domain is reconstructed, we can recover the shearmodulus µ , as shown in Figure 5. We reconstruct µ by minimizing the functional K andusing a gradient descent-type method. Note that gradient of K is computed with theadjoint state method, described previously. As it can be seen in Figure 5, the reconstruc-tion is very accurate but not so perfect on the boundaries of Ω, which is due to the poorestimation of u on ∂ Ω. 13 . . . . . . . . . . . . . . . . . .
91 101112131415161718
Figure 5.1: Optical index ε of the medium. . . . . . . . . . . . . . . . . . .
91 00 . . . . . . . . . Figure 5.2: Averaging kernel w δ . In this paper, we developed a novel algorithm which gives access not only to stiffnessquantitative information of biological tissues but also opens the way to other contrastssuch as mechanical anisotropy. In the heart, the muscle fibers have anisotropic mechanicalproperties. It would be very interesting to detect a change in fiber orientation using OCTelastographic tomography.
References [1] G. Alberti and C. Mantegazza, A note on the theory of SBV functions, Boll. Un.Mat. Ital., B 11 (1997), 375–382.[2] L. Ambrosio, N. Fusco, and D. Pallara,
Functions of Bounded Variation and FreeDiscontinuity Problems , Clarendon Press Oxford, 2000.14 . . . . . . . . . . . . . . . . . .
91 1 . . . . . . . . Figure 5.3: Conditioning of the matrix w δ (cid:63) ∇ ε ∇ ε T .0 0 . . . . . . . . u ∗ · e · − · − . . . . . . . . u ∗ · e · − · − . . . . . . . .
81 Initial guess u δ · e · − · − . . . . . . . .
81 Initial guess u δ · e · − · − . . . . . . . .
81 Reconstructed u · e · − · − . . . . . . . .
81 Reconstructed u · e · − · − Figure 5.4: Displacement field and its reconstruction.15 . . . . . . . . . . . . . . . . . .
91 Shear modulus distribution µ . . . .
55 0 0 . . . . . . . . . . . . . . . . . .
91 Reconstructed shear modulus distribution µ rec . . . . Figure 5.5: Shear modulus reconstruction.[3] H. Ammari,
An Introduction to Mathematics of Emerging Biomedical Imaging ,Math. Appl., Vol. 62, Springer-Verlag, Berlin, 2008.[4] H. Ammari, E. Bretin, J. Garnier, H. Kang, H. Lee, and A. Wahab,
MathematicalMethods in Elasticity Imaging , Princeton Series in Applied Mathematics, PrincetonUniversity Press, 2014.[5] H. Ammari, P. Garapon, H. Kang, and H. Lee, A method of biological tissues elas-ticity reconstruction using magnetic resonance elastography measurements, Quart.Appl. Math., 66 (2008), 139–175.[6] J.F. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal.Mach. Intell., 8 (1986), 679–697.[7] Y.Z. Chen and L.C. Wu,
Second Order Elliptic Equations and Elliptic Systems ,Translated from the 1991 Chinese original by Bei Hu. Translations of MathematicalMonographs, 174. American Mathematical Society, Providence, RI, 1998.[8] F.H. Clarke, Yu. S. Ledyaev, R.J. Stern, and P.R. Wolenski,
Nonsmooth Analysisand Control Theory , Graduate Texts in Mathematics, Springer-Verlag, New York,1998.[9] A. Dubois and A.C. Boccara, Full-field optical coherence tomography, in
OpticalCoherence Tomography , 565–591 , Biological and Medical Physics, Biomedical En-gineering, Springer, 2008.[10] A. Dubois, K. Grieve, G. Moneron, R. Lecaque, L. Vabre, and C. Boccara, Ultrahigh-resolution full-field optical coherence tomography, Appl. Optics, 43 (2004), 2874–2883.[11] P. Elbau, L. Mindrinos, and O. Scherzer, Mathematical modeling of optical coherencetomography, arXiv: 1403.0726. 1612] M. Giaquinta and L. Martinazzi,
An Introduction to the Regularity Theory for Ellip-tic Systems, Harmonic Maps and Minimal Graphs , Second edition. Appunti. ScuolaNormale Superiore di Pisa (Nuova Serie), 11. Edizioni della Normale, Pisa, 2012.[13] Y.Y. Li and L. Nirenberg, Estimates for elliptic systems from composite material.Dedicated to the memory of J¨urgen K. Moser, Comm. Pure Appl. Math., 56 (2003),892–925.[14] X. Liang, V. Crecea, and S. Boppart, Dynamic optical coherence elastography: Areview, J. Innov. Opt. Health Sci., 3 (2010), 221–233.[15] A. Manduca, T.E. Oliphant, M.A. Dresner, J.L. Mahowald, S.A. Kruse, E. Amromin,J.P. Felmlee, J.F. Greenleaf, and R.L. Ehman, Magnetic resonance elastography:Non-invasive mapping of tissue elasticity, Med. Imag. Anal., 5 (2001), 237–254.[16] R. Muthupillai and R.L. Ehman, Magnetic resonance elastography, Nat. Med., 2(1996), 601–603.[17] W. Naetar and O. Scherzer, Quantitative photoacoustic tomography with piecewiseconstant material parameters, arXiv:1403.2620.[18] A. Nahas, M. Bauer, S. Roux, and A.C. Boccara, 3D static elastography at themicrometer scale using Full Field OCT, Biomedical Opt. Expr., 4 (2013), 2138–2149.[19] M. Razami, A. Mariampillai, C. Sun, V.X.D. Yang, and M.C. Kolios, Biomechanicalproperties of soft tissue measurement using optical coherence elastography, Proc.SPIE, 8207 (2012), 820758.[20] J. Rogowska, N.A. Patel, J.G. Fujimoto, and M.E. Brezinski, Optical coherence tomo-graphic elastography technique for measuring deformation and strain of atheroscle-rotic tissues, Heart, 90 (2004), 556–562.[21] J.M. Schmitt, OCT elastography: imaging microscopic deformation and strain intissue, Opt. Express, 3 (1998), 199–211.[22] J.K. Seo and E.J. Woo,