Numerical code audit - #363
Closed
AlexanderFabisch wants to merge 7 commits into
Closed
Conversation
slerp_weights checked angle == 0.0 with exact equality. For antipodal quaternions (q and -q) angle_between_vectors returns pi, and sin(pi) ~= 1.2e-16 in float64, causing the ratio sin(t*pi)/sin(pi) to overflow to ~8e15. Replace the exact-zero guard with abs(sin(angle)) < eps so that both the near-zero and the near-pi (antipodal) degenerate cases fall back to a linear blend instead of producing +/-Inf. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
The previous threshold, np.finfo(float).eps (~2.2e-16), was far too low. For any theta in [eps, sqrt(6*eps)] the closed-form expressions (1 - sin(t)/t) and (1 - t/(2*tan(t/2))) suffer catastrophic cancellation: the correction is smaller than half a ULP (unit in the last place, i.e. the spacing between adjacent floats) of 1.0, so the subtraction rounds to exactly 0 and all significant bits are lost. The correct threshold comes from the tightest problematic term, the inverse Jacobian coefficient 1 - t/(2*tan(t/2)) ~= t^2/12 for small t. Requiring t^2/12 >= eps/2 (half a ULP of 1.0, where eps = ULP(1.0) = 2^-52 ~= 2.2e-16 for float64) gives t >= sqrt(6*eps) ~= 3.65e-8. This threshold also covers the forward Jacobian (1 - sin(t)/t ~= t^2/6). The 10-term Taylor series has truncation error O(t^20) at this threshold, which is ~1e-150, far below float64 precision. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Both estimate_gaussian_rotation_matrix_from_samples and estimate_gaussian_transform_from_samples called np.cov(..., bias=True), which divides by N rather than N-1. Their docstrings state the functions use an unbiased estimator as outlined by Eade (2017), which uses N-1. Changed bias=True to bias=False (the numpy default) to match the documented behaviour and Eade's formulation. The covariance was previously underestimated by the factor N/(N-1), most visible for small sample counts. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Two issues with the omega_norm validation: 1. Dead zone: line 126 used abs(omega_norm) < eps to decide "near zero" while line 120 used abs(omega_norm) > 10*eps to decide "not near zero". For omega_norm in (eps, 10*eps) neither condition was true, so the v_norm check was silently skipped and any v passed validation. Fixed by using the same threshold (10*eps) in both places. 2. Near-one tolerance: the checks for abs(omega_norm - 1) and abs(v_norm - 1) used 10*eps with no stated justification. IEEE 754 first-order analysis of np.linalg.norm on a 3-vector shows the accumulated error is at most 5u/2 = 5/4*eps (u = eps/2 unit roundoff). 2*eps gives a factor-of-1.6 safety margin over that bound and is used for all near-one checks. Both thresholds are extracted into named constants _NEAR_ZERO_TOL and _NEAR_ONE_TOL to make the distinction explicit and prevent future drift. Adds a regression test: omega_norm = 3*eps (in the former dead zone) with v_norm = 0.5 must raise ValueError. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
norm_angle maps to (-pi, pi], and the negation branch only flips angles below zero, leaving pi unchanged. The actual output range is [0, pi]; the pi endpoint is included and represents a valid 180-degree rotation. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
….., 3)
np.einsum('nii', ...) produces one scalar trace per matrix, so the
precomputed traces array has shape (...), not (..., 3).
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
clementperon
reviewed
Jun 9, 2026
| a : array, shape (..., 4) | ||
| Axis of rotation and rotation angle: (x, y, z, angle). The length | ||
| of the axis vector is 1 and the angle is in [0, pi). No rotation | ||
| of the axis vector is 1 and the angle is in [0, pi]. No rotation |
There was a problem hiding this comment.
Was reading the doc and about to open a MR for this issue 👍
There was a problem hiding this comment.
As there could be two representation for this maybe it would be nice to force a deterministic axis sign no?
Member
Author
There was a problem hiding this comment.
Hi, yes, this makes sense. Thanks for the suggestion. What would you suggest as a solution? Would you like to open a pull request to the develop branch? We also have to make sure that rotations.norm_axis_angle is consistent.
There was a problem hiding this comment.
The rule I have set in my repo is
"
At angle = 180°, axis and -axis represent the same rotation. To make the representation unique, the first non-zero component of axis SHALL be positive.
"
Member
Author
There was a problem hiding this comment.
Would you like to open a pull request for this? If not, I will do it when I have time.
There was a problem hiding this comment.
I would be glad to contribute but I have finally implemented it natively without using pytransform3d lib (that's the only feature I was looking for)
As I have zero knowledge in this repo, sure an AI agent will do better than me 😄
Member
Author
|
Merged with 1a9eb5b |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.