New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP / ENH: spatial.transform.Spline [GSOC 2018] #9024

Open
wants to merge 24 commits into
base: master
from

Conversation

Projects
None yet
5 participants
@adbugger
Copy link
Contributor

adbugger commented Jul 12, 2018

@nmayorov , preprocessing during initialization of class has been implemented.

@adbugger adbugger changed the title WIP / ENH: Quaternion cubic spline interpolation [GSOC 2018] WIP / ENH: spatial.transform.Spline [GSOC 2018] Jul 12, 2018

w_int, diff = _intermediate_step(np.zeros(len(rotations) - 2, 3),
rotvecs, dt_inv)

while diff >= 1e-6:

This comment has been minimized.

@rgommers

rgommers Jul 12, 2018

Member

Should this 1e-6 magic number be configurable? and why is 1e-6 the right number?

This comment has been minimized.

@adbugger

adbugger Jul 13, 2018

Author Contributor

This 1e-6 was mentioned in the paper which Nikolay provided. I'm not sure it suits all applications. Should I add a tolerance argument to the function and explain the convergence steps in the docstring?

This comment has been minimized.

@rgommers

rgommers Jul 13, 2018

Member

Not sure without having read the paper. At least I would suggest documenting in a comment here where this number came from, and why it's a "good enough" number (I hope the paper says something about that?).

@adbugger adbugger force-pushed the adbugger:rotation_spline branch from 77f2a21 to 016ddd3 Jul 20, 2018

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jul 20, 2018

Fixed a few bugs regarding intermediate angular velocity calculation and intermediate error calculation.

Added a tol variable to compute the error tolerance in angular velocity. Set default as 1e-12. Paper mentions 1e-6 but only as an estimate of number of iterations. In practice, the number of steps taken to reach convergence depends on the 'non-linearity' or the 'jaggedness' of the path which depends on the orientations, the times and the specified initial and final angular velocities.

On a jagged path generated by 10 random rotations, ~100 iterations yield an additional accuracy of 6 orders of magnitude. Have chosen default value of tol as 1e-12 to converge at ~200 iterations. Let me know if that seems acceptable.

@nmayorov
Copy link
Contributor

nmayorov left a comment

It looks OK, but we need to make the code easier to understand. Stronger tests are required.

@@ -15,13 +15,14 @@
Rotation
Slerp
Spline

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

Perhaps QSpline or QuaternionSplie or RotationSpline or something else to make it more specific.

@@ -1639,3 +1639,263 @@ def __call__(self, times):

return (self.rotations[ind] *
Rotation.from_rotvec(self.rotvecs[ind] * alpha[:, None]))


def _Rfunc(rotvec, omega):

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

I suggest to try and avoid names like this, because it's very hard to understand what this function does without reading the paper.

Can you come with some descriptive names.

This comment has been minimized.

@adbugger

adbugger Jul 23, 2018

Author Contributor

The thing is that the functions are quite simply replacements for long terms in the equations. I'm not sure if their purpose can be put in a succinct name. I'll go over the paper and try again.

r1 * e * dot_prod * np.cross(cross_prod, e))


def _Bfunc(rotvec, omega):

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

The same here.

rotations : `Rotation` instance
Rotations to perform the interpolation between. Must contain N
rotations.
omega_i : array_like, shape (3,), optional

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

I think we should provide an option to automatically determine some reasonable value for initial and final angular velocity.

One of simpler ways is to use "rotation vector / delta_time" for the first and last intervals. From experience, it creates some undesired boundary effects, but they are not that bad. I will be fine if you implement just that for now. For example omega_i='auto' or omega_i='const'.

omega_f : array_like, shape (3,), optional
Final angular velocity in radians at the time of the last orientation.
Default is `[0, 0, 0]`.
tol : float, optional

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

This is implementation detail which is not of interest for a user. If it works well with 1e-15 --- let's just use this value.

self.omega[-1] = omega_f
self.omega[1:-1] = w_int

# Formula for theta(t) from eqn 16, coefficients from page 5 of ref

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

Let's use http://scipy.github.io/devdocs/generated/scipy.interpolate.PPoly.html#scipy.interpolate.PPoly to store coefficients and do polynomial evaluations.

This comment has been minimized.

@adbugger

adbugger Jul 23, 2018

Author Contributor

This will require some re-formulation of the polynomial in the paper, I'll do that.

@@ -889,3 +889,80 @@ def test_slerp_call_time_out_of_range():
s([0, 1, 2])
with pytest.raises(ValueError, match="times must be within the range"):
s([1, 2, 6])


def test_spline_trivial():

This comment has been minimized.

@nmayorov

nmayorov Jul 23, 2018

Contributor

Add non-trivial tests.

Try to find a way to somehow check the continuity of the solution.

This comment has been minimized.

@adbugger

adbugger Jul 23, 2018

Author Contributor

How about calculating angular velocity (and angular acceleration?) close to the break points using a finite difference approximation?

This comment has been minimized.

@adbugger

adbugger Jul 23, 2018

Author Contributor

So I implemented this method and it looks like there are some bugs which still need to be ironed out. The angular velocities at key_times + delta_t agree to within 1e-3, but the velocities at key_times - delta_t do not. I'll keep looking into it.

This comment has been minimized.

@nmayorov

nmayorov Jul 24, 2018

Contributor

It might be not exactly a bug. What I suggest to do first is to implement angular velocity and acceleration computations (suggested in my email) and plot them as the function of time. It will be clear whether the behavior is correct. If it is correct, you will find a way to check it formally in tests.

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jul 25, 2018

So I've used scipy.interpolate.PPoly for various interpolations. One conceptual error I was making earlier was confusing theta_dot from the paper with the angular velocity vector. Further, calculating the angular velocity vector from quaternion rotations is not as straightforward as using the earlier finite difference approximation. For now I have just stored the coefficients for theta_dot and theta_dot2. Computing the angular velocity and angular acceleration from them requires some more work. I'm on it right now.

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jul 27, 2018

Figured out the bugs. Added a function for returning angular velocity and renamed to QSpline. Angular velocity continuity is now satisfied for smooth changes in orientation; for random rotations continuity at the right ends of the intervals is still sometimes an issue, however looking at the graphs in the paper it appears that problem cannot be completely removed.

@nmayorov the from scipy.interpolate import PPoly fails on CIs. What could be the issue there? The tests pass locally.

c2 = np.empty_like(c1)
c2[~small_theta] = (1 - np.cos(large_angle)) / large_angle
c2[small_theta] = small_angle / 2 - (small_angle ** 3) / 24

This comment has been minimized.

@ewmoore

ewmoore Jul 27, 2018

Member

The need to do this small vs. large angle calculation for c1 and c2 can probably be replaced by special.sinc (or np.sinc) and special.cosm1.

This comment has been minimized.

@nmayorov

nmayorov Jul 28, 2018

Contributor

@ewmoore nice suggestion!

Certainly can be used in a lot of places (including some previous code @adbugger wrote). In this particular case we should represent 1 - cos(x) as 2 sin^2 (x/2) and then sinc emerges.

One thing which is slightly annoying is that scipy/numpy defines sinc as sin(pi x) / (pi x), I would much prefer to have a clean definition sin x / x. I guess we may write a simple wrapper to alleviate that.

@adbugger investigate that idea. I think we can remove all small/large angle branches using sinc, here and in Rotation class as well.

This comment has been minimized.

@adbugger

adbugger Jul 30, 2018

Author Contributor

Will implement this in a separate PR.

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Jul 30, 2018

I think I figured out the source of the import error. The file scipy/interpolate/interpnd.pyx performs import scipy.spatial.qhull as qhull. This causes an import cycle when importing from scipy.interpolate within rotation.py.

Any ideas on how to work around it?

@rgommers

This comment has been minimized.

Copy link
Member

rgommers commented Jul 31, 2018

Any ideas on how to work around it?

Removing from . import transform in spatial/__init__.py may be enough. There's anyway a case that can be made for that.

Another option is to only import interpolate within QSpline. In general it's good if scipy submodules are not all dependent on each other; and here the use of interpolate within spatial is (a) new and (b) very limited.

return w, diff


class QSpline(object):

This comment has been minimized.

@rgommers

rgommers Jul 31, 2018

Member

How about calling this QuaternionSpline? Or QuaternionInterpolator? Or yet something else - just want to avoid giving the impression that this is a spline thingy that can be reused (it's very similar to B-spline now).

This comment has been minimized.

@adbugger

adbugger Jul 31, 2018

Author Contributor

How about renaming both functions to RotationSlerp and RotationSpline? In my mind this makes is explicit that these classes are for use only with Rotations, as opposed to general interpolators for quaternions.

This comment has been minimized.

@rgommers

rgommers Jul 31, 2018

Member

That sounds logical to me.

This comment has been minimized.

@nmayorov

nmayorov Jul 31, 2018

Contributor

RotationSlerp and RotationSpline sounds logical to me as well.

with pytest.raises(ValueError, match="times must be within the range"):
s([0, 1, 2])
with pytest.raises(ValueError, match="times must be within the range"):
s([1, 2, 6])

This comment has been minimized.

@tylerjereddy

tylerjereddy Jul 31, 2018

Contributor

Could reduce code duplication here with @pytest.mark.parametrize


vel = spline.omega[1]
assert_allclose(vel, omega_i)
assert_allclose(vel, omega_f)

This comment has been minimized.

@tylerjereddy

tylerjereddy Jul 31, 2018

Contributor

omega_i and omega_f are identical list objects in this unit test; why the assertion for vel compared to both as apposed to a single reference value? Are you trying to ensure that QSpline doesn't modify the omega argument objects passed in?

This comment has been minimized.

@adbugger

adbugger Jul 31, 2018

Author Contributor

You're right. Only one check is required here.

@adbugger adbugger force-pushed the adbugger:rotation_spline branch from b3c43a2 to ab2dc7e Jul 31, 2018

@adbugger adbugger force-pushed the adbugger:rotation_spline branch from 49ccc92 to 3a0db11 Jul 31, 2018

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Aug 1, 2018

@nmayorov , @rgommers is this mergable? I'll investigate the idea of using np.special.sinc in a separate PR.


return initial_rot * interpolating_quat

def angular_velocity(self, times):

This comment has been minimized.

@nmayorov

nmayorov Aug 2, 2018

Contributor

I think it should go to __call__ with some argument like order.

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Aug 2, 2018

@adbugger I'm not 100% satisfied with the code. I will try to review thoroughly on weekends. Meanwhile I think you should also implement computation of the angular acceleration (since we guarantee its continuity --- it will be interesting to at least see that).

I suggest the interface __call__(times, order=0)

@adbugger

This comment has been minimized.

Copy link
Contributor Author

adbugger commented Aug 2, 2018

@nmayorov I'll get on it. It's just that on 6 Aug the official coding period ends and the evaluation period starts (6 - 14 Aug). So, depending on the changes required, a review over the weekend might be cutting it a little too close.

I'll implement the __call__ interface and the acceleration functionality meanwhile.

@nmayorov

This comment has been minimized.

Copy link
Contributor

nmayorov commented Aug 5, 2018

Experimented a bit with the RotationSpline: https://gist.github.com/nmayorov/871380118428ae9e41bf6e4dd36cd35a

Computation of angular velocity and acceleration seems to be incorrect. There should not be any breaks in lines.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment