Determining gait cycle from heel marker motion capture data

I have motion capture data of participants walking, jogging and running on a treadmill. Participants had a marker on each heel. I need to determine each gait cycle from this data so that I can then calculate breast range of motion, velocity and acceleration for each gait cycle.

Answer 1

[recent works]

Obtain the position data (typically x, y, z coordinates) for each heel marker from the motion capture system. Assume you have the heel marker data in a matrix where each row corresponds to a frame and columns represent x, y, and z coordinates for the left and right heel markers.

% Load your heel marker data (assuming the data is in a variable called 'heelData')
% heelData should be an Nx6 matrix where N is the number of frames
% Columns are [Lx, Ly, Lz, Rx, Ry, Rz] for left and right heels

% Example: load data
% load('heelData.mat'); % Load your data file here

% Extract z-coordinate (vertical) data for left and right heels
Lz = heelData(:, 3); % Left heel vertical position
Rz = heelData(:, 6); % Right heel vertical position

% Define a threshold to detect heel strike (customize based on your data)
% This threshold might need to be adjusted depending on the data characteristics
threshold = 0.05; 

% Find heel strike (HS) events based on local minima and threshold crossing
[~, LHS] = findpeaks(-Lz, 'MinPeakHeight', threshold); % Inverted peaks
[~, RHS] = findpeaks(-Rz, 'MinPeakHeight', threshold);

% Find toe-off (TO) events if needed (local maxima or sharp increase)
% [~, LTO] = findpeaks(Lz, 'MinPeakHeight', threshold); % Example for toe-off detection

% Plotting to visualize
figure;
plot(Lz, 'b'); hold on;
plot(Rz, 'r');
plot(LHS, Lz(LHS), 'bo'); % Mark HS events for left heel
plot(RHS, Rz(RHS), 'ro'); % Mark HS events for right heel
xlabel('Frame');
ylabel('Vertical Position (z)');
title('Heel Marker Vertical Position and Gait Events');
legend('Left Heel', 'Right Heel', 'Left HS', 'Right HS');
hold off;

% Example output: displaying detected events
disp('Detected Heel Strikes for Left Heel:');
disp(LHS);
disp('Detected Heel Strikes for Right Heel:');
disp(RHS);

% Further analysis can include calculating parameters like velocity, acceleration, ROM, etc.
% Calculate velocity (simple finite difference approach)
L_velocity = diff(Lz) ./ diff((1:length(Lz))'); % Approximate velocity
R_velocity = diff(Rz) ./ diff((1:length(Rz))');

% Calculating acceleration (second derivative)
L_acceleration = diff(L_velocity) ./ diff((1:length(L_velocity))');
R_acceleration = diff(R_velocity) ./ diff((1:length(R_velocity))');

Answer 2

[William Rose]

@Katie,

Post some example data if you wish.

Here is a paper I co-authored on the topic of identifying gait cycle during teadmill running and overground running.

https://pubmed.ncbi.nlm.nih.gov/20478742/

The subjects ran on a treadmill that had a force plate, and they also ran over ground across a force plate. The vertical ground reaction force was used as the "gold standard" for determining foot strike and toe off. We compared that gold standard to five methods that did not use the force recording. The two methods that worked best for identifying foot strike accurately were the minimum height of the heel marker, and the time that heel marker vertical velocity changed from negative to positive. These methods were the best, and about equal in accuracy to one another, for both treadmill and overground. Toe off was best determined as time of peak knee extension, for treadmill and for overground.

Once you have the foot strike times identified accurately, you can scale each cycle to 100% and thn overage them, or do other analyses as dictated by your study gioals.

I have used the phrase foot strike instead of heel strike since some subjects may be mid- to forefoot strikers.

Good luck with your research.

[William Rose]

@Katie,

I have attached a sample mocap file from https://mocap.cs.sfu.ca/nusmocap/0005_Jogging001.txt downloaded 2024-07-30. I infer from the headers that columns 162-164 are X,Y,Z for left heel, and columns 198-200 are XYZ for right heel. The accompanying animation at the site above shows that the subject jogs in a circle with a diameter of a few meters during the capture, making about 1.25 revolutions. The .BVH file includes "Frame Time: 0.008333", i.e. 120 fps. There are 1377 frames, i.e. 11.5 seconds at 120 fps.

data=readmatrix('0005_Jogging001.txt');
HipL=data(:,142:144);
KneeL=data(:,146:148);
HeelL=data(:,162:164);
HipR=data(:,178:180);
KneeR=data(:,182:184);
HeelR=data(:,198:200);
disp([mean(HipL),mean(HipR)])

  410.1786  602.8888 -291.2815  462.2005  615.7610 -326.9488

disp([mean(KneeL),mean(KneeR)])

  421.4925  428.5886 -290.3981  431.1805  446.5158 -323.6278

disp([mean(HeelL),mean(HeelR)])

  479.3735  121.0343 -275.5647  496.6989  143.6351 -296.2394

I infer from the mean values that +X is forward, +Y is up, and +Z is to the right, and that the units are mm. Let's plot the left and right heel vertical coordinates verfsus time. The animation shows 15 right heel strikes and 14 or 15 left heel strikes.

N=length(HeelR); % number of frames
fps=120; % frames per second
t=(0:N-1)/fps; % time (s)
plot(t,HeelR(:,2),'-r',t,HeelL(:,2),'-g')
legend('Right','Left'); hold on
xlabel('Time (s)'); ylabel('Marker Height (mm)'); title('Heel Height')

The plot shows 15 minima for each, which is consistent with what we expect from the animation. Let's find the minima of the heel heights. Note that tsR and tsL are the heel strike times, in seconds.

[pksR,tsR]=findpeaks(-HeelR(:,2),fps,'MinPeakDistance',0.3);
pksR=-pksR;
[pksL,tsL]=findpeaks(-HeelL(:,2),fps,'MinPeakDistance',0.3);
pksL=-pksL;

Plot the minima

plot(tsR,pksR,'r*',tsL,pksL,'g*')
legend('Right','Left','R strike','L strike')

Now you have the heel strike times for this set of marker data. You can do likewise with your marker data.

I used findpeaks() which is in the Signal Processing Toolbox. I think findpeaks() is also in some other toolboxes.

[William Rose]

@Katie,

[Edit: Fix spelling errors. No changes to the code..]

You said "I need to determine each gait cycle from this data so that I can then calculate breast range of motion, velocity and acceleration for each gait cycle." I will illustrrate an approach, using the same marker data I posted earlier. This data set has X=forward,Y=up,Z=right for 53 markers at 120 frames per second for 11.5 seconds during jogging.

I showed above how to extract foot strike times from the heel markers. Left heel XYZ is columns 162-164. The left heel strike times (tsL) are found as follows (see previous post):

data=readmatrix('0005_Jogging001.txt');
HeelL=data(:,162:164);
N=length(HeelL);        % number of frames
fps=120;                % frames per second

Find the times of the minima of the heel vertical position

[pksL,tsL]=findpeaks(-HeelL(:,2),fps,'MinPeakDistance',0.3);
pksL=-pksL;

Plot left heel vertical position and the minima:

t=(0:N-1)/fps;          % time (s)
plot(t,HeelL(:,2),'-g',tsL,pksL,'g*')
legend('Left','L strike'); title('Left Heel Height')
xlabel('Time (s)'); ylabel('Marker Height (mm)')

Next: Use the left heel strike times (tsL) to analyze gait cycles. Remember this subject was jogging in a small circle, not running on a treadmill, so the whole body is moving around the lab in the horizontal plane (the X-Z plane), with a small amount of bouncing up-and-down (Y-axis) motion. We don't have breast markers in this data set. Let's use marker LUPA (probably left upper arm, but I'm not sure), columns 46-48, to illustrate vertical range of motion in each gait cycle.

UpaL=data(:,46:48);     % left upper arm XYZ (mm)

Compute M=number of complete cycles = number of left heel strikes - 1. Compute tc=vector with duration of each cycle (s), and tcI=vector with duration of each cycle (points).

M=length(tsL)-1;        % number of complete cycles
tc = diff(tsL);         % cycle durations (s)
tcI=round(tc*fps,0);    % cycle durations (Integer: points in each cycle)

Allocate an array of NaNs for left upper arm vertical position in each cycle. It will have one column for each cycle. The number of rows will equal the length of the longest cycle. The number of columns equals the number of complete cycles. I use NaNs because the cycles vary in length. Cycles with less-than-maximal length will have NaNs at the end of each column.

UpaLy=NaN(max(tcI),M);    % allocate array for Upa L vertical position 

Put the vertical (y) data for each cycle into the correct column of UpaLy:

for j=1:M
    UpaLy(1:tcI(j),j)=UpaL(t>=tsL(j) & t<tsL(j+1),2);
end

Make a set of M distinct colors for plotting:

colr=hsv2rgb([linspace(0,1,M)',ones(M,2)]);

Plot the vertical position versus cycle time, for each cycle:

figure
for j=1:M
    plot((0:tcI(j)-1)/fps,UpaLy(1:tcI(j),j),'Color',colr(j,:)); hold on
end
xlabel('Cycle Time (s)'); ylabel('UpaL vertical position (mm)'); grid on
title('Vertical Position for Each Stride'); legend(); hold off

If you want to compute the velocity or acceleration, differentiate the position signal once or twice. Be sure to use the appropriate value, 1/fps, for delta T, so that the units come out right. You may want to apply a gentle low pass filter each time you differentiate, since differentiation amplifies noise at high frequencies. There are many articles in the biomechanics literature on this subject.

If you want to get the mean signal for one stride cycle, then interpolate the signals from the different strides (use interp1()), so the signals from all strides have the same length. Then average the signals.

[Katie]

Hi William,

Thanks for this. I have got to the point of the script which produces the figure of the left heel height but I am now confused as what is happening next. Are you able to explain any further please?

[William Rose]

@Katie,

Once I have the 15 left heel strike times (in vector tsL), I use those times to divide the recording into M=14 compete stride cycles. Then I can look at any other marker (I chose left upper arm) and plot how that marker moves over each stride cycle. That is what is shown in the second plot: 14 traces of left upper arm marker height. The traces aer not all the same length because the stride cycle duraiton varies somewhat.

Details

I want to know the duration of each cycle, in seconds (for plotting) and in points (for extracting the data into separate cycles). Therefore I compute tc=vector of cycle durations (s), and tcI=vector of cycle durations (in points). For the duration in points, I multiply the time in seconds by the frames per seocnd, and round to the nearest integer. If I don't round, some of the numbers in tcI are slightly non-integer, which causes problems later.

M=length(tsL)-1;        % number of complete cycles
tc = diff(tsL);         % cycle durations (s)
tcI=round(tc*fps,0);    % cycle durations (Integer: points in each cycle)

Next , I want to make an array to contain the cycle-by-cycle data of interest. Arrays must be exactly rectangular, but the cycles have varying length. Therefore I allocate an array big enough to hold the longest cycle. Shorter cycles will have some unused array elements. I allocate an array, "UpaLy" (left upper arm, y-coordinate), of NaNs, with M columns, i.e. one column per cycle. The number of rows equals the length (in points) of the longest cycle. I use NaNs, because NaNs in the unused elements (for less than maximal-length cycles) cannot be confused with recorded data. Once the data for each cycle has been copied into array UpaLy, the columns for less-than-maximal length cycles will have NaNs at the ends.

The next thing is to copy the left upper arm data from the big array, data, into the array UpaLy. You may recall that data has 1377 rows by 213 columns. To do the copying, I use the heel strike times (vector tsL), and I use logical array indexing. See here for more on this useful concept. And see here. The basic idea is that cycle 1 is all the points from tsL(1) up to but not including tsl(2). Cycle 2 is all the points from tsl(2) up to but not including tsL(3). And so on.

for j=1:M
    UpaLy(1:tcI(j),j)=UpaL(t>=tsL(j) & t<tsL(j+1),2);
end

Then I plot the data from all M cycles on a single plot.

I recommend that you average the data of interest across the cycles, to get the average waveform for one cycle. Tis is standard in biomechanical analysis of data from walking or running. It is not comletely trivial, since the cyce lengths vary. That is why it is useful to interpolate to a common time base befofe averaging. If you want to do this and are unsure how to do it, please email.

If you would like to dicuss further, please email me securely by clicking on the envelope icon, which you will see at top right of the pop-up window that appears when you click on the "WR" circle next to my posts. Include your email address, so I can reply directly.

Answer 3

[Ayesha]

Is this dataset available for the generl public? I wish to use it for my research in marker based mocap systems.

[Katie]

No sorry it's not

[William Rose]

@Ayesha,

There is free-for-researchers mocap data at the Carnegie Mellon site:

http://mocap.cs.cmu.edu/