SEPIA Documentation¶

Welcome to the SEPIA documentation! Here you can find all the documents related to SEPIA.

SEPIA provides a graphical user interface in Matlab to build data processing pipelines related to quantitative susceptibility mapping (QSM).

SEPIA is designed to solve two issues I encountered during my research in QSM:

Using algorithms from toolboxes developed by different research groups,
Having a platform such that I can quickly adjust (and remember) parameters of each algorithm.

The purpose of SEPIA is to provide a platform for easy access to different QSM processing methods in the field. It provides interfaces to use the following three toolboxes freely available online for academic purposes, i.e.

Through SEPIA I hope researchers who are not coming from the field could also be able to use QSM for their research. For those more experienced researchers, I also hope this tool can simplify your workflow.

Table of Contents¶

Installation¶

Prerequisite¶

To use SEPIA you need three QSM toolboxes which are freely available for the research community, namely MEDI toolbox, STI Suite and FANSI toolbox. You can download these toolboxes with the following links:

MEDI toolbox (version 25 June 2018)
STI Suite (version 3.0)
FANSI toolbox (commit d33759b9)

If you encounter any difficulty when downloading these toolboxes please let me know.

Installing SEPIA¶

Once you have all the toolboxes in place, then you have to specify the directory of each toolbox in SpecifyToolboxesDirectory.m:

MEDI_dir = '/path/to/MEDI/toolbox/';
STISuite_dir = '/path/to/STISuite/toolbox/';
FANSI_dir = '/path/to/FASI/toolbox/';

For example, I have all my external toolboxes stored under the SEPIA home directory. Additionally, for each toolbox, I have different copies representing different versions when they were published

and here is an example of how is my SpecifyToolboxesDirectory.m defined:

% 1. sepcify the toolbox version you want to run
MEDI_version        = 'MEDI_toolbox_20180625';
FANSI_version       = 'FANSI-toolbox-d33759b970790cc8754adc9d0398cc3d07546074/';
STISuite_version    = 'STISuite_V3.0';

% 2. get the Sepia HOME directory from this script
fullName        = mfilename('fullpath');
SEPIA_HOME      = fileparts(fullName);

% 3. specify the top level of external dependency directory
external_dir    = [SEPIA_HOME filesep 'external' filesep];

% 4. specify the HOME directory of each toolbox
MEDI_HOME       = [external_dir 'MEDI_toolbox' filesep];
FANSI_HOME      = [external_dir 'FANSI_toolbox' filesep];
STISuite_HOME   = [external_dir 'STI_Suite' filesep];

% 5. sepcify the final destination of each toolbox you want to run in Sepia
MEDI_dir        = [MEDI_HOME MEDI_version filesep];
FANSI_dir       = [FANSI_HOME FANSI_version filesep];
STISuite_dir    = [STISuite_HOME STISuite_version filesep];

IMPORTANT: Please do not modify the original structure of these toolboxes, SEPIA searches the path of the related functions based on the original folder structure.

Finally, you have to add the directory containing sepia.m, i.e. the SEPIA HOME directory, to your MATLAB PATH

This can be done by: ‘Set Path’ -> ‘Add Folder’ -> /your/sepia/directory/ -> ‘Save’

(Caution: To ensure only the selected algorithm is used for the QSM processing, please do not manually add the paths to the external toolboxes you want to run in SEPIA to the Matlab PATH, the sepia_addpath function will do the job for you:).)

or

with MATLAB’s command: addpath('/your/sepia/directory')

Now you can start the GUI by entering sepia in the MATLAB’s command window.

Compatibility¶

SEPIA is developed mainly in MATLAB R2016b in Linux. In general, all methods should compatible with earlier MATLAB versions. Most of the methods should also compatible with MATLAB R2017a or later, and other OS, except the following functions/algorithms

DICOM reader
Laplacian Boundary Value (LBV) for background field removal
Graphcut for phase unwrapping

Data Preparation¶

SEPIA supports the following type of data input:

uncompressed or compressed NIfTI images (.nii and .nii.gz)

DICOM images¶

DICOM images input has been deprecated since v.0.7.0.

NIfTI images¶

NIfTI data is a handy way to work with SEPIA. SEPIA is mainly designed for multi-echo GRE (mGRE) data. Therefore, the main standalone of SEPIA can only work with 4D (row, column, slice, time) NIfTI data (3D for some standalone applications).

Data conversion¶

So far there are three data conversion tools tested in SEPIA:

One way to convert DICOM images into 4D NIfTI data is to use MRIConvert with the following setting checked: ‘Option’ -> ‘Save multivolumes series as 4D files’

In this way, your mGRE data will be stored as 4D FSL NIfTI data that is a valid input of SEPIA.

Alternatively, you can also use dcm2niix to convert your data. Please make sure the ‘merge’ option (-m) is set to ‘no’ for the conversion (i.e. dcm2niix -m n). In this way, multiple 3D volumes (the number of volumes depends on the number of echoes acquired) will be created together with the JSON files containing the TE of each echo (if you enable the merge option of dcm2niix you will only have one JSON file containing one TE). You can then merge the echo images into 4D using tools like fslmerge.
File formats generated by dicm2nii are also supported in SEPIA.

Data input¶

Once you have the NIfTI files ready, SEPIA provides two options to load the data:

select the required files directly, or
prepare the data with specific names and put all of them in a common directory, from which you can specify the input directory in SEPIA. The name requirement depends on the standalone you are working on. For more specific details please check the wiki pages of each standalone applications.

SEPIA header¶

To obtain a frequency shifted map or magnetic susceptibility map with correct units, it is important to inform the corresponding algorithm(s) with information including the main magnetic field strength (B0) of the scanner, the main magnetic field direction with respect to the field of view (FOV) and the echo times (TEs) used in the acquisition. Unfortunately, such information may not be (directly) available in the NIfTI header after the DICOM conversion. As a result, if you are using NIfTI files with SEPIA, you have to also supply a special header file alongside with your multi-echo GRE data.

SEPIA requires a specific MAT-file (.mat) that stores some header information of the input data. The following example shows the variables (case-sensitive) that must be stored in the header file:

B0 = 3;                 % magnetic field strength, in Tesla
B0_dir = [0;0;1];       % main magnetic field direction, [x,y,z]
CF = 3*42.58*1e6;       % imaging frequency, in Hz (B0*gyromagnetic_ratio*1e6)
TE = [te1,te2,te3,te4]; % echo time for each GRE image, in second
delta_TE = te2-te1;     % echo spacing, in second
matrixSize = [64,64,64] % image matrix size
voxelSize = [1,1,1]     % spatial resolution of the data, in mm

You can create this file manually in Matlab. Alternatively, you can use the Utility standalone application to generate the file. Please visit the corresponding wiki page for more information about how to generate the header file.

SEPIA (One-stop QSM processing)¶

What is SEPIA?¶

SEPIA is a pipeline analysis tool for quantitative susceptibility mapping (QSM) in Brain Imaging. It provides all the essential functions people would need to compute a susceptibility map from a 3D multi-echo GRE phase data, including phase unwrapping, background field contribution removal and dipole inversion. Incorporate with different toolboxes in SEPIA giving users the advantages of having a variety of options to build a pipeline that works the best for their data. When you use the SEPIA graphical user interface to process the data, a log file will be generated that contains also the setting and command that you’ve chosen in the pipeline. This log file will be particularly useful for batch processing.

Structure of the application¶

This standalone consists of 4 panels:

Input/Output(I/O) panel,
Total field recovery and phase unwrapping panel,
Background field removal panel, and
QSM panel.

The detailed description of each panel is given below:

I/O panel¶

Data input

This application accepts three types of data input method:
1. Specify a directory that contains all NIfTI images. Please specify the names of your data as in the following:
  - Phase: must contain the string ‘ph’ in the filename, e.g. phase.nii.gz;
  - Magnitude: must contain the string ‘mag’ in the filename, e.g. magn.nii.gz;
  - Header: must contain the string ‘header’ in the filename, e.g. header.mat;
  - (optional) Mask data: if provided, must contain string ‘mask’ in the filename, e.g. mask.nii.gz, or
2. Specify the required data separately using the GUI buttons.
Data output

You can specify the prefix of the data output name in the editable field ‘Output basename’. By default, the SEPIA output will be stored in a directory named ‘output’ under the input directory, i.e. ‘_/your/input/directory/output/_’ with prefix ‘sepia’. You can change the default output directory and basename to whatever you need. If the output directory does not exist, the application will create the directory.
Brain mask

QSM-related algorithms often require a mask that contains only brain tissue. If you already have the brain mask data in NIfTI format, you can select the file manually, or named it with a specific name (see Data input section) and put it in the directory with other NIfTI files. Alternatively, you can check the ‘FSL brain extraction’ checkbox to extract the brain mask. Enabling this option will run the Matlab implementation of FSL’s brain extraction tool (bet) implemented with MEDI toolbox.
Invert phase data

Due to the way of how the phase data is stored, values of the final local field map and/or magnetic susceptibility map might be inverted, i.e. positive frequency shift and paramagnetic magnetic susceptibility will show in negative, and vice versa. If it is the case, you can invert the values of the results by checking the ‘Invert phase data’ option. This will apply a conjugate operation to the phase data after the data being loaded.

Total field recovery and phase unwrapping panel¶

Echo phase combination

As the first step to process the multi-echo data, we need to recover the total frequency shift of the tissue across times. SEPIA provides two different ways to do this:
1. Optimum weights
This is a weighted combination of the phase difference between successive echoes, in which the weights are inversely proportional to the variance of the noise of the fieldmap estimated from the magnitude echo images.
1. MEDI nonlinear fit
This is a method in the MEDI toolbox
Phase unwrapping

There are 5 phase unwrapping method supported in SEPIA
1. Laplacian
Laplacian unwrapping implementation in MEDI toolbox
1. Laplacian STI suite
Laplacian unwrapping implementation in STI Suite v3.0
1. 3D best path
Robust region growing method yet only works in the DCCN cluster (recommended if you use this toolbox in the DCCN cluster)
1. Region growing
Region growing method in the MEDI toolbox
1. Graphcut
Graph-cut algorithm in the MEDI toolbox, sometimes uses with water-fat imaging.
Bipolar readout eddy current correction:

enable to correct the phase inconsistency between odd and even echoes, and a gradient-like field contribution by eddy current effect due to bipolar readout. If this option is enabled, the eddy current corrected data will be stored in the output directory with the following name:
- phase_eddy-correct.nii.gz (eddy current corrected phase data)
Exclude unreliable voxels, Threshold:

enable to exclude low SNR voxels that can create strong artefacts in susceptibility map (you may check with ‘relative-residual.nii.gz’ to adjust the threshold). Voxels that have relative fitting residual greater than the threshold will be weighted with zero in subsequent processes. Only available for region growing and 3D best path unwrapping methods.
Output

The output of this step are given below:
- total-field.nii.gz (unwrapped total (background+local) field, in Hz)
- fieldmap-sd.nii.gz (normalised field map standard deviation)
- mask.nii.gz (FSL’s bet brain mask, optional)
- mask-reliable.nii.gz (thresholded brain mask, optional)
- relative-residual.nii.gz (relative residual of fitting a mono-exponential decay function with a single frequnecy shift, depends on unwrapping method)

Background field removal panel¶

Method
1. LBV
  
  Laplacian boundary value approach to removal background field
2. PDF
  
  Projection onto dipole field
3. RESHARP
  
  regularisation enabled SHARP
4. SHARP
  
  Sophisticated harmonic artefact reduction for phase data
5. VSHARP STI suite
  
  STI suite v3.0 variable-kernel SHARP
6. VSHARP
7. iHARPERELLA
  
  (not optimised with SEPIA yet)
Refine local field by 4th order 3D polynomial fit

Enable to remove residual B1(+ & -) contribution in the local field
Output

The output of this step are given below:
- local-field.nii.gz (local (or tissue) field, in Hz)
- mask-qsm.nii.gz (brain mask where local field is reliable, might be eroded and depended on the background field removal algorithms and ‘exclude unreliable voxels’ threshold value)

QSM panel¶

Method:
1. TKD
  
  Thresholded k-space division
2. Closed-form solution
  
  closed-form solution with L2-norm regularisation
3. STI suite iLSQR
  
  STI suite v3.0 implementation of iterative LSQR approach
4. iLSQR
5. FANSI
  
  Fast algorithm for nonlinear susceptibility inversion
6. Star
  
  STI suite v3.0 Star-QSM (recommended)
7. MEDI
  
  Morphology enabled dipole inversion (MEDI+0)
Output

The output of this step is given below:
- QSM.nii.gz (quantitative susceptibility map, in ppm)

Phase Unwarpping Standalone¶

Why do we need phase unwrapping?¶

Phase wrapping occurs when continuous phase information is sampled in a discrete wrapped phase. The measured phase accumulation larger than one phase cycle is wrapped into the interval [- $\pi$ , $\pi$ ), causing the discontinuity in the phase data. To recover the true phase values, one must solve this ambiguity problem by adding the correct integer number of phase cycles to the phase data in order to recover the true phase revolution.

The objective of this standalone application is to recover the actual, total phase shift of the acquired data.

Structure of the application¶

This standalone consists of two panels:

I/O panel, and
Total field recovery and phase unwrapping panel.

The detailed description of each panel is given below:

I/O panel¶

Data input

This application accepts three types of data input method:
1. Specify a directory that contains all DICOM images,
2. Specify a directory that contains all NIfTI images. Please specify the names of your data as in the following: - Phase: must contain the string ‘ph’ in the filename, e.g. phase.nii.gz; - Magnitude: must contain the string ‘mag’ in the filename, e.g. magn.nii.gz; - Header: must contain the string ‘header’ in the filename, e.g. header.mat; - (optional) Mask data: if provided, must contain string ‘mask’ in the filename, e.g. mask.nii.gz, or
3. Specify the required data separately using the GUI buttons.
Data output

You can specify the prefix of the data output name in the editable field ‘Output basename’. By default, the SEPIA output will be stored in a directory named ‘output’ under the input directory, i.e. ‘_/your/input/directory/output/_’ with prefix ‘sepia’. You can change the default output directory and basename to whatever you need. If the output directory does not exist, the application will create the directory.
Brain mask

QSM related algorithms often require a mask that contains only brain tissue. If you already have the brain mask data in NIfTI format, you can select the file manually, or named it with a specific name (see Data input section) and put it in the directory with other NIfTI files. Alternatively, you can check the ‘FSL brain extraction’ checkbox to extract the brain mask. Enabling this option will run the Matlab implementation of FSL’s brain extraction tool (bet) implemented with MEDI toolbox.
Invert phase data

Due to the way of how the phase data is stored, values of the final local field map and/or magnetic susceptibility map might be inverted, i.e. positive frequency shift and paramagnetic magnetic susceptibility will show in negative, and vice versa. If it is the case, you can invert the values of the results by checking the ‘Invert phase data’ option. This will apply a conjugate operation to the phase data after the data being loaded.

Total field recovery and phase unwrapping panel¶

_images/phase-unwrapping_total-field.png

Echo phase combination

As the first step to process the multi-echo data, we need to recover the total frequency shift of the tissue across times. SEPIA provides two different ways to do this:
1. Optimum weights
This is a weighted combination of the phase difference between successive echoes, in which the weights are inversely proportional to the variance of the noise of the fieldmap estimated from the magnitude echo images.
1. MEDI nonlinear fit
This is a method in the MEDI toolbox
Phase unwrapping

There are 5 phase unwrapping method supported in SEPIA
1. Laplacian
Laplacian unwrapping implementation in MEDI toolbox
1. Laplacian STI suite
Laplacian unwrapping implementation in STI Suite v3.0
1. 3D best path
Robust region growing method yet only works in the DCCN cluster (recommended if you use this toolbox in the DCCN cluster)
1. Region growing
Region growing method in MEDI toolbox
1. Graphcut
Graph-cut algorithm in MEDI toolbox, usually uses in water-fat imaging.
Bipolar readout eddy current correction:

enable to correct the phase inconsistency between odd and even echoes, and a gradient-like field contribution by eddy current effect due to bipolar readout. If this option is enabled, the eddy current corrected data will be stored in the output directory with the following name:
- phase_eddy-correct.nii.gz (eddy current corrected phase data)
Exclude unreliable voxels, Threshold:

enable to exclude low SNR voxels that can create strong artefacts in susceptibility map (you may check with ‘relative-residual.nii.gz’ to adjust the threshold). Voxels that have relative fitting residual greater than the threshold will be weighted with zero in subsequent processes. Only available for region growing and 3D best path unwrapping methods.
Output

The output of the standalone application is given below:
- total-field.nii.gz (unwrapped total (background+local) field, in Hz)
- fieldmap-sd.nii.gz (normalised field map standard deviation)
- mask.nii.gz (FSL’s bet brain mask, optional)
- mask-reliable.nii.gz (thresholded brain mask, optional)
- relative-residual.nii.gz (relative residual of fitting a mono-exponential decay function with a single frequnecy shift, optional)

Background Magnetic Field Removal Standalone¶

Why do we need to remove background field?¶

The phase we measured in a GRE acquisition is affected by not only the brain tissue but also sources like B0 inhomogeneity and air sinus. In order to compute the susceptibility sources only contributed by the brain tissue, it is important to remove all the non-local field effect. Fortunately, the characteristic of the local field is different from that of the non-local field, it is possible to separate the two fields from the unwrapped total field. This standalone application is designed to solve the field separation problem.

Caution: It is crucial that the background field contribution is removed accurately in this stage. Otherwise, the remaining field due to background sources will be treated as part of the local field, degrading the quality of QSM result.

Structure of the application¶

This application consists of two panels:

I/O panel, and
Background field removal panel.

The detailed description of each panel is given below:

I/O panel¶

Data input

This application accepts two types of data input method:
1. Specify a directory that contains all NIfTI images. Please specify the names of your data as in the following:
  - Total field data: must contain the string ‘total-field’ in the filename, e.g. total-field.nii.gz;
  - Header: must contain the string ‘header’ in the filename, e.g. header.mat;
  - (optional) Fieldmap standard deviation data: must contain the string ‘fieldmap-sd’ in the filename, e.g. fieldmap-sd.nii.gz;
  - (optional) Mask data: if provided, must contain string ‘mask’ in the filename, e.g. mask.nii.gz, or
2. Specify the required data separately using the GUI buttons.
Data output

You can specify the prefix of the data output name in the editable field ‘Output basename’. By default, the Sepia output will be stored in a directory named ‘output’ under the input directory, i.e. ‘/your/input/directory/output/’ with prefix ‘sepia’. You can change the default output directory and basename to whatever you need. If the output directory does not exist, the application will create the directory.
Brain mask

QSM related algorithms often require a mask that contains only brain tissue. If you already have the brain mask data in NIfTI format, you can select the file manually, or named it with a specific name (see Data input section) and put it in the directory with other NIfTI files.

Background field removal panel¶

_images/background-field-removal_background-field.png

Method
1. LBV
Laplacian boundary value approach to removal background field
1. PDF
Projection onto dipole field
1. RESHARP
Regularisation enabled SHARP
1. SHARP
Sophisticated harmonic artefact reduction for phase data
1. VSHARP STI suite
STI suite v3.0 variable-kernel SHARP
1. VSHARP
2. iHARPERELLA
(not optimised with SEPIA yet)
Refine local field by 4th order 3D polynomial fit

Enable to remove residual B1(+ & -) contribution in the local field
Output

The output of the standalone application is given below:
- local-field.nii.gz (local (or tissue) field, in Hz)
- mask-qsm.nii.gz (brain mask where local field is reliable, might be eroded and depended on the background field removal algorithms and ‘exclude unreliable voxels’ threshold value)

QSM Standalone¶

How can we map the magnetic susceptibility sources from local field map?¶

A local field (or tissue field) is the field generated by the magnetic susceptibility sources from the brain tissue. The assumption of the field generated by a magnetic susceptibility point source is a dipole field. We can then try to compute the susceptibility values by deconvolving the local field with a dipole field. However, due to the fact that the dipole field in the k-space (i.e. the dipole kernel) has zero values on the cone surface at 54.7° and values close to this surface will also be close to zeros. As a result, the number of unknown parameters is more than the measurements we have in the data, leading to the problem ill-conditioned and the resulting QSM map has the so-called streaking artefact. To solve the ill-conditioned problem, most of the QSM algorithms try to add a regularisation term in the QMS dipole inversion problem, imposing spatial smoothness in the QSM map and/or the QSM map has similar anatomical features as shown in the magnitude images in order to reduce the streaking artefact in the QSM result.

Structure of the application¶

This standalone consists of two panels:

I/O panel, and
QSM panel.

The detailed description of each panel is given below:

I/O panel¶

Data input

This application accepts two types of data input method:
1. Specify a directory that contains all NIfTI images. Please specify the names of your data as in the following:
  - Local field data: must contain the string ‘local-field’ in the filename, e.g. local-field.nii.gz;
  - Header: must contain the string ‘header’ in the filename, e.g. header.mat;
  - Magnitude data: must contain the string ‘magn’ in the filename, e.g. magn.nii.gz;
  - (optional) Fieldmap standard deviation data: must contain the string ‘fieldmap-sd’ in the filename, e.g. fieldmap-sd.nii.gz;
  - (optional) Mask data: if provided, must contain string ‘mask’ in the filename, e.g. mask.nii.gz, or
2. Specify the required data separately using the GUI buttons.
Data output

You can specify the prefix of the data output name in the editable field ‘Output basename’. By default, the Sepia output will be stored in a directory named ‘output’ under the input directory, i.e. ‘/your/input/directory/output/’ with prefix ‘squirrel’. You can change the default output directory and basename to whatever you need. If the output directory does not exist, the application will create the directory.

QSM panel¶

Method:
1. TKD
Thresholded k-space division
1. Closed-form solution
closed-form solution with L2-norm regularisation
1. STI suite iLSQR
STI suite v3.0 implementation of iterative LSQR approach
1. iLSQR
2. FANSI
Fast algorithm for nonlinear susceptibility inversion
1. Star
STI suite v3.0 Star-QSM
1. MEDI
Morphology enabled dipole inversion (MEDI+0)
Output

The output of the standalone application is given below:
- QSM.nii.gz (quantitative susceptibility map, in ppm)

Choosing a right phase unwrapping method for your data¶

SEPIA provides a range of phase unwrapping methods from the supported toolboxes. In this demonstration, I will try to show the differences between two main phase unwrapping methods: (1) region growing and (2) Laplacian-based - an important processing step to reveal the (true) phase of MRI gradient echo data for QSM processing.

Data preparation¶

The in vivo test data was downloaded from Cornell MRI lab and converted to NIfTI format by dicm2nii. FSL’s BET was used to extract the brain mask to improve the processing result.

Processing¶

Test 1: Unwrapping the phase with a small amount of wrapping effect¶

To begin with, let have a look of the phase image from the 1st echo GRE data (TE=3.6ms)

Under a linear model assumption the phase of each voxel over time can be expressed by the following equation:

$phase = angular frequency \times time [1]$

However, this phase value is constrained in the interval of [-pi,pi) in the MRI phase data. Any value that is outside of this interval will be wrapped back to this interval (see red arrows). For QSM processing, it is important to unwrap the phase such that the true frequency shift inside the voxel can be computed using the above equation.

To unwrap the phase images, there are two main methods being used: (1) region growing and (2) Laplacian operator. In SEPIA, you can use the following commands to access the methods:

% region growing
unwrappedPhase = UnwrapPhaseMacro(wrappedPhase,matrixSize,voxelSize,'method','rg','magn',magn,'mask',mask);
% Laplacian
unwrappedPhase = UnwrapPhaseMacro(wrappedPhase,matrixSize,voxelSize,'method','laplacian');
unwrappedPhase = UnwrapPhaseMacro(wrappedPhase,matrixSize,voxelSize,'method','laplacian_stisuite');