Skip to content

DualQuaternion implementation and YawPitchRoll and euler addition to Quaternion.#60

Closed
nanjizal wants to merge 13 commits intotbrosman:masterfrom
nanjizal:master
Closed

DualQuaternion implementation and YawPitchRoll and euler addition to Quaternion.#60
nanjizal wants to merge 13 commits intotbrosman:masterfrom
nanjizal:master

Conversation

@nanjizal
Copy link

@nanjizal nanjizal commented May 7, 2020

Hi Thomas

I have implemented a DualQuaternion based on my version but following where I can your Quaternion api design and your code style, I added the import to the test so it compiles without errors on neko.

My geom.matrix.DualQuaternion is used for Keyboard rotation translation here:
https://nanjizal.github.io/trilateral2Setup/binWebGL/index.html ( arrows keys, tab, shift, return, delete, space... ).

I believe the code should work as it's very similar to mine. But not yet implemented a visual test, or added any unit tests yet, I normally use DualQuaternion.create and pass the Quaternion real in with YawPitchRoll so I have added that to your Quaternion class, but you can use Axis approach like this because that uses your previous Quaternion code:

Basic use

var degrees = 45;
var axis = new Vector( 0,1,1);// x, y, z
var translation = new Vector4( 1, 1, 1, 1 ); // s, x, y, z 
var arm = hxmath.math.DualQuaternion.fromAxisAngle( 45, axis, translation );

( I have am not too sure if slerp or lerp will work, they need to be tried if they do work that might be fun.)

Generally you can build IK type structures quite easily with DualQuaternions just by multiplying them in any order to build rotation and translations like

dq0 * dq1 * dq2 == dq2 * dq1 * dq0 == dq1 * dq2 * dq0;

So you can get the position of a finger tip by setting up the positions:

var fingerTip = arm*forarm*wrist*hand*finger;

I am unsure on details of setting up constraints, it's quite hard to find lots of good details on DualQuaternion especially ones that don't get too technical, but they are used for robot arms.

To get the Quaternion rotation out you just call dq0.real, I have added a clone for that to increase clarity for users. To get the Vector4 out you can use getTranslation.

Added my name to haxelib authors file but feel free to remove, and I left a comment at top of the DualQuaternion class, but that can be removed ( especially when I get code role, not working in code at moment ).

Hope the code looks good it took awhile, so can't promise I will look at porting the Complex numbers implementation soon, although that is much better tested.

Best Justin

@tbrosman
Copy link
Owner

tbrosman commented May 7, 2020

I skimmed through it, looks like it covers the basics. I like it. I'll review more thoroughly when I get more free time. Can you add some tests around transformations? I have a few you can adapt: https://github.com/tbrosman/hxmath/blob/master/test/Test3D.hx It might also be good to create a few comparison tests vs. Frame3 to verify that it is equivalent.

@nanjizal
Copy link
Author

nanjizal commented May 7, 2020

I have added a test and some code to test eval and cppia target.
https://github.com/nanjizal/hxmath/blob/master/test/Test3D.hx#L133

Confused it seems to pass on some targets not others?!
Do you have any thoughts on how that could be possible?

+ neko pass
+ cppia pass
- cpp fail
- eval fail 
! hashlink untested 
Class: test.Test2D 
Class: test.Test3D 
Class: test.TestMathUtil 
Class: test.TestConverters 
Class: test.TestFrames 
Class: test.TestIntMath 
Class: test.TestGeom 
Class: test.TestDataStructures 
OK 0 tests, 0 failed, 0 success

Class: test.TestStructures ...............
Class: test.Test2D ..............
Class: test.Test3D .......F........
Class: test.TestMathUtil ..
Class: test.TestConverters ..
Class: test.TestFrames .......
Class: test.TestIntMath .
Class: test.TestGeom ..........
Class: test.TestDataStructures ...........
* test.Test3D::testMatrixFrameDualQuaternionInverse()
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 3.53537713894499861
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 3.53537713894499772
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 3.53537713894499861
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 3.53537713894499772
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 2.00670505870921811
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 2.00670505870922
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 2.00670505870921811
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 2.006705058709219
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 5.35721468440821891
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 5.35721468440822157
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 5.35721468440821713
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 5.35721468440821891
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 5.93403787126897
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 5.93403787126897342
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 5.93403787126897
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 5.93403787126897253
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 2.74635540911803977
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 2.74635540911803977
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 2.74635540911803977
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 2.74635540911803933
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 7.16634592299498507
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 7.16634592299498507
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 7.16634592299498507
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 7.16634592299498419
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 2.91834413116484104
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 2.91834413116484237
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 2.91834413116484193
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 2.91834413116484148
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 6.50769214412353
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 6.50769214412353136
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 6.50769214412353048
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 6.50769214412353136
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 3.01771591918229598
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 3.01771591918229554
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 3.01771591918229687
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 3.01771591918229642
./test/Test3D.hx:166: Test failed : expected 0 +-1e-06 but was 3.2330929345293784
./test/Test3D.hx:167: Test failed : expected 0 +-1e-06 but was 3.23309293452938
./test/Test3D.hx:168: Test failed : expected 0 +-1e-06 but was 3.23309293452938
./test/Test3D.hx:169: Test failed : expected 0 +-1e-06 but was 3.2330929345293784

@nanjizal
Copy link
Author

nanjizal commented May 7, 2020

Not sure invert is the ideal test... since that method is not currently in my version, so I am not totally sure on that specific method.
Perhaps an alternate test would be more suitable at the moment and come back to invert - just thinking...

I started trying to add a js version, not sure the tests are running correctly for neko and cppia, this is my first time with this test framework.

OK 0 tests, 0 failed, 0 success

@nanjizal
Copy link
Author

nanjizal commented May 7, 2020

I think it is definitely failing this test :( , but I would like to understand why neko/cppia do not seem to apply testing properly.

@nanjizal
Copy link
Author

nanjizal commented May 7, 2020

perhaps I have setup the matrix incorrectly since I use a 4x3 matrix with different labels, time to find some food. edit: Quick check modifying matrix still failed. need to dig deeper.

- Invert formula was wrong
- Multiplication re-ordered terms (quaternion multiplication is not commutative)
- addWith subtracted

Also added more tests.
@tbrosman
Copy link
Owner

tbrosman commented May 8, 2020

I ran Hashlink, cppia, and eval. I got the same failures on all three targets. I think your hxml targets might not be running stuff correctly.

I fixed the math issues and opened a PR into your fork: nanjizal#1

(Bonus: in the process I found a case where Frame3's internal Matrix can be modified directly. I'll open a bug for that.)

tbrosman and others added 2 commits May 7, 2020 23:00
@nanjizal
Copy link
Author

nanjizal commented May 9, 2020

Thomas

DON'T PULL THIS YET, I THINK YOUR MODIFICATION IS NOT GOING TO WORK

My DualQuaternion implementation is based on this paper, your changes don't seem to agree with the paper.
http://wscg.zcu.cz/wscg2012/short/A29-full.pdf

The multiplication within the paper is defined:

    rhs.m_real*lhs.m_real,
     rhs.m_dual*lhs.m_real + rhs.m_real*lhs.m_dual

Haxe implementation

@:op(A * B) public static inline
    function multiplyQ( q1: DualQuaternion, q2: DualQuaternion ):DualQuaternion {
        return new DualQuaternion({ real: q2.real * q1.real
                                  , dual: q2.dual * q1.real + q2.real * q1.dual });
    }

This is very important since with the amends you made, when applied to my library geom breaks or rather disables the 3D translation in my trilateral2 tests, ie applying any translation has no effect ( instead it seems to scale it in the z-axis on initial load, I checked my original code against the paper and there was an error last two terms swapped over but did not seem to matter to translation. I have updated geom with the function above and not changed getTranslation code yet.

[ I am not sure if it's relevant but my geom library uses real matrix positions in normal maths, and transposes the transform for the WebGL shaders only at the end just before sending to shader, are you doing it this way, lots of online stuff is a little confusing in this regard, are you doing it this way? ]

To apply keyboard translations and rotations my code currently uses Axis3 effectively it is doing
dualQ_a * dualQ_b when ever the delta rotation/translation changes ( in theory dualQ_b * dualQ_a should be the same, I might be doing the former ). The code uses a three state Bool to apply the key changes so you can do orthogonal changes at the same time. I think it probably gets round lerp/slerp not being always linear so ideal for simple cube test since the delta is always the same.
https://github.com/nanjizal/geom/blob/master/src/geom/move/Axis3.hx#L49

I can have a go at setting up a WebGL demo with hxmath ( safer for a test than using OpenFL because you can be more sure of internal details from start to finish ) to try to double check physical results. But before doing so, as it can often be fiddly, I thought I should check in and get your thoughts, probably I should have mentioned the paper earlier, I think your likely much better at the maths and will appreciate the paper.

Prior to consideration on inverse it's important to resolve multiplication definition, I feel the one above is suitable - agreeing with the theoretical paper and seems to work practically with my WebGL experiments?

Best

Justin

@tbrosman
Copy link
Owner

tbrosman commented May 9, 2020

I'm not sure why their code swaps the order of multiplication. Section 5.2 of the same paper actually defines multiplication correctly. Since dual quaternion multiplication is non-commutative (like regular quaternions), the result would be different.

(p + e q) (r + e s) = p r + e (p s + q r)

Quaternion multiplication for p = [s, u] and q = [t, v] can be written using mixed vectors and scalars:

[s, u][t, v] = [s t - u * v, s v + t u + u X v]

(where * is a dot product and X is a cross product). Since u X v = - u X v, the quaternion product does not commute.

It might be that their multiplication is like "concat" and not really multiply, and then everything was multiplied in reverse. Kind of weird. This is an even further reach, but it might have been done intentionally to match up with matrix multiplication order.

@nanjizal
Copy link
Author

nanjizal commented May 9, 2020

This paper seems to have:
https://pdfs.semanticscholar.org/7bf0/ff0c2a6161f25f1f9d669da65ee896a8e99c.pdf?_ga=2.104136838.74812921.1589012781-191394131.1589012781

// Multiplication order - left to right
public static DualQuaternion_c operator *(DualQuaternion_c lhs, DualQuaternion_c rhs)
{
 lhs = DualQuaternion_c.Normalize( lhs );
 rhs = DualQuaternion_c.Normalize( rhs );
 return new DualQuaternion_c( rhs.m_real * lhs.m_real,
 rhs.m_dual * lhs.m_real + rhs.m_real * lhs.m_dual);
}

The original paper is the paper is the one the Luxe engine used.
https://github.com/underscorediscovery/luxe/blob/master/phoenix/DualQuaternion.hx

@nanjizal
Copy link
Author

nanjizal commented May 9, 2020

@tbrosman
Copy link
Owner

tbrosman commented May 9, 2020

I created a test comparing my Frame implementation to DualQuaternion and a short constructive walkthrough: https://gist.github.com/tbrosman/43c451f10cddd0826353ccf8af716898 They are equivalent in the current implementation. Test:

public function testDualQuaternionMult_ConsistentWithFrame3()
{
    var dualQuatA = DualQuaternion.fromAxisAngle(90, Vector3.yAxis, new Vector4(10, 0, 0, 1));
    var dualQuatB = DualQuaternion.fromAxisAngle(90, Vector3.xAxis, new Vector4(0, 0, 5, 1));
    var dualQuatC = dualQuatA * dualQuatB;

    var frameA = new Frame3(new Vector3(10, 0, 0), Quaternion.fromAxisAngle(90, Vector3.yAxis));
    var frameB = new Frame3(new Vector3(0, 0, 5), Quaternion.fromAxisAngle(90, Vector3.xAxis));
    var frameC = frameA.concat(frameB);

    var dualQuatCRotation = dualQuatC.getRotationQuat();

    // TODO: Make getTranslation return a Vector3
    var dualQuatCTranslation4 = dualQuatC.getTranslation();
    var dualQuatCTranslation = new Vector3(dualQuatCTranslation4.x, dualQuatCTranslation4.y, dualQuatCTranslation4.z);

    var frameCRotation = frameC.orientation;
    var frameCTranslation = frameC.offset;

    assertApproxEquals(0.0, (dualQuatCRotation - frameCRotation).length);
    assertApproxEquals(0.0, (dualQuatCTranslation - frameCTranslation).length);
}

The flipped multiplication just seems weird. An example of why it is weird: say you have two dual quaternions with only a rotation (labeled "real" currently) portion. What happens when you concatenate them? In the code you provided:

public static DualQuaternion_c operator *(DualQuaternion_c lhs, DualQuaternion_c rhs)
{
 lhs = DualQuaternion_c.Normalize( lhs );
 rhs = DualQuaternion_c.Normalize( rhs );
 return new DualQuaternion_c( rhs.m_real * lhs.m_real,
 rhs.m_dual * lhs.m_real + rhs.m_real * lhs.m_dual);
}

This is equivalent to some operator (X) that acts like this:

(p + e q) (X) (r + e s) = r*p + e(s*p + r*q)

Which means in a pure rotation case (q = 0, s = 0) you would have r * p. Since these are quaternions, r * p != p * r and (X) would give results that were inconsistent with just multiplying two regular quaternions (even though that's what the inputs were).

@nanjizal
Copy link
Author

Thomas
I thought I would enquire on 'bivector' discord. I am wondering if both should be supported, but not yet at a stage to know, I found some pdf's but not clear if any are the correct one referenced below, it's quite late to try reading them, I have yet to see how easy it is to use your Biquaternion/DualQuaternion multiplication in my code, correctness is less important to me than ease of use I guess. Sorry it seems there is more than one acceptable implementation, probably need to read more before understanding the reasons.

I'm probably out of my depth I have posted a link to multiplication in geom and in hxmath perhaps they will comment further. Here is snipit from conversation:

Nanjizal>
"Hi I have a question about DualQuaternions perhaps you can advise.
The code in this paper
https://cs.gmu.edu/~jmlien/teaching/cs451/uploads/Main/dual-quaternion.pdf
has multipication
real = rhs.m_reallhs.m_real
dual = rhs.m_duallhs.m_real + rhs.m_real*lhs.m_dual
which seems to be the wrong way round? With my webgl test it seems to work
https://nanjizal.github.io/dice/binWebGL/ but maybe it's down to how I am using it.
Also is dq0 * dq1 == dq1 * dq0
Have been making a pull to another maths library and needed more clarity on these."

chakravala>
'In my opinion, the "dual" basis element which squares to zero must be commutative, it's actually a different type of basis element that doesnt anti-commute, but a lot of other people insist on using a degenerate metric with a non-commuting basis (which i consider incorrect for the dual number basis). I assume you are asking about commutativity, but it's not clear'

mewertx0rz>
'@nanjizal it's worth reading Cliffords original writing on dual quaternions. "Preliminary Sketch of Biquaternions" Little tough to read ( 1870's language ) but it has some interesting insights. For example the non-degenerate dual is in there. There's more than one correct answer to how to implement dual quaternions. If yours is working, great.'

@tbrosman
Copy link
Owner

I think I can put some geometric intuition behind the non-commutativity. A dual quaternion can be thought of as a map that takes some point v and then applies (in order) rotation and then translation. In general, rigid body transformations do not commute. If you rotate a point 90 degrees counter-clockwise around the X axis and then add (0, 10, 0) to it the result is not the same as adding (10, 0, 0) first and then rotating.

Just rotation:

R(90) * (x, y, z) = (x, -z, y)

Just translation:

T(10, 0, 0) * (x, y, z) = (x, y + 10, z)

Translation, then rotation:

R(90) * T(10, 0, 0) * P = R(90) (x, y + 10, z) = (x, -z, y + 10)

Rotation, then translation:

T(10, 0, 0) * R(90) * P = T(10, 0, 0) (x, -z, y) = (x + 10, -z, y)

@nanjizal
Copy link
Author

Ok cool, I did a lot more research your approach seems the best option, sorry for my caution.

I have not yet created a visual test, I have been fighting with general setting up WebGL shaders in OpenFL, but if your happy to accept pull then go for it. I will look to add the 'complex number class' soon. It is likely a lot more robust as it has quite a few tests.

*
* @return The modified object.
*/
public inline function applyConjugate():DualQuaternion
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be good to skip this operation. There are two conjugates for DualQuaternion: the dual conjugate and the complex conjugate.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

adjusted name

* @return ~a
*/
@:op(~A)
public static inline function conjugate(a: DualQuaternion):DualQuaternion
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See comment below. Should either be skipped, or the documentation should make it clear this is the complex conjugate. I could also see this being implemented as the complex + dual conjugate. This is the way transformation by a dual quaternion is implemented: q * v * complexConjugate(dualConjugate(q)).

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you have lost me here.

* @param z
* @return self
*/
public inline function setReal(s:Float, x:Float, y:Float, z:Float):DualQuaternion
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure these helpers are useful. The fields are forwarded, so probably easier to set it that way.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed

* but set /gets internally as the Quaternion value
*
**/
public var euler(get, set): Quaternion;
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wouldn't expect this as a Quaternion. Either needs its own structure, or you could use Vector3.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added Euler abstract over Vector3

test_eval.hxml Outdated
-lib nanotest
-dce full
-cmd mkdir -p bin/eval
#-cmd mv test.cppia bin/cppia/test.cppia
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor: cleanup

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed commented line

*
* @return The rotational Quaternion
*/
public inline function getRotationQuat(): Quaternion
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about just "getRotation"?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

amended

@tbrosman
Copy link
Owner

tbrosman commented May 15, 2020

Cool, I'll try to add more documentation/references when I get time to work on this again. This is one of the more useful papers I found while refreshing my brain: https://www.cs.utah.edu/~ladislav/kavan06dual/kavan06dual.pdf It goes over how to set up the product for transforming points.

In the meantime there are a few things this PR needs:

  • Code style should match the rest of the library
  • Test for Quaternion.fromYawPitchRoll
  • Test for DualQuaternion.fromAxisAngle, fromTranslation (should match result of getTranslation())
  • Test for DualQuaternion.identity

One more thing I keep forgetting: if I include your name in contributors, haxelib associates the library with your username. It is easier for me to just credit you in the CHANGELOG.md (and class file). I found this out the hard way a year ago when someone updated my Haxe REST Client library and I tried to submit the change. Basically haxelib doesn't differentiate contributors and owners: https://lib.haxe.org/documentation/faq/

nanjizal added 2 commits May 15, 2020 08:56
Copy link
Author

@nanjizal nanjizal left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have adjusted the Quaternion stuff, it should be tested but internally the code agrees with Wikipedia.
I have left the new 'applyDualConjugate' method with a throw not implemented please advise.
Testing not yet addressed.

@nanjizal
Copy link
Author

I will look at testing hopefully later today.

@nanjizal
Copy link
Author

I have not added the Tests you asked for yet, sorry but I have instead implemented ( ported ) a Complex class. The Complex has coverage for unit tests ( for most aspects ).
I have not yet added Docs for all the methods, it already took lot of time to get formatting like yours.
Let me know what you think, I can flesh out the method docs.
I expect you will not want as many standard complex numbers.

Best Justin

@tbrosman
Copy link
Owner

Can you make the Complex class a separate PR? I think it's out of the scope of this change. Let me know if you need any suggestions/ideas for the tests.

@nanjizal
Copy link
Author

I have moved the code into two separate branches. Unfortunately it seems simplest to close this pull the Dual code is now in this pull.
#64

@nanjizal nanjizal closed this May 19, 2020
@nanjizal
Copy link
Author

Note that so far have only separated them. They are both WIP.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants