Definition and general description of the DPP Suite and modules
The process management system named Diffusion/Perfusion Project (DPP) described in the paper is an in-house developed suite written in Matlab (The MathWorks, Natick, MA, USA) and requiring the Image Processing Toolbox. The DPP is a collection of software modules all related to MRI data management, sharing the same GUI (Graphical User Interface) and exchanging data with each other. DPP Suite integrates commonly used software tools which are able to perform different types of data analysis and management. The primary aim of DPP Suite is to create a uniform environment where it can be possible to assess a considerably larger number of data from MS patients compared to the other current analysis methods, reducing as much as possible analysis complexity, time required and potential human errors. The most relevant procedures manageable through DPP Suite are schematically presented in Fig. 1, where data and operations are shown as a flowchart. In summary, DPP Suite operations include: 1) image pre-processing; 2) registration of T1-weighted images; 3) lesion segmentation and classification; 4) registration of lesion masks; 5) normal appearing tissue segmentation; 6) PWI and DWI map generation; 7) registration of tissue segmentation and quantitative MRI maps; 8) data analysis.
Pre-processing
Image pre-processing starts with the automatic anonymization of MRI sequences for each patient which are exported from a PACS database and copied in a local repository dynamically linked with the DPP Suite. This process aims at protecting patient identity during the whole procedure. Furthermore, demographic data, resulting analysis and patient identity are stored as a protected file controlled by the DPP Suite. The anonymization procedure is performed by Image Processing Toolbox Matlab function with some adaptations. The original embedded task for data anonymization, called dicomanon, was modified generating a function called dpp_dicomanon that can strengthen the patient data confidentiality level if needed. After anonymization, DICOM format [27] files are converted in NIFTI format [28] using the Statistical Parametric Mapping (SPM) 8 toolbox (http://www.fil.ion.ucl.ac.uk/spm/) [29]. Notably, pre-processing is performed on a limited dataset that includes the following MRI sequences: axial fluid-attenuated inversion recovery (FLAIR); two-dimensional (2D) axial non-enhanced T1-weighted spin-echo or three-dimensional (3D) non-enhanced T1 gradient-echo; axial Gadolinium (Gd)-enhanced T1-weighted spin-echo; axial DWI and axial PWI. Therefore, before data conversion, a DPP function selects patients having all MRI sequences requested for analysis, excluding those with incomplete radiological data. Concurrently, a set of report files containing a checklist of unavailable sequences for each non-conforming patient is generated. In addition, all MRI data are reorganized in a storage structure aimed at making a uniform data format, including sequences names, filenames or data storing and other details, that is independent of that given by different MRI scanners (e.g. Philips, Siemens or GE), without modifications in information content.
Registration of T1-weighted images
In this first step (Fig. 2), both axial 2D non-enhanced T1-weighted spin-echo or 3D non-enhanced T1 gradient-echo and axial Gd-enhanced T1-weighted spin-echo are registered in FLAIR space using an automated process based on the employment of a DPP module invoking one of the two following external software packages: SPM or FMRIB’s Linear Image Registration Tool (FLIRT) from FMRIB Software Library (FSL) suite (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/flirt) [30, 31]. Although T2-weighted images are generally considered the most sensitive in detecting infratentorial lesions [32], FLAIR images serve for the following lesion segmentation since they improve lesion/brain contrast due to the suppression of cerebrospinal fluid signal [33]. On the other hand, T1-weighted images allow the identification of MS T1 hypointense (black holes) and T1 Gd-enhancing lesions as well as the evaluation of brain atrophy [17]. In this regard, it is well-accepted that 3D are superior to 2D T1 images in the determination of brain volume [34]. However, we chose to process also 2D T1 images since, in clinical practice, not all centers are equipped to routinely perform 3D T1 sequences in MS patients. The choice between SPM and FSL FLIRT is user selectable.
Lesion segmentation and classification
The development of automated techniques for lesion detection is one of the most central challenges in MS research. Therefore, a number of approaches have recently been proposed [21–23]. However, it is generally accepted that there are no automatic lesion segmentation methods with 100 % reliability [18, 19]. This is the reason why a series of semi-automated algorithms were designed and tested in order to support and simplify lesion segmentation and classification. This step is one of the most critical and complex because it involves both human operation and automated processes. In fact, two basic operations are implemented: a) semi-automatic lesion segmentation; b) semi-automatic lesion classification.
Semi-automatic lesion segmentation
As depicted in Fig. 3, this step is performed on FLAIR images using Jim (Jim 6.0, Xinapse Systems, Leicester, UK; http://www.xinapse.com) [34] as an external software package. After visual identification, each lesion is semi-automatically segmented by a local threshold-based technique. In addition, a region of interest (ROI) of the NAWM is manually outlined. The output of this process is an ROI file corresponding to the delineated areas.
Semi-automatic lesion classification
Using the FLAIR-derived lesion ROI file, DDP Suite automatically masks the non-enhanced and Gd-enhanced T1-weighted images (Fig. 4). The DPP Suite then calculates the intensity of each lesion in both non-enhanced and Gd-enhanced T1-weighted sequence. Based on the recently proposed classification [20], the algorithm categorizes each T2-weighted hyperintense lesion, based on its intensity in four classes: 1) Gd-enhancing and T1-weighted isointense (C1); 2) Gd-enhancing and T1-weighted hypointense (C2); 3) non Gd-enhancing and T1-weighted isointense (C3); 4) non Gd-enhancing and T1-weighted hypointense, i.e. black holes (C4). As, in absence of serial MRI examinations, lesion activity is demonstrable only by Gd enhancement [35], C1 and C2 type lesions are considered as acute, whereas C3 and C4 lesions are judged as chronic [36, 37]. The output of the process is a recombined ROI file enriched with classification of each lesion given by a different color code for each lesion type and a descriptive text label where the type and characteristics of the single lesion are briefly reported (Fig. 4). Lesion classification is accomplished using a DPP algorithm based on a lesion intensity comparison to the NAWM ROI intensity level. The thresholds to identify the lesion classes are defined as a configurable multiplier of the standard deviation of the intensity values in the NAWM ROI.
$$ Thr{e}_{High}={\overline{I}}_{default}+Mf\cdot St-dev\left({I}_{NAWM}\right) $$
(1)
$$ Thr{e}_{Low}={\overline{I}}_{default}-Mf\cdot St-dev\left({I}_{NAWM}\right) $$
(2)
Where Ī is the mean intensity of the NAWM and Mf is the multiplier factor and St_dev is the standard deviation calculation function. The mean intensity of ach lesion is then calculated and compared with the two thresholds referred to in formulas (1) and (2) that are calculated for both non-enhanced and Gd-enhanced T1-weighted images. The automatic lesion classification algorithm performance is briefly presented in Fig. 5 to document the accuracy of the procedure. The error percentage of automatic classification was less than 1 % for C1 and C2 lesion classes and about 21 % and 23 % for C3 and C4 lesion classes, respectively. It is important to underline that the low error percentage found for C1 and C2 lesion classes could be related to the small number of this type of lesions occurred in the selected patients, as well as in all patients analyzed. Therefore, these results support the need of a visual correction for automatic lesion classification. The output of the automatic lesion classification process is a text file in Jim ROI Analysis Tool format, containing all the ROI definitions and classifications. Automatic lesion classifications are visually checked and revised by the operator. The new checked file is then automatically reloaded and converted by DDP Suite into different revised ROI files, including the lesions categorized as total lesions, lesion classes (C1, C2, C3 and C4) and single lesion, which are then used in the following steps.
Registration of lesion masks
This stage requires the creation of a lesion mask in the T1 space and is performed by SPM or FSL FLIRT external software packages invoked by a DPP module. Briefly, FLAIR images previously segmented are used to automatically produce a lesion mask that is coregistered in the non-enhanced T1-weighted space to obtain a registered lesion mask. When FLIRT is chosen, coregistration is done applying to the lesion mask in the FLAIR space the transformation matrix generated by the registration of FLAIR and T1. If SPM is selected, T1 is set as the reference image, FLAIR as source and the lesion mask as other. The result is a black and white lesion mask that is a binarization of the brain lesion area where lesions are white and the remaining brain tissue is black (Fig. 6).
Normal appearing tissue segmentation
As depicted in Fig. 7, DPP Suite invokes automatic segmentation of the basal ganglia and thalami via the external FSL tool FMRIB’s Integrated Registration and Segmentation Tool (FIRST) (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/first) [38]. The input is the non-enhanced T1-weighted sequence. The segmented basal ganglia structures include the right and left caudate, putamen and globus pallidus. The following step is represented by the automatic segmentation of grey and white matter for which DPP Suite invokes FSL SIENAX (Structural Image Evaluation using Normalization of Atrophy) (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/SIENA) [39]. In this case, the inputs are non-enhanced T1-weighted sequence and registered T1 lesion mask, whereas the output are two files including NAGM and NAWM images obtained after subtracting registered T1 lesion mask from non-enhanced T1-weighted. In addition, FSL SIENAX generates a report file including total brain tissue volume, as a whole and normalized according to skull size, and normalized NAGM and NAWM volumes.
PWI map generation
PWI studies are performed with a dynamic susceptibility contrast (DSC) MRI first-pass bolus-tracking technique using echo-planar gradient-echo T2* sequences. PWI analysis provides relative measurements of brain hemodynamic parameters such as Cerebral Blood Flow (CBF), Cerebral Blood Volume (CBV) and Mean Transit Time (MTT). CBF, CBV and MTT maps are generated by a singular value decomposition (SVD) deconvolution operation based on the measurement of an arterial input function (AIF) [40]. The calculation of PWI maps is automatically performed by Jim Perfusion Analysis tool invoked by DPP Suite. DSC images consist of a time-series of sequences, one volume for each time point, which monitor the concentration of an injected paramagnetic contrast agent transiting from the blood vessels to the brain tissue. Therefore, AIF determination implies the knowledge of the exact volume where the contrast agent is perceptible in the brain (Contrast Arrival Point). For that reason, the DPP Suite implements an algorithm that identifies the contrast arrival volume number using Jim Stats. This algorithm evaluates the mean intensity in volumes, using in subsequent steps the average intensity of the previous analyzed volumes, in order to reduce the noise affecting the images. It detects when the difference of mean intensity between two subsequent volumes is more than an adjustable threshold. The threshold is not an absolute value or an absolute percent, but is a fraction of the difference between maximum and minimum intensity value measured in all volumes of the sequence. Thus, the threshold is less dependent from the volume intensity absolute level.
DWI map generation
DWI studies are often performed with a single-shot echo-planar T2 spin-echo sequence according to the Stejskal-Tanner method [41]. The diffusion gradients are applied in three orthogonal directions (x, y, z) with two b-values (0 and 1000 s/mm2) to form the isotropic DWI images at b 1000 s/mm2. DWI analysis takes as input T2 images and extracts from these images apparent Diffusion Coefficient (ADC) maps, related to each Cartesian axis. More precisely, as reported elsewhere [42], ADC maps are generated using T2-weighted images at b 0s/mm2 and isotropic DWI images at b 1000 s/mm2 obtained in all three orthogonal directions. All calculations are automatically performed by Jim Image Algebra tool invoked by DPP Suite controlling the entire process. The process output is an average ADC map, obtained calculating the mean of the previously calculated ADC on the three Cartesian axes. After this process, previously generated ADC maps and PWI maps are ready to be registered in the next step (Fig. 8).
Registration of tissue segmentation and quantitative MRI maps
This step is aimed at automatically translating the results of all intermediate processes previously described in the FLAIR space to allow the final data analysis in which DWI and PWI values are measured in every lesion as well as in the basal ganglia and thalami, NAGM and NAWM. DPP Suite invokes SPM or FSL FLIRT external software tools which generate the following registered maps:
From T1 space:
-
1.
Basal ganglia and thalami maps;
-
2.
NAGM map;
-
3.
NAWM map;
from DWI space:
-
4.
Average ADC map;
from PWI space:
-
5.
CBV, CBF and MTT maps.
The different maps are registered into the FLAIR space and used in the next data analysis (Fig. 9). The coregistration is performed using a similar method as in registration of lesion masks.
Data analysis
The final step represents the key advance performed by the DPP Suite compared to existing semi-automated methods. Measurements are obtained by the integration of the entire set of images and maps previously generated. This result is obtained using distinct software modules. In this stage, two groups of automatic analysis are developed. In the first one (Fig. 10, panel a), NAGM, NAWM and basal ganglia maps are used to mask ADC and PWI maps owing to Jim Masker invoked by DPP Suite. At this point, DPP Suite activates the Jim Stats tool that is able to calculate ADC, CBF, CBV and MTT values in these brain areas. In addition, the Jim Stats tool also measures NAGM, NAWM and basal ganglia/thalami volumes and voxel number and DPP Suite parses previously saved FSL SIENAX report file including whole and normalized total brain tissue volume and normalized NAGM and NAWM volumes. In the second group of data analysis (Fig. 10, panel b), DPP Suite invokes Jim Masker to mask DWI and PWI maps with FLAIR segmented lesions which DPP Suite has previously saved and stored as total, single class (C1, C2, C3 and C4) and single lesion maps. In this way, it is possible to obtain ADC, CBF, CBV and MTT values in all types of lesions where, in addition, also volume and voxel number are measured. All of these data are stored in a comma separated values (csv) file.
Graphical User Interface
The Graphical user Interface (GUI) is shown in Fig. 11. The GUI keeps hidden all configuration parameters. Thus, the user is unable to modify the system parametrization and has a limited number of operations to perform. The configuration parameters are collected in a text configuration file. However, the configuration can be modified as needed by an expert operator, based on the needs of the study. These parameters are not shown in the GUI to facilitate the use of DPP suite by non-expert operators. The GUI configuration presents two sections: the first area in which all the process steps are singularly selectable and the second area where only three process macro steps can be activated. All these steps are performed automatically by DPP Suite, except lesion segmentation and classification which requires human intervention. The GUI presents also two directory browsing selection buttons for a more user friendly and flexible selection of input and output data directories. In order to control the process status and evolution, three status bars are also present in the lower part of the GUI window: a) a status bar defining the step under execution or the process execution status; b) a status-sub bar describing the sub-module under execution; c) a status-info bar indicating the patient number whose data are under analysis. The lower button is the step/process activation button, which starts the selected operations to be executed in an automated way by the DPP Suite.
Process safety and data security
All process steps, such as utilized parameters, the final results and the intermediate data are logged, saved and stored in the DPP Suite. This ensures complete environment preservation. For example, when an external tool is invoked by the DPP Suite both images and results and the tool’s standard output/error data log is read and saved. This approach allows post processing visual verification of e.g. co-registration (SPM or FSL FLIRT), tissue segmentation (SIENAX or FSL FIRST) or automatic detection of AIF-like voxels (Jim Perfusion tool). In Fig. 12 is presented an example of coregistration check performed using the Check Reg tool part of SPM module for functional MRI. Although processing parameters were tuned and optimized for every single step in DPP suite, in case of failure, configuration can be modified to correct processing errors. In this way, a clinician can check the results and make sure that they are accurate. Other security measures are integrated in DPP suite as well. For example, if a set of data represent both an output for the (n-1)th step and the input for the (n)th step, the entire data are copied before being modified by the (n)th step. Moreover, execution times and software version are preserved for every run of the suite (e.g. pre-processing data structure). This approach allows data reconstruction, process verification and software debugging, making it possible to reproduce every single step or the entire process and, thus, future validation of the DPP Suite.
