A common approach is to define values of friction and restitution for each individual material and then combine them to get the material-material values. Given individual material friction values \(x\) and \(y\), we need to define a combination function \(f(x,y)\). Similarly, \(r(x,y)\) for restitution (using \(x\) and \(y\) to mean the material restitution values this time).

The value of \(f(x,y)\) and \(r(x,y)\) should lie between \(x\) and \(y\). So \(f(x,x) = x\) and \(r(x,x) = x\).

For any constant \(c\), \(f(x,c)\) and \(r(x,c)\) should be monotonically increasing (it shouldn't ever decrease as \(x\) increases). It should also avoid long flat sections, so as to be able to discriminate bewteen materials. For instance, putting boxes of ice and rubber on an plane and increasing the slope, the ice should slip first, whether the plane is made of ice or rubber.

\(f(x,y)\) should favour slippy over grippy, i.e. \(f(0,1) \lt 0.5\). This corresponds to the intuitive notion that ice should have a greater effect on the combined restitution value than rubber, e.g. vehicle tyres on ice.

I decided that I wanted \(r(x,y)\) to favour the extremes, in a similar way that \(f(x,y)\) should favour values towards \(0\). So very bouncy or very energy-absorbing materials predominate over average ones. In a game world, you then have the option of making the world less or more bouncy for all objects within it by changing the world material, while maintaining a reasonably wide range of restitution values when the world is set to an averagely bouncy material.

Candidates for \(f\) that are often used are the minimum, the arithmetic average, and the geometric average.

\(f_{min}(x,y) = min(x,y)\)

This has too many flat sections and doesn't adequately discrimiate between materials, e.g. \(f(0.1,0.1) = f(0.1,1)\).

\(f_{aa}(x,y) = \frac{x + y}{2}\)

This doesn't favour slippy over grippy enough.

\(f_{ga}(x,y) = \sqrt{{x}{y}}\)

This is good, but the values aren't equally spaced. For this reason, I have found an alternative function - a weighted sum of \(x\) and \(y\).

\(f_{ws}(x,y) = \frac{{x}{w_{x}} + {y}{w_{y}}}{w_{x} + w_{y}}\) where \(w_{x} = \sqrt{2}(1 - x) + 1\)

\(f_{ws}\) has a more regular spacing between values than \(f_{ga}\). The trade-off is that it doesn't cover the full range of values for \(f(x,1)\) (\(f(0,1) \approx 0.3\)). However, they are approximately equal for \(0.2 \le x \le 1\).

As with friction, the minimum, the arithmetic average, and the geometric average are often used.

\(r_{min}\) has the same objections as \(f_{min}\).

\(r_{ga}\) is not suitable because it would require the world material to have a restiution near \(1\) to provide a wide range of combined values.

\(r_{aa}\) is better, but it doesn't give a very wide range of values for \(r(x,0.5)\). We can improve on this range by allowing some flatness and defining \(r\) as a piecewise min/max/sum function.

\(r_{pmms} = \begin{cases}min(x,y) & \text{if }x \lt 0.5 \land y \lt 0.5\\max(x,y) & \text{if }x \gt 0.5 \land y \gt 0.5\\x + y - 0.5 & \text{otherwise}\end{cases}\)

This has the same shortcomings as \(r_{min}\) at the corners where the min and max functions are used. But you can't avoid that if you want to have the full range of values for \(r(x,0.5)\) and still satisfy \(r(x,x) = x\). Similar to the friction, I have created a weighted sum function that I think is better.

\(r_{ws}(x,y) = \frac{{x}{w_{x}} + {y}{w_{y}}}{w_{x} + w_{y}}\) where \(w_{x} = \sqrt{2}\left\vert{2x - 1}\right\vert + 1\)

As with \(f_{ws}\), \(r_{ws}\) sacrifices some of the range for \(r(x,0.5)\) to provide better continuity in the corners.

It's the choice that maximizes the function range, while maintaining monotonicity. These graphs show what \(f_{ws}\) looks like with \(w_{x} = c(1 - x) + 1\), for \(c = 0.5\), \(c = 1\), \(c = \sqrt{2}\), \(c = 2\), and \(c = 4\). For \(c \lt \sqrt{2}\), \(f_{ws}(x,1)\) has less range. For \(c \gt \sqrt{2}\), \(f_{ws}\) is no longer monotonic - near \(f_{ws}(0.1,0.9)\) it starts to curve back, making a more grippy material have a lower combined friction! You can see this more clearly for higher values of \(c\).

The method to combine material friction and restitution values may seem like a minor detail, but sometimes the devil's in the details.

Since the whole concept of defining values for friction and restitution for a single material and then combining them isn't physically accurate, this is obviously just a case of finding a suitable function rather than attempting to model reality. You may have different requirements for your functions, but I hope this discussion of the alternatives is useful.

I made the graphs at the WolframAlpha web site. I find it's a pretty useful tool. Here's an example plot.

There are several variations of the problem, where perhaps the projectile is a missile that needs to accelerate, or where the shooter is a turret that is currently aiming in some direction and needs time to aim somewhere else... We'll cover the simple case first, and then we'll present a general template for solving variations of the problem.

Let's assume the shooter can aim and shoot anywhere instantly, the target is moving at a constant velocity and the projectile will travel at a constant velocity too. We are given as inputs the target's current position, its velocity and the speed of our projectile. We'll use coordinates where the shooter is at the origin and has zero velocity.

As mentioned before, if we naively aim at the target's current position, by the time the projectile gets there, the target will have moved. We can compute how long it will take for the projectile to get to the target's current position, compute where the target will be then and aim there instead.

Position compute_first_order_correction(Position target_position, Vector target_velocity, float projectile_speed) { float t = distance(Origin, target_position) / projectile_speed; return target_position + t * target_velocity; }

This simple piece of code is probably good enough in many cases (if the target is moving slowly compared to the projectile speed, if the target is moving perpendicularly to the shooter-to-target vector, or if we want to sometimes miss because a more precise solution would be detrimental to the fun of the game).

For a more precise solution, you could iterate this first-order correction until it converges.

Position iterative_approximation(Position target_position, Vector target_velocity, float projectile_speed) { float t = 0.0f; for (int iteration = 0; iteration < MAX_ITERATIONS; ++iteration) { float old_t = t; t = distance(Origin, target_position + t * target_velocity) / projectile_speed; if (t - old_t < EPSILON) break; } return target_position + t * target_velocity; }

In the iterative approximation, we would stop if we found a place where

Let's do some computations to try to solve it.

t^2 * projectile_speed^2 = dot_product(target_position + t * target_velocity, target_position + t * target_velocity)

t^2 * projectile_speed^2 = dot_product(target_position, target_position) + 2 * t * dot_product(target_position, target_velocity) + t^2 * dot_product(target_velocity, target_velocity)

This is a second-degree equation in

// a*x^2 + b*x + c = 0 float first_positive_solution_of_quadratic_equation(float a, float b, float c) { float discriminant = b*b - 4.0f*a*c; if (discriminant < 0.0f) return -1.0f; // Indicate there is no solution float s = std::sqrt(discriminant); float x1 = (-b-s) / (2.0f*a); if (x1 > 0.0f) return x1; float x2 = (-b+s) / (2.0f*a); if (x2 > 0.0f) return x2; return -1.0f; // Indicate there is no positive solution } Position direct_solution(Position target_position, Vector target_velocity, float projectile_speed) { float a = dot_product(target_velocity, target_velocity) - projectile_speed * projectile_speed; float b = 2.0f * dot_product(target_position, target_velocity); float c = dot_product(target_position, target_position); float t = first_positive_solution_to_quadratic_equation(a, b, c); if (t <= 0.0f) return Origin; // Indicate we failed to find a solution return target_position + t * target_velocity; }

There are many variations of the problem we could consider: Accelerating targets, accelerating projectiles, situations where it takes time to aim at a new direction... All of them can be solved following the same template.

The things that could change can be encoded in two functions:

`position_of_target_at(time)``time_to_hit(position)`

We then compute where the target will be at time

Just as before, we could do a first-order correction, use iterative approximation or solve the problem directly. It might be the case that an analytical solution can be found, like we did in the previous section, but things can get messy quickly and you may have to resort to a numerical solution.

This article covered three methods to implement deflection in your games: First-order correction, iterative approximation and directly finding the solution. You'll need to use some judgement to decide which one to use. Hopefully this article gives you enough to make an informed decision.

Click Here to view presentation.

Firefox, Chrome or Opera recommended as they support animations; see Requirements below for details. Press

Transformations are a general mathematical concept that is applicable to anyone working in:

- Computer Graphics
- Animation
- Game Programming (2D and 3D)
- Image Processing
- Motion Capture
- UI Design

- Computer Vision
- Robotics
- Aeronautics
- Parallel Computing

Concepts discussed are dimension-indepedant; it's just easier to explain and visualize things in 2D but they're applicable to higher dimensions without loss of generality.

Instead of dealing with only elementary transforms (rotate, scale, translate) on points, which most resources do, it also covers:

- Basic math behind transforms (without matrices)
- Matrix introduced past the basics since it's just a tool

- Composite transforms — concatenation of multiple transforms
- Anti-commutativity
- Transforms about an arbitrary origin

- Active and passive viewpoints of transforms
- Transformation of coordinate systems
- Mapping between multiple coordinate systems
- Hierarchical Transforms

The presentation was hand-crafted using HTML5, CSS3, JavaScript and SVG; so nothing more than a browser is required to view the presentation. Firefox, Chrome or Opera, even if you've an older version, is highly recommended since support for SVG animations isn't that great in other browsers; with browsers like Safari, Edge or IE you will not be able to see these animations — caveat lector.

You can also view the presentation on a mobile or a tablet; the content, without any loss of fidelity, should get resized and fit to given form factor, thanks to vector graphics and CSS3. Touch (tap and swipe left/right) is supported both for navigation and animation control.

Transformations are better understood visually than with just dry theory. Hence every transformation workout is accompanied, along with the math, by an animation. Images with a graph sheet background contain animations. Simply click on such an image to cycle through the animation. Every click changes the image's state, animating it into a new figure. When you see the original figure you started with, it denotes the end of all the animation the image contains.

If you're using a keyboard to navigate the presentation, clicking the image steals focus from the presentation and sets in on to the image; simply click anywhere outside the image to give the focus back to the presentation. Now press *space* / *->* to continue with the presentation.

The solution to problem of mapping the boy space to street space that is posed at the end of the presentation is in *map_boy_to_street.md*.

The demonstrations/interactive animations embedded within the presentation was done in SVG, the XML-based, open standard file format for vector graphics with support for animation. In English, all it has are just instructions like *move to*, *line to*, *rect*, *circle*, *ellipse*, etc. in a very readable, XML format and no unintelligible binary data. So a reader (yes, a human one too) can easily read and understand it; a renderer can render it at any resolution without any loss in fidelity. The presentation's first slide has a 3D-transformed background which too is an SVG — it should show up as something similar to this; a simple check to see how well your browser supports SVGs and CSS3 transforms.

It's highly recommended that you fiddle with the SVGs under *images* directory. SVG format is very similar to PostScript (which also has commands like *move to*, *line to*, etc.) and is an excellent *Hello World* test bed (Short, Self Contained, Correct, Example) for learning transformations or 2D graphics in general. Oh they also have a tag for groupings *<g>* which may be used to learn hierarchical transformations. An SVG is only more readable than a PS, PDF or XAML. Just open it in a (modern) browser to view it (no, not Edge, it doesn't do 3D CSS3 transforms in SVGs yet :), open it in your favourite text editor, muck around, save and refresh your browser to see the changes immediately; rinse and repeat.

- Computer Graphics using OpenGL,
*Francis Hill*and*Stephen Kelley* - 3-D Computer Graphics,
*Samuel R. Buss* - Essential Math for Games and Interactive Applications,
*James Van Verth*and*Lars Bishop* - 3D Math Primer,
*Fletcher Dunn*and*Ian Parberry* - reveal.js, the presentation framework, generously shared under the MIT licence by
*Hakim El Hattab* - MathJax for rendering beautiful math equations on any browser,
*American Mathematical Society* - Elementary affine transforms chart shared under CC 3.0, used as first slide background,
*CM Lee* - What is or isn't a linear transform, shared under CC 3.0,
*Ldo* - Reveal.js on Github pages,
*Vasko Zdravevski*

Note: the method presented in this article is described in the paper "Interactive Dynamics" and was written by Witkin, Gleicher, and Welch. It was published in ACM in 1990.

A posted PDF of the paper can be found here: http://www.cs.cmu.edu/~aw/pdf/interactive.pdf

A link to another article by Witkin on this subject can be found here: https://www.cs.cmu.edu/~baraff/pbm/constraints.pdf

Everyone's familiar with Newton's second law: \( F = ma\). It forms the basis of Newtonian mechanics. It looks very simple by itself, but usually it's very hard to deal with Newton's laws because of the number of equations involved. The number of ways a body can move in space is called the

We need to devise an easy way to build and compute this system of equations. For a point mass moving in 3D, we can set up the general equations as a matrix equation:

\[ \left [ \begin{matrix} m_1 & 0 & 0 \\ 0 & m_1 & 0 \\ 0 & 0 & m_1 \\ \end{matrix} \right ] \left [ \begin{matrix} a_{1x} \\ a_{1y} \\ a_{1z} \\ \end{matrix} \right ] = \left [ \begin{matrix} F_{1x} \\ F_{1y} \\ F_{1z} \\ \end{matrix} \right ]\]

This can obviously be extended to include accelerations and net forces for many particles as well. The abbreviated equation is:

\[ M \ddot{q} = F \]

where \(M\) is the mass matrix, \(\ddot{q}\) is acceleration (the second time derivative of position), and \(F\) is the sum of all the forces on the body.

One of the problems with computing with Newton-Euler methods is that we have to compute all the forces in the system to understand how the system will evolve, or in other words, how the bodies will move with respect to each other. Let's take a simple example of a pendulum.

Technically, we have a force on the wall by the string, and a force on the ball by the string. In this case, we can reduce it to the forces shown and solve for the motion of the ball. Here, we have to figure out how the string is pulling on the ball ( \( T = mg \cos{\theta}\) ), and then break it into components to get the forces in each direction. This yields the following:

\[ \left [ \begin{matrix} m & 0 \\ 0 & m \\ \end{matrix} \right ] \left [ \begin{matrix} \ddot{q}_x \\ \ddot{q}_y \\ \end{matrix} \right ] = \left [ \begin{matrix} -T\sin{\theta} \\ -mg+T\cos{\theta} \\ \end{matrix} \right ]\]

We can then model this motion without needing to use the string. The ball can simply exist in space and move according to this equation of motion.

One thing we've glossed over without highlighting in the pendulum example is the string. We were able to ignore the fact that the mass is attached to the string, so does the string actually do anything in this example? Well, yes and no. The string provides the tension to hold up the mass, but anything could do that. We could have had a rod or a beam hold it up. What it really does is define the possible paths the mass can travel on. The motion of the mass is dictated, or

If we rearrange the pendulum's equations of motion, we can illustrate a point:

\[ \left [ \begin{matrix} m & 0 \\ 0 & m \\ \end{matrix} \right ] \left [ \begin{matrix} \ddot{q}_x \\ \ddot{q}_y \\ \end{matrix} \right ] = \left [ \begin{matrix} 0 \\ -mg \\ \end{matrix} \right ] + \left [ \begin{matrix} -T\sin{\theta} \\ T\cos{\theta} \\ \end{matrix} \right ]\]

In our example, the only applied force on the mass is gravity. That is represented as the first term on the right hand side of the equation. So what's the second term? That is the

\[ M_{ij} \ddot{q}_j = Q_j + C_j \]

where \(M\) is the mass matrix, \(\ddot{q}\) is acceleration (the second time derivative of position), \(Q\) is the sum of all the

The notion of adding constraint forces can be a bit unsettling because we are adding more forces to the body, which you would think would change the energy in the system. However, if we take a closer look at the pendulum example, we can see that the tension in the string is acting perpendicular (orthogonal) to the motion of the mass. If the constraint force is orthogonal to the displacement, then there is no additional work being done on the system, meaning no energy is being added to the system. This is called d'Alembert's principle, or the principle of virtual work.

As illustrated above, as long as the constraint forces don't do any work on the system, we can use them. There are lots of combinations of constraint forces that can be used that satisfy d'Alembert's principle, but we can illustrate a simple way to get those forces.

Usually a simulation has a starting position \(q = \left [ \begin{matrix} x_1 & y_1 & z_1 & x_2 & y_2 & z_2 & \cdots \\ \end{matrix} \right ] \) and velocity \(\dot{q} = \left [ \begin{matrix} \dot{x}_1 & \dot{y}_1 & \dot{z}_1 & \dot{x}_2 & \dot{y}_2 & \dot{z}_2 & \cdots \\ \end{matrix} \right ] \). The general constraint function is based on the state \(q(t)\) and possibly on time as well: \(c(q(t))\). The constraints need to be implicit, meaning that the constraints should be an equation that equals zero. For example, the 2D implicit circle equation is \(x^2 + y^2 - R^2 = 0\).

Remember there are multiple constraints to take into consideration. The vector that stores them will be denoted in an index notation as \(c_i\). Taking the total derivative with respect to time is:

\[\dot{c}_i=\frac{d}{dt}c_i(q(t),t)=\frac{\partial c_i}{\partial q_j}\dot{q}_j+\frac{\partial c_i}{\partial t}\]

The first term in this equation is actually a matrix, the Jacobian of the constraint vector \(J = \partial c_i/\partial q_j\), left-multiplied to the velocity vector. The Jacobian is made of all the partial derivatives of the constraints. The second term is just the partial time derivative of the constraints.

Differentiating again to get \(\ddot{c_i}\) yields:

\[\ddot{c}_i = \frac{\partial c_i}{\partial q_j} \ddot{q}_j + \frac{\partial \dot{c}_i}{\partial q_j} \dot{q}_j + \frac{\partial^2 c_i}{\partial t^2}\]

Looking at the results, the first term is the Jacobian of the constraints multiplied by the acceleration vector. The second term is actually the Jacobian of the time derivative of the constraint. The third term is the second partial time derivative of the constraints.

The formulas for the complicated terms, like the Jacobians, can be calculated analytically ahead of time. As well, since the constraints are position constraints, the second time derivatives are accelerations.

If we solve the matrix Newton's Law equation for the accelerations, we get:

\[\ddot{q}_j = W_{jk}\left ( C_k + Q_k \right )\]

where \(W = M^{-1}\), the mass matrix inverse. If we were to replace this with the acceleration vector from our constraint equation, we would get the following:

\[\frac{\partial c_i}{\partial q_j} W_{jk}\left ( C_k + Q_k \right ) + \frac{\partial \dot{c}_i}{\partial q_j} \dot{q}_j + \frac{\partial^2 c_i}{\partial t^2} = 0\]

\[JW(C+Q) + \dot{J} \dot{q} + c_{tt} = 0\]

Here, the only unknowns are the constraint forces. From our discussion before, we know that the constraint forces must satisfy the principle of virtual work. As we said before, the forces need to be orthogonal to the displacements, or the legal paths. We will take the gradient of the constraint path to get vectors orthogonal to the path. The reason why this works will be explained later. Since the constraints are placed in a vector, the gradient of that vector would be the Jacobian matrix of the constraints: \(\partial c/\partial q\). Although the row vectors of the matrix have the proper directions to make the dot product with the displacements zero, they don't have the right magnitudes to force the masses to lie on the constraint. We can construct a vector of scalars that will multiply with the row vectors to make the magnitudes correct. These are known as

\[C_j = \lambda_i \frac{\partial c_i}{\partial q_j} = J^T \lambda_i\]

Plugging that equation back into the augmented equation for Newton's law:

\[ \left ( -JWJ^T \right ) \lambda = JWQ + \dot{J}\dot{q} + c_{tt}\]

Note that the only unknowns here are the Lagrange multipliers.

If you're confused at how Witkin got that equation for the constraint forces, that's normal. I'll attempt to relate it to something easier to visualize and understand: surfaces. Let's take a look at the equation of a quadric surface:

\[Ax^2+By^2+Cz^2+Dxy+Eyz+Fxz+Gx+Hy+Iz+J=0\]

The capital letters denote constants. Notice also the equation is implicit. We can see the equation for an ellipse is a quadric surface:

\[f(x,y,z) = (1/a^2)x^2+(1/b^2)y^2+(1/c^2)z^2-1=0\]

For a point (x,y,z) to be on the surface, it must satisfy this equation. To put it into more formal math terms, we could say the surface takes a point in \(\mathbb{R}^3\) and maps it to the zero vector in \(\mathbb{R}\), which is just 0. Any movement on this surface is "legal" because the new point will still satisfy the surface equation. If we were to take the gradient of this ellipse equation, we'd get:

\[ \left [ \begin{matrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z} \\ \end{matrix} \right ] = \left [ \begin{matrix} \frac{2x}{a^2} \\ \frac{2y}{b^2} \\ \frac{2z}{c^2} \\ \end{matrix} \right ] \]

This vector is the normal to the surface at the given (x,y,z) coordinates. If we were to visualize a plane tangent to the ellipse at that point, the dot product of the normal and the tangent plane would be zero by definition.

With this same type of thinking, we can see the constraints as a type of algebraic surface since they're also implicit equations (it's hard to visualize these surfaces geometrically since they can be n-dimensional). Just like with the geometry example, if we were to take the gradient of the constraints the resulting vector would be orthogonal to a tangent plane on the surface. In more formal math terms, the constraints/surfaces can be called the

The purpose of these math terms are to help generalize this concept for use in situations where the problems are not easy to visualize or intuitively understand.

With these equations in place, the process of calculating the system state can now be summed up as follows:

- Construct the \(W\), \(J\), \(\dot{J}\), \(c_{tt}\) matrices.
- Multiply and solve the \(Ax=b\) problem for the Lagrange multipliers \(\lambda\).
- Compute the constraint forces using \(C = J^T \lambda\).
- Compute the accelerations using \(\ddot{q} = W(C+Q)\).
- Integrate the accelerations to get the new velocities \( \dot{q}(t) \) and positions \( q(t) \).

This process can be optimized to take advantage of sparsity in the matrices. The Jacobian matrix wll generally be sparse since the constraints won't generally depend on a large number of particles. This can help with the matrix multiplication on both sides. The main challenge will be in building these matrices.

Due to numerical error, there will be a tendency to drift away from the proper solution. Recall that in order to derive the modified Newton's law equation above, we forced \( c_i = 0 \) and \( \dot{c_i} = 0 \). Numerical error can produce solutions that won't satisfy these equations within a certain tolerance, so we can use a method used in control systems engineering: the PID loop.

In most systems we want to control, there are inputs to the system (force, etc.) and measurements of a system's state (angle, position, etc.). A PID loop feeds the error in the actual and desired state back into the force inputs to drive the error to zero. For example, the human body has many different inputs to the force on the muscles when we stand upright. The brain measures many different things to see if we're centered on our feet or if we're teetering to one side or another. If we're falling or we're off-center, the brain makes adjustments to our muscles to stay upright. A PID loop does something similar by measuring the error and feeding that back into the inputs. If done correctly, the PID system will drive the error in the measured state to zero by changing the inputs to the system as needed.

Here, we use the error in the constraint and the constraint derivative to feedback into the system to better control the numerical drift. We augment the forces by adding terms that account for the \(c_i =0\) and \(\dot{c_i}=0\):

\[ F_j = Q_j + C_j + \alpha c_i \frac{\partial c_i}{\partial q_j} + \beta \dot{c_i} \frac{\partial c_i}{\partial q_j} \]

This isn't a true force being added since these extra terms will vanish when the forces are correct. This is just to inhibit numerical drift, so the constants \(\alpha\) and \(\beta\) are "magic", meaning that they are determined empirically to see what "fits" better.

Witkin's method for interactive simulations is pretty widely applicable. Although this can be obviously used for movable truss structures and models of Tinkertoys, he also talks about using them for deformable bodies, non-linear optimization and keyframe animation as well. There are lots of applications of this method. Hopefully this showcase of Witkin's method will help make this interesting solution more accessible to anyone doing any type of simulation engines.

BACKGROUND

Having worked on a physics-library-based space shooter, I was acutely aware of the difficulties. First of all, you need concave polygons. Our spaceships weren't, so we had to divide them into a bunch of little triangles using a fancy algorithm. When the ships collided, the triangles would fly apart and snap back together as if held together with elastic, sometimes with part of the other ship wedged in between - not the effect we had in mind. And then there were issues with "bullet tunneling", and the general pain of storing part of our objects' state in a black box - even if it is an open-source black box.

If I had to do it all over again, I'd use simple bounding-circle collisions. That's exactly what the MARS guys did.

After that ordeal, and with HTML5 coming into vogue a few years back, I decided to make a nice, simple, slow-paced, sailing ship battle game... in JavaScript. I tried a JS port of the Chipmunk physics engine, but it was just too slow. So, I looked up the algorithms I really needed, and implemented them in about 100 lines of JavaScript. And when HTML5 soured on me and I switched to C++ for my next game, porting the code was a piece of cake. I love working this way.

DESIGN

There are three kinds of objects in my naval combat game:

1. Ships - a few slow-moving

2. Cannonballs - many fast-moving point particles

3. Land masses - a few large static polygons

I only need to check for ship-vs-ship, ship-vs-cannonball, and ship-vs-land collisions. What are the algorithms for those?

1. Ship vs Cannonball - polygon vs line segment.

2. Ship vs Land - point vs polygon. Shorelines are fuzzy and imprecise, so for this purpose I treat ships as point particles.

3. Ship vs Ship - polygon vs polygon, specifically

There's no need for joints, polygon stacking, or polygon-vs-polygon continuous collision detection (Box2D creator Erin Catto wrote some great papers about that; it's rather advanced for DIYers but definitely worth reading).

I also use bounding boxes for first-pass collision detection. That's just an optimization, and it's easy, so I'll leave it as an exercise to the reader.

Potentially I might need a spatial partitioning scheme (e.g. quadtrees) if I create an MMO with huge sprawling maps. I haven't done that, and it could get complicated, so I'll leave it at that and move on.

ALGORITHMS

A few notes on the following code:

- This is working JavaScript code from my game, with only minor edits for clarity.

- Porting to C, C++, and similar languages should be very direct in most cases.

- Polygons are stored as flat interleaved arrays:

- Points are stored as

One piece of advice: if you intend to roll your own physics engine, allow yourself a few weeks to iron it out (maybe months, if you have ambitious goals). Study the algorithms. Code up some demos. Step through the algorithms in a debugger. Get comfortable with the math. Look for other algorithms that might suit your particular needs better than these ones.

LINE VS LINE

This is the standard line intersection formula, useful for all kinds of things.

function lineIntersection(a,b, c,d){ var v,u,d,r,s; v= [ d[0]-c[0], d[1]-c[1] ]; u= [ b[0]-a[0], b[1]-a[1] ]; d= u[0]*v[1] - u[1]*v[0]; r= ((a[1]-c[1])*v[0] - (a[0]-c[0])*v[1]) / d; s= ((a[1]-c[1])*u[0] - (a[0]-c[0])*u[1]) / d; return (0<=r && r<=1 && 0<=s && s<=1); }

LINE VS POLYGON

To detect when a cannonball hits a ship, trace a line from the cannonball's previous position (last frame) to its current position, and test to see if it crossed any line segment of the ship's hull. This is reasonably accurate even if the framerate drops to a miserable 10 fps. This is an easy form of Continuous Collision Detection that's ideal for CCD's main use cases: bullets and beam weapons.

// a, b = [x,y] // edges = [x1,y1, x2,y2, ...] // function linePolygonIntersection(a,b, edges){ var c,d,i; d= edges.slice(edges.length-2); for(i=0; i<edges.length; i++){ c= d; d= edges.slice(i, i+2); if(lineIntersection(a,b,c,d)) return true; } return false; } // How to use it: for each cannonball for each ship if(linePolygonIntersection(cannonball.oldpos, cannonball.pos, ship.hullVerts) BOOM

POINT VS POLYGON

This is one of the many elegant "even-odd rule" algorithms. Starting at the point in question, trace a line in any direction and count how many times it crosses the polygon. If it's an odd number, the point is inside the polygon.

It's really just a Line-vs-Polygon algorithm, optimized for a horizontal line.

// point = [x,y] // polygon = [x1,y1, x2,y2, ...] // function isPointInPolygon(point, polygon){ var tx= point[0], ty=point[1]; var a; var b= polygon.slice(polygon.length-2); var yflag0= (b[1] >= ty); var yflag1; var inside_flag=false; for(var j=0; j<polygon.length; j+=2){ a= b; b= polygon.slice(j, j+2); yflag1= (b[1] >= ty); if(yflag0 != yflag1){ // vertexes straddle our point's Y-coord: potential hit. // LERP to check X-coord... if(yflag1 == ((b[1]-ty) * (a[0]-b[0]) >= (b[0]-tx) * (a[1]-b[1]))) inside_flag= !inside_flag; } yflag0= yflag1; } return inside_flag; }

POLYGON VS POLYGON

I'm using the Separating Axis Theorem (SAT) algorithm. It's not "the best" but it's relatively simple and fast enough for most purposes.

Right: if we project along every line of both ship polygons, we always get an overlap (no separation), so the ships have collided. The axis with the smallest overlap (red) provides the 'separation vector' which tells us exactly how far we must move one ship to resolve the collision state. We could also move each ship half that amount, or whatever makes sense.

Left: there is at least one separating axis, therefore the ships have not collided.

// Algorithm: // For each side, compute its "axis" (normal vector), B. Normalize it. // Project each vertex, of both polygons, onto the axis. // Note: dot(A,B) is an adequate approximation of proj(A,B) for this purpose. // Keep track of min and max values for each polygon. // If there's an axis with no overlap - a SEPARATING AXIS - we can stop. // If the min/max values overlap on all sides, we have a collision. // Use the smallest overlap to separate the shapes. // // P1,P2 = two *convex* polygons in [x,y, x,y, ...] format // cx = canvas context for debug drawing (optional) // function collidePolyPoly(P1, P2, cx){ var collision= true; var dmin= -Infinity; var collaxis= null; function _projectVerts(P1, P2, flip){ var A=P1, B=P2; var v1; var v2= A.slice(A.length-2); for(var i=0; i<A.length; i+=2){ v1= v2; v2= A.slice(i,i+2); var edge= vsub(v2,v1); var axis= vnormal(edge); var midpt= vlerp(v1, v2, 0.5); if(cx) drawLine(cx, midpt, vadd(midpt,axis), 'red', 0.05); var minA= vdot(axis, A.slice(0,2)); var maxA= minA; for(var j=2; j<A.length; j+=2){ var t= vdot(axis, A.slice(j,j+2)); if(t<minA) minA=t; else if(t>maxA) maxA=t; } var minB= vdot(axis, B.slice(0,2)); var maxB= minB; for(var j=2; j<B.length; j+=2){ var t= vdot(axis, B.slice(j,j+2)); if(t<minB) minB=t; else if(t>maxB) maxB=t; } var d0= minA-maxB; var d1= minB-maxA; var d= Math.max(d0,d1); if(d>0){ collision=false; return null; } if(d>dmin){ // keep track of minimum separation vector dmin=d; collaxis={ vector: axis, dist: d, midpt: midpt, poly: (d0<d1) ? P1 : P2, flip: (d0>d1)^flip, }; } } } _projectVerts(P1,P2,false); if(!collision) return null; _projectVerts(P2,P1,true); if(!collision) return null; var sep= vscale(collaxis.vector, collaxis.dist); if(cx) drawLine(cx, collaxis.midpt, vadd(collaxis.midpt, sep), 'red', 0.2); return { vector: sep, poly: collaxis.poly, flip: collaxis.flip, }; } // vector math helpers function vadd(a,b){ return [a[0]+b[0], a[1]+b[1]]; } function vsub(a,b){ return [a[0]-b[0], a[1]-b[1]]; } function vscale(v,s){ return [v[0]*s, v[1]*s]; } function vdot(a,b) { return a[0] * b[0] + a[1] * b[1]; } function vcross(a,b) { return a[0] * b[1] - a[1] * b[0]; } function vlerp(a, b, lerp) { return [a[0] + lerp * (b[0] - a[0]), a[1] + lerp * (b[1] - a[1])]; } function vnormal(v){ // return a perpendicular unit vector var mag = Math.sqrt(v[0]*v[0] + v[1]*v[1]); return [-v[1]/mag, v[0]/mag]; // normalize & perpendicularize } // debug-draw helper function drawLine(cx, p1, p2, color, weight){ if(!cx) return; cx.strokeStyle=color; cx.lineWidth=weight; cx.beginPath(); cx.moveTo(p1[0],p1[1]); cx.lineTo(p2[0],p2[1]); cx.stroke(); }

NOTES

I mentioned that I'm using similar algorithms in my current C++ sidescroller project. One advantage over physics libraries is that I'm able to use the same structure-of-arrays data for collisions, skeleton animation, and rendering.

These old algorithms are heavy on Structures-of-Arrays, a staple of modern CPU/GPU/cache optimization and Data Oriented Design. That was a happy accident. The code felt questionable to me at first, I hadn't heard of DOD, and optimization would have been premature at the time. It was only when I started porting to C++/OpenGL that I realized how well these algorithms fit today's hardware.

While backporting recent improvements from C++ to JS, I used interleaved arrays to overcome JS's lack of pointers. For example, by storing splines as flat arrays

SOURCES

Line Intersection formula from the comp.graphics.algorithms FAQ (section 1.03)

Point vs Polygon

ptinpoly.c (by Eric Haines, in Graphic Gems IV) (I used the last algorithm, "CrossingsMultiplyTest")

SAT visual explanation + Flash demo

SAT tutorial with formulas + VB code

Related article on gamedev.net: Shooting At Stuff (how to aim at moving targets)

ARTICLE UPDATE LOG

The approach uses a Proportional-Integral-Derivative (PID) control loop on the angle of the body to apply torque to it in a well controlled and predictable manner. Two examples of using the approach are shown in the video, one for a missile that can only move in the direction it is facing and another for an "Entity" like a game character that moves independent of its facing direction.

The video below shows this in action. You don't have to watch the whole video to get the idea, but there is a lot in there to see...

Whether your game is in 2-D or 3-D, you often have the need to make an object "turn" to face another direction. This could be a character's walking direction as they are moving, the direction they are shooting while sitting crouched, the direction a missile is flying in, the direction a car is racing towards, etc. This is the job of the "character controller", the piece of code in your system responsible for the basic "movement" operations that a character must undergo (seek, turn, arrive, etc.).

Building games that use physics engines is a lot of fun and adds a level of realism to the game that can dramatically improve the gameplay experience. Objects collide, break, spin, bounce, and move in more realistic ways.

Moving in more realistic ways is

Most important of all, this is a physics engine. You should be using it as such to make bodies move as expected.

Our goal is to create a solution to change the facing direction of the body by applying turning torque to it.

If we decouple the problem of "how it turns" from "how it moves", we can use the same turning solution for other types of moving bodies where the facing direction needs to be controlled. For this article, we are considering a missile that is moving towards its target.

Here, the missile is moving in a direction and has a given velocity. The angle of the velocity is measured relative to the x-axis. The "facing direction" of the missile is directly down the nose and the missile can only move forward. We want to turn it so that it is facing towards the target, which is at a different angle. For the missile to hit the target, it has to be aiming at it. Note that if we are talking about an object that is not moving, we can just as easily use the angle of the body relative to the x-axis as the angle of interest.

The basic idea behind a control system is to take the difference of "what you want the value to be" and "what the value is" and adjust your input to the system so that, over time, the system converges to your desired value.

From this wikipedia article:

A familiar example of a control loop is the action taken when adjusting hot and cold faucets (valves) to maintain the water at a desired temperature. This typically involves the mixing of two process streams, the hot and cold water. The person touches the water to sense or measure its temperature. Based on this feedback they perform a control action to adjust the hot and cold water valves until the process temperature stabilizes at the desired value.

There is a huge body of knowledge in controls system theory. Polynomials, poles, zeros, time domain, frequency domain, state space, etc. It can seem daunting to the uninitiated.

Let's start with the basic variable we want to "control", the difference between the angle we want to be facing and the angle of the body/velocity:

\(e(t) = desired - actual\)

Here, \(e(t)\) is the "error". We want to drive the error to 0. We apply forces to the body to make it turn in a direction to make it turn so that it moves \(e(t)\) towards 0. To do this, we create a function \(f(.)\), feed \(e(t)\) into it, and apply torque to the body based on it. Torque makes bodies turn:

\(torque(t) = I * f(e(t)), I \equiv Angular Inertia \)

The first and most obvious choice is to apply a torque that is proportional to the \(e(t)\) itself. When the error is large, large force is applied. When the error is small, small force is applied. Something like:

\(f(e(t)) = K_p * e(t)\)

And applying the torque with just proporational feedback will work. Somewhat. The problem is that when the error is small, the corrective torque is small as well. So as the body is turning \(e(t)\) gets small (nearing the intended angle), the retarding torque is small also. So the body overshoots and goes past the angle intended. Then you start to swing back and eventually oscillate to a steady state. If the value \(K_p\) is not too large, the oscillation should settle into "damped" (exponentially decaying) sinusoid, dropping a large amount each oscillation (stable solution). It may also spiral off towards infinity (unstable solution), or just oscillate forever around the target point (marginally stable solution). If you reduce \(K_p\) so that it is not moving so fast, then when \(e(t)\) is large, you don't have a lot of driving force to get moving.

A pure proportional error also has a bias (Steady State Error) that keeps the final output different from the input. The error is a function of \(K_p\) of the form:

\(Steady State Error = \frac{desired}{ [1 + constant * K_p]} \)

So increasing the \(K_p\) value makes the bias smaller (good). But this will also make it oscillate more (bad).

The next term to add is the integral term, the "I" in PID:

\(f(e(t)) = K_p * e(t) + \int\limits_{-\infty}^{now} K_i * e(t) \)

For each time step, if \(e(t)\) has a constant value, the integral term will work to counter it:

- If direction to the target suddenly changes a small amount, then over each time step, this difference will build up and create turning torque.
- If there is a bias in the direction (e.g. Steady State Error), this will accumulate over the time steps and be countered.

We don't have to calculate the actual integral. We probably don't want to anyway since it stretches back to \(-\infty\) and an error back in the far past should have little effect on our near term decisions.

We can estimate the integral over the short term by summing the value of \(e(t)\) over the last several cycles and multiplying by the time step (Euler Integration) or some other numerical technique. In the code base, the Composite Simpson's Rule technique was used.

Most PID controllers stop at the "PI" version. The proportional part gets the output swinging towards the input and the integral part knocks out the bias or any steady external forces that might be countering the proportional control.

\(f(e(t)) = K_p * e(t) + \int\limits_{-\infty}^{now} K_i * e(t) dt + K_d * \frac{de(t)}{dt}\)

Consider what happens when \(e(t)\) is oscillating. Its behavior is like a sine function. The derivative of this is a cosine function and its maximum occurs when sin(e(t)) = 0. That is to say, the derivative is largest when \(e(t)\) is swinging through the position we want to achieve. Conversely, when the oscillation is at the edge, about to change direction, its rate of change switches from positive to negative (or vice versa), so the derivative is smallest (minimum).

Just like the integral, the derivative can be estimated numerically. This is done by taking differences over the last several \(e(t)\) values (see the code).

Because we are software minded, whatever algorithm we want to use for a PID controller, we want to wrap it into a convenient package, give it a clean interface, and hide everything except what the user needs. This needs to be "owned" by the entity that is doing the turning.

The

The interface itself is generic so that the

The

/******************************************************************** * File : PIDController.h * Project: Interpolator * ******************************************************************** * Created on 10/13/13 By Nonlinear Ideas Inc. * Copyright (c) 2013 Nonlinear Ideas Inc. All rights reserved. ******************************************************************** * This software is provided 'as-is', without any express or implied * warranty. In no event will the authors be held liable for any * damages arising from the use of this software. * * Permission is granted to anyone to use this software for any * purpose, including commercial applications, and to alter it and * redistribute it freely, subject to the following restrictions: * * 1. The origin of this software must not be misrepresented; you must * not claim that you wrote the original software. If you use this * software in a product, an acknowledgment in the product * documentation would be appreciated but is not required. * 2. Altered source versions must be plainly marked as such, and * must not be misrepresented as being the original software. * 3. This notice may not be removed or altered from any source * distribution. */ #ifndef __Interpolator__PIDController__ #define __Interpolator__PIDController__ #include "CommonSTL.h" #include "MathUtilities.h" /* This class is used to model a Proportional- * Integral-Derivative (PID) Controller. This * is a mathemtical/control system approach * to driving the state of a measured value * towards an expected value. * */ class PIDController { private: double _dt; uint32 _maxHistory; double _kIntegral; double _kProportional; double _kDerivative; double _kPlant; vector<double> _errors; vector<double> _outputs; enum { MIN_SAMPLES = 3 }; /* Given two sample outputs and * the corresponding inputs, make * a linear pridiction a time step * into the future. */ double SingleStepPredictor( double x0, double y0, double x1, double y1, double dt) const { /* Given y0 = m*x0 + b * y1 = m*x1 + b * * Sovle for m, b * * => m = (y1-y0)/(x1-x0) * b = y1-m*x1 */ assert(!MathUtilities::IsNearZero(x1-x0)); double m = (y1-y0)/(x1-x0); double b = y1 - m*x1; double result = m*(x1 + dt) + b; return result; } /* This funciton is called whenever * a new input record is added. */ void CalculateNextOutput() { if(_errors.size() < MIN_SAMPLES) { // We need a certain number of samples // before we can do ANYTHING at all. _outputs.push_back(0.0); } else { // Estimate each part. size_t errorSize = _errors.size(); // Proportional double prop = _kProportional * _errors[errorSize-1]; // Integral - Use Extended Simpson's Rule double integral = 0; for(uint32 idx = 1; idx < errorSize-1; idx+=2) { integral += 4*_errors[idx]; } for(uint32 idx = 2; idx < errorSize-1; idx+=2) { integral += 2*_errors[idx]; } integral += _errors[0]; integral += _errors[errorSize-1]; integral /= (3*_dt); integral *= _kIntegral; // Derivative double deriv = _kDerivative * (_errors[errorSize-1]-_errors[errorSize-2]) / _dt; // Total P+I+D double result = _kPlant * (prop + integral + deriv); _outputs.push_back(result); } } public: void ResetHistory() { _errors.clear(); _outputs.clear(); } void ResetConstants() { _kIntegral = 0.0; _kDerivative = 0.0; _kProportional = 0.0; _kPlant = 1.0; } PIDController() : _dt(1.0/100), _maxHistory(7) { ResetConstants(); ResetHistory(); } void SetKIntegral(double kIntegral) { _kIntegral = kIntegral; } double GetKIntegral() { return _kIntegral; } void SetKProportional(double kProportional) { _kProportional = kProportional; } double GetKProportional() { return _kProportional; } void SetKDerivative(double kDerivative) { _kDerivative = kDerivative; } double GetKDerivative() { return _kDerivative; } void SetKPlant(double kPlant) { _kPlant = kPlant; } double GetKPlant() { return _kPlant; } void SetTimeStep(double dt) { _dt = dt; assert(_dt > 100*numeric_limits<double>::epsilon());} double GetTimeStep() { return _dt; } void SetMaxHistory(uint32 maxHistory) { _maxHistory = maxHistory; assert(_maxHistory >= MIN_SAMPLES); } uint32 GetMaxHistory() { return _maxHistory; } void AddSample(double error) { _errors.push_back(error); while(_errors.size() > _maxHistory) { // If we got too big, remove the history. // NOTE: This is not terribly efficient. We // could keep all this in a fixed size array // and then do the math using the offset from // the beginning and module math. But this // gets complicated fast. KISS. _errors.erase(_errors.begin()); } CalculateNextOutput(); } double GetLastError() { size_t es = _errors.size(); if(es == 0) return 0.0; return _errors[es-1]; } double GetLastOutput() { size_t os = _outputs.size(); if(os == 0) return 0.0; return _outputs[os-1]; } virtual ~PIDController() { } };

This is a very simple class to use. You set it up calling the

Looking at the

void ApplyTurnTorque() { Vec2 toTarget = GetTargetPos() - GetBody()->GetPosition(); float32 angleBodyRads = MathUtilities::AdjustAngle(GetBody()->GetAngle()); if(GetBody()->GetLinearVelocity().LengthSquared() > 0) { // Body is moving Vec2 vel = GetBody()->GetLinearVelocity(); angleBodyRads = MathUtilities::AdjustAngle(atan2f(vel.y,vel.x)); } float32 angleTargetRads = MathUtilities::AdjustAngle(atan2f(toTarget.y, toTarget.x)); float32 angleError = MathUtilities::AdjustAngle(angleBodyRads - angleTargetRads); _turnController.AddSample(angleError); // Negative Feedback float32 angAcc = -_turnController.GetLastOutput(); // This is as much turn acceleration as this // "motor" can generate. if(angAcc > GetMaxAngularAcceleration()) angAcc = GetMaxAngularAcceleration(); if(angAcc < -GetMaxAngularAcceleration()) angAcc = -GetMaxAngularAcceleration(); float32 torque = angAcc * GetBody()->GetInertia(); GetBody()->ApplyTorque(torque); }

If you look carefully at the video, there is a distinct difference in the way path following works for the missile vs. the character (called the

The

I have also, quite deliberately, left out a bit of key information on how to tune the constants for the PID controller. There are numerous articles on Google for how to tune a PID control loop, and I have to leave something for you to do, after all.

You will also note that the default value for

The source code for this, written in cocos2d-x/C++, can be found on github here. The

There are two cases for "shooting" at something:

- Non-Rotating Shooter
- Rotating Shooter

This article only covers the case of the non-rotating shooter. The second case is a bit more complicated. I've posted a second article "Throwing a Winning Pass" to cover that.

Your goal is to predict the future position of the target so that you can launch something at it and have them collide. It is given that you know the following:

- The position of the shooter when the projectile will be launched: \(\vec{P_s}\)
- The position of the target when the shooter will launch the projectile (i.e. now): \(\vec{P_T^0}\)
- The speed at which your projectiles travel: \(S_b\)
- The velocity of the target, \(\vec{v_T}\)

What you would like to find is \(\vec{P_T^1}\), the final position where the projectile will intercept the target.

We are going to have to assume that the target will travel at a constant velocity while the projectile is moving. Without this, the projectile will have to adapt after it has been launched (because it cannot predict the future position of the target if it is changing velocity in an unpredictable way). It should be pointed out that this is not limiting and in fact, still works out pretty well as long as the "edges" are handled well (target right in front of you or headed straight at you).

Consider the image below:

This image shows the situation for the shooter. The target is moving with some constant velocity and the intercept position for the projectile, when launched, is unknown until the calculation is done. When the projectile is launched, it will intercept the target after \(t_B\) seconds. The projectile is launched with a constant speed. We don't know its direction yet, but we know the magnitude of its velocity, its speed will be \(S_b\).

If the target is at position \(\vec{P_T^0}\) now, it will be at position \(\vec{P_T^1}\), given by:

(1) \(\vec{P_T^1} = \vec{P_T^0} + \vec{v_T} * t_B\)

Since the projectile will have traveled for the same amount of time, it will have moved from \(\vec{P_s}\) to \(\vec{P_T^1}\) as well. In that time, it will have moved a distance of \(S_B x t_B\). Since we are talking about vector quantities here, we can write this as:

\(\mid\vec{P_T^1}-\vec{P_s}\mid = S_b * t_B\)

If we square both sides and break it into components to get rid of the absolute value:

(2) \((P_{Tx}^1 - P_{Sx})^2 +(P_{Ty}^1 - P_{Sy})^2 = S_b^2 * t_B^2\)

Breaking (1) into components as well and substituting back into (2) for the value of \(P_{Tx}^1\) and \(P_{Ty}^1\), we get the following:

\((P_{T0x} - P_{Sx} + v_{Tx}t_B)^2 + (P_{T0y} - P_{Sy} + v_{Ty}t_B)^2 = S_b^2 * t_B^2\)

For the sake of simplicity, we going to redefine:

\(P_T^0 - P_s = R \)(this is a constant)

After some algebra, this gives us the **final** equation:

\(t_B^2(v_{Tx}^2 + v_{Ty}^2-S_B^2) + t_B(2*R_x*v_{Tx} + 2*R_y*v_{Ty}) + (R_x^2 + R_y^2) = 0\)

This is a quadratic in \(t_B\):

\(t_b = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}\)

where:

\( a =v_{Tx}^2 + v_y^2-S_B^2\)

\( b =2(R_x*v_{Tx} + R_y*v_{Ty})\)

\( c = R_x^2 + R_y^2\)

You can test the discriminant, \(b^2-4ac\):

< 0 \(\Rightarrow\) No Solution.

= 0 \(\Rightarrow\) One solution.

> 0 \(\Rightarrow\) Two solutions, pick the lowest positive value of \(t_B\).

Once you have solved the quadratic for \(t_B\), you can then substitute it back into (1) and calculate the intercept position, \(\vec{P_T^1}\).

Putting this together and covering some edge cases:

/* Calculate the future position of a moving target so that * a projectile launched immediately can intercept (collide) * with it. * * Some situations where this might be useful for an AI to * make this calculation. * * 1. Shooting a projectile at a moving target. * 2. Launching a football or soccer ball to a player. * 3. Figuring out the best position to jump towards in * a platform game. * * * The output value, solution, is the position that the * intercept will occur at and the location that the * projectile should be launched towards. * * The function will return false if a solution cannot * be found. Consider the case of a target moving away * from the shooter faster than the speed of the * projectile and you will see at least one case where * this calculation may fail. */ bool CalculateInterceptShotPosition(const Vec2& pShooter, const Vec2& pTarget0, const Vec2& vTarget, float64 sProjectile, Vec2& solution ) { // This formulation uses the quadratic equation to solve // the intercept position. Vec2 R = pTarget0 - pShooter; float64 a = vTarget.x*vTarget.x + vTarget.y*vTarget.y - sProjectile*sProjectile; float64 b = 2*(R.x*vTarget.x + R.y*vTarget.y); float64 c = R.x*R.x + R.y*R.y; float64 tBullet = 0; // If the target and the shooter have already collided, don't bother. if(R.LengthSquared() < 2*DBL_MIN) { return false; } // If the squared velocity of the target and the bullet are the same, the equation // collapses to tBullet*b = -c. If they are REALLY close to each other (float tol), // you could get some weirdness here. Do some "is it close" checking? if(fabs(a) < 2*DBL_MIN) { // If the b value is 0, we can't get a solution. if(fabs(b) < 2*DBL_MIN) { return false; } tBullet = -c/b; } else { // Calculate the discriminant to figure out how many solutions there are. float64 discriminant = b*b - 4 * a * c; if(discriminant < 0) { // All solutions are complex. return false; } if (discriminant > 0) { // Two solutions. Pick the smaller one. // Calculate the quadratic. float64 quad = sqrt(discriminant); float64 tBullet1 = (-b + quad)/(2*a); float64 tBullet2 = (-b - quad)/(2*a); if((tBullet1 < 0) && (tBullet2 < 0)) { // This would be really odd. // Both times are negative. return false; } else if(tBullet2 < 0 && tBullet1 >= 0) { // One negative, one positive. tBullet = tBullet1; } else if(tBullet1 < 0 && tBullet2 >= 0) { // One negative, one positive. tBullet = tBullet2; } else if(tBullet1 < tBullet2) { // First less than second tBullet = tBullet1; } else { // Only choice left tBullet = tBullet2; } } else { tBullet = -b / (2*a); } } // If the time is negative, we can't get there from here. if(tBullet < 0) { return false; } // Calculate the intercept position. solution = pTarget0 + tBullet*vTarget; return true; }

I have posted a working solution, with a simulation of using the above function and which you can tinker with, on github.

It is assumed that the reader has insight in theory behind inertia tensors and/or has an existing physics system that can make use of such a matrix. Knowledge in multivariate calculus is also assumed.

Firstly, what is meant by saying capsule? Well, simply put, it is a cylinder that (instead of flat ends) has hemispherical ends. By cylinder I mean a circular one (not elliptical). Hemispheres have the same radius. It should also be mentioned that the given body is solid (not hollow).

Where

**Figure 2.** Capsule decomposition.

The approach that we will be using to calculate the inertia tensor can be seen from Figure 2. The body will be decomposed into three parts: upper hemisphere, lower hemisphere and cylinder. Inertia tensors will be calculated for each of them and then just summed together to result in a complete capsule tensor. Here are given, in Cartesian coordinates, the well-known equations for moments and products of inertia needed for the inertia tensor in general.

(**eq. 1)**

(**eq. 2)**

Note that in eq. 2 products of inertia are equal to zero. There exists a theorem that says that it is possible to find three (mutually perpendicular) axes for every body so the tensor produced from those axes has zero products of inertia. Those axes are called the principal axes. Furthermore, for bodies with constant density an axis of rotational symmetry is a principal axis. We have three perpendicular axes of rotational symmetry for our capsule. First one is the y axis of continuous rotational symmetry and then there are x and z axes of discrete rotational symmetry. I will not go into detail about what axes of rotational symmetry are, but I recommend to go and check it out.

(**eq. 3)**

Since we will integrate over space eq. 3 will be used. If density is required to be some other constant than 1, the end result can simply be multiplied by that constant.

Before we dive into integration here is the Cartesian coordiante system we will be using for cylinder:

**Figure 3.** Cylinder and its position relative to coordinate system.

Let's start with easiest to compute moment of inertia, the one around the y axis. As previously mentioned, integrand is just a squared distance from the axis. With that in mind we can leave the Cartesian coordinate system and use the cylindrical one. The y axis remains unchanged. Also note that solutions will be written immediately without showing any steps.

(**eq. 4)**

First we add up (integrate) over 2*pi*r interval to get the moment of inertia of a circular curve around the y axis. After that we integrate that curve over the interval R to get the moment of inertia of a disc surface perpendicular to the y axis. Lastly those discs are added up across [-H/2, H/2] interval to get the moment of inertia of a whole cylinder. The integral for the mass of a cylinder is very similar and follows the same intuition. It is given by:

(**eq. 5)**

Now, slightly more complicated integrals are the ones for moments of inertia around x and z axis. Those are equal, so only one calculation is needed. Note that we are back in Cartesian coordinates for this one and the integral is set up like we are calculating

(**eq. 6)**

This is it for the cylinder as I will not go into detail about intuition for this integral since there is more work to be done for hemisphere caps.

The following is how the hemisphere is defined relative to Cartesian coordinate system for the purpose of integration:

**Figure 4.** Hemisphere and its position relative to coordinate system.

The easiest to calculate moment of inertia of hemispheres is also the one around y axis.

(**eq. 7)**

We are in the cylindrical coordinate system and the intuition is similar to eq. 4. The only difference is that now the radius of a disc we need to integrate is a variable of the y position of a disc (sqrt(R^2 - y^2)). Then we just integrate those discs from 0 to R and we end up with the whole hemisphere region. As for the cylinder, the integral for the mass of a hemisphere follows the same intuition and is given by:

(**eq. 8)**

Now, when we try to calculate the moments of inertia around x and z axis (which are equal), things get tricky. Note that integral is set up like we are calculating

(**eq. 9)**

**Figure 5.** Center of mass of a hemisphere.

As it can be seen in figure 5, center of mass is not in the origin of the coordinate system, but the axis for which we calculated the moment of inertia is. In order to work with this moment of inertia (to apply Steiner's rule) we need it to be in the center of mass. This is where we calculate the distance

(**eq. 10)**

(**eq. 11)**

Where

(**eq. 12)**

Moments of inertia for the hemisphere cap on the bottom of the capsule do not differ to those here.

(**eq. 13)**

Equations in eq. 13 give us separated masses of cylinder and hemisphere as well as the mass of the whole body. Notice how mass of the hemisphere (

(**eq. 14)**

Although previously shown formulas for moments of inertia are not given by

The following code implementation is given in C-like language and is fairly optimized.

#define PI 3.141592654f #define PI_TIMES2 6.283185307f const float oneDiv3 = (float)(1.0 / 3.0); const float oneDiv8 = (float)(1.0 / 8.0); const float oneDiv12 = (float)(1.0 / 12.0); void ComputeRigidBodyProperties_Capsule(float capsuleHeight, float capsuleRadius, float density, float &mass, float3 ¢erOfMass, float3x3 &inertia) { float cM; // cylinder mass float hsM; // mass of hemispheres float rSq = capsuleRadius*capsuleRadius; cM = PI*capsuleHeight*rSq*density; hsM = PI_TIMES2*oneDiv3*rSq*capsuleRadius*density; // from cylinder inertia._22 = rSq*cM*0.5f; inertia._11 = inertia._33 = inertia._22*0.5f + cM*capsuleHeight*capsuleHeight*oneDiv12; // from hemispheres float temp0 = hsM*2.0f*rSq / 5.0f; inertia._22 += temp0 * 2.0f; float temp1 = capsuleHeight*0.5f; float temp2 = temp0 + hsM*(temp1*temp1 + 3.0f*oneDiv8*capsuleHeight*capsuleRadius); inertia._11 += temp2 * 2.0f; inertia._33 += temp2 * 2.0f; inertia._12 = inertia._13 = inertia._21 = inertia._23 = inertia._31 = inertia._32 = 0.0f; mass = cM + hsM * 2.0f; centerOfMass = {0.0f, 0.0f, 0.0f}; }

It would seem that

Capsules defined in a way which was presented here are very commonly used in games and physical simulations, so this tensor should be of use to anyone building their own 3D game engine or similar software.

Although, right-handed Cartesian coordinate system was used in this article, you can also use the code in left-handed system as long as y axis is the longitudinal axis (the axis that goes through caps).

The somewhat tedious processes of integration were skipped and only the end results were shown. This is mostly because it would take too much space. If you are curious whether those expressions are correct, you can try them on Wolfram or try to integrate them yourself.

Also, you can follow me on Twitter.

The forums are replete with people trying to learn how to efficiently find inverses of matrices or solve the classic \(Ax=b\) problem. Here's a decent method that is fairly easy to learn and implement. Hopefully it might also serve as a stepping stone to learning some of the more advanced matrix factorization methods, like Cholesky, QR, or SVD.

In 1948, Alan Turing came up with LU decomposition, a way to factor a matrix and solve \(Ax=b\) with numerical stability. Although there are many different schemes to factor matrices, LU decomposition is one of the more commonly-used algorithms. Interestingly enough, Gauss elimination can be implemented as LU decomposition. The computational effort expended is about the same as well.

So why would anyone want to use this method? First of all, the time-consuming elimination step can be formulated so that only the entries of \(A\) are required. As well, if there are several matrices of \(b\) that need to be computed for \(Ax=b\), then we only need to do the decomposition once. This last benefit helps us compute the inverse of \(A\) very quickly.

Let's start with an \(Ax=b\) problem, where you have an \(n\)x\(n\) matrix \(A\) and \(n\)x\(1\) column vector \(b\) and you're trying to solve for a \(n\)x\(1\) column vector \(x\). The idea is that we break up \(A\) into two matrices: \(L\), which is a lower triangular matrix, and \(U\), which is an upper triangular matrix. It would look something like this:

\[

\left [

\begin{matrix}

a_{11} & a_{12} & a_{13} \\

a_{21} & a_{22} & a_{23} \\

a_{31} & a_{32} & a_{33} \\

\end{matrix} \right ] = \left [

\begin{matrix}

l_{11} & 0 & 0 \\

l_{21} & l_{22} & 0 \\

l_{31} & l_{32} & l_{33} \\

\end{matrix} \right ] \left [

\begin{matrix}

u_{11} & u_{12} & u_{13} \\

0 & u_{22} & u_{23} \\

0 & 0 & u_{33} \\

\end{matrix} \right ]

\]

This seems like a lot more work, but allows us to use some cool math tricks. Our original equation, substituted with our decomposition is now \((LU)x=b\):

\[

\left (

\left [ \begin{matrix}

l_{11} & 0 & 0 \\

l_{21} & l_{22} & 0 \\

l_{31} & l_{32} & l_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

u_{11} & u_{12} & u_{13} \\

0 & u_{22} & u_{23} \\

0 & 0 & u_{33} \\

\end{matrix} \right ]

\right )

\left [ \begin{matrix}

x_1 \\

x_2 \\

x_3 \\

\end{matrix} \right ] =

\left [ \begin{matrix}

b_1 \\

b_2 \\

b_3 \\

\end{matrix} \right ]

\]

If we shift the parentheses, we get the following:

\[

\left [ \begin{matrix}

l_{11} & 0 & 0 \\

l_{21} & l_{22} & 0 \\

l_{31} & l_{32} & l_{33} \\

\end{matrix} \right ]

\left (

\left [ \begin{matrix}

u_{11} & u_{12} & u_{13} \\

0 & u_{22} & u_{23} \\

0 & 0 & u_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

x_1 \\

x_2 \\

x_3 \\

\end{matrix} \right ]

\right )

=

\left [ \begin{matrix}

b_1 \\

b_2 \\

b_3 \\

\end{matrix} \right ]

\]

Looking just inside the parentheses, we can see another \(Ax=b\) type problem. If we say that \(Ux=d\), where \(d\) is a different column vector from \(b\), we have 2 separate \(Ax=b\) type problems:

\[

\left [ \begin{matrix}

l_{11} & 0 & 0 \\

l_{21} & l_{22} & 0 \\

l_{31} & l_{32} & l_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

d_1 \\

d_2 \\

d_3 \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

b_1 \\

b_2 \\

b_3 \\

\end{matrix} \right ]

\\

\left [ \begin{matrix}

u_{11} & u_{12} & u_{13} \\

0 & u_{22} & u_{23} \\

0 & 0 & u_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

x_1 \\

x_2 \\

x_3 \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

d_1 \\

d_2 \\

d_3 \\

\end{matrix} \right ]

\]

It looks like we just created more work for ourselves, but we actually made it easier. If we look at the problems, the matrices on the left are in a form like row echelon form. Using forward and back substitutions, we can solve easily for all the elements of the \(x\) and \(d\) matrices. We first solve \(Ld=b\) for \(d\), then we substitute into \(Ux=d\) to solve for \(x\). The real trick to this method is the decomposition of \(A\).

There are a couple catches to this method:

- The matrix \(A\) must be square to use LU factorization. Other factorization schemes will be necessary if \(A\) is rectangular.
- We have to be sure that \(A\) is a nonsingular (i.e. invertible) matrix. If it can't be inverted, then the decomposition will produce an \(L\) or \(U\) that is singular and the method will fail because there is no unique solution.
- To avoid division by zero or by really small numbers, we have to implement a pivoting scheme just like with Gaussian elimination. This makes the problem take the form \(PA=LU\), where P is a permutation matrix that allows us to swap the rows of A. P is usually the identity matrix with rows swapped such that \(PA\) produces the \(A\) matrix with the same rows swapped as P. Then the \(Ax=b\) problem takes the form \(LUx=Pb\) since \(PA=LU\).
- Some of the entries in the \(L\) and \(U\) matrices must be known before the decomposition, or else the system has too many unknowns and not enough equations to solve for all the entries of both matrices. For what's formally known as Doolittle decomposition, the diagonal entries of the \(L\) matrix are all 1. If we use Crout decomposition, the diagonals of the \(U\) matrix are all 1.

For a square matrix \(A\) with entries \(a_{ij},\,i=1,\cdots,n,\,j=1,\cdots,n\), the Crout decomposition is as follows:

First:

\[

l_{i1} = a_{i1},\,\textrm{for}\,i=1,2,\cdots,n\\

u_{1j} = \frac{a_{1j}}{l_{11}},\,\textrm{for}\,j=2,3,\cdots,n

\]

For \(j=2,3,\cdots,n-1\):

\[

l_{ij} = a_{ij}-\sum_{k=1}^{j-1}l_{ik}u_{kj},\,\textrm{for}\,i=j,j+1,\cdots,n \\

u_{jk} = \frac{a_{jk}-\sum_{i=1}^{j-1}l_{ji}u_{ik}}{l_{jj}},\,\textrm{for}\,k=j+1,j+2,\cdots,n

\]

Finally:

\[

l_{nn}=a_{nn}-\sum_{k=1}^{n-1}l_{nk}u_{kn}

\]

That's it! If you notice, \(A\) is being traversed in a specific way to build \(L\) and \(U\). To start building \(L\), we start with \(a_{11}\) and then traverse down the rows in the same column until we hit \(a_{n1}\). To start building \(U\), we start at \(a_{12}\) and traverse along the same row until we hit \(a_{1n}\). We then build \(L\) further by starting at \(a_{22}\) and traverse along the column until we hit \(a_{n2}\), then we build \(U\) further by starting at \(a_{23}\) and traverse along the row until \(a_{2n}\), and so on.

Since the 1's on the diagonal of \(U\) are given as well as the 0's in both the \(L\) and \(U\) matrices, there's no need to store anything else. In fact, there's a storage-saving scheme where all the calculated entries of both matrices can be stored in a single matrix as large as \(A\).

From here, \(d\) can be solved very easily with forward-substitution:

\[

d_1 = \frac{b_1}{l_{11}} \\

d_i = \frac{b_i - \sum_{j=1}^{i-1}l_{ij}d_j}{l_{ii}},\,\textrm{for}\,i=2,3,\cdots,n

\]

As well, \(x\) can be solved very easily with back-substitution:

\[

x_n = d_n \\

x_i = d_i - \sum_{j=i+1}^n u_{ij}x_j,\,\textrm{for}\,i=n-1,n-2,\cdots,1

\]

Let's try to solve the following problem using LU decomposition:

\[

\left [ \begin{matrix}

3 & -0.1 & -0.2 \\

0.1 & 7 & -0.3 \\

0.3 & -0.2 & 10 \\

\end{matrix} \right ]

\left [ \begin{matrix}

x_1 \\ x_2 \\ x_3 \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

7.85 \\ -19.3 \\ 71.4 \\

\end{matrix} \right ]

\]

First, we copy the first column to the \(L\) matrix and the scaled first row except the first element to the \(U\) matrix:

\[

L =

\left [ \begin{matrix}

3 & 0 & 0 \\

0.1 & - & 0 \\

0.3 & - & - \\

\end{matrix} \right ],\,\,\,

U =

\left [ \begin{matrix}

1 & -0.0333 & -0.0667 \\

0 & 1 & - \\

0 & 0 & 1 \\

\end{matrix} \right ]

\]

Then we compute the following entries:

\[

l_{22} = a_{22} - l_{21}u_{12} = 7-(0.1)(-0.0333) = 7.00333 \\

l_{32} = a_{32} - l_{31}u_{12} = -0.2-(0.3)(-0.0333) = -0.19 \\

u_{23} = \frac{a_{23}-l_{21}u_{13}}{l_{22}} = \frac{-0.3-(0.1)(-0.0667)}{7} = -0.0419 \\

l_{33} = a_{33}-l_{31}u_{13}-l_{32}u_{23} = 10-(0.3)(-0.0667)-(-0.19)(-0.0419) = 10.012 \\

\]

Inserting them into our matrices:

\[

L =

\left [ \begin{matrix}

3 & 0 & 0 \\

0.1 & 7.00333 & 0 \\

0.3 & -0.19 & 10.012 \\

\end{matrix} \right ],\,\,\,

U =

\left [ \begin{matrix}

1 & -0.0333 & -0.0667 \\

0 & 1 & -0.0419 \\

0 & 0 & 1 \\

\end{matrix} \right ]

\]

This is the LU factorization of the matrix. You can check if \(LU=A\). Now, we use back- and forward-substitution to solve the problem:

\[

d_1 = \frac{b_1}{l_{11}} = 7.85/3 = 2.6167 \\

d_2 = \frac{b_2-l_{21}d_1}{l_{22}} = (-19.3-(0.1)(2.6167))/7.00333 = -2.7932 \\

d_3 = \frac{b_3-l_{31}d_1-l_{32}d_2}{l_{33}} = (71.4-(0.3)(2.6167)-(-0.19)(-2.7932))/10.012 = 7 \\

\]

\[

x_3 = d_3 = 7 \\

x_2 = d_2 - u_{23}x_3 = -2.7932-(-0.0419)(7) = -2.5 \\

x_1 = d_1 - u_{12}x_2 - u_{13}x_3 = 2.6167-(-0.0333)(-2.5)-(-0.0667)(7) = 3 \\

\]

So the solution to the problem is:

\[

\left [ \begin{matrix}

x_1 \\ x_2 \\ x_3 \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

3 \\ -2.5 \\ 7 \\

\end{matrix} \right ]

\]

This solution can easily be verified by plugging it back into \(Ax=b\).

Here's some quick MATLAB code for LU decomposition:

function [L,U] = lucrout(A) [~,n] = size(A); L = zeros(n,n); U = eye(n,n); L(1,1) = A(1,1); for j=2:n L(j,1) = A(j,1); U(1,j) = A(1,j) / L(1,1); end for j=2:n-1 for i=j:n L(i,j) = A(i,j); for k=1:j-1 L(i,j) = L(i,j) - L(i,k)*U(k,j); end end for k=j+1:n U(j,k) = A(j,k); for i=1:j-1 U(j,k) = U(j,k) - L(j,i)*U(i,k); end U(j,k) = U(j,k) / L(j,j); end end L(n,n) = A(n,n); for k=1:n-1 L(n,n) = L(n,n) - L(n,k)*U(k,n); end end

LU decomposition is nice for solving a series of \(Ax=b\) problems with the same \(A\) matrix and different \(b\) matrices. This is advantageous for computing the inverse of \(A\) because only one decomposition is required. The definition of the inverse of a matrix \(A^{-1}\) is a matrix such that \(AA^{-1}=I\), where \(I\) is the identity matrix. This looks like this for a general 3x3 matrix:

\[

\left [ \begin{matrix}

A_{11} & A_{12} & A_{13} \\

A_{21} & A_{22} & A_{23} \\

A_{31} & A_{32} & A_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

a_{11}^\prime & a_{12}^\prime & a_{13}^\prime \\

a_{21}^\prime & a_{22}^\prime & a_{23}^\prime \\

a_{31}^\prime & a_{32}^\prime & a_{33}^\prime \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

1 & 0 & 0 \\

0 & 1 & 0 \\

0 & 0 & 1 \\

\end{matrix} \right ]

\]

This can be set up like \(n\) \(Ax=b\) problems with different \(b\) matrices and the \(x\) matrix becomes the nth column in the inverse matrix. For example, the first problem is:

\[

\left [ \begin{matrix}

A_{11} & A_{12} & A_{13} \\

A_{21} & A_{22} & A_{23} \\

A_{31} & A_{32} & A_{33} \\

\end{matrix} \right ]

\left [ \begin{matrix}

a_{11}^\prime \\

a_{21}^\prime \\

a_{31}^\prime \\

\end{matrix} \right ]

=

\left [ \begin{matrix}

1 \\

0 \\

0 \\

\end{matrix} \right ]

\]

The second column of the inverse can be computed by changing \(b\) to \([0,1,0]^T\), the third column with \([0,0,1]^T\), and so on. This method is quick because only back- and forward-substitution is required to solve for the column vectors after the initial LU decomposition.

There are a lot of other matrix factorization schemes besides LU, like Cholesky or QR factorization, but the general idea of decomposing a matrix into other matrices is roughly the same. The real key to computational savings comes from knowing beforehand what kind of matrix you're factoring and choosing the appropriate algorithm. For example, in structural finite element analysis, the matrix being decomposed is always symmetric positive definite. Cholesky decomposition is way more efficient and quicker than LU for those kinds of matrices, so it's preferred.

LU decomposition is a great tool for anyone working with matrices. Hopefully you can make use of this simple, yet powerful method.

Making a Game Engine: Core Design Principles

Moving objects around in a scene is not hard to do. You simply draw objects at different locations each frame. Adding rotation, and resizing on top of that makes the code a little trickier but still not hard to manage. But then start trying to move objects relative to another such as putting a sword in a character's hand and things can get really tricky. The sword should stay in the character's hand regardless of what direction they rotate or where they move. This is where using matrices to represent these transformations is extremely valuable.

In order for a coordinate to have any meaning there needs to be a point of reference to base that coordinate off of. Imagine you have a dot on a piece of paper. The dot does not have a coordinate unless you come up with a way to measure its location. One way would be to determine the distance the point is from the left side of the page and its distance from the bottom. These two measurements give the coordinate meaning since it has a point of reference, the bottom left corner of the page. The point of reference is known as a basis and any point measured relative to the same basis are said to be in the same coordinate system. If a worm were on a ball you could devise a way to measure the location of the worm on the ball. This coordinate would be the worm's location in the ball's coordinate space.

The ball itself also has a location but its location is measured in world coordinates. Let's assume that the location of the ball is measured from the bottom left corner of the image to the center of the ball. If you were to throw that ball the worm would move with the ball. The position of the worm in world coordinates would move as the ball moves, even though its position in ball coordinates stay the same. Transformations give us a way to convert between these coordinate systems relatively easy using matrices allowing us to easily keep the worm on the ball as it travels through the air.

Matrices have a lot of useful properties that lend themselves to solving many complex problems. We will be using them to transform objects on the screen. This article will focus on three simple transformations: translation, rotation, and scaling. I will be using 2D transformations for the sake of simplicity. The principles in the article also apply to 3D transformations. If you want to do 3D transformations, search around on the web for 3D transformation matrices. The jump from 2D to 3D is trivial with the exception of rotations. I am assuming you already know what matrices are and how to multiply them, what the inverse and transpose of a matrix are, and other basic matrix operations. (For help with matrices, try this video series)

So far this talk of coordinate spaces and matrices may seem to be complicated things when all you need to do is move or rotate objects. How exactly do these things help move objects on the screen anyway? Well, imagine you define a character sprite with a width of 50 and a height of 100. You could define four points that represent the corners of the sprite image (-25, 0), (-25, 0), (25, 100), (25, 0). These points represent the corners of the sprite relative to the character. So what whappens when the character moves? The transform of the player converts the points of the sprite to world space. Then the camera transform converts from world space to screen space to be drawn. If the player clicks on the screen you can use the inverse of these transforms to determine where on the character they clicked.

A translation simply changes an object's location. The matrix used to represent this is pretty simple too.

T

Rotates points around the origin.

The above matrix will rotate counter clockwise by θ radians.

Scaling an object changes its size. Scalar values in the range (0, 1) make objects smaller. A scale of 1 keeps it the same size, a scale greater than 1 makes it bigger. Really all the scale is doing is taking the current size and multiplying by the scale value.

Scales along the x and y directions by S

Now that we have a few different transformation matrices we need a way to transform points using them. To do this you simply treat the point as a 1x3 matrix padding the extra slots with a 1 and multiplying it with the transform matrix.

The result of that multiplication is a another 1x3 matrix. Just taking the top two elements of the resulting matrix gives a new x, y pair of the point in the new coordinate system.

Transforming a direction is slightly different than transforming a point. Take a look at the example below.

The ball has been squashed. This is achieved by scaling the y direction down while keeping the x the same. If directions were transformed using the same matrix that were used to transform points, then you would get the red arrow as a result, but notice how the red arrow is no longer perpendicular to the ball. We want the green arrow that still is. To get that result, use the following equation.

As shown above, to get the proper matrix for transforming direction you need the inverse transpose of the transform matrix. You can avoid transposing a 3x3 matrix if you simply do a row vector multiplication in front of the matrix. Also take note that the third element is a 0 instead of a 1. This makes transforming a direction uneffected by any translations. This is proper behavior for transforming a direction since directions still point the same way regaurdless of where they are moved to.

The three basic tranformations listed above aren't that useful on their own. To make things interesting we want to be able to join them together allowing objects to be translated, rotated, and scaled simultaneously. To do this we simply multiply the matrices together.

The resulting matrix,

One thing to note is that roatations pivot around the origin. This means if you move an object up ten units then rotate it clockwise 90 degrees it will end up 10 units to the right. If you want the object to pivot around its own center, rotate it first then move it.

Since matrices can be multiplied to concatenate transforms you can easily keep an object at the same position relative to another, like the example of the worm on a ball described above. You first create a transform using the position and orientation of the worm on the ball. This transform represents the conversion of the worm's local coordinates to the ball coordinates. You then create the transform of the ball relative to world coordinates. Since the worm is on the ball, to find the transform from worm coordinates to world coordinates you simply multiply them together.

Matrix insectToWorldTransform = ballToWorldTransform * insectToBallTransform;

Before moving into the next section I should mention the difference between column vectors and row vectors. I have been using column vectors in this article. This means vectors are represented as a 3x1 matrix that is multiplied on the right of the transform.

The alternative are row vectors. A row vector is a 1x3 matrix and is multiplied on the left of the matrix.

Mathematically, these two methods are the same. Either one works. Some prefer how compact column vectors are to write out on paper. Others find row vectors to be more intiuative. The important thing is to pick one and be consistent with it. However, the graphics package you are using may require matrices for row vectors or column vectors. If you are using the other method you will need to convert your matrices before passing them in. To do so, simply transpose the matrix. As an example, the translate matrix using column vectors is

The same matrix for row vectors is

Also, keep in mind that when using row vectors the order to multiply matrices in is reversed.

So using row vectors, the left most matrix is applied first to the point and the matrix on the far right is applied last. An easy way to remember the difference is the matrix closest to the vector is applied first. This means the order is left to right using row matrices and right to left using column matrices.

Now that we have covered some basics of how transforms work let's build a transform class. Since we don't want to deal with matrices at a high level, we will allow the transform to be manipulated using a position, rotation, and scale. The transform class will also have children to allow for a heirarchy. When the parent transform moves all of the children move with it and the children's position, rotation, and scale will be defined relative to its parent. Finally the transform will extend the component class that is attached to a game object. Game objects and components are touched upon in this article.

class Transform extends Component { // the parent transform of this transform // if it is null then the parent transform // is the world coordinate system private Transform parent; // all of the transforms that have this // transform set as their parent private Transform[] children; // the position relative to the parent transform private Vector2 localPosition = new Vector2(0.0f, 0.0f); // rotation relative to the parent private float localRotation = 0.0f; // scale relative to the parent private Vector2 localScale = new Vector2(1.0f, 1.0f); // specifies if the localToWorldTransform // needs to be recalulated private bool isDirty = false; // the transform that converts local coordinates // to world coordinates private Matrix localToWorldMatrix = Matrix.identity; // specifies if the worldToLocalMatrix // needs to be recalculated private bool isInverseDirty = false; // the transform that converts world cooridnates // to local coordinates private Matrix worldToLocalMatrix = Matrix.identity; /* * Whenever any change happens that changes the localToWorldMatrix * this should be called. That way the next time localToWorldMatrix * is requested it will be recalculated */ private void setDirty() { // only update dirty boolean if it isn't already dirty if (!isDirty) { isDirty = true; isInverseDirty = true; // set all children to be dirty since any modification // of a parent transform also effects its children's // localToWorldTransform foreach (Transform child in children) { child.setDirty(); } } } // change the parent transform. // setting it to null makes the // transform a child of world coordinates public void setParent(Transform value) { // remove this from the previous parent if (parent != null) { parent.children.remove(this); } // assign new parent parent = value; // add this to new parent if (parent) { parent.children.add(this); } // changes parents effects // the world position setDirty(); } public Transform getParent() { return parent; } // calculates the transform matrix that converts // from local coordinates to the coordinate space // of the parent transform public Matrix calculateLocalToParentMatrix() { // Matrix.translate creates a translation matrix // that shifts by (localPosition.x, localPosition.y) // Matrix.rotate rotates by localRotation radians // Matrix.scale scales by a factor of (localScale.x, localScale.y) // These are the basic transforms that are described previously // in this article return Matrix.translate(localPosition) * Matrix.rotate(localRotation) * Matrix.scale(localScale); } // gets the matrix that converts from local // coordinates to world coordinates public Matrix getLocalToWorldMatrix() { // if the dirty flag is set, the the // localToWorldMatrix is out of date // and needs to be reclaculated if (isDirty) { if (parent == null) { // if the parent is null then the parent is // the world so the localToWorldMatrix // is the same as local to parent matrix localToWorldMatrix = calculateLocalToParentMatrix(); } else { // if there is a parent, then the localToWorldMatrix // is calcualted recursively using the parent's localToWorldMatrix // concatenated with the local to parent matrix localToWorldMatrix = parent.getLocalToWorldMatrix() * calculateLocalToParentMatrix(); } // clear the dirty flag since the // matrix is now up to date isDirty = false; } return localToWorldMatrix; } public Matrix getWorldToLocalMatrix() { if (isInverseDirty) { // the inverse is out of date // so it needs to be updated // the worldToLocalMatrix is the inverse of // the localToWorldMatrix worldToLocalMatrix = getLocalToWorldMatrix().inverse(); // clear the dirty flag since the // matrix is now up to date isInverseDirty = false; } return worldToLocalMatrix; } // transforms a point from local coordinates to world coordinates public Vector2 transformPoint(Vector2 point) { // matrix multiply padding the extra element with a 1 Matrix1x3 transformResult = getLocalToWorldMatrix() * Matrix1x3(point.x, point.y, 1); return new Vector2(transformResult[1,1], transformResult[1,2], transformResult[1,3]); } // transforms a direction from local coordinates to world coordinates public Vector2 transformDirection(Vector2 point) { // matrix multiply padding the extra element with a 0 // notice that the worldToLocalMatrix is used here // and the point is multiplied as a row matrix before the // transform matrix. This is the proper way to transform // directions as described before in this article Matrix3x1 transformResult = Matrix3x1(point.x, point.y, 0) * getWorldToLocalMatrix(); return new Vector2(transformResult[1,1], transformResult[2,1], transformResult[3,1]); } // transforms a point from world coordinates to local coordinates public Vector2 inverseTransformPoint(Vector2 point) { // same logic as transformPoint only with the inverse matrix Matrix1x3 transformResult = getWorldToLocalMatrix() * Matrix1x3(point.x, point.y, 1); return new Vector2(transformResult[1,1], transformResult[1,2], transformResult[1,3]); } // transforms a direction from world coordinates to local coordinates public Vector2 inverseTransformDirection(Vector2 point) { // same logic as transformDirection only with the inverse of the // inverse localToWorldMatrix which is just the localToWorldMatrix Matrix3x1 transformResult = Matrix3x1(point.x, point.y, 0) * getLocalToWorldMatrix(); return new Vector2(transformResult[1,1], transformResult[2,1], transformResult[3,1]); } public Vector2 getLocalPosition() { return localPosition; } // sets the position relative to the parent // and marks the transform as dirty public void setLocalPosition(Vector2 value) { localPosition = value; // set the dirty flag since the localToWorldMatrix needs to be updated setDirty(); } /* localRoation and localScale should also have getters and setters * like the local position does. Be sure to call setDirty in the * setters for each of them */ }

There is a lot of code there so let me explain the key points. First of all, I am using a dirty flag to indicate when the transform matrices are out of date. That allows the position, rotation, and scale to be changed multiple times and the matrix is only actually recalculated once it is requested, reducing the amount of uneeded recalculations. Also take note that the dirty flag is only set if the transform is not already dirty. This is to keep the

The inverse matrix is also stored in the transform to keep from having to calculate it everytime it is needed. It has its own dirty flag so it isn't calculated if it isn't needed.

Now that we have a transform class, let's see how it can be used. First of all, the

class Renderer extends Component { void render(graphics) { graphics.setTransformMatrix(this.gameObject.transform.getLocalToWorldMatrix()); graphics.drawSprite(this.frame); } // whatever else the renderer does // would go here ... }

Keep in mind the above code doesn't account for the view transform. Think of the view transform as the camera for the scene. If you wanted to be able to scroll the view of the scene you should specify a camera object and multiply all transforms by the

graphics.setTransformMatrix( camera.transform.getWorldToLocalMatrix() * gameObject.transform.getLocalToWorldMatrix());

The transform could be used in game logic code

class Mine extends Behavoir { void update(float deltaTime, InputState input) { // transforming the zero vector gets that // transform's origin in world coordinates Vector2 playerWorldPosition = player.transform.transformPoint(new Vector2(0.0, 0.0)); // take the players world position and convert it the mine's // local coordinates Vector2 playerLocalPosition = this.transform.inverseTransformPoint(playerWorldPosition); // since playerLocalPosition is the players position relative to the mine // the magnitude of the position is the distance from the mine to the player if (playerLocalPosition.magnitudeSquared() < sensorRange * sensorRange) { this.detonate(); } } }

Transforms may take a while to fully grasp but once understood help simplify problems that would otherwise be a huge challenge to tackle. Hopefully this article helped you understand how to use them in your games better.

For this article, I'm assuming you know all about vectors, points, dot and cross products. We'll cover some quick properties of polynomials, some things about some basic curves, and then go over intersection tests.

The Stone-Weierstrass theorem states that any continuous function defined on a closed interval can be uniformly approximated as closely as desired with a polynomial function. For graphics and game development, that means that most 2D objects we deal with can be expressed as polynomials. The key to these intersection test methods is using that to our advantage.

In order to find all our intersection points, we need to know how many intersections there can be between 2 polynomials. Bezout's theorem states that for a planar algebraic curve of degree

Bezout's theorem also extends to surfaces. A surface of degree

The key here is to identify how to count intersections. For example, if 2 curves are tangent, they intersect twice. If they have the same curvature, then they intersect 3 times. Simple intersections (not tangent and not self-intersecting) are counted once. So how do we count intersections at infinity? Using homogeneous coordinates, of course!

Counting intersections at infinity sounds hard, but we can use

This helps us with define infinities with finite numbers. Let's say we have the point \((2,3,1)\) in homogeneous coordinates. If we changed the weight coordinate W so that we had \((2,3,0.5)\), the 2D point would be \((4,6)\). Smaller weights mean larger projected coordinates. As the weight approaches zero, the 2D point gets larger and larger, until at W = 0, the projected point is at infinity. Thus, any homogeneous point \((X,Y,0)\) is a point at infinity.

Usually, polynomials have terms that differ in their algebraic degrees. Some are quadratic, some cubic, some constant, etc. The polynomial takes the following form:

\[f(x,y) = \sum_{i+j\le n} a_{ij}x^iy^j = 0 \]

However, the same curve can be expressed in homogeneous form by adding a homogenizing variable

\[f(X,Y,W) = \sum_{i+j+k = n} a_{ij}X^iY^jW^k = 0 \]

Here, every term in the polynomial is of degree

To use polynomials effectively, we need to be familiar with how they can be expressed. There are basically 3 types of equations that can be used to describe planar curves: parametric, implicit, and explicit. If you've only had high-school math, you've probably dealt with explicit curves a lot and not so much with the others.

Let's start with a very basic curve: the line. It's a degree-1 polynomial. There are many definitions of this kind of curve. Some are helpful for games and others...not so much. Some you may have seen in algebra class and others you may be seeing for the first time.

In secondary schools, they use the first 2 methods almost exclusively, but these turn out to be the least helpful for computing. Here, I've tried to order the definitions from least helpful to most helpful for our needs.

Sometimes it's helpful to illustrate how useful these forms can be with a motivating example. This will hopefully help to connect formulas like dot product to traditional algebra and polynomials so you can see the connections.

In 3D, we can't represent a line as an implicit equation like ax+by+c=0. It becomes a parametric equation:

\[ L(x(t),y(t),z(t)) = \begin{cases} x&=x_0+ut \\ y&=y_0+vt \\ z&=z_0+wt \end{cases}\]

where \(\vec{a}=(u,v,w)\) defines the vector along the line and \(P_0=(x_0,y_0,z_0)\) defines a point on the line. These two things define the line completely.

What we want is an expression for the closest point on that line to another given point, which we'll name \(Q\). If we drew a vector from the closest point on the line to \(Q\) we would see that this vector would be perpendicular to \(\vec{a}\), or in other words, the dot product will be zero. This point will be unique on the line, so we should be able to solve for it algebraically. Remember the equation for the general point is \(P=P_0+\vec{a}t\).

The dot product equation would look like this:

\[ \begin{aligned} \vec{a} \cdot (Q-P) &= 0 \\ (u,v,w) \cdot (dx-ut,dy-vt,dz-wt) &= 0 \end{aligned} \]

where \( (dx,dy,dz) = (x_Q-x_{P_0},y_Q-y_{P_0},z_Q-z_{P_0} ) \). If we expand this out and solve for \(t\), we get the following expression:

\[ t = \frac{\vec{a} \cdot (Q-P_0)}{\vec{a} \cdot \vec{a}} \]

If we restrict \(\vec{a}\) to be a unit vector, the denominator will become 1 and vanish. This gives us a very simple expression for the point on the line, \( P = P_0 + \left [ \vec{a} \cdot (Q-P_0) \right ] \vec{a} \). The parameter multiplied to the vector is simply the dot product of the vector \(\vec{a}\) with the vector from the known point on the line to the given point. This serves to show 2 things:

- The dot product helps us find components of vectors in certain directions. We can see that what we get back from the dot product is the component of the \((Q-P_0)\) vector in the direction of \(\vec{a}\) vector. This is consistent with the linear algebra definition.
- If we're clever, we can use parametric and implicit equations to build polynomials and solve for the parameter values using algebra. These parameter values will help us find the points they're related to, if we need that information.

Probably one of the most useful forms, the implicit form is nice for a few reasons. First, we can express it in homogeneous coordinates: \(aX+bY+cW=0\). From here, we can define the coordinates as an ordered triple \((a,b,c)\) and then modify the equation to use the dot product so we can use a simple test if a point lies on a given line in implicit form:

\[L(a,b,c) \cdot P(X,Y,W) = 0\]

bool pointOnLine = fabs(dot(L,P)) < epsilon

Although this isn't strictly an intersection test, this method is valuable enough to mention here. Although points that lie on a line satisfy the equation \(L(a,b,c) \cdot P(X,Y,W) = 0\), if the point is not on the line, the scalar value can let us know on which side of the line the point lies by the sign. If \(L(a,b,c) \cdot P(X,Y,W) > 0\), the point lies on one side of the line; if \(L(a,b,c) \cdot P(X,Y,W) < 0\) it lies on the other. This is a convenient way of dividing a plane into 2 separate regions.

This scalar value is also a scaled distance from the line. If only the relative position of a point from the closest point on a line is desired, then this scalar can suffice. However, if true distance is required, then the full signed distance formula for a point \((x,y)\) from an implicit line defined by the ordered triple \((a,b,c)\) is given by the following:

\[D = \frac{ax+by+c}{\sqrt{a^2+b^2}}\]

\[\begin{aligned}A&=-v\\B&=u\\C&=vx_0-uy_0\\\end{aligned}\]

We want to use this notion of direction to apply to the signed distance formula so we can tell what side of the line we are on when looking in the direction of the line. Using this convention, we can say that if the sign of the dot product is greater than zero, it lies on the left side of the line. If the dot product is less than zero, it lies on the right side.

The reader should note that if this "right is less than zero, left is greater than zero" convention bothers anyone, you can flip this by multiplying the A, B, C coefficients by -1 to make "right is > 0, left is < 0".

This ordered triple form a convenient way to define coordinate axes. For example, the x-axis can be defined as \((0,1,0)\) and the y-axis can be defined as \((1,0,0)\). The ordered triple is also really nice because we can compute it really easily if we have 2 points. Let's go back to our vector math to see how this might work.

In the above picture, two vectors \(\vec{a}\) and \(\vec{b}\) are crossed to produce a 3rd vector \(\vec{c} = \vec{a}\times\vec{b}\). This vector \(\vec{c}\) is orthogonal to the vectors \(\vec{a}\) and \(\vec{b}\), meaning that \(\vec{c}\cdot\vec{a} = 0\) and \(\vec{c}\cdot\vec{b} = 0\). This is important for this next neat trick.

Imagine we have points \(P_1\) and \(P_2\). To tie back into the example above, let's let \(\vec{a} = P_1 = (x_1,y_1,w_1)\) and \(\vec{b} = P_2 = (x_2,y_2,w_2)\). If we cross these vectors, we get a vector \(\vec{c} = \vec{a} \times\vec{b}\) for which \(\vec{c}\cdot\vec{a} = 0\) and \(\vec{c}\cdot\vec{b} = 0\). Since \(\vec{a}\) and \(\vec{b}\) are points, we could say that \(\vec{c} = (a,b,c)\), an implicit ordered triple that passes through both points. We can check that the line given by \((a,b,c)\) passes through both points by the dot product method.

\[L(a,b,c) = P(X_1,Y_1,W_1) \times P(X_2,Y_2,W_2)\]

Using the same kind of logic, we can get the point at which 2 lines intersect. To use the cross product example again, let's let the vectors be the ordered triples of 2 implicit lines, \(\vec{a} = (a_1,b_1,c_1)\) and \(\vec{b} = (a_2,b_2,c_2)\). Again, crossing these vectors yields a vector \(\vec{c} = \vec{a} \times\vec{b}\) for which \(\vec{c}\cdot\vec{a} = 0\) and \(\vec{c}\cdot\vec{b} = 0\). Since \(\vec{a}\) and \(\vec{b}\) are lines, we could say that \(\vec{c} = (X,Y,W)\), the intersection point of both lines. This point has to be on both lines, and we can verify that using the dot-product method.

\[P(X,Y,W) = L(a_1,b_1,c_1) \times L(a_2,b_2,c_2)\]

Lines are degree-1 algebraic curves. Bezout's theorem states that 2 lines must intersect at exactly 1 point. So what happens with parallel lines? Well, if we use homogeneous coordinates, the intersection point will take the form \((X,Y,0)\), a point at infinity. Bezout's theorem still holds for that case.

Vector P = cross(new Vector(a1,b1,c1), new Vector(a2,b2,c2)); // check to see if the lines are parallel if(P[2] != 0) { // this is to transform back to (x,y) from (X,Y,W) x = P[0] / P[2]; y = P[1] / P[2]; }

Sometimes it's advantageous to define some curves as parametric and some as implicit to solve for intersections. Most times, it's better to define the simpler curve as parametric and the more complex curve as implicit, if possible. This method solves for all algebraic intersections, which may or may not be "real" intersections.

There are cases where there might be multiple roots, in which case we have to reconcile this with our intersection counting methods above (i.e. we have to count some of them differently). Let's take a fairly common case: line-circle intersection. Bezout's theorem says we will get 2 intersections since a circle is a degree-2 algebraic curve. The parametric line at point \((x_0,y_0)\) with direction \((c,d)\) and implicit circle centered at \((a,b)\) with radius \(r\) are defined as follows:

\[ \begin{aligned} C(x,y) &= (x-a)^2+(y-b)^2-r^2 = 0 \\

L(x(t),y(t)) &= \begin{cases} x &= x_0+ct \\ y&=y_0+dt \end{cases} \end{aligned}\]

We can substitute the parametric equations into the implicit equation to get a polynomial in terms of the parameter \(t\). The roots of the polynomial are the parameter values at which the line intersects the circle. Substituting, we get:

\[

\begin{aligned}

0 &= (x_0+ct-a)^2+(y_0+dt-b)^2-r^2 \\

&= [c^2t^2+2c(x_0-a)t+(x_0-a)^2]+[d^2t^2+2d(y_0-b)t+(y_0-b)^2]-r^2 \\

&= (c^2+d^2)t^2+2[c(x_0-a)+d(y_0-b)]t + [(x_0-a)^2+(y_0-b)^2-r^2] \\

&= At^2+Bt+C \\

\end{aligned}

\]

This looks promising. The quadratic formula can solve really quickly for the roots of the polynomial:

\[t=\frac{-B\pm\sqrt{B^2-4AC}}{2A}\]

The result of this will either be 2 complex roots, 1 real root, or 2 real roots, depending on the discriminant \(B^2-4AC\). The 1 real root case means the line is tangent to the circle, which according to our intersection counting rules above, we count twice.

That might be nice mathematically, but what does that mean for us? Well, the parameter \(t\) has to be a real number, so the real roots are values of \(t\) that the line intersects, so plugging the roots back into the parametric line gives us the Cartesian points of the intersection. What about the complex roots? Well, the restriction on \(t\) is that it has to be a real number, so in the case of 2 complex roots, we say the line doesn't intersect the circle.

A = c*c+d*d; B = 2*(c(x0-a)+d(y0-b)); C = (x0-a)*(x0-a)+(y0-b)*(y0-b)-r*r; disc = B*B-4*A*C; if(disc == 0) { t = -B/(2*A); x1 = x0 + c*t; y1 = y0 + d*t; return [Point(x1,y1)]; } else if(disc > 0) { t1 = (-B+sqrt(disc))/(2*A); t2 = (-B-sqrt(disc))/(2*A); x1 = x0 + c*t1; y1 = y0 + d*t1; x2 = x0 + c*t2; y2 = y0 + d*t2; return [Point(x1,y1), Point(x2,y2)]; } else return [];

Line-sphere intersection is just like line-circle intersection. Algebraically, a sphere is also a degree-2 polynomial, and Bezout's theorem says we will still only get 2 intersections. The parametric line at point \((x_0,y_0,z_0)\) with direction \((p,q,r)\) and implicit sphere centered at \((a,b,c)\) with radius \(R\) are defined as follows:

\[ \begin{aligned} S(x,y,z) &= (x-a)^2+(y-b)^2+(z-c)^2-R^2 = 0 \\

L(x(t),y(t),z(t)) &= \begin{cases} x &= x_0+pt \\ y&=y_0+qt \\ z&=z_0+rt \end{cases} \end{aligned}\]

We can substitute the parametric equations into the implicit equation to get a polynomial in terms of the parameter \(t\). The roots of the polynomial are the parameter values at which the line intersects the circle. Substituting, we get:

\[

\begin{aligned}

0 &= (x_0+pt-a)^2+(y_0+qt-b)^2+(z_0+rt-c)^2-R^2 \\

&= [p^2t^2+2p(x_0-a)t+(x_0-a)^2]+[q^2t^2+2q(y_0-b)t+(y_0-b)^2]+[r^2t^2+2r(z_0-c)t+(z_0-c)^2]-R^2 \\

&= (p^2+q^2+r^2)t^2+2[p(x_0-a)+q(y_0-b)+r(z_0-c)]t + [(x_0-a)^2+(y_0-b)^2+(z_0-c)^2-R^2] \\

&= At^2+Bt+C \\

\end{aligned}

\]

The solution method is exactly the same as for line-circle intersections.

This is a simpler case than the line-circle intersection, although it involves a surface and a curve. A plane is algebraically a degree-1 polynomial in implicit form, so according to Bezout's theorem, they should intersect at exactly 1 point. We take the plane in implicit form and the line in parametric form and apply our method as above:

\[

\begin{aligned}

P(x,y) &= ax+by+cz+d=0 \\

L(x(t),y(t),z(t)) &= \begin{cases} x &= x_0+ut \\ y&=y_0+vt \\ z &= z_0+wt \\ \end{cases} \end{aligned}\]

We substitute the parametric equations into the implicit form and solve for \(t\) as before:

\[

\begin{aligned}

P(x(t),y(t),z(t)) = 0 &= a(x_0+ut)+b(y_0+vt)+c(z_0+wt)+d \\

&= (ax_0+by_0+cy_0+d) + (au+bv+cw)t \\

t &= -\frac{ax_0+by_0+cy_0+d}{au+bv+cw} \\

&= -\frac{(a,b,c,d)\cdot(x_0,y_0,z_0,1)}{(a,b,c,d)\cdot(u,v,w,0)}

\end{aligned}

\]

By substituting the value for \(t\) into the parametric line, we get the intersection point of the line and the plane. If the line runs parallel to the plane, the dot product of the plane normal and the line direction (which is the denominator) will be zero.

Vector4 PL = new Vector4(a,b,c,d); Vector4 P = new Vector4(x0,y0,z0,1); Vector4 D = new Vector4(u,v,w,0); den = dot(PL,D); if(den != 0) // check for div by zero { t = -dot(PL,P) / den; x = x0 + u*t; y = y0 + v*t; z = z0 + w*t; }

One method of ray-triangle intersection is basically a two-step process, the first step being computing the point on the plane that intersects the ray (line). The next is to determine if the point lies inside the triangle. Barycentric coordinates is a common method of doing this, but another method would be to create lines in the plane with directions such that a point inside the triangle would be on the same side of each line (i.e. the point would be on the left or right sides of each line). The side of the line can be computed using the scaled signed distance as detailed above.

Just to illustrate how general this method is, let's take a more advanced example:

\[ \begin{aligned} C(x,y) &= x^2+y^2-1 = 0 \\

L(x(t),y(t)) &= \begin{cases} x &= 2t-1 \\ y&=8t^2-9t+1 \end{cases} \end{aligned}\]

The picture above shows the curves. The red curve is the parametric curve and the black curve is the implicit curve. As you can see, there are 4 real intersections, and since both curves are degree-2 we should end up with all real roots.

Inserting definitions of the parametric equations into the implicit form, we get:

\[

\begin{aligned}

f(x(t),y(t)) &= (2t-1)^2+(8t^2-9t+1)^2-1 \\

&= 64t^4-144t^3+101t^2-22t+1 \\

&= 0

\end{aligned}

\]

This is a quartic polynomial, so a more advanced numeric root finding method needs to be used, like bisection, regula falsi, or Newton's method. The roots of the polynomial are t = 0.06118, 0.28147, 0.90735, and 1.0. We do have all real roots, so Bezout's theorem is satisfied.

If the parametric curves are rational, then this method needs to be slightly modified. Rational parametric curves are usually of the form:

\[x = \frac{a(t)}{c(t)}, \, y = \frac{b(t)}{c(t)}\]

We can use homogeneous coordinates here really well. Since we have the mapping \((x,y) = (X/W,Y/W)\), we can define each homogeneous coordinate as \(X = a(t),\,Y = b(t),\,W=c(t)\). We also need to modify the implicit curve to handle homogeneous coordinates too, but this isn't very hard. Let's illustrate with an example:

\[

\begin{aligned}

S(x,y) &= x^2+2x+y^2+1 = 0 \\

P(x(t),y(t)) &= \begin{cases}

x &= \frac{t^2+t}{t^2+1} \\

y &= \frac{2t}{t^2+1} \\

\end{cases}

\end{aligned}\]

Our parametric curve in homogeneous coordinates is \(X = t^2+t\),\(Y=2t\),and \(W={t^2+1}\). Our implicit curve can be changed to use homogeneous coordinates by substituting in the mapping \((x,y) = (X/W,Y/W)\):

\[

\begin{aligned}

S &= \left(\frac{X}{W}\right)^2+2\left(\frac{X}{W}\right)+\left(\frac{Y}{W}\right)^2+1 \\

&= X^2+2XW+Y^2+W^2 \\

\end{aligned}

\]

We can substitute our homogeneous coordinate equation into our modified implicit function:

\[

\begin{aligned}

S &= (t^2+t)^2 + 2(t^2+t)(t^2+1) + (2t)^2 + (t^2+1)^2 \\

&= 4t^4+6t^3+5t^2+4t+1 \\

&= 0 \\

\end{aligned}

\]

We get 2 real roots at t = 0.3576 and t = 1, and 2 complex roots. Both curves are degree-2, so we expect 4 intersections. We have 4 intersections total, so Bezout's theorem is satisfied.

The way to solve intersections of 2 parametric curves is via

A Bezier curve is just a degree-n Bernstein polynomial, which means it's just a regular polynomial of a different form. The above methods work well for Bezier curves, but they are more efficient and more numerically stable if modified slightly to take advantage of the Bernstein form. We won't cover this here, but there is literature out there on this as well.

These common methods of intersection testing can greatly aid any game programmer, whether purely for graphics or other game-specific logic. Hopefully you can make some use of these simple, yet powerful techniques.

A lot of the information here was taught to me by Dr. Thomas Sederberg (associate dean at Brigham Young University, inventor of the T-splines technology, and recipient of the 2006 ACM SIGGRAPH Computer Graphics Achievement Award). The rest came from my own study of relevant literature. The article image is from his self-published book,

Written by Michael Schmidt Nissen

Version 0.5.2, December 18th 2014

The damped spring based on Hooke's law of elasticity is often regarded as a black sheep in the field of game physics, and some papers and websites even warn the reader against implementing them. There are a couple of good reasons for this.

The most straight-forward reason is that springs are notoriously difficult to tune. It can be very time consuming to adjust the spring stiffness and damping coefficients in order to make a simulation behave just right. For large systems of interconnected particles and springs the task can be daunting. If a variable is changed, for instance a spring coefficient, the mass of a particle, or the simulation time-step, you might have to re-tune the entire system from scratch.

A more profound reason for avoiding springs is that they are potentially unstable. When increasing the spring coefficients above some unknown threshold, the simulation will start to oscillate chaotically and explode in all directions. Everyone who ever worked with springs has experienced this frustrating phenomenon. The problem is that there’s no good method to determine in advance if a given spring system will show this behavour or not. Thus, the programmers or game designers ability to balance a spring system often rely on experience and gut feeling, which makes it look more like black magic than actual physics.

In this article I’ll introduce a modified spring equation that is easy to tune, displays maximum stiffness and damping, and is guaranteed to keep the spring system stable under all conditions. The math behind the equation is surprisingly simple, which makes me wonder why it doesn’t appear to have been mentioned or implemented earlier.

The improved spring equation stems from trying to answer a simple question: Is it possible to make a spring that reaches its rest state in just one simulation loop? It turns out that the answer is yes! This wouldn’t be possible in the real world, since it would require infinitely high stiffness and damping coefficients or zero mass. But in the not-quite real world of discrete, time-stepping physics simulation, this is achievable.

To explain how this is done, let's look at a one dimensional harmonic oscillator with a damped spring attached to a fixed point in one end and to a moving body in the other. Assume that the body is displaced distance

Distance can be expressed as velocity multiplied by time

\[ \large x = -{v} {\Delta t} \]

The minus sign simply means that we're moving in the opposite direction of the displacement. Velocity can be expressed as acceleration multiplied by time

\[ \large x = -{a} {\Delta t^{2}} \]

Newton's 2nd law of motion states that acceleration equals force over mass

\[ \large x = -\dfrac{f}{m} {\Delta t^{2}} \]

Now we can isolate force

\[ \large f= -\dfrac{m}{\Delta t^{2}} x \]

And since the spring force in Hooke's law of elasticity is defined as

\[ \large f= -k x \]

It becomes immediately clear that

\[ \large k = \dfrac{m}{\Delta t^{2}} \]

Which is the exact spring stiffness coefficient value required to move the particle back to its rest position in one simulation loop. However, since we are not doing anything to stop the particle, it will keep oscillating back and forth through the rest position. We need to add damping, which is done with the second part of the spring equation.

Velocity is equal to acceleration multiplied by time step

\[ \large v = -{a} {\Delta t} \]

Which is equal to force over mass

\[ \large v = -\dfrac{f}{m} {\Delta t} \]

When isolating force we get

\[ \large f = -\dfrac{m}{\Delta t} v \]

And since the damping force is defined

\[ \large f = -d v \]

We immediately see by inspection that the damping coefficient is

\[ \large d = \dfrac{m}{\Delta t} \]

This is the exact damping coefficient value needed to stop the particle in one simulation loop. Now we can write the complete damped spring equation.

\[ \large f= -\dfrac{m}{\Delta t^{2}} x -\dfrac{m}{\Delta t} v \]

This spring equation has some very interesting properties. The first thing we notice is the lack of coefficients. We've simply replaced

Now we have a really powerful spring equation. It is easy to implement, very stable, and reaches its rest state in just one loop. But in our quest towards a better spring it has lost its springiness. We need to somehow get the softness and bounciness back again, and for this purpose we will re-introduce the spring and damping coefficients. To avoid confusion, the new coefficients are named

\[ \large f= -\dfrac{m}{\Delta t^{2}} C_{k} x -\dfrac{m}{\Delta t} C_{d} v \qquad [ 0 \leq C_{k}, C_{d} \leq 1 ] \]

As we can see, the new coefficients are simply the fraction of completely rigid behavior we would like the spring to display. Soft, bouncy springs would usually have values in the range of 0.001 – 0.00001. In the other end of the scale, values of just 0.1 – 0.01 is enough to display rigid behavior. Setting both values to 1.0 would of course still satisfy the constraint in one loop.

Please notice that spring behavior is determined exclusively by these two coefficients. Particle mass or simulation time-step has no influence on how the spring behaves, and changing them wouldn't make it necessary to tune the system!

Interestingly, the spring will get less rigid and less stable if we increase

\[ \large k_{max} = \dfrac{m}{\Delta t^{2}} \]

\[ \large d_{max} = \dfrac{m}{\Delta t} \]

This allows us to simplify the spring equation

\[ \large f= -k_{max} C_{k} x -d_{max} C_{d} v \]

There is an important conclusion to be drawn from this. It is a misunderstanding to think that spring stiffness and damping can be increased towards infinity by simply increasing

It is only slightly more complicated to constrain two free-moving particles with the improved spring. To do this, we need to introduce the concept of

\[ \large m_{red} = \dfrac{m_1 m_2}{(m_1 + m_2)} \]

Since the inverse mass quantity is often already pre-computed for other purposes, it can also be useful to define reduced mass as

\[ \large m_{red}= \dfrac{1}{ ( \dfrac{1}{m_1} + \dfrac{1}{m_2} ) } \]

For two connected particles, the maximum coefficient values are

\[ \large k_{max} = \dfrac{m_{red}}{\Delta t^{2}} \]

\[ \large d_{max} = \dfrac{m_{red}}{\Delta t} \]

When replacing mass with reduced mass we get

\[ \large f= -\dfrac{m_{red}}{\Delta t^{2}} C_{k} x -\dfrac{m_{red}}{\Delta t} C_{d} v \]

This spring equation will bring the two connected particles to equilibrium distance in one loop and make them stand still relative to each other. However, since energy and momentum are conserved, the particles may rotate around each other, which will make the bodies come to rest at a larger distance depending on how fast they rotate.

The spring equation can also be used for angular springs, also commonly named rotational or torsional springs. Rather than keeping two particles at a fixed distance by applying opposite equal forces, the angular spring will try to keep two rotating bodies at a fixed angle by applying opposite equal torques. The equation introduces the concept

\[ \large I_{red} = \dfrac{I_1 I_2}{(I_1 + I_2)} \]

Or alternatively

\[ \large I_{red}= \dfrac{1}{ ( \dfrac{1}{I_1} + \dfrac{1}{I_2} ) } \]

The maximum coefficient values for angular springs are

\[ \large k_{max} = \dfrac{I_{red}}{\Delta t^{2}} \]

\[ \large d_{max} = \dfrac{I_{red}}{\Delta t} \]

When replacing the variables of linear motion with those of angular motion we get

\[ \large \tau = -\dfrac{I_{red}}{\Delta t^{2}} C_{k} \theta -\dfrac{I_{red}}{\Delta t} C_{d} \omega \]

This spring equation will keep two rotating bodies at any given rest angle. If both coefficients are set to 1, the constraint will be solved in one loop. The spring allows for angular displacements larger than one full turn, making it possible to “wind up” bodies like the coil spring of a mechanical clock.

Today a lot of physics engines and simulation methods are based on impulses - direct changes in velocities - rather than forces and acceleration. The linear and angular spring equations described above works equally well if we redesign them to work with impulses. Here are the equations without further ado

\[ \large J_{linear}= -\dfrac{m_{red}}{\Delta t} C_{k} x - C_{d} v \]

\[ \large J_{angular} = -\dfrac{I_{red}}{\Delta t} C_{k} \theta - C_{d} \omega \]

In this particular case there is no difference between the force-based and impulse-based spring. They should return the exact same results.

When connecting a multitude of springs and particles into larger bodies, we run into the same trouble as any other type of distance and velocity constraint. Rather than cooperating, the constraints tend compete against each other, and this spells trouble. When a spring moves two particles to satisfy distance and velocity, it usually means dissatisfying one or several other springs. It is outside the scope of this article to dig deeply into this problem, but the author would like to provide a bit of advice on how to prevent the worst disasters.

When two or more springs are connected to the same particle, which is the case in any kind of rope, mesh, or deformable body, setting the coefficients to the maximum value of 1.0 will lead to stability problems. Although the spring equation is stable when used with a single or two particles, this is sadly not the case for higher number of particles. The author of this article has after some lengthy tampering worked out that a safe upper limit for both the stiffness and damping coefficient is

\[ \large c_{max} \approx \dfrac{1}{n + 1} \]

Where

http://en.wikipedia.org/wiki/Reduced_mass]]>