---
# System prepended metadata

title: Implementación de MRDS
tags: ["MRtrix\_", MRDS, dMRI]

---

Paulina J. Villaseñor
Laboratorio de Conectividad Cerebral C13
LANIREM

## MULTI-RESOLUTION DISCRETE-SEARCH

En este manual repasaremos paso por paso como implementar el método multi-tensor MRDS desarrollado por el Dr. Ricardo Coronado-Leija, para literatura mas extensa de como funciona el método, porfavor consulta el artículo [aquí](https://www.sciencedirect.com/science/article/abs/pii/S1361841517300956).

---

MRDS sigue una serie de pasos para poder ser implementado, primero antes que nada vamos a cargar los módulos necesarios:

```python=
module load mrtrix/3.0.4
module load mrds/1.0
```
MRDS lo puedes ejecutar en toooodo el cerebro o en una rebanada o en solo la corteza, eso dependerá de tu análisis. Para fines didácticos, vamos a hacer todo el cerebro de rata de una sola rebanada.

Antes de empezar, necesitamos preparar algunos datos que MRDS nos pedirá para su ejecución:

#### 1. Scheme
El scheme no es nada mas que la información de los gradientes de difusión de la imágen, es decir, tus archivos `.bvec` y `.bval`. Sin embargo, MRDS pide que estos esten juntos en un archivo de texto organizados de manera específica con 4 columnas en el orden *X Y Z b*, donde *x y x* es la dirección de los gradientes y *b* es el b-value asociado. 

Por ejemplo, si vemos como esta ordenado tu archivo `.bvec` original, vemos que esta de esta manera:

![image](https://hackmd.io/_uploads/H1OlbbogA.png)


Lo mismo sucede para el archivo `.bval`. Lo que necesitamos hacer es juntar ambos archivos y acomodar los datos en formato columna/tabla. Para esto, podemos usar un script que forma parte de los `inb_tools` que hace el trabajo!

```python=
# Unimos los archivos .bvec y .bval:

cat 64A_dwi.bvec 64A_dwi.bval > bvec_bval

# Usamos el script para modificar la estructura

transpose_table.sh bvec_bval > scheme
```
Teniendo como output lo siguiente:

![image](https://hackmd.io/_uploads/rJHMQbogR.png)

Vemos que las columnas siguen sin alinearse adecuadamente y esto es debido a que tenemos numeros muy largos despues del punto. Esto lo corregimos con el siguiente código:

```python=
cut -d' ' -f1-4 scheme | while read a b c d; do printf "%.5f %.5f %.5f %.4f\n" "$a" "$b" "$c" "$d"; done > scheme_nuevo
```
> Aprende [aquí](https://linuxize.com/post/linux-cut-command/) mas acerca de como usar el comando `cut`.


Ahora vemos que ahora las columnas estan bien alineadas y en orden:

![image](https://hackmd.io/_uploads/SyuFFWilA.png)

Como ultimo paso para el scheme, hay que remover los primeros bvals (b=21.01.08) y convertirlos en ceros, ya que MRDS no puede ser ejecutado sin volúmenes b=0, en este caso puedes tomar los mas bajos y convertilos. Aquí un código para poderlo hacer:


```python=
shells=`mrinfo -bvalue_scaling false -grad scheme_nuevo 64A_dwi.nii -shell_bvalues`
firstbval=`echo $shells | awk '{print $1}'`
if (( $(echo "$firstbval > 0 " | bc -l)  ))
then
sed -i -e "s/${firstbval}/0.0000/g" $scheme
fi
```
![image](https://hackmd.io/_uploads/BJH9he2xR.png)



#### 2. Máscara
Primero necesitamos una máscara del cerebro o del área de interes donde quieres que se implemente MRDS. Sin embargo, dado que MRDS se ejecuta utilizando una *response function*, asegurá que tu máscara incluya una estructura de alta anisotropía, como por ejemplo el cuerpo calloso. 
> La *response funcion* es definida como la señal esperada/ideal en un voxel que contiene una organización altamente coherente y perfecta de un solo manojo de fibras. Un gran ejemplo es el cuerpo calloso en su porción medial. 

Aquí hay de dos sopas, o creas una máscara general del cerebro utilizando `dwi2mask` o creas tu propio ROI/máscara binaria de la estructura y rebanada que quieras. Esta ultima opción requere que el ROI se dibuje manualmente (al menos que tengas un pipeline que lo haga por ti). Con fines didácticos y de tiempo, voy a crear un ROI utilizando `mrview`:

```python=
mrview 64A_dwi.nii.gz 
## Aquí haces tu ROI y lo guardas
```
Alternativamente puedes simplemente utilizar: 
```python=
dwi2mask -grad scheme_nuevo 64A_dwi.nii.gz mascara_64A.nii.gz
```
![image](https://hackmd.io/_uploads/S18fmM2xC.png)


### Ajuste de MRDS primera parte:

Una vez que tenemos estos preparativos, estamos listos para correr la primera parte de MRDS que es ajustando el DTI con un método no-linear. También puedes ajustarlo de manera más tradicional, etre otros, para esto hay que revisar la documentación y ajustar `dti` a tus necesidades y objetivos específicos:

```
dti

dti [ options ] dwi_input scheme tensor_output

        dwi_input
                name of the input dwi file (.nii).
        scheme
                name of the scheme (.txt) corresponding to the dwi image. Each line of the file must be:
                         x1 y1 z1 b1
                         x2 y2 z2 b2
                         x3 y3 z3 b3
                         .  .  .  .
                         .  .  .  .
                         .  .  .  .
                         xn yn zn bn
        tensor_output
                name of the output diffusion tensor (.nii).

Compute Diffusion Tensor from Diffusion Weighted Magnetic Resonance Images.

Options:

         -save option
                Which files to save: (0) tensors; (1) diff parameters; (2) both. Default: 2.
         -mask file(.nii)
                mask for selecting the voxels that will be processed.
         -nonorm
                Work with original signals Si (not recommended). By default the algorithm works with A_i = S_i/S_0.
         -beta num
                For the spatial bilateral estimation of S0^{s}_{snr}, S0^{s}_{std} and S0^{s}_{mean}.
                The estimation includes close voxels inside the range S0^{v}_{mean} +/- beta*S0^{v}_{std}.
                (s) -> spatial, (v) -> voxel.
                Default: beta = 0.5.
         -method name
                select the method used for estimate the diffusion tensor coefficients.
                the choices are:
                        -method naive      -> Log Linear Least Squares. Very simple implementation.
                        -method linear     -> Log Linear Least Squares. GSL Implementation.
                        -method nonlinear  -> Non Linear Least Squares. GSL Implementation.
                        Default: -nonlinear.
         -response nvoxels
                Estimate and compute the mean of the diffusion tensor eigenvalues using the voxels
                in the volume where there is high probability only one fiber bundle exist.
                This could be used for starting point on other multi-fiber algorithms.
                The algorithm will use the nvoxels with maximum FA in a confidence interval.
                Setting nvoxels = 0 will use all voxels in the confidence interval.
         -adc
                compute the apparent diffusion coefficient on each orientation of the scheme.
         -lps
                compute the linear, planar and spherical coeficients of the diffusion tensor.
         -fa
                compute the fractional anisotropy of the diffusion tensor.
         -md
                compute the mean diffusivity of the diffusion tensor.
         -dec
                compute the direction encoding color of the diffusion tensor DEC = FA*PDD.
         -mse
                compute the mse of the signals for the estimated diffusion tensors and the measured signal.
         -correction num
                Apply a transformation on the diffusion data in order to correct the Rician bias.
                This is not a denoising, it is just a bias correction. The choices are:
                        -correction  0, 1, 2 -> No correction: y = x.
                        -correction  3, 4, 5 -> Gudbjartsson correction: y = sqrt(fabs(x*x - sigma*sigma)).
                        -correction  6, 7, 8 -> Gudbjartsson correction modified: y = sqrt(max(x*x - sigma*sigma,0.0)).
                        -correction  9,10,11 -> Manjon correction: y = sqrt(max(x*x - 2.0*sigma*sigma,0.0)).
                        -correction 12,13,14 -> Maximum Likelihood correction y_ML = max_y(p(x;y,sigma)).
                                                Exhaustive evaluation.
                        -correction 15,16,17 -> Mean Posterior correction y = sum x*p(x).
                                                Exhaustive evaluation.
                        -correction 18,19,20 -> Mean Posterior correction adaptive according to each measurement.
                        Aditional consideration.
                        num % 3 => 0 - min value = 0.
                        num % 3 => 1 - min value = estimated from data.
                        num % 3 => 2 - min value = physically plausible exp(-b*3e-3).
                        Default: -correction 8.
         -help
                show this help

```
Entonces, nosotros podemos ajustar `dti` de la siguiente manera:

```python=
dti 
-mask ROI.nii \
-response 0 \
-correction 0 \
-fa -md \
-64A_dwi.nii \
-scheme_nuevo.txt \
-dwi.nii
```
> Ojo, si tus archivos los tienes en `.nii.gz`, utiliza `mrconvert` para cambiar a `.nii`.

Checamos nuestros outputs:

```python=
ls dwi_DTInolin_*

dwi_DTInolin_COMP_SIZE.nii    
dwi_DTInolin_ISOTROPIC.nii  
dwi_DTInolin_PDDs_CARTESIAN.nii           
dwi_DTInolin_ResponseIsotropicMask.nii
dwi_DTInolin_EIGENVALUES.nii  
dwi_DTInolin_MD.nii        
dwi_DTInolin_ResponseAnisotropicMask.nii  
dwi_DTInolin_ResponseIsotropic.txt
dwi_DTInolin_FA.nii           
dwi_DTInolin_NUM_COMP.nii   
dwi_DTInolin_ResponseAnisotropic.txt      
dwi_DTInolin_Tensor.nii
```
Vemos que tenemos de regreso una serie de archivos que nos serviran para ejecutar el segundo segmento de MRDS que es precisamente ajustar el multi-tensor. Como práctica, te recomiendo que abras todos los archivos con `mrview` y explores que es cada output. Pero por lo pronto, vemos que el comando `dti` nos generó los mapas de fracción de anisotropía (FA) y de difusividad media (MD), asi como el archivo de la *response function* `dwi_DTInolin_ResponseAsotropic.txt` que necesitaremos para el siguiente paso.


### Ajuste de MRDS segunda parte:

En esta segunda parte vamos a ajustar los multi-tensores voxel por voxel tomando algunos de los outputs de la primera parte. Esta segunda parte es un poco mas compleja en cuanto a los parámetros y también mucho mas tardado en correr. Veamos el manual del comando `mdtmrds`:

```
mdtmrds
 
mdtmrds [ options ] dwi_input scheme mt_output

        dwi_input
                name of the input dwi file (.nii).
        scheme
                name of the scheme (.txt) corresponding to the dwi image. Each line of the file must be:
                         x1 y1 z1 b1
                         x2 y2 z2 b2
                         x3 y3 z3 b3
                         .  .  .  .
                         .  .  .  .
                         .  .  .  .
                         xn yn zn bn
        mt_output
                name of the output multi-diffusion tensor file (.nii).

Compute the Multiple Radially Symmetryc Diffusion Tensor from Diffusion Weighted Magnetic Resonance Images.
Multi-Resolution Discrete-Search method is used.

Options:

         -mask file(.nii)
                mask for selecting the voxels that will be processed.
         -response l1,l2(,Diso)
                By default the eigenvalues l1 = 1.5e-3, l2 = 0.3e-3 are used as initial response function.
                Using this option, they could be specified so the eigenvalues of the tensor be [l1,l2,l2].
                If isotropic compartment flag is set, Diso can be supplied as the third value. Default 0.8e-3.
         -nonorm
                Work with original signals Si (not recommended). By default the algorithm works with A_i = S_i/S_0.
         -beta
                For the spatial bilateral estimation of S0^{s}_{snr}, S0^{s}_{std} and S0^{s}_{mean}.
                The estimation includes close voxels inside the range S0^{v}_{mean} +/- beta*S0^{v}_{std}.
                (s) -> spatial, (v) -> voxel.
                Default: beta = 0.5.
         -modsel name
                select the criteria used for estimate the number of fiber bundles (#param k) on each voxel.
                the choices are:
                        -modsel bic   -> Bayesian Information Critearia: BIC = -2 log L + n * log(k).
                        -modsel aic   -> Akaike Information Critearia:   AIC = -2 log L + 2 * k.
                        -modsel aicc  -> Akaike Information Critearia corrected for finite samples.
                        -modsel hqic  -> Hanan-Quinn Information Critearia: HQIC = -2 log L + 2k log(log(n))
                        -modsel l0    -> pseudo-l0 norm. n*MSE + lambda*sigma*sigma*log(n)*k (BIC: lambda = 1)
                        -modsel ftest -> F-test: F = [(RSS1-RSS2)/(p2-p1)] / [(RSS2)/(n-p2)] < pvalue.
                        -modsel bhq   -> compute BIC and HQIC.
                        -modsel all   -> compute all: bic, aic, aicc and hqic (sometimes needed).
                        Default: -ftest.
         -reg value
                If the option modsel is ftest, this option sets the pvalue for performing the f-test.
                If the option modsel is any other, this option will be ignored.
                Defaults: 0.001.
         -alt
                For the model selection, by default, the sigma dependent Gaussian Log Likelihood is used:
                'log L(sigma) = -(n/2)[ log(2pi) + log(sigma^2) + (1/sigma^2) MSE ]'.

                With this option a sigma free approximation of the Gaussian Log Likelihood is used instead:
                'log L = -(n/2) ln(MSE).

                If the option modsel is ftest, this option will be ignored.
         -method name
                Select the method to use (based on the estimation of the eigenvalues).
                        -method fixed: The eigenvalues are not estimated, they are kept fixed during the process.
                        -method equal: The eigenvalues are estimated equal for all the bundles inside the voxel.
                        -method diff:  The eigenvalues are estimated different for each bundle inside the voxel.
                Default: -method equal (is more stable).
         -iso
                Adding the isotropic compartment to the estimation (Still not working properly).
         -each
                By default, only the multi-tensor with the selected number of bundles is saved.
                With this option the multi-tensors with N = 1,2,3,... are also saved.
         -intermediate
                By default, only the multi-tensor(s) of the final stage are saved.
                With this option, the multi-tensor(s) of the intermediate stages are also saved.
         -stages np1,np2,np3,...
                Define the number of stages (resolutions of the orientation sets) used in the method.
                npi defines the number of orientations for the orientation set used in the ith stage.
         -fa
                compute the fractional anisotropy of the multi diffusion tensors on each voxel.
         -md
                compute the mean diffusivity of the multi diffusion tensors on each voxel.
         -mse
                compute the mse of the signals for the estimated multi diffusion tensors and the measured signal.
         -correction num
                Apply a transformation on the diffusion data in order to correct the Rician bias.
                This is not a denoising, it is just a bias correction. The choices are:
                        -correction  0, 1, 2 -> No correction: y = x.
                        -correction  3, 4, 5 -> Gudbjartsson correction: y = sqrt(fabs(x*x - sigma*sigma)).
                        -correction  6, 7, 8 -> Gudbjartsson correction modified: y = sqrt(max(x*x - sigma*sigma,0.0)).
                        -correction  9,10,11 -> Manjon correction: y = sqrt(max(x*x - 2.0*sigma*sigma,0.0)).
                        -correction 12,13,14 -> Maximum Likelihood correction y_ML = max_y(p(x;y,sigma)).
                                                Exhaustive evaluation.
                        -correction 15,16,17 -> Mean Posterior correction y = sum x*p(x).
                                                Exhaustive evaluation.
                        -correction 18,19,20 -> Mean Posterior correction adaptive according to each measurement.
                        Aditional consideration.
                        num % 3 => 0 - min value = 0.
                        num % 3 => 1 - min value = estimated from data.
                        num % 3 => 2 - min value = physically plausible exp(-b*3e-3).
                        Default: -correction 8.
         -help
                show this help

```


Como indica el manual en la opción de `response` (leer arriba), debemos de especificar los dos primeros lamdas. Esta información esta contenida en el archivo `dwi_DTInolin_ResponseAnisotropic.txt`. Para extraer esos datos y guardarlos en una variable puedes utilizar el siguiente código:

```python=
responsef=`cat dwi_DTInolin_ResponseAnisotropic.txt | cut -d" " -f1,2 --output-delimiter=,`
```
Hasta aquí ya tenemos los datos suficientes para correr el multi-tensor. Como puedes ver, hay diferentes opciones para correr la función, recuerda que hay que adaptarlo de acuerdo a lo que necesites. 
Sin embargo, algunos de los puntos son claves:
* La selección del modelo (`-modsel`), donde eliges como se resuelve el número de poblaciones de fibras por voxel. Donde BIC (Criterio de Información Bayesiano) es de los mas robustos,  
* El método para estimar los eigenvalores (`-method`), este punto es crucial si lo que buscas es que tus tensores sean completamente **independientes** entre sí, si es así, `diff` es la opción. 


El siguiente código es un ejemplo de como puedes ajustarlo:
```python=
mdtmrds
-mask ROI.nii \
-response $responsef \
-correction 0 \
-modsel bic \
-fa -md \
method diff 1 \
each \
64A_dw.nii \
scheme_nuevo.txt \
dwi.nii
```
> Ten mucha paciencia por que MRDS tarda horas en correr, incluso días, dependiendo de que tanto cerebro estes procesando....

Ahora vemos nuestros outputs:

```python=
 ls dwi_MRDS_Diff*
 
dwi_MRDS_Diff_BIC_COMP_SIZE.nii
dwi_MRDS_Diff_BIC_EIGENVALUES.nii
dwi_MRDS_Diff_BIC_FA.nii
dwi_MRDS_Diff_BIC_ISOTROPIC.nii
dwi_MRDS_Diff_BIC_MD.nii
dwi_MRDS_Diff_BIC_NUM_COMP.nii
dwi_MRDS_Diff_BIC_PDDs_CARTESIAN.nii
dwi_MRDS_Diff_V1_COMP_SIZE.nii
dwi_MRDS_Diff_V1_EIGENVALUES.nii
dwi_MRDS_Diff_V1_FA.nii
dwi_MRDS_Diff_V1_ISOTROPIC.nii
dwi_MRDS_Diff_V1_MD.nii
dwi_MRDS_Diff_V1_NUM_COMP.nii
dwi_MRDS_Diff_V1_PDDs_CARTESIAN.nii
dwi_MRDS_Diff_V2_COMP_SIZE.nii
dwi_MRDS_Diff_V2_EIGENVALUES.nii
dwi_MRDS_Diff_V2_FA.nii
dwi_MRDS_Diff_V2_ISOTROPIC.nii
dwi_MRDS_Diff_V2_MD.nii
dwi_MRDS_Diff_V2_NUM_COMP.nii
dwi_MRDS_Diff_V2_PDDs_CARTESIAN.nii
dwi_MRDS_Diff_V3_COMP_SIZE.nii
dwi_MRDS_Diff_V3_EIGENVALUES.nii
dwi_MRDS_Diff_V3_FA.nii
dwi_MRDS_Diff_V3_ISOTROPIC.nii
dwi_MRDS_Diff_V3_MD.nii
dwi_MRDS_Diff_V3_NUM_COMP.nii
dwi_MRDS_Diff_V3_PDDs_CARTESIAN.nii
```

Nuevamente te recomiendo mucho que explores tus outputs y veas que significan. Pero por lo pronto, podemos observar que hay cuatro sets de datos: `V1`, `V2`, `V3` y `BIC`. Donde `BIC` son el resultdado final una vez aplicado el criterio de información bayesiana y son los que al final puedes utilizar el análisis.

Entonces por el momento nos concentraremos en estos outputs:

```
ls dwi_MRDS_Diff_BIC*
 
dwi_MRDS_Diff_BIC_COMP_SIZE.nii
dwi_MRDS_Diff_BIC_EIGENVALUES.nii
dwi_MRDS_Diff_BIC_FA.nii
dwi_MRDS_Diff_BIC_ISOTROPIC.nii
dwi_MRDS_Diff_BIC_MD.nii
dwi_MRDS_Diff_BIC_NUM_COMP.nii
dwi_MRDS_Diff_BIC_PDDs_CARTESIAN.nii
```



