Non-rigid or deformable registration is an important tool for assessing spatial and temporal changes between images crum2004non; yu2019learning. It allows a non-uniform mapping between the source and target images. Among these non-rigid registration techniques are the diffeomorphic ones which have presented some interesting properties such as invertibility, differentiability, and smoothness of the resulting deformation fields. Furthermore, these techniques are suitable for parameterizing the non-linear deformations of anatomical structures or shapes which include the aspect of change over time arsigny2006log; durrleman2014morphometry. The Log Euclidean Polyrigid Framework (LEPF) is a parametric registration approach that composes a set of rigid transformations into a global diffeomorphism arsigny2006log; arsigny2009fast
. The main idea is to parameterize a non-linear geometrical deformation with a small number of flexible degrees of freedom. Hence, it can be used to estimate human joint deformations by fusing a set of local bone rigid transformations, composing the joint of interestmakki2019vivo. This approach, which is based on solving a first order autonomous ODE modeling the temporal evolution of a dynamical system, offers the possibility to estimate infinitesimal diffeomorphisms, parameterized by a dynamical time scale, via the exponential map from Lie Algebra
to the corresponding Lie group. This framework maps the interpolation of a deformation field into the tangent space, and then transforms the resulting velocity-vector field back onto the manifold. Which makes it possible to estimate smooth deformation fields in continuous time where theinvertibility is guaranteed. The LEPF can be used to estimate a smooth and continuous curves of moving articulated structures by interpolating the estimated polyrigid joint transformations in the domain of matrix logarithms. Which is in fact a generalization of the parameterization of rigid body transformations via the exponential map, as presented in vzefran1998interpolation; belta2002computation; beltanew; vemulapalli2014human.
Since the LEPF relies on the computationally heavy solution of an ODE, the efficient computation of the exponential map over a regular grid requires the implementation of a fast algorithm to deal with the high-number of point trajectories to be estimated simultaneously. In this context, Arsigny et al. have proposed to use a fast algorithm employing the scaling and squaring method arsigny2009fast. However, this method have a high memory requirement to store all matrices in the main computer memory during the repeated squarings.
Towards the end of the 1970s, Moler and Van Loan have synthetized a study to present the different ways of computing the exponential of a square matrix moler1978nineteen. This study has been revisited more than two decades later in moler2003nineteen
to take advantage of the advances of computational resources and to classify the existing methods and algorithms in terms of generality, accuracy, storage requirements, and efficiency. This study has pointed out that the algorithms which avoid use of the eigenvalues are more time-consuming for any particular problem.
Inspired by these conclusions, we propose to generate trajectories and motion interpolants from bone rigid transformations by scaling the transformation eigenvalues using a matrix-diagonalization based method.
The matrix diagonalization method has been already used to interpolate symmetric positive-definite matrices (i.e.tensor matrices) in arsigny2007geometric. However, this computational approach has never been used to interpolate homogeneous rigid transformations and the diagonalizability of such transformations has never been discussed in the context of the LEPF before.
In this work, we investigate the applicability of this technique for interpolating homogeneous bone rigid transformations, expressed in the image coordinate system.
The diagonalizability of bone transformations under joint mechanical contraints is justified in section 2.2.1. Then, we employ this technique to propose a fast algorithm for evaluating a Log-Euclidean polyrigid transformations on a regular grid. The main advantage of our algorithm is its capacity to drastically reduce memory requirements.
To demonstrate the effectiveness of the proposed algorithm on both synthetic data and real dynamic MRI data, we performed an objective comparison between the two methods in terms of accuracy, computation time, and memory requirements. And we concluded that the eigendecomposition method is more capable of balancing the trade-off between them.
Finally, we have applied the eigendecomposition technique to estimate a smooth continuous curve of the ankle joint trajectory from a discrete set of bone transformations. This allowed to interpolate intermediate time frames in-betweens acquired time frames and thus to enhance the temporal resolution of the dynamic MRI sequence.
2.1 Smooth interpolation on
2.1.1 Formulation of the interpolation problem
Consider two elements , . The interpolation function to estimate a smooth trajectory between them according to a time-varying parameter , is given by:
The group element that takes to can be obtained as follows:
The function transforms the interpolation into the tangent space for linearization purposes, and then transforms the resulting velocity vector back onto the manifold. A scaled map of to the Lie Algebra (where is the scaling factor) is given by:
Finally, the exponential map allows for transforming back into the manifold, yielding the following interpolation function:
2.1.2 An ODE approach for the interpolation on
Suppose now that we have the initial condition , then the displacement of the rigid body can be modeled by the following autonomous first order ODE beltanew:
To simplify notations, we denote . The solution of this ODE gives a smooth continuous curve of the rigid body trajectory:
The Equation (8) shows that the exponential map takes the linear tangential trajectory , into a one-parameter subgroup of , satisfying , for all . This gives a local parameterization for the geometric transformation .
2.1.3 Smooth interpolation of polyrigid transformations
To move from simply rigid to polyrigid interpolation of geometric transformations, Arsigny et al. proposed a registration framework to deal with autonomous continuous-time dynamical systems arsigny2006log
. The main idea is to parameterize a non-linear geometrical deformation with a small number of intuitive parameters by fusing a set of linear transformations into a global diffeomorphism, according to a predefined positive weighting functions.
Let be a set of transformations in , then the associated Log-Euclidean Fréchet mean, representing the Riemannian equivalent of the Euclidean arithmetic mean is obtained by minimizing the following metric dispersion:
where d(., .) is the distance associated to the Riemannian metric.
The resulting Log-Euclidean Fréchet mean is given by:
Finally, the polyrigid interpolation problem can be formulated as follows:
For a given point in the source image , the continuous-time trajectory of from to the target image can be computed according to:
where ; ; is the time parameter; , is the rigid transformation of the component from to ; is the total number of rigid components; is the infinitesimal diffeomorphism from source to target; is a normalized weight function (i.e., ) which defines the local influence of the component displacement on the resulting transformation at (see section 3.2.2).
2.1.4 Numerical computation of the exponential map
The matrix exponential of an matrix can be defined by the following convergent power series:
where is the identity matrix.
The computation of the exponential map for solving ODEs requires the development of a fast and efficient algorithm to deal with the repeated evaluation of the matrix exponential over a regular grid. In this context, Arsigny et al. have proposed a fast algorithm for the estimation of exponential map to parameterize polyrigid transformations using the scaling and squaring method arsigny2009fast. The scaling and squaring method al2009new; li2011matrix
is a recursive technique that exploits the fact that the matrix exponential can be easily estimated for matrices close to zero using the Padé approximants. This method is based on the relation. The scaling step consists of evaluating while the squaring step consists of squaring the approximant times to finally obtain an estimation of . The overall algorithm is composed of the three following steps:
Scaling step: so that , and .
Padé approximant to : , with: ;
; and , where represents the degree of the Padé approximant. The resulting error satisfies:
Repeated squaring: .
In arsigny2009fast, the ODE integration is performed using the scaling and squaring method, based on the following equations:
where ; is called the scaling factor, and .
The efficiency of this technique over a regular grid is somehow comparable to that of the Fast Fourier Transform (FFT) according to Arsigny et alarsigny2009fast. However, this method increases the computation time dramatically with an important increase of grid size because of the high memory requirements for storing all matrices in the main computer memory during the repeated squaring.
For example, if in (14), matrix multiplications are needed to compute the exponential map at each grid point, as follows:
until we obtain:
In general, the total number of matrix multiplications required for obtaining the complete trajectory of the articulated system is given by so that the number increases exponentially with the scaling factor. From another angle, the increase of the scaling factor has been recommended in higham2005scaling for enhancing the numerical accuracy of matrix exponentials and their inverses. A review the Higham’s method is proposed later in higham2009scaling for assessing the same accuracy with fewer matrix multiplications by increasing the degree of the Padé approximant. In arsigny2009fast, Arsigny et al. recommended the use of a scaling factor to obtain accurate and stable results.
2.2 Proposed approach
In this section, we propose to use a matrix diagonalizisation based method for fast computation of the exponential map with reduced memory consumption.
2.2.1 Diagonalizability of the homogeneous transformation matrices
In general, a rigid transformation matrix is diagonalizable. With the exception of screw transformations (i.e., when the motion consists of a rotation about a single axis with a non-zero displacement along the same axis). If we take the example of screw motion through the , the screw transformation will be:
where , is the translation vector, and is the rotation angle about the . This matrix is diagonalizable if and only if the translation along the , ghosal2006note. Otherwise, the eigenvalues of in the complex domain are:
. So that the eigenvectors corresponding to the repeated eigenvalueare linearly dependant, and thus the matrix is non-diagonalizable. For more details, the reader is referred to the work of Ghosal et al. ghosal2006note.
In the context of articulated motion, the human bones may never be able to perform screw motions in normal conditions due to the natural joint mechanical constraints. This suggests that rigid bone transformations are always diagonalizable.
2.2.2 Smooth interpolation of polyrigid transformations via matrix eigendecomposition
We propose to employ an eigenvalue-based method for smooth interpolation on , assuming that the map from the Lie Algebra to the corresponding Lie group is performed according to the equation (8).
Let be the rigid transform to be interpolated, then there exist an invertible matrix and a diagonal matrix such that . The diagonal elements are the eigenvalues of which satisfy the characteristic equation , while the columns of are the corresponding eigenvectors which satisfy the linear equation , also known as the eigenvalue problem. In the case of bone rigid transformations, the eigenvalue is very close to but still not exactly equal to . Based on this matrix eigendecomposition, one can compute the rigid body trajectory with respect to (8) as follows:
Thanks to the property , one can reformulate the interpolation function as follows:
Changing continuously from to will infinitesimally change the function from the identity to the matrix . This allows for interpolating between two rigid body poses from an homogeneous geometric transformation, by directly scaling its eigenvalues where is the scaling factor here. Figure 1 illustrates the linear interpolation of one simulated 3D rigid transformation using (17).
This technique can be also applied to interpolate any linear combination of bone rigid transformations and can thereby be embedded in the LEPF in this context.
3 Applications and results
3.1 Exponential map and parameterization of pure rotations: example of counterclockwise rotations
Let be the following matrix,
letting ; ; and , then the matrix is the homogeneous rigid transformation which generates a counterclockwise rotation by
radians about the origin of a 2D cartesian coordinate system. The logarithm of the matrix, which maps to its tangent space, is unique for and is equal to its principal matrix logarithm:
a one-parameter subgroup of the Lie group , which is parametrized with the parameter , is given by:
For and for different values of , Table 1 reports the error in the computation of the exponential map from the Lie algebra to the corresponding Lie group , using each of the two methods. This error is defined by: .
Table 1 shows that both methods are prone to neglectable rounding errors: occurring during the matrix inversion and multiplications for the eigendecomposition method; and during the repeated squarings for the scaling and squaring method. This confirms that both methods have nearly the same precision.
|Method||Error per rotation angle|
|2. Scaling and squaring|
In order to compare the performances of the two methods in terms of computation times, we have generated a set of regular grids with a size of , with . Each element of each grid is equal to , where is the homogeneous rigid transformation that encodes a counterclockwise rotation by radians around the . Our code is implemented in Python using the LAPACK (Linear Algebra Package) routines for computing the exponential map: we have used the routine tensorflow.linalg.expm which uses a combination of the scaling and squaring method and the Padé approximation (this routine optimizes the scaling factor automatically in order to increase the computational accuracy, depending on the matrix norm; for example for matrices very close to zero, the minimal possible scaling factor is used), and the routine numpy.linalg.eig to compute the matrix eigenvalues and eigenvectors. All experiments are performed on an Intel Xeon Processor E3-1271 v3 3.60 GHz, with a physical memory of 32GB.
Both methods are fast and numerically stable for sparse regular grids (i.e., for ). As illustrated in Figure 2, the eigendecomposition method is twice slower than the scaling and squaring method for sparse grids because of the additional time required for matrix inversion and multiplications, after solving the eigenvalue problem. However, the scaling and squaring method increased the computation time dramatically with an important increase of grid size (i.e., for ) because of the high memory requirements for storing all matrices in the main computer memory during the repeated squarings.
3.2 Application to the enhancement of temporal resolution of dynamic MRI sequences
Dynamic MRI is a non-invasive imaging technique that allows to qualitatively assess the behavior of the human joints in vivo. Real-time Fast Field Echo (FFE) sequences have the potential to reduce the effect of motion artifacts by acquiring the image data within a few milliseconds schaeffter2001projection. However, the short acquisition times affect the temporal resolution of the acquired sequences (i.e. the scanning duration is short relative to the joint motion). In this work, we propose to reconstruct the missing frames of the sequence given the reduced amount of acquired data, by applying the proposed motion interpolation technique which is based on scaling the transformation eigenvalues. This leads to estimate a continuous-time trajectory of the joint from a discrete set of rigid-body poses. To do this: we first estimate the rigid motion of each bone from the acquired discrete data using a robust intensity-based registration algorithm. Then, we fuse these local transformations into a global diffeomorphisms to compute the non-linear joint deformations between successive images according to (11). The idea is to reconstruct the missing time frames by interpolating the joint velocity vector field and then transforming the resulting infinitesimal velocity vector field back onto the manifold via the exponential map. The algorithm has been applied and validated using dynamic data from five subjects performing passive ankle dorsi-plantar flexion.
3.2.2 Estimation of temporal joint deformation field
The joint deformation field between each successive acquired time frames is estimated by fusing a set of rigid body transformations (corresponding to the bones of interest, where is the length of the acquired discrete sequence), according to Eq (11). A continuous trajectory of the joint is then obtained by varying the time parameter from to .
In commowick2006efficient; commowick2008efficient, Commowick et al. proposed an inverse-distance based weighting function that preserves bone shapes after registration (18). However, such weighting functions yield inaccurate deformation outside the segmented bones.
where , is the binary mask of the component, and is the corresponding Euclidean distance map at time (see Figure 3).
As illustrated in Figure 6, experiments show that an increase in in Eq (18) increases the accuracy of the estimated deformation field. Inspired by the relation (, when ), we have redefined the associated weight functions as follows makki2019temp:
To take into account the effects of all dependee components during fusion, each component weight function is normalized as follows:
Figure 4 illustrates each bone normalized weight function as defined in (19).
Finally, the source image intensities are mapped to new coordinates in the target image space by spline interpolation. The robustness of the method and the accuracy of the results have been evaluated using a local leave-one-out cross-validation technique. All results on joint motion estimation and interpolation are validated and reported in makki2019temp.
The exponential map of Eq (11) is computed within seconds on a grid using the scaling and squaring method (which required 1.8 GB of RAM), and within nearly seconds using the proposed eigendecomposition method (which required 0.8 GB of RAM). Figure 5 illustrates the interpolation quality of one 3D+time MRI sequence from our dataset.
If is the total number of intermediate points chosen to discretize the continuous trajectory of each point , the scaling and squaring method is the most commonly used technique for computing the matrix exponential in previous works arsigny2009fast; commowick2008efficient; ehrhardt2019temporal; porras2018locally. However, the use of this method results in excessive memory requirements as it requires multiple matrix multiplications (see Figure 2) and thus it is not suitable for estimating or interpolating very dense deformation fields (e.g. when registering very high-resolution MRI data makki2018high). To overcome this limitation, we have proposed to compute the matrix exponential using matrix eigendecomposition. The grid is interpreted as a tensor of matrices, each of size (with , , and are the image dimensions), and the trajectories of all the points of the regular grid are computed simultaneously. Furthermore, computations are performed in the complex domain: for all real transformation matrices expressed in homogeneous coordinates, the associated eigenvalues occur in complex-conjugate pairs where the imaginary part of the eigenvalue is the frequency of voxel rotation.
Both methods lead to nearly the same level of accuracy and they can also reach the machine precision since the error is always of order .
For sparse regular grids, the eigendecomposition method is a little slower than the scaling and squaring method, but it offers much faster computations in the case of dense regular grids. Hence, this method is more capable of balancing the trade-off between computation times and memory requirements. For example, the optimization of large storage requirements could be very useful when estimating a dense deformation field from a high-resolution static MRI scan to a low-resolution dynamic time frame makki2018high; makki2019vivo. In these previous works, a deformation vector field of millions deformation vectors is computed in 15 min using the eigendecompostion method while consuming the entire physical memory (32GB). However, the use of the scaling and squaring method for estimating the same deformation vector field results severe memory swapping which drastically slows down computations until totally stopping the computation process.
The eigendecomposition method allows for interpolating the desirable number of intermediate points and thus for obtaining a smooth and continuous trajectory of the dynamical system by scaling the transformation eigenvalues, without the need to perform a large number of matrix multiplications. In fact, once the eigenvalue problem is solved over the regular grid, the transformation eigenvalues and the sub-matrices containing the associated eigenvectors (with their inverses) can be stored in the main computer memory, and the interpolation of an intermediate point requires only two matrix multiplications.
The proposed interpolation technique can be also extended to diagonalizable polyaffine transformations, for which local transformations to be fused offers additional degrees of freedom like the scalings and shearings . Figure 1 illustrates the consistent interpolation of a simulated affine transformation using matrix eigendecomposition.
In most applications such as human skeleton tracking during motion, the amount of local rotation of each rigid component is strictly below radians, so that the matrix logarithm of a bone transformation exists always and is equal to the principal matrix logarithm. This suggests that the property , is verified for these transformations. Concerning the diagonalizability of homogeneous rigid transformations, we argue that the joint bones may never be able to perform screw motions in normal conditions due to the effect of mechanical constraints on the joint motion. Hence, the estimated bone transformations expressed in homogeneous coordinates possess 4 distinct eigenvalues and thus are always diagonalizable.
In this work, a motion-interpolation-based method for increasing the temporal resolution of anatomical dynamic MRI sequences is presented. We have generalized the LEPF to dynamic articulated structures and we have also proposed new weight functions to increase the local effect of bone transformations on the deformations of the surrounding tissues. The proposed inverse-distance-based weight functions are more suitable than the Gaussian weight functions seiler2012capturing for articulated registration. Since the Gaussian weight functions induces some smoothing which affects the bone shapes by affecting sharp peaks particularly. Both our new weight functions and the weight functions defined in commowick2006efficient; commowick2008efficient conserve the bone topologies. However, our weight functions yield more accurate deformation vectors outside the skeleton as illustrated in Figure 6 for the Achilles tendon deformations.
The proposed post processing pipeline aims to overcome the physical limitations related to real-time dynamic MR imaging algorithms which are generally based on compressed sensing theory lustig2008compressed, 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.
In this paper, we proposed a fast and efficient algorithm to compute the exponential map for solving ODEs modelizing the evolution of dynamical systems over a regular grid. An objective comparison between the eigendecomposition method and the scaling and squaring technique used in previous studies is performed using simulated data and four dimensional (3D + time) dynamic MRI data. Experiments show that the scaling and squaring method is more suitable for sparse grids since computational burden and memory requirements grow exponentially with grid size. The obtained results demonstrate also the potential of the eigendecomposition method to optimally balance the trade-off between accuracies, computation times, and memory requirements for the computation of the exponential map.
Acknowledgements.This work was supported by the French region of Brittany and the chaire d’excellence INSERM, IMT Atlantique.
The source code is available on github: https://github.com/rousseau/dynMRI/blob/master/Exponential_map.py