Motion tracking III: Triangulation via Pose2sim, Inverse kinematics & dynamics via OpenSim

Overview

In the previous script, we have extracted 2D keypoints for each trial (per each video, 3 in total). In this script, we will use Pose2sim (Pagnon et al., 2022) to calibrate the 3 cameras present in the lab setup, and triangulate the 3D position of the keypoints.

Demo of this pipeline has been published on EnvisionBOX

While in preregistration we used version 0.1.0, the final processing uses 0.10.39. The new version allows one to also perform inverse kinematics (i.e., joint angles) using OpenSim API (Seth et al., 2018). As it improves and simplifies the process and the scripting workflow, we deem this a legitimate deviation from the previous script.

For more information on Pose2sim, please refer to the Pose2sim documentation

Lastly, using OpenSim (Seth et al., 2018), we finish the workflow with calculating inverse dynamics (i.e., joint forces). In the preregistration, this was done in a separate script but given that we moved the inverse kinematics here, we perform the inverse dynamics here as well.

The documentation of OpenSim project is available here

from Pose2Sim import Pose2Sim
import os
import glob
import pandas as pd
from trc import TRCData
import pandas as pd
import shutil
import cv2
import numpy as np
import toml
import matplotlib.pyplot as plt
import re
from pathlib import Path
import scipy.signal
import opensim
from scipy.signal import savgol_filter
import xml.etree.ElementTree as ET

curfolder = os.getcwd()
print('Current folder is: ' + curfolder)

# Here is our config.file
pose2simprjfolder = os.path.join(curfolder, 'Pose2Sim', 'FLESH_setup')

# Here we store the data
inputfolder = os.path.join(curfolder, 'projectdata')
folderstotrack = glob.glob(os.path.join(inputfolder, '*'))
#print(folderstotrack)

# Here we store mass information (weight, height) about participants
META = pd.read_csv(os.path.join(curfolder, '..', '00_raw', 'Demographics_all', 'all_demodata.csv')) 

# Initiate empty list
pcnfolders = []

# Get all the folders per session, per participant
for i in folderstotrack:
    pcn1folders = glob.glob(os.path.join(i, '*'))
    pcn2folders = glob.glob(os.path.join(i, '*'))
    pcnfolders_in_session = pcn1folders + pcn2folders
    pcnfolders = pcnfolders + pcnfolders_in_session


# Get rid of all pontetially confusing files/folders
pcnfolders = [x for x in pcnfolders if 'Config' not in x]
pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
pcnfolders = [x for x in pcnfolders if 'xml' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
pcnfolders = [x for x in pcnfolders if 'sto' not in x]
pcnfolders = [x for x in pcnfolders if 'txt' not in x]
pcnfolders = [x for x in pcnfolders if 'calibration' not in x]
print(pcnfolders[0:10])


# Prepare special log txt
error_log = os.path.join(curfolder, 'errors', 'pose2sim_error_log.txt')
if not os.path.exists(error_log):
    with open(error_log, 'w') as f:
        f.write('Pose2Sim processing error log\n')
        f.write('============================\n\n')

Current folder is: f:\FLESH_ContinuousBodilyEffort\02_MotionTracking_processing
['f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_11_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_12_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_13_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_14_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_15_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_16_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_17_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_20_p0', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_21_p0', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_22_p0']

import copy
import logging
import traceback

def load_toml(file_path):
    """
    Loads a TOML file and returns its contents as a dictionary.

    Args:
        file_path (str): Path to the TOML file to be loaded.

    Returns:
        dict: The contents of the TOML file as a dictionary.
    """
    with open(file_path, 'r') as file:
        return toml.load(file)

def save_toml(data, file_path):
    """
    Saves a dictionary to a TOML file.

    Args:
        data (dict): The dictionary to be saved as a TOML file.
        file_path (str): Path to the TOML file to be created or overwritten.
    """
    with open(file_path, 'w') as file:
        toml.dump(data, file)

def update_participant_info(toml_data, height, mass):
    """
    Updates the participant information in the TOML data with the given height and mass.

    Args:
        toml_data (dict): The TOML data to be updated.
        height (float): The height of the participant.
        mass (float): The mass of the participant.

    Returns:
        dict: The updated TOML data with the new participant information.
    """

    # Make copy of toml data
    toml_data_new = copy.deepcopy(toml_data)
    if 'project' in toml_data:
        toml_data_new['project']['participant_height'] = height
        toml_data_new['project']['participant_mass'] = mass
        toml_data_new['kinematics']['default_height'] = height
    else:
        raise KeyError("The key 'project' is not present in the TOML data.")
    return toml_data_new

def update_framerate(toml_data, framerate):
    """
    Updates the frame rate in the TOML data with the given value.

    Args:
        toml_data (dict): The TOML data to be updated.
        framerate (float): The new frame rate to be set.

    Returns:
        dict: The updated TOML data with the new frame rate.
    """

    # Make copy of toml data
    toml_data_new = copy.deepcopy(toml_data)
    if 'project' in toml_data:
        toml_data_new['project']['frame_rate'] = framerate
    else:
        raise KeyError("The key 'project' is not present in the TOML data.")
    return toml_data_new

def select_and_save_frame_every_x(video_paths, output_dirs, x_interval=3):
    """
    Display every x-th frame from camera 1, let the user choose one to save,
    and save that same frame index from camera 2 and 3 into their respective output_dirs.

    Args:
        video_paths (list): List of paths to the three videos (cam1, cam2, cam3).
        output_dirs (list): List of output directories for the three cameras.
        x_interval (int): Interval to skip frames (every x-th frame is displayed).
    """
    assert len(video_paths) == 3 and len(output_dirs) == 3, "Expect 3 videos and 3 output folders (cam1, cam2, cam3)"
    
    # Read first video to choose frame
    cap1 = cv2.VideoCapture(video_paths[0])
    frame_index = 0
    selected_index = None

    while True:
        ret = cap1.grab()  # Skip non-x frames
        if not ret:
            break

        if frame_index % x_interval == 0:
            ret, frame = cap1.retrieve()
            if not ret:
                break

            cv2.imshow('Select Frame - Press s to save, n to skip', frame)
            key = cv2.waitKey(0)

            if key == ord('s'):
                selected_index = frame_index
                print(f'Selected frame: {selected_index}')
                break
            elif key == ord('q'):
                print('Quitting selection')
                break

        frame_index += 1

    cap1.release()
    cv2.destroyAllWindows()

    if selected_index is not None:
        for cam_index, (video_path, output_dir) in enumerate(zip(video_paths, output_dirs)):
            cap = cv2.VideoCapture(video_path)
            cap.set(cv2.CAP_PROP_POS_FRAMES, selected_index)
            ret, frame = cap.read()
            if ret:
                os.makedirs(output_dir, exist_ok=True)
                save_path = os.path.join(output_dir, f'selected_frame.png')
                cv2.imwrite(save_path, frame)
                print(f'Saved frame {selected_index} from cam{cam_index+1} to {save_path}')
            else:
                print(f'Could not read frame {selected_index} from {video_path}')
            cap.release()


def log_step_error(step_name, session_id, folder, error_log_path):
    """
    Logs an error that occurred during a processing step, including a human-readable message and full traceback.

    Args:
        step_name (str): Name of the processing step that failed.
        session_id (str): Identifier for the session where the error occurred.
        folder (str): Path to the folder associated with the session.
        error_log_path (str): Path to the file where the error details should be logged.

    The function performs the following actions:
    1. Creates error message combining the step name, session ID, and folder path
    2. Prints the message to the console
    3. Logs the message using Python's logging module
    4. Appends the message and full stack trace to the specified error log file
    """
    
    # Create a human-readable message
    msg = f'{step_name} failed for session {session_id} in folder {folder}'
    print(msg)
    logging.error(msg)

    # Write message + full traceback into a file
    with open(error_log_path, 'a') as f:
        f.write(msg + '\n')
        traceback.print_exc(file=f)  # full stack trace

Pose2sim triangulation

The Pose2sim pipeline comes in three steps: - calibration - triangulation - filtering

In calibration, we will use the calibration videos with checkerboard to calibrate the intrinsic and extrinsic parameters of the cameras. Note that we calibrate intrinsic parameters only once, and copy the file to the rest of the sessions. Extrinsic parameters are calibrated for each session (in part 1, and copied into part 2)

Fig. 1: Points of detection (determined by number of inner corners on the checkerboard, see Config.toml)

As noted in the Pose2sim documentation, intrinsic error should be below 0.5 pixels, and extrinsic error should be below 1 cm (but acceptable until 2.5 cm)

Note that extrinsic are sometimes not automatically detected so the corners need to be annotated manually.

Fig. 2: Automatically detected checkerboard corners

In triangulation, we will use the keypoints extracted in the previous script to triangulate the 3D position of the keypoints. The output will be a 3D position for each keypoint in each frame.

In filtering, we will filter the 3D position of the keypoints to remove noise and outliers with the in-build Butterworth filter (order 4, cut-off frequency 10 Hz).

Refer to the Config.toml file in Pose2Sim/FLESH_setup gfolder for the configuration of the pipeline.

There are additional three steps available in Pose2sim that we will not utilize in this script - synchronization, person association, and marker augmentation.

Copy intrinsic calibration file

# Find the toml file
toml_file = glob.glob(os.path.join(curfolder, 'projectdata', 'Session_1_1', 'calibration', '*.toml'), recursive=True)

# Copy the toml file to all folders
sessionfolders = glob.glob(os.path.join(curfolder, 'projectdata', '*'))

for session in sessionfolders:
    # Copy it to calibration folder
    calibration_folder = os.path.join(session, 'calibration')

    if not os.path.exists(calibration_folder):
        os.makedirs(calibration_folder)

    # Copy if there is no toml file in the calibration folder
    if not os.path.exists(os.path.join(calibration_folder, 'Calib_board.toml')):
        print('Intrinsic Calibration file not found')
        shutil.copy(toml_file, calibration_folder)
    else:
        print('Intrinsic Calibration file already exists')

Copy Config files and opensim folder

import copy

# Copy a folder in pose2simprjfolder and its contents to folders
source1 = os.path.join(pose2simprjfolder, 'Config.toml')
source2 = os.path.join(pose2simprjfolder, 'opensim')


for i in folderstotrack:
    print('working on folder: ' + i)

    sessionID = i.split(os.sep)[-1].split('_')[1]

    # First we need to prepare Config.file to all levels of folders
    # Copy to session folder
    shutil.copy(source1, i)
    toml_path = os.path.join(i, 'Config.toml')
    input_toml = load_toml(toml_path)

    # Update the p0 info
    mass_p0 = float(META.loc[(META['pcn_ID'] == sessionID + '_0'), 'weight'].values[0]) # this is new
    height_p0 = float(META.loc[(META['pcn_ID'] == sessionID + '_0'), 'height'].values[0])
    updated_toml_p0 = update_participant_info(copy.deepcopy(input_toml), height_p0, mass_p0)

    # Update p1 info
    mass_p1 = float(META.loc[(META['pcn_ID'] == sessionID + '_1'), 'weight'].values[0]) # this is new
    height_p1 = float(META.loc[(META['pcn_ID'] == sessionID + '_1'), 'height'].values[0])
    updated_toml_p1 = update_participant_info(copy.deepcopy(input_toml), height_p1, mass_p1)

    print('Participant p0 of session ' + sessionID + ', height = ' + str(height_p0) + ', mass = ' + str(mass_p0))
    print('Participant p1 of session ' + sessionID + ', height = ' + str(height_p1) + ', mass = ' + str(mass_p1))

    # Prepare all trial folders
    pcnfolders = glob.glob(os.path.join(i, '*'))
    pcnfolders = [x for x in pcnfolders if 'Config' not in x]
    pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
    pcnfolders = [x for x in pcnfolders if 'xml' not in x]
    pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
    pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
    pcnfolders = [x for x in pcnfolders if 'sto' not in x]
    pcnfolders = [x for x in pcnfolders if 'txt' not in x]
    pcnfolders = [x for x in pcnfolders if 'calibration' not in x]

    # Copy necessary files 
    for j in pcnfolders:
        print('working on subfolder: ' + j)

        # We also need to update framerate in the Config file such that it fits to each video
        # Find folder raw-2d and find avi file inside, read the fps
        raw2d_folder = os.path.join(j, 'raw-2d')
        if os.path.exists(raw2d_folder):
            avi_files = glob.glob(os.path.join(raw2d_folder, '*.avi'))
            if len(avi_files) > 0:
                # Read FPS from the first AVI file found
                cap = cv2.VideoCapture(avi_files[0])
                fps = cap.get(cv2.CAP_PROP_FPS)
                cap.release()

        # Update the framerate in the toml files
        updated_toml_p0 = update_framerate(copy.deepcopy(updated_toml_p0), fps)
        updated_toml_p1 = update_framerate(copy.deepcopy(updated_toml_p1), fps)

        if '_p0' in j:
            
            # Update toml and save in j 
            new_toml_p0_path = os.path.join(j, 'Config.toml')
            save_toml(updated_toml_p0, new_toml_p0_path)

        if '_p1' in j:
            new_toml_p1_path = os.path.join(j, 'Config.toml')
            save_toml(updated_toml_p1, new_toml_p1_path)

    opensim_path = os.path.join(i, 'opensim')

    if not os.path.exists(opensim_path):
        shutil.copytree(source2, opensim_path)
        print('source = ' + source2 + ' to destination: ' + os.path.join(i, 'opensim'))

Calibration, triangulation, filtering

for i in folderstotrack:
    os.chdir(i)
    print('working on folder: ' + i)

    sessionID = i.split(os.sep)[-1].split('_')[1]

    # Now we calibrate
    print('Step: Calibration')

    # Calibrate only if there is no toml file in the calibration folder
    if not os.path.exists(os.path.join(i, 'calibration', 'Object_points.trc')): 
        print('Extrensic Calibration file not found')
        
        # Now we prepare images from calibration videos
        calib_folders = glob.glob(os.path.join(i, 'calibration', '*', '*'))

        if not os.path.exists(os.path.join(i, 'calibration', 'extrinsics', 'cam1', 'selected_frame.png')):
            input_videos = [
                glob.glob(os.path.join(i, 'calibration', 'extrinsics', 'cam1', '*cam1.avi'))[0],
                glob.glob(os.path.join(i, 'calibration', 'extrinsics', 'cam2', '*cam2.avi'))[0],
                glob.glob(os.path.join(i, 'calibration', 'extrinsics', 'cam3', '*cam3.avi'))[0]
            ]

            output_dirs = [
                os.path.join(i, 'calibration', 'extrinsics', 'cam1'),
                os.path.join(i, 'calibration', 'extrinsics', 'cam2'),
                os.path.join(i, 'calibration', 'extrinsics', 'cam3')
            ]
            
            # Call the function to select and save frames
            select_and_save_frame_every_x(input_videos, output_dirs, x_interval=10)

        print('Calibration file does not exist, calibrating...')
        Pose2Sim.calibration() 
    
        # Get the last element of the i
        split = i.split(os.path.sep)
        parts = split[-1].split('_')

        # Get the sessionID
        session_id = parts[1]
        session_part = parts[-1]

        # If session_part is 1, we copy trc and calib file to the session that has some id, but part 2
        if session_part == '1':
            copy_to_part = '2'
            new_folder = 'Session_'+session_id+'_'+copy_to_part
            new_folder_path = os.path.join(inputfolder, new_folder)
            if not os.path.exists(os.path.join(new_folder_path, 'calibration')):
                os.makedirs(os.path.join(new_folder_path, 'calibration'))
            
            # Get the trc and toml files
            trc_file = os.path.join(i, 'calibration', 'Object_points.trc')
            toml_file = os.path.join(i, 'calibration', 'Calib_board.toml')
            
            # Copy the files to the new folder
            shutil.copy(trc_file, os.path.join(new_folder_path, 'calibration'))
            shutil.copy(toml_file, os.path.join(new_folder_path, 'calibration'))
        
        # Part 2 does not need to be calibrated so we can just proceed
        else:
            continue

    # If calibration file exists, then we can skip calibration
    else:
        print('Calibration file found, no need to calibrate')

    try:
        print('Step: triangulation')
        Pose2Sim.triangulation()
    except:
        print('Triangulation failed')
        log_step_error('Triangulation', sessionID, i, error_log)
        continue

    try:
        print('Step: filtering')
        Pose2Sim.filtering()
    except:
        print('Filtering failed')
        log_step_error('Filtering', sessionID, i, error_log)
        continue

Check if all folders have 3d-pose folder

# Get all the folders per session, per participant
pcnfolders = []

for i in folderstotrack:
    pcnfolders_in_session = glob.glob(os.path.join(i, '*'))   
    pcnfolders = pcnfolders + pcnfolders_in_session 

# Get rid of all pontetially confusing files/folders
pcnfolders = [x for x in pcnfolders if 'Config' not in x]
pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
pcnfolders = [x for x in pcnfolders if 'xml' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
pcnfolders = [x for x in pcnfolders if 'sto' not in x]
pcnfolders = [x for x in pcnfolders if 'txt' not in x]
pcnfolders = [x for x in pcnfolders if 'calibration' not in x]
print(pcnfolders[0:10])

print('length of pcnfolders: ' + str(len(pcnfolders)))

# Check if there is folder pose-3d that contains 2 trc files (original + filtered)
for j in pcnfolders:
    pose3dfolder = os.path.join(j, 'pose-3d')
    trcfiles = glob.glob(os.path.join(pose3dfolder, '*.trc'))
    if len(trcfiles) < 2: # one filtered and one unfiltered
        print('Warning: Less than 2 TRC files found in ' + pose3dfolder + ' Found ' + str(len(trcfiles)) + ' TRC files.')

Inspecting measurement error

Now let’s get information about the intrinsic and extrinsic error for each trial. The instrinsic error is identical to the one in preregistration. The extrinsic error is saved in pose2sim_triangulation_error_log.txt. To get more comprehensive overview of the error for each file, we will extract the relevant information for each trial by acessing regular expressions in the logs.txt file.

import re

# Function to load text file
def load_log_file(filename):
    """Loads the contents of a text file and returns it as a string.

    Args:
        filename (str): Path to the file to be loaded.

    Returns:
        str: The contents of the file as a single string.
    """
    with open(filename, 'r') as file:
        text = file.read()
    return text

# This is the reprojection/triangulation error file
log = os.path.join(curfolder, 'errors', 'pose2sim_triangulation_error_final.txt')

# Load the log file text
log_data = load_log_file(log)

# Pattern to match the header line for each trial
header_pattern = r"Triangulation of 2D points for (\S+),\s*for all frames\."
headers = list(re.finditer(header_pattern, log_data))
trial_results = []

# Get the extrinsic error for each session
for idx, h in enumerate(headers):
    trial_name_raw = h.group(1)
    start = h.end()
    end = headers[idx+1].start() if idx+1 < len(headers) else len(log_data)
    chunk = log_data[start:end]

    # Pattern to match the summary line
    m = re.search(
        r"Mean reprojection error for all points on .*? is ([0-9]+(?:\.[0-9]+)?) px,?\s*which roughly corresponds to ([0-9]+(?:\.[0-9]+)?) mm",
        chunk,
        flags=re.DOTALL
    )

    if m:
        px = m.group(1)
        mm = m.group(2)
    else:
        # fallback: more permissive, no "all frames" requirement
        m2 = re.search(
            r"Mean reprojection error for all points on .*? is ([0-9]+(?:\.[0-9]+)?) px",
            chunk
        )
        if m2:
            px = m2.group(1)
            mmsearch = re.search(
                r"which roughly corresponds to ([0-9]+(?:\.[0-9]+)?) mm",
                chunk
            )
            mm = mmsearch.group(1) if mmsearch else ''
        else:
            px = ''
            mm = ''

    # normalize trial name: remove trailing punctuation or whitespace
    trial_name = re.sub(r"[^\w\-_]+$", "", trial_name_raw)
    trial_results.append({
        'trial_name': trial_name,
        'mean_reprojection_error_px': px,
        'mean_reprojection_error_mm': mm
    })

df_tr = pd.DataFrame(trial_results)
df_tr.head(15)

	trial_name	mean_reprojection_error_px	mean_reprojection_error_mm
0	10_1_11_p1	3.4	14.9
1	10_1_12_p1	3.8	16.7
2	10_1_13_p1	3.5	15.4
3	10_1_14_p1	3.5	15.4
4	10_1_15_p1	3.5	15.4
5	10_1_16_p1	3.7	16.3
6	10_1_17_p1	3.4	14.9
7	10_1_20_p0	3.9	17.1
8	10_1_21_p0	3.7	16.3
9	10_1_22_p0	3.8	16.7
10	10_1_23_p0	3.8	16.7
11	10_1_24_p0	3.5	15.4
12	10_1_25_p0	3.8	16.7
13	10_1_26_p0	3.3	14.5
14	10_1_29_p1	3.5	15.4

What is the mean reprojection error across trials in mm?

Mean reprojection error for all trials is  14.943116820164809  mm

Is there any trial that has an error above 2 cm?

There are trials with reprojection error larger than 20 mm

	trial_name	mean_reprojection_error_px	mean_reprojection_error_mm
21	10_1_35_p1	4.6	20.2
402	13_1_21_p0	5.1	20.8
679	15_1_48_p1	5.2	20.2
1073	18_1_48_p1	4.5	20.1
1103	18_2_31_p1	5.2	23.3
...	...	...	...
8015	7_2_96_p1	8.0	32.7
8016	7_2_97_p1	7.8	31.9
8017	7_2_98_p1	8.0	32.7
8018	7_2_99_p1	8.1	33.1
8019	7_2_9_p0	8.3	34.0

302 rows × 3 columns

Some trials have error above 2 cm but not above 2.3 cm. Since acceptable error is up to 2.5 cm, we will keep these trials as is. Session 7 has after repeated re-calibrations and triangulations still trials with error above 2.5 cm. We will therefore exclude session 7 from further analysis.

Final error documentation

Intrinsic calibration error

For intrinsic calibration, we use the calibration file from preregistration with an error of 0.24 pixels for each camera (recommended below 0.5 pixels)

Extrinsic calibration error

ex_log = os.path.join(curfolder, 'errors', 'pose2sim_extrinsic_error.txt')
text = Path(ex_log).read_text(encoding="utf-8")

# Pattern to find information about extrinsic error
pattern = re.compile(
    r"(?m)^\s*(\d+)[^\n]*\n\s*--> Residual.*?\[([^\]]+)\]\s*px,.*?\[([^\]]+)\]\s*mm",
    re.S
)

rows = []
for m in pattern.finditer(text):
    session = int(m.group(1))
    # Extract numbers robustly, ignoring stray characters
    px_vals = [float(x) for x in re.findall(r"[-+]?\d*\.?\d+", m.group(2))][:3]
    mm_vals = [float(x) for x in re.findall(r"[-+]?\d*\.?\d+", m.group(3))][:3]
    if len(px_vals) == 3 and len(mm_vals) == 3:
        rows.append({
            "session": session,
            "cam1_px": px_vals[0], "cam2_px": px_vals[1], "cam3_px": px_vals[2],
            "cam1_mm": mm_vals[0], "cam2_mm": mm_vals[1], "cam3_mm": mm_vals[2],
        })

# Show in a dataframe
df_ex = pd.DataFrame(rows).sort_values("session").reset_index(drop=True)
print(df_ex.head(20))

# Calculate mean across all cameras in px
px_list = df_ex[['cam1_px', 'cam2_px', 'cam3_px']].values.flatten().tolist()
mm_list = df_ex[['cam1_mm', 'cam2_mm', 'cam3_mm']].values.flatten().tolist()
mean_px = np.mean(px_list).round(3)
mean_mm = np.mean(mm_list).round(3)
print(f'Mean reprojection/extrinsic error across all cameras is {mean_px} pixels')
print(f'Mean reprojection/extrinsic error across all cameras is {mean_mm} mm')

    session  cam1_px  cam2_px  cam3_px  cam1_mm  cam2_mm  cam3_mm
0         1    2.014    1.583    2.273    9.293    6.099    9.092
1         2    0.819    2.099    0.776    3.699    7.852    3.041
2         3    3.011    2.980    3.064   12.898   10.593   11.604
3         4    2.288    1.349    2.382    9.603    4.600    8.552
4         5    2.712    2.358    2.493   11.820    8.372    9.281
5         6    2.900    2.283    2.271   12.096    7.869    8.448
6         7    1.527    4.085    1.455    6.248   13.599    5.179
7         8    3.437    2.486    3.220   14.702    8.692   11.978
8         9    3.388    1.931    2.345   14.494    6.727    8.659
9        10    3.255    2.389    2.441   14.312    8.660    9.359
10       11    3.436    2.768    3.034   14.421    9.249   10.565
11       12    2.752    3.101    2.204   11.209   10.073    7.634
12       13    1.079    2.214    2.555    4.401    7.237    9.007
13       14    3.663    2.063    4.938   15.095    6.937   17.849
14       15    2.195    3.364    2.728    8.516   10.762    9.527
15       16    3.778    2.650    4.217   15.749    9.141   15.551
16       17    4.103    3.589    3.201   17.436   12.434   11.715
17       18    2.800    2.697    3.516   12.519   10.004   13.834
18       19    5.123    3.848    4.440   21.116   12.673   15.405
19       20    2.952    2.676    3.630   12.870    9.509   13.600
Mean reprojection/extrinsic error across all cameras is 2.515 pixels
Mean reprojection/extrinsic error across all cameras is 9.552 mm

The only session that is relatively high (above 20 mm) is session 19. Given that it does not exceed the acceptable threshold of 25 mm and the triangulation also does not exceed this threshold (see below), we will keep this session for further analysis.

Triangulation error

# create a mean reprojection error per session

df_tr['session'] = df_tr['trial_name'].apply(lambda x: '_'.join(x.split('_')[:2]))
session_means = df_tr.groupby('session')['mean_reprojection_error_mm'].mean().reset_index()
session_means['mean_reprojection_error_mm'] = session_means['mean_reprojection_error_mm'].round(3)
session_means.head(5)

	session	mean_reprojection_error_mm
0	10_1	15.786
1	10_2	15.648
2	11_1	12.843
3	11_2	12.559
4	12_1	12.052

Inspecting through animation

def trc_to_csv(file, framerate=60):
    """
    Convert TRC motion capture data to a CSV format DataFrame.

    Args:
        file (str): Path to the TRC file to be converted.
        framerate (int, optional): Framerate of the motion capture data. Defaults to 60.

    Returns:
        pd.DataFrame: A DataFrame containing the motion capture data with separated x, y, z columns for each marker and a time column.
    """
    
    mocap_data = TRCData()
    mocap_data.load(os.path.abspath(file))
    num_frames = mocap_data['NumFrames']
    markernames = mocap_data['Markers'] # the marker names are not

    # Convert mocap_data to pandas dataframe
    mocap_data_df = pd.DataFrame(mocap_data, columns=mocap_data['Markers'])
    # Each value within the dataframe consists a list of x,y,z coordinates, we want to seperate these out so that each marker and dimension has its own column
    colnames = []
    for marker in markernames:
        colnames.append(marker + '_x')
        colnames.append(marker + '_y')
        colnames.append(marker + '_z')

    # Create a new DataFrame to store separated values
    new_df = pd.DataFrame()

    # Iterate through each column in the original DataFrame
    for column in mocap_data_df.columns:
        # Extract the x, y, z values from each cell
        xyz = mocap_data_df[column].tolist()
        # Create a new DataFrame with the values in the cell separated into their own columns
        xyz_df = pd.DataFrame(xyz, columns=[column + '_x', column + '_y', column + '_z'])
        # Add the new columns to the new DataFrame
        new_df = pd.concat([new_df, xyz_df], axis=1)

    # Add a new time column to the new dataframe assuming the framerate was 60 fps
    time = []
    ts = 0
    for i in range(0, int(num_frames)):
        ts = ts + 1/framerate
        time.append(ts)

    # Add the time column to the new dataframe
    new_df['Time'] = time

    return new_df

This is the video that it corresponds to

We can see that there are some issues with the tracking that we need to correct for:

• the orientation of the Z-axis is inverted - this is probably an artifact of calibration where we ticked corners from top to bottom. Because this is important for working with forces (as gravity needs to be correctly oriented), we will invert the Z-axis immediately in this script. We use code that is adapted from here
• as we see in the animation, Y and Z axis are flipped. Again, for correct estimation of forces, we need to correct for this.
• lastly, we can see that the skeleton does not align with the ground plane, i.e., the feet are not on the ground (Z=0). We will therefore also correct for this by shifting the whole skeleton up or down such that the lowest point (foot) is at Z=0. This is done using part of Pose2sim augmentation script
• we will further smooth the data using Savitzky-Golay filter to remove high-frequency noise and prevent its propagation into inverse kinematics and dynamics.

This is how resulting time series look like for a single trial (keypoints: wrist, knee)

Comparing it to the video above, we can see that the tracking is still a bit jittery despite the Butterworth filter. Since this noise is likely going to propagate when solving inverse kinematics using OpenSim, we will apply Savitzky-Golay filter to smooth the coordinates further.

Inversion, offset, smoothing plus csv-convert on all files

upperbodycols = [
    'Head', 'Neck', 'RShoulder', 'RElbow', 'RWrist',
    'LShoulder', 'LElbow', 'LWrist', 'Nose',
    'RIndex', 'LIndex'
]

lowerbodycols = [
    'RHip', 'RKnee', 'RAnkle', 'RHeel',
    'LHip', 'LKnee', 'LAnkle', 'LHeel'
]

def remove_global_pitch_safe(coords_out, marker_indices, max_correction_deg=30):
    """
    Corrects global pitch based on pelvis orientation only.
    Limits correction to +/- max_correction_deg to avoid flips.
    Applies correction ONLY based on the first frame.

    Args:
        coords_out (pd.DataFrame): DataFrame containing marker coordinates (columns: X1,Y1,Z1,...)
        marker_indices (dict): Mapping of marker names to their 1-based indices
        max_correction_deg (float): Maximum allowed correction in degrees (default: 30)
    """

    pelvis_markers = ['LHip', 'RHip']
    for pm in pelvis_markers:
        if pm not in marker_indices:
            print("[WARN] Missing pelvis markers; skipping pitch correction")
            return coords_out

    def get_first_frame(marker):
        idx = marker_indices[marker]
        base = 3 * (idx - 1)
        return coords_out.iloc[0, base:base+3].values

    # Pelvis center
    L = get_first_frame('LHip')
    R = get_first_frame('RHip')
    pelvis = (L + R) / 2

    # Estimate pitch using pelvis-to-shoulder vector IF shoulders exist
    if all(m in marker_indices for m in ['LShoulder', 'RShoulder']):
        LS = get_first_frame('LShoulder')
        RS = get_first_frame('RShoulder')
        shoulders = (LS + RS) / 2
        vec = shoulders - pelvis
    else:
        # fallback: use hip-to-knee
        LK = get_first_frame('LKnee')
        RK = get_first_frame('RKnee')
        vec = ((LK + RK) / 2) - pelvis

    dy = vec[1]
    dz = vec[2]

    raw_pitch = np.degrees(np.arctan2(dz, dy))

    # Clamp correction
    correction = np.clip(raw_pitch, -max_correction_deg, max_correction_deg)

    rad = np.radians(-correction)

    cos_t = np.cos(rad)
    sin_t = np.sin(rad)

    # Rotate all frames around X-axis
    for marker, idx in marker_indices.items():
        base = 3 * (idx - 1)

        Y = coords_out.iloc[:, base + 1].values
        Z = coords_out.iloc[:, base + 2].values

        Y_new = cos_t * Y - sin_t * Z
        Z_new = sin_t * Y + cos_t * Z

        coords_out.iloc[:, base + 1] = Y_new
        coords_out.iloc[:, base + 2] = Z_new

    print(f"[OK] Applied pitch correction of {correction:.2f} degrees")
    return coords_out

def remove_global_roll_safe(coords_out, marker_indices, max_correction_deg=15):
    """
    Corrects global roll (left/right tilt) using pelvis markers (LHip, RHip).
    Applies rotation around Z-axis.
    Limits correction to +/- max_correction_deg.

    Args:
            coords_out (pd.DataFrame): DataFrame containing marker coordinates (columns: X1,Y1,Z1,...)
            marker_indices (dict): Mapping of marker names to their 1-based indices
            max_correction_deg (float): Maximum allowed correction in degrees (default: 15)

    """

    if not all(m in marker_indices for m in ['LHip', 'RHip']):
        print("[WARN] Missing pelvis markers; skipping roll correction")
        return coords_out

    def ff(marker):
        idx = marker_indices[marker]
        base = 3 * (idx - 1)
        return coords_out.iloc[0, base:base+3].values

    L = ff('LHip')
    R = ff('RHip')

    # Vertical difference between hips
    dy = L[1] - R[1]
    dx = L[0] - R[0]

    # Roll angle (hip height difference)
    raw_roll = np.degrees(np.arctan2(dy, dx))

    # Clamp
    correction = np.clip(raw_roll, -max_correction_deg, max_correction_deg)

    rad = np.radians(+correction)
    cos_t = np.cos(rad)
    sin_t = np.sin(rad)

    # Rotate all points around Z axis
    for marker, idx in marker_indices.items():
        base = 3 * (idx - 1)

        X = coords_out.iloc[:, base + 0].values
        Y = coords_out.iloc[:, base + 1].values

        X_new = cos_t * X - sin_t * Y
        Y_new = sin_t * X + cos_t * Y

        coords_out.iloc[:, base + 0] = X_new
        coords_out.iloc[:, base + 1] = Y_new

    print(f"[OK] Applied roll correction of {correction:.2f} degrees")
    return coords_out

def process_trc(trc_path):
    """
    Pipeline for a single TRC:

      1. load TRC
      2. apply axis transform per marker:
           X_new = X_old
           Y_new = -Z_old
           Z_new =  Y_old
      3. smooth (upper body: 11/3, lower body: 25/3, with fallbacks)
      4. shift vertically so feet are on the floor (min Y of feet ≈ 0)
      5. save ONE final TRC: <name>_processed.trc
      6. save ONE CSV with columns Marker_x, Marker_y, Marker_z, ..., Time

    Returns:
        final_trc_path, final_csv_path
    """

    base, ext = os.path.splitext(trc_path)
    final_trc_path = base + "_processed.trc"
    final_csv_path = base + "_processed.csv"

    # Read header and marker names
    with open(trc_path, "r") as f:
        header = [next(f) for _ in range(5)]

    marker_header = header[3].rstrip("\n")
    parts = marker_header.split("\t")

    # marker_name -> index (1-based, matching X1/Y1/Z1,...)
    marker_indices = {}
    marker_num = 0
    for name in parts[2:]:   # skip "Frame#" and "Time"
        if name != "":
            marker_num += 1
            marker_indices[name] = marker_num

    # Read numeric block
    df = pd.read_csv(trc_path, sep="\t", skiprows=4, encoding="utf-8")

    # df columns: ['Frame#','Time','X1','Y1','Z1',..., 'Xn','Yn','Zn']
    frames = df.iloc[:, 0]
    time_col = df.iloc[:, 1]
    coords = df.drop(df.columns[[0, 1]], axis=1)

    # this will hold transformed (invertZ + swap Y/Z) and smoothed coords
    coords_out = coords.copy()

    # Axis transform
    for marker_name, marker_idx in marker_indices.items():
        base_idx = 3 * (marker_idx - 1)
        x_col = coords.columns[base_idx]
        y_col = coords.columns[base_idx + 1]
        z_col = coords.columns[base_idx + 2]

        X_old = coords[x_col]
        Y_old = coords[y_col]
        Z_old = coords[z_col]

        coords_out[x_col] = X_old           # X_new
        coords_out[y_col] = -Z_old          # Y_new = -Z_old (invert original Z)
        coords_out[z_col] = Y_old           # Z_new =  Y_old

    # Smooth
    for marker_name, marker_idx in marker_indices.items():

        # choose smoothing parameters
        if marker_name in upperbodycols:
            primary = (11, 3)
        elif marker_name in lowerbodycols:
            primary = (25, 3)
        else:
            primary = (11, 3)

        fallbacks = [(10, 2), (5, 2)]

        base_idx = 3 * (marker_idx - 1)
        col_names = [
            coords_out.columns[base_idx],
            coords_out.columns[base_idx + 1],
            coords_out.columns[base_idx + 2],
        ]

        for col_name in col_names:
            signal = coords_out[col_name].values
            applied = False

            for win, ord_ in (primary,) + tuple(fallbacks):
                try:
                    coords_out[col_name] = scipy.signal.savgol_filter(
                        signal, win, ord_
                    )
                    applied = True
                    break
                except ValueError:
                    continue

            if not applied:
                coords_out[col_name] = signal

    # Feet on floor
    ground_markers = ['RHeel', 'LHeel', 'RAnkle', 'LAnkle']

    ground_y_cols = []
    for m in ground_markers:
        if m in marker_indices:
            idx = marker_indices[m]
            base_idx = 3 * (idx - 1)
            y_col = coords_out.columns[base_idx + 1]
            ground_y_cols.append(y_col)

    if ground_y_cols:
        min_y = coords_out[ground_y_cols].min().min()

        for marker_name, marker_idx in marker_indices.items():
            base_idx = 3 * (marker_idx - 1)
            y_col = coords_out.columns[base_idx + 1]
            coords_out[y_col] = coords_out[y_col] - min_y

    if 'Session_52_' in trc_path:
        print("[INFO] Applying special global roll correction for trial starting with '52_'")
        coords_out = remove_global_roll_safe(coords_out, marker_indices)

    coords_out = remove_global_pitch_safe(coords_out, marker_indices)

    # Save final TRC
    with open(final_trc_path, "w") as out:
        for line in header:
            out.write(line)

        df_out = coords_out.copy()
        df_out.insert(0, "Time", time_col)
        df_out.insert(0, "Frame#", frames)

        df_out.to_csv(out, sep="\t", index=False, header=None, lineterminator="\n")

    print(f"[OK] Final TRC saved → {final_trc_path}")

    # Save CSV using real timestamps from TRC (not synthetic)
    new_df = pd.DataFrame()

    for marker_name, marker_idx in marker_indices.items():
        base_idx = 3 * (marker_idx - 1)
        x_col = coords_out.columns[base_idx]
        y_col = coords_out.columns[base_idx + 1]
        z_col = coords_out.columns[base_idx + 2]

        new_df[f"{marker_name}_x"] = coords_out[x_col].values
        new_df[f"{marker_name}_y"] = coords_out[y_col].values
        new_df[f"{marker_name}_z"] = coords_out[z_col].values

    new_df["Time"] = time_col.values  # real timestamps from TRC

    new_df.to_csv(final_csv_path, index=False)
    print(f"[OK] Final CSV saved → {final_csv_path}")

    return final_trc_path, final_csv_path

trctoconvert = []

for j in pcnfolders:
    # Here we store the 3D pose data
    posefolder = 'pose-3d'
    trcfiles = glob.glob(os.path.join(j, posefolder, '*.trc'))
    
    # keep only butterworth smoothed files
    trcfiles = [file for file in trcfiles if 'butterworth' in file]
    trcfiles = [file for file in trcfiles if 'processed.trc' not in file] # get rid of the processed ones already
    
    trctoconvert = trctoconvert + trcfiles

# Loop through files and convert to csv
for file in trctoconvert:
    print(file)
    
    # Process each TRC through the full pipeline
    final_trc, final_csv = process_trc(file)

Check the animation again

Total processed TRC files: 16504

Plot as before

Inverse kinematics with Pose2sim

Currently, Pose2sim implements two steps of the OpenSim pipeline - scaling and inverse kinematics.

In scaling, we scale the model to match the anthropometry of the subject. We use the Pose2sim model with 135 keypoints (BODY_135). While OpenSim requires t-pose for scaling the participants, Pose2Sim uses every trial’s coordinates to estimate the participant’s anthropometry. We end up with a model scaled for each participant.

In inverse kinematics, we use the scaled model to estimate the joint angles of the participant. We use the motion tracking data to estimate the joint angles. Joint angles are saved as .mot files.

# Prepare special log txt
error_log = os.path.join(curfolder, 'pose2sim_error_IK_log.txt')
if not os.path.exists(error_log):
    with open(error_log, 'w') as f:
        f.write('Pose2Sim processing error log\n')
        f.write('============================\n\n')

for i in folderstotrack:
    os.chdir(i)
    print('working on folder: ' + i)

    sessionID = i.split(os.sep)[-1].split('_')[1]

    try:
        print('Step: inverse kinematics')
        Pose2Sim.kinematics()

    except:
        print('IK failed')
        log_step_error('Kinematics', sessionID, i, error_log)
        continue

First, check that every folder has kinematics with osim and mot files

# Get all the folders per session, per participant
pcnfolders = []

for i in folderstotrack:
    pcnfolders_in_session = glob.glob(os.path.join(i, '*'))   
    pcnfolders = pcnfolders + pcnfolders_in_session 

# Get rid of all pontetially confusing files/folders
pcnfolders = [x for x in pcnfolders if 'Config' not in x]
pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
pcnfolders = [x for x in pcnfolders if 'xml' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
pcnfolders = [x for x in pcnfolders if 'sto' not in x]
pcnfolders = [x for x in pcnfolders if 'txt' not in x]
pcnfolders = [x for x in pcnfolders if 'calibration' not in x]
print(pcnfolders[0:10])


print('Length of pcnfolders: ' + str(len(pcnfolders)))

for j in pcnfolders:
    ikfolder = os.path.join(j, 'kinematics')
    osimfiles = glob.glob(os.path.join(ikfolder, '*.osim'))
    if len(osimfiles) == 0:
        print('Warning: No osim files found in ' + ikfolder)
    motfiles = glob.glob(os.path.join(ikfolder, '*.mot'))
    if len(motfiles) == 0:
        print('Warning: No mot files found in ' + ikfolder)

['f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_11_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_12_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_13_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_14_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_15_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_16_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_17_p1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_20_p0', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_21_p0', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1\\10_1_22_p0']
Length of pcnfolders: 8252

Then, for each file get marker_error_RMS and max_error from the errors.sto files

ikerrors = []

for j in pcnfolders:

    # Skip poor tracking sessions
    if os.path.basename(j).startswith('7_'):
        continue
    
    ikfolder = os.path.join(j, 'kinematics')
    errorfiles = glob.glob(os.path.join(ikfolder, '*_ik_marker_errors.sto'))
    if len(errorfiles) == 0:
        print('Warning: No ik error files found in ' + ikfolder)

    error = pd.read_csv(errorfiles[0], sep="\t", skiprows=6)
    # get mean and sd of marker_error_RMS
    mean_error = error['marker_error_RMS'].mean()
    sd_error = error['marker_error_RMS'].std()

    # get mean and sd of marker_error_max
    max_error = error['marker_error_max'].mean()
    sd_max_error = error['marker_error_max'].std()

    filename = os.path.basename(j)

    ikerrors.append({
        'pcnfolder': filename,
        'mean_marker_error_RMS': mean_error,
        'sd_marker_error_RMS': sd_error,
        'mean_marker_error_max': max_error,
        'sd_marker_error_max': sd_max_error
    })

# create dataframe
df_ikerrors = pd.DataFrame(ikerrors)

# round all to three decimal places
df_ikerrors = df_ikerrors.round(3)
df_ikerrors.head(15)

# save it
df_ikerrors.to_csv(os.path.join('errors', 'ik_errors_all.csv'), index=False)

Top 10 worst remaining trials:
    pcnfolder  mean_marker_error_RMS  mean_marker_error_max
0  40_1_17_p1                  0.115                  0.339
1  27_2_39_p0                  0.099                  0.236
2  64_2_29_p1                  0.098                  0.243
3  26_2_69_p1                  0.092                  0.209
4  26_2_94_p0                  0.090                  0.206
5  70_1_24_p0                  0.090                  0.202
6   8_2_27_p1                  0.090                  0.187
7  10_1_20_p0                  0.087                  0.218
8   8_2_29_p1                  0.085                  0.173
9  35_2_41_p0                  0.085                  0.218

861 trials above 4cm threshold (10.6%)
['10_1_12_p1', '10_1_15_p1', '10_1_20_p0', '10_1_22_p0', '10_1_23_p0', '10_1_35_p1', '10_1_49_p1', '10_1_50_p1', '10_1_51_p1', '10_1_6_p0', '10_2_11_p0', '10_2_16_p0', '10_2_31_p1', '10_2_32_p1', '10_2_34_p1', '10_2_39_p0', '10_2_41_p0', '10_2_42_p0', '10_2_46_p0', '10_2_47_p0', '10_2_50_p1', '10_2_57_p1', '10_2_5_p0', '10_2_6_p0', '10_2_91_p1', '10_2_92_p1', '10_2_97_p1', '10_2_98_p1', '11_1_3_p0', '11_1_6_p0', '11_1_7_p0', '11_1_8_p0', '11_2_100_p1', '11_2_47_p0', '11_2_48_p0', '11_2_86_p0', '12_1_16_p1', '12_2_12_p0', '12_2_14_p0', '12_2_15_p0', '12_2_27_p1', '12_2_31_p1', '12_2_32_p1', '12_2_33_p1', '12_2_88_p0', '12_2_89_p0', '12_2_97_p0', '13_1_23_p0', '13_1_34_p1', '13_1_38_p0', '13_1_44_p0', '13_1_49_p1', '13_1_50_p1', '13_1_51_p1', '13_1_52_p1', '13_1_53_p1', '13_2_106_p1', '13_2_107_p1', '13_2_53_p0', '13_2_56_p0', '13_2_66_p1', '13_2_72_p1', '13_2_73_p1', '13_2_74_p1', '13_2_90_p0', '13_2_92_p0', '13_2_97_p1', '14_1_39_p0', '14_1_49_p1', '14_1_5_p0', '14_2_11_p0', '14_2_12_p0', '14_2_6_p0', '15_1_12_p1', '15_1_15_p1', '15_1_17_p1', '15_1_47_p1', '15_1_48_p1', '15_2_100_p1', '15_2_101_p1', '15_2_103_p1', '15_2_104_p1', '15_2_10_p0', '15_2_20_p1', '15_2_22_p1', '15_2_23_p1', '15_2_24_p1', '15_2_26_p1', '15_2_27_p1', '15_2_28_p1', '15_2_30_p1', '15_2_47_p0', '15_2_62_p1', '15_2_63_p1', '15_2_6_p0', '15_2_79_p0', '15_2_96_p1', '15_2_98_p1', '15_2_99_p1', '16_2_19_p1', '16_2_24_p1', '16_2_25_p1', '16_2_45_p0', '16_2_56_p1', '16_2_5_p0', '16_2_61_p1', '16_2_62_p1', '16_2_63_p1', '17_1_39_p0', '17_1_42_p0', '17_1_44_p0', '17_1_6_p0', '17_2_94_p0', '17_2_99_p0', '18_1_24_p0', '18_1_39_p0', '18_1_40_p0', '18_1_48_p1', '18_2_15_p0', '18_2_31_p1', '18_2_4_p0', '18_2_5_p0', '18_2_6_p0', '18_2_8_p0', '18_2_9_p0', '18_2_rein_44_p0', '18_2_rein_72_p1', '18_2_rein_80_p0', '18_2_rein_84_p0', '18_2_rein_85_p0', '18_2_rein_89_p0', '19_1_21_p0', '19_1_38_p0', '19_2_47_p0', '19_2_49_p0', '19_2_54_p0', '19_2_63_p1', '19_2_66_p1', '19_2_68_p1', '19_2_85_p0', '19_2_86_p0', '19_2_8_p0', '19_2_98_p1', '1_1_21_p0', '1_1_26_p0', '1_2_rein_100_p0', '1_2_rein_101_p0', '1_2_rein_102_p0', '1_2_rein_112_p1', '1_2_rein_121_p1', '1_2_rein_60_p0', '1_2_rein_61_p0', '1_2_rein_63_p0', '1_2_rein_64_p0', '1_2_rein_68_p0', '1_2_rein_79_p1', '1_2_rein_99_p0', '20_1_42_p0', '20_1_43_p0', '20_1_51_p1', '20_2_11_p0', '20_2_48_p0', '20_2_49_p0', '20_2_50_p0', '20_2_69_p1', '20_2_84_p0', '20_2_85_p0', '20_2_90_p0', '20_2_96_p1', '20_2_99_p1', '20_2_9_p0', '21_1_30_p1', '21_1_53_p1', '22_1_20_p0', '22_1_50_p1', '22_2_111_p1', '22_2_22_p1', '22_2_32_p1', '23_1_2_p0', '23_1_38_p0', '23_1_41_p0', '23_1_44_p0', '23_1_48_p1', '23_1_8_p0', '24_1_14_p1', '24_1_15_p1', '24_1_16_p1', '24_1_42_p0', '24_2_105_p0', '24_2_110_p1', '24_2_111_p1', '24_2_113_p1', '24_2_114_p1', '24_2_21_p0', '24_2_27_p1', '24_2_32_p1', '24_2_73_p1', '24_2_76_p1', '24_2_77_p1', '26_1_11_p1', '26_1_12_p1', '26_1_17_p1', '26_1_35_p1', '26_1_47_p1', '26_1_48_p1', '26_1_49_p1', '26_1_4_p0', '26_1_50_p1', '26_1_51_p1', '26_1_53_p1', '26_1_5_p0', '26_1_7_p0', '26_2_106_p1', '26_2_110_p1', '26_2_111_p1', '26_2_19_p0', '26_2_20_p0', '26_2_29_p1', '26_2_31_p1', '26_2_33_p1', '26_2_45_p0', '26_2_54_p0', '26_2_55_p0', '26_2_61_p0', '26_2_66_p1', '26_2_69_p1', '26_2_84_p0', '26_2_85_p0', '26_2_86_p0', '26_2_88_p0', '26_2_90_p0', '26_2_92_p0', '26_2_94_p0', '27_1_20_p0', '27_1_22_p0', '27_1_34_p1', '27_1_6_p0', '27_2_20_p1', '27_2_27_p1', '27_2_28_p1', '27_2_39_p0', '27_2_3_p0', '27_2_52_p1', '27_2_53_p1', '27_2_55_p1', '27_2_71_p0', '27_2_80_p0', '27_2_81_p0', '28_1_25_p0', '28_1_33_p1', '28_1_39_p0', '28_1_42_p0', '28_2_19_p1', '28_2_34_p0', '28_2_36_p0', '28_2_78_p0', '28_2_79_p0', '28_2_80_p0', '28_2_92_p1', '29_1_25_p0', '29_1_29_p1', '29_1_40_p0', '29_1_49_p1', '29_2_19_p1', '29_2_20_p1', '29_2_75_p0', '29_2_84_p0', '29_2_88_p1', '2_1_12_p1', '2_1_13_p1', '2_1_14_p1', '2_1_16_p1', '2_2_102_p1', '2_2_103_p1', '2_2_104_p1', '2_2_25_p1', '2_2_30_p1', '2_2_95_p1', '2_2_97_p1', '2_2_98_p1', '2_2_99_p1', '30_1_21_p0', '30_1_25_p0', '30_1_33_p1', '30_1_39_p0', '30_2_103_p1', '30_2_46_p0', '30_2_47_p0', '30_2_63_p1', '30_2_65_p1', '30_2_74_p1', '30_2_76_p1', '30_2_77_p1', '30_2_78_p1', '30_2_91_p0', '30_2_97_p1', '30_2_98_p1', '31_1_14_p1', '31_1_20_p0', '31_1_42_p0', '31_1_44_p0', '31_2_rein_100_p1', '31_2_rein_107_p1', '31_2_rein_32_p1', '32_2_30_p1', '32_2_58_p1', '32_2_59_p1', '32_2_61_p1', '33_1_35_p1', '33_1_52_p1', '33_2_106_p1', '33_2_107_p1', '33_2_11_p0', '33_2_14_p0', '33_2_15_p0', '33_2_18_p0', '33_2_23_p1', '33_2_24_p1', '33_2_26_p1', '33_2_63_p1', '33_2_66_p1', '33_2_80_p0', '33_2_81_p0', '33_2_82_p0', '33_2_86_p0', '33_2_8_p0', '33_2_9_p0', '34_1_48_p1', '34_1_50_p1', '34_2_50_p0', '34_2_52_p0', '34_2_55_p0', '34_2_70_p1', '34_2_73_p0', '34_2_74_p0', '34_2_7_p0', '34_2_92_p1', '35_2_16_p0', '35_2_22_p1', '35_2_23_p1', '35_2_32_p1', '35_2_41_p0', '35_2_44_p0', '35_2_46_p0', '35_2_54_p1', '35_2_62_p1', '36_1_5_p0', '36_2_13_p0', '36_2_rein_50_p0', '36_2_rein_62_p0', '36_2_rein_65_p0', '36_2_rein_90_p0', '36_2_rein_91_p0', '36_2_rein_93_p0', '36_2_rein_94_p0', '36_2_rein_99_p0', '37_1_12_p1', '37_1_14_p1', '37_1_17_p1', '37_1_39_p0', '37_1_3_p0', '37_1_43_p0', '37_2_106_p1', '37_2_15_p0', '37_2_19_p0', '37_2_25_p1', '37_2_26_p1', '37_2_30_p1', '37_2_31_p1', '37_2_6_p0', '37_2_82_p1', '37_2_91_p0', '37_2_92_p0', '37_2_9_p0', '38_1_12_p1', '38_2_12_p0', '38_2_23_p0', '38_2_46_p0', '38_2_54_p0', '38_2_62_p1', '38_2_66_p1', '39_1_32_p1', '39_2_102_p1', '39_2_103_p1', '39_2_107_p1', '39_2_113_p1', '39_2_71_p1', '39_2_72_p1', '39_2_73_p1', '39_2_76_p1', '39_2_81_p1', '39_2_88_p0', '39_2_97_p0', '3_1_25_p0', '3_1_33_p1', '3_1_38_p0', '3_1_39_p0', '3_1_40_p0', '3_1_42_p0', '3_1_43_p0', '3_1_52_p1', '3_2_14_p0', '3_2_43_p0', '3_2_44_p0', '3_2_45_p0', '3_2_46_p0', '3_2_49_p0', '3_2_53_p0', '3_2_58_p1', '3_2_65_p1', '3_2_67_p1', '3_2_71_p0', '3_2_73_p0', '3_2_74_p0', '3_2_75_p0', '3_2_76_p0', '3_2_79_p0', '3_2_80_p0', '3_2_92_p1', '3_2_93_p1', '3_2_9_p0', '40_1_17_p1', '40_1_23_p0', '40_1_29_p1', '40_1_3_p0', '40_2_103_p1', '40_2_105_p1', '40_2_12_p0', '40_2_13_p0', '40_2_15_p0', '40_2_34_p1', '40_2_37_p1', '40_2_46_p0', '40_2_54_p0', '40_2_88_p0', '41_1_25_p0', '41_2_6_p0', '41_2_89_p1', '41_2_90_p1', '45_1_11_p1', '45_1_13_p1', '45_1_3_p0', '45_2_106_p1', '45_2_111_p1', '45_2_51_p0', '45_2_52_p0', '45_2_56_p0', '46_1_11_p1', '46_1_15_p1', '46_1_16_p1', '46_1_25_p0', '46_1_29_p1', '46_1_30_p1', '46_1_33_p1', '46_1_49_p1', '46_1_4_p0', '46_1_50_p1', '46_2_101_p1', '46_2_16_p1', '46_2_17_p1', '46_2_25_p1', '46_2_60_p1', '46_2_82_p0', '46_2_86_p0', '46_2_90_p1', '46_2_91_p1', '46_2_92_p1', '46_2_94_p1', '47_1_23_p0', '47_1_31_p1', '47_1_40_p0', '47_1_42_p0', '47_1_51_p1', '47_1_52_p1', '48_1_26_p0', '48_1_31_p1', '48_2_37_p0', '48_2_40_p0', '48_2_43_p0', '48_2_44_p0', '48_2_45_p0', '48_2_4_p0', '48_2_53_p1', '48_2_58_p1', '48_2_8_p0', '49_2_rein_74_p0', '49_2_rein_76_p0', '49_2_rein_94_p1', '50_1_rein_49_p1', '50_2_105_p1', '50_2_109_p1', '50_2_45_p0', '50_2_47_p0', '50_2_48_p0', '50_2_54_p0', '50_2_71_p1', '50_2_72_p1', '50_2_76_p1', '50_2_83_p0', '50_2_89_p0', '51_1_22_p0', '51_1_24_p0', '51_1_33_p1', '51_1_39_p0', '51_1_43_p0', '51_1_50_p1', '51_1_52_p1', '51_2_48_p0', '51_2_50_p0', '51_2_53_p0', '51_2_57_p0', '51_2_62_p1', '51_2_65_p1', '51_2_68_p1', '51_2_82_p0', '51_2_83_p0', '51_2_86_p0', '51_2_87_p0', '51_2_93_p1', '51_2_95_p1', '51_2_99_p1', '52_1_15_p1', '52_1_21_p0', '52_1_23_p0', '52_1_2_p0', '52_1_35_p1', '52_2_4_p0', '52_2_5_p0', '52_2_69_p0', '52_2_70_p0', '52_2_75_p0', '52_2_76_p0', '52_2_84_p1', '52_2_85_p1', '52_2_8_p0', '52_2_91_p1', '52_2_95_p1', '53_1_14_p1', '53_2_12_p0', '53_2_16_p0', '53_2_22_p1', '53_2_30_p1', '53_2_80_p0', '53_2_8_p0', '54_1_13_p1', '54_1_14_p1', '54_1_38_p0', '54_1_39_p0', '54_1_40_p0', '54_1_44_p0', '54_1_47_p1', '54_1_48_p1', '54_1_49_p1', '54_1_52_p1', '54_1_8_p0', '54_2_10_p0', '54_2_11_p0', '54_2_15_p1', '54_2_17_p1', '54_2_18_p1', '54_2_20_p1', '54_2_22_p1', '54_2_23_p1', '54_2_24_p1', '54_2_25_p1', '54_2_26_p1', '54_2_27_p1', '54_2_28_p1', '54_2_4_p0', '54_2_9_p0', '54_2_rein_102_p1', '54_2_rein_103_p1', '54_2_rein_48_p1', '54_2_rein_50_p1', '54_2_rein_52_p1', '54_2_rein_53_p1', '54_2_rein_55_p1', '54_2_rein_57_p1', '54_2_rein_81_p0', '54_2_rein_92_p1', '54_2_rein_94_p1', '55_2_108_p1', '55_2_109_p1', '55_2_21_p1', '55_2_23_p1', '55_2_25_p1', '55_2_29_p1', '55_2_2_p0', '55_2_32_p1', '55_2_88_p0', '56_1_34_p1', '56_1_4_p0', '56_2_107_p1', '56_2_15_p0', '56_2_63_p0', '56_2_68_p1', '56_2_69_p1', '56_2_73_p1', '56_2_74_p1', '56_2_77_p1', '56_2_83_p0', '56_2_88_p0', '56_2_90_p0', '57_1_11_p1', '57_2_100_p1', '57_2_101_p1', '57_2_6_p0', '57_2_91_p0', '57_2_92_p0', '57_2_93_p0', '58_1_12_p1', '58_1_17_p1', '58_1_44_p0', '58_2_23_p1', '58_2_24_p1', '58_2_29_p1', '58_2_30_p1', '58_2_31_p1', '58_2_39_p0', '58_2_45_p0', '58_2_60_p1', '58_2_7_p0', '58_2_88_p0', '59_2_80_p0', '59_2_91_p1', '5_1_11_p1', '5_1_13_p1', '5_1_14_p1', '5_1_17_p1', '5_1_47_p1', '5_1_48_p1', '5_1_51_p1', '5_1_52_p1', '5_1_53_p1', '5_2_100_p1', '5_2_13_p0', '5_2_20_p1', '5_2_21_p1', '5_2_23_p1', '5_2_24_p1', '5_2_27_p1', '5_2_29_p1', '5_2_30_p1', '5_2_31_p1', '5_2_32_p1', '5_2_33_p1', '5_2_34_p1', '5_2_54_p1', '5_2_56_p1', '5_2_57_p1', '5_2_58_p1', '5_2_59_p1', '5_2_61_p1', '5_2_62_p1', '5_2_63_p1', '5_2_69_p0', '60_1_20_p0', '60_1_23_p0', '60_1_26_p0', '60_2_12_p0', '60_2_27_p1', '60_2_28_p1', '60_2_5_p0', '60_2_78_p0', '60_2_7_p0', '60_2_80_p0', '60_2_82_p0', '60_2_85_p0', '60_2_89_p0', '60_2_94_p0', '60_2_9_p0', '61_1_20_p0', '61_1_22_p0', '61_1_23_p0', '61_1_24_p0', '61_1_25_p0', '61_1_26_p0', '61_1_29_p1', '61_2_102_p1', '61_2_12_p0', '61_2_13_p0', '61_2_14_p0', '61_2_24_p1', '61_2_30_p1', '61_2_33_p1', '61_2_43_p0', '61_2_44_p0', '61_2_45_p0', '61_2_48_p0', '61_2_51_p0', '61_2_5_p0', '61_2_64_p1', '61_2_6_p0', '61_2_78_p0', '61_2_79_p0', '62_1_16_p1', '62_1_20_p0', '62_1_2_p0', '62_1_3_p0', '62_1_7_p0', '62_2_15_p0', '62_2_23_p1', '62_2_24_p1', '62_2_27_p1', '62_2_33_p1', '62_2_3_p0', '62_2_44_p0', '62_2_48_p0', '62_2_5_p0', '62_2_61_p1', '62_2_67_p1', '62_2_79_p0', '62_2_7_p0', '62_2_81_p0', '63_1_15_p1', '63_1_20_p0', '63_1_21_p0', '63_1_2_p0', '63_1_31_p1', '63_1_34_p1', '63_1_4_p0', '63_1_51_p1', '63_1_6_p0', '63_2_15_p1', '63_2_17_p1', '63_2_19_p1', '63_2_23_p0', '63_2_25_p0', '63_2_35_p1', '63_2_37_p1', '63_2_43_p1', '63_2_46_p1', '63_2_5_p0', '63_2_77_p1', '64_1_23_p0', '64_1_25_p0', '64_1_26_p0', '64_1_29_p1', '64_1_35_p1', '64_1_42_p0', '64_1_48_p1', '64_1_4_p0', '64_1_5_p0', '64_1_7_p0', '64_2_103_p1', '64_2_108_p1', '64_2_13_p0', '64_2_26_p1', '64_2_29_p1', '64_2_39_p0', '64_2_46_p0', '64_2_81_p0', '64_2_83_p0', '64_2_86_p0', '64_2_88_p0', '64_2_89_p0', '65_1_13_p1', '65_2_24_p1', '65_2_25_p1', '65_2_36_p0', '65_2_81_p0', '65_2_87_p0', '66_1_15_p1', '66_2_100_p1', '66_2_58_p1', '66_2_61_p1', '66_2_63_p1', '67_1_21_p0', '67_2_34_p1', '67_2_46_p1', '67_2_53_p0', '67_2_63_p1', '67_2_65_p1', '67_2_81_p0', '67_2_85_p0', '67_2_92_p1', '69_1_30_p1', '69_1_40_p0', '69_1_42_p0', '69_2_20_p0', '69_2_67_p1', '69_2_69_p1', '69_2_71_p1', '69_2_72_p1', '69_2_7_p0', '69_2_8_p0', '6_1_25_p0', '6_1_33_p1', '6_1_43_p0', '6_2_105_p1', '6_2_109_p1', '6_2_17_p0', '6_2_21_p0', '6_2_53_p0', '6_2_63_p0', '6_2_83_p0', '6_2_84_p0', '6_2_85_p0', '6_2_86_p0', '6_2_91_p0', '70_1_21_p0', '70_1_24_p0', '70_1_39_p0', '70_1_47_p1', '70_2_10_p0', '70_2_21_p1', '70_2_43_p0', '70_2_4_p0', '70_2_50_p0', '70_2_58_p1', '70_2_7_p0', '71_2_103_p1', '71_2_10_p0', '71_2_23_p1', '71_2_27_p1', '71_2_70_p1', '71_2_84_p0', '72_1_48_p1', '72_2_102_p0', '72_2_108_p1', '72_2_22_p1', '72_2_34_p1', '72_2_90_p0', '8_2_104_p1', '8_2_12_p0', '8_2_13_p0', '8_2_18_p1', '8_2_19_p1', '8_2_20_p1', '8_2_21_p1', '8_2_22_p1', '8_2_24_p1', '8_2_25_p1', '8_2_26_p1', '8_2_27_p1', '8_2_28_p1', '8_2_29_p1', '8_2_32_p1', '8_2_33_p1', '8_2_37_p0', '8_2_39_p0', '8_2_41_p0', '8_2_42_p0', '8_2_46_p0', '8_2_48_p0', '8_2_52_p1', '8_2_53_p1', '8_2_55_p1', '8_2_56_p1', '8_2_57_p1', '8_2_58_p1', '8_2_59_p1', '8_2_8_p0', '9_2_104_p1', '9_2_105_p1', '9_2_107_p1', '9_2_109_p1', '9_2_38_p1', '9_2_64_p1', '9_2_65_p1', '9_2_69_p1']

We will save these error trials and exclude them during feature extraction

This is the final IK error report when we exclude those that have mean RMS above 4 cm.

RMS error between 0.02–0.04 m is considered acceptable

	session	mean_marker_error_RMS	sd_marker_error_RMS	mean_marker_error_max	sd_marker_error_max
0	10_1	0.033	0.004	0.065	0.013
1	10_2	0.033	0.004	0.066	0.016
2	11_1	0.028	0.004	0.054	0.011
3	11_2	0.029	0.004	0.056	0.013
4	12_1	0.028	0.003	0.056	0.010
5	12_2	0.027	0.004	0.054	0.012
6	13_1	0.031	0.003	0.061	0.013
7	13_2	0.029	0.005	0.059	0.016
8	14_1	0.028	0.003	0.055	0.010
9	14_2	0.028	0.003	0.057	0.011
10	15_1	0.029	0.005	0.059	0.016
11	15_2	0.030	0.008	0.064	0.024
12	16_1	0.031	0.003	0.060	0.010
13	16_2	0.030	0.003	0.060	0.011
14	17_1	0.030	0.003	0.061	0.011

Inverse Dynamics with OpenSim

# This is where we store the data
projectdata = os.path.join(curfolder, 'projectdata')
sessionstotrack = glob.glob(os.path.join(projectdata, 'Session*_*'))
print(sessionstotrack)

# Here we store mass information (weight, height) about participants
META = pd.read_csv(os.path.join(curfolder, '..', '00_raw', 'all_demodata.csv'))

# Get sessionIDs
sessionIDs = []
for session in sessionstotrack:
    sessionIDs.append(session.split(os.sep)[-1])
    sessionIDs[-1] = sessionIDs[-1].split('_')[1]
    # Keep only unique values
    sessionIDs = list(set(sessionIDs))

print('Tracking sessions to process: ' + str(sessionIDs))

['f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_10_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_11_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_11_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_12_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_12_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_13_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_13_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_14_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_14_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_15_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_15_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_16_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_16_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_17_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_17_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_18_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_18_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_19_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_19_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_1_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_1_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_20_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_20_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_21_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_22_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_22_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_23_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_24_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_24_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_26_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_26_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_27_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_27_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_28_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_28_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_29_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_29_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_2_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_2_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_30_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_30_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_31_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_31_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_32_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_32_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_33_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_33_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_34_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_34_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_35_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_35_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_36_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_36_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_37_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_37_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_38_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_38_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_39_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_39_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_3_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_3_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_40_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_40_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_41_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_41_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_43_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_44_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_45_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_45_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_46_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_46_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_47_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_48_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_48_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_49_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_49_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_4_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_50_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_50_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_51_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_51_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_52_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_52_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_53_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_53_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_54_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_54_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_55_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_55_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_56_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_56_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_57_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_57_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_58_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_58_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_59_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_59_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_5_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_5_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_60_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_60_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_61_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_61_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_62_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_62_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_63_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_63_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_64_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_64_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_65_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_65_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_66_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_66_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_67_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_67_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_69_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_69_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_6_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_6_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_70_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_70_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_71_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_71_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_72_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_72_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_7_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_7_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_8_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_8_2', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_9_1', 'f:\\FLESH_ContinuousBodilyEffort\\02_MotionTracking_processing\\projectdata\\Session_9_2']
Tracking sessions to process: ['61', '29', '14', '46', '31', '51', '21', '28', '45', '9', '13', '3', '17', '5', '41', '63', '7', '64', '57', '36', '43', '62', '39', '53', '48', '58', '70', '67', '47', '33', '59', '22', '60', '15', '72', '8', '30', '40', '49', '54', '66', '56', '71', '19', '27', '20', '12', '6', '35', '44', '16', '4', '34', '11', '18', '50', '10', '24', '32', '55', '65', '1', '23', '26', '52', '38', '2', '69', '37']

# Function to update XML file
def update_xml_file(dir, input_file, output_file, new_mass=None, new_model_file=None, new_marker_file=None, new_timerange=None, new_output_model_file=None, new_coord_file=None):
    """
    Updates an XML file with new parameters for model file, time range, coordinate file, and output model file.

    Args:
        dir (str): Directory path where the files are located.
        input_file (str): Path to the input XML file to be updated.
        output_file (str): Path where the updated XML file will be saved.
        new_mass (float, optional): New mass value to update in the XML. Defaults to None.
        new_model_file (str, optional): New path for the model file. Defaults to None.
        new_marker_file (str, optional): New path for the marker file. Defaults to None.
        new_timerange (str, optional): New time range string. Defaults to None.
        new_output_model_file (str, optional): New path for the output model file. Defaults to None.
        new_coord_file (str, optional): New path for the coordinate file. Defaults to None.

    Returns:
        None: The function writes the updated XML to the specified output file.
    """

    # Load the XML document
    tree = ET.parse(input_file)
    root = tree.getroot()

    if 'ID' in input_file:
        # We update model_file
        if new_model_file is not None:
            model_file_element = root.find('.//model_file')
            if model_file_element is not None:
                model_file_element.text = new_model_file

        # We need to update time range
        if new_timerange is not None:
            timerange_elements = root.findall('.//time_range')
            for timerange_element in timerange_elements:
                timerange_element.text = new_timerange

        # And we need to update the path to mot file
        if new_coord_file is not None:
            coord_file_element = root.find('.//coordinates_file')
            if coord_file_element is not None:
                coord_file_element.text = new_coord_file

        # And output_motion_file
        if new_output_model_file is not None:
            output_model_file_element = root.find('.//output_gen_force_file')
            if output_model_file_element is not None:
                output_model_file_element.text = new_output_model_file

    tree.write(output_file, encoding='UTF-8', xml_declaration=True)

# Function to extract time range from a trial
def extract_first_and_last_time(file_path):
    """
    Extracts the first and last time values from a .mot file.

    Args:
        file_path (str): Path to the .mot file to be processed.

    Returns:
        tuple: A tuple containing the first and last time values (first_time, last_time).
    """

    df = pd.read_csv(file_path, sep='\t', skiprows=10)
    
    # Extract the first and last time values, time is the second column
    first_time = df.iloc[0, 0]
    last_time = df.iloc[-1, 0]
    
    return first_time, last_time

# Function to smooth .mot file
def smooth_data(input_path, output_path, smoothing_params, plot=False, plot_column=None):
    """
    Smooths the data in a .mot file using Savitzky-Golay filter.

    Args:
        input_path (str): Path to the input .mot file to be smoothed.
        output_path (str): Path where the smoothed .mot file will be saved.
        smoothing_params (dict): Dictionary containing smoothing parameters:
            - window_length (int): The length of the filter window (must be odd).
            - polyorder (int): The order of the polynomial used to fit the samples.
        plot (bool, optional): If True, plots the unsmoothed vs smoothed data for the specified column. Defaults to False.
        plot_column (str, optional): The column name to plot if plot=True. Defaults to None.

    Returns:
        None: The function writes the smoothed data to the specified output file.
    """

    # Read the entire file
    with open(input_path, 'r') as file:
        lines = file.readlines()

    # Identify header and data section
    header_lines = []
    data_start_index = 0

    for i, line in enumerate(lines):
        if line.strip() == 'endheader':
            header_lines = lines[:i + 1]
            data_start_index = i + 1
            break

    # Extract the column headers and numerical data
    column_headers = lines[data_start_index].split()
    data_lines = lines[data_start_index + 1:]
    data = np.array([list(map(float, line.split())) for line in data_lines])

    # Identify the column index for plotting (if applicable)
    plot_column_idx = column_headers.index(plot_column) if plot_column else None

    # Apply smoothing to each column except 'time'
    smoothed_data = data.copy()
    for col_idx in range(1, data.shape[1]):  # Skip 'time' (assumed to be the first column)
        smoothed_data[:, col_idx] = savgol_filter(
            data[:, col_idx], 
            window_length=smoothing_params['window_length'], 
            polyorder=smoothing_params['polyorder']
        )

    if plot == True:
        #Plot unsmoothed vs smoothed data for the specified column
        if plot_column and plot_column_idx is not None:
            plt.figure(figsize=(10, 6))
            plt.plot(data[:, 0], data[:, plot_column_idx], label='Unsmoothed', alpha=0.7)
            plt.plot(data[:, 0], smoothed_data[:, plot_column_idx], label='Smoothed', alpha=0.7)
            plt.xlabel('Time')
            plt.ylabel(plot_column)
            plt.title(f'Unsmoothed vs Smoothed: {plot_column}')
            plt.legend()
            plt.grid(True)
            plt.show()

    # Write back the original structure with smoothed data
    with open(output_path, 'w') as output_file:
        # Write the header
        output_file.writelines(header_lines)
        # Write the column headers
        output_file.write('\t'.join(column_headers) + '\n')
        # Write the smoothed data row by row - this is necessary to maintain the same formatting, otherwise inverse dynamics will fail
        for row in smoothed_data:
            output_file.write('\t'.join(f'{x:.6f}' for x in row) + '\n')
    

The opensim pipeline has three steps, two of which were already covered within Pose2sim: - scaling - inverse kinematics

The last step is inverse dynamics. In this step, we will use the joint angles estimated in inverse kinematics to estimate the joint forces.

idfile = os.path.join(curfolder, 'Pose2Sim', 'OpenSim_Setup', 'ID_Setup_Pose2Sim_Body135.xml')
id_errors = []

for i in folderstotrack:
    print('working on folder: ' + i)

    sessionID = i.split(os.sep)[-1].split('_')[1]
    ID = sessionID + '_' + i.split(os.sep)[-1].split('_')[2]

    if ID not in sessionIDs:
        print('Skipping sessionID: ' + sessionID)
        continue

    pcnfolders = glob.glob(os.path.join(i, '*'))
    pcnfolders = [x for x in pcnfolders if 'Config' not in x]
    pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
    pcnfolders = [x for x in pcnfolders if 'xml' not in x]
    pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
    pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
    pcnfolders = [x for x in pcnfolders if 'sto' not in x]
    pcnfolders = [x for x in pcnfolders if 'txt' not in x]
    pcnfolders = [x for x in pcnfolders if 'calibration' not in x]

    # Copy necessary files 
    for j in pcnfolders:
        print('working on subfolder: ' + j)

        basename = os.path.basename(j)

        scaled_model = glob.glob(os.path.join(j, 'kinematics', '*.osim'))[0]
        if scaled_model == '':
            print('No scaled model found in ' + j)
            id_errors.append(('No scaled model', j))
            continue

        motfile = glob.glob(os.path.join(j, 'kinematics', '*.mot'))[0]
        if motfile == '':
            print('No mot file found in ' + j)
            id_errors.append(('No mot file', j))
            continue

        # open the file and get time range
        first_time, last_time = extract_first_and_last_time(motfile)
        new_timerange = str(first_time) + ' ' + str(last_time)

        newmotfile = motfile.replace('.mot', '_smoothed.mot')

        # Smooth angles before ID
        smoothing_params = {'window_length': 21, 'polyorder': 3}
        try:
            smooth_data(motfile, newmotfile, smoothing_params, plot=False, plot_column='arm_flex_r')
        except ValueError as e:
            smoothing_params = {'window_length': 5, 'polyorder': 3}
            smooth_data(motfile, newmotfile, smoothing_params, plot=False, plot_column='arm_flex_r')

        # Output force file
        newidfile = newmotfile.replace('.mot', '.sto')

        # Update the XML file
        new_idfile = os.path.join(j, 'kinematics', 'ID_Setup_Pose2Sim_Body135_FLESH_' + basename + '.xml')
        update_xml_file(j, idfile, new_idfile, new_model_file=scaled_model, new_timerange=new_timerange, new_coord_file=newmotfile, new_output_model_file=newidfile)

        print('Inverse Dynamics...')
        try:
            opensim.InverseDynamicsTool(new_idfile).run()
        except:
            print('Error in ID')
            id_errors.append(('ID', j))
            continue

# Get all the folders per session, per participant
pcnfolders = []

for i in folderstotrack:
    pcnfolders_in_session = glob.glob(os.path.join(i, '*'))   
    pcnfolders = pcnfolders + pcnfolders_in_session 

# Get rid of all pontetially confusing files/folders
pcnfolders = [x for x in pcnfolders if 'Config' not in x]
pcnfolders = [x for x in pcnfolders if 'opensim' not in x]
pcnfolders = [x for x in pcnfolders if 'xml' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseDynamics' not in x]
pcnfolders = [x for x in pcnfolders if 'ResultsInverseKinematics' not in x]
pcnfolders = [x for x in pcnfolders if 'sto' not in x]
pcnfolders = [x for x in pcnfolders if 'txt' not in x]
pcnfolders = [x for x in pcnfolders if 'calibration' not in x]
print(f"There are {len(pcnfolders)} pcn folders")

There are 8252 pcn folders

This is an example of visualized time series

References

Pagnon, D., Domalain, M., & Reveret, L. (2022). Pose2Sim: An End-to-End Workflow for 3D Markerless Sports Kinematics—Part 2: Accuracy. Sensors, 22(7, 7), 2712. https://doi.org/10.3390/s22072712

Seth, A., Hicks, J. L., Uchida, T. K., Habib, A., Dembia, C. L., Dunne, J. J., Ong, C. F., DeMers, M. S., Rajagopal, A., Millard, M., Hamner, S. R., Arnold, E. M., Yong, J. R., Lakshmikanth, S. K., Sherman, M. A., Ku, J. P., & Delp, S. L. (2018). OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. PLOS Computational Biology, 14(7), e1006223. https://doi.org/10.1371/journal.pcbi.1006223