Added customisable percentages to AP_rise_time.

efel/DependencyV5.txt

@@ -38,7 +38,7 @@ LibV5:AP_end_indices        #LibV5:peak_indices #LibV1:interpolate
 LibV2:AP_fall_indices       #LibV5:peak_indices     #LibV5:AP_begin_indices     #LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_duration	        #LibV5:AP_begin_indices #LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_duration_half_width	#LibV2:AP_rise_indices	#LibV2:AP_fall_indices #LibV1:interpolate 
-LibV2:AP_rise_time	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
+LibV2:AP_rise_time	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:AP_amplitude #LibV1:interpolate 
 LibV2:AP_fall_time	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 
 LibV2:AP_rise_rate	#LibV5:AP_begin_indices	#LibV5:peak_indices #LibV1:interpolate 
 LibV2:AP_fall_rate	#LibV5:peak_indices	#LibV5:AP_end_indices #LibV1:interpolate 

efel/api.py

@@ -86,6 +86,8 @@ def reset():
     setDoubleSetting('voltage_base_end_perc', 1.0)
     setDoubleSetting('current_base_start_perc', 0.9)
     setDoubleSetting('current_base_end_perc', 1.0)
+    setDoubleSetting('rise_start_perc', 0.0)
+    setDoubleSetting('rise_end_perc', 1.0)
     setDoubleSetting("initial_perc", 0.1)
     setDoubleSetting("min_spike_height", 20.0)
     setIntSetting("strict_stiminterval", 0)

efel/cppcore/LibV2.cpp

@@ -369,13 +369,42 @@ int LibV2::AP_duration_half_width(mapStr2intVec& IntFeatureData,
 // end of AP_duration_half_width
 
 // *** AP_rise_time according to E9 and E17 ***
-static int __AP_rise_time(const vector<double>& t,
+static int __AP_rise_time(const vector<double>& t, 
+                          const vector<double>& v,
                           const vector<int>& apbeginindices,
                           const vector<int>& peakindices,
+                          const vector<double>& apamplitude,
+                          double beginperc,
+                          double endperc,
                           vector<double>& aprisetime) {
   aprisetime.resize(std::min(apbeginindices.size(), peakindices.size()));
+  double begin_v;
+  double end_v;
+  double begin_indice;
+  double end_indice;
   for (size_t i = 0; i < aprisetime.size(); i++) {
-    aprisetime[i] = t[peakindices[i]] - t[apbeginindices[i]];
+    begin_v = v[apbeginindices[i]] + beginperc * apamplitude[i];
+    end_v = v[apbeginindices[i]] + endperc * apamplitude[i];
+
+    // Get begin indice
+    size_t j=apbeginindices[i];
+    // change slightly begin_v for almost equal case
+    // truncature error can change begin_v even when beginperc == 0.0
+    while (j<peakindices[i] && v[j] < begin_v - 0.0000000000001){
+      j++;
+    }
+    begin_indice = j;
+
+    // Get end indice
+    j=peakindices[i];
+    // change slightly end_v for almost equal case
+    // truncature error can change end_v even when beginperc == 0.0
+    while (j>apbeginindices[i] && v[j] > end_v + 0.0000000000001){
+      j--;
+    }
+    end_indice = j;
+
+    aprisetime[i] = t[end_indice] - t[begin_indice];
   }
   return aprisetime.size();
 }
@@ -400,8 +429,36 @@ int LibV2::AP_rise_time(mapStr2intVec& IntFeatureData,
   retval = getVec(IntFeatureData, StringData, "peak_indices",
                      peakindices);
   if (retval < 0) return -1;
+  vector<double> v;
+  retval = getVec(DoubleFeatureData, StringData, "V", v);
+  if (retval < 0) return -1;
+  vector<double> AP_amplitude;
+  retval =
+      getVec(DoubleFeatureData, StringData, "AP_amplitude", AP_amplitude);
+  if (retval < 0) {
+    GErrorStr += "Error calculating AP_amplitude for mean_AP_amplitude";
+    return -1;
+  } else if (retval == 0) {
+    GErrorStr += "No spikes found when calculating mean_AP_amplitude";
+    return -1;
+  } else if (AP_amplitude.size() == 0) {
+    GErrorStr += "No spikes found when calculating mean_AP_amplitude";
+    return -1;
+  }
+  // Get rise begin percentage
+  vector<double> risebeginperc;
+  retval = getVec(DoubleFeatureData, StringData, "rise_start_perc", risebeginperc);
+  if (retval <= 0) {
+    risebeginperc.push_back(0.0);
+  }
+  // Get rise end percentage
+  vector<double> riseendperc;
+  retval = getVec(DoubleFeatureData, StringData, "rise_end_perc", riseendperc);
+  if (retval <= 0) {
+    riseendperc.push_back(1.0);
+  }
   vector<double> aprisetime;
-  retval = __AP_rise_time(t, apbeginindices, peakindices, aprisetime);
+  retval = __AP_rise_time(t, v, apbeginindices, peakindices, AP_amplitude, risebeginperc[0], riseendperc[0], aprisetime);
   if (retval >= 0) {
     setVec(DoubleFeatureData, StringData, "AP_rise_time", aprisetime);
   }

efel/tests/test_basic.py

@@ -1963,3 +1963,62 @@ def test_unfinished_peak():
     spikecount = traces_results[0]['Spikecount'][0]
 
     nt.assert_equal(spikecount, 3)
+
+
+def rise_time_perc(
+    time, voltage,
+    AP_begin_indices,
+    peak_indices,
+    rise_start_perc,
+    rise_end_perc
+):
+    """AP_rise_time numpy implementation with percentages"""
+    rise_times = []
+    AP_amp = voltage[peak_indices] - voltage[AP_begin_indices]
+    begin_voltages = AP_amp * rise_start_perc + voltage[AP_begin_indices]
+    end_voltages = AP_amp * rise_end_perc + voltage[AP_begin_indices]
+
+    for AP_begin_indice, peak_indice, begin_v, end_v in zip(
+        AP_begin_indices, peak_indices, begin_voltages, end_voltages
+    ):
+        voltage_window = voltage[AP_begin_indice:peak_indice]
+
+        new_begin_indice = AP_begin_indice + numpy.min(
+            numpy.where(voltage_window >= begin_v)[0]
+        )
+        new_end_indice = AP_begin_indice + numpy.max(
+            numpy.where(voltage_window <= end_v)[0]
+        )
+
+        rise_times.append(time[new_end_indice] - time[new_begin_indice])
+
+    return numpy.array(rise_times)
+
+
+def test_rise_time_perc():
+    """basic: Test AP rise time percentage"""
+
+    import efel
+    efel.reset()
+    trace, time, voltage, stim_start, stim_end = load_data(
+        'mean_frequency1', interp=True
+    )
+
+    trace['rise_start_perc'] = [0.2]
+    trace['rise_end_perc'] = [0.8]
+
+    features = ['AP_rise_time', 'AP_begin_indices', 'peak_indices']
+
+    feature_values = efel.getFeatureValues(
+        [trace], features, raise_warnings=False
+    )
+    ap_rise_time = feature_values[0]['AP_rise_time']
+    AP_begin_indices = feature_values[0]['AP_begin_indices']
+    peak_indices = feature_values[0]['peak_indices']
+
+    expected = rise_time_perc(
+        time, voltage, AP_begin_indices, peak_indices, 0.2, 0.8
+    )
+
+    for exp, rise_time in zip(expected, ap_rise_time):
+        nt.assert_almost_equal(exp, rise_time)