# Tag: math

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

Uwaga! Informacje na tej stronie mają ponad 5 lat. Nadal je udostępniam, ale prawdopodobnie nie odzwierciedlają one mojej aktualnej wiedzy ani przekonań.

# Further Experiments with Fractals

Sun
13
Dec 2009

Today I woke up with an idea of inventing my own simple IFS fractal with four pieces sheared or rotated to form a hole in the center. I must admin that effects are not impressive :) What I've learned is that it's difficult to imagine how overall fractal will generally look like when you design affine transforms for a single step that will be processed recursively. But anyway I'll show my results, with single step on the left and final render on the right:

Shear transform which unexpectedly formed round shapes:  Rotation transform with boundaries exceeding -1..1 range:  Rotation transform bound in range (-1..1, -1..1):  # Random Choice with Given Probabilities

Fri
11
Dec 2009

When we want to randomly choose one option from some possible items, we usually do something as simple as (rand() % numberOfItems). But what if we want to vary probabilities of these items? Here is an algorithm for this:

As input data, we have number of items and an array of their probabilities. Each probability is a nonnegative real number. I like to allow these numbers to be as big as I want, knowing that an item with probability 2 (or 200) is just two times more probable than the item with probability 1 (or 100). I will normalize them later.

```// == Preprocessing step ==

// Input:
std::vector<float> ItemProbabilities;
// Output:
std::vector<float> ProcessedProbabilities;

// Algorithm:

size_t count = ItemProbabilities.size();
assert(count > 0);

// First pass: summation
float sum = 0.f;
for (size_t i = count; i--; ) // count-1, count-2, ..., 0
{
float p = ItemProbabilities[i];
assert(p >= 0.f);
// Optional: Delete items with zero probability
if (p == 0.f)
{
Items.RemoveByIndex(i);
// That's how we delete an element from std::vector by index :P
ItemProbabilities.erase(ItemProbabilities.begin()+i);
}
sum += p;
}
assert(sum > 0.f);
float sumInv = 1.f / sum;

// Second pass: building preprocessed array
ProcessedProbabilities.resize(count);
ProcessedProbabilities = ItemProbabilities * sumInv;
for (size_t i = 1; i < count; i++)
ProcessedProbabilities[i] =
ProcessedProbabilities[i-1] + ItemProbabilities[i] * sumInv;```

The algorithm above does two things: It prepares ProcessedProbabilities array so that each probability is the item's probability plus the sum of probabilities of all previous items. It also normalizes them to range 0..1. As the result, the output array forms an ascending sequence, where the difference to the previous element is proportional to the item's probability in respect to full 0..1 range. For example:

```Input: 100, 100, 200
Temporary data: count=3, sum=400, sumInv=0.0025
Output: 0.25, 0.5, 1```

Now as we have these preprocessed data, generating random choices with certain probabilities is fast and simple. We just generate a random real number in range 0..1 and choose the first item for which ProcessedProbabilities[i] is greater than this. Additional optimization would be to use binary search here.

```// == Random choice ==

// Input:
std::vector<float> ProcessedProbabilities;
// Output:
size_t choice;

// Algorithm:
choice = 0;
float randNum = RandFloat(); // 0..1
while (choice < ProcessedProbabilities.size()-1
&& ProcessedProbabilities[choice] < randNum)
{
choice++;
}```

# Rendering Mandelbrot and Julia Fractal

Thu
10
Dec 2009

I've been playing recently with Mandelbrot and Julia sets. These are fractals rendered in the completely different way that my last experiments with IFS (Iterated Function Systems). In IFS I had a set of points with given (x,y) positions and I've been transforming their positions iteratively by some functions before I rendered them. I'll blog more about IFS another time. Rendering Mandelbrot and Julia fractals is more like writing a pixel shader - for every pixel on the screen with its (x,y) coordinates we do some computations to calculate final color. I didn't extend my framework to render such things yet, but instead I've been doing my experiments in EvalDraw (where I can freely pan and zoom the view) and AMD RenderMonkey (where I can write pixel shaders executed on GPU).

Complex Numbers for Game Developers :) Fractal is a mathematical concept and so there is lots of hardcore math involved in the foundation of these beautiful objects. But fortunately we don't have to understand all that stuff to be able to render them. Mandelbrot and Julia fractals are created using complex numbers, but for us they can be thought of just as float2/D3DXVECTOR2 2-component (x,y) vectors with some operations defined:

• Addition: (x1, y1) + (x2, y2) = (x1+x2, y1+y2)
• Muliplication: (x1, y1) * (x2, y2) = (x1*x2 - y1*y2, y1*x2 + x1*y2)
• Square: (x, y)^2 = (x, y) * (x, y) = (x*x - y*y, 2*x*y)

That's actually all we need to render Mandelbrot and Julia fractals. The general formula for both of them is: z = z^2 + c where z and c are complex numbers, c is a constant and z is recalculated many times through this function. It was a mystery for me for a long time how to apply this formula. Now it turned out that when we calculate it many times for each pixel with staring z_start = (0, 0) and c dependant on pixel position in range (-2..2, -2..2) and draw all the pixels for which z is bounded and do not go to infinity, we get the shape of the Mandelbrot fractal. Of course we are not going to iterate this function inifinite times in our computers like mathematicians do in their heads, so instead we do it several times and treat all pixels for which length(z_final) < 2 as being "inside" the shape. Here is a draft made in EvalDraw. Unfortunately we meet some limitations of floating point numbers here, like limited precision and infinites and that's why the space outside is white. To render Julia set we do something similar, but this time we use pixel position as starting z and set c to some constant, like (-0.1, 0.7). Different constants make different images.

A question arises about how to render a fractal - how to map a complex number to pixel color? The algorithm I like the most is called Escape Time Algorith and it goes like this: recalculate z for several iterations, but only while the length of z=(x,y) vector is less than 2 (that is zx*zx + zy*zy < 2*2). If the loop reached the iteration count limit, this pixel is inside the fractal - draw in in black. If not, then color of this pixel depends on number of iterations done.

Here are my implementations of these fractals in EvalDraw and RenderMonkey. They are very simple so the whole source code is visible on the screenshots. You can also download the code here: Fractals_for_RenderMonkey.rfx, Mandelbrot_for_EvalDraw_1.kc, Julia_for_EvalDraw_1.kc, Multibrot_for_EvalDraw_1.kc. You can find more screenshots in my fractals gallery.

Some variations of these fractals have also been discovered. For example when we do absolute value of z components before every iteration in Mandelbrot fractal, we get an interesting shape called Burning Ship.

Alternatively, if we raise z to any real power instead of 2, we make Mandelbrot fractal "unfolding", having more and more "bulbs" as the exponent grows. It is called Multibrot set. Raising a complex number to power with any real exponent is a little bit more complicated, as it has to use polar form (HLSL-like pseudocode here):

```float2 ComplexPow(float2 v, float n)
{
float len = sqrt(x*x + y*y); // Standard 2D vector length
float n_fi = n * atan2(y, x);
return pow(len, n) * float2( cos(n_fi), sin(n_fi) ); // scalar * vector
}```

# New Fractal Renders

Sun
06
Dec 2009

Continuing playing with fractals, I've improved my experiment framework so now I can render any IFS (Iterated Function System) fractals described by a set of affine transformations with subpixel quality. Here are my today renders (these with colors are processed with GIMP). You can find more in my Gallery / Fractals. Where do I get information from? There is a great website with ready formulas for each of these fractals: Classic Iterated Function Systems by Larry Riddle from Agnes Scott College.

# Properties of Pascal Triangle

Thu
12
Nov 2009

Pascal Triangle is a mathematical object that looks like triangle with numbers arranged the way like bricks in the wall. Each next row has one more number, ones on both sides and every inner number is the sum of two numbers above it. It can span infinitely. Despite simple algorithm this triangle has some interesting properties. First, if you draw only odd numbers, you get a fractal - Sierpinski Triangle. Second, graph of numbers in each row resembles Gaussian distribution function. The lower row you take from the triangle (containing more numbers), the more precise graph you get. I've made these images with Octave - a free Matlab alternative. I like it and I think concepts behind this Matlab programming language can be quite useful for doing computations, but doing graphics is painfully slow. Whole seconds to draw several dozens of texts or rectangles? WTF?!

# Calculating Linear and Quadratic Equation Coefficients

Sun
30
Aug 2009

Some time ago I've written about formulas to calculate coefficients of linear and quadratic equation [pl] when having 2 or 3 given points (x,y). Yesterday I suddenly noticed that my code needs to do it every frame so I need functions to calculate these coefficients. Below you can find the code of my functions. They are templates, so they work with many types including float, double, as well as vectors of any dimmension and other types which have addition, subtraction and assignment operators, as well as multiplication and divistion by float scalar. Read full entry | Comments | #rendering #math #algorithms

# Reflect and Refract Functions

Mon
24
Aug 2009

Every game developer must have some knowledge about maths and especially geometry. Math library for game developers always contains structures like vector, matrix and functions like vector normalization or matrix inversion. But high level shader languages, designed to do geometrical computations, have many more interesting intrinsic functions. Some of more complex ones are reflect and refract, which calculate vector reflection and refraction in relation to the normal of a surface. How to implement them in C++? Reflection works according the simple rule that "angle of incidence equals angle of reflection". Although there are no angles in the formula - all the necessary work is done by dot product.

```void Reflect(VEC3 &out, const VEC3 &incidentVec, const VEC3 &normal)
{
out = incidentVec - 2.f * Dot(incidentVec, normal) * normal;
}```

For reflect function, the formula is given in th HLSL documentation, but for refract function it is not. It took me some time to find the implementation, but finally I've found it in the GLSL Specification (thanks KriS!). So here is the code converted to C++. Eta is the refraction index (e.g. 1.5f).

```inline void Refract(
VEC3 &out, const VEC3 &incidentVec, const VEC3 &normal, float eta)
{
float N_dot_I = Dot(normal, incidentVec);
float k = 1.f - eta * eta * (1.f - N_dot_I * N_dot_I);
if (k < 0.f)
out = VEC3(0.f, 0.f, 0.f);
else
out = eta * incidentVec - (eta * N_dot_I + sqrtf(k)) * normal;
}```

These functions are part of my library: CommonLib.

# Three Ways to Calculate Mean and Variance

Mon
20
Jul 2009

There is a free e-book on DSP (Digital Signal Processing) called The Scientist and Engineer's Guide to Digital Signal Processing (by Steven W. Smith, Ph.D.). As I have a holiday now, I've started reading it and I've found some interesting information even in introductory chapters. For example, there are three algorithms to calculate mean and variance. Let's say we have some data in a vector of bytes and we want to calculate its statistics.

```std::vector<unsigned char> Bytes;
// Fill Bytes...
uint N = Bytes.size();```

```float Mean = 0.f, Variance = 0.f;
for (uint i = 0; i < N; i++)
Mean += (float)Bytes[i];
Mean /= (float)N;
for (uint i = 0; i < N; i++)
Variance += ((float)Bytes[i] - Mean) * ((float)Bytes[i] - Mean);
Variance /= (float)(N-1);```

Now the incremental one, which can return current mean and variance at any time while you can also post next piece of data. All we need to implement it is to keep track of current sample number, sum of samples and sum of squared samples. I've created template class for that, which can be parametrized by float, double, int or other numeric types.

```template <typename T>
class IncrementalMeanAndVarianceCalc
{
public:
IncrementalMeanAndVarianceCalc() : m_Sum(T()), m_SqSum(T()), m_Count(0) { }
void Clear() { m_Sum = m_SqSum = T(); m_Count = 0; }

{
m_Sum += v;
m_SqSum += v*v;
m_Count++;
}
void Add(const T *values, uint valCount)
{
for (uint i = 0; i < valCount; i++)
{
m_Sum += values[i];
m_SqSum += values[i]*values[i];
}
m_Count += valCount;
}

bool IsEmpty() { return m_Count == 0; }
uint GetCount() { return m_Count; }
T GetSum() { return m_Sum; }
T GetSqSum() { return m_SqSum; }

T GetMean()
{
assert(!IsEmpty());
return m_Sum / m_Count;
}
T GetVariance(bool varianceBiased = true)
{
if (varianceBiased)
return (m_SqSum - m_Sum*m_Sum/m_Count) / (m_Count-1);
else
return (m_SqSum - m_Sum*m_Sum/m_Count) / m_Count;
}
void GetMeanAndVariance(T &outMean, T &outVariance, bool varianceBiased = true)
{
assert(!IsEmpty());
outMean = m_Sum / m_Count;
if (varianceBiased)
outVariance = (m_SqSum - m_Sum*m_Sum/m_Count) / (m_Count - 1);
else
outVariance = (m_SqSum - m_Sum*m_Sum/m_Count) / m_Count;
}

private:
T m_Sum, m_SqSum;
uint m_Count;
};```

Finally, there is an algorithm to calculate mean and variance from histogram. It is very efficient method. If we have only 256 possible values for each sample, we don't have to do binning and calculating histogram is very simple:

```uint Histogram = { 0 };
for (uint i = 0; i < N; i++)
Histogram[Bytes[i]]++;```

Now, the mean and variance can be calculated this way:

```float Mean = 0.f, Variance = 0.f;
for (uint i = 0; i < 256; i++)
Mean += (float)i * Histogram[i];
Mean /= (float)N;
for (uint i = 0; i < 256; i++)
Variance += (i - Mean)*(i - Mean) * Histogram[i];
Variance /= N - 1;```