Merge pull request #1089 from AntoineTheb/atheb/s0_div_0

FIX: Division by 0 in `scil_dti_metrics`
scilus · Dec 17, 2024 · 0817790 · 0817790
2 parents cbb3501 + 3d29717
commit 0817790
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 3 deletions.
diff --git a/scripts/scil_dti_metrics.py b/scripts/scil_dti_metrics.py
@@ -340,6 +340,7 @@ def main():
                          "mask was provided.")
         # Mean residual image
         S0 = np.mean(data[..., gtab.b0s_mask], axis=-1)
+
         tenfit2_predict = np.zeros(data.shape, dtype=np.float32)
 
         for i in range(data.shape[0]):
@@ -348,7 +349,8 @@ def main():
             else:
                 tenfit2 = tenmodel.fit(data[i, :, :, :])
 
-            tenfit2_predict[i, :, :, :] = tenfit2.predict(gtab, S0[i, :, :])
+            S0_i = np.maximum(S0[i, :, :], tenfit2.model.min_signal)
+            tenfit2_predict[i, :, :, :] = tenfit2.predict(gtab, S0_i)
 
         R, data_diff = compute_residuals(
             predicted_data=tenfit2_predict.astype(np.float32),
@@ -370,7 +372,7 @@ def main():
 
         # Plotting and saving figure
         plot_residuals(data_diff, mask, R_k, q1, q3, iqr,
-                        residual_basename)
+                       residual_basename)
 
 
 if __name__ == "__main__":

diff --git a/scripts/tests/test_dti_metrics.py b/scripts/tests/test_dti_metrics.py
@@ -17,7 +17,7 @@ def test_help_option(script_runner):
     assert ret.success
 
 
-def test_execution_processing(script_runner, monkeypatch):
+def test_execution_processing_diff_metrics(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
     in_dwi = os.path.join(SCILPY_HOME, 'processing', 'dwi_crop_1000.nii.gz')
     in_bval = os.path.join(SCILPY_HOME, 'processing', '1000.bval')
@@ -39,11 +39,37 @@ def test_execution_processing(script_runner, monkeypatch):
                             '--mask', mask_uint8)
     assert ret.success
 
+
+def test_execution_processing_b0_threshold(script_runner, monkeypatch):
+    monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi = os.path.join(SCILPY_HOME, 'processing', 'dwi_crop_1000.nii.gz')
+    in_bval = os.path.join(SCILPY_HOME, 'processing', '1000.bval')
+    in_bvec = os.path.join(SCILPY_HOME, 'processing', '1000.bvec')
+
+    # No mask fitting with this data? Creating our own.
+    mask = os.path.join(SCILPY_HOME, 'processing', 'ad.nii.gz')
+    mask_uint8 = os.path.join('mask_uint8.nii.gz')
+    script_runner.run('scil_volume_math.py', 'convert',
+                      mask, mask_uint8, '--data_type', 'uint8')
+
     ret = script_runner.run('scil_dti_metrics.py', in_dwi,
                             in_bval, in_bvec, '--not_all',
                             '--fa', 'fa.nii.gz', '--b0_threshold', '1', '-f')
     assert not ret.success
 
+
+def test_execution_processing_rgb(script_runner, monkeypatch):
+    monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
+    in_dwi = os.path.join(SCILPY_HOME, 'processing', 'dwi_crop_1000.nii.gz')
+    in_bval = os.path.join(SCILPY_HOME, 'processing', '1000.bval')
+    in_bvec = os.path.join(SCILPY_HOME, 'processing', '1000.bvec')
+
+    # No mask fitting with this data? Creating our own.
+    mask = os.path.join(SCILPY_HOME, 'processing', 'ad.nii.gz')
+    mask_uint8 = os.path.join('mask_uint8.nii.gz')
+    script_runner.run('scil_volume_math.py', 'convert',
+                      mask, mask_uint8, '--data_type', 'uint8')
+
     ret = script_runner.run('scil_dti_metrics.py', in_dwi,
                             in_bval, in_bvec, '--not_all',
                             '--rgb', 'rgb.nii.gz', '-f')