Merge pull request #78 from oesteban/enh/NativeSpace

[ENH] Compute metrics in their appropriate image, use original image whenever possible

build/circle_tests.sh

@@ -9,13 +9,14 @@ declare -a MRIQC_TESTS=(
 # -i data/ds003_sub-01/ -o outputs/ss-all/out -w outputs/ss-all/work
 )
 
+if [[ $CIRCLE_NODE_INDEX -eq 0 ]]; then
+	docker run -i -v /etc/localtime:/etc/localtime:ro --entrypoint="/usr/bin/run_tests" oesteban/mriqc
+fi
+
 i=0
 for test_params in "${MRIQC_TESTS[@]}"; do
-	if [[ $i -eq 0 ]]; then
-		docker run -i -e TZ=PST --entrypoint="/usr/bin/run_tests" oesteban/mriqc
-	fi
 	if [ $(($i % $CIRCLE_NODE_TOTAL)) -eq $CIRCLE_NODE_INDEX ]; then
-	    docker run -i -e TZ=PST -v ~/scratch:/scratch -w /scratch oesteban/mriqc ${test_params}
+	    docker run -i -v /etc/localtime:/etc/localtime:ro -v ~/scratch:/scratch -w /scratch oesteban/mriqc ${test_params}
 	fi
 	((i=i+1))
 done

mriqc/__init__.py

@@ -31,7 +31,7 @@
 
 """
 
-__version__ = '0.4.0a1'
+__version__ = '0.5.0a2'
 __author__ = 'Oscar Esteban'
 __email__ = 'code@oscaresteban.es'
 __maintainer__ = 'Oscar Esteban'

mriqc/interfaces/qc.py

@@ -8,7 +8,7 @@
 # @Date:   2016-01-05 11:29:40
 # @Email:  code@oscaresteban.es
 # @Last modified by:   oesteban
-# @Last Modified time: 2016-04-13 15:35:14
+# @Last Modified time: 2016-04-15 15:39:39
 """ Nipype interfaces to quality control measures """
 
 import numpy as np
@@ -23,7 +23,8 @@
 IFLOGGER = logging.getLogger('interface')
 
 class StructuralQCInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc='File to be plotted')
+    in_file = File(exists=True, mandatory=True, desc='file to be plotted')
+    in_noinu = File(exists=True, mandatory=True, desc='image after INU correction')
     in_segm = File(exists=True, mandatory=True, desc='segmentation file from FSL FAST')
     in_bias = File(exists=True, mandatory=True, desc='bias file')
     air_msk = File(exists=True, mandatory=True, desc='air mask')
@@ -76,6 +77,10 @@ def _run_interface(self, runtime):
         # Remove negative values
         imdata[imdata < 0] = 0
 
+        # Load image corrected for INU
+        inudata = np.nan_to_num(nb.load(self.inputs.in_noinu).get_data())
+        inudata[inudata < 0] = 0
+
         segnii = nb.load(self.inputs.in_segm)
         segdata = segnii.get_data().astype(np.uint8)
 
@@ -86,25 +91,25 @@ def _run_interface(self, runtime):
         snrvals = []
         self._results['snr'] = {}
         for tlabel in ['csf', 'wm', 'gm']:
-            snrvals.append(snr(imdata, segdata, airdata, fglabel=tlabel))
+            snrvals.append(snr(inudata, segdata, airdata, fglabel=tlabel))
             self._results['snr'][tlabel] = snrvals[-1]
         self._results['snr']['total'] = np.mean(snrvals)
 
         # CNR
-        self._results['cnr'] = cnr(imdata, segdata)
+        self._results['cnr'] = cnr(inudata, segdata)
 
         # FBER
-        self._results['fber'] = fber(imdata, segdata, airdata)
+        self._results['fber'] = fber(inudata, segdata, airdata)
 
         # EFC
-        self._results['efc'] = efc(imdata)
+        self._results['efc'] = efc(inudata)
 
         # Artifacts
         self._results['qi1'] = art_qi1(airdata, artdata)
         self._results['qi2'] = art_qi2(imdata, airdata, artdata)
 
         # CJV
-        self._results['cjv'] = cjv(imdata, segdata)
+        self._results['cjv'] = cjv(inudata, segdata)
 
         pvmdata = []
         for fname in self.inputs.in_pvms:

mriqc/workflows/anatomical.py

@@ -7,7 +7,7 @@
 # @Date:   2016-01-05 11:24:05
 # @Email:  code@oscaresteban.es
 # @Last modified by:   oesteban
-# @Last Modified time: 2016-04-12 10:58:14
+# @Last Modified time: 2016-04-15 15:33:31
 """ A QC workflow for anatomical MRI """
 import os.path as op
 from nipype.pipeline import engine as pe
@@ -64,12 +64,12 @@ def anat_qc_workflow(name='aMRIQC', settings=None, sub_list=None):
         output_names=['out_qc'], function=_merge_dicts), name='merge_qc')
 
     arw = mri_reorient_wf()  # 1. Reorient anatomical image
-    n4itk = pe.Node(ants.N4BiasFieldCorrection(dimension=3, save_bias=True,
-                    bspline_fitting_distance=30.), name='Bias')
+    n4itk = pe.Node(ants.N4BiasFieldCorrection(dimension=3, save_bias=True),
+                    name='Bias')
     asw = skullstrip_wf()    # 2. Skull-stripping (afni)
     mask = pe.Node(fsl.ApplyMask(), name='MaskAnatomical')
     hmsk = headmsk_wf()
-    amw = airmsk_wf(save_memory=settings['save_memory'])
+    amw = airmsk_wf(save_memory=settings.get('save_memory', False))
 
     # Brain tissue segmentation
     segment = pe.Node(fsl.FAST(
@@ -87,17 +87,19 @@ def anat_qc_workflow(name='aMRIQC', settings=None, sub_list=None):
         (datasource, arw, [('anatomical_scan', 'inputnode.in_file')]),
         (arw, asw, [('outputnode.out_file', 'inputnode.in_file')]),
         (arw, n4itk, [('outputnode.out_file', 'input_image')]),
-        (asw, n4itk, [('outputnode.out_mask', 'weight_image')]),
+        # (asw, n4itk, [('outputnode.out_mask', 'mask_image')]),
         (n4itk, mask, [('output_image', 'in_file')]),
         (asw, mask, [('outputnode.out_mask', 'mask_file')]),
         (mask, segment, [('out_file', 'in_files')]),
-        (n4itk, hmsk, [('output_image', 'inputnode.in_file')]),
+        (arw, hmsk, [('outputnode.out_file', 'inputnode.in_file')]),
         (segment, hmsk, [('tissue_class_map', 'inputnode.in_segm')]),
+        (n4itk, measures, [('output_image', 'in_noinu')]),
         (arw, measures, [('outputnode.out_file', 'in_file')]),
-        (n4itk, fwhm, [('output_image', 'in_file')]),
+        (arw, fwhm, [('outputnode.out_file', 'in_file')]),
         (asw, fwhm, [('outputnode.out_mask', 'mask')]),
 
-        (n4itk, amw, [('output_image', 'inputnode.in_file')]),
+        (arw, amw, [('outputnode.out_file', 'inputnode.in_file')]),
+        (n4itk, amw, [('output_image', 'inputnode.in_noinu')]),
         (asw, amw, [('outputnode.out_mask', 'inputnode.in_mask')]),
         (hmsk, amw, [('outputnode.out_file', 'inputnode.head_mask')]),
 
@@ -234,7 +236,7 @@ def airmsk_wf(name='AirMaskWorkflow', save_memory=False):
     workflow = pe.Workflow(name=name)
 
     inputnode = pe.Node(niu.IdentityInterface(
-        fields=['in_file', 'in_mask', 'head_mask']), name='inputnode')
+        fields=['in_file', 'in_noinu', 'in_mask', 'head_mask']), name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'artifact_msk']),
                          name='outputnode')
 
@@ -244,23 +246,43 @@ def _invt_flags(transforms):
     # Spatial normalization, using ANTs
     norm = pe.Node(ants.Registration(dimension=3), name='normalize')
     norm.inputs.initial_moving_transform_com = 1
+    norm.inputs.winsorize_lower_quantile = 0.05
+    norm.inputs.winsorize_upper_quantile = 0.98
+    norm.inputs.float = True
+
     norm.inputs.transforms = ['Rigid', 'Affine']
     norm.inputs.transform_parameters = [(2.0,), (1.0,)]
-    norm.inputs.number_of_iterations = [[500], [1000, 200]]
+    norm.inputs.number_of_iterations = [[500], [200]]
     norm.inputs.convergence_window_size = [50, 20]
     norm.inputs.metric = ['Mattes', 'GC']
-    norm.inputs.metric_weight = [1] * 2
+    norm.inputs.metric_weight = [1] * 3
     norm.inputs.radius_or_number_of_bins = [64, 3]
     norm.inputs.sampling_strategy = ['Random', None]
     norm.inputs.sampling_percentage = [0.2, 1.]
-    norm.inputs.convergence_threshold = [1.e-8, 1.e-9]
-    norm.inputs.smoothing_sigmas = [[8], [4, 2]]
+    norm.inputs.smoothing_sigmas = [[8], [4]]
+    norm.inputs.shrink_factors = [[3], [2]]
+    norm.inputs.convergence_threshold = [1.e-8] * 2
     norm.inputs.sigma_units = ['mm'] * 2
-    norm.inputs.shrink_factors = [[3], [2, 1]]
     norm.inputs.use_estimate_learning_rate_once = [True] * 2
     norm.inputs.use_histogram_matching = [True] * 2
-    norm.inputs.winsorize_lower_quantile = 0.001
-    norm.inputs.winsorize_upper_quantile = 0.999
+
+#    norm.inputs.transforms = ['Rigid', 'Affine', 'SyN']
+#    norm.inputs.transform_parameters = [(2.0,), (1.0,), (.2, 3, 0)]
+#    norm.inputs.number_of_iterations = [[500], [200], [100]]
+#    norm.inputs.convergence_window_size = [50, 20, 10]
+#    norm.inputs.metric = ['Mattes', 'GC', 'Mattes']
+#    norm.inputs.metric_weight = [1] * 3
+#    norm.inputs.radius_or_number_of_bins = [64, 3, 64]
+#    norm.inputs.sampling_strategy = ['Random', None, 'Random']
+#    norm.inputs.sampling_percentage = [0.2, 1., 0.1]
+#    norm.inputs.convergence_threshold = [1.e-8] * 3
+#    norm.inputs.smoothing_sigmas = [[8], [4], [2]]
+#    norm.inputs.sigma_units = ['mm'] * 3
+#    norm.inputs.shrink_factors = [[3], [2], [2]]
+#    norm.inputs.use_estimate_learning_rate_once = [True] * 3
+#    norm.inputs.use_histogram_matching = [True] * 3
+
+
     if save_memory:
         norm.inputs.fixed_image = op.join(get_mni_template(), 'MNI152_T1_2mm.nii.gz')
         norm.inputs.fixed_image_mask = op.join(get_mni_template(),
@@ -283,7 +305,7 @@ def _invt_flags(transforms):
 
     workflow.connect([
         (inputnode, qi1, [('in_file', 'in_file')]),
-        (inputnode, norm, [('in_file', 'moving_image'),
+        (inputnode, norm, [('in_noinu', 'moving_image'),
                            ('in_mask', 'moving_image_mask')]),
         (norm, invt, [('forward_transforms', 'transforms'),
                       (('forward_transforms', _invt_flags), 'invert_transform_flags')]),