From bd7de807e20c8f7f98900f5e270c3c81d01dad0e Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Wed, 27 Nov 2019 15:56:00 -0500 Subject: [PATCH 01/19] remove scale clipping, suggestion by @suheyla2 --- lib/buildTemplate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index 3bb0474..5de173b 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -183,7 +183,7 @@ def stat_calc(ref, target, mask): np.nan_to_num(per_diff).clip(max=100., min=-100., out= per_diff) per_diff_smooth= smooth(per_diff) scale= ref/(target+eps) - scale.clip(max=10., min= 0., out= scale) + # scale.clip(max=10., min= 0., out= scale) return (delta, per_diff, per_diff_smooth, scale) From e70f31bbf3bafd45159ba7f6236680ed20aa52e5 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Wed, 27 Nov 2019 16:12:02 -0500 Subject: [PATCH 02/19] use Popen instead of check_call to launch spm_bspline_exec --- lib/resampling.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/resampling.py b/lib/resampling.py index 019055b..fdda527 100644 --- a/lib/resampling.py +++ b/lib/resampling.py @@ -18,6 +18,7 @@ from scipy.io import loadmat, savemat from normalize import normalize_data, find_b0 from util import * +from subprocess import Popen FILEDIR= os.path.dirname(os.path.abspath(__file__)) @@ -28,7 +29,9 @@ def resize_spm(lowResImg, inPrefix): savemat(dataFile, {'lowResImg': lowResImg}) # call MATLAB_Runtime based spm bspline interpolation - check_call([os.path.join(FILEDIR,'spm_bspline_exec', 'bspline')+' '+inPrefix], shell= True) + cmd= (' ').join([os.path.join(FILEDIR,'spm_bspline_exec', 'bspline'), inPrefix]) + p= Popen(cmd, shell= True) + p.wait() highResImg= np.nan_to_num(loadmat(inPrefix+'_resampled.mat')['highResImg']) From 993d66c28f21fcd8d1e90b206b90423f0b015f62 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Thu, 28 Nov 2019 15:04:49 -0500 Subject: [PATCH 03/19] use Popen instead of check_call --- lib/reconstSignal.py | 6 ++++-- lib/resampling.py | 8 ++++---- lib/tests/fa_skeleton_test.py | 6 ++++-- lib/tests/test_util.py | 2 +- lib/util.py | 2 +- 5 files changed, 14 insertions(+), 10 deletions(-) diff --git a/lib/reconstSignal.py b/lib/reconstSignal.py index aca377a..d7ab6be 100644 --- a/lib/reconstSignal.py +++ b/lib/reconstSignal.py @@ -32,20 +32,22 @@ def antsReg(img, mask, mov, outPrefix): if mask: - check_call((' ').join(['antsRegistrationSyNQuick.sh', + p= Popen((' ').join(['antsRegistrationSyNQuick.sh', '-d', '3', '-f', img, '-x', mask, '-m', mov, '-o', outPrefix, '-e', '123456']), shell= True) + p.wait() else: - check_call((' ').join(['antsRegistrationSyNQuick.sh', + p= Popen((' ').join(['antsRegistrationSyNQuick.sh', '-d', '3', '-f', img, '-m', mov, '-o', outPrefix, '-e', '123456']), shell= True) + p.wait() def antsApply(templatePath, directory, prefix): diff --git a/lib/resampling.py b/lib/resampling.py index fdda527..0125e89 100644 --- a/lib/resampling.py +++ b/lib/resampling.py @@ -18,7 +18,6 @@ from scipy.io import loadmat, savemat from normalize import normalize_data, find_b0 from util import * -from subprocess import Popen FILEDIR= os.path.dirname(os.path.abspath(__file__)) @@ -29,8 +28,7 @@ def resize_spm(lowResImg, inPrefix): savemat(dataFile, {'lowResImg': lowResImg}) # call MATLAB_Runtime based spm bspline interpolation - cmd= (' ').join([os.path.join(FILEDIR,'spm_bspline_exec', 'bspline'), inPrefix]) - p= Popen(cmd, shell= True) + p= Popen((' ').join([os.path.join(FILEDIR,'spm_bspline_exec', 'bspline'), inPrefix]), shell=True) p.wait() highResImg= np.nan_to_num(loadmat(inPrefix+'_resampled.mat')['highResImg']) @@ -136,7 +134,9 @@ def resampling(lowResImgPath, lowResMaskPath, lowResImg, lowResImgHdr, lowResMas # unring the b0 highResB0Path = lowResImgPath.split('.')[0] + '_resampled_bse.nii.gz' - check_call(['unring.a64', highResB0PathTmp, highResB0Path]) + p= Popen((' ').join(['unring.a64', highResB0PathTmp, highResB0Path]), shell= True) + p.wait() + check_call(['rm', highResB0PathTmp]) b0_gibs = load(highResB0Path).get_data() np.nan_to_num(b0_gibs).clip(min= 0., out= b0_gibs) # using min= 1. is unnecessary diff --git a/lib/tests/fa_skeleton_test.py b/lib/tests/fa_skeleton_test.py index 58671a6..d24d82b 100755 --- a/lib/tests/fa_skeleton_test.py +++ b/lib/tests/fa_skeleton_test.py @@ -42,20 +42,22 @@ def printStat(ref_mean, imgs): def antsReg(img, mask, mov, outPrefix, n_thread=1): if mask: - check_call((' ').join(['antsRegistrationSyNQuick.sh', + p= Popen((' ').join(['antsRegistrationSyNQuick.sh', '-d', '3', '-f', img, '-x', mask, '-m', mov, '-n', str(n_thread), '-o', outPrefix]), shell= True) + p.wait() else: - check_call((' ').join(['antsRegistrationSyNQuick.sh', + p= Popen((' ').join(['antsRegistrationSyNQuick.sh', '-d', '3', '-f', img, '-m', mov, '-n', str(n_thread), '-o', outPrefix]), shell= True) + p.wait() def register_subject(imgPath, warp2mni, trans2mni, templatePath, siteName): diff --git a/lib/tests/test_util.py b/lib/tests/test_util.py index 90f8d77..bc7e90a 100644 --- a/lib/tests/test_util.py +++ b/lib/tests/test_util.py @@ -16,7 +16,7 @@ import sys from shutil import rmtree, copyfile import unittest -from subprocess import check_call +from subprocess import check_call, Popen FILEDIR= abspath(dirname(__file__)) LIBDIR= dirname(FILEDIR) diff --git a/lib/util.py b/lib/util.py index b9c4336..1f9494b 100644 --- a/lib/util.py +++ b/lib/util.py @@ -13,7 +13,7 @@ import os, shutil, configparser, sys import numpy as np -from subprocess import check_call +from subprocess import check_call, Popen import warnings with warnings.catch_warnings(): From 3c61a4d899d20f509f7cbd93eb11e3437c8c1ed7 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Thu, 28 Nov 2019 17:30:29 -0500 Subject: [PATCH 04/19] write statistics in a file for future --- lib/harmonization.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/lib/harmonization.py b/lib/harmonization.py index 3bbf871..2996fd8 100755 --- a/lib/harmonization.py +++ b/lib/harmonization.py @@ -361,9 +361,17 @@ def post_debug(self): def showStat(self): from debug_fa import analyzeStat + from datetime import datetime print('\n\nPrinting statistics :\n\n') - + + # save statistics for future + statFile= os.path.join(self.templatePath, 'meanFAstat.txt') + f= open(statFile,'w') + stdout= sys.stdout + sys.stdout= f + + print(datetime.now().strftime('%c'),'\n') print(f'{self.reference} site: ') ref_mean = analyzeStat(self.ref_csv, self.templatePath) printStat(ref_mean, self.ref_csv) @@ -375,6 +383,15 @@ def showStat(self): print(f'{self.target} site after harmonization: ') target_mean_after = analyzeStat(self.harm_csv, self.templatePath) printStat(target_mean_after, self.harm_csv) + + f.close() + sys.stdout= stdout + + # print statistics on console + with open(statFile) as f: + print(f.read()) + + print('\nThe statistics are also saved in ', statFile) def sanityCheck(self): From a6b66f39faf1b76e0c092c376dc081795c5f5622 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 13:15:53 -0500 Subject: [PATCH 05/19] decouple rish compute from spm_bspline interp, use hard-coded 4 threads for rish compute, introduce env var for num of cores to use by antsMult --- lib/buildTemplate.py | 6 +++--- lib/preprocess.py | 48 ++++++++++++++++++++++++++++---------------- lib/reconstSignal.py | 2 +- 3 files changed, 35 insertions(+), 21 deletions(-) diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index 5de173b..860db78 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -22,7 +22,6 @@ eps= 2.2204e-16 SCRIPTDIR= os.path.dirname(__file__) config = configparser.ConfigParser() -# config.read(os.path.join(SCRIPTDIR,'config.ini')) config.read(f'/tmp/harm_config_{os.getpid()}.ini') N_shm = int(config['DEFAULT']['N_shm']) N_proc = int(config['DEFAULT']['N_proc']) @@ -83,7 +82,8 @@ def createAntsCaselist(imgs, file): def antsMult(caselist, outPrefix): - + + N_core=os.getenv('TEMPLATE_CONSTRUCT_CORE') check_call((' ').join([os.path.join(SCRIPTDIR, 'antsMultivariateTemplateConstruction2_fixed_random_seed.sh'), '-d', '3', '-g', '0.2', @@ -91,7 +91,7 @@ def antsMult(caselist, outPrefix): '-t', "BSplineSyN[0.1,26,0]", '-r', '1', '-c', '2', - '-j', str(N_proc), + '-j', str(N_core) if N_core else str(N_proc), '-f', '8x4x2x1', '-o', outPrefix, caselist]), shell= True) diff --git a/lib/preprocess.py b/lib/preprocess.py index a1e4c60..8996b12 100644 --- a/lib/preprocess.py +++ b/lib/preprocess.py @@ -55,7 +55,7 @@ def read_caselist(file): return (imgs, masks) -def dti_harm(imgPath, maskPath): +def dti_harm(imgPath, maskPath, nargout=0): directory = os.path.dirname(imgPath) inPrefix = imgPath.split('.')[0] @@ -69,15 +69,10 @@ def dti_harm(imgPath, maskPath): outPrefix = os.path.join(directory, 'harm', prefix) b0, shm_coeff, qb_model= rish(imgPath, maskPath, inPrefix, outPrefix, N_shm) + + if nargout==3: + return (b0, shm_coeff, qb_model) - return (b0, shm_coeff, qb_model) - - -# def pre_dti_harm(imgPath, maskPath): -def pre_dti_harm(itr): - imgPath, maskPath = preprocessing(itr[0], itr[1]) - dti_harm(imgPath, maskPath) - return (imgPath, maskPath) # convert NRRD to NIFTI on the fly def nrrd2nifti(imgPath): @@ -163,21 +158,40 @@ def preprocessing(imgPath, maskPath): def common_processing(caselist): - imgs, masks = read_caselist(caselist) - f = open(caselist + '.modified', 'w') - pool = multiprocessing.Pool(N_proc) # Use all available cores, otherwise specify the number you want as an argument + imgs, masks = read_caselist(caselist) + + # to avoid MemoryError, decouple preprocessing (spm_bspline) and dti_harm (rish) + res=[] + pool = multiprocessing.Pool(N_proc) + for imgPath,maskPath in zip(imgs,masks): + res.append(pool.apply_async(func= preprocessing, args= (imgPath,maskPath))) + + attributes= [r.get() for r in res] + + pool.close() + pool.join() - res = pool.map_async(pre_dti_harm, np.hstack((np.reshape(imgs, (len(imgs), 1)), np.reshape(masks, (len(masks), 1))))) - attributes = res.get() + + f = open(caselist + '.modified', 'w') for i in range(len(imgs)): imgs[i] = attributes[i][0] masks[i] = attributes[i][1] f.write(f'{imgs[i]},{masks[i]}\n') - + f.close() + + + # the following imgs, masks is for diagnosing MemoryError i.e. computing rish w/o preprocessing + # to diagnose, comment all the above and uncomment the following + # imgs, masks = read_caselist(caselist+'.modified') + + # experimentally found ncpu=4 to be memroy optimal + pool = multiprocessing.Pool(4) + for imgPath,maskPath in zip(imgs,masks): + pool.apply_async(func= dti_harm, args= (imgPath,maskPath)) pool.close() pool.join() + - f.close() + return (imgs, masks) - return (imgs, masks) \ No newline at end of file diff --git a/lib/reconstSignal.py b/lib/reconstSignal.py index d7ab6be..7b875e8 100644 --- a/lib/reconstSignal.py +++ b/lib/reconstSignal.py @@ -168,7 +168,7 @@ def reconst(imgPath, maskPath, moving, templatePath, preFlag): # provide full sampled shm_coeff, qb_model.B # provide imgPath header img = load(imgPath) - b0, shm_coeff, qb_model = dti_harm(imgPath, maskPath) + b0, shm_coeff, qb_model = dti_harm(imgPath, maskPath, nargout=3) directory = os.path.dirname(imgPath) inPrefix = imgPath.split('.')[0] From 7c00f0916b70cbcdeeb7d461028d5a73879690ea Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 14:19:55 -0500 Subject: [PATCH 06/19] print a msg when travelHeads is true --- lib/buildTemplate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index 860db78..6e947c5 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -211,6 +211,7 @@ def difference_calc(refSite, targetSite, refImgs, targetImgs, per_diff_smooth= [] scale= [] if travelHeads: + print('Using travelHeads for computing scale maps') for refImg, targetImg in zip(refImgs, targetImgs): prefix = os.path.basename(refImg).split('.')[0] ref= load_nifti(os.path.join(templatePath, f'{prefix}_Warped{dm}.nii.gz'))[0] From 460bfd1603deb30010c7899c54205e2eebe1a3f9 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 15:23:22 -0500 Subject: [PATCH 07/19] clip scale map to min 0 --- lib/buildTemplate.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index 6e947c5..a99f569 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -183,6 +183,7 @@ def stat_calc(ref, target, mask): np.nan_to_num(per_diff).clip(max=100., min=-100., out= per_diff) per_diff_smooth= smooth(per_diff) scale= ref/(target+eps) + scale.clip(min=0., out=scale) # scale.clip(max=10., min= 0., out= scale) return (delta, per_diff, per_diff_smooth, scale) @@ -211,7 +212,7 @@ def difference_calc(refSite, targetSite, refImgs, targetImgs, per_diff_smooth= [] scale= [] if travelHeads: - print('Using travelHeads for computing scale maps') + print('Using travelHeads for computing scale maps of',dm) for refImg, targetImg in zip(refImgs, targetImgs): prefix = os.path.basename(refImg).split('.')[0] ref= load_nifti(os.path.join(templatePath, f'{prefix}_Warped{dm}.nii.gz'))[0] From d41abc281dbf765b2374ffb6d00168ed7c46c764 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 15:34:12 -0500 Subject: [PATCH 08/19] use --travelHeads in one of the tests --- lib/tests/pipeline_test.sh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lib/tests/pipeline_test.sh b/lib/tests/pipeline_test.sh index cf0b8c5..10cf5ea 100755 --- a/lib/tests/pipeline_test.sh +++ b/lib/tests/pipeline_test.sh @@ -55,7 +55,6 @@ write_list prisma.txt ### Run pipeline and obtain statistics when same number of matched reference and target images are used in ### tempalate creation and harmonization - # run test ../../harmonization.py \ --bvalMap 1000 \ @@ -87,7 +86,6 @@ rm -rf template && \ ### Run pipeline and obtain statistics when small set of matched reference and target images are used in template creation ### and a larger set (does not have to be mutually exclusive from the former) of target images are used in harmonization - # --create and --debug block ../../harmonization.py \ --bvalMap 1000 \ @@ -97,6 +95,7 @@ rm -rf template && \ --tar_list prisma.txt \ --ref_name CONNECTOM \ --tar_name PRISMA \ +--travelHeads \ --nproc -1 \ --create --debug || EXIT 'harmonization.py with --create --debug failed' From ca4a036f9c32ad8611eb00562f3c065175148de3 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 15:34:41 -0500 Subject: [PATCH 09/19] minor correction --- lib/harmonization.py | 6 +++--- lib/reconstSignal.py | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/lib/harmonization.py b/lib/harmonization.py index 2996fd8..d03eec8 100755 --- a/lib/harmonization.py +++ b/lib/harmonization.py @@ -232,7 +232,7 @@ def createTemplate(self): # ATTN: antsMultivariateTemplateConstruction2.sh requires absolute path for caselist antsMult(os.path.abspath(antsMultCaselist), self.templatePath) - # # load templateHdr + # load templateHdr templateHdr= load(os.path.join(self.templatePath, 'template0.nii.gz')).header @@ -244,9 +244,9 @@ def createTemplate(self): pool.close() pool.join() - print('dti statistics: mean, std(FA, MD) calculation of reference site') + print('dti statistics: mean, std calculation of reference site') refMaskPath= dti_stat(self.reference, refImgs, refMasks, self.templatePath, templateHdr) - print('dti statistics: mean, std(FA, MD) calculation of target site') + print('dti statistics: mean, std calculation of target site') targetMaskPath= dti_stat(self.target, targetImgs, targetMasks, self.templatePath, templateHdr) print('masking dti statistics of reference site') diff --git a/lib/reconstSignal.py b/lib/reconstSignal.py index 7b875e8..4f6688a 100644 --- a/lib/reconstSignal.py +++ b/lib/reconstSignal.py @@ -22,7 +22,6 @@ eps= 2.2204e-16 SCRIPTDIR= os.path.dirname(__file__) config = configparser.ConfigParser() -# config.read(os.path.join(SCRIPTDIR,'config.ini')) config.read(f'/tmp/harm_config_{os.getpid()}.ini') N_shm = int(config['DEFAULT']['N_shm']) N_proc = int(config['DEFAULT']['N_proc']) From e9b87198e84ca23d57ee4abe2bd2b1d2d8f6c1eb Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Mon, 2 Dec 2019 15:45:39 -0500 Subject: [PATCH 10/19] document TEMPLATE_CONSTRUCT_CORES env var --- README.md | 6 ++++++ lib/buildTemplate.py | 2 +- lib/tests/pipeline_test.sh | 4 ++++ 3 files changed, 11 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 4efa130..20ea014 100644 --- a/README.md +++ b/README.md @@ -449,6 +449,12 @@ is advisable to leave out at least two cores for other processes to run smoothly **NOTE** See [Caveats/Issues](#caveatsissues) that may occur while using many processors in parallel. +Furthermore, you can define the environment variable `TEMPLATE_CONSTRUCT_CORES` to use a different number of processors +for `antsMultivariateTemplateConstruction2.sh` independent of `--nproc` used for rest of the processes in *dMRIharmonization*: + + export TEMPLATE_CONSTRUCT_CORES=32 + + # Order of spherical harmonics RISH features are derived from spherical harmonic coefficients. The order of spherical harmonic coefficients you can use diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index a99f569..0031524 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -83,7 +83,7 @@ def createAntsCaselist(imgs, file): def antsMult(caselist, outPrefix): - N_core=os.getenv('TEMPLATE_CONSTRUCT_CORE') + N_core=os.getenv('TEMPLATE_CONSTRUCT_CORES') check_call((' ').join([os.path.join(SCRIPTDIR, 'antsMultivariateTemplateConstruction2_fixed_random_seed.sh'), '-d', '3', '-g', '0.2', diff --git a/lib/tests/pipeline_test.sh b/lib/tests/pipeline_test.sh index 10cf5ea..ce299f5 100755 --- a/lib/tests/pipeline_test.sh +++ b/lib/tests/pipeline_test.sh @@ -86,6 +86,10 @@ rm -rf template && \ ### Run pipeline and obtain statistics when small set of matched reference and target images are used in template creation ### and a larger set (does not have to be mutually exclusive from the former) of target images are used in harmonization + +# test the following advanced parameter +export TEMPLATE_CONSTRUCT_CORES=32 + # --create and --debug block ../../harmonization.py \ --bvalMap 1000 \ From 5d4174cf17ba5563ef9795441c5d70ebdcbc0ce8 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Wed, 4 Dec 2019 12:19:46 -0500 Subject: [PATCH 11/19] remove dead code --- lib/harmonization.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/lib/harmonization.py b/lib/harmonization.py index d03eec8..acd1968 100755 --- a/lib/harmonization.py +++ b/lib/harmonization.py @@ -342,9 +342,6 @@ def post_debug(self): from debug_fa import sub2tmp2mni - refImgs, _ = read_imgs_masks(self.ref_csv) - targetImgs, _= read_imgs_masks(self.target_csv) - print('\n\n Reference site') sub2tmp2mni(self.templatePath, self.reference, self.ref_csv, ref= True) From 4458fe8ac14e1147f506f54cec20b489d25c92fb Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Wed, 4 Dec 2019 16:10:00 -0500 Subject: [PATCH 12/19] add w/space --- lib/tests/fa_skeleton_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/tests/fa_skeleton_test.py b/lib/tests/fa_skeleton_test.py index d24d82b..b632b1e 100755 --- a/lib/tests/fa_skeleton_test.py +++ b/lib/tests/fa_skeleton_test.py @@ -143,8 +143,8 @@ def main(): and then to MNI space. Finally, calculates mean FA over IITmean_FA_skeleton.nii.gz''') parser.add_argument('-i', '--input', type=str, required=True, help='a .txt/.csv file having one column for FA imgs, ' - 'or two columns for (img,mask) pair, the latter list is what you used in/obtained from harmonization.py' - 'see documentation for more details') + 'or two columns for (img,mask) pair, the latter list is what you used in/obtained from harmonization.py. ' + 'See documentation for more details') parser.add_argument('-s', '--site', type= str, required=True, help='site name for locating template FA and mask in tempalte directory') parser.add_argument('-t', '--template', type=str, required=True, From 8577f3ea3eae61c0d4d59d15c9c058c585f814c8 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Wed, 4 Dec 2019 16:40:31 -0500 Subject: [PATCH 13/19] always append to statFile --- lib/harmonization.py | 2 +- lib/tests/fa_skeleton_test.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/harmonization.py b/lib/harmonization.py index acd1968..18dfa3c 100755 --- a/lib/harmonization.py +++ b/lib/harmonization.py @@ -364,7 +364,7 @@ def showStat(self): # save statistics for future statFile= os.path.join(self.templatePath, 'meanFAstat.txt') - f= open(statFile,'w') + f= open(statFile,'a') stdout= sys.stdout sys.stdout= f diff --git a/lib/tests/fa_skeleton_test.py b/lib/tests/fa_skeleton_test.py index b632b1e..66411a1 100755 --- a/lib/tests/fa_skeleton_test.py +++ b/lib/tests/fa_skeleton_test.py @@ -108,8 +108,8 @@ def sub2tmp2mni(templatePath, siteName, faImgs, N_proc): pool= multiprocessing.Pool(N_proc) res=[] for imgPath in faImgs: - res.append(pool.apply_async(func= register_subject, - args= (imgPath, warp2mni, trans2mni, templatePath, siteName, ))) + res.append(pool.apply_async(func= register_subject, + args= (imgPath, warp2mni, trans2mni, templatePath, siteName, ))) mniFAimgs= [r.get() for r in res] From 464f728848a54896a9f0061ec986345bf461f742 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 16:10:31 -0500 Subject: [PATCH 14/19] remove scale clipping, compute scale for rish only --- lib/buildTemplate.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/lib/buildTemplate.py b/lib/buildTemplate.py index 0031524..b97fe8d 100644 --- a/lib/buildTemplate.py +++ b/lib/buildTemplate.py @@ -183,8 +183,6 @@ def stat_calc(ref, target, mask): np.nan_to_num(per_diff).clip(max=100., min=-100., out= per_diff) per_diff_smooth= smooth(per_diff) scale= ref/(target+eps) - scale.clip(min=0., out=scale) - # scale.clip(max=10., min= 0., out= scale) return (delta, per_diff, per_diff_smooth, scale) @@ -212,7 +210,7 @@ def difference_calc(refSite, targetSite, refImgs, targetImgs, per_diff_smooth= [] scale= [] if travelHeads: - print('Using travelHeads for computing scale maps of',dm) + print('Using travelHeads for computing templates of',dm) for refImg, targetImg in zip(refImgs, targetImgs): prefix = os.path.basename(refImg).split('.')[0] ref= load_nifti(os.path.join(templatePath, f'{prefix}_Warped{dm}.nii.gz'))[0] @@ -248,8 +246,9 @@ def difference_calc(refSite, targetSite, refImgs, targetImgs, save_nifti(os.path.join(templatePath, f'PercentageDiff_{dm}smooth.nii.gz'), np.mean(per_diff_smooth, axis= 0), templateAffine, templateHdr) - save_nifti(os.path.join(templatePath, f'Scale_{dm}.nii.gz'), - np.sqrt(np.mean(scale, axis= 0)), templateAffine, templateHdr) + if 'L' in dm: + save_nifti(os.path.join(templatePath, f'Scale_{dm}.nii.gz'), + np.sqrt(np.mean(scale, axis= 0)), templateAffine, templateHdr) From 5ada6daee3e5701672973488245f3dc0745f0597 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 16:11:15 -0500 Subject: [PATCH 15/19] remove nan in GFA image --- lib/dti.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/dti.py b/lib/dti.py index 47b4dba..44c9fd9 100644 --- a/lib/dti.py +++ b/lib/dti.py @@ -47,7 +47,7 @@ def dti(imgPath, maskPath, inPrefix, outPrefix, tool='FSL'): ] & FG - gfa_vol = gfa(masked_vol) + gfa_vol = np.nan_to_num(gfa(masked_vol)) save_nifti(outPrefix + '_GFA.nii.gz', gfa_vol, vol.affine, vol.header) From e873835be4fd579abb539c6c0bbe3604734ddc12 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 16:12:22 -0500 Subject: [PATCH 16/19] use 6 cores which is all needed for total 6 images: FA and L0 --- lib/tests/pipeline_test.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/tests/pipeline_test.sh b/lib/tests/pipeline_test.sh index ce299f5..5259520 100755 --- a/lib/tests/pipeline_test.sh +++ b/lib/tests/pipeline_test.sh @@ -88,7 +88,7 @@ rm -rf template && \ ### and a larger set (does not have to be mutually exclusive from the former) of target images are used in harmonization # test the following advanced parameter -export TEMPLATE_CONSTRUCT_CORES=32 +export TEMPLATE_CONSTRUCT_CORES=6 # --create and --debug block ../../harmonization.py \ From 54ea7d6e70531ec5c135e6575b2d56338c991c4c Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 16:15:46 -0500 Subject: [PATCH 17/19] update status msgs --- lib/harmonization.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/harmonization.py b/lib/harmonization.py index 18dfa3c..1dbdcf3 100755 --- a/lib/harmonization.py +++ b/lib/harmonization.py @@ -244,9 +244,9 @@ def createTemplate(self): pool.close() pool.join() - print('dti statistics: mean, std calculation of reference site') + print('calculating dti statistics i.e. mean, std for reference site') refMaskPath= dti_stat(self.reference, refImgs, refMasks, self.templatePath, templateHdr) - print('dti statistics: mean, std calculation of target site') + print('calculating dti statistics i.e. mean, std for target site') targetMaskPath= dti_stat(self.target, targetImgs, targetMasks, self.templatePath, templateHdr) print('masking dti statistics of reference site') @@ -254,16 +254,16 @@ def createTemplate(self): print('masking dti statistics of target site') templateMask= template_masking(refMaskPath, targetMaskPath, self.templatePath, self.target) - print('rish_statistics mean, std(L{i}) calculation of reference site') + print('calculating rish_statistics i.e. mean, std calculation for reference site') rish_stat(self.reference, refImgs, self.templatePath, templateHdr) - print('rish_statistics mean, std(L{i}) calculation of target site') + print('calculating rish_statistics i.e. mean, std calculation for target site') rish_stat(self.target, targetImgs, self.templatePath, templateHdr) - print('calculating scale map for diffusionMeasures') + print('calculating templates for diffusionMeasures') difference_calc(self.reference, self.target, refImgs, targetImgs, self.templatePath, templateHdr, templateMask, self.diffusionMeasures) - print('calculating scale map for rishFeatures') + print('calculating templates for rishFeatures') difference_calc(self.reference, self.target, refImgs, targetImgs, self.templatePath, templateHdr, templateMask, [f'L{i}' for i in range(0, self.N_shm+1, 2)]) From 7f60d0a4bb62225d96086fbff2184056f9eb0dc6 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 17:54:40 -0500 Subject: [PATCH 18/19] decouple dti and rish computation --- lib/preprocess.py | 13 +++---------- lib/reconstSignal.py | 9 +++++---- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/lib/preprocess.py b/lib/preprocess.py index 8996b12..7fc87c6 100644 --- a/lib/preprocess.py +++ b/lib/preprocess.py @@ -23,7 +23,6 @@ SCRIPTDIR= os.path.dirname(__file__) config = configparser.ConfigParser() -# config.read(os.path.join(SCRIPTDIR,'config.ini')) config.read(f'/tmp/harm_config_{os.getpid()}.ini') N_shm = int(config['DEFAULT']['N_shm']) @@ -55,24 +54,18 @@ def read_caselist(file): return (imgs, masks) -def dti_harm(imgPath, maskPath, nargout=0): +def dti_harm(imgPath, maskPath): directory = os.path.dirname(imgPath) inPrefix = imgPath.split('.')[0] prefix = os.path.split(inPrefix)[-1] outPrefix = os.path.join(directory, 'dti', prefix) - - # if the dti output exists with the same prefix, don't dtifit again - if not os.path.exists(outPrefix+'_FA.nii.gz'): - dti(imgPath, maskPath, inPrefix, outPrefix) + dti(imgPath, maskPath, inPrefix, outPrefix) outPrefix = os.path.join(directory, 'harm', prefix) - b0, shm_coeff, qb_model= rish(imgPath, maskPath, inPrefix, outPrefix, N_shm) + rish(imgPath, maskPath, inPrefix, outPrefix, N_shm) - if nargout==3: - return (b0, shm_coeff, qb_model) - # convert NRRD to NIFTI on the fly def nrrd2nifti(imgPath): diff --git a/lib/reconstSignal.py b/lib/reconstSignal.py index 4f6688a..0038a3c 100644 --- a/lib/reconstSignal.py +++ b/lib/reconstSignal.py @@ -18,6 +18,7 @@ from buildTemplate import applyXform from local_med_filter import local_med_filter from preprocess import dti_harm, preprocessing +from rish import rish eps= 2.2204e-16 SCRIPTDIR= os.path.dirname(__file__) @@ -164,14 +165,14 @@ def reconst(imgPath, maskPath, moving, templatePath, preFlag): if preFlag: imgPath, maskPath = preprocessing(imgPath, maskPath) - # provide full sampled shm_coeff, qb_model.B - # provide imgPath header img = load(imgPath) - b0, shm_coeff, qb_model = dti_harm(imgPath, maskPath, nargout=3) directory = os.path.dirname(imgPath) inPrefix = imgPath.split('.')[0] - prefix = os.path.split(inPrefix)[-1] + prefix = os.path.split(inPrefix)[-1] + outPrefix = os.path.join(directory, 'harm', prefix) + b0, shm_coeff, qb_model = rish(imgPath, maskPath, inPrefix, outPrefix, N_shm) + print(f'Registering template FA to {imgPath} space ...') outPrefix = os.path.join(directory, 'harm', 'ToSubjectSpace_' + prefix) From 38b4915348482e123cde193b0bb46b9ecc721b12 Mon Sep 17 00:00:00 2001 From: "Billah, Tashrif" Date: Sat, 7 Dec 2019 17:55:04 -0500 Subject: [PATCH 19/19] save statistics for future --- lib/tests/fa_skeleton_test.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/lib/tests/fa_skeleton_test.py b/lib/tests/fa_skeleton_test.py index 66411a1..02866d2 100755 --- a/lib/tests/fa_skeleton_test.py +++ b/lib/tests/fa_skeleton_test.py @@ -180,12 +180,30 @@ def main(): # register and obtain *_InMNI_FA.nii.gz mniFAimgs= sub2tmp2mni(templatePath, siteName, faImgs, N_proc) + + + # save statistics for future + statFile= os.path.join(self.templatePath, 'meanFAstat.txt') + f= open(statFile,'a') + stdout= sys.stdout + sys.stdout= f + + print(datetime.now().strftime('%c'),'\n') # pass *_InMNI_FA.nii.gz list to analyzeStat site_means= analyzeStat(mniFAimgs) print(f'{siteName} site: ') printStat(site_means, mniFAimgs) + f.close() + sys.stdout= stdout + + # print statistics on console + print('') + with open(statFile) as f: + print(f.read()) + + print('\nThe statistics are also saved in ', statFile) if __name__ == '__main__': main()