Image reconstruction through metamorphosis
IImage reconstruction through metamorphosis
Barbara Gris ∗ Chong Chen † Ozan ¨Oktem ‡ Abstract
This article adapts the framework of metamorphosis to solve inverseproblems in imaging that includes joint reconstruction and image reg-istration. The deformations in question have two components, one thatis a geometric deformation moving intensities and the other a defor-mation of intensity values itself, which, e.g., allows for appearance ofa new structure. The idea developed here is to reconstruct an imagefrom noisy and indirect observations by registering, via metamorpho-sis, a template to the observed data. Unlike a registration with onlygeometrical changes, this framework gives good results when intensi-ties of the template are poorly chosen. We show that this method isa well-defined regularisation method (proving existence, stability andconvergence) and present several numerical examples. In shape based reconstruction or spatiotemporal image reconstruction , a keydifficulty is to match an image against an indirectly observed target ( indirectimage registration ). This paper provides theory and algorithms for indirectimage registration applicable to general inverse problems. Before proceed-ing, we give a brief overview of these notions along with a short survey ofexisting results. Shape based reconstruction
The goal is to recover shapes of interiorsub-structures of an object whereas variations within these is of less im-portance. Examples of such imaging studies are nano-characterisation ofspecimens by means of electron microscopy or x-ray phase contrast imag-ing, e.g., nano-characterisation of materials by electron electron tomography(ET) primarily focuses on the morphology of sub-structures [5]. Another ex-ample is quantification of sub-resolution porosity in materials by means ofx-ray phase contrast imaging. ∗ LJLL - Laboratoire Jacques-Louis Lions, UPMC, Paris, France. † LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academyof Sciences, Beijing 100190, People’s Republic of China ‡ Department of Mathematics, KTH – Royal Institute of Technology, Stockholm, Swe-den. a r X i v : . [ c s . C V ] J un n these imaging applications it makes sense to account for qualitativeprior shape information during the reconstruction. Enforcing an exact spa-tial match between a template and the reconstruction is often too strongsince realistic shape information is almost always approximate, so the natu-ral approach is to perform reconstruction assuming the structures are ‘shapewise similar’ to a template. Spatiotemporal imaging
Imaging an object that undergoes temporalvariation leads to a spatiotemporal reconstruction problem where both theobject and its time variation needs to be recovered from noisy time seriesof measured data. An important case is when the only time dependency isthat of the object.Spatiotemporal imaging occurs in medical imaging, see, e.g., [23] for asurvey of organ motion models. It is particular relevant for techniques likepositron emission tomography (PET) and single photon emission computedtomography (SPECT), which are used for visualising the distribution of in-jected radiopharmaceuticals (activity map). The latter is an inherently dy-namic quantity, e.g., anatomical structures undergo motion, like the motionof the heart and respiratory motion of the lungs and thoracic wall, duringthe data acquisition. Not accounting for organ motion is known to degradethe spatial localisation of the radiotracer, leading to spatially blurred im-ages. Furthermore, even when organ motion can be neglected, there areother dynamic processes, such as the uptake and wash-out of radiotracersfrom body organs. Visualising such kinetics of the radiotracers can actuallybe a goal in itself, as in pre-clinical imaging studies related to drug discov-ery/development. The term ‘dynamic’ in PET and SPECT imaging oftenrefers to such temporal variation due to radiotracers kinetics rather thanorgan movement [13].To exemplify the above mentioned issues, consider SPECT based cardiacperfusion studies and F-FDG-PET imaging of lung nodules/tumours. Theformer needs to account for the beating heart and the latter for respiratorymotion of the lungs and thoracic wall. Studies show a maximal displacementof 23 mm (average 15–20 mm) due to respiratory motion [21] and 42 mm(average 8–23 mm) due to cardiac motion in thoracic PET [25].
Indirect image registration (matching)
In image registration the aimis to deform a template image so that it matches a target image, whichbecomes challenging when the template is allowed to undergo non-rigid de-formations.A well developed framework is diffeomorphic image registration wherethe image registration is recast as the problem of finding a suitable diffeo-morphism that deforms the template into the target image [26, 2]. Theunderlying assumption is that the target image is contained in the orbit of2he template under the group action of diffeomorphisms. This can be statedin a very general setting where diffeomorphisms act on various structures,like landmark points, curves, surfaces, scalar images, or even vector/tensorvalued images.The registration problem becomes more challenging when the target isonly known indirectly through measured data. This is referred to as indirectimage registration , see [19] for using small diffeomorphic deformations and[9, 14] for adapting the large deformation diffeomorphic metric mapping(LDDMM) framework to indirect image registration.
The paper adapts the metamorphosis framework [24] to the indirect im-age registration setting. Metamorphosis is an extension of the LDDMMframework (diffeomorphometry) [26, 16] where not only the geometry of thetemplate, but also the grey-scale values undergo diffeomorphic changes.We start by recalling necessary theory from LDDMM-based indirect reg-istration (section 3). Using the notions from section 3, we adapt the meta-morphosis framework to the indirect setting (section 4). We show how thisframework allows to define a regularization method for inverse problems,satisfying properties of existence, stability and convergence (section 4.3).The numerical implementation is outlined in section 4.4. We present severalnumerical examples from 2D tomography, and in particular give a prelimi-nary result for motion reconstruction when the acquisition is done at severaltime points. We also study the robustness of our methods with respects tothe parameters (section 5).
We recall here the notion of large diffeomorphic deformations defined byflows of time-varying vector fields, as formalized in [1].Let Ω ⊂ R d be a fixed bounded domain and let X := L (Ω , R ) representgrey scale images on Ω. Next, let V denote a fixed Hilbert space of vectorfields on R d . We will assume V ⊂ C p (Ω), i.e., the vector fields are supportedon Ω and p times continuously differentiable. Finally, L ([0 , , V ) denotesthe space of time-dependent V -vector fields that are integrable, i.e., ν ( t, · ) ∈ V and t (cid:55)→ (cid:13)(cid:13) ν ( t, · ) (cid:13)(cid:13) C p is integrable on [0 , (cid:107) ν (cid:107) p := (cid:16)(cid:90) (cid:13)(cid:13) ν ( t, · ) (cid:13)(cid:13) pV d t (cid:17) /p (cid:107) · (cid:107) V is the naturally defined norm based upon the inner product ofthe Hilbert space V of vector fields.The following proposition allows one to consider flows of elements in L ([0 , , V ) and ensures that these flows belong to Diff p (Ω) (set of p -diffeomorphisms that are supported in Ω ⊂ R d , and if Ω is unbounded,tend to zero towards infinity). Proposition 1.
Let ν ∈ L ([0 , , V ) and consider the ordinary differentialequation (flow equation): dd t φ ( t, x ) = ν (cid:0) t, φ ( t, x ) (cid:1) φ (0 , x ) = x for any x ∈ Ω and t ∈ [0 , . (1) Then, (1) has a unique absolutely continuous solution φ ( t, · ) ∈ Diff p ( R d ) . The above result is proved in [1] and the unique solution of (1) is hence-forth called the flow of ν . We also introduce to notation ϕ ν s,t : R d → R d thatrefers to ϕ ν s,t := φ ( t, · ) ◦ φ ( s, · ) − for s, t ∈ [0 ,
1] (2)where φ : Ω → R d denotes the unique solution to (1).As stated next, the set of diffeomorphisms that are given as flows formsa group that is a complete metric space [1]. Proposition 2.
Let V ⊂ C p (Ω) ( p ≥ ) be an admissible reproducing kernelHilbert space (RKHS) and define G V := (cid:110) φ : R d → R d | φ = ϕ ν , for some ν ∈ L ([0 , , V ) (cid:111) . Then G V forms a sub-group of Diff p ( R d ) that is a complete metric spaceunder the metric d G ( φ , φ ) := inf (cid:110) (cid:107) ν (cid:107) : ν ∈ L ([0 , , V ) and φ = φ ◦ ϕ ν , (cid:111) = inf (cid:110) (cid:107) ν (cid:107) : ν ∈ L ([0 , , V ) and φ = φ ◦ ϕ ν , (cid:111) . The elements of G V are called large diffeomorphic deformations and G V acts on X via the geometric group action that is defined by the operator W : G V × X → X where W ( φ, I ) := I ◦ φ − . (3)We conclude by stating regularity properties of flows of velocity fields aswell as the group action in (3), these will play an important role in what isto follow. The proof is given in [6]. 4 roposition 3. Assume V ⊂ C p (Ω) ( p ≥ ) is a fixed admissible Hilbertspace of vector fields on Ω and { ν n } n ⊂ L ([0 , , V ) a sequence that con-verges weakly to ν ∈ L ([0 , , V ) . Then, the following holds with ϕ nt := ϕ ν n ,t :1. ( ϕ nt ) − converges to ( ϕ ν ,t ) − uniformly w.r.t. t ∈ [0 , and uniformlyon compact subsets of Ω ⊂ R d .2. lim n →∞ (cid:13)(cid:13)(cid:13) W ( ϕ nt , I ) − W ( ϕ ν ,t , I ) (cid:13)(cid:13)(cid:13) X = 0 for any f ∈ X . Image registration (matching) refers to the task of deforming a given tem-plate image I ∈ X so that it matches a given target image I ∗ ∈ X .The above task can also be stated in an indirect setting, which refersto the case when the template I ∈ X is to be registered against a target I ∗ ∈ X that is only indirectly known through data g ∈ Y where g = A ( I ∗ ) + e. (4)In the above, A : X → Y (forward operator) is known and assumed to bedifferentiable and e ∈ Y is a single sample of a Y -valued random elementthat denotes the measurement noise in the data.A further development requires specifying what is meant by deforminga template image, and we will henceforth consider diffeomorphic (non-rigid)deformations, i.e., diffeomorphisms that deform images by actin g on themthrough a group action. LDDMM-based registration
An example of using large diffeomorphic(non-rigid) deformations for image registration is to minimize the followingfunctional: G V (cid:51) φ (cid:55)→ γ G (Id , φ ) + (cid:13)(cid:13) W ( φ, I ) − I ∗ (cid:13)(cid:13) X given γ > V is admissible, then minimizing the above functional on G V amounts tominimizing the following functional on L ([0 , , V ) [26, Theorem 11.2 andLemma 11.3]: L ([0 , , V ) (cid:51) ν (cid:55)→ γ (cid:107) ν (cid:107) + (cid:13)(cid:13) W ( ϕ ν , , I ) − f (cid:13)(cid:13) X given γ > L ([0 , , V ) is a vector space,whereas G V is not, so it is easier to minimize a functional over L ([0 , , V )rather than over G V .The above can be extended to the indirect setting as shown in [9], whichwe henceforth refer to as LDDMM-based indirect registration . More pre-cisely, the corresponding indirect registration problem can be adressed by5inimising the functional L ([0 , , V ) (cid:51) ν (cid:55)→ γ (cid:107) ν (cid:107) + L (cid:0) ( A ◦ W )( ϕ ν ,t , I ) , g (cid:1) . Here, L : Y × Y → R is typically given by an appropriate affine transformof the data negative log-likelihood [4], so minimizing f (cid:55)→ L (cid:0) A ( f ) , g (cid:1) corre-sponds to seeking a maximum likelihood solution of (4).An interpretation of the above is that the template image I , which isassumed to be given a priori, acts as a shape prior when solving the inverseproblem in (4) and γ > min ν ∈ L ([0 , ,V ) (cid:20) γ (cid:107) ν (cid:107) + L (cid:16) ( A ◦ W ) (cid:0) φ (1 , · ) , I (cid:1) , g (cid:17)(cid:21) ddt φ ( t, x ) = ν (cid:0) t, φ ( t, x ) (cid:1) ( t, x ) ∈ Ω × [0 , ,φ (0 , x ) = x x ∈ Ω . (5) As shown in [9], access to a template that can act as a shape prior canhave profound effect in solving challenging inverse problem in imaging. Asan example, tomographic imaging problems that are otherwise intractable(highly noisy and sparsely sampled data) can be successfully addressed usingindirect registration even when using a template is far from the target imageused for generating the data.When template has correct topology and intensity levels, then LDDMM-based indirect registration with geometric group action is remarkably stableas shown in [9]. Using a geometric group action, however, makes it impos-sible to create or remove intensity, e.g., it is not possible to start out from atemplate with a single isolated structure and deform it to a image with twoisolated structures. This severely limits the usefulness of LDDMM-based in-direct registration, e.g., spatiotemporal images (moves) are likely to involvechanges in both geometry (objects appear or disappear) and intensity. Seefig. 1 for an example of how wrong intensity influences the registration.As noted in [9], one approach is to replace the geometric group actionwith one that alters intensities, e.g., a mass preserving group action. An-other is to keep the geometric group action, but replace LDDMM witha framework for diffeomorphic deformations that acts on both geometryand intensities, e.g., metamorphosis. This latter approach is the essence ofmetamorphosis-based indirect registration.6 a) Template. (b) Target. (c) Reconstruction. (d) Data.
Figure 1: Reconstruction by LDDMM-based indirect registration (c) usinga template (a) with a geometry that matches the target (b), but with in-correct background intensity values. Target is observed indirectly throughtomographic data (d), which is 2D parallel beam Radon transform with100 evenly distributed directions (see section 5.1 for details). The artefactsin the reconstruction are due to incorrect background intensity in template.
In metamorphosis diffeomorphisms are still generated by flows as in LD-DMM, but the difference is that they now act with a geometric group actionon both intensities and underlying points. As such, metamorphosis extendsLDDMM. The abstract definition of a metamorphosis reads as follows.
Definition 1 (Metamorphosis [24]) . Let V ⊂ C p (Ω) be an admissible Hilbertspace and “ . ” denotes some group action of G V on X . A Metamorphosis is a curve t (cid:55)→ ( φ t , J t ) in G V × X . The curve t (cid:55)→ f t := φ t .J t is called the image part, t (cid:55)→ φ t is the deformation part, and t (cid:55)→ f t is the template part. The image part represents the temporal evolution that is not relatedto intensity changes, i.e., evolution of underlying geometry, whereas thetemplate part is the evolution of the intensity. Both evolutions, which arecombined in metamorphosis, are driven by the same underlying flow of dif-feomorphisms in G V .A important case is when the metamorphosis t (cid:55)→ ( φ t , f t ) has a defor-mation part that solves the flow equation (1) and a template part is C intime. More precisely, L ([0 , , X ) denotes the space of functions in X thatare square integrable, i.e., ζ ( t, · ) ∈ X and t (cid:55)→ (cid:13)(cid:13) ζ ( t, · ) (cid:13)(cid:13) X is in L ([0 , , R ).The norm on L ([0 , , X ) is then (cid:107) ζ (cid:107) := (cid:16)(cid:90) (cid:13)(cid:13) ζ ( t, · ) (cid:13)(cid:13) X d t (cid:17) / . We will also use the notation L ([0 , , V × X ) := L ([0 , , V ) × L ([0 , , X ) . ν , ζ ) ∈ L ([0 , , V × X ) and I ∈ X , define the curve t (cid:55)→ I ν ,ζt , which is absolutely continuous on [0 , dd t I ν ,ζt ( x ) = ζ (cid:0) t, ϕ ν ,t ( x ) (cid:1) I ν ,ζ ( x ) = I ( x ) with ϕ ν ,t ∈ G V as in (2). (6)The metamorphosis can now be parametrised as t (cid:55)→ ( ϕ ν ,t , I ν ,ζt ). Indirect registration
The indirect registration problem in section 3.2 canbe approached by metamorphosis instead of LDDMM. Similar to LDDMM-based indirect image registration in [9], we define metamorphosis-based in-direct image registration as the minimization of the objective functional J γ,τ ( · ; g ) : L ([0 , , V × X ) → R defined as J γ,τ ( ν , ζ ; g ) := γ (cid:107) ν (cid:107) + τ (cid:107) ζ (cid:107) + L (cid:16) A (cid:0) W ( ϕ ν , , I ν ,ζ ) (cid:1) , g (cid:17) (7)for given regularization parameters γ, τ ≥
0, measured data g ∈ Y , andinitial template I ∈ X that sets the initial condition I ν ,ζ ( x ) := I ( x ).Hence, performing metamorphosis-based indirect image registration of atemplate I against a target indirectly observed through data g amounts tosolving ( (cid:98) ν , (cid:98) ζ ) ∈ arg min ( ν ,ζ ) J γ,τ ( ν , ζ ; g ) . (8)The above always has a solution assuming the data discrepancy and theforward operator fulfills some weak requirements (see proposition 4). Froma solution we then obtain the following: • Initial template: I ∈ X such that I ν ,ζ := I . • Reconstruction:
Final registered template f (cid:98) ν , (cid:98) ζ = W (cid:0) ϕ (cid:98) ν , , I (cid:98) ν , (cid:98) ζ (cid:1) ∈ X . • Image trajectory:
The evolution of both geometry and intensity of thetemplate, given by t (cid:55)→ W (cid:0) ϕ (cid:98) ν ,t , I (cid:98) ν , (cid:98) ζt (cid:1) . • Template trajectory:
The evolution of intensities of the template, i.e.,the part that does not include evolution of geometry: t (cid:55)→ I (cid:98) ν , (cid:98) ζt . • Deformation trajectory:
The geometric evolution of the template, i.e.,the part that does not include evolution of intensity: t (cid:55)→ W ( ϕ (cid:98) ν ,t , I ).8 .3 Regularising properties In the following we prove several properties (existence, stability and conver-gence) of metamorphosis-based indirect image registration, which are nec-essary if the approach is to constitute a well defined regularisation method (notion defined in [12]). We set X := L (Ω , R ) and Y a Hilbert space. Proposition 4 (Existence) . Assume A : X → Y is continuous and the datadiscrepancy L ( · , g ) : Y → R is weakly lower semi-continuous for any g ∈ Y .Then, J γ,τ ( · , g ) : L ([0 , , V × X ) → R defined through (6) and (7) has aminimizer in L ([0 , , V × X ) for any I ∈ L (Ω , R ) .Proof. We follow here the strategy to prove existence of minimal trajecto-ries for metamorphosis (as in [8] for instance). One considers a minimiz-ing sequence of J γ,τ ( · ; g ), i.e., a sequence that converges to the infimumof J γ,τ ( · ; g ) (such a sequence always exists). The idea is to prove thatsuch a minimizing sequence has a sub-sequence that converges to a point in L ([0 , , V × X ), i.e., the infimum is contained in L ([0 , , V × X ) whichproves existence of a minima.Bearing in mind the above, we start by considering a minimizing se-quence (cid:8) ( ν n , ζ n ) (cid:9) n ⊂ L ([0 , , V × X ) to J γ,τ ( · ; g ), i.e.,lim n →∞ J γ,τ ( ν n , ζ n ; g ) = inf ν ,ζ J γ,τ ( ν , ζ ; g ) . Since (cid:8) ν n (cid:9) n ⊂ L ([0 , , V ) is bounded, it has a sub-sequence that convergesto an element ν ∞ ∈ L ([0 , , V ). Likewise, (cid:8) ζ n (cid:9) n ⊂ L ([0 , , X ) has a sub-sequence that converges to an element ζ ∞ ∈ L ([0 , , X ). Hence, with aslight abuse of notation, we conclude that ν n (cid:42) ν ∞ and ζ n (cid:42) ζ ∞ as n → ∞ .The aim is now to prove existence of minimizers by showing that ( ν ∞ , ζ ∞ )is a minimizer to J γ,τ ( · ; g ) : L ([0 , , V × X ) → R .Before proceeding, we introduce some notation in order to simplify theexpressions. Define I nt := I ν n ,ζ n t and ϕ ns,t := ϕ ν n s,t for n ∈ N (cid:91) {∞} . (9)Hence, assuming geometric group action (3) and using (2), we can write J γ,τ ( ν n , ζ n ; g ) = γ (cid:107) ν n (cid:107) + τ (cid:107) ζ n (cid:107) + L (cid:16) A (cid:0) I n ◦ ϕ n , (cid:1) , g (cid:17) for n ∈ N (cid:83) {∞} . Assume next that the following holds: I n ◦ ϕ n , (cid:42) I ∞ ◦ ϕ ∞ , as n → ∞ . (10)9he data discrepancy term L ( · , g ) : Y → R is weakly lower semi continuousand the forward operator A : X → Y is continuous, so L ( · , g ) ◦ A is alsoweakly lower semi continuous and then (10) implies L (cid:0) A ( I ∞ ◦ ϕ ∞ , ) , g (cid:1) ≤ lim inf n →∞ L ( A ( I n ◦ ϕ n , ) , g ) . (11)Furthermore, from the weak convergences of ν n and ζ n , we get γ (cid:107) ν ∞ (cid:107) + τ (cid:107) ζ ∞ (cid:107) ≤ lim inf n →∞ (cid:104) γ (cid:107) ν n (cid:107) + τ (cid:107) ζ n (cid:107) (cid:105) . (12)Hence, combining (11) and (12) we obtain J γ,τ ( ν ∞ , ζ ∞ ; g ) ≤ lim n →∞ J γ,τ ( ν n , ζ n ; g ) . Since (cid:8) ( ν n , ζ n ) (cid:9) n ⊂ L ([0 , , V × X ) is a minimizing sequence, this yields J γ,τ ( ν ∞ , ζ ∞ ; g ) = inf ( ν ,ζ ) ∈ L ([0 , ,V × X ) J γ,τ ( ν , ζ ; g ) , which proves ( ν ∞ , ζ ∞ ) ∈ L ([0 , , V × X ) is a minimizer to J γ,τ ( · ; g ).Hence, to finalize the proof we need to show that (10) holds. We startby observing that the solution of (6) can be written as I nt := I n ( x ) + (cid:90) t ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) d s for n ∈ N ∪ {∞} , (13)and note that ( t, x ) (cid:55)→ I nt ( x ) ∈ C ([0 , × Ω , R ). Next, we claim that I n (cid:42) I ∞ for some I ∞ ∈ X ,which is equivalent tolim n →∞ (cid:104) I n − I ∞ , J (cid:105) = 0 for any J ∈ L (Ω , R ). (14)To prove (14), note first that since continuous functions are dense in L , itis enough to show (14) holds for J ∈ C (Ω , R ). Next, (cid:104) I n − I ∞ , J (cid:105) = (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) − ζ ∞ (cid:0) s, ϕ ∞ ,s ( x ) (cid:1)(cid:17) J ( x )d s d x (15)= (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) − ζ n (cid:0) s, ϕ ∞ ,s ( x ) (cid:1)(cid:17) J ( x )d s d x (16)+ (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) − ζ ∞ (cid:0) s, ϕ ∞ ,s ( x ) (cid:1)(cid:17) J ( x )d s d x. (17)10et us now take a closer look at the term in (16): (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) − ζ n (cid:0) s, ϕ ∞ ,s ( x ) (cid:1)(cid:17) J ( x )d s d x = (cid:90) Ω (cid:90) t ζ n ( s, x ) J (cid:0) ϕ n ,s ( x ) (cid:1)(cid:12)(cid:12) Dϕ n ,s ( x ) (cid:12)(cid:12) d s d x − (cid:90) Ω (cid:90) t ζ ∞ ( s, x ) J (cid:0) ϕ ∞ ,s ( x ) (cid:1)(cid:12)(cid:12) Dϕ ∞ ,s ( x ) (cid:12)(cid:12) d s d x = (cid:90) Ω (cid:90) t ζ n ( s, x ) (cid:16) J (cid:0) ϕ n ,s ( x ) (cid:1)(cid:12)(cid:12) Dϕ n ,s ( x ) (cid:12)(cid:12) − J (cid:0) ϕ ∞ ,s ( x ) (cid:1)(cid:12)(cid:12) Dϕ ∞ ,s ( x ) (cid:12)(cid:12)(cid:17) d s d x − (cid:90) Ω (cid:90) t (cid:16) ζ ∞ ( s, x ) − ζ n ( s, x ) (cid:17) J (cid:0) ϕ ∞ ,s ( x ) (cid:1)(cid:12)(cid:12) Dϕ ∞ ,s ( x ) (cid:12)(cid:12) d s d x = (cid:104) ζ n , J n − J ∞ (cid:105) − (cid:104) ζ ∞ − ζ n , J ∞ (cid:105) where J n ∈ L ([0 , , X ) is defined as J n ( s, x ) := J (cid:0) ϕ ns, ( x ) (cid:1)(cid:12)(cid:12) Dϕ ns, ( x ) (cid:12)(cid:12) for n ∈ N (cid:91) {∞} . (18)By proposition 3 we know that ϕ ns, → ϕ ∞ s, and Dϕ ns, → Dϕ ∞ s, uniformlyon Ω. Since J is continuous on Ω, we conclude that (cid:107) J n − J ∞ (cid:107) tends to 0.Since ζ n is bounded, we conclude that (cid:104) ζ n , J n − J ∞ (cid:105) ≤ (cid:107) ζ n (cid:107) · (cid:107) J n − J ∞ (cid:107) → . Furthermore, since ζ n (cid:42) ζ ∞ , we also get (cid:104) ζ ∞ − ζ n , J ∞ (cid:105) →
0. Hence, wehave shown that (16) tends to zero, i.e.,lim n →∞ (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ n ,s ( x ) (cid:1) − ζ n (cid:0) s, ϕ ∞ ,s ( x ) (cid:1)(cid:17) J ( x )d s d x = 0 . Finally, we consider the term in (17). Since ζ n (cid:42) ζ ∞ , we immediatelyobtain (cid:90) Ω (cid:90) t (cid:16) ζ n (cid:0) s, ϕ ∞ s ( x ) (cid:1) − ζ ∞ (cid:0) s, ϕ ∞ s ( x ) (cid:1)(cid:17) J ( x )d s d x = (cid:10) ζ n − ζ ∞ , J ∞ (cid:11) → . To summarise, we have just proved that both terms (16) and (17) tend to 0as n → ∞ , which implies that (14) holds, i.e., I n (cid:42) I ∞ .To prove (10), i.e., I n ◦ ϕ n , (cid:42) I ∞ ◦ ϕ ∞ , , we need to show thatlim n →∞ (cid:10) I n ◦ ϕ n , − I ∞ ◦ ϕ ∞ , , J (cid:11) = 0 for any J ∈ L (Ω , R ), (19)and as before, we may assume J ∈ C (Ω , R ). Using (18) we can express theterm in (19) whose limit we seek as (cid:12)(cid:12) (cid:104) I n ◦ ϕ n , − I ∞ ◦ ϕ ∞ , , J (cid:105) (cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:10) I n , J n (1 , · ) − J ∞ (1 , · ) (cid:11)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:10) I n − I ∞ , J ∞ (1 , · ) (cid:11)(cid:12)(cid:12)(cid:12) ≤ (cid:107) I n (cid:107) · (cid:13)(cid:13) J n (1 , · ) − J ∞ (1 , · ) (cid:13)(cid:13) + (cid:12)(cid:12) (cid:104) I n − I ∞ , J ∞ (1 , · ) (cid:105) (cid:12)(cid:12) . (cid:107) I n (cid:107) is bounded (because (cid:107) ζ n (cid:107) is bounded) and since I n (cid:42) I ∞ (whichwe shoed before), all terms above tend to 0 as n → ∞ , i.e., (19) holds.This concludes the proof of (10), which in turn implies the existence ofa minimizer of J γ,τ ( · ; g ).Our next result shows that the solution to the indirect registration prob-lem is (weakly) continuous w.r.t. variations in data, and as such, it is a kindof stability result. Proposition 5 (Stability) . Let { g k } k ⊂ Y and assume this sequence con-verges (in norm) to some g ∈ Y . Next, for each γ, τ > and each k , define ( ν k , ζ k ) ∈ L ([0 , , V × X ) as ( ν k , ζ k ) = arg min ( ν ,ζ ) J γ,τ ( ν , ζ ; g k ) . Then there exists a sub sequence of ( ν k , ζ k ) that converges weakly to a min-imizer of J γ,τ ( · ; g ) in (7) .Proof. J γ,τ ( · ; g k ) has a minimizer ( ν k , ζ k ) ∈ L ([0 , , V × X ) for any g k ∈ Y (proposition 4). The idea is first to show that the sequences ( ν k ) k and( ζ k ) k are bounded. Next, we show that there exists a weakly convergingsubsequence of ( ν k , ζ k ) that converges to a minimizer ( ν , ζ ) of J γ,τ ( · ; g ),which also exists due to proposition 4.Since ( ν k , ζ k ) minimizes J γ,τ ( · ; g k ), by (7) we have (cid:107) ν k (cid:107) ≤ γ J γ,τ ( · ; g k )( ν k , ζ k ) ≤ γ J γ,τ ( · ; g k )( ,
0) for each k . (20)Observe now that if ν = and ζ = 0, then ϕ ν , = Id by (1) and I ν ,ζ = I by (6), so in particular W (cid:0) ϕ ν , , I ν ,ζ ) (cid:1) = I whenever ν = and ζ = 0,Hence, J γ,τ ( · ; g k )( ,
0) = L (cid:0) A ( I ) , g k (cid:1) and, in addition, (cid:107) ν (cid:107) = 0 and (cid:107) ζ (cid:107) = 0, so (20) becomes (cid:107) ν k (cid:107) ≤ γ L (cid:0) A ( I ) , g k (cid:1) → L ( A ( I ) , g ) as k → ∞ . (21)In conclusion, the sequence ( ν k ) k ⊂ L ([0 , , V ) is bounded. In a similarway, we can show that ( ζ k ) k ⊂ L ([0 , , X ) is bounded.The boundedness of both sequences implies that there are sub sequencesto these that converge weakly to some element ν ∞ ∈ L ([0 , , V ) and ζ ∞ ∈ L ([0 , , X ), respectively. Thus, to complete the proof, we need to showthat ( ν ∞ , ζ ∞ ) ∈ L ([0 , , V × X ) minimizes J γ,τ ( · ; g ), i.e., that J γ,τ ( ν ∞ , ζ ∞ ; g ) ≤ J γ,τ ( ν , ζ ; g ) holds for any ( ν , ζ ) ∈ L ([0 , , V × X ).12rom the weak convergences, we obtain γ (cid:107) ν ∞ (cid:107) + τ (cid:107) ζ ∞ (cid:107) ≤ γ k (cid:107) ν k (cid:107) + τ k (cid:107) ζ k (cid:107) ≤
12 lim inf k (cid:104) γ (cid:107) ν k (cid:107) + τ (cid:107) ζ k (cid:107) (cid:105) . (22)The weak convergence also implies (see proof of proposition 4) that W (cid:0) ϕ k , , I ∞ (cid:1) (cid:42) W (cid:0) ϕ ∞ , , I ∞ (cid:1) in X .In the above, we have used the notational convention introduced in (9). Bythe lower semi-continuity of L , we get L (cid:16) A (cid:0) W ( ϕ ∞ , , I ∞ ) (cid:1) , g (cid:17) ≤ lim inf k L (cid:16) A (cid:0) W ( ϕ k , , I k ) (cid:1) , g k (cid:17) . (23)Hence, J γ,τ ( ν ∞ , ζ ∞ ; g ) = γ (cid:107) ν ∞ (cid:107) + τ (cid:107) ζ ∞ (cid:107) + L (cid:16) A (cid:0) W ( ϕ ∞ , , I ∞ ) (cid:1) , g (cid:17) . ≤
12 lim inf k (cid:104) γ (cid:107) ν k (cid:107) + τ (cid:107) ζ k (cid:107) (cid:105) + lim inf k L (cid:16) A (cid:0) W ( ϕ k , , I k ) (cid:1) , g k (cid:17) ≤ lim inf k J γ,τ ( ν k , ζ k ; g k ) . (24)Next, since ( ν k , ζ k ) ∈ L ([0 , , V × X ) minimizes J γ,τ ( · ; g k ), we get J γ,τ ( ν ∞ , ζ ∞ ; g ) ≤ lim inf k J γ,τ ( ν , ζ ; g k )for any ( ν , ζ ) ∈ L ([0 , , V × X ). Furthermore, J γ,τ ( ν , ζ ; g k ) → J γ,τ ( ν , ζ ; g ) , so J γ,τ ( ν ∞ , ζ ∞ ; g ) ≤ J γ,τ ( ν , ζ ; g ) for all ( ν , ζ ) ∈ L ([0 , , V × X ).In particular, we have shown that ( ν ∞ , ζ ∞ ) minimises J γ,τ ( · ; g ).Our final results concerns convergence, which investigates the behaviourof the solution as data error tends to zero and regularization parameters areadapted accordingly through a parameter choice rule against the data error. Proposition 6 (Convergence) . Let g ∈ Y and assume A (cid:0) W ( ϕ ν , , I ν ,ζ ) (cid:1) = g for some ( ν , ζ ) ∈ L ([0 , , V × X ) .Next, for parameter choice rules δ (cid:55)→ γ ( δ ) and δ (cid:55)→ τ ( δ ) with δ > , define ( ν δ , ζ δ ) ∈ arg min ( ν ,ζ ) J γ ( δ ) ,τ ( δ ) ( ν , ζ ; g + e δ )13 here e δ ∈ Y (data error) has magnitude (cid:107) e δ (cid:107) = δ . Finally, assume that δ (cid:55)→ γ ( δ ) /τ ( δ ) and δ (cid:55)→ τ ( δ ) /γ ( δ ) are bounded, and lim δ → γ ( δ ) = lim δ → τ ( δ ) = lim δ → δ γ ( δ ) = lim δ → δ τ ( δ ) = 0 . Then, for any sequence δ k → there exists a sub-sequence δ k (cid:48) such that ( ν δ k (cid:48) , ζ δ k (cid:48) ) converges weakly to a ( ν ∗ , ζ ∗ ) satisfying A (cid:0) W ( ϕ ν ∗ , , I ν ∗ ,ζ ∗ ) (cid:1) = g .Proof. Let ( δ k ) be a sequence converging to 0 and, for each k , let us denote g k := g + e δ k , ν k := ν δ k , and ζ k := ζ δ k . Similarly to previous proofs, we will show that the sequences ( ν k ) and ( ζ k )are bounded, and then that the weakly converging subsequence that can beextracted from ( ν k , ζ k ) converges to a suitable solution.Define γ k := γ ( δ k ) and τ k := γ ( δ k ). Then, for each k we have (cid:107) ν k (cid:107) ≤ γ k J γ k ,τ k ,g k ( ν k , ζ k ) ≤ γ k J γ k ,τ k ,g k ( (cid:98) ν , ˆ ζ )= 1 γ k (cid:16) γ k (cid:107) (cid:98) ν (cid:107) + τ k (cid:107) ˆ ζ (cid:107) + L ( g, g k ) (cid:17) ≤ (cid:107) (cid:98) ν (cid:107) + τ k γ k (cid:107) ˆ ζ (cid:107) + δ k γ k . From the assumptions on the parameter choice rules, we conclude that( ν k ) ⊂ L ([0 , , V ) is bounded. Similarly, one can show that ( ζ k ) ⊂ L ([0 , , X ) is bounded.From the above, we conclude that there is a subsequence of ( ν k , ζ k ) thatconverges weakly to ( (cid:101) ν , (cid:101) ζ ) in L ([0 , , V ) × L ([0 , , V ). Then (see proofof proposition 4) L (cid:16) A (cid:0) W ( ϕ (cid:101) ν , , I (cid:101) ν , (cid:101) ζ ) (cid:1) , g (cid:17) ≤ lim inf k L (cid:16) A ( W (cid:0) ϕ ν k , , I (cid:101) ν , (cid:101) ζ ) (cid:1) , g k (cid:17) . Furthermore, the above quantity converges to 0 since L (cid:16) A (cid:0) W ( ϕ ν k , , I ν k ,ζ k ) (cid:1) , g k (cid:17) ≤ J γ k ,τ k ,g k ( ν k , ζ k ) ≤ J γ k ,τ k ,g k ( (cid:98) ν , ˆ ζ )= γ k (cid:107) (cid:98) ν (cid:107) + τ k (cid:107) ˆ ζ (cid:107) + L ( g, g k ) → k → ∞ .Hence, A (cid:0) W ( ϕ (cid:101) ν , , I (cid:101) ν , (cid:101) ζ ) (cid:1) = g . In order to solve (8), we use a gradient descent scheme on the variable( ν , ζ ) ∈ L ([0 , , V × X ) with a uniform discretization of the interval [0 , N parts, i.e., t i = 1 /N for i = 0 , . . . , N and the gradient descent is per-formed on ν ( t i , · ), ζ ( t i , · ), for i = 0 , , . . . , N . An alternative approach14eveloped in [18] extends the time discrete path method in [11] to the indi-rect setting.In order to compute numerical integrations, we use a Euler scheme onthis discretization. The flow equation (1) is computed using the followingapproximation with small deformations: ϕ ν t i , ≈ ϕ ν t i − , ◦ (cid:16) Id − N ν ( t i − , · ) (cid:17) .Algorithm 1 presents the implementation for computing the gradient of J and it is based on expressions from appendix A). The computation of theJacobian determinant (cid:12)(cid:12)(cid:12) Det(d ϕ ν t i , ( x )) (cid:12)(cid:12)(cid:12) at each time point is based on thefollowing approximation similar to [9]: (cid:12)(cid:12)(cid:12) Det (cid:0) d ϕ ν t i , ( x ) (cid:1)(cid:12)(cid:12)(cid:12) ≈ (cid:16) N div ν ( t i , · ) (cid:17)(cid:12)(cid:12) D ϕ ν t i +1 , (cid:12)(cid:12) ◦ (cid:16) Id + 1 N ν ( t i , · ) (cid:17) . The forward operator
Let X = L (Ω , R ) whose elements represent 2Dimages on a fixed bounded domain Ω ⊂ R . In the application shown here,diffeomorphisms act on X through a geometric group action in (3) and thegoal is to register a given differentiable template image I ∈ X against atarget that observed indirectly as in (4).The forward operator A : X → Y in 2D tomographic imaging is the 2Dray/Radon transform, i.e., A ( f )( ω, x ) = (cid:90) R f ( x + sω )d s for ω ∈ S and x ∈ ω ⊥ .Here, S is the unit circle, so ( ω, x ) encodes the line s (cid:55)→ x + sω in R withdirection ω through x . The data manifold M is the set of such lines that areincluded in the measurements, i.e., M is given by the experimental set-up.We will consider parallel lines in R (parallel beam data), i.e., tomographicdata are noisy digitized values of an L -function on this manifold so g ∈ Y = L ( M, R ). The forward operator is linear, so it is particular Gateauxdifferentiable, and the adjoint of its derivative is given by the backprojection,see [17, 15] for further details.If data is corrupted by additive Gaussian noise, so a suitable data like-lihood is the 2-norm, i.e., L : Y × Y → R with L ( g, h ) = (cid:107) g − h (cid:107) . The noise level in data is specified by the peak signal-to-noise ratio (PSNR),which is defined asPSNR( g ) = 10 log (cid:18) (cid:107) g − g (cid:107) (cid:107) e − e (cid:107) (cid:19) for g = g + e . https://github.com/bgris/odl/tree/IndirectMatching/examples/Metamorphosis lgorithm 1 Computation of ∇ J ( ν , ζ ). Require: ν ( t i , · ) and ζ ( t i , · ) with t i ← i/N for i = 0 , , . . . , N . for i = 1 , . . . , N do (cid:46) Compute ζ ( t i , · ) ◦ ϕ ν ,t i temp ← ζ ( t i , · ) for j = i − , . . . , do temp ← temp ◦ (cid:16) Id + N ν ( t j , · ) (cid:17) end for ζ ( t i , · ) ◦ ϕ ν ,t i ← temp end for for i = 1 , . . . , N do (cid:46) Compute f ν ,ζ ( t i , · ) := I ν ,ζ ( t i , · ) ◦ ϕ ν t i , I ν ,ζ ( t i , · ) ← I + (cid:80) i − j =0 I ν ,ζ ( t j , · ) + N ζ ( t j , · ) ◦ ϕ ν ,t j I ν ,ζ ( t i , · ) ◦ ϕ ν , ← I ν ,ζ ( t i , · ) for j = 1 , . . . , i do I ν ,ζ ( t i , · ) ◦ ϕ ν t j , ← (cid:0) I ν ,ζ ( t i , · ) ◦ ϕ ν ,t j − (cid:1) ◦ (cid:16) Id − N ν ( t j − , · ) (cid:17) end for end for for i = 1 , . . . , N do (cid:46) Compute I ◦ ϕ ν t i , I ◦ ϕ ν , ← I ◦ ϕ ν , = I I ◦ ϕ ν t i , ← (cid:0) I ◦ ϕ ν t i − , (cid:1) ◦ (cid:16) Id − N ν ( t i − , · ) (cid:17) end for for i = 1 , . . . , N do G ( t i , · ) ← ∇ ( I ◦ ϕ ν t i , ) + (cid:80) t i − j =0 N ∇ ( ζ ( t j , · ) ◦ ϕ ν t i ,t j ) end for (cid:12)(cid:12) D ϕ ν t N , (cid:12)(cid:12) = (cid:12)(cid:12) D ϕ ν , (cid:12)(cid:12) = 1 (cid:46) Compute (cid:12)(cid:12) D ϕ ν t i , (cid:12)(cid:12) for i = N − , . . . , do (cid:12)(cid:12) D ϕ ν t i , (cid:12)(cid:12) ← (cid:16) N div ν ( t i , · ) (cid:17)(cid:12)(cid:12) D ϕ ν t i +1 , (cid:12)(cid:12) ◦ (cid:16) Id + N ν ( t i , · ) (cid:17) end for ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t N , (cid:1) ← ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1) for i = N − , . . . , do (cid:46) Compute ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t i , (cid:1) ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t i , (cid:1) ← ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t i +1 , (cid:1) ◦ (cid:16) Id + N ν ( t i , · ) (cid:17) end for for i = 1 , . . . , N do (cid:46) Compute ∇ J ( ν , ζ ) ∇ ν J γ,τ ( ν , ζ, g )( t i , , · ) ← γ ν ( t i , , · ) − (cid:90) Ω K ( x, · ) (cid:12)(cid:12)(cid:12) Det(d ϕ ν t i , ( x )) (cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t i , ( x ) (cid:1) G ( t i , x )d x ∇ ζ J γ,τ ( ν , ζ )( t i , , · ) ← τ ζ ( t i , , · )+ | Det(d ϕ ν t i , ) (cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζ (1 , · ) , g (cid:1)(cid:0) ϕ ν t i , ( x ) (cid:1) G ( t i , · ) end for return ∇ J ( ν )( t i , · ), ∇ J ( ζ )( t i , · ) for i = 1 , . . . , N .
16n the above, g is the noise-free part and e is the noise component of datawith g and e denoting the mean of g and e , respectively. The PSNR isexpressed in terms of dB. Joint tomographic reconstruction and registration
Under the ge-ometric group action (3), metamorphosis based-indirect registration readsas f (cid:98) ν , (cid:98) ζ = W (cid:0) ϕ (cid:98) ν , , I (cid:98) ν , (cid:98) ζ (cid:1) = I (cid:98) ν , (cid:98) ζ ◦ ϕ (cid:98) ν , where ( (cid:98) ν , (cid:98) ζ ) ∈ L ([0 , , V × X ) minimizes (7), i.e., given regularizationparameters γ, τ ≥ I ∈ X we solvemin ( ν ,ζ ) (cid:20) γ (cid:107) ν (cid:107) + τ (cid:107) ζ (cid:107) + (cid:13)(cid:13)(cid:13) A (cid:16) f (cid:0) , φ (1 , · ) − (cid:1)(cid:17) − g (cid:13)(cid:13)(cid:13) (cid:21) dd t f ( t, x ) = ζ (cid:0) t, φ ( t, x ) (cid:1) f (0 , x ) = I ( x ) ddt φ ( t, x ) = ν (cid:0) t, φ ( t, x ) (cid:1) φ (0 , x ) = x. (25)We will consider a set V of vector fields that is an RKHS with a reproducingkernel represented by symmetric and positive definite Gaussian. Then V isadmissible and is continuously embedded in L (Ω , R ). The kernel we pickis K σ : Ω × Ω → R × K σ ( x, y ) := exp (cid:16) − σ (cid:107) x − y (cid:107) (cid:17) (cid:18) (cid:19) for x, y ∈ R and σ >
0. (26)The kernel-size σ also acts as a regularization parameter. In the following we perform a number of experiments that tests variousaspects of using metamorphoses based indirect registration for joint tomo-graphic reconstruction and registration. The tomographic inverse problemalong with characteristics of the data are outlined in section 5.1. The resultsare obtained by solving (25) via a gradient descent, see appendix A for thecomputation of the gradient of the objective. For each reconstruction, welist the the number of angles of the parallel beam ray transform, the kernel-size σ in (26), and the two regularisation parameters γ, τ > ×
256 pixels in the imagedomain Ω = [ − , × [ − , , π ] with 362 lines/angle. Data is corruptedwith additive Gaussian noise with differing noise levels. Here, topology of the template is consistent with that of the target, butintensities differ. The template, which is shown in fig. 2(a), is registeredagainst tomographic data shown in fig. 2(c). The (unknown) target usedto generate data is shown in fig. 2(b). Also, data has a noise level corre-sponding to a PSNR of 15 . σ = 2, which should becompared to the size of the image domain Ω = [ − , × [ − , Here, both topology and intensities of the template differ from those in thetarget. The template, which is shown in fig. 3(a), is registered against tomo-graphic data shown in fig. 3(c). The (unknown) target used for generatingthe data is shown in fig. 3(b). Also, data has a noise level corresponding to aPSNR of 10 . σ = 2, which should be compared to thesize of the image domain Ω = [ − , × [ − , γ − − − − − − − − σ = 3 forseveral regularisation parameters.is shown in fig. 3(h), which is to be compared against the target in fig. 3(b).Figure 3 also shows image, deformation and template trajectories.We clearly see that metamorphosis based indirect registration can handlea template where both intensities and the topology are wrong. In particular,we can see follow both the deformation of the template and the appearanceof the white disc. Metamorphosis based indirect registration, which amounts to solving (25),requires selecting three parameters: the kernel-size σ and the two regularisa-tion parameters γ and τ . Here we study the influence of these parameters onthe final registered image (reconstruction) based on the setup in section 5.4.The reconstruction along with its template and deformation parts arenot that sensitive to the specific choice the two regularisation parameters γ and τ , see table 1 that shows the structural similarity (SSIM) and PSNRvalues for various values of γ and τ when σ = 3. The reconstruction is onthe other hand more sensitive to the choice of the kernel size, see table 2for a table of SSIM and PSNR values corresponding to different choices ofkernel size. Figure 4 also shows reconstructed image with the correspondingfinal template and deformation parts for various values of σ . Interestingly,even if the reconstruction is satisfying for the various values of the kernelsize σ , its template part and deformation parts are really different. Thegeometric deformation and the change in intensity values seem to balancein an non-intuitive way in order to produce a reasonable final image.19 . . σ and fixed regularisation parameters γ = τ = 10 − . The goal here is to recover the unknown temporal evolution of a templatematched against (gated) parallel beam 2D ray transform data acquired at10 different time points (from t = 0 . t = 1), so the target undergoesa temporal evolution. At each of the 10 time points, we only have limitedtomographic data in the sense that i :th acquisition corresponds to samplingthe parallel beam ray transform of the target at time t i using 10 angles ran-domly distributed in [( i − π/ , iπ/
10] using 362 lines/angles. Similarly toprevious experiments, the reconstruction space is Ω = [ − , × [ − , ×
256 pixel grey scale images.The registration of the template I against the temporal series of data g i ,1 ≤ i ≤
10 at the 10 time points t i is performed by minimizing the followingfunctional with respect to one trajectory ( ν , ζ ) ∈ L ([0 , , V × X ): J γ,τ ( ν , ζ ; g , . . . , g ) := γ (cid:107) ν (cid:107) + τ (cid:107) ζ (cid:107) + (cid:88) i =1 L (cid:16) A (cid:0) W ( ϕ ν ,t i , I ν ,ζt i ) (cid:1) , g i (cid:17) where t (cid:55)→ I ν ,ζt , is the absolutely continuous solution to dd t I ν ,ζt ( x ) = ζ (cid:0) t, ϕ ν ,t ( x ) (cid:1) I ν ,ζ ( x ) = I ( x ) with ϕ ν ,t ∈ G V as in (2).The target, the gated tomographic data, and the three trajectories (im-age, deformation and template) resulting from the metamorphosis basedindirect registration are shown in fig. 5. We see that metamorphosis basedindirect registration can be used for spatio-temporal reconstruction evenwhen (gated) data is highly under sampled. In particular, we can recoverthe evolution (both the geometric deformation and the appearance of thewhite disc) of the target. As a comparison, fig. 5(f) presents reconstructionsobtained from filtered back projection (FBP) and total variation (TV). Here,data is a concatenation of the 10 gated data sets, thereby corresponding thensampling the ray transform using 100 angles in [0 , π ]. Note however that thetemporal evolution of the target is not accounted for in these reconstruc-tions. 20 Conclusions and discussion
We introduced a metamorphosis-based framework for indirect registrationand showed that this corresponds to a well-defined regularization method.We also present several numerical examples from tomography.In particular, section 5.6 illustrates that this framework enables to re-cover the temporal evolution of a template from temporal data, even whendata are very limited for each time point. This approach assumes that onehas access to an initial template. In spatio-temporal reconstruction, suchan initial template is unknown and it needs to be recovered as well. Oneapproach for doing this is by an intertwined scheme that alternates betweento steps (similarly to [14]): (i) given a template, estimate its evolution thatis consistent with times series of data using the metamorphosis frameworkfor indirect registration, and (ii) estimate the initial template from timesseries of data given its evolution. The approach in section 5.6 solves thefirst step, which is the more difficult one.Another topic is the choice of hyperparameters. Our metamorphosis-based framework for indirect registration relies on three parameters, butas shown in section 5.5, the most important one is the kernel-size σ . Thelatter has a strong influence on the way the reconstructed image trajectorydecomposes into a deformation and a template part. Clearly it acts asa regularisation parameter and a natural problem is to devise a scheme forchoosing it depending on the size of features (scale) undergoing deformation.Unfortunately, similarly to direct registration using the LDDMM framework,the choice of this parameter (and more generally choice of kernel for theRKHS V ) is still an open problem [3, 9, 10]. One way is to use a multi-scale approach [7, 20, 22] but a general method for selecting an appropriatekernel-size remains to be determined. The work by Ozan ¨Oktem, Barbara Gris and Chong Chen was supportedby the Swedish Foundation for Strategic Research grant AM13-0049.21 a) Template. (b) Target. (c) Data (sinogram). I m ag e (d) t = 0. (e) t = 0 .
2. (f) t = 0 .
5. (g) t = 0 .
7. (h) t = 1. D e f o r m a t i o n (i) t = 0. (j) t = 0 .
2. (k) t = 0 .
5. (l) t = 0 .
7. (m) t = 1. I n t e n s i t y (n) t = 0. (o) t = 0 .
2. (p) t = 0 .
5. (q) t = 0 .
7. (r) t = 1. Figure 2: Metamorphosis based indirect-matching of template in (a) againstdata in (c), which represents 2D ray transform of target in (b) (100 uniformlydistributed angles in [0 , π ]). The second row (d)–(h) shows the image tra-jectory t (cid:55)→ W ( ϕ ν ,t , f t ( ν , ζ )), so the final registered template is in (h). Thethird row (i)–(m) shows the deformation trajectory t (cid:55)→ W ( ϕ ν ,t , I ), likewisethe fourth row (n)–(r) shows the intensity trajectory t (cid:55)→ f t ( ν , ζ ).22 a) Template I . (b) Target. (c) Data (sinogram). I m ag e (d) t = 0. (e) t = 0 .
2. (f) t = 0 .
5. (g) t = 0 .
7. (h) t = 1. D e f o r m a t i o n (i) t = 0. (j) t = 0 .
2. (k) t = 0 .
5. (l) t = 0 .
7. (m) t = 1. I n t e n s i t y (n) t = 0. (o) t = 0 .
2. (p) t = 0 .
5. (q) t = 0 .
7. (r) t = 1. Figure 3: Metamorphosis based indirect-matching of template in (a) againstdata in (c), which represents 2D ray transform of target in (b) (100 uniformlydistributed angles in [0 , π ]). The second row (d)–(h) shows the image tra-jectory t (cid:55)→ W ( ϕ ν ,t , f t ( ν , ζ )), so the final registered template is in (h). Thethird row (i)–(m) shows the deformation trajectory t (cid:55)→ W ( ϕ ν ,t , I ), likewisethe fourth row (n)–(r) shows the intensity trajectory t (cid:55)→ f t ( ν , ζ ).23mage Deformation Template σ = . σ = . σ = σ = σ = σ = σ = Figure 4: Reconstruction results and their recomposition in template partand deformation part for various kernel size σ .24 a) The temporal evolution of the target.(b) The (gated) tomographic data. Each data set is highly incomplete (limited angle).(c) Image trajectory (reconstruction), combines deformation and template trajectories.(d) Deformation trajectory, models mainly geometric changes.(e) Template trajectory, models mainly intensity changes.(f) FBP (left) and TV (middle) reconstructions from concatenated data(right). Figure 5: Reconstructing the temporal evolution of a template using meta-morphosis. Target (a), data (b), and results (c)–(e), are shown at selectedtime points t = 0 . , . , . , .
8, and 1 .
0. As a comparison we show re-constructions assuming static target obtained from concatenating the gatedtomographic data (f). 25
Gradient computation
This section presents the computation of the gradient of J γ,τ ( · ; g ), whichis useful for any first order optimisation metod for minimising the functional J γ,τ ( · ; g ) in (7). The computations assume I ∈ X ∩ C (Ω , R ) and ( ν , ζ ) ∈ L ([0 , , V × X ) with X = L (Ω , R ).Furthermore, for each t ∈ [0 ,
1] we also assume t (cid:55)→ ζ ( t, · ) ∈ C (Ω , R ). Innumerical implementations, we consider digitized images and considerationsof the above type are not that restrictive.Let us first compute the differential of the data discrepancy term withrespect to ζ using the notation f ν ,ζt := W ( ϕ ν t, , I ν ,ζt ) = I ν ,ζt ◦ ϕ ν t, . As notedin (13), we have f ν ,ζt ( x ) = ( I ν ,ζt ◦ ϕ ν t, )( x ) = I (cid:0) ϕ ν t, ( x ) (cid:1) + (cid:90) t ζ (cid:0) τ, ϕ ν t,τ ( x ) (cid:1) d τ. (27)Then ∂ ζ (cid:104) L (cid:0) f ν ,ζt , g (cid:1)(cid:105) ( ζ )( η ) = (cid:68) ∇ L ( f ν ,ζt , g ) , ∂ ζ f ν ,ζt ( ζ )( η ) (cid:69) = (cid:90) Ω (cid:90) t ∇ L (cid:0) f ν ,ζ , g (cid:1) η ( τ, ϕ ν t,τ ( x ))d τ d x = (cid:90) Ω (cid:90) τ ≤ t | Det(d ϕ ν τ,t ( x )) |∇ L ( f ν ,ζt , g )( ϕ ν τ,t ( x )) η ( τ, x )d τ d x = (cid:68) · ≤ t | Det(d ϕ ν · ,t ) |∇ L ( f ν ,ζt , g ))( ϕ ν · ,t ) , η (cid:69) L ([0 , ,L (Ω , R )) . In order to compute the differential of the discrepancy term with respectto ν , we start by computing the differential of f ν ,ζ with respect to ν . Hence,let µ ∈ L ([0 , , V ) and x ∈ Ω. Thendd (cid:15) f ν + (cid:15) µ ,ζt ( x ) (cid:12)(cid:12)(cid:12) (cid:15) =0 = (cid:68) ∇ I (cid:0) ϕ ν t, ( x ) (cid:1) , dd (cid:15) ϕ ν + (cid:15) µ t, ( x ) (cid:12)(cid:12) (cid:15) =0 (cid:69) + (cid:90) t (cid:68) ∇ ζ ( τ, ϕ ν t,τ ( x )) , dd (cid:15) ϕ ν + (cid:15) µ t,τ ( x ) (cid:12)(cid:12) (cid:15) =0 (cid:69) d τ = − (cid:90) t (cid:68) ∇ I ( ϕ ν t, ( x )) , d ϕ ν s, (cid:0) ϕ ν t,s ( x ) (cid:1)(cid:16) µ ( s, ϕ ν t,s ( x ) (cid:1)(cid:17)(cid:69) d s (cid:90) t (cid:28) ∇ ζ (cid:0) τ, ϕ ν t,τ ( x ) (cid:1) , (cid:90) τt d ϕ ν s,τ (cid:0) ϕ ν t,s ( x ) (cid:1)(cid:16) µ (cid:0) s, ϕ ν t,s ( x ) (cid:1)(cid:17) d s (cid:29) d τ = − (cid:90) t (cid:28) ∇ I (cid:0) ϕ ν t, ( x ) (cid:1) , d ϕ ν s, (cid:0) ϕ ν t,s ( x ) (cid:1)(cid:16) ( µ (cid:0) ( s, ϕ ν t,s ( x ) (cid:1)(cid:17)(cid:29) d s − (cid:90) t (cid:90) s (cid:28) ∇ ζ ( τ, · ) ◦ ϕ ν t,τ ( x ) , d ϕ ν s,τ (cid:0) ϕ ν t,s ( x ) (cid:1)(cid:16) µ (cid:0) s, ϕ ν t,s ( x ) (cid:1)(cid:17)(cid:29) d τ d s. (28)26sing (28), we can compute the derivative of (cid:15) (cid:55)→ L (cid:0) W ( ϕ ν + (cid:15) µ t , I ν + (cid:15) µ ,ζt ) (cid:1) at (cid:15) = 0:dd (cid:15) L (cid:0) W (cid:0) ϕ ν + (cid:15) µ t , I ν + (cid:15) µ ,ζt (cid:1)(cid:1)(cid:12)(cid:12) (cid:15) =0 = (cid:68) ∇ L (cid:0) f ν ,ζt , g (cid:1) , dd (cid:15) f ν + (cid:15) µ ,ζt (cid:12)(cid:12) (cid:15) =0 (cid:69) = − (cid:90) Ω (cid:40) (cid:90) t ∇ L (cid:0) f ν ,ζt , g (cid:1) ( x ) · (cid:34)(cid:68) ∇ I (cid:0) ϕ ν t, ( x ) (cid:1) , d ϕ ν s, (cid:0) ϕ ν t,s ( x ) (cid:1) ( µ ( s, ϕ ν t,s ( x ))) (cid:69) + (cid:90) s (cid:28) ∇ ζ ( τ, · ) ◦ ϕ ν t,τ ( x ) , d ϕ ν s,τ (cid:0) ϕ ν t,s ( x ) (cid:1)(cid:16) µ (cid:0) s, ϕ ν t,s ( x ) (cid:1)(cid:17)(cid:29) d τ (cid:35) d s (cid:41) d x = − (cid:90) Ω (cid:40) (cid:90) t (cid:12)(cid:12)(cid:12) Det (cid:0) d ϕ ν s,t ( x ) (cid:1)(cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζt , g (cid:1)(cid:0) ϕ ν s,t ( x ) (cid:1) · (cid:20)(cid:68) ∇ I (cid:0) ϕ ν s, ( x ) (cid:1) , d ϕ ν s, ( x ) (cid:0) µ ( s, x ) (cid:1)(cid:69) + (cid:90) s (cid:68) ∇ ζ ( τ, · ) ◦ ϕ ν s,τ ( x ) , d ϕ ν s,τ ( x ) (cid:0) µ ( s, x ) (cid:1)(cid:69) d τ (cid:21) d s (cid:41) d x = − (cid:90) Ω (cid:90) t (cid:12)(cid:12)(cid:12) Det (cid:0) d ϕ ν s,t ( x ) (cid:1)(cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζt , g (cid:1)(cid:0) ϕ ν s,t ( x ) (cid:1) · (cid:20)(cid:10) ∇ ( I ◦ ϕ ν s, )( x ) , µ ( s, x ) (cid:11) + (cid:90) s (cid:68) ∇ ( ζ ( τ, · ) ◦ ϕ ν s,τ )( x ) , µ ( s, x ) (cid:69) d τ (cid:35) d s d x = − (cid:90) Ω (cid:90) (cid:28) s ≤ t (cid:12)(cid:12)(cid:12) Det (cid:0) d ϕ ν s,t ( x ) (cid:1)(cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζt , g (cid:1)(cid:0) ϕ ν s,t ( x ) (cid:1)(cid:20) ∇ ( I ◦ ϕ ν s, )( x ) + (cid:90) s ∇ ( ζ ( τ, · ) ◦ ϕ ν s,τ )( x )d τ (cid:21) , µ ( s, x ) (cid:29) d s d x = − (cid:28) · ≤ t (cid:12)(cid:12)(cid:12) Det(d ϕ ν · ,t ) (cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζt , g (cid:1) ◦ ϕ ν · ,t (cid:20) ∇ ( I ◦ ϕ ν · , )( · )+ (cid:90) · ∇ ( ζ ( τ, · ) ◦ ϕ ν · ,τ )( · )d τ (cid:21) , µ (cid:29) L ([0 , ,L (Ω , R d )) = − (cid:28)(cid:90) Ω K ( x, · )1 · ≤ t (cid:12)(cid:12)(cid:12) Det(d ϕ ν · ,t ( x )) (cid:12)(cid:12)(cid:12) ∇ L (cid:0) f ν ,ζt , g (cid:1)(cid:0) ϕ ν · ,t ( x ) (cid:1)(cid:20) ∇ ( I ◦ ϕ ν · , )( x )+ (cid:90) · ∇ ( ζ ( τ, · ) ◦ ϕ ν · ,τ )( x )d τ (cid:21) , µ (cid:29) L ([0 , ,V ) . References [1] S. Arguillere, E. Tr´elat, A. Trouv´e, and L. Younes. Shape deformationanalysis from the optimal control viewpoint.
Journal de Math´ematiquesPures et Appliqu´es , 104(1):139–178, 2015.272] M. Bauer, M. Bruveris, and P. W. Michor. Overview of the geometriesof shape spaces and diffeomorphism groups.
Journal of MathematicalImaging and Vision , 50(1–2):60–97, 2014.[3] M. F. Beg, M. I. Miller, A. Trouv´e, and L. Younes. Computing largedeformation metric mappings via geodesic flows of diffeomorphisms.
International journal of computer vision , 61(2):139–157, 2005.[4] M. Bertero, H. Lant´eri, and L. Zanni. Iterative image reconstruction:a point of view. In Y. Censor, M. Jiang, and A. K. Louis, editors,
Pro-ceedings of the Interdisciplinary Workshop on Mathematical Methods inBiomedical Imaging and Intensity-Modulated Radiation (IMRT), Pisa,Italy , pages 37–63, 2008.[5] E. Bladt, D. M. Pelt, S. Bals, and K. J. Batenburg. Electron tomogra-phy based on highly limited data using a neural network reconstructiontechnique.
Ultramicroscopy , 158:81–88, 2015.[6] M. Bruveris and D. D. Holm. Geometry of image registration: Thediffeomorphism group and momentum maps. In D. E. Chang, D. D.Holm, G. Patrick, and T. Ratiu, editors,
Geometry, Mechanics, andDynamics , pages 19–56. Springer, 2015.[7] M. Bruveris, L. Risser, and F.-X. Vialard. Mixture of kernels and iter-ated semidirect product of diffeomorphisms groups.
Multiscale Modeling& Simulation , 10(4):1344–1368, 2012.[8] N. Charon, B. Charlier, and A. Trouv´e. Metamorphoses of functionalshapes in Sobolev spaces.
Foundations of Computational Mathematics ,pages 1–62, 2016.[9] C. Chen and O. ¨Oktem. Indirect image registration with large diffeo-morphic deformations.
SIAM Journal of Imaging Sciences , 11(1):575–617, 2018.[10] S. Durrleman, M. Prastawa, N. Charon, J. R. Korenberg, S. Joshi,G. Gerig, and A. Trouv´e. Morphometry of anatomical shape complexeswith dense deformations and sparse parameters.
NeuroImage , 101:35–49, 2014.[11] A. Effland, M. Rumpf, and F. Sch¨afer. Image extrapolation for thetime discrete metamorphosis model: Existence and applications.
SIAMJournal of Imaging Sciences , 11(1):834–862, 2018.[12] M. Grasmair. Generalized Bregman distances and convergence rates fornon-convex regularization methods.
Inverse Problems , 26(11):115014,2010. 2813] G. T. Gullberg, B. W. Reutter, A. Sitek, J. S. Maltz, and T. F.Budinger. Dynamic single photon emission computed tomography –basic principles and cardiac applications.
Physics in Medicine and Bi-ology , 55:R111–R191, 2010.[14] J. Hinkle, M. Szegedi, B. Wang, B. Salter, and S. Joshi. 4D CT im-age reconstruction with diffeomorphic motion model.
Medical imageanalysis , 16(6):1307–1316, 2012.[15] A. Markoe.
Analytic Tomography , volume 106 of
Encyclopedia of math-ematics and its applications . Cambridge University Press, 2006.[16] M. I. Miller, L. Younes, and A. Trouv´e. Diffeomorphometry andgeodesic positioning systems for human anatomy.
Technology , 2(1),2014.[17] F. Natterer and F. W¨ubbeling.
Mathematical Methods in Image Re-construction . Mathematical Modeling and Computation. Society forIndustrial and Applied Mathematics, 2001.[18] S. Neumayer, J. Persch, and G. Steidl. Regularization of inverseproblems via time discrete geodesics in image spaces. arXiv preprintarXiv:1805.06362 , 2018.[19] O. ¨Oktem, C. Chen, N. O. Domani¸c, P. Ravikumar, and C. Bajaj.Shape based image reconstruction using linearised deformations.
In-verse Problems , 33(3):035004, 2017.[20] L. Risser, F.-X. Vialard, R. Wolz, M. Murgasova, D. D. Holm, andD. Rueckert. Simultaneous multi-scale registration using large defor-mation diffeomorphic metric mapping.
IEEE transactions on medicalimaging , 30(10):1746–1759, 2011.[21] A. Schwarz and M. Leach. Implications of respiratory motion for thequantification of 2D MR spectroscopic imaging data in the abdomen.
Physics in Medicine and Biology , 45(8):2105—2116, 2000.[22] S. Sommer, M. Nielsen, F. Lauze, and X. Pennec. A multi-scale ker-nel bundle for LDDMM: towards sparse deformation description acrossspace and scales. In
Biennial International Conference on InformationProcessing in Medical Imaging , pages 624–635. Springer Verlag, 2011.[23] A. Sotiras, C. Davatzikos, and N. Paragios. Deformable medical im-age registration: A survey.
IEEE Transactions on Medical Imaging ,32(7):1153–1190, 2013.[24] A. Trouv´e and L. Younes. Metamorphoses through Lie group action.
Foundations of Computational Mathematics , 5(2):173–198, 2005.2925] Y. Wang, E. Vidan, and G. Bergman. Cardiac motion of coronaryarteries: Variability in the rest period and implications for coronaryMR angiography.
Radiology , 213(3):751—758, 1999.[26] L. Younes.
Shapes and Diffeomorphisms , volume 171 of