diff --git a/.coveragerc b/.coveragerc index 21b501e1c..7b98ae4cf 100644 --- a/.coveragerc +++ b/.coveragerc @@ -2,7 +2,7 @@ branch = True concurrency = multiprocessing data_file = .coverage -source_pkgs = +source_pkgs = scilpy scripts relative_files = True @@ -17,6 +17,9 @@ omit = [report] skip_empty = True skip_covered = True +exclude_also = + if __name__ == "__main__": + (? 1e-5])) + assert np.all(np.asarray(smallests) > min_expected_angle) + + +def test_get_new_gtab_order(): + # Using N=4 vectors + philips_table = np.asarray([[1, 1, 1, 1], + [2, 2, 2, 2], + [3, 3, 3, 3], + [4, 4, 4, 4]]) + dwi = np.ones((10, 10, 10, 4)) + bvecs = np.asarray([[3, 3, 3], + [4, 4, 4], + [2, 2, 2], + [1, 1, 1]]) + bvals = np.asarray([3, 4, 2, 1]) + order = get_new_gtab_order(philips_table, dwi, bvals, bvecs) -def test_get_new_order_philips(): - # Needs dwi files - philips version - pass + assert np.array_equal(bvecs[order, :], philips_table[:, 0:3]) + assert np.array_equal(bvals[order], philips_table[:, 3]) diff --git a/scilpy/gradients/utils.py b/scilpy/gradients/utils.py index 66923434c..f03f051c2 100644 --- a/scilpy/gradients/utils.py +++ b/scilpy/gradients/utils.py @@ -6,7 +6,11 @@ def random_uniform_on_sphere(nb_vectors): """ Creates a set of K pseudo-random unit vectors, following a uniform - distribution on the sphere. + distribution on the sphere. Reference: Emmanuel Caruyer's code + (https://github.com/ecaruyer). + + This is not intended to create a perfect result. It's usually the + initialization step of a repulsion strategy. Parameters ---------- @@ -15,14 +19,21 @@ def random_uniform_on_sphere(nb_vectors): Returns ------- - bvecs: nd.array - pseudo-random unit vector + bvecs: nd.array of shape (nb_vectors, 3) + Pseudo-random unit vectors """ + # Note. Caruyer's docstring says it's a uniform on the half-sphere, but + # plotted a few results: it is one the whole sphere. phi = 2 * np.pi * np.random.rand(nb_vectors) r = 2 * np.sqrt(np.random.rand(nb_vectors)) theta = 2 * np.arcsin(r / 2) + # See here: + # https://www.bogotobogo.com/Algorithms/uniform_distribution_sphere.php + # They seem to be using something like this instead: + # theta = np.arccos(2 * np.random.rand(nb_vectors) - 1.0) + bvecs = np.zeros((nb_vectors, 3)) bvecs[:, 0] = np.sin(theta) * np.cos(phi) bvecs[:, 1] = np.sin(theta) * np.sin(phi) @@ -31,54 +42,59 @@ def random_uniform_on_sphere(nb_vectors): return bvecs -def get_new_order_philips(philips_table, dwi, bvals, bvecs): +def get_new_gtab_order(ref_gradient_table, dwi, bvals, bvecs): """ - Reorder bval and bvec files based on the philips table. + Find the sorting order that could be applied to the bval and bvec files to + obtain the same order as in the reference gradient table. + + This is mostly useful to reorder bval and bvec files in the order they were + acquired by the Philips scanner (before version 5.6). Parameters ---------- - philips_table: nd.array - Philips gradient table - dwis: nibabel - dwis + ref_gradient_table: nd.array + Gradient table, of shape (N, 4). It will use as reference for the + ordering of b-vectors. + Ex: Could be the result of scil_gradients_generate_sampling.py + dwi: nibabel image + dwi of shape (x, y, z, N). Only used to confirm the dwi's shape. bvals : array, (N,) - bvals + bvals that need to be reordered. bvecs : array, (N, 3) - bvecs + bvecs that need to be reordered. Returns ------- new_index: nd.array New index to reorder bvals/bvec """ - # Check number of gradients, bvecs, bvals, dwi and oTable - if len(bvecs) != dwi.shape[3] or len(bvals) != len(philips_table): - raise ValueError('bvec/bval/dwi and original table \ - does not contain the same number of gradients') - - # Check bvals - philips_bval = np.unique(philips_table[:, 3]) + if not (len(bvecs) == dwi.shape[3] == len(bvals) == + len(ref_gradient_table)): + raise ValueError('bvec/bval/dwi and reference table do not contain ' + 'the same number of gradients.') - philips_dwi_shells = philips_bval[philips_bval > 1] - philips_b0s = philips_bval[philips_bval < 1] + ref_bval = np.unique(ref_gradient_table[:, 3]) + ref_dwi_shells = ref_bval[ref_bval > 1] + ref_b0s = ref_bval[ref_bval < 1] dwi_shells = np.unique(bvals[bvals > 1]) b0s = np.unique(bvals[bvals < 1]) - if len(philips_dwi_shells) != len(dwi_shells) or\ - len(philips_b0s) != len(b0s): - raise ValueError('bvec/bval/dwi and original table\ - does not contain the same shells') + if len(ref_dwi_shells) != len(dwi_shells) or \ + len(ref_b0s) != len(b0s): + raise ValueError('bvec/bval/dwi and reference table do not contain ' + 'the same shells.') new_index = np.zeros(bvals.shape) - for nbval in philips_bval: + for nbval in ref_bval: curr_bval = np.where(bvals == nbval)[0] - curr_bval_table = np.where(philips_table[:, 3] == nbval)[0] + curr_bval_table = np.where(ref_gradient_table[:, 3] == nbval)[0] if len(curr_bval) != len(curr_bval_table): - raise ValueError('bval/bvec and orginal table does not contain \ - the same number of gradients for shell {0}'.format(curr_bval)) + raise ValueError('bval/bvec and orginal table do not contain ' + 'the same number of gradients for shell {}.' + .format(curr_bval)) new_index[curr_bval_table] = curr_bval diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index 2de36e347..727e78c04 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -242,7 +242,7 @@ def register_image(static, static_grid2world, moving, moving_grid2world, def compute_snr(dwi, bval, bvec, b0_thr, mask, noise_mask=None, noise_map=None, split_shells=False, - basename=None, verbose=False): + basename=None): """ Compute snr @@ -265,16 +265,10 @@ def compute_snr(dwi, bval, bvec, b0_thr, mask, basename: string Basename used for naming all output files. - verbose: boolean - Set to use logging - Return ------ Dictionary of values (bvec, bval, mean, std, snr) for all volumes. """ - if verbose: - logging.getLogger().setLevel(logging.INFO) - data = dwi.get_fdata(dtype=np.float32) affine = dwi.affine mask = get_data_as_mask(mask, dtype=bool) @@ -417,17 +411,17 @@ def resample_volume(img, ref=None, res=None, iso_min=False, zoom=None, if interp not in interp_choices: raise ValueError("interp must be one of 'nn', 'lin', 'quad', 'cubic'.") - logging.debug('Data shape: %s', data.shape) - logging.debug('Data affine: %s', affine) - logging.debug('Data affine setup: %s', nib.aff2axcodes(affine)) - logging.debug('Resampling data to %s with mode %s', new_zooms, interp) + logging.info('Data shape: %s', data.shape) + logging.info('Data affine: %s', affine) + logging.info('Data affine setup: %s', nib.aff2axcodes(affine)) + logging.info('Resampling data to %s with mode %s', new_zooms, interp) data2, affine2 = reslice(data, affine, original_zooms, new_zooms, _interp_code_to_order(interp)) - logging.debug('Resampled data shape: %s', data2.shape) - logging.debug('Resampled data affine: %s', affine2) - logging.debug('Resampled data affine setup: %s', nib.aff2axcodes(affine2)) + logging.info('Resampled data shape: %s', data2.shape) + logging.info('Resampled data affine: %s', affine2) + logging.info('Resampled data affine setup: %s', nib.aff2axcodes(affine2)) if enforce_dimensions: if ref is None: diff --git a/scilpy/io/utils.py b/scilpy/io/utils.py index cd63a32fd..ca0a4485a 100644 --- a/scilpy/io/utils.py +++ b/scilpy/io/utils.py @@ -225,13 +225,17 @@ def add_force_b0_arg(parser): def add_verbose_arg(parser): - parser.add_argument('-v', action='store_true', dest='verbose', - help='If set, produces verbose output.') + parser.add_argument('-v', default="WARNING", const='INFO', nargs='?', + choices=['DEBUG', 'INFO', 'WARNING'], dest='verbose', + help='Produces verbose output depending on ' + 'the provided level. \nDefault level is warning, ' + 'default when using -v is info.') version = importlib.metadata.version('scilpy') logging.getLogger().setLevel(logging.INFO) logging.info("Scilpy version: {}".format(version)) + logging.getLogger().setLevel(logging.WARNING) def add_bbox_arg(parser): diff --git a/scilpy/reconst/fodf.py b/scilpy/reconst/fodf.py index 3cb1b614e..a64206a0c 100644 --- a/scilpy/reconst/fodf.py +++ b/scilpy/reconst/fodf.py @@ -91,13 +91,13 @@ def get_ventricles_max_fodf(data, fa, md, zoom, args): count += 1 mask[i, j, k] = 1 - logging.debug('Number of voxels detected: {}'.format(count)) + logging.info('Number of voxels detected: {}'.format(count)) if count == 0: logging.warning('No voxels found for evaluation! Change your fa ' 'and/or md thresholds') return 0, mask - logging.debug('Average max fodf value: {}'.format(sum_of_max / count)) + logging.info('Average max fodf value: {}'.format(sum_of_max / count)) return sum_of_max / count, mask diff --git a/scilpy/reconst/frf.py b/scilpy/reconst/frf.py index d40524926..f9c2a3137 100644 --- a/scilpy/reconst/frf.py +++ b/scilpy/reconst/frf.py @@ -107,7 +107,7 @@ def compute_ssst_frf(data, bvals, bvecs, mask=None, mask_wm=None, nvox = np.sum(mask) response, ratio = response_from_mask_ssst(gtab, data, mask) - logging.debug( + logging.info( "Number of indices is {:d} with threshold of {:.2f}".format( nvox, fa_thresh)) fa_thresh -= 0.05 @@ -117,14 +117,14 @@ def compute_ssst_frf(data, bvals, bvecs, mask=None, mask_wm=None, "Could not find at least {:d} voxels with sufficient FA " "to estimate the FRF!".format(min_nvox)) - logging.debug( + logging.info( "Found {:d} voxels with FA threshold {:.2f} for " "FRF estimation".format(nvox, fa_thresh + 0.05)) - logging.debug("FRF eigenvalues: {}".format(str(response[0]))) - logging.debug("Ratio for smallest to largest eigen value " - "is {:.3f}".format(ratio)) - logging.debug("Mean of the b=0 signal for voxels used " - "for FRF: {}".format(response[1])) + logging.info("FRF eigenvalues: {}".format(str(response[0]))) + logging.info("Ratio for smallest to largest eigen value " + "is {:.3f}".format(ratio)) + logging.info("Mean of the b=0 signal for voxels used " + "for FRF: {}".format(response[1])) full_response = np.array([response[0][0], response[0][1], response[0][2], response[1]]) diff --git a/scilpy/reconst/mti.py b/scilpy/reconst/mti.py index bf08679ed..e00492406 100644 --- a/scilpy/reconst/mti.py +++ b/scilpy/reconst/mti.py @@ -4,6 +4,7 @@ import nibabel as nib import numpy as np +import scipy.io import scipy.ndimage from scilpy.io.image import get_data_as_mask @@ -18,9 +19,10 @@ def py_fspecial_gauss(shape, sigma): Parameters ---------- - shape Vector to specify the size of square matrix (rows, columns). - Ex.: (3, 3) - sigma Value for standard deviation (Sigma value of 0.5 is recommended). + shape: Vector to specify the size of square matrix (rows, columns). + Ex.: (3, 3) + sigma: Value for standard deviation (Sigma value of 0.5 is + recommended). Returns ------- @@ -36,16 +38,16 @@ def py_fspecial_gauss(shape, sigma): return h -def compute_contrasts_maps(merged_images, single_echo=False, +def process_contrast_map(merged_images, single_echo=False, filtering=False): """ - Load echoes and compute corresponding contrast map. + Average echoes of a contrast map and apply gaussian filter. Parameters ---------- - 4d_image 4D images - single_echo Apply when only one echo is used to compute contrast maps - filtering Apply Gaussian filtering to remove Gibbs ringing + 4d_image: 4D image, containing every echoes of the contrast map + single_echo: Apply when only one echo is used to compute contrast map + filtering: Apply Gaussian filtering to remove Gibbs ringing (default is False). Returns @@ -68,13 +70,15 @@ def compute_contrasts_maps(merged_images, single_echo=False, return contrast_map -def compute_saturation(cPD1, cPD2, cT1, acq_parameters): +def compute_saturation_map(MT, PD, T1, a, TR): """ - Compute saturation of given contrast. - Saturation is computed from apparent longitudinal relaxation rate - (R1app) and apparent signal amplitude (Aapp). The estimation of the + Compute saturation of given contrast (MT). + Saturation is computed from apparent longitudinal relaxation time + (T1app) and apparent signal amplitude (Aapp), as the difference of + longitudinal relaxation times of bound to macromolecules pool + from free water pool saturation. The estimation of the saturation includes correction for the effects of excitation flip angle - and longitudinal relaxation rate, and remove the effect of T1-weighted + and longitudinal relaxation time, and remove the effect of T1-weighted image. see Helms et al., 2008 @@ -82,259 +86,293 @@ def compute_saturation(cPD1, cPD2, cT1, acq_parameters): Parameters ---------- - cPD1: Reference proton density (MT-off) - cPD2: Contrast of choice (can be a single contrast or the - mean of positive and negative for example) - cT1: Contrast T1-weighted - acq_parameters: List of TR and Flipangle for MT contrast and T1w images - [TR, Flipangle] + MT: Contrast of choice (can be a single contrast or the + mean of positive and negative for example) + PD: Reference proton density weighted image (MToff PD) + T1: Reference T1 weighted image (MToff T1) + a: List of two flip angles corresponding to PD and T1. If B1 + correction is on, this should be two 3D-array of flip + angles varying spatially with respect to the B1 map. + TR: List of two repetition times corresponding to PD and T1. + Returns ------- - MT ratio and MT saturation matrice in 3D-array. + MT saturation matrice in 3D-array (sat), computed apparent T1 map (T1app). """ - Aapp_num = ((2*acq_parameters[0][0] / (acq_parameters[0][1]**2)) - - (2*acq_parameters[1][0] / (acq_parameters[1][1]**2))) - Aapp_den = (((2*acq_parameters[0][0]) / (acq_parameters[0][1]*cPD1)) - - ((2*acq_parameters[1][0]) / (acq_parameters[1][1]*cT1))) + Aapp_num = (TR[0] / (a[0] ** 2)) - (TR[1] / (a[1] ** 2)) + Aapp_den = (TR[0] / (a[0] * PD)) - (TR[1] / (a[1] * T1)) Aapp = Aapp_num / Aapp_den - R1app_num = ((cPD1 / acq_parameters[0][1]) - (cT1 / acq_parameters[1][1])) - R1app_den = ((cT1*acq_parameters[1][1]) / (2*acq_parameters[1][0]) - - (cPD1*acq_parameters[0][1]) / (2*acq_parameters[0][0])) - R1app = R1app_num / R1app_den + T1app_num = (PD / a[0]) - (T1 / a[1]) + T1app_den = 0.5 * ((T1 * a[1]) / TR[1] - (PD * a[0]) / TR[0]) + T1app = T1app_num / T1app_den - sat = 100*(((Aapp*acq_parameters[0][1]*acq_parameters[0][0] / R1app - ) / cPD2) - (acq_parameters[0][0] / R1app) - - (acq_parameters[0][1]**2) / 2) + sat = Aapp * a[0] * TR[0] / T1app / MT - (TR[0] / T1app) - (a[0] ** 2) / 2 + sat *= 100 - return sat + return sat, T1app -def compute_ihMT_maps(contrasts_maps, acq_parameters): +def compute_ratio_map(mt_on_single, mt_off, mt_on_dual=None): """ - Compute Inhomogenous Magnetization transfer ratio and saturation maps. + Compute magnetization transfer ratio (MTR), and inhomogenous magnetization + transfer ratio (ihMTR) if mt_on_dual is given. + - MT ratio (MTR) is computed as the percentage difference of single + frequency mean image from the reference. - ihMT ratio (ihMTR) is computed as the percentage difference of dual from single frequency. - - ihMT saturation (ihMTsat) is computed as the difference of longitudinal - relaxation rates of bound to macromolecules pool from free water pool - saturation. The estimation of the ihMT saturation includes correction for - the effects of excitation flip angle and longitudinal relaxation rate, and - remove the effect of T1-weighted image. - cPD : contrast proton density - a : mean of positive and negative proton density - b : mean of two dual proton density - cT1 : contrast T1-weighted - num : numberator - den : denumerator - - see Varma et al., 2015 + + For ihMTR, see Varma et al., 2015 www.sciencedirect.com/science/article/pii/S1090780715001998 - see Manning et al., 2017 - www.sciencedirect.com/science/article/pii/S1090780716302488?via%3Dihub Parameters ---------- - contrasts_maps: List of matrices : list of all contrast maps computed - with compute_contrasts_maps - acq_parameters: List of TR and Flipangle values for ihMT and T1w images - [TR, Flipangle] + mt_on_single: MT-on image with single frequency pulse: can be one + single positive/negative frequency MT image or the + average of both images (positive and negative). + mt_off: MT-off image: the reference image without MT + preparation. + mt_on_dual: MT-on image with dual frequency pulse: can be one + dual alternating positive/negative frequency MT image + or the average of both images. Optional. If given, will + compute the ihMTR also. + Returns ------- - ihMT ratio (ihMTR) and ihMT saturation (ihMTsat) matrices in 3D-array. - + MT ratio (MTR), ihMT ratio (ihMTR). """ - # Compute ihMT ratio map - ihMTR = 100*(contrasts_maps[4] + contrasts_maps[3] - - contrasts_maps[1] - contrasts_maps[0]) / contrasts_maps[2] - - # Compute MT saturation maps - cPD1 = contrasts_maps[2] - cPDa = (contrasts_maps[4] + contrasts_maps[3]) / 2 - cPDb = (contrasts_maps[1] + contrasts_maps[0]) / 2 - cT1 = contrasts_maps[5] + # Compute MT Ratio map + MTR = 100 * ((mt_off - mt_on_single) / mt_off) - MT_sat_single = compute_saturation(cPD1, cPDa, cT1, acq_parameters) - MT_sat_dual = compute_saturation(cPD1, cPDb, cT1, acq_parameters) - ihMTsat = MT_sat_dual - MT_sat_single + # Compute ihMT ratio map + if mt_on_dual is not None: + # The factor 2 is there to account for the /2 in mt_on mean images. + ihMTR = 2 * 100 * (mt_on_single - mt_on_dual) / mt_off + return MTR, ihMTR - return ihMTR, ihMTsat + return MTR -def compute_MT_maps_from_ihMT(contrasts_maps, acq_parameters): +def threshold_map(computed_map, in_mask, + lower_threshold, upper_threshold, + idx_contrast_list=None, contrast_maps=None): """ - Compute Magnetization transfer ratio and saturation maps. - MT ratio is computed as the percentage difference of two images, one - acquired with off-resonance saturation (MT-on) and one without (MT-off). - MT saturation is computed from apparent longitudinal relaxation rate - (R1app) and apparent signal amplitude (Aapp). The estimation of the MT - saturation includes correction for the effects of excitation flip angle - and longitudinal relaxation rate, and remove the effect of T1-weighted - image. - cPD : contrast proton density - 1 : reference proton density (MT-off) - 2 : mean of positive and negative proton density (MT-on) - cT1 : contrast T1-weighted - num : numberator - den : denumerator - - see Helms et al., 2008 - https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.21732 + Remove NaN and apply different threshold based on + - maximum and minimum threshold value + - T1 mask + - combination of specific contrast maps + idx_contrast_list and contrast_maps are required for + thresholding of ihMT images. Parameters ---------- - contrasts_maps: List of 3D-array constrats matrices : list of all - contrast maps computed with compute_ihMT_contrasts - acq_parameters: List of TR and Flipangle for ihMT and T1w images - [TR, Flipangle] + computed_map: 3D-Array data. + Myelin map (ihMT or non-ihMT maps) + in_mask: Path to binary T1 mask from T1 segmentation. + Must be the sum of GM+WM+CSF. + lower_threshold: Value for low thresold + upper_thresold: Value for up thresold + idx_contrast_list: List of indexes of contrast maps corresponding to + that of input contrast_maps ex.: [0, 2, 4] + Altnp = 0; Atlpn = 1; Negative = 2; Positive = 3; + PD = 4; T1 = 5 + contrast_maps: List of 3D-Array. File must containing the + 5 or 6 contrast maps. Returns - ------- - MT ratio and MT saturation matrice in 3D-array. + ---------- + Thresholded matrix in 3D-array. """ - # Compute MT Ratio map - MTR = 100*((contrasts_maps[2] - - (contrasts_maps[4] + contrasts_maps[3]) / 2) / - contrasts_maps[2]) + # Remove NaN and apply thresold based on lower and upper value + computed_map[np.isnan(computed_map)] = 0 + computed_map[np.isinf(computed_map)] = 0 + computed_map[computed_map < lower_threshold] = 0 + computed_map[computed_map > upper_threshold] = 0 - # Compute MT saturation maps - cPD1 = contrasts_maps[2] - cPD2 = (contrasts_maps[4] + contrasts_maps[3]) / 2 - cT1 = contrasts_maps[5] + # Load and apply sum of T1 probability maps on myelin maps + if in_mask is not None: + mask_image = nib.load(in_mask) + mask_data = get_data_as_mask(mask_image) + computed_map[np.where(mask_data == 0)] = 0 - MTsat = compute_saturation(cPD1, cPD2, cT1, acq_parameters) + # Apply threshold based on combination of specific contrast maps + if idx_contrast_list and contrast_maps: + for idx in idx_contrast_list: + computed_map[contrast_maps[idx] == 0] = 0 - return MTR, MTsat + return computed_map -def compute_MT_maps(contrasts_maps, acq_parameters): +def adjust_B1_map_intensities(B1_map, nominal=100): """ - Compute Magnetization transfer ratio and saturation maps. - MT ratio is computed as the percentage difference of two images, one - acquired with off-resonance saturation (MT-on) and one without (MT-off). - MT saturation is computed from apparent longitudinal relaxation rate - (R1app) and apparent signal amplitude (Aapp). The estimation of the MT - saturation includes correction for the effects of excitation flip angle - and longitudinal relaxation rate, and remove the effect of T1-weighted - image. - cPD : contrast proton density - 1 : reference proton density (MT-off) - 2 : mean of positive and negative proton density (MT-on) - cT1 : contrast T1-weighted - num : numberator - den : denumerator + Adjust and verify the values of the B1 map. We want the B1 map to have + values around 1.0, so we divide by the nominal B1 value if it is not + already the case. Sometimes the scaling of the B1 map is wrong, leading + to weird values after this process, which what is verified here. - see Helms et al., 2008 - https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.21732 + Parameters + ---------- + B1_map: B1 coregister map. + nominal: Nominal values of B1. For Philips, should be 1. + + Returns + ---------- + Ajusted B1 map in 3d-array. + """ + B1_map /= nominal + med_bB = np.nanmedian(np.where(B1_map == 0, np.nan, B1_map)) + if not np.isclose(med_bB, 1.0, atol=0.2): + raise ValueError("Intensities of the B1 map are wrong.") + B1_map = np.clip(B1_map, 0.5, 1.5) + return B1_map + + +def smooth_B1_map(B1_map, wdims=5): + """ + Apply a smoothing to the B1 map. Parameters ---------- - contrasts_maps: List of 3D-array constrats matrices : list of all - contrast maps computed with compute_ihMT_contrasts - acq_parameters: List of TR and Flipangle for ihMT and T1w images - [TR, Flipangle] + B1_map: B1 coregister map. + wdims: Window dimension (in voxels) for the smoothing. + Returns - ------- - MT ratio and MT saturation matrice in 3D-array. + ---------- + Smoothed B1 map in 3d-array. """ - # Compute MT Ratio map - MTR = 100*(contrasts_maps[0] - contrasts_maps[1]) / contrasts_maps[0] + h = np.ones((wdims, wdims, 1)) / wdims ** 2 + B1_map_smooth = scipy.ndimage.convolve(B1_map, h).astype(np.float32) + return B1_map_smooth - # Compute MT saturation maps: cPD1 = mt-off; cPD2 = mt-on - cPD1 = contrasts_maps[0] - cPD2 = contrasts_maps[1] - cT1 = contrasts_maps[2] - Aapp_num = ((2*acq_parameters[0][0] / (acq_parameters[0][1]**2)) - - (2*acq_parameters[1][0] / (acq_parameters[1][1]**2))) - Aapp_den = (((2*acq_parameters[0][0]) / (acq_parameters[0][1]*cPD1)) - - ((2*acq_parameters[1][0]) / (acq_parameters[1][1]*cT1))) - Aapp = Aapp_num / Aapp_den +def apply_B1_corr_empiric(MT_map, B1_map): + """ + Apply an empiric B1 correction on MT map. - R1app_num = ((cPD1 / acq_parameters[0][1]) - (cT1 / acq_parameters[1][1])) - R1app_den = ((cT1*acq_parameters[1][1]) / (2*acq_parameters[1][0]) - - (cPD1*acq_parameters[0][1]) / (2*acq_parameters[0][0])) - R1app = R1app_num / R1app_den + see Weiskopf et al., 2013 + https://www.frontiersin.org/articles/10.3389/fnins.2013.00095/full - MTsat = 100*(((Aapp*acq_parameters[0][1]*acq_parameters[0][0]/R1app) / - cPD2) - (acq_parameters[0][0]/R1app) - - (acq_parameters[0][1]**2)/2) + Parameters + ---------- + MT_map: 3D-Array MT map. + B1_map: B1 coregister map. - return MTR, MTsat + Returns + ---------- + Corrected MT matrix in 3D-array. + """ + MT_map_B1_corrected = MT_map*(1.0-0.4)/(1-0.4*(B1_map)) + return MT_map_B1_corrected -def threshold_maps(computed_map, in_mask, - lower_threshold, upper_threshold, - idx_contrast_list=None, contrasts_maps=None): +def apply_B1_corr_model_based(MTsat, B1_map, R1app, fitvalues_file, B1_ref=1): """ - Remove NaN and apply different threshold based on - - maximum and minimum threshold value - - T1 mask - - combination of specific contrasts maps - idx_contrast_list and contrasts_maps are required for - thresholding of ihMT images. + Apply a model-based B1 correction on MT map. + + see Rowley et al., 2021 + https://onlinelibrary.wiley.com/doi/10.1002/mrm.28831 Parameters ---------- - computed_map 3D-Array data. - Myelin map (ihMT or non-ihMT maps) - in_mask Path to binary T1 mask from T1 segmentation. - Must be the sum of GM+WM+CSF. - lower_threshold Value for low thresold - upper_thresold Value for up thresold - idx_contrast_list List of indexes of contrast maps corresponding to - that of input contrasts_maps ex.: [0, 2, 5] - Altnp = 0; Atlpn = 1; Reference = 2; Negative = 3; - Positive = 4; T1weighted = 5 - contrasts_maps List of 3D-Array. File must containing the - 6 contrasts maps. + MTsat: 3D-Array MTsat map. + B1_map: B1 coregister map. + R1app: Apparent R1 map, obtained from compute_saturation_map. + fitvalues_file: Path to the fitValues_*.mat file corresponding to the + MTsat map. This file is computed with Christopher + Rowley's Matlab library. + B1_ref: Reference B1 value used for the fit (usually 1). + Returns ---------- - Thresholded matrix in 3D-array. + Corrected MTsat matrix in 3D-array. """ - # Remove NaN and apply thresold based on lower and upper value - computed_map[np.isnan(computed_map)] = 0 - computed_map[np.isinf(computed_map)] = 0 - computed_map[computed_map < lower_threshold] = 0 - computed_map[computed_map > upper_threshold] = 0 + cf_eq, R1_to_M0b = _read_fitvalues_from_file(fitvalues_file) + cf_map = _compute_B1_corr_factor_map(B1_map, R1app, cf_eq, + R1_to_M0b, B1_ref=B1_ref) + MTsat_corr = MTsat + MTsat * cf_map + return MTsat_corr - # Load and apply sum of T1 probability maps on myelin maps - mask_image = nib.load(in_mask) - mask_data = get_data_as_mask(mask_image) - computed_map[np.where(mask_data == 0)] = 0 - # Apply threshold based on combination of specific contrasts maps - if idx_contrast_list and contrasts_maps: - for idx in idx_contrast_list: - computed_map[contrasts_maps[idx] == 0] = 0 +def _read_fitvalues_from_file(fitvalues_file): + """ + Extract equations from fitValues_*.mat file. - return computed_map + see Rowley et al., 2021 + https://onlinelibrary.wiley.com/doi/10.1002/mrm.28831 + Parameters + ---------- + fitvalues_file: Path to the fitValues_*.mat file corresponding to the + MTsat map. This file is computed with Christopher + Rowley's Matlab library. -def apply_B1_correction(MT_map, B1_map): + Returns + ---------- + Correction factor equation (cf_eq), R1 to M0b equation (R1_to_M0b) + """ + fitvalues = scipy.io.loadmat(fitvalues_file) + fit_SS_eq = fitvalues['fitValues']['fit_SS_eqn'][0][0][0] + est_M0b_from_R1 = fitvalues['fitValues']['Est_M0b_from_R1'][0][0][0] + cf_eq = fit_SS_eq.replace('^', + '**').replace('Raobs.', + 'Raobs').replace('M0b.', + 'M0b').replace('b1.', + 'B1') + R1_to_M0b = est_M0b_from_R1.replace('^', + '**').replace('Raobs.', + 'Raobs').replace('M0b.', + 'M0b') + return cf_eq, R1_to_M0b + + +def _compute_B1_corr_factor_map(B1_map, Raobs, cf_eq, R1_to_M0b, B1_ref=1): """ - Function to apply an empiric B1 correction. + Compute the correction factor map for B1 correction. - see Weiskopf et al., 2013 - https://www.frontiersin.org/articles/10.3389/fnins.2013.00095/full + see Rowley et al., 2021 + https://onlinelibrary.wiley.com/doi/10.1002/mrm.28831 Parameters ---------- - MT_map 3D-Array Myelin map. - B1_map Path to B1 coregister map. + B1_map: B1 coregister map. + Raobs: R1 map, obtained from compute_saturation_map. + cf_eq: Correction factor equation extracted with + _read_fitvalues_from_file. + R1_to_M0b: Conversion equation from R1 to M0b extracted with + _read_fitvalues_from_file. + B1_ref: Reference B1 value used for the fit (usually 1). Returns ---------- - Corrected MT matrix in 3D-array. + Correction factor map. """ - # Load B1 image - B1_img = nib.load(B1_map) - B1_img_data = B1_img.get_fdata(dtype=np.float32) + M0b = eval(R1_to_M0b, {"Raobs": Raobs}) - # Apply a light smoothing to the B1 map - h = np.ones((5, 5, 1))/25 - B1_smooth_map = scipy.ndimage.convolve(B1_img_data, - h).astype(np.float32) + B1 = B1_map * B1_ref + cf_act = eval(cf_eq, {"Raobs": Raobs, "B1": B1, "M0b": M0b}) + B1 = B1_ref + cf_nom = eval(cf_eq, {"Raobs": Raobs, "B1": B1, "M0b": M0b}) - # Apply an empiric B1 correction via B1 smooth data on MT data - MT_map_B1_corrected = MT_map*(1.0-0.4)/(1-0.4*(B1_smooth_map/100)) + cf_map = (cf_nom - cf_act) / cf_act + return cf_map - return MT_map_B1_corrected + +def adjust_B1_map_header(B1_img, slope): + """ + Fixe issue in B1 map header, by applying the scaling (slope) and setting + the slope to 1. + + Parameters + ---------- + B1_img: B1 nifti image object. + slope: Slope value, obtained from the image header. + + Returns + ---------- + Adjusted B1 nifti image object (new_b1_img) + """ + B1_map = B1_img.get_fdata() + B1_map /= slope + B1_img.header.set_slope_inter(1, 0) + new_b1_img = nib.nifti1.Nifti1Image(B1_map, B1_img.affine, + header=B1_img.header) + return new_b1_img diff --git a/scilpy/reconst/tests/test_mti.py b/scilpy/reconst/tests/test_mti.py index 50c06499d..60165666b 100644 --- a/scilpy/reconst/tests/test_mti.py +++ b/scilpy/reconst/tests/test_mti.py @@ -6,36 +6,46 @@ def test_py_fspecial_gauss(): pass -def test_compute_contrasts_maps(): +def test_process_contrast_map(): # toDO pass -def test_compute_saturation(): +def test_compute_saturation_map(): # toDO pass -def test_compute_ihMT_maps(): +def test_compute_ratio_map(): # toDO pass -def test_compute_MT_maps_from_ihMT(): +def test_threshold_map(): # toDO pass -def test_compute_MT_maps(): +def test_adjust_B1_map_intensities(): # toDO pass -def test_threshold_maps(): +def test_smooth_B1_map(): # toDO pass -def test_apply_B1_correction(): +def test_apply_B1_corr_empiric(): + # toDO + pass + + +def test_apply_B1_corr_model_based(): + # toDO + pass + + +def test_adjust_B1_map_header(): # toDO pass diff --git a/scilpy/stats/matrix_stats.py b/scilpy/stats/matrix_stats.py index 8f1767d1e..d77b5e1fc 100644 --- a/scilpy/stats/matrix_stats.py +++ b/scilpy/stats/matrix_stats.py @@ -67,8 +67,8 @@ def ttest_two_matrices(matrices_g1, matrices_g2, paired, tail, fdr, sum_both_groups = np.sum(matrices_g1, axis=2) + np.sum(matrices_g2, axis=2) nbr_non_zeros = np.count_nonzero(np.triu(sum_both_groups)) - logging.debug('The provided matrices contain {} non zeros elements.' - .format(nbr_non_zeros)) + logging.info('The provided matrices contain {} non zeros elements.' + .format(nbr_non_zeros)) matrices_g1 = matrices_g1.reshape((np.prod(matrix_shape), nb_group_g1)) matrices_g2 = matrices_g2.reshape((np.prod(matrix_shape), nb_group_g2)) @@ -76,10 +76,10 @@ def ttest_two_matrices(matrices_g1, matrices_g2, paired, tail, fdr, matrix_pval = np.ones(np.prod(matrix_shape)) * -0.000001 text = ' paired' if paired else '' - logging.debug('Performing{} t-test with "{}" hypothesis.' - .format(text, tail)) - logging.debug('Data has dimensions {}x{} with {} and {} observations.' - .format(matrix_shape[0], matrix_shape[1], + logging.info('Performing{} t-test with "{}" hypothesis.' + .format(text, tail)) + logging.info('Data has dimensions {}x{} with {} and {} observations.' + .format(matrix_shape[0], matrix_shape[1], nb_group_g1, nb_group_g2)) # For conversion to p-values @@ -105,7 +105,7 @@ def ttest_two_matrices(matrices_g1, matrices_g2, paired, tail, fdr, corr_matrix_pval = matrix_pval.reshape(matrix_shape) if fdr: - logging.debug('Using FDR, the results will be q-values.') + logging.info('Using FDR, the results will be q-values.') corr_matrix_pval = np.triu(corr_matrix_pval) corr_matrix_pval[corr_matrix_pval > 0] = multipletests( corr_matrix_pval[corr_matrix_pval > 0], 0, method='fdr_bh')[1] @@ -158,8 +158,8 @@ def omega_sigma(matrix): transitivity_latt_list = [] path_length_rand_list = [] for i in range(10): - logging.debug('Generating random and lattice matrices, ' - 'iteration #{}.'.format(i)) + logging.info('Generating random and lattice matrices, ' + 'iteration #{}.'.format(i)) random = bct.randmio_und(matrix, 10)[0] lattice = bct.latmio_und(matrix, 10)[1] diff --git a/scilpy/stats/stats.py b/scilpy/stats/stats.py index 50ff83faa..f02246664 100644 --- a/scilpy/stats/stats.py +++ b/scilpy/stats/stats.py @@ -32,10 +32,10 @@ def verify_normality(data, alpha=0.05): # First, we verify if sample pass Shapiro-Wilk test W, p_value = scipy.stats.shapiro(data) if p_value < alpha and len(data) < 30: - logging.debug('The data sample can not be considered normal') + logging.info('The data sample can not be considered normal') normality = False else: - logging.debug('The data sample pass the normality assumption.') + logging.info('The data sample pass the normality assumption.') normality = True return normality, p_value @@ -76,12 +76,12 @@ def verify_homoscedasticity(data_by_group, normality=False, alpha=0.05): else: test = 'Levene' W, p_value = scipy.stats.levene(*data_by_group) - logging.debug('Test name: {}'.format(test)) + logging.info('Test name: {}'.format(test)) if p_value < alpha and mean_nb < 30: - logging.debug('The sample didnt pass the equal variance assumption') + logging.info('The sample didnt pass the equal variance assumption') homoscedasticity = False else: - logging.debug('The sample pass the equal variance assumption') + logging.info('The sample pass the equal variance assumption') homoscedasticity = True return test, homoscedasticity, p_value @@ -145,12 +145,12 @@ def verify_group_difference(data_by_group, normality=False, test = 'Kruskalwallis' T, p_value = scipy.stats.kruskal(*data_by_group) - logging.debug('Test name: {}'.format(test)) + logging.info('Test name: {}'.format(test)) if p_value < alpha: - logging.debug('There is a difference between groups') + logging.info('There is a difference between groups') difference = True else: - logging.debug('We are not able to detect difference between the groups.') + logging.info('We are not able to detect difference between the groups.') difference = False return test, difference, p_value @@ -191,9 +191,9 @@ def verify_post_hoc(data_by_group, groups_list, test, test : string Name of the test done to verify group difference """ - logging.debug('We need to do a post-hoc analysis since ' - 'there is a difference') - logging.debug('Post-hoc: {} pairwise'.format(test)) + logging.info('We need to do a post-hoc analysis since ' + 'there is a difference') + logging.info('Post-hoc: {} pairwise'.format(test)) differences = [] nb_group = len(groups_list) @@ -214,7 +214,7 @@ def verify_post_hoc(data_by_group, groups_list, test, data_by_group[x], data_by_group[y]) differences.append((groups_list[x], groups_list[y], p_value < alpha, p_value)) - logging.debug('Result:') - logging.debug(differences) + logging.info('Result:') + logging.info(differences) return test, differences diff --git a/scilpy/stats/utils.py b/scilpy/stats/utils.py index 6c6c837c5..973128a6a 100644 --- a/scilpy/stats/utils.py +++ b/scilpy/stats/utils.py @@ -51,8 +51,8 @@ def __init__(self, json_file, participants): self.data_dictionnary[participant['participant_id']]\ [variable] = participant[variable] - logging.debug('Data_dictionnary') - logging.debug(self.data_dictionnary[self.get_first_participant()]) + logging.info('Data_dictionnary') + logging.info(self.data_dictionnary[self.get_first_participant()]) with open('data.json', 'w') as fp: json.dump(self.data_dictionnary, fp, indent=4) @@ -64,15 +64,15 @@ def validation_participant_id(self, json_info, participants_info): # Create the list of participants id from the json dictionnary participants_from_json = list(json_info.keys()) - logging.debug('participant list from json dictionnary:') - logging.debug(participants_from_json) + logging.info('participant list from json dictionnary:') + logging.info(participants_from_json) # Create the list of participants id from the tsv list of dictionnary participants_from_tsv = [] for participant in participants_info: participants_from_tsv.append(participant['participant_id']) - logging.debug('participant list from tsv file:') - logging.debug(participants_from_tsv) + logging.info('participant list from tsv file:') + logging.info(participants_from_tsv) # Compare the two list participants_from_json.sort() @@ -80,15 +80,17 @@ def validation_participant_id(self, json_info, participants_info): if not participants_from_json == participants_from_tsv: if not len(participants_from_json) == len(participants_from_tsv): - logging.debug('The number of participants from json file is not the same ' - 'as the one in the tsv file.') + logging.info('The number of participants from json file is ' + 'not the same as the one in the tsv file.') is_in_tsv = np.in1d(participants_from_json, participants_from_tsv) is_in_json = np.in1d(participants_from_tsv, participants_from_json) - logging.debug('participants list from json file missing in tsv file :') - logging.debug(np.asarray(participants_from_json)[~is_in_tsv]) - logging.debug('participants list from tsv file missing in json file :') - logging.debug(np.asarray(participants_from_tsv)[~is_in_json]) + logging.info('participants list from json file missing in tsv ' + 'file :') + logging.info(np.asarray(participants_from_json)[~is_in_tsv]) + logging.info('participants list from tsv file missing in json ' + 'file :') + logging.info(np.asarray(participants_from_tsv)[~is_in_json]) logging.error('The subjects from the json file does not fit ' 'with the subjects of the tsv file. ' @@ -97,7 +99,7 @@ def validation_participant_id(self, json_info, participants_info): 'with the subjects of the tsv file. ' 'Impossible to build the data_for_stat object') else: - logging.debug('The json and the tsv are compatible') + logging.info('The json and the tsv are compatible') def get_participants_list(self): # Construct the list of participant_id from the data_dictionnary @@ -492,6 +494,6 @@ def visualise_distribution(data_by_group, participants_id, bundle, metric, fig.savefig(os.path.join(oFolder, 'Graph', bundle, metric)) - logging.debug('outliers:[(id, group)]') - logging.debug(outliers) + logging.info('outliers:[(id, group)]') + logging.info(outliers) return outliers diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py index d81a7c902..208a33fc6 100644 --- a/scilpy/tracking/tracker.py +++ b/scilpy/tracking/tracker.py @@ -179,8 +179,8 @@ def _set_nbr_processes(self, nbr_processes): if nbr_processes > self.nbr_seeds: nbr_processes = self.nbr_seeds - logging.debug("Setting number of processes to {} since there were " - "less seeds than processes.".format(nbr_processes)) + logging.info("Setting number of processes to {} since there were " + "less seeds than processes.".format(nbr_processes)) return nbr_processes def _prepare_multiprocessing_pool(self, tmpdir): diff --git a/scilpy/tractanalysis/features.py b/scilpy/tractanalysis/features.py index a49724548..267f0611c 100644 --- a/scilpy/tractanalysis/features.py +++ b/scilpy/tractanalysis/features.py @@ -115,9 +115,9 @@ def remove_loops_and_sharp_turns(streamlines, if tm.mean_curvature(clusters.centroids[i]) <= mean_curvature: ids.extend(clusters[i].indices) else: - logging.debug("Impossible to use the use_qb option because " + - "not more than one streamline left from the\n" + - "input file.") + logging.info("Impossible to use the use_qb option because " + + "not more than one streamline left from the\n" + + "input file.") return ids diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 77f8b59ae..fe1c13716 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -176,19 +176,19 @@ def filter_streamlines_by_total_length_per_dim( total_per_orientation = np.abs(np.asarray( [np.sum(d, axis=0) for d in all_dirs])) - logging.debug("Total length per orientation is:\n" - "Average: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" - "Min: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" - "Max: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" - .format(np.mean(total_per_orientation[:, 0]), - np.mean(total_per_orientation[:, 1]), - np.mean(total_per_orientation[:, 2]), - np.min(total_per_orientation[:, 0]), - np.min(total_per_orientation[:, 1]), - np.min(total_per_orientation[:, 2]), - np.max(total_per_orientation[:, 0]), - np.max(total_per_orientation[:, 1]), - np.max(total_per_orientation[:, 2]))) + logging.info("Total length per orientation is:\n" + "Average: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" + "Min: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" + "Max: x: {:.2f}, y: {:.2f}, z: {:.2f} \n" + .format(np.mean(total_per_orientation[:, 0]), + np.mean(total_per_orientation[:, 1]), + np.mean(total_per_orientation[:, 2]), + np.min(total_per_orientation[:, 0]), + np.min(total_per_orientation[:, 1]), + np.min(total_per_orientation[:, 2]), + np.max(total_per_orientation[:, 0]), + np.max(total_per_orientation[:, 1]), + np.max(total_per_orientation[:, 2]))) # Find good ids mask_good_x = np.logical_and(limits_x[0] < total_per_orientation[:, 0], @@ -264,11 +264,11 @@ def resample_streamlines_step_size(sft, step_size): if step_size == 0: raise ValueError("Step size can't be 0!") elif step_size < 0.1: - logging.debug("The value of your step size seems suspiciously low. " - "Please check.") + logging.info("The value of your step size seems suspiciously low. " + "Please check.") elif step_size > np.max(sft.voxel_sizes): - logging.debug("The value of your step size seems suspiciously high. " - "Please check.") + logging.info("The value of your step size seems suspiciously high. " + "Please check.") # Make sure we are in world space orig_space = sft.space @@ -299,9 +299,9 @@ def _warn_and_save(new_streamlines, sft): Warn that we loose data_per_point, then create resampled SFT.""" if sft.data_per_point is not None and sft.data_per_point.keys(): - logging.debug("Initial StatefulTractogram contained data_per_point. " - "This information will not be carried in the final " - "tractogram.") + logging.info("Initial StatefulTractogram contained data_per_point. " + "This information will not be carried in the final " + "tractogram.") new_sft = StatefulTractogram.from_sft( new_streamlines, sft, data_per_streamline=sft.data_per_streamline) @@ -328,7 +328,7 @@ def smooth_line_gaussian(streamline, sigma): raise ValueError('Cant have a 0 sigma with gaussian.') if length(streamline) < 1: - logging.debug('Streamline shorter than 1mm, corner cases possible.') + logging.info('Streamline shorter than 1mm, corner cases possible.') # Smooth each dimension separately x, y, z = streamline.T @@ -367,7 +367,7 @@ def smooth_line_spline(streamline, smoothing_parameter, nb_ctrl_points): raise ValueError('Cant have a 0 sigma with spline.') if length(streamline) < 1: - logging.debug('Streamline shorter than 1mm, corner cases possible.') + logging.info('Streamline shorter than 1mm, corner cases possible.') if nb_ctrl_points < 3: nb_ctrl_points = 3 diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 55b80bd82..64c794d42 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -809,12 +809,12 @@ def split_sft_randomly_per_cluster(orig_sft, chunk_sizes, seed, thresholds): nb_chunks = len(chunk_sizes) percent_kept_per_chunk = [nb / len(orig_sft) for nb in chunk_sizes] - logging.debug("Computing QBx") + logging.info("Computing QBx") clusters = qbx_and_merge(orig_sft.streamlines, thresholds, nb_pts=20, verbose=False) - logging.debug("Done. Now getting list of indices in each of the {} " - "cluster.".format(len(clusters))) + logging.info("Done. Now getting list of indices in each of the {} " + "cluster.".format(len(clusters))) total_indices = [[] for _ in range(nb_chunks + 1)] for cluster in clusters: if len(cluster.indices) > 1: diff --git a/scilpy/utils/metrics_tools.py b/scilpy/utils/metrics_tools.py index 809d2cc18..66860549d 100644 --- a/scilpy/utils/metrics_tools.py +++ b/scilpy/utils/metrics_tools.py @@ -157,9 +157,10 @@ def weighted_mean_std(weights, data): a tuple containing the mean and standard deviation of the data """ - masked_data = np.ma.masked_array(data, np.isnan(data)) - mean = np.average(masked_data, weights=weights) - variance = np.average((masked_data-mean)**2, weights=weights) + masked_data = np.ma.masked_array(data, np.logical_or(np.isnan(data), + np.isinf(data))) + mean = np.ma.average(masked_data, weights=weights) + variance = np.ma.average((masked_data-mean)**2, weights=weights) return mean, np.sqrt(variance) diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index 82b12e40a..b5eb9853d 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -91,6 +91,15 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + # COMMIT has some c-level stdout and non-logging print that cannot + # be easily stopped. Manual redirection of all printed output + if args.verbose == "WARNING": + f = io.StringIO() + redirected_stdout = redirect_stdout(f) + redirect_stdout_c() + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + redirected_stdout = redirect_stdout(sys.stdout) if args.compute_only and not args.save_kernels: parser.error('--compute_only must be used with --save_kernels.') @@ -102,16 +111,6 @@ def main(): args.out_dir, optional=args.save_kernels) - # COMMIT has some c-level stdout and non-logging print that cannot - # be easily stopped. Manual redirection of all printed output - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - redirected_stdout = redirect_stdout(sys.stdout) - else: - f = io.StringIO() - redirected_stdout = redirect_stdout(f) - redirect_stdout_c() - # Generage a scheme file from the bvals and bvecs files tmp_dir = tempfile.TemporaryDirectory() tmp_scheme_filename = os.path.join(tmp_dir.name, 'gradients.b') @@ -123,8 +122,8 @@ def main(): np.savetxt(tmp_bval_filename, shells_centroids[indices_shells], newline=' ', fmt='%i') fsl2mrtrix(tmp_bval_filename, args.in_bvec, tmp_scheme_filename) - logging.debug('Compute NODDI with AMICO on {} shells at found ' - 'at {}.'.format(len(shells_centroids), shells_centroids)) + logging.info('Compute NODDI with AMICO on {} shells at found ' + 'at {}.'.format(len(shells_centroids), shells_centroids)) with redirected_stdout: # Load the data diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index 30058d66b..822d30e7d 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -97,6 +97,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_AD, args.in_FA, args.in_MD, args.in_RD]) @@ -109,9 +110,6 @@ def main(): assert_same_resolution([args.in_AD, args.in_FA, args.in_MD, args.in_RD]) - log_level = logging.DEBUG if args.verbose else logging.INFO - logging.getLogger().setLevel(log_level) - fa_img = nib.load(args.in_FA) fa_data = fa_img.get_fdata(dtype=np.float32) affine = fa_img.affine @@ -149,13 +147,13 @@ def main(): max(int(cj - w), 0): min(int(cj + w), fa_shape[1]), max(int(ck - w), 0): min(int(ck + w), fa_shape[2])] - logging.debug('fa_min, fa_max, md_min: {}, {}, {}'.format( + logging.info('fa_min, fa_max, md_min: {}, {}, {}'.format( args.fa_min, args.fa_max, args.md_min)) indices = np.where((roi_fa > args.fa_min) & (roi_fa < 0.95)) N = roi_ad[indices].shape[0] - logging.debug('Number of voxels found in single fiber area: {}'.format(N)) + logging.info('Number of voxels found in single fiber area: {}'.format(N)) cc_avg_para = np.mean(roi_ad[indices]) cc_std_para = np.std(roi_ad[indices]) @@ -171,7 +169,7 @@ def main(): indices = np.where((roi_md > args.md_min) & (roi_fa < args.fa_max)) N = roi_md[indices].shape[0] - logging.debug('Number of voxels found in ventricles: {}'.format(N)) + logging.info('Number of voxels found in ventricles: {}'.format(N)) vent_avg = np.mean(roi_md[indices]) vent_std = np.std(roi_md[indices]) diff --git a/scripts/scil_aodf_metrics.py b/scripts/scil_aodf_metrics.py index 7961403c1..cdde59580 100755 --- a/scripts/scil_aodf_metrics.py +++ b/scripts/scil_aodf_metrics.py @@ -28,6 +28,7 @@ import argparse +import logging import nibabel as nib import numpy as np @@ -117,6 +118,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.asi_map = args.asi_map or 'asi_map.nii.gz' diff --git a/scripts/scil_bids_validate.py b/scripts/scil_bids_validate.py index 43c22ea51..375ae4699 100755 --- a/scripts/scil_bids_validate.py +++ b/scripts/scil_bids_validate.py @@ -453,14 +453,12 @@ def associate_dwis(layout, nSub): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + coloredlogs.install(level=logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [], args.bids_ignore) assert_outputs_exist(parser, args, args.out_json) - if args.verbose: - logging.getLogger().setLevel(logging.INFO) - coloredlogs.install(level=logging.INFO) - data = [] bids_indexer = BIDSLayoutIndexer(validate=False, ignore=_load_bidsignore_(os.path.abspath(args.in_bids), diff --git a/scripts/scil_bingham_metrics.py b/scripts/scil_bingham_metrics.py index 1d4b5829e..d4c2eabf0 100755 --- a/scripts/scil_bingham_metrics.py +++ b/scripts/scil_bingham_metrics.py @@ -78,8 +78,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.out_fd = args.out_fd or 'fd.nii.gz' diff --git a/scripts/scil_btensor_metrics.py b/scripts/scil_btensor_metrics.py index 3a1c5963f..bc341da44 100755 --- a/scripts/scil_btensor_metrics.py +++ b/scripts/scil_btensor_metrics.py @@ -144,11 +144,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.basicConfig(level=logging.DEBUG) - else: - logging.basicConfig(level=logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.md = args.md or 'md.nii.gz' diff --git a/scripts/scil_bundle_compute_centroid.py b/scripts/scil_bundle_compute_centroid.py index 4daa6b322..19adb053c 100755 --- a/scripts/scil_bundle_compute_centroid.py +++ b/scripts/scil_bundle_compute_centroid.py @@ -8,6 +8,7 @@ """ import argparse +import logging from dipy.io.stateful_tractogram import StatefulTractogram from dipy.io.streamline import save_tractogram @@ -43,6 +44,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundle) assert_outputs_exist(parser, args, args.out_centroid) diff --git a/scripts/scil_bundle_compute_endpoints_map.py b/scripts/scil_bundle_compute_endpoints_map.py index 725d5845f..f42f0af7b 100755 --- a/scripts/scil_bundle_compute_endpoints_map.py +++ b/scripts/scil_bundle_compute_endpoints_map.py @@ -66,6 +66,7 @@ def main(): parser = _build_arg_parser() args = parser.parse_args() swap = args.swap + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundle, args.reference) assert_outputs_exist(parser, args, [args.endpoints_map_head, diff --git a/scripts/scil_bundle_diameter.py b/scripts/scil_bundle_diameter.py index 37c8ebbbb..9be91fed4 100755 --- a/scripts/scil_bundle_diameter.py +++ b/scripts/scil_bundle_diameter.py @@ -29,6 +29,7 @@ import argparse import json +import logging import os from dipy.io.utils import is_header_compatible @@ -206,6 +207,7 @@ def fit_circle_in_space(positions, directions, dist_w=None): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # The number of labels maps must be equal to the number of bundles tmp = args.in_bundles + args.in_labels diff --git a/scripts/scil_bundle_filter_by_occurence.py b/scripts/scil_bundle_filter_by_occurence.py index 015d83b72..f7116578f 100755 --- a/scripts/scil_bundle_filter_by_occurence.py +++ b/scripts/scil_bundle_filter_by_occurence.py @@ -15,6 +15,7 @@ import argparse +import logging from dipy.io.stateful_tractogram import StatefulTractogram from dipy.io.streamline import save_tractogram @@ -67,6 +68,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundles) output_streamlines_filename = '{}streamlines.trk'.format( diff --git a/scripts/scil_bundle_generate_priors.py b/scripts/scil_bundle_generate_priors.py index 136d28f24..e9636188e 100755 --- a/scripts/scil_bundle_generate_priors.py +++ b/scripts/scil_bundle_generate_priors.py @@ -71,9 +71,9 @@ def _build_arg_parser(): def main(): - logging.getLogger().setLevel(logging.INFO) parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) required = [args.in_bundle, args.in_fodf, args.in_mask] assert_inputs_exist(parser, required) diff --git a/scripts/scil_bundle_label_map.py b/scripts/scil_bundle_label_map.py index 476a9a2ef..b803fb7a5 100755 --- a/scripts/scil_bundle_label_map.py +++ b/scripts/scil_bundle_label_map.py @@ -15,12 +15,12 @@ """ import argparse +import logging import os from dipy.align.streamlinear import StreamlineLinearRegistration from dipy.io.streamline import save_tractogram -from dipy.io.stateful_tractogram import (StatefulTractogram, - set_sft_logger_level) +from dipy.io.stateful_tractogram import StatefulTractogram from dipy.io.utils import is_header_compatible import matplotlib.pyplot as plt import nibabel as nib @@ -85,7 +85,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - set_sft_logger_level('ERROR') + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + assert_inputs_exist(parser, args.in_bundles + [args.in_centroid], optional=args.reference) assert_output_dirs_exist_and_empty(parser, args, args.out_dir) diff --git a/scripts/scil_bundle_mean_fixel_afd.py b/scripts/scil_bundle_mean_fixel_afd.py index 7e59f14b0..1a12ebd50 100755 --- a/scripts/scil_bundle_mean_fixel_afd.py +++ b/scripts/scil_bundle_mean_fixel_afd.py @@ -14,6 +14,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -59,6 +60,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bundle, args.in_fodf]) assert_outputs_exist(parser, args, [args.afd_mean_map]) diff --git a/scripts/scil_bundle_mean_fixel_afd_from_hdf5.py b/scripts/scil_bundle_mean_fixel_afd_from_hdf5.py index 2aef3f28a..5512a861e 100755 --- a/scripts/scil_bundle_mean_fixel_afd_from_hdf5.py +++ b/scripts/scil_bundle_mean_fixel_afd_from_hdf5.py @@ -15,8 +15,9 @@ import argparse import itertools -import os +import logging import multiprocessing +import os import shutil from dipy.io.stateful_tractogram import Space, Origin, StatefulTractogram @@ -97,6 +98,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_hdf5, args.in_fodf]) assert_outputs_exist(parser, args, [args.out_hdf5]) diff --git a/scripts/scil_bundle_mean_fixel_bingham_metric.py b/scripts/scil_bundle_mean_fixel_bingham_metric.py index efe3021da..469dc3abc 100755 --- a/scripts/scil_bundle_mean_fixel_bingham_metric.py +++ b/scripts/scil_bundle_mean_fixel_bingham_metric.py @@ -26,6 +26,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -70,6 +71,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bundle, args.in_bingham, diff --git a/scripts/scil_bundle_mean_std.py b/scripts/scil_bundle_mean_std.py index 12f1ee4b1..91a453072 100755 --- a/scripts/scil_bundle_mean_std.py +++ b/scripts/scil_bundle_mean_std.py @@ -81,6 +81,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bundle] + args.in_metrics, optional=[args.distance_weighting, args.in_labels, diff --git a/scripts/scil_bundle_pairwise_comparison.py b/scripts/scil_bundle_pairwise_comparison.py index 93091de5f..dc4961ad8 100755 --- a/scripts/scil_bundle_pairwise_comparison.py +++ b/scripts/scil_bundle_pairwise_comparison.py @@ -284,6 +284,7 @@ def compute_all_measures(args): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundles) assert_outputs_exist(parser, args, [args.out_json]) diff --git a/scripts/scil_bundle_score_many_bundles_one_tractogram.py b/scripts/scil_bundle_score_many_bundles_one_tractogram.py index 24672bdd2..cf29f4505 100755 --- a/scripts/scil_bundle_score_many_bundles_one_tractogram.py +++ b/scripts/scil_bundle_score_many_bundles_one_tractogram.py @@ -234,9 +234,7 @@ def read_config_file(args): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) (bundle_names, gt_masks, dimensions, vb_sft_list, wpc_sft_list, ib_sft_list, nc_sft, diff --git a/scripts/scil_bundle_score_same_bundle_many_segmentations.py b/scripts/scil_bundle_score_same_bundle_many_segmentations.py index 77a9de24b..1ea87417a 100755 --- a/scripts/scil_bundle_score_same_bundle_many_segmentations.py +++ b/scripts/scil_bundle_score_same_bundle_many_segmentations.py @@ -157,9 +157,7 @@ def compute_streamlines_measures(args): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundles) assert_outputs_exist(parser, args, args.out_json) diff --git a/scripts/scil_bundle_shape_measures.py b/scripts/scil_bundle_shape_measures.py index 2feb10896..1f9928a4c 100755 --- a/scripts/scil_bundle_shape_measures.py +++ b/scripts/scil_bundle_shape_measures.py @@ -197,6 +197,8 @@ def compute_span(streamline_coords): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + assert_inputs_exist(parser, args.in_bundles) assert_outputs_exist(parser, args, [], args.out_json) diff --git a/scripts/scil_bundle_volume_per_label.py b/scripts/scil_bundle_volume_per_label.py index fe946cb14..9b114525d 100755 --- a/scripts/scil_bundle_volume_per_label.py +++ b/scripts/scil_bundle_volume_per_label.py @@ -18,6 +18,7 @@ import argparse import json +import logging import nibabel as nib import numpy as np @@ -48,6 +49,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.voxel_label_map) diff --git a/scripts/scil_clean_qbx_clusters.py b/scripts/scil_clean_qbx_clusters.py index 24afc6570..5a973d78a 100755 --- a/scripts/scil_clean_qbx_clusters.py +++ b/scripts/scil_clean_qbx_clusters.py @@ -149,6 +149,7 @@ def keypress_callback(obj, _): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundles) assert_outputs_exist(parser, args, [args.out_accepted, args.out_rejected]) @@ -162,9 +163,6 @@ def keypress_callback(obj, _): args.out_rejected_dir, create_dir=True) - if args.verbose: - logging.getLogger().setLevel(logging.INFO) - if args.min_cluster_size < 1: parser.error('Minimum cluster size must be at least 1.') diff --git a/scripts/scil_compute_pca.py b/scripts/scil_compute_pca.py index 32313749a..d9f16b2f7 100755 --- a/scripts/scil_compute_pca.py +++ b/scripts/scil_compute_pca.py @@ -180,9 +180,7 @@ def apply_binary_mask(dictionary, mask): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_output_dirs_exist_and_empty(parser, args, args.out_folder, create_dir=True) diff --git a/scripts/scil_connectivity_compare_populations.py b/scripts/scil_connectivity_compare_populations.py index d91484764..74a9bcced 100755 --- a/scripts/scil_connectivity_compare_populations.py +++ b/scripts/scil_connectivity_compare_populations.py @@ -92,14 +92,12 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_g1+args.in_g2, args.filtering_mask) assert_outputs_exist(parser, args, args.out_pval_matrix) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - if args.filtering_mask: filtering_mask = load_matrix_in_any_format(args.filtering_mask) else: @@ -131,7 +129,7 @@ def main(): p_thresh = float(args.p_threshold[0]) matrix_shape = matrices_g1.shape[0:2] masked_pval_matrix = np.zeros(matrix_shape) - logging.debug('Threshold the p-values at {}'.format(p_thresh)) + logging.info('Threshold the p-values at {}'.format(p_thresh)) masked_pval_matrix[matrix_pval < p_thresh] = 1 masked_pval_matrix[matrix_pval < 0] = 0 diff --git a/scripts/scil_connectivity_compute_matrices.py b/scripts/scil_connectivity_compute_matrices.py index 27c920fc0..1f866ff67 100755 --- a/scripts/scil_connectivity_compute_matrices.py +++ b/scripts/scil_connectivity_compute_matrices.py @@ -43,8 +43,8 @@ import argparse import copy import itertools -import multiprocessing import logging +import multiprocessing import os import coloredlogs @@ -285,14 +285,12 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + coloredlogs.install(level=logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_hdf5, args.in_labels], args.force_labels_list) - log_level = logging.INFO if args.verbose else logging.WARNING - logging.getLogger().setLevel(log_level) - coloredlogs.install(level=log_level) - # Summarizing all options chosen by user in measures_to_compute. measures_to_compute = [] measures_output_filename = [] diff --git a/scripts/scil_connectivity_filter.py b/scripts/scil_connectivity_filter.py index 1900495b3..024501402 100755 --- a/scripts/scil_connectivity_filter.py +++ b/scripts/scil_connectivity_filter.py @@ -81,15 +81,13 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_outputs_exist(parser, args, args.out_matrix_mask) if not args.lower_than and not args.greater_than: parser.error('At least one of the two options is required.') - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - conditions_list = [] if args.lower_than: for input_list in args.lower_than: @@ -120,8 +118,8 @@ def main(): population_score = np.sum(empty_matrices, axis=2) - logging.debug('Condition {}_than (#{}) resulted in {} filtered ' - 'elements out of {}.'.format( + logging.info('Condition {}_than (#{}) resulted in {} filtered ' + 'elements out of {}.'.format( condition, condition_counter, len(np.where(population_score < @@ -149,9 +147,9 @@ def main(): 'filtering.\nApply threshold manually to binarize the ' 'output matrix.') else: - logging.debug('All condition resulted in {} filtered ' - 'elements out of {}.'.format(filtered_elem, - np.prod(shape))) + logging.info('All condition resulted in {} filtered ' + 'elements out of {}.'.format(filtered_elem, + np.prod(shape))) save_matrix_in_any_format(args.out_matrix_mask, output_mask.astype(np.uint8)) diff --git a/scripts/scil_connectivity_graph_measures.py b/scripts/scil_connectivity_graph_measures.py index 4a2122e01..de3e691ec 100755 --- a/scripts/scil_connectivity_graph_measures.py +++ b/scripts/scil_connectivity_graph_measures.py @@ -85,19 +85,17 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_length_matrix, args.in_conn_matrix]) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - if not args.append_json: assert_outputs_exist(parser, args, args.out_json) else: - logging.debug('Using --append_json, make sure to delete {} ' - 'before re-launching a group analysis.'.format( - args.out_json)) + logging.info('Using --append_json, make sure to delete {} ' + 'before re-launching a group analysis.'.format( + args.out_json)) if args.append_json and args.overwrite: parser.error('Cannot use the append option at the same time as ' diff --git a/scripts/scil_connectivity_hdf5_average_density_map.py b/scripts/scil_connectivity_hdf5_average_density_map.py index 758650e03..a343cd825 100755 --- a/scripts/scil_connectivity_hdf5_average_density_map.py +++ b/scripts/scil_connectivity_hdf5_average_density_map.py @@ -21,6 +21,7 @@ import argparse import itertools +import logging import multiprocessing import os @@ -95,6 +96,7 @@ def _average_wrapper(args): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_hdf5) assert_output_dirs_exist_and_empty(parser, args, args.out_dir, diff --git a/scripts/scil_connectivity_math.py b/scripts/scil_connectivity_math.py index 73db6dd49..85ccbb3d6 100755 --- a/scripts/scil_connectivity_math.py +++ b/scripts/scil_connectivity_math.py @@ -85,9 +85,7 @@ def load_matrix(arg): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_outputs_exist(parser, args, args.out_matrix) diff --git a/scripts/scil_connectivity_normalize.py b/scripts/scil_connectivity_normalize.py index 9998f20e6..fd5654267 100755 --- a/scripts/scil_connectivity_normalize.py +++ b/scripts/scil_connectivity_normalize.py @@ -45,6 +45,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -108,6 +109,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_matrix, [args.length, args.inverse_length, diff --git a/scripts/scil_connectivity_pairwise_agreement.py b/scripts/scil_connectivity_pairwise_agreement.py index b29eef6c0..41e62687d 100755 --- a/scripts/scil_connectivity_pairwise_agreement.py +++ b/scripts/scil_connectivity_pairwise_agreement.py @@ -13,6 +13,7 @@ import argparse import itertools import json +import logging import numpy as np @@ -50,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_matrices) assert_outputs_exist(parser, args, args.out_json) diff --git a/scripts/scil_connectivity_print_filenames.py b/scripts/scil_connectivity_print_filenames.py index c0836771c..e19222bb5 100755 --- a/scripts/scil_connectivity_print_filenames.py +++ b/scripts/scil_connectivity_print_filenames.py @@ -20,6 +20,7 @@ """ import argparse +import logging import numpy as np @@ -51,6 +52,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_matrix) assert_outputs_exist(parser, args, args.out_txt) diff --git a/scripts/scil_connectivity_reorder_rois.py b/scripts/scil_connectivity_reorder_rois.py index baa0f8886..3a92fefbc 100755 --- a/scripts/scil_connectivity_reorder_rois.py +++ b/scripts/scil_connectivity_reorder_rois.py @@ -24,6 +24,7 @@ """ import argparse +import logging import os import numpy as np @@ -97,6 +98,7 @@ def parse_ordering(in_ordering_file, labels_list=None): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_matrices, [args.labels_list, args.in_ordering]) diff --git a/scripts/scil_convert_tensors.py b/scripts/scil_convert_tensors.py index ea3f42051..95af331c8 100755 --- a/scripts/scil_convert_tensors.py +++ b/scripts/scil_convert_tensors.py @@ -7,6 +7,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -44,6 +45,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_file) assert_outputs_exist(parser, args, args.out_file) diff --git a/scripts/scil_denoising_nlmeans.py b/scripts/scil_denoising_nlmeans.py index 8588cfce4..f027344d7 100755 --- a/scripts/scil_denoising_nlmeans.py +++ b/scripts/scil_denoising_nlmeans.py @@ -75,13 +75,11 @@ def _get_basic_sigma(data): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) assert_outputs_exist(parser, args, args.out_image, args.logfile) - log_level = logging.INFO if args.verbose else logging.WARNING - logging.getLogger().setLevel(log_level) - if args.logfile is not None: logging.getLogger().addHandler(logging.FileHandler(args.logfile, mode='w')) diff --git a/scripts/scil_dki_metrics.py b/scripts/scil_dki_metrics.py index 109b54889..a2daf2127 100755 --- a/scripts/scil_dki_metrics.py +++ b/scripts/scil_dki_metrics.py @@ -146,11 +146,9 @@ def _build_arg_parser(): def main(): - logger = logging.getLogger("Compute_DKI_Metrics") - logger.setLevel(logging.INFO) - parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.dki_fa = args.dki_fa or 'dki_fa.nii.gz' diff --git a/scripts/scil_dti_metrics.py b/scripts/scil_dti_metrics.py index 7b496f6b6..723b447fd 100755 --- a/scripts/scil_dti_metrics.py +++ b/scripts/scil_dti_metrics.py @@ -148,6 +148,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.fa = args.fa or 'fa.nii.gz' diff --git a/scripts/scil_dwi_apply_bias_field.py b/scripts/scil_dwi_apply_bias_field.py index d4a7af611..a12c78b6e 100755 --- a/scripts/scil_dwi_apply_bias_field.py +++ b/scripts/scil_dwi_apply_bias_field.py @@ -10,6 +10,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -48,6 +49,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bias_field], args.mask) assert_outputs_exist(parser, args, args.out_name) diff --git a/scripts/scil_dwi_compute_snr.py b/scripts/scil_dwi_compute_snr.py index 46a6ff176..fe5c84d3f 100755 --- a/scripts/scil_dwi_compute_snr.py +++ b/scripts/scil_dwi_compute_snr.py @@ -90,12 +90,12 @@ def _build_arg_parser(): def main(): - parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: + if args.verbose == "WARNING": logging.getLogger().setLevel(logging.INFO) + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec, args.in_mask], @@ -125,8 +125,7 @@ def main(): noise_mask=noise_mask, noise_map=noise_map, split_shells=args.split_shells, - basename=basename, - verbose=args.verbose) + basename=basename) df = pd.DataFrame.from_dict(values).T diff --git a/scripts/scil_dwi_concatenate.py b/scripts/scil_dwi_concatenate.py index 58f9e835a..1ae54bc5c 100755 --- a/scripts/scil_dwi_concatenate.py +++ b/scripts/scil_dwi_concatenate.py @@ -9,6 +9,7 @@ """ import argparse +import logging from dipy.io.gradients import read_bvals_bvecs from dipy.io.utils import is_header_compatible @@ -53,6 +54,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if len(args.in_dwis) != len(args.in_bvals) \ or len(args.in_dwis) != len(args.in_bvecs): diff --git a/scripts/scil_dwi_convert_FDF.py b/scripts/scil_dwi_convert_FDF.py index ece390e94..538ad222e 100755 --- a/scripts/scil_dwi_convert_FDF.py +++ b/scripts/scil_dwi_convert_FDF.py @@ -13,6 +13,7 @@ """ import argparse +import logging from scilpy.io.varian_fdf import (correct_procpar_intensity, load_fdf, save_babel) @@ -53,6 +54,7 @@ def build_arg_parser(): def main(): parser = build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_outputs_exist(parser, args, args.out_path, optional=[args.bval, args.bvec]) diff --git a/scripts/scil_dwi_detect_volume_outliers.py b/scripts/scil_dwi_detect_volume_outliers.py index c5af9833c..0e9d8df00 100755 --- a/scripts/scil_dwi_detect_volume_outliers.py +++ b/scripts/scil_dwi_detect_volume_outliers.py @@ -15,6 +15,7 @@ """ import argparse +import logging from dipy.io.gradients import read_bvals_bvecs import nibabel as nib @@ -56,6 +57,10 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + if args.verbose == "WARNING": + logging.getLogger().setLevel(logging.INFO) + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) @@ -66,8 +71,7 @@ def main(): bvals.min(), args.b0_thr) bvecs = normalize_bvecs(bvecs) - detect_volume_outliers(data, bvecs, bvals, args.std_scale, - args.verbose, b0_thr) + detect_volume_outliers(data, bvecs, bvals, args.std_scale, b0_thr) if __name__ == "__main__": diff --git a/scripts/scil_dwi_extract_b0.py b/scripts/scil_dwi_extract_b0.py index 472a91521..9e6cf2f70 100755 --- a/scripts/scil_dwi_extract_b0.py +++ b/scripts/scil_dwi_extract_b0.py @@ -88,8 +88,7 @@ def _split_time_steps(b0, affine, header, output): def main(): parser = _build_arg_parser() args = parser.parse_args() - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) diff --git a/scripts/scil_dwi_extract_shell.py b/scripts/scil_dwi_extract_shell.py index f5ee73f76..eb173a335 100755 --- a/scripts/scil_dwi_extract_shell.py +++ b/scripts/scil_dwi_extract_shell.py @@ -71,9 +71,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) assert_outputs_exist(parser, args, [args.out_dwi, args.out_bval, diff --git a/scripts/scil_dwi_powder_average.py b/scripts/scil_dwi_powder_average.py index 9b8627fe6..b3d9bdf14 100755 --- a/scripts/scil_dwi_powder_average.py +++ b/scripts/scil_dwi_powder_average.py @@ -30,9 +30,6 @@ from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, assert_outputs_exist, add_verbose_arg) -logger = logging.getLogger("Compute_Powder_Average") -logger.setLevel(logging.INFO) - def _build_arg_parser(): p = argparse.ArgumentParser( @@ -75,6 +72,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + inputs = [args.in_dwi, args.in_bval] if args.mask: inputs.append(args.mask) @@ -82,9 +81,6 @@ def main(): assert_inputs_exist(parser, inputs) assert_outputs_exist(parser, args, args.out_avg) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - img = nib.load(args.in_dwi) data = img.get_fdata(dtype=np.float32) affine = img.affine diff --git a/scripts/scil_dwi_prepare_eddy_command.py b/scripts/scil_dwi_prepare_eddy_command.py index ee187e30c..0ccb32324 100755 --- a/scripts/scil_dwi_prepare_eddy_command.py +++ b/scripts/scil_dwi_prepare_eddy_command.py @@ -117,6 +117,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) try: devnull = open(os.devnull) @@ -127,9 +128,6 @@ def main(): "the command from the FSL library and make sure it is " "available in your path.".format(args.eddy_cmd)) - if args.verbose: - logging.getLogger().setLevel(logging.INFO) - required_args = [args.in_dwi, args.in_bvals, args.in_bvecs, args.in_mask] assert_inputs_exist(parser, required_args) diff --git a/scripts/scil_dwi_prepare_topup_command.py b/scripts/scil_dwi_prepare_topup_command.py index 4c4bad364..b1ce316a0 100755 --- a/scripts/scil_dwi_prepare_topup_command.py +++ b/scripts/scil_dwi_prepare_topup_command.py @@ -75,6 +75,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) try: devnull = open(os.devnull) @@ -85,9 +86,6 @@ def main(): "the command from the FSL library and make sure it is " "available in your path.") - if args.verbose: - logging.getLogger().setLevel(logging.INFO) - required_args = [args.in_forward_b0, args.in_reverse_b0] assert_inputs_exist(parser, required_args) diff --git a/scripts/scil_dwi_reorder_philips.py b/scripts/scil_dwi_reorder_philips.py index 5ac74d715..717c5f4eb 100755 --- a/scripts/scil_dwi_reorder_philips.py +++ b/scripts/scil_dwi_reorder_philips.py @@ -18,7 +18,7 @@ import nibabel as nib import numpy as np -from scilpy.gradients.utils import get_new_order_philips +from scilpy.gradients.utils import get_new_gtab_order from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg, assert_inputs_exist, @@ -56,9 +56,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) required_args = [args.in_dwi, args.in_bvec, args.in_bval, args.in_table] @@ -90,7 +88,7 @@ def main(): bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec) dwi = nib.load(args.in_dwi) - new_index = get_new_order_philips(philips_table, dwi, bvals, bvecs) + new_index = get_new_gtab_order(philips_table, dwi, bvals, bvecs) bvecs = bvecs[new_index] bvals = bvals[new_index] diff --git a/scripts/scil_dwi_split_by_indices.py b/scripts/scil_dwi_split_by_indices.py index 0aee08ded..d9cff7b16 100755 --- a/scripts/scil_dwi_split_by_indices.py +++ b/scripts/scil_dwi_split_by_indices.py @@ -60,9 +60,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) diff --git a/scripts/scil_dwi_to_sh.py b/scripts/scil_dwi_to_sh.py index 3f9a822c0..c2ec20564 100755 --- a/scripts/scil_dwi_to_sh.py +++ b/scripts/scil_dwi_to_sh.py @@ -8,6 +8,7 @@ """ import argparse +import logging from dipy.core.gradients import gradient_table from dipy.io.gradients import read_bvals_bvecs @@ -56,6 +57,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) assert_outputs_exist(parser, args, args.out_sh) diff --git a/scripts/scil_fodf_max_in_ventricles.py b/scripts/scil_fodf_max_in_ventricles.py index cbe22545c..40a070333 100755 --- a/scripts/scil_fodf_max_in_ventricles.py +++ b/scripts/scil_fodf_max_in_ventricles.py @@ -70,14 +70,12 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_fodfs, args.in_fa, args.in_md]) assert_outputs_exist(parser, args, [], [args.max_value_output, args.mask_output]) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - # Load input image img_fODFs = nib.load(args.in_fodfs) fodf = img_fODFs.get_fdata(dtype=np.float32) diff --git a/scripts/scil_fodf_memsmt.py b/scripts/scil_fodf_memsmt.py index c1f0e23de..543a925f5 100755 --- a/scripts/scil_fodf_memsmt.py +++ b/scripts/scil_fodf_memsmt.py @@ -127,11 +127,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.basicConfig(level=logging.DEBUG) - else: - logging.basicConfig(level=logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.wm_out_fODF = args.wm_out_fODF or 'wm_fodf.nii.gz' diff --git a/scripts/scil_fodf_metrics.py b/scripts/scil_fodf_metrics.py index b2002f14d..2aa032f57 100755 --- a/scripts/scil_fodf_metrics.py +++ b/scripts/scil_fodf_metrics.py @@ -34,6 +34,7 @@ """ import argparse +import logging import numpy as np import nibabel as nib @@ -111,6 +112,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.afd_max = args.afd_max or 'afd_max.nii.gz' diff --git a/scripts/scil_fodf_msmt.py b/scripts/scil_fodf_msmt.py index 82aaa8d1b..8143fd0d3 100755 --- a/scripts/scil_fodf_msmt.py +++ b/scripts/scil_fodf_msmt.py @@ -106,7 +106,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.wm_out_fODF = args.wm_out_fODF or 'wm_fodf.nii.gz' diff --git a/scripts/scil_fodf_ssst.py b/scripts/scil_fodf_ssst.py index 1b758b329..9c276b990 100755 --- a/scripts/scil_fodf_ssst.py +++ b/scripts/scil_fodf_ssst.py @@ -66,7 +66,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec, args.frf_file]) diff --git a/scripts/scil_fodf_to_bingham.py b/scripts/scil_fodf_to_bingham.py index 9a9b7663a..d1aa4004c 100755 --- a/scripts/scil_fodf_to_bingham.py +++ b/scripts/scil_fodf_to_bingham.py @@ -82,8 +82,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_sh, args.mask) assert_outputs_exist(parser, args, args.out_bingham) diff --git a/scripts/scil_freewater_maps.py b/scripts/scil_freewater_maps.py index 7e108148a..50a71caf0 100755 --- a/scripts/scil_freewater_maps.py +++ b/scripts/scil_freewater_maps.py @@ -98,6 +98,15 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + # COMMIT has some c-level stdout and non-logging print that cannot + # be easily stopped. Manual redirection of all printed output + if args.verbose == "WARNING": + f = io.StringIO() + redirected_stdout = redirect_stdout(f) + redirect_stdout_c() + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + redirected_stdout = redirect_stdout(sys.stdout) if args.compute_only and not args.save_kernels: parser.error('--compute_only must be used with --save_kernels.') @@ -109,16 +118,6 @@ def main(): args.out_dir, optional=args.save_kernels) - # COMMIT has some c-level stdout and non-logging print that cannot - # be easily stopped. Manual redirection of all printed output - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - redirected_stdout = redirect_stdout(sys.stdout) - else: - f = io.StringIO() - redirected_stdout = redirect_stdout(f) - redirect_stdout_c() - # Generage a scheme file from the bvals and bvecs files tmp_dir = tempfile.TemporaryDirectory() tmp_scheme_filename = os.path.join(tmp_dir.name, 'gradients.b') @@ -130,7 +129,7 @@ def main(): np.savetxt(tmp_bval_filename, shells_centroids[indices_shells], newline=' ', fmt='%i') fsl2mrtrix(tmp_bval_filename, args.in_bvec, tmp_scheme_filename) - logging.debug( + logging.info( 'Compute FreeWater with AMICO on {} shells at found at {}.'.format( len(shells_centroids), shells_centroids)) diff --git a/scripts/scil_frf_mean.py b/scripts/scil_frf_mean.py index 00e6586f1..2b62f62d7 100755 --- a/scripts/scil_frf_mean.py +++ b/scripts/scil_frf_mean.py @@ -38,7 +38,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.frf_files) assert_outputs_exist(parser, args, args.mean_frf) diff --git a/scripts/scil_frf_memsmt.py b/scripts/scil_frf_memsmt.py index 2484b2375..6cb24370b 100755 --- a/scripts/scil_frf_memsmt.py +++ b/scripts/scil_frf_memsmt.py @@ -179,11 +179,7 @@ def buildArgsParser(): def main(): parser = buildArgsParser() args = parser.parse_args() - - if args.verbose: - logging.basicConfig(level=logging.DEBUG) - else: - logging.basicConfig(level=logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [], optional=list(np.concatenate((args.in_dwis, diff --git a/scripts/scil_frf_msmt.py b/scripts/scil_frf_msmt.py index afeb5bf3c..c55093811 100755 --- a/scripts/scil_frf_msmt.py +++ b/scripts/scil_frf_msmt.py @@ -157,12 +157,9 @@ def buildArgsParser(): def main(): - parser = buildArgsParser() args = parser.parse_args() - - log_level = logging.DEBUG if args.verbose else logging.INFO - logging.getLogger().setLevel(log_level) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) assert_outputs_exist(parser, args, [args.out_wm_frf, args.out_gm_frf, diff --git a/scripts/scil_frf_set_diffusivities.py b/scripts/scil_frf_set_diffusivities.py index be564e3af..6b51e483b 100755 --- a/scripts/scil_frf_set_diffusivities.py +++ b/scripts/scil_frf_set_diffusivities.py @@ -14,6 +14,7 @@ import argparse from ast import literal_eval +import logging import numpy as np from scilpy.io.utils import (add_overwrite_arg, @@ -50,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.frf_file) assert_outputs_exist(parser, args, args.output_frf_file) diff --git a/scripts/scil_frf_ssst.py b/scripts/scil_frf_ssst.py index 5ab64117c..2ede8bf3d 100755 --- a/scripts/scil_frf_ssst.py +++ b/scripts/scil_frf_ssst.py @@ -91,9 +91,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - log_level = logging.DEBUG if args.verbose else logging.INFO - logging.getLogger().setLevel(log_level) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) assert_outputs_exist(parser, args, args.frf_file) diff --git a/scripts/scil_gradients_apply_transform.py b/scripts/scil_gradients_apply_transform.py index 4ede2bd7c..36f519fc6 100755 --- a/scripts/scil_gradients_apply_transform.py +++ b/scripts/scil_gradients_apply_transform.py @@ -8,6 +8,7 @@ """ import argparse +import logging from dipy.io.gradients import read_bvals_bvecs import numpy as np @@ -41,6 +42,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bvecs, args.in_transfo]) assert_outputs_exist(parser, args, args.out_bvecs) diff --git a/scripts/scil_gradients_convert.py b/scripts/scil_gradients_convert.py index f2c2d2766..c2209a74d 100755 --- a/scripts/scil_gradients_convert.py +++ b/scripts/scil_gradients_convert.py @@ -45,9 +45,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) input_is_fsl = args.input_fsl diff --git a/scripts/scil_gradients_generate_sampling.py b/scripts/scil_gradients_generate_sampling.py index 14c9d939b..850188fca 100755 --- a/scripts/scil_gradients_generate_sampling.py +++ b/scripts/scil_gradients_generate_sampling.py @@ -104,6 +104,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # ---- Checks out_basename, _ = os.path.splitext(args.out_basename) @@ -124,9 +125,6 @@ def main(): if args.b0_every is not None and args.b0_every <= 0: parser.error("--b0_every must be an integer > 0.") - log_level = logging.DEBUG if args.verbose else logging.INFO - logging.getLogger().setLevel(log_level) - # ---- b-vectors generation # Non-b0 samples: gradient sampling generation # Generates the b-vectors only. We will set their b-values after. diff --git a/scripts/scil_gradients_modify_axes.py b/scripts/scil_gradients_modify_axes.py index a642725c6..760e63526 100755 --- a/scripts/scil_gradients_modify_axes.py +++ b/scripts/scil_gradients_modify_axes.py @@ -8,6 +8,7 @@ Formerly: scil_flip_gradients.py or scil_swap_gradient_axis.py """ import argparse +import logging import os import numpy as np @@ -52,6 +53,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_gradient_sampling_file) assert_outputs_exist(parser, args, args.out_gradient_sampling_file) diff --git a/scripts/scil_gradients_round_bvals.py b/scripts/scil_gradients_round_bvals.py index ea9402e03..848195f2c 100755 --- a/scripts/scil_gradients_round_bvals.py +++ b/scripts/scil_gradients_round_bvals.py @@ -57,9 +57,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bval) assert_outputs_exist(parser, args, args.out_bval) diff --git a/scripts/scil_gradients_validate_correct.py b/scripts/scil_gradients_validate_correct.py index 86e46a42a..41e092911 100755 --- a/scripts/scil_gradients_validate_correct.py +++ b/scripts/scil_gradients_validate_correct.py @@ -76,8 +76,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) inputs = [args.in_bvec, args.in_peaks, args.in_FA] optional = [args.mask, args.peaks_vals] diff --git a/scripts/scil_gradients_validate_correct_eddy.py b/scripts/scil_gradients_validate_correct_eddy.py index 279e9f559..121e68513 100755 --- a/scripts/scil_gradients_validate_correct_eddy.py +++ b/scripts/scil_gradients_validate_correct_eddy.py @@ -10,6 +10,7 @@ """ import argparse +import logging import numpy as np @@ -41,6 +42,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bvec, args.in_bval]) assert_outputs_exist(parser, args, [args.out_bval, args.out_bvec]) diff --git a/scripts/scil_header_print_info.py b/scripts/scil_header_print_info.py index b6fcce713..ab685220b 100755 --- a/scripts/scil_header_print_info.py +++ b/scripts/scil_header_print_info.py @@ -9,6 +9,7 @@ """ import argparse +import logging import pprint import nibabel as nib @@ -34,6 +35,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_file) diff --git a/scripts/scil_header_validate_compatibility.py b/scripts/scil_header_validate_compatibility.py index 3c7cab613..c0c66d16a 100755 --- a/scripts/scil_header_validate_compatibility.py +++ b/scripts/scil_header_validate_compatibility.py @@ -11,7 +11,7 @@ """ import argparse - +import logging from scilpy.io.utils import ( add_reference_arg, @@ -36,6 +36,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_files) is_header_compatible_multiple_files(parser, args.in_files, diff --git a/scripts/scil_json_convert_entries_to_xlsx.py b/scripts/scil_json_convert_entries_to_xlsx.py index 7a84615e7..8d38a2e55 100755 --- a/scripts/scil_json_convert_entries_to_xlsx.py +++ b/scripts/scil_json_convert_entries_to_xlsx.py @@ -10,6 +10,7 @@ import argparse import json +import logging import numpy as np import pandas as pd @@ -475,6 +476,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_json) assert_outputs_exist(parser, args, args.out_xlsx) diff --git a/scripts/scil_json_harmonize_entries.py b/scripts/scil_json_harmonize_entries.py index 6f67c5637..f0e415176 100755 --- a/scripts/scil_json_harmonize_entries.py +++ b/scripts/scil_json_harmonize_entries.py @@ -18,6 +18,7 @@ import argparse from copy import deepcopy import json +import logging from deepdiff import DeepDiff @@ -48,6 +49,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_file) assert_outputs_exist(parser, args, args.out_file) diff --git a/scripts/scil_json_merge_entries.py b/scripts/scil_json_merge_entries.py index dd362cf98..834ec6a38 100755 --- a/scripts/scil_json_merge_entries.py +++ b/scripts/scil_json_merge_entries.py @@ -30,6 +30,7 @@ import argparse import json +import logging import os from scilpy.io.utils import (add_overwrite_arg, add_json_args, @@ -70,6 +71,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_json) assert_outputs_exist(parser, args, args.out_json) diff --git a/scripts/scil_labels_combine.py b/scripts/scil_labels_combine.py index 35e6f029c..161457ba4 100755 --- a/scripts/scil_labels_combine.py +++ b/scripts/scil_labels_combine.py @@ -18,6 +18,7 @@ import argparse +import logging from dipy.io.utils import is_header_compatible import nibabel as nib @@ -76,6 +77,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) image_files = [] indices_per_input_volume = [] diff --git a/scripts/scil_labels_dilate.py b/scripts/scil_labels_dilate.py index 3d44b7c9e..002b4ed33 100755 --- a/scripts/scil_labels_dilate.py +++ b/scripts/scil_labels_dilate.py @@ -70,7 +70,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_file, optional=args.mask) assert_outputs_exist(parser, args, args.out_file) diff --git a/scripts/scil_labels_remove.py b/scripts/scil_labels_remove.py index 54fc86eba..d69b922a5 100755 --- a/scripts/scil_labels_remove.py +++ b/scripts/scil_labels_remove.py @@ -11,6 +11,7 @@ import argparse +import logging import nibabel as nib @@ -48,6 +49,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_labels) assert_outputs_exist(parser, args, args.out_labels) diff --git a/scripts/scil_labels_split_volume_by_ids.py b/scripts/scil_labels_split_volume_by_ids.py index 1894d1624..cc34a3acd 100755 --- a/scripts/scil_labels_split_volume_by_ids.py +++ b/scripts/scil_labels_split_volume_by_ids.py @@ -12,6 +12,7 @@ """ import argparse +import logging import os import nibabel as nib @@ -52,6 +53,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_labels) diff --git a/scripts/scil_labels_split_volume_from_lut.py b/scripts/scil_labels_split_volume_from_lut.py index aed49464e..a1b2b103f 100755 --- a/scripts/scil_labels_split_volume_from_lut.py +++ b/scripts/scil_labels_split_volume_from_lut.py @@ -14,6 +14,7 @@ import argparse import json +import logging import os import nibabel as nib @@ -56,6 +57,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_label) diff --git a/scripts/scil_lesions_info.py b/scripts/scil_lesions_info.py index be8a5e2a9..dcc66356a 100755 --- a/scripts/scil_lesions_info.py +++ b/scripts/scil_lesions_info.py @@ -15,6 +15,7 @@ import argparse import json +import logging import os import nibabel as nib @@ -72,6 +73,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if (not args.bundle) and (not args.bundle_mask) \ and (not args.bundle_labels_map): diff --git a/scripts/scil_mti_adjust_B1_header.py b/scripts/scil_mti_adjust_B1_header.py new file mode 100644 index 000000000..85812a40c --- /dev/null +++ b/scripts/scil_mti_adjust_B1_header.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Correct B1 map header problem. + +""" + +import argparse +import json +import logging + +import nibabel as nib + +from scilpy.io.utils import (add_overwrite_arg, + add_verbose_arg, + assert_inputs_exist, + assert_outputs_exist) +from scilpy.reconst.mti import (adjust_B1_map_header) + + +def _build_arg_parser(): + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) + p.add_argument('in_B1_map', + help='Path to input B1 map file.') + p.add_argument('out_B1_map', + help='Path to output B1 map file.') + p.add_argument('in_B1_json', + help='Json file of the B1 map.') + + add_verbose_arg(p) + add_overwrite_arg(p) + + return p + + +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + + assert_outputs_exist(parser, args, args.out_B1_map) + assert_inputs_exist(parser, (args.in_B1_map, args.in_B1_json)) + + with open(args.in_B1_json) as curr_json: + b1_json = json.load(curr_json) + if 'PhilipsRescaleSlope' in b1_json.keys(): + slope = b1_json['PhilipsRescaleSlope'] + else: + raise ValueError('Rescale slope not in Json file.') + + b1_img = nib.load(args.in_B1_map) + + new_b1_img = adjust_B1_map_header(b1_img, slope) + + nib.save(new_b1_img, args.out_B1_map) + + +if __name__ == '__main__': + main() diff --git a/scripts/scil_mti_maps_MT.py b/scripts/scil_mti_maps_MT.py index 0bcb9819e..2c257280c 100755 --- a/scripts/scil_mti_maps_MT.py +++ b/scripts/scil_mti_maps_MT.py @@ -13,13 +13,12 @@ Different contrasts can be done with an off-resonance pulse to saturating the protons on non-aqueous molecules a frequency irradiation. The MT maps are -obtained using three contrasts: single frequency irradiation (MT-on, saturated -images) and an unsaturated contrast (MT-off); and a T1weighted image as -reference. +obtained using three or four contrasts: a single positive frequency image +and/or a single negative frequency image, and two unsaturated contrasts as +reference. These two references should be acquired with predominant PD +(proton density) and T1 weighting at different excitation flip angles +(a_PD, a_T1) and repetition times (TR_PD, TR_T1). -The output consist in two types of images: - Three contrasts images : MT-off, MT-on and T1weighted images. - MT maps corrected or not for an empiric B1 correction maps. Input Data recommendation: - it is recommended to use dcm2niix (v1.0.20200331) to convert data @@ -33,30 +32,54 @@ - ANTs can be used for the registration steps (http://stnava.github.io/ANTs/) -The output consist in two types of images in two folders : - 1. Contrasts_MT_maps which contains the 2 contrast images - - MT-off.nii.gz : pulses applied at positive frequency - - MT-on.nii.gz : pulses applied at negative frequency - - T1w.nii.gz : anatomical T1 reference images - - - 2. MT_native_maps which contains the 4 myelin maps - - MTR.nii.gz : Magnetization Transfer Ratio map - The MT ratio is a measure reflecting the amount of bound protons. - - - MTsat.nii.gz : Magnetization Transfer saturation map - The MT saturation is a pseudo-quantitative maps representing - the signal change between the bound and free water pools. - ->>> scil_mti_maps_MT.py path/to/output/directory path/to/mask_bin.nii.gz - --in_mtoff path/to/echo*mtoff.nii.gz --in_mton path/to/echo*mton.nii.gz - --in_t1w path/to/echo*T1w.nii.gz - -Formerly: scil_compute_MT_maps.py +The output consists of a MT_native_maps folder containing the 2 myelin maps: + - MTR.nii.gz : Magnetization Transfer Ratio map + The MT ratio is a measure reflecting the amount of bound protons. + - MTsat.nii.gz : Magnetization Transfer saturation map + The MT saturation is a pseudo-quantitative maps representing + the signal change between the bound and free water pools. + +As an option, the Complementary_maps folder contains the following images: + - positive.nii.gz : single positive frequency image + - negative.nii.gz : single negative frequency image + - mtoff_PD.nii.gz : unsaturated proton density weighted image + - mtoff_T1.nii.gz : unsaturated T1 weighted image + - MTsat_sp.nii.gz : MTsat computed from the single positive frequency image + - MTsat_sn.nii.gz : MTsat computed from the single negative frequency image + - R1app.nii.gz : Apparent R1 map computed for MTsat. + - B1_map.nii.gz : B1 map after correction and smoothing (if given). + + +The final maps from MT_native_maps can be corrected for B1+ field + inhomogeneity, using either an empiric method with + --in_B1_map option, suffix *B1_corrected is added for each map. + --B1_correction_method empiric + or a model-based method with + --in_B1_map option, suffix *B1_corrected is added for each map. + --B1_correction_method model_based + --B1_fitValues 1 or 2 .mat files, obtained externally from + https://github.com/TardifLab/OptimizeIHMTimaging/tree/master/b1Correction, + and given in this order: positive frequency saturation, negative frequency + saturation. +For both methods, the nominal value of the B1 map can be set with + --B1_nominal value + + +>>> scil_mti_maps_MT.py path/to/output/directory + --in_mtoff_pd path/to/echo*mtoff.nii.gz + --in_positive path/to/echo*pos.nii.gz --in_negative path/to/echo*neg.nii.gz + --in_mtoff_t1 path/to/echo*T1w.nii.gz --mask path/to/mask_bin.nii.gz + --in_jsons path/to/echo*mtoff.json path/to/echo*T1w.json + +By default, the script uses all the echoes available in the input folder. +If you want to use a single echo, replace the * with the specific number of +the echo. """ import argparse +import logging import os +import sys import nibabel as nib import numpy as np @@ -66,9 +89,14 @@ assert_output_dirs_exist_and_empty) from scilpy.io.image import load_img from scilpy.image.volume_math import concatenate -from scilpy.reconst.mti import (compute_contrasts_maps, - compute_MT_maps, threshold_maps, - apply_B1_correction) +from scilpy.reconst.mti import (adjust_B1_map_intensities, + apply_B1_corr_empiric, + apply_B1_corr_model_based, + compute_ratio_map, + compute_saturation_map, + process_contrast_map, + threshold_map, + smooth_B1_map) EPILOG = """ Helms G, Dathe H, Kallenberg K, Dechent P. High-resolution maps of @@ -83,32 +111,89 @@ def _build_arg_parser(): formatter_class=argparse.RawTextHelpFormatter) p.add_argument('out_dir', help='Path to output folder.') - p.add_argument('in_mask', - help='Path to the T1 binary brain mask. Must be the sum ' - 'of the three tissue probability maps from ' - 'T1 segmentation (GM+WM+CSF).') p.add_argument('--out_prefix', help='Prefix to be used for each output image.') - p.add_argument('--in_B1_map', - help='Path to B1 coregister map to MT contrasts.') + p.add_argument('--mask', + help='Path to the binary brain mask.') + p.add_argument('--extended', action='store_true', + help='If set, outputs the folder Complementary_maps.') p.add_argument('--filtering', action='store_true', help='Gaussian filtering to remove Gibbs ringing. ' 'Not recommended.') - g = p.add_argument_group(title='MT contrasts', description='Path to ' - 'echoes corresponding to contrasts images. All ' - 'constrasts must have the same number of echoes ' - 'and coregistered between them.' - 'Use * to include all echoes.') - g.add_argument("--in_mtoff", nargs='+', - help='Path to all echoes corresponding to the ' - 'no frequency saturation pulse (reference image).') - g.add_argument("--in_mton", nargs='+', + g = p.add_argument_group(title='Contrast maps', description='Path to ' + 'echoes corresponding to contrast images. All ' + 'constrasts must have \nthe same number of ' + 'echoes and coregistered between them. ' + 'Use * to include all echoes. \nThe in_mtoff_pd ' + 'input and at least one of in_positive or ' + 'in_negative are required.') + g.add_argument("--in_positive", nargs='+', + required='--in_negative' not in sys.argv, help='Path to all echoes corresponding to the ' - 'Positive frequency saturation pulse.') - g.add_argument("--in_t1w", nargs='+', + 'positive frequency \nsaturation pulse.') + g.add_argument("--in_negative", nargs='+', + required='--in_positive' not in sys.argv, help='Path to all echoes corresponding to the ' - 'T1-weigthed.') + 'negative frequency \nsaturation pulse.') + g.add_argument("--in_mtoff_pd", nargs='+', required=True, + help='Path to all echoes corresponding to the predominant ' + 'PD \n(proton density) weighting images with no ' + 'saturation pulse.') + g.add_argument("--in_mtoff_t1", nargs='+', + help='Path to all echoes corresponding to the predominant ' + 'T1 \nweighting images with no saturation pulse. This ' + 'one is optional, \nsince it is only needed for the ' + 'calculation of MTsat. \nAcquisition ' + 'parameters should also be set with this image.') + + a = p.add_argument_group(title='Acquisition parameters', + description='Acquisition parameters required ' + 'for MTsat and ihMTsat ' + 'calculation. \nThese are the ' + 'excitation flip angles ' + '(a_PD, a_T1), in DEGREES, and \n' + 'repetition times (TR_PD, TR_T1) of ' + 'the PD and T1 images, in SECONDS. ' + '\nCan be given through json files ' + '(--in_jsons) or directly ' + '(--in_acq_parameters).') + a1 = a.add_mutually_exclusive_group(required='--in_mtoff_t1' in sys.argv) + a1.add_argument('--in_jsons', nargs=2, + metavar=('PD_json', 'T1_json'), + help='Path to MToff PD json file and MToff T1 json file, ' + 'in that order. \nThe acquisition parameters will be ' + 'extracted from these files. \nMust come from a ' + 'Philips acquisition, otherwise, use ' + 'in_acq_parameters.') + a1.add_argument('--in_acq_parameters', nargs=4, type=float, + metavar=('PD flip angle', 'T1 flip angle', + 'PD repetition time', 'T1 repetition time'), + help='Acquisition parameters in that order: flip angle of ' + 'mtoff_PD, \nflip angle of mtoff_T1, repetition time ' + 'of mtoff_PD, \nrepetition time of mtoff_T1') + + b = p.add_argument_group(title='B1 correction') + b.add_argument('--in_B1_map', + help='Path to B1 coregister map to MT contrasts.') + b.add_argument('--B1_correction_method', + choices=['empiric', 'model_based'], default='empiric', + help='Choice of B1 correction method. Choose between ' + 'empiric and model-based. \nNote that the model-based ' + 'method requires a B1 fitvalues file. \nBoth method ' + 'will only correct the saturation measures. ' + '[%(default)s]') + b.add_argument('--B1_fitvalues', nargs='+', + help='Path to B1 fitvalues files obtained externally. ' + 'Should be one .mat \nfile per input MT-on image, ' + 'given in this specific order: \npositive frequency ' + 'saturation, negative frequency saturation.') + b.add_argument('--B1_nominal', default=100, type=float, + help='Nominal value for the B1 map. For Philips, should be ' + '100. [%(default)s]') + b.add_argument('--B1_smooth_dims', default=5, type=int, + help='Dimension of the squared window used for B1 ' + 'smoothing, in number of voxels. [%(default)s]') add_verbose_arg(p) add_overwrite_arg(p) @@ -119,92 +204,201 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - assert_output_dirs_exist_and_empty(parser, args, - os.path.join(args.out_dir, - 'Contrasts_MT_maps'), - os.path.join(args.out_dir, - 'MT_native_maps'), - create_dir=True) + output_dir = os.path.join(args.out_dir, 'MT_native_maps') + if args.extended: + extended_dir = os.path.join(args.out_dir, 'Complementary_maps') + assert_output_dirs_exist_and_empty(parser, args, extended_dir, + output_dir, create_dir=True) + else: + assert_output_dirs_exist_and_empty(parser, args, output_dir, + create_dir=True) # Merge all echos path into a list - maps = [args.in_mtoff, args.in_mton, args.in_t1w] - - maps_flat = (args.in_mtoff + args.in_mton + args.in_t1w) - - jsons = [curr_map.replace('.nii.gz', '.json') - for curr_map in maps_flat] + input_maps = [] + contrast_names = [] + if args.in_positive: + input_maps.append(args.in_positive) + contrast_names.append('positive') + if args.in_negative: + input_maps.append(args.in_negative) + contrast_names.append('negative') + input_maps.append(args.in_mtoff_pd) + contrast_names.append('mtoff_PD') + if args.in_mtoff_t1: + input_maps.append(args.in_mtoff_t1) + contrast_names.append('mtoff_T1') # check data - assert_inputs_exist(parser, jsons + maps_flat) - for curr_map in maps[1:]: - if len(curr_map) != len(maps[0]): + assert_inputs_exist(parser, args.in_mtoff_pd) # Problem with maps_flat... + # cannot verify the not required input. Somehow it breaks the input_maps... + # even if it is not linked at all. WTF. + for curr_map in input_maps[1:]: + if len(curr_map) != len(input_maps[0]): parser.error('Not the same number of echoes per contrast') - - # Set TR and FlipAngle parameters for MT (mtoff contrast) - # and T1w images - parameters = [] - for curr_map in maps[0][0], maps[2][0]: - acq_parameter = get_acq_parameters(curr_map.replace('.nii.gz', - '.json'), - ['RepetitionTime', 'FlipAngle']) - acq_parameter = acq_parameter[0]*1000, acq_parameter[1]*np.pi/180 - parameters.append(acq_parameter) + if len(input_maps[0]) == 1: + single_echo = True + else: + single_echo = False + + if args.in_B1_map and not args.in_mtoff_t1: + logging.warning('No B1 correction was applied because no MTsat or ' + 'ihMTsat can be computed without the in_mtoff_t1.') + + if args.B1_correction_method == 'model_based' and not args.B1_fitvalues: + parser.error('Fitvalues files must be given when choosing the ' + 'model-based B1 correction method. Please use ' + '--B1_fitvalues.') + + # Set TR and FlipAngle parameters + if args.in_acq_parameters: + flip_angles = np.asarray(args.in_acq_parameters[:2]) * np.pi / 180. + rep_times = np.asarray(args.in_acq_parameters[2:]) * 1000 + if rep_times[0] > 10000 or rep_times[1] > 10000: + logging.warning('Given repetition times do not seem to be given ' + 'in seconds. MTsat results might be affected.') + elif args.in_jsons: + rep_times = [] + flip_angles = [] + for curr_json in args.in_jsons: + acq_parameter = get_acq_parameters(curr_json, + ['RepetitionTime', 'FlipAngle']) + if acq_parameter[0] > 10: + logging.warning('Repetition time found in {} does not seem to ' + 'be given in seconds. MTsat results might be ' + 'affected.'.format(curr_json)) + rep_times.append(acq_parameter[0] * 1000) # convert ms. + flip_angles.append(acq_parameter[1] * np.pi / 180.) # convert rad. # Fix issue from the presence of invalide value and division by zero np.seterr(divide='ignore', invalid='ignore') # Define reference image for saving maps - ref_img = nib.load(maps[0][0]) + affine = nib.load(input_maps[0][0]).affine + + # Load B1 image + if args.in_B1_map and args.in_mtoff_t1: + B1_img = nib.load(args.in_B1_map) + B1_map = B1_img.get_fdata(dtype=np.float32) + B1_map = adjust_B1_map_intensities(B1_map, nominal=args.B1_nominal) + B1_map = smooth_B1_map(B1_map, wdims=args.B1_smooth_dims) + if args.B1_correction_method == 'model_based': + # Apply the B1 map to the flip angles for model-based correction + flip_angles[0] *= B1_map + flip_angles[1] *= B1_map + if args.extended: + nib.save(nib.Nifti1Image(B1_map, affine), + os.path.join(extended_dir, "B1_map.nii.gz")) # Define contrasts maps names - contrasts_name = ['mt_off', 'mt_on', 'T1w'] - + contrast_names_og = contrast_names + if args.filtering: + contrast_names = [curr_name + '_filter' + for curr_name in contrast_names] + if single_echo: + contrast_names = [curr_name + '_single_echo' + for curr_name in contrast_names] if args.out_prefix: - contrasts_name = [args.out_prefix + '_' + curr_name - for curr_name in contrasts_name] + contrast_names = [args.out_prefix + '_' + curr_name + for curr_name in contrast_names] # Compute contrasts maps - computed_contrasts = [] - for idx, curr_map in enumerate(maps): + contrast_maps = [] + for idx, curr_map in enumerate(input_maps): input_images = [] for image in curr_map: img, _ = load_img(image) input_images.append(img) merged_curr_map = concatenate(input_images, input_images[0]) - - computed_contrasts.append(compute_contrasts_maps( - merged_curr_map, filtering=args.filtering)) - - nib.save(nib.Nifti1Image(computed_contrasts[idx].astype(np.float32), - ref_img.affine), - os.path.join(args.out_dir, 'Contrasts_MT_maps', - contrasts_name[idx] + '.nii.gz')) - - # Compute and thresold MT maps - MTR, MTsat = compute_MT_maps(computed_contrasts, parameters) - for curr_map in MTR, MTsat: - curr_map = threshold_maps(curr_map, args.in_mask, 0, 100) - if args.in_B1_map: - curr_map = apply_B1_correction(curr_map, args.in_B1_map) - - # Save MT maps - img_name = ['MTR', 'MTsat'] - + contrast_maps.append(process_contrast_map(merged_curr_map, + filtering=args.filtering, + single_echo=single_echo)) + if args.extended: + nib.save(nib.Nifti1Image(contrast_maps[idx].astype(np.float32), + affine), + os.path.join(extended_dir, + contrast_names[idx] + '.nii.gz')) + + # Compute MTR + if 'positive' in contrast_names_og and 'negative' in contrast_names_og: + MTR = compute_ratio_map((contrast_maps[0] + contrast_maps[1]) / 2, + contrast_maps[2]) + else: + MTR = compute_ratio_map(contrast_maps[0], contrast_maps[1]) + + img_name = ['MTR'] + img_data = [MTR] + + # Compute MTsat + if args.in_mtoff_t1: + MTsat_maps = [] + if 'positive' in contrast_names_og: + MTsat_sp, T1app = compute_saturation_map(contrast_maps[0], + contrast_maps[-2], + contrast_maps[-1], + flip_angles, rep_times) + MTsat_maps.append(MTsat_sp) + if 'negative' in contrast_names_og: + MTsat_sn, T1app = compute_saturation_map(contrast_maps[-3], + contrast_maps[-2], + contrast_maps[-1], + flip_angles, rep_times) + MTsat_maps.append(MTsat_sn) + R1app = 1000 / T1app # convert 1/ms to 1/s + if args.extended: + if 'positive' in contrast_names_og: + nib.save(nib.Nifti1Image(MTsat_sp, affine), + os.path.join(extended_dir, "MTsat_positive.nii.gz")) + if 'negative' in contrast_names_og: + nib.save(nib.Nifti1Image(MTsat_sn, affine), + os.path.join(extended_dir, "MTsat_negative.nii.gz")) + nib.save(nib.Nifti1Image(R1app, affine), + os.path.join(extended_dir, "apparent_R1.nii.gz")) + + # Apply model-based B1 correction + if args.in_B1_map and args.B1_correction_method == 'model_based': + for i, MTsat_map in enumerate(MTsat_maps): + MTsat_maps[i] = apply_B1_corr_model_based(MTsat_map, B1_map, + R1app, + args.B1_fitvalues[i]) + + # Compute MTsat and ihMTsat from saturations + if 'positive' in contrast_names_og and 'negative' in contrast_names_og: + MTsat = (MTsat_maps[0] + MTsat_maps[1]) / 2 + else: + MTsat = MTsat_maps[0] + + # Apply empiric B1 correction + if args.in_B1_map and args.B1_correction_method == 'empiric': + # MTR = apply_B1_correction_empiric(MTR, B1_map) + MTsat = apply_B1_corr_empiric(MTsat, B1_map) + + img_name.append('MTsat') + img_data.append(MTsat) + + # Apply thresholds on maps + for i, map in enumerate(img_data): + img_data[i] = threshold_map(map, args.mask, 0, 100) + + # Save ihMT and MT images + if args.filtering: + img_name = [curr_name + '_filter' + for curr_name in img_name] + if single_echo: + img_name = [curr_name + '_single_echo' + for curr_name in img_name] if args.in_B1_map: img_name = [curr_name + '_B1_corrected' for curr_name in img_name] - if args.out_prefix: img_name = [args.out_prefix + '_' + curr_name for curr_name in img_name] - img_data = MTR, MTsat for img_to_save, name in zip(img_data, img_name): nib.save(nib.Nifti1Image(img_to_save.astype(np.float32), - ref_img.affine), - os.path.join(args.out_dir, 'MT_native_maps', - name + '.nii.gz')) + affine), + os.path.join(output_dir, name + '.nii.gz')) if __name__ == '__main__': diff --git a/scripts/scil_mti_maps_ihMT.py b/scripts/scil_mti_maps_ihMT.py index b7af1f67d..72320dca3 100755 --- a/scripts/scil_mti_maps_ihMT.py +++ b/scripts/scil_mti_maps_ihMT.py @@ -14,9 +14,12 @@ Different contrasts can be done with an off-resonance pulse prior to image acquisition (a prepulse), saturating the protons on non-aqueous molecules, by applying different frequency irradiation. The two MT maps and two ihMT maps -are obtained using five contrasts: single frequency positive or negative and -dual frequency with an alternation of both positive and negative frequency -(saturated images); and one unsaturated contrast as reference (T1weighted). +are obtained using six contrasts: single positive frequency image, single +negative frequency image, dual alternating positive/negative frequency image, +dual alternating negative/positive frequency image (saturated images); +and two unsaturated contrasts as reference. These two references should be +acquired with predominant PD (proton density) and T1 weighting at different +excitation flip angles (a_PD, a_T1) and repetition times (TR_PD, TR_T1). Input Data recommendation: @@ -24,48 +27,66 @@ https://github.com/rordenlab/dcm2niix/releases/tag/v1.0.20200331 - dcm2niix conversion will create all echo files for each contrast and corresponding json files - - all input must have a matching json file with the same filename - all contrasts must have a same number of echoes and coregistered - between them before running the script. + between them before running the script - Mask must be coregistered to the echo images - ANTs can be used for the registration steps (http://stnava.github.io/ANTs/) -The output consist in two types of images in two folders : - 1. Contrasts_ihMT_maps which contains the 5 contrast images - - altnp.nii.gz : alternating negative and positive frequency pulses - - altpn.nii.gz : alternating positive and negative frequency pulses - - positive.nii.gz : pulses applied at positive frequency - - negative.nii.gz : pulses applied at negative frequency - - reference.nii.gz : no pulse - - 2. ihMT_native_maps which contains the 4 myelin maps - - MTR.nii.gz : Magnetization Transfer Ratio map - - ihMTR.nii.gz : inhomogeneous Magnetization Transfer Ratio map - The (ih)MT ratio is a measure reflecting the amount of bound protons. - - - MTsat.nii.gz : Magnetization Transfer saturation map - - ihMTsat.nii.gz : inhomogeneous Magnetization Transfer saturation map - The (ih)MT saturation is a pseudo-quantitative maps representing - the signal change between the bound and free water pools. - - These final maps can be corrected by an empiric B1 correction with +The output consists of a ihMT_native_maps folder containing the 4 myelin maps: + - MTR.nii.gz : Magnetization Transfer Ratio map + - ihMTR.nii.gz : inhomogeneous Magnetization Transfer Ratio map + The (ih)MT ratio is a measure reflecting the amount of bound protons. + - MTsat.nii.gz : Magnetization Transfer saturation map + - ihMTsat.nii.gz : inhomogeneous Magnetization Transfer saturation map + The (ih)MT saturation is a pseudo-quantitative maps representing + the signal change between the bound and free water pools. + +As an option, the Complementary_maps folder contains the following images: + - altnp.nii.gz : dual alternating negative and positive frequency image + - altpn.nii.gz : dual alternating positive and negative frequency image + - positive.nii.gz : single positive frequency image + - negative.nii.gz : single negative frequency image + - mtoff_PD.nii.gz : unsaturated proton density weighted image + - mtoff_T1.nii.gz : unsaturated T1 weighted image + - MTsat_d.nii.gz : MTsat computed from the mean dual frequency images + - MTsat_sp.nii.gz : MTsat computed from the single positive frequency image + - MTsat_sn.nii.gz : MTsat computed from the single negative frequency image + - R1app.nii.gz : Apparent R1 map computed for MTsat. + - B1_map.nii.gz : B1 map after correction and smoothing (if given). + + +The final maps from ihMT_native_maps can be corrected for B1+ field + inhomogeneity, using either an empiric method with --in_B1_map option, suffix *B1_corrected is added for each map. + --B1_correction_method empiric + or a model-based method with + --in_B1_map option, suffix *B1_corrected is added for each map. + --B1_correction_method model_based + --B1_fitValues 3 .mat files, obtained externally from + https://github.com/TardifLab/OptimizeIHMTimaging/tree/master/b1Correction, + and given in this order: positive frequency saturation, negative frequency + saturation, dual frequency saturation. +For both methods, the nominal value of the B1 map can be set with + --B1_nominal value + ->>> scil_mti_maps_ihMT.py path/to/output/directory path/to/mask_bin.nii.gz +>>> scil_mti_maps_ihMT.py path/to/output/directory --in_altnp path/to/echo*altnp.nii.gz --in_altpn path/to/echo*altpn.nii.gz - --in_mtoff path/to/echo*mtoff.nii.gz --in_negative path/to/echo*neg.nii.gz - --in_positive path/to/echo*pos.nii.gz --in_t1w path/to/echo*T1w.nii.gz + --in_mtoff_pd path/to/echo*mtoff.nii.gz + --in_negative path/to/echo*neg.nii.gz --in_positive path/to/echo*pos.nii.gz + --in_mtoff_t1 path/to/echo*T1w.nii.gz --mask path/to/mask_bin.nii.gz + --in_jsons path/to/echo*mtoff.json path/to/echo*T1w.json By default, the script uses all the echoes available in the input folder. -If you want to use a single echo add --single_echo to the command line and -replace the * with the specific number of the echo. - -Formerly: scil_compute_ihMT_maps.py +If you want to use a single echo, replace the * with the specific number of +the echo. """ import argparse +import logging import os +import sys import nibabel as nib import numpy as np @@ -75,10 +96,14 @@ assert_output_dirs_exist_and_empty) from scilpy.io.image import load_img from scilpy.image.volume_math import concatenate -from scilpy.reconst.mti import (compute_contrasts_maps, - compute_ihMT_maps, threshold_maps, - compute_MT_maps_from_ihMT, - apply_B1_correction) +from scilpy.reconst.mti import (adjust_B1_map_intensities, + apply_B1_corr_empiric, + apply_B1_corr_model_based, + compute_ratio_map, + compute_saturation_map, + process_contrast_map, + threshold_map, + smooth_B1_map) EPILOG = """ Varma G, Girard OM, Prevost VH, Grant AK, Duhamel G, Alsop DC. @@ -102,45 +127,94 @@ def _build_arg_parser(): formatter_class=argparse.RawTextHelpFormatter) p.add_argument('out_dir', help='Path to output folder.') - p.add_argument('in_mask', - help='Path to the T1 binary brain mask. Must be the sum ' - 'of the three tissue probability maps from ' - 'T1 segmentation (GM+WM+CSF).') + p.add_argument('--out_prefix', help='Prefix to be used for each output image.') - p.add_argument('--in_B1_map', - help='Path to B1 coregister map to MT contrasts.') + p.add_argument('--mask', + help='Path to the binary brain mask.') + p.add_argument('--extended', action='store_true', + help='If set, outputs the folder Complementary_maps.') p.add_argument('--filtering', action='store_true', help='Gaussian filtering to remove Gibbs ringing. ' 'Not recommended.') - p.add_argument('--single_echo', action='store_true', - help='Use this option when there is only one echo.') - g = p.add_argument_group(title='ihMT contrasts', description='Path to ' - 'echoes corresponding to contrasts images. All ' - 'constrasts must have the same number of echoes ' - 'and coregistered between them. ' + g = p.add_argument_group(title='Contrast maps', description='Path to ' + 'echoes corresponding to contrast images. All ' + 'constrasts must have \nthe same number of ' + 'echoes and coregistered between them. ' 'Use * to include all echoes.') g.add_argument('--in_altnp', nargs='+', required=True, help='Path to all echoes corresponding to the ' - 'alternation of Negative and Positive' + 'alternation of \nnegative and positive ' 'frequency saturation pulse.') g.add_argument('--in_altpn', nargs='+', required=True, help='Path to all echoes corresponding to the ' - 'alternation of Positive and Negative ' + 'alternation of \npositive and negative ' 'frequency saturation pulse.') - g.add_argument("--in_mtoff", nargs='+', required=True, - help='Path to all echoes corresponding to the ' - 'no frequency saturation pulse (reference image).') g.add_argument("--in_negative", nargs='+', required=True, help='Path to all echoes corresponding to the ' - 'Negative frequency saturation pulse.') + 'negative frequency \nsaturation pulse.') g.add_argument("--in_positive", nargs='+', required=True, help='Path to all echoes corresponding to the ' - 'Positive frequency saturation pulse.') - g.add_argument("--in_t1w", nargs='+', required=True, - help='Path to all echoes corresponding to the ' - 'T1-weigthed.') + 'positive frequency \nsaturation pulse.') + g.add_argument("--in_mtoff_pd", nargs='+', required=True, + help='Path to all echoes corresponding to the predominant ' + 'PD \n(proton density) weighting images with no ' + 'saturation pulse.') + g.add_argument("--in_mtoff_t1", nargs='+', + help='Path to all echoes corresponding to the predominant ' + 'T1 \nweighting images with no saturation pulse. This ' + 'one is optional, \nsince it is only needed for the ' + 'calculation of MTsat and ihMTsat. \nAcquisition ' + 'parameters should also be set with this image.') + + a = p.add_argument_group(title='Acquisition parameters', + description='Acquisition parameters required ' + 'for MTsat and ihMTsat ' + 'calculation. \nThese are the ' + 'excitation flip angles ' + '(a_PD, a_T1), in DEGREES, and \n' + 'repetition times (TR_PD, TR_T1) of ' + 'the PD and T1 images, in SECONDS. ' + '\nCan be given through json files ' + '(--in_jsons) or directly ' + '(--in_acq_parameters).') + a1 = a.add_mutually_exclusive_group(required='--in_mtoff_t1' in sys.argv) + a1.add_argument('--in_jsons', nargs=2, + metavar=('PD_json', 'T1_json'), + help='Path to MToff PD json file and MToff T1 json file, ' + 'in that order. \nThe acquisition parameters will be ' + 'extracted from these files. \nMust come from a ' + 'Philips acquisition, otherwise, use ' + 'in_acq_parameters.') + a1.add_argument('--in_acq_parameters', nargs=4, type=float, + metavar=('PD flip angle', 'T1 flip angle', + 'PD repetition time', 'T1 repetition time'), + help='Acquisition parameters in that order: flip angle of ' + 'mtoff_PD, \nflip angle of mtoff_T1, repetition time ' + 'of mtoff_PD, \nrepetition time of mtoff_T1') + + b = p.add_argument_group(title='B1 correction') + b.add_argument('--in_B1_map', + help='Path to B1 coregister map to MT contrasts.') + b.add_argument('--B1_correction_method', + choices=['empiric', 'model_based'], default='empiric', + help='Choice of B1 correction method. Choose between ' + 'empiric and model-based. \nNote that the model-based ' + 'method requires a B1 fitvalues file. \nBoth method ' + 'will only correct the saturation measures. ' + '[%(default)s]') + b.add_argument('--B1_fitvalues', nargs=3, + help='Path to B1 fitvalues files obtained externally. ' + 'Should be three .mat \nfiles given in this specific ' + 'order: positive frequency saturation, \nnegative ' + 'frequency saturation, dual frequency saturation.') + b.add_argument('--B1_nominal', default=100, type=float, + help='Nominal value for the B1 map. For Philips, should be ' + '100. [%(default)s]') + b.add_argument('--B1_smooth_dims', default=5, type=int, + help='Dimension of the squared window used for B1 ' + 'smoothing, in number of voxels. [%(default)s]') add_verbose_arg(p) add_overwrite_arg(p) @@ -151,123 +225,209 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - assert_output_dirs_exist_and_empty(parser, args, - os.path.join(args.out_dir, - 'Contrasts_ihMT_maps'), - os.path.join(args.out_dir, - 'ihMT_native_maps'), - create_dir=True) + output_dir = os.path.join(args.out_dir, 'ihMT_native_maps') + if args.extended: + extended_dir = os.path.join(args.out_dir, 'Complementary_maps') + assert_output_dirs_exist_and_empty(parser, args, extended_dir, + output_dir, create_dir=True) + else: + assert_output_dirs_exist_and_empty(parser, args, output_dir, + create_dir=True) - # Merge all echos path into a list - maps = [args.in_altnp, args.in_altpn, args.in_mtoff, args.in_negative, - args.in_positive, args.in_t1w] - - maps_flat = (args.in_altnp + args.in_altpn + args.in_mtoff + - args.in_negative + args.in_positive + args.in_t1w) + if args.out_prefix: + out_prefix = args.out_prefix + "_" + else: + out_prefix = "" - jsons = [curr_map.replace('.nii.gz', '.json') - for curr_map in maps_flat] + # Merge all echos path into a list + input_maps = [args.in_altnp, args.in_altpn, args.in_negative, + args.in_positive, args.in_mtoff_pd] + maps_flat = (args.in_altnp + args.in_altpn + args.in_negative + + args.in_positive + args.in_mtoff_pd) + if args.in_mtoff_t1: + input_maps.append(args.in_mtoff_t1) # check echoes number and jsons - assert_inputs_exist(parser, jsons + maps_flat) - for curr_map in maps[1:]: - if len(curr_map) != len(maps[0]): + assert_inputs_exist(parser, maps_flat, optional=args.in_mtoff_t1) + for curr_map in input_maps[1:]: + if len(curr_map) != len(input_maps[0]): parser.error('Not the same number of echoes per contrast') - - # Set TR and FlipAngle parameters for ihMT (positive contrast) - # and T1w images - parameters = [] - for curr_map in maps[4][0], maps[5][0]: - acq_parameter = get_acq_parameters(curr_map.replace('.nii.gz', - '.json'), - ['RepetitionTime', 'FlipAngle']) - acq_parameter = acq_parameter[0]*1000, acq_parameter[1]*np.pi/180 - parameters.append(acq_parameter) + if len(input_maps[0]) == 1: + single_echo = True + else: + single_echo = False + + if args.in_B1_map and not args.in_mtoff_t1: + logging.warning('No B1 correction was applied because no MTsat or ' + 'ihMTsat can be computed without the in_mtoff_t1.') + + if args.B1_correction_method == 'model_based' and not args.B1_fitvalues: + parser.error('Fitvalues files must be given when choosing the ' + 'model-based B1 correction method. Please use ' + '--B1_fitvalues.') + + # Set TR and FlipAngle parameters + if args.in_acq_parameters: + flip_angles = np.asarray(args.in_acq_parameters[:2]) * np.pi / 180. + rep_times = np.asarray(args.in_acq_parameters[2:]) * 1000 + if rep_times[0] > 10000 or rep_times[1] > 10000: + logging.warning('Given repetition times do not seem to be given ' + 'in seconds. MTsat results might be affected.') + elif args.in_jsons: + rep_times = [] + flip_angles = [] + for curr_json in args.in_jsons: + acq_parameter = get_acq_parameters(curr_json, + ['RepetitionTime', 'FlipAngle']) + if acq_parameter[0] > 10: + logging.warning('Repetition time found in {} does not seem to ' + 'be given in seconds. MTsat and ihMTsat ' + 'results might be affected.'.format(curr_json)) + rep_times.append(acq_parameter[0] * 1000) # convert ms. + flip_angles.append(acq_parameter[1] * np.pi / 180.) # convert rad. # Fix issue from the presence of invalide value and division by zero np.seterr(divide='ignore', invalid='ignore') - # Define reference image for saving maps - ref_img = nib.load(maps[4][0]) + # Define affine + affine = nib.load(input_maps[4][0]).affine + + # Load B1 image + if args.in_B1_map and args.in_mtoff_t1: + B1_img = nib.load(args.in_B1_map) + B1_map = B1_img.get_fdata(dtype=np.float32) + B1_map = adjust_B1_map_intensities(B1_map, nominal=args.B1_nominal) + B1_map = smooth_B1_map(B1_map, wdims=args.B1_smooth_dims) + if args.B1_correction_method == 'model_based': + # Apply shift to the B1 map for better correction + # shift = 0.05 + # expt = 1.3 + # B1_map = (B1_map + shift) ** expt + # Apply the B1 map to the flip angles for model-based correction + flip_angles[0] *= B1_map + flip_angles[1] *= B1_map + if args.extended: + nib.save(nib.Nifti1Image(B1_map, affine), + os.path.join(extended_dir, out_prefix + "B1_map.nii.gz")) # Define contrasts maps names - contrasts_name = ['altnp', 'altpn', 'reference', 'negative', 'positive', - 'T1w'] + contrast_names = ['altnp', 'altpn', 'negative', 'positive', 'mtoff_PD', + 'mtoff_T1'] if args.filtering: - contrasts_name = [curr_name + '_filter' - for curr_name in contrasts_name] - if args.single_echo: - contrasts_name = [curr_name + '_single_echo' - for curr_name in contrasts_name] - + contrast_names = [curr_name + '_filter' + for curr_name in contrast_names] + if single_echo: + contrast_names = [curr_name + '_single_echo' + for curr_name in contrast_names] if args.out_prefix: - contrasts_name = [args.out_prefix + '_' + curr_name - for curr_name in contrasts_name] + contrast_names = [out_prefix + curr_name + for curr_name in contrast_names] # Compute contrasts maps - computed_contrasts = [] - for idx, curr_map in enumerate(maps): + contrast_maps = [] + for idx, curr_map in enumerate(input_maps): input_images = [] for image in curr_map: img, _ = load_img(image) input_images.append(img) merged_curr_map = concatenate(input_images, input_images[0]) - computed_contrasts.append(compute_contrasts_maps( - merged_curr_map, filtering=args.filtering, - single_echo=args.single_echo)) - - nib.save(nib.Nifti1Image(computed_contrasts[idx].astype(np.float32), - ref_img.affine), - os.path.join(args.out_dir, 'Contrasts_ihMT_maps', - contrasts_name[idx] + '.nii.gz')) - - # Compute and thresold ihMT maps - ihMTR, ihMTsat = compute_ihMT_maps(computed_contrasts, parameters) - ihMTR = threshold_maps(ihMTR, args.in_mask, 0, 100, - idx_contrast_list=[4, 3, 1, 0, 2], - contrasts_maps=computed_contrasts) - ihMTsat = threshold_maps(ihMTsat, args.in_mask, 0, 10, - idx_contrast_list=[4, 3, 1, 0], - contrasts_maps=computed_contrasts) - if args.in_B1_map: - ihMTR = apply_B1_correction(ihMTR, args.in_B1_map) - ihMTsat = apply_B1_correction(ihMTsat, args.in_B1_map) - - # Compute and thresold non-ihMT maps - MTR, MTsat = compute_MT_maps_from_ihMT(computed_contrasts, parameters) - for curr_map in MTR, MTsat: - curr_map = threshold_maps(curr_map, args.in_mask, 0, 100, - idx_contrast_list=[4, 2], - contrasts_maps=computed_contrasts) - if args.in_B1_map: - curr_map = apply_B1_correction(curr_map, args.in_B1_map) + contrast_maps.append(process_contrast_map(merged_curr_map, + filtering=args.filtering, + single_echo=single_echo)) + if args.extended: + nib.save(nib.Nifti1Image(contrast_maps[idx].astype(np.float32), + affine), + os.path.join(extended_dir, + contrast_names[idx] + '.nii.gz')) + + # Compute ratio maps + MTR, ihMTR = compute_ratio_map((contrast_maps[2] + contrast_maps[3]) / 2, + contrast_maps[4], + mt_on_dual=(contrast_maps[0] + + contrast_maps[1]) / 2) + img_name = ['ihMTR', 'MTR'] + img_data = [ihMTR, MTR] + + # Compute saturation maps + if args.in_mtoff_t1: + MTsat_sp, T1app = compute_saturation_map(contrast_maps[3], + contrast_maps[4], + contrast_maps[5], + flip_angles, rep_times) + MTsat_sn, _ = compute_saturation_map(contrast_maps[2], + contrast_maps[4], + contrast_maps[5], + flip_angles, rep_times) + MTsat_d, _ = compute_saturation_map((contrast_maps[0] + + contrast_maps[1]) / 2, + contrast_maps[4], contrast_maps[5], + flip_angles, rep_times) + R1app = 1000 / T1app # convert 1/ms to 1/s + if args.extended: + nib.save(nib.Nifti1Image(MTsat_sp, affine), + os.path.join(extended_dir, + out_prefix + "MTsat_single_positive.nii.gz")) + nib.save(nib.Nifti1Image(MTsat_sn, affine), + os.path.join(extended_dir, + out_prefix + "MTsat_single_negative.nii.gz")) + nib.save(nib.Nifti1Image(MTsat_d, affine), + os.path.join(extended_dir, + out_prefix + "MTsat_dual.nii.gz")) + nib.save(nib.Nifti1Image(R1app, affine), + os.path.join(extended_dir, + out_prefix + "apparent_R1.nii.gz")) + + MTsat_maps = [MTsat_sp, MTsat_sn, MTsat_d] + + # Apply model-based B1 correction + if args.in_B1_map and args.B1_correction_method == 'model_based': + for i, MTsat_map in enumerate(MTsat_maps): + MTsat_maps[i] = apply_B1_corr_model_based(MTsat_map, B1_map, + R1app, + args.B1_fitvalues[i]) + + # Compute MTsat and ihMTsat from saturations + MTsat = (MTsat_maps[0] + MTsat_maps[1]) / 2 + ihMTsat = MTsat_maps[2] - MTsat + + # Apply empiric B1 correction + if args.in_B1_map and args.B1_correction_method == 'empiric': + # MTR = apply_B1_correction_empiric(MTR, B1_map) + # ihMTR = apply_B1_correction_empiric(ihMTR, B1_map) + MTsat = apply_B1_corr_empiric(MTsat, B1_map) + ihMTsat = apply_B1_corr_empiric(ihMTsat, B1_map) + + img_name.extend(('ihMTsat', 'MTsat')) + img_data.extend((ihMTsat, MTsat)) + + # Apply thresholds on maps + upper_thresholds = [100, 100, 10, 10] + idx_contrast_lists = [[0, 1, 2, 3, 4], [3, 4], [0, 1, 2, 3], [3, 4]] + for i, map in enumerate(img_data): + img_data[i] = threshold_map(map, args.mask, 0, upper_thresholds[i], + idx_contrast_list=idx_contrast_lists[i], + contrast_maps=contrast_maps) # Save ihMT and MT images - img_name = ['ihMTR', 'ihMTsat', 'MTR', 'MTsat'] - if args.filtering: img_name = [curr_name + '_filter' for curr_name in img_name] - - if args.single_echo: + if single_echo: img_name = [curr_name + '_single_echo' for curr_name in img_name] - if args.in_B1_map: img_name = [curr_name + '_B1_corrected' for curr_name in img_name] - if args.out_prefix: img_name = [args.out_prefix + '_' + curr_name for curr_name in img_name] - img_data = ihMTR, ihMTsat, MTR, MTsat for img_to_save, name in zip(img_data, img_name): nib.save(nib.Nifti1Image(img_to_save.astype(np.float32), - ref_img.affine), - os.path.join(args.out_dir, 'ihMT_native_maps', - name + '.nii.gz')) + affine), + os.path.join(output_dir, name + '.nii.gz')) if __name__ == '__main__': diff --git a/scripts/scil_outlier_rejection.py b/scripts/scil_outlier_rejection.py index 1004a00c2..977df8ad6 100755 --- a/scripts/scil_outlier_rejection.py +++ b/scripts/scil_outlier_rejection.py @@ -51,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundle) assert_outputs_exist(parser, args, args.out_bundle, args.remaining_bundle) diff --git a/scripts/scil_plot_stats_per_point.py b/scripts/scil_plot_stats_per_point.py index d2f230c66..4f7a6e321 100755 --- a/scripts/scil_plot_stats_per_point.py +++ b/scripts/scil_plot_stats_per_point.py @@ -13,6 +13,7 @@ import argparse import itertools import json +import logging import os import numpy as np @@ -62,6 +63,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_json) assert_output_dirs_exist_and_empty(parser, args, args.out_dir, diff --git a/scripts/scil_project_streamlines_to_map.py b/scripts/scil_project_streamlines_to_map.py index f4f96c33f..24b1adaf3 100755 --- a/scripts/scil_project_streamlines_to_map.py +++ b/scripts/scil_project_streamlines_to_map.py @@ -171,6 +171,7 @@ def _pick_data(args, sft): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_bundle], args.in_metrics + args.load_dps + args.load_dpp) diff --git a/scripts/scil_qball_metrics.py b/scripts/scil_qball_metrics.py index 22d36d0c6..817ffd565 100755 --- a/scripts/scil_qball_metrics.py +++ b/scripts/scil_qball_metrics.py @@ -96,6 +96,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not args.not_all: args.gfa = args.gfa or 'gfa.nii.gz' diff --git a/scripts/scil_rgb_convert.py b/scripts/scil_rgb_convert.py index 0eb201318..10b6921fe 100755 --- a/scripts/scil_rgb_convert.py +++ b/scripts/scil_rgb_convert.py @@ -22,6 +22,7 @@ """ import argparse +import logging from dipy.io.utils import decfa, decfa_to_float import nibabel as nib @@ -52,6 +53,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) assert_outputs_exist(parser, args, args.out_image) diff --git a/scripts/scil_score_tractogram.py b/scripts/scil_score_tractogram.py index 9e3d8569e..32d0dc7ba 100755 --- a/scripts/scil_score_tractogram.py +++ b/scripts/scil_score_tractogram.py @@ -401,9 +401,7 @@ def read_config_file(args): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # Load (gt_tails, gt_heads, sft, bundle_names, list_rois, bundle_lengths, angles, diff --git a/scripts/scil_screenshot_bundle.py b/scripts/scil_screenshot_bundle.py index 137bc593f..809740d81 100755 --- a/scripts/scil_screenshot_bundle.py +++ b/scripts/scil_screenshot_bundle.py @@ -99,12 +99,12 @@ def prepare_data_for_actors(bundle_filename, reference_filename, target_template_affine = target_template_img.affine # Register the DWI data to the template - logging.debug('Starting registration...') + logging.info('Starting registration...') transformed_reference, transformation = register_image(target_template_data, target_template_affine, reference_data, reference_affine) - logging.debug('Transforming streamlines...') + logging.info('Transforming streamlines...') streamlines = transform_streamlines(streamlines, np.linalg.inv(transformation), in_place=True) @@ -154,12 +154,11 @@ def plot_glass_brain(args, sft, img, output_filenames): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + required = [args.in_bundle, args.in_anat] assert_inputs_exist(parser, required, args.target_template) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - output_filenames_3d = [] output_filenames_glass = [] for axis_name in ['sagittal', 'coronal', 'axial']: diff --git a/scripts/scil_screenshot_dti.py b/scripts/scil_screenshot_dti.py index a37c7e2d2..cbfbf9833 100755 --- a/scripts/scil_screenshot_dti.py +++ b/scripts/scil_screenshot_dti.py @@ -10,6 +10,7 @@ """ import argparse +import logging import os from dipy.core.gradients import gradient_table, get_bval_indices @@ -143,6 +144,8 @@ def prepare_slices_mask(mask_data, x_slice, y_slice, z_slice): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + required = [args.in_dwi, args.in_bval, args.in_bvec, args.in_template] assert_inputs_exist(parser, required) diff --git a/scripts/scil_screenshot_volume.py b/scripts/scil_screenshot_volume.py index 983c19ada..646a40c39 100755 --- a/scripts/scil_screenshot_volume.py +++ b/scripts/scil_screenshot_volume.py @@ -62,6 +62,7 @@ """ import argparse +import logging from itertools import zip_longest import numpy as np @@ -124,6 +125,7 @@ def _parse_args(parser): def main(): parser = _build_arg_parser() args = _parse_args(parser) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) vol_img, t_mask_img, labelmap_img, mask_imgs, mask_colors = \ get_default_screenshotting_data(args) diff --git a/scripts/scil_screenshot_volume_mosaic_overlap.py b/scripts/scil_screenshot_volume_mosaic_overlap.py index 4227a391b..c1ce9edf6 100755 --- a/scripts/scil_screenshot_volume_mosaic_overlap.py +++ b/scripts/scil_screenshot_volume_mosaic_overlap.py @@ -63,6 +63,7 @@ """ import argparse +import logging import numpy as np @@ -147,6 +148,7 @@ def _parse_args(parser): def main(): parser = _build_arg_parser() args = _parse_args(parser) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) vol_img, t_mask_img, labelmap_img, mask_imgs, mask_colors = \ get_default_screenshotting_data(args) diff --git a/scripts/scil_search_keywords.py b/scripts/scil_search_keywords.py index e9a23d9a6..d65d58f2e 100755 --- a/scripts/scil_search_keywords.py +++ b/scripts/scil_search_keywords.py @@ -50,10 +50,10 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - # Use INFO as default log level, switch to DEBUG if verbose - log_level = logging.DEBUG if args.verbose else logging.INFO - logging.getLogger().setLevel(log_level) + if args.verbose == "WARNING": + logging.getLogger().setLevel(logging.INFO) + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # Use directory of this script, should work with most installation setups script_dir = pathlib.Path(__file__).parent diff --git a/scripts/scil_sh_convert.py b/scripts/scil_sh_convert.py index 006805bff..ca4bbc56a 100755 --- a/scripts/scil_sh_convert.py +++ b/scripts/scil_sh_convert.py @@ -12,6 +12,7 @@ """ import argparse +import logging from dipy.data import get_sphere import nibabel as nib @@ -50,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_sh) assert_outputs_exist(parser, args, args.out_sh) diff --git a/scripts/scil_sh_fusion.py b/scripts/scil_sh_fusion.py index cc2860f9d..748835754 100755 --- a/scripts/scil_sh_fusion.py +++ b/scripts/scil_sh_fusion.py @@ -15,6 +15,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -56,6 +57,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_shs) assert_outputs_exist(parser, args, args.out_sh) diff --git a/scripts/scil_sh_to_aodf.py b/scripts/scil_sh_to_aodf.py old mode 100644 new mode 100755 index ffda00ddf..e383b3426 --- a/scripts/scil_sh_to_aodf.py +++ b/scripts/scil_sh_to_aodf.py @@ -97,9 +97,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if args.use_gpu and args.method == 'cosine': parser.error('Option --use_gpu is not supported for cosine filtering.') diff --git a/scripts/scil_sh_to_rish.py b/scripts/scil_sh_to_rish.py index 9463edae5..efccc079a 100755 --- a/scripts/scil_sh_to_rish.py +++ b/scripts/scil_sh_to_rish.py @@ -21,6 +21,7 @@ Formerly: scil_compute_rish_from_sh.py """ import argparse +import logging from dipy.reconst.shm import order_from_ncoef, sph_harm_ind_list import nibabel as nib @@ -54,6 +55,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_sh, optional=args.mask) diff --git a/scripts/scil_sh_to_sf.py b/scripts/scil_sh_to_sf.py index 52840fe0e..8a55dcd1e 100755 --- a/scripts/scil_sh_to_sf.py +++ b/scripts/scil_sh_to_sf.py @@ -14,6 +14,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -85,6 +86,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_sh, optional=[args.in_bvec, args.in_bval, args.in_b0]) diff --git a/scripts/scil_stats_group_comparison.py b/scripts/scil_stats_group_comparison.py index 79ca36900..919f0c8fc 100755 --- a/scripts/scil_stats_group_comparison.py +++ b/scripts/scil_stats_group_comparison.py @@ -117,9 +117,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) required_args = [args.in_json, args.in_participants] assert_inputs_exist(parser, required_args) @@ -154,8 +152,8 @@ def main(): curr_comparison_measure = ('_').join([b, m, v]) - logging.debug('______________________') - logging.debug('Measure to compare: {}'.format(curr_comparison_measure)) + logging.info('______________________') + logging.info('Measure to compare: {}'.format(curr_comparison_measure)) # Check normality of that metric across all groups @@ -164,19 +162,19 @@ def main(): groups_array = [] for group in my_group_dict: curr_sample = get_group_data_sample(my_group_dict, group, b, m, v) - logging.debug('Group {}'.format(group)) + logging.info('Group {}'.format(group)) current_normality[group] = verify_normality( curr_sample, alpha_error) if not current_normality[group][0]: overall_normality = False groups_array.append(curr_sample) - logging.debug('Normality result:') - logging.debug(current_normality) - logging.debug('Overall Normality:') - logging.debug(overall_normality) - logging.debug('Groups array:') - logging.debug(groups_array) + logging.info('Normality result:') + logging.info(current_normality) + logging.info('Overall Normality:') + logging.info(overall_normality) + logging.info('Groups array:') + logging.info(groups_array) # Generate graph of the metric if args.generate_graph: @@ -198,8 +196,8 @@ def main(): groups_array, normality=overall_normality, alpha=alpha_error) - logging.debug('Equality of variance result:') - logging.debug(variance_equality) + logging.info('Equality of variance result:') + logging.info(variance_equality) # Now we compare the groups population @@ -208,8 +206,8 @@ def main(): normality=overall_normality, homoscedasticity=variance_equality[1], alpha=alpha_error) - logging.debug('Main test result:') - logging.debug(difference) + logging.info('Main test result:') + logging.info(difference) # Finally if we have more than 2 groups and found a difference # We do a post hoc analysis to explore where is this difference @@ -234,8 +232,8 @@ def main(): else: diff_2_by_2 = [] - logging.debug('Summary of difference 2 by 2:') - logging.debug(diff_2_by_2) + logging.info('Summary of difference 2 by 2:') + logging.info(diff_2_by_2) # Write the current metric in the report curr_dict = write_current_dictionnary(curr_comparison_measure, diff --git a/scripts/scil_surface_apply_transform.py b/scripts/scil_surface_apply_transform.py index 405fbe6e4..a572dcb67 100755 --- a/scripts/scil_surface_apply_transform.py +++ b/scripts/scil_surface_apply_transform.py @@ -21,6 +21,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -65,6 +66,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_surface, args.ants_affine], args.ants_warp) diff --git a/scripts/scil_surface_convert.py b/scripts/scil_surface_convert.py index 181145b39..bebb0265a 100755 --- a/scripts/scil_surface_convert.py +++ b/scripts/scil_surface_convert.py @@ -13,6 +13,7 @@ Formerly: scil_convert_surface.py """ import argparse +import logging import os from trimeshpy.vtk_util import (load_polydata, @@ -59,9 +60,9 @@ def _build_arg_parser(): def main(): - parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_surface) assert_outputs_exist(parser, args, args.out_surface) diff --git a/scripts/scil_surface_flip.py b/scripts/scil_surface_flip.py index 6b742b6f0..9a8a2f440 100755 --- a/scripts/scil_surface_flip.py +++ b/scripts/scil_surface_flip.py @@ -11,6 +11,7 @@ """ import argparse +import logging from trimeshpy.io import load_mesh_from_file @@ -52,6 +53,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_surface) assert_outputs_exist(parser, args, args.out_surface) diff --git a/scripts/scil_surface_smooth.py b/scripts/scil_surface_smooth.py index ced6de253..de55814d6 100755 --- a/scripts/scil_surface_smooth.py +++ b/scripts/scil_surface_smooth.py @@ -62,6 +62,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_surface, args.vts_mask) assert_outputs_exist(parser, args, args.out_surface) @@ -73,9 +74,6 @@ def main(): if args.step_size <= 0.0: parser.error("Step size should be positive") - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - # Step size (zero for masked vertices) if args.vts_mask: mask = np.load(args.vts_mask) diff --git a/scripts/scil_tracking_local.py b/scripts/scil_tracking_local.py index e5f26da2e..7ab4b5631 100755 --- a/scripts/scil_tracking_local.py +++ b/scripts/scil_tracking_local.py @@ -134,9 +134,7 @@ def main(): t_init = perf_counter() parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if args.use_gpu: batch_size = args.batch_size or DEFAULT_BATCH_SIZE diff --git a/scripts/scil_tracking_local_dev.py b/scripts/scil_tracking_local_dev.py index f8c17fb24..37bcdcf86 100755 --- a/scripts/scil_tracking_local_dev.py +++ b/scripts/scil_tracking_local_dev.py @@ -50,8 +50,7 @@ import nibabel as nib import numpy as np -from dipy.io.stateful_tractogram import StatefulTractogram, Space, \ - set_sft_logger_level +from dipy.io.stateful_tractogram import StatefulTractogram, Space from dipy.io.stateful_tractogram import Origin from dipy.io.streamline import save_tractogram from nibabel.streamlines import detect_format, TrkFile @@ -156,9 +155,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if not nib.streamlines.is_supported(args.out_tractogram): parser.error('Invalid output streamline file format (must be trk or ' + @@ -194,7 +191,7 @@ def main(): our_origin = Origin('center') # Preparing everything - logging.debug("Loading seeding mask.") + logging.info("Loading seeding mask.") seed_img = nib.load(args.in_seed) seed_data = seed_img.get_fdata(caching='unchanged', dtype=float) if np.count_nonzero(seed_data) == 0: @@ -219,19 +216,19 @@ def main(): parser.error('Seed mask "{}" does not have any voxel with value > 0.' .format(args.in_seed)) - logging.debug("Loading tracking mask.") + logging.info("Loading tracking mask.") mask_img = nib.load(args.in_mask) mask_data = mask_img.get_fdata(caching='unchanged', dtype=float) mask_res = mask_img.header.get_zooms()[:3] mask = DataVolume(mask_data, mask_res, args.mask_interp) - logging.debug("Loading ODF SH data.") + logging.info("Loading ODF SH data.") odf_sh_img = nib.load(args.in_odf) odf_sh_data = odf_sh_img.get_fdata(caching='unchanged', dtype=float) odf_sh_res = odf_sh_img.header.get_zooms()[:3] dataset = DataVolume(odf_sh_data, odf_sh_res, args.sh_interp) - logging.debug("Instantiating propagator.") + logging.info("Instantiating propagator.") # Converting step size to vox space # We only support iso vox for now. assert odf_sh_res[0] == odf_sh_res[1] == odf_sh_res[2] @@ -245,7 +242,7 @@ def main(): args.sf_threshold, args.sf_threshold_init, theta, args.sphere, space=our_space, origin=our_origin) - logging.debug("Instantiating tracker.") + logging.info("Instantiating tracker.") tracker = Tracker(propagator, mask, seed_generator, nbr_seeds, min_nbr_pts, max_nbr_pts, args.max_invalid_nb_points, compression_th=args.compress, @@ -258,13 +255,13 @@ def main(): verbose=args.verbose) start = time.time() - logging.debug("Tracking...") + logging.info("Tracking...") streamlines, seeds = tracker.track() str_time = "%.2f" % (time.time() - start) - logging.debug("Tracked {} streamlines (out of {} seeds), in {} seconds.\n" - "Now saving..." - .format(len(streamlines), nbr_seeds, str_time)) + logging.info("Tracked {} streamlines (out of {} seeds), in {} seconds.\n" + "Now saving..." + .format(len(streamlines), nbr_seeds, str_time)) # save seeds if args.save_seeds is given # We seeded (and tracked) in vox, center, which is what is expected for @@ -274,10 +271,6 @@ def main(): else: data_per_streamline = {} - # Silencing SFT's logger if our logging is in DEBUG mode, because it - # typically produces a lot of outputs! - set_sft_logger_level('WARNING') - # Compared with scil_tracking_local, using sft rather than # LazyTractogram to deal with space. # Contrary to scilpy or dipy, where space after tracking is vox, here diff --git a/scripts/scil_tracking_pft.py b/scripts/scil_tracking_pft.py index dafae15cf..58e2dc098 100755 --- a/scripts/scil_tracking_pft.py +++ b/scripts/scil_tracking_pft.py @@ -138,8 +138,7 @@ def _build_arg_parser(): help='If set, save the seeds used for the tracking \n ' 'in the data_per_streamline property.') - log_g = p.add_argument_group('Logging options') - add_verbose_arg(log_g) + add_verbose_arg(p) return p @@ -147,9 +146,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_sh, args.in_seed, args.in_map_include, diff --git a/scripts/scil_tracking_pft_maps.py b/scripts/scil_tracking_pft_maps.py index 6d2dda730..c01206169 100755 --- a/scripts/scil_tracking_pft_maps.py +++ b/scripts/scil_tracking_pft_maps.py @@ -61,9 +61,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_wm, args.in_gm, args.in_csf]) assert_outputs_exist(parser, args, diff --git a/scripts/scil_tracking_pft_maps_edit.py b/scripts/scil_tracking_pft_maps_edit.py index f45ec1cb1..551e0c45d 100755 --- a/scripts/scil_tracking_pft_maps_edit.py +++ b/scripts/scil_tracking_pft_maps_edit.py @@ -8,6 +8,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -43,6 +44,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.map_include, args.map_exclude, args.additional_mask]) diff --git a/scripts/scil_tractogram_apply_transform.py b/scripts/scil_tractogram_apply_transform.py index d1e3da044..fd2b764a0 100755 --- a/scripts/scil_tractogram_apply_transform.py +++ b/scripts/scil_tractogram_apply_transform.py @@ -101,9 +101,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_moving_tractogram, args.in_target_file, @@ -130,8 +128,8 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written " - "(0 streamline).".format(args.out_tractogram)) + logging.info("The file {} won't be written " + "(0 streamline).".format(args.out_tractogram)) return if args.keep_invalid: diff --git a/scripts/scil_tractogram_apply_transform_to_hdf5.py b/scripts/scil_tractogram_apply_transform_to_hdf5.py index 382cac05e..9affc33c5 100755 --- a/scripts/scil_tractogram_apply_transform_to_hdf5.py +++ b/scripts/scil_tractogram_apply_transform_to_hdf5.py @@ -12,6 +12,7 @@ """ import argparse +import logging import os from dipy.io.stateful_tractogram import Space, Origin, StatefulTractogram @@ -67,6 +68,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_hdf5, args.in_target_file, args.in_transfo], args.in_deformation) diff --git a/scripts/scil_tractogram_assign_custom_color.py b/scripts/scil_tractogram_assign_custom_color.py index b49cbbb45..74b936f12 100755 --- a/scripts/scil_tractogram_assign_custom_color.py +++ b/scripts/scil_tractogram_assign_custom_color.py @@ -170,7 +170,7 @@ def save_colorbar(cmap, lbound, ubound, args): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.WARNING) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram, diff --git a/scripts/scil_tractogram_assign_uniform_color.py b/scripts/scil_tractogram_assign_uniform_color.py index edb37251c..420ee645e 100755 --- a/scripts/scil_tractogram_assign_uniform_color.py +++ b/scripts/scil_tractogram_assign_uniform_color.py @@ -65,7 +65,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.WARNING) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if len(args.in_tractograms) > 1 and args.out_tractogram: parser.error('Using multiple inputs, use --out_suffix.') diff --git a/scripts/scil_tractogram_commit.py b/scripts/scil_tractogram_commit.py index f03b9aba2..2f7e16792 100755 --- a/scripts/scil_tractogram_commit.py +++ b/scripts/scil_tractogram_commit.py @@ -218,7 +218,7 @@ def _save_results_wrapper(args, tmp_dir, ext, hdf5_file, offsets_list, new_hdf5_file.attrs['voxel_order'] = sft.voxel_order # Assign the weights into the hdf5, while respecting # the ordering of connections/streamlines - logging.debug('Adding commit weights to {}.'.format(new_filename)) + logging.info('Adding commit weights to {}.'.format(new_filename)) for i, key in enumerate(list(hdf5_file.keys())): new_group = new_hdf5_file.create_group(key) old_group = hdf5_file[key] @@ -275,9 +275,9 @@ def _save_results_wrapper(args, tmp_dir, ext, hdf5_file, offsets_list, essential_ind = np.where(streamline_weights > 0)[0] nonessential_ind = np.where(streamline_weights <= 0)[0] - logging.debug('{} essential streamlines were kept at'.format( + logging.info('{} essential streamlines were kept at'.format( len(essential_ind))) - logging.debug('{} nonessential streamlines were kept'.format( + logging.info('{} nonessential streamlines were kept'.format( len(nonessential_ind))) save_tractogram(sft[essential_ind], @@ -288,7 +288,7 @@ def _save_results_wrapper(args, tmp_dir, ext, hdf5_file, offsets_list, 'nonessential_tractogram.trk')) if args.keep_whole_tractogram: output_filename = os.path.join(out_dir, 'tractogram.trk') - logging.debug('Saving tractogram with weights as {}'.format( + logging.info('Saving tractogram with weights as {}'.format( output_filename)) save_tractogram(sft, output_filename) @@ -296,6 +296,15 @@ def _save_results_wrapper(args, tmp_dir, ext, hdf5_file, offsets_list, def main(): parser = _build_arg_parser() args = parser.parse_args() + # COMMIT has some c-level stdout and non-logging print that cannot + # be easily stopped. Manual redirection of all printed output + if args.verbose == "WARNING": + f = io.StringIO() + redirected_stdout = redirect_stdout(f) + redirect_stdout_c() + else: + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + redirected_stdout = redirect_stdout(sys.stdout) assert_inputs_exist(parser, [args.in_tractogram, args.in_dwi, args.in_bval, args.in_bvec], @@ -335,21 +344,11 @@ def main(): parser.error('{} does not have a compatible header with {}'.format( args.in_tractogram, args.in_dwi)) - # COMMIT has some c-level stdout and non-logging print that cannot - # be easily stopped. Manual redirection of all printed output - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - redirected_stdout = redirect_stdout(sys.stdout) - else: - f = io.StringIO() - redirected_stdout = redirect_stdout(f) - redirect_stdout_c() - tmp_dir = tempfile.TemporaryDirectory() hdf5_file = None offsets_list = None if ext == '.h5': - logging.debug('Reconstructing {} into a tractogram for COMMIT.'.format( + logging.info('Reconstructing {} into a tractogram for COMMIT.'.format( args.in_tractogram)) hdf5_file = h5py.File(args.in_tractogram, 'r') @@ -389,7 +388,7 @@ def main(): np.savetxt(tmp_bval_filename, shells_centroids[indices_shells], newline=' ', fmt='%i') fsl2mrtrix(tmp_bval_filename, args.in_bvec, tmp_scheme_filename) - logging.debug('Lauching COMMIT on {} shells at found at {}.'.format( + logging.info('Lauching COMMIT on {} shells at found at {}.'.format( len(shells_centroids), shells_centroids)) @@ -424,13 +423,13 @@ def main(): mit.set_model('StickZeppelinBall') if args.ball_stick: - logging.debug('Disabled zeppelin, using the Ball & Stick model.') + logging.info('Disabled zeppelin, using the Ball & Stick model.') para_diff = args.para_diff or 1.7E-3 perp_diff = [] isotropc_diff = args.iso_diff or [2.0E-3] mit.model.set(para_diff, perp_diff, isotropc_diff) else: - logging.debug('Using the Stick Zeppelin Ball model.') + logging.info('Using the Stick Zeppelin Ball model.') para_diff = args.para_diff or 1.7E-3 perp_diff = args.perp_diff or [0.85E-3, 0.51E-3] isotropc_diff = args.iso_diff or [1.7E-3, 3.0E-3] diff --git a/scripts/scil_tractogram_compress.py b/scripts/scil_tractogram_compress.py index ff7678bf7..5a3b75d5f 100755 --- a/scripts/scil_tractogram_compress.py +++ b/scripts/scil_tractogram_compress.py @@ -50,6 +50,7 @@ def compress_streamlines_wrapper(tractogram, error_rate): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) diff --git a/scripts/scil_tractogram_compute_TODI.py b/scripts/scil_tractogram_compute_TODI.py index 652d1b4fc..7612ffc1a 100755 --- a/scripts/scil_tractogram_compute_TODI.py +++ b/scripts/scil_tractogram_compute_TODI.py @@ -87,7 +87,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram, [args.mask, args.reference]) diff --git a/scripts/scil_tractogram_compute_density_map.py b/scripts/scil_tractogram_compute_density_map.py index c17486c41..b86effc95 100755 --- a/scripts/scil_tractogram_compute_density_map.py +++ b/scripts/scil_tractogram_compute_density_map.py @@ -11,6 +11,7 @@ Formerly: scil_compute_streamlines_density_map.py """ import argparse +import logging import numpy as np import nibabel as nib @@ -46,6 +47,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundle, optional=args.reference) assert_outputs_exist(parser, args, args.out_img) diff --git a/scripts/scil_tractogram_convert.py b/scripts/scil_tractogram_convert.py index 49bb5952a..2a161fae5 100755 --- a/scripts/scil_tractogram_convert.py +++ b/scripts/scil_tractogram_convert.py @@ -9,6 +9,7 @@ """ import argparse +import logging import os from dipy.io.streamline import save_tractogram @@ -42,6 +43,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram, args.reference) diff --git a/scripts/scil_tractogram_convert_hdf5_to_trk.py b/scripts/scil_tractogram_convert_hdf5_to_trk.py index 1a00a662c..39d962972 100755 --- a/scripts/scil_tractogram_convert_hdf5_to_trk.py +++ b/scripts/scil_tractogram_convert_hdf5_to_trk.py @@ -23,6 +23,7 @@ """ import argparse +import logging import os from dipy.io.stateful_tractogram import Space, Origin, StatefulTractogram @@ -75,6 +76,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_hdf5) assert_output_dirs_exist_and_empty(parser, args, args.out_dir, diff --git a/scripts/scil_tractogram_count_streamlines.py b/scripts/scil_tractogram_count_streamlines.py index 9918b891d..b42005993 100755 --- a/scripts/scil_tractogram_count_streamlines.py +++ b/scripts/scil_tractogram_count_streamlines.py @@ -10,6 +10,7 @@ import argparse import json +import logging import os from scilpy.io.utils import (add_json_args, @@ -38,6 +39,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) diff --git a/scripts/scil_tractogram_cut_streamlines.py b/scripts/scil_tractogram_cut_streamlines.py index ab2f8b20f..ad331eddb 100755 --- a/scripts/scil_tractogram_cut_streamlines.py +++ b/scripts/scil_tractogram_cut_streamlines.py @@ -70,9 +70,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_tractogram, args.in_mask]) assert_outputs_exist(parser, args, args.out_tractogram) diff --git a/scripts/scil_tractogram_detect_loops.py b/scripts/scil_tractogram_detect_loops.py index 279e81f3d..8dd119c4b 100755 --- a/scripts/scil_tractogram_detect_loops.py +++ b/scripts/scil_tractogram_detect_loops.py @@ -78,6 +78,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram, diff --git a/scripts/scil_tractogram_extract_ushape.py b/scripts/scil_tractogram_extract_ushape.py index f3156009a..54b7b9ee2 100755 --- a/scripts/scil_tractogram_extract_ushape.py +++ b/scripts/scil_tractogram_extract_ushape.py @@ -63,6 +63,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram, @@ -81,11 +82,11 @@ def main(): if len(ids_c) == 0: if args.no_empty: - logging.debug("The file {} won't be written " - "(0 streamline).".format(args.out_tractogram)) + logging.info("The file {} won't be written " + "(0 streamline).".format(args.out_tractogram)) return - logging.debug('The file {} contains 0 streamline.'.format( + logging.info('The file {} contains 0 streamline.'.format( args.out_tractogram)) save_tractogram(sft[ids_c], args.out_tractogram) @@ -100,8 +101,8 @@ def main(): if args.remaining_tractogram: if len(ids_l) == 0: if args.no_empty: - logging.debug("The file {} won't be written (0 streamline" - ").".format(args.remaining_tractogram)) + logging.info("The file {} won't be written (0 streamline" + ").".format(args.remaining_tractogram)) return logging.warning('No remaining streamlines.') diff --git a/scripts/scil_tractogram_filter_by_anatomy.py b/scripts/scil_tractogram_filter_by_anatomy.py index d2e2625f1..8bf23a1bb 100755 --- a/scripts/scil_tractogram_filter_by_anatomy.py +++ b/scripts/scil_tractogram_filter_by_anatomy.py @@ -212,16 +212,16 @@ def save_intermediate_sft(sft, outliers_sft, new_path, in_sft_name, if len(sft.streamlines) == 0: if no_empty: - logging.debug("The file" + sft_name + - " won't be written (0 streamlines)") + logging.info("The file" + sft_name + + " won't be written (0 streamlines)") save_tractogram(sft, sft_name) else: save_tractogram(sft, sft_name) if len(outliers_sft.streamlines): if no_empty: - logging.debug("The file" + outliers_sft_name + - " won't be written (0 streamlines)") + logging.info("The file" + outliers_sft_name + + " won't be written (0 streamlines)") save_tractogram(outliers_sft, outliers_sft_name) else: save_tractogram(outliers_sft, outliers_sft_name) @@ -248,8 +248,8 @@ def save_rejected(sft, new_sft, rejected_sft_name, no_empty): if len(rejected_sft.streamlines) == 0: if no_empty: - logging.debug("The file" + rejected_sft_name + - " won't be written (0 streamlines)") + logging.info("The file" + rejected_sft_name + + " won't be written (0 streamlines)") return save_tractogram(rejected_sft, rejected_sft_name) @@ -260,7 +260,7 @@ def display_count(o_dict, indent, sort_keys): Display the streamline count. """ o_dict_str = json.dumps(o_dict, indent=indent, sort_keys=sort_keys) - logging.debug("Streamline count:\n{}".format(o_dict_str)) + logging.info("Streamline count:\n{}".format(o_dict_str)) def save_count(o_dict, out_path, indent, sort_keys): @@ -275,6 +275,7 @@ def save_count(o_dict, out_path, indent, sort_keys): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_inputs_exist(parser, args.in_wmparc) @@ -284,9 +285,6 @@ def main(): nbr_cpu = validate_nbr_processes(parser, args) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - if args.angle <= 0: parser.error('Angle "{}" '.format(args.angle) + 'must be greater than or equal to 0') @@ -308,14 +306,14 @@ def main(): 'not compatible.') if args.minL == 0 and np.isinf(args.maxL): - logging.debug("You have not specified minL nor maxL. Output will " - "not be filtered according to length!") + logging.info("You have not specified minL nor maxL. Output will " + "not be filtered according to length!") if np.isinf(args.angle): - logging.debug("You have not specified the angle. Loops will " - "not be filtered!") + logging.info("You have not specified the angle. Loops will " + "not be filtered!") if args.ctx_dilation_radius == 0: - logging.debug("You have not specified the cortex dilation radius. " - "The wmparc atlas will not be dilated!") + logging.info("You have not specified the cortex dilation radius. " + "The wmparc atlas will not be dilated!") o_dict = {} step_dict = ['length', 'no_csf', 'end_in_atlas', 'no_loops'] @@ -353,9 +351,9 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written".format( - out_sft_name) + "(0 streamlines after " - + step + " filtering).") + logging.info("The file {} won't be written".format( + out_sft_name) + "(0 streamlines after " + + step + " filtering).") if args.verbose: display_count(o_dict, args.indent, args.sort_keys) @@ -365,8 +363,8 @@ def main(): save_tractogram(initial_sft, rejected_sft_name) return - logging.debug('The file {} contains 0 streamlines after '.format( - out_sft_name) + step + ' filtering') + logging.info('The file {} contains 0 streamlines after '.format( + out_sft_name) + step + ' filtering') save_tractogram(new_sft, out_sft_name) if args.save_rejected: @@ -414,9 +412,9 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written".format( - out_sft_name) + "(0 streamlines after " - + step + " filtering).") + logging.info("The file {} won't be written".format( + out_sft_name) + "(0 streamlines after " + + step + " filtering).") if args.verbose: display_count(o_dict, args.indent, args.sort_keys) @@ -426,8 +424,8 @@ def main(): save_tractogram(sft, rejected_sft_name) return - logging.debug('The file {} contains 0 streamlines after '.format( - out_sft_name) + step + ' filtering') + logging.info('The file {} contains 0 streamlines after '.format( + out_sft_name) + step + ' filtering') save_tractogram(new_sft, out_sft_name) if args.save_rejected: @@ -487,9 +485,9 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written".format( - out_sft_name) + "(0 streamlines after " - + step + " filtering).") + logging.info("The file {} won't be written".format( + out_sft_name) + "(0 streamlines after " + + step + " filtering).") if args.verbose: display_count(o_dict, args.indent, args.sort_keys) @@ -499,8 +497,8 @@ def main(): save_tractogram(sft, rejected_sft_name) return - logging.debug('The file {} contains 0 streamlines after '.format( - out_sft_name) + step + ' filtering') + logging.info('The file {} contains 0 streamlines after '.format( + out_sft_name) + step + ' filtering') save_tractogram(new_sft, out_sft_name) if args.save_rejected: @@ -539,9 +537,9 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written".format( - out_sft_name) + "(0 streamlines after " - + step + " filtering).") + logging.info("The file {} won't be written".format( + out_sft_name) + "(0 streamlines after " + + step + " filtering).") if args.verbose: display_count(o_dict, args.indent, args.sort_keys) @@ -551,11 +549,11 @@ def main(): save_tractogram(sft, rejected_sft_name) return - logging.debug('The file {} contains 0 streamlines after '.format( - out_sft_name) + step + ' filtering') + logging.info('The file {} contains 0 streamlines after '.format( + out_sft_name) + step + ' filtering') save_tractogram(new_sft, out_sft_name) - if args.verbose: + if args.verbose == "INFO" or args.verbose == "DEBUG": display_count(o_dict, args.indent, args.sort_keys) if args.save_counts: save_count(o_dict, args.out_path, args.indent, args.sort_keys) diff --git a/scripts/scil_tractogram_filter_by_length.py b/scripts/scil_tractogram_filter_by_length.py index 729026fa8..0ba977460 100755 --- a/scripts/scil_tractogram_filter_by_length.py +++ b/scripts/scil_tractogram_filter_by_length.py @@ -54,16 +54,14 @@ def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - if args.minL == 0 and np.isinf(args.maxL): - logging.debug("You have not specified minL nor maxL. Output will " - "simply be a copy of your input!") + logging.info("You have not specified minL nor maxL. Output will " + "simply be a copy of your input!") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) @@ -78,13 +76,13 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written " - "(0 streamline).".format(args.out_tractogram)) + logging.info("The file {} won't be written " + "(0 streamline).".format(args.out_tractogram)) return - logging.debug('The file {} contains 0 streamline'.format( - args.out_tractogram)) + logging.info('The file {} contains 0 streamline'.format( + args.out_tractogram)) save_tractogram(new_sft, args.out_tractogram) diff --git a/scripts/scil_tractogram_filter_by_orientation.py b/scripts/scil_tractogram_filter_by_orientation.py index e45066f65..0af0b652e 100755 --- a/scripts/scil_tractogram_filter_by_orientation.py +++ b/scripts/scil_tractogram_filter_by_orientation.py @@ -22,7 +22,6 @@ import json import logging -from dipy.io.stateful_tractogram import set_sft_logger_level from dipy.io.streamline import save_tractogram import numpy as np @@ -85,19 +84,13 @@ def _build_arg_parser(): def main(): - parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram, args.save_rejected) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - # Silencing SFT's logger if our logging is in DEBUG mode, because it - # typically produces a lot of outputs! - set_sft_logger_level('WARNING') - if args.min_x == 0 and np.isinf(args.max_x) and \ args.min_y == 0 and np.isinf(args.max_y) and \ args.min_z == 0 and np.isinf(args.max_z): @@ -120,11 +113,11 @@ def main(): if len(new_sft.streamlines) == 0: if args.no_empty: - logging.debug("The file {} won't be written " - "(0 streamline).".format(args.out_tractogram)) + logging.info("The file {} won't be written " + "(0 streamline).".format(args.out_tractogram)) else: - logging.debug('The file {} contains 0 streamline'.format( - args.out_tractogram)) + logging.info('The file {} contains 0 streamline'.format( + args.out_tractogram)) save_tractogram(new_sft, args.out_tractogram) diff --git a/scripts/scil_tractogram_filter_by_roi.py b/scripts/scil_tractogram_filter_by_roi.py index ff5ae9b2d..13dff1c2f 100755 --- a/scripts/scil_tractogram_filter_by_roi.py +++ b/scripts/scil_tractogram_filter_by_roi.py @@ -42,7 +42,6 @@ import os from copy import deepcopy -from dipy.io.stateful_tractogram import set_sft_logger_level from dipy.io.streamline import save_tractogram from dipy.io.utils import is_header_compatible import nibabel as nib @@ -237,24 +236,21 @@ def check_overwrite_distance(parser, args): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) overwrite_distance = check_overwrite_distance(parser, args) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram, args.save_rejected) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - set_sft_logger_level('WARNING') - if overwrite_distance: - logging.debug('Overwrite distance dictionnary {}'.format( + logging.info('Overwrite distance dictionnary {}'.format( overwrite_distance)) roi_opt_list, only_filtering_list = prepare_filtering_list(parser, args) o_dict = {} - logging.debug("Loading the tractogram...") + logging.info("Loading the tractogram...") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) if args.save_rejected: initial_sft = deepcopy(sft) @@ -266,7 +262,7 @@ def main(): total_kept_ids = np.arange(len(sft.streamlines)) for i, roi_opt in enumerate(roi_opt_list): - logging.debug("Preparing filtering from option: {}".format(roi_opt)) + logging.info("Preparing filtering from option: {}".format(roi_opt)) curr_dict = {} # Atlas needs an extra argument (value in the LUT) if roi_opt[0] == 'atlas_roi': @@ -365,8 +361,8 @@ def main(): else: raise ValueError("Unexpected filter type.") - logging.debug('The filtering options {0} resulted in ' - '{1} streamlines'.format(roi_opt, len(filtered_sft))) + logging.info('The filtering options {0} resulted in ' + '{1} streamlines'.format(roi_opt, len(filtered_sft))) sft = filtered_sft @@ -385,12 +381,12 @@ def main(): if not filtered_sft: if args.no_empty: - logging.debug("The file {} won't be written (0 streamline)".format( + logging.info("The file {} won't be written (0 streamline)".format( args.out_tractogram)) return - logging.debug('The file {} contains 0 streamline'.format( + logging.info('The file {} contains 0 streamline'.format( args.out_tractogram)) save_tractogram(sft, args.out_tractogram) @@ -400,8 +396,8 @@ def main(): total_kept_ids) if len(rejected_ids) == 0 and args.no_empty: - logging.debug("Rejected streamlines file won't be written (0 " - "streamline).") + logging.info("Rejected streamlines file won't be written (0 " + "streamline).") return sft = initial_sft[rejected_ids] diff --git a/scripts/scil_tractogram_fix_trk.py b/scripts/scil_tractogram_fix_trk.py index 404159f68..c4302bbb8 100755 --- a/scripts/scil_tractogram_fix_trk.py +++ b/scripts/scil_tractogram_fix_trk.py @@ -126,8 +126,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - logging.getLogger().setLevel(logging.INFO) warning_msg = """ # This script is still experimental, DSI-Studio and Startrack # evolve quickly and results may vary depending on the data itself diff --git a/scripts/scil_tractogram_flip.py b/scripts/scil_tractogram_flip.py index 18ed8c739..a95449fd1 100755 --- a/scripts/scil_tractogram_flip.py +++ b/scripts/scil_tractogram_flip.py @@ -12,6 +12,7 @@ """ import argparse +import logging from dipy.io.streamline import save_tractogram @@ -48,6 +49,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) diff --git a/scripts/scil_tractogram_math.py b/scripts/scil_tractogram_math.py index 508ad1121..0e5bf0f65 100755 --- a/scripts/scil_tractogram_math.py +++ b/scripts/scil_tractogram_math.py @@ -111,9 +111,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractograms) assert_outputs_exist(parser, args, args.out_tractogram, diff --git a/scripts/scil_tractogram_print_info.py b/scripts/scil_tractogram_print_info.py index 84cbcac74..48e604d29 100755 --- a/scripts/scil_tractogram_print_info.py +++ b/scripts/scil_tractogram_print_info.py @@ -17,6 +17,7 @@ import argparse import json +import logging from dipy.tracking.streamlinespeed import length import numpy as np @@ -24,6 +25,7 @@ from scilpy.io.streamlines import load_tractogram_with_reference from scilpy.io.utils import (add_json_args, add_reference_arg, + add_verbose_arg, assert_inputs_exist) @@ -34,6 +36,7 @@ def _build_arg_parser(): p.add_argument('in_tractogram', help='Tractogram file.') add_reference_arg(p) + add_verbose_arg(p) add_json_args(p) return p @@ -42,6 +45,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) diff --git a/scripts/scil_tractogram_qbx.py b/scripts/scil_tractogram_qbx.py index d27f9873c..fdc8da4a0 100755 --- a/scripts/scil_tractogram_qbx.py +++ b/scripts/scil_tractogram_qbx.py @@ -56,6 +56,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, [], optional=args.out_centroids) @@ -63,9 +64,6 @@ def main(): args.out_clusters_dir, create_dir=True) - if args.verbose: - logging.getLogger().setLevel(logging.INFO) - sft = load_tractogram_with_reference(parser, args, args.in_tractogram) streamlines = sft.streamlines thresholds = [40, 30, 20, args.dist_thresh] diff --git a/scripts/scil_tractogram_register.py b/scripts/scil_tractogram_register.py index f3ce41806..c8e467a54 100755 --- a/scripts/scil_tractogram_register.py +++ b/scripts/scil_tractogram_register.py @@ -13,6 +13,7 @@ """ import argparse +import logging import os from dipy.align.streamlinear import whole_brain_slr @@ -59,6 +60,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.moving_tractogram, args.static_tractogram]) diff --git a/scripts/scil_tractogram_remove_invalid.py b/scripts/scil_tractogram_remove_invalid.py index 8bb3a4bae..db72c2f70 100755 --- a/scripts/scil_tractogram_remove_invalid.py +++ b/scripts/scil_tractogram_remove_invalid.py @@ -61,6 +61,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # Equivalent of add_bbox_arg(p): always ignoring invalid streamlines for # this script. diff --git a/scripts/scil_tractogram_resample.py b/scripts/scil_tractogram_resample.py index 606b6cfdf..f00c01eb2 100755 --- a/scripts/scil_tractogram_resample.py +++ b/scripts/scil_tractogram_resample.py @@ -32,7 +32,6 @@ import argparse import logging -from dipy.io.stateful_tractogram import set_sft_logger_level from dipy.io.streamline import save_tractogram from scilpy.io.streamlines import load_tractogram_with_reference @@ -116,6 +115,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if (args.point_wise_std is not None and args.point_wise_std <= 0) or \ (args.streamline_wise_std is not None and @@ -125,21 +125,15 @@ def main(): assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) - log_level = logging.WARNING - if args.verbose: - log_level = logging.DEBUG - set_sft_logger_level('INFO') - logging.getLogger().setLevel(log_level) - - logging.debug("Loading sft.") + logging.info("Loading sft.") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) original_number = len(sft.streamlines) if args.never_upsample and args.nb_streamlines > original_number: args.nb_streamlines = original_number - logging.debug("Done. Now getting {} streamlines." - .format(args.nb_streamlines)) + logging.info("Done. Now getting {} streamlines." + .format(args.nb_streamlines)) if args.nb_streamlines > original_number: # Check is done here because it is not required if downsampling @@ -155,8 +149,8 @@ def main(): # output contains rejected streamlines, we don't use them. sft, _ = split_sft_randomly_per_cluster( sft, [args.nb_streamlines], args.seed, args.qbx_thresholds) - logging.debug("Kept {} out of {} expected streamlines." - .format(len(sft), args.nb_streamlines)) + logging.info("Kept {} out of {} expected streamlines." + .format(len(sft), args.nb_streamlines)) else: # output is a list of two: kept and rejected. sft = split_sft_randomly(sft, args.nb_streamlines, args.seed)[0] diff --git a/scripts/scil_tractogram_resample_nb_points.py b/scripts/scil_tractogram_resample_nb_points.py index 0bf7d7f89..5cd3e5a2b 100755 --- a/scripts/scil_tractogram_resample_nb_points.py +++ b/scripts/scil_tractogram_resample_nb_points.py @@ -8,6 +8,7 @@ Formerly: scil_resample_streamlines.py """ import argparse +import logging from dipy.io.streamline import save_tractogram @@ -44,9 +45,9 @@ def _build_arg_parser(): def main(): - parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) diff --git a/scripts/scil_tractogram_seed_density_map.py b/scripts/scil_tractogram_seed_density_map.py index 6ae570088..25954f5ec 100755 --- a/scripts/scil_tractogram_seed_density_map.py +++ b/scripts/scil_tractogram_seed_density_map.py @@ -8,6 +8,7 @@ """ import argparse +import logging from dipy.io.streamline import load_tractogram from nibabel import Nifti1Image @@ -49,6 +50,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.tractogram_filename]) assert_outputs_exist(parser, args, [args.seed_density_filename]) diff --git a/scripts/scil_tractogram_segment_bundles.py b/scripts/scil_tractogram_segment_bundles.py index e23731834..790004cd9 100755 --- a/scripts/scil_tractogram_segment_bundles.py +++ b/scripts/scil_tractogram_segment_bundles.py @@ -73,9 +73,6 @@ def _build_arg_parser(): p.add_argument('--out_dir', default='voting_results', help='Path for the output directory [%(default)s].') - p.add_argument('--log_level', default='INFO', - choices=['DEBUG', 'INFO', 'WARNING', 'ERROR'], - help='Log level of the logging class.') p.add_argument('--minimal_vote_ratio', type=float, default=0.5, help='Streamlines will only be considered for saving if\n' 'recognized often enough [%(default)s].') @@ -96,6 +93,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + args.in_models_directories = [os.path.join(args.in_directory, x) for x in os.listdir(args.in_directory) if os.path.isdir(os.path.join( @@ -122,9 +121,8 @@ def main(): fmt='%(asctime)s, %(name)s %(levelname)s %(message)s', datefmt='%H:%M:%S') file_handler.setFormatter(formatter) - logging.getLogger().setLevel(args.log_level) logging.getLogger().addHandler(file_handler) - coloredlogs.install(level=args.log_level) + coloredlogs.install(level=logging.getLevelName(args.verbose)) transfo = load_matrix_in_any_format(args.in_transfo) if args.inverse: diff --git a/scripts/scil_tractogram_segment_bundles_for_connectivity.py b/scripts/scil_tractogram_segment_bundles_for_connectivity.py index fab98419f..11755ad78 100755 --- a/scripts/scil_tractogram_segment_bundles_for_connectivity.py +++ b/scripts/scil_tractogram_segment_bundles_for_connectivity.py @@ -37,8 +37,7 @@ import time import coloredlogs -from dipy.io.stateful_tractogram import (StatefulTractogram, - set_sft_logger_level) +from dipy.io.stateful_tractogram import StatefulTractogram from dipy.io.streamline import save_tractogram from dipy.io.utils import get_reference_info, is_header_compatible from dipy.tracking.streamlinespeed import length @@ -243,6 +242,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + coloredlogs.install(level=logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractograms+[args.in_labels], args.reference) @@ -264,11 +265,6 @@ def main(): assert_output_dirs_exist_and_empty(parser, args, args.out_dir, create_dir=True) - log_level = logging.INFO if args.verbose else logging.WARNING - logging.getLogger().setLevel(log_level) - coloredlogs.install(level=log_level) - set_sft_logger_level('WARNING') - # Load everything img_labels = nib.load(args.in_labels) data_labels = get_data_as_labels(img_labels) diff --git a/scripts/scil_tractogram_segment_one_bundles.py b/scripts/scil_tractogram_segment_one_bundles.py index 7fd4274ef..ae41b209f 100755 --- a/scripts/scil_tractogram_segment_one_bundles.py +++ b/scripts/scil_tractogram_segment_one_bundles.py @@ -96,13 +96,11 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_tractogram, args.in_transfo]) assert_outputs_exist(parser, args, args.out_tractogram) - log_level = logging.INFO if args.verbose else logging.WARNING - logging.getLogger().setLevel(log_level) - wb_file = load_tractogram_with_reference(parser, args, args.in_tractogram) wb_streamlines = wb_file.streamlines model_file = load_tractogram_with_reference(parser, args, args.in_model) diff --git a/scripts/scil_tractogram_shuffle.py b/scripts/scil_tractogram_shuffle.py index 9cf4430f8..48a8fbbb6 100755 --- a/scripts/scil_tractogram_shuffle.py +++ b/scripts/scil_tractogram_shuffle.py @@ -8,6 +8,7 @@ """ import argparse +import logging from dipy.io.streamline import save_tractogram @@ -40,6 +41,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) diff --git a/scripts/scil_tractogram_smooth.py b/scripts/scil_tractogram_smooth.py index c25e44fdb..6cc5f07c5 100755 --- a/scripts/scil_tractogram_smooth.py +++ b/scripts/scil_tractogram_smooth.py @@ -77,13 +77,11 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) assert_outputs_exist(parser, args, args.out_tractogram) - log_level = logging.INFO if args.verbose else logging.WARNING - logging.getLogger().setLevel(log_level) - sft = load_tractogram_with_reference(parser, args, args.in_tractogram) smoothed_streamlines = [] for streamline in sft.streamlines: diff --git a/scripts/scil_tractogram_split.py b/scripts/scil_tractogram_split.py index 5220acf9d..5a2a4049e 100755 --- a/scripts/scil_tractogram_split.py +++ b/scripts/scil_tractogram_split.py @@ -17,7 +17,6 @@ import logging import os -from dipy.io.stateful_tractogram import set_sft_logger_level from dipy.io.streamline import save_tractogram import numpy as np @@ -83,6 +82,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_tractogram) _, out_extension = os.path.splitext(args.in_tractogram) @@ -93,13 +93,7 @@ def main(): assert_outputs_exist(parser, args, os.path.join( args.out_dir, '{}_0{}'.format(args.out_prefix, out_extension))) - log_level = logging.WARNING - if args.verbose: - log_level = logging.DEBUG - set_sft_logger_level('INFO') - logging.getLogger().setLevel(log_level) - - logging.debug("Loading sft.") + logging.info("Loading sft.") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) streamlines_count = len(sft.streamlines) diff --git a/scripts/scil_tractogram_uniformize_endpoints.py b/scripts/scil_tractogram_uniformize_endpoints.py index 3a165a374..f8617b594 100755 --- a/scripts/scil_tractogram_uniformize_endpoints.py +++ b/scripts/scil_tractogram_uniformize_endpoints.py @@ -75,9 +75,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_bundle) assert_outputs_exist(parser, args, args.out_bundle) diff --git a/scripts/scil_visualize_bingham_fit.py b/scripts/scil_visualize_bingham_fit.py index ad4e16cad..8e272029a 100755 --- a/scripts/scil_visualize_bingham_fit.py +++ b/scripts/scil_visualize_bingham_fit.py @@ -9,6 +9,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -128,6 +129,7 @@ def main(): args = _parse_args(parser) data = _get_data_from_inputs(args) sph = get_sphere(args.sphere) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) actors = create_bingham_slicer(data, args.axis_name, args.slice_index, sph, diff --git a/scripts/scil_visualize_bundles.py b/scripts/scil_visualize_bundles.py index a29d5b130..17422871a 100755 --- a/scripts/scil_visualize_bundles.py +++ b/scripts/scil_visualize_bundles.py @@ -24,6 +24,7 @@ import glob import json import itertools +import logging import nibabel as nib import numpy as np import os @@ -106,6 +107,7 @@ def random_rgb(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # Make sure the colors are consistent between executions if args.random_coloring is not None: diff --git a/scripts/scil_visualize_bundles_mosaic.py b/scripts/scil_visualize_bundles_mosaic.py index 5f323c2a0..bfec4b29b 100755 --- a/scripts/scil_visualize_bundles_mosaic.py +++ b/scripts/scil_visualize_bundles_mosaic.py @@ -174,6 +174,7 @@ def random_rgb(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_volume) assert_outputs_exist(parser, args, args.out_image) diff --git a/scripts/scil_visualize_connectivity.py b/scripts/scil_visualize_connectivity.py index e16dd06c5..c9b3168b5 100755 --- a/scripts/scil_visualize_connectivity.py +++ b/scripts/scil_visualize_connectivity.py @@ -148,6 +148,7 @@ def write_values(ax, matrix, properties): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_matrix) if not args.show_only: diff --git a/scripts/scil_visualize_fodf.py b/scripts/scil_visualize_fodf.py index dfaf24f00..422b0c843 100755 --- a/scripts/scil_visualize_fodf.py +++ b/scripts/scil_visualize_fodf.py @@ -14,6 +14,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -265,6 +266,7 @@ def main(): data = _get_data_from_inputs(args) sph = get_sphere(args.sphere) sh_order, full_basis = get_sh_order_and_fullness(data['fodf'].shape[-1]) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) actors = [] diff --git a/scripts/scil_visualize_gradients.py b/scripts/scil_visualize_gradients.py index 5062b735b..155c5efd0 100755 --- a/scripts/scil_visualize_gradients.py +++ b/scripts/scil_visualize_gradients.py @@ -80,8 +80,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # -- Perform checks if args.in_gradient_scheme is not None: diff --git a/scripts/scil_visualize_histogram.py b/scripts/scil_visualize_histogram.py index e8b79cc86..ca490fb56 100755 --- a/scripts/scil_visualize_histogram.py +++ b/scripts/scil_visualize_histogram.py @@ -12,6 +12,7 @@ """ import argparse +import logging import matplotlib.pyplot as plt import nibabel as nib @@ -57,6 +58,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_metric, args.in_mask]) assert_outputs_exist(parser, args, args.out_png) diff --git a/scripts/scil_visualize_scatterplot.py b/scripts/scil_visualize_scatterplot.py index 566d3d54c..04977c401 100755 --- a/scripts/scil_visualize_scatterplot.py +++ b/scripts/scil_visualize_scatterplot.py @@ -41,6 +41,7 @@ import argparse import copy import json +import logging import os import matplotlib.pyplot as plt @@ -140,6 +141,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_x_map, args.in_y_map]) diff --git a/scripts/scil_visualize_seeds.py b/scripts/scil_visualize_seeds.py index 76daf1a7e..de1822283 100755 --- a/scripts/scil_visualize_seeds.py +++ b/scripts/scil_visualize_seeds.py @@ -11,6 +11,7 @@ """ import argparse +import logging from dipy.io.streamline import load_tractogram from fury import window, actor @@ -41,6 +42,8 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + assert_inputs_exist(parser, [args.tractogram]) assert_outputs_exist(parser, args, [], [args.save]) diff --git a/scripts/scil_visualize_seeds_3d.py b/scripts/scil_visualize_seeds_3d.py index fe527559d..3cc0c66f1 100755 --- a/scripts/scil_visualize_seeds_3d.py +++ b/scripts/scil_visualize_seeds_3d.py @@ -11,6 +11,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -65,6 +66,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_seed_map, [args.tractogram]) diff --git a/scripts/scil_volume_apply_transform.py b/scripts/scil_volume_apply_transform.py index a940ac2ea..6a6085dde 100755 --- a/scripts/scil_volume_apply_transform.py +++ b/scripts/scil_volume_apply_transform.py @@ -11,7 +11,7 @@ """ import argparse -import warnings +import logging import nibabel as nib import numpy as np @@ -51,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_file, args.in_target_file, args.in_transfo]) @@ -72,9 +73,9 @@ def main(): moving = nib.load(args.in_file) if moving.get_fdata().ndim == 4: - warnings.warn('You are applying a transform to a 4D dwi volume, ' - 'make sure to rotate your bvecs with ' - 'scil_gradients_apply_transform.py') + logging.warning('You are applying a transform to a 4D dwi volume, ' + 'make sure to rotate your bvecs with ' + 'scil_gradients_apply_transform.py') reference = nib.load(args.in_target_file) diff --git a/scripts/scil_volume_count_non_zero_voxels.py b/scripts/scil_volume_count_non_zero_voxels.py index a72eafc26..c8b984df9 100755 --- a/scripts/scil_volume_count_non_zero_voxels.py +++ b/scripts/scil_volume_count_non_zero_voxels.py @@ -13,6 +13,7 @@ """ import argparse +import logging import os import nibabel as nib @@ -51,6 +52,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) # out_filename can exist or not diff --git a/scripts/scil_volume_crop.py b/scripts/scil_volume_crop.py index e9a61a204..71230fc08 100755 --- a/scripts/scil_volume_crop.py +++ b/scripts/scil_volume_crop.py @@ -14,6 +14,7 @@ """ import argparse +import logging import pickle import nibabel as nib @@ -56,6 +57,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image, args.input_bbox) assert_outputs_exist(parser, args, args.out_image, args.output_bbox) diff --git a/scripts/scil_volume_flip.py b/scripts/scil_volume_flip.py index 524f95897..a9bc9897d 100755 --- a/scripts/scil_volume_flip.py +++ b/scripts/scil_volume_flip.py @@ -7,6 +7,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -36,9 +37,9 @@ def _build_arg_parser(): def main(): - parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) assert_outputs_exist(parser, args, args.out_image) diff --git a/scripts/scil_volume_math.py b/scripts/scil_volume_math.py index 601d8a0bf..ae99a93a2 100755 --- a/scripts/scil_volume_math.py +++ b/scripts/scil_volume_math.py @@ -65,9 +65,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - - if args.verbose: - logging.getLogger().setLevel(logging.INFO) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_outputs_exist(parser, args, args.out_image) diff --git a/scripts/scil_volume_remove_outliers_ransac.py b/scripts/scil_volume_remove_outliers_ransac.py index 8635b3568..7c3377b35 100755 --- a/scripts/scil_volume_remove_outliers_ransac.py +++ b/scripts/scil_volume_remove_outliers_ransac.py @@ -51,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) assert_outputs_exist(parser, args, args.out_image) @@ -65,9 +66,6 @@ def main(): parser.error('--fit_thr should be greater than 0. Current value: {}' .format(args.fit_thr)) - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - in_img = nib.load(args.in_image) in_data = in_img.get_fdata(dtype=np.float32) diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index a037206dc..42ed7d88a 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -61,6 +61,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) # Checking args assert_inputs_exist(parser, args.in_image, args.ref) @@ -68,9 +69,6 @@ def main(): if args.enforce_dimensions and not args.ref: parser.error("Cannot enforce dimensions without a reference image") - if args.verbose: - logging.getLogger().setLevel(logging.DEBUG) - if args.volume_size and (not len(args.volume_size) == 1 and not len(args.volume_size) == 3): parser.error('Invalid dimensions for --volume_size.') @@ -79,7 +77,7 @@ def main(): not len(args.voxel_size) == 3): parser.error('Invalid dimensions for --voxel_size.') - logging.debug('Loading raw data from %s', args.in_image) + logging.info('Loading raw data from %s', args.in_image) img = nib.load(args.in_image) @@ -90,7 +88,7 @@ def main(): enforce_dimensions=args.enforce_dimensions) # Saving results - logging.debug('Saving resampled data to %s', args.out_image) + logging.info('Saving resampled data to %s', args.out_image) nib.save(resampled_img, args.out_image) diff --git a/scripts/scil_volume_reshape_to_reference.py b/scripts/scil_volume_reshape_to_reference.py index ffeddc10a..bb96a0775 100755 --- a/scripts/scil_volume_reshape_to_reference.py +++ b/scripts/scil_volume_reshape_to_reference.py @@ -13,6 +13,7 @@ """ import argparse +import logging import nibabel as nib import numpy as np @@ -50,6 +51,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, [args.in_file, args.in_ref_file]) assert_outputs_exist(parser, args, args.out_file) diff --git a/scripts/scil_volume_stats_in_ROI.py b/scripts/scil_volume_stats_in_ROI.py index c57a0b520..e0009be4c 100755 --- a/scripts/scil_volume_stats_in_ROI.py +++ b/scripts/scil_volume_stats_in_ROI.py @@ -67,6 +67,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) if args.metrics_dir and os.path.exists(args.metrics_dir): list_metrics_files = glob.glob(os.path.join(args.metrics_dir, diff --git a/scripts/scil_volume_stats_in_labels.py b/scripts/scil_volume_stats_in_labels.py index 6471597e3..c4b926afb 100755 --- a/scripts/scil_volume_stats_in_labels.py +++ b/scripts/scil_volume_stats_in_labels.py @@ -11,6 +11,7 @@ import argparse import json +import logging import nibabel as nib import numpy as np @@ -42,6 +43,7 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) required = args.in_labels, args.in_seed_maps, args.in_labels_lut assert_inputs_exist(parser, required) diff --git a/scripts/tests/test_dwi_detect_volume_outliers.py b/scripts/tests/test_dwi_detect_volume_outliers.py index d2782e865..4f86d7761 100644 --- a/scripts/tests/test_dwi_detect_volume_outliers.py +++ b/scripts/tests/test_dwi_detect_volume_outliers.py @@ -25,5 +25,5 @@ def test_execution(script_runner): in_bvec = os.path.join(get_home(), 'processing', 'dwi.bvec') ret = script_runner.run('scil_dwi_detect_volume_outliers.py', in_dwi, - in_bval, in_bvec) + in_bval, in_bvec, '-v') assert ret.success diff --git a/scripts/tests/test_mti_adjust_B1_header.py b/scripts/tests/test_mti_adjust_B1_header.py new file mode 100644 index 000000000..cf432d190 --- /dev/null +++ b/scripts/tests/test_mti_adjust_B1_header.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import tempfile + +from scilpy.io.fetcher import fetch_data, get_home, get_testing_files_dict + +# If they already exist, this only takes 5 seconds (check md5sum) +fetch_data(get_testing_files_dict(), keys=['ihMT.zip']) +tmp_dir = tempfile.TemporaryDirectory() + + +def test_help_option(script_runner): + ret = script_runner.run('scil_mti_adjust_B1_header.py', '--help') + assert ret.success + + +def test_execution_ihMT_no_option(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_b1_map = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + + # no option + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + tmp_dir.name, in_b1_json, '-f') + assert ret.success diff --git a/scripts/tests/test_mti_maps_MT.py b/scripts/tests/test_mti_maps_MT.py index 2658984c5..f8b77eb76 100644 --- a/scripts/tests/test_mti_maps_MT.py +++ b/scripts/tests/test_mti_maps_MT.py @@ -21,6 +21,11 @@ def test_execution_MT_no_option(script_runner): in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + in_e1_mtoff = os.path.join(get_home(), 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') in_e2_mtoff = os.path.join(get_home(), @@ -56,13 +61,14 @@ def test_execution_MT_no_option(script_runner): # no option ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, - in_mask, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, - '--in_mton', in_e1_mton, in_e2_mton, in_e3_mton, - in_e4_mton, in_e5_mton, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, in_e4_t1w, in_e5_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, '-f') assert ret.success @@ -72,6 +78,11 @@ def test_execution_MT_prefix(script_runner): in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + in_e1_mtoff = os.path.join(get_home(), 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') in_e2_mtoff = os.path.join(get_home(), @@ -107,23 +118,145 @@ def test_execution_MT_prefix(script_runner): # --out_prefix ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, - in_mask, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, - '--in_mton', in_e1_mton, in_e2_mton, in_e3_mton, - in_e4_mton, in_e5_mton, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, in_e4_t1w, in_e5_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, '--out_prefix', 'sub_01', '-f') assert ret.success +def test_execution_MT_extended(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + in_e2_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mtoff_mtsat.nii.gz') + in_e3_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mtoff_mtsat.nii.gz') + in_e4_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mtoff_mtsat.nii.gz') + in_e5_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + in_e2_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mton_mtsat.nii.gz') + in_e3_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mton_mtsat.nii.gz') + in_e4_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mton_mtsat.nii.gz') + in_e5_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + in_e2_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-t1w_mtsat.nii.gz') + in_e3_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-t1w_mtsat.nii.gz') + in_e4_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-t1w_mtsat.nii.gz') + in_e5_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-t1w_mtsat.nii.gz') + + # --extended + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, + in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, + in_e4_t1w, in_e5_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, + '--extended', + '-f') + assert ret.success + + +def test_execution_MT_filtering(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + in_e2_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mtoff_mtsat.nii.gz') + in_e3_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mtoff_mtsat.nii.gz') + in_e4_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mtoff_mtsat.nii.gz') + in_e5_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + in_e2_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mton_mtsat.nii.gz') + in_e3_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mton_mtsat.nii.gz') + in_e4_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mton_mtsat.nii.gz') + in_e5_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + in_e2_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-t1w_mtsat.nii.gz') + in_e3_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-t1w_mtsat.nii.gz') + in_e4_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-t1w_mtsat.nii.gz') + in_e5_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-t1w_mtsat.nii.gz') + + # --filtering + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, + in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, + in_e4_t1w, in_e5_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, + '--filtering', + '-f') + assert ret.success + + def test_execution_MT_B1_map(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + in_e1_mtoff = os.path.join(get_home(), 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') in_e2_mtoff = os.path.join(get_home(), @@ -156,19 +289,234 @@ def test_execution_MT_B1_map(script_runner): 'MT', 'sub-001_echo-4_acq-t1w_mtsat.nii.gz') in_e5_t1w = os.path.join(get_home(), 'MT', 'sub-001_echo-5_acq-t1w_mtsat.nii.gz') + in_b1_map = os.path.join(get_home(), 'MT', 'sub-001_run-01_B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') # --in_B1_map ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, - in_mask, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, - '--in_mton', in_e1_mton, in_e2_mton, in_e3_mton, - in_e4_mton, in_e5_mton, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_negative', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, in_e4_t1w, in_e5_t1w, - '--in_B1_map', in_b1_map, + '--in_jsons', in_mtoff_json, in_t1w_json, + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'empiric', '--out_prefix', 'sub-01', '-f') assert ret.success + + +def test_execution_MT_wrong_echoes(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + in_e2_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mtoff_mtsat.nii.gz') + in_e3_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mtoff_mtsat.nii.gz') + in_e4_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mtoff_mtsat.nii.gz') + in_e5_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + in_e2_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-mton_mtsat.nii.gz') + in_e3_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-mton_mtsat.nii.gz') + in_e4_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-mton_mtsat.nii.gz') + in_e5_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + in_e2_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-2_acq-t1w_mtsat.nii.gz') + in_e3_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-3_acq-t1w_mtsat.nii.gz') + in_e4_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-4_acq-t1w_mtsat.nii.gz') + in_e5_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-5_acq-t1w_mtsat.nii.gz') + + # Wrong number of echoes for negative + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, in_e2_mtoff, + in_e3_mtoff, in_e4_mtoff, in_e5_mtoff, + '--in_positive', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, in_e5_mton, + '--in_negative', in_e1_mton, in_e2_mton, + in_e3_mton, in_e4_mton, + '--in_mtoff_t1', in_e1_t1w, in_e2_t1w, in_e3_t1w, + in_e4_t1w, in_e5_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, + '-f') + assert (not ret.success) + + +def test_execution_MT_single_echoe(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + + # Single echoe + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, + '--in_positive', in_e1_mton, + '--in_negative', in_e1_mton, + '--in_mtoff_t1', in_e1_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, + '-f') + assert ret.success + + +def test_execution_MT_B1_not_T1(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + + in_b1_map = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') + + # B1 no T1 should raise warning. + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, + '--in_positive', in_e1_mton, + '--in_negative', in_e1_mton, + '--in_jsons', in_mtoff_json, in_t1w_json, + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'empiric', + '-f') + assert ret.success + + +def test_execution_MT_B1_no_fit(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + + in_b1_map = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') + + # B1 model_based but no fit values + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, + '--in_positive', in_e1_mton, + '--in_negative', in_e1_mton, + '--in_mtoff_t1', in_e1_t1w, + '--in_jsons', in_mtoff_json, in_t1w_json, + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'model_based', + '-f') + assert (not ret.success) + + +def test_execution_MT_acq_params(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'MT', 'mask.nii.gz') + + in_mtoff_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.json') + in_t1w_json = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.json') + + in_e1_mtoff = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mtoff_mtsat.nii.gz') + + in_e1_mton = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-mton_mtsat.nii.gz') + + in_e1_t1w = os.path.join(get_home(), + 'MT', 'sub-001_echo-1_acq-t1w_mtsat.nii.gz') + + # Acquisition parameters + ret = script_runner.run('scil_mti_maps_MT.py', tmp_dir.name, + '--mask', in_mask, + '--in_mtoff_pd', in_e1_mtoff, + '--in_positive', in_e1_mton, + '--in_negative', in_e1_mton, + '--in_mtoff_t1', in_e1_t1w, + '--in_acq_parameters', "15", "15", "0.1", "0.1", + '-f') + assert ret.success diff --git a/scripts/tests/test_mti_maps_ihMT.py b/scripts/tests/test_mti_maps_ihMT.py index 529d0daeb..960a39bbb 100644 --- a/scripts/tests/test_mti_maps_ihMT.py +++ b/scripts/tests/test_mti_maps_ihMT.py @@ -21,6 +21,11 @@ def test_execution_ihMT_no_option(script_runner): in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + in_e1_altnp = os.path.join(get_home(), 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') in_e2_altnp = os.path.join(get_home(), @@ -35,12 +40,12 @@ def test_execution_ihMT_no_option(script_runner): in_e3_altpn = os.path.join(get_home(), 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') - in_e1_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') - in_e2_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') - in_e3_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') in_e1_neg = os.path.join(get_home(), 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') @@ -56,25 +61,28 @@ def test_execution_ihMT_no_option(script_runner): in_e3_pos = os.path.join(get_home(), 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') - in_e1_t1w = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') - in_e2_t1w = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') - in_e3_t1w = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') # no option ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, - in_mask, + '--mask', in_mask, '--in_altnp', in_e1_altnp, in_e2_altnp, in_e3_altnp, '--in_altpn', in_e1_altpn, in_e2_altpn, in_e3_altpn, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, - in_e3_mtoff, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, '-f') assert ret.success @@ -84,6 +92,11 @@ def test_execution_ihMT_prefix(script_runner): in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + in_e1_altnp = os.path.join(get_home(), 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') in_e2_altnp = os.path.join(get_home(), @@ -98,12 +111,12 @@ def test_execution_ihMT_prefix(script_runner): in_e3_altpn = os.path.join(get_home(), 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') - in_e1_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') - in_e2_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') - in_e3_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') in_e1_neg = os.path.join(get_home(), 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') @@ -119,37 +132,119 @@ def test_execution_ihMT_prefix(script_runner): in_e3_pos = os.path.join(get_home(), 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') - in_e1_t1w = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') - in_e2_t1w = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') - in_e3_t1w = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') # --out_prefix ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, - in_mask, + '--mask', in_mask, '--in_altnp', in_e1_altnp, in_e2_altnp, in_e3_altnp, '--in_altpn', in_e1_altpn, in_e2_altpn, in_e3_altpn, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, - in_e3_mtoff, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, '--out_prefix', 'sub_01', '-f') assert ret.success +def test_execution_ihMT_extended(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + + in_e1_altnp = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') + in_e2_altnp = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altnp_ihmt.nii.gz') + in_e3_altnp = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altnp_ihmt.nii.gz') + + in_e1_altpn = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altpn_ihmt.nii.gz') + in_e2_altpn = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altpn_ihmt.nii.gz') + in_e3_altpn = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') + + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + + in_e1_neg = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') + in_e2_neg = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-neg_ihmt.nii.gz') + in_e3_neg = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-neg_ihmt.nii.gz') + + in_e1_pos = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-pos_ihmt.nii.gz') + in_e2_pos = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-pos_ihmt.nii.gz') + in_e3_pos = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') + + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + + # --extended + ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, + '--mask', in_mask, + '--in_altnp', in_e1_altnp, in_e2_altnp, + in_e3_altnp, + '--in_altpn', in_e1_altpn, in_e2_altpn, + in_e3_altpn, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, + '--in_negative', in_e1_neg, in_e2_neg, + in_e3_neg, + '--in_positive', in_e1_pos, in_e2_pos, + in_e3_pos, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, + '--extended', + '-f') + assert ret.success + + def test_execution_ihMT_filtering(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + in_e1_altnp = os.path.join(get_home(), 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') in_e2_altnp = os.path.join(get_home(), @@ -164,12 +259,12 @@ def test_execution_ihMT_filtering(script_runner): in_e3_altpn = os.path.join(get_home(), 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') - in_e1_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') - in_e2_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') - in_e3_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') in_e1_neg = os.path.join(get_home(), 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') @@ -185,25 +280,28 @@ def test_execution_ihMT_filtering(script_runner): in_e3_pos = os.path.join(get_home(), 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') - in_e1_t1w = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') - in_e2_t1w = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') - in_e3_t1w = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') # --filtering ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, - in_mask, + '--mask', in_mask, '--in_altnp', in_e1_altnp, in_e2_altnp, in_e3_altnp, '--in_altpn', in_e1_altpn, in_e2_altpn, in_e3_altpn, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, - in_e3_mtoff, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, '--out_prefix', 'sub-01', '--filtering', '-f') @@ -215,6 +313,11 @@ def test_execution_ihMT_B1_map(script_runner): in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + in_e1_altnp = os.path.join(get_home(), 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') in_e2_altnp = os.path.join(get_home(), @@ -229,12 +332,12 @@ def test_execution_ihMT_B1_map(script_runner): in_e3_altpn = os.path.join(get_home(), 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') - in_e1_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') - in_e2_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') - in_e3_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') in_e1_neg = os.path.join(get_home(), 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') @@ -250,38 +353,54 @@ def test_execution_ihMT_B1_map(script_runner): in_e3_pos = os.path.join(get_home(), 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') - in_e1_t1w = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') - in_e2_t1w = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') - in_e3_t1w = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_b1_map = os.path.join(get_home(), 'ihMT', 'B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') - # --filtering ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, - in_mask, + '--mask', in_mask, '--in_altnp', in_e1_altnp, in_e2_altnp, in_e3_altnp, '--in_altpn', in_e1_altpn, in_e2_altpn, in_e3_altpn, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, - in_e3_mtoff, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, '--out_prefix', 'sub-01', - '--in_B1_map', in_b1_map, + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'empiric', + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, '-f') assert ret.success -def test_execution_ihMT_single_echo(script_runner): +def test_execution_ihMT_B1_no_T1(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + in_e1_altnp = os.path.join(get_home(), 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') in_e2_altnp = os.path.join(get_home(), @@ -296,12 +415,12 @@ def test_execution_ihMT_single_echo(script_runner): in_e3_altpn = os.path.join(get_home(), 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') - in_e1_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') - in_e2_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') - in_e3_mtoff = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') in_e1_neg = os.path.join(get_home(), 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') @@ -317,27 +436,258 @@ def test_execution_ihMT_single_echo(script_runner): in_e3_pos = os.path.join(get_home(), 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') - in_e1_t1w = os.path.join(get_home(), - 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') - in_e2_t1w = os.path.join(get_home(), - 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') - in_e3_t1w = os.path.join(get_home(), - 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + in_b1_map = os.path.join(get_home(), + 'ihMT', 'B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') - # --out_prefix ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, - in_mask, + '--mask', in_mask, '--in_altnp', in_e1_altnp, in_e2_altnp, in_e3_altnp, '--in_altpn', in_e1_altpn, in_e2_altpn, in_e3_altpn, - '--in_mtoff', in_e1_mtoff, in_e2_mtoff, - in_e3_mtoff, - '--in_negative', in_e1_neg, in_e2_neg, - in_e3_neg, - '--in_positive', in_e1_pos, in_e2_pos, - in_e3_pos, - '--in_t1w', in_e1_t1w, in_e2_t1w, in_e3_t1w, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, + '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, + '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'empiric', + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, + '-f') + assert ret.success + + +def test_execution_ihMT_wrong_echoes(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + + in_e1_altnp = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') + in_e2_altnp = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altnp_ihmt.nii.gz') + in_e3_altnp = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altnp_ihmt.nii.gz') + + in_e1_altpn = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altpn_ihmt.nii.gz') + in_e2_altpn = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altpn_ihmt.nii.gz') + in_e3_altpn = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') + + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + + in_e1_neg = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') + in_e2_neg = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-neg_ihmt.nii.gz') + in_e3_neg = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-neg_ihmt.nii.gz') + + in_e1_pos = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-pos_ihmt.nii.gz') + in_e2_pos = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-pos_ihmt.nii.gz') + in_e3_pos = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') + + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + + ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, + '--mask', in_mask, + '--in_altnp', in_e1_altnp, in_e2_altnp, + in_e3_altnp, + '--in_altpn', in_e1_altpn, in_e2_altpn, + in_e3_altpn, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, + '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, + '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, + '-f') + assert (not ret.success) + + +def test_execution_ihMT_B1_no_fit(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + + in_e1_altnp = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') + in_e2_altnp = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altnp_ihmt.nii.gz') + in_e3_altnp = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altnp_ihmt.nii.gz') + + in_e1_altpn = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altpn_ihmt.nii.gz') + in_e2_altpn = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-altpn_ihmt.nii.gz') + in_e3_altpn = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-altpn_ihmt.nii.gz') + + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + in_e2_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-mtoff_ihmt.nii.gz') + in_e3_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-mtoff_ihmt.nii.gz') + + in_e1_neg = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') + in_e2_neg = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-neg_ihmt.nii.gz') + in_e3_neg = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-neg_ihmt.nii.gz') + + in_e1_pos = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-pos_ihmt.nii.gz') + in_e2_pos = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-pos_ihmt.nii.gz') + in_e3_pos = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-pos_ihmt.nii.gz') + + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + in_e2_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-2_acq-T1w_ihmt.nii.gz') + in_e3_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-3_acq-T1w_ihmt.nii.gz') + + in_b1_map = os.path.join(get_home(), + 'ihMT', 'B1map.nii.gz') + in_b1_json = os.path.join(get_home(), + 'MT', 'sub-001_run-01_B1map.json') + out_b1_map = tmp_dir.name + '/B1map.nii.gz' + + # Temporary trick to have the B1 map with proper header. + ret = script_runner.run('scil_mti_adjust_B1_header.py', in_b1_map, + out_b1_map, in_b1_json, '-f') + + ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, + '--mask', in_mask, + '--in_altnp', in_e1_altnp, in_e2_altnp, + in_e3_altnp, + '--in_altpn', in_e1_altpn, in_e2_altpn, + in_e3_altpn, + '--in_mtoff_pd', in_e1_mtoff_pd, in_e2_mtoff_pd, + in_e3_mtoff_pd, + '--in_negative', in_e1_neg, in_e2_neg, in_e3_neg, + '--in_positive', in_e1_pos, in_e2_pos, in_e3_pos, + '--in_mtoff_t1', in_e1_mtoff_t1, in_e2_mtoff_t1, + in_e3_mtoff_t1, + '--out_prefix', 'sub-01', + '--in_B1_map', out_b1_map, + '--B1_correction_method', 'model_based', + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, + '-f') + assert (not ret.success) + + +def test_execution_ihMT_single_echo(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + + in_mtoff_pd_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.json') + in_mtoff_t1_json = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.json') + + in_e1_altnp = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') + + in_e1_altpn = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altpn_ihmt.nii.gz') + + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + + in_e1_neg = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') + + in_e1_pos = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-pos_ihmt.nii.gz') + + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + + ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, + '--mask', in_mask, + '--in_altnp', in_e1_altnp, + '--in_altpn', in_e1_altpn, + '--in_mtoff_pd', in_e1_mtoff_pd, + '--in_negative', in_e1_neg, + '--in_positive', in_e1_pos, + '--in_mtoff_t1', in_e1_mtoff_t1, + '--out_prefix', 'sub_01', + '--in_jsons', in_mtoff_pd_json, + in_mtoff_t1_json, '-f') + assert ret.success + + +def test_execution_ihMT_acq_params(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + + in_mask = os.path.join(get_home(), 'ihMT', 'mask_resample.nii.gz') + + in_e1_altnp = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altnp_ihmt.nii.gz') + + in_e1_altpn = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-altpn_ihmt.nii.gz') + + in_e1_mtoff_pd = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-mtoff_ihmt.nii.gz') + + in_e1_neg = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-neg_ihmt.nii.gz') + + in_e1_pos = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-pos_ihmt.nii.gz') + + in_e1_mtoff_t1 = os.path.join(get_home(), + 'ihMT', 'echo-1_acq-T1w_ihmt.nii.gz') + + ret = script_runner.run('scil_mti_maps_ihMT.py', tmp_dir.name, + '--mask', in_mask, + '--in_altnp', in_e1_altnp, + '--in_altpn', in_e1_altpn, + '--in_mtoff_pd', in_e1_mtoff_pd, + '--in_negative', in_e1_neg, + '--in_positive', in_e1_pos, + '--in_mtoff_t1', in_e1_mtoff_t1, '--out_prefix', 'sub_01', - '--single_echo', '-f') + '--in_acq_parameters', '15', '15', '0.1', '0.1', + '-f') assert ret.success diff --git a/scripts/tests/test_search_keywords.py b/scripts/tests/test_search_keywords.py index 600bd3402..5eda2b4fc 100644 --- a/scripts/tests/test_search_keywords.py +++ b/scripts/tests/test_search_keywords.py @@ -5,3 +5,18 @@ def test_help_option(script_runner): ret = script_runner.run('scil_search_keywords.py', '--help') assert ret.success + + +def test_no_verbose(script_runner): + ret = script_runner.run('scil_search_keywords.py', 'mti') + assert ret.success + + +def test_verbose_option(script_runner): + ret = script_runner.run('scil_search_keywords.py', 'mti', '-v') + assert ret.success + + +def test_not_find(script_runner): + ret = script_runner.run('scil_search_keywords.py', 'toto') + assert ret.success diff --git a/scripts/tests/test_tractogram_segment_bundles.py b/scripts/tests/test_tractogram_segment_bundles.py index 4cdbce030..bc7ace46e 100644 --- a/scripts/tests/test_tractogram_segment_bundles.py +++ b/scripts/tests/test_tractogram_segment_bundles.py @@ -39,5 +39,5 @@ def test_execution_bundles(script_runner): in_tractogram, 'config.json', in_models, in_aff, '--inverse', - '--processes', '1', '--log', 'WARNING') + '--processes', '1', '-v', 'WARNING') assert ret.success