This week I was exploring fractals - specifically computing fractals using Julia Set
quadratic equation. I considered using SYCL DPC++ for implementation, which I've been
working on for few weeks now. I was interested in using image interface of DPC++, which
just same as buffer interface creates an abstraction over raw memory address based
data read/ write ( read pointer arithmetic ), while providing some image specific functionality for data interpretation,
reading & writing. I find buffer & image interfaces easier to work with rather than raw pointers,
allocated using Unified Shared Memory, due to their elegance in design; accessor based
data read/ write on device/ host, providing implicit data movement/ synchronization/ visibility, which in-turn
lets me implicitly define task graph dependencies. So it's good, nice, powerful abstraction.
I'll now spend some time explaining how to sequentially compute Julia Set fractal.
z = z ** 2 + c, where z, c ∈ ℂ && c is a constant
I sequentially apply aforementioned equation on each pixel of target fractal image as long as certain conditions satisfy. Let us assume, I'm working with an image of dimension 4 x 4. This is how each cell ( read pixel ) is addressed in 4 x 4 image. Notice the origin of reference frame in use. These indices can easily be generated using two nested for loops.
(uint, uint) A[4][4];
for i in 0..4 {
for j in 0..4 {
A[i][j] = (i, j);
}
}
But for computing fractal, I've to reinterpret indices in a way where origin of reference frame should move to center of image. And those points on 2D cartesian plane to be interpreted as complex numbers ( read a + ib, where a, b are components of 2D point along X & Y-axis respectively ). It can be easily generated using meshgrid method.
$ python3
>>> w = np.linspace(-2, 1, 4); w // along X-axis
array([-2., -1., 0., 1.])
>>> h = np.linspace(2, -1, 4); h // along Y-axis
array([ 2., 1., 0., -1.])
>>> a, b = np.meshgrid(w, h)
>>> a // visualise 4 X-axes placed along rows of matrix
array([[-2., -1., 0., 1.],
[-2., -1., 0., 1.],
[-2., -1., 0., 1.],
[-2., -1., 0., 1.]])
>>> b // visualise 4 Y-axes placed along columns of matrix
array([[ 2., 2., 2., 2.],
[ 1., 1., 1., 1.],
[ 0., 0., 0., 0.],
[-1., -1., -1., -1.]])
>>> a + b*1j // complex numbers, required for fractal computation
array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j],
[-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j],
[-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j],
[-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j]])
Notice along X-axis points are in range [-2, 2), while along Y-axis they're in range [2, -2). For scaling them to S, I first take them to (-1, 1) interval & then multiply by S. Now as I've these complex numbers in matrix form properly scaled, I'd like to go to actual fractal computation part.
>>> W, H = 4, 4 // image dimensions
>>> factor = ((W/2) ** 2 + (H/2) ** 2) ** 0.5; factor // √((w/2) ** 2 + (h/2) ** 2)
2.8284271247461903
>>> a / factor // all values in interval (-1, 1)
array([[-0.70710678, -0.35355339, 0. , 0.35355339],
[-0.70710678, -0.35355339, 0. , 0.35355339],
[-0.70710678, -0.35355339, 0. , 0.35355339],
[-0.70710678, -0.35355339, 0. , 0.35355339]])
>>> # do same for b i.e. `b / factor`
>>> a / factor * S + (b/factor * S) * 1j // S = scale, here 2.f
array([[-1.41421356+1.41421356j, -0.70710678+1.41421356j, 0.+1.41421356j, 0.70710678+1.41421356j],
[-1.41421356+0.70710678j, -0.70710678+0.70710678j, 0.+0.70710678j, 0.70710678+0.70710678j],
[-1.41421356+0.j , -0.70710678+0.j , 0. +0.j, 0.70710678+0.j ],
[-1.41421356-0.70710678j, -0.70710678-0.70710678j, 0.-0.70710678j, 0.70710678-0.70710678j]])
Quadratic iteration based Julia Set computation can be encapsulated in following function. As long as z.real()2 + z.imag()2 < S2, I keep applying z = z2 + c. There is another bound on number of iterations that can be performed.
const uint MAX_ITR = 1 << 10; // 1024
const float S = 2.f; // scale
func length(complex z) -> float {
return (z.real() ** 2 + z.imag() ** 2) ** .5f; // z = a + ib; √(a**2 + b**2)
}
func quadratic_iteration(complex z, c) {}
uint i = 0;
for ; i < MAX_ITR; i++ {
z = z ** 2 + c;
// z ** 2 = (a+ib) * (a+ib) = (a**2-b**2) + 2ab*i
if (length(z) > S) {
break;
}
}
return i;
}
I choose S ( = scale ) such that S2 - S >= abs(c), where c is a constant complex number. Say, I choose to work with c = (-0.8+0.156j), in that case setting S = 2.f, satisfies aforementioned constraint.
>>> S = 2.; S // scale
2.0
>>> c = -0.8 + 0.156j; c // constant
(-0.8+0.156j)
>>> S ** 2 - S >= abs(c)
True ✅
From quadratic_iteration(...), I obtain #-of evaluations it takes for z
to determine whether it's in set with radius S or not. Depending upon #-of evaluations for certain
z, I assign different color to respective image pixel. #-of iterations denote how many of it
is required by z to get outside of Julia Set with radius S, which is alongside capped
at maximum number of iterations ( = MAX_ITR ). After every evaluation of z = z2 + c, I check
whether distance of updated z from center (0.+0j) is still less than radius of Julia Set defined on
complex plane. This way I keep increasing z by (z2 + c) in each iteration & in some time z will reach
some value in complex plane, which goes outside of set. I take that certain iteration number & use it to represent
color.
Recently I discovered a nice way to assign
beautiful colors to each pixel of image using cosine function. It simply requires me to
evaluate color(t) = a + b * cos(2 * pi * (c * t + d)), where I can variate
a, b, c, d ( 3-element floating point tuples ) to get different palettes.
In this case t = quadratic_iteration(...) / MAX_ITR.
I recommend reading more about cosine palette based color generation technique
here.
Today I'm using following method for colorizing each pixel of fractal.
[float; 4] colorize(uint itr) {
[float; 3] d = {0.5, 0.5, 0.5};
[float; 3] e = {0.5, 0.5, 0.5};
[float; 3] f = {2.0, 1.0, 0.0};
[float; 3] g = {0.3, 0.2, 0.4};
// last 1.f dictates α channel value of RGBA image
[float; 4] pixel_val = {d + e * cos(6.28318f * (f * ((float)itr / (float)(MAX_ITR)) + g)), 1.f};
return pixel_val;
}
Now I've one sequential program which computes Julia Set Fractal of dimension 512x512, by evaluating quadratic equation z = z2 + c. Following is the computed fractal with c = (-0.8 + 0.156j).
It's good time to move to parallel implementation of Julia Set computation on heterogeneous accelerators ( read CPU, GPGPU, FPGA etc. ).
Being a good candidate for SIMD model, Julia Set computation fits well in data parallel programming domain, where each pixel
of resulting image is processed independently & RGBA pixel value is written in respective memory location.
The only challenge I see --- how to compute complex number (a + ib) located in pixel location (r, c) ?
In previous iteration, I was using meshgrid method to do so, this time I've to implement it manually.
Each kernel invocation processes single ( unique ) pixel of image, so certain invocation has knowledge
of which pixel it is responsible of processing. Let me assume certain kernel invocation knows, it has to process
image pixel (r, c). Again my interest is to move reference frame to center of image.
r ∈ [0, H) && H = image height
c ∈ [0, W) && W = image width
Notice each (x, y) point on 2D cartesian plane, is seperated by distance 1, so I can compute what should be (x, y) in current pixel location (r, c), with following formula.
w = h = 4 // image dimension
func f(r, c) {
step = 1.;
// a + s * k
return -(w/2) + step * c, (h/2) - step * r;
}
>>> assert f(0, 0) == (-2., 2.) ✔ // top left
>>> assert f(0, 3) == (1., 2.) ✔ // top right
>>> assert f(3, 0) == (-2., -1.) ✔ // bottom left
>>> assert f(3, 3) == (1., -1.) ✔ // bottom right
Rest is same as I've already done while implementing sequential algorithm for fractal computation. I need to scale (x, y) point by factor S, so I rewrite aforementioned function as below, which I make use of for converting (r, c) image pixel to complex number (a+ib) which is to be operated on during fractal computation i.e. checking whether that complex number is part of Julia Set or not.
// check https://gist.github.com/itzmeanjan/0a13447c634afbc34f64d1f22751011d#file-julia_set-cpp-L16-L23
// for more info
[float; 2] convert_to_complex_number(uint r, c) {
float step = 1.;
float factor = ((w/2) ** 2 + (h/2) ** 2) ** .5;
// S = scale factor
return {( ( -(w/2) + step * c ) / factor ) * S, ( ( (h/2) - step * r ) / factor ) * S}; // complex number (a+ib) as (a, b)
}
I play around with complex constat c, with c = -0.123 + 0.745j, it generates following fractal.
Here's a dendrite shaped fractal, with c = 0.f + 1.0j.
With c = -0.75 + -0.2j, beautiful galaxy shaped fractal is generated.
Good thing about cosine based color palette is just changing to following starts generating fractals with very different color.
{
// check https://gist.github.com/itzmeanjan/0a13447c634afbc34f64d1f22751011d#file-julia_set-cpp-L38-L41
// for more info
[float; 3] d, e, f = /* unchanged */ ;
[float; 3] g = {0.5, 0.2, 0.25}; // updated
// ...
}
I keep whole DPC++ Julia Set computation implementation
here
for future reference. In coming days, I plan to compute Julia Set with sine/ cosine function, as this time I did with quadratic
function. I'm also interested in Newton Fractal, which I'll consider implementing next week in DPC++.
Have a great time !