search

jdelauney / SIMD-VectorMath-UnitTest

For testing asm SIMD (SSE/SSE 2/SSE 3/SSE 4.x / AVX /AVX 2) vector math library (2f, 4f, matrix, quaternion...) with Lazarus and FreePascal Compiler
Mozilla Public License 2.0
8 stars 0 forks source link

Quaternions #27

Open dicepd opened 6 years ago

dicepd commented 6 years ago

What fun this has been so far.

Just to let you know the changes.

Made mul conform to first second as there is a MulFirst (not tested yet) which does non Math order. Commented out large parts of CreateEuler, trying to hit the exceptions that they are guarding against, but not managed yet. Rotations conform to existing rotations. Euler create gives correct numbers.

Default create rpy not done as we will probably have to decide on world up for the system as this is meant to be the quick quat that defaults to world.

dicepd commented 6 years ago

I found this about quats https://blog.molecular-matters.com/2013/05/24/a-faster-quaternion-vector-multiplication/ I amy code this up and test when I have finished the quat tests and added to from Matrix

jdelauney commented 6 years ago

It have also some interesting articles. It will be no difficult to implement

dicepd commented 6 years ago

Right when I get to Slerp ordinary Slerp works fine. How ever Slerp with spin is fundamentally flawed from the first line of code onward. It may give increasing angles but not in a defined manner.

cost:=VectorAngleCosine(QStart.ImagPart, QEnd.ImagPart);  <--- Imag part of quat holds half angles

Again no one has used this in Scene, another for the bin? I will check in as is, so you can play with the behaviour to see if I am 'doing it wrong'

dicepd commented 6 years ago

Ok maybe you can scale with quat, something I had not seen before

You don’t seem to find it in writing anywhere, but quaternions support uniform scaling just fine: p’ = q p q (simplified rotation form)…multiply q by non-zero scalar ‘s’: (qs) p (qs) = ss (q p q*) = ss p’. So sqrt(s) is a uniform scale factor and these obviously naturally compose.

jdelauney commented 6 years ago

Many things don't work here with Quaternion :

What the change you did exactly ? Normalize and mul worked before

dicepd commented 6 years ago

Are you normalising 3 element or 4 with Quats you have to normalise all four not just the inaginary vector

"CreateTwoUnitAffine:Sub3 Z failed " expected: <0,7071> but was: <1>

this result would point to that, Quats are 4 dimensional beasts so are normalised as such. It is not a vector + flag or Vector + length. So if you copied Vector normalise then this would be the cause.

jdelauney commented 6 years ago

I' m doing nothing just ran the test

jdelauney commented 6 years ago

I found article in french http://mecaspa.cannes-aero-patrimoine.net/SCAO/QUATERN/complements/som_quat.htm#note_calcul_perso but don't understand well all yet. I need little more time

dicepd commented 6 years ago

Ok I just looked at the win64 SSE and the above is exactly the reason for most of the errors. only error in win64 native is slerp spin

jdelauney commented 6 years ago

by comment andps xmm2, [RIP+cSSE_MASK_NO_W] in normailize, now the SlerpSingle func is working but not others

dicepd commented 6 years ago

Now onto the gimbal lock code in Euler Create. What I think has happened here is the original code in the sub function was the first written, someone read about gimbal lock and decided to 'protect' original code with some code found somewhere, not realising that the original code was gimbal lock safe as it creates 3 quats and returns the product of the three, This is inherently safe from gimbal lock as it is combining quats not euler angles.

dicepd commented 6 years ago

maybe this

Made mul conform to first second as there is a MulAsSecond (not tested yet) which does non Math order.

This cured a lot of ills in the pascal. So tests reflect what it naturally broke. Try swapping A and B registers

jdelauney commented 6 years ago

Try swapping A and B registers

In which function ???? i'm lost, normailze ? or Mul

dicepd commented 6 years ago

I have checked in the fixes

dicepd commented 6 years ago

Normalize could do with reworking and simplifying but it does work

jdelauney commented 6 years ago

Ok, i'll compare with old code, i don't see what you're doing. So now only 1 fail :

dicepd commented 6 years ago

That one is to be expected that is the routine that makes no sense to me.

dicepd commented 6 years ago

Thinking ahead a bit while I finish Quat and Matrix what functionallity are we missing from fustrum. That I think, is the last piece of the jigsaw to a simple software world line render. We have plane and the ability to test for in plane. Do we need a create from Quat/EyeVector or similar.

jdelauney commented 6 years ago

Yes, you're right missing some functions, but not many. For a 3D word render. After some learning somethings about quaternion. I'm thinking it will good if we'll add RotateWithQuaternion functions in TGLZVectorHelper. From what I read, we would have several advantages to using quaternions in rotations rather than matrices

dicepd commented 6 years ago

The one issue I had withquat was needing to scale. Take boidz as an example with a simple 3D mesh . Place master boidz at center fustrum at suitable scale. For each boid work out where center of boid[i] is in relation to master boid locally orient boid work out scale near to far add Boid orient + boid scale + boid pos to a matrix and apply

Without quats that is a lot of Mat manipulation for each item with quat we have delta orient, and if scale works then a quat move and scale from two vectors center of viewplane (fixed for all) and point project center boid to viewplane, so two quat create and two quat mul array. No idea if this is more efficient for low poly but may be, for high poly maybe not. But then it maybe more efficent to create the two quats convert to Mat and Add two mat.

jdelauney commented 6 years ago

I found this for unity (it's in french) https://www.youtube.com/watch?v=EWQCsjp3NvY it describe very well how Unity work with quaternion and why. I understand much better quaternion with this and making the link with maths formulas. In terms of performance computing quaternion is less expensive than matrices some infos (near twice less operations) I'm also falling on that and this. in an article. What do you think about quaternion structure and functions ?

jdelauney commented 6 years ago

I'm trying to fix slerpSpin

Like we can see the result is in W

jdelauney commented 6 years ago

Ok by accessing var in with aq4.Realpart instead aq4.Z get :

jdelauney commented 6 years ago

After some mino changes, now i've :

Have you swap components results in functionnal test ? or it is something wrong

jdelauney commented 6 years ago

After some minors change, the problems come always from the Z

Now the result don't match. It's a pain

jdelauney commented 6 years ago

if i normailze the result i get :

dicepd commented 6 years ago

sorry just back from shopping, tests may may not be correct as I had no idea if spin was 180 or 360 thats what I was trying to determine when nothing worked in a sane manner

jdelauney commented 6 years ago

now if i negate result.z : ( and i change slerpsingle to SlerpSpin, for not to be confused)

-"SlerpSpin:Sub7 Z failed (ImagePart.X: 0,00000 ,ImagePart.Y: 0,00000 ,ImagePart.Z: -0,17365 , RealPart.W: -0,98481)" expected: <0,766> but was: <-0,1736>

I don't understand why result is false at this stage

this is the change i've made for testing : 

function TGLZQuaternion.Slerp(const QEnd: TGLZQuaternion; Spin: Integer; t: Single): TGLZQuaternion;
var
   to1: array[0..4] of Single;
   phi,omega, cosom, sinom, scale0, scale1: Extended;
// t goes from 0 to 1
// absolute rotations
begin
   // calc cosine
   cosom:= Self.ImagePart.X*QEnd.ImagePart.X
          +Self.ImagePart.Y*QEnd.ImagePart.Y
          +Self.ImagePart.Z*QEnd.ImagePart.Z
      +Self.RealPart*QEnd.RealPart;
   // adjust signs (if necessary)
   if cosom<0 then
   begin
      cosom := -cosom;
      to1[0] := - QEnd.ImagePart.X;
      to1[1] := - QEnd.ImagePart.Y;
      to1[2] := - QEnd.ImagePart.Z;
      to1[3] := - QEnd.RealPart;
   end
   else
   begin
      to1[0] := QEnd.ImagePart.X;
      to1[1] := QEnd.ImagePart.Y;
      to1[2] := QEnd.ImagePart.Z;
      to1[3] := QEnd.RealPart;
   end;
   // calculate coefficients
   if ((1.0-cosom)>cEpsilon) then
   begin
      //Cosom := Cosom
      omega:=GLZMath.ArcCos(cosom);
      Omega :=(Omega + DegToRadian(Spin*360));
      sinom:=1/Sin(omega);
      phi:=omega;// + Spin * cPi;
      scale0:=sin(omega - t * phi) * sinom;
      scale1:=sin(t * phi) * sinom;
      //scale0:=Sin((1.0-t)*omega)*sinom;
      //scale1:=Sin(t*omega)*sinom;
   end
   else  // "from" and "to" quaternions are very close
   begin
      //  ... so we can do a linear interpolation
      scale0:=1.0-t;
      scale1:=t;
   end;
   // calculate final values
   Result.ImagePart.V[0] := scale0 * Self.ImagePart.V[0] + scale1 * to1[0];
   Result.ImagePart.V[1] := scale0 * Self.ImagePart.V[1] + scale1 * to1[1];
   Result.ImagePart.V[2] := -(scale0 * Self.ImagePart.V[2] + scale1 * to1[2]); // Don't understand why we need negate here
   Result.RealPart := (scale0 * Self.RealPart + scale1 * to1[3]);
   Result.Normalize;
end;
jdelauney commented 6 years ago

sorry just back from shopping

No, problems Peter :) 👍

jdelauney commented 6 years ago

Ok changed Omega :=(Omega + DegToRadian(Spin*360)); by Omega :=(Omega + DegToRadian(Spin*180)); it's remove the negation of Z but now it's failed on W'sign :'(

jdelauney commented 6 years ago

Here we need determine if is -180 or 180, i think. It seem depend of the quarter part where is the "Spin". (sorry don't say how to explain, clearly, what is in my mind)

dicepd commented 6 years ago

ok Plugged that in, get a quat of [ 0, 0, -0.3826834, 0.9238795 ] which in converted to euler angle is [ x: 0, y: 0, z: -44.9999976 ] now to work out what quad that is half of to me it looks like -45 which would be 315 , now how that relates to (90 + x*pi) / 2 will have to find out. X being unkown as yet.

Using this site as checks https://www.andre-gaschler.com/rotationconverter/

jdelauney commented 6 years ago 
     omega:=GLZMath.ArcCos(cosom);
      //Omega :=(Omega + DegToRadian((Spin*180)+180));
      sinom:=1/Sin(omega);
      phi:=omega;// + Spin * cPi; // Try to negate 
      scale0:=sin(omega - t * phi) * sinom;
      scale1:=sin(t * phi) * sinom;

something is not clear here with the sign of omega or the range formed by Scale0 and Scale1

jdelauney commented 6 years ago

if i make` Omega :=(Omega + ((90 + Spin cPi) 0.5));

if i make Omega :=(Omega + ((-90 + Spin cPi) 0.5));

dicepd commented 6 years ago

Try the online calc quaternions can have differing signs but be the same rotation. Due to the half angle nature of quat they can make unique rotations in the -2pi to +2pi. Only when you convert back to euler can you see the equivalence. which make the things so dammed hard to test.;)

jdelauney commented 6 years ago

By spining we must find the direction for know on what quarter we are. not ?

jdelauney commented 6 years ago

ok if i understand we must negate value of axis x,y or z if the angle is less than zero and thats it. i'm right ? Your SlerpSpin make me confused

Axis x 1 y 1 z 0 Angle (degrees) -45 result = [ -0.2705981, -0.2705981, 0, 0.9238795 ]

Axis x 1 y 1 z 0 Angle (degrees) 45 result = [0.2705981, 0.2705981, 0, 0.9238795 ]

Axis x 0 y 0 z 1 Angle (degrees) -45 = [ 0, 0, -0.3826834, 0.9238795 ] Axis x 0 y 0 z 1 Angle (degrees) 45 = [ 0, 0, 0.3826834, 0.9238795 ]

dicepd commented 6 years ago

(720 - 90) / 2 = 315 so perhaps we are doubling the spin rotation and negating the initial rotation. this is still from the unchanged code you posted, will plug some more numbers in to test

dicepd commented 6 years ago

Ok that is it from changing 90 to 20 we get 720 - 20 = 700 : 700 / 2 = 350 and I get a euler answer of -10

dicepd commented 6 years ago

so need to half spin value and negate initial somehow.

dicepd commented 6 years ago

Quat to Mat4 done and checked in. Now for Mat4 to Quat.

jdelauney commented 6 years ago

Always don't understand SlerpSpin.

I don't understand "(720 - 90) / 2 = 315" = ((360*2) - 90) / 2

Why 360 our destination is at 90°

I don't understand slerpSpin in the same maner like you.

why in

if ((1.0-cosom)>cEpsilon) then
   begin
      //Cosom := Cosom
      omega:=GLZMath.ArcCos(cosom);

RadiansToDeg(Omega) = 45 and not 90 like the initial aqt2.Create(90,ZVector);

at firts i'm understading SlerpSpin like this :

so normally we have made a rotation of 90° = (90° 2) 0.5 or something escape me

And for Matrix test

jdelauney commented 6 years ago

if this is the case an Angular lerp is needed here

jdelauney commented 6 years ago

Ok by making like i think SlerpSpin is ok

First change : 

procedure TGLZQuaternion.Create(const angle  : Single; const axis : TGLZAffineVector);
//procedure TGLZQuaternion.Create(const angle  : Single; const axis : TGLZVector);
var
   f, s, c : Single;
   vaxis : TGLZVector;
begin
   GLZMath.SinCos(DegToRadian(angle*cOneHalf), s, c); //--> Angle div by 2 here
   Self.RealPart:=c;
   vaxis.AsVector3f := axis;
   vAxis.w :=1; // -----> Need set as affine
   f:=s/vAxis.Length;
   Self.ImagePart.V[0]:=axis.V[0]*f;
   Self.ImagePart.V[1]:=axis.V[1]*f;
   Self.ImagePart.V[2]:=axis.V[2]*f;
end;

after

function TGLZQuaternion.Slerp(const QEnd: TGLZQuaternion; Spin: Integer; t: Single): TGLZQuaternion; begin Result := Self.Slerp(Qend,t*spin); end;

and for finish : 

procedure TQuaternionFunctionalTestCase.TestSlerpSpin;
begin
//   aqt1.AsVector4f := WHmgVector;  // null rotation as start point.
   aqt1.create(1e-14,ZVector);  // null rotation as start point.
   aqt2.Create(90,ZVector); // 90  = 90
   aqt4 := aqt1.Slerp(aqt2, 2, 0.5); //  90 [ 0, 0, 0.7071068, 0.7071068 ]
   AssertEquals('SlerpSpin:Sub1 X failed ', 0.0, aqt4.X);
   AssertEquals('SlerpSpin:Sub2 Y failed ', 0.0, aqt4.Y);
   AssertEquals('SlerpSpin:Sub3 Z failed ', 0.7071068, aqt4.Z);
   AssertEquals('SlerpSpin:Sub4 W failed ', 0.7071068, aqt4.W);
   aqt4 := aqt1.Slerp(aqt2,2,2/9); // 40  [ 0, 0, 0.3420185, 0.9396932 ]
   AssertEquals('SlerpSpin:Sub5 X failed ', 0.0, aqt4.X);
   AssertEquals('SlerpSpin:Sub6 Y failed ', 0.0, aqt4.Y);
   AssertEquals('SlerpSpin:Sub7 Z failed ',  0.3420185, aqt4.Z);
   AssertEquals('SlerpSpin:Sub8 W failed ', 0.9396932, aqt4.W);
   aqt4 := aqt1.Slerp(aqt2,2,8/9); // 160  [ 0, 0, 0.9848078, 0.1736482 ]
   AssertEquals('SlerpSpin:Sub9 X failed ',   0.0, aqt4.X);
   AssertEquals('SlerpSpin:Sub10 Y failed ',  0.0, aqt4.Y);
   AssertEquals('SlerpSpin:Sub11 Z failed ', 0.9848078, aqt4.Z);
   AssertEquals('SlerpSpin:Sub12 W failed ', 0.1736482, aqt4.W);
end;

This is the behaviour you want ?

dicepd commented 6 years ago

Ok I just did the last two parts and it works, you have changed the spin to mean rotation * spin which is fine by me. previous code looked like it was adding some multiple of pi.

vAxis.w :=1; // -----> Need set as affine

This has no effect, with it or without it. Best not putting it in as it may make us add some asembler code where not required.

dicepd commented 6 years ago

And for Matrix test "ConvertToMatrix:Sub1 m33 failed " expected: <0,8936> but was: <1,0026>

Fixed

dicepd commented 6 years ago

thinking about it need to test for large number of 'spirals' and see what the limits are, if any.

jdelauney commented 6 years ago

If it not working just go on Bin, not really necessary. If you see others functions like this don't hesitate to trash it

dicepd commented 6 years ago

have you seen these quad_tess4x3 tri_tess3

patterns before for creating tri-strips without degenerates?

jdelauney commented 6 years ago

I see the problem, but I do not understand what you want to talk about and about what. Can you be more precise