Skip to content

Commit

Permalink
skip dki/wdki/smi if less than 2 shells and less than 21 unique dir. …
Browse files Browse the repository at this point in the history 
…don't require rpeb0 bval and bvec
  • Loading branch information
Jenny Chen authored and Jenny Chen committed Oct 22, 2024
1 parent 605c117 commit 5e1ff62
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 18 deletions.
47 changes: 37 additions & 10 deletions designer2/tmi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,13 +265,22 @@ def execute(): #pylint: disable=unused-variable
bvec_dki = bvec[:,bval < maxb]
bval_dki = bval[bval < maxb]

if np.max(bval_dki) < 2:
logger.warning("Max b-value is <2. Skipping DKI and WDKI fit.")
bvec_dki_temp=bvec_dki[:,bval_dki > 0.1]
ndir=np.shape(np.unique(bvec_dki_temp,axis=1))[1]

if len(set(np.round(bval_dki, 2))) <= 2:
logger.warning(f"Fewer than 2 nonzero shells found for DKI fitting. Skipping DKI/WDKI.", extra={"bval_dki": bval_dki.tolist()})
app.ARGS.DKI=False
app.ARGS.WDKI=False
app.ARGS.akc_outliers=False
elif ndir < 21:
logger.warning("Fewer than 21 unique directions detected for DKI fitting. Skipping DKI/WDKI.")
app.ARGS.DKI=False
app.ARGS.WDKI=False
app.ARGS.akc_outliers=False
else:
if len(set(np.round(bval_dki, 2))) <= 2:
logger.warning("Fewer than 2 nonzero shells found for DKI fitting.", extra={"bval_dki": bval_dki.tolist()})
if np.max(bval_dki) < 2:
logger.warning("Max b-value is <2 for DKI fitting.")

grad_dki = np.hstack((bvec_dki.T, bval_dki[None,...].T))
dki = tensor.TensorFitting(grad_dki, int(app.ARGS.n_cores))
Expand Down Expand Up @@ -425,13 +434,22 @@ def execute(): #pylint: disable=unused-variable
bvec_dki = bvec[:,bval < maxb]
bval_dki = bval[bval < maxb]

if np.max(bval_dki) < 2:
logger.warning("Max b-value is <2. Skipping DKI and WDKI fit.")
bvec_dki_temp=bvec_dki[:,bval_dki > 0.1]
ndir=np.shape(np.unique(bvec_dki_temp,axis=1))[1]

if len(set(np.round(bval_dki, 2))) <= 2:
logger.warning(f"Fewer than 2 nonzero shells found for DKI fitting at TE={te}. Skipping DKI/WDKI.", extra={"bval_dki": bval_dki.tolist()})
app.ARGS.DKI=False
app.ARGS.WDKI=False
app.ARGS.akc_outliers=False
elif ndir < 21:
logger.warning("Fewer than 21 unique directions detected for DKI fitting. Skipping DKI/WDKI.")
app.ARGS.DKI=False
app.ARGS.WDKI=False
app.ARGS.akc_outliers=False
else:
if len(set(np.round(bval_dki, 2))) <= 2:
logger.warning(f"Fewer than 2 nonzero shells found for DKI fitting at TE={te}.", extra={"bval_dki": bval_dki.tolist()})
if np.max(bval_dki) < 2:
logger.warning("Max b-value is <2 for DKI fitting.")

grad = np.hstack((bvec_dki.T, bval_dki[None,...].T))
dki = tensor.TensorFitting(grad, int(app.ARGS.n_cores))
Expand Down Expand Up @@ -560,10 +578,19 @@ def execute(): #pylint: disable=unused-variable
warnings.simplefilter('always', UserWarning)

logger.info("Starting SMI fitting process...")
if np.max(bval) < 2:
logger.warning("Max b-value is <2. Skipping SMI fit.")
bvec_temp=bvec[:,bval > 0.1]
ndir=np.shape(np.unique(bvec_temp,axis=1))[1]

if len(set(np.round(bval, 2))) <= 2:
logger.warning(f"Fewer than 2 nonzero shells found for SMI fitting. Skipping SMI.")
app.ARGS.SMI=False
elif ndir < 21:
logger.warning("Fewer than 21 unique directions detected for SMI fitting. Skipping SMI.")
app.ARGS.SMI=False
else:
if np.max(bval) < 2:
logger.warning("Max b-value is <2 for SMI fitting.")

if not app.ARGS.sigma:
logger.warning("No sigma map provided. SMI may be poorly conditioned.")
sigma = None
Expand Down
16 changes: 8 additions & 8 deletions lib/designer_func_wrappers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -490,11 +490,11 @@ def run_eddy(shell_table, dwi_metadata):

rpe_size = [ int(s) for s in image.Header(app.ARGS.rpe_pair).size() ]
if len(rpe_size) == 4:
run.command('mrconvert -coord 3 0 -strides -1,+2,+3 -fslgrad "%s" "%s" -json_import "%s" "%s" "%s"' %
(rpe_bvec_path, rpe_bvals_path, rpe_bids_path, app.ARGS.rpe_pair, 'rpe_b0.mif'))
run.command('mrconvert -coord 3 0 -strides -1,+2,+3 -json_import "%s" "%s" "%s"' %
(rpe_bids_path, app.ARGS.rpe_pair, 'rpe_b0.mif'))
else:
run.command('mrconvert -strides -1,+2,+3 -fslgrad "%s" "%s" -json_import "%s" "%s" "%s"' %
(rpe_bvec_path, rpe_bvals_path, rpe_bids_path, app.ARGS.rpe_pair, 'rpe_b0.mif'))
run.command('mrconvert -strides -1,+2,+3 -json_import "%s" "%s" "%s"' %
(rpe_bids_path, app.ARGS.rpe_pair, 'rpe_b0.mif'))
run.command('mrconvert pe_original_meanb0.mif pe_original_meanb0.nii')
run.command('mrconvert rpe_b0.mif rpe_b0.nii')
else:
Expand Down Expand Up @@ -692,11 +692,11 @@ def run_eddy(shell_table, dwi_metadata):
run.command('dwiextract -bzero dwi.mif - | mrconvert -coord 3 0 - b0pe.mif')
rpe_size = [ int(s) for s in image.Header(app.ARGS.rpe_pair).size() ]
if len(rpe_size) == 4:
run.command('mrconvert "%s" -coord 3 0 -strides -1,+2,+3 -fslgrad "%s" "%s" -json_import "%s" b0rpe.mif' %
(app.ARGS.rpe_pair,rpe_bvec_path, rpe_bvals_path, rpe_bids_path))
run.command('mrconvert "%s" -coord 3 0 -strides -1,+2,+3 -json_import "%s" b0rpe.mif' %
(app.ARGS.rpe_pair, rpe_bids_path))
else:
run.command('mrconvert "%s" -fslgrad "%s" "%s" -strides -1,+2,+3 -json_import "%s" b0rpe.mif' %
(app.ARGS.rpe_pair, rpe_bvec_path, rpe_bvals_path, rpe_bids_path))
run.command('mrconvert "%s" -strides -1,+2,+3 -json_import "%s" b0rpe.mif' %
(app.ARGS.rpe_pair, rpe_bids_path))

run.command('mrinfo b0pe.mif -export_pe_eddy topup_config_1.txt topup_indicies_1.txt')
run.command('mrinfo b0rpe.mif -export_pe_eddy topup_config_2.txt topup_indicies_2.txt')
Expand Down

0 comments on commit 5e1ff62

Please sign in to comment.