Improve step detection (#122)

BlueBrain · Jan 18, 2023 · bcdf986 · bcdf986
1 parent 3cc23d6
commit bcdf986
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 25 deletions.
diff --git a/bluepyefe/ecode/step.py b/bluepyefe/ecode/step.py
@@ -79,8 +79,9 @@ def interpret(self, t, current, config_data, reader_data):
         # Smooth the current
         smooth_current = None
 
-        # Set the threshold to look for steps
-        epsilon = 0.005
+        # Set the threshold to detect the step
+        noise_level = numpy.std(numpy.concatenate((self.current[:50], self.current[-50:])))
+        step_threshold = numpy.clip(2. * noise_level, 1e-5, numpy.max(self.current))
 
         # The buffer prevent miss-detection of the step when artifacts are
         # present at the very start or very end of the current trace
@@ -124,7 +125,7 @@ def interpret(self, t, current, config_data, reader_data):
             if smooth_current is None:
                 smooth_current = scipy_signal2d(current, 85)
             _ = numpy.abs(smooth_current[idx_buffer:] - self.hypamp)
-            self.ton = idx_buffer + numpy.argmax(_ > epsilon)
+            self.ton = idx_buffer + numpy.argmax(_ > step_threshold)
 
         else:
             # Infer the base current hypamp
@@ -137,7 +138,7 @@ def interpret(self, t, current, config_data, reader_data):
                 numpy.abs(smooth_current[:-idx_buffer] - self.hypamp)
             )
             self.toff = (
-                (len(current) - numpy.argmax(_ > epsilon)) - 1 - idx_buffer
+                (len(current) - numpy.argmax(_ > step_threshold)) - 1 - idx_buffer
             )
 
         # Get the amplitude of the step current (relative to hypamp)

diff --git a/bluepyefe/recording.py b/bluepyefe/recording.py
@@ -299,7 +299,10 @@ def compute_spikecount(self, efel_settings=None, offset_voltage=20.):
             idx_toff = self.ms_to_index(self.toff)
             step_voltage = numpy.median(self.voltage[idx_ton:idx_toff])
             base_voltage = numpy.median(self.voltage[:idx_ton])
-            thresh = step_voltage + offset_voltage
+            if base_voltage > step_voltage:
+                thresh = base_voltage + offset_voltage
+            else:
+                thresh = step_voltage + offset_voltage
             tmp_settings["Threshold"] = thresh
             efel_vals = self.call_efel(['peak_time'], tmp_settings)
             self.peak_time = efel_vals[0]['peak_time']