diff --git a/designer2/designer.py b/designer2/designer.py index ee828f6..7cc5cb3 100644 --- a/designer2/designer.py +++ b/designer2/designer.py @@ -126,8 +126,6 @@ def execute(): #pylint: disable=unused-variable index = shell_rows) print('input DWI data has properties:') print(shell_df) - - import pdb; pdb.set_trace() app.goto_scratch_dir() diff --git a/designer2/tmi.py b/designer2/tmi.py index 4ad2011..a1d2658 100644 --- a/designer2/tmi.py +++ b/designer2/tmi.py @@ -109,7 +109,7 @@ def usage(cmdline): #pylint: disable=unused-variable options.add_argument('-n_cores',metavar=(''),help='specify the number of cores to use in parallel tasts, by default designer will use available cores - 2', default=-3) options.add_argument('-echo_time',metavar=(''),help='specify the echo time used in the acquisition (comma separated list the same length as number of inputs)') options.add_argument('-bshape',metavar=(''),help='specify the b-shape used in the acquisition (comma separated list the same length as number of inputs)') - + dki_options = cmdline.add_argument_group('tensor options for the TMI script') dki_options.add_argument('-DKI', action='store_true', help='Include DKI parameters in output folder (mk,ak,rk)') dki_options.add_argument('-DTI', action='store_true', help='Include DTI parameters in output folder (md,ad,rd,fa,eigenvalues, eigenvectors') @@ -581,13 +581,19 @@ def execute(): #pylint: disable=unused-variable logger.info("No lmax specified for SMI. Using default.") logger.info("Initializing SMI fitting...") + echo_times = dwi_metadata['echo_time_per_volume'] + if np.max(echo_times) > 1.0: + logger.info("Echo times in ms, converting to seconds.") + echo_times /= 1000 + if multi_te_beta: smi = SMI(bval=bval_orig, bvec=bvec_orig, rotinv_lmax=lmax) smi.set_compartments(compartments) - smi.set_echotime(dwi_metadata['echo_time_per_volume']) + + smi.set_echotime(echo_times) smi.set_bshape(dwi_metadata['bshape_per_volume']) logger.info("SMI model initialized for multi-TE/beta data.") - + params_smi = smi.fit(dwi_orig, mask=mask, sigma=sigma) logger.info("SMI fitting completed for multi-TE/beta data.", extra={"params_smi_shape": {key: value.shape for key, value in params_smi.items()}}) @@ -596,7 +602,7 @@ def execute(): #pylint: disable=unused-variable else: smi = SMI(bval=bval, bvec=bvec, rotinv_lmax=lmax) smi.set_compartments(compartments) - smi.set_echotime(dwi_metadata['echo_time_per_volume']) + smi.set_echotime(echo_times) smi.set_bshape(dwi_metadata['bshape_per_volume']) logger.info("SMI model initialized for single-TE/beta data.") diff --git a/lib/designer_input_utils.py b/lib/designer_input_utils.py index fa7d101..135a94e 100644 --- a/lib/designer_input_utils.py +++ b/lib/designer_input_utils.py @@ -172,6 +172,7 @@ def assert_inputs(dwi_metadata, args_pe_dir, args_pf): import json from mrtrix3 import app, MRtrixError import numpy as np + bidslist = dwi_metadata['bidslist'] # if bids files are provided, check that the user inputs match the bids files @@ -283,6 +284,9 @@ def convert_input_data(dwi_metadata): import numpy as np import os import warnings + import inspect + + caller = os.path.basename(inspect.stack()[-1].filename) miflist = [] phaselist = [] @@ -337,51 +341,53 @@ def convert_input_data(dwi_metadata): dwi_header = image.Header('%s/dwi.mif' % (app.SCRATCH_DIR)) dwi_ind_size.append([ int(s) for s in dwi_header.size() ]) - if app.ARGS.pf: - pf_per_series_app = convert_to_float(app.ARGS.pf) - try: - pf_per_series_bids = dwi_header.keyval()['PartialFourier'] - if pf_per_series_app != pf_per_series_bids: - warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value') - except: - pass - pf_per_series.append(pf_per_series_app) - else: - try: - pf_per_series.append(dwi_header.keyval()['PartialFourier']) - except: - warnings.warn('No partial fourier factor found in header, assuming full sampling') - pf_per_series.append(1) + if caller == 'designer': + if app.ARGS.pf: + pf_per_series_app = convert_to_float(app.ARGS.pf) + try: + pf_per_series_bids = dwi_header.keyval()['PartialFourier'] + if pf_per_series_app != pf_per_series_bids: + warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value') + except: + pass + pf_per_series.append(pf_per_series_app) + else: + try: + pf_per_series.append(dwi_header.keyval()['PartialFourier']) + except: + warnings.warn('No partial fourier factor found in header, assuming full sampling') + pf_per_series.append(1) - if app.ARGS.pe_dir: - pe_dir_app = app.ARGS.pe_dir - pe_dir = convert_pe_dir_to_ijk(pe_dir_app) - try: - ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection'] - if ped_per_series_bids != pe_dir_app: - warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value') - except: - pass - ped_per_series.append(pe_dir) - else: - try: - ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection']) - except: - raise MRtrixError('No phase encoding direction found in header, please specify manually') + if app.ARGS.pe_dir: + pe_dir_app = app.ARGS.pe_dir + pe_dir = convert_pe_dir_to_ijk(pe_dir_app) + try: + ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection'] + if ped_per_series_bids != pe_dir_app: + warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value') + except: + pass + ped_per_series.append(pe_dir) + else: + try: + ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection']) + except: + raise MRtrixError('No phase encoding direction found in header, please specify manually') if app.ARGS.echo_time: - te_per_series.append(app.ARGS.echo_time) + te_app = np.round(float(app.ARGS.echo_time), 3) + te_per_series.append(te_app) try: - te_per_series_bids = dwi_header.keyval()['EchoTime'] - if te_per_series_bids != app.ARGS.echo_time: + te_per_series_bids = np.round(float(dwi_header.keyval()['EchoTime']), 3) + if te_per_series_bids != te_app: warnings.warn('User defined echo time does not match that found in bids json. Using user defined value') except: pass else: try: - te_per_series.append(dwi_header.keyval()['EchoTime']) + te_per_series.append(np.round(float(dwi_header.keyval()['EchoTime']), 3)) except: - warnings.warn('No echo time found in header, assuming 0') + warnings.warn('No echo time found in header, assuming 0 unless specifified by a .echotime file') te_per_series.append(0) else: @@ -417,51 +423,53 @@ def convert_input_data(dwi_metadata): dwi_ind_size.append([ int(s) for s in dwi_header.size() ]) miflist.append('%s/dwi%s.mif' % (app.SCRATCH_DIR, str(idx))) - if app.ARGS.pf: - pf_per_series_app = convert_to_float(app.ARGS.pf) - try: - pf_per_series_bids = dwi_header.keyval()['PartialFourier'] - if pf_per_series_app != pf_per_series_bids: - warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value') - except: - pass - pf_per_series.append(pf_per_series_app) - else: - try: - pf_per_series.append(dwi_header.keyval()['PartialFourier']) - except: - warnings.warn('No partial fourier factor found in header, assuming full sampling') - pf_per_series.append(1) - - if app.ARGS.pe_dir: - pe_dir_app = app.ARGS.pe_dir - pe_dir = convert_pe_dir_to_ijk(pe_dir_app) - try: - ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection'] - if ped_per_series_bids != pe_dir_app: - warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value') - except: - pass - ped_per_series.append(pe_dir) - else: - try: - ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection']) - except: - raise MRtrixError('No phase encoding direction found in header, please specify manually') + if caller == 'designer': + if app.ARGS.pf: + pf_per_series_app = convert_to_float(app.ARGS.pf) + try: + pf_per_series_bids = dwi_header.keyval()['PartialFourier'] + if pf_per_series_app != pf_per_series_bids: + warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value') + except: + pass + pf_per_series.append(pf_per_series_app) + else: + try: + pf_per_series.append(dwi_header.keyval()['PartialFourier']) + except: + warnings.warn('No partial fourier factor found in header, assuming full sampling') + pf_per_series.append(1) + + if app.ARGS.pe_dir: + pe_dir_app = app.ARGS.pe_dir + pe_dir = convert_pe_dir_to_ijk(pe_dir_app) + try: + ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection'] + if ped_per_series_bids != pe_dir_app: + warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value') + except: + pass + ped_per_series.append(pe_dir) + else: + try: + ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection']) + except: + raise MRtrixError('No phase encoding direction found in header, please specify manually') if app.ARGS.echo_time: - te_per_series.append(app.ARGS.echo_time) + te_app = [np.round(float(i), 3) for i in app.ARGS.echo_time.rsplit(',')] + te_per_series.append(te_app[idx]) try: - te_per_series_bids = dwi_header.keyval()['EchoTime'] - if te_per_series_bids != app.ARGS.echo_time: + te_per_series_bids = np.round(float(dwi_header.keyval()['EchoTime']), 3) + if te_per_series_bids != te_app[idx]: warnings.warn('User defined echo time does not match that found in bids json. Using user defined value') except: pass else: try: - te_per_series.append(dwi_header.keyval()['EchoTime']) + te_per_series.append(np.round(float(dwi_header.keyval()['EchoTime']), 3)) except: - warnings.warn('No echo time found in header, assuming 0') + warnings.warn('No echo time found in header, assuming 0 unless specified by a .echotime file') te_per_series.append(0) DWImif = ' '.join(miflist) @@ -476,7 +484,7 @@ def convert_input_data(dwi_metadata): else: for idx,i in enumerate(phase_n_list): run.command('mrconvert %s%s %s/phase%s.nii' % - (i, dwi_ext[idx], app.SCRATCH_DIR, str(idx))) + (i, phase_ext[idx], app.SCRATCH_DIR, str(idx))) phaselist.append('%s/phase%s.nii' % (app.SCRATCH_DIR, str(idx))) run.command('mrcat -axis 3 %s %s/phase.nii' % (' '.join(phaselist), app.SCRATCH_DIR)) except: @@ -494,11 +502,8 @@ def convert_input_data(dwi_metadata): dwi_metadata['bshape'] = bshape bshape_per_series = dwi_metadata['bshape'] - try: - if (len(set(te_per_series)) > 1) and (not app.ARGS.rpe_te): - raise MRtrixError('If data has variable echo time and no RPE TE is specified, please use the -rpe_te flag to specify the RPE TE') - except: - pass + if (len(set(te_per_series)) > 1) and (not app.ARGS.rpe_te): + raise MRtrixError('If data has variable echo time and no RPE TE is specified, please use the -rpe_te flag to specify the RPE TE') if not all(x == ped_per_series[0] for x in ped_per_series): raise MRtrixError('input series have different phase encoding directions, series should be processed separately') @@ -524,8 +529,6 @@ def convert_input_data(dwi_metadata): dwi_ind_size = [i + [1] if len(i)==3 else i for i in dwi_ind_size] - import pdb; pdb.set_trace() - nvols = [i[3] for i in dwi_ind_size] for idx,i in enumerate(dwi_n_list): if len(dwi_n_list) == 1: diff --git a/lib/smi.py b/lib/smi.py index 4b1b940..3388880 100644 --- a/lib/smi.py +++ b/lib/smi.py @@ -860,7 +860,7 @@ def generate_sm_wfw_b_beta_te_ws0_training_data(self): if not self.flag_compartments[2]: f_fw = np.zeros((n_training, 1)) else: - f_fw = self.prior[:,[5]] + f_fw = self.prior[:,[4]] if not self.flag_compartments[1]: f = 1 - f_fw @@ -890,7 +890,6 @@ def generate_sm_wfw_b_beta_te_ws0_training_data(self): kernel_params[:,[8]] * abs(k2_all), kernel_params[:,[9]] * abs(k4_all), kernel_params[:,[10]] * abs(k6_all))) - return rot_invs def standard_model_mlfit_rot_invs(self, rot_invs, sigma_norm_limits): diff --git a/rpg_cpp/fftw-3.3.10/api/.libs/libapi.a b/rpg_cpp/fftw-3.3.10/api/.libs/libapi.a index 063967f..fc2ccb4 100644 Binary files a/rpg_cpp/fftw-3.3.10/api/.libs/libapi.a and b/rpg_cpp/fftw-3.3.10/api/.libs/libapi.a differ diff --git a/rpg_cpp/fftw-3.3.10/dft/.libs/libdft.a b/rpg_cpp/fftw-3.3.10/dft/.libs/libdft.a index 7f89829..dcaf04c 100644 Binary files a/rpg_cpp/fftw-3.3.10/dft/.libs/libdft.a and b/rpg_cpp/fftw-3.3.10/dft/.libs/libdft.a differ diff --git a/rpg_cpp/fftw-3.3.10/dft/scalar/.libs/libdft_scalar.a b/rpg_cpp/fftw-3.3.10/dft/scalar/.libs/libdft_scalar.a index da8ead8..79304f8 100644 Binary files a/rpg_cpp/fftw-3.3.10/dft/scalar/.libs/libdft_scalar.a and b/rpg_cpp/fftw-3.3.10/dft/scalar/.libs/libdft_scalar.a differ diff --git a/rpg_cpp/fftw-3.3.10/dft/scalar/codelets/.libs/libdft_scalar_codelets.a b/rpg_cpp/fftw-3.3.10/dft/scalar/codelets/.libs/libdft_scalar_codelets.a index 75f8a1c..ff5bf56 100644 Binary files a/rpg_cpp/fftw-3.3.10/dft/scalar/codelets/.libs/libdft_scalar_codelets.a and b/rpg_cpp/fftw-3.3.10/dft/scalar/codelets/.libs/libdft_scalar_codelets.a differ diff --git a/rpg_cpp/fftw-3.3.10/kernel/.libs/libkernel.a b/rpg_cpp/fftw-3.3.10/kernel/.libs/libkernel.a index 6614e80..980f136 100644 Binary files a/rpg_cpp/fftw-3.3.10/kernel/.libs/libkernel.a and b/rpg_cpp/fftw-3.3.10/kernel/.libs/libkernel.a differ diff --git a/rpg_cpp/fftw-3.3.10/libbench2/libbench2.a b/rpg_cpp/fftw-3.3.10/libbench2/libbench2.a index c453cd3..544ca48 100644 Binary files a/rpg_cpp/fftw-3.3.10/libbench2/libbench2.a and b/rpg_cpp/fftw-3.3.10/libbench2/libbench2.a differ diff --git a/rpg_cpp/fftw-3.3.10/rdft/.libs/librdft.a b/rpg_cpp/fftw-3.3.10/rdft/.libs/librdft.a index d1680e7..8bbea7c 100644 Binary files a/rpg_cpp/fftw-3.3.10/rdft/.libs/librdft.a and b/rpg_cpp/fftw-3.3.10/rdft/.libs/librdft.a differ diff --git a/rpg_cpp/fftw-3.3.10/rdft/scalar/.libs/librdft_scalar.a b/rpg_cpp/fftw-3.3.10/rdft/scalar/.libs/librdft_scalar.a index 787ffae..fa9f1b9 100644 Binary files a/rpg_cpp/fftw-3.3.10/rdft/scalar/.libs/librdft_scalar.a and b/rpg_cpp/fftw-3.3.10/rdft/scalar/.libs/librdft_scalar.a differ diff --git a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cb/.libs/librdft_scalar_r2cb.a b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cb/.libs/librdft_scalar_r2cb.a index 31a70d9..80f65b9 100644 Binary files a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cb/.libs/librdft_scalar_r2cb.a and b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cb/.libs/librdft_scalar_r2cb.a differ diff --git a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cf/.libs/librdft_scalar_r2cf.a b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cf/.libs/librdft_scalar_r2cf.a index ddf8a3c..c76dc5e 100644 Binary files a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cf/.libs/librdft_scalar_r2cf.a and b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2cf/.libs/librdft_scalar_r2cf.a differ diff --git a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2r/.libs/librdft_scalar_r2r.a b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2r/.libs/librdft_scalar_r2r.a index 6523f72..8cf04e0 100644 Binary files a/rpg_cpp/fftw-3.3.10/rdft/scalar/r2r/.libs/librdft_scalar_r2r.a and b/rpg_cpp/fftw-3.3.10/rdft/scalar/r2r/.libs/librdft_scalar_r2r.a differ diff --git a/rpg_cpp/fftw-3.3.10/reodft/.libs/libreodft.a b/rpg_cpp/fftw-3.3.10/reodft/.libs/libreodft.a index 63b040b..0e920ad 100644 Binary files a/rpg_cpp/fftw-3.3.10/reodft/.libs/libreodft.a and b/rpg_cpp/fftw-3.3.10/reodft/.libs/libreodft.a differ diff --git a/rpg_cpp/fftw-3.3.10/simd-support/.libs/libsimd_support.a b/rpg_cpp/fftw-3.3.10/simd-support/.libs/libsimd_support.a index fe1e1e0..9a3663e 100644 Binary files a/rpg_cpp/fftw-3.3.10/simd-support/.libs/libsimd_support.a and b/rpg_cpp/fftw-3.3.10/simd-support/.libs/libsimd_support.a differ diff --git a/setup.py b/setup.py index 91234ff..91cc406 100644 --- a/setup.py +++ b/setup.py @@ -76,7 +76,7 @@ def __str__(self): setup( name ='designer2', - version ='2.0.8', + version ='2.0.9', author ='Benjamin Ades-Aron', author_email ='benjamin.ades-aron@nyulangone.org', url ='https://github.com/NYU-DiffusionMRI/DESIGNER-v2',