ಠ_ಠ: Quasicrystals as a sum of waves on a GPU

Table of Contents

2.1. Fixed wave, 0.0 angle
2.2. Fixed wave, arbitrary angle
2.3. Adding phase to waves
2.4. Summing the waves

1. Introduction

In this short post, I will describe how to implement "quasicrystals" as described in this blog post, using GPU ARB_fragment_program extension.

I will concentrate mostly on the shader implementation, rather than on setting up the OpenGL environment, passing program parameters, setting viewport and so on. Interested can find the latter explained to some detail in one of my other posts on the topic. One (not very important) difference is that since this is not a cellular automaton, there is no need to copy resulting texture back on each step of simulation.

2. Building the waves

So, to generate nice picture as shown above, we will need to:

Generate 7 rectangular "wave" pictures, for different rotation angles: 0, π/7, 2π/7, ... , 6π/7. Let's say the possible value for each pixel is in [-1.0, 1.0] range.
Sum them together, using some normalization, so the resulting values are in [0.0, 1.0] range, which we then can directly display as a color.
Change phase of the waves on each frame, to generate different quasicrystals.

Figure 2. 7 waves for different rotation angles

As an experiment, I'll try to go through the implementation step by step, adding new features on each iteration. Or you can jump to the complete shader code directly.

2.1. Fixed wave, 0.0 angle

Drawing fixed wave is simple, we basically need to implement this pseudocode:

result_color = cos(pixel.x * frequency)

Which in ARB assembly language becomes:

!!ARBfp1.0
OPTION ARB_precision_hint_nicest;

TEMP wave_freq;
MOV  wave_freq, 300.0;                            # let's define constant wave frequency for convenience

TEMP angle;                                       # define a couple of variables
TEMP cur_clr;

MUL  angle, fragment.texcoord, wave_freq;         # get current pixel coordinates to 'angle', multiplied by frequency
COS  cur_clr, angle.x;                            # calculate scalar cos(angle.x), all 4 components are the same
MOV  result.color, cur_clr;                       # draw, only [0.0, 1.0] values are drawn, but we don't care for now

END

2.2. Fixed wave, arbitrary angle

Now let's rotate the picture, the pseudocode would be:

pixel.x' = pixel.x * sin(rotate) + pixel.y * cos(rotate)
result_color = cos(pixel.x' * frequency)

Or, equivalent, after changing rotation and multiplication ordering:

pixel.x' = pixel.x * frequency
pixel.x'' = pixel.x' * sin(rotate) + pixel.y' * cos(rotate)
result_color = cos(pixel.x'')

Which gives us the following ARB code (some boilerplate code omitted this time):

TEMP rotate;
MOV  rotate, 0.448798571;                         # define scalar rotation angle, π/7
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;                          # rotate = { sin(π/7), cos(π/7), π/7, π/7 }

MUL  angle, fragment.texcoord, wave_freq;         # angle = {pixel.x * frequency, pixel.y * frequency, ..., ... }
DP4  angle, angle, rotate;                        # rotate
COS  cur_clr, angle.x;
MOV  result.color, cur_clr;

2.3. Adding phase to waves

Now, to have these waves moving, let's introduce the phase parameter, which will grow on every frame. Here is pseudocode from the previous section with added phase parameter:

pixel.x' = pixel.x * frequency
pixel.x'' = pixel.x' * sin(rotate) + pixel.y' * cos(rotate)
result_color = cos(pixel.x'' + phase)

And the corresponding ARB code:

TEMP phase;
MOV  phase, program.env[0].xxxx;                  # phase parameter passed to the shader as program environment  

MUL  angle, fragment.texcoord, wave_freq;
DP4  angle, angle, rotate;
ADD  angle, angle, phase;                         # add the phase before taking cos()
COS  cur_clr, angle.x;
MOV  result.color, cur_clr;

2.4. Summing the waves

Getting the sum of waves which will look nice on the screen is a bit tricky. Basically, we have 7 values in [-1.0, 1.0] range, and we need to sum them, bounding result nicely to [0.0, 1.0] range, so all values are mapped to some color.

Original article suggests wrapping around 1.0, while changing wrapping direction at each integer. I.e. the value of 1.3 wraps to 0.7, and 2.3 wraps to 0.3. That's exactly what we're going to do. In pseudocode:

if (floor(sum) is even)
  wrapped = fraction(sum);
else
  wrapped = 1.0 - fraction(sum);

There is no conditional jumps in ARB assembly, only conditional moves, so we'll have to execute both branches and select the needed result:

TEMP sum_clr;                                     # let's say we have sum of 7 waves in sum_clr scalar, so the value is in [-7.0, 7.0] range

ADD  sum_clr, sum_clr, 7.0;                       # sum_clr is in [0.0, 14.0] range
MUL  sum_clr, sum_clr, 0.5;                       # sum_clr is in [0.0, 7.0] range

FRC  sum_clr.x, sum_clr.aaaa;                     # sum_clr.x = fraction(sum)
SUB  sum_clr.z, 1.0, sum_clr.xxxx;                # sum_clr.z = 1.0 - fraction(sum)
FLR  sum_clr.y, sum_clr.aaaa;                     # sum_clr = { frac(sum), floor(sum), 1.0-frac(sum), ... }

# so, here we just need to choose from sum_clr.x and sum_clr.z based on whether sum_clr.y is odd or even
TEMP is_odd;
MUL  is_odd, sum_clr.yyyy, 0.5;
FRC  is_odd, is_odd;                              # is_odd != 0 when floor(sum) is odd

CMP  sum_clr, -is_odd, sum_clr.zzzz, sum_clr.xxxx;
MOV  result.color, sum_clr;

3. Complete shader program

First, for completeness, this is how external phase parameter is passed to the shader (using jogl OpenGL bindings):

// 200.0f is an arbitrary constant defining the morphing speed
gl.glProgramEnvParameter4fvARB(GL.GL_FRAGMENT_PROGRAM_ARB, 0, new float[] { _num_step / 200.0f, 0.0f, 0.0f, 0.0f }, 0);

And here is a complete ARB code, generating 7 waves and summing them:

!!ARBfp1.0
OPTION ARB_precision_hint_nicest;

TEMP wave_freq;
TEMP phase;
TEMP coord;

MOV  wave_freq, 250.0;                            # let's define constant wave frequency for convenience
MOV  phase, program.env[0].xxxx;                  # phase passed as program environment
MUL  coord, fragment.texcoord, wave_freq;         # get pixel coordinates, multiplied by frequency

TEMP angle;
TEMP cur_clr;
TEMP sum_clr;
TEMP rotate;

MOV  rotate, 0.0;                                 # 0.0 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;                          # rotate = { sin(π/7), cos(π/7), π/7, π/7 }
DP4  angle, coord, rotate;                        # angle = {x_rotated, y_rotated, z', w' }
ADD  angle, angle, phase;                         # add phase to the cos() argument
COS  cur_clr.r, angle.x;                          # store result wave's value to 1^st cur_clr's component

MOV  rotate, 0.448798571;                         # π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.g, angle.x;                          # store result wave's value to 2^nd cur_clr's component

MOV  rotate, 0.897597143;                         # 2π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.b, angle.x;                          # store result wave's value to 3^rd cur_clr's component

MOV  rotate, 1.346395713;                         # 3π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.a, angle.x;                          # store result wave's value to 4^th cur_clr's component

# first 4 components used, add them together
DP4  cur_clr, cur_clr, 1.0;
MOV  sum_clr, cur_clr;                            # sum_clr = cur_clr.r + cur_clr.g + cur_clr.b + cur_clr.a

MOV  rotate, 1.795194284;                         # 4π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.r, angle.x;                          # store result wave's value to 1^st cur_clr's component

MOV  rotate, 2.243992855;                         # 5π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.g, angle.x;                          # store result wave's value to 2^nd cur_clr's component

MOV  rotate, 2.692791426;                         # 6π/7 wave
SIN  rotate.x, rotate.x;
COS  rotate.y, rotate.y;
DP4  angle, coord, rotate;
ADD  angle, angle, phase;
COS  cur_clr.b, angle.x;                          # store result wave's value to 3^rd cur_clr's component

# remaining 3 components used, summate
DP4  cur_clr, cur_clr, {1.0, 1.0, 1.0, 0.0};
ADD  sum_clr, sum_clr, cur_clr;                   # sum_clr = sum_clr + cur_clr.r + cur_clr.g + cur_clr.b

ADD  sum_clr, sum_clr, 7.0;
MUL  sum_clr, sum_clr, 0.5;                       # transform sum_clr to [0.0, 7.0] range

# wrap the sum, so it goes from 0 to 1 and back (on even steps)
FRC  sum_clr.x, sum_clr.aaaa;
SUB  sum_clr.z, 1.0, sum_clr.xxxx;
FLR  sum_clr.y, sum_clr.aaaa;

ALIAS is_odd = angle;                             # re-use one of the variables
MUL  is_odd, sum_clr.yyyy, 0.5;
FRC  is_odd, is_odd;

CMP  sum_clr, -is_odd, sum_clr.zzzz, sum_clr.xxxx;
MOV  result.color, sum_clr;

END

Table 1. A couple of examples how all this looks in the end

Figure 3. Quasicrystal on GPU, wavelength 150 (HD)

Figure 4. Quasicrystal on GPU, wavelength 50 (HD)

4. References

"Quasicrystals as sums of waves in the plane", original blog post with the idea.
More on quasicrystals, a link from the post above.
Nice ARB_fragment_program extension specification.
All GPU programming posts from this blog.

ಠ_ಠ

Jan 1, 2012

Quasicrystals as a sum of waves on a GPU