Tuesday, December 11, 2007

Creating fractal images with C# 3.0 features and Parallel Extensions

In this post I'm going to show a small example of using the Parallel Extensions Library to improve a program that renders a escape time fractal image.

A couple of months ago I wrote about a small program that creates an image with the escape time algorithm using C# 3.0 features . For this post I'm going to change the implementation to use the Parallel Extensions Library.

Along the the library documentation and MSDN articles, the blog entries from the Parallel Extensions Team blog provided a nice introduction and guidance . The original post was inspired by Luke Hoban's blog entry A Ray Tracer in C#3.0 .

The escape time algorithm is very simple and also considered a embarrassingly parallelproblem. For previous posts I did a similar experiment of using Fortress parallel for loops to render the Mandelbrot Set fractal.

The code for the original program was 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);
}

}



We can change this code in several ways to take advantage of the library. Here I'm going to present two of them.

Write one big LINQ expression

The code for the generation of pixel data is coded in a single LINQ expression.

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).AsParallel()
from xi in Enumerable.Range(0, imageWidth)
let mappedX = xF(xi)
let mappedY = yF(yi)
let p0 = new Complex(xF(xi), yF(yi))
let function = functionConstructor(p0)
select new
{
x = xi,
y = yi,
xD = mappedX,
yD = mappedY,
i = ApplyFunction(function, p0)
.TakeWhile(
(x, j) => j < maxIteration && x.NormSquared() < 4.0)
.Count()
})
{
HandlePixel(p.x, p.y, p.i);
}


Elements from the body of the foreach loop were moved to the LINQ expression by using the let keyword.

The AsParallel extension method called does the magic of distributing the work for calculating the lines of the image. A nice discussion on were to put this call can be found in Parallelizing a query with multiple “from” clauses and Chunk partitioning vs range partitioning in PLINQ.

Although this problem is much more simpler than Ray tracing, this alternative tries to follow the same approach as: Taking LINQ to Objects to Extremes: A fully LINQified RayTracer.



Add a Parallel.For loop

The other alternative is to replace one of the from elements with a Parallel.For loop.


Func<double, double> xF = MathUtils.InterpFunc(0, minX, imageWidth, maxX);
Func<double, double> yF = MathUtils.InterpFunc(0, minY, imageHeight, maxY);

Parallel.For(0, imageHeight , delegate(int yi)
{
foreach (var p in from xi in Enumerable.Range(0, imageWidth)
let mappedX = xF(xi)
let mappedY = yF(yi)
let p0 = new Complex(xF(xi), yF(yi))
let function = functionConstructor(p0)
select new
{
x = xi,
y = yi,
xD = mappedX,
yD = mappedY,
i = ApplyFunction(function, p0)
.TakeWhile(
(x, j) => j < maxIteration && x.NormSquared() < 4.0)
.Count()
})
{
lock (this)
{
HandlePixel(p.x, p.y, p.i);
}
}
});



A lock needed to be added because HandlePixel calls Bitmap.SetPixel which needs to be protected.

Results

I tested this code by generating a 2000 by 2000 image of a small section of the Mandelbrot fractal located between (-1.1752491998171,0.223337905807042) and (-1.17342021033379,0.225166895290352) with a escape time of 512. A reduced image of this location looks like this:

Calculated Fractal Image

By using the library the program took about 40% less of the time of original version on my dual core machine. The performance for the two approaches presented above was very similar.

The nicest thing about this experiment is that the code required just a couple of modifications in order to take advantage of the other CPU.

Another approach for rendering this fractal using Parallel extensions is from Jon Skeet's Coding Blog : LINQ to Silliness: Generating a Mandelbrot with parallel potential and A cautionary parallel tale: ordering isn't simple.

Code for this post can be found here.