diff --git a/designer2/tmi.py b/designer2/tmi.py index 24b1071a..5171bb41 100644 --- a/designer2/tmi.py +++ b/designer2/tmi.py @@ -5,7 +5,7 @@ from logging import StreamHandler, FileHandler from lib.designer_input_utils import get_input_info, convert_input_data, create_shell_table, assert_inputs -from lib.designer_fit_wrappers import refit_or_smooth, save_params +from lib.designer_fit_wrappers import refit_or_smooth, save_params, b0restore_slope, akc_out from lib.io import load_mrtrix from mrtrix3 import run, path @@ -122,6 +122,15 @@ def usage(cmdline): #pylint: disable=unused-variable dki_options.add_argument('-polyreg',action='store_true',help='polynomial regression based DKI estimation') dki_options.add_argument('-maxb', metavar=(''),help='maximum b-value for DKI fitting, default=3.') + #=======================================tmi black voxel====================================== + dki_options.add_argument('-akc_lowerlim', metavar=(''),help='akc lower threshold, default=-1') + dki_options.add_argument('-akc_upperlim', metavar=(''),help='akc upper threshold, default=10') + dki_options.add_argument('-b0restore', action='store_true',help='b0-restore dki outlier correction') + dki_options.add_argument('-kernal', metavar=(''),help='kernal/patch size for b0-restore correction, default=5') + dki_options.add_argument('-percentile', metavar=(''),help='percentile of patch to include, default=10') + dki_options.add_argument('-thresh_criteria', metavar=(''),help='outlier percent improvement iteration threshold, default 0.05') + #=======================================tmi black voxel====================================== + smi_options = cmdline.add_argument_group('tensor options for the TMI script') smi_options.add_argument('-SMI', action='store_true',help='Perform estimation of SMI (standard model of Diffusion in White Matter). Please use in conjunction with the -bshape, -echo_time, -sigma, and -compartments options.') smi_options.add_argument('-compartments', metavar=(''),help='SMI compartments (IAS, EAS, and FW), default=IAS,EAS') @@ -431,25 +440,257 @@ def execute(): #pylint: disable=unused-variable dt_poly_dti = dti.train_rotated_bayes_fit(dwi_dti, dt_dti, s0_dti, b_dti, mask, True) logger.info("Polyreg DTI fit completed.", extra={"dt_poly_dti_shape": dt_poly_dti.shape}) + + #=======================================tmi black voxel====================================== + if app.ARGS.b0restore: + from lib.mpunits import vectorize + import scipy.io as sio + + #DETECTING OUTLIERS + logger.info("Starting AKC outlier detection...") + dwd = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + mat = sio.loadmat(os.path.join(dwd,'constant','dirs256.mat')) + dir = mat['dirs'] + # mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) + # dir = mat['dir'] + + if not app.ARGS.akc_lowerlim: + akc_lowerlim=0 + else: + akc_lowerlim=int(app.ARGS.akc_lowerlim) + + if not app.ARGS.akc_upperlim: + akc_upperlim=10 + else: + akc_upperlim=int(app.ARGS.akc_upperlim) + + if not (app.ARGS.DKI or app.ARGS.WDKI): + logger.error("AKC Outlier detection must be accompanied by DKI option") + raise MRtrixError("AKC Outlier detection must be accompanied by DKI option") + else: + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir,akc_lowerlim, akc_upperlim) + akc_mask_copy = akc_mask.copy() + + logger.info("Outlier detection completed.", extra={"num_outliers": str(np.sum(akc_mask))}) + + outlier_mask = {} + akc_mask_tmp = vectorize(akc_mask, mask) + outlier_mask['akc'] = akc_mask_tmp + save_params(outlier_mask, mif, model='dki', outdir=outdir) + + x=np.shape(mask)[0] + y=np.shape(mask)[1] + z=np.shape(mask)[2] + akc_dirs=np.zeros((x,y,z,np.shape(akc_d)[0])) + # print(np.shape(akc_d)) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + + + # ==================================================================================== + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + # first akc outlier mask + print("============AKC mask============") + outlier_inds=np.array(np.where(mask>0)) + # n=np.shape(outlier_inds)[-1] + dwi_norm = abs(dwi_dki) / np.amax(dwi_dki, axis=(0,1,2)) + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + + outlier_mask = {} + akc_mask[akc_mask!=0]=1 + akc_mask=vectorize(akc_mask, mask) + akc_mask[akc_mask_new==1]=1 + outlier_mask['outlier_mask'] = akc_mask + save_params(outlier_mask, mif, model='dki', outdir=outdir) + akc_mask = akc_mask.astype(bool) + noutlier0=np.sum(akc_mask) + + # =======================================b0-restore======================================= + print("==============b0-restore==============") + if not app.ARGS.percentile: + percentile=10 + else: + percentile=int(app.ARGS.percentile) + + if not app.ARGS.kernal: + kernal=5 + else: + kernal=int(app.ARGS.kernal) + + #new dwi with restored b0 + dwi_new = b0restore_slope(akc_mask, dwi_dki, bval_dki, kernal,percentile,fa, md,mask=None, n_cores=-3) + + dt_new,s0_new,b_dki = dki.dki_fit(dwi_new, akc_mask) + dtishell = (bval_dki <= 0.1) | ((bval_dki > .5) & (bval_dki <= 1.5)) + + x,y,z = np.where(akc_mask == 1) + DT = vectorize(dt_dki, mask) + DT[x,y,z,:] = dt_new.T + dt_dki = vectorize(DT, mask) + + # Detect Outlier + print("============Detect Outlier now with conservative outlier detection (AKC < -2)============") + akc_lowerlim=-2 + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + akc_mask_copy = akc_mask.copy() + akc_mask = vectorize(akc_mask, mask) + akc_mask_tmp = akc_mask + + #extract new rk, md, fa + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + + #new akc outlier mask based on rk and md + xx=np.shape(mask)[0] + yy=np.shape(mask)[1] + zz=np.shape(mask)[2] + dwi_norm = abs(dwi_new) / np.amax(dwi_new, axis=(0,1,2)) + akc_dirs=np.zeros((xx,yy,zz,np.shape(akc_d)[0])) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + akc_mask[akc_mask!=0]=1 + akc_mask[akc_mask_new==1]=1 + # akc_dirs_mask= np.repeat(np.reshape(akc_mask,(xx,yy,zz,1)),np.shape(akc_d)[0],axis=3) + akc_mask = akc_mask.astype(bool) + + # # print('dir shape: {}'.format(np.shape(dir))) + # # print('dir type: {}'.format(type(dir))) + # # print('bvec shape: {}'.format(np.shape(np.reshape(bvec_dki,(-1,3))))) + # _,akc_d_temp = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + # akc_dirs_temp=np.zeros((xx,yy,zz,np.shape(akc_d_temp)[0])) + # for i in range(np.shape(akc_d_temp)[0]): + # akc_dirs_temp[:,:,:,i]=vectorize(akc_d_temp[i,:],mask) + # akc_dirs_mask[np.where(akc_dirs_mask>0) and np.where(akc_dirs_temp>1)]=0 + + noutlier=np.sum(akc_mask) + c=np.sum(akc_mask)+1 + count=0 + improve=(noutlier0-noutlier)/noutlier0 + print('Number of outliers: {}'.format(np.sum(mask))) + + # ==========================iteration=========================== + if app.ARGS.thresh_criteria: + thresh=float(app.ARGS.thresh_criteria) + else: + thresh=0.05 + + if improve>thresh: + print(f"{improve} > {thresh}") + + while True: + count=count+1 + print('iteration {}'.format(count)) + noutlier0=np.sum(akc_mask) + + # iterate b0-restore + print('Start correction {}'.format(count)) + dwi_new = b0restore_slope(akc_mask, dwi_new, bval_dki, kernal,percentile,fa,md,mask=None,n_cores=-3) + + # detect new outliers + dt_new,s0_new,b_dki = dki.dki_fit(dwi_new, akc_mask) + x,y,z = np.where(akc_mask == 1) + # DT = vectorize(dt_dki, mask) + DT[x,y,z,:] = dt_new.T + dt_dki = vectorize(DT, mask) + + print('detecting AKC outliers') + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + akc_mask_copy = akc_mask.copy() + + #akc mask + dwi_norm = abs(dwi_new) / np.amax(dwi_new, axis=(0,1,2)) + akc_dirs=np.zeros((xx,yy,zz,np.shape(akc_d)[0])) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + akc_mask = vectorize(akc_mask, mask).astype(bool) + akc_mask_tmp = akc_mask + + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + + print('detecting RK outliers') + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + akc_mask[akc_mask!=0]=1 + akc_mask[akc_mask_new==1]=1 + akc_mask = akc_mask.astype(bool) + noutlier=np.sum(akc_mask) + improve=(noutlier0-noutlier)/noutlier0 + print("nOutlier {}".format(np.sum(akc_mask))) + # print('{} > {}'.format(improve, thresh)) + + logger.info("Outlier correction iteration {}".format(count), extra={"num_outliers": str(np.sum(akc_mask))}) + if improve < thresh: + break + # ==========================iteration=========================== + + #save new b0 + newdwi = {} + b0_idx = np.where(bval_dki < 0.01)[0] + np.savetxt('{}/b0_indices.txt'.format(outdir),b0_idx,fmt='%d') + newdwi['b0'] = dwi_new[:,:,:,b0_idx] + save_params(newdwi, mif, model='dki', outdir=outdir) + + #save new dt + dt_ = {} + dt_['dt'] = DT + save_params(dt_, mif, model='dki_b0restore', outdir=outdir) + logger.info("DKT with b0-restore saved.") + + #save new outlier mask + outlier_mask = {} + akc_mask_copy[akc_mask_copy!=0]=1 + akc_mask_copy=vectorize(akc_mask_copy, mask) + akc_mask_copy[akc_mask_new==1]=1 + outlier_mask['outliermask_final'] = akc_mask_copy + save_params(outlier_mask, mif, model='dki', outdir=outdir) + + logger.info("AKC outlier post-processing completed.", extra={"num_outliers": str(np.sum(akc_mask))}) + else: + akc_mask = np.zeros_like(mask) + + #=======================================tmi black voxel====================================== + + if app.ARGS.akc_outliers: from lib.mpunits import vectorize import scipy.io as sio logger.info("Starting AKC outlier detection...") dwd = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) - mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) - dir = mat['dir'] + # mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) + mat = sio.loadmat(os.path.join(dwd,'constant','dirs256.mat')) + dir = mat['dirs'] if not (app.ARGS.DKI or app.ARGS.WDKI): logger.error("AKC Outlier detection must be accompanied by DKI option") raise MRtrixError("AKC Outlier detection must be accompanied by DKI option") else: - akc_mask = dki.outlierdetection(dt_dki, mask, dir) + if not app.ARGS.akc_lowerlim: + akc_lowerlim=0 + else: + akc_lowerlim=int(app.ARGS.akc_lowerlim) + + if not app.ARGS.akc_upperlim: + akc_upperlim=10 + else: + akc_upperlim=int(app.ARGS.akc_upperlim) + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir,akc_lowerlim,akc_upperlim) akc_mask = vectorize(akc_mask, mask).astype(bool) logger.info("Outlier detection completed.", extra={"num_outliers": str(np.sum(akc_mask))}) - dwi_new = refit_or_smooth(akc_mask, dwi_dki, n_cores=int(app.ARGS.n_cores)) + if app.ARGS.b0restore: + dwi_new = refit_or_smooth(akc_mask, dwi_new, n_cores=int(app.ARGS.n_cores)) + else: + dwi_new = refit_or_smooth(akc_mask, dwi_dki, n_cores=int(app.ARGS.n_cores)) dt_new,_,_ = dki.dki_fit(dwi_new, akc_mask) x,y,z = np.where(akc_mask == 1) @@ -462,7 +703,7 @@ def execute(): #pylint: disable=unused-variable logger.info("DKT with AKC saved.") dt_dki = vectorize(DT, mask) - akc_mask = dki.outlierdetection(dt_dki, mask, dir) + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) akc_mask = vectorize(akc_mask, mask).astype(bool) logger.info("AKC outlier post-processing completed.", extra={"num_outliers": str(np.sum(akc_mask))}) else: @@ -482,7 +723,10 @@ def execute(): #pylint: disable=unused-variable logger.info("DTI fit after smoothing completed.", extra={"dt_dti_shape": dt_dti.shape}) if (app.ARGS.DKI or app.ARGS.WDKI): - dwi_new = refit_or_smooth(akc_mask, dwi_dki, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) + if app.ARGS.b0restore: + dwi_new = refit_or_smooth(akc_mask, dwi_new, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) + else: + dwi_new = refit_or_smooth(akc_mask, dwi_dki, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) dt_dki,_,_ = dki.dki_fit(dwi_new, mask) DT = vectorize(dt_dki, mask) @@ -616,32 +860,263 @@ def execute(): #pylint: disable=unused-variable dt_poly_dti = dti.train_rotated_bayes_fit(dwi_dti, dt_dti, s0_dti, b_dti, mask, True) logger.info(f"Polyreg DTI fit completed for TE={te}.", extra={"dt_poly_dti_shape": dt_poly_dti.shape}) + #=======================================tmi black voxel====================================== + if app.ARGS.b0restore: + from lib.mpunits import vectorize + import scipy.io as sio + + #DETECTING OUTLIERS + logger.info("Starting AKC outlier detection...") + dwd = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + mat = sio.loadmat(os.path.join(dwd,'constant','dirs256.mat')) + dir = mat['dirs'] + # mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) + # dir = mat['dir'] + + if not app.ARGS.akc_lowerlim: + akc_lowerlim=0 + else: + akc_lowerlim=int(app.ARGS.akc_lowerlim) + + if not app.ARGS.akc_upperlim: + akc_upperlim=10 + else: + akc_upperlim=int(app.ARGS.akc_upperlim) + + if not (app.ARGS.DKI or app.ARGS.WDKI): + logger.error("AKC Outlier detection must be accompanied by DKI option") + raise MRtrixError("AKC Outlier detection must be accompanied by DKI option") + else: + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir,akc_lowerlim, akc_upperlim) + akc_mask_copy = akc_mask.copy() + + logger.info("Outlier detection completed.", extra={"num_outliers": str(np.sum(akc_mask))}) + + outlier_mask = {} + akc_mask_tmp = vectorize(akc_mask, mask) + outlier_mask['akc'] = akc_mask_tmp + save_params(outlier_mask, mif, model='dki', outdir=outdir) + + x=np.shape(mask)[0] + y=np.shape(mask)[1] + z=np.shape(mask)[2] + akc_dirs=np.zeros((x,y,z,np.shape(akc_d)[0])) + # print(np.shape(akc_d)) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + + + # ==================================================================================== + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + # first akc outlier mask + print("============AKC mask============") + outlier_inds=np.array(np.where(mask>0)) + # n=np.shape(outlier_inds)[-1] + dwi_norm = abs(dwi_dki) / np.amax(dwi_dki, axis=(0,1,2)) + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + + outlier_mask = {} + akc_mask[akc_mask!=0]=1 + akc_mask=vectorize(akc_mask, mask) + akc_mask[akc_mask_new==1]=1 + outlier_mask['outlier_mask'] = akc_mask + save_params(outlier_mask, mif, model='dki', outdir=outdir) + akc_mask = akc_mask.astype(bool) + noutlier0=np.sum(akc_mask) + + # =======================================b0-restore======================================= + print("==============b0-restore==============") + if not app.ARGS.percentile: + percentile=10 + else: + percentile=int(app.ARGS.percentile) + + if not app.ARGS.kernal: + kernal=5 + else: + kernal=int(app.ARGS.kernal) + + #new dwi with restored b0 + dwi_new = b0restore_slope(akc_mask, dwi_dki, bval_dki, kernal,percentile,fa, md,mask=None, n_cores=-3) + + dt_new,s0_new,b_dki = dki.dki_fit(dwi_new, akc_mask) + dtishell = (bval_dki <= 0.1) | ((bval_dki > .5) & (bval_dki <= 1.5)) + + x,y,z = np.where(akc_mask == 1) + DT = vectorize(dt_dki, mask) + DT[x,y,z,:] = dt_new.T + dt_dki = vectorize(DT, mask) + + # Detect Outlier + print("============Detect Outlier now with conservative AKC outlier detection (AKC < -2)============") + akc_lowerlim=-2 + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + akc_mask_copy = akc_mask.copy() + akc_mask = vectorize(akc_mask, mask) + akc_mask_tmp = akc_mask + + #extract new rk, md, fa + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + + #new akc outlier mask based on rk and md + xx=np.shape(mask)[0] + yy=np.shape(mask)[1] + zz=np.shape(mask)[2] + dwi_norm = abs(dwi_new) / np.amax(dwi_new, axis=(0,1,2)) + akc_dirs=np.zeros((xx,yy,zz,np.shape(akc_d)[0])) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + akc_mask[akc_mask!=0]=1 + akc_mask[akc_mask_new==1]=1 + # akc_dirs_mask= np.repeat(np.reshape(akc_mask,(xx,yy,zz,1)),np.shape(akc_d)[0],axis=3) + akc_mask = akc_mask.astype(bool) + + # # print('dir shape: {}'.format(np.shape(dir))) + # # print('dir type: {}'.format(type(dir))) + # # print('bvec shape: {}'.format(np.shape(np.reshape(bvec_dki,(-1,3))))) + # _,akc_d_temp = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + # akc_dirs_temp=np.zeros((xx,yy,zz,np.shape(akc_d_temp)[0])) + # for i in range(np.shape(akc_d_temp)[0]): + # akc_dirs_temp[:,:,:,i]=vectorize(akc_d_temp[i,:],mask) + # akc_dirs_mask[np.where(akc_dirs_mask>0) and np.where(akc_dirs_temp>1)]=0 + + noutlier=np.sum(akc_mask) + c=np.sum(akc_mask)+1 + count=0 + improve=(noutlier0-noutlier)/noutlier0 + print('Number of outliers: {}'.format(np.sum(mask))) + + # # ====================================iteration=================================== + # iteration + if app.ARGS.thresh_criteria: + thresh=float(app.ARGS.thresh_criteria) + else: + thresh=0.05 + + if improve>thresh: + print('{} > {}'.format(improve, thresh)) + + while True: + count=count+1 + print('iteration {}'.format(count)) + noutlier0=np.sum(akc_mask) + + # iterate b0-restore + print('Start correction {}'.format(count)) + dwi_new = b0restore_slope(akc_mask, dwi_new, bval_dki, kernal,percentile,fa,md,mask=None,n_cores=-3) + + # detect new outliers + dt_new,s0_new,b_dki = dki.dki_fit(dwi_new, akc_mask) + x,y,z = np.where(akc_mask == 1) + # DT = vectorize(dt_dki, mask) + DT[x,y,z,:] = dt_new.T + dt_dki = vectorize(DT, mask) + print('detecting outliers') + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) + akc_mask_copy = akc_mask.copy() + + #akc mask + dwi_norm = abs(dwi_new) / np.amax(dwi_new, axis=(0,1,2)) + akc_dirs=np.zeros((xx,yy,zz,np.shape(akc_d)[0])) + for i in range(np.shape(akc_d)[0]): + akc_dirs[:,:,:,i]=vectorize(akc_d[i,:],mask) + akc_mask = vectorize(akc_mask, mask).astype(bool) + akc_mask_tmp = akc_mask + params_dki = dki.extract_parameters(dt_dki, b_dki, mask, extract_dti=True, extract_dki=True, fit_w=False) + rk = params_dki['rk'] + md = params_dki['md'] + fa = params_dki['fa'] + print('detecting outliers') + akc_mask_new=akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3) + akc_mask[akc_mask!=0]=1 + akc_mask[akc_mask_new==1]=1 + akc_mask = akc_mask.astype(bool) + noutlier=np.sum(akc_mask) + improve=(noutlier0-noutlier)/noutlier0 + print("nOutlier {}".format(np.sum(akc_mask))) + # print('{} > {}'.format(improve, thresh)) + + logger.info("Outlier correction iteration {} for TE={}".format(count,te), extra={"num_outliers": str(np.sum(akc_mask))}) + if improve < thresh: + break + # ====================================iteration=================================== + + #save new b0 + newdwi = {} + b0_idx = np.where(bval_dki < 0.01)[0] + np.savetxt('{}/b0_indices.txt'.format(outdir),b0_idx,fmt='%d') + newdwi['b0'] = dwi_new[:,:,:,b0_idx] + save_params(newdwi, mif, model='dki', outdir=outdir) + + #save new dt + dt_ = {} + dt_['dt'] = DT + save_params(dt_, mif, model='dki_b0restore', outdir=outdir) + logger.info("DKT with b0-restore saved.") + + #save new outlier mask + outlier_mask = {} + akc_mask_copy[akc_mask_copy!=0]=1 + akc_mask_copy=vectorize(akc_mask_copy, mask) + akc_mask_copy[akc_mask_new==1]=1 + outlier_mask['outliermask_final'] = akc_mask_copy + save_params(outlier_mask, mif, model='dki', outdir=outdir) + + logger.info("b0-restore AKC outlier post-processing completed for TE={}.".format(te), extra={"num_outliers": str(np.sum(akc_mask))}) + else: + akc_mask = np.zeros_like(mask) + + #=======================================tmi black voxel====================================== + + + if app.ARGS.akc_outliers: from lib.mpunits import vectorize import scipy.io as sio logger.info(f"Starting AKC outlier detection for TE={te}...") dwd = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) - mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) - dir = mat['dir'] + # mat = sio.loadmat(os.path.join(dwd,'constant','dirs10000.mat')) + mat = sio.loadmat(os.path.join(dwd,'constant','dirs256.mat')) + dir = mat['dirs'] + + if not app.ARGS.akc_lowerlim: + akc_lowerlim=0 + else: + akc_lowerlim=int(app.ARGS.akc_lowerlim) + + if not app.ARGS.akc_upperlim: + akc_upperlim=10 + else: + akc_upperlim=int(app.ARGS.akc_upperlim) if not (app.ARGS.DKI or app.ARGS.WDKI): logger.error(f"AKC Outlier detection for TE={te} must be accompanied by DKI option.") raise MRtrixError("AKC Outlier detection must be accompanied by DKI option") else: - akc_mask = dki.outlierdetection(dt_dki, mask, dir) + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) akc_mask = vectorize(akc_mask, mask).astype(bool) logger.info(f"Outlier detection completed for TE={te}.", extra={"num_outliers": str(np.sum(akc_mask))}) - dwi_new = refit_or_smooth(akc_mask, dwi_dki, n_cores=int(app.ARGS.n_cores)) + if app.ARGS.b0restore: + dwi_new = refit_or_smooth(akc_mask, dwi_new, n_cores=int(app.ARGS.n_cores)) + else: + dwi_new = refit_or_smooth(akc_mask, dwi_dki, n_cores=int(app.ARGS.n_cores)) dt_new,_,_ = dki.dki_fit(dwi_new, akc_mask) x,y,z = np.where(akc_mask == 1) DT = vectorize(dt_dki, mask) DT[x,y,z,:] = dt_new.T dt_dki = vectorize(DT, mask) - akc_mask = dki.outlierdetection(dt_dki, mask, dir) + akc_mask,akc_d = dki.outlierdetection(dt_dki, mask, dir, akc_lowerlim,akc_upperlim) akc_mask = vectorize(akc_mask, mask).astype(bool) logger.info(f"AKC outlier post-processing completed for TE={te}.", extra={"num_outliers": str(np.sum(akc_mask))}) else: @@ -654,7 +1129,10 @@ def execute(): #pylint: disable=unused-variable dt_dti,_,_ = dti.dti_fit(dwi_new, mask) logger.info(f"DTI fit after smoothing completed for TE={te}.", extra={"dt_dti_shape": dt_dti.shape}) if (app.ARGS.DKI or app.ARGS.WDKI): - dwi_new = refit_or_smooth(akc_mask, dwi_dki, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) + if app.ARGS.b0restore: + dwi_new = refit_or_smooth(akc_mask, dwi_new, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) + else: + dwi_new = refit_or_smooth(akc_mask, dwi_dki, mask=mask, smoothlevel=int(app.ARGS.fit_smoothing)) dt_dki,_,_ = dki.dki_fit(dwi_new, mask) logger.info(f"DKI fit after smoothing completed for TE={te}.", extra={"dt_dki_shape": dt_dki.shape}) diff --git a/docs/docs/TMI/usage.md b/docs/docs/TMI/usage.md index 81dd336b..5858dfd4 100644 --- a/docs/docs/TMI/usage.md +++ b/docs/docs/TMI/usage.md @@ -15,7 +15,7 @@ Main usage: `tmi ` ### `input` -Input to `tmi` can be any diffusion MRI image file that is compatible with mrtrix3 as an input. Ideally the input to `tmi` is the output from the `deisgner` preprocessing pipeline. +Input to `tmi` can be any diffusion MRI image file that is compatible with mrtrix3 as an input. Ideally the input to `tmi` is the output from the `designer` preprocessing pipeline. ### `output` Name of the folder which will contain parameter images. @@ -55,8 +55,15 @@ By default, if none of the below options are used, TMI will not estimate paramet ## Options for outlier replacement ### `-akc_outliers` -- Brute force K tensor outlier detection and filtering. The kurtosis tensor is projected onto a sphere with 10,000 directions. Voxels where mean kurtosis is less than -1 or greater than 10 are labelled as outliers. +- Brute force K tensor outlier detection and filtering. The kurtosis tensor is projected onto a sphere with 256 directions. Voxels where mean kurtosis is less than -1 or greater than 10 are labelled as outliers. +- Use `-akc_lowerlim ` and `akc_upperlim ` to set your own lower and upper AKC limit respectively. - Outliers are replaced by the median value of neighboring (26 adjacent) non-outlier voxels in all diffusion weighted images and the DKI fit is done again. + ### `-fit_smoothing ` - Windowed adaptive nonlocal means filter on the dwi. diff --git a/lib/designer_fit_wrappers.py b/lib/designer_fit_wrappers.py index 9f4c5763..e2e6f6b3 100644 --- a/lib/designer_fit_wrappers.py +++ b/lib/designer_fit_wrappers.py @@ -76,6 +76,174 @@ def refit_or_smooth(outlier_locations, dwi, mask=None, smoothlevel=None, n_cores return dwi_new +#=======================================tmi black voxel====================================== +def parallel_outlier_akc(inds, akc_mask_tmp, akc_dirs,rk,md,fa): + import numpy as np + k = 7 // 2 + x = inds[0] + y = inds[1] + z = inds[2] + O=akc_mask_tmp[x,y,z] + val_rk=rk[x,y,z] + val_md=md[x,y,z] + val_fa=fa[x,y,z] + max_akc=np.max(akc_dirs[x,y,z,:]) + min_akc=np.min(akc_dirs[x,y,z,:]) + # if max_akc>10 or min_akc<0: + # O=1 + + if val_rk<0.5 and val_md<2: + O=1 + return int(O) + +def akc_out(outlier_inds, akc_mask_tmp, akc_dirs,rk,md,fa, n_cores=-3): + from joblib import Parallel, delayed + import numpy as np + + new_akc_mask=akc_mask_tmp.copy() + + wval = (Parallel(n_jobs=n_cores, prefer='processes') + (delayed(parallel_outlier_akc)( + outlier_inds[:,i], akc_mask_tmp, akc_dirs,rk,md,fa + ) for i in range(len(outlier_inds[0])))) + + wval=np.asarray(wval) + wval=wval.astype(int) + new_akc_mask[outlier_inds[0,:],outlier_inds[1,:],outlier_inds[2,:]] = wval[:] + + return new_akc_mask + + +def parallel_outlier_slope(inds, kernel, outlier_locations, bval, dwi_norm, dwi, fa,md, md_mask, smoothlevel): + import numpy as np + import warnings + + warnings.filterwarnings("ignore", + message="Mean of empty slice", + category=RuntimeWarning) + + x, y, z = inds + k = kernel // 2 + # --- precomputed shells --- + bval_rounded = np.round(bval, 2) + nonb0shell = bval > 0.05 + lowbval = np.sort(np.unique(bval_rounded[bval_rounded < 1.2])) + bx = lowbval[-1] + bxx = lowbval[0] + bxshell = bval_rounded == bx + bxxshell = bval_rounded == bxx + + # --- grow patch --- + while True: + xmin = max(x - k - 1, 0) + xmax = min(x + k, dwi.shape[0]) + ymin = max(y - k - 1, 0) + ymax = min(y + k, dwi.shape[1]) + zmin = max(z - k - 1, 0) + zmax = min(z + k, dwi.shape[2]) + + akcpatch = outlier_locations[xmin:xmax, ymin:ymax, zmin:zmax] + + if not akcpatch.all(): + break + k += 2 + + # --- flatten patches --- + psize = akcpatch.size + + fapatch = fa[xmin:xmax, ymin:ymax, zmin:zmax].ravel() + mdpatch = md[xmin:xmax, ymin:ymax, zmin:zmax].ravel() + csfpatch = md_mask[xmin:xmax, ymin:ymax, zmin:zmax].ravel() + + # --- similarity intensities --- + patch = dwi_norm[xmin:xmax, ymin:ymax, zmin:zmax, nonb0shell] + patch = patch.reshape(psize, -1) + + ref = dwi_norm[x, y, z, nonb0shell][None, :] + diff = patch - ref + intensities = np.sqrt((diff * diff).sum(axis=1)) / patch.shape[1] + + # --- FA / MD rejection mask --- + omit = ( + (fapatch > fa[x,y,z] + fapatch.std()) | + (fapatch < fa[x,y,z] - fapatch.std()) | + (mdpatch > md[x,y,z] + mdpatch.std()) | + (mdpatch < md[x,y,z] - mdpatch.std()) + ) + + # --- rank + exclude --- + min_wgs = intensities.copy() + wgs_max = min_wgs.max() + min_wgs[akcpatch.ravel()] = wgs_max + min_wgs[csfpatch.ravel()] = wgs_max + min_wgs[omit] = wgs_max + + if not smoothlevel: + goodidx = min_wgs <= min_wgs.mean() + else: + thr = np.percentile(min_wgs, smoothlevel) + goodidx = (min_wgs <= thr) & (min_wgs != wgs_max) + + # --- slopes --- + bx_bv = np.log(dwi[x,y,z,bxshell].mean()) + bxx_bv = np.log(dwi[x,y,z,bxxshell]) + + slope_bv = np.abs((bx_bv - bxx_bv) / (bx - bxx)) + + patch_bx = np.log(dwi[xmin:xmax,ymin:ymax,zmin:zmax][:,:,:,bxshell]).reshape(psize, -1) + patch_bxx = np.log(dwi[xmin:xmax,ymin:ymax,zmin:zmax][:,:,:,bxxshell]).reshape(psize, -1) + + bx_ok = patch_bx[goodidx].mean(axis=1) + bxx_ok = patch_bxx[goodidx] + + slope_ok = np.abs((bx_ok[:,None] - bxx_ok) / (bx - bxx)) + + mask = slope_ok > slope_bv + masked = np.where(mask, slope_ok, np.nan) + mean_vals = np.nanmean(masked, axis=0) + slope_ok2 = np.where(np.isnan(mean_vals), slope_bv, mean_vals) + + wval_log = np.abs(slope_ok2) * bx + bx_bv + return np.exp(wval_log).squeeze() + +def b0restore_slope(outlier_locations, dwi, bval,k,perc,fa,md, mask=None, n_cores=-3): + from joblib import Parallel, delayed + import numpy as np + from tqdm import tqdm + + if mask is None: + outinds = np.array(np.where(outlier_locations == 1)) + else: + outinds = np.array(np.where(mask == 1)) + + dwi_norm = abs(dwi) / np.amax(dwi, axis=(0,1,2)) + dwi_new = dwi.copy() + kernel = k + + csf_t=1.2 + md_mask = md.copy() + md_mask[md>=csf_t] = int(1) + md_mask[md 10), axis=0) + return np.any((akc < akc_lowerlim) | (akc > akc_uplim), axis=0) - def outlierdetection(self, dt, mask, dir): + def akc_dirs(self, dt, dir): + akc = self.kurtosis_coeff(dt, dir) + # print(np.shape(akc)) + # print(np.sum(np.any((akc < akc_lim) | (akc > 10), axis=0))) + return akc + + def outlierdetection(self, dt, mask, dir, akc_lowerlim, akc_uplim): nvxls = dt.shape[1] akc_mask = np.zeros((nvxls)) - nblocks = 200 + # nblocks = 200 + nblocks = 4 + + N=np.shape(dir)[0] try: - N = 10000 + N=np.shape(dir)[0] + # N = 10000 akc_mask = Parallel(n_jobs=self.n_cores, prefer='threads') \ (delayed(self.compute_outliers) - (dt, dir[int(N / nblocks * (i - 1)):int(N / nblocks * i), :]) for i in range(1, nblocks + 1) + (dt, dir[int(N / nblocks * (i - 1)):int(N / nblocks * i), :], akc_lowerlim, akc_uplim) for i in range(1, nblocks + 1) ) except: N = 1000 akc_mask = Parallel(n_jobs=8) \ (delayed(self.compute_outliers) - (dt, dir[int(N / nblocks * (i - 1)):int(N / nblocks * i), :]) for i in range(1, nblocks + 1) + (dt, dir[int(N / nblocks * (i - 1)):int(N / nblocks * i), :], akc_lowerlim, akc_uplim) for i in range(1, nblocks + 1) ) - return np.sum(akc_mask, axis=0) - + try: + N=np.shape(dir)[0] + # N=10000 + akc_d = Parallel(n_jobs=self.n_cores, prefer='threads')\ + (delayed(self.akc_dirs) + (dt, dir[int(N/nblocks*(i-1)):int(N/nblocks*i),:]) for i in range(1, nblocks + 1) + ) + except: + N = 1000 + akc_d = Parallel(n_jobs=8)\ + (delayed(self.akc_dirs) + (dt, dir[int(N/nblocks*(i-1)):int(N/nblocks*i),:]) for i in range(1, nblocks + 1) + ) + # print(np.shape(akc_d)) + akc_d=np.asarray(akc_d) + akc_d = np.reshape(akc_d,(-1,np.sum(mask))) + + return np.sum(akc_mask, axis=0),akc_d + #==================================black voxel=============================== def identify_outliers(self, dt, percentiles): low_thresh = np.percentile(dt, percentiles[0], method='median_unbiased')