PythonPlasma

plasma.py

# Python PIL (pillow) code to generate frame images out of a plasma algorithm
# Author: 2015 Thomas O. Robinson, cyranix@cyranix.net 
# Open Source, no particular rights reserved. Would love to hear from anyone who reads, uses or improves it

# Certain research was done using http://www.bidouille.org/prog/plasma and http://lodev.org/cgtutor/plasma.html as references.

# I started this project because of a need to generate some random, psychedelicish grayscale animations for use in another project I am working on, and old school plasma screensavers seemed like a logical solution. Since there really wasn't any good tutorials on the net for implementing this, and literally nothing for python (which I am only using because the libraries for my other project only have bindings for C++ or python), which has some pretty poor documentation for working with images and graphics and pixel data and such, I set out on my own, and this is the result of about 3 days of inconsistent work. I've documented it the best I could for now, so I hope this either proves useful to someone else, or else someone improves on my own formulas and gives me a better method of achieving my images. Either way, I donate this work to the public domain. Enjoy!

import PIL, math
import numpy as np
from scipy import ndimage as nd # Might not need this, but included just incase

# This command just creates a new PIL image, 640x480 pixels with black background, and simultaneously converts it into a float32 array in the form [y][x][r, g, b], which I do for a lot of reasons, primarily because I just really hate python, and PIL documentation sucks, and really, its just easier to work with arrays (in fact, I usually treat them as a matrix and use numpy/scipy to manipulate them, but such sourcery (yes, I did that on purpose, sorry) was not strictly necessary for this exercise. 
img = np.float32(PIL.Image.new("RGB", (640, 400), 0))

# I think the normal iteration here in math conventions goes from 0 to 2pi (lol radians), but I wanted my first sinusoid to start on the left side of the image rather than the top when I started this project, so thats why start and end go from -pi to pi. Technically, you can use any multiple of 2pi from any starting point here to get that many complete cycles of any of the lissajous curves used below (the path which the center of the moving sines follow)
start = -math.pi
end = math.pi

# Frames per second, and number of seconds to run for.
fps = 30
s = 30

# These arrays are storing the "key" values for making some lissajous curves a few lines down from here... The normal formula for a Lissajous curve is:
# x = A sin(a*t + Delta), y = B sin(b*t)
# But the real concept here is just that you have x and y sinusoids not in phase, so a simple parametric equation using a sin and cos (since their phases are off by fraction of a period of each other):
# x = A cos(kx*t), y = B sin(ky*t)
# does the job very elegantly (you can invert the sin/cos or swap them if you like and the effect will still work). I didn't need multiple phases, so we just assume '1' for A and B (you can add multipliers if you like), and now kx becomes the number of "peaks" we will see on the sides, and ky becomes the number of peaks we see across the top and bottom. For further visualization, the first iteration, kx=1, ky=2, gives us a kind of flattened out infinity symbol, whereas the next iteration, kx=2, ky=1, is like a figure 8. kx=3, ky=3 would be like a hash/pound symbol shape. Oh yeah, I just made this an array for simplicity of calling below.
kx = [1, 2, 2, 3]
ky = [2, 1, 3, 4]

i = 0 # Incremential counter for numbering save files

# Okay, so this creates a full loop, from -pi to pi (remember above?), and does so in increments of (pi - -pi)/(30fps * 30s), so in total (assuming you didn't modify these values above), this will generate 900 frames per complete cycle of one lissajous curve.
for t in [round(step, 2) for step in np.arange(start, end, (end - start)/(fps * s))]:
    # Okay, so later on, I want to move 4 different circular sine waves around my lissajous curves. These two arrays each loop are storing the x and y locations for that spot at that "time" (t) on the curve (so if I were to place a point on the image at each of these locations, it would draw the shape of the curve, but I 'm not interested in drawing the curve, I just want to use it as a "path" from which to draw out later). Also, because the array measures from (0,0) in the upper left corner, and I want to use all 4 quadrants of a "graph" for the points, I move (0,0) to the center of our "virtual" graph by calculating the cos/sin (which gives us a value between -1 and 1) and multiplying that by half the width or height of our graph (so figure, if its 640 wide, -1 * (640 / 2) == -320, 1 * (640 / 2) == 320, giving us a range of -320 to 320, works the same way for height), and I add that value to half the width or height (thats actually the first parameter, not the second one, but I had to explain it in this order to make sense), and finally subtract 1 at the end, because we have to calculate from (0, 0) to start, not (1, 1).
    xc = [
      ((w / 2) + (math.cos(kx[0] * t) * (w / 2))) - 1,
      ((w / 2) + (math.cos(kx[1] * t) * (w / 2))) - 1,
      ((w / 2) + (math.cos(kx[2] * t) * (w / 2))) - 1,
      ((w / 2) + (math.cos(kx[3] * t) * (w / 2))) - 1
    ]
    yc = [
      ((h / 2) + (math.sin(ky[0] * t) * (h / 2))) - 1,
      ((h / 2) + (math.sin(ky[1] * t) * (h / 2))) - 1,
      ((h / 2) + (math.sin(ky[2] * t) * (h / 2))) - 1,
      ((h / 2) + (math.sin(ky[3] * t) * (h / 2))) - 1
    ]

    # r was originally for radius, although its really more of a modifier, technically, I should rename it 'p' because its the period of the function used below, but you can call it fontandiddyhopper if you want, its just a name. For reference, the idea on this current equation is that it will always be at LEAST 32 (not exact pixels, but rough estimate in a single linear direction), plus a value which increases as a percentage of our pi based time function above (in laymans terms, as our loop t approaches pi, our period approaches the value 64, or whatever you change the two coefficients to, I used 16 and 32 for testing). You can make it a static value if you don't want the period to grow/shrink
    r = 32 + (((math.pi + t) / (2 * math.pi)) * 32)

    # Finally, a basic loop which just draws from left to right, top to bottom. 
    for y in xrange(h):
        for x in xrange(w):
            # So, our np.float32 array is in the form of y, x, [R, G, B]. I just wanted grayscale, so I affect RGB all together, hence the [0:3]. Then basically, at whatever point we're at on our lissajous curves in this loop, I am calculating the sum of radial sinewaves (formula: sin(sqrt(x^2 + y^2) / period) is the simplest way to do this based on a radian circle formula of x^2=y^2, although I might change this to an ellipse instead for some more interesting geometry), you'll note here again that for each sine, I add 1 at the beginning so we get a value between 0 and 2, which is why in the second to last term, I divide by 8 (4 sines * 2.0 maximum value) to get my percent value, which I then multiply by 255 to get our shade of color, which is what we're assigning to that "pixel" in our array. Its a lot simpler than it sounds. 
            img[y][x][0:3] = (
                (
                    # A note on the math here, the reason for the x - xc, y - yc factors is because we need to calculate the value of the pixel we want based on its offset from the "virtual" centers we calculated above, not where 'x' and 'y' are currently drawing, which is based off our drawing loop above, so without this parameter, we'd be calculating our offset from (0, 0). Also, in the final term, I multiplied the periods by slightly different amounts because at the last second this was easier for me to fine tune the size/shape of each radial sine period than modifying 'r' above.
                    (1 + math.sin(math.sqrt((x - xc[0])**2 + (y - yc[0])**2) / (2.0 * r))) +
                    (1 + math.sin(math.sqrt((x - xc[1])**2 + (y - yc[1])**2) / (1.8 * r))) +
                    (1 + math.sin(math.sqrt((x - xc[2])**2 + (y - yc[2])**2) / (1.8 * r))) +
                    (1 + math.sin(math.sqrt((x - xc[3])**2 + (y - yc[3])**2) / (1.5 * r)))
                ) / float(8)
            ) * 255

    # So that was all the guts of the code above. So now, at each iteration of 't', we're going to save a frame. Depending on what shell you run this from, you will probably want to specify a directory to save these images in. To save space, I convert our array back to PIL image and save it all in one line. Because python. And it would be stupid to divide this out into several commands, but to make this as easy to understand as possible, the PIL.Image.fromarray() can't read float32 values, so the first thing I do, using np.clip, to clip all of our pixel values to between 0 and 255 (they should already be between 0 and 255, otherwise numbers outside of that bound are automatically "clipped" at that value, which by the way corresponds to completely black and completely white, if I failed to mention the meaning for those numbers above), and then round those values off to whole numbers and convert them to 8 bit integer values using np.uint8 (which is what fromarray() wants to read). Of course, the %03d and %i terms number the images for us in standard pr1ntf format :). Alternatively, you can use .show() instead of .save() if you just want to output the image to a window and view it, however, keep in mind, this is going to open a window for every frame of output, which is 900 frames at 30fps for 30s ;). Its useful if you just want to see how a single frame of output might look at a certain point though (I'll leave it to you to write the debug code for all that though). I output a print statement just so we can know that our code is working in the background, and then increment our frame number and then the t loop iterates and runs again!
    PIL.Image.fromarray(np.uint8(np.clip(alpha_layer, 0, 255))).save("lissajoussin%03d.jpg" %i, fmt="jpeg")
    print "Saved: lisajoussin%05d.jpg" %i
    i += 1

# Once all the images are saved, you can animate them easily using ffmpeg (I will leave it to you to read the man page for that, its not that complicated ;). If I have more time later, I'll post updated code that does some color tricks, or you could try to implement some of the tricks listed in my references up top! Have fun, let me know if you have any questions about it. :)
# --TR 16 Sep 2015
# cyranix@cyranix.net