misoSR: Medical Image Isotropic Super-Resolution Reconstruction

Isotropic volumetric acquisition in Magnetic Resonance Imaging (MRI) is often challenging. A large number of factors such as patient or physiological motion, signal-to-noise ratio (SNR), available static magnetic ﬁeld B 0 and total scanning time limit the acquired resolution. Super-resolution (SR) is a post-processing technique that optimally combines several anisotropic scans into a single isotropic volume that was not - or could not be - acquired in practice. If necessary conditions are met, the resulting isotropic volume oﬀers clear improvements over the initial acquisitions such as reduced partial volume eﬀect, oblique visualization and improved sharpness and SNR. This paper details the misoSR implementation of a SR isotropic reconstruction algorithm using the Insight Toolkit Library (ITK) library. It is a generic implementation that reconstructs an isotropic volume from any number of anisotropic volumes acquired from any orientation. The algorithm takes advantage of the inputs header information to handle the diﬀerent scans properties such as ﬁeld of view (FOV), resolution parameters and orientation. Step by step details on the implementation are given, parameters are individually detailed, and results are shown on diﬀerent applications as an example of SR reconstruction. The algorithm is hosted on the Creatis Virtual Imaging Platform (VIP), which allows users to run misoSR without having to install the software on their system.


Introduction
In magnetic resonance imaging (MRI), a large number of factors such as patient or physiological motion, signal-to-noise ratio (SNR), available static magnetic field B 0 and total scanning time limit the acquired resolution.In order to reduce the acquisition time and increase the SNR, highly elongated voxels are usually acquired by increasing the slice thickness.This inevitably limits the accurate assessment of anatomical structures in the through-plane dimension, increases the partial volume effect (less contrast and more blur), penalizes oblique reconstructions and further processing such as structure segmentation, 3D modelling or registration.Super-resolution (SR) is a post-processing technique that optimally combines several anisotropic acquisitions of the same scene into a single isotropic volume that was not -or could not be -acquired in practice.If each acquisition provides independent information, which is the case when patient motion occurs or when different scanning planes are used, resolution enhancement is possible.
Unlike classic interpolation techniques (linear, cubic, zero-padding-interpolation), SR is able to add significant high-frequency content into the reconstructed volume and efficiently remove some aliasing that might be present in the initial images if the acquisition sampling does not meet the Nyquist criterion.Note that two types of resolution enhancement can theoretically be performed using SR, namely i) in-plane and ii) through-plane resolution enhancement.The feasibility of in-plane improvement in MRI using SR remains a subject of controversy because of i) the inherent lack of high-frequency components acquired in the in-plane dimensions and ii) the fact that a shift in the in-plane dimension corresponds to a simple phase modulation in the acquired k-space [1].This means that for two spatially shifted volumes with similar field-of-views (FOV), similar k-space points are sampled which does not meet the necessary observation independence criterion for resolution enhancement.Later studies have been performed to clarify this point [2,3,4].On the other hand, through-plane resolution improvement, and in particular isotropic reconstruction, has been extensively and successfully applied to MRI in various contexts such as functional imaging [5,6], volumetric imaging [7,8], brain MRI [9,10,11].Review articles can be found in [12,13].Various acquisition strategies have been compared in [14] to optimize the reconstruction.Combining orthogonal acquisitions provides an optimal sampling of the k-space and offers a good trade-off between resolution, signalto-noise ratio, and acquisition time [15,8].This paper presents a generic implementation of a SR algorithm -misoSR -that reconstruct an isotropic volume from several observations.It can handle various acquisition protocols (number of inputs, volume orientation, different FOVs) since the implementation is based on the header information.The proposed algorithm is written in C++ and based on the Insight Toolkit (ITK) library [16].It is uploaded on the Virtual Imaging Platform (VIP) [17] hosted at Creatis (vip.creatis.insa-lyon.fr),which allows users to run misoSR without system compatibility constraints.This article starts by presenting the theory of SR, then details the implementation step by step and presents each individual parameter.Different working examples including the set of used parameters are given and limitations and perspectives are listed to conclude.

Image Acquisition Model
SR algorithms are based on an acquisition model given in Eq. 1 that links the observations Y k to the high-resolution image X. Observations are considered displaced, blurred, down-sampled and noisy versions of the ideal and unknown high-resolution image: In this equation, N is the total number of observations, G k is the geometric transformation operator that encompasses both the transformation due to the plane acquisition orientation and the motion between acquisitions.B k is a space invariant 3-D point-spread function (PSF), D k is the down-sampling operator from the high-resolution space to the observation space, and V k is an additive zero mean Gaussian noise.It is assumed that the PSF (B k ) is separable and has a Gaussian profile in each dimension.The choice of a Gaussian PSF is motivated by the study of Greenspan et al. in which it is shown that such a model is able to produce high-quality reconstructions [9].In the through-plane dimension, the variance (σ t ) of the PSF is usually chosen so that the full width at half maximum (FWHM) equals to the slice thickness (t), such as: σ t = t 2.35 .The PSF profile in this dimension is thus defined by: The in-plane PSF is constructed as a convolution of two identical Gaussian functions which variance is user defined.

Regularization
The general model given in Eq. 1 is reformulated as a regularized inverse problem because of its ill-posed nature.Let Ŷk = D k B k G k X be the observations estimator, and Γ T T V ( X) the regularization term: where λ balances the influence of the regularization term and the data fidelity term.The chosen regularization term prevents the apparition of noise by penalizing the L 1 norm of the magnitude of the image gradient, which favours the reconstruction of smooth images with sharp transitions [18].Farsiu et al. introduced in [19] the bilateral-TV that proved to produce high-quality results and to be computationally efficient, which is an interesting consideration when processing large dataset.In misoSR, the bilateral-TV regularizer is extended to 3-D to match the dimension of the isotropic reconstruction problem.It is defined as: In this equation, P is the size of the regularizer kernel, S lmn xyz is the global shift operator and regroups respectively the horizontal (S l x ), vertical (S m y ) and depth (S n z ) circular shift operators.For example, (S l x ) shifts X by l pixels in the x axis and allows the calculation of the gradient in this dimension.The scalar parameter (0 < α < 1) controls the influence of higher scales derivatives in the regularization term.In a 3D problem the choice of P has a significant impact on the computation time.If P is set to 1, the TTV operator can be viewed as the classical TV operator and α does not impact the regularization process.

Convergence
To iteratively converge towards the optimal estimation, a gradient descent algorithm is used as explicitly shown in Eq. 3. The step size γ is decreased at each iteration by the parameter β, (0 < β ≤ 1).Note that the step size is also normalized with respect to the number of inputs to facilitate convergence in case a large number of inputs are used.In Eq. 3, T represents the transpose operator, sgn the sign function and Id the identity matrix.A stable solution is considered to be reached when the difference between two successive iterations is below a pre-determined threshold which is detailed in section 3.3.

Implementation
In this section, the successive steps of the algorithm and their associated parameters are detailed.

Registration
Registration is a key step in the SR process and should be optimized with great care.It corresponds to the deformation operator G k of Eq. 1 that estimates the displacement between successive acquisitions.As mentioned earlier, the first input volume specified in the option file is used as the rigid image and all subsequent inputs are registered to it.The output image will thus be reconstructed in the same position as Y 1 .Note that the choice of the rigid input might have an significant impact in some applications.A multi-resolution scheme is used to improve speed, accuracy and robustness.For more details on the multi-resolution scheme, the reader is referred to [16].The number of multi-resolution levels can be set.Setting it to one is equivalent to a normal registration strategy.To avoid divergence, the maximum step length can be set, together with the maximum number of iterations per resolution level.Convergence is considered reached when the difference between successive iterations is below a minimum step length or when the maximum iteration number is reached.Iteration details are output to monitor the registration convergence with the selected parameters.Having a look at it is strongly advised at least for first uses since it is hard to set default registration parameters robust to any applications.The registration used in misoSR is a 6 degree-of-freedom rigid transformation, allowing 3D translations and 3D rotations since it is the most common transformation found in practice.This might be a limitation for some applications but deformable registration often requires fine tuning which we decided to let the user handle separately.

Image Modelling
This section explains how the different steps of imaging model of Eq. 1 are implemented, and details the parameters associated to each process.

Initialization
A first guess of the output image, X0 , is performed by averaging all the inputs after they have been registered to the rigid image and linearly interpolated to the isotropic sample grid.We have noted that the choice of the initialization does not dramatically influence the reconstructed image but a clever choice might reduce the required iteration number.

Deformation
This step uses the result of the registration step where the deformation operators were calculated.Note that in the registration step, the displacement from position k to position 1 is estimated, which actually corresponds to the operator G T (k): Y k → Y 1 for k = {2, . . ., N }.The current iteration estimation (in position 1), Xn , is successively displaced with the operator G k to match the k-th position.Note that G 1 is the identity operator.

PSF Filtering
The next step corresponds to the PSF filtering of each previously deformed image.The 3D Gaussian PSF function is defined in part 2.1.The PSF variance in all three dimensions can be set.It should be noted that when no in-plane resampling is performed, i.e. when all inputs have similar in-plane resolutions, a low variance (σ IP = 0.75 px) is usually a good choice.The variance of the through-plane PSF is set by default so that the full width at half maximum (FWHM) equals the slice thickness.However, user can

Down-sampling
The last operator of the imaging model consists of down-sampling the current estimation to each of the observation sampling space.This operation is performed using a simple linear interpolation.

Error Computation
Successively applying the previously described operators allows to create the estimated observations: Ŷk = D k B k G k X that lie in the same sampling space as the actual observations Y k .In order to improve the current estimation at next iteration, the error, e k = Y k − Ŷk , is computed and back-projected into the high-resolution space as can be seen in Eq. 3.

Error Back-projection
The error back-projection step consists of applying the transform operators to the estimation error.It reverses the image acquisition process.The transform of the downsampling operator is performed by upsampling e k back to the isotropic sampling grid.Zero values are inserted in between known values as illustrated in Figure 1.For noninteger sampling factor, voxel values are linearly interpolated to match the isotropic grid positions.The transform of the Gaussian blur operator is the blur operator itself.The same filtering that was applied for the error computation is applied to the upsampled error.Note that after upsampling, the number of zero-values present in the upsampled volume increases with the upsampling ratio.For this matter, the intensity values obtained after PSF filtering are normalized with respect to the through-plane to in-plane ratio.Finally the displacement operator is inverted to move each error image back to the reference position (X) by applying:

Stopping criterion
To stop the optimization process, a stopping criterion is set by the user (tolerance).
When the difference between two successive iterations is less than the fixed threshold, the optimization process stops.The stopping criterion is normalized with respect to the difference between the first and second iteration.So if set to 0.03, process stops when the difference between the current and previous iterations is less than 3% of the difference between the initialization and the second iteration.In practice, this has proved to be a reliable method for most applications.At the end of each iteration, the value of this criterion is displayed for the user to monitor the convergence process.In parallel, the maximum number of iterations can also be set to avoid excessive running time.

Examples
This sections presents two examples of SR reconstruction based on misoSR.Both examples with inputs and parameter values can be found on VIP.

Data
A first application of misoSR is shown in the context of T2-weighted rat brain images.Three orthogonal images are acquired in each anatomical plane on an horizontal MRI 4.7T Bruker Biospec with a surface quadrature receiving coil adapted to rat brain (Rapid Biomed).They are illustrated in Figure 2. The acquired voxel size is 0.1×0.1×0.5 mm 3 .A 2D multi-slice RARE sequence is used without slice gaps (TE=72.9ms, TR>4.8 s depending on the number of slices needed to cover the whole brain).The TR is considered long enough for the longitudinal magnetization to recover.Acquisitions are triggered on the respiratory motion, and the animal is anaesthetized with isoflurane during the imaging process.

Reconstruction
An isotropic reconstruction was performed at a voxel size of 0.1 × 0.1 × 0.1 mm 3 and is shown in Figure 2. The coronal image was used as the reference for the registration process.Note that a field inhomogeneity corrector is used as the only pre-processing step [22].The parameters used for the reconstructions are listed in Table 1.Total computation was around 7 minutes on a 3.7Ghz, 8 cores machine.13 iterations were needed for convergence.The full volume is 256×128×273 pixels and can be downloaded from VIP for accurate evaluation.
The coronal plane is taken as reference image for registration.All inputs are considered being at a similar phase of the breathing cycle.Orthogonal inputs and SR reconstruction are shown in Figure 3.The parameters used for registration and super-resolution are listed in Table 2. Total computation time was around 1 minute on 3.7 GHz, 8 cores machine to output a 128 × 97 × 120 pixels volume in 10 iterations.

Inputs -Outputs
Inputs are required to have correct header information.Standard image formats (Nifti, Dicom) are expected.Each input must be a single volumetric image.
The number of inputs used is up to the user.Using a single input is equivalent to single image interpolation with TV regularization and can be compared to classical interpolation methods (linear, cubic).A large input number will induce longer processing time but generally increases the reconstruction quality.Note that the first input is taken as the reference image for the subsequent registration steps.Depending on the application, this can significantly influence the reconstruction.Finally, note that the orientation of the output file is set to the following matrix : (1 0 0 ; 0 1 0 ; 0 0 1).This usually causes the orientation of the output volume to be different from the orientation of the inputs.

Parameters summary I/O:
Reference Image : The complete name (name + extension) of the image that will be used as a reference for the registration step.
Input Images : A single archive file (tar.gz)containing all volumes that will be used for the SR reconstruction.It must at least contain the file previously set as the reference image.

Registration:
Registration multi-res nb : Number of resolution levels used for registration Registration max step size : Maximum step size (in mm) allowed for registration Registration max iteration nb : Maximum iteration number allowed for registration at each resolution level SR reconstruction:

Conclusion
This article describes misoSR, a new tool for isotropic reconstruction of medical images.One or more volumes representing the same scene are combined to produce a single highresolution volume.Inputs are iteratively combined using an image acquisition model that comprises a rigid 6 degrees-of-freedom registration, a Gaussian PSF filtering and a down-sampling step.The overall reconstruction quality mainly depends on various factors : • the image acquisition model accuracy (registration accuracy, PSF choice).Default parameters are advised in the option file but should be optimized for each application.
• the resampling factor.The reconstruction quality usually decreases when the initial voxels anisotropy increases.A study performed in [24,8] illustrates this point and advises that anisotropy ratios higher than 5 should be avoided, i.e. slice thickness should not be more than 5 times than the in-plane pixel spacing.For high anisotropy ratios, it is usually hard to reduce the partial volume effect present in the initial scans.
• the initial scan orientations.Reconstruction quality increases when the redundancy between initial scans decreases.A denser coverage of the k-space logically leads to an optimal isotropic reconstruction, which explains why many studies combine orthogonal scans.
• the input homogeneity.It is important to keep in mind that the SR model assumes that all inputs provide the same information at similar spatial locations.Contrast in MRI can vary due to acquisition factors or field inhomogeneity.Important variations can thus alter the reconstruction procedure.Users might consider preprocessing such as field inhomogeneity corrections [22] before using misoSR.

Figure 1 :
Figure 1: Illustration of the upsampling process used for error back-projection

Table 1 :Figure 2 :
Figure 2: SR Reconstruction example on T2-weighted rat brain data.Boxed images show the original acquisitions ; the others are orthogonal reconstructions.

Axial ( 2 . 1 ×Figure 3 :
Figure 3: SR Reconstruction example on lung data.Boxed images show the original acquisitions ; the others are orthogonal reconstructions.

Table 2 :
Parameters used for the lung exampleThe second example of misoSR is on TWIST images of a lung cancer patient.3D TWIST images are obtained in three orthogonal planes on a whole body Siemens Trio 3T MR scanner (Siemens Trio MRI scanner, Siemens Healthcare, Erlangen, Germany).The inplane resolution of images is 2.1 × 2.1 mm 2 and the slice thickness is 7 mm with no slice gaps (TE= 0.8 ms, TR = 2.05 ms, FOV= 269 mm; Acquisition matrix = 128 x 128).