Math: add APIs for quaternion reflection.

It isn't really useful for anything in practice, but I'm putting it
there for learning purposes because the topic comes up every now and
then.

doc/changelog.dox

@@ -207,6 +207,10 @@ See also:
     @relativeref{Math::Intersection,pointSphere()}, which are just wrappers
     over trivial code but easier to discover
 -   Added an unary @cpp operator+() @ce to all @ref Math classes
+-   Added @ref Math::Quaternion::reflection() and
+    @ref Math::Quaternion::reflectVector(), but mainly just for documentation
+    purposes as reflections cannot be combined with rotations and thus are
+    mostly useless in practice
 
 @subsubsection changelog-latest-new-materialtools MaterialTools library
 

src/Magnum/Math/Matrix4.h

@@ -244,8 +244,8 @@ template<class T> class Matrix4: public Matrix4x4<T> {
          * @cpp Matrix4::scaling(Vector3::yScale(-1.0f)) @ce. @f[
          *      \boldsymbol{A} = \boldsymbol{I} - 2 \boldsymbol{NN}^T ~~~~~ \boldsymbol{N} = \begin{pmatrix} n_x \\ n_y \\ n_z \end{pmatrix}
          * @f]
-         * @see @ref Matrix3::reflection(), @ref Vector::isNormalized(),
+         * @see @ref Quaternion::reflection(), @ref Matrix3::reflection(),
-         *      @ref reflect()
+         *      @ref Vector::isNormalized(), @ref reflect()
          */
         static Matrix4<T> reflection(const Vector3<T>& normal);
 

src/Magnum/Math/Quaternion.h

@@ -290,6 +290,25 @@ template<class T> class Quaternion {
          */
         static Quaternion<T> rotation(Rad<T> angle, const Vector3<T>& normalizedAxis);
 
+        /**
+         * @brief Reflection quaternion
+         * @param normal        Normal of the plane through which to reflect
+         * @m_since_latest
+         *
+         * Expects that the normal is normalized. @f[
+         *      q = [\boldsymbol n, 0]
+         * @f]
+         * Note that reflection quaternions behave differently from usual
+         * rotations, in particular they *can't* be concatenated together with
+         * usual quaternion multiplication, @ref toMatrix() will *not* create a
+         * reflection matrix out of them and @ref transformVector() will *not*
+         * do a proper reflection either, you have to use @ref reflectVector()
+         * instead. See its documentation for more information.
+         * @see @ref Matrix4::reflection(), @ref Vector::isNormalized(),
+         *      @ref reflect()
+         */
+        static Quaternion<T> reflection(const Vector3<T>& normal);
+
         /**
          * @brief Create a quaternion from a rotation matrix
          *
@@ -659,8 +678,11 @@ template<class T> class Quaternion {
          * quaternions. @f[
          *      v' = qvq^{-1} = q [\boldsymbol v, 0] q^{-1}
          * @f]
+         * Note that this function will not give the correct result for
+         * quaternions created with @ref reflection(), for those use
+         * @ref reflectVector() instead.
          * @see @ref Quaternion(const Vector3<T>&), @ref vector(),
-         *      @ref Matrix4::transformVector(),
+         *      @ref reflectVector(), @ref Matrix4::transformVector(),
          *      @ref DualQuaternion::transformPoint(),
          *      @ref Complex::transformVector()
          */
@@ -689,6 +711,31 @@ template<class T> class Quaternion {
          */
         Vector3<T> transformVectorNormalized(const Vector3<T>& vector) const;
 
+        /**
+         * @brief Reflect a vector with a reflection quaternion
+         * @m_since_latest
+         *
+         * Compared to the usual vector transformation performed with
+         * rotation quaternions and @ref transformVector(), the reflection is
+         * done like this: @f[
+         *      v' = qvq = q [\boldsymbol v, 0] q
+         * @f]
+         * You can use @ref reflection() to create a quaternion reflecting
+         * along given normal. Note that it's **not possible to combine
+         * reflections and rotations with the usual quaternion multiplication.
+         * Assuming a (normalized) rotation quaternion @f$ r @f$, a combined
+         * rotation and reflection of vector @f$ v @f$ would look like this
+         * instead: @f[
+         *      v' = rqvqr^{-1} = rqvqr^* = rq [\boldsymbol v, 0] qr^*
+         * @f]
+         * See also [quaternion reflection at Euclidean Space](https://www.euclideanspace.com/maths/geometry/affine/reflection/quaternion/index.htm).
+         * @see @ref Quaternion(const Vector3<T>&), @ref vector(),
+         *      @ref Matrix4::transformVector()
+         */
+        Vector3<T> reflectVector(const Vector3<T>& vector) const {
+            return ((*this)*Quaternion<T>{vector}*(*this)).vector();
+        }
+
     private:
         #ifndef DOXYGEN_GENERATING_OUTPUT
         /* Doxygen copies the description from Magnum::Quaternion here. Ugh. */
@@ -786,6 +833,12 @@ template<class T> inline Quaternion<T> Quaternion<T>::rotation(const Rad<T> angl
     return {normalizedAxis*std::sin(T(angle)/2), std::cos(T(angle)/2)};
 }
 
+template<class T> inline Quaternion<T> Quaternion<T>::reflection(const Vector3<T>& normal) {
+    CORRADE_DEBUG_ASSERT(normal.isNormalized(),
+        "Math::Quaternion::reflection(): normal" << normal << "is not normalized", {});
+    return {normal, 0.0f};
+}
+
 template<class T> inline Quaternion<T> Quaternion<T>::fromMatrix(const Matrix3x3<T>& matrix) {
     /* Checking for determinant equal to 1 ensures we have a pure rotation
        without shear or reflections.

src/Magnum/Math/Test/QuaternionTest.cpp

@@ -94,6 +94,8 @@ struct QuaternionTest: Corrade::TestSuite::Tester {
 
     void rotation();
     void rotationNotNormalized();
+    void reflection();
+    void reflectionNotNormalized();
     void angle();
     void angleNormalizedButOver1();
     void angleNotNormalized();
@@ -121,6 +123,7 @@ struct QuaternionTest: Corrade::TestSuite::Tester {
     void transformVector();
     void transformVectorNormalized();
     void transformVectorNormalizedNotNormalized();
+    void reflectVector();
 
     void strictWeakOrdering();
 
@@ -177,6 +180,8 @@ QuaternionTest::QuaternionTest() {
 
               &QuaternionTest::rotation,
               &QuaternionTest::rotationNotNormalized,
+              &QuaternionTest::reflection,
+              &QuaternionTest::reflectionNotNormalized,
               &QuaternionTest::angle,
               &QuaternionTest::angleNormalizedButOver1,
               &QuaternionTest::angleNotNormalized,
@@ -206,6 +211,7 @@ QuaternionTest::QuaternionTest() {
               &QuaternionTest::transformVector,
               &QuaternionTest::transformVectorNormalized,
               &QuaternionTest::transformVectorNormalizedNotNormalized,
+              &QuaternionTest::reflectVector,
 
               &QuaternionTest::strictWeakOrdering,
 
@@ -517,6 +523,22 @@ void QuaternionTest::rotationNotNormalized() {
     CORRADE_COMPARE(out.str(), "Math::Quaternion::rotation(): axis Vector(-1, 2, 2) is not normalized\n");
 }
 
+void QuaternionTest::reflection() {
+    Vector3 axis(1.0f/Constants<Float>::sqrt3());
+    Quaternion q = Quaternion::reflection(axis);
+    CORRADE_COMPARE(q.vector(), axis);
+    CORRADE_COMPARE(q.scalar(), 0.0f);
+}
+
+void QuaternionTest::reflectionNotNormalized() {
+    CORRADE_SKIP_IF_NO_DEBUG_ASSERT();
+
+    std::ostringstream out;
+    Error redirectError{&out};
+    Quaternion::reflection({-1.0f, 2.0f, 2.0f});
+    CORRADE_COMPARE(out.str(), "Math::Quaternion::reflection(): normal Vector(-1, 2, 2) is not normalized\n");
+}
+
 void QuaternionTest::angle() {
     auto a = Quaternion({1.0f, 2.0f, -3.0f}, -4.0f).normalized();
     auto b = Quaternion({4.0f, -3.0f, 2.0f}, -1.0f).normalized();
@@ -916,6 +938,25 @@ void QuaternionTest::transformVectorNormalizedNotNormalized() {
     CORRADE_COMPARE(out.str(), "Math::Quaternion::transformVectorNormalized(): Quaternion({0.398736, 0, 0}, 1.95985) is not normalized\n");
 }
 
+void QuaternionTest::reflectVector() {
+    Vector3 normal = Vector3{-1.0f, 0.5f, -0.5f}.normalized();
+    Quaternion reflection = Quaternion::reflection(normal);
+    Matrix4 reflectionMatrix = Matrix4::reflection(normal);
+    Vector3 v{1.0f, 2.0f, 3.0f};
+
+    Vector3 reflected = reflection.reflectVector(v);
+    CORRADE_COMPARE(reflected, reflectionMatrix.transformVector(v));
+    CORRADE_COMPARE(reflected, (Vector3{-1.0f, 3.0f, 2.0f}));
+
+    /* Combining with rotations is ... involved */
+    Quaternion rotation = Quaternion::rotation(35.0_degf, Vector3{0.5f, 0.7f, 0.1f}.normalized());
+    Matrix4 rotationMatrix = Matrix4::rotation(35.0_degf, Vector3{0.5f, 0.7f, 0.1f}.normalized());
+    Vector3 transformed = (rotation*reflection*Quaternion{v}*reflection*rotation.conjugated()).vector();
+    CORRADE_COMPARE(transformed, rotation.transformVector(reflection.reflectVector(v)));
+    CORRADE_COMPARE(transformed, (rotationMatrix*reflectionMatrix).transformVector(v));
+    CORRADE_COMPARE(transformed, (Vector3{0.126405f, 2.03274f, 3.13879f}));
+}
+
 void QuaternionTest::strictWeakOrdering() {
     StrictWeakOrdering o;
     const Quaternion a{{1.0f, 2.0f, 3.0f}, 4.0f};