Issues with extracting COP data

[ClayZhang999]

Hi,

First, thanks for developing this amazing tool!

I'm using ezc3d(version: 1.5.9) in Python to extract GRFs and COPs from the Qualisys mocap system. I have attached the example c3d file (walk_1.zip) and code below to replicate the issue.

#%%
from ezc3d import c3d
import pandas as pd
import numpy as np

def read_c3d(path):
    c = c3d(path, extract_forceplat_data=True)
    analog_fs = c['header']['analogs']['frame_rate']
    analog_headers = c['parameters']['ANALOG']['LABELS']['value']
    analog_data = c['data']['analogs']

    # Reshape the analog_data to 2D shape
    num_columns = analog_data.shape[1]
    num_rows = analog_data.shape[2]
    analog_data_reshaped = analog_data[0].reshape(num_rows, num_columns)

    # Create DataFrame for analog data and drop unwanted columns
    unwanted_columns = ['Channel_01', 'Channel_02', 'Channel_03', 'Channel_04',
                        'Channel_05', 'Channel_06', 'Channel_07', 'Channel_08',
                        'Channel_09', 'Channel_10', 'Channel_11', 'Channel_12',
                        'Channel_13', 'Channel_14', 'Channel_15', 'Channel_16', 'Sync']
    df_analog = pd.DataFrame(analog_data_reshaped, columns=analog_headers)
    df_analog.drop(columns=unwanted_columns, inplace=True)
    df_analog.insert(0, 'time(s)', np.arange(0, len(df_analog) / analog_fs, 1 / analog_fs))

    # Extract force plate data
    fp_1 = c['data']['platform'][0]['force']
    fp_1_reshaped = fp_1.T
    df_fp_1 = pd.DataFrame(fp_1_reshaped, columns=['fp1_X', 'fp1_Y', 'fp1_Z'])

    # Extract the center of pressure data
    cop_1 = c['data']['platform'][0]['center_of_pressure']
    cop_1_reshaped = cop_1.T
    df_cop_1 = pd.DataFrame(cop_1_reshaped, columns=['cop1_X', 'cop1_Y', 'cop1_Z'])

    # Merge the analog and force plate data
    df = pd.concat([df_analog, df_fp_1, df_cop_1], axis=1)
    return df, c

#%%
import matplotlib.pyplot as plt
df, c = read_c3d("...\walk_1.c3d") # Set file path here
plt.plot(df['cop1_X'], label='cop1_X')
plt.plot(df['cop1_Y'], label='cop1_Y')
plt.legend()
plt.show()

Then I got COPs on AP and ML directions from ezc3d like this (which I cannot make sense of them):

However, if I load this c3d file in modelling software like Visual3D and plot AP and ML COPs, they would look more sensible like below:

I hope this describes the question and please feel free to let me know if you need me to clarify.

Thanks again for your work and time.

Kind regards, Clay

[pariterre]

Hi @ClayZhang999

Thanks for reporting! I'll have a look today or tomorrow, but at first glance, it looks like a division by a near 0 number as the values are very high but for what it seems to be a single frame. I'll try to figure out why and I come back to you!

[ClayZhang999]

Hi,

Thank you very much for getting back to me and your effort.

Just for your reference, I dug out some similar problems reported by other users. Please see below.

https://github.com/pyomeca/ezc3d/issues/194#issue-840835545 https://github.com/pyomeca/ezc3d/issues/193#issue-840143965

In issue 194, you suggested removing data where forces are too small so that the trial can be trimmed to the stance phase only and then CoP can be calculated correctly. However, in my particular case, I cannot trim the data because I need to keep accelerometer signals in analog channels as a whole for synchronising with the external instrumented insole with IMUs.

So maybe this issue is worth a thorough debugging?

Kind regards, Clay

[pariterre]

Hi again!

As suspected, the main culprit is this line of code https://github.com/pyomeca/ezc3d/blob/2c2adf4ee78360851985572a5fdbe0a31892895b/src/modules/ForcePlatforms.cpp#L437 from where the center of pressure is computed.

When the force on the platform is too small, the center of pressure cannot be properly computed. I could have pre-removed these extreme values when they are detected, but I prefer giving the data as less process as possible to the user. The main reason is there is nothing I could do in background that the end user won't be able to do (i.e. detecting when the force is almost 0 and then setting the cop to 0 at these points). Worst, I could end up deleting actual data if the cut-off is not set properly...

For your specific case, I suggest to compute the norm of the force and removing all the cop data where the force are too small to be meaningful. Based on your data, a cutoff of 50N seems reasonable, but can obviously be adjusted. Slightly changing your read_c3d function, it gives:

# Extract the center of pressure data
force_cutoff = 50
cop_1 = c['data']['platform'][0]['center_of_pressure']
cop_1[:, np.linalg.norm(fp_1, axis=0) < force_cutoff] = 0
cop_1_reshaped = cop_1.T
df_cop_1 = pd.DataFrame(cop_1_reshaped, columns=['cop1_X', 'cop1_Y', 'cop1_Z'])

[ClayZhang999]

Hi,

Thank you very much for your suggestion. This workaround should be good for my specific case.

All the best, Clay