pyomeca / ezc3d

Easy to use C3D reader/writer for C++, Python and Matlab
https://pyomeca.github.io/Documentation/ezc3d/index.html
MIT License
142 stars 44 forks source link

Possible bug in relation to the new rotation matrices feature? #327

Closed felixchenier closed 2 months ago

felixchenier commented 2 months ago

Hi @pariterre

I have this bug in KTK https://github.com/kineticstoolkit/kineticstoolkit/issues/229#issuecomment-2118586183 with C3D files with rotations. It's probably related to #316

All using ezc3d, if I create a C3D file with rotations saved at twice the point rate, and then load it back, I get rotations at 4x the point rate. I most certainly see a dumb mix between * and / somewhere, and I don't think it's in my code. Would you mind checking if it could be on ezc3d side?

Tx!

pariterre commented 2 months ago

Hi @felixchenier I am currently in vacation, I can have a look when I am back home (about 2 weeks)! Hopefully, this is okay for you!

felixchenier commented 2 months ago

@pariterre Absolutely no problem! It's for a new feature in KTK that's unreleased yet. Take your time and enjoy your vacation! 🍺

pariterre commented 2 months ago

Hi again! Was in a café and had some time haha

Using an ezc3d implementation of what you provided in the As stated in https://github.com/kineticstoolkit/kineticstoolkit/issues/229, I noticed a bug in the data dispatching. It made sense since Analog were not dispatch the same. Now it seems coherent.

That said, I could not reproduce the 2x error based on this implementation.

For reference the following code now pass

import ezc3d
import numpy as np
np.random.seed(42)
n_frames = 2
ratio = 4

class Dummy:
    pass

points = Dummy()
points.names = ()
points.data = np.zeros((4, 0, n_frames))
points.data[3, :] = 1

rotations = Dummy()
rotations.names = ('pelvis_4X4', 'pelvis2_4X4', 'pelvis3_4X4', 'pelvis4_4X4', 'pelvis5_4X4', 'pelvis6_4X4')
rotations.data = np.random.rand(4, 4, 6, n_frames * ratio)
rotations.data[:, :, :, 0] = 1
rotations.data[:, :, :, 1] = 2
rotations.data[:, :, :, 2] = 3
rotations.data[:, :, :, 3] = 4
rotations.data[:, :, :, 4] = 5
rotations.data[:, :, :, 5] = 6
rotations.data[:, :, :, 6] = 7
rotations.data[:, :, :, 7] = 8

c3d = ezc3d.c3d()

c3d['parameters']['POINT']['UNITS']['value'] = ['m']
c3d['parameters']['POINT']['RATE']['value'] = [n_frames]
c3d['parameters']['POINT']['LABELS']['value'] = points.names
c3d['data']['points'] = points.data

c3d.add_parameter("ROTATION", "RATE", [n_frames * ratio])
c3d.add_parameter("ROTATION", "LABELS", rotations.names)
c3d['data']['rotations'] = rotations.data

c3d.write("test.c3d")

c3d_2 = ezc3d.c3d("test.c3d")

assert np.allclose(c3d_2.data["rotations"], rotations.data)
assert np.allclose(c3d_2.data["points"].data, points.data)

Could you provide an example that I could run using vanilla ezc3d that reproces your bug?

felixchenier commented 2 months ago

Thanks @pariterre

The test you wrote in https://github.com/kineticstoolkit/kineticstoolkit/issues/229 passes. And I added this to this unit test:

assert np.allclose(
    c3d_2["parameters"]["ROTATION"]["RATE"]["value"],
    c3d["parameters"]["ROTATION"]["RATE"]["value"],
)

which is important because this is what I thought was the issue. But it's not, it passes.

The problem must be in my code then, back to mystery...

Thanks for checking!

felixchenier commented 2 months ago

@pariterre OK after inspecting my code I couldn't find anything, but now that I think more, the subframe inversion was sure problematic since all I do with the data is to permute it, and then I create a time array based on the sampling frequency, which is completely unrelated to the data. This would sure explain why I'd have twice the sampling frequency (not 4X though) but I'll try again with the future release of ezc3d (no stress for that!).

pariterre commented 2 months ago

@felixchenier Can you confirm that I should close this issue? :)

felixchenier commented 2 months ago

@pariterre

Félix: I guess you can :-) I didn't try it because I'll wait for a release (much simpler for me), but I'm confident you found the problem, and I'll reopen if it doesn't!

Benjamin: NO, I don't wanna see this issue raised again!

Félix: But it will work.

Benjamin: It will.