## Moire Patterns with PyCUDA

Before the semester ended, two of my good friends – Dave and Adrian – procrastinated a bit and developed a very cool Matlab program. The program generates Moire patterns by creating arrangements of concentring circles and having them decay across the screen into one another. It’s a very clever program, and I’ve decided to one-up them a little bit. I just finished prototyping my own version of the program, but there’s a catch – rather than use Matlab, I decided to tackle the problem with a different set of tools.

For starters, I’m addicted to Python. I think it’s the go-to language for a myriad of tasks as diverse as scirpting, software prototyping, and even applicaiton development. Although the climate models I study are coded mostly in Fortran, I’ve developed my own set of tools for processing their outpat data through Python – mostly through the help of RPy, NumPy, and SciPy (with the help of NCO for working with netCDF files, of course, although SciPy has a great I/O functionality for them). I decided to port Dave and Adrian’s code to Python.

However, there were issues. For starters, I wasn’t sure how to do 2D graphics with Python. I first tried TKinter but I couldn’t optimize the drawing very well. So, I settled with PyGame, which turned out to be *very *easy to use! Even with PyGame, though, the prototype program ran very sluggishly. So, I decided to up the ante and try something a bit esoteric – PyCUDA. PyCUDA is a convenient module for utilizing CUDA code within a Python script or program. For those of you who don’t know, CUDA is nVidia’s solution for bringing GPGPU to the masses – that is, general purpose programming on graphics processing units. In a nutshell, modern GPU’s are massively parallel constrcuts; they have hundreds of cores and have theoretical computing powerof over a teraflop! They are perfectly suited for SIMD situations.

The code for my prototype program follows after the break, and serves as a practice attempt at utilizing PyCUDA. Now, this isn’t real hardcore CUDA programming here – I’m using very convenient abstractions available in PyCUDA because my translation algorithm is so simple. As far as I can tell, there isn’t even any real speed increase in this situation (although I plan on modifying things a bit here in a second and doing some time trials). However, this is a neat project for seeing how PyCUDA can be effortlessly thrown into a script or small program. Again, this is only a prototype program to explore PyGame and PyCUDA; hopefully in the coming days I’ll flesh out the same simulation that Dave and Adrian ran, and then I’ll post both their Matlab code and my Python code for comparison and experimentation.

Some of this is cut off on the margin, so I published it from my Google Docs here.

################################################ ## moire_v1.py - Daniel Rothenberg ## ## 5/27/2009 ## ## ## ## Performs a simple translation of ## ## concentric circles to create a Moire ## ## interference pattern. ## ################################################ ## Import Modules import math, numpy, pygame import pycuda.autoinit import pycuda.driver as drv import pycuda.gpuarray as gpuarray ## Function Definitions def translateCircle(circle, delta_X): """ Given a 'circle', utlizes the GPU to rapidly translate that circle delta_X units horizontally. circle = a circle delta_X = number of units to translate """ num_circles = len(circle) num_points = len(circle[0][0]) for c in range(0, num_circles): coords = numpy.array(circle[c][0]).astype(numpy.float32) c_gpu = gpuarray.to_gpu(coords) c_trans = (delta_X + c_gpu).get() circle[c][0] = c_trans def drawCircle(screen, circle, color): """ Draws a circle on a given pygame screen. screen = a pygame scren circle = a circle color = a tuple representing a color """ num_circles = len(circle) num_points = len(circle[0][0]) for c in range(0, num_circles): for i in range(0, num_points): pygame.draw.line(screen, color, (circle[c][0][i-1], circle[c][1][i-1]), (circle[c][0][i], circle[c][1][i])) def createCircle(c_x, c_y, radii, theta): """ Creates a circle from a given set of parameters. c_x = x coordinate of center c_y = y coordinate of center radii = radii of the concentric circles theta = list of all the theta values for calc vectors """ num_circles = len(radii) num_points = len(theta) circle = [] for c in range(0, num_circles): x = [] y = [] for i in range(0, num_points): radius = radii[c] x.append(c_x + radius*math.cos(theta[i])) y.append(c_y + radius*math.sin(theta[i])) circle.append([x, y]) return circle ## Main Program if __name__ == '__main__': ## Constants ## Note - all these should be ints, except for max_radius, ## which should be a float. num_circles = 100 max_radius = 300.0 num_points = 150 DIM_X = 800 DIM_Y = 800 dX = 1 timesteps = 250 BLACK = (0,0,0) WHITE = (255,255,255) ## Generate two lists - one to hold the radii of each of ## the concentric circles in our larger circles, and one to ## hold the steps of theta used in calculating the vector ## coordinates of those circles R = [(max_radius/num_circles)*i for i in range(0, num_circles+1)] t = [(2.0*math.pi/num_points)*i for i in range(0, num_points+1)] ## Create our circles. You can customize c_x and c_y ## depending on the dimensions of our screen circle1 = createCircle(300, 400, R, t) circle2 = createCircle(500, 400, R, t) ## Initialize a pygame screen with our custom dimensions screen = pygame.display.set_mode((DIM_X, DIM_Y)) running = True ## Main animation loop. Need to keep track of the timesteps ## whic have passed to auto exit the program at its end. counter = 0 while running: event = pygame.event.poll() counter = counter+1 ## Clear the screen and paint the current iteration ## of our circles. screen.fill(WHITE) drawCircle(screen, circle1, BLACK) drawCircle(screen, circle2, BLACK) pygame.display.flip() ## Behind the scenes, translate the circles. translateCircle(circle1, dX) translateCircle(circle2, -dX) ## Check exit conditions if ((event.type == pygame.QUIT) or (counter > timesteps)): running = False

A screenshot of the script in action: (once I figure out how to have PyGame save a gif, I’ll post the entire animation)

################################################

## moire_v1.py – Daniel Rothenberg ##

## 5/27/2009 ##

## ##

## Performs a simple translation of ##

## concentric circles to create a Moire ##

## interference pattern. ##

################################################## Import Modules

import math, numpy, pygame

import pycuda.autoinit

import pycuda.driver as drv

import pycuda.gpuarray as gpuarray## Function Definitions

def translateCircle(circle, delta_X):

“””

Given a ‘circle’, utlizes the GPU to rapidly translate

that circle delta_X units horizontally.circle = a circle

delta_X = number of units to translate

“””

num_circles = len(circle)

num_points = len(circle[0][0])

for c in range(0, num_circles):

coords = numpy.array(circle[c][0]).astype(numpy.float32)

c_gpu = gpuarray.to_gpu(coords)

c_trans = (delta_X + c_gpu).get()

circle[c][0] = c_transdef drawCircle(screen, circle, color):

“””

Draws a circle on a given pygame screen.screen = a pygame scren

circle = a circle

color = a tuple representing a color

“””

num_circles = len(circle)

num_points = len(circle[0][0])

for c in range(0, num_circles):

for i in range(0, num_points):

pygame.draw.line(screen, color,

(circle[c][0][i-1], circle[c][1][i-1]),

(circle[c][0][i], circle[c][1][i]))def createCircle(c_x, c_y, radii, theta):

“””

Creates a circle from a given set of parameters.c_x = x coordinate of center

c_y = y coordinate of center

radii = radii of the concentric circles

theta = list of all the theta values for calc vectors

“””

num_circles = len(radii)

num_points = len(theta)

circle = []

for c in range(0, num_circles):

x = []

y = []

for i in range(0, num_points):

radius = radii[c]

x.append(c_x + radius*math.cos(theta[i]))

y.append(c_y + radius*math.sin(theta[i]))

circle.append([x, y])

return circle## Main Program

if __name__ == ‘__main__’:## Constants

## Note – all these should be ints, except for max_radius,

## which should be a float.

num_circles = 100

max_radius = 300.0

num_points = 150

DIM_X = 800

DIM_Y = 800

dX = 1

timesteps = 250

BLACK = (0,0,0)

WHITE = (255,255,255)## Generate two lists – one to hold the radii of each of

## the concentric circles in our larger circles, and one to

## hold the steps of theta used in calculating the vector

## coordinates of those circles

R = [(max_radius/num_circles)*i for i in range(0, num_circles+1)]

t = [(2.0*math.pi/num_points)*i for i in range(0, num_points+1)]## Create our circles. You can customize c_x and c_y

## depending on the dimensions of our screen

circle1 = createCircle(300, 400, R, t)

circle2 = createCircle(500, 400, R, t)## Initialize a pygame screen with our custom dimensions

screen = pygame.display.set_mode((DIM_X, DIM_Y))

running = True## Main animation loop. Need to keep track of the timesteps

## whic have passed to auto exit the program at its end.

counter = 0

while running:

event = pygame.event.poll()

counter = counter+1## Clear the screen and paint the current iteration

## of our circles.

screen.fill(WHITE)

drawCircle(screen, circle1, BLACK)

drawCircle(screen, circle2, BLACK)

pygame.display.flip()## Behind the scenes, translate the circles.

translateCircle(circle1, dX)

translateCircle(circle2, -dX)## Check exit conditions

if ((event.type == pygame.QUIT) or

(counter > timesteps)):

running = False

Hi

Can I have the MATLAB code for moire pattern please?

Thank you

Shadi said this on September 23, 2010 at 2:06 am |

You can definitely see your enthusiasm within the article you write.

The world hopes for even more passionate writers such

as you who are not afraid to say how they believe.

At all times go after your heart.

Jule Z. Clausing said this on September 27, 2014 at 7:48 am |