Yu M Shukrinov

A Mazanik

I R Rahmonov

A E Botha

A Buzdin

Re-orientation of the easy axis in ϕ 0 -junction

Keywords: 74, 50, +r -Tunneling phenomena; Josephson effects 85, 25, -j -Superconducting devices 85, 25, Cp -Josephson devices

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A ferromagnet can exert a strong influence on a Josephson junction (JJ), particularly, the superconductorferromagnet-superconductor (SFS) JJ oscillates between 0-and π-junctions as the thickness of the ferromagnet increases [1][2][3][4][5][6]. In a closed superconducting loop with π-junction a spontaneous circulating current appears without applied magnetic flux [7]. A very special situation is possible when the weak link is a non-centrosymmetric magnetic metal with broken inversion symmetry like MnSi or FeGe. In this case, Rashba-type spin-orbit coupling leads to a general current-phase relation j = j c sin(ϕ-ϕ 0 ), where ϕ 0 is proportional to the strength of the broken inversion symmetry magnetic interaction [8][9][10]. An important issue is that the ϕ 0 -junction provides a mechanism for direct coupling between the suppercurrent (superconducting phase) and the magnetic moment.

The mechanical analog of a simple pendulum with applied torque is very useful to gain insight into the physics of Josephson devices [START_REF] Barone | Physics and Applications of the Josephson Effect[END_REF]. Here we introduce another analog between the magnetic ϕ 0 -junction and a pendulum with an oscillating point of suspension. A particle moving simultaneously in a permanent field and in a field oscillating with a high frequency, demonstrates a very interesting feature: new stability points appear at some parameters of the particle and fields [START_REF] Kapitza | [END_REF]13]. The system in which this feature was first studied is called the Kapitza pendulum. Particularly, in a pendulum with a vibrating point of suspension (pivot point), the external sinusoidal force can invert the stability position of the pendulum [START_REF] Kapitza | [END_REF]14,15]. In this case the bottom equilibrium position is no longer stable. Any tiny deviation from the vertical increases in amplitude with time. In the presence of a periodic drive, the unstable fixed point can become dynamically stable. Kapitza provided an analytical insight into the reasons for stability, by splitting the motion into "fast" and "slow" variables, and by introducing an effective potential. His result was first obtained in the high-frequency limit. By averaging the classical equations of motion over the fast oscillations of the drive, Kapitza found that the "upper" equilibrium point becomes stable for large enough driving amplitudes, while the "lower" point becomes unstable. Kapitza's pioneering work initiated the field of vibrational mechanics, and his method is used to describe periodic processes in a variety of different physical systems, ranging from atomic physics, plasma physics, to cybernetical physics (see [16][17][18][19][20][21] and references therein). In the nonlinear control theory the Kapitza pendulum is used as an example of a parametric oscillator that demonstrates the concept of "dynamic stabilization".

In this paper we study the dynamics of a ϕ 0 -junction with direct coupling between the magnetic moment and Josephson current [22,23]. We find that the application of an external voltage may lead to different stability regions for the magnetization, i.e., Kapitza pendulum features in the ϕ 0 -junction. Developing our previous approach [24], we investigate the effect of superconducting current on the dynamics of the magnetic moment. We show that, starting with magnetization along the z-axis, the character of the m dynamics changes crucially as a stable position of m develops between the z-and y-axes, depending on the parameters of the system.

The ac Josephson effect provides an ideal tool to study magnetic dynamics in a ϕ 0 -junction. To realize the ac Josephson effect, we apply a constant voltage V to the ϕ 0 -junction. In this case the superconducting phase varies with time, as ϕ(t) = ω J t, where ω J = 2 eV/h is the Josephson frequency [23,[START_REF] Josephson | Superconductivity[END_REF]. The dynamics of the considered system is then described by the Landau-Lifshitz-Gilbert equation

dM dt = γH eff × M + α M 0 M × dM dt , (1)
with effective magnetic field H eff in the form [23]

H eff = K M 0 Γ sin ωt -ϕ 0 y + M z M 0 z , ( 2
)
where γ is the gyromagnetic ratio, α is the phenomenological damping constant, M 0 = M , and M i are the components of M. Below we use the normalized components m i = Mi M0 . Furthermore, ϕ 0 = r My M0 , r = lυ so /υ F , l = 4hL/hυ F , L is the length of the F -layer, h is the exchange field of the F layer, Γ = Gr, G = E J /(KV), and E J = Φ 0 I c /2π is the Josephson energy. Here Φ 0 is the flux quantum, I c is the critical current, υ F is the Fermi velocity, and the parameter υ so /υ F characterizes a relative strength of spin-orbit interaction, while K is the anisotropic constant, and V is the volume of the F layer. To investigate the dynamics of considered system numerically, we write eq. ( 1) in the dimensionless form (see system of equations ( 1) in the Supplementary Material Supplementarymaterial.pdf (SM)). In that system the time is normalized to the inverse ferromagnetic resonance frequency ω F = γK/M 0 : (t → tω F ), so ω is normalized to the ω F .

First we present results of numerical simulations of the system (1) and (2). Figure 1 shows the dynamics of the magnetization components m z and m y at different values of the parameter G, demonstrating the re-orientation of the oscillations around the z-axis to the oscillations around the y-axis. With increase in G, the component m y goes from zero to m y = 1. and in the middle of time interval is shown in the insets. We see that the average value of m z deviates from one. Figure 1(b) shows the corresponding oscillation of m y . Figure 1(c) shows a stabilization of m y oscillations around some average value of m y , between the z-and y-directions with increase in G, i.e., it shows the re-orientation of the oscillations at three values of G = 10π, 20π, 50π. With an increase in G, the time of re-orientation from the z-direction to the y-direction decreases substantially. At a sufficiently large value of G, the average value of m y becomes close to 1, as we see in fig. 1(c We note that the effect of r leads to similar features as we observe for changes in G. So, the situation is reminiscent of the Kapitza pendulum, i.e., a pendulum whose point of suspension vibrates, leading to the inversion of the stability positions of the pendulum [START_REF] Kapitza | [END_REF].

Variation of average m y as a function of frequency ω and G is shown in fig. 2. We see that an increase in G makes orientation of m y along the y-axis stable, but that the frequency dependence differs from the characteristic Kapitza pendulum behavior, as we can see in fig. 4. The system of equations ( 1), describing the dynamics of the magnetic moment in the angular variables m z = cos θ, m x = sin θ cos φ, m y = sin θ sin φ, can be written as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ φ = cos θ 1 + α 2 - Gr 1+α 2 1 sin θ [cos θ sin φ -α cos φ] sin [ωτ -r sin θ sin φ], θ = - α sin 2θ 2(1 + α 2 ) + Gr 1 + α 2 [α cos θ sin φ + cos φ] sin [ωτ -r sin θ sin φ].
(3) In the case of no coupling between the magnetic moment and the Josephson junction, at G = 0, we have

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ φ = cos θ 1 + α 2 , θ = - α sin 2θ 2(1 + α 2 ) . ( 4
)
The solution to eq. ( 4) has the form

θ(τ ) = ArcTan exp log tan θ 0 - ατ 1 + α 2 , ( 5
)
where tan θ 0 is considered to be positive for simplicity. So, the characteristic scale in time to order the magnetic moment along the easy axes is τ = 1+α 2 α . To realize the Kapitza pendulum we should investigate the case ω 2π τ = 2πα 1+α 2 , where ω is the frequency of the external fastvarying field.

Let us make the replacements

θ → Θ + ξ, φ → Φ + η. ( 6
)
Here Θ and Φ describe slow movement, while ξ and η are coordinates for fast-varying movement. Conditions for ξ and η are discussed in the SM. We consider θ = Θ, φ = Φ, where averaging is taken over the period of the fast-varying force T = 2π ω . The system of equations for slow movement has a form (see SM)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Φ = cos Θ 1 + α 2 - (Gr) 2 rα 2ω(1 + α 2 ) 2 • 1 sin Θ [cos Θ sin Φ -α cos Φ] 1 -sin 2 Θ sin 2 Φ , Θ = - α sin 2Θ 2(1 + α 2 ) + (Gr) 2 rα 2ω(1 + α 2 ) 2 [α cos Θ sin Φ + cos Φ] 1 -sin 2 Θ sin 2 Φ . (7)
There are some equilibrium points of the system following from Θ = 0, Φ = 0.

The first pair is located on the equator

Θ 0 = π/2, Φ 0 = π/2 and Θ 0 = π/2, Φ 0 = 3π/2. ( 8
)
Linearization of ( 7) after substitution Φ = Φ 0 + δφ,

Θ = Θ 0 + δθ gives ⎧ ⎪ ⎨ ⎪ ⎩ δ φ = 1 1 + α 2 δθ, δθ = α 1 + α 2 δθ. (9)
It is clear that δθ ∼ e ατ 1+α 2 , so these points are unstable. Using (7) we find that the equilibrium points at Φ 0 = π 2 are described by equation

sin Θ 0 = -1 + 1 + 4β 2 2β , ( 10
)
where β = (Gr) 2 rα 2ω(1+α 2 ) , which can be approximated as sin Θ 0 = β at small β. It has two solutions in the interval 0 ≤ Θ 0 ≤ π. Note that at Φ 0 = 3π 2 we have sin

Θ 0 = 1- √ 1+4β 2 2β
which leads to the negative sin Θ 0 , but Θ is always positive. So, there are not any stable points at Φ 0 = 3π 2 . To find out if these points are stable we have linearized eq. ( 7) near the points, and tested the eigenvalues of the linearized system. This straightforward calculation shows that the real parts of eigenvalues of the corresponding system are always negative. It means that the second pair of points are always stable (see SM). Figure 3 shows the location of the unstable and stable pairs of equilibrium points.

Results of numerical calculations of the averaged m y as a function G at different frequencies are presented in fig. 4. We see that system reminiscent the Kapitza pendulum behavior: the averaged m y component characterizing changes of stability direction is growing with G. The character of behavior depends essentially on the frequency of the fast-varying field getting very sharp at small ω. Figure 4 also compares the m y (G) dependence obtained by analytical and numerical calculations. We see results at three frequencies: ω = 0.5, 20 and 70. The analytical dependence is calculated according to formula (10). Two methods are used for numerics: numerics-1 is based on the standard program "Mathematica" (shown by triangles, calculated for ω = 70) and numerics-2 presents results of the direct solution of system (1) based on the Runge-Kutta fourth order. Both numerics give practically the same results. An excellent agreement between analytical and numerical results is obtained at low G, depending on the frequency of the fast-varying field. For high frequencies the coincidence is still good at rather large G.

Table 1 shows the values of the averaged m y obtained by analytic and numerical calculations at different G. The difference in values of the averaged m y obtained by all the three used methods is rather small. There is an essential difference between the original Kapitza pendulum and our system. In the Kapitza pendulum the stability of a new equilibrium point is determined by the amplitude and frequency of the external force. In our case, two new points are always stable and their positions on the unit sphere is determined by the parameters of the system. Figure 5 shows the vector fields in the plane 0 ≤ Φ ≤ 2π, 0 < Θ < π according to eq. ( 7) at two different values of G: 5π and 50π. This demonstrates that with an increase in G the two equilibrium positions approach the unstable point at Φ = π 2 , Θ = 0. This result is in agreement with (10), which shows that an increase in G leads to an increase in β, and to the convergence of sin Θ 0 to one.

There is another interesting phenomenon which is realized in our case. With changing parameters of the system, the stable equilibrium positions of the "slow" movement are approaching the unstable one at Φ = π 2 , Θ = π/2. When the distance between them becomes comparable to the amplitude of the "fast" motion, all the three special points of the "slow" system effectively merge.

Finally, we specify the set of ferromagnetic layer parameters and junction geometry that make it possible to reach the values used in our numerical calculations, for the possible experimental observation of the predicted effect. In our model the interaction between the Josephson current and the magnetization is determined by the parameter G = E J /K * V , which describes the ratio between the Josephson energy and the magnetic anisotropy energy. To have a strong interaction, one needs to choose a ferromagnet with weak magnetic anisotropy, like, for example, the permalloy with K ∼ 4 * 10 -5 K Å-3 [START_REF] Rusanov | [END_REF], and a junction with a relatively high critical current. Recently, a high critical current density of (3 * 10 5 -5 * 10 6 ) A/cm 2 was reported for a Nb/NiGdNi/Nb junction [27]. For the volume V of the ferromagnetic film V = 10 -15 cm 3 (10 -5 cm × 10 -4 cm × 10 -6 cm), we find G ∼ 20-200, for these parameters.

The strength of the specific Rashba-type spin-orbit interaction in our case is controlled by the parameter Γ = Gr, which should be stronger in metals with large atomic number Z. A suitable candidate may be a permalloy doped with Pt [28]. In Pt the spin-orbit interaction plays a crucial role in electronic band formation and the parameter v so /v F ∼ 1. The doping of the permalloy by Pt, up to 10 percent, did do not change the magnetic properties significantly, and therefore we may expect to have v so /v F ∼ 0.1, in this case. If the length of the ferromagnetic layer is of the order of the magnetic decaying length, i.e., l ∼ 1, we have r ∼ 0.1, the value of the parameter Γ is in the range 2-20. Note that ferromagnets without inversion symmetry, like MnSi or FeGe, may also be interesting candidates for the observation of the discussed phenomena.

In conclusion, we have demonstrated that the direct coupling between the magnetic moment and Josephson current in a ϕ 0 -junction may effectively lead to strong re-orientation of the magnetic easy axis in response to an applied voltage. This serves as a manifestation of Kapitzalike pendulum behavior and opens a new way to control the magnetization via a superconducting current.

Figure 1 (Fig. 1 :

 11 Fig. 1: (Color online) (a) Dynamics of the mz component at G = 5π, r = 0.1; (b) the same for my; (c) effect of G at r = 0.5. The numbers show the G value.

  ), for case G = 50π. The oscillations show periodically splashing related to the Josephson frequency. The inset demonstrates the jumps of m y from the value m y = 1 with Josephson frequency. The amplitudes of the jumps decrease in time, becoming smaller and shorter in time, with an increase in G.

Fig. 2 :

 2 Fig. 2: (Color online) The ω-G diagram for averaged my shown by the color scale.

Fig. 3 :Fig. 4 :

 34 Fig. 3: (Color online) Demonstration of the equilibrium points of eq. (7). Arrows show stable points, the other two points are unstable.

Fig. 5 :

 5 Fig. 5: (Color online) Phase planes (7) at G = 10π (left) and G = 50π (right). Here r = 0.5, α = 1, ω = 70. Red points are stable equilibrium points which are calculated from (10), black points are unstable(8). With increase in G the stable red points are approaching unstable ones Φ0 = π/2, Θ0 = π/2.

Table 1 :

 1 m y from the theory and the numerical results.

	G	Θ 0	analytic m y numerics-1 numerics-2
	15.7 1.462	0.109	0.108	0.109
	31.4 0.387	0.369	0.376	0.378
	47.1 0.663	0.577	0.611	0.616
	62.8 0.859	0.689	0.748	0.756
	157 1.272	0.817	0.945	0.956

The authors thank K. Sengupta, K. Kulikov, I. Bobkova, and A. Bobkov for useful discussions. The study was partially funded by the RFBR (research project 18-02-00318), and the SA-JINR Collaborations. YuMS and AEB gratefully acknowledge support from the University of South Africa's visiting researcher program. AB acknowledges the French ANR projects "SUPERTRON-ICS" and Leverhulme Trust for supporting his stay in Cambridge.


Julian Krebs

Tommaso Mansi

Boris Mailhé

Nicholas Ayache

Hervé Delingette

Unsupervised Probabilistic Deformation Modeling for Robust Diffeomorphic Registration

We propose a deformable registration algorithm based on unsupervised learning of a low-dimensional probabilistic parameterization of deformations. We model registration in a probabilistic and generative fashion, by applying a conditional variational autoencoder (CVAE) network. This model enables to also generate normal or pathological deformations of any new image based on the probabilistic latent space. Most recent learning-based registration algorithms use supervised labels or deformation models, that miss important properties such as diffeomorphism and sufficiently regular deformation fields. In this work, we constrain transformations to be diffeomorphic by using a differentiable exponentiation layer with a symmetric loss function. We evaluated our method on 330 cardiac MR sequences and demonstrate robust intrasubject registration results comparable to two state-of-the-art methods but with more regular deformation fields compared to a recent learningbased algorithm. Our method reached a mean DICE score of 78.3% and a mean Hausdorff distance of 7.9mm. In two preliminary experiments, we illustrate the model's abilities to transport pathological deformations to healthy subjects and to cluster five diseases in the unsupervised deformation encoding space with a classification performance of 70%.

Introduction

Deformable registration is an essential task in medical image analysis. It describes the process of finding voxel correspondences in a pair of images [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Traditional registration approaches aim to optimize a local similarity metric between deformed and target image, while being regularized by various energies [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. In order to retrieve important properties such as invertible deformation fields, diffeomorphic registration was introduced. Among other parametrizations, one way to parametrize diffeomorphisms are stationary velocity fields (SVF) [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF].

In recent years, major drawbacks of these approaches like high computational costs and long execution times have led to an increasing popularity of learningbased algorithms -notably deep learning (DL). One can classify these algorithms as supervised or unsupervised. Due to the difficulty of finding ground truth voxel correspondences, supervised methods need to rely on predictions from existing algorithms [START_REF] Yang | Fast predictive image registration[END_REF], simulations [START_REF] Sokooti | Nonrigid image registration using multi-scale 3d convolutional neural networks[END_REF] or both [START_REF] Krebs | Robust non-rigid registration through agent-based action learning[END_REF]. These methods are either limited by the performance of the used existing algorithms or the realism of simulations. On the other hand, unsupervised approaches make use of spatial transformer layers (STN [START_REF] Jaderberg | Spatial transformer networks[END_REF]) to warp the moving image in a differentiable way such that loss functions can operate on the warped image (similarity metric) and on the transformation itself (regularization) [START_REF] Balakrishnan | An unsupervised learning model for deformable medical image registration[END_REF][START_REF] Jason | Back to basics: Unsupervised learning of optical flow via brightness constancy and motion smoothness[END_REF][START_REF] De Vos | End-to-end unsupervised deformable image registration with a convolutional neural network[END_REF]. While unsupervised approaches perform well in minimizing a similarity metric, it remains unclear if the retrieved deformation fields are sufficiently regular which is of high interest for intra-subject registration. Furthermore, important properties like symmetry or diffeormorphisms [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] are still missing in DL-based approaches.

In this paper, we suggest to learn a low-dimensional probabilistic parameterization of deformations which is restricted to follow a prescribed distribution. This stochastic encoding is defined by a latent code vector of an encoder-decoder neural network and it restricts the space of plausible deformations with respect to the training data. By using a conditional variational autoencoder (CVAE [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]), our generative network constrains encoder and decoder on the moving image. After training, the probabilistic encoding can be potentially used for deformation analysis tasks such as clustering of deformations or the generation of new deformations for a given image -similar to the deformations seen during training. Furthermore, we include a generic vector field exponentiation layer to generate diffeomorphic transformations. Our framework contains an STN and can be trained with a choice of similarity metrics. To avoid asymmetry, we use a symmetric local cross correlation criterion. The main contributions are:

• A probabilistic formulation of the registration problem through unsupervised learning of an encoded deformation model. • A differentiable exponentiation and an user-adjustable smoothness layer that ensure the outputs of neural networks to be regular and diffeomorphic. • As a proof of concept, first experiments on deformation transport and disease clustering.

Methods

The goal of image registration is to find the spatial transformation T z : R 3 → R 3 , parametrized by a d-dimensional vector z ∈ R d , which best warps the moving image M to match the fixed image F. Both images are defined in the spatial domain Ω ∈ R 3 . Typically, this is done by minimizing an objective function of the form: arg min z F(z, M, F) = D (F, M • T z ) + R(T z ) with the image similarity D of the fixed F and the warped moving image M • T z and a spatial regularizer R. Recent unsupervised DL-based approaches (e.g. [START_REF] Balakrishnan | An unsupervised learning model for deformable medical image registration[END_REF][START_REF] Jason | Back to basics: Unsupervised learning of optical flow via brightness constancy and motion smoothness[END_REF]) mimic the optimization of such an objective function.

Instead, we propose to model the registration probabilisitcally by parametrizing the deformation as a vector z to follow a prior p(z). To learn this probabilistic space, we define the latent vector of dimensionality d in an encoder-decoder neural network as this z. Given the moving and the fixed image as input, a variational inference method (CVAE [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]) is used to reconstruct the fixed by warping the moving image. An exponentiation layer interprets the network's output as velocities v (an SVF) and returns a diffeomorphism φ which is used by a dense STN to retrieve the warped image M * . To enforce an user-adjustable level of deformation smoothness (comparable to [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF]), a convolutional Gaussian layer is added before the exponentiation with Gaussian weights according to the variance σ 2 S . During training, the network parameters are updated through back-propagation of the gradients. The network architecture can be seen in Fig. 1a. Finally, registration is done in a single forward path. The trained probabilistic framework can be also used for the sampling of deformations as shown in Fig. 1b.

Learning a Probabilistic Deformation Encoding Learning a generative model typically involves a latent variable model (as in VAE)

, where an encoder maps an image to its z-code -a low-dimensional latent vector, from which a decoder aims to reconstruct the original image. Typically, the encoder and decoder are defined as distributions q ω and p γ with trainable network parameters ω and γ. The network is trained by maximizing a lower bound on the data likelihood with respect to a prior distribution p(z). We define the prior as multivariate unit Gaussians p(z) = N (0, I) with the identity matrix I. In CVAE [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF], encoder q ω and decoder p γ distributions are additionally conditioned on extra information (e.g. classes). We propose to frame image registration as a reconstruction problem in which the moving image M acts as the conditioning data and is warped to reconstruct or to match the fixed image F. Thus, the decoder reconstructs F given z and M: p γ (F | z, M). To have z, the encoder serves as an approximation of the intractable true posterior probability of z given F and M and is denoted as q ω (z | F, M). Since the prior p(z) is defined as multivariate unit Gaussians, the encoder network predicts the mean µ ∈ R d and diagonal covariance σ ∈ R d , from which z is drawn:

q ω (z | F, M) = N (µ(F, M), σ(F, M)).
Both distributions can be combined in a two-term loss function [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF] where the first term describes the reconstruction loss as the expected negative loglikelihood of p γ (F | z, M). In other words, the reconstruction loss represents a similarity metric between input F and output M * . The second term acts as a regularization term on the deformation latent space by forcing the encoded distribution q ω (z | F, M) to be close to the prior probability distribution p(z) using a Kullback-Leibler (KL) divergence. The loss function results in:

l(ω, γ, F, M) = -E z∼qω(•|F,M) [logp γ (F | z, M)]+KL [q ω (z | F, M) p(z)] , (1)
where the KL-divergence can be computed in closed form [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]. Assuming a Gaussian log-likelihood term of p γ is equivalent to minimizing a weighted SSD criterion (cf. [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]). We propose instead to use a symmetric local cross-correlation (LCC) criterion due to its favorable properties for registration [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF] and assume a LCC Boltzmann distribution p γ (F | z, M) ∼ exp(-λD LCC (F, M, v)) with the LCC criterion D LCC and the weighting factor λ. Using the velocities v and a small constant , which is added for numerical stability, we define:

D LCC (F, M, v) = 1 P x∈Ω F x • exp -vx 2 M x • exp vx 2 2 F x • exp -vx 2 2 M x • exp vx 2 2 + , (2)
with a total number of P pixels x ∈ Ω and where • symbolizes the local mean image derived by Gaussian smoothing with a strength of σ G and kernel size k. To help the reconstruction task, we introduce conditioning by involving M not only as the image to be warped in the STN, but also in the first decoder layers by concatenating down-sampled versions of M with the filter maps on each scale. The hypothesis is that in order to better optimize the reconstruction loss, the network makes use of the provided extra information of M such that less anatomical but more deformation information are conveyed by the low-dimensional latent layer, which would make the encoding more geometry-invariant.

Exponentiation Layer: Generating Diffeomorphisms In the SVF setting, the transformation φ is defined as the Lie group exponential map with respect to the velocities v: φ(x) = exp(v). For efficient computation, the scaling and squaring algorithm is typically used [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF]. In order to generate diffeomorphic transformations φ in a neural network, we propose an exponentiation layer that implements this algorithm in a fully differentiable way. To this end, the layer expects a vector field as input (the velocities v) which is scaled with a factor N which we precompute on a subset of the training data according to the formulations in [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF]. In the squaring step, the approximated φ 0 ≈ id + v * 2 -N (with id as a regular grid) is recursively squared, N -times, from k = 1 to N :

φ k = φ k-1 • φ k-1 .
The result is the diffeomorphism φ N ≡ φ [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF]. The squaring step requires the composition of two vector fields on regular grids which we realized by linear interpolation. All these computations consist of standard operations that can be added to the computational graph and are auto-differentiable in modern deep learning libraries. This differentiable layer can be added to any neural network which predicts (stationary) velocity fields.

Experiments

We evaluate our framework on an intra-subject task of cardiac MRI cine registration where end-diastole frames are registered to end-systole frames (ED-ES) -a very large deformation. Furthermore, we show preliminary experiments evaluating the learned deformation encoding: its potentials for transporting encoded deformations from one subject to another and showing the clustering of diseases in the encoding space. All experiments are in 3-D. We used 184 short-axis datasets acquired from different hospitals and 150 cases from the Automatic Cardiac Diagnosis Challenge (ACDC) at STACOM 2017 3 , mixing congenital heart diseases with images from adults. We used 234 cases for training and for testing the remaining 100 cases from ACDC, that contain segmentation and disease label information from five cardiac diseases. Both information were only used for evaluation purposes. All images were sampled with a spacing of 1.5 × 1.5 × 3.15 mm and cropped to a size of 128 × 128 × 32 voxels. These dimensions were chosen to save computation time and are not a limitation of the framework (validated on different image sizes).

Un Dem VM Our x10 -1 x10 -2 x10 -4
Implementation Details The encoder of our neural network consisted of four convolutional layers with strides (2, 2, 2, 1) (Fig. 1a). The bottleneck layers (µ, σ, z) were fully-connected. The decoder had one fully-connected and three deconvolutional layers, where the outputs at each layer were concatenated with sub-sampled versions of M. Two convolutional layers and a convolutional Gaussian layer with σ S = 3 (kernel size 15) were placed in front of the exponentiation and transformer layer. The latent code size d was set to 16 as a trade off between registration quality and generalizability. This leads to a total of ∼267k trainable parameters. L2 weight decay with a factor of 0.0001 was applied. The numbers of iterations in the exponentiation layer was set to N = 4 in all experiments. In training, the strength of the Gaussians for computing the LCC was set to σ G = 2 with a kernel size k = 9. The loss balancing factor λ = 5000 was empirically chosen such that encoded training samples roughly had zero means and variances of 1 and the reconstruction loss was optimized. We used the Adam optimizer with a learning rate of 0.0005 and a batch size of one. We performed online data augmentation by randomly shifting, rotating, scaling and mirroring training images. The framework has been implemented using Keras with Tensorflow. Training took 24 hours on a NVIDIA GTX TITAN X GPU.

Registration Results

We compare our registration algorithm with the LCCdemons [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF] with manually tuned parameters (on training images) and the nondiffeomorphic DL-based method VoxelMorph-2 [2] (VM) with a regularization weighting parameter of 1.5, as recommended. As a surrogate measure of registration performance, we used the intensity root mean square error (RMSE), mean DICE score and 95%-tile Hausdorff distance (HD) in mm on the following anatomical structures: myocardium and epicardium of the left ventricle (LV-Epi, LV-Myo), left bloodpool (BP) and heart (Heart). The LCC-demons showed better mean DICE scores (averaged over the five structures, in %) with 79.9 compared to our algorithm with 78.3 and VM with 77.5 (cf. Fig. 2). The Voxel-Morph algorithm reached a very low RMSE of 0.025 compared to ours (0.031) and the demons (0.034), but could not reach the other algorithms in terms of HD with a mean score of 9.4mm compared to ours with 7.9mm and the demons with 8.2mm. Besides these metrics, VM produced very irregular and highly nondiffeomorphic deformation fields since 2.2% of the displacements had a negative Jacobian determinant (cf. in Fig. 3). In general, our approach led to deformation fields with both smaller amplitudes and smaller gradients than the demons and the VM algorithm. Furthermore, our results were more robust as variances were lower for all metrics compared to the demons and lower or comparable to VM. This is also visible in Fig. 3 and further shown by the fact that HD scores are the smallest experienced in the experiments. Average execution time per test case was 0.32s using the mentioned GPU and an Intel Xeon CPU E5, compared to 108s for the demons on CPU. Deformation Encoding For evaluating the learned deformation encoding, we show geometry-invariance by transporting a deformation from one subject to another. Therefore, we take a z-code from a pathological subject and condition the decoder on the ED image of healthy subjects (Fig. 1b). More precisely, in Fig. 4 we transported a cardiomyopathy (DCM) and hypertrophy (HCM) deformation to two healthy cases (Normal). One can see the disease-specific deformation (DCM: reduced cardiac contraction) which are different from the healthy transformations. The resulting deformation fields are adapted to the anatomy of the conditioning image and they are translation-invariant. In a second experiment, we used the encoded z-codes and disease information of our cardiac test set to visualize the structure of the learned space. Therefore, we linearly projected the 16-D z-codes to a 2-D space by using the two most discriminative CCA components (canonical correlation analysis). We used the ACDC classes: dilated cardiomyopathy DCM, hypertrophic cardiomyopathy HCM, myocardial infarction MNF, abnormal right ventricle RV and normal NOR. In Fig. 5, one can see that the classes of the 100 test sets are clustered in the projected space. The five class classification accuracy reaches 70% with 10fold cross-validation, by using the six most discriminative CCA components and applying support vector machine (SVM) on-top. These results which are solely based on unsupervised deformation z-codes suggest that similar deformations are close to each other in the deformation encoding space.

Conclusion

We presented an unsupervised deformable registration approach that learns a probabilistic deformation encoding. This encoding constrains the registration and leads to robust and accurate registration results on a large dataset of cardiac images. Furthermore, an exponentiation layer has been introduced that creates diffeomorphic transformations. The performance of the proposed method was comparable and partially superior to two state-of-the-art algorithms. Our approach produced more regular deformation fields than a DL-based algorithm. Furthermore, first results show, that the probabilistic encoding could potentially be used for deformation transport and clustering tasks. In future work, we plan to further explore the deformation encoding to evaluate these tasks more deeply.

Fig. 1 :

 1 Fig. 1: (a) CVAE registration network during training and registration including diffeomorphic layer (exponentiation). Deformations are encoded in z from which velocities are decoded while being conditioned on the moving image. (b) Decoder network for sampling and deformation transport: Apply z-code conditioned on any new image M.

Fig. 2 :

 2 Fig. 2: Comparing registration performance: unregistered (Un), LCC-Demons (Dem), VoxelMorph (VM) and our method in terms of RMSE and mean deformation magnitude and gradient, DICE and 95%-tile Hausdorff distances (HD).

Fig. 3 :

 3 Fig. 3: Two random examples of end-diastole to end-systole registration: (Row 1) original images. The LCC-demons (Dem, Row 2) and VoxelMorph (VM, Row 3) versus our method (Row 4), showing the warped moving image, the deformation field and the Jacobian determinants. All results are in 3-D, showing the central short-axis slices.

Fig. 4 :

 4 Fig. 4: Transport the z-code of pathological deformations (top row: cardiomyopathy DCM and hypertrophy HCM) to two healthy subjects (bottom rows: Normal). The simulated deformation fields are similar compared to the pathological deformations but are adapted to the geometry of the healthy image (e.g. translated).

Fig. 5 :

 5 Fig. 5: Distribution of cardiac diseases after projecting 100 z-codes of test images on 2 CCA components.

https://www.creatis.insa-lyon.fr/Challenge/acdc/index.html

Acknowledgements: Data used in this article were obtained from the EU FP7-funded project MD-Paedigree and the ACDC STACOM challenge 2017. Disclaimer: This feature is based on research, and is not commercially available. Due to regulatory reasons its future availability cannot be guaranteed.