K Makki

B Borotikar

M Garetier

O Acosta

S Brochard

D Ben Salem

F Rousseau

4D in vivo quantification of ankle joint space width using dynamic MRI

Keywords: intensity-based registration, joint space width (JSW), partial differential equations (PDEs). I

Spatio-temporal evolution of joint space width (JSW) during motion is of great importance to help with making early treatment plans for degenerative joint diseases like osteoarthritis (OA). These diseases can affect people of all ages leading to an acceleration of joint degeneration and to limitations in the activities of daily living. However, only a few studies have attempted to quantify the JSW from moving joints. In this paper, we present a generic pipeline to accurately determine the changes of the JSW during the joint motion cycle. The key idea is to combine spatial information of static MRI with temporal information of low-resolution (LR) dynamic MRI sequences via an intensity-based registration framework, leading to a high-resolution (HR) temporal reconstruction of the joint. This allows the temporal JSW to be measured in the HR domain using an Eulerian approach for solving partial differential equations (PDEs) inside a deforming inter-bone area where the HR reconstructed bone segmentations are considered as temporal Dirichlet boundaries. The proposed approach has been applied and evaluated on in vivo MRI data of five healthy children to non-invasively quantify the spatio-temporal evolution of the JSW of the ankle (tibiotalar joint) during the entire dorsi-plantar flexion motion cycle. Promising results were obtained, showing that this pipeline can be useful to perform large-scale studies containing subjects with OA for different joints like ankle and knee.

Index Terms-Dynamic MRI,

overlap between the two contact cartilage layers under the assumption that each cartilage layer should undergo the same rigid transform as that of the attached bone. However, it was needed to manually segment the dynamic data as they used mesh-based registration for reconstructing HR temporal bone models. Moreover, this manual intervention is time consuming and prone to segmentation bias because of noise and low-resolution of dynamic data. Another limitation of this method is that it offers limited spatial information (only overlapped regions are highlighted). In summary, tracking the entire joint space area and performing accurate measures of the JSW throughout a specific joint trajectory are still challenging because articular cartilages undergo non-linear deformations during motion. In the context of soft-tissue thickness quantification, Yezzi et al. [START_REF] Yezzi | An Eulerian PDE approach for computing tissue thickness[END_REF] presented an Eulerian framework for computing the thickness of tissues based on the solution of Laplace's equation between two simply connected boundaries. These techniques have been successfully employed to measure the cortical thickness [START_REF] Acosta | Automated voxel-based 3D cortical thickness measurement in a combined Lagrangian-Eulerian PDE approach using partial volume maps[END_REF] and to segment the articular cartilages from conventional MRI of the knee joint [START_REF] Fripp | Automatic segmentation and quantitative analysis of the articular cartilages from magnetic resonance images of the knee[END_REF]. In this paper, we extend these techniques to the measurement of the JSW during the entire joint motion cycle using dynamic MRI sequences, which will certainly provide new temporal motion features. In our previous work [START_REF] Makki | High-resolution temporal reconstruction of ankle joint from dynamic MRI[END_REF], we have proposed a motion-based algorithm for estimating HR dynamic MRI sequences via a log-euclidean polyrigid framework (LEPF) [START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF] using both static and dynamic MRI. In the current work, we use the associated motion estimation algorithm to track the tibiotalar joint with high accuracy and without the need for manual segmentation of dynamic data. The main goal is to identify the spatial location of each bone and of the joint space section at each time frame, but in the HR domain of the static scan (i.e. image resolution of 0.26 × 0.26 × 0.8mm). So that the temporal JSW measures are performed in this HR domain.

To summarize, we present a complete pipeline for : 1) estimating spatio-temporal bone rigid transformations (Section III-A), 2) computing spatio-temporal non-rigid deformations of the joint space area using the LEPF (Section III-B), 3) measuring dynamic, in vivo JSW via an iterative relaxation method (3) using the high-resolution-reconstructed bone segmentations as Dirichlet boundaries (Section III-C).

II. MRI SCANNING

Data were acquired as part of a pilot study for five healthy children (range 10-14 years) at CHRU of Brest, France, under a protocol approved by the Institutional Review Board (IRB). Children were selected with no contraindications to MRI and with no history of lower limb musculo-skeletal injury or surgery in the past six months in order to establish normative contact mechanics before focusing on pathomechanics studies. HR images of the ankle were obtained using a 3D T1-weighted GE sequence (TR= 7.81ms, TE= 2.75ms, FA= 10 • , FOV 150mm × 150mm, voxel size of 0.26 × 0.26 × 0.8mm, and matrix 576 × 576), using a 3T clinical MR imaging scanner (Achieva dStream, Philips Medical Systems, Best, NL). Ankle passive movements were standardized for all subjects using a specifically designed orthotic fixture in order to perform coherent averaging of the results.

III. METHODS

A. Bone motion estimation

The proposed algorithm for bone motion tracking consists in estimating the rigid transformation of each bone between the HR static image and each LR dynamic image. This leads to a high-resolution temporal reconstruction of the ankle joint from LR dynamic MRI sequence. This algorithm require only the manual segmentations of bones of interest in the high-resolution static image as inputs. For more details about this algorithm, the reader is referred to our previous works [START_REF] Makki | High-resolution temporal reconstruction of ankle joint from dynamic MRI[END_REF], [START_REF] Makki | In vivo ankle joint kinematics from dynamic magnetic resonance imaging using a registration-based framework[END_REF].

B. Temporal tracking of joint space area

Following previous works [START_REF] Goker | Quantification of the radiographic joint space width of the ankle[END_REF], [START_REF] Imai | Joint space width of the tibiotalar joint in the healthy foot[END_REF], we define the tibiotalar joint space in healthy feet as the space between the upper surface of the talus and the inferior articular surface of the tibia. The joint space area is manually segmented in the static image (figure 1) and then automatically tracked throughout the dynamic sequence using non-rigid registration. In this work, we propose to compute the local spatio-temporal deformations of this area using the LEPF. The segmentations of the two bones, composing the joint of interest, are used to compute the spatial weighting functions {w i (x)} i∈{1,2} . The use of this bijective framework ensure the estimation of a smooth and invertible transformation at each voxel x inside the joint space section R s by integrating its trajectories (from the static image S to each k th time-frame D k ) in the domain of matrix logarithms, via the following equation :

T S→D k (x) = exp w1 S (x). log(T 1 k ) + w2 S (x). log(T 2 k )).x (1
) where : T i k is the rigid transform from S to D k for the bone i ; wi S (x) is the normalized weighting function for the bone i in S (i.e. w1 S (x) + w2 S (x) = 1∀x ∈ R s , R s ⊂ S). To compute the exponential mapping of Equation 1 in regular grids with high-accuracy, we propose to use the matrix eigendecomposition * instead of the classical use of the scaling and squaring method involving approximation theory where the accuracy always remains debatable [START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF].

C. Computation of temporal joint space width

Given the position of each bone (B 1 k and B 2 k ) and each joint space area (R k ) at each time frame D k , the timedependent width of the area separating two neighboring bones can be computed with high accuracy using an Eulerian PDE approach [START_REF] Yezzi | An Eulerian PDE approach for computing tissue thickness[END_REF]. The JSW is computed by first solving Laplace's equation ( u k = 0) inside R k . For this, we have used the following Jacobi iterative method which is simple to implement and numerically robust :

u i+1 k (x, y, z) = 1 6 [u i k (x + 1, y, z) + u i k (x -1, y, z) + u i k (x, y +1, z) + u i k (x, y -1, z) + u i k (x, y, z + 1) + u i k (x, y, z -1)] (2)
where u i k (x, y, z) is the potential of the point p = (x, y, z) ∈ R k during the i th iteration. Bones are considered as Dirichlet boundaries so that they are set to fixed potentials (u k (B 1 k ) = 0V for the inner boundary and u k (B 2 k ) = 100V for the outer boundary). The potential is set to zero at all grid points outside the joint segmentation (including R k ) to connect the two boundaries as shown in figure 5.a in such a way that these values will serve as inner boundary condition. Solving this second order PDE gives the harmonic interpolant u k which divides R k into a set of equipotential sub-layers (figure 5.a). Figure 2 shows the solution of Laplace's equation inside a synthetic sphere using the Jacobi iterative method (2). The normalized gradient vector field N k = ∇u k ∇u k is then computed from u k with 3D finite differences. From this normalized gradient vector flow, it is then possible to compute the two correspondence trajectories between p and each of the two boundaries, that we note L 0 k and L 1 k , by solving a pair of PDEs in a parallel way : ∇L 0 . N = -∇L 1 . N = 1.

To simplify notations, we denote L 0 = L 0 k , L 1 = L 1 k , and N = N k . In this work, L 0 and L 1 are first initialized to -(dx+dy+dz) 6 at all grid points and then iteratively updated only over R k using a fast Gauss-Seidel method * :

L i+1 0 [x, y, z] = 1 |Nx|+|Ny|+|Nz| .(1 + |N x | L i 0 [x ∓ 1, y, z]+ |N y | L i 0 [x, y ∓ 1, z] + |N z | L i 0 [x, y, z ∓ 1]) L i+1 1 [x, y, z] = 1 |Nx|+|Ny|+|Nz| .(1 + |N x | L i 1 [x ± 1, y, z]+ |N y | L i 1 [x, y ± 1, z] + |N z | L i 1 [x, y, z ± 1]) (3)
where : dx, dy, and dz are the voxel spacing values in the x, y and z directions, respectively ; m ± 1 = m + sgn(N m ), Fig. 2: 3D solution of the Laplace equation u inside a synthetic sphere S ( i.e. between its center c and its surfaces) using the Jacobi iterative method. The initial Dirichlet boundary conditions are set as follows : u(c) = 0 and u( S) = 100. Three coronal slices are presented to show the harmonic interpolant inside S. From another physical point of view, this can be interpreted as the heat propagation from the outer surface to the sphere center as the Laplace equation is a particular case of the Heat equation.

m ∓ 1 = m -sgn(N m ) for m ∈ {x, y, z} ; sgn(.) is the sign function ; and i is the iteration index. Since L 0 and L 1 may never intersect (as they start from opposite endpoints in opposite bones), the JSW W k (p) can be computed from the tangent field N by adding these two optimal lengths yielding a unique value at each voxel p inside R k (i.e., W k (p) = L 0 (p) + L 1 (p)).

IV. VALIDATION

The accuracies of bone motion estimation have been presented in our previous work using DICE score [START_REF] Makki | High-resolution temporal reconstruction of ankle joint from dynamic MRI[END_REF]. In this work, we have also quantified the accuracy of joint space area tracking. To do this, we have computed the overlap between automatic (or propagated) and ground-truth (or manual) segmentations of the joint space section at each time frame using the Dice coefficient. The temporal Dice scores across all subjects was around 0.9 ± 0.03, indicating that the joint space area was accurately tracked over the entire ankle trajectory. To validate the accuracy of the implemented Eulerian PDE approach, we have generated a sphere S with center c = (64, 64, 64) and radius r = 50 voxels in a regular grid of 128×128×128 voxels (figure 2) and with an isotropic resolution (voxel size of 1×1×1mm), under the assumption that the thickness values should be equal to the radius over the sphere. The obtained thickness values inside S using our implementation are effectively around 50.2 ± 0.4mm.

V. RESULTS

The tibiotalar JSW is computed in 22sec on an Intel ® Xeon ® Processor E3-1271 v3 3.60 GHz on a 576×576×210 grid using our optimized implementation of the Eulerian framework (figure 5). The volume of the tracked joint space area increased in the order of maximum dorsiflexion (183 ± 10cm 3 ), neutral position (190±12cm 3 ), and maximum plantar flexion (196±13cm 3 ). Figure 4 shows the estimated JSW values projected on tibial contact surface for several time

The source code is available on GitHub : https://github.com/ rousseau/dynMRI. frames for one subject. It illustrates the temporal evolution of 3D JSW during continuous passive motion for one subject. An advantage of the use of the Eulerian framework is that it yields width values in a bijective way, as shown in the second row in figure 4. From these temporal maps, we can extract temporal features for each subject to provide an average behavior of the JSW. Assuming the central area to be the center of rotation for the tibiotalar joint in healthy feet [START_REF] Imai | Joint space width of the tibiotalar joint in the healthy foot[END_REF], we define the landmark P 1 shown in figure 3 as the centroid of the joint-space-mesh. For each time-frame for each subject, we then compute the mean average of the JSW inside a sphere centered at P 1 (with radius r = 20voxels 5mm and by taking into account only the voxels inside the joint space area). The medial JSW was around 2.7±0.60mm across all healthy ankles throughout the whole motion cycle. The lateral JSW was around 3 ± 0.60mm. Obtained measures demonstrate correspondance with the JSW measures presented in [START_REF] Imai | Joint space width of the tibiotalar joint in the healthy foot[END_REF] for 10 healthy adults from stationary scans. In this work, the authors did not provide any details about the techniques they have used to measure the JSW making it difficult to perform objective comparisons in terms of measurement errors. Figure 3 shows the temporal mean feature of the JSW for all children during the whole ankle trajectory.

In our study, the JSW in the central part (area centred around point P 1) was not significantly different in any of the positions of the ankle joint while the normal ankle is moving (see Figure 4, thinner orange-yellow areas). In accordance with these results, the ankle moves in this middle position as a centre of rotation for the tibiotalar joint in healthy subjects. This also confirms the results of Imai et al. [START_REF] Imai | Joint space width of the tibiotalar joint in the healthy foot[END_REF] for healthy adults, after dividing the tibial plafond into nine areas. In summary, our methodologies offer many advantages when compared to previous studies such as : reducing the manual interventions as much as possible, tracking the deforming joint space area with high-accuracy, performing 3D width measures over the entire joint space area, in addition to the non-invasive nature of MR imaging.

VI. CONCLUSION

In this paper, we have presented a complete pipeline for measuring spatio-temporal JSW, providing qualitative visual supports and allowing quantitative data analysis. The proposed approach relying on HR reconstructed dynamic sequences can provide 3D+t morphometric information throughout the whole motion. This work is the first attempt to track the entire joint space area and to evaluate JSWs without having to perform time-consuming manual segmentations. The Eulerian PDE approach provides an efficient and fast solution for quantifying widths of segmented areas, bounded by two bone masks. To conclude, dynamic MRI and dedicated image processing techniques can be applied to study both children and adults with degenerative joint diseases. Moreover, they can eventually be applied to other joints and related muskuloskeletal disorders like the knee OA.

  Real-time FFE (Fast Field Echo) sequences have been acquired to reduce the dynamic-scanning time to 18 sec for a single cycle of dorsi-plantar flexion where motion was consistent and continuous during scanning. With TR= 20.61ms, TE= 1.8ms, FA= 15 • , resolution of 0.57×0.57× 8mm, FOV 200mm × 200mm, 6 slices per time frame (15 time frames per sequence) and matrix 352 × 352.

Fig. 1 :

 1 Fig. 1: Segmentation of the tibiotalar joint. Regions of interest are : tibia (green), talus (red), tibiotalar joint space (cyan). From left to right : Mid-coronal image from the HR static image ; 3D rendering of the joint.

Fig. 3 :

 3 Fig. 3: Temporal evolution of the JSW, from left to right : Location of the centroid of the joint space mesh points (P 1 ) ; temporal evolution of the JSW around P 1. The average across subjects is shown as a solid line with one standard deviation above and below this line (grey shaded area).

Fig. 4 :

 4 Fig. 4: Spatio-temporal evolution of the tibiotalar JSW during motion. The first row is composed of a set of dynamic anatomical images from one healthy subject : bone contours in neutral position show the joint motion across time (first row). Width values were projected on tibial contact surface (second row).

Fig. 5 :

 5 Fig. 5: Calculation of tibiotalar JSW : (a) Harmonic function, (b) Length L 0 , (c) Length L 1 , (d) width (L 0 + L 1 ).

*This work was supported by Région Bretagne, Chaire


K Makki

B Borotikar

M Garetier

S Brochard

Ben Salem

François Rousseau

Temporal resolution enhancement of dynamic MRI sequences within a motion-based framework

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Temporal resolution enhancement of dynamic MRI sequences within a motion-based framework K. Makki 1, , B. Borotikar 2 , M. Garetier 3 , S. Brochard 2 , D. Ben Salem 2 , F. Rousseau 1

Abstract-Dynamic MRI has made it possible to noninvasively capture the moving human joints in vivo. Realtime Fast Field Echo (FFE) sequences have the potential to reduce the effect of motion artifacts by acquiring the image data within a few milliseconds. However, the short acquisition times affect the temporal resolution of the acquired sequences. In this paper, we propose a post-processing technique to reconstruct the missing frames of the sequence given the reduced amount of acquired data, which leads to recover the entire joint trajectory outside the MR scanner. To do this, we generalize the Log-Euclidean polyrigid registration framework to deal with dynamic three-dimensional articulated structures by adding the time as fourth dimension : we first estimate the rigid motion of each bone from the acquired data using linear intensity-based registration. Then, we fuse these local transformations to compute the non-linear joint deformations between successive images using a spatio-temporal log-euclidean polyrigid framework. The idea is to reconstruct the missing time frames by interpolating the realistic joint deformation fields in the domain of matrix logarithms assuming the motion to be consistent over a short period of time. The algorithm has been applied and validated using dynamic data from five children performing passive ankle dorsi-plantar flexion.

Index Terms-Dynamic MRI, tracking, intensity-based registration, diffeomorphisms, matrix analysis, eigenvalue

I. INTRODUCTION

Dynamic MRI was originally developed for cardiovascular imaging to quantify blood flow to study heart valve functions towards the end of the 1980s [START_REF] Pettigrew | Dynamic cardiac MR imaging. Techniques and applications[END_REF]. Since then, MRI acquisition and reconstruction techniques have improved by leaps and bounds with current sequences like fast-PC MRI, ultrafast MRI, and real-time FFE. These techniques have been successfully employed to quantify in vivo musculoskeletal biomechanics such as the ankle kinematics during dorsiplantar flexion. Sequences such as Cine-PC MRI or Fast-PC MRI [START_REF] Behnam | Assessing the accuracy and precision of musculoskeletal motion tracking using Cine-PC MRI on a 3.0 Tesla platform[END_REF] can provide velocity changes in the three directions of the space. However, these sequences require repeatable motions of the joint for up to two minutes (e.g. motion rate = 30 cycles/minute for one dorsi plantarflexion cycle) and the dynamic data are built using these repeatable motions. Hence, velocity data does not meet the requirements of pediatric studies. Ultrafast Contrast-Enhanced MRI is another dynamic imaging technique which allows faster scanning speed by exploiting spatio-temporal redundancy and can acquire images during a single dorsi-plantar flexion cycle within 50 seconds [START_REF] Clarke | A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF]. However, these sequences require slow motion of the joint (mean rotation speed of 1 • /s) as they are sensitive to motion artifacts so the motion itself is reduced to almost static or quasi-static nature where no real functionality can be assessed. Real-time FFE sequences [START_REF] Schaeffter | Projection reconstruction balanced FFE for interactive real-time cardiac imaging[END_REF] are less sensitive to motion artifacts and they can reduce the scanning time to only 18 seconds to acquire a single cycle of dorsi-plantar flexion. However, the short acquisition times reduce the number of acquired time frames.To resume, balancing the trade-off between acquisition speed and temporal resolution of dynamic sequences is still a challenging task. In this work, we propose to use FFE sequences, allowing the assessment of real functionalities (of muscles and tendons for example) with more comfortable and rapid acquisitions. The idea is to fastly acquire a reduced amount of data and then to estimate the missing amount with the help of image post-processing techniques. To this end, we propose to recover the missing data by continuously interpolating the temporal deformation field of the joint, which can be estimated using diffeomorphic registrations from the acquired data assuming the motion to be continuous and consistent during scanning. In this context, we have extended the Logeuclidean polyrigid framework proposed by Arsigny et al. [START_REF] Arsigny | Polyrigid and polyaffine transformations : a novel geometrical tool to deal with non-rigid deformations-application to the registration of histological slices[END_REF] to deal with non-stationary dynamic MRI data by adding the time dimension into the stationary spatial framework. In the literature [START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF], [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF], the exponential mapping for solving the polyrigid Ordinary Differential Equation (ODE) is recursively computed using the scaling and squaring method involving approximation theory for which the accuracy always remains debatable. In the work of Moler et al. [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF], the authors have synthesized a study to compare the different existing techniques for calculating the matrix exponential (like methods involving approximation theory, differential equations, or matrix eigenvalues).Hence, we propose to use matrix eigendecomposition for exactly computing the exponential mapping since computations are performed in nondense regular grids due to the low-resolution of dynamic data, so that all the matrices can be stored in the main computer memory without the need of high-capacity RAM.

II. DYNAMIC MRI SCANNING

Five healthy children volunteers (between 10 and 14 years) participated in this pilot study. Each child was laid prone inside an Achieva dStream 3T MRI scanner (Philips Medical System, NL) with fully-extended knee, and two flex coils were placed on the ankle. The rate of ankle motion was standardized for all subjects using an MRI safe orthotic fixture (i.e. mean rotation speed of 4 • /s, each sequence captured the motion in 15 time frames). Real-time FFE sequences were used to acquire dynamic data during a single cycle of ankle plantal-dorsi flexion-extension (TR/TE=20.61/1.8, FOV=200x200, FA=15°, voxel=0.57x0.57x8mm, acquisition time = 18 sec). For each child, two sequences during passive joint motion were acquired for repeatability measures.

III. METHODS

A. Skeleton-tracking algorithm

In this section, we present a robust skeleton-tracking algorithm employing intensity-based registration. For a given joint trajectory, each data-sequence is first splitted to K 3D time frames (352 x 352 x 6 voxels each one). Bones of interest (calcaneus, talus, and tibia) are manually segmented in the first or last time-frame and then idividually and automatically tracked throughout the sequence, forward or backward, using rigid registration. Our tracking method within the dynamic sequence is composed of two distinct steps : 1) Estimating bone rigid transforms by locally maximizing the correlation ratio between neighboring time frames using voxel-wise weighting based on the bone segmentations.

2) Propagating the bone segmentations forward or backward by applying the estimated rigid transforms using nearest neighbor interpolation. Contrary to other interpolation methods like trilinear which induces smoothing of the bone mask, the use of nearest neighbor interpolation allows for preservation of bone shapes. To go forward in time, only one image interpolation is needed to propagate segmentations from the first time frame into each other. Since the composition of mappings in the Lie group SE(3) corresponds to the multiplication of elementary transformations (i.e., SE(3) × . . . × SE(3)), the direct path between images is computed by composing a set of estimated successive rigid transforms in order to avoid the propagation of interpolation errors as much as possible. The same for backward propagation for which the last time frame becomes the starting point for the tracking process. To improve the robustness of the algorithm, rotations are estimated based on their quaternion representations. Quaternions or versors are the most widely used in computer graphics and robotics, precisely in motion interpolation algorithms [START_REF] Mukundan | Quaternions : From classical mechanics to computer graphics, and beyond[END_REF] as they offer a compact representation of rotations.

B. Linear interpolation of bone transformation matrices

Let T i,k be the rigid transform from time frame D k to time frame D k+1 for the bone i. Assuming that this transformation is diagonalizable (which is the case for most transformations in an orthonormal basis, such as the image coordinate system), then there exist an orthogonal matrix P and a diagonal matrix D such that T i,k = P DP -1 . The nonzero elements of D are the eigenvalues of T i,k (λ 1 , λ 2 , λ 3 , and λ 4 = 1), while the columns of P are the corresponding eigenvectors. Based on this matrix eigendecomposition, one can define T i,k raised to its δk th power as follows :

T δk i,k = P D δk P -1 = P.diag(λ δk 1 , λ δk 2 , λ δk 3 , 1).P -1 (1)
Integrating δk continuously from 0 to 1 will infinitesimally change the matrix T δk i,k from the identity to the matrix T i,k . This allows for interpolating between two bone poses from a realistic transformation matrix, expressed in homogeneous coordinates (i.e., in the Special Euclidean group SE(3)).

C. Spatio-temporal Log-euclidean polyrigid framework

The Log-euclidean polyrigid framework provides an efficient way to synthesize joint deformation fields with nice properties like invertibility and smoothness, independently of the way local bone rigid transformations are first estimated. Based on the algebraic properties of this framework (the fact that it outputs diagonalziable velocity components in the Lie Algebra se(3)), it is then possible to add the time dimension to the stationary framework (since log(T δk ) = δk.log(T )). Given a point x in the k th time frame D k , the target location of this point in the interpolated time frame D k+δk (located between D k and D k+1 ) can be computed as follows :

Φ D k →D k+δk (x) = exp δk N i=1 wi D k (x) log(T i,k ) .x (2)
where Φ δk=0 (x) = x ; δk ∈ [0, 1] is the time-term added to the stationary polyrigid fusion formula ; T i,k is an element of the Lie group SE(3) ; N is the total number of rigid components ; Φ D k →D k+δk is the infinitesimal deformation field from D k to the time frame to be interpolated D k+δk ; wi In the literature, the exponential map of ( 2) is approximated using the scaling and squaring method [START_REF] Higham | The scaling and squaring method for the matrix exponential revisited[END_REF]. In this work, this exponential map is computed exactly using eigendecomposition. Similarly to the matrix power defined in (1), the matrix exponential is obtained by exponentiating the associated eigenvalues (i.e., the non-zero elements of the diagonal matrix D). Thus, we just needed to cast the exponentiation of the complex eigenvalues of the transformation matrices in regular grids leading to smooth dynamic deformation fields that preserve the topology of the bones. Finally, the floating image intensities are mapped to new coordinates in the target image space by spline interpolation.

D. Definition of weighting functions

In the context of articulated registration, the choice of weight functions is very important as the main goal is to accurately estimate the non-linear joint deformations without affecting bone shapes. It is important to note that the Gaussian weight functions are more suitable for polyaffine fusion than the polyrigid fusion, as they do not guarantee the rigidity and accuracy of the transformations in the bonny skeleton (by affecting sharp peaks particularly while smoothing). In [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF], Commowick et al. proposed an inverse-distance weighting function that preserves bone shapes after registration (3). However, such weighting functions yield inaccurate deformation outside the segmented bones.

w i D k (x) = 1 1 + αdist(x, B i D k ) β (3)
where B i D k is the binary mask of bone i in time frame D k , and dist(x, B i D k ) is the Euclidean distance between x and B i D k . Experiments show that an increase in β leads to an increase in the deformation accuracy for non-rigid structures. Thus, we propose new weighting functions suited for articulated registration. These new weight functions, which are inversely proportional to a rate of the distance exponential can be used within the log-euclidean framework in order to obtain more realistic deformations outside the segmented bones :

w i D k (x) = 2 1 + exp(γdist(x, B i D k )) (4)
Note that these new weight functions yield more accurate transformations especially in the case of non-large deformations, which is always the case between two successive frames. γ ∈ [0.4, 0.8] to ensure smooth interpolation.

IV. RESULTS

The proposed pipeline has been applied on five dynamic MRI sequences. The exponential map of equation ( 2) is computed within 3 s on an Intel ® Xeon ® Processor E3-1271 v3 3.60 GHz on a 352×352×6 regular grid using the proposed eigendecomposition method. Our code is implemented in Python T M using the LAPACK (Linear Algebra Package) routines for computing eigenvalues and eigenvectors. The RMSE between Achilles tendon contours was around 2mm when using the weigthing functions defined in in equation (3) with β < 3. Thus, the new definition of weighting functions played a crucial role in the enhancement of deformation accuracy outside the segmented bones while protecting skeleton edges, as reported in table I. Figure 1 illustrates, for one subject, such high-temporal-resolution reconstructed data. The robustness of the method and the accuracy of the results have been evaluated using a local leave-one-out crossvalidation technique. This was done by reconstructing each acquired time frame D k based on the deformation field halfway between D k-1 and D k+1 for k ∈ {2 . . . K -1}. Meaning each image D k has been reconstructed from the image-intensities of D k-1 , where T i,k is the rigid transformation from time frame D k-1 to time frame D k+1 for the bone i and δk is set to 0.5 in (2) (in the case of forward tracking). This gave a set of reconstructed time frames {D k } k∈2...K-1 . Results on bone motion interpolation are reported in table I. The accuracy of estimated deformation fields was then validated by computing the DICE overlap between the manually delineated masks of bones of interest on each acquired time frame D k and the bone segmentations, automatically propagated onto each reconstructed time frame D k . For all sequences, the DICE coefficient was greater than 0.82 indicating accurate propagations of bone masks over each dynamic sequence. In order to provide a meaningful information about the accuracy of estimated deformation fields outside the segmented bones, we have computed the RMSE (Root Mean Square Error) between contour-points of Achilles tendon on each acquired time frame D k and on each reconstructed time frameD k . In all cases, the mean error was less than 1mm (using the new weighting functions), indicating the accuracy of the estimated deformation fields for the non-rigid structures. Figure 2 illustrates the importance of the weighting functions when estimating and interpolating the temporal joint deformation field.

The proposed design for bone motion tracking is still robust in either direction (i.e., when going forward or backward in time). Note that this interpolation technique is still available in the case of polyaffine transformations, in where local transformations are affine (i.e. include additional degrees of freedom like the scaling and shearing parameters).

V. CONCLUSION

In this work, we have presented a motion-interpolationbased method for increasing the temporal resolution of anatomical dynamic MRI sequences. We have generalized the Log-Euclidean polyrigid registration framework to dynamic articulated structures and we have also proposed new weight functions which are well adapted to our context. The exponential map of the Log euclidean polyrigid framework is computed in an efficient and elegant way using matrix diagonalization-based techniques. To conclude, the proposed post processing technique aims to overcome the physical limitations related to real-time dynamic MR imaging algorithms which are generally based on compressed sensing theory [START_REF] Gamper | Compressed sensing in dynamic MRI[END_REF], for which it is hard to fastly acquire the entire or nearly the entire joint trajectory inside the MR scanner because of the limited k-space sampling.

  weighting function (i.e., N i=1 wi D k (x) = 1, ∀x ∈ D k ) which defines the local influence of the bone i on the final voxel transformation (see section III-D).

This work was supported by Région Bretagne, Chaire d'excellence INSERM-IMT Atlantique, Fondation de l'Avenir, and Fondation Motrice. Code: https://github.com/rousseau/dynMRI

 4) with γ = 0.4. The contours of bones (magenta) and of Achilles tendon (white) have been drawn in the target image to show the reconstruction accuracies when changing the weighting functions. Despite the fact that all the used weighting functions conserve bone topologies (this was validated by checking that the Jacobian map associated to the joint deformation field is equal to 1 over bone segmentations), the reconstructed image is more accurate for non-rigid structures such as the Achilles tendon especially when using our proposed weighting functions (d).


domain: sdv spi

text:
Guillaume Auzias

François Rousseau

Sylvain Takerkart

Nadine Girard

Christine Deruelle

Olivier A Coulon

Julien Lefèvre

Quantifying perinatal brain maturation using anatomical MRI

à 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.

I. INTRODUCTION

Continuous, sustained methodological advances in image acquisitions and processing (fast acquisition, movement correction, super-resolution/reconstruction) allow to address new problems that were not affordable few years ago (e.g. RS-fMRI, tractography) [START_REF] Dubois | The early development of brain white matter: a review of imaging studies in fetuses, newborns and infants[END_REF], [START_REF] Rousseau | In Vivo Human Fetal Brain Analysis Using MR Imaging[END_REF]. More specifically, dedicated image processing pipelines allow to extract a 3D mesh representation of the cortical surface from an anatomical MRI acquired either before or after birth [START_REF] Makropoulos | A review on automatic fetal and neonatal brain MRI segmentation[END_REF]. Global shape descriptors such as the gyrification index can then be computed for each individual and compared across groups. Early disturbance in cortical development at birth in preterm babies can be quantified using this approach, and the extracted measures are predictors of infants'outcome at term equivalent age [START_REF] Dubois | Primary cortical folding in the human newborn: an early marker of later functional development[END_REF]. Global gyrification measures were also used in [START_REF] Lefèvre | Are Developmental Trajectories of Cortical Folding Comparable Between Cross-sectional Datasets of Fetuses and Preterm Newborns?[END_REF] to compare intra uterine maturation in fetuses and extra uterine maturation in age-matched preterm infants and showed that the folding dynamics is different.

II. FINER SCALE DESCRIPTORS

In order to get access to more local shape descriptors, we recently adapted a surface registration technique developed initially for adult brain to address the problem of defining spatio-temporal correspondences across cortical surfaces from a group of fetuses [START_REF] Auzias | Model-driven parameterization of fetal cortical surfaces[END_REF]. The analysis of local curvature showed that brain folding is a progressive process, and that large folds observed in adults do emerge from smaller indivisible entities. Based on this observation, we developed a technique to decompose the sulci into smaller entities corresponding to putative first sulcal folds, denoted as sulcal pits [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. Using a closely related approach [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF], Im and colleagues have shown that sulcal pits are valid descriptors of brain folding in fetuses, allowing to detect abnormal folding patterns [START_REF] Im | Quantitative Folding Pattern Analysis of Early Primary Sulci in Human Fetuses with Brain Abnormalities[END_REF]. On the other hand we have shown that sulcal pits are also suitable for characterizing postnatal brain maturation and abnormal folding in a group of children with Autism (Corresponding author e-mail: guillaume.auzias@univ-amu.fr) spectrum disorder between 1.5 and 10 years old [START_REF] Brun | Localized Misfolding Within Broca's Area as a Distinctive Feature of Autistic Disorder[END_REF]. Finally, sulcal pits and corresponding shape descriptors can also be extracted from synthetic folded surfaces resulting from bio-mechanical modeling such as proposed in [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF] and compared to real data. In the future, these predictive models will be decisive for enhancing our understanding of neural mechanisms resulting in abnormal cortical maturation in developmental brain disorders.

III. CONCLUSION

Quantitative local/regional descriptors of the geometry of the cortical surface have a great potential to highlight abnormalities in cortical folding pattern and to detect departure from typical developmental trajectories in the perinatal period, both before and after birth. They will be instrumental to uncover the biological mechanisms underlying the cortical maturation.

  G.A., S.T., C.D., O.C. and J.L. are with Aix Marseille Univ, CNRS, INT, Institut de Neurosciences de la Timone, Marseille, France. F.R. is with IMT Atlantique, LaTIM U1101 INSERM, UBL, Brest, France. N.G. is with CRMBM UMR 7339, Aix-Marseille Université, CNRS, France and APHM Timone, Service de Neuroradiologie Diagnostique et Interventionnelle, Marseille, France.