{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "accelerator": "GPU", "colab": { "name": "OpenColab.ipynb", "provenance": [], "collapsed_sections": [], "toc_visible": true }, "kernelspec": { "display_name": "Python 3", "name": "python3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "J4eKRuEwayuB" }, "source": [ "

\"OpenSim

\n", "\n", "

What is OpenColab?

\n", "\n", "[**Open**Sim](https://opensim.stanford.edu/) in [Google **Colab**](https://colab.research.google.com/notebooks/intro.ipynb): **OpenColab** Project\n", "\n", "This project aims to help biomechancial modeling in [OpenSim](https://opensim.stanford.edu/) easily accessible to all via Google Colab.\n", "\n", "We use the power of cloud computing (Google Colab) and install OpenSim via Conda. Following installation of relevant packages, we run several examples from OpenSim resources for validations and compare the outcomes with OpenSim GUI. Most of the scripts, images come from OpenSim [GitHub](https://github.com/opensim-org) or [website](https://opensim.stanford.edu/) for comparison purposes. \n", "\n", "To run the following codes, you only need a **Gmail account** and **Internet access** plus passion to learn/enjoy biomechanical modeling. \n", "\n", "For more information, please visit us here:\n", "https://simtk.org/projects/opencolab\n", "\n", "Project Lead: Dr. Hossein Mokhtarzadeh\n", "\n", "Other contributors: Fangwei Jiang;\n", "Andy Shengzhe Zhao;\n", "Fatemeh Malekipour\n", "\n", "Contact: mokhtarzadeh.hossein@gmail.com\n", "\n", "\n", "

\"OpenColab\"

\n", "\n", "**Note:** If you run all the cells at once, it may take somewhere between *25min to 1hr.*\n" ] }, { "cell_type": "markdown", "metadata": { "id": "DY4Dw32kA6K-" }, "source": [ "#Brief Introduction on how to use this notebook!\n", "## We also created several video tutorials and how to use OpenColab in [this link](https://www.tinyurl.com/xukhmnez)." ] }, { "cell_type": "code", "metadata": { "id": "y_0pU9LRA5qu", "cellView": "form" }, "source": [ "#@title How to use this notebook? A video tutorial. For quality watch it as full screen! We also created several video tutorials and how to use OpenColab in this link: www.tinyurl.com/xukhmnez\n", "from IPython.display import YouTubeVideo\n", "YouTubeVideo('iEjd7OSOitg',1000,700) #version 4.3\n", "\n", "# version 4.2\n", "# YouTubeVideo('HR3zLuUfsUo',1000,700)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "h-P7YrMuZmBp" }, "source": [ "#What is OpenSim?\n", "Check the video below from [OpenSim website](https://opensim.stanford.edu/) to learn more about how this tool can be used to develop Biomechancial models." ] }, { "cell_type": "code", "metadata": { "id": "4IO6S0lZZjty", "cellView": "form" }, "source": [ "#@title Run and then watch to learn more about OpenSim \n", "from IPython.display import YouTubeVideo\n", "YouTubeVideo('ScaZP-PZxpI')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "5fCEDCU_qrC0" }, "source": [ "

\"Colaboratory

\n", "\n", "

What is Colaboratory?

\n", "\n", "Colaboratory, or \"Colab\" for short, allows you to write and execute Python in your browser, with \n", "- Zero configuration required\n", "- Free access to GPUs\n", "- Easy sharing\n", "\n", "Whether you're a **student**, a **data scientist** or an **AI researcher**, Colab can make your work easier. Watch [Introduction to Colab](https://www.youtube.com/watch?v=inN8seMm7UI) to learn more, or just get started [here](https://colab.research.google.com/notebooks/intro.ipynb)!" ] }, { "cell_type": "markdown", "metadata": { "id": "NGaSPESzJj3e" }, "source": [ "#What is OpenColab? Scripting in OpenColab\n", "As mentioned by OpenSim team, \"[Scripting](https://simtk-confluence.stanford.edu/display/OpenSim/Scripting) allows you to access OpenSim's functionality through the following programming languages:\n", "\n", "* The scripting shell in the OpenSim GUI (which is a Jython interpreter embedded in the application)\n", "* Matlab\n", "* Python\n", "\n", "In other words, you can access OpenSim's Application Programming Interface without compiling your code in C++.\"\n", "\n", "***OpenColab*** introduces a versatile, & web-based (cloud) scripting platform with Python, zero configuration, free access to GPUs and easy to share." ] }, { "cell_type": "markdown", "metadata": { "id": "ADViBTkvWkuq" }, "source": [ "#Comparison between OpenSim APIs\n", "0: Less established 1: Half established 2: Well established\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
C++Matlab/PythonJython GUIOpenColabNote
Installation on computer2222linux ok, windows for old version
Installation on cloud0102Internet connection needed
Collaboration1112
3D Visualization1221cannot directly run in colab hostruntime, but we can use the local runtime instead
Open Source2222
Pre-/Post-processing1222
Use friendliness0122
Parallel Computing e.g. GPU/TPU1202
Machine Learning1202
Other Opensim modules (e.g. MoCo)2211
Large Scale Simulations1201
Zero Configuration on PC1112
Sensitive Data Sharing (Security)2221
Learning Curve1122
Internet Connection required2121
Sum18241926
\n", " \n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "kDme-AxLFjdx" }, "source": [ "**Summary:** Initially, OpenColab is a free, new and promising software that researchers can easily use and work collaboratively. Although OpenColab and colab have limitations, they are much more extensible and likely to be improved in the future. Besides, Matlab, Python and Jython GUI have implemented some features, but it may be hard to collaborate since they are designed locally. Moreover, Matlab is expensive as it requires a licence to use. Finally, C++ is not an interpreted language. The compatibility of C++ in the cloud is not as high as that of Python and Matlab.\n", "\n", "Please note that we compare them from the end user prespective e.g. a clinician, a researcher, etc. Indeed development is still being done in C.\n", "\n", "**OpenColab > Matlab/Python on PC > Jython GUI> C++**" ] }, { "cell_type": "markdown", "metadata": { "id": "nVyfValqgVFu" }, "source": [ "#Run the following steps to install OpenSim package (mandatory)\n", "Step one may take up to 7 min. [Opensim V4.3 Colab Conda Package](https://anaconda.org/ember123/opencolab)\n", "\n", "\n", "**How to use this notebook?**\n", "---\n", "You have two options: a) Running all the cells at once or b) Step by step.\n", "\n", "Before running any part, please do the following for a quality experience:\n", "\n", "1. Please terminate any session by manage session-->terminate\n", "2. Edit-->Clear All Outputs\n", "3. Connect to hosted runtime\n", "4. If you want to run all, just do Ctrl+F9 or Runtime-->Run all (get a coffee/tea as this may take up to 1hr depending on your connection)\n", "5. If you want to run cells one by one, keep doing them as you move on the Table of Contents\n", "6. Check your network setting and make sure your Chrome version is latest.\n", "\n", "**Note:** If you run up to Static Optimization section, you may do it much faster. [CMC](https://simtk-confluence.stanford.edu/display/OpenSim/How+CMC+Works) tool takes time as we run it for ~2.5 s of gait." ] }, { "cell_type": "code", "metadata": { "id": "FbXIRuVBE7K5", "cellView": "form" }, "source": [ "#@title Step 1: Install OpenSim package from link above & other dependency packages\n", "import time\n", "start_time = time.time()\n", "!pip uninstall -y -q plotly; pip install -q plotly\n", "!pip uninstall -y -q pandas; pip install -q pandas\n", "import pandas as pd \n", "import plotly\n", "pd.set_option('plotting.backend','plotly')\n", "def enable_plotly_in_cell():\n", " import IPython\n", " from plotly.offline import init_notebook_mode\n", " display(IPython.core.display.HTML(''''''))\n", " init_notebook_mode(connected=False)\n", "get_ipython().events.register('pre_run_cell', enable_plotly_in_cell)\n", "!wget -c https://repo.anaconda.com/miniconda/Miniconda3-py37_4.8.3-Linux-x86_64.sh\n", "!chmod +x Miniconda3-py37_4.8.3-Linux-x86_64.sh\n", "!bash ./Miniconda3-py37_4.8.3-Linux-x86_64.sh -b -f -p /usr/local\n", "import sys\n", "sys.path.append('/usr/local/lib/python3.7/site-packages')\n", "!conda install -y --prefix /usr/local -c ember123 opencolab\n", "!apt-get update -y\n", "!apt-get install -y x11-apps\n", "!apt install mesa-utils \n", "!apt-get install xvfb x11-utils\n", "!pip install pyvirtualdisplay\n", "from pyvirtualdisplay import Display\n", "Display(visible=0, size=(1400, 900)).start()\n", "!pip install c3d\n", "!pip install numpy --upgrade\n", "!pip install -q --upgrade ipython\n", "!pip install -q --upgrade ipykernel\n", "!pip install tensorflow==2.0.0.alpha0\n", "!conda install -c plotly -y plotly-orca\n", "!pip install svgutils\n", "!apt install libadolc2\n", "!apt-get install coinor-libipopt-dev \n", "import opensim as osim\n", "print('OpenSim Version Installed is version:',osim.__version__)\n", "print(f'The execution time of OpenSim Package Installation is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "KOsY_BZfZjaX", "cellView": "form" }, "source": [ "#@title Step 2: Let's import dataset from 5 GitHub URLs\n", "# OpenSim models etc from github\n", " \n", "!git clone https://github.com/opensim-org/opensim-models.git # models needed\n", "!git clone https://github.com/HernandezVincent/OpenSim.git # models needed\n", "!cp /content/opensim-models/Models/Gait2354_Simbody/* /content/opensim-models/Pipelines/Gait2354_Simbody\n", "!git clone https://github.com/ESJiang/2354_xml.git\n", "!cp -rf /content/2354_xml/* /content/opensim-models/Pipelines/Gait2354_Simbody\n", "!git clone https://github.com/ESJiang/2354_result.git\n", "!git clone https://github.com/ESJiang/testc3d.git" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "JYnLWfJtXb86", "cellView": "form" }, "source": [ "#@title **`7 useful functions for plotting `** (Mandatory)\n", "#- by Fangwei: feel free to reuse my functions for other opensim model analysis` (mandatory for generating all the data and do plotting in this notebook)**\n", "#this is right foot off timing!\n", "from linecache import getline\n", "import plotly.graph_objects as go\n", "from plotly.subplots import make_subplots\n", "from matplotlib import pyplot\n", "import numpy as np \n", "import cv2\n", "\n", "#dash='dashdot'\n", "#mode=\"lines\"\n", "\n", "def find_nearest_element(time, start, end):\n", " print(min(time, key=lambda x: abs(x-Decimal(start))))\n", " print(time.index(min(time, key=lambda x: abs(x-Decimal(start)))))\n", " print(min(time, key=lambda x: abs(x-Decimal(end))))\n", " print(time.index(min(time, key=lambda x: abs(x-Decimal(end)))))\n", " print((min(time, key=lambda x: abs(x-Decimal(1.840)))+min(time, key=lambda x: abs(x-Decimal(0.600))))/2)\n", " print(time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end)))+min(time, key=lambda x: abs(x-Decimal(start))))/2))))\n", " \n", "def judge_line_number(path):\n", " count = 1\n", " f = open(path,\"r\")\n", " line = f.readline()\n", " while line!=\"\":\n", " s=line[:4]\n", " if s.upper()==\"time\".upper():\n", " return count\n", " else:\n", " count = count + 1\n", " line = f.readline()\n", " \n", "def generate_dict(path):\n", " def sto_Getline(path):\n", " return getline(path, judge_line_number(path)).strip('\\n').split()\n", " for each in dict(zip(sto_Getline(path), range(0,len(sto_Getline(path))))):\n", " print(each, ':', dict(zip(sto_Getline(path), range(0,len(sto_Getline(path)))))[each])\n", "\n", "def three(fig, y1, y1_name, y2, y2_name, y3, y3_name, title):\n", " fig.add_trace(go.Scatter(x=time_5, y=y1, name=y1_name, mode=\"lines\",\n", " line=dict(color='red', width=0.5, dash='solid')))\n", " fig.add_trace(go.Scatter(x=time_10, y=y2, name=y2_name, mode=\"lines\",\n", " line=dict(color='blue', width=0.5, dash='solid')))\n", " fig.add_trace(go.Scatter(x=time_20, y=y3, name=y3_name, mode=\"lines\",\n", " line=dict(color='green', width=0.5, dash='solid')))\n", " fig.update_layout(title=''+title+'', xaxis_title='time', yaxis_title='value',font=dict(\n", " family=\"Courier New, monospace\",\n", " size=24,\n", " #color=\"RebeccaPurple\"\n", " ))\n", " fig.show()\n", "\n", "def four(fig, x, y1, y1_name, y2, y2_name, y3, y3_name, y4, y4_name, title, fig_name):\n", " fig.add_trace(go.Scatter(x=x, y=y1, name=y1_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y2, name=y2_name, mode=\"lines\",\n", " line=dict(color='grey', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y3, name=y3_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='dash')))\n", " fig.add_trace(go.Scatter(x=x, y=y4, name=y4_name, mode=\"lines\",\n", " line = dict(color='grey', width=2, dash='dash')))\n", " fig.update_layout(title={\n", " 'text': ''+title+'',\n", " 'y': 0.9,\n", " 'x': 0.5,\n", " 'xanchor': 'center',\n", " 'yanchor': 'top'}, xaxis_title='time', yaxis_title='value',font=dict(\n", " family=\"Courier New, monospace\",\n", " size=24), template=\"simple_white\"\n", " )\n", " \n", "r_FO = 1.43 #@param {type:\"number\"}\n", "\n", "def four_IDTool(fig, x, y1, y1_name, y2, y2_name, y3, y3_name, y4, y4_name, title, fig_name):\n", " fig.add_trace(go.Scatter(x=x, y=y1, name=y1_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y2, name=y2_name, mode=\"lines\",\n", " line=dict(color='grey', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y3, name=y3_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='dash')))\n", " fig.add_trace(go.Scatter(x=x, y=y4, name=y4_name, mode=\"lines\",\n", " line = dict(color='grey', width=2, dash='dash')))\n", " fig.add_vline(x=r_FO, annotation_position=\"top\",annotation_text=\"r_FO\", line_width=2, line_dash=\"dash\", line_color=\"red\")\n", " #fig.add_vline(x=0.308, annotation_position=\"bottom\",annotation_text=\"l_FO\", line_width=1, line_dash=\"dash\", line_color=\"blue\")\n", " #fig.add_vline(x=0.888, annotation_position=\"top\",annotation_text=\"l-FS\", line_width=1, line_dash=\"dash\", line_color=\"red\")\n", " #fig.add_vline(x=0.918, annotation_position=\"bottom\",annotation_text=\"r_FO\", line_width=1, line_dash=\"dash\", line_color=\"blue\")\n", " fig.update_layout(title={\n", " 'text': ''+title+'',\n", " 'y': 0.95,\n", " 'x': 0.5,\n", " 'xanchor': 'center',\n", " 'yanchor': 'top'}, yaxis_title='value',font=dict(\n", " family=\"Courier New, monospace\",\n", " size=24), template=\"simple_white\"\n", " )\n", "\n", "def six(fig, x, y1, y1_name, y2, y2_name, y3, y3_name, y4, y4_name, y5, y5_name, y6, y6_name, title, fig_name):\n", " fig.add_trace(go.Scatter(x=x, y=y1, name=y1_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y2, name=y2_name, mode=\"lines\",\n", " line=dict(color='grey', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y3, name=y3_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y4, name=y4_name, mode=\"lines\",\n", " line = dict(color='grey', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y5, name=y5_name, mode=\"lines\",\n", " line=dict(color='black', width=2, dash='solid')))\n", " fig.add_trace(go.Scatter(x=x, y=y6, name=y6_name, mode=\"lines\",\n", " line = dict(color='grey', width=2, dash='solid')))\n", " fig.update_layout(title=''+title+'', xaxis_title='time', yaxis_title='value',font=dict(\n", " family=\"Courier New, monospace\",\n", " size=24,\n", " #color=\"RebeccaPurple\"\n", " ))\n", " fig.show()\n", " fig.write_image(fig_name)\n", "\n", "print('finished')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Tj8RnW3VAjRI" }, "source": [ "#Optional (C3D Files): The [3D Standard](https://www.c3d.org/) in Biomechanics.\n", "\n", "C3D file formats are used to collect variety of data in e.g. experiments and record \"*synchronized 3D and analog data*\". Several methods have been developed to read/write c3d files. OpenSim seems to integarete with the following package (ezc3d)\n", "Michaud, B., & Begon, M. (2021). [ezc3d](https://joss.theoj.org/papers/10.21105/joss.02911.pdf): An easy C3D file I/O cross-platform solution for C++, Python and MATLAB. Journal of Open Source Software, 6(58), 2911. " ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "F2XXG9qup7x7" }, "source": [ "#@title From c3d import Reader, an example from https://pypi.org/project/c3d/\n", "from c3d import Reader\n", "r = Reader(open('/content/testc3d/walking5.c3d', 'rb'))\n", "for frame_no, points, analog in r.read_frames():\n", " print('{0.shape} points in this frame'.format(points))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "DFiLoUTcdfxe" }, "source": [ "#Tutorial - Scaling, Inverse Kinematics, and Inverse Dynamics \n", "\n", "---\n", "**Note:** the working directory is:\n", "*/content/opensim-models/Pipelines/Gait2354_Simbody/*\n", "\n", "*Some visualisation functions only can be completed in OpenSim GUI*\n", "\n", "The focus in this notebook is on the following [tutorial](https://simtk-confluence.stanford.edu/display/OpenSim/Tutorial+3+-+Scaling%2C+Inverse+Kinematics%2C+and+Inverse+Dynamics) in OpenSim." ] }, { "cell_type": "markdown", "metadata": { "id": "apN4qZ5Y9wLl" }, "source": [ "## OpenSim Workflow, Model Introduction and Inverse Problem:\n", "OpenSim models consist of several elemets including actuators, joints, ligaments and bones. We use one of the simple models (Gait 2354) that is used for educational and demo purposes since it runs faster. For furthur information please visit the following [link](https://simtk-confluence.stanford.edu/display/OpenSim/Tutorial+3+-+Scaling%2C+Inverse+Kinematics%2C+and+Inverse+Dynamics).\n", "\n", "This model includes \"54 musculotendon actuators to represent 46 muscles in the lower extremities and torso\". \n", "\n", "It is highly recommended to familiarize yourself with OpenSim's basic [tutorials](https://simtk-confluence.stanford.edu/display/OpenSim/Examples+and+Tutorials) on GUI and then run these scripts on Colab to compare and enjoy simplicity of OpenColab. \n", "\n", "In this notebook, we will run inverse problme in Biomechanics using OpenSim gait2354 generic model. Inverse problem uses human motion (joint motions using e.g. external markers on body's anatomical points) and ground reaction forces to calcualate or predict joint angles (e.g. ankle, knee, etc), joint moments, and muscle function (activation/forces) and joint reaction forces during any movements. We will do scaling, inverse kinematics (IK), inverse dynamics (ID), residual reduction algorithm (RRA) and finally computed muscle control (cmc) in OpenSim.\n", "\n", "If you are not familiar with these concepts in Biomechanics, please refer to any basic biomechanics book (e.g. [Biomechanics and Motor Control of Human Movement](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470549148)) or check out [OpenSim website](https://simtk-confluence.stanford.edu/display/OpenSim/Tutorial+3+-+Scaling%2C+Inverse+Kinematics%2C+and+Inverse+Dynamics) for furthur infomration.\n", "\n", "First let's explore an OpenSim model. Please see images below for better undestanding of the workflow in OpenSim from Scaling all the way to muscle function. " ] }, { "cell_type": "markdown", "metadata": { "id": "oyUie_j8AYlL" }, "source": [ "##Explore the model - test \"from opensim import Model\"\n", "\n", "This simple script loads an OpenSim model and prints bodies, jointsets, markerset, probset, etc and number of muscles, coordinates." ] }, { "cell_type": "code", "metadata": { "id": "Z84ohVVL-qTZ", "cellView": "form" }, "source": [ "#@title OpenSim model details e.g. joints, bodies, muscles\n", "from opensim import Model\n", "a = Model(\"/content/opensim-models/Pipelines/Gait2354_Simbody/gait2354_simbody.osim\")\n", "print(\"bodyset:\")\n", "for d in a.getBodySet():\n", " print(\" \" + d.getName())\n", "print()\n", "print(\"Jointset:\")\n", "for d in a.getJointSet():\n", " print(\" \" + d.getName())\n", "print()\n", "print(\"Forceset:\")\n", "for d in a.get_ForceSet():\n", " print(\" \" + d.getName())\n", "print()\n", "print(\"Markerset:\")\n", "for d in a.getMarkerSet():\n", " print(\" \" + d.getName())\n", "print()\n", "print(\"Probeset:\")\n", "for d in a.get_ProbeSet():\n", " print(\" \" + d.getName())\n", "print()\n", "print(\"FrameList:\")\n", "for d in a.getFrameList():\n", " print(d.getName())\n", "print()\n", " \n", "w = Model(\"/content/opensim-models/Pipelines/Gait2354_Simbody/OutputReference/subject01_simbody.osim\")\n", "state = w.initSystem()\n", "count = 0\n", "for d in w.getMuscleList():\n", " count = count + 1\n", "print(f'muscles: {count} \\nbodies: {w.getNumBodies()} \\ndegree: {w.getNumCoordinates()}')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "HGKnoWSMhWp6" }, "source": [ "#[OpenSim](https://journals.plos.org/ploscompbiol/article/figures?id=10.1371/journal.pcbi.1006223): Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement\n", "We will focus on the Inverse Problem in Biomechanics from official [website](https://simtk-confluence.stanford.edu/display/OpenSim/Overview+of+OpenSim+Workflows)\n", "\n", "

\"Inverse

\n", "\n", "
\n", "\n", "#Biomechanical modeling \n", "\n", "\n", "

\"MSK

\n", "\n", "\n", "
\n", "\n", "\n", "#OpenSim Pipeline from Hossein's PhD: [link](https://minerva-access.unimelb.edu.au/handle/11343/38446)\n", "

\"Inverse

\n", "\n", "

Mokhtarzadeh, H. (2013). Anterior cruciate ligament injury mechanism during impact load. PhD thesis, Department of Mechanical Engineering, Melbourne School of Engineering, The University of Melbourne.

\n", "\n", "

Dorn T.W. Computational modelling of lower-limb muscle function in human running. The University of Melbourne, 2012.

\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "1BSp7CeBBSls" }, "source": [ "#Scaling A Generic Musculoskeletal Model\n", "To develop a subject-specific modeling, we scale a generic model in OpenSim. Please see more details [here](https://simtk-confluence.stanford.edu/display/OpenSim/Tutorial+3+-+Scaling%2C+Inverse+Kinematics%2C+and+Inverse+Dynamics).\n", "\n", "Files needed for Scaling from official website:\n", "\n", "

\"Files

\n", "\n", "
\n", "\n", "#Scaling shown in GUI from OpensimGUI\n", "After perforing the scaling in GUI, please note that marker errors will be reported in the messages. We will also report them here in OpenColab in validation part of scaling.\n", "

\"Scaling\"

\n" ] }, { "cell_type": "code", "metadata": { "id": "jHYW1uu09gne", "cellView": "form" }, "source": [ "#@title How to scale a generic model in OpenColab?\n", "from opensim import ScaleTool\n", "import time\n", "start_time = time.time()\n", "ScaleTool(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_Setup_Scale.xml\").run()\n", "print(f'The execution time of ScaleTool is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "bbT5DOJ_HeBj" }, "source": [ "The above code will generate subject01_scaledOnly.osim, subject01_static_output.mot and subject01_scaleSet_applied.xml for you. Now we can load the mot file and save as a new model if need be." ] }, { "cell_type": "code", "metadata": { "id": "CMyT3PARwJ3p", "cellView": "form" }, "source": [ "#@title Scaling Validation\n", "#You can compare these values below and the ones from GUI above.\n", "import os\n", "if os.path.exists('/content/opensim-models/Pipelines/Gait2354_Simbody/opensim.log'):\n", " with open('/content/opensim-models/Pipelines/Gait2354_Simbody/opensim.log','r') as f:\n", " print([line for line in f if 'Frame at (t = 1.0):\t total squared error = ' in line])\n", "else:\n", " with open('/content/opensim.log','r') as f:\n", " print([line for line in f if 'Frame at (t = 1.0):\t total squared error = ' in line])" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "vaF3SbpPKC_V" }, "source": [ "#Inverse Kinematics\n", "Inverse kinematics (IK) aims to find the most accurate joint motions using an optimization method. For this, IK minimizes the distance between experimental and virtual markers on anatomical landmarks. \n", "\n", "Full Body Musculoskeletal Model for Muscle-Driven Simulation of Human OpenSim. \n", "\n", "\n", "\n", "
\n", "\n", "Inverse Kinematics Tool from official website\n", "\n", "\n", "\n", "\n", "
\n", "\n", "#IK shown in GUI from OpensimGUI\n", "After performing scaling, we do Inverse kinematics to calcualte (predict via optimization) the joint motions.\n", "

\"Scaling\"

\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "Mk3Z5XLyKFY0", "cellView": "form" }, "source": [ "#@title How to perform Inverse Kinematics (IK)?\n", "from opensim import InverseKinematicsTool\n", "import time\n", "start_time = time.time()\n", "InverseKinematicsTool(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_Setup_IK.xml\").run()\n", "print(f'The execution time of InverseKinematicsTool is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "ifjtbSehLQEo" }, "source": [ "The above code will generate subject01_walk1_ik.mot, subject01_ik_marker_errors.sto for you. Although subject01_walk1_ik.mot can help us do the visualization, we currently cannot watch it in the colab." ] }, { "cell_type": "markdown", "metadata": { "id": "nZ6OunLgnNJ-" }, "source": [ "##Plotting Inverse Kinematics Results & Validation with GUI\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "_PnAr9l5TKp9" }, "source": [ "#@title subject01_walk1_ik.mot generated by Colab\n", "generate_dict(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_walk1_ik.mot\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "Y8J_m5izTafx" }, "source": [ "#@title subject01_walk1_ik.mot generated by OpensimGUI\n", "generate_dict(\"/content/2354_result/subject01_walk1_ik.mot\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "QxtbxuJuTs0F" }, "source": [ "#@title Read data from subject01_walk1_ik.mot (Colab)\n", "from decimal import Decimal\n", "time, ankle_angle_r, ankle_angle_l, knee_angle_r,knee_angle_l= [],[],[],[],[]\n", "hip_flexion_r, hip_flexion_l=[],[]\n", "def load_file_1(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time.append(Decimal(m[0])) \n", " ankle_angle_r.append(Decimal(m[11]))\n", " ankle_angle_l.append(Decimal(m[18]))\n", " knee_angle_r.append(Decimal(m[10]))\n", " knee_angle_l.append(Decimal(m[17]))\n", " hip_flexion_r.append(Decimal(m[7]))\n", " hip_flexion_l.append(Decimal(m[14]))\n", " c.flush()\n", " c.close()\n", "load_file_1(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_walk1_ik.mot\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "5Zam1lc0Tsl9" }, "source": [ "#@title Read data from subject01_walk1_ik.mot (OpenSimGUI)\n", "from decimal import Decimal\n", "time1, ankle_angle_r1, ankle_angle_l1, knee_angle_r1,knee_angle_l1= [],[],[],[],[]\n", "hip_flexion_r1, hip_flexion_l1 = [],[]\n", "def load_file_1(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time1.append(Decimal(m[0])) \n", " ankle_angle_r1.append(Decimal(m[11]))\n", " ankle_angle_l1.append(Decimal(m[18]))\n", " knee_angle_r1.append(Decimal(m[10]))\n", " knee_angle_l1.append(Decimal(m[17]))\n", " hip_flexion_r1.append(Decimal(m[7]))\n", " hip_flexion_l1.append(Decimal(m[14]))\n", " c.flush()\n", " c.close()\n", "load_file_1(\"/content/2354_result/subject01_walk1_ik.mot\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "p8pVSNfFWd3-" }, "source": [ "#@title Test IKTools panda GUI quick draw full time (0 to 2.5s)\n", "#(code is very simple and execution is short) better for a short demo - fangweij\n", "df1 = pd.DataFrame({'ankle_angle_r_Colab': ankle_angle_r, 'ankle_angle_l_Colab': ankle_angle_l, 'ankle_angle_r_GUI': ankle_angle_r1, 'ankle_angle_l_GUI': ankle_angle_l1}, index = time)\n", "df2 = pd.DataFrame({'knee_angle_r_Colab': knee_angle_r, 'knee_angle_l_Colab': knee_angle_l, 'knee_angle_r_GUI': knee_angle_r1, 'knee_angle_l_GUI': knee_angle_l1}, index = time)\n", "df3 = pd.DataFrame({'hip_flexion_r_Colab': hip_flexion_r, 'hip_flexion_l_Colab': hip_flexion_l, 'hip_flexion_r_GUI': hip_flexion_r1, 'hip_flexion_l_GUI': hip_flexion_l1}, index = time)\n", "df1.index.name = df2.index.name = df3.index.name = 'time'\n", "df1.plot(title='ankle_angle_IK').show()\n", "df2.plot(title='knee_angle_IK').show()\n", "df3.plot(title='hip_flexion_IK').show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "gkGHX7HURljm" }, "source": [ "#@title Find nearest element (gait events, optional)\n", "#- by Fangwei\n", "#This is to find gait events: 0.6s and 1.84s are right side foot strike (i.e. 100% gait cycle). \n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "find_nearest_element(time, start_time, end_time)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "gQ1dEeoT6XZc" }, "source": [ "#@title Test IK panda GUI gait cycle\n", "# - Fangwei (already simplified - prepare for article) (original data starting from 0)\n", "#Don't write them in one line (fig1 = fig2 = fig3 = fig4 = fig5 = fig6 = go.Figure()), we don't want a shadow copy.\n", "\n", "start_time = 0.600#@param {type:\"number\"}\n", "end_time = 1.840#@param {type:\"number\"}\n", "\n", "fig1 = go.Figure()\n", "fig4 = go.Figure()\n", "fig7 = go.Figure()\n", "\n", "four_IDTool(fig7, time, ankle_angle_r, 'r_Colab', ankle_angle_l, 'l_Colab', ankle_angle_r1, 'r_GUI', ankle_angle_l1, 'l_GUI', 'ankle', 'fig1.png')\n", "fig7.add_annotation(textangle=270, text='dorsi', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-5, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-10, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig7.add_annotation(textangle=270, text='plantar', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=5, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=10, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig7.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'ankle flexion (deg)', 'standoff': 35})\n", "fig7.show()\n", "#fig7.write_image('fig7.webp')\n", "fig7.update_layout(title=\"\")\n", "fig7.write_image('fig7.svg')\n", "\n", "four_IDTool(fig4, time, knee_angle_r, 'r_Colab', knee_angle_l, 'l_Colab', knee_angle_r1, 'r_GUI', knee_angle_l1, 'l_GUI', 'knee', 'fig2.png')\n", "fig4.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-20, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-80, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig4.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=0, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig4.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'knee flexion (deg)', 'standoff': 35})\n", "fig4.show()\n", "#fig4.write_image('fig4.webp')\n", "fig4.update_layout(title=\"\")\n", "fig4.write_image('fig4.svg')\n", "\n", "four_IDTool(fig1, time, hip_flexion_r, 'r_Colab', hip_flexion_l, 'l_Colab', hip_flexion_r1, 'r_GUI', hip_flexion_l1, 'l_GUI', 'hip', 'fig3.png')\n", "fig1.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-20, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-60, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig1.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=0, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig1.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'hip flexion (deg)', 'standoff': 35})\n", "fig1.show()\n", "#fig1.write_image('fig1.webp')\n", "fig1.update_layout(title=\"\")\n", "fig1.write_image('fig1.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "L0hAexXDRqTC" }, "source": [ "#Inverse Dynamics\n", "Inverse Dynamics (ID) analysis is used in Biomechanics (inclduing OpenSim) to calculate joint moments (generalized loadings). Figures below show how IK and force data are used to calcualte joint moments. Then these moments can be used to predict muscle activations (or forces). Please continue this tutorial to learn how to quantify muscle function using some optimization methods. OpenSim uses Static Optimization (SO) and Computed Muscle Control (CMC) or EMG-Informed Methods to predict muscle activations. For more information, please visit [Overview of OpenSim Workflows.](https://simtk-confluence.stanford.edu/display/OpenSim/Overview+of+OpenSim+Workflows)\n", "\n", "\n", "ID process in OpenSim from official website\n", "

\"Scaling\"

\n", "
\n", "Files needed for Inverse dynamic from official website\n", "

\"Files

\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "HRb91TXHR7hA" }, "source": [ "#@title Scripts to perform Inverse Dynamics (ID)\n", "from opensim import InverseDynamicsTool\n", "import time\n", "start_time = time.time()\n", "InverseDynamicsTool(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_Setup_InverseDynamics.xml\").run()\n", "print(f'The execution time of InverseDynamicsTool is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "vYLekyopSJcv" }, "source": [ "The above code will save the sto data(inverse_dynamics.sto) to ResultsInverseDynamics folder" ] }, { "cell_type": "markdown", "metadata": { "id": "rPy0K3MrYfWF" }, "source": [ "##Plotting Inverse Dynamics Results & Validation with GUI" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "uRd3lB7WSh5Z" }, "source": [ "#@title Inverse_dynamics.sto generated by Colab\n", "generate_dict(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsInverseDynamics/inverse_dynamics.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "W7L5PGLTj385" }, "source": [ "#@title Inverse_dynamics.sto generated by OpensimGUI\n", "generate_dict(\"/content/2354_result/inverse_dynamics.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "t2Gj4MWfTEIW" }, "source": [ "#@title Read data from inverse_dynamics.sto (Colab)\n", "from decimal import Decimal\n", "time, ankle_angle_r_moment, ankle_angle_l_moment,pelvis_tx_force,pelvis_ty_force,pelvis_tz_force,knee_angle_r_moment,knee_angle_l_moment,mtp_angle_r_moment,mtp_angle_l_moment,subtalar_angle_r_moment,subtalar_angle_l_moment = [],[],[],[],[],[],[],[],[],[],[],[]\n", "hip_flexion_r_moment,hip_flexion_l_moment,hip_adduction_r_moment,hip_adduction_l_moment,hip_rotation_r_moment,hip_rotation_l_moment=[],[],[],[],[],[]\n", "def load_file_1(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time.append(Decimal(m[0])) \n", " ankle_angle_r_moment.append(Decimal(m[18]))\n", " ankle_angle_l_moment.append(Decimal(m[19]))\n", " knee_angle_r_moment.append(Decimal(m[16]))\n", " knee_angle_l_moment.append(Decimal(m[17]))\n", " hip_flexion_r_moment.append(Decimal(m[7]))\n", " hip_flexion_l_moment.append(Decimal(m[10]))\n", " c.flush()\n", " c.close()\n", "load_file_1(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsInverseDynamics/inverse_dynamics.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "_1Z2nAf3mN6z" }, "source": [ "#@title Read data from inverse_dynamics.sto (OpensimGUI)\n", "from decimal import Decimal\n", "time1, ankle_angle_r_moment1, ankle_angle_l_moment1,pelvis_tx_force1,pelvis_ty_force1,pelvis_tz_force1,knee_angle_r_moment1,knee_angle_l_moment1,mtp_angle_r_moment1,mtp_angle_l_moment1,subtalar_angle_r_moment1,subtalar_angle_l_moment1 = [],[],[],[],[],[],[],[],[],[],[],[]\n", "hip_flexion_r_moment1,hip_flexion_l_moment1,hip_adduction_r_moment1,hip_adduction_l_moment1,hip_rotation_r_moment1,hip_rotation_l_moment1=[],[],[],[],[],[]\n", "def load_file_11(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time1.append(Decimal(m[0])) \n", " ankle_angle_r_moment1.append(Decimal(m[18]))\n", " ankle_angle_l_moment1.append(Decimal(m[19]))\n", " knee_angle_r_moment1.append(Decimal(m[16]))\n", " knee_angle_l_moment1.append(Decimal(m[17]))\n", " hip_flexion_r_moment1.append(Decimal(m[7]))\n", " hip_flexion_l_moment1.append(Decimal(m[10]))\n", " c.flush()\n", " c.close()\n", "load_file_11(\"/content/2354_result/inverse_dynamics.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "10zLyviXBK8X" }, "source": [ "#@title X-axis shifting to start from zero\n", "time[:]=[x-time[0] for x in time]" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "6cmuoVKfHfhA" }, "source": [ "#@title Test inversedynamics panda GUI quick draw full time (0 to 2.5s)\n", "# (code is very simple and execution is short) better for a short demo - Fangwei\n", "df1 = pd.DataFrame({'ankle_angle_r_Colab': ankle_angle_r_moment, 'ankle_angle_l_Colab': ankle_angle_l_moment, 'ankle_angle_r_GUI': ankle_angle_r_moment1, 'ankle_angle_l_GUI': ankle_angle_l_moment1}, index = time)\n", "df2 = pd.DataFrame({'knee_angle_r_Colab': knee_angle_r_moment, 'knee_angle_l_Colab': knee_angle_l_moment, 'knee_angle_r_GUI': knee_angle_r_moment1, 'knee_angle_l_GUI': knee_angle_l_moment1}, index = time)\n", "df3 = pd.DataFrame({'hip_flexion_r_Colab': hip_flexion_r_moment, 'hip_flexion_l_Colab': hip_flexion_l_moment, 'hip_flexion_r_GUI': hip_flexion_r_moment1, 'hip_flexion_l_GUI': hip_flexion_l_moment1}, index = time)\n", "df1.index.name = df2.index.name = df3.index.name ='time'\n", "df1.plot(title='ankle_angle_ID').show()\n", "df2.plot(title='knee_angle_ID').show()\n", "df3.plot(title='hip_flexion_ID').show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "Dmv4-jpHU2Ma" }, "source": [ "#@title Find nearest element (gait events, optional)\n", "#- by Fangwei\n", "#This is to find gait events: 0.6s and 1.84s are right side foot strike (i.e. 100% gait cycle). \n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "find_nearest_element(time, start_time, end_time)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "1mmY_pQj_nh8" }, "source": [ "#@title Test ID gait cycle\n", "\n", "# test inverse dynamics panda GUI subplot & slow draw - Fangwei (already simplified - prepare for artical) (original data starting from 0)\n", "#Don't write them in one line (fig1 = fig2 = fig3 = fig4 = fig5 = fig6 = go.Figure()), we don't want a shadow copy.\n", "\n", "start_time = 0.600#@param {type:\"number\"}\n", "end_time = 1.840#@param {type:\"number\"}\n", "\n", "fig2 = go.Figure()\n", "fig5 = go.Figure()\n", "fig8 = go.Figure()\n", "\n", "four_IDTool(fig8, time, ankle_angle_r_moment, 'r_Colab', ankle_angle_l_moment, 'l_Colab', ankle_angle_r_moment1, 'r_GUI', ankle_angle_l_moment1, 'l_GUI', 'ankle', 'fig1.png')\n", "fig8.add_annotation(textangle=270, text='dorsi', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-60, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-120, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig8.add_annotation(textangle=270, text='plantar', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-40, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=20, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig8.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'ankle flex. moment (N.m)', 'standoff': 35})\n", "fig8.show()\n", "#fig8.write_image('fig8.webp')\n", "fig8.update_layout(title=\"\")\n", "fig8.write_image('fig8.svg')\n", "\n", "four_IDTool(fig5, time, knee_angle_r_moment, 'r_Colab', knee_angle_l_moment, 'l_Colab', knee_angle_r_moment1, 'r_GUI', knee_angle_l_moment1, 'l_GUI', 'knee', 'fig2.png')\n", "fig5.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-20, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-80, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig5.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=0, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig5.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'knee flex. moment (N.m)', 'standoff': 35})\n", "fig5.show()\n", "#fig5.write_image('fig5.webp')\n", "fig5.update_layout(title=\"\")\n", "fig5.write_image('fig5.svg')\n", "\n", "four_IDTool(fig2, time, hip_flexion_r_moment, 'r_Colab', hip_flexion_l_moment, 'l_Colab', hip_flexion_r_moment1, 'r_GUI', hip_flexion_l_moment1, 'l_GUI', 'hip', 'fig3.png')\n", "fig2.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-20, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-60, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig2.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=0, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig2.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'hip flex. moment (N.m)', 'standoff': 35})\n", "fig2.show()\n", "#fig2.write_image('fig2.webp')\n", "fig2.update_layout(title=\"\")\n", "fig2.write_image('fig2.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "z2PSxVPNTppb" }, "source": [ "#Residual Reduction Algorithm (RRA)\n", "Toward dynamically consistant models using RRATool. The following inputs are required for RRA to work.\n", "\n", "

\"MSK

\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "tKPb1ZWmTt3J", "cellView": "form" }, "source": [ "#@title Residual Reduction Algorithm (RRA) in OpenSim\n", "from opensim import RRATool\n", "import time\n", "start_time = time.time()\n", "RRATool(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_Setup_RRA.xml\").run()\n", "print(f'The execution time of RRATool is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "sPLp5gWNbwWz" }, "source": [ "##Plotting RRA results & Validation with GUI" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "cZCvZg9dT6A2" }, "source": [ "#@title subject01_walk1_RRA_Actuation_force.sto generated by Colab\n", "generate_dict(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsRRA/subject01_walk1_RRA_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "Gb-yhLTKrtuS" }, "source": [ "#@title subject01_walk1_RRA_Actuation_force.sto generated by OpensimGUI\n", "generate_dict(\"/content/2354_result/subject01_walk1_RRA_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "P1u3oG94RMd_" }, "source": [ "#@title Read data from subject01_walk1_RRA_Actuation_force.sto (Colab)\n", "from decimal import Decimal\n", "time, pelvis_tx, pelvis_ty, pelvis_tz, ankle_angle_r, ankle_angle_l, knee_angle_l, knee_angle_r, hip_flexion_l, hip_flexion_r = [],[],[],[],[],[],[],[],[],[]\n", "hip_adduction_r, hip_adduction_l, hip_rotation_r, hip_rotation_l = [],[],[],[]\n", "def load_file_2(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time.append(Decimal(m[0])) \n", " ankle_angle_r.append(Decimal(m[11]))\n", " ankle_angle_l.append(Decimal(m[18]))\n", " knee_angle_r.append(Decimal(m[10]))\n", " knee_angle_l.append(Decimal(m[17]))\n", " hip_flexion_l.append(Decimal(m[14]))\n", " hip_flexion_r.append(Decimal(m[7]))\n", " c.flush()\n", " c.close()\n", "load_file_2(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsRRA/subject01_walk1_RRA_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "KSLMm4ZiuBYO" }, "source": [ "#@title Read data from subject01_walk1_RRA_Actuation_force.sto (OpensimGUI)\n", "from decimal import Decimal\n", "time1, pelvis_tx1, pelvis_ty1, pelvis_tz1, ankle_angle_r1, ankle_angle_l1, knee_angle_l1, knee_angle_r1, hip_flexion_l1, hip_flexion_r1 = [],[],[],[],[],[],[],[],[],[]\n", "hip_adduction_r1, hip_adduction_l1, hip_rotation_r1, hip_rotation_l1 = [],[],[],[]\n", "def load_file_22(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time1.append(Decimal(m[0])) \n", " ankle_angle_r1.append(Decimal(m[11]))\n", " ankle_angle_l1.append(Decimal(m[18]))\n", " knee_angle_r1.append(Decimal(m[10]))\n", " knee_angle_l1.append(Decimal(m[17]))\n", " hip_flexion_l1.append(Decimal(m[14]))\n", " hip_flexion_r1.append(Decimal(m[7]))\n", " c.flush()\n", " c.close()\n", "load_file_22(\"/content/2354_result/subject01_walk1_RRA_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "mjEGYBsaK5GR" }, "source": [ "#@title Test RRATools panda GUI quick draw full time (0 to 2.5s)\n", "# (code is very simple and execution is short) better for a short demo - Fangwei\n", "df1 = pd.DataFrame({'ankle_angle_r_Colab': ankle_angle_r, 'ankle_angle_l_Colab': ankle_angle_l, 'ankle_angle_r_GUI': ankle_angle_r1, 'ankle_angle_l_GUI': ankle_angle_l1}, index = time)\n", "df2 = pd.DataFrame({'knee_angle_r_Colab': knee_angle_r, 'knee_angle_l_Colab': knee_angle_l, 'knee_angle_r_GUI': knee_angle_r1, 'knee_angle_l_GUI': knee_angle_l1}, index = time)\n", "df3 = pd.DataFrame({'hip_flexion_r_Colab': hip_flexion_r, 'hip_flexion_l_Colab': hip_flexion_l, 'hip_flexion_r_GUI': hip_flexion_r1, 'hip_flexion_l_GUI': hip_flexion_l1}, index = time)\n", "df1.index.name = df2.index.name = df3.index.name = 'time'\n", "df1.plot(title='ankle_angle_RRA').show()\n", "df2.plot(title='knee_angle_RRA').show()\n", "df3.plot(title='hip_flexion_RRA').show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "7okmz7TkVz5F" }, "source": [ "#@title Find nearest element (gait events, optional)\n", "#- by Fangwei\n", "#This is to find gait events: 0.6s and 1.84s are right side foot strike (i.e. 100% gait cycle). \n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "find_nearest_element(time, start_time, end_time)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "ge7srVQSWtLu" }, "source": [ "#@title Test RRATool gait cycle\n", "# - Fangwei (already simplified - prepare for article) (original data starting from 0)(gait cycle)\n", "#Don't write them in one line (fig1 = fig2 = fig3 = fig4 = fig5 = fig6 = go.Figure()), we don't want a shadow copy.\n", "\n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "\n", "fig3 = go.Figure()\n", "fig6 = go.Figure()\n", "fig9 = go.Figure()\n", "\n", "four_IDTool(fig9, time, ankle_angle_r, 'r_Colab', ankle_angle_l, 'l_Colab', ankle_angle_r1, 'r_GUI', ankle_angle_l1, 'l_GUI', 'ankle', 'fig1.png')\n", "fig9.add_annotation(textangle=270, text='dorsi', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-50, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-100, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig9.add_annotation(textangle=270, text='plantar', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=0, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=50, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig9.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(1.840)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'ankle flex. moment (N.m)', 'standoff': 35})\n", "fig9.show()\n", "#fig9.write_image('fig9.webp')\n", "fig9.update_layout(title=\"\")\n", "fig9.write_image('fig9.svg')\n", "\n", "four_IDTool(fig6, time, knee_angle_r, 'r_Colab', knee_angle_l, 'l_Colab', knee_angle_r1, 'r_GUI', knee_angle_l1, 'l_GUI', 'knee', 'fig2.png')\n", "fig6.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-40, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-70, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig6.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-30, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=0, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig6.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(1.840)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'knee flex. moment (N.m)', 'standoff': 35})\n", "fig6.show()\n", "#fig6.write_image('fig6.webp')\n", "fig6.update_layout(title=\"\")\n", "fig6.write_image('fig6.svg')\n", "\n", "four_IDTool(fig3, time, hip_flexion_r, 'r_Colab', hip_flexion_l, 'l_Colab', hip_flexion_r1, 'r_GUI', hip_flexion_l1, 'l_GUI', 'hip', 'fig3.png')\n", "fig3.add_annotation(textangle=270, text='extension', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=-10, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=-40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='right', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "\n", "fig3.add_annotation(textangle=270, text='flexion', ax=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], axref='x', ay=10, ayref='y',\n", " x=time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], xref='x', y=40, yref='y', xshift=-60, arrowcolor=\"black\", xanchor='left', yanchor='bottom',\n", " arrowwidth=2, arrowhead=2, font=dict(size=30, color='black'))\n", "fig3.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(1.840)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(1.840))))]],\n", " ticktext=['0','50','100']\n", "), xaxis_title='Gait cycle (%)', yaxis_title={'text': 'hip flex. moment (N.m)', 'standoff': 35})\n", "fig3.show()\n", "#fig3.write_image('fig3.webp')\n", "fig3.update_layout(title=\"\")\n", "fig3.write_image('fig3.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "9FqB3lvAoLTq" }, "source": [ "#@title save subplot (ID,IK,RRA) to svg\n", "# - Fangwei 04/06/2021\n", "import svgutils.transform as sg\n", "\n", "fig = sg.SVGFigure(\"16cm\", \"16cm\")\n", "fig1 = sg.fromfile('/content/fig1.svg')\n", "fig2 = sg.fromfile('/content/fig2.svg')\n", "fig3 = sg.fromfile('/content/fig3.svg')\n", "fig4 = sg.fromfile('/content/fig4.svg')\n", "fig5 = sg.fromfile('/content/fig5.svg')\n", "fig6 = sg.fromfile('/content/fig6.svg')\n", "fig7 = sg.fromfile('/content/fig7.svg')\n", "fig8 = sg.fromfile('/content/fig8.svg')\n", "fig9 = sg.fromfile('/content/fig9.svg')\n", "\n", "plot1 = fig1.getroot()\n", "plot2 = fig2.getroot()\n", "plot3 = fig3.getroot()\n", "plot4 = fig4.getroot()\n", "plot5 = fig5.getroot()\n", "plot6 = fig6.getroot()\n", "plot7 = fig7.getroot()\n", "plot8 = fig8.getroot()\n", "plot9 = fig9.getroot()\n", "\n", "plot1.moveto(100, 0, 1, None)\n", "plot2.moveto(800, 0, 1, None)\n", "plot3.moveto(1500, 0, 1, None)\n", "plot4.moveto(100, 500, 1, None)\n", "plot5.moveto(800, 500, 1, None)\n", "plot6.moveto(1500, 500, 1, None)\n", "plot7.moveto(100, 1000, 1, None)\n", "plot8.moveto(800, 1000, 1, None)\n", "plot9.moveto(1500, 1000, 1, None)\n", "\n", "txt1 = sg.TextElement(300, 30, \"Inverse kinematics\", size=40, weight=\"bold\")\n", "txt2 = sg.TextElement(1000, 30, \"Inverse dynamics\", size=40, weight=\"bold\")\n", "txt3 = sg.TextElement(1500, 30, \"Residual reduction algorithm\", size=40, weight=\"bold\")\n", "txt4 = sg.TextElement(-150, -20, \"Hip\", size=40, weight=\"bold\")\n", "txt4.rotate(-90, 100)\n", "txt5 = sg.TextElement(-650, -20, \"Knee\", size=40, weight=\"bold\")\n", "txt5.rotate(-90, 100)\n", "txt6 = sg.TextElement(-1150, -20, \"Ankle\", size=40, weight=\"bold\")\n", "txt6.rotate(-90, 100)\n", "\n", "fig.append([plot1, plot2, plot3, plot4, plot5, plot6, plot7, plot8, plot9])\n", "fig.append([txt1, txt2, txt3, txt4, txt5, txt6])\n", "fig.set_size([\"58cm\", \"40cm\"])\n", "\n", "fig.save(\"/content/fig_subplot_ik_id_rra.svg\")\n", "print('Finished, please check fig_subplot_ik_id_rra.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "GD7JArnwkrYc" }, "source": [ "#Static Optimization (SO) in OpenSim\n", "Static optimization is an extension to inverse dynamics that further resolves the net joint moments into individual muscle forces at each instant in time. The muscle forces are resolved by minimizing the sum of squared (or other power) muscle activations.\n", "\n", "\n", "Files required to perform Static Optimization (SO).\n", "

\"Files

" ] }, { "cell_type": "code", "metadata": { "id": "ud4PFXMOmAy6", "cellView": "form" }, "source": [ "#@title Static Optimization (SO) in OpenSim 30s\n", "# from this link: https://www.gitmemory.com/issue/opensim-org/opensim-core/2076/615401510\n", "import os \n", "def perform_so(model_file, ik_file, grf_file, grf_xml, reserve_actuators, results_dir):\n", " \"\"\"Performs Static Optimization using OpenSim.\n", "\n", " Parameters\n", " ----------\n", " model_file: str\n", " OpenSim model (.osim)\n", " ik_file: str\n", " kinematics calculated from Inverse Kinematics\n", " grf_file: str\n", " the ground reaction forces\n", " grf_xml: str\n", " xml description containing how to apply the GRF forces\n", " reserve_actuators: str\n", " path to the reserve actuator .xml file\n", " results_dir: str\n", " directory to store the results\n", " \"\"\"\n", " # model\n", " model = osim.Model(model_file)\n", "\n", " # prepare external forces xml file\n", " name = os.path.basename(grf_file)[:-8]\n", " external_loads = osim.ExternalLoads(grf_xml, True)\n", " # external_loads.setExternalLoadsModelKinematicsFileName(ik_file) \n", " external_loads.setDataFileName(grf_file)\n", " # external_loads.setLowpassCutoffFrequencyForLoadKinematics(6)\n", " external_loads.printToXML(results_dir + name + '.xml')\n", "\n", " # add reserve actuators\n", " force_set = osim.SetForces(reserve_actuators, True)\n", " force_set.setMemoryOwner(False) # model will be the owner\n", " for i in range(0, force_set.getSize()):\n", " model.updForceSet().append(force_set.get(i))\n", "\n", " # construct static optimization\n", " motion = osim.Storage(ik_file)\n", " static_optimization = osim.StaticOptimization()\n", " static_optimization.setStartTime(motion.getFirstTime())\n", " static_optimization.setEndTime(motion.getLastTime())\n", " static_optimization.setUseModelForceSet(True)\n", " static_optimization.setUseMusclePhysiology(True)\n", " static_optimization.setActivationExponent(2)\n", " static_optimization.setConvergenceCriterion(0.0001)\n", " static_optimization.setMaxIterations(100)\n", " model.addAnalysis(static_optimization)\n", "\n", " # analysis\n", " analysis = osim.AnalyzeTool(model)\n", " analysis.setName(name)\n", " analysis.setModel(model)\n", " analysis.setInitialTime(motion.getFirstTime())\n", " analysis.setFinalTime(motion.getLastTime())\n", " analysis.setLowpassCutoffFrequency(6)\n", " analysis.setCoordinatesFileName(ik_file) \n", " analysis.setExternalLoadsFileName(results_dir + name + '.xml')\n", " analysis.setLoadModelAndInput(True)\n", " analysis.setResultsDir(results_dir)\n", " analysis.run()\n", " so_force_file = results_dir + name + '_so_forces.sto'\n", " so_activations_file = results_dir + name + '_so_activations.sto'\n", " return (so_force_file, so_activations_file)\n", "\n", "#these may need to change for our final paper (I just did it quickly to run it once):\n", "model_file = '/content/opensim-models/Pipelines/Gait2354_Simbody/OutputReference/subject01_simbody.osim'\n", "ik_file = '/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_walk1_ik.mot'\n", "grf_file = '/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_walk1_grf.mot'\n", "grf_xml = '/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_walk1_grf.xml'\n", "reserve_actuators = '/content/2354_result/Reserve_Actuators.xml'\n", "results_dir = '/content/opensim-models/Pipelines/Gait2354_Simbody/OutputReference/ResultsSO1'\n", "perform_so(model_file, ik_file, grf_file, grf_xml, reserve_actuators, results_dir)\n", "print(\"SO finished\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "WY1tdUpfnc2j" }, "source": [ "##Plotting SO results & Validation with GUI" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "bkFaECZ9mnZH" }, "source": [ "#@title Read data from subject01_walk1_StaticOptimization_force.sto (Colab)\n", "from decimal import Decimal\n", "time, bifemlh_r, bifemsh_r, bifemlh_l, bifemsh_l=[],[],[],[],[]\n", "rect_fem_r,vas_int_r,rect_fem_l,vas_int_l=[],[],[],[]\n", "glut_med1_r,glut_med2_r,glut_med3_r,glut_med1_l,glut_med2_l,glut_med3_l=[],[],[],[],[],[]\n", "med_gas_r, med_gas_l, soleus_r, soleus_l, tib_ant_r, tib_ant_l=[],[],[],[],[],[]\n", "def load_file_2(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time.append(Decimal(m[0]))\n", " # rhams\n", " bifemlh_r.append(Decimal(m[4]))\n", " bifemsh_r.append(Decimal(m[5]))\n", " # lhams\n", " bifemlh_l.append(Decimal(m[28]))\n", " bifemsh_l.append(Decimal(m[29]))\n", " # rquad\n", " rect_fem_r.append(Decimal(m[19]))\n", " vas_int_r.append(Decimal(m[20]))\n", " # lquad\n", " rect_fem_l.append(Decimal(m[43]))\n", " vas_int_l.append(Decimal(m[44]))\n", " # rgluts\n", " glut_med1_r.append(Decimal(m[1]))\n", " glut_med2_r.append(Decimal(m[2]))\n", " glut_med3_r.append(Decimal(m[3]))\n", " # lgluts\n", " glut_med1_l.append(Decimal(m[25]))\n", " glut_med2_l.append(Decimal(m[26]))\n", " glut_med3_l.append(Decimal(m[27]))\n", " # rgast\n", " med_gas_r.append(Decimal(m[21]))\n", " # lgast\n", " med_gas_l.append(Decimal(m[45]))\n", " # rsole\n", " soleus_r.append(Decimal(m[22]))\n", " # lsole\n", " soleus_l.append(Decimal(m[46]))\n", " # rTAs\n", " tib_ant_r.append(Decimal(m[24]))\n", " # lTAs\n", " tib_ant_l.append(Decimal(m[48]))\n", " c.flush()\n", " c.close()\n", "load_file_2(\"/content/opensim-models/Pipelines/Gait2354_Simbody/OutputReference/ResultsSO1/subject01_walk1_StaticOptimization_force.sto\")\n", "rhams,lhams,rquad,lquad,rgluts,lgluts,rgast,lgast,rsole,lsole,rTAs,lTAs=[],[],[],[],[],[],[],[],[],[],[],[]\n", "mid_r,mid_l=[],[]\n", "rhams[:] = [x + y for x, y in zip(bifemlh_r, bifemsh_r)]\n", "lhams[:] = [x + y for x, y in zip(bifemlh_l, bifemsh_l)]\n", "rquad[:] = [x + y for x, y in zip(rect_fem_r,vas_int_r)]\n", "lquad[:] = [x + y for x, y in zip(rect_fem_l,vas_int_l)]\n", "mid_r[:] = [x + y for x, y in zip(glut_med1_r, glut_med2_r)]\n", "mid_l[:] = [x + y for x, y in zip(glut_med1_l, glut_med2_l)]\n", "rgluts[:] = [x + y for x, y in zip(mid_r, glut_med3_r)]\n", "lgluts[:] = [x + y for x, y in zip(mid_l, glut_med3_l)]\n", "rgast = med_gas_r\n", "lgast = med_gas_l\n", "rsole = soleus_r\n", "lsole = soleus_l\n", "rTAs = tib_ant_r\n", "lTAs = tib_ant_l" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "2stw43sAmzVa", "cellView": "form" }, "source": [ "#@title Read data from 3DGaitModel2354_StaticOptimization_force.sto (OpenSimGUI)\n", "from decimal import Decimal\n", "time1, bifemlh_r1, bifemsh_r1, bifemlh_l1, bifemsh_l1=[],[],[],[],[]\n", "rect_fem_r1,vas_int_r1,rect_fem_l1,vas_int_l1=[],[],[],[]\n", "glut_med1_r1,glut_med2_r1,glut_med3_r1,glut_med1_l1,glut_med2_l1,glut_med3_l1=[],[],[],[],[],[]\n", "med_gas_r1, med_gas_l1, soleus_r1, soleus_l1, tib_ant_r1, tib_ant_l1=[],[],[],[],[],[]\n", "def load_file_2(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time1.append(Decimal(m[0]))\n", " # rhams\n", " bifemlh_r1.append(Decimal(m[4]))\n", " bifemsh_r1.append(Decimal(m[5]))\n", " # lhams\n", " bifemlh_l1.append(Decimal(m[28]))\n", " bifemsh_l1.append(Decimal(m[29]))\n", " # rquad\n", " rect_fem_r1.append(Decimal(m[19]))\n", " vas_int_r1.append(Decimal(m[20]))\n", " # lquad\n", " rect_fem_l1.append(Decimal(m[43]))\n", " vas_int_l1.append(Decimal(m[44]))\n", " # rgluts\n", " glut_med1_r1.append(Decimal(m[1]))\n", " glut_med2_r1.append(Decimal(m[2]))\n", " glut_med3_r1.append(Decimal(m[3]))\n", " # lgluts\n", " glut_med1_l1.append(Decimal(m[25]))\n", " glut_med2_l1.append(Decimal(m[26]))\n", " glut_med3_l1.append(Decimal(m[27]))\n", " # rgast\n", " med_gas_r1.append(Decimal(m[21]))\n", " # lgast\n", " med_gas_l1.append(Decimal(m[45]))\n", " # rsole\n", " soleus_r1.append(Decimal(m[22]))\n", " # lsole\n", " soleus_l1.append(Decimal(m[46]))\n", " # rTAs\n", " tib_ant_r1.append(Decimal(m[24]))\n", " # lTAs\n", " tib_ant_l1.append(Decimal(m[48]))\n", " c.flush()\n", " c.close()\n", "load_file_2(\"/content/2354_result/ResultsSO_1/3DGaitModel2354_StaticOptimization_force.sto\")\n", "rhams1,lhams1,rquad1,lquad1,rgluts1,lgluts1,rgast1,lgast1,rsole1,lsole1,rTAs1,lTAs1=[],[],[],[],[],[],[],[],[],[],[],[]\n", "mid_r,mid_l=[],[]\n", "rhams1[:] = [x + y for x, y in zip(bifemlh_r1, bifemsh_r1)]\n", "lhams1[:] = [x + y for x, y in zip(bifemlh_l1, bifemsh_l1)]\n", "rquad1[:] = [x + y for x, y in zip(rect_fem_r1,vas_int_r1)]\n", "lquad1[:] = [x + y for x, y in zip(rect_fem_l1,vas_int_l1)]\n", "mid_r[:] = [x + y for x, y in zip(glut_med1_r1, glut_med2_r1)]\n", "mid_l[:] = [x + y for x, y in zip(glut_med1_l1, glut_med2_l1)]\n", "rgluts1[:] = [x + y for x, y in zip(mid_r, glut_med3_r1)]\n", "lgluts1[:] = [x + y for x, y in zip(mid_l, glut_med3_l1)]\n", "rgast1 = med_gas_r1\n", "lgast1 = med_gas_l1\n", "rsole1 = soleus_r1\n", "lsole1 = soleus_l1\n", "rTAs1 = tib_ant_r1\n", "lTAs1 = tib_ant_l1" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "li0Fvb4-mSzd" }, "source": [ "#@title x-axis shifting\n", "time[:]=[x-time[0] for x in time]" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "6ESgrPBQngid" }, "source": [ "#@title Test SO panda GUI quick draw full time (0 to 2.5s)\n", "# (code is very simple and execution is short) better for a short demo - Fangwei\n", "df1 = pd.DataFrame({'rhams_Colab': rhams, 'lhams_Colab': lhams, 'rhams_GUI': rhams1, 'lhams_GUI': lhams1}, index = time)\n", "df2 = pd.DataFrame({'rquad_Colab': rquad, 'lquad_Colab': lquad, 'rquad_GUI': rquad1, 'lquad_GUI': lquad1}, index = time)\n", "df3 = pd.DataFrame({'rgluts_Colab': rgluts, 'lgluts_Colab': lgluts, 'rgluts_GUI': rgluts1, 'lgluts_GUI': lgluts1}, index = time)\n", "df4 = pd.DataFrame({'rgast_Colab': rgast, 'lgast_Colab': lgast, 'rgast_GUI': rgast1, 'lgast_GUI': lgast1}, index = time)\n", "df5 = pd.DataFrame({'rsole_Colab': rsole, 'lsole_Colab': lsole, 'rsole_GUI': rsole1, 'lsole_GUI': lsole1}, index = time)\n", "df6 = pd.DataFrame({'rTAs_Colab': rTAs, 'lTAs_Colab': lTAs, 'rTAs_GUI': rTAs1, 'lTAs_GUI': lTAs1}, index = time)\n", "df1.index.name = df2.index.name = df3.index.name = df4.index.name = df5.index.name = df6.index.name = 'time'\n", "df1.plot(title='rhams/lhams_SO').show()\n", "df2.plot(title='rquad/lquad_SO').show()\n", "df3.plot(title='rgluts/lgluts_SO').show()\n", "df4.plot(title='rgast/lgast_SO').show()\n", "df5.plot(title='rsole/lsole_SO').show()\n", "df6.plot(title='rTAs/lTAs_SO').show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "JffU64L9aaey" }, "source": [ "#@title Find nearest element (gait events, optional)\n", "#- by Fangwei\n", "#This is to find gait events: 0.6s and 1.84s are right side foot strike (i.e. 100% gait cycle). \n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "find_nearest_element(time, start_time, end_time)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "3dutPgjogvC0" }, "source": [ "#@title Test SO gait cycle\n", "# - Fangwei (already simplified - prepare for artical)(gait cycle) \n", "# Don't write them in one line (fig1 = fig2 = fig3 = fig4 = fig5 = fig6 = go.Figure()), we don't want a shadow copy.\n", "\n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "\n", "fig1 = go.Figure()\n", "fig2 = go.Figure()\n", "fig3 = go.Figure()\n", "fig4 = go.Figure()\n", "fig5 = go.Figure()\n", "fig6 = go.Figure()\n", "\n", "four_IDTool(fig1, time, rhams, 'r_Colab', lhams, 'l_Colab', rhams1, 'r_GUI', lhams1, 'l_GUI', 'rhams/lhams', 'fig1.png')\n", "fig1.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig1.show()\n", "#fig1.write_image('fig1.webp')\n", "\n", "four_IDTool(fig2, time, rquad, 'r_Colab', lquad, 'l_Colab', rquad1, 'r_GUI', lquad1, 'l_GUI', 'rquad/lquad', 'fig1.png')\n", "fig2.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig2.show()\n", "#fig2.write_image('fig2.webp')\n", "\n", "four_IDTool(fig3, time, rgluts, 'r_Colab', lgluts, 'l_Colab', rgluts1, 'r_GUI', lgluts1, 'l_GUI', 'rgluts/lgluts', 'fig1.png')\n", "fig3.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig3.show()\n", "#fig3.write_image('fig3.webp')\n", "\n", "four_IDTool(fig4, time, rgast, 'r_Colab', lgast, 'l_Colab', rgast1, 'r_GUI', lgast1, 'l_GUI', 'rgast/lgast', 'fig1.png')\n", "fig4.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig4.show()\n", "#fig4.write_image('fig4.webp')\n", "\n", "four_IDTool(fig5, time, rsole, 'r_Colab', lsole, 'l_Colab', rsole1, 'r_GUI', lsole1, 'l_GUI', 'rsole/lsole', 'fig1.png')\n", "fig5.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig5.show()\n", "#fig5.write_image('fig5.webp')\n", "\n", "four_IDTool(fig6, time, rTAs, 'r_Colab', lTAs, 'l_Colab', rTAs1, 'r_GUI', lTAs1, 'l_GUI', 'rTAs/lTAs', 'fig1.png')\n", "fig6.update_layout(xaxis_range=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time[time.index(min(time, key=lambda x: abs(x-Decimal(start_time))))], time[time.index(min(time, key=lambda x: abs(x-(min(time, key=lambda x: abs(x-Decimal(end_time)))+min(time, key=lambda x: abs(x-Decimal(start_time))))/2)))], time[time.index(min(time, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig6.show()\n", "#fig6.write_image('fig6.webp')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "0cFepT3RevMG" }, "source": [ "#Computed Muslce Control (CMC) in OpenSim\n", "In OpenSim CMC is used to predict muscle activities using a combination of SO and forward dynamics (FD) as shown below:\n", "\n", "

\"Files

\n", "\n", "
\n", "\n", "Files required to perform Computed Muscle Control (CMC).\n", "\n", "

\"Files

\n", "\n" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "QHsjGU6vcQj-" }, "source": [ "#@title Computed Muscle Control (CMC) in Colab 15 mins\n", "from opensim import CMCTool\n", "import time\n", "start_time = time.time()\n", "CMCTool(\"/content/opensim-models/Pipelines/Gait2354_Simbody/subject01_Setup_CMC.xml\").run()\n", "print(f'The execution time of CMCTool is {(time.time() - start_time)} sec')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "qH4nP9CqdyYV" }, "source": [ "##Plotting CMC results & Validation with GUI" ] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "w8NC7Ae2I9VR" }, "source": [ "#@title subject01_walk1_Actuation_force.sto generated by Colab\n", "generate_dict(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsCMC/subject01_walk1_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "xQCf4926JDZv" }, "source": [ "#@title subject01_walk1_Actuation_force.sto generated by OpensimGUI\n", "generate_dict(\"/content/2354_result/ResultsCMC/subject01_walk1_Actuation_force.sto\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "dPhEzbtEbfIg" }, "source": [ "#@title Read data from subject01_walk1_Actuation_force.sto (Colab)\n", "from decimal import Decimal\n", "time, bifemlh_r, bifemsh_r, bifemlh_l, bifemsh_l=[],[],[],[],[]\n", "rect_fem_r,vas_int_r,rect_fem_l,vas_int_l=[],[],[],[]\n", "glut_med1_r,glut_med2_r,glut_med3_r,glut_med1_l,glut_med2_l,glut_med3_l=[],[],[],[],[],[]\n", "med_gas_r, med_gas_l, soleus_r, soleus_l, tib_ant_r, tib_ant_l=[],[],[],[],[],[]\n", "def load_file_2(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time.append(Decimal(m[0]))\n", " # rhams\n", " bifemlh_r.append(Decimal(m[4]))\n", " bifemsh_r.append(Decimal(m[5]))\n", " # lhams\n", " bifemlh_l.append(Decimal(m[28]))\n", " bifemsh_l.append(Decimal(m[29]))\n", " # rquad\n", " rect_fem_r.append(Decimal(m[19]))\n", " vas_int_r.append(Decimal(m[20]))\n", " # lquad\n", " rect_fem_l.append(Decimal(m[43]))\n", " vas_int_l.append(Decimal(m[44]))\n", " # rgluts\n", " glut_med1_r.append(Decimal(m[1]))\n", " glut_med2_r.append(Decimal(m[2]))\n", " glut_med3_r.append(Decimal(m[3]))\n", " # lgluts\n", " glut_med1_l.append(Decimal(m[25]))\n", " glut_med2_l.append(Decimal(m[26]))\n", " glut_med3_l.append(Decimal(m[27]))\n", " # rgast\n", " med_gas_r.append(Decimal(m[21]))\n", " # lgast\n", " med_gas_l.append(Decimal(m[45]))\n", " # rsole\n", " soleus_r.append(Decimal(m[22]))\n", " # lsole\n", " soleus_l.append(Decimal(m[46]))\n", " # rTAs\n", " tib_ant_r.append(Decimal(m[24]))\n", " # lTAs\n", " tib_ant_l.append(Decimal(m[48]))\n", " c.flush()\n", " c.close()\n", "load_file_2(\"/content/opensim-models/Pipelines/Gait2354_Simbody/ResultsCMC/subject01_walk1_Actuation_force.sto\")\n", "rhams,lhams,rquad,lquad,rgluts,lgluts,rgast,lgast,rsole,lsole,rTAs,lTAs=[],[],[],[],[],[],[],[],[],[],[],[]\n", "mid_r,mid_l=[],[]\n", "rhams[:] = [x + y for x, y in zip(bifemlh_r, bifemsh_r)]\n", "lhams[:] = [x + y for x, y in zip(bifemlh_l, bifemsh_l)]\n", "rquad[:] = [x + y for x, y in zip(rect_fem_r,vas_int_r)]\n", "lquad[:] = [x + y for x, y in zip(rect_fem_l,vas_int_l)]\n", "mid_r[:] = [x + y for x, y in zip(glut_med1_r, glut_med2_r)]\n", "mid_l[:] = [x + y for x, y in zip(glut_med1_l, glut_med2_l)]\n", "rgluts[:] = [x + y for x, y in zip(mid_r, glut_med3_r)]\n", "lgluts[:] = [x + y for x, y in zip(mid_l, glut_med3_l)]\n", "rgast = med_gas_r\n", "lgast = med_gas_l\n", "rsole = soleus_r\n", "lsole = soleus_l\n", "rTAs = tib_ant_r\n", "lTAs = tib_ant_l\n", "print('finished')" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "5KBfIO-4j8_n" }, "source": [ "#@title Read data from subject01_walk1_Actuation_force.sto (OpenSimGUI)\n", "from decimal import Decimal\n", "time1, bifemlh_r1, bifemsh_r1, bifemlh_l1, bifemsh_l1=[],[],[],[],[]\n", "rect_fem_r1,vas_int_r1,rect_fem_l1,vas_int_l1=[],[],[],[]\n", "glut_med1_r1,glut_med2_r1,glut_med3_r1,glut_med1_l1,glut_med2_l1,glut_med3_l1=[],[],[],[],[],[]\n", "med_gas_r1, med_gas_l1, soleus_r1, soleus_l1, tib_ant_r1, tib_ant_l1=[],[],[],[],[],[]\n", "def load_file_2(file_name):\n", " c = open(file_name)\n", " for x in range(judge_line_number(file_name)):\n", " next(c)\n", " for i in c.readlines():\n", " m = i.strip('\\n').split()\n", " time1.append(Decimal(m[0]))\n", " # rhams\n", " bifemlh_r1.append(Decimal(m[4]))\n", " bifemsh_r1.append(Decimal(m[5]))\n", " # lhams\n", " bifemlh_l1.append(Decimal(m[28]))\n", " bifemsh_l1.append(Decimal(m[29]))\n", " # rquad\n", " rect_fem_r1.append(Decimal(m[19]))\n", " vas_int_r1.append(Decimal(m[20]))\n", " # lquad\n", " rect_fem_l1.append(Decimal(m[43]))\n", " vas_int_l1.append(Decimal(m[44]))\n", " # rgluts\n", " glut_med1_r1.append(Decimal(m[1]))\n", " glut_med2_r1.append(Decimal(m[2]))\n", " glut_med3_r1.append(Decimal(m[3]))\n", " # lgluts\n", " glut_med1_l1.append(Decimal(m[25]))\n", " glut_med2_l1.append(Decimal(m[26]))\n", " glut_med3_l1.append(Decimal(m[27]))\n", " # rgast\n", " med_gas_r1.append(Decimal(m[21]))\n", " # lgast\n", " med_gas_l1.append(Decimal(m[45]))\n", " # rsole\n", " soleus_r1.append(Decimal(m[22]))\n", " # lsole\n", " soleus_l1.append(Decimal(m[46]))\n", " # rTAs\n", " tib_ant_r1.append(Decimal(m[24]))\n", " # lTAs\n", " tib_ant_l1.append(Decimal(m[48]))\n", " c.flush()\n", " c.close()\n", "load_file_2(\"/content/2354_result/ResultsCMC/subject01_walk1_Actuation_force.sto\")\n", "rhams1,lhams1,rquad1,lquad1,rgluts1,lgluts1,rgast1,lgast1,rsole1,lsole1,rTAs1,lTAs1=[],[],[],[],[],[],[],[],[],[],[],[]\n", "mid_r,mid_l=[],[]\n", "rhams1[:] = [x + y for x, y in zip(bifemlh_r1, bifemsh_r1)]\n", "lhams1[:] = [x + y for x, y in zip(bifemlh_l1, bifemsh_l1)]\n", "rquad1[:] = [x + y for x, y in zip(rect_fem_r1,vas_int_r1)]\n", "lquad1[:] = [x + y for x, y in zip(rect_fem_l1,vas_int_l1)]\n", "mid_r[:] = [x + y for x, y in zip(glut_med1_r1, glut_med2_r1)]\n", "mid_l[:] = [x + y for x, y in zip(glut_med1_l1, glut_med2_l1)]\n", "rgluts1[:] = [x + y for x, y in zip(mid_r, glut_med3_r1)]\n", "lgluts1[:] = [x + y for x, y in zip(mid_l, glut_med3_l1)]\n", "rgast1 = med_gas_r1\n", "lgast1 = med_gas_l1\n", "rsole1 = soleus_r1\n", "lsole1 = soleus_l1\n", "rTAs1 = tib_ant_r1\n", "lTAs1 = tib_ant_l1\n", "print('finished')" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "459NjzUwyfDO", "cellView": "form" }, "source": [ "#@title Offset x-axis value\n", "def xxx(array, qq):\n", " for i in range(qq):\n", " array.pop(0)\n", " array.pop()\n", "\n", "def xxx1(array,qq1,qq2):\n", " for i in range(qq1):\n", " array.pop(0)\n", " for i in range(qq2):\n", " array.pop() \n", "\n", "if ( abs(len(time)-len(time1)) % 2 == 0 ):\n", " qq = int(abs(len(time)-len(time1)) / 2)\n", " xxx(rhams,qq)\n", " xxx(lhams,qq)\n", " xxx(rquad,qq)\n", " xxx(lquad,qq)\n", " xxx(rgluts,qq)\n", " xxx(lgluts,qq)\n", " xxx(rgast,qq)\n", " xxx(lgast,qq)\n", " xxx(rsole,qq)\n", " xxx(lsole,qq)\n", " xxx(rTAs,qq)\n", " xxx(lTAs,qq)\n", "else:\n", " qq1 = int((abs(len(time)-len(time1)) + 1) / 2)\n", " qq2 = int((abs(len(time)-len(time1)) - 1) / 2)\n", " xxx1(rhams,qq1,qq2)\n", " xxx1(lhams,qq1,qq2)\n", " xxx1(rquad,qq1,qq2)\n", " xxx1(lquad,qq1,qq2)\n", " xxx1(rgluts,qq1,qq2)\n", " xxx1(lgluts,qq1,qq2)\n", " xxx1(rgast,qq1,qq2)\n", " xxx1(lgast,qq1,qq2)\n", " xxx1(rsole,qq1,qq2)\n", " xxx1(lsole,qq1,qq2)\n", " xxx1(rTAs,qq1,qq2)\n", " xxx1(lTAs,qq1,qq2)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "nNvrI6gJQsum", "cellView": "form" }, "source": [ "#@title x-axis shifting\n", "time1[:]=[x-time1[0] for x in time1]" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "aGlyP4JdvkIi", "cellView": "form" }, "source": [ "#@title Test CMCTool panda GUI quick draw full time (0 to 2.5s)\n", "# (code is very simple and execution is short) better for a short demo - Fangwei\n", "df1 = pd.DataFrame({'rhams_Colab': rhams, 'lhams_Colab': lhams, 'rhams_GUI': rhams1, 'lhams_GUI': lhams1}, index = time1)\n", "df2 = pd.DataFrame({'rquad_Colab': rquad, 'lquad_Colab': lquad, 'rquad_GUI': rquad1, 'lquad_GUI': lquad1}, index = time1)\n", "df3 = pd.DataFrame({'rgluts_Colab': rgluts, 'lgluts_Colab': lgluts, 'rgluts_GUI': rgluts1, 'lgluts_GUI': lgluts1}, index = time1)\n", "df4 = pd.DataFrame({'rgast_Colab': rgast, 'lgast_Colab': lgast, 'rgast_GUI': rgast1, 'lgast_GUI': lgast1}, index = time1)\n", "df5 = pd.DataFrame({'rsole_Colab': rsole, 'lsole_Colab': lsole, 'rsole_GUI': rsole1, 'lsole_GUI': lsole1}, index = time1)\n", "df6 = pd.DataFrame({'rTAs_Colab': rTAs, 'lTAs_Colab': lTAs, 'rTAs_GUI': rTAs1, 'lTAs_GUI': lTAs1}, index = time1)\n", "df1.index.name = df2.index.name = df3.index.name = df4.index.name = df5.index.name = df6.index.name = 'time'\n", "df1.plot(title='rhams/lhams_CMC').show()\n", "df2.plot(title='rquad/lquad_CMC').show()\n", "df3.plot(title='rgluts/lgluts_CMC').show()\n", "df4.plot(title='rgast/lgast_CMC').show()\n", "df5.plot(title='rsole/lsole_CMC').show()\n", "df6.plot(title='rTAs/lTAs_CMC').show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "dHvn_OPhjZ-E" }, "source": [ "#@title Find nearest element (gait events, optional)\n", "#- by Fangwei\n", "#This is to find gait events: 0.6s and 1.84s are right side foot strike (i.e. 100% gait cycle). \n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "find_nearest_element(time1, start_time, end_time)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "cellView": "form", "id": "iEv5X6jVSiue" }, "source": [ "#@title Test CMCTool gait cycle\n", "# - Fangwei (already simplified - prepare for article)(gait cycle) \n", "# Don't write them in one line (fig1 = fig2 = fig3 = fig4 = fig5 = fig6 = go.Figure()), we don't want a shadow copy.\n", "\n", "start_time = 0.600 #@param {type:\"number\"}\n", "end_time = 1.840 #@param {type:\"number\"}\n", "\n", "fig1 = go.Figure()\n", "fig2 = go.Figure()\n", "fig3 = go.Figure()\n", "fig4 = go.Figure()\n", "fig5 = go.Figure()\n", "fig6 = go.Figure()\n", "\n", "four_IDTool(fig1, time1, rhams, 'r_Colab', lhams, 'l_Colab', rhams1, 'r_GUI', lhams1, 'l_GUI', 'rhams/lhams', 'fig1.png')\n", "fig1.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig1.show()\n", "#fig1.write_image('fig1.webp')\n", "fig1.update_layout(title=\"\",yaxis_title='')\n", "fig1.write_image('fig1.svg')\n", "\n", "four_IDTool(fig2, time1, rquad, 'r_Colab', lquad, 'l_Colab', rquad1, 'r_GUI', lquad1, 'l_GUI', 'rquad/lquad', 'fig1.png')\n", "fig2.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig2.show()\n", "#fig2.write_image('fig2.webp')\n", "fig2.update_layout(title=\"\",yaxis_title='')\n", "fig2.write_image('fig2.svg')\n", "\n", "four_IDTool(fig3, time1, rgluts, 'r_Colab', lgluts, 'l_Colab', rgluts1, 'r_GUI', lgluts1, 'l_GUI', 'rgluts/lgluts', 'fig1.png')\n", "fig3.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig3.show()\n", "#fig3.write_image('fig3.webp')\n", "fig3.update_layout(title=\"\",yaxis_title='')\n", "fig3.write_image('fig3.svg')\n", "\n", "four_IDTool(fig4, time1, rgast, 'r_Colab', lgast, 'l_Colab', rgast1, 'r_GUI', lgast1, 'l_GUI', 'rgast/lgast', 'fig1.png')\n", "fig4.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig4.show()\n", "#fig4.write_image('fig4.webp')\n", "fig4.update_layout(title=\"\",yaxis_title='')\n", "fig4.write_image('fig4.svg')\n", "\n", "four_IDTool(fig5, time1, rsole, 'r_Colab', lsole, 'l_Colab', rsole1, 'r_GUI', lsole1, 'l_GUI', 'rsole/lsole', 'fig1.png')\n", "fig5.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig5.show()\n", "#fig5.write_image('fig5.webp')\n", "fig5.update_layout(title=\"\",yaxis_title='')\n", "fig5.write_image('fig5.svg')\n", "\n", "four_IDTool(fig6, time1, rTAs, 'r_Colab', lTAs, 'l_Colab', rTAs1, 'r_GUI', lTAs1, 'l_GUI', 'rTAs/lTAs', 'fig1.png')\n", "fig6.update_layout(xaxis_range=[time[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]], xaxis=dict(\n", " tickmode='array',\n", " tickvals=[time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(start_time))))], time1[time1.index(min(time1, key=lambda x: abs(x-(min(time1, key=lambda x: abs(x-Decimal(end_time)))+min(time1, key=lambda x: abs(x-Decimal(start_time))))/2)))], time1[time1.index(min(time1, key=lambda x: abs(x-Decimal(end_time))))]],\n", " ticktext=['0','50','100']\n", "),xaxis_title='Gait cycle (%)')\n", "fig6.show()\n", "#fig6.write_image('fig6.webp')\n", "fig6.update_layout(title=\"\",yaxis_title='')\n", "fig6.write_image('fig6.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "FnsIToeypRVA", "cellView": "form" }, "source": [ "#@title save subplot (CMC) to svg\n", "# - Fangwei 04/06/2021\n", "import svgutils.transform as sg\n", "\n", "fig = sg.SVGFigure(\"16cm\", \"16cm\")\n", "fig1 = sg.fromfile('/content/fig1.svg')\n", "fig2 = sg.fromfile('/content/fig2.svg')\n", "fig3 = sg.fromfile('/content/fig3.svg')\n", "fig4 = sg.fromfile('/content/fig4.svg')\n", "fig5 = sg.fromfile('/content/fig5.svg')\n", "fig6 = sg.fromfile('/content/fig6.svg')\n", "\n", "plot1 = fig1.getroot()\n", "plot2 = fig2.getroot()\n", "plot3 = fig3.getroot()\n", "plot4 = fig4.getroot()\n", "plot5 = fig5.getroot()\n", "plot6 = fig6.getroot()\n", "\n", "plot1.moveto(100, 0, 1, None)\n", "plot2.moveto(800, 0, 1, None)\n", "plot3.moveto(1500, 0, 1, None)\n", "plot4.moveto(100, 500, 1, None)\n", "plot5.moveto(800, 500, 1, None)\n", "plot6.moveto(1500, 500, 1, None)\n", "\n", "txt1 = sg.TextElement(400, 30, \"Hamstrings\", size=30, weight=\"bold\")\n", "txt2 = sg.TextElement(1130, 30, \"Quadriceps\", size=30, weight=\"bold\")\n", "txt3 = sg.TextElement(1800, 30, \"Gluteus maximus\", size=30, weight=\"bold\")\n", "txt4 = sg.TextElement(-300, 20, \"Muscle forces(N)\", size=30, weight=\"bold\")\n", "txt4.rotate(-90, 100)\n", "txt5 = sg.TextElement(-800, 20, \"Muscle forces(N)\", size=30, weight=\"bold\")\n", "txt5.rotate(-90, 100)\n", "txt6 = sg.TextElement(380, 550, \"Gastrocnemius\", size=30, weight=\"bold\")\n", "txt7 = sg.TextElement(1150, 550, \"Solues\", size=30, weight=\"bold\")\n", "txt8 = sg.TextElement(1800, 550, \"Tibiali Anterior\", size=30, weight=\"bold\")\n", "\n", "fig.append([plot1, plot2, plot3, plot4, plot5, plot6])\n", "fig.append([txt1, txt2, txt3, txt4, txt5, txt6, txt7, txt8])\n", "fig.set_size([\"58cm\", \"26.5cm\"])\n", "\n", "fig.save(\"/content/fig_subplot_cmc.svg\")\n", "print('Finished, please check fig_subplot_cmc.svg')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "A759Wsgzcyrw" }, "source": [ "# Future work in OpenColab\n", "\n", "1. Automate the whole process and OpenSim workflow from c3d to muscle/joint function (our next study)\n", "2. Installation of other packages of OpenSim e.g. MOCO, OpenSense, etc (some has been done but may need to easily implemented in Colab)\n", "3. 3D Visualization similar to GUI and a Colab GUI for the whole process\n", "4. Link to TF and how to run the models in OpenSim and utilize the ML/AL in Colab (quite easy when #1 done)\n", "5. How to speed up the whole process? Now if Internet is down, we lose all the functions and need to reinstall OpenSim and other packages. [Colab issue]\n", "6. Can we install OpenSim faster? It takes now up to 6 min each time.\n", "7. Build a win-64 conda package to avoid manual installation.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "35AfnhqgXX5G" }, "source": [ "# References\n", "\n", "1. OpenSim Tutorials, https://simtk-confluence.stanford.edu/display/OpenSim/Examples+and+Tutorials\n", "\n", "2. Delp, S.L., Loan, J.P., Hoy, M.G., Zajac, F.E., Topp E.L., Rosen, J.M. An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures. IEEE Transactions on Biomedical Engineering, vol. 37, pp. 757-767, 1990.\n", "\n", "2. Anderson, F.C., Pandy, M.G. A dynamic optimization solution for vertical jumping in three dimensions. Computer Methods in Biomechanical and Biomedical Engineering, vol. 2, pp. 201-231, 1999.\n", "\n", "3. Kuo, A.D. A least squares estimation approach to improving the precision of inverse dynamics computations, Journal of Biomechanical Engineering, vol. 120, pp. 148-159, 1998.\n", "\n", "4. Winter, D.A. Biomechanics and Motor Control of Human Movement, Wiley and Sons, pp. 77-79, 1990.\n", "\n", "5. Thelen, D.G., Anderson, F.C. Using computed muscle control to generate forward dynamic simulations of human walking from experimental data, Journal of Biomechanics, vol. 39, pp. 1107-1115, 2006.\n", "\n", "6. John, C.T., Anderson, F.C., Guendelman, E., Arnold, A.S., Delp, S.L. An algorithm for generating muscle-actuated simulations of long-duration movements, Biomedical Computation at Stanford (BCATS) Symposium, Stanford University, 21 October 2006, Poster Presentation.\n", "\n", "7. Delp, S.L., Anderson, F.C., Arnold, A.S., Loan, P., Habib, A., John, C.T., Guendelman, E., Thelen, D.G. OpenSim: Open-source software to create and analyze dynamic simulations of movement. IEEE Transactions on Biomedical Engineering, vol. 55, pp. 1940-1950, 2007.\n", "\n", "8. Chand T. John, Frank C. Anderson, Jill S. Higginson & Scott L. Delp (2012): Stabilisation of walking by intrinsic muscle properties revealed in a three-dimensional muscle-driven simulation, Computer Methods in Biomechanics and Biomedical Engineering.\n", "\n", "9. Mokhtarzadeh H, Anderson DE, Allaire BT, Bouxsein ML. Patterns of load‐to‐strength ratios along the spine in a population‐based cohort to evaluate the contribution of spinal loading to vertebral fractures. J Bone Miner Res 2020. \n", "\n", "10. Mokhtarzadeh H, Yeow CH, Hong Goh JC, Oetomo D, Malekipour F, Lee PV-S. Contributions of the Soleus and Gastrocnemius muscles to the anterior cruciate ligament loading during single-leg landing. J Biomech 2013;46:1913–20. https://doi.org/10.1016/j.jbiomech.2013.04.010.\n", "\n", "11. Dorn TW. Computational modeling of lower-limb muscle function in human running. The University of Melbourne, 2011. \n", "\n", "12. Winter DA. Biomechanics and motor control of human movement. John Wiley & Sons, 2009. \n", "\n", "13. Mokhtarzadeh Hossein, Fangwei Jiang; Andy Shengzhe Zhao; Fatemeh Malekipour. *OpenColab Project*. https://simtk.org/projects/opencolab\n", "\n", "14. Mokhtarzadeh, H. (2013). Anterior cruciate ligament injury mechanism during impact load. PhD thesis, Department of Mechanical Engineering, Melbourne School of Engineering, The University of Melbourne.\n", "\n", "15. Bartosz Telenczuk, svgutils 0.3.4. https://pypi.org/project/svgutils/\n", "\n", "\n", "\n", "\n", "\n" ] } ] }