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
    quadstrips = numpy.empty((nr-1,nc*4),'f')
    quadstrips[:,0::4] = x[:-1,:]
    quadstrips[:,2::4] = x[1:,:]
    quadstrips[:,1::4] = y[:-1,:]
    quadstrips[:,3::4] = y[1:,:]
    return quadstrips

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

def draw_quads(x,y,z,clim=[0,1],cmap=None):
    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)
    glMultiDrawArrays(GL_QUAD_STRIP, #starts, counts, nr)
                      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)
    glLoadIdentity()
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()

    glShadeModel(GL_FLAT)
    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)

    draw_quads(x,y,z,clim=[z.min(),z.max()],cmap=cmap)
    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 :slight_smile: ) 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
Graduate Research Assistant
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 :slight_smile: ) 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
Graduate Research Assistant
School of Meteorology
University of Oklahoma