On the suggestion of a friend, I rendered this animation:
This quasicrystal is full of emergent patterns, but it can be described in a simple way. Imagine that every point in the plane is shaded according to the cosine of its y coordinate. The result would look like this:
Now we can rotate this image to get other waves, like these:
Each frame of the animation is a summation of such waves at evenlyspaced rotations. The animation occurs as each wave moves forward.
I recommend viewing it up close, and then from a few feet back. There are different patterns at each spatial scale.
The code
To render this animation I wrote a Haskell program, using the Repa array library. For my purposes, the advantages of Repa are:
 Immutable arrays, supporting clean, expressive code
 A fast implementation, including automatic parallelization
 Easy output to image files, via
repadevil
Here is a simplified (but complete!) program, which renders a single still image.
import Data.Array.Repa ( Array, DIM2, DIM3, Z(..), (:.)(..) )
import qualified Data.Array.Repa as R
import qualified Data.Array.Repa.IO.DevIL as D
import Data.Word ( Word8 )
import Data.Fixed ( divMod' )
For clarity, we define a few type synonyms:
type R = Float
type R2 = (R, R)
type Angle = R
We'll convert pixel indices to coordinates in the real plane, with origin at the image center. We have to decide how many pixels to draw, and how much of the plane to show.
pixels :: Int
pixels = 800
scale :: R
scale = 128
Repa's array indices are "snoc lists" of the form (Z :. x :. y)
. By contrast, our planar coordinates are conventional tuples.
point :: DIM2 > R2
point = \(Z :. x :. y) > (adj x, adj y) where
adj n = scale * ((2 * fromIntegral n / denom)  1)
denom = fromIntegral pixels  1
A single wave is a cosine depending on x and y coordinates in some proportion, determined by the wave's orientation angle.
wave :: Angle > R2 > R
wave th = f where
(cth, sth) = (cos th, sin th)
f (x,y) = (cos (cth*x + sth*y) + 1) / 2
To combine several functions, we sum their outputs, and wrap to produce a result between 0 and 1. As n
increases, (wrap n)
will rise to 1, fall back to 0, rise again, and so on. sequence
converts a list of functions to a function returning a list, using the monad instance for ((>) r)
.
combine :: [R2 > R] > (R2 > R)
combine xs = wrap . sum . sequence xs where
wrap n = case divMod' n 1 of
(k, v)  odd k > 1v
 otherwise > v
To draw the quasicrystal, we combine waves at 7 angles evenly spaced between 0 and π.
angles :: Int > [Angle]
angles n = take n $ enumFromThen 0 (pi / fromIntegral n)
quasicrystal :: DIM2 > R
quasicrystal = combine (map wave (angles 7)) . point
We convert an array of floatingpoint values to an image in two steps. First, we map floats in [0,1] to bytes in [0,255]. Then we copy this to every color channel. The result is a 3dimensional array, indexed by (row, column, channel). repadevil
takes such an array and outputs a PNG image file.
toImage :: Array DIM2 R > Array DIM3 Word8
toImage arr = R.traverse arr8 (:. 4) chans where
arr8 = R.map (floor . (*255) . min 1 . max 0) arr
chans _ (Z :. _ :. _ :. 3) = 255  alpha channel
chans a (Z :. x :. y :. _) = a (Z :. x :. y)
main :: IO ()
main = do
let arr = R.fromFunction (Z :. pixels :. pixels) quasicrystal
D.runIL $ D.writeImage "out.png" (toImage arr)
Running it
Repa's array operations automatically run in parallel. We just need to enable GHC's threaded runtime.
$ ghc O2 rtsopts threaded quasicrystal.lhs
$ ./quasicrystal +RTS N
$ xview out.png
And it looks like this:
Note that repadevil
silently refuses to overwrite an existing file, so you may need to rm out.png
first.
On my 6core machine, this parallel code ran in 3.72 seconds of wallclock time, at a CPU utilization of 474%. The same code compiled without threaded
took 14.20 seconds, so the net efficiency of parallelization is 382%. This is a good result; what's better is how little work it required on my part. Cutting a mere 10 seconds from a single run is not a big deal. But it starts to matter when rendering many frames of animation, and trying out variations on the algorithm.
As a side note, switching from Float
to Double
increased the run time by about 30%. I suspect this is due to increased demand for memory bandwidth and cache space.
You can grab the Literate Haskell source and try it out on your own machine. This is my first Repa program ever, so I'd much appreciate feedback on improving the code.
Be sure to check out Michael Rule's work on animating quasicrystals.
How can you be doing both
ReplyDelete> import qualified Data.Array.Repa as R
and
> type R = Float
at the same time? My Haskell is very rusty, but I can't see how that would work...
Should be fine... "R" will refer to the type, since "R" by itself as a module qualifier makes no sense.
ReplyDeletePorted to HTML5 canvas: http://www.jasondavies.com/animatedquasicrystals/
ReplyDeletenice idea! reminds me of the old plasmas in the demoscene...
ReplyDeleteported to webgl via shadertoy  realtime in your browser!
https://gist.github.com/f448ba84e94c61ab5924
enjoy.
The resultant image isn't a true structure of a quasicrystal, is it? Doesn't the pattern repeat beyond the bounds of the image?
ReplyDeleteIt has 7fold symmetry about the origin, so it can't be periodic.
ReplyDeleteThese patterns never truly repeat, try it for yourself. You can render as much of the crystal as you like, and drive yourself crazy looking for alignment, but all you'll ever find are local regions that approximately align.
ReplyDeleteAlso ported to WebGL via PhiloGL. Can toggle fullscreen :)
ReplyDeletehttp://senchalabs.github.com/philogl/PhiloGL/examples/quasicrystal/
Does not work on nvidia because you define "float v" twice, just rename it to some other value and it works fine.
ReplyDelete// paste this into http://www.iquilezles.org/apps/shadertoy/ for pictures!
ReplyDelete// this is a gpu/webgl version of http://mainisusuallyafunction.blogspot.com/2011/10/quasicrystalsassumsofwavesinplane.html
// 'implementation' by @mmalex for http://news.ycombinator.com/item?id=3153835
// fixed up by colski
#ifdef GL_ES
precision highp float;
#endif
uniform float time;
uniform vec2 resolution;
void main(void)
{
vec2 p = gl_FragCoord.xy / resolution.xy;
float t=time*5.0; // change this for more speed
vec2 up=vec2(256.0,0.0); // change this for finer/coarser stripes
const float c = 0.90096886790241, s = 0.43388373911755; // sin and cos of pi / 7
const vec2 rx=vec2(c, s), ry=vec2(s, c);
float sum = 0.0;
for (int i = 0; i < 7; ++i)
{
sum += cos(dot(p,up)+t);
up = vec2(dot(up,rx), dot(up,ry));
}
float a1=(cos(sum+2.0)+1.0) * 0.5; // the important wrap around term. play with this for different effects
gl_FragColor = vec4(a1,a1,a1,1.0);
}
Ported to Processing: http://openprocessing.org/visuals/?visualID=43954
ReplyDeleteWhat about an R version?
ReplyDeleteframe < function(N,Scale,Sym,FS,time) {
x=matrix(rep((1:(2*N+1))N,(2*N+1)),nrow=(2*N+1))*Scale
y=t(x)
S=sin(pi/Sym)
C=cos(pi/Sym)
total=matrix(rep(0,(2*N+1)^2),nrow=(2*N+1))
Upx=FS; Upy=0;
for (i in 1:Sym) {
total=total+cos(x*Upx+y*Upy+time)
oldUpx=Upx
Upx=Upx*C+Upy*S
Upy=oldUpx*S+Upy*C
}
(cos(total+2.0)+1.0)*0.5
}
library(png)
for (i in 0:49) {
fname=paste("quasi",i,".png",sep="")
writePNG(frame(512,1/32,7,8,2*pi*i/50),fname)
}
I found some code on the HaskellWiki for writing PNGs. I used it here:
ReplyDeletehttps://github.com/dagit/haray/blob/master/src/Data/PNG.hs
https://github.com/dagit/haray/blob/master/src/Data/Bitmap.hs
To use my code you create a BMP, write to it (it's just a mutable array under the hood), and then write that to PNG when you're done. No need for 3rd party libraries that are hard to install as it's just pure Haskell.
The code is BSD3 so feel free to use it if you like it.
Nice !! Same logic can be applied for CG water reflection using multiple level of normal maps moving in different direction and velocities.
ReplyDeleteInteractive, color mapped version implemented in Processing.
ReplyDeletehttp://tryptamine.net/quasicrystal/
wrong broken or removed?
DeleteGreat work. Can's say that I care for the Haskell though. Here's some Matlab code to recreate the image at the end of your post (which, by the way, I think is rotated by 90 degrees):
ReplyDeleteN=800; scale=1/2.5; angles=7; t=0;
x=(0.5*(1N):0.5*(N1))*scale;
x=x(ones(1,N),:);
y=x';
rx=cos(0:pi/angles:pi);
ry=sin(0:pi/angles:pi);
qc=cos(x+t);
for i=2:angles
qc=qc+cos(x*rx(i)+y*ry(i)+t);
end
imagesc(0.5*(cos(qc+2)+1));
colormap(gcf,gray);
set(gcf,'Units','pixels','Position',[0 0 N N]);
set(gca,'Position',[0 0 1 1]);
The 800x800 image takes 0.2 sec. on my old 2.4 GHz Intel Core 2 Duo MacBook Pro. The R version from @fermin above was helpful. However, it used repeated multiplication of rotation matrices by each other in order to get the sequence of rotated waves, which can lead to numeric issues. Here's a more general Matlab function that outputs an arbitrarilysized movie (takes under 3 sec. for a 512x512 25frame movie like the animated GIF above):
function quasicrystal(N,scale,angles,frames)
narginchk(3,4);
if ~isvector(N)  isempty(N)  ~isnumeric(N)  length(N) > 2
error('First input must be one or two element numeric vector.');
end
if N(1) < 1  mod(N(1),1)  (length(N) > 2 && (N(2) < 1  mod(N(2),1)))
error('First input must be vector >= 1.');
end
if ~isscalar(scale)  isempty(scale)  ~isnumeric(scale)  scale <0
error('Second input must be positive scalar numeric value.');
end
if ~isscalar(angles)  isempty(angles)  ~isnumeric(angles)
error('Third input must be scalar numeric value.');
end
if angles < 1  mod(angles,1)
error('Third input must be integer >= 1.');
end
if nargin == 3
frames=1;
else
if ~isscalar(frames)  isempty(frames)  ~isnumeric(frames)
error('Fourth input must be scalar numeric value.');
end
if frames < 1  mod(frames,1)
error('Fourth input must be integer >= 1.');
end
end
x=(0.5*(1N(1)):0.5*(N(1)1))*scale;
if length(N) == 2 && N(1) ~= N(2)
x=x(ones(1,N(2)),:);
y=(0.5*(1N(2)):0.5*(N(2)1))'*scale;
y=y(:,ones(1,N(1)));
else
x=x(ones(1,N(1)),:);
y=x';
end
rx=cos(0:pi/angles:pi);
ry=sin(0:pi/angles:pi);
vid=VideoWriter('qc.avi');
open(vid);
for t=2*pi/frames*(0:frames1)
qc=cos(x+t);
for i=2:angles
qc=qc+cos(x*rx(i)+y*ry(i)+t);
end
writeVideo(vid,0.5*(cos(qc+2)+1));
end
close(vid);
Here's an updated and improved version of my Matlab code that can output movies, animated GIFs, or PNGs in Cartesian or logpolar coordinates (inspired by Michael Rule) with the ability to adjust contrast and apply a colormap:
ReplyDeletehttp://biorobots.case.edu/personnel/adh/quasicrystal/
There are also some example animated GIFs and AVIs. The code will be particularly useful to anyone working in R or numPy/sciPy.
Also @madan's comment above about water reflection made me think of when I tried outputting in just polar coordinates, as opposed to logpolar. When the scale is 1, the resulting animation is a very pleasing circular ripple pattern. It's especially nice if you have a full color colormap.
Cheers,  Andy
Joining in the choir of reimplementations, here it is in GLSL (WebGL): http://glsl.heroku.com/e#5384.1
ReplyDeleteIt lets you muck about with the code interactively, so that's neat :)
Hi, you make mind blowing ideas and a spectacular article hereanimated video production
ReplyDeleteIs there a java equivalent for this code? :(
ReplyDeletecheck the rikrd's version. Processing is Java based, so this shouldn't be hard to port.
Deletehttp://openprocessing.org/visuals/?visualID=43954
see also this Gerasimov fractal
ReplyDeleteHey Keegan, I absolutely love this quasicrystal image, would you be interested in it being used as a vinyl LP album cover? I'm in the band Moths & Locusts from Vancouver Island, Canada, we're currently working on our 2nd full album. Pls drop me a line at david@vinylrecordguru.com. Hope it's OK to post this here, looking forward to hearing back from you! ~Dave
ReplyDeleteThanks for posting! I really like what you've acquired here; You should keep it up forever! Best of luck
ReplyDeletepuin afvoeren
Nice article, thanks for the information. It's very complete information. I will bookmark for next reference
ReplyDeletejaring futsal  jaring golf  jaring pengaman proyek 
jaring pengaman bangunan  jaring pengaman gedung
http://www.jualjaring.blogspot.com/
http://www.agenjaring.blogspot.com/
http://www.pancasamuderasafetynet.blogspot.com/
http://www.tokojaring.blogspot.com/
http://www.pusatjaring.blogspot.com/
http://jualjaringpengaman.blogspot.com/
https://pancasamudera.wordpress.com/
https://pasangjaringfutsal.wordpress.com/
https://jualtambangmurah.wordpress.com/
https://tokojaring.wordpress.com/
https://jualjaringfutsal.wordpress.com/
https://jaringfutsal.wordpress.com/
Thank you for the hard work you have made in writing this post. In the future I am hopeful the same best work from you as well. Actually your creative writing abilities has motivated me
ReplyDeleteYou can't really say what is beautiful about a place, but the image of the place will remain vividly with you.
ReplyDeletelevelشركة تسليك مجارى بالرياض
ReplyDeleteشركة تنظيف بالرياض
شركة تنظيف شقق بالرياض
شركة تنظيف منازل بالرياض
شركة تنظيف خزنات بالرياض
شركة مكافحة حشرات بالرياض
شركة رش مبيدات بالرياض
شركة تخزين اثاث بالرياض
شركة تنظيف مجالس بالرياض
شركة تنظيف فلل بالرياض
great tutorial :)
ReplyDeletebonus bagging login
thanks for share this with us :)
ReplyDeleteDoes Vo Genesis Work
thanks for your times
ReplyDeleteSalehoo Review 2016
ReplyDelete;شركة مكافحة حشرات بجدة
شركة تنظيف منازل بجدة
شركة مكافحة النمل الابيض بالرياض شركة تنظيف فلل بالدمام
شركة تنظيف موكيت بالدمام
شركة تنظيف فلل بالجبيل
شركة تنظيف بالقطيف شركة تنظيف فلل بالرياض
شركة تنظيف منازل بجدة
شركة تنظيف خزانات بالدمام
تنظيف الستائر بسهوله

ReplyDelete

شركة مكافحة حشرات بالدمام
شركة مكافحة النمل الابيض بالدمام
شركة مكافحة حشرات بالاحساء
شركة تنظيف منازل بالرياض
شركة كشف تسربات المياه بالرياض upgrading
شركة مكافحة حشرات بالرياض

شركة مكافحة حشرات بالاحساء
شركة مكافحة النمل الابيض بالاحساء
شركة مكافحة حشرات بالدمام
شركة مكافحة النمل الابيض بجدة
التخلص من الحشرات