Merge pull request #3969 from opensim-org/trim_filtered_table_fix

Trim filtered table to original time range in TabOpLowPassFilter

CHANGELOG.md

@@ -47,6 +47,8 @@ v4.6
 - Added the templatized `MocoStudy::analyze<T>()` and equivalent scripting counterparts: `analyzeVec3`, `analyzeSpatialVec`, `analyzeRotation`. (#3940)
 - Added `ConstantCurvatureJoint` to the SWIG bindings; it is now available in Matlab and Python (#3957). 
 - Added methods and `Output`s for calculating the angular momentum of a `Body`. (#3962)
+- Updated `TabOpLowPassFilter` so that by default the processed table is trimmed to the original time range after padding and filtering. 
+  The original behavior (no trimming) can be enabled via the new property `trim_to_original_time_range`. (#3969)
 
 
 v4.5.1

OpenSim/Moco/Test/testMocoInverse.cpp

@@ -222,23 +222,20 @@ TEST_CASE("MocoInverse Rajagopal2016, 18 muscles", "[casadi]") {
                 "/jointset/back/lumbar_extension/speed"));
     }
 }
-// Next test_case fails on linux while parsing .sto file, disabling for now 
-#ifdef _WIN32
+
 TEST_CASE("Test IMUDataReporter for gait") {
 
     // Compute accelerometer signals from MocoInverse solution.
     ModelProcessor modelProcessor =
-        ModelProcessor("subject_walk_armless_18musc.osim") |
-        ModOpReplaceJointsWithWelds({"subtalar_r", "subtalar_l",
-            "mtp_r", "mtp_l"}) |
-        ModOpReplaceMusclesWithDeGrooteFregly2016() |
-        ModOpIgnorePassiveFiberForcesDGF() |
-        ModOpTendonComplianceDynamicsModeDGF("implicit") |
-        ModOpAddExternalLoads("subject_walk_armless_external_loads.xml");
-    std::vector<std::string> paths = {"/bodyset/pelvis",
-                                       "/bodyset/femur_r",
-                                       "/bodyset/tibia_r",
-                                       "/bodyset/calcn_r"};
+            ModelProcessor("subject_walk_armless_18musc.osim") |
+            ModOpReplaceJointsWithWelds({"subtalar_r", "subtalar_l",
+                "mtp_r", "mtp_l"}) |
+            ModOpReplaceMusclesWithDeGrooteFregly2016() |
+            ModOpIgnorePassiveFiberForcesDGF() |
+            ModOpTendonComplianceDynamicsModeDGF("implicit") |
+            ModOpAddExternalLoads("subject_walk_armless_external_loads.xml");
+    std::vector<std::string> paths = {"/bodyset/pelvis", "/bodyset/femur_r",
+            "/bodyset/tibia_r",  "/bodyset/calcn_r"};
     Model model = modelProcessor.process();
     OpenSenseUtilities().addModelIMUs(model, paths);
     model.initSystem();
@@ -247,15 +244,17 @@ TEST_CASE("Test IMUDataReporter for gait") {
             analyzeMocoTrajectory<SimTK::Vec3>(model, std,
                                                {".*accelerometer_signal"});
     STOFileAdapter_<SimTK::Vec3>::write(accelSignals,
-        "testMocoInverse_accelerometer_signals.sto");
+            "testMocoInverse_accelerometer_signals.sto");
 
     // Compute accelerometer signals with no forces and a PositionMotion created
     // from the model coordinates.
 
     // Create the coordinates table (in radians).
     auto tableProcessor =
             TableProcessor("subject_walk_armless_coordinates.mot") |
-            TabOpLowPassFilter(6) |
+            // Trimming to the original time range increases the acceleration
+            // error slightly, so we disable it here.
+            TabOpLowPassFilter(6, false) |
             TabOpUseAbsoluteStateNames();
     auto coordinatesRadians = tableProcessor.processAndConvertToRadians(model);
     for (const auto& label : coordinatesRadians.getColumnLabels()) {
@@ -272,11 +271,11 @@ TEST_CASE("Test IMUDataReporter for gait") {
                           "subject_walk_armless_coordinates_radians.sto");
     // Create a model with no muscles (or other forces) and add IMU components.
     ModelProcessor modelProcessorNoMuscles =
-        ModelProcessor("subject_walk_armless_18musc.osim") |
-        ModOpReplaceJointsWithWelds({"subtalar_r", "subtalar_l",
-            "mtp_r", "mtp_l"}) |
-        ModOpRemoveMuscles() |
-        ModOpAddExternalLoads("subject_walk_armless_external_loads.xml");
+            ModelProcessor("subject_walk_armless_18musc.osim") |
+            ModOpReplaceJointsWithWelds({"subtalar_r", "subtalar_l",
+                                         "mtp_r", "mtp_l"}) |
+            ModOpRemoveMuscles() |
+            ModOpAddExternalLoads("subject_walk_armless_external_loads.xml");
     Model modelNoMuscles = modelProcessorNoMuscles.process();
     modelNoMuscles.updForceSet().clearAndDestroy();
     OpenSenseUtilities().addModelIMUs(modelNoMuscles, paths);
@@ -285,7 +284,7 @@ TEST_CASE("Test IMUDataReporter for gait") {
 
     // Add a PositionMotion to the model based on the coordinate trajectories.
     auto statesTraj = StatesTrajectory::createFromStatesTable(model,
-          coordinatesRadians, true, true, true);
+            coordinatesRadians, true, true, true);
     auto posmot = PositionMotion::createFromStatesTrajectory(model, statesTraj);
     posmot->setName("position_motion");
     modelNoMuscles.addComponent(posmot.release());
@@ -319,12 +318,12 @@ TEST_CASE("Test IMUDataReporter for gait") {
     // Compute error in accelerometer signals.
     auto accelSignalsFlat = accelSignals.flatten();
     auto accelBlock = accelSignalsFlat.getMatrixBlock(0, 0,
-        accelSignalsFlat.getNumRows(),
-        accelSignalsFlat.getNumColumns());
+            accelSignalsFlat.getNumRows(),
+            accelSignalsFlat.getNumColumns());
     auto accelBlockNoMuscles =
-        accelSignalsNoMusclesDownSampledFlat.getMatrixBlock(0, 0,
-        accelSignalsNoMusclesDownSampledFlat.getNumRows(),
-        accelSignalsNoMusclesDownSampledFlat.getNumColumns());
+            accelSignalsNoMusclesDownSampledFlat.getMatrixBlock(0, 0,
+            accelSignalsNoMusclesDownSampledFlat.getNumRows(),
+            accelSignalsNoMusclesDownSampledFlat.getNumColumns());
     auto diff = accelBlock - accelBlockNoMuscles;
     auto diffSqr = diff.elementwiseMultiply(diff);
     auto sumSquaredError = diffSqr.rowSum();
@@ -391,4 +390,3 @@ TEST_CASE("Test IMUDataReporter for gait") {
     // TODO The error when using AnalyzeTool is slightly larger.
     SimTK_TEST_EQ_TOL(integratedSumSquaredErrorAnalyze, 0.0, 1e-3);
 }
-#endif

OpenSim/Simulation/TableProcessor.h

@@ -183,19 +183,36 @@ class OSIMSIMULATION_API TabOpLowPassFilter : public TableOperator {
     OpenSim_DECLARE_PROPERTY(cutoff_frequency, double,
             "Low-pass cutoff frequency (Hz) (default is -1, which means no "
             "filtering).");
-    TabOpLowPassFilter() { constructProperty_cutoff_frequency(-1); }
+    OpenSim_DECLARE_PROPERTY(trim_to_original_time_range, bool,
+            "Trim the rows of the output table to match the original table's "
+            "time range after filtering (default: true).");
+    TabOpLowPassFilter() { 
+        constructProperty_cutoff_frequency(-1); 
+        constructProperty_trim_to_original_time_range(true);
+    }
     TabOpLowPassFilter(double cutoffFrequency) : TabOpLowPassFilter() {
         set_cutoff_frequency(cutoffFrequency);
     }
+    TabOpLowPassFilter(double cutoffFrequency, bool trimToOriginalTimeRange) : 
+            TabOpLowPassFilter() {
+        set_cutoff_frequency(cutoffFrequency);
+        set_trim_to_original_time_range(trimToOriginalTimeRange);
+    }
     void operate(TimeSeriesTable& table, const Model* model = nullptr)
             const override {
         if (get_cutoff_frequency() != -1) {
             OPENSIM_THROW_IF(get_cutoff_frequency() <= 0, Exception,
                     "Expected cutoff frequency to be positive, but got {}.",
                     get_cutoff_frequency());
 
+            const auto& times = table.getIndependentColumn();
+            double initialTime = times.front();
+            double finalTime = times.back();
             TableUtilities::filterLowpass(
                     table, get_cutoff_frequency(), true);
+            if (get_trim_to_original_time_range()) {
+                table.trim(initialTime, finalTime);
+            }
         }
     }
 };

OpenSim/Simulation/Test/testTableProcessor.cpp

@@ -76,3 +76,13 @@ TEST_CASE("TableProcessor") {
         }
     }
 }
+
+TEST_CASE("TabOpLowPassFilter retains original time range") {
+    TimeSeriesTable table(std::vector<double>{1, 2, 3, 4, 5, 6},
+            SimTK::Test::randMatrix(6, 2), std::vector<std::string>{"a", "b"});
+    TableProcessor proc = TableProcessor(table) | TabOpLowPassFilter(6);
+    TimeSeriesTable out = proc.process();
+    CHECK(out.getNumRows() == 6);
+    CHECK(out.getIndependentColumn().front() == 1);
+    CHECK(out.getIndependentColumn().back() == 6);
+}