Merge pull request #9 from pnlbwh/fix-warp

do not glob transform files with *ToMNI* pattern during --create, structure outputs into csv, provide demonstrative plots

Author: tashrifbillah

README.md

@@ -29,6 +29,8 @@ Table of Contents
       * [Create template](#create-template)
       * [Harmonize data](#harmonize-data)
       * [Debug](#debug)
+         * [1. Same target list](#1-same-target-list)
+         * [2. Different target list](#2-different-target-list)
    * [Tests](#tests)
       * [1. pipeline](#1-pipeline)
       * [2. unittest](#2-unittest)
@@ -339,29 +341,103 @@ or provided(--nshm 4), data can be harmonized.
 ## Harmonize data
 
     multi-shell-dMRIharmonization/lib/multi-shell-harmonization.py --tar_list tar_list.txt 
-    --ref_name REF --tar_name TAR --template template
+    --tar_name TAR --template template
     --process
 
     multi-shell-dMRIharmonization/lib/multi-shell-harmonization.py --tar_list tar_list.txt 
-    --ref_name REF --tar_name TAR --template template
+    --tar_name TAR --template template
     --process --debug
 
 
 ## Debug
 
-Runs together with `--create` and `--process`
+### 1. Same target list
+
+Run together with `--create` and `--process`:
 
     multi-shell-dMRIharmonization/lib/multi-shell-harmonization.py --ref_list ref_list.txt --tar_list tar_list.txt 
     --ref_name REF --tar_name TAR --template template
     --create --process --debug
 
 
+### 2. Different target list
+
+In theory, you want to create a template with small number of data from each sites and then use the template to 
+harmonize all of your data in the target site. Since the data for template creation and harmonization are not same, 
+we don't have luxury of using `--create --process --debug` altogether. Instead, you would use `--debug` with each of 
+`--create` and `--process`. The `--debug` flag creates some files that are used to obtain statistics later.
+The steps are described below:
+
+(i) Create template with `--debug` enabled:
+    
+    multi-shell-dMRIharmonization/lib/multi-shell-harmonization.py --ref_list ref_list.txt --tar_list tar_small_list.txt 
+    --ref_name REF --tar_name TAR --template template
+    --create --debug
+
+(ii) Harmonize data with `--debug` enabled:
+
+    multi-shell-dMRIharmonization/lib/multi-shell-harmonization.py --tar_list tar_list.txt 
+    --tar_name TAR --template template
+    --process --debug
+
+(iii) Obtain statistics
+    
+    ## reference site ##
+    
+    # start with highest bvalue shell
+    
+    multi-shell-dMRIharmonization/lib/tests/fa_skeleton_test.py 
+    -i ref_list_b3000.csv.modified -s REF 
+    -t template/ --bshell_b 3000    
+    
+    # repeat for other non-zero bvalues
+    
+    multi-shell-dMRIharmonization/lib/tests/fa_skeleton_test.py 
+    -i ref_list_b2000.csv.modified -s REF 
+    -t template/ --bshell_b 2000
+    
+    multi-shell-dMRIharmonization/lib/tests/fa_skeleton_test.py 
+    -i ref_list_b1000.csv.modified -s REF 
+    -t template/ --bshell_b 1000
+
+    ...
+    ...
+    
+    
+    
+    ## target site before harmonization ##
+    
+    # again, start with highest bvalue shell, notice the absence of ".modified" at the end of -i
+    
+    multi-shell-dMRIharmonization/lib/tests/fa_skeleton_test.py 
+    -i tar_list_b3000.csv -s TAR 
+    -t template/ --bshell_b 3000
+    
+    # repeat for other non-zero bvalues
+    
+    ...
+    ...
+    
+    
+    
+    ## target site after harmonization ##
+    
+    # once again, start with highest bvalue shell, notice the presence of ".modified.harmonized" at the end of -i
+    
+    multi-shell-dMRIharmonization/lib/tests/fa_skeleton_test.py 
+    -i tar_list_b3000.csv.modified.harmonized -s TAR 
+    -t template/ --bshell_b 3000
+    
+    # repeat for other non-zero bvalues
+    ...
+    ...
+
 
-**TBD** Tests are incomplete as of now but will be completed soon.
 
 # Tests
 
 A small test data is provided with each [release](https://github.com/pnlbwh/multi-shell-dMRIharmonization/releases). 
+**TBD** Tests are incomplete as of now but will be completed soon.
 
 
 ## 1. pipeline
@@ -442,7 +518,10 @@ each subject FA (reference, target before harmonization, and after harmonization
     optional arguments:
       -h, --help            show this help message and exit
       -i INPUT, --input INPUT
-                            input list of FA images
+                            a .txt/.csv file that you used in/obtained from
+                            harmonization.py having two columns for (img,mask)
+                            pair. See pnlbwh/dMRIharmonization documentation for
+                            more details
       -s SITE, --site SITE  site name for locating template FA and mask in
                             tempalte directory
       -t TEMPLATE, --template TEMPLATE
@@ -466,7 +545,7 @@ in the same space, registration is performed only once. The script is intelligen
 files if registration was performed before.
 
 In multi-shell-dMRIharmonization approach, registration is performed with highest b-shell. Obtained transform files 
-are used to warp rest of the b-shells.
+are used to warp rest of the b-shells. See [Different target list](#2-different-target-list) for details.
 
 
 # Caveats/Issues

lib/buildTemplate.py

@@ -47,8 +47,8 @@ def warp_bands(imgPath, maskPath, templatePath):
     prefix= basename(imgPath).split('.nii')[0]
     transPrefix= prefix.replace(f'_b{bshell_b}','')
     directory= dirname(imgPath)
-    warp = glob(pjoin(templatePath, transPrefix + f'*_FA*[!Inverse]Warp.nii.gz'))
-    trans = glob(pjoin(templatePath, transPrefix + f'*_FA*GenericAffine.mat'))
+    warp = glob(pjoin(templatePath, transPrefix + f'*_FA*[!ToMNI]1Warp.nii.gz'))
+    trans = glob(pjoin(templatePath, transPrefix + f'*_FA*[!ToMNI]0GenericAffine.mat'))
 
     # warping the mask
     applyXform(maskPath,

lib/harmonization.py

@@ -20,26 +20,11 @@
 from determineNshm import verifyNshmForAll, determineNshm
 from util import *
 from fileUtil import read_caselist, check_dir, check_csv
-from conversion import read_imgs_masks
 
 N_CPU= psutil.cpu_count()
 SCRIPTDIR= dirname(__file__)
 
 
-def printStat(ref_mean, csvFile):
-
-    print('mean FA over IIT_mean_FA_skeleton.nii.gz for all cases: ')
-    imgs, _ = read_imgs_masks(csvFile)
-    for i, imgPath in enumerate(imgs):
-        print(basename(imgPath), ref_mean[i])
-
-    print('')
-    print('mean meanFA: ', np.mean(ref_mean))
-    print('std meanFA: ', np.std(ref_mean))
-    print('')
-
-
-
 class pipeline(cli.Application):
 
     """
@@ -62,19 +47,19 @@ class pipeline(cli.Application):
     ref_csv = cli.SwitchAttr(
         ['--ref_list'],
         cli.ExistingFile,
-        help='reference csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\ndwi2,mask2\n...',
+        help='reference csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\\n dwi2,mask2\\n...',
         mandatory=False)
 
     target_csv = cli.SwitchAttr(
         ['--tar_list'],
         cli.ExistingFile,
-        help='target csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\ndwi2,mask2\n...',
+        help='target csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\\n dwi2,mask2\n...',
         mandatory=False)
 
     harm_csv = cli.SwitchAttr(
         ['--harm_list'],
         cli.ExistingFile,
-        help='harmonized csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\ndwi2,mask2\n...',
+        help='harmonized csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\\n dwi2,mask2\n...',
         mandatory=False)
 
     templatePath = cli.SwitchAttr(
@@ -145,7 +130,7 @@ class pipeline(cli.Application):
     reference = cli.SwitchAttr(
         '--ref_name',
         help= 'reference site name',
-        mandatory= True)
+        mandatory= False)
 
     target = cli.SwitchAttr(
         '--tar_name',
@@ -263,6 +248,7 @@ def createTemplate(self):
     def harmonizeData(self):
 
         from reconstSignal import reconst
+        from preprocess import dti_harm
 
         # check the templatePath
         if not exists(self.templatePath):
@@ -324,7 +310,7 @@ def harmonizeData(self):
         pool.close()
         pool.join()
 
-        
+       
         # loop for debugging
         # res= []
         # for imgPath, maskPath in zip(imgs, masks):
@@ -341,6 +327,17 @@ def harmonizeData(self):
         if preFlag:
             fm.close()
         fh.close()
+
+        
+        if self.debug:
+            harmImgs, harmMasks= read_caselist(self.harm_csv)
+            pool = multiprocessing.Pool(self.N_proc)
+            for imgPath,maskPath in zip(harmImgs,harmMasks):
+                pool.apply_async(func= dti_harm, args= (imgPath,maskPath))
+            pool.close()
+            pool.join()
+
+
         print('\n\nHarmonization completed\n\n')
 
 
@@ -365,42 +362,56 @@ def showStat(self):
 
         from debug_fa import analyzeStat
         from datetime import datetime
+        from harm_plot import generate_csv, harm_plot
+        import pandas as pd
 
-        print('Printing statistics :\n\n')
-        
-        # save statistics for future
-        statFile= pjoin(self.templatePath, 'meanFAstat.txt') 
-        f= open(statFile,'a')
-        stdout= sys.stdout
-        sys.stdout= f
-        
-        print(datetime.now().strftime('%c'),'\n')
-        
-        print('b-shell', self.bshell_b, '\n')
+        print('\n\nComputing statistics\n\n')
 
         print(f'{self.reference} site: ')
         ref_mean = analyzeStat(self.ref_csv, self.templatePath)
-        printStat(ref_mean, self.ref_csv)
+        generate_csv(self.ref_csv, ref_mean, pjoin(self.templatePath, self.reference), self.bshell_b)
 
         print(f'{self.target} site before harmonization: ')
-        target_mean_before = analyzeStat(self.tar_unproc_csv, self.templatePath)        
-        printStat(target_mean_before, self.tar_unproc_csv)
+        target_mean_before = analyzeStat(self.tar_unproc_csv, self.templatePath)
+        generate_csv(self.tar_unproc_csv, target_mean_before, pjoin(self.templatePath, self.target)+'_before', self.bshell_b)
 
         print(f'{self.target} site after harmonization: ')
-        target_mean_after = analyzeStat(self.harm_csv, self.templatePath) 
-        printStat(target_mean_after, self.harm_csv)
+        target_mean_after = analyzeStat(self.harm_csv, self.templatePath)
+        generate_csv(self.harm_csv, target_mean_after, pjoin(self.templatePath, self.target)+'_after', self.bshell_b)
+        
+        
+        print('\n\nPrinting statistics\n\n')
+        # save statistics for future
+        statFile= pjoin(self.templatePath, 'meanFAstat.csv')
+        if isfile(statFile):
+            df= pd.read_csv(statFile)
+        else:
+            timestamp= datetime.now().strftime('%m/%d/%y %H:%M')
+            sites= [f'{self.reference}',f'{self.target}_before',f'{self.target}_after']
+            df= pd.DataFrame({timestamp:sites})
+        
+        header= f'mean meanFA b{self.bshell_b}'
+        value= [np.mean(x) for x in [ref_mean, target_mean_before, target_mean_after]]
+        df= df.assign(**{header:value})
+        
+        # print an empty line so future results, if appended, are visually separate
+        # df=df.append(pd.Series(),ignore_index=True)
+        
+        df.to_csv(statFile, index=False)
 
-        print('\n\n')
-
-        f.close()
-        sys.stdout= stdout
-
         # print statistics on console        
         with open(statFile) as f:
-            print(f.read())
+            print(f.read())            
 
-        print('\nThe statistics are also saved in ', statFile)
-    
+        # generate graph
+        ebar= harm_plot([ref_mean, target_mean_before, target_mean_after],
+                         labels=[self.reference, self.target+'_before', self.target+'_after'],
+                         outPrefix=pjoin(self.templatePath,'meanFAstat'), bshell_b=self.bshell_b)
+
+        print(f'\nDetailed statistics, summary results, and demonstrative plots are saved in:\n\n{self.templatePath}/*_stat.csv'
+              f'\n{statFile}\n{ebar}\n')
+              
+   
 
 
     def sanityCheck(self):
@@ -447,9 +458,9 @@ def main(self):
         # determine N_shm in default mode during template creation
         if self.N_shm==-1 and self.create:
             if self.ref_csv:
-                ref_nshm_img = read_imgs_masks(self.ref_csv)[0][0]
+                ref_nshm_img = read_caselist(self.ref_csv)[0][0]
             elif self.target_csv:
-                ref_nshm_img = read_imgs_masks(self.target_csv)[0][0]
+                ref_nshm_img = read_caselist(self.target_csv)[0][0]
 
             directory= dirname(ref_nshm_img)
             prefix= basename(ref_nshm_img).split('.nii')[0]

lib/joinBshells.py

@@ -72,11 +72,7 @@ def joinBshells(imgPath, ref_bvals_file=None, ref_bvals=None, sep_prefix=None):
 
             joinedDwi[:,:,:,ind] = b0_bshell[:,:,:,1:]
 
-    if not isfile(harmPrefix + '.nii.gz'):
-        save_nifti(harmPrefix + '.nii.gz', joinedDwi, b0Img.affine, b0Img.header)
-    else:
-        print(harmPrefix + '.nii.gz', 'already exists, not overwritten.')
-
+    save_nifti(harmPrefix + '.nii.gz', joinedDwi, b0Img.affine, b0Img.header)
 
 
 def joinAllBshells(tar_csv, ref_bvals_file, separatedPrefix=None, ncpu=4):

lib/multi-shell-harmonization.py

@@ -54,7 +54,6 @@ class multi_shell_pipeline(cli.Application):
         help='target csv/txt file with first column for dwi and 2nd column for mask: dwi1,mask1\\ndwi2,mask2\\n...',
         mandatory=False)
 
-
     templatePath = cli.SwitchAttr(
         ['--template'],
         help='template directory',
@@ -107,7 +106,7 @@ class multi_shell_pipeline(cli.Application):
     reference = cli.SwitchAttr(
         '--ref_name',
         help= 'reference site name',
-        mandatory= True)
+        mandatory= False)
 
     target = cli.SwitchAttr(
         '--tar_name',
@@ -122,6 +121,9 @@ class multi_shell_pipeline(cli.Application):
 
     def main(self):
 
+        # FIXME https://github.com/pnlbwh/multi-shell-dMRIharmonization/issues/12
+        self.force=False
+
         if self.N_proc=='-1':
             self.N_proc= N_CPU
 

lib/preprocess.py

@@ -20,6 +20,7 @@
 from resampling import resampling
 from dti import dti
 from rish import rish
+from glob import glob
 
 SCRIPTDIR= os.path.dirname(__file__)
 config = ConfigParser()
@@ -44,10 +45,12 @@ def dti_harm(imgPath, maskPath):
     prefix = basename(inPrefix)
 
     outPrefix = os.path.join(directory, 'dti', prefix)
-    dti(imgPath, maskPath, inPrefix, outPrefix)
+    if not glob(outPrefix+'_FA.nii.gz'):
+       dti(imgPath, maskPath, inPrefix, outPrefix)
 
     outPrefix = os.path.join(directory, 'harm', prefix)
-    rish(imgPath, maskPath, inPrefix, outPrefix, N_shm)
+    if not glob(outPrefix+'_L0.nii.gz'):
+       rish(imgPath, maskPath, inPrefix, outPrefix, N_shm)