Apache Commons Math3 bug in Rotation / mat2quat

package org.example.bug.mat2quat;

import com.sun.javafx.geom.Matrix3f;
import com.sun.javafx.geom.Quat4f;
import org.apache.commons.lang.NotImplementedException;
import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;
import org.junit.Test;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import javax.vecmath.Matrix4d;
import javax.vecmath.Quat4d;
import java.util.Locale;

public class MatrixToQuaternionUnitTest
{
    private static final Logger log = LoggerFactory.getLogger(MatrixToQuaternionUnitTest.class);

    @Test
    public void commonsMathTest()
    {
        testCommonsMath(Math.PI / 2, 0);      // branch 1 in Rotation#mat2quat(final double[][])
        testCommonsMath(Math.PI / 10 * 9, 2); // branch 2
        testCommonsMath(Math.PI / 10 * 9, 1); // branch 3
        testCommonsMath(Math.PI / 10 * 9, 0); // branch 4

        /*
            Output:
                [ 0.7  0.0  0.0 -0.7] (Commons Math)
                [ 0.7  0.0  0.0  0.7] (JavaFX)
                [ 0.7  0.0  0.0  0.7] (JavaX Vecmath)
                [ 0.7  0.0  0.0  0.7] (NumPy + transformations.py)

                [-0.2  1.0  0.0  0.0] (Commons Math)
                [ 0.2  1.0  0.0  0.0] (JavaFX)
                [ 0.2  1.0  0.0  0.0] (JavaX Vecmath)
                [ 0.2  1.0  0.0  0.0] (NumPy + transformations.py)

                [ 0.2  0.0  1.0  0.0] (Commons Math)
                [ 0.2  0.0 -1.0  0.0] (JavaFX)
                [ 0.2  0.0 -1.0  0.0] (JavaX Vecmath)
                [-0.2  0.0  1.0  0.0] (NumPy + transformations.py)

                [-0.2  0.0  0.0  1.0] (Commons Math)
                [ 0.2  0.0  0.0  1.0] (JavaFX)
                [ 0.2  0.0  0.0  1.0] (JavaX Vecmath)
                [ 0.2  0.0  0.0  1.0] (NumPy + transformations.py)
         */
    }

    private static void testCommonsMath(double angle, int id)
    {
        Rotation rot = new Rotation(makeMatrix(angle, id).getData(), 1e-10);
        log.info(String.format((Locale)null, "[%4.1f %4.1f %4.1f %4.1f] (Commons Math)", rot.getQ0(), rot.getQ1(), rot.getQ2(), rot.getQ3()));

        Quat4f quatFx = new Quat4f();
        quatFx.set(makeJavaFxMatrix(angle, id));
        log.info(String.format((Locale)null, "[%4.1f %4.1f %4.1f %4.1f] (JavaFX)", quatFx.w, quatFx.x, quatFx.y, quatFx.z));

        Quat4d quatX = new Quat4d();
        makeJavaXMatrix(angle, id).get(quatX);
        log.info(String.format((Locale)null, "[%4.1f %4.1f %4.1f %4.1f] (JavaX Vecmath)", quatX.w, quatX.x, quatX.y, quatX.z));
    }

    private static RealMatrix makeMatrix(double angle, int id)
    {
        double[][] data;
        double c = Math.cos(angle);
        double s = Math.sin(angle);

        switch(id)
        {
            case 0: data = new double[][] { {c, -s, 0}, {s, c,  0}, {0, 0, 1} }; break;
            case 1: data = new double[][] { {c, 0, -s}, {0, 1,  0}, {s, 0, c} }; break;
            case 2: data = new double[][] { {1, 0,  0}, {0, c, -s}, {0, s, c} }; break;
            default: throw new NotImplementedException();
        }
        return new Array2DRowRealMatrix(data);
    }

    private static Matrix3f makeJavaFxMatrix(double angle, int id)
    {
        float[] data;
        float c = (float)Math.cos(angle);
        float s = (float)Math.sin(angle);

        switch(id)
        {
            case 0: data = new float[] {c, -s, 0,  s, c,  0,  0, 0, 1}; break;
            case 1: data = new float[] {c, 0, -s,  0, 1,  0,  s, 0, c}; break;
            case 2: data = new float[] {1, 0,  0,  0, c, -s,  0, s, c}; break;
            default: throw new NotImplementedException();
        }
        return new Matrix3f(data);
    }

    private static Matrix4d makeJavaXMatrix(double angle, int id)
    {
        double[] data;
        double c = Math.cos(angle);
        double s = Math.sin(angle);

        switch(id)
        {
            case 0: data = new double[] {c, -s, 0, 0, s,  c,  0, 0,  0, 0, 1, 0,  0, 0, 0, 1}; break;
            case 1: data = new double[] {c, 0, -s, 0, 0,  1,  0, 0,  s, 0, c, 0,  0, 0, 0, 1}; break;
            case 2: data = new double[] {1, 0,  0, 0, 0,  c, -s, 0,  0, s, c, 0,  0, 0, 0, 1}; break;
            default: throw new NotImplementedException();
        }
        return new Matrix4d(data);
    }
}