Merge pull request #16 from pnlbwh/fix-force

Fix force

Author: tashrifbillah

lib/buildTemplate.py

@@ -95,7 +95,7 @@ def antsMult(caselist, outPrefix):
         print(f'See {logFile} for details of template construction')
 
 
-    N_core=os.getenv('TEMPLATE_CONSTRUCT_CORE')
+    N_core=getenv('TEMPLATE_CONSTRUCT_CORE')
     check_call((' ').join([pjoin(SCRIPTDIR, 'antsMultivariateTemplateConstruction2_fixed_random_seed.sh'),
                            '-d', '3',
                            '-g', '0.2',

lib/debug_fa.py

@@ -18,30 +18,30 @@
 import multiprocessing
 from util import *
 
-SCRIPTDIR= os.path.dirname(__file__)
-ROOTDIR= os.path.abspath(os.path.join(SCRIPTDIR, '..'))
-mniTmp = os.path.join(ROOTDIR, 'IITAtlas', 'IITmean_FA.nii.gz')
+SCRIPTDIR= dirname(__file__)
+ROOTDIR= abspath(pjoin(SCRIPTDIR, '..'))
+mniTmp = pjoin(ROOTDIR, 'IITAtlas', 'IITmean_FA.nii.gz')
 
 config = ConfigParser()
-config.read(f'/tmp/harm_config_{os.getpid()}.ini')
+config.read(f'/tmp/harm_config_{getpid()}.ini')
 N_proc = int(config['DEFAULT']['N_proc'])
 diffusionMeasures = [x for x in config['DEFAULT']['diffusionMeasures'].split(',')]
 bshell_b = int(config['DEFAULT']['bshell_b'])
 
 def register_reference(imgPath, warp2mni, trans2mni, templatePath):
 
     print(f'Warping {imgPath} diffusion measures to standard space')
-    directory = os.path.dirname(imgPath)
+    directory = dirname(imgPath)
     inPrefix = imgPath.split('.nii')[0]
     prefix = basename(inPrefix)
 
     for dm in diffusionMeasures:
 
-        output = os.path.join(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
+        output = pjoin(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
 
         # reference site have been already warped to reference template space in buildTemplate.py: warp_bands()
-        # warped data are os.path.join(templatePath, prefix, prefix + f'_WarpedFA.nii.gz')
-        moving = os.path.join(templatePath, prefix + f'_Warped{dm}.nii.gz')
+        # warped data are pjoin(templatePath, prefix, prefix + f'_WarpedFA.nii.gz')
+        moving = pjoin(templatePath, prefix + f'_Warped{dm}.nii.gz')
 
         # so warp diffusion measure to MNI space directly
         antsApplyTransforms[
@@ -56,22 +56,22 @@ def register_reference(imgPath, warp2mni, trans2mni, templatePath):
 def register_target(imgPath, templatePath):
 
     print(f'Warping {imgPath} diffusion measures to standard space')
-    directory = os.path.dirname(imgPath)
+    directory = dirname(imgPath)
     inPrefix = imgPath.split('.nii')[0]
     prefix = basename(inPrefix)
 
-    dmImg = os.path.join(directory, 'dti', prefix + f'_FA.nii.gz')
-    outPrefix = os.path.join(templatePath, prefix.replace(f'_b{bshell_b}','') + '_FA_ToMNI')
+    dmImg = pjoin(directory, 'dti', prefix + f'_FA.nii.gz')
+    outPrefix = pjoin(templatePath, prefix.replace(f'_b{bshell_b}','') + '_FA_ToMNI')
     warp2mni = outPrefix + '1Warp.nii.gz'
     trans2mni = outPrefix + '0GenericAffine.mat'
     # unprocessed target data is given, so in case multiple debug is needed, pass the registration
-    if not os.path.exists(warp2mni):
+    if not exists(warp2mni):
         antsReg(mniTmp, None, dmImg, outPrefix)
 
     for dm in diffusionMeasures:
-        output = os.path.join(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
+        output = pjoin(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
 
-        moving = os.path.join(directory, 'dti', prefix + f'_{dm}.nii.gz')
+        moving = pjoin(directory, 'dti', prefix + f'_{dm}.nii.gz')
         # warp diffusion measure to template space first, then to MNI space
         antsApplyTransforms[
             '-d', '3',
@@ -85,25 +85,25 @@ def register_target(imgPath, templatePath):
 def register_harmonized(imgPath, warp2mni, trans2mni, templatePath, siteName):
 
     print(f'Warping {imgPath} diffusion measures to standard space')
-    directory = os.path.dirname(imgPath)
+    directory = dirname(imgPath)
     inPrefix = imgPath.split('.nii')[0]
     prefix = basename(inPrefix)
 
-    dmImg = os.path.join(directory, 'dti', prefix + f'_FA.nii.gz')
-    dmTmp = os.path.join(templatePath, f'Mean_{siteName}_FA_b{bshell_b}.nii.gz')
-    maskTmp = os.path.join(templatePath, f'{siteName}_Mask.nii.gz')
-    outPrefix = os.path.join(templatePath, prefix.replace(f'_b{bshell_b}','') + '_FA')
+    dmImg = pjoin(directory, 'dti', prefix + f'_FA.nii.gz')
+    dmTmp = pjoin(templatePath, f'Mean_{siteName}_FA_b{bshell_b}.nii.gz')
+    maskTmp = pjoin(templatePath, f'{siteName}_Mask.nii.gz')
+    outPrefix = pjoin(templatePath, prefix.replace(f'_b{bshell_b}','') + '_FA')
     warp2tmp = outPrefix + '1Warp.nii.gz'
     trans2tmp = outPrefix + '0GenericAffine.mat'
 
     # check existence of transforms created with _b{bmax}
-    if not os.path.exists(warp2tmp):
+    if not exists(warp2tmp):
         antsReg(dmTmp, maskTmp, dmImg, outPrefix)
 
     for dm in diffusionMeasures:
-        output = os.path.join(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
+        output = pjoin(templatePath, prefix + f'_InMNI_{dm}.nii.gz')
 
-        moving = os.path.join(directory, 'dti', prefix + f'_{dm}.nii.gz')
+        moving = pjoin(directory, 'dti', prefix + f'_{dm}.nii.gz')
         # warp diffusion measure to template space first, then to MNI space
         antsApplyTransforms[
             '-d', '3',
@@ -117,14 +117,14 @@ def register_harmonized(imgPath, warp2mni, trans2mni, templatePath, siteName):
 def sub2tmp2mni(templatePath, siteName, caselist, ref= False, tar_unproc= False, tar_harm= False):
 
     # obtain the transform
-    moving = os.path.join(templatePath, f'Mean_{siteName}_FA_b{bshell_b}.nii.gz')
+    moving = pjoin(templatePath, f'Mean_{siteName}_FA_b{bshell_b}.nii.gz')
 
-    outPrefix= os.path.join(templatePath, f'TemplateToMNI_{siteName}')
+    outPrefix= pjoin(templatePath, f'TemplateToMNI_{siteName}')
     warp2mni= outPrefix+'1Warp.nii.gz'
     trans2mni= outPrefix+'0GenericAffine.mat'
 
     # check existence of transforms created with _b{bmax}
-    if not os.path.exists(warp2mni):
+    if not exists(warp2mni):
         antsReg(mniTmp, None, moving, outPrefix)
 
     imgs, _= read_caselist(caselist)
@@ -160,7 +160,7 @@ def analyzeStat(file, templatePath):
     :return: mean of the images
     '''
 
-    skel= load(os.path.join(ROOTDIR, 'IITAtlas', 'IITmean_FA_skeleton.nii.gz'))
+    skel= load(pjoin(ROOTDIR, 'IITAtlas', 'IITmean_FA_skeleton.nii.gz'))
     skel_mask= (skel.get_data()>0)*1.
 
     imgs, _ = read_caselist(file)
@@ -170,7 +170,7 @@ def analyzeStat(file, templatePath):
         inPrefix = imgPath.split('.nii')[0]
         prefix = basename(inPrefix)
 
-        faImg= os.path.join(templatePath, prefix + f'_InMNI_FA.nii.gz')
+        faImg= pjoin(templatePath, prefix + f'_InMNI_FA.nii.gz')
         data= load(faImg).get_data()
         temp= data*skel_mask
         meanAttr.append(temp[temp>0].mean())

lib/harm_plot.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import pandas as pd
+from conversion import read_imgs_masks
+from os.path import isfile, basename
+
+XERR=0.1
+ELINEWIDTH=3
+CAPSIZE=5
+CAPTHICK=3
+FMT='cD'
+
+def harm_plot(ydata, labels, outPrefix, bshell_b):
+    '''
+    :param ydata: list of [y1, y2, y3, ...] where each yi is a list
+    :param labels: list of strings
+    :param outPrefix:
+    :return:
+    '''
+    outPrefix += f'_b{bshell_b}'
+    labels= list(labels)
+
+    num_series, num_sub= np.shape(ydata)
+
+    iter_obj= [i for i in range(num_series)]
+
+    # errorbar plot
+    plt.figure(1)
+    plt.grid(True)
+    for i in iter_obj:
+        x= list(i*np.ones((num_sub,)))
+        y= ydata[i]
+        plt.plot(x, y, 'r*')
+
+        plt.errorbar([i], np.mean(y), xerr=XERR, yerr=np.std(y),
+                     ecolor='k', capsize=CAPSIZE, capthick=CAPTHICK, elinewidth=ELINEWIDTH, fmt=FMT)
+
+    plt.xticks(iter_obj, labels)
+    plt.title('Comparison of meanFA before and after harmonization')
+    plt.ylabel('meanFA over IIT_mean_FA_skeleton')
+    plt.savefig(outPrefix+'_ebarplot.png')
+    # plt.show()
+
+    # box plot
+    # plt.figure(2)
+    # plt.grid(True)
+    # for i in iter_obj:
+    #     x = list(i * np.ones((num_sub,)))
+    #     y = ydata[i]
+    #     plt.plot(x, y, 'r*')
+    #
+    # plt.boxplot(ydata, labels=labels, positions=iter_obj,
+    #             boxprops=dict(linewidth=4),
+    #             medianprops=dict(linewidth=4),
+    #             whiskerprops=dict(linewidth=2))
+    #
+    #
+    # plt.title(f'Comparison of boxplot before and after harmonization for b{bshell_b}')
+    # plt.ylabel('meanFA over IIT_mean_FA_skeleton')
+    # plt.savefig(outPrefix+'_boxplot.png')
+    # plt.show()
+    # return (outPrefix+'_ebarplot.png', outPrefix+'_boxplot.png')
+
+    return outPrefix+'_ebarplot.png'
+
+
+def generate_csv(imgs, site_means, outPrefix, bshell_b):
+
+    try:
+        imgs, _= read_imgs_masks(imgs)
+    except:
+        pass
+
+    statFile = outPrefix + '_stat.csv'
+    
+    if isfile(statFile):
+        df= pd.read_csv(statFile)
+        df= df.assign(**{f'meanFA b{bshell_b}':site_means})
+    else:
+        stat = {'subject': [basename(f) for f in imgs], f'meanFA b{bshell_b}': site_means}
+        df = pd.DataFrame(stat)
+    
+    df.to_csv(statFile, index=False)
+
+
+if __name__=='__main__':
+
+    sub=['hi','hello','go','come']
+    ref_mean= [0.46, 0.49, 0.44, 0.40]
+    target_mean_before= [0.42, 0.58, 0.43, 0.66]
+    target_mean_after= [0.5 , 0.45, 0.40, 0.55]
+    labels=['Reference','Target before','Target after']
+    bshell_b= 2000
+    
+    harm_plot([ref_mean, target_mean_before, target_mean_after], labels, '/tmp/abc', bshell_b)
+    # harm_plot([ref_mean], ['Reference'], '/tmp/abc', bshell_b)
+    
+    generate_csv(sub, ref_mean, '/tmp/abc', bshell_b)
+    
+    

lib/harmonization.py

@@ -16,6 +16,7 @@
 from plumbum import cli
 from distutils.spawn import find_executable
 import multiprocessing, psutil
+import io
 
 from determineNshm import verifyNshmForAll, determineNshm
 from util import *
@@ -254,13 +255,12 @@ def harmonizeData(self):
         if not exists(self.templatePath):
             raise NotADirectoryError(f'{self.templatePath} does not exist')
         else:
-            if not os.listdir(self.templatePath):
+            if not listdir(self.templatePath):
                 raise ValueError(f'{self.templatePath} is empty')
 
         # go through each file listed in csv, check their existence, create dti and harm directories
         check_csv(self.target_csv, self.force)
 
-
         # target data is not manipulated in multi-shell-dMRIharmonization i.e. bvalMapped, resampled, nor denoised
         # this block may be uncommented in a future design
         # from preprocess import dti_harm
@@ -281,13 +281,10 @@ def harmonizeData(self):
 
         # read target image list
         moving= pjoin(self.templatePath, f'Mean_{self.target}_FA_b{self.bshell_b}.nii.gz')
-        imgs, masks= read_caselist(self.target_csv)
+        imgs, masks= read_caselist(self.tar_unproc_csv)
 
-        preFlag= 1 # omit preprocessing of target data again
-        if self.target_csv.endswith('.modified'):
-            preFlag= 0
-        else:
-            # this file will be used later for debugging
+        fm= None
+        if not self.target_csv.endswith('.modified'):
             self.target_csv += '.modified'
             fm = open(self.target_csv, 'w')
 
@@ -297,12 +294,12 @@ def harmonizeData(self):
         pool = multiprocessing.Pool(self.N_proc)
         res= []
         for imgPath, maskPath in zip(imgs, masks):
-            res.append(pool.apply_async(func= reconst, args= (imgPath, maskPath, moving, self.templatePath, preFlag, )))
+            res.append(pool.apply_async(func= reconst, args= (imgPath, maskPath, moving, self.templatePath,)))
 
         for r in res:
             imgPath, maskPath, harmImg, harmMask= r.get()
 
-            if preFlag:
+            if isinstance(fm, io.IOBase):
                 fm.write(imgPath + ',' + maskPath + '\n')
             fh.write(harmImg + ',' + harmMask + '\n')
 
@@ -314,17 +311,17 @@ def harmonizeData(self):
         # loop for debugging
         # res= []
         # for imgPath, maskPath in zip(imgs, masks):
-        #     res.append(reconst(imgPath, maskPath, moving, self.templatePath, preFlag))
+        #     res.append(reconst(imgPath, maskPath, moving, self.templatePath))
         #
         # for r in res:
         #     imgPath, maskPath, harmImg, harmMask= r
         #
-        #     if preFlag:
+        #     if isinstance(fm, io.IOBase):
         #         fm.write(imgPath + ',' + maskPath + '\n')
         #     fh.write(harmImg + ',' + harmMask + '\n')
 
 
-        if preFlag:
+        if isinstance(fm, io.IOBase):
             fm.close()
         fh.close()
 
@@ -486,7 +483,7 @@ def main(self):
 
 
         # write config file to temporary directory
-        configFile= f'/tmp/harm_config_{os.getpid()}.ini'
+        configFile= f'/tmp/harm_config_{getpid()}.ini'
         with open(configFile,'w') as f:
             f.write('[DEFAULT]\n')
             f.write(f'N_shm = {self.N_shm}\n')
@@ -498,6 +495,7 @@ def main(self):
             f.write(f'denoise = {1 if self.denoise else 0}\n')
             f.write(f'travelHeads = {1 if self.travelHeads else 0}\n')
             f.write(f'debug = {1 if self.debug else 0}\n')
+            f.write(f'force = {1 if self.force else 0}\n')
             f.write(f'verbose = {1 if self.verbose else 0}\n')
             f.write('diffusionMeasures = {}\n'.format((',').join(self.diffusionMeasures)))
 
@@ -506,6 +504,13 @@ def main(self):
 
         if self.create:
             self.createTemplate()
+            import fileinput
+            for line in fileinput.input(configFile, inplace=True):
+                if 'force' in line:
+                    print('force = 0')
+                else:
+                    print(line)
+            self.force= False
 
         if self.process:
             self.harmonizeData()
@@ -514,7 +519,7 @@ def main(self):
             self.post_debug()
 
 
-        os.remove(configFile)
+        remove(configFile)
 
 
 if __name__ == '__main__':