This post is inspired by the very nice blog post A Ray Tracer in C#3.0.
The Escape time algorithm is very simple, but depending on its input it can generate very complex images.
Although not mandatory, the first thing to create is a class from complex numbers. The complex number class contains the following members:
public class Complex
{
public Complex(double real, double img)
public Complex Multiply(Complex c)
public Complex Add(Complex c)
public double Norm()
public double NormSquared()
public static Complex operator *(Complex c1,Complex c2)
public static Complex operator +(Complex c1, Complex c2)
public double Real
public double Img
}
Now we need something to convert from the image coordinate system to the real coordinate system. To do this the following function was created:
public static Func<double, double> InterpFunc(double x1,double y1,double x2,double y2)
{
double m = (y2 - y1) / (x2 - x1);
double b = y1 - (m * x1);
return (x => m * x + b);
}
The
InterpFunc method creates a function that maps values from one coordinate system to the other.The code that creates the image is the following:
void CreateFractal(double minX, double minY,double maxX, double maxY,int imageWidth, int imageHeight)
{
Func<double, double> xF = MathUtils.InterpFunc(0, minX, imageWidth, maxX);
Func<double, double> yF = MathUtils.InterpFunc(0, minY, imageHeight, maxY);
foreach (var p in from yi in Enumerable.Range(0, imageHeight)
from xi in Enumerable.Range(0, imageWidth)
select new
{
x = xi,
y = yi,
xD = xF(xi),
yD = yF(yi)
})
{
Complex p0 = new Complex(p.xD, p.yD);
Func<Complex, Complex> function = functionConstructor(p0);
int i = ApplyFunction(function, p0)
.TakeWhile(
(x, j) => j < maxIteration && x.NormSquared() < 4.0)
.Count();
HandlePixel(p.x, p.y, i);
}
}
Basically
select expression generates the coordinates for all the points in the image . The functionConstructor function is used to create the function that will be applied to generate the fractal. One example of this functions is the one used for the Mandelbrot fractal f(x) = x^2 + p0 . The escape time is calculated by counting the number of elements in the generated sequence of function applications before the value escaped.The
ApplyFunction method is interesting since is the one that creates the sequence of recursive function applications required for this algorithm. The method looks like this:
IEnumerable<Complex> ApplyFunction(Func<Complex, Complex> function,Complex initial)
{
Complex last = initial;
while (true)
{
last = function(last);
yield return last;
}
}
By running this algorithm using the Mandelbrot formula:
Func<Complex,Func<Complex,Complex>> fc =
p0 => (c => c.Multiply(c).Add(p0));
EscapeTimeFractal p = new EscapeTimeFractal(300, 300, fc);
p.MinX = -0.03;
p.MinY = 0.68;
p.MaxX = 0.03;
p.MaxY = 0.62;
p.Run();
The generated image is the following:
By running it using the Szegedi Butterfly 1 formula:
Func<Complex,Func<Complex,Complex>> fc =
p0 =>
(c => (new Complex(
((c.Img * c.Img) - Math.Sqrt(Math.Abs(c.Real))),
((c.Real * c.Real) - Math.Sqrt(Math.Abs(c.Img)))))
+ p0);
EscapeTimeFractal p = new EscapeTimeFractal(300, 300,fc);
p.MinX = -2;
p.MinY = 2;
p.MaxX =2;
p.MaxY = -2;
p.MaxIteration = 127;
p.OutputFileName = @"c:\temp\output.bmp";
p.Run();
The generated image is the following:
The code for this experiment can be found here.
