Applying positivity constraints to q-space trajectory imaging: The QTI+ implementation

Diffusion MRI is a powerful technique sensitive to the microstructure of heterogeneous media. By relating the dMRI signal obtained via general gradient waveforms to the moments of an underlying diffusion tensor distribution, q-space trajectory imaging (QTI) provides several quantities indicative of the structural composition of the medium. Substantial improvements in the reliability of the produced estimates has been achieved via incorporating necessary positivity constraints in the estimation by employing Semidefinite Programming. Here we present the Matlab code implementing said constraints, provide a simple example showing the main functionalities of the package


Introduction
Diffusion magnetic resonance imaging (dMRI) is a powerful imaging modality with exquisite sensitivity to tissue microstructure.By combining diffusion sensitizing MR sequences and models that capture the relevant descriptors of complex environments, it is possible to obtain valuable information about the scanned specimen with potential diagnostic utility.By utilizing general gradient waveforms [1] and a diffusion tensor distribution (DTD) model for the tissue microstructure [2,3], q-space trajectory imaging (QTI) [4] relates the signal to the statistical moments of the quantities describing the local characteristics [5].In QTI, the first statistical moments are the mean and covariance tensors of the underlying DTD, which can be obtained from the data and used to compute meaningful biomarkers.Several recent studies showed the sensitivity of metrics obtainable through the QTI analysis to various diseases [4,[6][7][8][9], thus prompting interest towards this technique.
A key component of such endeavor is to ensure proper estimation of the model's parameters.Noise and artifacts often affect the measured quantities of interest, thus limiting the reliability of such measurements or of parameters derived from them.To alleviate such confounding effects, it is desirable to constrain the fitting so that the estimated parameters fulfill conditions dictated by the underlying physics of the problem.In the here considered QTI framework, this translates into certain positivity conditions to be applied during the estimation of the mean and covariance tensors.Following the approach taken by Dela Haije et al. [10] for other dMRI models, we recently studied the problem of enforcing relevant constraints in QTI using Semidefinite Programming (SDP) and proposed to incorporate such constraints in a framework called QTI+ [11].The results showed that by forcing the two estimated tensors to respect mathematical and physical conditions, it is possible to obtain more stable and reliable quantities compared to the original QTI implementation [12] (available at https://github.com/markus-nilsson/md-dmri).Here we focus on providing an overview of the blocks building the Matlab (The Mathworks Inc, Natick, Massachussets) QTI+ software and its functionalities.
An in-depth comparison between this software and the original implementation is outside the scope of this work, but is available at [11].The results obtained in [11] can be reproduced using this software (available at https://github.com/DenebBoito/qtiplus), the unconstrained QTI implementation (available at https://github.com/markus-nilsson/md-dmri),and the publicly available data [13] used in that study.Using only this software and the data available at https://github.com/filip-szczepankiewicz/Szczepankiewicz_DIB_2019/tree/master/DATA/brain/MD-dMR I , it is possible to reproduce part of the results presented in [11] by following the instructions provided in the software's User Guide and collected in a script available in the package's Github repository.

Problem formulation
We keep the notation consistent with that of Herberthson et al. [11].Briefly, boldface characters are reserved for matrices and second order tensors, while italic characters denote scalars.As an example, A ij is a tensor, while A ij is the ij th element of the tensor A ij .Double struck characters indicate 4th order objects, and the ones considered in this can also be represented using 6 × 6 matrices.To distinguish between the two representations, i, j, k, and ℓ indexes ranging from 1 to 3 are used for 4th order objects, while α and β indexes ranging from 1 to 6 are used for the 6×6 matrix representation.As an example, C ijkℓ and C αβ are the same tensor in the 4th order and 6×6 matrix representations, respectively.The Einstein summation convention is used to represent products between tensors.For example, the tensor product, indicated with '':'', The QTI framework implies that, for a number of diffusion measurements whose experimental parameters can be condensed in the measurement tensor B ij , the MR signal is given by [4]: where S 0 is the signal without diffusion weighting, Dij is the mean diffusion tensor, and C ijkl is the covariance tensor.The problem then consists of the estimation of S 0 , Dij , and C ijkl from a series of signals Given the symmetries of Dij and C ijkl , there are 28 parameters that need to be estimated in each voxel: one for the value of S 0 , 6 for the unique components of Dij , and 21 for the unique components of C ijkl .
The problem formulation can be written as argmin which can be linearized and written as a weighted linear least squares problem (accounting for heteroscedasticity introduced by the logarithm [14,15]) In QTI+, the goal is to compute (2) or (3) with nonlinear routines or SDP while enforcing appropriate positivity conditions on the tensors Dij , C ijkl , and M ijkl , where

⟨
• ⟩ indicates the mean (expectation) value.The set of constraints is summarized as follows, with the conditions taking the name of the tensor they are enforced upon: For details, we refer to [11].

Software framework and architecture
The flowchart in Fig. 1 shows the QTI+ framework and the adopted nomenclature for the various parts.Each step is presented briefly in the following subsections.We note here that it is not necessary to follow all the steps in the framework.The first step is able to provide results which are, in our experience, only mildly improved upon in the subsequent steps in most voxels.Thus, as it can become time-consuming to perform all the steps, by default the software performs only SDPdc and computes the QTI accessible invariants on the estimates obtained this way.However, we provide other pipelines as options with which it is possible to pass the input through either all or parts of the framework.
We also shall remark that there are more constraints which could be enforced.Additional conditions could be easily incorporated into the software by, for example, adding blocks to its structure.

CVX
Before introducing the various steps, we briefly mention the library used to solve all the Semidefinite Programming problem.CVX is a Matlab package for specifying and solving convex problems [16,17].It is freely available for academic research and provides an interface to various SDP solvers, some of which are free to use while others require a license.By default in QTI+, and to obtain the results in [11], we used Mosek (version 9.1.9)(MOSEK ApS, Denmark).This software requires a license to be used, which is free for academia.The QTI+ software checks if Mosek is available before starting the computation.In case it is not, it uses SDPT3 [18] which is shipped and installed with CVX.Choosing one or the other solver should not affect the final result.From our experience, we found relevant differences only in the computational times, Mosek being the faster alternative.In the user guide we provide a comparison between the computational time achievable with both solvers on the different steps of the QTI+ framework.

Data and B-tensors
The required input to the framework are the diffusion MRI data in arbitrary units, and the measurement tensors in SI units (s/m 2 ).The data are expected to be a 4D Matlab matrix with shape As the first step, the problem described in ( 3) is solved with a Semidefinite Programming formulation while enforcing conditions (d) and (c), see [11] for details.SDP naturally incorporates positive semidefinite matrices in the problem as constraints.
Therefore, it is sufficient to pass Dij and C αβ to the solver to ensure the positivity of these estimates.

NLLSdc
This can be viewed as the second step of the algorithm, where we seek to improve the residuals by solving the original nonlinear problem in (2).The solution found with SDPdc is used as initial guess for the non-linear solver.Here the software uses the lsqcurvefit function in Matlab, with Levenberg-Marquardt as algorithm.To enforce conditions (d) and (c), the two tensors are expressed in their Cholesky decomposition throughout the fitting.In case this step does not provide an improvement to the residual, the initial guess is retained as solution.

m-check and SDPdcm
These are the last two steps where the algorithm checks whether condition (m) is satisfied, and if necessary imposes it.SDP is employed in both steps.The detailed description of how this is achieved can be found in [11].Note that here we need to fix the estimates of Dij to be able to impose condition (m).

Output
The output consists of the model parameters and the QTI derived scalars.The model parameters are returned in a matrix with shape [nx, ny, nz, 28], where the fourth dimension contains the estimates of S 0 and of the unique elements of Dij and C ijkl .
In addition, a number of scalar invariants derived from the estimated quantities is returned in a structure where each field corresponds to the maps listed in [4].

Software call and functionalities
The software is intended to be interacted with using the function qtiplus_fit.This is the main function controlling the different settings, and eventually calling the subroutines that perform the fit to the data.Through this function, the user can control the behavior of the software via key-value pairs.It is for example possible to select which of the framework steps to perform using the keyword pipeline, and which solver to use via the keyword solver.These and more options are discussed in-depth in the documentation accompanying the software.An example and a simple synthetic dataset are provided to get the user acquainted with the software and the different functionalities.These are discussed in the next section.
The most basic call to this function is made by only inputting the data and the measurements tensors as follows:

≫ [model, invariants] = qtiplus_fit(data, btensors)
When called as such, the software uses all the default options to compute the model parameters and the scalar invariants.briefly, • SDPdc is the only step performed • A mask is created to avoid computations on voxels located outside the area of interest.The algorithm used for this purpose is a modified implementation of the median_otsu segmentation method from Dipy [19] • Computations are performed in parallel exploiting multicore processing, if available • 50 voxels are sent at a time to the SDP solver • All inserted diffusion volumes are used to get the estimates We would like to note that, while the software is intended to be easily used through interacting with the function qtiplus_fit, the routines implemented for SDPdc, NLLSdc, and SDPdcm could be used as stand-alone as well as integrated in other processing pipelines.

Illustrative examples
In the software package we provide two examples illustrating the main functionalities of the software.The first example uses synthetic data available within the software.The second example shows how to reproduce some of the results in [11] on experimental data presented in [13], which can be downloaded from https://github.com/filip-szczepankiewicz/Szczepankiewicz_DIB_2019/tree/master/DATA/brain/MD-dMRI.This last example can be found in the software package in a script called exam-ple_qtiplus_fit_experimental_data.Both examples are described in detail in the User Guide.
Here we briefly describe the example using synthetic data which follows the simulations performed in [11].This example is available in the software's package in the script exam-ple_qtiplus_fit_synthetic_data.
For a family of 217 B-tensors employed in [13], the diffusion MR signal is generated using a non-central Wishart distribution [20] with known mean and covariance.As in [11], two signals are generated, one for a DTD with isotropic Dij , and one for a DTD with anisotropic Dij .Rician distributed signal is then obtained  Fig. 2 showcases how options can be passed to the fitting functions by means of key-value pairs.For this example, we decide to limit the computations to a central region of the dataset delimited by a [11 × 11] voxels window.Moreover, we exclude some of the measurements from the computations by passing a vector of indices indicating which volumes to use for the fitting.The other optional parameters are left to their default settings.The qtiplus_fit function is then called and the model parameters and the computed invariants are returned.The obtained results are displayed in Fig. 3 where the returned estimated values of one of the invariants, namely µFA, are shown.

Impact
Enforcing positivity constraints in different diffusion MRI models and in QTI in particular, was recently investigated [10,11].The results in both articles clearly showed a marked improvement as conditions are imposed during the estimation.In the case of QTI, it was also shown that applying positivity constraints not only increases resilience to noise, but also helps in retaining a higher degree of information when data are scarcely sampled [11,21].Therefore, an added benefit of employing QTI+ concerns the possibility of reducing the time spent acquiring data, a true limitation in clinical settings.Moreover, QTI+ is found to be robust to acquisition imperfections, which makes it suitable for frequently encountered situations such as incomplete data sampling or certain images having severe artifacts.Overall, the QTI+ approach is a step forward in providing more reliable and accurate parameters to be used in clinical studies.

Conclusions
In this paper we described the software implementing QTI+, a Matlab package for imposing positivity constraints in q-space trajectory imaging.The software design allows the user to easily control which of the available constraints to impose, and opens for possible implementation of additional constraints in future releases.The software was recently employed in a proof-of-concept study, the results of which demonstrated increased accuracy in the estimated metrics when compared to the original implementation [11].Parts of these results can be directly reproduced using

Fig. 1 .
Fig. 1.Flowchart of the QTI+ framework and nomenclature of the various components.

[
nx, ny, nz, nd], where n denotes the number of voxels, x, y, z indicate the respective spatial dimensions, and d indicates the ''diffusion" dimension.As an example, nx is the number of voxels along the x dimension, and nd is the number of diffusion measurements.The family of B-tensors which generated the signals are input as a [nd, 6] matrix, where each row contains one of the tensors in Voigt format.The following Voigt convention is used to convert a tensor from a [3 × 3] matrix representation to a [1 × 6] vector:

Fig. 2 .
Fig. 2. Snippet showing how to pass options to the software using key-value pairs.Additional information can be found in the examples provided with the software package, and in the User Guide.
by adding Gaussian noise to the real and imaginary parts of the noiseless signal, creating 4096 signal values at SNR= 1/σ = 30 and SNR= 15 for both cases, where σ is the standard deviation of the underlying Gaussian distribution.The data are then organized in a [64, 64, 6, 217] matrix, where each of the 6 slices [:,:,k,:], k = 1, . . .,6, represents the following: 1. Noiseless signals for a non-central Wishart distribution with isotropic Dij 2. Signals from 1. with SNR = 30 3. Signals from 1. with SNR = 15 4. Noiseless signals for a non-central Wishart distribution with anisotropic Dij

Fig. 3 .
Fig. 3. Results obtained on the example dataset at two levels of noise.Displayed are the µFA maps derived from the model parameters and accessible in the invariants structure.The first column indicates the ideal value obtained by fitting the QTI model to the noiseless signals, here considered as ground truth (GT), while the two middle columns show the results obtained on the noisy signals.Results are displayed in the [0, 1] range.The last column shows the value distribution around the ground truth.5. Signals from 4. with SNR = 30 6. Signals from 4. with SNR = 15 The interested user can find in the software package a function called create_new_example_dataset with which the same dataset with possibly different levels of noise can be created.This function allows the user to select what SNR the data should have in slice pairs 2,5 and 3,6.An example call to this function is provided in the example_qtiplus_fit_synthetic_data script.Fig.2showcaseshow options can be passed to the fitting functions by means of key-value pairs.For this example, we decide to limit the computations to a central region of the dataset