haegarr

Please notice that updating an entity by simply iterating its components may not be sufficient. The overall architecture may need that specific (by type) components need to be updates before others, and that some components may need more than a single update step (look e.g. at this excerpt from Jason Gregory's Game Engine Architecture). Hence I've gone the way of sub-systems where belonging components are to be registered, and updating is left to be run by the sub-system. This doesn't exclude components to have their own update functions if necessary. Furthermore, there may be a dependency in the update order. Think e.g. of parenting the placement of models. This kind of problem can be solved using a dependency graph inside a sub-system where appropriate.

Error is in cross-product: Replace this line

                retVec.y = v1.z*v2.z - v1.x*v2.z;

with this line

                retVec.y = v1.z*v2.x - v1.x*v2.z;

BTW: Your approach will fail if dv and tv become very close, because the axis vanishes but will be enforced to unit length. You then will get more or less arbitrary axes. You have to check for this and simply set the new direction vector directly to tv in that case.

Another thing. when you build up your matrix from the base componants you should probely convert the rotation into 3 separate rotations (1 for x, 1 for y and 1 for z) and convert  each of them into quaternions and then multiply the 3 resulting queternions together to get 1 rotation again that does all 3 rotations in one instead of first x then y then z. This should prevent the anoying gimbal lock effect.

But it doesn't: Gimbal lock is an issue of concatenating rotations so that a rotation axis gets coincident with an already used rotation axis. It is not an issue of whether one uses matrices or quaternions for rotation representation. Hence the above suggestion doesn't help against gimbal lock but makes things needlessly more complex.

Furthermore, there is no difference (besides perhaps efficiency and rounding errors) whether one applies particular rotations one by one or else computes the composition matrix first and applies that:
( R₁ * R₂ ) * v == R₁ * ( R₂ * v )
a.k.a. associativity

The problems with timing are these:
* You want a minimum rendering rate to give the player the illusion of a smooth flow.
* You want to have a fixed timestep for physic simulations, because the underlying math works fine with fixed timesteps (only few methods work with variable timestep).
* You don't want to waste time; e.g. AI isn't required to decide 50 times per second what to do.
* You don't want input to be missed, because that would frustrate the player soon.

To handle all this in a single loop with unknown rate, you have to measure the time that has elapsed. E.g.

while ( running ) {
    now = GetSystemTime() ;
    elapsed = now - justNow;
    // do something
    justNow = now;
}

So, knowing the current moment in time "now", the previous moment in time "justNow", and hence the elapsed period "elapsed", you can drive your simulation. A fast computer may "do something" at 0.2, 0.4, 0.6, ... while a slow computer way "do something" at 0.3, 0.6, 0.9 ... seconds. However, both will reach e.g. 0.6 seconds, the former one within 3 steps and the latter one within 2 steps.

Now think of your simulations being defined in real time units. E.g. an animation with its timeline given in seconds. Assume that the animation has a keyframe with position P1 at 0s, and another keyframe with position P2 at 6s. The fast computer asks for a position at times 0.2s and 0.4s which need to be interpolated as usual, and then it hits P2 at 0.6s directly. The slow computer asks for the position at 0.3s which needs to be interpolated, and then hits P2 at 0.6s directly. Notice that both computers use P2 at the same time. Although the computers run differently fast, the animation is played in the same rate.

In reality the "elapsed" time will differ from loop pass to pass, but that still works fine, assuming that a minimum rate is reached at all.

Now coming to the fix timestep of physic simulation. The principle is simple: You use an inner loop that passes as often as the fixed timestep fits into the current "elapsed".

delta = 0;
while ( running ) {
    now = GetSystemTime() ;
    elapsed = now - justNow;

    delta += elapsed;
    while ( delta > 0 ) {
        // do physic simulation step
        delta -= physicsTimestep;
    }

    justNow = now;
}

The devil is in the details, of course. E.g. it is a sheer fortuity if delta will be reduced to 0. Instead, it will be slightly less than zero in general. Then you can do an interpolation of the current physic state and the previous one in order to compute the active state.

It computes the product
M := M * M_f
where M is the matrix on the top of the current matrix stack, and M_f is the matrix given as argument to the function. Notice that the argument matrix is on the right like any matrix generated by glTranslate, glRotate, glScale, ..., too. In other words: It multiplies the argument matrix with the matrix on the top of the stack and replaces the top matrix with the product.

IMHO the shown adopted code snippet is missing a normalization of both n2 and t2 (assuming that "using V,N,T as regular" doesn't include it). The code snippet probably only works well if Bones.z == 1 and Bones.y < 0. In all other cases the lengths of n2 and t2 may be anything but are later expected to be 1.

 bMat *= oMat;
 bMat += oMat2;

Why are you adding? All matrix chaining transformations should be multiplications...

Please let me clarify what
M_T := R_R * T_T(-c) + T_T(c)
actually means. As mentioned my post above, R_R is a 3x3 matrix and T_T is a 1x3 matrix (a.k.a. column vector). Multiplying a 3x3 matrix on the left and a 1x3 matrix on the right gives you a 1x3 matrix. Adding a 1x3 matrix onto a 1x3 matrix gives you a 1x3 matrix.

So notice that the result M_T is a 1x3 matrix (and that M_R is a 3x3 matrix), while M itself is a usual homogeneous 4x4 matrix. The correct assembly then looks like

bMat.make_identity();
// MR
bMat.element(0, 0) = cos(rx) * cos(ry);
bMat.element(1, 0) = -cos(rz) * sin(rx) - cos(rx)*sin(ry)*sin(rz);
bMat.element(2, 0) = cos(rx) * cos(rz) * sin(ry) + sin(rx) * sin(rz);;
bMat.element(0, 1) = cos(ry) * sin(rx);
bMat.element(1, 1) = cos(rx) * cos(rz) - sin(rx)*sin(ry)*sin(rz);
bMat.element(2, 1) = cos(rz) * sin(rx) * sin(ry) - cos(rx) * sin(rz);
bMat.element(0, 2) = -sin(ry);
bMat.element(1, 2) = cos(ry) * sin(rz);
bMat.element(2, 2) = cos(ry) * cos(rz);
// MT = MR * TT(-c) + TT(c)
bMat.element(0, 3) = bMat.element(0, 0) * (-ox) + bMat.element(0, 1) * (-oy) + bMat.element(0, 2) * (-oz) + ox;
bMat.element(1, 3) = bMat.element(1, 0) * (-ox) + bMat.element(1, 1) * (-oy) + bMat.element(1, 2) * (-oz) + oy;
bMat.element(2, 3) = bMat.element(2, 0) * (-ox) + bMat.element(2, 1) * (-oy) + bMat.element(2, 2) * (-oz) + oz;

if I have interpreted the indexing scheme correctly.


EDIT: It is for sure possible to compose the desired rotation simply by computing T(c) * R * T(-c). The above way just shows (as mentioned) the minimal computational effort to do; it avoids all that nasty scalar products with 0 and 1. However, this kind of optimization will probably not be noticeable.

Assuming that you use column vectors, if T(c) is the translation matrix to the center of rotation and R is the rotation matrix, then the product
M := T(c) * R * T(-c)
defines the desired new rotation matrix. It translates by -c, rotates, and translates back by c like you do now separately.

EDIT:
You can take advantage from knowing the structures of the matrices when you compute the above product. Think of M having a 3x3 sub-marix M_R on the upper left and a 1x3 sub-matrix M_T on the upper right, and doing so with R and T as well, then
M_R := R_R
M_T := R_R * T_T(-c) + T_T(c)
means the minimum computations to be done for yielding in M.

First notice that there are point vectors P and Q_i in use, as well as difference vectors p₁, p₂, and q_i. The correspondence between point and difference vectors is that a difference vector is the difference of 2 point vectors. E.g.
d := D₂ - D₁
what means that the difference vector d is the difference from point D₁ to point D₂. One can obviously rewrite the above equation to yield in
D₂ = D₁ + d
what means that point D₂ is reached when adding the difference vector d onto point D₁.

But how to reach the point D₃ exactly on the half way between D1 and D2? Well, we have to add the half of vector d:
D₃ = D₁ + 0.5 * d
In general we can add any elongation of d to D₁ by weighting d with a scalar s.
D( s ) := D₁ + s * d
with the special cases
D( s=0 ) = D₁ + 0 * d = D₁
D( s=1 ) = D₁ + 1 * d = D₂

The above are ray equations. They describe a directed line starting from D₁ in direction of D₂, where D can be any point on that line in dependence of s. That is exactly what the right side of the equation in the paper is (but using its own variable names).

Now, think of not moving on a 1-dimensional line but on a 2-dimensional plane. You can reach any point of a plane if you 1st wander along a ray as done above, and 2nd wander on another ray but now starting from the point you've reached during your previous wandering and, of course, wander into another direction. This would look like
D( s, t ) := D( s ) + t * d₂ = D₁ + s * d₁ + t * d₂
and that is exactly what is written on the left side of the article's equation (but using its own variable names).

In summary, the equation means a point that is both (a) on the plane in which triangle A lies and is (b) on a ray in which the edge of a 2nd triangle B lies.

I thought that in 4x4 matrices a tranlation was possible through multiplication... but in truth I haven't even got an negation function.... Nor had I thought about them....

Well, things are more complicated. Assume a 2D space in which to do transformations. Assume that you have a 2x2 matrix R for rotation and a 2 vector t for translation. Using row vectors, a transformed point P' then is computed from its original point P as
P' := P * R + t
Writing this as set of scalar equations gives
P'_x := P_x * R_xx + P_y * R_yx + t_x
P'_y := P_x * R_xy + P_y * R_yy + t_y

If we want to express the translation within the matrix product, the t_x and t_y must become part of R, but then we need to multiply it with another component of the point P. This is because In a matrix product the count of columns in the left matrix (the vector P in this case) and the count of rows in the right matrix (R in this case) must be equal. This additional component need to be 1, because e.g. t_x * 1 == t_x doesn't change anything:
P'_x := P_x * R_xx + P_y * R_yx + t_x * 1
P'_y := P_x * R_xy + P_y * R_yy + t_y * 1

But now R is no longer square and P' has 2 components but P has 3. To overcome this and to allow the resulting P' to be used in a further transformation, the same extension under consideration of the mechanics of the matrix product shows us that
P'_x := P_x * R_xx + P_y * R_yx + t_x * 1
P'_y := P_x * R_xy + P_y * R_yy + t_y * 1
P'_w  := 1 =  P_x * 0 + P_y * 0 + 1 * 1
will do the job. Hence, using such an extension into a 3rd dimension allows us to write
P'_h := P_h * M_h = P_h * R_h * T_h

This additional co-ordinate is named the homogeneous co-ordinate, and vectors or matrices using this are named homogeneous vectors or matrices, resp.

So: It is not only an additional dimension that allows to incorporate a translation as matrix product. It is further a special use of this co-ordinate. You'll see during your journey that there are more things with this homogeneous stuff like differing between point and direction vectors, normalized and unnormalized points, and its use for perspective projection.

Here are the methods for multiplication and inversion I wrote (just in case there are obviouse errors)...
The multiplication function is just a standard row dot column loop,
the inversion is using  Gauss-Jordan Elimination (well; my understanding of it).
Neither are written for speed... just trying to put together a working class.

It is just so that the unary minus is commonly used for negation, and I've interpreted its occurrence just that way. Using it as matrix inversion is confusing most people, and that is a Bad Thing for a library.

Some more things:

1. It is Good Practice to make non-member functions for operations that are not directly related to a class. E.g.
float dot( vec2f const&, vec2f const&) instead of float vec2f::dot( vec2f const& ) const
vec3f cross( vec3f const&, vec3f const& ) instead of vec3f vec3f::cross( vec3f const& ) const

2. The cross product is not commutative. The rule is
v1 x v2 = -( v2 x v1 )
Now, with the usual order X Y Z of axes, this means that
X x Y = Z
Y x Z = X
Z x X = Y
but especially
Z x Y = -X
X x Z = -Y
Please check your code snippet w.r.t. this issue.

--Where are the text mark-up options gone to? --


Don't confuse inversion with negation. Inversion is meant w.r.t. an operation, and its outcome is also dependent on the operation. E.g. negation is the inversion for addition, so that
f + -f == 0
while the reziprocal is the inversion of the multiplication, so that
f * f^-1 == 1
(I'm not a mathematician but I hope you understand what I mean.)

For matrices both kinds exists, and "inversion" is usually meant in the context of multiplication:
M * M^-1 == I


The order of multiplication of matrices depends on whether you use column or row matrices, because of
( M1 * M2 )^t == M2^t * M1^t
where t denotes the transpose operator. As you can see, the order reverses if you choose "the other" convention.

Besides this, the relative order of transformations depends on the desires of the user. In general one wants to transform by using the local axes. This implies the common order
S * R * T
(when using row vectors), where S means scaling, R menas rotation, and T means translation.

Often R is composed using Euler angles. There is a couple of rotations known as Euler angles. A typical usage is pitching, then heading, then banking. Because you want to do so around local axes, the formula using row vectors would be (I neglect the details here)
R := B * H * P

Sometimes you want to use arbitrary origins or axes. In such cases you need to multiply on both sides of the original transformation matrix (I neglect the details here once again). Eg. to use an arbitrary center of rotation, you need to use
C^-1 * R * C

In summary: A library must provide both multiplication on the left and on the right side. And it has to be absolutely clear whether column or row vectors are in use!

Forget about Game Engine and APIs. Its almost clear now.
Just answer my two questions :-

Excuse me coming back to this ... I think it is not clear. AFAIK:

An API is an interface. I.e. it is a definition of how to invoke the functionality of, e.g., a game engine. As that, an API is never a game engine itself. A game engine may have several APIs. Further, an API can exist without any backing software.

OpenGL, strictly speaking, is an API for a graphics rendering engine. It defines how an application can interact with and what it can expect from a piece of software that e.g. a graphics card manufacturer provides. The said implementation of the software is not OpenGL but a software for that the OpenGL API is implemented. One can use the same engine and put another API in front of it; however, that would no longer be OpenGL.

Is XNA a game engine?

I'm not familiar with XNA, so I don't contribute to answer this.

Are drivers the *.dll files which interact with the hardware?

Drivers interact with hardware, that's true. They came usually as DLL files (under Windows; other OS'es use other formats and naming), because such kind of external modules allow a separation of the functionality from the using application. E.g. a driver is implemented by the hardware manufacturer, while the using applications are implemented by several other parties.

I have to hint at this post in the 2nd of the threads cited in the OP ;)

Please notice also the message of the linked post: nlerp may be suited better than slerp. Its a question of the use case.

Maybe you want to read this article on lerp, nlerp, and slerp.

You'll be surprised how often a single set of vertex data may be copied until it is used for actual rendering. Remember that OpenGL enqueues its commands. What happens if you try to alter a VBO before its content was rendered? Your invocation may stall until the VBO is free to use.

To avoid such pitfalls, you have to use advanced techniques:

a. for GL3 is available or the extension GL_map_buffer_range: ring buffers with MapBufferRange
b. for GL2 on MacOS: ring buffers with extension APPLE_flush_buffer_range
c. else, e.g. for GLES: VBO double buffering

There is some in-depth discussion on ring buffering. Some thoughts about VBO double buffering (a 2 part article). But be warned: Those stuff is advanced... I don't suggest you to implement such a solution now, but reading that stuff may help to understand what happens under the hood.