Potential bug when writing c3d file with rotations

[Alek050]

Hi,

I am using python 3.12 and ezc3d 1.5.9 on a macOS system. I was working on writing c3d files using this great package and was thrilled to see that you added support for c3d files with rotation matrices last year. Reading in the files is really easy and works as expected.

However, when writing c3d files with rotations I noticed some unexpected behaviour (which might be because I am not fully aware of the proper usecase of the package). To exemplify I used this example file from the c-motion site. When reading in the data, everything works as expected:

from ezc3d import c3d
read_file = c3d("C3DRotationExample.c3d")
read_file["data"]["rotations"] # np.ndarray of shape (4, 4, 21, 340)
read_file["header"]["rotations"] # {'size': 21, 'frame_rate': 85.0, 'first_frame': 0.0, 'last_frame': 28899.0}

When writing files, it seems like the "rotations" are not expected:

from ezc3d import c3d
write_file = c3d()
write_file["data"].keys() # dict_keys(['points', 'meta_points', 'analogs'])
write_file["header"].keys() # dict_keys(['points', 'analogs', 'events'])

# manually try to add header info and data
import numpy as np
trans_mat = np.random.rand(4, 4, 100)
write_file["data"]["rotations"] =  trans_mat[:, :, np.newaxis, :] # works fine, just adding a new key
write_file["header"]["rotations"] = {'size': 1, 'frame_rate': 85.0, 'first_frame': 0.0, 'last_frame':100.0} # raises a TypeError `*** TypeError: 'Header' object does not support item assignment`

As you can see, the "rotations" is not pre-added to the data and header objects, whereas all other possibilities are. For the data that is not a big problem since you can just add a new key, however, that is not allowed for the headers. Is there a workaround for this?

ADDED LATER:
I think I may have stumbled on another unexpected behaviour when writing rotations. When reading rotation matrixes where there are missing values, there are no real problems:

import ezc3d
file = ezc3d.c3d("C3DRotationExample.c3d")
print(file["data"]["rotations"].shape)
print(file["data"]["rotations"][:, :, 1, 0])

4, 4, 21, 340)
[[nan nan nan nan]
 [nan nan nan nan]
 [nan nan nan nan]
 [ 0.  0.  0.  1.]]

If I now want to write that same rotation matrix, and read it again, I get the following:

time = np.linspace(0, 1, 10, endpoint=False)
frame_rate = 10

# create empty points
c3d = ezc3d.c3d()
c3d["header"]["points"]["size"] = 0
c3d["header"]["points"]["frame_rate"] = frame_rate
c3d["header"]["points"]["first_frame"] = 0
c3d["header"]["points"]["last_frame"] = 9
c3d.add_parameter("POINT", "RATE", [frame_rate])
c3d.add_parameter("POINT", "LABELS", [])
c3d.add_parameter("POINT", "UNITS", ["m"])
c3d.add_parameter("POINT", "USED", [0])
c3d["data"]["points"] = np.empty((4, 0, len(time)))

# create rotations
c3d["data"]["rotations"] = file["data"]["rotations"][:, :, 1:2, 0:10]

c3d.add_parameter("ROTATION", "RATE", [frame_rate])
c3d.add_parameter("ROTATION", "LABELS", ["pelvis_4X4"])
c3d.add_parameter("ROTATION", "UNITS", [])
c3d.add_parameter("ROTATION", "USED", [1])

# before writing and reading
print(c3d["data"]["rotations"][:, :, 0, 0])

c3d.write("test.c3d")
file2 = ezc3d.c3d("test.c3d")

# after writing and reading
print(file2["data"]["rotations"].shape)
print(file2["data"]["rotations"][:, :, 0, 0])

[[nan nan nan nan]
 [nan nan nan nan]
 [nan nan nan nan]
 [ 0.  0.  0.  1.]]
(4, 4, 1, 10)
[[nan nan nan nan]
 [nan nan nan nan]
 [nan nan nan nan]
 [nan nan nan nan]]

As you can see, after writing and reading the full matrix is filled with nans. Although it is not a huge problem in analysis, since you can not read the transformations from the initial transoframtion matrices either. I do think it is unexpected behaviour. In my situation, it creates problems when writing tests for the kineticstoolkit where I want to test if we correctly parse the writing and reading of ezc3d commands.

[pariterre]

Thanks for reporting! I am not surprised that there are still bugs in rotations support as not much people use them! I'll have a look soon :)

[pariterre]

Dear @Alek050

I'll slowly answer the comments in different posts as I am going through the problems, sorry if it spams your email box!

As you can see, the "rotations" is not pre-added to the data and header objects, whereas all other possibilities are. For the data that is not a big problem since you can just add a new key, however, that is not allowed for the headers. Is there a workaround for this?

This is actually expected behavior. I made the header is read-only because it is enterely determined by the parameters section (it is redundant information). Hence, the user is not expected to fill the header.rotation sections, it will be filled automatically at write time, if the the ROTATION parameter if properly provided.

Another way to say it is that any information the user would provide in the header section (if they changed it) is simply ignored when writting a new C3D based on the parameters section, this is true for rotations, but also points and analogs

[pariterre]

As you can see, after writing and reading the full matrix is filled with nans. Although it is not a huge problem in analysis, since you can not read the transformations from the initial transoframtion matrices either.

Do you mean: if there are data it works properly, but if there are missing data, instead of getting the last row being 0 0 0 1, it now gets to nan nan nan nan. If that is, the reason if the company who created the original data actually stored 0 0 0 1 in the file while ezc3d would fill the matrices with nan and only fill the last line when there is data to fill.

If my understanding is correct, it is expected behavior. The reason I chose to do that is that filling with nan is a first initialisation and filling the last line with 0 0 0 1 is a second initialisation. I figured it was loosing time to do that while there isn't any information anyway.

One way to fix that for the tests could be to ignore the last line as it is always a line equals to 0 0 0 1 when the rotation/translation parts of the matrix are relevant.

[Alek050]

Hi @pariterre

Thanks for both explenations! I was not aware that the headers object was written automatically based on the parameters, but then it makes sense.

I understand your second point as well, and from a efficiency and functionality perspective you are totally right. From a users perspective I still think it is a bit unexpected to write an array with partially floats, and a get back fully nans, but I am also not aware of how much extra time/resources it would take to still do that. On top of that, you created a great package and your suggestion solves my current problem.

Thanks for the clear explenation and I think the issue can be closed!

[Alek050]

Sorry to bother you again @pariterre, but I still have some questions regarding writing c3d files with rotations:

QUESTION 1:
Firstly, there might be something with the generation of the last_frame and first_frame in the header of the rotations:

>>> import ezc3d
>>> c3d_rotations= ezc3d.c3d("data/C3DRotationExample.c3d")
>>> print(c3d_rotations["header"]["rotations"])
{'size': 21, 'frame_rate': 85.0, 'first_frame': 0.0, 'last_frame': 28899.0}
>>> print(c3d_rotations["data"]["rotations"].shape)
(4, 4, 21, 340)

Here, the last_frame in the header is 28899 (the product of the frame_rate and the actual number of frames (340)). In similar data where I have point data this is different:

>>> c3d_points = ezc3d.c3d("data/my_points_data.c3d")
>>> c3d_points["header"]["points"]
{'size': 142, 'frame_rate': 200.0, 'first_frame': 0, 'last_frame': 1232}
>>> c3d_points["data"]["points"].shape
(4, 142, 1233)

Here the last_frame corresponds with the number of data_points. This is similar for other c3d files I have with points and rotations. I would at least expect both the refer to the same value, am I overlooking something?

QUESTION2
How do I specify the the start_frame if I can not change the header?

>>> import kineticstoolkit as ktk
>>> import ezc3d
>>> filename = ktk.doc.download("walk.c3d")
>>> c3d = ezc3d.c3d(filename)
>>> c3d["header"]["points"]
{'size': 96, 'frame_rate': 100.0, 'first_frame': 152, 'last_frame': 372}
>>> c3d["header"]["analogs"]
{'size': 248, 'frame_rate': 2000.0, 'first_frame': 3040, 'last_frame': 7459}
>>> c3d["parameters"]["EVENT"]["TIMES"]["value"][1]
array([2.22350001, 2.89199996, 1.66050005, 2.36400008, 3.32999992,
       1.80999994, 2.74000001, 3.5       ])

As you can see, both the points and the analogs start after 1.52 seconds, which is important since the events are linked to the data based on the time. If I were to write this to a new c3d file:

import kineticstoolkit as ktk
import ezc3d
filename = ktk.doc.download("walk.c3d")
c3d = ezc3d.c3d(filename)

c3d_writer = ezc3d.c3d()
c3d_writer.add_parameter("POINT", "RATE", [100.])
c3d_writer.add_parameter("POINT", "LABELS", c3d["parameters"]["POINT"]["LABELS"]["value"])
c3d_writer.add_parameter("POINT", "UNITS", ["mm"])
c3d_writer.add_parameter("POINT", "USED", [len(c3d["parameters"]["POINT"]["LABELS"]["value"])])
c3d_writer.add_parameter("POINT", "DATA_START", 1)
c3d_writer.add_parameter("TRIAL", "ACTUAL_START_FIELD", [c3d["header"]["points"]["first_frame"] + 1, 0])
c3d_writer.add_parameter("TRIAL", "ACTUAL_END_FIELD", [c3d["header"]["points"]["last_frame"] + 1, 0])
c3d_writer["data"]["points"] = c3d["data"]["points"]

c3d_writer.add_parameter("ANALOG", "RATE", [2000.])
c3d_writer.add_parameter("ANALOG", "LABELS", c3d["parameters"]["ANALOG"]["LABELS"]["value"])
c3d_writer.add_parameter("ANALOG", "USED", (len(c3d["parameters"]["ANALOG"]["LABELS"]["value"])))
c3d_writer["data"]["analogs"] = c3d["data"]["analogs"]

for event_time, event_label, event_context in zip(c3d["parameters"]["EVENT"]["TIMES"]["value"].T, c3d["parameters"]["EVENT"]["LABELS"]["value"], c3d["parameters"]["EVENT"]["CONTEXTS"]["value"]):
    c3d_writer.add_event(event_time, event_context, event_label)

c3d_writer.write("test.c3d")

new_c3d = ezc3d.c3d("test.c3d")
print(new_c3d["header"]["points"])
print(new_c3d["header"]["analogs"])

{'size': 96, 'frame_rate': 100.0, 'first_frame': 0, 'last_frame': 220}
{'size': 248, 'frame_rate': 2000.0, 'first_frame': 0, 'last_frame': 4419}

The header is fully filled, except for the first_frame variable. I played around with the TRIAL group and the ACTUAL_START_FIELD and ACTUAL_END_FIELD, as they represent the start and end frames of the actual data. For intance, c3d["parameter"]["TRIAL"]["ACTUAL_START_FIELD"] returns [153, 0]. But I have not been able to get the header object to be changed accordingly.

PS, its fine if you answer in different comments! I appreciate your expertise.

[pariterre]

Hi again!

Firstly, there might be something with the generation of the last_frame and first_frame in the header of the rotations:

Can I see the parameter "ROTATION:RATIO" you set for the rotations? It may seem counter intuitive, but the number of frames for the rotations (as well as the analogs) is based on the number of points (even though there are no points in the C3D). The computation is:

self.header.frameRate() * rotation_info.ratio() * (self.header.lastFrame() + 1) - 1

where frameRate is the point frame rate, ratio should probably be 1 in your case and last frame is the last frame of the points. The reason for this confusing way of computing the last frame is that historically C3D was designed for points with other data).

[pariterre]

How do I specify the the start_frame if I can not change the header?

As of now, it seems that it is indeed not possible to do that from the Python interface! I'll have to add a setter... sorry about that!

EDIT: If you want a dirty workaround for the short term, you can tap into the core of the c3d by accessing the c3d.c3d_swig structure where you can use setFirstFrame(value). (not tested)

[pariterre]

RE-EDIT:
As I am skimming the code, I think I know why I did not previously allowed to change the first_frame value. The "number of frames" is computed by using the last_frame and the first_frame (subtracting them), which make sense. But, if last_frame and first_frame can be changed freely, then ezc3d gets confused on where to start and stop reading the data for the ANALOGs and the ROTATIONs. I may have to dig a bit deaper into the code to update the code properly if first_frame is not 0 (and last_frame is not the last_frame).

[Alek050]

Can I see the parameter "ROTATION:RATIO" you set for the rotations? It may seem counter intuitive, but the number of frames for the rotations (as well as the analogs) is based on the number of points (even though there are no points in the C3D). The computation is:

self.header.frameRate() * rotation_info.ratio() * (self.header.lastFrame() + 1) - 1

where frameRate is the point frame rate, ratio should probably be 1 in your case and last frame is the last frame of the points. The reason for this confusing way of computing the last frame is that historically C3D was designed for points with other data).

Certainly:

>>> c3d_rotations["header"]["points"]
{'size': 0, 'frame_rate': 85.0, 'first_frame': 0, 'last_frame': 339}
>>> c3d_rotations["parameters"]["ROTATION"]["RATIO"]
{'type': 2, 'description': '', 'is_locked': False, 'value': array([1])}

Assuming that the self.header refers to the points header, the calculation would be 85 * 1 * (339 + 1) - 1 which is indeed 28899. However, last_frame in the rotations then does not refer to the index of the last frame availabe (there are only 340 frames) but to number of frames times * the frame rate. This is in contrast to the value of last_frame in the points where it refers to the actual index of the last frame, which would make more sense to me. Is this just something which is counterintuitive to me but actual logical as you shortly explained here?

The reason for this confusing way of computing the last frame is that historically C3D was designed for points with other data).

[Alek050]

RE-EDIT:
As I am skimming the code, I think I know why I did not previously allowed to change the first_frame value. The "number of frames" is computed by using the last_frame and the first_frame (subtracting them), which make sense. But, if last_frame and first_frame can be changed freely, then ezc3d gets confused on where to start and stop reading the data for the ANALOGs and the ROTATIONs. I may have to dig a bit deaper into the code to update the code properly if first_frame is not 0 (and last_frame is not the last_frame).

I understand your line if thought, and the previous solution for us was to just change the header manually, which worked fine:

# previous code ...
c3d_writer["header"]["points"]["first_frame"] = 152
c3d_writer["header"]["analogs"]["first_frame"] = 3040
print(c3d_writer["header"]["points"])
print(c3d_writer["header"]["analogs"])

c3d_writer.write("test.c3d")

new_c3d = ezc3d.c3d("test.c3d")
print(new_c3d["header"]["points"])
print(new_c3d["header"]["analogs"])

{'size': 0, 'frame_rate': 0.0, 'first_frame': 152, 'last_frame': 0}
{'size': 0, 'frame_rate': 0.0, 'first_frame': 3040, 'last_frame': -1}
{'size': 96, 'frame_rate': 100.0, 'first_frame': 152, 'last_frame': 372}
{'size': 248, 'frame_rate': 2000.0, 'first_frame': 3040, 'last_frame': 7459}

As you can see, the last_frame is automatically updated an now the first frame is kept in place. This may not be the right way to do it officially, but it does work for the points and analogs. However, as mentioned before (my first message): the rotations key is not assigned in headers and I am not allowed to create one manually, nor can I change the start_frame using the parameters. For me a solution would be if the headers also initialised an empty rotatations, although I could understand that you might see more problems emerging from this than that it solves in the backend.

[pariterre]

Is this just something which is counterintuitive to me but actual logical as you shortly explained here?

Nah, that is just me being dumb. Obviously the computation should not be self.header.frameRate() * rotation_info.ratio() * (self.header.lastFrame() + 1) - 1 but simply rotation_info.ratio() * (self.header.lastFrame() + 1) - 1.

Fixed in #317 !

[pariterre]

the rotations key is not assigned in headers and I am not allowed to create one manually

I think it makes sense to have a standard output, regardless of the internal data (that is returning a structure with "rotations" set to 0 instead of skipping it). I've added it in the same PR.

I also added easy accessors to the dictionary keys, so now, for instance, header can be accessed using:

import ezc3d
c = ezc3d.c3d("path")

print(c.header)  # Equivalent to c["header"]
print(c.parameters.POINT.USED)  # Equivalent to c["parameters"]["POINT"]["USED"]
print(c.data.points)  # Equivalent to c["data"]["points"]

[Alek050]

That's a great update @pariterre ! Thanks for the quick and clear replies and solves!

[pariterre]

Please let me know if it fixes your issues!

[Alek050]

I think all issues are fixed:

• The last_frame in the header rotations object now represents the index of the last frame
• The header object now has a initialized "rotations" key
• If I change the start_frame in header["rotations"] it works similar as with the points and analogs

Thanks again!