Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@omegahm
Last active April 28, 2021 01:44
Show Gist options
  • Star 12 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save omegahm/e823a68c201406d32a94 to your computer and use it in GitHub Desktop.
Save omegahm/e823a68c201406d32a94 to your computer and use it in GitHub Desktop.
Quasicrystals in Python
# -*- coding: utf-8 -*-
import bohrium as np
from math import pi
import matplotlib.pyplot as plt
import matplotlib.colors as colors
fig = plt.figure()
plt.xticks([])
plt.yticks([])
k = 5 # number of plane waves
stripes = 37 # number of stripes per wave
N = 512 # image size in pixels
ite = 30 # iterations
phases = np.arange(0, 2*pi, 2*pi/ite)
image = np.empty((N, N))
d = np.arange(-N/2, N/2, dtype=np.float64)
xv, yv = np.meshgrid(d, d)
theta = np.arctan2(yv, xv)
r = np.log(np.sqrt(xv*xv + yv*yv))
r[np.isinf(r) == True] = 0
tcos = theta * np.cos(np.arange(0, pi, pi/k))[:, np.newaxis, np.newaxis]
rsin = r * np.sin(np.arange(0, pi, pi/k))[:, np.newaxis, np.newaxis]
inner = (tcos - rsin) * stripes
cinner = np.cos(inner)
sinner = np.sin(inner)
i = 0
for phase in phases:
image[:] = np.sum(cinner * np.cos(phase) - sinner * np.sin(phase), axis=0) + k
plt.imshow(image.copy2numpy(), cmap="RdGy")
fig.savefig("quasi-{:03d}.png".format(i), bbox_inches='tight', pad_inches=0)
i += 1
@omegahm
Copy link
Author

omegahm commented Mar 7, 2016

Convert the 30 pngs to an animated gif, by running:

$ convert -delay 8 *.png quasi.gif

@Zulko
Copy link

Zulko commented Mar 8, 2016

Your pipeline is a little heavy (plot with matplotlib, export to many pngs, then animate separately with convert). As a note you can create GIFs directly from the numpy arrays using MoviePy, it is simpler and much faster:

import numpy as np
import matplotlib.cm as cm
import moviepy.editor as mpe

k       = 5    # number of plane waves
stripes = 37   # number of stripes per wave
N       = 512  # image size in pixels
duration = 2 # the final gif's duration in seconds

d = np.arange(-N/2, N/2, dtype=np.float64)
xv, yv = np.meshgrid(d, d)
theta  = np.arctan2(yv, xv)
r      = np.log(np.sqrt(xv*xv + yv*yv))
r[np.isinf(r) == True] = 0

tcos   = theta * np.cos(np.arange(0, np.pi, np.pi/k))[:, np.newaxis, np.newaxis]
rsin   = r * np.sin(np.arange(0, np.pi, np.pi/k))[:, np.newaxis, np.newaxis]
inner  = (tcos - rsin) * stripes

cinner = np.cos(inner)
sinner = np.sin(inner)

def make_frame(t):
    phase = 2*np.pi*t / duration 
    array = np.sum(cinner * np.cos(phase) - sinner * np.sin(phase), axis=0) + k
    normalized_array = (array - array.min()) / (array.max() - array.min()) 
    colored_array = cm.RdGy(normalized_array, bytes=True)[:,:,:3] 
    return colored_array

clip = mpe.VideoClip(make_frame, duration=duration)
clip.write_gif("animation.gif", fps=10)

This script takes under one second in total and produces this GIF.

@omegahm
Copy link
Author

omegahm commented Mar 8, 2016

Cool. I'm still learning a lot of Python.

Some other guy recommended to dump the raw NumPy arrays and using ffmpeg to convert it to gif, which also seemed like a long journey for the data.

Your solution is quite nice, however I am not sure I like the fact that you have to state gif-duration instead of defining how many frames there should be in total (which in this case will end up being 20, right?)

@Zulko
Copy link

Zulko commented Mar 9, 2016

It's actually quite practical. The duration basically sets the speed of your gif (since it is looping) so it makes sense. As long as you are testing the animation, you make the rendering very quick by setting the fps to 5, and when you are happy with the result you set the fps to 30 and get a final, much fluider gif.

If you want to control the number of frames, is see these solutions:

  • Set the fps automatically to n_frames / duration
  • Use MoviePy's ImageSequenceClip, which lets you make a clip from a series of images or numpy arrays.
  • Use Imageio, a library simpler than MoviePy (and used by MoviePy) which let you work frame by frame.

You may also be interested by these other examples.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment