BTK Matlab Wrapper  0.3dev.0
Matlab bindings for the Biomechanical ToolKit library (BTK)
Tutorial

The following list gives you the steps to use BTK in Matlab.


Build Matlab support from the sources of BTK

This section is only necessary if you downloaded the sources of BTK. If you use a binary version, go directly to the tutorial Getting Started.

To simplify the way to build the Matlab support, the tool EasyInstallWithRedistributableMatlabToolbox.jar is proposed in the root of the sources.

To use this tool, we need to download the following software:

When these software are installed, you only need to double click on the file EasyInstallWithRedistributableMatlabToolbox.jar. This tool will configure the sources, compile them and install the toolbox btk for Matlab.

Note: If your OS doesn't have Java installed (and you don't want or cannot install it), then you have again the possibility to use the tool EasyInstallWithRedistributableMatlabToolbox. For that, you have to go in the folder Batch and double clic on the file corresponding to your OS:

Getting Started

First of all, you need to add the path of the BTK toolbox in the list of the directories loaded by Matlab. Depending your OS, you will have to add the following path:

To be sure that BTK is loaded in Matlab, type the command help btk in the Matlab's command window. You will see the the main documentation for the BTK toolbox. If it not the case, be sure of the path that you gave to Matlab to register the BTK toolbox.

Now, that BTK is recognized in Matlab, it comes easy to open your acquisition files. You have only to use the function btkReadAcquisition. Based on the content of your file, this function tries to open the given file and return a handle which will help you to access to the content of the acquisition.

Note: Due to the use of C++ object hidden by the handle, Matlab doesn't know the memory allocated for each acquisition. In lots of case, you don't need to manage the memory used by BTK in Matlab (for example, if you want to check the content of an acquistion or read a dozen of files to analyze gait data) as the command clear all will clear all the memory. But if you process hundred of files, then you have to use the command btkDeleteAcquisition. This command free the memory associated with the acquisition. This function can also be understood as a 'close' function.

Acquisition data explained

This section presents you some possibilities to access or modify data stored in an acquisition.

An acquisition contains points (parameter with 3 components), analog channels (parameter with 1 component), events and metadata.

Points

A point in BTK is a generic type which contains 3D data. Then markers, forces, moments, angles, powers are regrouped in this category. For example, to extract markers' trajectories you need to use the function btkGetMarkers like this:

acq = btkReadAcquisition('myAcquisitionFile.c3d');
markers = btkGetMarkers(acq);
% markers is a scructure containing the 3D trajectory of each markers.
% For example, the movement of the center of the knee (KNEE) can be plotted by using:
plot(markers.KNEE);

If you load a C3D file already processed which contains kinematics or inverse dynamic (for example, by using PluginGait for the Vicon users), it is very simple to extract angles, forces, moments, powers, etc. You only have to use the corresponding function. All these functions return a structure where each field correspond to the label of the requested points.

angles = btkGetAngles(acq);
forces = btkGetForces(acq);
moments = btkGetMoments(acq);
powers = btkGetPowers(acq);
% ...

You can also want to extract only few points from an acquisition and return then as a matrix. This can be realized by using a little script based on BTK functions.

%
filename = 'myAcquisition.c3d';
pointsLabel = {'C7', 'CLAV', 'CenterOfMass','KneeJointForce'};
%
acq = btkReadAcquisition(filename);
num = length(pointsLabel);
points = btkGetPoints(acq);
extractedPoints = zeros(btkGetPointFrameNumber(acq), num*3);
for i = 1:num
extractedPoints(:, 1+(i-1)*3:i*3) = points.(pointsLabel{i});
end

Note: With BTK, it is not necessary to use the metadata (named group and parameter in a C3D file) TRIAL:ACTUAL_START_FIELD and TRIAL:ACTUAL_END_FIELD to determine what is the index of the first frame and the last frame. There is already some functions for this. The function btkGetFirstFrame returns the index of the first frame. The last frame will be determined by the function btkGetLastFrame which use the formula: btkGetPointFrameNumber(handle) - btkGetFirstFrame(handle) + 1.

Note: The frames extracted from the C3D correspond to the number of frames between ACTUAL_START_FIELD and ACTUAL_END_FIELD. If your acquisition start at the frame 511 and finish at the frame 8450, and you want the 1450th you only have to type extractedPoints(1450,:). However, if you want the frame #1450, then you have to substract the index of the first frame (here 511): extractedPoints(1450-btkGetFirstFrame(acq)+1,:).

If you can access to the data, you can also edit them and finally write them in a new acquisition file. It could be useful if you have your own algorithm to compute kinematic and inverse dynamic (or to compute the mean of each parameter over a gait cycle, or ...) and would like to save them in a C3D file to be used in another software like Polygon or Mokka.

% We can use a new acquisition for that
acq = btkNewAcquisition(20, 2000); % 20 points and 2000 frames
% Acquisition frequency was 200 Hz
btkSetFrequency(acq, 200);
% Values setting
for i = 1:btkGetPointNumber(acq)
btkSetPointLabel(acq, i, labels{i}):
btkSetPoint(acq, i, values(i,:,:)): % values: array of 20 by 3 by 2000
end
btkWriteAcquisition(acq, 'myNewAcquisitionFile.c3d');

The output is then a new C3D file containing the standard parameters required to be used in others softwares.

Analog channels

The analog channels contains all the 1D measures. All the sensors measuring voltage, angular velocity, acceleration, etc., should be stored in analog channels. From the hardware point of view, it corresponds to the sensors plugged into the analog to digital converter (ADC) card.

To extract the data related to the analog channels, the simplest is the command btkGetAnalogs. As the function btkGetPoints, it returns a structure where each field correspond to a analog channel. If you use a analog sample rate greater than the video frequency, then the number of analog frames is greater than the number of point frames. For example, using the functions btkGetPointFrequency and btkGetAnalogFrequency, you find that the video and analog frequencies are 100 Hz and 1000 Hz respectively. So it means that for 1 video frame you have 10 samples for each analog channels. This can be confirmed by the function btkGetAnalogSampleNumberPerFrame. So, for an acquisition of 10 seconds, you will have 1000 frames for the points and 10000 frames for the analog channels.

An easy way to extract the analog frame corresponding to the video frame can be realized with the following code:

% acq = btkReadAcquisition ...
analogs = btkGetAnalogs(acq);
analogsDownsampled = [];
labels = fieldnames(analogs);
for i = 1:btkGetAnalogNumber(acq)
analogsDownsampled.(labels{i}) = analogs.(labels{i})(1:ratio:end);
% ...
end

An analog channel contains also of a numerical offset and a scale factor used for the analog to digital conversion or vice-versa. You can find them in the second output of the function btkGetAnalogs. These parameters are important if you want to add analog channels in your acquisition and then save it in a file. These informations should be found in your hardware system or in the configuration of your acquisition software.

Events

To extract events, you only have to use the function btkGetEvents which sorts the events by label and by time.

events = btkGetEvents(acq)

events = 

     Left_Foot_Strike: [5.6300 6.6900 7.8100]
    Right_Foot_Strike: [6.1600 7.3000 8.3000]
        Left_Foot_Off: [6.2500 7.4500]
       Right_Foot_Off: [6.9400 7.9800] 

Metadata

A metadata is a generic container to store the acquisition configuration. Or said differently, a metadata contains every informations which cannot be set into markers' trajectories nor analog channels' measures, nor events. For example, if you want to include subject's informations (sex, weight, height, ...) or force platform's configuration (corner's coordinates, analog channel used, ...), then the metadata are the right place to do this. In the C3D file format, the metadata are known as groups and parameters.

To access to the acquisition's metadata, you have to use the function btkGetMetaData which will all the content of the metadata from its root. The output of the function btkGetMetaData is a structure which gives you a tree with children and values. For example, if you want to access to the value of the metadata SUBJECTS:WEIGHT, you can use the following code:

% acq = btkReadAcquisition ...
md = btkGetMetaData(acq);
weight = md.children.SUBJECTS.children.WEIGHT.info.values;

If you are not sure that the metadata SUBJECTS:WEIGHT exists you can use the function btkFindMetaData.

% acq = btkReadAcquisition ...
md = btkFindMetaData(acq, 'SUBJECTS', 'WEIGHT');
if (md ~= 0) % exists
weight = md.info.values; % the returned structure is directly the content of the metadata SUBJECTS:WEIGHT.
end

Add force platform type 2 (AMTI or equivalent) into an acquisition

In some case, your acquisition setup cannot record all the data by the same system and then, you have several files for one trial. This section explain you how to add a force platform into an acquistion.

A force platform is represented in BTK as a set of analog channels combined with a configuration stored in the metadata FORCE_PLATFORM and its children. This configuration is the same than the one proposed in the C3D file format (http://www.c3d.org/HTML/theforceplatformgroup.htm). To easily integrate a force platfrom, it exists the Matlab function btkAppendForcePlatformType2. This function appends an AMTI (or equivalent with FxFyFzMxMyMz measurements) force platform by adding required analog channels and the creation (or update) of its configuration.

% Create with an acquisition with 19 markers, 1000 frames, no analog channel but sampled 10 times more than video frame.
acq = btkNewAcqusition(19,1000,0,10);
% Set the frequency: Markers: 100Hz, Analog channels: 1000Hz (due to the use of the factor 10 in the function btkNewAcquisition)
% Set 3D coordinates for markers
% - markersV is a matrix wit the dimensions (19*3,1000)
btkSetMarkersValues(acq, markersV);
% Append a force platform type 2:
% - forcesV is a matrix (3,10000) with reaction forces ;
% - momentsV is a matrix (3,10000) with reaction moments ;
% - corners is a matrix (4,3) with the coordinates of the 4 corners expressed in the global (laboratory) frame ;
% - origin is a vector (1,3) containing the location of the sensor related to the center of the working surface and expressed in the global frame ;
% - 1 is a flag to say that the reaction forces and moments are measured in the force platform (local) frame.
btkAppendForcePlatformType2(acq, forcesV, momentsV, corners, origin, 1);
% Write the new acquisition
btkWriteAcquisition(acq, 'myAcquisitionWithFP.c3d');

Processing data

One of the goal of BTK is also to give some tools to compute standard biomechanical parameters (kinematic, dynamic inverse, ...). All the existing tools are stored in the BasicFilters section in the help. For example, to determine the ground reaction forces (wrenches) from an acquisition containing a force platform, you only have to use the function btkGetGroundReactionWrenches. This function take into account the type of force platform (AMTI, Kistler, ...) and compute the group reaction wrench (force, moment and position). The point of application of each wrench is calculated from the formula of Shimba (1984) which also was described in Zatsiorsky (2002). Contrary to the center of pressure (COP) the point of wrench application (PWA) is calculated by using also the horizontal moments. The result is a structure with a number of elements equal to the number of force platforms.

% Import acquisition
acq = btkReadAcquisition('myAcquisition.c3d');
% Extract analogs values as a matrix.
% Filter data
% The parameters must be multiplied due to the use of filtfilt.
% The final order is divided by 2 and the cut-off frequency is multiplied by 1.5162
% These values are computed from the article:
% "Design and responses of Butterworth and critically damped digital filters",
% Robertson & Dowling, Journal of Electromyography and Kinesiology, 2003.
% Warning: the factor 1.5162 is dependant of the filter order and must be
% adapted in consequence for other order using the formula given in the paper.
[b,a] = butter(4/2, 2 * 30 * 1.5162 / btkGetAnalogFrequency(acq));
av = filtfilt(b,a, av);
% Compute ground reaction wrenches with a threshold of 3 newtons
figure; plot(grws(1).P) % Position
figure; plot(grws(1).F) % Force
figure; plot(grws(1).M) % Moment

Batch processing

Even if BTK is well adapted to open hundred of acquisition files and process them, there is a limitation in the mechanism used to manage memory between the C++ objects and Matlab. If you open too many acquisitions without releasing the memory, Matlab can send you an "Out of Memory" error and may crash... (see FAQ). One way to correct this problem is to use the function btkDeleteAcquisition (or btkCloseAcquisiton) which will release the memory associated to the C++ object.

% acquisitionsFilename: Large list of filenames
for 1:length(acquisitionsFilename)
acq = btkReadAcquistion(acquisitionsFilename{i});
% ... All of your code processing the acquisition data
btkCloseAcquisition(acq); % Release the memory of the C++ object associated with the handle 'acq'
end

References

  1. Shimba T., An estimation of center of gravity from force platform data, Journal of Biomechanics, 1984, 17(1), 53–60.
  2. Zatsiorsky V.M., Kinetics of Human Motion, Human Kinetics Publishers, Champaign, IL, 2002.