Open C3D files

Marcos Duarte

Text and binary files
There are two kinds of computer files: text and binary files. Text files are structured as a sequence of lines of electronic text. The most common formats of a text file are ASCII (with 128 ($2^7$) different characters and UTF-8 (which includes non-English characters). A binary file is a file that is not structured as a text file. Because in fact everything in a computer is stored in binary format (a sequence of zeros and ones), text files are binary files that store text codes.

To open and read a text file is simple and straightforward. A text file doesn't need additional information to be read, and can be openned by any text-processing software. This is not the case of a binary file, we need to have extra information about how the file is structured to be able to read it. However, binary files can store more information per file size than text files and we can read and write binary files faster than text files. This is one of the reasons why software developers would choose a binary format.

C3D format

The C3D format is a public domain, binary file format that has been used in Biomechanics, Animation and Gait Analysis laboratories to record synchronized 3D and analog data since the mid 1980's. It is supported by all 3D major Motion Capture System manufacturers, as well as other companies in the Biomechanics, Motion Capture and Animation Industries (http://www.c3d.org/).

There is a very detailed technical implementation manual of the C3D file format.

The C3D file format has three basic components:

  • Data: at this level the C3D file is simply a binary file that stores raw 3D and analog information.
  • Standard parameters: default information about the raw 3D and analog data that is required to access the data.
  • Custom parameters: information specific to a particular manufacturers’ software application or test subject.

Regarding to what is useful to the analysis, a C3D file normally has four types of information:

  • 3D point: 3D coordinates of markers or any related biomechanical quantities such as angles, joint forces and moments, etc.
  • Analaog: analog data acquired with an A/D converter.
  • Event: specific frames as events of the acquisition.
  • Metadata: information about the subject, the system configuration, etc.

All C3D files contain a minimum of three sections of information:

  • A single, 512 byte, header section.
  • A parameter section consisting of one, or more, 512-byte blocks.
  • 3D point/analog data section consisting of one, or more, 512-byte blocks.

To demonstrate the process of reading a binary file, let's read only part of the header of a C3D file:

In [1]:
from __future__ import division, print_function
from struct import unpack       # convert C structs represented as Python strings
from cStringIO import StringIO  # reads and writes a string buffer

def getFloat(floatStr,processor):
    "16-bit float string to number conversion"
    if processor > 1: #DEC or SGI
        floatNumber = unpack('f',floatStr[2:] + floatStr[0:2])[0]/4
    else:
        floatNumber = unpack('f',floatStr)[0]
    return floatNumber   

filename = './../data/Gait.c3d'  # data from sample03.zip @ http://www.c3d.org/sampledata.html
fid = open(filename, 'rb')       # open file for reading in binary format
#Header section:
bytes = fid.read(512)
buf = StringIO(bytes)
firstBlockParameterSection, fmt = unpack('BB', buf.read(2))
if fmt != 80:
    print('This file is not a valid C3D file.')
    fid.close()
# First block of parameter section:
firstBlockByte = 512*(firstBlockParameterSection - 1) + 2
fid.seek(firstBlockByte)
nparameter_blocks, processor = unpack('BB', fid.read(2))
processor = processor - 83
processors = ['unknown', 'Intel', 'DEC', 'SGI']
#Back to the header section:
n3Dpoints, = unpack('H', buf.read(2))
nTotalAnalogDataPerFrame, = unpack('H', buf.read(2))
nFirstFrame3D, = unpack('H', buf.read(2))
nLastFrame3D, = unpack('H', buf.read(2))
maxInterpGap, = unpack('H', buf.read(2))
scaleFactor = getFloat(buf.read(4), processor)
dataStartBlock, = unpack('H', buf.read(2))
nAnalogDataPerFrame, = unpack('H', buf.read(2))
frameRate = getFloat(buf.read(4), processor)

print('File "%s":' %filename)
print('Number of 3D points: %d' % n3Dpoints)
print('Number of frames: %d' %(nLastFrame3D - nFirstFrame3D + 1))
print('Frame rate: %.2f' % frameRate)
print('Analog data per frame: %d' % nAnalogDataPerFrame)   
File "data/Gait.c3d":
Number of 3D points: 33
Number of frames: 487
Frame rate: 100.00
Analog data per frame: 10

To read a complete C3D file, and fast enough, is not simple. Fortunately, there is a Python package to work with C3D files: the BTK library. The same developer of BTK, Arnaud Barré, also wrote another nice software: Mokka, an 'open-source and cross-platform software to easily analyze biomechanical data'.

See this page on how to install the BTK library for Python.

Let's use BTK to read a C3D file.
There is a different version of BTK for each OS; here is a workaround to import the right BTK library according to your OS:

In [2]:
import platform, sys

if platform.system() == 'Windows':
    if sys.maxsize > 2**32:
        sys.path.insert(1, r'./../functions/btk/win64')
    else:
        sys.path.insert(1, r'./../functions/btk/win32')
elif platform.system() == 'Linux':
    if sys.maxsize > 2**32:
        sys.path.insert(1, r'./../functions/btk/linux64')
    else:
        sys.path.insert(1, r'./../functions/btk/linux32')
elif platform.system() == 'Darwin':
    sys.path.insert(1, r'./../functions/btk/mac')

The BTK files above were taken from https://code.google.com/p/b-tk/wiki/PythonBinaries with the exception of the files for the Mac OS. I compiled these files for the latest Mac OS version (10.9) and Python 2.7.6. Use the Mac files from the BTK website if you have older Mac OS and Python versions.

With the BTK files and the path to them in the Python path, to use them it's easy:

In [3]:
import btk

Read a C3D file:

In [4]:
reader = btk.btkAcquisitionFileReader()  # build a btk reader object 
reader.SetFilename("./../data/Gait.c3d") # set a filename to the reader
acq = reader.GetOutput()                 # btk aquisition object
acq.Update()                             # Update ProcessObject associated with DataObject

Get information about the C3D file content:

In [5]:
print('Acquisition duration: %.2f s' %acq.GetDuration()) 
print('Point frequency: %.2f Hz' %acq.GetPointFrequency())
print('Number of frames: %d' %acq.GetPointFrameNumber())
print('Point unit: %s' %acq.GetPointUnit())
print('Analog frequency: %.2f Hz' %acq.GetAnalogFrequency())
print('Number of analog channels: %d' %acq.GetAnalogNumber()) 
print('Number of events: %d' %acq.GetEventNumber())
Acquisition duration: 4.87 s
Point frequency: 100.00 Hz
Number of frames: 487
Point unit: mm
Analog frequency: 1000.00 Hz
Number of analog channels: 28
Number of events: 0

Get the labels of the 3D points and analog channels:

In [6]:
print('Marker labels:')
for i in range(0, acq.GetPoints().GetItemNumber()):
    print(acq.GetPoint(i).GetLabel(), end='  ')   
print('\n\nAnalog channels:')
for i in range(0, acq.GetAnalogs().GetItemNumber()):
    print(acq.GetAnalog(i).GetLabel(), end='  ')   
Marker labels:
R.Shoulder  R.Offset  R.Elbow  R.Wrist  L.Shoulder  L.Elbow  L.Wrist  R.ASIS  L.ASIS  V.Sacral  R.Thigh  R.Knee  R.Shank  R.Ankle  R.Heel  R.Toe  L.Thigh  L.Knee  L.Shank  L.Ankle  L.Heel  L.Toe  R.Knee.Medial  R.Ankle.Medial  L.Knee.Medial  L.Ankle.Medial  R.Foot.Medial  R.Foot.Lateral  L.Foot.Medial  L.Foot.Lateral  C7  IJ  PX  

Analog channels:
f1x  f1y  f1z  m1x  m1y  m1z  f2x  f2y  f2z  m2x  m2y  m2z  L Rectus  R Rectus  L Vaste Ext  R Vaste Ext  L GRF  R GRF  L Ischio  R Ischio  L Biceps  R Biceps  L Tibialis  R Tibialis  L Sol  R Sol  L Jum Int  R Jum Int  

Get all events:

In [7]:
for i in range(0, acq.GetEvents().GetItemNumber()):
    print(acq.GetEvent(i).GetLabel() + ' at frame %d' %acq.GetEvent(i).GetFrame())   

This file doesn't have any event.

Get all metadata:

In [8]:
for i in range(acq.GetMetaData().GetChildNumber()):
    print(acq.GetMetaData().GetChild(i).GetLabel() + ':')
    for j in range(acq.GetMetaData().GetChild(i).GetChildNumber()):
        print(acq.GetMetaData().GetChild(i).GetChild(j).GetLabel(), end='  ')
    print('\n')
POINT:
USED  FRAMES  DESCRIPTIONS  LABELS  INITIAL_COMMAND  X_SCREEN  Y_SCREEN  LAB_ROT  RATE  UNITS  SCALE  DATA_START  

ANALOG:
USED  DESCRIPTIONS  LABELS  RATE  BITS  FORMAT  GEN_SCALE  OFFSET  UNITS  SCALE  

FORCE_PLATFORM:
USED  TYPE  ORIGIN  CORNERS  CAL_MATRIX  ZERO  CHANNEL  

Get the 3D position data of all markers as a Numpy array:

In [9]:
import numpy as np

data = np.empty((3, acq.GetPointFrameNumber(), 1))
for i in range(0, acq.GetPoints().GetItemNumber()):
    label = acq.GetPoint(i).GetLabel()
    data = np.dstack((data, acq.GetPoint(label).GetValues().T))
data = data.T
data = np.delete(data, 0, axis=0)  # first marker is noisy for this file
data[data==0] = np.NaN             # handle missing data (zeros)
In [10]:
data.shape
Out[10]:
(33L, 487L, 3L)
In [11]:
print(np.nanmin(data), np.nanmax(data))
-753.301107883 2438.03477883

There are 33 markers (with x, y, z coordinates) and 487 frames.
Let's visualize these data.

First, the vertical ground reaction force:

In [12]:
import matplotlib.pyplot as plt
%matplotlib inline

ana = acq.GetAnalog("f1z")
fz1 = ana.GetValues()/ana.GetScale()
ana = acq.GetAnalog("f2z")
fz2 = ana.GetValues()/ana.GetScale()
freq2 = acq.GetAnalogFrequency()
t = np.linspace(1, len(fz1), num=len(fz1))/freq2
plt.figure(figsize=(10, 4))
plt.plot(t, fz1, label='Fz1')
plt.plot(t, fz2, label='Fz2')
plt.legend()
plt.xlabel('Time [s]')
plt.ylabel('GRF vertical [N]')
plt.show()

Let's visualize the markers positions in a 3D animation:

In [13]:
# Use the IPython magic %matplotlib qt to plot a figure in a separate window
#  because the Ipython Notebook doesn't play inline matplotlib animations yet.
%matplotlib qt
In [14]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

dat = data[:, 130:340, :]
freq = acq.GetPointFrequency()

fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
ax.view_init(10, 150)
pts = []
for i in range(dat.shape[0]):
    pts += ax.plot([], [], [], 'o')

ax.set_xlim3d([np.nanmin(dat[:, :, 0]), np.nanmax(dat[:, :, 0])])
ax.set_ylim3d([np.nanmin(dat[:, :, 1])-400, np.nanmax(dat[:, :, 1])+400])
ax.set_zlim3d([np.nanmin(dat[:, :, 2]), np.nanmax(dat[:, :, 2])])
ax.set_xlabel('X [mm]')
ax.set_ylabel('Y [mm]')
ax.set_zlabel('Z [mm]')

# animation function
def animate(i):
    for pt, xi in zip(pts, dat):
        x, y, z = xi[:i].T
        pt.set_data(x[-1:], y[-1:])
        pt.set_3d_properties(z[-1:])   
    return pts

# Animation object
anim = animation.FuncAnimation(fig, func=animate, frames=dat.shape[1], interval=1000/freq, blit=True)

plt.show()
In [15]:
# get back the inline plot
%matplotlib inline

This is the animation you see in a separate window when you run the code above:

walking

To save the animation generated above as a video file (you need to have FFmpeg or ImageMagick installed):

In [18]:
anim.save('walking.mp4', writer='ffmpeg', fps=50)
# or
anim.save('walking.gif', writer='imagemagick', fps=50)

You can also write to a C3D file using BTK; look at its webpage.