| Subject: | creating Quaternions based on two vectors |
Thanks for creating Math::Quaternion. Very handy.
Versions: Math::Quaternion 0.03
Perl and OS: irrelevant
When creating a new Quaternion based on two vectors, the code crashed on
a divide-by-zero when the two vectors are parallel. The underlying cause
is that the rotation vector obtained from two parallel vectors through
their cross product is the null vector, which can not be normalised.
my $quat = Math::Quaternion->new({ v1 => [ 1,1,1],
v2 => [ 2,2,2],
});
# Illegal division by zero at Math/Quaternion_orig.pm line 685.
However, it is possible to generate a valid Quaternion corresponding to
a zero-angle rotation, thereby avoiding a hard crash, and allowing the
calling program to continue (applying such a Quaternion to do a rotation
has the desired effect of not rotating the vectors).
In the attached diff file, two solutions are proposed, and the most
elaborate one is implemented. This generates a Quaternion with a zero
degree rotation angle, and a rotation vector perpendicular to the
original vectors (but in a random direction, as that can not be defined)
The attached diff file also proposed a few extra checks for invalid
input data.
I marked the issue "important", as it crashes the library for a
situation that arises frequently in general code where a vector can be
mapped on itself, or where the vectors were already aligned because
that's how the data happens to be.
| Subject: | Quaternion.diff |
--- Quaternion_orig.pm 2011-02-03 00:09:18.000000000 +0100
+++ Quaternion.pm 2011-02-03 22:51:09.000000000 +0100
@@ -647,6 +647,11 @@
# Both args references to vectors
my ($ax,$ay,$az)=@{$_[0]};
my ($bx,$by,$bz)=@{$_[1]};
+ if ( (($ax == 0) and ($ay == 0) and ($az == 0)) or
+ (($bx == 0) and ($by == 0) and ($bz == 0)) ) {
+ croak("Math::Quaternion::rotation() passed zero-length vector");
+ }
+
# Find cross product. This is a vector
# perpendicular to both argument vectors,
# and is therefore the axis of rotation.
@@ -659,6 +664,38 @@
my $mod2 = sqrt($bx*$bx+$by*$by+$bz*$bz);
# Find the angle of rotation.
$theta=Math::Trig::acos($dotprod/($mod1*$mod2));
+
+ # Check for parallel vectors (cross product is zero)
+ if (($x == 0) and ($y == 0) and ($z == 0)) {
+ # Vectors a and b are parallel, such that rotation vector is the zero-length vector (0,0,0), with theta 0.
+ # To remove round-off errors in theta, set it to 0
+ $theta = 0;
+
+ # Such a zero-length rotation vector is annoying (e.g. division by 0 on normalization, and problems combining rotations)
+ # Simple solution would be to set a random rotation vector (e.g. 1,0,0) to go with the zero rotation angle. This
+ # would satisfy as a quaternion that rotates the first vector on the second, by a zero degree rotation.
+ #
+ # A more elegant solution is to select a random rotation vector that is also perpendicular to
+ # both parallel vectors a and b. This satisfies the rotation requirement, and helps programs relying on the logic
+ # that the rotation vector has to be perpendicular to both vectors given (even if there are an infinite amount of
+ # rotation vectors that would satisfy that condition).
+ # Algorithm: Find a random vector b at any non-zero angle to vector a. One of the main axis will do. To reduce round-off errors, make
+ # b as perpendicular as possible to a by selecting one of the smallest components of vector a as the main component of b. This
+ # also avoid accidentally selecting a vector parallel to a
+ if ( (abs($ax) <= abs($ay)) and (abs($ax) <= abs($az)) ) {
+ ($bx,$by,$bz)=(1,0,0);
+ } elsif ( (abs($ay) <= abs($ax)) and (abs($ay) <= abs($az)) ) {
+ ($bx,$by,$bz)=(0,1,0);
+ } else {
+ ($bx,$by,$bz)=(0,0,1);
+ }
+ # Then, take the cross product between vector a and the new vector b, to generate some vector exactly perpendicular to vector a
+ # and hence also perpendicular to the original vector b (i.e. @{$_[1]})
+ $x = $ay*$bz-$az*$by;
+ $y = $az*$bx-$ax*$bz;
+ $z = $ax*$by-$ay*$bx;
+ # ($x,$y,$z) is now a random yet valid rotation vector perpendicular to the two original vectors. theta is still 0.
+ }
} else {
# 0 is a ref, 1 is not.
$theta = $_[1];
@@ -682,6 +719,9 @@
}
my $modulus = sqrt($x*$x+$y*$y+$z*$z); # Make it a unit vector
+ if ($modulus == 0) {
+ croak("Math::Quaternion::rotation() passed zero-length rotation vector");
+ }
$x /= $modulus;
$y /= $modulus;
$z /= $modulus;