dMRIharmonization repository is developed by Tashrif Billah, Sylvain Bouix, Suheyla Cetin Karayumak, and Yogesh Rathi, Brigham and Women's Hospital (Harvard Medical School).

• Table of Contents
• Multi-site dMRI harmonization
• Citation
• Dependencies
• Installation
  • 1. Install prerequisites
    • Check system architecture
    • Python 3
    • unringing
  • 2. Install pipeline
  • 3. Download IIT templates
  • 4. Configure your environment
• Running
• Tests
  • 1. pipeline
  • 2. unittest
• List of images
• Site names
• Multi threading
• Order of spherical harmonics
• Number of zero padding
• NRRD support
• Preprocessing
  • 1. Denoising
  • 2. Bvalue mapping
  • 3. Resampling
• Config
• Template
• List of outputs
  • 1. Folders
  • 2. Files
• Template creation
• Data harmonization
• Debugging
  • 1. With the pipeline
  • 2. Use seperately
  • 3. With a list of FA images
• Travel heads
• Caveats/Issues
  • 1. Template path
  • 2. Multi-processing
  • 3. Tracker
• Reference

Table of Contents created by gh-md-toc

Integrated study of multi-site diffusion MRI (dMRI) data can enable diagnosis, monitoring, and treatment of several brain diseases. However, data acquired on a variety of scanners cannot be integrated due to differences in acquisition parameters and scanner artifacts. Therefore, dMRI data has to be harmonized for joint analysis by removing scanner-specific differences.

dMRIharmonization is a Python command line module that implements a method capable of removing scanner-specific effects. The proposed method eliminates inter-site variability in acquisition parameters, while preserves inter-subject anatomical variability.

See Reference for more details

The method leverages on rotation invariant spherical harmonics (RISH) features derived from spherical harmonic coefficients. In brief, the method constructs a scale map for each pair of reference and target sites. Applying the scale map to the target site, it aims to remove inter-site variability.

If this repository is useful in your research, please cite all of the following:

• Billah T*, Cetin Karayumak S*, Bouix S, Rathi Y. Multi-site Diffusion MRI Harmonization, https://github.com/pnlbwh/dMRIharmoniziation, 2019, doi: 10.5281/zenodo.2584275

  * denotes equal first authroship

• Cetin Karayumak S, Bouix S, Ning L, James A, Crow T, Shenton M, Kubicki M, Rathi Y. Retrospective harmonization of multi-site diffusion MRI data acquired with different acquisition parameters. Neuroimage. 2019 Jan 1;184:180-200. doi: 10.1016/j.neuroimage.2018.08.073. Epub 2018 Sep 8. PubMed PMID: 30205206; PubMed Central PMCID: PMC6230479.

• Mirzaalian H, Ning L, Savadjiev P, Pasternak O, Bouix S, Michailovich O, Karmacharya S, Grant G, Marx CE, Morey RA, Flashman LA, George MS, McAllister TW, Andaluz N, Shutter L, Coimbra R, Zafonte RD, Coleman MJ, Kubicki M, Westin CF, Stein MB, Shenton ME, Rathi Y. Multi-site harmonization of diffusion MRI data in a registration framework. Brain Imaging Behav. 2018 Feb;12(1):284-295. doi:10.1007/s11682-016-9670-y. PubMed PMID: 28176263.

• ANTs = 2.2.0
• reisert/unring
• FSL = 5.0.11
• numpy = 1.16.2
• scipy = 1.2.1
• scikit-image = 0.15.0
• dipy = 0.16.0
• nibabel = 2.3.0
• pynrrd = 0.3.6
• conversion = 2.0

NOTE The above versions were used for developing the repository. However, dMRIharmonization should work on any advanced version.

You may ignore installation instruction for any software module that you have already.

uname -a # check if 32 or 64 bit

Download Miniconda Python 3.6 bash installer (32/64-bit based on your environment):

sh Miniconda3-latest-Linux-x86_64.sh -b # -b flag is for license agreement

Activate the conda environment:

source ~/miniconda3/bin/activate # should introduce '(base)' in front of each line

git clone https://bitbucket.org/reisert/unring.git
cd unring/fsl
export PATH=$PATH:`pwd`
unring --help

Now that you have installed the prerequisite software, you are ready to install the pipeline:

git clone https://github.com/pnlbwh/dMRIharmonization.git && cd dMRIharmonization
conda env create -f environment.yml    # you may comment out any existing package from environment.yml
conda activate harmonization           # should introduce '(harmonization)' in front of each line

Alternatively, if you already have ANTs, you can continue using your python environment by directly installing the prerequisite libraries:

pip install -r requirements.txt --upgrade

dMRIharmonization toolbox is provided with a debugging capability to test how good has been the harmonization. For debug to work and tests to run, download the following data from IIT HUMAN BRAIN ATLAS and place them in IITAtlas/ directory:

• IITmean_FA.nii.gz
• IITmean_FA_skeleton.nii.gz

Make sure the following executables are in your path:

antsMultivariateTemplateConstruction2.sh
antsApplyTransforms
antsRegistrationSyNQuick.sh
unring.a64

You can check them as follows:

which dtifit

If any of them does not exist, add that to your path:

export PATH=$PATH:/directory/of/executable

conda activate harmonization should already put the ANTs scripts in your path. Yet, you should set ANTSPATH as follows:

export ANTSPATH=~/miniconda3/envs/harmonization/bin/

However, if you choose to use pre-installed ANTs scripts, you can define ANTSPATH according to this instruction.

Upon successful installation, you should be able to see the help message

$ lib/harmonization.py --help

Template creation, harmonization, and debugging

Usage:
    harmonization.py [SWITCHES] 

Meta-switches:
    -h, --help                          Prints this help message and quits
    --help-all                          Prints help messages of all sub-commands and quits
    -v, --version                       Prints the program's version and quits

Switches:
    --bvalMap VALUE:str                 specify a bmax to scale bvalues into
    --create                            turn on this flag to create template
    --debug                             turn on this flag to debug harmonized data (valid only with --process)
    --denoise                           turn on this flag to denoise voxel data
    --force                             turn on this flag to overwrite existing data
    --harm_list VALUE:ExistingFile      harmonized csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1
                                        dwi2,mask2 ...
    --nproc VALUE:str                   number of processes/threads to use (-1 for all available, may slow down your system);
                                        the default is 8
    --nshm VALUE:str                    spherical harmonic order; the default is 6
    --nzero VALUE:str                   number of zero padding for denoising skull region during signal reconstruction; the default is 10
    --process                           turn on this flag to harmonize
    --ref_list VALUE:ExistingFile       reference csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1
                                        dwi2,mask2 ...
    --ref_name VALUE:str                reference site name; required
    --resample VALUE:str                voxel size MxNxO to resample into
    --tar_list VALUE:ExistingFile       target csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1
                                        dwi2,mask2 ...
    --tar_name VALUE:str                target site name; required
    --template VALUE:str                template directory; required
    --travelHeads                       travelling heads

The harmonization.py cli takes in the following arguments that are explained below.

• --bvalMap VALUE:str
• --create
• --debug
• --denoise
• --force
• --harm_list VALUE:ExistingFile
• --nproc VALUE:str
• --nshm VALUE:str
• --nzero VALUE:str
• --process
• --ref_list VALUE:ExistingFile
• --ref_name VALUE:str
• --resample VALUE:str
• --tar_list VALUE:ExistingFile
• --tar_name VALUE:str
• --template VALUE:str
• --travelHeads

A small test data is provided with each release.

You may test the whole pipeline as follows*:

cd lib/tests
./test_pipeline.sh

NOTE: running the above test should take roughly an hour.

./pipeline_test.sh will download test data*, and run the whole processing pipeline on them. If the test is successful and complete, you should see the following output on the command line.

CONNECTOM mean FA:  0.675917931134666
PRISMA mean FA before harmonization:  0.8008729791383536
PRISMA mean FA after harmonization:  0.6366103024088171

* If there is any problem downloading test data, try manually downloading and unzipping it to lib/tests/ folder.

You may run smaller and faster unittest as follows:

python -m unittest discover -v lib/tests/    

Two image lists (.txt or .csv) are required for the reference and target sites:

--ref_list cidar_imgs_masks.csv
--tar_list mrc_imgs_masks.txt

Each \n delimited line in the list file should have , seperated dwi and mask:

case001_dwi.nii.gz,case001_mask.nii.gz
case002_dwi.nrrd,case002_mask.nii.gz
case003_dwi.nii,case003_mask.nhdr
...
...

It is assumed that .bval and .bvec are in the same location and have the same prefix as that of the dwi.nii.gz image.

During template construction, you want to provide a small list of images to --ref_list and --tar_list. However, provide a complete list of images you want to harmonize to --tar_list during data harmonization

Site names are used to properly identify files in the template:

 --ref_name CIDAR
 --tar_name MRC

Processing can be multi-threaded over the cases. Besides, antsMultivariateTemplateConstruction2.sh utilizes multiple threads to speed-up template construction.

--nproc 8 # default is 4, use -1 for all available

However, multi-threading comes with a price of slowing down other processes that may be running in your system. So, it is advisable to leave out at least two cores for other processes to run smoothly.

NOTE See Caveats/Issues that may occur while using many processors in parallel.

RISH features are derived from spherical harmonic coefficients. The order of spherical harmonic coefficients you can use is limited by the lowest number of gradients present in the diffusion images between reference and target sites.

During signal reconstruction, skull region can be noisy. So, each RISH feature is zero padded, median filtered in the skull region, and then used for signal reconstruction. The more noisy data you have, the more number of zeros you want to pad. The value should be 5-20.

--nzero 10

The pipeline is written for NIFTI image format. However, NRRD support is incorporated through NIFTI --> NRRD conversion on the fly.

See Tashrif Billah, Sylvain Bouix and Yogesh Rathi, Various MRI Conversion Tools, https://github.com/pnlbwh/conversion, 2019, DOI: 10.5281/zenodo.2584003 for more details on the conversion method.

There are three preprocessing steps. More than one step can be specified at a time. However, they are performed in the following order.

--denoise        # turn on this flag to denoise voxel data

--bvalMap VALUE  # specify a bmax to scale bvalues into    

--resample VALUE # voxel size MxNxO to resample into

After preprocessing, the image lists are saved with .modified extension in the same location of provided lists, and used for further processing.

With all the arguments passed to harmonization.py cli, a config.ini file is created in lib/. This file is used by submodules of the harmonization pipeline during executation. A sample config.ini file looks like the following:

[DEFAULT]
N_shm = 6
N_proc = 8
N_zero = 10
resample = 1.5x1.5x1.5
bvalMap = 1000
denoise = 0
travelHeads = 0
debug = 1
diffusionMeasures = MD,FA,GFA

Template directory is provided with --template argument. This directory is primarily used by antsMultivariateTemplateConstruction2.sh to save transform files. Then, various scales maps (Scale_*.nii.gz) and mean templates (Mean_*.nii.gz) are created in this directory. There are other template files: Delta_*.nii.gz, PercentageDiff_*.nii.gz, and PercentageDiff_*smooth.nii.gz. You may look at them to know how good are the created templates.

Several files are created down the pipeline. They are organized with proper folder hierarchy and naming:

In each input image directory, two folders are created: dti and rish. The dti folder stores diffusion measures and corresponding transform files. On the other hand rish folder stores RISH features and corresponding transform files.

The harmonization pipeline will generate a bunch of diffusion measures and RISH features. All the measures are saved as .nii.gz format with associated .bval and .bvec where necessary.

Firstly, see the files with prefix harmonized_ in target site image directories. They are the harmonized diffusion weighted images.

Secondly, you can look at files with suffix _denoised.nii.gz, _bmapped.nii.gz, and _resampled.nii.gz. Denoising, bvalue mapping, and resampling are performed in the above order and files are saved with the last suffix unless --debug option is used. In the latter case, files after each preprocessing step are saved.

Finally, there are other files:

RISH features and diffusion measures are computed from all the input images. Images of two modalites: FA and L0 (RISH order 0) are provided as input to antsMultivariateTemplateConstruction2.sh. A template0.nii.gz is created. Afterwards, various scales maps (Scale_*.nii.gz) and mean templates (Mean_*.nii.gz) are created in --template directory. See Reference for details on the template construction method and Template for list of outputs. A sample template construction command is given below:

lib/harmonization.py \
--ref_list test_data/ref_caselist.txt \
--tar_list test_data/target_caselist.txt \
--ref_name REF \
--tar_name TAR \
--template test_data/template/ \
--nshm 6 \
--bvalMap 1000 \
--resample 1.5x1.5x1.5 \
--nproc 8 \
--create

Note: Replace each of the above paths with absolute path as needed.

During template construction, you want to provide a small list of images to --ref_list and --tar_list. Section 2.4 of our paper gives a concise analysis about the number of samples to be chosen for template creation. The paper posits 20 as the "Gold Standard" for template creation while 16 as the minimum number. The rule of thumb should be to choose 16-20 subjects matched across sites for age, gender and IQ. However, provide a complete list of images you want to harmonize to --tar_list during data harmonization

After template creation is completed, there should be scales maps (Scale_*.nii.gz) in the template directory. The scale maps are used to scale RISH features of the reference site and then reconstruct diffusion weighted images. See Reference for details on the signal reconstruction.

lib/harmonization.py \
--tar_list test_data/target_caselist.txt \
--tar_name TAR \
--template test_data/template/ \
--nshm 6 \
--bvalMap 1000 \
--resample 1.5x1.5x1.5 \
--nproc 8 \
--process

Note: Replace each of the above paths with absolute path as needed.

During template construction, you want to provide a small list of images to --ref_list and --tar_list. However, provide a complete list of images you want to harmonize to --tar_list during data harmonization

dMRIharmonization toolbox is provided with a debugging capability to test how good has been the harmonization. --debug can be run with any (or all) of --create and --process options or separately.

Once harmonization is done, there are three types of data:

• reference site
• target site
• harmonized target

If the three data sets are brought into a common space, then looking at the mean FA over the whole brain, we can comment on the goodness of harmonization. If data is properly harmonized, then mean FA of the harmonized target should be almost equal to the mean FA of that of reference site.

In details:

• reference data is proprocessed, and registered to reference template space and then to MNI space (IITmean_FA.nii.gz)
• unprocessed target data is directly registered to MNI space
• harmonized target data is registered to target template space and then to MNI space
• once the data are in MNI space, we calculate mean FA over the IITmean_FA_skeleton.nii.gz

NOTE: Download the above data from IIT HUMAN BRAIN ATLAS and place them in IITAtlas/ directory.

The numbers should like like below:

Printing statistics :
CIDAR mean FA:  0.5217237675408243
BSNIP mean FA before harmonization:  0.5072286796848892
BSNIP mean FA after harmonization:  0.5321998242139347

As we see above, BSNIP (target) mean FA over IITmean_FA_skeleton.nii.gz for all the site images after harmonization increased to be almost equal to that of CIDAR (reference) mean FA.

Now there are two ways to debug:

Use --debug flag with any (or all) of --create and --process

If you would like to debug at a later time, you need to specify three images lists:

• --ref_list: use the reference list with .modified extension

• --tar_list: use the unprocessed target list without the .modified extension

• --harm_list: use the harmonized target list that has .harmonized extension

  lib/harmonization.py --ref_list ref.txt.modified --tar_list target.csv --harm_list target.csv.modified.harmonized

NOTE: You should run the pipeline first before debugging separately because --debug makes use of files created in the pipeline.

The repository provides a more discrete discrete script for finding the goodness of harmonization.

$ lib/tests/fa_skeleton_test.py --help

usage: fa_skeleton_test.py [-h] -i INPUT -s SITE -t TEMPLATE

Warps diffusion measures (FA, MD, GFA) to template space and then to subject
space. Finally, calculates mean FA over IITmean_FA_skeleton.nii.gz

optional arguments:
  -h, --help            show this help message and exit
  -i INPUT, --input INPUT
                        input list of FA images
  -s SITE, --site SITE  site name for locating template FA and mask in
                        tempalte directory
  -t TEMPLATE, --template TEMPLATE
                        template directory where Mean_{site}_FA.nii.gz and
                        {site}_Mask.nii.gz is located

This script does not depend of registration performed during the harmonization process. Rather, it performs all the steps mentioned above (Debugging: In details) and computes mean FA over skeleton across all subjects in a site.

Traveling heads mean same set of subjects scanned with different scanners. For example, scanning same set of people in Brigham and Women's Hospital and Massachusetts General Hospital. Therefore, users should set --travelHeads if they have same set of subjects across sites.

antsMultivariateTemplateConstruction2.sh: all the images need to have unique prefix because transform files are created in the same --template ./template/ directory. The value of --template should have / at the end. The pipeline appends one if there is not, but it is good to include it when specifying.

Multi threading requires memory and processor availability. If pipeline does not continue past unringing or shm_coeff computation, your machine likely ran out of memory. Reducing --nproc to lower number of processors (i.e. 1-4) or swithcing to a powerful machine should help in this regard.

In any case, feel free to submit an issue here. We shall get back to you as soon as possible.

Zhang S, Arfanakis K. Evaluation of standardized and study-specific diffusion tensor imaging templates of the adult human brain: Template characteristics, spatial normalization accuracy, and detection of small inter-group FA differences. Neuroimage 2018;172:40-50.

Billah, Tashrif; Bouix, Sylvain; Rathi, Yogesh; Various MRI Conversion Tools, https://github.com/pnlbwh/conversion, 2019, DOI: 10.5281/zenodo.2584003.