Read/write transforms in read_c3d and write_c3d

[felixchenier]

As discussed here: https://github.com/orgs/kineticstoolkit/discussions/116#discussioncomment-8644282

It's possible to store full transform matrices in c3d files. These transforms are supported by ezc3d. It would be logical to be able to read/write those transforms as requested by @Alek050.

[Alek050]

Hi @felixchenier,

Just started working on the change. The read_c3d is pretty simple since the rotations only hold rotation matrices and the labels (no issues with unit conversions etc.). However, I wanted to add some changes to the write_c3d but I think the function will change a bit, so before making any changes I wanted to discuss what is the best option.

The issue mainly lies in the assumption of the function that there must be point data in the c3d file. In files where there are rotation matrices, there is often no point data available. In the current format, that is not possible since the first check in the function is points._check_not_empty_data() which will raise an value error.

I was thinking about making points also default to None and raise an error if points, analogs, and rotations are all None:

if points is None and analgos is None and rotations is None:
    raise ValueError("Message")

if points is not None:
   # add points

if analogs is not None:
   # add analogs, may need some changes since it now checks wheter the frame rate related to the points frame rate

if rotations is not None:
   # add rotations

# wride c3d and add events

The main function change will that there is a possibility for a c3d file with analogs and without points. I am not sure if that makes sense, but I think this could be the case when you want to measure jump height using a force plate without any marker data?

I am asking since I am probably missing something and am not sure what your intentions and ideas have been thusfar.

[felixchenier]

Hi @Alek050

This is an interesting question. I agree that

def write_c3d(
    filename: str,
    points: TimeSeries,
    analogs: TimeSeries | None = None,
) -> None:

should be changed to:

def write_c3d(
    filename: str,
    points: TimeSeries | None = None,
    analogs: TimeSeries | None = None,
    rotations: TimeSeries | None = None,
) -> None:

with None by default for the points. It just makes sense.

Now for the implementation, I read on https://www.c-motion.com/v3dwiki/index.php?title=ROTATION_DATA_TYPE :

---

ROTATION:RATE and ROTATION:RATIO

RATE = Sampling Rate

RATIO = Ratio of ROTATION sample rate to POINT sample rate

The C3D file format is designed to store synchronized 3D Point data, analog data, and Rotation data; Rotation measurements are recorded at fixed intervals (set by the ROTATION:RATE parameter).

For consistency with C3D file prior to the introduction of the ROTATION group, the POINT:RATE is considered the base rate of a C3D file. If the ROTATION rate is stored at a higher, integer multiple ratio of the POINT:RATE, the ANALOG:RATIO parameter defines this multiple. Storing the RATE and RATIO is redundant, but for historical reasons both are stored.

---

Therefore, we must have a POINT:RATE, even if we have no points. My idea, in case there are no points, would be to create, at the very beginning of the function, a dummy points TimeSeries that would contain no data. It would however have a time attribute that we may copy from rotations if it exists, or from analogs if it doesn't. By setting this dummy "base rate" TimeSeries, we won't have edge cases anymore: the data rate for both analogs and rotations will be based on the points data rate, as expected.

As you mention, it would also solve the hypothetical problem of having only analogs in the C3D, which I effectively overlooked.

[Alek050]

Hi @felixchenier

Thanks for you clear response. I just updated the code to also include writing files. I have not tested it elaborately but it seems like it gives the expected output. There are two issues though. First, the c3d.write(filename) in write_c3d() raises an error:

miniconda3/envs/ktk_env/lib/python3.12/site-packages/ezc3d/ezc3d.py", line 1874, in __init__
    _ezc3d.c3d_swiginit(self, _ezc3d.new_c3d(*args))
                              ^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Could not open the c3d file: unspecified iostream_category error

I've played around a bit but the error pursues. Funny enough, the output c3d is exactly as expected and there is no problem reading in this file. Have you seen this error earlier?

Secondly, I could be wrong, but I think ezc3d might not be fully up to date regarding the transformations data. In the c3d you can add headers info about points and analogs (c3d["header"]["points"]["start_frame"]=0.0), but it does not recognise transformations yet, even though in the example c3d I pointed out in the discussion earlier there is a rotations key in the c3d["header"] dictionary. I do not think this is a big problem for now, but wanted to mention it for completeness.

I have opened a work in progress (WIP) pull request at #232, could you have a look at the code, provide some feedback and shine your light on the above mentioned issues? Much appreciated!

PS I have yet to write the tests for the write_c3d, but first want to make sure that the code is actually what you want it to be.

[felixchenier]

Hi @Alek050

Great work, really :-)

First, the c3d.write(filename) in write_c3d() raises an error

No, I don't remember to have encountered this error before. Does it happen with any file, or just when you try to write rotations?

In the c3d you can add headers info about points and analogs (c3d["header"]["points"]["start_frame"]=0.0), but it does not recognise transformations yet, even though in the example c3d I pointed out in the discussion earlier there is a rotations key in the c3d["header"] dictionary.

It may become a problem if data doesn't start at 0s. I understand that c3d["header"]["rotations"]["start_frame"] exist when we read the c3d, but we cannot create it when we write the c3d? If it's the case, then I think we should file it as a bug in ezc3d.

In general I agree with everything you wrote. This review also allowed me to spot a potential problem in writing events: at the moment, only the events from the points TimeSeries are written. The better way to define the events would be to merge all events from all TimeSeries:

points_events = points.copy(copy_data=False)
analogs_events = analogs.copy(copy_data=False) if analogs is not None else points.copy(copy_data=False)
rotations_events = rotations.copy(copy_data=False) if rotations is not None else points.copy(copy_data=False)

points_events.merge(analogs_events, in_place=True)
points_events.merge(rotations_events, in_place=True)
points_events.remove_duplicate_events(in_place=True)

Then use points_events instead of points to set the c3d events.

[Alek050]

Hi @felixchenier,

Just opened an issue in ezc3d, as I suspect it is a problem from their side. Thanks for providing the code for the events, I am not working with events within c3d files for now so this is new for me, but I will update the code based on your recommendations.

Regarding the error, I get it with writing any file, for instance when reading in the kinematics_racing_full.c3d provided by your package:

import kineticstoolkit as ktk
ts = ktk.read_c3d("data/kinematics_racing_full.c3d")
ktk.write_c3d("test_file", points=ts["Points"])

File "~/miniconda3/envs/ktk_env/lib/python3.12/site-packages/ezc3d/ezc3d.py", line 1874, in __init__
    _ezc3d.c3d_swiginit(self, _ezc3d.new_c3d(*args))
                              ^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Could not open the c3d file: unspecified iostream_category error

So I think it has something to do with my environment in combination with ezc3d, will look a bit further if I can find why its not working.

UPDATED:
Interstingly, I do not get an error when using the example provided by ezc3d:

import numpy as np

import ezc3d

# Load an empty c3d structure
c3d = ezc3d.c3d()

# Adjust some mandatory values of the parameters and fill the data with random values
c3d["parameters"]["POINT"]["RATE"]["value"] = [100]
c3d["parameters"]["POINT"]["LABELS"]["value"] = ("point1", "point2", "point3", "point4", "point5")
c3d["data"]["points"] = np.random.rand(3, 5, 100)

c3d["parameters"]["ANALOG"]["RATE"]["value"] = [1000]
c3d["parameters"]["ANALOG"]["LABELS"]["value"] = ("analog1", "analog2", "analog3", "analog4", "analog5", "analog6")
c3d["data"]["analogs"] = np.random.rand(1, 6, 1000)

# Create a custom parameter to the POINT group
c3d.add_parameter("POINT", "NewParam", [1, 2, 3])

# Create a custom parameter a new group
c3d.add_parameter("NewGroup", "NewParam", ["MyParam1", "MyParam2"])

# Write a new modified C3D and read back the data
c3d.write("temporary.c3d")

Similarly, the example of ktk also seems to work (even when reading in only the analogs data):

import kineticstoolkit.lab as ktk
import numpy as np

points = ktk.TimeSeries()
points.time = np.linspace(0, 10, 10*240, endpoint=False)
points.data["Marker1"] = np.ones((2400, 4))
points.data["Marker2"] = np.ones((2400, 4))

analogs = ktk.TimeSeries()
analogs.time = np.linspace(0, 10, 10*2400, endpoint=False)
analogs.data["Signal1"] = np.sin(analogs.time)
analogs.data["Signal2"] = np.cos(analogs.time)

rotations = ktk.TimeSeries()
rotations.time = np.linspace(0, 10, 10*480, endpoint=False)
rotations.data["Rotation1"] = np.ones((4800, 4, 4))
rotations.data["Rotation2"] = np.ones((4800, 4, 4))

ktk.write_c3d("testfile.c3d", rotations=rotations)
ktk.write_c3d("testfile.c3d", points=points)
ktk.write_c3d("testfile.c3d", analogs=analogs)
ktk.write_c3d("testfile.c3d", points=points, analogs=analogs, rotations=rotations)

I am not really sure why it does not work since type(ts["Points"]) and type(points)) both give are <class 'kineticstoolkit.timeseries.TimeSeries'>. (ts as difined in the first code block and points as defined in the last code block)

[felixchenier]

Hi @Alek050

It's weird, but then it reminded my something, an issue that I reported myself several months ago on ezc3d:

pyomeca/ezc3d#252

Is it possible that:

1. You're on Windows, and
2. In one case, you were trying to write in a folder that contained accented characters, and another time you were not?

I did a workaround for this issue on Windows in read_c3d, but not in write_c3d, and most probably I should do this also for write_c3d. The trick would be to write to the temp folder and then move the file back to its final place. I created a new issue for it: #233. Could you please fill the information about your OS and OS version please?

Thanks

[Alek050]

Hi @felixchenier,

Thanks for looking into this, I replied in the new opened issue. I was looking into your earlier mentioned suggestion regarding the events and found some potential issues with this approach:

First, merging TimeSeries only seem to work when the number of frames/ time vecotors are equal:

points = ktk.TimeSeries(time=np.linspace(0, 10, 100,), data={"point1": np.random.rand(100)})
analogs = ktk.TimeSeries(time=np.linspace(0, 10, 1000), data={"force1": np.random.rand(1000)})

points.merge(analogs)

Traceback (most recent call last):
  File "/~/test.py", line 9, in <module>
    points.merge(analogs)
  File "/~/kineticstoolkit/kineticstoolkit/timeseries.py", line 3503, in merge
    raise ValueError(
ValueError: Time vectors do not match, resampling is required.

Obviously, this makes sense. We could of course simply resample the analogs, merge the TimeSeries, combine and add the events. I think this should work since I guess the events are primarily coupled to the points. However, I am not sure of this since I've not worked with events yet.

Secondly, I found in the resample attribute of TimeSeries in the docstring a part that said:

"""
Caution
        -------
        Attempting to resample a series of homogeneous matrices would likely
        produce non-homogeneous matrices, and as a result, transforms would not
        be rigid anymore. This function can't detect if you attempt to resample
        series of homogeneous matrices, and therefore won't generate an
        error or warning.
"""

When resampling rotations to points, this applies. However, as long as we do not set in_place=True I think the data change shouldn't matter since we already added the rotations data to the c3d writer. But I just wanted to check if you agree with this line of thought.

[felixchenier]

@Alek050 You know you rock, right? ;-)

Yeah, I was thinking the same, this merging stuff is running from problems. The way better way is to do something like:

points = points.copy()

if analogs is not None:
    points.events.extend(analogs.events)
if rotations is not None:
    points.events.extend(rotations.events)
points.remove_duplicate_events(in_place=True)

# Then add the events as in the original function

Much better isn't it?

[Alek050]

Hi @felixchenier, Thanks for the appreciation!

To be fair, it is truly insane that you have methods like merge, resample and extend implemented and lying around for anyone to use, this makes the package really convenient. Happy I can contribute!

I wanted to finalize the tests, but stumbled on another unexpected behaviour in the writing of ezc3d when using rotations, so I extended my issue there, lets hope it can be fixed easily so that we can finish this issue.

[Alek050]

Hi @felixchenier

Just wanted to give a quick update, I had some conversations with pariterre regarding how to change the headers section in ezc3d. He mentioned that it is not the prefered way, since it will be overwritten based on the parameters. However, this has been the way the kineticstoolkit has done it so far (c3d["header"]["points"]["first_frame"] = round(points.time[0] * point_rate)). This solution is not possible for the rotations since they are not initialised in the the headers and I can not add them manually.

For the kineticstoolkit a solution would be to ignore the start frame when writing files and there are rotations, but this can only be done if different systems (points, analogs) always start at the same time, which is not an assumption we can make I guess. On top of that it could create problems with backward compatability since we need to change the times of the events, points, and analogs will change based on what is expected.

I think we will have to wait if things change on the ezc3d side, but let me know if you have any other ideas.

[felixchenier]

Hi @Alek050
Thanks for the work and for the update. I've been following your issue on ezc3d. It's difficult to answer definitely until we get a clear expected behavior in ezc3d. At least we now have:

• full read_c3d support
• partial write_c3d support

And I guess that documented and robust partial support is preferred to no support at all. For now, maybe we should make a check condition where as soon as there are rotations involved, time[0] must be 0s. for everything? I think it's the safest direction that we can take, since it shouldn't break, and it will allow users to write rotations (if there are, really! 😉) even if they need to offset their time.

This would let us finish this issue, and then we could open a new, related improvement issue to allow writing rotations that start at other than 0s?

[felixchenier]

Done in branch PR232

[felixchenier]

Hi @Alek050

A new version of ezc3d has been released today with the fix for the rotations. However, this test that you wrote in large part fails and I can't find why because it should pass. In this test, rotations have twice the sample rate of the points. When we save and load back, they are 4 times as fast as they were recorded. This problem doesn't happen when points and rotations have the same sample rate, everything works well.

    rotations = ktk.TimeSeries()
    rotations.time = np.linspace(0, 1, 240, endpoint=False)
    rotations.data["pelvis_4X4"] = ktk.geometry.create_transforms(
        "x", np.linspace(0, 2*np.pi, 240, endpoint=False)
    )
    rotations.add_event(0.5, "TestEvent")

    # add some point data
    points = ktk.TimeSeries()
    points.time = np.linspace(0, 1, 120, endpoint=False)
    points.data["point1"] = np.random.rand(120, 4)
    points.data["point1"][:, 3] = 1

    ktk.write_c3d("test.c3d", rotations=rotations, points=points)
    c3d = ktk.read_c3d("test.c3d")
    assert np.allclose(
        c3d["Rotations"].data["pelvis_4X4"], rotations.data["pelvis_4X4"]
    )
    assert c3d["Points"].events[0].name == "TestEvent"
    assert np.allclose(c3d["Points"].data["point1"], points.data["point1"])
    os.remove("test.c3d")

I wonder if it's something in my code or in ezc3d, and since you work with a custom version of ezc3d, I'd like to know if this code works or not on your side.

By the way I'm on the PR232 branch of the main kineticstoolkit repository.

Tx