# Tag: math

Entries for tag "math", ordered from most recent. Entry count: 58.

Warning! Some information on this page is older than 5 years now. I keep it for reference, but it probably doesn't reflect my current knowledge and beliefs.

Pages: 1 2 3 4 ... 8

# Vector Register Size - Diagram

Thu
29
Sep 2011

It may be hard to imagine and remember what is the exact number of bits, bytes, words or floats in some piece of data, like a SIMD register. So today I've made following diagram/cheatsheet:

Here you can find its "source" in OpenOffice Draw format: Vector_register_size.odg.

# DevMeeting and Demoscene Night

Sun
24
Jul 2011

Yesterday I attended two events related to my interests. First was about collision detection in JavaScript games, organized by Marek Paw³owski from devmeetings.pl. During this almost 12-hour workshop we could learn about JavaScript, as well as some theory of 2D collision detection and space partitioning techniques. Nothing new for me, but I liked formula of this workshop. We could learn some theory from slides, but most of the time we had some tasks to code on our laptops - first to implement a simple library for aspects in JavaScript and then to code Asteroids game, including precise collision detection and QuadTree. Most of the code for this game was ready - we were given a framework, as well as libraries with math and a class for QuadTree. Our task was to just to connect it together. I think it's a good way of teaching programming in practice.

Web technologies are not my main interest, but I like JavaScript. I even think  it's the most beautiful scripting language in terms of syntax. Being able to freely draw 2D (using HTML Canvas) as well as 3D graphics (using WebGL) makes it a good technology for learning and prototyping geometry or gameplay algorithms. Marek also pointed on the DevMeeting an interesting conclusion that it's very easy to port algorithms from C/C++ to JavaScript. But on the other hand, Dab reminded me today that Lua is a scripting language with very similar features, but faster and more lightweight interpreter, so it's still better choice as a scripting language embedded in high-performance software like games.

Second event was Noc z Demoscen± - Demoscene Night. It took place in Fabryka Kot³ów club (former No Mercy). During that night we could see some olschool, as well as newschool demos presented on a a big screen and, what is most important, meet some nice people there. Thanks Voyager for organizing this event! It was great especially because the 3-day demoscene party RiverWash, which took place here in Warsaw in last 2 years, this year is in £ód¼ and I can't attend it. This night was very inspiring and now I feel like developing a technology to make a demo :)

Comments | #demoscene #webdev #javascript #math #events

# Bit Tricks with Modulo

Sat
26
Feb 2011

I like to emphasize that programming is not so similar to maths as some people say. We are used to thinking about numbers in decimal numeral system and the operations that seem most natural for us are arithmetic computations like additions, multiplication and division. On the other hand, computer works with numbers stored in binary system and bitwise operations are most efficient and natural for it. Another difference is modulo operation (the remainder after division) - a function not so frequently used in math classes, but very important in programming. That's what I'd like to talk about here.

Modulo in C++ can be done using % operator, but just as division, it is quite slow. So when doing a common operation of cycling between subsequent numbers from a limited range:

x2 = (x1 + 1) % n;

few alternatives have been developed for some special cases, using bitwise operations.

First, when the divisor is a power of 2, we can use bitwise and (operator &) to mask only least significant bits of the addition result. But this optimization can be done automatically by the compiler (at least my Visual C++ 2010).

// 7 = 0b111, so it masks 3 least significant bits.
// Equal to (x1 + 1) % 8
x2 = (x1 + 1) & 7;

When we only cycle between two values - 0 and 1 - we can use XOR (operator ^):

// Equal to (x1 + 1) % 2, for values 0, 1.
x2 = x1 ^ 1;

Finally, when we cycle between three values - 0, 1, 2 (like when indexing x, y, z components of a 3D vector), we can use the trick invented by Christer Ericson (at the end of the linked page). I came across it recently on Twitter and that's what inspired me to write this article.

// Equal to (x1 + 1) % 3, for values 0, 1, 2.
x2 = (1 << x1) & 3;

# XNA Math and Access Violation

Wed
16
Feb 2011

XNA Math is a great math library bundled with DirectX SDK. I've written about it here: Introduction to XNA Math, XNA Math Cheat Sheet. Recently I've started coding some small raytracer using this library (especially XMVECTOR type) and I came across an error:

Google didn't display many results when asked about this error in connection with "XNA Math" or "XMVECTOR", but I found a topic on our Polish gamdev forum: xnamath i dziwne crashe. Directed to other web places from there, I discovered that, just as I thought, this bug comes from misaligned address of the vector. XMVECTOR on PC is an alias for SSE __m128 type, which has to be always aligned to 16-byte boundary. Compiler seems to know about this alignment requirement and respects it whenever XMVECTOR is defined on the stack, passed as function parameter or returned from a function, but not when the vector is a member of some dynamically allocated class or struct. So this won't work:

struct Scene {
XMVECTOR dir_to_light;
};

Scene *scene = new Scene();
scene->dir_to_light = XMVector3Normalize(
XMVectorSet(1.5f, 1.f, -2.f, 0.f)); // Crash!

Adding special alignment declaration supported by Visual C++ - __declspec(align(16)) - doesn't help. it looks like operator new always aligns allocated objects to 8 bytes and there is no way to change this behavior other than redefining it to use your own custom memory allocator.

So I suppose the most simple and convenient way to walk around this problem is to use more "casual" types for storing vectors - the ones that do not use SSE and do not require alignment, like XMFLAOT3 (which is just struct with float x, y, z). Anyway XMVECTOR is designed to optimize computations when compiler is able to place it in the CPU registers. So my solution looks like this:

struct Scene {
XMFLOAT3 dir_to_light;
};

Scene *scene = new Scene();
XMVECTOR dir_to_light = XMVector3Normalize(
XMVectorSet(1.5f, 1.f, -2.f, 0.f));
XMStoreFloat3(&scene->dir_to_light, dir_to_light);

...

XMVECTOR n_dot_l = XMVectorSaturate(XMVector3Dot(normal, dir_to_light));

# Finding Polygon of Plane-AABB Intersection

Tue
15
Feb 2011

While doing volumetric rendering according to the article Volume Rendering Techniques (Milan Ikits, Joe Kniss, Aaron Lefohn, Charles Hansen, Chapter 39 of "GPU Gems" book), I came across a problem of calculating intersection between a plane and an AABB (axis-aligned bounding box), which forms a 3- to 6-vertex polygon. The solution is actually described in the article, but here I'd like to share my implementation (using D3DX).

The solution is based on finding points of intersection between the plane and every edge of the box. So we first need a ray-to-plane intersection code. My function for this is based on: Ray-Plane Intersection.

// OutVD > 0 means ray is back-facing the plane
// returns false if there is no intersection because ray is perpedicular to plane
bool ray_to_plane(const D3DXVECTOR3 &RayOrig, const D3DXVECTOR3 &RayDir, const D3DXPLANE &Plane, float *OutT, float *OutVD)
{
*OutVD = Plane.a * RayDir.x + Plane.b * RayDir.y + Plane.c * RayDir.z;
if (*OutVD == 0.0f)
return false;
*OutT = - (Plane.a * RayOrig.x + Plane.b * RayOrig.y + Plane.c * RayOrig.z + Plane.d) / *OutVD;
return true;
}

Now having a segment A-B between two points, we can calculate ray_to_plane intersection where RayOrig=A, RayDir=B-A and we get parameter T. The segment is a subset of the infinite ray. The intersection point lies on this segment when 0 <= T <= 1. The position of this point is A+T*B. Now it's time for plane-box intersection code:

// Maximum out_point_count == 6, so out_points must point to 6-element array.
// out_point_count == 0 mean no intersection.
// out_points are not sorted.
void calc_plane_aabb_intersection_points(const D3DXPLANE &plane,
const D3DXVECTOR3 &aabb_min, const D3DXVECTOR3 &aabb_max,
D3DXVECTOR3 *out_points, unsigned &out_point_count)
{
out_point_count = 0;
float vd, t;

// Test edges along X axis, pointing right.
D3DXVECTOR3 dir = D3DXVECTOR3(aabb_max.x - aabb_min.x, 0.f, 0.f);
D3DXVECTOR3 orig = aabb_min;
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_max.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_max.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;

// Test edges along Y axis, pointing up.
dir = D3DXVECTOR3(0.f, aabb_max.y - aabb_min.y, 0.f);
orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_max.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_max.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;

// Test edges along Z axis, pointing forward.
dir = D3DXVECTOR3(0.f, 0.f, aabb_max.z - aabb_min.z);
orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
orig = D3DXVECTOR3(aabb_max.x, aabb_max.y, aabb_min.z);
if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
out_points[out_point_count++] = orig + dir * t;
}

But that's not all. Returned points are not sorted! If we don't sort them, we may in some cases end up with a mess like this:

Here comes another clever algorithm - the one for sorting vertices (counter?)clockwise. I came up with it on my own and it's a bit similar to the algorithm for finding a convex hull. We are given a set of points that all lie in a single plane and already form a convex polygon. We want to sort them clockwise "around some center point". Beginner could probably think here about calculating such center point, then for each vertex calculating and comparing some angle using atan2, of course in degrees ;)

But that's not the solution. What we need here is a way to tell whether one vector is more "to the right" then the other, relative to some reference point. Whether one vector lies on the right-hand side of the other can be easily told in 2D using the sign of a scalar returned by 2D cross-product: a.x*b.y - b.x*a.y. In 3D it's not that easy. My idea is to use 3D cross-product. It returns 3D vector, but has similar property - it is anticommutative, which means: cross(a,b) == -cross(b,a). Now if only I could tell whether the direction of the resulting vector is "right" or "opposite"... That's what the sign of dot product means! I'm just going to use dot product between the cross product result and the normal vector of the plane. Eureka! Now I have all necessary parts. I choose first vertex as a reference. I also use lambda expression from C++0x (supported in Visual C++ 2010) to more conveniently pass custom comparator to std::sort STL algorithm (otherwise I would have to write a functor). Here is the final solution:

#include <algorithm>

void sort_points(D3DXVECTOR3 *points, unsigned point_count, const D3DXPLANE &plane)
{
if (point_count == 0) return;

const D3DXVECTOR3 plane_normal = D3DXVECTOR3(plane.a, plane.b, plane.c);
const D3DXVECTOR3 origin = points[0];

std::sort(points, points + point_count, [&](const D3DXVECTOR3 &lhs, const D3DXVECTOR3 &rhs) -> bool {
D3DXVECTOR3 v;
D3DXVec3Cross(&v, &(lhs - origin), &(rhs - origin));
return D3DXVec3Dot(&v, &plane_normal) < 0;
} );
}

# Technology for Data Processing

Mon
04
Oct 2010

When doing some engineering work on a computer, no matter if gamedev or any other field, there is sometimes a need to process or visualize some tabular data, especially numbers, e.g. statistics about performance or something gathered during program execution. What technology is best for this purpose? At the moment I know about following solutions:

Spreadsheet software, like Microsoft Excel or OpenOffice Calc. They can import CSV files and draw great variety of plots, but for more advanced data processing it would be useful to have something more like a programming language.

At the other end of the spectrum, we can write normal C++ or C# programs to do the work. It can be hard though as we have to code everything on our own, including data structures, loading files and drawing plots. It would be nice to use some higher-level, scripting language with rich standard library, built-in data structures and the like.

Any scripting language can do this. For example, PHP (which I know the best) can be used as a normal scripting language, not only in connection with a web server. All in all, its name means "PHP Hypertext Prerocessor", because it's suited for processing strings and text files.

A technology designed specially for the purpose of crunching numbers is Matlab and its free alternative - Scilab. It provides its own programming language with convenient built-in data types like matrix and is also able to draw plots.

Python is something in between - a normal scripting language with weird but nice syntax and a language-level support for operations like array slicing and complex numbers. It looks like many scientists and engineers use it for computation and data analysis, because there are Python libraries designed for this, like Numpy and Scipy.

There is also The R Project metioned at the Nick Darnells' Blog. It looks like another environment with its own, different programming language, designed especially for statistical computing. It can also load data from files and draw plots.

And finally, some tasks can be accomplished in environments that allow playing around with computer graphics, like EvalDraw (with its own, C-like simple programming language) or Processing (with a language based on Java).

So there are many possibilities in this subject. I've played a bit with all of them at some time, but obviously learning chosen one thoroughly would require much time and effort. So maybe you can help me decide? Which solution would you recommend? Personally I feel a little more convinced to Python, because it is a general purpose language that I can use in different fields too, like coding Blender plugins.

# Music Analysis - Spectrogram

Wed
14
Jul 2010

I've started learning about sound analysis. I have some deficiencies in education when it comes to digital signal processing (greetings for the professor who taught this subject at our university ;) but Wikipedia comes to the rescue. As a starting point, here is a spectrogram I've made from one of my recent favourite songs: Sarge Devant feat. Emma Hewitt - Take Me With You.

Now I'm going to exaplain in details how I've done this by showing some C++ code. First I had to figure out how to decode an MP3, OGG or other compressed sound format. FMOD is my favourite sound library and I knew it can play many file formats. It took me some time though to find functions for fast decoding uncompressed PCM data from a song without actually playing it for all 3 minutes. I've found on the FMOD forum that Sound::seekData and Sound::readData can do the job. Finally I've finished with this code (all code shown here is stripped from error checking which I actually do everywhere):

Read full entry | Comments | #math #libraries #music #rendering #dsp #algorithms

# Calculating AABB of a Rotated 2D Rectangle

Wed
02
Jun 2010

Today I've solved an interesting geometrical problem of calculating an AABB rectangle (Axis-Aligned Bounding Box) around an oriented rectangle (rotated by an angle) in 2D.

It seems simple but it's easy to come up with suboptimal solution like calculating all 4 corners of the rotated rectangle or calling sin and cos functions multiple times. So I've coded a simple demonstration in ActionScript 3 (Flash). Here is my algorithm:

Input data:

// Center position of the rectangle.
private const m_PosX : Number, m_PosY : Number;
// Half-width and half-height of the rectangle.
private const m_HalfSizeX : Number, m_HalfSizeY : Number;
private var m_Orientation : Number;

Algorithm:

// corner_1 is right-top corner of unrotated rectangle, relative to m_Pos.
// corner_2 is right-bottom corner of unrotated rectangle, relative to m_Pos.
var corner_1_x : Number = m_HalfSizeX;
var corner_2_x : Number = m_HalfSizeX;
var corner_1_y : Number = -m_HalfSizeY;
var corner_2_y : Number =  m_HalfSizeY;

var sin_o : Number = Math.sin(m_Orientation);
var cos_o : Number = Math.cos(m_Orientation);

// xformed_corner_1, xformed_corner_2 are points corner_1, corner_2 rotated by angle m_Orientation.
var xformed_corner_1_x : Number = corner_1_x * cos_o - corner_1_y * sin_o;
var xformed_corner_1_y : Number = corner_1_x * sin_o + corner_1_y * cos_o;
var xformed_corner_2_x : Number = corner_2_x * cos_o - corner_2_y * sin_o;
var xformed_corner_2_y : Number = corner_2_x * sin_o + corner_2_y * cos_o;

// ex, ey are extents (half-sizes) of the final AABB.
var ex : Number = Math.max(Math.abs(xformed_corner_1_x), Math.abs(xformed_corner_2_x));
var ey : Number = Math.max(Math.abs(xformed_corner_1_y), Math.abs(xformed_corner_2_y));

var aabb_min_x : Number = m_PosX - ex;
var aabb_max_x : Number = m_PosX + ex;
var aabb_min_y : Number = m_PosY - ey;
var aabb_max_y : Number = m_PosY + ey;