Index: lib/matplotlib/image.py
===================================================================
--- lib/matplotlib/image.py	(revision 6039)
+++ lib/matplotlib/image.py	(working copy)
@@ -394,6 +394,7 @@
                 ):
         AxesImage.__init__(self, ax,
                            **kwargs)
+        AxesImage.set_interpolation(self, 'nearest')
 
     def make_image(self, magnification=1.0):
         if self._A is None:
@@ -407,7 +408,9 @@
         height *= magnification
         im = _image.pcolor(self._Ax, self._Ay, self._A,
                            height, width,
-                           (x0, x0+v_width, y0, y0+v_height))
+                           (x0, x0+v_width, y0, y0+v_height),
+                           self._interpd[self._interpolation])
+
         fc = self.axes.patch.get_facecolor()
         bg = mcolors.colorConverter.to_rgba(fc, 0)
         im.set_bg(*bg)
@@ -450,8 +453,8 @@
         raise NotImplementedError('Method not supported')
 
     def set_interpolation(self, s):
-        if s != None and s != 'nearest':
-            raise NotImplementedError('Only nearest neighbor supported')
+        if s != None and not s in ('nearest','bilinear'):
+            raise NotImplementedError('Only nearest neighbor and bilinear interpolations are supported')
         AxesImage.set_interpolation(self, s)
 
     def get_extent(self):
Index: src/_image.cpp
===================================================================
--- src/_image.cpp	(revision 6039)
+++ src/_image.cpp	(working copy)
@@ -1138,21 +1138,211 @@
   return Py::asObject(imo);
 }
 
+// utilities for irregular grids
+void _bin_indices_middle(unsigned int *irows, int nrows, float *ys1, int ny,float dy, float y_min)
+{ int  i,j, j_last;
+   unsigned  int * rowstart = irows;
+  float *ys2 = ys1+1;
+  float *yl = ys1 + ny ;
+  float yo = y_min + dy/2.0;
+  float ym = 0.5f*(*ys1 + *ys2);
+  // y/rows
+   j = 0;
+   j_last = j;
+  for (i=0;i<nrows;i++,yo+=dy,rowstart++) {
+      while(ys2 != yl && yo > ym) {
+          ys1 = ys2;
+          ys2 = ys1+1;
+          ym = 0.5f*(*ys1 + *ys2);
+          j++;
+      }
+      *rowstart = j - j_last;
+      j_last = j;
+  }
+}
 
+void _bin_indices_middle_linear(float *arows, unsigned int *irows, int nrows, float *y, int ny,float dy, float y_min)
+{ int i;
+    int ii = 0;
+    int iilast = ny-1;
+    float sc = 1/dy;
+    int iy0 = (int)floor(sc * (y[ii]  - y_min) );
+    int iy1 = (int)floor(sc * (y[ii+1]  - y_min) );
+    float invgap=1.0f/(iy1-iy0);
+    for (i=0; i<nrows && i<=iy0; i++) {
+        irows[i] = 0;
+        arows[i] = 1.0;
+        //std::cerr<<"i="<<i<<"  ii="<<0<<" a="<< arows[i]<< std::endl;
+    }
+    for (; i<nrows; i++) {
+        while (i > iy1 && ii < iilast) {
+            ii++;
+            iy0 = iy1;
+            iy1 = (int)floor(sc * (y[ii+1] - y_min));
+            invgap=1.0f/(iy1-iy0);
+        }
+        if (i >= iy0 && i <= iy1) {
+            irows[i] = ii;
+            arows[i]=(iy1-i)*invgap;
+            //std::cerr<<"i="<<i<<"  ii="<<ii<<" a="<< arows[i]<< std::endl;
+        }
+        else break;
+    }
+    for (; i<nrows; i++) {
+        irows[i] =iilast-1;
+        arows[i] = 0.0;
+        //std::cerr<<"i="<<i<<"  ii="<<iilast-1<<" a="<< arows[i]<< std::endl;
+    }
+}
+
+void _bin_indices(int *irows, int nrows, double *y, int ny,
+                                            double sc, double offs)
+{
+    int i;
+    if (sc*(y[ny-1] - y[0]) > 0)
+    {
+        int ii = 0;
+        int iilast = ny-1;
+        int iy0 = (int)floor(sc * (y[ii]  - offs));
+        int iy1 = (int)floor(sc * (y[ii+1]  - offs));
+        for (i=0; i<nrows && i<iy0; i++) {
+            irows[i] = -1;
+        }
+        for (; i<nrows; i++) {
+            while (i > iy1 && ii < iilast) {
+                ii++;
+                iy0 = iy1;
+                iy1 = (int)floor(sc * (y[ii+1] - offs));
+            }
+            if (i >= iy0 && i <= iy1) irows[i] = ii;
+            else break;
+        }
+        for (; i<nrows; i++) {
+            irows[i] = -1;
+        }
+    }
+    else
+    {
+        int iilast = ny-1;
+        int ii = iilast;
+        int iy0 = (int)floor(sc * (y[ii]  - offs));
+        int iy1 = (int)floor(sc * (y[ii-1]  - offs));
+        for (i=0; i<nrows && i<iy0; i++) {
+            irows[i] = -1;
+        }
+        for (; i<nrows; i++) {
+            while (i > iy1 && ii > 1) {
+                ii--;
+                iy0 = iy1;
+                iy1 = (int)floor(sc * (y[ii-1] - offs));
+            }
+            if (i >= iy0 && i <= iy1) irows[i] = ii-1;
+            else break;
+        }
+        for (; i<nrows; i++) {
+            irows[i] = -1;
+        }
+    }
+}
+
+void _bin_indices_linear(float *arows, int *irows, int nrows, double *y, int ny,
+                                            double sc, double offs)
+{
+    int i;
+    if (sc*(y[ny-1] - y[0]) > 0)
+    {
+        int ii = 0;
+        int iilast = ny-1;
+        int iy0 = (int)floor(sc * (y[ii]  - offs));
+        int iy1 = (int)floor(sc * (y[ii+1]  - offs));
+        float invgap=1.0/(iy1-iy0);
+        for (i=0; i<nrows && i<iy0; i++) {
+            irows[i] = -1;
+        }
+        for (; i<nrows; i++) {
+            while (i > iy1 && ii < iilast) {
+                ii++;
+                iy0 = iy1;
+                iy1 = (int)floor(sc * (y[ii+1] - offs));
+                invgap=1.0/(iy1-iy0);
+            }
+            if (i >= iy0 && i <= iy1) {
+                irows[i] = ii;
+                arows[i]=(iy1-i)*invgap;
+            }
+            else break;
+        }
+        for (; i<nrows; i++) {
+            irows[i] = -1;
+        }
+    }
+    else
+    {
+        int iilast = ny-1;
+        int ii = iilast;
+        int iy0 = (int)floor(sc * (y[ii]  - offs));
+        int iy1 = (int)floor(sc * (y[ii-1]  - offs));
+        float invgap=1.0/(iy1-iy0);
+        for (i=0; i<nrows && i<iy0; i++) {
+            irows[i] = -1;
+        }
+        for (; i<nrows; i++) {
+            while (i > iy1 && ii > 1) {
+                ii--;
+                iy0 = iy1;
+                iy1 = (int)floor(sc * (y[ii-1] - offs));
+                invgap=1.0/(iy1-iy0);
+            }
+            if (i >= iy0 && i <= iy1) {
+                irows[i] = ii-1;
+                arows[i]=(i-iy0)*invgap;
+            }
+            else break;
+        }
+        for (; i<nrows; i++) {
+            irows[i] = -1;
+        }
+    }
+}
+
+
+
 char __image_module_pcolor__doc__[] =
 "pcolor(x, y, data, rows, cols, bounds)\n"
 "\n"
-"Generate a psudo-color image from data on a non-univorm grid using\n"
-"nearest neighbour interpolation.\n"
+"Generate a pseudo-color image from data on a non-uniform grid using\n"
+"nearest neighbour or linear interpolation.\n"
 "bounds = (x_min, x_max, y_min, y_max)\n"
+"interpolation = NEAREST or BILINEAR \n"
 ;
+
+void _pcolor_cleanup(PyArrayObject* x, PyArrayObject* y,  PyArrayObject *d,
+                                        unsigned int * rowstarts ,unsigned int*colstarts ,
+                                        float *acols , float *arows) {
+    if (x)
+      Py_XDECREF(x);
+    if (y)
+      Py_XDECREF(y);
+    if(d)
+      Py_XDECREF(d);
+    if(rowstarts)
+      PyMem_Free(rowstarts);
+    if(colstarts)
+      PyMem_Free(colstarts);
+    if(acols)
+      PyMem_Free(acols);
+    if(arows)
+      PyMem_Free(arows);
+    return;
+}
+    
 Py::Object
 _image_module::pcolor(const Py::Tuple& args) {
   _VERBOSE("_image_module::pcolor");
 
 
-  if (args.length() != 6)
-      throw Py::TypeError("Incorrect number of arguments (6 expected)");
+  if (args.length() != 7)
+      throw Py::TypeError("Incorrect number of arguments (7 expected)");
 
   Py::Object xp = args[0];
   Py::Object yp = args[1];
@@ -1160,6 +1350,7 @@
   unsigned int rows = Py::Int(args[3]);
   unsigned int cols = Py::Int(args[4]);
   Py::Tuple bounds = args[5];
+  unsigned int interpolation = Py::Int(args[6]);
 
   if (rows > 1 << 15 || cols > 1 << 15) {
     throw Py::ValueError("rows and cols must both be less than 32768");
@@ -1179,26 +1370,29 @@
   // Check we have something to output to
   if (rows == 0 || cols ==0)
       throw Py::ValueError("Cannot scale to zero size");
-
+  
+  PyArrayObject *x = NULL; PyArrayObject *y = NULL; PyArrayObject *d = NULL;
+  unsigned int * rowstarts = NULL; unsigned int*colstarts = NULL;
+  float *acols = NULL; float *arows = NULL;
+  
   // Get numpy arrays
-  PyArrayObject *x = (PyArrayObject *) PyArray_ContiguousFromObject(xp.ptr(), PyArray_FLOAT, 1, 1);
-  if (x == NULL)
+  x = (PyArrayObject *) PyArray_ContiguousFromObject(xp.ptr(), PyArray_FLOAT, 1, 1);
+  if (x == NULL) {
+      _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::ValueError("x is of incorrect type (wanted 1D float)");
-  PyArrayObject *y = (PyArrayObject *) PyArray_ContiguousFromObject(yp.ptr(), PyArray_FLOAT, 1, 1);
+  }
+  y = (PyArrayObject *) PyArray_ContiguousFromObject(yp.ptr(), PyArray_FLOAT, 1, 1);
   if (y == NULL) {
-      Py_XDECREF(x);
-      throw Py::ValueError("y is of incorrect type (wanted 1D float)");
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
+     throw Py::ValueError("y is of incorrect type (wanted 1D float)");
   }
-  PyArrayObject *d = (PyArrayObject *) PyArray_ContiguousFromObject(dp.ptr(), PyArray_UBYTE, 3, 3);
+  d = (PyArrayObject *) PyArray_ContiguousFromObject(dp.ptr(), PyArray_UBYTE, 3, 3);
   if (d == NULL) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::ValueError("data is of incorrect type (wanted 3D UInt8)");
   }
   if (d->dimensions[2] != 4) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
-      Py_XDECREF(d);
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::ValueError("data must be in RGBA format");
   }
 
@@ -1206,26 +1400,21 @@
   int nx = x->dimensions[0];
   int ny = y->dimensions[0];
   if (nx != d->dimensions[1] || ny != d->dimensions[0]) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
-      Py_XDECREF(d);
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::ValueError("data and axis dimensions do not match");
   }
 
   // Allocate memory for pointer arrays
-  unsigned int * rowstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int)*rows));
-  if (rowstarts == NULL) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
-      Py_XDECREF(d);
+  rowstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int)*rows));
+  arows = reinterpret_cast<float *>(PyMem_Malloc(sizeof(float)*rows));
+  if (rowstarts == NULL || arows == NULL ) {
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::MemoryError("Cannot allocate memory for lookup table");
   }
-  unsigned int * colstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int*)*cols));
-  if (colstarts == NULL) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
-      Py_XDECREF(d);
-      PyMem_Free(rowstarts);
+  colstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int)*cols));
+  acols = reinterpret_cast<float*>(PyMem_Malloc(sizeof(float)*cols));
+  if (colstarts == NULL || acols == NULL) {
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::MemoryError("Cannot allocate memory for lookup table");
   }
 
@@ -1238,54 +1427,17 @@
   size_t NUMBYTES(rows * cols * 4);
   agg::int8u *buffer = new agg::int8u[NUMBYTES];
   if (buffer == NULL) {
-      Py_XDECREF(x);
-      Py_XDECREF(y);
-      Py_XDECREF(d);
-      PyMem_Free(rowstarts);
-      PyMem_Free(colstarts);
+     _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
       throw Py::MemoryError("Could not allocate memory for image");
   }
+  
 
   // Calculate the pointer arrays to map input x to output x
-  unsigned int i, j, j_last;
+  unsigned int i, j;
   unsigned int * colstart = colstarts;
   unsigned int * rowstart = rowstarts;
   float *xs1 = reinterpret_cast<float*>(x->data);
   float *ys1 = reinterpret_cast<float*>(y->data);
-  float *xs2 = xs1+1;
-  float *ys2 = ys1+1;
-  float *xl = xs1 + nx - 1;
-  float *yl = ys1 + ny - 1;
-  float xo = x_min + dx/2.0;
-  float yo = y_min + dy/2.0;
-  float xm = 0.5*(*xs1 + *xs2);
-  float ym = 0.5*(*ys1 + *ys2);
-  // x/cols
-  j = 0;
-  j_last = j;
-  for (i=0;i<cols;i++,xo+=dx,colstart++) {
-      while(xs2 != xl && xo > xm) {
-          xs1 = xs2;
-          xs2 = xs1+1;
-          xm = 0.5f*(*xs1 + *xs2);
-          j++;
-      }
-      *colstart = j - j_last;
-      j_last = j;
-  }
-  // y/rows
-  j = 0;
-  j_last = j;
-  for (i=0;i<rows;i++,yo+=dy,rowstart++) {
-      while(ys2 != yl && yo > ym) {
-          ys1 = ys2;
-          ys2 = ys1+1;
-          ym = 0.5f*(*ys1 + *ys2);
-          j++;
-      }
-      *rowstart = j - j_last;
-      j_last = j;
-  }
 
 
   // Copy data to output buffer
@@ -1297,87 +1449,82 @@
   agg::int8u * position = buffer;
   agg::int8u * oldposition = NULL;
   start = reinterpret_cast<unsigned char*>(d->data);
-  for(i=0;i<rows;i++,rowstart++)
-  {
-      if (i > 0 && *rowstart == 0) {
-          memcpy(position, oldposition, rowsize*sizeof(agg::int8u));
-          oldposition = position;
-          position += rowsize;
-      } else {
-          oldposition = position;
-          start += *rowstart * inrowsize;
-          inposition = start;
-          for(j=0,colstart=colstarts;j<cols;j++,position+=4,colstart++) {
-              inposition += *colstart * 4;
-              memcpy(position, inposition, 4*sizeof(agg::int8u));
+    int s0 = d->strides[0];
+    int s1 = d->strides[1];
+  
+  if(interpolation == Image::NEAREST) {
+      _bin_indices_middle(colstart, cols, xs1,  nx,dx,x_min);
+      _bin_indices_middle(rowstart, rows, ys1,  ny, dy,y_min);
+      for(i=0;i<rows;i++,rowstart++)
+      {
+          if (i > 0 && *rowstart == 0) {
+              memcpy(position, oldposition, rowsize*sizeof(agg::int8u));
+              oldposition = position;
+              position += rowsize;
+          } else {
+              oldposition = position;
+              start += *rowstart * inrowsize;
+              inposition = start;
+              for(j=0,colstart=colstarts;j<cols;j++,position+=4,colstart++) {
+                  inposition += *colstart * 4;
+                  memcpy(position, inposition, 4*sizeof(agg::int8u));
+              }
           }
       }
   }
+  else if(interpolation == Image::BILINEAR) {
+      _bin_indices_middle_linear(acols, colstart, cols, xs1,  nx,dx,x_min);
+      _bin_indices_middle_linear(arows, rowstart, rows, ys1,  ny, dy,y_min);
+      double a00,a01,a10,a11,alpha,beta;
+      
+      
+        agg::int8u * start00;
+        agg::int8u * start01;
+        agg::int8u * start10;
+        agg::int8u * start11;
+      // Copy data to output buffer
+      for (i=0; i<rows; i++)
+        {
+            for (j=0; j<cols; j++)
+            {
+                alpha=arows[i];
+                beta=acols[j];
+                
+                a00=alpha*beta;
+                a01=alpha*(1.0-beta);
+                a10=(1.0-alpha)*beta;
+                a11=1.0-a00-a01-a10;
+                
+                start00=(agg::int8u *)(start + s0*rowstart[i] + s1*colstart[j]);
+                start01=start00+s1;
+                start10=start00+s0;
+                start11=start10+s1;
+                position[0] =(agg::int8u)(start00[0]*a00+start01[0]*a01+start10[0]*a10+start11[0]*a11);
+                position[1] =(agg::int8u)(start00[1]*a00+start01[1]*a01+start10[1]*a10+start11[1]*a11);
+                position[2] =(agg::int8u)(start00[2]*a00+start01[2]*a01+start10[2]*a10+start11[2]*a11);
+                position[3] =(agg::int8u)(start00[3]*a00+start01[3]*a01+start10[3]*a10+start11[3]*a11);
+                position += 4;
+            }
+        }
+        
+    }
 
+
   // Attatch output buffer to output buffer
   imo->rbufOut = new agg::rendering_buffer;
   imo->bufferOut = buffer;
   imo->rbufOut->attach(imo->bufferOut, imo->colsOut, imo->rowsOut, imo->colsOut * imo->BPP);
 
-  Py_XDECREF(x);
-  Py_XDECREF(y);
-  Py_XDECREF(d);
-  PyMem_Free(rowstarts);
-  PyMem_Free(colstarts);
+  
+   _pcolor_cleanup(x,y,d,rowstarts,colstarts,acols,arows);
 
   return Py::asObject(imo);
-}
+    
 
-void _bin_indices(int *irows, int nrows, double *y, int ny,
-                                            double sc, double offs)
-{
-    int i;
-    if (sc*(y[ny-1] - y[0]) > 0)
-    {
-        int ii = 0;
-        int iilast = ny-1;
-        int iy0 = (int)floor(sc * (y[ii]  - offs));
-        int iy1 = (int)floor(sc * (y[ii+1]  - offs));
-        for (i=0; i<nrows && i<iy0; i++) {
-            irows[i] = -1;
-        }
-        for (; i<nrows; i++) {
-            while (i > iy1 && ii < iilast) {
-                ii++;
-                iy0 = iy1;
-                iy1 = (int)floor(sc * (y[ii+1] - offs));
-            }
-            if (i >= iy0 && i <= iy1) irows[i] = ii;
-            else break;
-        }
-        for (; i<nrows; i++) {
-            irows[i] = -1;
-        }
-    }
-    else
-    {
-        int iilast = ny-1;
-        int ii = iilast;
-        int iy0 = (int)floor(sc * (y[ii]  - offs));
-        int iy1 = (int)floor(sc * (y[ii-1]  - offs));
-        for (i=0; i<nrows && i<iy0; i++) {
-            irows[i] = -1;
-        }
-        for (; i<nrows; i++) {
-            while (i > iy1 && ii > 1) {
-                ii--;
-                iy0 = iy1;
-                iy1 = (int)floor(sc * (y[ii-1] - offs));
-            }
-            if (i >= iy0 && i <= iy1) irows[i] = ii-1;
-            else break;
-        }
-        for (; i<nrows; i++) {
-            irows[i] = -1;
-        }
-    }
+  
 }
 
+
 char __image_module_pcolor2__doc__[] =
 "pcolor2(x, y, data, rows, cols, bounds, bg)\n"
 "\n"
@@ -1477,7 +1624,7 @@
         Py_XDECREF(bg);
         throw Py::MemoryError("Cannot allocate memory for lookup table");
     }
-    int * jcols = reinterpret_cast<int*>(PyMem_Malloc(sizeof(int*)*cols));
+    int * jcols = reinterpret_cast<int*>(PyMem_Malloc(sizeof(int)*cols));
     if (jcols == NULL) {
         Py_XDECREF(x);
         Py_XDECREF(y);
