# Tag: math

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

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.

# 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 planebool 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;    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;```

AabbTest01.swf - online Flash demo
AabbTest01.zip - full project source code

I believe it can be optimized further (not only by rewriting it in C++ :) Anyone knows how?

# Introduction to XNA Math

Tue
04
May 2010

I've recently shared my XNA Math Cheat Sheet Today I'd like to post some more information about XNA Math. It is a new C++ math library for gamedev distributed with DirectX SDK, intended to replace math support from D3DX. Its online documentation can be found here: XNA Math Library Reference and the offline one is in DXSDK\Documentation\DirectX9\directx_sdk.chm file. It's worth noticing that the documentation in recent DX SDK versions is split across two separate CHM files. windows_graphics.chm describes only graphics programming (Direct3D 9, 10 and 11, D3DX, DXGI, HLSL language, effects format, DDS format, X format etc.), while directx_sdk.chm is about all other stuff (DirectSetup, DXUT, XNA Math, XAudio2, XInput etc.).

XNA Math is different from old D3DX and I'm afraid also harder to use. But I believe it's better. It's your decision whether you start to use it or not, because it's completely separate library that can be used with DirectX 11, 10, 9 or without anyone of them. On the other hand, you can still link with old D3DX while coding in new DirectX 11, so it's all your choice.

By the way, I don't know why is this library called "XNA Math". It looks like it has nothing to do with the XNA technology. XNA is a .NET library with its own vector and matrix classes while XNA Math is pure native C++ library distributed with DirectX SDK. It looks like a misnomer.

Beginning with XNA Math is simple. You just #include <xnamath.h> and that's all. There is even no need to link with any LIB file - all code is inline in this header or INL files included by it (which sum up to 26530 lines :) It provides support for vectors (up to 4D), matrices (up to 4x4), colors, planes and quaternions - all made of 32-bit floats, just as we all use in gamedev.

But design ideas of this API are totally different from D3DX. There is now one main data type called XMVECTOR - a vector of 4 floats. When the library uses SSE (which is the default on PC), this type is simply a typedef to __m128 (if you know what I mean). So it must be always aligned to 16-byte boundary. As a consequence, it's safe to cast this type to any other one interpreted as a vector of floats (like float or D3DVECTOR), but not the other way around (because the pointer you cast might not be property aligned). This single type is used in all computations and interpreted as a 2D, 3D or 4D vector, plane, quaternion, part of 4x4 matrix or even a vector of unsigned integers, depending on the context. It may also feel strange at the beginning that this not-so-small object is usually passed and returned by value, not by reference.

Old DirectX versions used a D3DVECTOR and D3DMATRIX structures, which were extended by inheritance in D3DX library as D3DXVECTOR and D3DXMATRIX - structures with lots of overloaded operators and methods to operate on. The idea behind new DirectX 11 is that the graphics API is not dependent on any particular math library - it only uses float* or even void* and it's you job to use some library that would help you with math computations (like D3DX or XNA Math).

If you ever coded in SSE or other vector instruction set, you might know the concept of loading and storing. Here is an example that sets and gets vector as well as its components in many different ways.

```XMVECTOR v = XMVectorZero();
// v is now (0, 0, 0, 0)
v = XMVectorSplatOne();
// v is now (1, 1, 1, 1)
v = XMVectorReplicate(5.f);
// v is now (5, 5, 5, 5)
v = XMVectorSet(1.f, 2.f, 3.f, 4.f);
// v is now (1, 2, 3, 4)
v = XMVectorSetX(v, 10.f);
// v is now (10, 2, 3, 4)
std::cout << XMVectorGetW(v);
// Prints 4```

There are lots of structures defined in XNA Math that represent different ways vectors can be stored and packet in the memory. Many of them are very weird, but there is good reason for their existence - they are CPU equivalents for some of DXGI_FORMAT-s available in DirectX 10/11. For example, the XMUDECN4 type, representing DXGI_FORMAT_R10G10B10A2_UNORM, is defined as follows:

```union {
struct { UINT x : 10; UINT y : 10; UINT z : 10; UINT w : 2; };
UINT v;
};```

It is a 32-bit value that stores 4D vector in 10+10+10+2 bit format, where float values 0.0 .. 1.0 are mapped to integers 0..1023 for x, y, z components and 0..3 for w component (just a bit smaller precision, nothing more ;) XNA Math provides us with functions to load and store XMVECTOR from/to all these strange types, so you can just do something like this:

```XMVECTOR v = XMVectorSet(0.f, 1.f, 0.5f, 0.5f);
XMUDECN4 packedV;
XMStoreUDecN4(&packedV, v);
// packedV is 0x5ffffc00
// v2 is (0, 1, 0.499511, 0.333333)```

XNA Math contains lots of functions that we all know from old good D3DX. (The only missing thing are SH* functions for spherical harmonics, but well - I have never used them anyway :) For example, a function for normalizing 3D vector now has a signature:

`XMVECTOR XMVector3Normalize(XMVECTOR V)`

XNA Math also supports matrices. XMMATRIX type is a 4x4 matrix. It contains 4 XMVECTOR-s inside. Functions for manipulating matrices include XMMatrixTranslation, XMMatrixRotationRollPitchYaw, XMMatrixPerspectiveFovLH and so on. You can also do computations on 3D planes if only you agree to treat (x, y, z, w) vectors as (A, B, C, D) factors of plane equation. Same applies to colors - functions like XMColorAdjustSaturation also operate on XMVECTOR, but treat its components as (R, G, B, A) values.

XNA Math wouldn't be complete without some basic functions, constants and macros for scalars, such as XM_PI, XMMax, XMScalarSin. Functions with -Est suffix (from "estimated") are faster but less precise versions, e.g. XMScalarSinEst.

If you want to make full use of the performance of vector arithmetic, you should also use swizzling, permutation and comparison functionality. Swizzling is about permutation of components of a vector. It's elegant and intuitive in languages designed with gamedev in mind (like HLSL), where you just write:

`float4 v2 = v1.xyxy;`

In C++ with XNA Math, same operation can be written as:

`XMVECTOR v2 = XMVectorSwizzle(v1, 0, 1, 0, 1);`

When it comes to comparisons, functions with -R suffix return additional UINT value, which is a bitfield designed to be tested with special functions. For example:

```XMVECTOR v1 = XMVectorSet(0.0f, 1.0f, 0.5f, 0.5f);
XMVECTOR v2 = XMVectorSet(0.0f, 0.0f, 1.0f, 1.0f);
UINT cr; XMVECTORU32 res;
res.v = XMVectorGreaterR(&cr, v1, v2);
// res contains (0, 0xffffffff, 0, 0)
std::cout << XMComparisonAnyTrue(cr);
// Prints 1.```

# XNA Math Cheat Sheet

Wed
28
Apr 2010

Today I've prepared a...

XNA Math Cheatsheet (PDF)
XNA Math Cheatsheet (ODT)

(NEW) Updated links to version 1.1 for DX SDK June 2010.

# Coloring of Fractal Flame

Wed
23
Dec 2009

My next step in Fractal Flames rendering was to enhance the way I give colors to pixels on the final texture. As positions of the points iterated through the function system are discretized to a 3D matrix, but before it's processed into the final texture, I've introduced a choice of matrix format (inspired by texture formats in DirectX).

The simplest one is MATRIX_FORMAT_DENSITY, as it has only one component per cell - a density, equal to number of points that hit this particular cell. No coloring is done here, so each pixel is just white and there is only density influencing its brightness through tone mapping.

Second format is MATRIX_FORMAT_COLOR_DENSITY. It works similar to what is described in The Fractal Flame Algorithm paper. Each processed point carries a "color" with it, but this "color" is just a scalar value in range 0..1. Every possible transform (transforms are chosen in random manner according to their probabilities) pushes the color towards its "desired" value, defined for each transform. Current point color is added to the matrix cell so that final color is an average of colors of all points that hit this cell.

To visualize this scalar "color", it's nice to process it through some lookup table, gradient, palette or however you call it. I seems easy to generate or hardcode one, but I prefer to rely on existing gradients prepared by more talented people, so I've implemented reading of two gradient formats so far: gradient from a single-row texture or from a "cmap.UGR" file shipped with Apophysis (nice, free Windows application that render fractal flames according to the paper mentioned above). I'm also planning to support "ggr" files with gradients from GIMP, but it's far more sophisticated and powerful format. And finally there is the third way of coloring fractal flames that seems most intuitive to me - MATRIX_FORMAT_RGB_DENSITY. Here, full three color components are carried with points and saved to the matrix. Each transform pushes the point color towards its own "desired" RGB value. Thanks to this, each transform influences not only shape, but also color of the fractal in direct and observable way (like the scaling transform here, that shrinks points to the center, which is red). # I Render Fractal Flames

Sun
20
Dec 2009

Continuing my learning/research of fractals, today I've finally managed to render something that looks like Fractal Flame. AFAIK the most important paper on this subject is The Fractal Flame Algorithm by Scott Draves and Erik Reckase. According to this document, I've first added to my program the ability to use colors, so each point, as it has the position iteratively processed by some functions, also carries current color. Now I am able to render colourful IFS fractals:

Implementing logarithmic tone mapping and gamma correction significantly improved quality. Colors can help visualising internal structure of a fractal:

Next I've coded (almost) all functions mentioned in this paper, so now not only affine transforms are possible in my program. This allows me to mix traditional fractals with new functionality:

And finally to render something that looks like famous Fractal Flames. But it's not easy to achieve this. Modifying numeric parameters in a text file is not convenient, and even if I had an editor with real-time preview, it's not obvious what parameters give what results and whether they will degenerate to something like a single point. There are simply too many degrees of freedom :) If we had a continuous function Beauty(Params), we could use some optimization algorithms to find and render its local extrema :) Anyway, here are my first renders (more in Fractal Flames Gallery).