From 9158312faa20dca46858869f2cca382f683424a1 Mon Sep 17 00:00:00 2001
From: "Billah, Tashrif" <tbillah@bwh.harvard.edu>
Date: Thu, 31 Dec 2020 17:48:36 -0500
Subject: [PATCH] determine bzero based on absolute s/mm^2

---
 ukf/dwi_normalize.cc | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/ukf/dwi_normalize.cc b/ukf/dwi_normalize.cc
index 7fe04c8..826dcbc 100644
--- a/ukf/dwi_normalize.cc
+++ b/ukf/dwi_normalize.cc
@@ -56,6 +56,17 @@ void dwiNormalize(const Nrrd *raw, Nrrd *& normalized)
     char *value;
     nrrdKeyValueIndex(raw, &key, &value, i);
     std::string keyStr(key);
+
+    // Obtain DWMRI_b-value
+    int bmax;
+    if( keyStr.length() == 13 && !keyStr.compare("DWMRI_b-value") )
+      {
+        sscanf(value, "%d", &bmax);
+        std::cout << "DWMRI_b-value " << bmax << std::endl;
+
+      }
+
+
     if( keyStr.length() > 14 && !keyStr.substr(0, 14).compare("DWMRI_gradient") )
       {
       float gx, gy, gz;
@@ -64,8 +75,12 @@ void dwiNormalize(const Nrrd *raw, Nrrd *& normalized)
         std::cout << "The gradients must have 3 components!" << std::endl;
         throw;
         }
+
       const float gradient_magnitude = sqrt(gx * gx + gy * gy + gz * gz);
-      nonZeroGradientFlag.push_back( gradient_magnitude > 0.05F );
+
+      // See https://github.com/pnlbwh/ukftractography/issues/136#issuecomment-753173654
+      // for details about the following criteria
+      nonZeroGradientFlag.push_back( bmax * gradient_magnitude * gradient_magnitude > 50 );
       }
     keyValuePairsOfRaw.push_back(std::make_pair(std::string(key), std::string(value)) );
     free(key);   //Key and value generated by Nrrd by malloc, so free needs to be used.