1 Introduction
In mathematical imaging, image morphing is the problem of computing a visually appealing transition of two images such that semantically corresponding regions are mapped onto each other. A wellknown approach for image morphing is the metamorphosis model originally introduced by Miller, Trouvé, and Younes [MY01, TY05b, TY05a], which generalizes the flow of diffeomorphism model and the large deformation diffeomorphic metric mapping (LDDMM) which dates back to the pioneering work of Arnold [Arn66] with its exploration and extension in imaging by Dupuis, Grenander and others [DGM98, BMTY05, JM00, MTY02, VS09, VRRC12]. From the perspective of the flow of diffeomorphism model, each point of the reference image is transported to the target image in an energetically optimal way such that the image intensity is preserved along the trajectories of the pixels. Here, the energy measures the total dissipation of the underlying flow. The metamorphosis model additionally allows for image intensity modulations along the trajectories by incorporating the magnitude of these modulations, which is reflected by the integrated squared material derivative of image trajectories as a penalization term in the energy functional. Recently, the metamorphosis model has been extended to images on Hadamard manifolds [NPS18, ENR19], to reproducing kernel Hilbert spaces [RY16], to functional shapes [CCT18] and to discrete measures [RY13]. For a more detailed exposition of these models we refer the reader to [You10, MTY15] and the references therein.
Starting from the general framework for variational time discretization in geodesic calculus [RW15], a variational time discretization of the metamorphosis model for squareintegrable images was proposed in [BER15]. Moreover, the existence of discrete geodesic paths as well as the Mosco–convergence of the time discrete to the time continuous metamorphosis model was proven. However, the classical metamorphosis model, its time discrete counterpart and the spatial discretization based on finite elements in [BER15] exhibit several drawbacks:

The comparison of images in their original gray or color space is not invariant to natural radiometric transformations caused by lighting or material changes, shadows etc. and hence might lead to a blending along the discrete geodesic path instead of flowinduced geometric transformations.

Texture patterns, which are important for a natural appearance of images, are often destroyed along the geodesic path due to the colorbased matching.

Sharp interfaces such as object boundaries, which frequently coincide with depth discontinuities of a scene, are in general not preserved along a geodesic path because of the strong smoothness implied by the homogeneous and isotropic variational prior for the deformation fields.
To overcome these problems originating from the intensitybased matching, we propose a multiscale feature space approach incorporating the deep convolution neural network introduced in
[SZ14]. In detail, this convolutional neural network, which was trained to classify the ImageNet dataset
[KSH12], extracts semantic features using weight layers, each composed of smallconvolution filters with subsequent nonlinear ReLU activation functions. This network defines a feature extraction operator, where each feature map is considered as a continuous map into some higherdimensional feature space consisting of vectors in
, where ranges from to depending on the considered scale associated with a certain network layer. Compared to the original time discrete metamorphosis model [BER15] we advocate a metamorphosis model in a deep feature space, which amounts to replacing the input images by feature vectors combining image intensities and semantic information generated by the feature extraction operator. To explicitly allow for discontinuities in the deformation fields, we introduce an anisotropic regularization of the time discrete deformation sequence. Since motion discontinuities and object interfaces in images commonly coincide, the considered anisotropy solely depends on the magnitude of image gradients.We prove the existence of discrete geodesic paths for the deep feature metamorphosis model and discuss its Mosco–convergence to the appropriate time continuous metamorphosis model in deep feature space. This in particular implies the convergence of time discrete to time continuous geodesic paths and establishes the existence of time continuous geodesics as minimizers of the time continuous metamorphosis model.
We propose a finite difference/third order Bspline discretization for the fully discrete feature space metamorphosis model and use the iPALM algorithm [PS16] for the optimization, which leads to an efficient and robust computation of morphing sequences that visually outperform the prior intensitybased finite element discretization discussed in [BER15]. This scheme is significantly less sensitive to intensity modulations due to the exploitation of semantic information.
Note that this publication is an extended version of the conference proceeding [EKPR19], in which the model is adapted and in addition a rigorous mathematical analysis of this novel model is presented. In fact, the morphing sequence is no longer retrieved in a postprocessing step. Instead, the color values are part of the feature vector. Different from the prior proceedings article, we prove the existence of time discrete geodesics in feature space, present a time continuous model and discuss the issue of convergence of the discrete functionals.
Notation.
Throughout this paper, we assume that the image domain for is bounded and strongly Lipschitz. We use standard notation for Lebesgue and Sobolev spaces from the image domain to a Banach space , i.e. and and omit if the space is clear from the context. The associated norms are denoted by and , respectively, and the seminorm in is given by , i.e.
for . We use the notation for Hölder spaces of order with regularity , the corresponding (semi)norm is
The symmetric part of a matrix is denoted by , i.e. and the symmetrized Jacobian of a differentiable function by . We denote by the elements of with positive determinant, and by
both the identity map and the identity matrix. Finally,
refers to the temporal derivative of a differentiable function .Organization.
This paper is structured as follows: in section 2, we review the classical metamorphosis model and present its extension to deep feature spaces. Then, in section 3 we introduce the time discrete deep feature metamorphosis model and prove the existence of geodesic paths. In section 4, we present a time continuous metamorphosis model and comment on the Mosco–convergence in deep feature space. The fully discrete model and the optimization scheme using the iPALM algorithm are presented in section 5. Finally, in section 6 several examples demonstrate the applicability of the proposed methods to real image data.
2 Metamorphosis model
In this section, we briefly review the classical flow of diffeomorphism model and the metamorphosis model as its generalization. Then, we extend the metamorphosis model to the space of deep features, where we additionally incorporate an anisotropic regularization.
2.1 Flow of diffeomorphism
In what follows, we present a very short exposition of the flow of diffeomorphism model and refer the reader to [DGM98, BMTY05, JM00, MTY02] for further details. In the flow of diffeomorphism model, the temporal change of image intensities is determined by a family of diffeomorphisms describing a flow transporting image intensities along particle paths. The main assumption of this model is the brightness constancy assumption, which is equivalent to a vanishing material derivative along a path in the space of images, where denotes the timedependent Eulerian velocity. The Riemannian space of images is endowed with the following metric and path energy
Note that we use as a shortcut for the function . Here, the quadratic form is the higher order elliptic operator
where and . Physically, the metric describes the viscous dissipation in a multipolar fluid model as investigated by Nečas and Šilhavý [Nv91]. The first two terms of the integrand represent the dissipation density in a Newtonian fluid and the third term can be regarded as a higher order measure for friction. Following [DGM98, Theorem 2.5], paths with a finite energy, which connect two diffeomorphisms and , are actually oneparameter families of diffeomorphisms. Given two image intensity functions , an associated geodesic path is a family of images with and , which minimizes the path energy. The resulting flow of images is given by .
2.2 Metamorphosis model in image space
The metamorphosis approach originally proposed by Miller, Trouvé, Younes and coworkers in [MY01, TY05b, TY05a] generalizes the flow of diffeomorphism model by allowing for image intensity variations along motion paths and penalizing the squared material derivative in the metric. Under the assumption that the image path is sufficiently smooth, the metric and the path energy read as
for a penalization parameter , where denotes the material derivative of . The Lagrangian formulation of this variation of the image intensity along motion trajectories can be phrased as follows: for all we have
(1) 
Hence, the flow of diffeomorphism model is the limit case of the metamorphosis model for . This definition of the metric has two major drawbacks: In general, paths in the space of images do not exhibit any smoothness properties (neither in space nor time), and therefore the evaluation of the material derivative is not welldefined. Moreover, since different pairs of velocity fields and material derivatives can imply the same time derivative of the image path , the restriction to equivalence classes of pairs is required, where two pairs are equivalent if and only if they induce the same temporal change of the image path .
To tackle both problems, Trouvé and Younes [TY05a] proposed a nonlinear geometric structure in the space of RGB images . In detail, for a given image path and an associated velocity field , where denotes the velocity space, the weak material derivative is incorporated in the model, which is implicitly given by
(2) 
for a smooth test function . We consider as a tangent vector in the tangent space of at the image and write defined by (2). Indeed, represents a variation of the image via transport and change of intensity. This (weak) formulation and the consideration of equivalence classes of motion fields and material derivatives inducing the same temporal change of the image intensity gives rise to the notion for regular paths in the space of images. For details we refer the reader to [TY05a]. The path energy in the metamorphosis model for a regular path is then defined as
(3) 
Then, image morphing of two input images amounts to computing a shortest geodesic path in the metamorphosis model, which is defined as a minimizer of the path energy in the class of regular curves such that and . The existence of a shortest geodesic is proven in [TY05a, Theorem 6]. Note that the infimum in (3) is attained, which is shown in [TY05a, Proposition 1 & Theorem 2].
2.3 Metamorphosis model in deep feature space
In this subsection, we extend the metamorphosis model to images as maps into a deep feature space with the aim to increase the reliability and robustness of the resulting morphing. To further improve the quality of the deformations, we incorporate an anisotropic regularization of the deformation field. We will compute geodesics path in the feature space ) for . Here, the first part of a feature vector encodes the RGB image intensity values, the remaining component represents deep features, which are highdimensional local image patterns describing the local structure of the image as a superposition on different levels of a multiscale image approximation. Let us denote by the projection onto the image component of a feature, i.e. . To compute the geodesic sequence in the deep feature space, we extract the features from the fixed input images and define for a fixed (small)
The computation of the VGG features is composed of convolution operators and nonlinear ReLU activation functions which are both continuous mappings. Hence, it is reasonable to assume in our mathematical model that the mapping is continuous. Following [SZ14], we define for the fully discrete model discussed in section 5 a discrete feature operator to incorporate semantic information in image morphing based on convolutional neural networks, where ranges from 64 to 512. The parameter
is used to scale down the RGB component mainly needed to compute the anisotropy (see below) and to primarily focus on the actual VGG features when estimating the transport.
Next, we include an anisotropic elliptic operator in our model to properly account for image structures such as sharp edges or corners. To this end, we consider an anisotropy operator fulfilling the following assumptions:

[leftmargin=5.6ex,itemsep=.2parsep=.1label=(A0)]

boundedness and coercivity: for and all and a.e. ,

compactness: in implies in ,

Lipschitz continuity: for all neighborhoods there exists such that for all .
In the numerical experiments, we use the operator [PM90]
(4) 
for fixed , where
is the Gaussian kernel with standard deviation
. Note that (4) satisfies 1–3.Now we are in the position to introduce the variational model for deep feature metamorphosis. Instead of generalizing the definition of regular paths and adapting the notion of a weak material derivative (2) originally proposed by Trouvé and Younes, we follow the relaxed material derivative approach proposed in [ENR19], in which the material derivative quantity is retrieved from a variational inequality. In [ENR19, Section 3], the equivalence of this energy functional and (3) in the isotropic case has been shown. Let as above denote the Lagrangian flow map induced by the Eulerian motion field with and . Then, we replace the equality (1) (rephrased for the feature map as with being the weak material derivative) by the inequality
(5) 
for a.e. and all , where formally the scalar valued replaces the actually vectorvalued material derivative. In fact, this inequality defines a set of admissible pairs given a path in . This relaxed approach will turn out to be very natural when it comes to lower semicontinuity of the path energy in the context of the existence proof for geodesic paths. For more details we refer the reader to section 4.
Definition 1 (Continuous path energy).
We consider the anisotropic elliptic operator
for an anisotropy weight , a velocity field and . Then, we define the path energy
(6) 
for a path , where
denotes the set of admissible pairs of the velocity and a scalar quantity fulfilling (5).
Let us stress that the anisotropy solely takes into account local RGB values and not the actual VGG features with their discriminative multiscale characteristics.
Geodesic curves in the deep feature space joining are defined as minimizers of the path energy among all curves with the fixed boundary conditions and .
Remark 1.
One observes that a path in feature space with finite energy exhibits additional smoothness properties. Indeed, the boundedness of in implies that the flow is in and, by using Sobolev embedding arguments, in with . The same observation holds for by noting that is the flow associated with the backward motion field . This together with the variational inequality (5) and in ensures that . Using approximation by smooth functions one shows that is uniformly continuous, and by using 3 the mapping is welldefined and in .
3 Variational time discretization
In this section, we develop a variational time discretization of the deep feature space metamorphosis model taking into account the approach presented in [RW15, BER15].
We define the time discrete pairwise energy for two feature maps by
where is given by
(7) 
Here, the set of admissible deformations is
Note that the anisotropy weight only depends on the image component of the second feature in the pairwise energy. We make the following assumptions with respect to the energy density function :

[leftmargin=6.5ex,itemsep=.2parsep=.1label=(W0)]

and is polyconvex and , ,

there exist constants such that for all the growth estimates
hold true,

for all the relation
holds true.
The first two assumptions ensure existence of a minimizing deformation in and the third is a consistency assumption with respect to the differential operator required to guarantee that the below defined discrete path energy is consistent with the time continuous path energy (6).
The particular energy density function
(8) 
used for all numerical experiments satisfies 1–3. The first term enforces the positivity of the determinant of the Jacobian matrix of a deformation and favors a balance of shrinkage and growth as advocated in [DR04, BMR13], while the second term penalizes large deviations of the deformation from the identity. Here, the positivity constraint of the determinant of the Jacobian of the deformations prohibits interpenetration of matter [Bal81].
We proceed with the definition of the discrete path energy and the discrete geodesic between two features .
Definition 2 (Discrete path energy).
Let and . The discrete path energy for a discrete path is defined as
(9) 
A discrete geodesic path morphing into is a discrete tuple that minimizes over all discrete paths with .
For arbitrary vectors and we set
(10) 
In what follows, we will investigate the existence of discrete geodesic curves in the time discrete deep feature space metamorphosis model. To this end, we combine the proofs of the local wellposedness of the pairwise energy with the existence result of a feature vector minimizing for a fixed vector of deformations. We remark that the structure of all proofs is similar to the corresponding proofs in [BER15, Eff18] and we focus on the adaptations necessitated by the anisotropic regularization.
The following lemma, which provides an estimate for the norm of the displacement, is crucial for the wellposedness of the energy.
Lemma 1.
Proof.
Set . An application of the Gagliardo–Nirenberg inequality [Nir66] yields
(11) 
The last term in (11) is bounded by
(12) 
By using the embedding of into and the uniform boundedness of the minimizing sequence in we get To control the lower order term appearing on the righthand side of (11), we define and use 1 and 2 to obtain
which implies . Hence, by the embedding we infer
(13) 
We remark that the inequality
(14) 
holds true, which follows from Korn’s inequality and the Poincaré inequality. Thus, the lemma follows by combining (11), (12), (13) and (14). ∎
Proposition 1 (Wellposedness of ).
Proof.
For fixed , let be a feature vector satisfying (15) for a constant specified below. Let be a minimizing sequence that converges to . Since we can deduce using 1 that
for all . Using again the Gagliardo–Nirenberg inequality we infer that is uniformly bounded in because of the estimate . Due to the reflexivity of there exists a weakly convergent subsequence (not relabeled) such that in . By using the Sobolev embedding theorem as well as the Arzelà–Ascoli theorem we can additionally infer that for a subsequence (again not relabeled) in for holds true. Then, Lemma 1 implies
Thus, by choosing sufficiently small and taking into account the Lipschitz continuity of the determinant we obtain for a constant and all , which implies for a constant . Note that all estimates remain valid for the limit deformation . By [Cia88, Theorem 5.52] the deformations and are diffeomorphisms. Finally, 1 and the lower semicontinuity of the seminorm imply
It remains to verify that
(16) 
as . To this end, we approximate by smooth functions with . Then, using the transformation formula we obtain
where and are pointwise estimated by Finally, by first choosing and then we obtain (16) and thereby verify the claim. ∎
This proposition guarantees the existence of an admissible vector of deformations such that provided that each pair of features contained in satisfies (15).
In what follows, we prove the existence of an energy minimizing vector of features for a fixed vector of deformations.
Proposition 2.
Proof.
We consider a minimizing sequence of features , , for the energy . Then,
A straightforward computation reveals
where we used 1, (17) and the transformation formula. Furthermore, again by (17) one obtains
(18) 
Thus, an induction argument (starting from ) shows that is uniformly bounded in independently of , which implies for a subsequence (not relabeled) in .
In what follows, we prove the weak lower semicontinuity of the discrete path energy along the minimizing sequence. We observe that 2 implies in , which yields
for every . It remains to verify the weak lower semicontinuity of the matching functional, i.e.
(19) 
for every . To this end, we first show in . For every the transformation formula yields
which converges to since due to (17). Hence, in , which readily implies (19). Hence,
which proves the proposition. ∎
We can now combine both previous propositions to prove the existence of discrete geodesics for the deep feature space metamorphosis model.
Comments
There are no comments yet.