Author Topic: Euler angles from transformation matrix?  (Read 2423 times)

0 Members and 1 Guest are viewing this topic.

kaefer

• Swamp Rat
• Posts: 572
Euler angles from transformation matrix?
« on: September 02, 2011, 05:18:00 am »
There's been a contribution to the Autodesk AutoCAD .NET Discussion Group which is based on an Inventor matrix, using the arccos function. I'm suspecting that Inventor, missing the translational part of the matrix, keeps additional constraints that make such calculations possible. Whatever the reason, I wasn't able to port that approach to AutoCAD.

Now here's a naive alternative, which works at least for cartesian rotation (along one axis). Has anybody a better idea?

Code: [Select]
`        let ent = tr.GetObject(per.ObjectId, OpenMode.ForRead) :?> Entity        let cs = ent.CompoundObjectTransform.CoordinateSystem3d        let a = cs.Yaxis.GetAngleTo(Vector3d.ZAxis, Vector3d.XAxis) - pi2        let b = -cs.Zaxis.GetAngleTo(Vector3d.XAxis, Vector3d.YAxis) + pi2        let c = -cs.Xaxis.GetAngleTo(Vector3d.YAxis, Vector3d.ZAxis) + pi2        ed.WriteMessage("\na: {0}, b: {1}, c: {2} ", a, b, c )`

kaefer

• Swamp Rat
• Posts: 572
Re: Euler angles from transformation matrix?
« Reply #1 on: September 02, 2011, 06:01:09 pm »
I'm suspecting that Inventor, missing the translational part of the matrix, keeps additional constraints

More trivially, Inventor.Matrix indexer is 1-based. Besides, I've been advised to apply the fishmonger's treatment: Remove scales.

Here's a loose translation of the original VB, with test command, in C#.

Code: [Select]
`        [CommandMethod("CalculateRotationAngles")]        public void CalculateRotationAngleCommand()         {            Document doc = AcadApp.DocumentManager.MdiActiveDocument;            Editor ed = doc.Editor;            PromptEntityOptions peo = new PromptEntityOptions("Select BlockReference");            peo.SetRejectMessage("BlockReference only");            peo.AddAllowedClass(typeof(BlockReference), false);            PromptEntityResult per = ed.GetEntity(peo);            if(per.Status != PromptStatus.OK) return;            using(Transaction tr = doc.Database.TransactionManager.StartTransaction())            {                BlockReference bref = (BlockReference)tr.GetObject(per.ObjectId, OpenMode.ForRead);                double[] rot = CalculateRotationAngles(bref.BlockTransform);                ed.WriteMessage("\na: {0}, b: {1}, c: {2} ", rot[0], rot[1], rot[2]);                tr.Commit();            }        }        double Acos(double value)        {            return value >= 1 ? 0 : value <= -1 ? Math.PI : Math.Acos(value);        }        double[] CalculateRotationAngles(Matrix3d oMatrix)        {            oMatrix = Scale3d.RemoveScale(oMatrix).Inverse().GetMatrix() * oMatrix;            // Original VB for Inventor, posted by philippe.leefsma; found here:            // http://forums.autodesk.com/t5/NET/Getting-blocs-amp-solids-rotation-angle-in-3D-space/td-p/3140470            double[] aRotAngles = new double[3];            // Choose aRotAngles[0] about x which transforms axes[2] onto the x-z plane            Vector2d v = new Vector2d(oMatrix[1, 2], oMatrix[2, 2]);            aRotAngles[0] = v.Length <= 0.000001 ? 0 : Math.Sign(v.X) * Acos(v.Y / v.Length);            oMatrix = Matrix3d.Rotation(aRotAngles[0], Vector3d.XAxis, Point3d.Origin) * oMatrix;             // Choose aRotAngles[1] about y which transforms axes[3] onto the z axis            aRotAngles[1] = Math.Sign(-oMatrix[0, 2]) * Acos(oMatrix[2, 2]);            oMatrix = Matrix3d.Rotation(aRotAngles[1], Vector3d.YAxis, Point3d.Origin) * oMatrix;                // Choose aRotAngles[2] about z which transforms axes[0] onto the x axis            aRotAngles[2] = Math.Sign(-oMatrix[1, 0]) * Acos(oMatrix[0, 0]);             // if you want to get the result in degrees            const double toDegrees = 180 / Math.PI;            aRotAngles[0] = aRotAngles[0] * toDegrees;            aRotAngles[1] = aRotAngles[1] * toDegrees;            aRotAngles[2] = aRotAngles[2] * toDegrees;                        return aRotAngles;        }`

LE3

• Guest
Re: Euler angles from transformation matrix?
« Reply #2 on: September 02, 2011, 08:19:06 pm »
some links that might help:

http://geometrictools.com/
&
http://cmldev.net/ [The CML (Configurable Math Library) is a free C++ math library for games and graphics.]
&
http://tog.acm.org/resources/GraphicsGems/gems.html
&
http://www.flipcode.com/documents/matrfaq.html
« Last Edit: September 02, 2011, 10:53:56 pm by DeVo »

kaefer

• Swamp Rat
• Posts: 572
Re: Euler angles from transformation matrix?
« Reply #3 on: September 05, 2011, 11:48:38 am »
http://www.flipcode.com/documents/matrfaq.html

Many thanks! Those links sure do help:
Code: [Select]
`    // Implementation lifted from    // http://www.flipcode.com/documents/matrfaq.html#Q37    let getAngles (mat: Matrix3d) =        let mat =  (Scale3d.RemoveScale mat).Inverse().GetMatrix() * mat        let ay = acos(mat.ElementAt(0, 2)) - System.Math.PI / 2.        let c = cos ay        let (ax, az) =            if abs c > 0.005 then                let xtrx = mat.ElementAt(2, 2) / c                let xtry = -mat.ElementAt(1, 2) / c                let ztrx = mat.ElementAt(0, 0) / c                let ztry = -mat.ElementAt(0, 1) / c                atan2 xtry xtrx, atan2 ztry ztrx            else                let ztrx = mat.ElementAt(1, 1)                let ztry = -mat.ElementAt(1, 0)                0., atan2 ztry ztrx        // if you want to get the result in degrees        let toDegrees = 180. / System.Math.PI        ax * toDegrees,        ay * toDegrees,        az * toDegrees`

LE3

• Guest
Re: Euler angles from transformation matrix?
« Reply #4 on: September 05, 2011, 11:55:12 am »
^
Great.
« Last Edit: September 05, 2011, 12:16:57 pm by DeVo »

gile

• Water Moccasin
• Posts: 2094
• Marseille, France
Re: Euler angles from transformation matrix?
« Reply #5 on: December 18, 2016, 10:48:34 am »
Hi,

I know this is an old post but a similar question in the Autodesk Discussion group and I wanted to share here my reply.
It implements a base 'EulerAngles' for common purposes and two classes for proper Euler angles and Tait-Bryan angles which allows conversions from Matrix3d to angles and vice-versa.

A base 'EulerAngles' abstract class
Code - C#: [Select]
`using Autodesk.AutoCAD.Geometry;using static System.Math; namespace Gile.AutoCAD.Geometry{    /// <summary>    /// Base class for Euler angles defintions.    /// </summary>    public abstract class EulerAngles    {        protected double psi, theta, phi;        protected Matrix3d xform;         /// <summary>        /// Gets the angle alpha (or psi).        /// </summary>        public double Alpha => psi;         /// <summary>        /// Gets the angle beta (or theta).        /// </summary>        public double Beta => theta;         /// <summary>        /// Gets the angle gamma (or phi).        /// </summary>        public double Gamma => phi;         /// <summary>        /// Get the transformation matrix.        /// </summary>        public Matrix3d Transform => xform;         /// <summary>        /// Base constructor.        /// </summary>        /// <param name="transform">Transformation matrix.</param>        public EulerAngles(Matrix3d transform)        {            if (!transform.IsUniscaledOrtho())                throw new System.ArgumentException("Non uniscaled ortho matrix.");            xform = transform;        }         /// <summary>        /// Base constructor.        /// </summary>        /// <param name="alpha">Precession angle.</param>        /// <param name="beta">Nutation angle.</param>        /// <param name="gamma">Intrinsic rotation.</param>        public EulerAngles(double alpha, double beta, double gamma)        {            psi = Wrap(alpha);            theta = Wrap(beta);            phi = Wrap(gamma);        }         /// <summary>        /// Equality operator.        /// </summary>        /// <param name="x">Left operand.</param>        /// <param name="y">Right operand.</param>        /// <returns>true if the values of its operands are equal, false otherwise.</returns>        public static bool operator ==(EulerAngles x, EulerAngles y) => x.Equals(y);         /// <summary>        /// Inequality operator.        /// </summary>        /// <param name="x">Left operand.</param>        /// <param name="y">Right operand.</param>        /// <returns>true if the values of its operands are not equal, false otherwise.</returns>        public static bool operator !=(EulerAngles x, EulerAngles y) => !x.Equals(y);         /// <summary>        /// Determines whether the specified object is equal to the current object.        /// </summary>        /// <param name="obj">The object to compare with the current object.</param>        /// <returns>true if the specified object is equal to the current object; otherwise, false.</returns>        public override bool Equals(object obj)        {            var other = obj as EulerAngles;            return psi == other?.psi && theta == other.theta && phi == other.phi;        }         /// <summary>        /// Serves as hash function.        /// </summary>        /// <returns>A hash code for the current object.</returns>        public override int GetHashCode() => psi.GetHashCode() ^ theta.GetHashCode() ^ phi.GetHashCode();         /// <summary>        /// Returns a string that represents the current object.        /// </summary>        /// <returns>A string that represents the current object.</returns>        public override string ToString() => \$"({Alpha}, {Beta}, {Gamma})";         /// <summary>        /// Return the angle in [0,2*PI) range.        /// </summary>        private double Wrap(double angle) =>             0.0 <= angle && angle < 2 * PI ? angle : Asin(Sin(angle));    }}`

The 'ProperEuler' class (can also be used with Normal and Rotation)
Transform = Rotation_Z(Alpha) * Rotation_X(Beta) * Rotation_Z(Gamma)
Code - C#: [Select]
`using Autodesk.AutoCAD.Geometry;using static System.Math; namespace Gile.AutoCAD.Geometry{    /// <summary>    /// Defines the relations between a transformation matrix and Euler angles    /// (proper Euler angles using z-x'-z" convention).    /// </summary>    public class ProperEuler : EulerAngles    {        Matrix3d planeToWorld;         /// <summary>        /// Gets the normal of the plane.        /// </summary>        public Vector3d Normal => planeToWorld.CoordinateSystem3d.Zaxis;         /// <summary>        /// Gets the rotation on the plane;        /// </summary>        public double Rotation => phi;         /// <summary>        /// Create a new intance of ProperEuler.        /// </summary>        /// <param name="transform">Transformation matrix.</param>        public ProperEuler(Matrix3d transform) : base(transform)        {            theta = Acos(xform[2, 2] / xform.GetScale());            if (Abs(theta) < 1e-7)            {                theta = 0.0;                psi = Atan2(xform[1, 0], xform[1, 1]);                phi = 0.0;            }            else            {                psi = Atan2(xform[0, 2], -xform[1, 2]);                phi = Atan2(xform[2, 0], xform[2, 1]);            }            planeToWorld =                 Matrix3d.Rotation(psi, Vector3d.ZAxis, Point3d.Origin) *                Matrix3d.Rotation(theta, Vector3d.XAxis, Point3d.Origin);        }         /// <summary>        /// Create a new intance of ProperEuler.        /// </summary>        /// <param name="alpha">Rotation angle around the Z axis.</param>        /// <param name="beta">Rotation angle around the X' axis.</param>        /// <param name="gamma">Rotation angle around the Z" axis.</param>        public ProperEuler(double alpha, double beta, double gamma) : base(alpha, beta, gamma)        {            planeToWorld =                Matrix3d.Rotation(alpha, Vector3d.ZAxis, Point3d.Origin) *                Matrix3d.Rotation(beta, Vector3d.XAxis, Point3d.Origin);            xform = planeToWorld *                 Matrix3d.Rotation(gamma, Vector3d.ZAxis, Point3d.Origin);        }         /// <summary>        /// Create a new intance of ProperEuler.        /// </summary>        /// <param name="normal">Plane normal.</param>        /// <param name="rotation">Proper rotation.</param>        public ProperEuler(Vector3d normal, double rotation)            : this(Matrix3d.PlaneToWorld(normal) *                   Matrix3d.Rotation(rotation, Vector3d.ZAxis, Point3d.Origin))        { }     /// <summary>    /// Determines whether the specified object is equal to the current object.    /// </summary>    /// <param name="obj">The object to compare with the current object.</param>    /// <returns>true if the specified object is equal to the current object; otherwise, false.</returns>    public override bool Equals(object obj) => obj is ProperEuler && base.Equals(obj);         /// <summary>        /// Serves as hash function.        /// </summary>        /// <returns>A hash code for the current object.</returns>        public override int GetHashCode() => base.GetHashCode();    }}`

The TaitBryan class (yaw,pitch, roll)
Transform = Rotation_Z(Alpha) * Rotation_Y(Beta) * Rotation_X(Gamma)
Code - C#: [Select]
`using Autodesk.AutoCAD.Geometry;using static System.Math; namespace Gile.AutoCAD.Geometry{    /// <summary>    /// Defines the relations between a transformation matrix and Euler angles    /// (specific Euler angles called Tait-Ryan angles using z-y'-x" convention).    /// </summary>    public class TaitBryan : EulerAngles    {        /// <summary>        /// Create a new intance of TaitBryan.        /// </summary>        /// <param name="transform">Transformation matrix.</param>        public TaitBryan(Matrix3d transform) : base(transform)        {            theta = -Asin(xform[2, 0] / xform.GetScale());            if (Abs(theta - PI * 0.5) < 1e-7)            {                theta = PI * 0.5;                psi = Atan2(xform[1, 2], xform[1, 1]);                phi = 0.0;            }            else if (Abs(Beta + PI * 0.5) < 1e-7)            {                theta = -PI * 0.5;                psi = Atan2(-xform[1, 2], xform[1, 1]);                phi = 0.0;            }            else            {                psi = Atan2(xform[1, 0], xform[0, 0]);                phi = Atan2(xform[2, 1], xform[2, 2]);            }        }         /// <summary>        /// Create a new intance of TaitBryan.        /// </summary>        /// <param name="alpha">Rotation angle around the Z axis.</param>        /// <param name="beta">Rotation angle around the Y' axis.</param>        /// <param name="gamma">Rotation angle around the X" axis.</param>        public TaitBryan(double alpha, double beta, double gamma) : base(alpha, beta, gamma)        {            xform =                 Matrix3d.Rotation(alpha, Vector3d.ZAxis, Point3d.Origin) *                Matrix3d.Rotation(beta, Vector3d.YAxis, Point3d.Origin) *                Matrix3d.Rotation(gamma, Vector3d.XAxis, Point3d.Origin);        }         /// <summary>        /// Determines whether the specified object is equal to the current object.        /// </summary>        /// <param name="obj">The object to compare with the current object.</param>        /// <returns>true if the specified object is equal to the current object; otherwise, false.</returns>        public override bool Equals(object obj) => obj is TaitBryan && base.Equals(obj);         /// <summary>        /// Serves as hash function.        /// </summary>        /// <returns>A hash code for the current object.</returns>        public override int GetHashCode() => base.GetHashCode();    }}`
« Last Edit: December 18, 2016, 06:43:32 pm by gile »
Speaking English as a French Frog