This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 560 - Allow standard linear algebra operations for quaternions
Allow standard linear algebra operations for quaternions
 Status: DECISIONNEEDED None Eigen Unclassified Geometry (show other bugs) unspecified All All Normal Feature Request Nobody 758 458 459 601 719 3.4 Show dependency tree / graph

 Reported: 2013-03-04 01:32 UTC by Boris D. 2019-12-04 12:07 UTC (History) 6 users (show) chtz gael.guennebaud hauke.heibel jacob.benoit.1 Keenan.Crane minorlogic

Attachments

 Boris D. 2013-03-04 01:32:09 UTC ```It seems it is currently not allowed by Eigen to multiply a quaternion by a scalar, or add quaternions together (simple component-wise operations). Indeed, if a quaternion is seen as a 3D rotation, these operations do not make sense, but quaternions could be used for other purposes, and even if it is used to represent rotations, algorithms might use such operations followed by normalization (my application case was an integration step in a physics simulator). I guess it should be easy to implement these: Quaternion & operator+(Quaternion & q); Quaternion & operator+=(Quaternion & q); Quaternion & operator*(Scalar s); Quaternion & operator*(Scalar s, Quaternion & q); Quaternion & operator*=(Scalar s); But I could not figure out how to do it myself, not being familiar with the Eigen's code source structures and patterns, hence directly used code like: q.w() += other.w(); q.x() += other.x(); q.y() += other.y(); q.z() += other.z(); That unfortunately does not take advantage of vectorization.``` Boris D. 2013-03-04 02:48:33 UTC ```I do not know how to create a patch, but these two methods added in the body defintion of QuaternionBase in Eigen/src/Geometry/Quaternion.h seem to work: ________________________________________________________________ /** \returns the addition of \c *this and \a other * \warning This operation has no direct geometric interpretation, and in the * general case result in a non unit quaternion. */ template inline Quaternion operator+(const QuaternionBase& other) const { return Quaternion(coeffs() + other.coeffs()); } /** \returns the multiplication of \c *this by the scalar \a s * \warning This operation has no direct geometric interpretation, and in the * general case result in a non unit quaternion. */ inline Quaternion operator*(Scalar s) const { return Quaternion(s * coeffs()); } ________________________________________________________________ that allows to write code like: q1 = q1 * s1 + q2 * s2; which is already much more readable than what I had previously. I do not know where I should write the implementation for s * q (instead of q * s) and didn't manage to implement operator+= ( coeffs() returns a const Coefficients and then can't be modified ).``` Gael Guennebaud 2013-03-04 09:35:36 UTC ```This has been discussed many times, and the problem to add such methods is that our Quaternion class represents unit quaternions. As a workaround, you can use the .coeffs() method, e.g.: q.coeffs() +=other.coeffs();``` Christoph Hertzberg 2013-03-04 09:50:16 UTC ```How about introducing a new class (e.g. GeneralQuaternion -- although eventually it might be better to call the current class UnitQuaternion and the new class just Quaternion)? Then UnitQuaternion might inherit from RotationBase and essentially hold a GeneralQuaternion inside, but only provide unit quaternion operations, while GeneralQuaterion also supports addition, exponentiation, etc.``` Boris D. 2013-03-04 14:27:56 UTC ```> q.coeffs() += other.coeffs(); Oh, this works indeed, thanks! This workaround is really fine for my particular usage (non unit operations are rare), but I tend to agree with the proposition of Christoph: it makes a lot of sense and would probably be useful for several users. By the way, I think exponentiation is already supported for unit quaternions in Eigen through slerp: q^t = q.slerp( t, <1,0,0,0> ). For instance, sqrt(q) represents "half" the rotation between a no-rotation and q, with the ambiguity in the case of a 180 degree angle (as for slerp).``` Christoph Hertzberg 2013-03-17 14:45:02 UTC ```(In reply to comment #4) > By the way, I think exponentiation is already supported for unit quaternions in > Eigen through slerp: q^t = q.slerp( t, <1,0,0,0> ). For instance, sqrt(q) > represents "half" the rotation between a no-rotation and q, with the ambiguity > in the case of a 180 degree angle (as for slerp). That would work for unit quaternions only, and I guess it is not the most efficient way. Actually, what I also had in mind was doing exp(q), as defined by the standard power series (of course implemented using a closed form solution). Further functions such as log, sin, cos, tan, etc would be possible. I'm not sure if there are actual use-cases for most of them, though.``` Boris D. 2013-03-17 23:24:46 UTC ```I may be wrong, but I think slerp is actually a pretty efficient way to compute the exponential of a (at least unit) quaternion. It involves computing one 4D dot product, then one acos and 3 sin of real numbers (the dot product, acos and one sin can be precomputed if slerp is used for different t but same vectors). I'm not fluent in numerical methods to compute power series, but my guess is just that computing 4 trigonometric functions of real numbers (well known and optimized power series using lookout tables) is faster than directly computing a power series of a 4D vector. Though, it has indeed the advantage of being more generic to allow other power series and non-unit quaternions, which I think is good and would always find some use-cases we don't think of.``` Christoph Hertzberg 2013-03-18 00:23:58 UTC ```(In reply to comment #6) > [...] is faster than directly computing a power series of a 4D vector. As I said, there are closed form solutions for these. E.g. quaternion power (assuming unit-quaternions for fair comparison) [imag; real]: [v; w]^t = [v*sin(t*acos(w))/sqrt(1-w*w); cos(t*acos(w))] Involving 1 acos, sin, cos and sqrt [sqrt(1-w*w) can be replaced by sin(acos(w))] and some scalar multiplications. When only t changes one sin and one cos (of the same angle) need to be calculated). Of course similar safe-guards as for SLERP need to be done for w>1-eps. For non-unit quaternions it basically requires an atan or atan2 instead of acos and an additional real power.``` Boris D. 2013-03-18 00:47:09 UTC ```Oh sorry, I indeed missed "of course implemented using a closed form solution"! Then this closed form is exactly the way slerp computes [v; w]^t, expressed as slerp between [0; 1] and [v, w] (with very slight optimizations due to the known [0; 1] component that of course would have been done), then we are talking about the same thing. To me, a closed form for non-unit quaternions would be called a generalization of slerp, but it's just a matter of terminology ;-)``` Christoph Hertzberg 2014-06-15 04:11:26 UTC ```It is currently not obvious if Quaternion allows non-unit quaternions. I'd like to bring up the discussion of introducing new classes GeneralQuaternion and UnitQuaternion to make the distinction clear. The latter can still exploit all current optimizations, the former can provide additional functionality, such as quaternion addition or _rotateAndScaleVector as suggested in Bug 459.``` Gael Guennebaud 2014-06-16 13:56:40 UTC `I'd rather generalize the current Quaternion class, remove the unitary assumptions, mark the functions related to unitary-rotation obsolete (and lacking efficiency) and then add a UnitQuaternion class or maybe something more general relying on Quaternion to represent rotations or rotations+uniform_positive_scaling. The Quaternion would then be very low level, more like std::complex.` Boris D. 2014-06-16 20:07:07 UTC `For what my opinion is worth: I personally prefer the option suggested by Gael, the class names seem more intuitive this way. :-)` Christoph Hertzberg 2014-06-17 10:34:46 UTC `I'm totally ok with letting our current Quaternion class mean "generalized quaternion". I was merely afraid that we introduce API incompatibilities by that.` Michael Norel 2014-06-20 12:24:08 UTC ```My small comment on quaternion unity. Here i assume, quaternion is used to represent rotation. From my experience, introducing the "unit" contract leads to lot of problems. Let's look at example: 1) We keeps rotation transformation in quaternion. Making some modeling of physics or visualization, we update this orientation with other "deltaRotation" quaternion. Every multiplication add some noise to quaternions and its length comes away from 1.0. 2) Than we call rotationMatrix method that assumes the length 1.0 , and our rotation matrix is distorted (cause assumes unit length). 3) making some manipulation with matrix we convert it back to quaternion , assuming that it is ortonormal. And have quaternion length even farer from 1.0. 4) This is very hard to debug bug with current code. We don't check the invariants of input data (matrix) the quaternion unity. So , to avoid such kind of errors we Must keep quaternion unity invariant after every modification (multiplication). It can be done inside Quaternion class, or outside manually. This solution is more expensive than use of nonunit quaternions. Nonunit quaternion require some careful modifications to code. Main is matrix extraction and vector rotation, that should be weighted by inverse squaredNorm of quaternion. Also quaternion interpolation need additional normalization (but i prefer separate objects "interpolators" that cashes all pre calculations). My main use of quaternion class in computer vision software, and i prefer to use nonunit. But can't give general recommendation.``` Christoph Hertzberg 2014-06-20 15:54:09 UTC ```We could avoid adding another class, by providing 3 variants for each rotation-related method: rotate: res = q*[0;v]*q.inverse(); // what Michael and many others want rotateAndScale: res = q*[0;v]*q.conj(); // cf Bug 459 rotateAssumeUnity: // What we do know: Undefined behavior for non-unit quaternions! // To get the expected behavior, normalize() must be called once in a while Then we'd only need to decide what shall be the default behavior, or perhaps make it configurable using an Options-flag. @Michael, about the "interpolators": If you want this, please open a new bug.``` Michael Norel 2014-06-20 17:02:25 UTC ```Agree with Christoph Hertzberg I prefer nonunit versions for default rotation-related method. My explanation. Default , common functions should be save with minor performance penalty, and if user is friendly with quaternions and knows that quaternion is unit, he can write optimized code with rotateAssumeUnity. Also solution can come from Eigen design, default prefer speed or robustness. Another tip is to implement runtime, debug checks for quaternion unity in danger methods.``` Gael Guennebaud 2014-06-20 17:07:06 UTC ```I would make > rotate: > res = q*[0;v]*q.inverse(); the default, and offer the other options through subclasses (or a template option like Transform?). It's important to keep a unified way to apply a transformation to a vector.``` Christoph Hertzberg 2014-07-03 18:24:04 UTC ```I'm afraid allowing Scalar*Quaternion will lead to unexpected behavior because multiplication wouldn't be associative anymore: Quateriond q; Vector3d v; Vector3d w1= (-q)*v, w2= -(q*v), w3 = (-2*q)*v, w4=-2*(q*v); Note that in w1 and w3 the parentheses only emphasize the default precedence, but w2 and w4 are likely the intended behavior. The current behavior for Scalar*Quaternion is definitely buggy. See Bug 758.``` Keenan Crane 2016-06-10 15:10:47 UTC `I will voice another vote for including standard vector space operations in the Quaternion class. There are of course many uses of quaternions beyond rotations. One that may not be immediately obvious is quaternionic linear algebra, i.e., matrices with quaternion-valued entries (as used in, e.g., 3D geometry processing, conformal geometry, and particle physics). It is difficult to add this functionality to traditional linear algebra libraries, which typically support only real and complex. Since Eigen is more flexible about data types it could potentially provide some unique value here---if quaternions were given vector space operations. (Note that representing quaternionic matrices as real matrices is not a terrific solution, since this artificially inflates the storage/bandwidth cost by 4x.)` Michael Norel 2016-06-18 20:49:36 UTC ```Also can be introduced classes: "QuaternionRotation" "UnitQuaternionRotation" just for quaternions representing rotation, and restricted some math that brokes rotation invariants and "Quaternion" have all math``` Christoph Hertzberg 2016-06-21 13:43:29 UTC ```I'm fine with either name suggestion for unit Quaternions. And I think, eventually, the best name for the generic class is simply "Quaternion". The question is how do we transition to that. Most importantly, do we allow Quaternion * Vector3? This will be problematic when the Quaternion is scalable (cf comment 17). If we disallow Quaternion * Vector3, we will break code. We could deprecate that operator and suggest using the UnitQuaternion class instead, or provide dedicated rotate/rotateAndScale methods (cf Bug 459) -- this would still influence a lot of existing code. Also, shall we do the first step of transition in 3.3? I think, if we want to deprecate methods, we should do that in 3.3.0, if possible. The easiest solution would be just to introduce new Quaternion classes under new names and keep the current Quaternion class more or less as it is (and suggest in the documentation to transition to QuaternionRotation or to GeneralizedQuaternion) -- but we would waste the name Quaternion then. Another alternative would be to use several #ifdef switches to choose the wanted behaviour (this will be bad for linking with different libraries, however).``` Christoph Hertzberg 2016-06-21 13:47:58 UTC `Regarding comment 18: Quaternionic linear algebra sounds interesting -- we'd need to refactor quite some code for that, I guess. For example, we make some optimizations for (scalar * matrix * scalar * matrix * scalar) products, which only work as long as multiplication is commutative.` Keenan Crane 2016-06-21 16:10:51 UTC ```If I understand comment 17 correctly, w3 and w4 generate different results because the semantics of Quaterniond::operator*(Vector3d) are to treat the action of a quaternion on a vector as a rotation. In other words, 2(q*v) means "use q to rotate v, then scale the result by 2" whereas (2*q)*v means "use 2q to rotate v", which is likely the same as 4*(q*v) (since q appears twice in the formula for the rotation of a vector by a quaternion). For those of us who work a lot with quaternions, this interpretation of "q*v" is confusing, and violates the "principal of least surprise." It also seems to be at the root of the problems described in comment 17 and elsewhere. If one is going to allow quaternion-vector multiplication, a more standard interpretation is that any vector v=(a,b,c) is naturally represented as the quaternion 0 + ai + bj + ck. With this interpretation, the expression q*v is just a multiplication of two quaternions; to represent a rotation, one would write q*v*conj(q). (The operation q*v without the trailing conj(q) also has a standard interpretation in geometry/physics: for unit quaternions it is the action of Spin(3) on R^3.) Perhaps a safer way to handle all this in Eigen is to have a method Vector3d Quaterniond::rotate(Vector3d v) that is responsible for rotation. Direct multiplication of quaternions with vectors could be forbidden, except via explicit conversion. E.g., there could be an explicit type conversion method like Vector3d::toQuat() that converts a 3-vector to an imaginary quaternion, as described above. In this way, q * v.toQuat() has clear semantics, as does q.rotate(v).``` Christoph Hertzberg 2016-06-23 14:41:53 UTC ```Yes, the problem is that the current Quaternion class essentially was introduced only as representation of a SO(3) [3d-rotation] object. For this use-case it makes perfect sense to directly allow multiplication with 3d vectors (or with other Transformation objects). Requiring to write q.rotate(v) instead of q*v will make some code unnecessarily verbose. However, the naming is admittedly confusing and something like QuaternionRotation probably would be have been better. Ideally, Scalar*QuaternionRotation will return an expression which is not assignable to a Quaternion, but assignable to a Matrix3 or multipliable by another Vector3 (and resulting in s*q.rotate(v)) We can definitely not suddenly remove (or even worse: re-define) the Quaternion*Vector3 operator, because we would break a lot of existing code. At the moment, I see basically only the possibility to introduce one or two new Quaternion classes (or Option flags for the existing class) which have the expected behaviour of either generalized Quaternions or SO(3) objects and at some point deprecate the existing methods and operator*(Vector3) overload.``` Michael Norel 2016-06-23 15:02:19 UTC ```(In reply to Christoph Hertzberg from comment #23) SO(3) - is very useful entity for 3dmath and computer vision, but quaternion is only one of possible way to represent rotation space? Do you want to use only UNIT quaternions? (because a lot of operations can be performed faster with NONunit quaternions)``` Keenan Crane 2016-06-23 16:00:33 UTC ```In response to comment #23: Makes sense. One quick path forward might be: - Retain Quaternion*Vector3 as a way of expressing rotations. - Add vector space operations to Quaternion, including scalar multiplication. - Add a run-time check in Quaternion*Vector3 that the quaternion has unit norm (only in debug mode, perhaps). - Add the warning below (which can be removed in a later version): Quaternion operator*( const Scalar&, const Quaternion& ) { #warning "Please note that Eigen now supports scalar-quaternion multiplication. For quaternion rotation, expressions of the form Scalar*Quaternion*Vector should be made explicit via Scalar*(Quaternion*Vector)" } This path would still break some code, though. A slower but safer alternative might be: 1. Add a QuaternionRotation class that is basically typedef'd to the current Quaternion class. 2. Add a GeneralQuaternion class that implements vector space operations. 3. At some future date, deprecate Quaternion. 4. At some even later future date, deprecate GeneralQuaternion, replacing it with Quaternion. > Requiring to write q.rotate(v) instead of q*v will make some code unnecessarily verbose. Perhaps, though I notice that dot and cross products in Eigen are expressed as u.dot(v) and u.cross(v) rather than, say, u*v and u^v. So perhaps q.rotate(v) would provide some uniformity with the rest of the library - at very least, it could be good to include this as one possible syntax, i.e., users can write either q*v or q.rotate(v), if they want to be more explicit.``` Gael Guennebaud 2016-06-23 21:31:42 UTC ```Ok, let's resolve this issue for 3.3, it's standing for too long. In my opinion Christoph's suggestion about adding a template parameter controlling the semantic of Quaternion should allow us to reach a good compromise: - By default, it would be set to "RotationSemantic", and only the current set of methods would be exposed (for backward compatibility). - To get a "general" quaternion, one would pass something like "GeneralSemantic", (in which case all rotation-oriented specificities would be dropped). - The key advantage of this approach is that it permits to make the default behavior configurable through a preprocessor directive, just like Matrix can default to row-major. There won't be linking issues, because the types will remain different. - Finally, in c++11, one can easily make convenient template typedefs like: Quat <-> Quaternion QuatRot <-> Quaternion``` Gael Guennebaud 2019-02-22 16:38:54 UTC ```See PR595: https://bitbucket.org/eigen/eigen/pull-requests/595/add-a-general-quaternionnumber/diff``` Gael Guennebaud 2019-03-02 12:45:53 UTC ```The aforementioned PR generated extended discussions that yielded to a simpler design and a new PR: https://bitbucket.org/eigen/eigen/pull-requests/600 So far, this new PR simply extends the current Quaternion and seamlessly disable dangerous expressions: s*q*v, q*v*s, -q*v with s a scalar, q, a quaternion, and v a 3D vector.``` Gael Guennebaud 2019-03-09 12:56:46 UTC ```Still remains to conclude on the default behavior of q*v, to toRotationMatrix, and slerp. Let's start with q*v. Making it a fast alias for q*(0,v)*conj(q) was appealing to me at first but: - q*v comes from RotationBase, so making it applies scaling makes little sense - Using quaternion for rotation+scaling is not very handy because you have to scale q by sqrt(scale). The unique user requesting for rotation+scale rather wanted q.norm()==scale... (this has the advantage of making s*q*v==s*(q*v) though). Overall, I would not even attempt to leverage any facilities for quaternion as rotation+scale. So it remains to keep it as is or to perform implicit normalization, and I would strongly advocate for the later: - This is the opinion expressed in this thread. - q*v with implicit normalization is faster than q.normalized()*v (no sqrt) - The overhead is not that high. - We could still optimize existing q.normalized()*v expression if needed (to avoid double normalization). - This choice makes the most sense with q*v being defined in RotationBase. Then, toRotationMatrix should really do implicit normalization as well. What is less clear to me is how to leverage the existing fast paths. For toRotationMatrix we could simply add an option as a second parameter, but this cannot be done the same way for q*v. For``` Michael Norel 2019-03-11 08:38:05 UTC ```Regarding "q*v with implicit normalization" vs "q.normalized()*v" and scale policy. Sorry, i can have leak of current develop information. 0. I strongly prefer use of "SCALE_POLICY" for quaternions. It can have both advantages from unit and nonunit quaternions and operations. 1. There was idea that quaternion modification can be done less often than rotation transform with vector. So implicit normalization can lose performance. Since several years, i changed opinion about "implicit normalization". 2. Consider idea to keep quaternion unit scale as class invariant. Having different normalization policy. For me , it is legacy to use quaternion for rotation as unit , but without internal invariant of unity (we can see a lot of hard to discover bugs). And , as for now , without even runtime check for unity. 3. Having scale policy, we can use nonunit quaternions with implicit normalization during rotation and a lot of different operations faster than unit. It is creation from different representation of rotation, slerp, from 2 vectors and many others. A lot of quaternion creation operations faster for nonunit. But we still have correct functionality for rotation (with implicit normalization). 4. We can have simple conversions between Unit and Nonunit quaternions. Also it prevents a lot of potential bugs. By default, developers can use Unit quaternions, with strong unity invariants, And nonunit with all of its benefits for advanced use. For example, all the temporal calculations can be done with NONunit quaternions, but result is stored with unit. 5. UNIT_SCALE_POLICY pro - auto normalization, remove huge pack of hard bugs. - maximum speed of vector transformation - possible some tricks for normalization ONCE per expression template collapse - some operations can be implemented as operation of nonunit and than cast to unit(Slerp, conversions, multiplications). So it is easy to implement neg - performance loose for implicit normalization of all quaternion modifications. NONUNIT_SCALE_POLICY pro - a lot of faster operations neg - not save for beginner developers - vector rotation (matrix creation) operation require additional implicit scale``` Gael Guennebaud 2019-03-13 10:23:36 UTC ```Thank you for your feedback, however I'm not sure to clearly follow you. Maybe my question was not clear enough. In the PR linked in comment 28, Quaternion is generalized to a general purpose non-unit quaternion. There is no plan at all for unit norm invariant, nor to have two quaternions classes with different invariants (too complicated API-wise, raise more questions than it solve, very confusing, etc.). So to move forward, let's consider this to be established. The open questions are only about what to do by default for operations that currently assumes a quaternion with a unit-norm as input: q*v, toRotationMatrix, slerp. 1) q*v, three options: 1.a) no change, i.e., undefined behavior if q.norm()!=1 1.b) pure rotation, i.e., implicit normalization, i.e.: make q*v <=> q*(0,v)*conj(q)/q.squaredNorm() 1.c) make q*v <=> q*(0,v)*conj(q) Performance-wise, the difference between 1.b) and 1.c) is small. Option 1.b) is also faster than 1.a) plus calls to q.normalized() (no sqrt). For the reasons listed in comment 29, I suggested 1.b), but from you comment I'm not sure which option you would suggest! If that's a concern, we could optimize 1.b) to fall-back to 1.a) for people who have already written stuff like: q.normalized()*v. If we go with 1.b), then it remains to find a way to expose 1.a, (and if we go with others, we also need a way to expose 1.b ), possibilities: q.assumeUnit()*v; assumeUnit(q)*v; q.apply(v, AssumeUnit); All approaches can be applied to toRotationMatrix and slerp. The first two could be more easily applied to other functions and ctors assuming unit vectors as input. One we converge on q*v, handling toRotationMatrix and slerp should be easier.``` Michael Norel 2019-03-13 12:04:13 UTC ```"There is no plan at all for unit norm invariant, nor to have two quaternions classes with different invariants (too complicated API-wise, raise more questions than it solve, very confusing, etc.)." It is sad to read. Cause for me, additional parametrization of quaternion class solves a lot of problems. Quaternion 1. safety 2. performance 3. compatibility 4. easy to implement I can share my vision if you ready to talk about. For me it looks as best way to deal with quaternion representing rotation. //======================================================== "1) q*v, three options: 1.a) no change, i.e., undefined behavior if q.norm()!=1 1.b) pure rotation, i.e., implicit normalization, i.e.: make q*v <=> q*(0,v)*conj(q)/q.squaredNorm() 1.c) make q*v <=> q*(0,v)*conj(q)" 1.a) Is a worst way, cause it is NOT safe. If you will go this way, i prefer to add , runtime checks for unity. Cause this operation can produce invalid matrix. Even more, all operation assuming unit magnitude, should do runtime checks for unity (in debug). 1.b) vs 1.c) Current "q*v" is transformation of 3d point by rotation represented as quaternion. If this operator will stay as part of interface , it should perform as q.toRotationMatrix() * v. And off cause, toRotationMatrix() must do 1.b) I prefer 1.b). 1. It is more safe than 1.a) and 1.c) in case of scale drift. 2. It is rotation by quaternion, and have a lot legacy code and understanding. It is more compatible than 1.c).``` Nobody 2019-12-04 12:07:57 UTC ```-- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/560.```

 Note You need to log in before you can comment on or make changes to this bug.