From 4e3a0381565bbd544006ec080b3da3e803962ba6 Mon Sep 17 00:00:00 2001 From: barma7 Date: Tue, 15 Oct 2024 15:29:39 -0700 Subject: [PATCH 1/4] added b1 correction to qDESS T2 mapping (beta) --- dosma/core/numpy_routines.py | 2 +- dosma/scan_sequences/mri/qdess.py | 50 ++++++++++++++++++++++--------- 2 files changed, 37 insertions(+), 15 deletions(-) diff --git a/dosma/core/numpy_routines.py b/dosma/core/numpy_routines.py index a96e9a4..b887d64 100644 --- a/dosma/core/numpy_routines.py +++ b/dosma/core/numpy_routines.py @@ -161,7 +161,7 @@ def nan_to_num(x, copy=True, nan=0.0, posinf=None, neginf=None): return x._partial_clone(volume=vol) -@implements(np.around, np.round, np.round_) +@implements(np.around, np.round, np.round) def around(x, decimals=0, affine=False): """Round medical image pixel data (and optionally affine) to the given number of decimals. diff --git a/dosma/scan_sequences/mri/qdess.py b/dosma/scan_sequences/mri/qdess.py index a308a39..9481860 100644 --- a/dosma/scan_sequences/mri/qdess.py +++ b/dosma/scan_sequences/mri/qdess.py @@ -107,6 +107,7 @@ def generate_t2_map( tissue: Tissue = None, suppress_fat: bool = False, suppress_fluid: bool = False, + b1map: MedicalVolume = None, beta: float = 1.2, gl_area: float = None, tg: float = None, @@ -211,9 +212,15 @@ def generate_t2_map( # Flip Angle (degree -> radians) alpha = float(ref_dicom.FlipAngle) if alpha is None else alpha - alpha = math.radians(alpha) - if np.allclose(math.sin(alpha / 2), 0): - warnings.warn("sin(flip angle) is close to 0 - t2 map may fail.", stacklevel=2) + + if b1map: + alpha = np.multiply(self.b1map.volume, math.radians(alpha)) + if np.allclose(np.sin(alpha / 2), 0): + warnings.warn("sin(flip angle) is close to 0 - t2 map may fail.") + else: + alpha = math.radians(alpha) + if np.allclose(math.sin(alpha / 2), 0): + warnings.warn("sin(flip angle) is close to 0 - t2 map may fail.") mask = xp.ones([r, c, num_slices]) @@ -228,22 +235,37 @@ def generate_t2_map( dkL = gamma * Gl * Tg # Simply math - k = ( - xp.power((xp.sin(alpha / 2)), 2) - * (1 + xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) - / (1 - xp.cos(alpha) * xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) - ) + if b1map: + # treat k as an array + k = np.square((np.sin(alpha / 2))) * ( + 1 + np.exp(-TR / T1 - TR * np.square(dkL) * diffusivity)) / ( + 1 - np.cos(alpha) * np.exp(-TR / T1 - TR * np.square(dkL) * diffusivity)) + else: + k = ( + xp.power((xp.sin(alpha / 2)), 2) + * (1 + xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) + / (1 - xp.cos(alpha) * xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) + ) + ## END MODIFICATION c1 = (TR - Tg / 3) * (xp.power(dkL, 2)) * diffusivity else: # T2 fit - k = ( - xp.power((xp.sin(alpha / 2)), 2) - * (1 + xp.exp(-TR / T1)) - / (1 - xp.cos(alpha) * xp.exp(-TR / T1)) - ) - c1 = 0 + if b1map: + k = ( + np.square((np.sin(alpha / 2))) + * (1 + np.exp(-TR / T1)) + / (1 - np.cos(alpha) * np.exp(-TR / T1)) + ) + c1 = 0 + else: + k = ( + xp.power((xp.sin(alpha / 2)), 2) + * (1 + xp.exp(-TR / T1)) + / (1 - xp.cos(alpha) * xp.exp(-TR / T1)) + ) + c1 = 0 # have to divide division into steps to avoid overflow error t2map = -2000 * (TR - TE) / (xp.log(abs(ratio) / k) + c1) From d7e3273ca89b123deda8d6427ae340910af4ab0c Mon Sep 17 00:00:00 2001 From: barma7 Date: Thu, 24 Oct 2024 11:53:46 -0700 Subject: [PATCH 2/4] modified model_loading_util.py: yaml.load(f) --> yaml.load(f, Loader=yaml.FullLoader) --- dosma/models/model_loading_util.py | 2 +- dosma/scan_sequences/mri/qdess.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dosma/models/model_loading_util.py b/dosma/models/model_loading_util.py index 55c7310..faad8a7 100644 --- a/dosma/models/model_loading_util.py +++ b/dosma/models/model_loading_util.py @@ -72,7 +72,7 @@ def _gen_mask(func, *_args, **_kwargs): if isinstance(cfg_file_or_dict, str): with open(cfg_file_or_dict, "r") as f: - cfg = yaml.load(f) + cfg = yaml.load(f, Loader=yaml.FullLoader) else: cfg = cfg_file_or_dict diff --git a/dosma/scan_sequences/mri/qdess.py b/dosma/scan_sequences/mri/qdess.py index 9481860..ee8a857 100644 --- a/dosma/scan_sequences/mri/qdess.py +++ b/dosma/scan_sequences/mri/qdess.py @@ -214,7 +214,7 @@ def generate_t2_map( alpha = float(ref_dicom.FlipAngle) if alpha is None else alpha if b1map: - alpha = np.multiply(self.b1map.volume, math.radians(alpha)) + alpha = np.multiply(b1map.volume, math.radians(alpha)) if np.allclose(np.sin(alpha / 2), 0): warnings.warn("sin(flip angle) is close to 0 - t2 map may fail.") else: From e47bcdcc765cdc49e1c96b26bb82fd3f2672a4c3 Mon Sep 17 00:00:00 2001 From: barma7 Date: Thu, 31 Oct 2024 14:59:20 -0700 Subject: [PATCH 3/4] assigned unique names to oaiunet2d.py conv2d layers as newr Keras versions require unique name for layers. --- dosma/models/oaiunet2d.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dosma/models/oaiunet2d.py b/dosma/models/oaiunet2d.py index 5a919f7..05fba69 100644 --- a/dosma/models/oaiunet2d.py +++ b/dosma/models/oaiunet2d.py @@ -79,6 +79,7 @@ def __load_keras_model__(self, input_shape, force_weights=False): padding="same", activation="relu", kernel_initializer="he_normal", + name=f"conv2d_down_{depth_cnt}_1" # Unique name for each layer )(pool) conv = Conv2D( nfeatures[depth_cnt], @@ -86,6 +87,7 @@ def __load_keras_model__(self, input_shape, force_weights=False): padding="same", activation="relu", kernel_initializer="he_normal", + name=f"conv2d_down_{depth_cnt}_2" # Unique name for each layer )(conv) conv = BN(axis=-1, momentum=0.95, epsilon=0.001)(conv) @@ -134,6 +136,7 @@ def __load_keras_model__(self, input_shape, force_weights=False): padding="same", activation="relu", kernel_initializer="he_normal", + name=f"conv2d_up_{depth_cnt}_1" # Unique name for each layer )(up) conv = Conv2D( nfeatures[depth_cnt], @@ -141,6 +144,7 @@ def __load_keras_model__(self, input_shape, force_weights=False): padding="same", activation="relu", kernel_initializer="he_normal", + name=f"conv2d_up_{depth_cnt}_2" # Unique name for each layer )(conv) conv = BN(axis=-1, momentum=0.95, epsilon=0.001)(conv) From 05cd8b685d104c2cc2cce5ebb6268aad18f2521b Mon Sep 17 00:00:00 2001 From: barma7 Date: Thu, 31 Oct 2024 15:43:21 -0700 Subject: [PATCH 4/4] removed unecessary if statements in handling b1 correction --- dosma/scan_sequences/mri/qdess.py | 42 +++++++++++-------------------- 1 file changed, 14 insertions(+), 28 deletions(-) diff --git a/dosma/scan_sequences/mri/qdess.py b/dosma/scan_sequences/mri/qdess.py index ee8a857..341c864 100644 --- a/dosma/scan_sequences/mri/qdess.py +++ b/dosma/scan_sequences/mri/qdess.py @@ -214,8 +214,9 @@ def generate_t2_map( alpha = float(ref_dicom.FlipAngle) if alpha is None else alpha if b1map: + alpha_nominal = deepcopy(alpha) alpha = np.multiply(b1map.volume, math.radians(alpha)) - if np.allclose(np.sin(alpha / 2), 0): + if np.allclose(np.sin(alpha_nominal / 2), 0): warnings.warn("sin(flip angle) is close to 0 - t2 map may fail.") else: alpha = math.radians(alpha) @@ -235,38 +236,23 @@ def generate_t2_map( dkL = gamma * Gl * Tg # Simply math - if b1map: - # treat k as an array - k = np.square((np.sin(alpha / 2))) * ( - 1 + np.exp(-TR / T1 - TR * np.square(dkL) * diffusivity)) / ( - 1 - np.cos(alpha) * np.exp(-TR / T1 - TR * np.square(dkL) * diffusivity)) - else: - k = ( - xp.power((xp.sin(alpha / 2)), 2) - * (1 + xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) - / (1 - xp.cos(alpha) * xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) - ) - ## END MODIFICATION + k = ( + xp.power((xp.sin(alpha / 2)), 2) + * (1 + xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) + / (1 - xp.cos(alpha) * xp.exp(-TR / T1 - TR * xp.power(dkL, 2) * diffusivity)) + ) c1 = (TR - Tg / 3) * (xp.power(dkL, 2)) * diffusivity else: # T2 fit - if b1map: - k = ( - np.square((np.sin(alpha / 2))) - * (1 + np.exp(-TR / T1)) - / (1 - np.cos(alpha) * np.exp(-TR / T1)) - ) - c1 = 0 - else: - k = ( - xp.power((xp.sin(alpha / 2)), 2) - * (1 + xp.exp(-TR / T1)) - / (1 - xp.cos(alpha) * xp.exp(-TR / T1)) - ) - c1 = 0 - + k = ( + xp.power((xp.sin(alpha / 2)), 2) + * (1 + xp.exp(-TR / T1)) + / (1 - xp.cos(alpha) * xp.exp(-TR / T1)) + ) + c1 = 0 + # have to divide division into steps to avoid overflow error t2map = -2000 * (TR - TE) / (xp.log(abs(ratio) / k) + c1) t2map = xp.nan_to_num(t2map)