Table of Contents
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.
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.
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.
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
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;
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;
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.0frac(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;
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; # reuse 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

"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.
No comments:
Post a Comment