Problem when interpolate direction vectors with rotation interpolation method

[sysuyl]

Hello,

I am tring to interpolate a set of unit direction vectors' endpoints with spline. It means that the spline will interpolate/cross each direction vector's endpoint on a unit sphere surface. The result will look like this picture,

I used the following process,

• convert the direction vectors to unit quanterions/rotations
• do rotation interpolation with this library's method(e.g Squad or CatmullRom)
• and finnaly convert the interpolated samples back to vectors

The problem is that the interpolated result seems not correct with some inputs. Below is my code and test data,

import numpy as np

def interpolate_catmull_rom_centripetal(qs, num_samples, closed=False):
    """
    Interpolate quanterions/roations

      - https://github.com/AudioSceneDescriptionFormat/splines
      - https://splines.readthedocs.io/en/0.3.0/rotation/squad.html#Non-Uniform-Parameterization
    """
    from splines.quaternion import UnitQuaternion, Squad, CatmullRom

    endconditions = "closed" if closed else "natural"

    xyzw_list = [[q[1], q[2], q[3], q[0]] for q in qs]
    rotations = [UnitQuaternion.from_unit_xyzw(xyzw) for xyzw in xyzw_list]

    # note, alpha=0.5, centripetal catmul-rom spline
    spline = CatmullRom(rotations, alpha=0.5, endconditions=endconditions)
    # spline = Squad(rotations)

    ts = np.linspace(spline.grid[0], spline.grid[-1], num_samples, endpoint=False)
    sample_rotations = spline.evaluate(ts)

    sample_xyzw_list = [q.xyzw for q in sample_rotations]
    sample_qs = [[q[3], q[0], q[1], q[2]] for q in sample_xyzw_list]

    return sample_qs

def norm_vector(v):
    v = np.array(v)
    v /= np.linalg.norm(v)
    return v.tolist()

def quanterions_from_direction_vectors(vs):
    vs = [norm_vector(v) for v in vs]
    qs = [[0, *v] for v in vs]
    return qs

def quanterions_to_direction_vectors(qs):
    vs = [q[1:] for q in qs]
    vs = [norm_vector(v) for v in vs]
    return vs

def interpolate_vectors(vectors, sample_num, closed=True):
    """
    Interpolate vectors by converting vectors to quanterions/rotations and use quanterion interpolation.
    """
    # convert vector to quanterion
    qs = quanterions_from_direction_vectors(vectors)

    # interpolate quanterion/rotation
    sample_qs = interpolate_catmull_rom_centripetal(qs, sample_num, closed=closed)

    # convert quanterion back to vector
    sample_vs = quanterions_to_direction_vectors(sample_qs)

    return sample_vs

def interpolate_and_check_if_cross_control_point(vs):
    """
    Interpolate a spline and check if the sample points across control points.
    For each control point, compute the maximum dot value between it and all sample points. And the dot
    should be close to 1.0 if the samples cross it.
    """
    sample_vs = interpolate_vectors(vs, sample_num=1000, closed=True)

    dots = np.array(vs).dot(np.array(sample_vs).T)
    dot_max = np.max(dots, axis=1)

    print("dot max: ", dot_max)

def test():
    # data1, seems ok
    vs = [[-0.82537162, -0.5640741,  -0.02412931],
          [-0.61665761, -0.75438958, -0.22501047],
          [-0.5571503,  -0.8302387,   0.01694991],
          [-0.66785687, -0.71488637,  0.20713463]]
    vs = [norm_vector(v) for v in vs]
    print("data1: ", np.array(vs))
    interpolate_and_check_if_cross_control_point(vs)

    # data2, seems bad
    vs = [[-0.94211662, -0.32190701, -0.09376714],
          [-0.36887884, -0.92903411,  0.0287057 ],
          [ 0.52999741,  0.5104602,  -0.67715067],
          [-0.04777308,  0.8693673,  -0.49185181]]
    vs = [norm_vector(v) for v in vs]
    print("\ndata2: ", np.array(vs))
    interpolate_and_check_if_cross_control_point(vs)

if __name__ == "__main__":
    test()

The output is,

data1:  [[-0.82537157 -0.56407406 -0.02412931]
 [-0.61665762 -0.7543896  -0.22501047]
 [-0.55715028 -0.83023868  0.01694991]
 [-0.66785684 -0.71488634  0.20713462]]
dot max:  [1.         0.99999995 0.99999992 1.        ]

data2:  [[-0.94211656 -0.32190699 -0.09376713]
 [-0.36887884 -0.92903411  0.0287057 ]
 [ 0.52999744  0.51046023 -0.6771507 ]
 [-0.04777308  0.86936731 -0.49185182]]
dot max:  [ 1.          1.         -0.59510805 -0.18848043]

When I plot the original direction vectors and sampled vectors, it shows that the problem more clearly, test-data1, test-data2,

Besides, I'm not that confident about the direction vector and quanterion conversion, and maybe also use this library's interpolation in wrong way. But it worked for part of test data, and I also tested other implmentation of Squad method and interpolated correctly when this library fails.

[mgeier]

Your unit direction vectors have only 2 degrees of freedom, while 3D rotations and unit quaternions have 3. I don't think the splines.quaternion module is the right tool for your situation.

I'm sure that most of what I derived for unit quaternions can similarly be derived for a "normal" sphere. I don't know the details, though.

You can have a look at #5, which mentions an application with 2 degrees of freedom and it references this paper: https://doi.org/10.1016/S0010-4485(00)00049-X

I also tested other implmentation of Squad method and interpolated correctly when this library fails.

Out of curiosity, which implementation was that?

[sysuyl]

Hello,

First, I updated the code when checking if interpolated samples cross the control points, use interpolate_and_check_if_cross_control_point instead of the old one '''interpolate_and_compare_last''' which is incorrect. The new one will show that if the samples cross a control point, the dot product will near to be 1.0.

And I want to ask is there any problem in function interpolate_catmull_rom_centripetal when do quanterion interpolation? Because I'm not familiar with the concept of grid in code,

    ts = np.linspace(spline.grid[0], spline.grid[-1], num_samples, endpoint=False)
    sample_rotations = spline.evaluate(ts)

Out of curiosity, which implementation was that?

I tested the Squad(ShoemakeC1Spline) in interpolation-methods library. It seems just a directly implemenation of squad interpolation.

By the way, I'm interested in your library because it discussed the centripetal parameterization, which has the very nice property that it guarantees no cusps and no self-intersections.

Finally, thanks for your work!

[mgeier]

I'm not familiar with the concept of grid in code

spline.grid is the list of parameter values (a.k.a. time values) corresponding to the list of quaternions used to construct spline (there is one more grid value if closed=True was chosen).

Those time values could be given when constructing a spline, but when using the alpha parameter, they are instead chosen automatically (starting with zero).

spline.evaluate(spline.grid[0]) gives you the first quaternion in the spline.

spline.evaluate(spline.grid[-1]) gives you the last quaternion in the spline (if closed=False was given).