# opengl backend

Hi,

There was a recent discussion about opengl and matplotlib in the
context of matplotlib rendering speeds.

At the scipy sprints we put together a proof of concept renderer
for quad meshes using the opengl frame buffer object, which we
then render as a matplotlib image. Once the data is massaged
to the correct form we were seeing rendering speeds of about
10 million quads per second. See below.

Using this technique we can get the advantage of opengl without
having to write a full backend. Please let me know if you have
time to contribute to writing a more complete backend.

- Paul

-- meshgl.py --
# This program is public domain
import numpy
from pyglet.gl import *
import pyglet
from ctypes import POINTER,c_long,c_ubyte,c_int

def buildmesh(x,y):
"""
x,y are corners of the mesh quadrants
"""
# Quadlist is a list of pairs of vertices, left edge, right edge, right
# edge, one for each row. The resulting map will look like:
# [[p00 p10 p01 p11 p02 p12 ...] [p10 p20 p11 p21 ...] ...]
# Note that the x,y shape will be one more than the number of data values.
nr,nc = x.shape

def colorindex(z,clim,cmap):
# Calc colour indices
ncolors = cmap.shape[0]
step = float(clim[1]-clim[0])/(ncolors-1)
idx = numpy.asarray((z-clim[0])/step, dtype='i')
idx = numpy.clip(idx, 0, ncolors-1)
return idx

def buildcolor(z,clim=[0,1],cmap=None):
"""
z are centers of the mesh quadrants.

Note there size of z are less than size of x,y by 1 in each dim.
"""
# Generate colours
idx = colorindex(z, clim, cmap)

nr,nc = z.shape
# Quadlist is a list of pairs of vertices, left edge, right edge, right
# edge. The resulting map will look like:
# [[p00 p10 p01 p11 p02 p12 ...] [p10 p20 p11 p21 ...] ...]
# The z values correspond to the bottom right corner of each quad. This
# means we do not need the first two quads of each row and every odd
# quad after that. The resulting color list looks like the following:
# [[0 0 0 z00 0 z01 0 z02 ...] [0 0 0 z10 0 z11 0 z12 ...] ...]
# Each i,j has a four byte colour value.
colors = numpy.zeros((nr,2*nc+2,4),'B')
colors[:,3::2,:] = cmap[idx,:]
return colors

strips = buildmesh(x,y)
colors = buildcolor(z,clim,cmap)
nr,nc = strips.shape
nc = nc/2 # pair of vertices per point
starts = numpy.arange(0,nr*nc,nc,'i')
counts = numpy.ones(nr,'i')*nc

glEnableClientState(GL_VERTEX_ARRAY)
glEnableClientState(GL_COLOR_ARRAY)
glVertexPointer(2,GL_FLOAT,0,strips.ctypes.data)
glColorPointer(4,GL_UNSIGNED_BYTE,0,colors.ctypes.data)
starts.ctypes.data_as(POINTER(c_int)),
counts.ctypes.data_as(POINTER(c_int)),
nr)
glDisableClientState(GL_VERTEX_ARRAY)
glDisableClientState(GL_COLOR_ARRAY)

def set_scene(w,h):
"""Set up the openGL transform matrices"""
glViewport(0, 0, w, h)
glMatrixMode(GL_PROJECTION)
glMatrixMode(GL_MODELVIEW)

glEnable(GL_BLEND)
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

# Clear scene
glClearColor(1, 1, 1, 0)
glClear(GL_COLOR_BUFFER_BIT)

def renderbuffer(w,h,draw):
"""Render buffer to off screen framebuffer object"""

# Make sure we support frame buffers
if not pyglet.gl.gl_info.have_extension('GL_EXT_framebuffer_object'):
raise RuntimeError("This OpenGL does not support framebuffer objects")

# Create a frame buffer context
framebuffers = (GLuint*1)(0)
glGenFramebuffersEXT(1, framebuffers)
glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, framebuffers[0])

# Create a render buffer
renderbuffers = (GLuint*1)(0)
glGenRenderbuffersEXT(1, renderbuffers)
glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, renderbuffers[0])
glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT, GL_RGBA8, w, h)
glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT,
GL_RENDERBUFFER_EXT, renderbuffers[0])

status = glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT)
if status != GL_FRAMEBUFFER_COMPLETE_EXT:
raise RuntimeError("Could not create GL framebuffer")

# Render scene
set_scene(w,h)
draw()
glFlush()

# Retrieve buffer
buffer = numpy.empty((h,w,4),'B')
glReadPixels(0, 0, w, h, GL_RGBA, GL_UNSIGNED_BYTE, buffer.ctypes.data)
# View the buffer in the window
#glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0)
# Release buffer
glDeleteRenderbuffersEXT(1, renderbuffers)

return buffer

# === DEMO CODE BELOW ===

def testmesh(shape=(3,2)):
"""sample scene"""
# Quad mesh vertices (nr x nc)
nr,nc = shape
x = numpy.arange(nr)[None,:] + numpy.arange(nc)[:,None]
y = numpy.arange(nr)[None,:] + numpy.zeros((nc,1))
# Quad mesh values (nr-1 x nc-1)
z = x[:-1,:-1]

# Grayscale RGBA colormap
R = numpy.linspace(0,255,64)
alpha = numpy.ones(64,'B')*255
cmap = numpy.array([R,R,R,alpha],'B').T

if True:
# Green box in center using default -1,1 coords
glColor3f(0., 1., 0.)
glRectf(-0.5, 0.5, 0.5, -0.5)

# Set data coordinates
glPushMatrix()
glOrtho(0,x.max()*2,0,y.max()*2,-1,1)

if True:
# Blue box unit size at 1,1 in data coords
glColor4f(0., 0., 1., 1.)
glRectf(1, 2, 2, 1)

glPopMatrix()

class Window(pyglet.window.Window):
"""Pyglet driver"""
def __init__(self, w, h):
super(Window, self).__init__(width=w, height=h, resizable=True)
glClearColor(1.,1.,1.,1.)
self._dirty = True

def on_resize(self, w, h):
print "resize"
set_scene(w,h);
self._dirty = True

def draw(self):
if self._dirty:
print "draw"
#window.clear()
testmesh()
self._dirty = False
self.flip()

def run(self):
counter = 0
while not self.has_exit and counter < 1000:
self.dispatch_events()
self.draw()
counter += 1

def plot():
"""Matplotlib driver"""
import pylab
print "start rendering"
b = renderbuffer(400,500,lambda: testmesh((1000,1000)))
print "done rendering"
pylab.imshow(b,origin='lower')
pylab.show()

if __name__ == "__main__":
#Window(400,500).run()
plot()

Paul Kienzle wrote:

Hi,

There was a recent discussion about opengl and matplotlib in the
context of matplotlib rendering speeds.

At the scipy sprints we put together a proof of concept renderer
for quad meshes using the opengl frame buffer object, which we
then render as a matplotlib image. Once the data is massaged
to the correct form we were seeing rendering speeds of about
10 million quads per second. See below.

Using this technique we can get the advantage of opengl without
having to write a full backend. Please let me know if you have
time to contribute to writing a more complete backend.

I know that I (as the other contributer ) plan on making this into a full backend provided:

1) The quadmesh code can be shown yield the gains when integrated into matplotlib more fully

2) I find the time (which is a matter of *when*, not if)

I'm certainly finding thus far that pyglet makes it a lot easier to do a full backend than some of the other python->opengl methods I'd explored in the past.

Ryan

···

--
Ryan May
School of Meteorology
University of Oklahoma

On RHEL4, I get the following message:

This OpenGL does not support framebuffer objects

I understand that I'm probably just suffering from a relatively old Mesa/OpenGL stack here. And obviously, just because some older systems won't support this is not a reason to not include it as an optional backend.

Can you possibly send a screen shot to the list for those of us who can't see it? Is the quadmesh interpolated? That would be a huge improvement over anything we can do currently.

As for developing a full OpenGL backend, note that the process of writing a backend has been greatly simplified as of 0.98.x (there's only about four rendering methods to implement now, with three more optional ones to re-implement for performance reasons). You can start by looking at backend_bases.py. I'm happy to help with any questions you have as you go along, but as I don't have a working pyglet/OpenGL, I can't be of much help for testing.

Cheers,
Mike

Ryan May wrote:

···

Paul Kienzle wrote:

Hi,

There was a recent discussion about opengl and matplotlib in the
context of matplotlib rendering speeds.

At the scipy sprints we put together a proof of concept renderer
for quad meshes using the opengl frame buffer object, which we
then render as a matplotlib image. Once the data is massaged
to the correct form we were seeing rendering speeds of about
10 million quads per second. See below.

Using this technique we can get the advantage of opengl without
having to write a full backend. Please let me know if you have
time to contribute to writing a more complete backend.

I know that I (as the other contributer ) plan on making this into a full backend provided:

1) The quadmesh code can be shown yield the gains when integrated into matplotlib more fully

2) I find the time (which is a matter of *when*, not if)

I'm certainly finding thus far that pyglet makes it a lot easier to do a full backend than some of the other python->opengl methods I'd explored in the past.

Ryan

--
Michael Droettboom
Science Software Branch
Operations and Engineering Division
Space Telescope Science Institute
Operated by AURA for NASA

Michael Droettboom wrote:

On RHEL4, I get the following message:

This OpenGL does not support framebuffer objects

I understand that I'm probably just suffering from a relatively old Mesa/OpenGL stack here. And obviously, just because some older systems won't support this is not a reason to not include it as an optional backend.

Can you possibly send a screen shot to the list for those of us who can't see it? Is the quadmesh interpolated? That would be a huge improvement over anything we can do currently.

As for developing a full OpenGL backend, note that the process of writing a backend has been greatly simplified as of 0.98.x (there's only about four rendering methods to implement now, with three more optional ones to re-implement for performance reasons). You can start by looking at backend_bases.py. I'm happy to help with any questions you have as you go along, but as I don't have a working pyglet/OpenGL, I can't be of much help for testing.

I'll let you know if I need any help. 0.98 definitely seems simpler and I think Pyglet will handle a lot of the context stuff to make things simpler as well.

Your error on RHEL4 indicates more that it doesn't support our particular approach, which uses offscreen rendering to render to an image which is then blitted in. It's an expedient approach to get a drop in quadmesh replacement. A full pyglet backend would not suffer from this limitation.

I don't have this hooked up fully within matplotlib yet, so right now it just generates a test image with a 2000x2000 mesh of quadrilaterals which is displayed using imshow(). (simple example attached.) I have not looked at interpolation or even turned on antialiasing (should just be one or two function calls). I can tell you that most of the time is spent generating the list of vertices in a form suitable for OpenGL QUAD_LISTS. The actual rendering and handing off to the graphics card is *really* fast.

I have a couple other things on my plate ATM, but hopefully next week some time I will get back to this, as I really, really, really want/need to pursue the opengl stuff further.

Ryan

···

--
Ryan May