Merge pull request #121 from edickie/enh_more_tests

[ENC] various added tests and bug fixes

ciftify/bin/ciftify_PINT_vertices.py

@@ -263,7 +263,7 @@ def roi_surf_data(df, vertex_colname, surf, hemisphere, roi_radius):
     with ciftify.utils.TempDir() as lil_tmpdir:
         ## write a temp vertex list file
         vertex_list = os.path.join(lil_tmpdir, 'vertex_list.txt')
-        df.loc[df.hemi == hemisphere, vertex_colname].to_csv(vertex_list,sep='\n',index=False)
+        df.loc[df.hemi == hemisphere, vertex_colname].to_csv(vertex_list,sep='\n',index=False, header=False)
 
         ## from the temp text build - func masks and target masks
         roi_surf = os.path.join(lil_tmpdir,'roi_surf.func.gii')
@@ -329,7 +329,7 @@ def linalg_calc_residulals(X, Y):
     ----------------
     Residuals as 1D array
     '''
-    betas = np.linalg.lstsq(X, Y)[0]
+    betas = np.linalg.lstsq(X, Y, rcond = -1)[0]
     residuals = Y - X.dot(betas)
     return(residuals)
 
@@ -407,7 +407,7 @@ def pint_move_vertex(df, idx, vertex_incol, vertex_outcol,
             o_networks = set(netmeants.columns.tolist()) - set([network])
             seed_corrs[idx_mask] = mass_partial_corr(meants,
                                       func_data[idx_mask, :],
-                                      netmeants.loc[:,o_networks].as_matrix())
+                                      netmeants.loc[:,o_networks].values)
         else:
             seed_corrs[idx_mask] = np.corrcoef(meants,
                                                   func_data[idx_mask, :])[0, 1:]

ciftify/bin/ciftify_dlabel_to_vol

@@ -0,0 +1 @@
+ciftify_dlabel_to_vol.py

ciftify/bin/ciftify_dlabel_to_vol.py

@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+"""
+Takes a atlas or cluster map ('dlabel.nii') and output a nifti image in a particular participants space
+
+Usage:
+    ciftify_dlabel_to_vol [options] --input-dlabel <input.dlabel.nii> --left-mid-surface <input.L.surf.gii> --volume-template <voltemplate.nii.gz> --output-nifti <output.nii.gz>
+
+Arguments:
+    --input-dlabel <input.dlabel.nii>       Input atlas or cluster map.
+    --left-mid-surface <input.L.surf.gii>   Left midthickness surface file (note the other 5 needed surfaces are inferred from this)
+    --volume-template <voltemplate.nii.gz>  Volume that defines the output space and resolution
+    --output-nifti <output.nii.gz>          Output nifti atlas
+
+Options:
+    --map-number INT         The map number [default: 1] within the dlabel file to report on.
+    --use-nearest-vertex MM  Use wb_commands nearest vertex method (instead 
+                             of ribbon contrained) where voxels nearest 
+                             to the midthickness surface are labeled 
+                             within the specified MM distance. 
+
+    --debug                Debug logging
+    -n,--dry-run           Dry run
+    -h, --help             Prints this message
+
+DETAILS:
+
+Useful for building regions of interest in volume space for diffusion analysis, or as an intermediate for checking volume-to-surface mapping.
+
+"""
+from docopt import docopt
+import os, sys
+import ciftify.utils
+from ciftify.meants import NibInput
+import logging
+from ciftify.niio import cifti_info, voxel_spacing
+from nilearn.image import resample_to_img
+
+config_path = os.path.join(os.path.dirname(ciftify.config.find_ciftify_global()), 'bin', "logging.conf")
+logging.config.fileConfig(config_path, disable_existing_loggers=False)
+logger = logging.getLogger(os.path.basename(__file__))
+
+class UserSettings(object):
+    def __init__(self, arguments):
+        self.dlabel_in = NibInput(arguments['--input-dlabel'])
+        self.output_nifti = self.__get_output_path(arguments['--output-nifti'])
+        self.volume_in = NibInput(arguments['--volume-template'])
+        self.use_nearest = self.__get_surf_method(arguments['--use-nearest-vertex'])
+        self.surfs = self.__get_surfs(arguments['--left-mid-surface'])
+        self.map_number = str(arguments['--map-number'])
+        #returns a bool of wether or not there are volume space labels to deal with
+        self.dlabel_in.has_volume = cifti_info(self.dlabel_in.path)['maps_to_volume']
+
+    def __get_surfs(self, L_mid_path):
+
+        surfs = {
+            'L': {
+                'pial': L_mid_path.replace('midthickness', 'pial'),
+                'midthickness': L_mid_path,
+                'white': L_mid_path.replace('midthickness', 'white'),
+            },
+            'R' : {
+                'pial': L_mid_path.replace('midthickness', 'pial').replace('L.', 'R.').replace('hemi-L','hemi-R'),
+                'midthickness': L_mid_path.replace('L.', 'R.').replace('hemi-L','hemi-R'),
+                'white': L_mid_path.replace('midthickness', 'white').replace('L.', 'R.').replace('hemi-L','hemi-R'),
+            }
+        }
+        return surfs
+
+    def __get_output_path(self, output_arg):
+        ciftify.utils.check_output_writable(output_arg)
+        if not output_arg.endswith('.nii.gz'):
+            logger.warning('Recommended output extension is ".nii.gz", path {} given'.format(output_arg))
+        return output_arg
+
+    def __get_surf_method(self, nearest_arg):
+        if not nearest_arg:
+            return None
+        else:
+            return str(nearest_arg)
+
+def dlabel_number_maps(dlabel_file_path):
+    '''make sure that we get the correct settings'''
+    c_info = ciftify.utils.get_stdout(['wb_command', '-file-information',
+                dlabel_file_path, '-no-map-info'])
+    for line in c_info.split(os.linesep):
+        if "Number of Columns:" in line:
+            number_of_cols_str = line.replace("Number of Columns:",'')
+            number_of_cols_int = int(number_of_cols_str)
+    return number_of_cols_int
+
+def run_ciftify_dlabel_to_vol(arguments, tmpdir):
+
+    settings = UserSettings(arguments)
+
+    if dlabel_number_maps(settings.dlabel_in.path) > 1:
+        dlabel_file = os.path.join(tmpdir, 'split_in.dlabel.nii')
+        ciftify.utils.run(['wb_command', '-cifti-merge',
+            dlabel_file,
+            '-cifti', settings.dlabel_in.path,
+            '-column', settings.map_number])
+    else:
+        dlabel_file = settings.dlabel_in.path
+
+    ## split the files into two label.gii files then map them to the volume
+    ciftify.utils.run(['wb_command', '-cifti-separate',
+        dlabel_file, 'COLUMN',
+        '-label', 'CORTEX_LEFT', os.path.join(tmpdir, 'atlas.L.label.gii'),
+        '-label', 'CORTEX_RIGHT', os.path.join(tmpdir, 'atlas.R.label.gii')])
+
+    for hemi in ['L', 'R']:
+        if not settings.use_nearest:
+            ciftify.utils.run(['wb_command', '-label-to-volume-mapping',
+                os.path.join(tmpdir, 'atlas.{}.label.gii'.format(hemi)),
+                settings.surfs[hemi]['midthickness'],
+                settings.volume_in.path,
+                os.path.join(tmpdir, '{}.nii.gz'.format(hemi)),
+                '-ribbon-constrained',
+                settings.surfs[hemi]['white'],
+                settings.surfs[hemi]['pial']])
+        else:
+            ciftify.utils.run(['wb_command', '-label-to-volume-mapping',
+                os.path.join(tmpdir, 'atlas.{}.label.gii'.format(hemi)),
+                settings.surfs[hemi]['midthickness'],
+                settings.volume_in.path,
+                os.path.join(tmpdir, '{}.nii.gz'.format(hemi)),
+                '-nearest-vertex', settings.use_nearest])
+
+    # if there is a volume component - then separate that out to nii files
+    if not settings.dlabel_in.has_volume:
+        ciftify.utils.run(['wb_command', '-volume-math',
+            '"Lsurf + Rsurf"',
+            settings.output_nifti,
+            '-var Lsurf',os.path.join(tmpdir, 'L.nii.gz'),
+            '-var Rsurf', os.path.join(tmpdir, 'R.nii.gz')])
+    else:
+        subcort_first_labels = os.path.join(tmpdir, 'subcort1.nii.gz')
+        ciftify.utils.run(['wb_command', '-cifti-separate',
+            dlabel_file, 'COLUMN',
+            '-volume-all', subcort_first_labels])
+
+        ## add bit here that checks if the two match in voxel resolution
+        if ciftify.niio.voxel_spacing(subcort_first_labels) != ciftify.niio.voxel_spacing(settings.volume_in.path):
+        ## if not then resample
+            subcort_final_labels = os.path.join(tmpdir, 'subcort2.nii.gz')
+            resampled_ref_vol1 = resample_to_img(
+                source_img = subcort_first_labels,
+                target_img = settings.volume_in.path)
+            resampled_ref_vol1.to_filename(subcort_final_labels)
+        else:
+            subcort_final_labels = subcort_labels
+
+        ## then add all three together
+        ciftify.utils.run(['wb_command', '-volume-math',
+            '"Lsurf + Rsurf + subcort"',
+            settings.output_nifti,
+            '-var Lsurf', os.path.join(tmpdir, 'L.nii.gz'),
+            '-var Rsurf', os.path.join(tmpdir, 'R.nii.gz'),
+            '-var subcort', subcort_final_labels])
+
+    ## get the labels dictionary out of the dlabel files and write it into the outputs
+    ciftify.utils.run(['wb_command', '-cifti-label-export-table',
+        dlabel_file, settings.map_number,
+        os.path.join(tmpdir, 'label_table.txt')])
+
+    ciftify.utils.run(['wb_command', '-volume-label-import',
+       settings.output_nifti,
+       os.path.join(tmpdir, 'label_table.txt'),
+       settings.output_nifti])
+
+def main():
+    arguments = docopt(__doc__)
+
+    logger.setLevel(logging.WARNING)
+
+    if arguments['--debug']:
+        logger.setLevel(logging.DEBUG)
+        logging.getLogger('ciftify').setLevel(logging.DEBUG)
+
+    ## set up the top of the log
+    logger.info('{}{}'.format(ciftify.utils.ciftify_logo(),
+        ciftify.utils.section_header('Starting ciftify_dlabel_to_vol')))
+
+    ciftify.utils.log_arguments(arguments)
+
+    with ciftify.utils.TempDir() as tmpdir:
+        logger.info('Creating tempdir:{} on host:{}'.format(tmpdir,
+                    os.uname()[1]))
+        ret = run_ciftify_dlabel_to_vol(arguments, tmpdir)
+
+if __name__ == '__main__':
+    main()

ciftify/bin/ciftify_meants.py

@@ -99,7 +99,7 @@ def cifti_parcellate_to_meants(settings):
             tmp_parcelated = os.path.join(tempdir, 'parcellated.pscalar.nii')
         ciftify.utils.run(['wb_command', '-cifti-parcellate',
             settings.func.path, settings.seed.path,
-            'COLUMN', tmp_parcelated])
+            'COLUMN', tmp_parcelated, '-include-empty'])
         ciftify.utils.run(['wb_command', '-cifti-convert', '-to-text',
             tmp_parcelated, settings.outputcsv,'-col-delim ","'])
         if settings.outputlabels:

ciftify/bin/ciftify_seed_corr.py

@@ -93,7 +93,7 @@ def get_output_prefix(self, outputname):
     def get_outputcsv(self, output_ts):
         '''set outputcsv name if this is asked for'''
         if output_ts:
-            outputcsv = '_meants.csv'.format(self.output_prefix)
+            outputcsv = '{}_meants.csv'.format(self.output_prefix)
         else:
             outputcsv = None
         return(outputcsv)
@@ -153,7 +153,7 @@ def run_ciftify_seed_corr(settings, tempdir):
     seed_ts = seed_ts.reshape(seed_ts.shape[0]*seed_ts.shape[1])
     logger.debug('seed_ts shape after reshaping {}'.format(seed_ts.shape))
     logger.debug('Writing output with prefix: {}'.format(settings.output_prefix))
-
+    logger.debug('Writing meants: {}'.format(settings.outputcsv))
     logger.info('Using numpy to calculate seed-correlation')
 
     ## convert to nifti

ciftify/bin/ciftify_surface_rois.py

@@ -109,7 +109,7 @@ def run_ciftify_surface_rois(arguments, tmpdir):
         vertices = df.loc[df[hemi_col] == hemisphere, vertex_col]
         logger.info('{} vertices are: {}'.format(hemisphere, vertices))
         if len(vertices) > 0:
-            vertices.to_csv(vertex_list,sep='\n',index=False)
+            vertices.to_csv(vertex_list,sep='\n',index=False, header = False)
 
             if gaussian:
                 run(['wb_command', '-surface-geodesic-rois', surf,
@@ -133,7 +133,7 @@ def run_ciftify_surface_rois(arguments, tmpdir):
                       rois_2D, 'SUM', rois_1D])
 
         else:
-            pd.Series([1]).to_csv(vertex_list,sep='\n',index=False)
+            pd.Series([1]).to_csv(vertex_list,sep='\n',index=False, header = False)
             run(['wb_command', '-surface-geodesic-rois', surf,
                 str(radius),  vertex_list, rois_2D])
             run(['wb_command -metric-math "x*0"', rois_1D,

ciftify/meants.py

@@ -172,8 +172,7 @@ def calc_meants_with_numpy(settings, outputlabels = None):
     with ciftify.utils.TempDir() as tempdir:
         func_data, seed_data, mask_data = load_data_as_numpy_arrays(settings, tempdir)
 
-    std_array = np.std(func_data, axis=1)
-    mask_indices = np.where(std_array > 0)[0]
+    mask_indices = np.where(np.isfinite(func_data[:,0]))[0]
 
     if settings.mask:
         # attempt to mask out non-brain regions in ROIs

tests/data/PINT/pint_clean_sm0_sm8_pint.log

@@ -0,0 +1,48 @@
+
+'||'''|, |''||''| '||\   ||` |''||''|
+ ||   ||    ||     ||\  ||     ||
+ ||...|'    ||     || \ ||     ||
+ ||         ||     ||  \||     ||
+.||      |..||..| .||   \||.   .||.
+
+-------------------------------------------------------------
+2019-07-25 17:01:34.407253 : Personalized Intrinsic Network Topolography (PINT)
+-------------------------------------------------------------
+
+Arguments: 
+    functional data: /scratch/edickie/test_PINT_bug_2019/clean_img/output_clean_s0.dtseries.nii
+    left surface: /projects/edickie/code/ciftify/tests/data/sub-50005.L.midthickness.32k_fs_LR.surf.gii
+    right surface: /projects/edickie/code/ciftify/tests/data/sub-50005.R.midthickness.32k_fs_LR.surf.gii
+    pint template csv: /projects/edickie/code/ciftify/ciftify/data/PINT/Yeo7_2011_80verts.csv
+    output prefix: /scratch/edickie/test_PINT_bug_2019/ciftify_PINT_newenv_rcondup/testsub_clean_s0_s8
+Pre-smoothing FWHM: 8 (Sigma 3.3972872011520763)
+    Sampling ROI radius (mm): 6
+    Search ROI radius (mm): 6
+    Paddding ROI radius (mm): 12
+    Maximizing partial correlation
+
+-------------------------------------------------------------
+2019-07-25 17:01:34.408238 : Starting PINT
+-------------------------------------------------------------
+
+Iteration 0 	max distance: 5.994324207305908	Vertices Moved: 78
+Iteration 1 	max distance: 5.987587928771973	Vertices Moved: 69
+Iteration 2 	max distance: 5.98065185546875	Vertices Moved: 49
+Iteration 3 	max distance: 5.938047409057617	Vertices Moved: 40
+Iteration 4 	max distance: 5.992596626281738	Vertices Moved: 23
+Iteration 5 	max distance: 5.946619987487793	Vertices Moved: 22
+Iteration 6 	max distance: 5.197894096374512	Vertices Moved: 13
+Iteration 7 	max distance: 5.907351016998291	Vertices Moved: 12
+Iteration 8 	max distance: 5.334118843078613	Vertices Moved: 8
+Iteration 9 	max distance: 4.705900192260742	Vertices Moved: 9
+Iteration 10 	max distance: 4.10476541519165	Vertices Moved: 9
+Iteration 11 	max distance: 5.178483486175537	Vertices Moved: 8
+Iteration 12 	max distance: 5.8631439208984375	Vertices Moved: 8
+Iteration 13 	max distance: 2.76534366607666	Vertices Moved: 4
+Iteration 14 	max distance: 1.7977086305618286	Vertices Moved: 3
+Iteration 15 	max distance: 0.0	Vertices Moved: 0
+
+-------------------------------------------------------------
+2019-07-25 17:05:25.707749 : Done PINT
+-------------------------------------------------------------
+