Martin Laprise

June 18, 2010

githubGraph + ubigraph server = pure awesomeness

Example of a GitHub user’s social graph. This undirected graph of users which mutually follow each other, contains 7243 nodes and 14484 edges. The video shows the convergence of the Force-based algorithm used to visualize the graph.

April 15, 2010

Weekend project: compute a 68 gigapixels fractal image: Part II

So, as promised, here is the source code for the gigapixels fractal project: gmpano.tar.gz. The archive contain the following files:

mandelbrot.cpp: The actual C++ code for the computation
numpy.i and PyMandelbrot.i: The SWIG descriptions files used for the function wrapper (C++ to Python)
gigaMandelbrot.py: This one contain the Python class used for the image tile rendering
gmPano.py: The final script that putting it all together and write the final image in the Deep Zoom Format. You will need to write the xml file by hand . You can use HDView SL or HDView for viewing the image.

You will also need SWIG and the following python package: numpy, matplotlib, python parallel.

Here is a example of compilation process on a 64-bits linux machine:
swig -c++ -python PyMandelbrot.i
g++ -c -O3 -fPIC -fomit-frame-pointer -funroll-loops -fstrict-aliasing -m64 mandelbrot.cpp -I/usr/include/python2.6
g++ -c -O3 -fPIC -fomit-frame-pointer -funroll-loops -fstrict-aliasing -m64 PyMandelbrot_wrap.cxx -I/usr/include/python2.6
g++ -shared mandelbrot.o PyMandelbrot_wrap.o -o _PyMandelbrot.so -lm -lstdc+

March 27, 2010

Weekend project: compute a 68 gigapixels fractal image

Here is a project I made in order to test HDView SL. I’m not a big fan of the Flash/Silverlight thing, but I have to admit that the Seadragon zooming technique include in Silverlight is pretty neat. Put it in Fullscreen and have fun !

Quick facts:
- The final image size is around 3.2 GB.
- A physical print @ 300 dpi would be around 72.75 feet x 72.75 feet ! Which is equivalent to a 8-story building
- It took around 8 hrs to compute on a 4 cores CPU @ 2.4 Ghz.
- The main loop of the code was done in C.
- The final multithreaded code was done in Python with the help of Parallel Python
- The function wrapper (C -> Python) was done with SWIG and numpy.i.
- I’ll publish the code in the next posts

March 2, 2010

Fractal & Python: Part II

After the Lorenz attractor, here is a little code for exploring the Mandelbrot Set. Of course, you could compute the same thing in CUDA with a 1000X speedup… but it wouldn’t be pythonish as this !

from numpy import *
import matplotlib.pylab as pl
 
maxIteration = 128
z_min = -2-1j
z_max = 1+1j
 
# Set the image size here
imageSize = [512,512]
image = zeros(imageSize,int)
xValue = linspace(z_min.real, z_max.real, imageSize[0])
yValue = linspace(z_min.imag, z_max.imag, imageSize[1])
 
for i in arange(imageSize[0]):
	for k in arange(imageSize[1]):
		iteration = 0
		x = xValue[i]
		y = yValue[k]
		while (x*x+y*y < (2*2)) & (iteration < maxIteration):
			xtemp = x*x - y*y + xValue[i]
 			y = 2*x*y + yValue[k]
			x = xtemp
			iteration = iteration + 1
		if iteration == maxIteration:
			image[i,k] = 0
		else:
			image[i,k] = 255-iteration
 
pl.imshow(image)
pl.prism()
pl.show()

February 28, 2010

Visualizing the Lorenz attractor in 3D with VTK

I just finished my nth reading of the James Gleick’s classic Chaos: Making a New Science and I realized that the famous Lorenz attractor should be much more exciting to explore in 3D. So here is a very simple Python code that resolve the nonlinear differential equation describing Lorenz oscillator and show the solution in a 3D phase space with VTK.

$\frac{dx}{dt} = \sigma (y - x)$
$\frac{dy}{dt} = x (\rho - z) - y$
$\frac{dz}{dt} = x y - \beta z$

The coolest think about the Lorenz attractor is the simplicity of it. You can retrieve it simply by resolving those three coupled equations which can be done in 10 line of code ! The rest is just the VTK kitchen sink.

Required packages: python-vtk, python-scipy, python-matplotlib

#!/usr/bin/env python
 
import vtk
from scipy import *
from scipy import integrate
from pylab import *
from vtk.util.colors import tomato, banana
 
nbrPoints = 5000
 
def Lorenz(w, t, S, R, B):
    x, y, z = w
    return array([S*(y-x), R*x-y-x*z, x*y-B*z])
 
w0 = array([0.0, 1.0, 0.0])
time = linspace(0.0, 100.0, nbrPoints)
S = 10.0; R = 28.0; B = 8.0/3.0
 
# Conditions initiales
y0=array([-7.5, -3.6, 30.0])
 
# On resous le tout avec odeint de scipy qui est
# lui-meme un wrap de ODEPACK en fortran
trajectory = integrate.odeint(Lorenz, w0, time, args=(S, R, B))
 
# This will be used later to get random numbers.
math = vtk.vtkMath()
 
# Total number of points.
numberOfInputPoints = nbrPoints
 
# One spline for each direction.
aSplineX = vtk.vtkCardinalSpline()
aSplineY = vtk.vtkCardinalSpline()
aSplineZ = vtk.vtkCardinalSpline()
 
inputPoints = vtk.vtkPoints()
for i in range(0, numberOfInputPoints):
    x = trajectory[i,0]
    y = trajectory[i,1]
    z = trajectory[i,2]
    aSplineX.AddPoint(i, trajectory[i,0])
    aSplineY.AddPoint(i, trajectory[i,1])
    aSplineZ.AddPoint(i, trajectory[i,2])
    inputPoints.InsertPoint(i, x, y, z)
 
# The following section will create glyphs for the pivot points
# in order to make the effect of the spline more clear.
 
# Create a polydata to be glyphed.
inputData = vtk.vtkPolyData()
inputData.SetPoints(inputPoints)
 
# Use sphere as glyph source.
balls = vtk.vtkSphereSource()
balls.SetRadius(.01)
balls.SetPhiResolution(10)
balls.SetThetaResolution(10)
 
glyphPoints = vtk.vtkGlyph3D()
glyphPoints.SetInput(inputData)
glyphPoints.SetSource(balls.GetOutput())
 
glyphMapper = vtk.vtkPolyDataMapper()
glyphMapper.SetInputConnection(glyphPoints.GetOutputPort())
 
glyph = vtk.vtkActor()
glyph.SetMapper(glyphMapper)
glyph.GetProperty().SetDiffuseColor(tomato)
glyph.GetProperty().SetSpecular(.3)
glyph.GetProperty().SetSpecularPower(30)
 
# Generate the polyline for the spline.
points = vtk.vtkPoints()
profileData = vtk.vtkPolyData()
 
# Number of points on the spline
numberOfOutputPoints = nbrPoints*10
 
# Interpolate x, y and z by using the three spline filters and
# create new points
for i in range(0, numberOfOutputPoints):
    t = (numberOfInputPoints-1.0)/(numberOfOutputPoints-1.0)*i
    points.InsertPoint(i, aSplineX.Evaluate(t), aSplineY.Evaluate(t),
                       aSplineZ.Evaluate(t))
 
# Create the polyline.
lines = vtk.vtkCellArray()
lines.InsertNextCell(numberOfOutputPoints)
for i in range(0, numberOfOutputPoints):
    lines.InsertCellPoint(i)
 
profileData.SetPoints(points)
profileData.SetLines(lines)
 
# Add thickness to the resulting line.
profileTubes = vtk.vtkTubeFilter()
profileTubes.SetNumberOfSides(8)
profileTubes.SetInput(profileData)
profileTubes.SetRadius(.05)
 
profileMapper = vtk.vtkPolyDataMapper()
profileMapper.SetInputConnection(profileTubes.GetOutputPort())
 
profile = vtk.vtkActor()
profile.SetMapper(profileMapper)
profile.GetProperty().SetDiffuseColor(banana)
profile.GetProperty().SetSpecular(.3)
profile.GetProperty().SetSpecularPower(30)
 
# Now create the RenderWindow, Renderer and Interactor
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
 
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
 
# Add the actors
ren.AddActor(glyph)
ren.AddActor(profile)
 
renWin.SetSize(500, 500)
 
iren.Initialize()
renWin.Render()
iren.Start()

« Older Entries