Case Studies

Game Mechanics
#1

In current version of JOML, class Vector3f


    public float angleCos(Vector3f v) {
        double length1 = Math.sqrt(x * x + y * y + z * z);
        double length2 = Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
        double dot = x * v.x + y * v.y + z * v.z;
        return (float) (dot / (length1 * length2));
    }

   public float angle(Vector3f v) {
        float cos = angleCos(v);
        // This is because sometimes cos goes above 1 or below -1 because of lost precision
        cos = Math.min(cos, 1);
        cos = Math.max(cos, -1);
        return (float) Math.acos(cos);
    }

And Vector4f looks the same as the 3D (just extended). The problem is 2D regardless of the number of dimensions, so you want to try to work it as a 2D problem. The version is slower and has more error than is needed.

This is the exact sub-problem is the starting point of formulating slerp and the analysis of quaternion lerp. The only difference between 2D problems in 2D vs. more than 2 is you have to find the plane. (and put it back in place if needed)

#2

You forgot your point.

#3

This solution is lost in the equation. I’m asking someone to reason about the meaning which is geometric.

#4

Likewise for the divide and conquer nlerps. 2d and yet they are manipulating 4d values. Also bisected an angle is reducible. Say think a PBR half vector.

#5

I never knew about this vector-replacement for finding the cosine of an angle. Great stuff. I couldn’t find a version for Math.sin on google which is a pity.

I think the additional optimisation that @Roquen is suggesting is that the two square roots can be replaced by one.
Since sqrt(x)sqrt(y) = sqrt(xy)

Old:


    public float angleCos(Vector3f v) {
        double length1 = Math.sqrt(x * x + y * y + z * z);
        double length2 = Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
        double dot = x * v.x + y * v.y + z * v.z;
        return (float) (dot / (length1 * length2));
    }

New:


    public float angleCos(Vector3f v) {
        double lengthSq1 = x * x + y * y + z * z;
        double lengthSq2 = v.x * v.x + v.y * v.y + v.z * v.z;
        double dot = x * v.x + y * v.y + z * v.z;
        return (float) (dot / Math.sqrt(lengthSq1 * lengthSq2));
    }
#6

Useful tip, CommanderKeith, I’ll certainly be using that in the future, thanks.

Would it be taking it too literally to be only using X, and Y, values of any size vector, in an attempt to make the problem 2D?

Edit: Removed previous post because it was wrong and just taking up space.

#7

Sin can be calculated either using the cross product between two vectors if you need the sign as well, or if you don’t need the sign by simply using the fact that (cos^2 + sin^2 = 1), which translates to [icode]double sin = Math.sqrt(1 - cos*cos);[/icode], but this is less useful.

The cos = dot-product can be used for lots of useful things. A common appliance is field-of-view testing for AI. Given an NPC that is facing a certain angle, you can precompute its normalized view vector each frame. If his field-of-view is 30 60 degrees, you can calculate cos(30). To check if someone is in the NPC’s field-of-view:


//Precomputed:
Vector3 normalizedViewVector = ...;
float cosCutoff = (float)Math.cos(30);

//For each enemy:
Vector3 vectorToEnemy = target.position - me.position;
boolean canSee = dot(me.normalizedViewVector, vectorToEnemy) / vectorToEnemy.length() > cosCutoff;
if(canSee){
    ...
}

This has extremely good performance as it only uses conventional math plus a single square root to calculate length(). The precomputation is generally cheap as it only needs to be done N times per frame, while the enemy testing may be done up to N^2 times.

EDIT:
The dot-product calculates the cosine of the angle from the view vector, so for a 90 degree field of view, you calculate cos(45). This even works for values over 90 degrees, i.e to give the NPC a blind spot directly behind it.

#8

I had in mind geometric reasoning…but this is all headed in the right direction. Gosh we have |a||b|cos(t), if we could cheaply find |a||b|sin(t) we could use atan2.

vectorToEnemy.length() is positive.

#9

Why is atan2() so much better than acos()?

#10

Short answer atan & atan2 are (vs. acos and asin):

  1. faster to compute
  2. don’t have input domain issues
  3. are more robust.

I can give some details if you want.

#11

Is it possible to make arc-cos, arc-sin and arc-tan using vectors too?

The alternative cos and sin implementations with vectors are very interesting. I suppose that tan can be calculated as sinAngle/cosAngle which reduces to the cross/dot.

By the way, one disadvantage of the new angleCos method that does just one square root is a higher chance of numerical overflow if the x and y coordinates are very large. This is because they get multiplied together one extra time before being square rooted. I believe that many robust topographical API’s such as JTS (https://en.wikipedia.org/wiki/JTS_Topology_Suite) first ‘normalise’ the coordinates by shrinking or growing the vectors to manageable distances first. Given that JOML probably favours performance over edge-case stability, perhaps it’s an acceptable trade-off.

PS: thanks for the detailed explanation theagentd

#12

Hmm, interesting. 2 and 3 I can understand, but why is it faster to compute?

#13

Short answer: http://www.wolframalpha.com/input/?i=Plot[{ArcCos(x)%2C+ArcTan(x)}%2C+{x%2C0%2C1}]

I can expand on this if desired.

#14

Pseudo-code for more robust computation (but we’re not dropping real missiles into real teacups):


{
  ma = magnitude(a);
  mb = magnitude(b);
  y  = magnitude(a*mb - b*ma);
  x  = mangitude(a*mb + b*ma);
  return 2 * atan2(y, x)
}
#15

OK: we’re at we could use atan2 and compute the ‘y’ with |cross(a,b)| and ‘x’ with dot(a,b). ‘y’ seems reducible. The inside of the root looks like: dot( cross(a,b), cross(a,b) )

#16

I’ll answer my own question. If we look at a table of basic vector identities one will be:

(a×b) · (c×d) = (a·c)(b·d) - (b·c)(a·d)

Neat the cross product disappeared. And we have:

(a×b) · (a×b) = (a·a)(b·b) - (b·a)(a·b) = (a·a)(b·b) - (a·b)2

tying all of this together into pseudo-code:


px = dot(a,b);                       // scaled projection in 'x'
py = sqrt(dot(a,a)*dot(b,b)-px*px);  // scaled projection in 'y'
return atan2(py, px);

the update JOML version (with CommanderKeith’s sqrt comment) in pseudo-code:


t = sqrt(dot(a,a)*dot(b,b));  // composed scale of a and b
t = dot(a,b)/t;               // cos(angle)
t = clamp(t,-1,1);            // clamp to legal range to account for compounding errors
return acos(t);

One thing to note here is that both of these work in any number of dimension(1). WHAT? Notice how the cross product above disappeared? The reason is simple:

THIS IS A 2D PROBLEM: SO WHY AREN’T WE WORKING A 2D PROBLEM???

(1) Ok…there are a couple of assumptions here, but not important for this monologue.

(EDIT: forgot the acos return for the JOML version)

#17

Interesting, thanks for explaining.
An additional computation saving might be making an atan2Sq which gives the square atan2 and therefore avoids the square root, and adds just a multiplication of the dot product in the numerator.
As agentd explained, this can be compared to the atan2 squared of the NPC view angle

#18

Something I hinted at but didn’t say was:


boolean canSee = dot(me.normalizedViewVector, vectorToEnemy) / vectorToEnemy.length() > cosCutoff;

which I said: “vectorToEnemy.length() is positive.”

This is because is the length is positive, you can multiply both sides by it without changing the comparison. Then you can loose the root.

#19

I’m going to now very quickly repeat everything above but directly in 2D. If someone wants an informal visualization aid then let me know. We have two vectors ‘a’ & ‘b’ (directions or coordinates) in N dimensions (but let’s think 3…the number doesn’t matter). We visualize the vectors in the standard way: tails starting from some origin and their magnitude as a length with little arrows drawn at the end. Assuming that ‘a’ and ‘b’ are not pointing in the same or opposite directions (in the visualization they don’t fall on the same line), then there is a plane which contains both of them. We’re using ‘a’ as our reference direction so it falls on the positive x-axis. For our positive ‘y’ axis we’re going to choose the direction that contains ‘b’. We’re forcing the visualization of ‘b’ to be above the x-axis (in the first or second quadrant). From the head of ‘b’ (the tip of the arrow), we can draw a vertical line to the x-axis and label the point it cross as ‘px’. Likewise we can draw a horizontal line and label the point it cross the y-axis as ‘py’. These values are: px = parallel projection of ‘b’ onto the x-axis and parallel projection of ‘b’ onto the y-axis. These two values give us a local 2D coordinate of ‘b’.

So in this local 2D coordinate system we have (assuming b is unit, |b|=1):

a = |a|(1,0)
b = |b|(px,py)

if ‘b’ isn’t unit then the |b| is already contained in our px and py. Showing it this way because we can also express ‘b’ as:

b = |b|(cos(t), sin(t))

From this point of view what the trig operations are really saying is: cos(t) = how much is in ‘x’ and sin(t) = how much is in ‘y’. Now I’m really still making this overly complicated because I attempt to give a picture of how the pieces fit together. Now we’re ready for the punch-line:

Let X=|b|px, Y=|b|py and H=|b|

X2+Y2=H2 or (s cos(t))2+(s sin(t))2 = s2

So we’ve done a long dog-and-pony show to derive: Pythagorean theorem

The Pythagorean theorem doesn’t care about the scale of things. The relation is independent of uniform scaling. Given two of values we can find the third. The only higher order math notion required is the dot product which is typically shown as:

(a·b) = |a||b|cos(t)

which is a parallel projection(1) of b onto a , which is a scaled parallel projection onto the local x-axis. This extra scale factor of |a| doesn’t matter since the Pythagorean theorem doesn’t care about it. That gives us a ‘how much in x’ and we know the magnitude of ‘b’ so our unknown is ‘how much in y’. Solve for Y:

X2+Y2 = H2 -> Y2 = H2 - X2

X2 =(a·b)2
H2 = (a·a)(b·b)

Y2 = (a·a)(b·b) - (a·b)2

H may not be immediately clear. (b·b) is the magnitude squared of ‘b’. The magnitude of ‘a’ squared (a·a) is being multiplied in to account for the extra scaling on X.

Now for some other parts. If we wanted to know the unit vector in the direction of the positive ‘x’ axis then it’s simply: a/|a|. Now what about the same for ‘y’? Let’s do the dog-and-pony vector version.

(a×b) would give us a vector orthogonal to the local plane. Then if we cross that again with ‘a’ would give us the direction of ‘y’. So: (a×b)×a.

This is a form of the vector triple product: (a×b)×c = (a·c)b - (b·c)a

So our ‘y’ direction is: (a×b)×a = (a·a)b - (b·a)a = (a·a)b - (a·b)a

if |a| = 1, then (a·a)=1 and it would be: b - (a·b)a

the (a·b) is “how much in X” which we multiply by ‘a’ to get it back into a direction instead of a scalar. And we’re subtracting that from ‘b’. The end result is we’re taking ‘b’ and removing the part of it that’s in the ‘x’ direction. Or in other words we’re zeroing out the x-component in the local coordinate system. When |a| isn’t 1 then the extra (a·a) is accounting for the scaling introduced by (a·b). The magnitude of the computed Y direction is Pythagorean result above.

Yeah I could have made this easier to understand, but it’s already long enough. If something isn’t clear and you care then ask.

(1) Or ‘a’ onto ‘b’…the projection is the same.

#20

Ha Ha! Jerkface. Made some dumb-ass mistakes trying to bang that out quickly. Each mistake found before I correct the post get’s a medal!