Hi all,
here is a patch which implement bilinear interpolation on irregular grid
( i.e. it allows NonUniformImage
to accept both 'nearest' and 'bilinear' interpoaltion, instead of only
'nearest'.)
It is not perfect, given the current architecture of the image module I
think there is not simple way
to specify other interpolations (except for 'nearest', all
interpolations are implemented at the AGG level, not in
matplotlib itself).
For the same reason, it is not possible to interpolate before colormap
lookup instead of after (and this
is usually where the biggest effect is, when dealing with coarse grid).
However, I think it is already quite useful and the best one ca do
without a -very- extensive rewrite
of the matrix map modules....
BTW, I think it also corrects a small bug in the 'nearest'
interpolation: the last intervals was ignored
in the CVS version, now it is taken into account.
BOTH nearest and bilinear interpolation do the same for values oustside
the data grid: they just use boundary values (constant extrapolation).
I avoided linear extrapolation for 'bilinear', as extrapolation is
usually dangerous or meaningless (we can go outside the colormap very
fast)...
I have included a small example showing how both interpolation works....
Any remarks, could this be added before the next release? ;-)
Greg.
On Mon, 2008-08-11 at 16:50 +0200, Grégory Lielens wrote:
> On Fri, 2008-08-08 at 16:05 +0200, Grégory Lielens wrote:
> > Hello everybody,
> >
> > I have sent this message to the user group, but thinking of it, it may be
> > more
> > relevant to the development mailing list...so here it is again.
> >
> >
> >
> > We are looking for the best way to plot a waterfall diagram in
> > Matplotlib. The 2 functions which could be used
> > to do that are (as far as I have found) imshow and pcolormesh. Here is a
> > small script that use both to compare the output:
> >
> > -----------------
> >
> > from pylab import *
> >
> >
> > delta = 0.2
> > x = arange(-3.0, 3.0, delta)
> > y = arange(-2.0, 2.0, delta)
> > X, Y = meshgrid(x, y)
> > Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
> > Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
> > # difference of Gaussians
> > Z = 10.0 * (Z2 - Z1)
> > figure(1)
> > im = imshow(Z,extent=(-3,3,-2,2))
> > CS = contour(X, -Y, Z, 6,
> > colors='k', # negative contours will be dashed by default
> > )
> > clabel(CS, fontsize=9, inline=1)
> > title('Using imshow')
> > figure(2)
> > im = pcolormesh(X,-Y,Z)
> > CS = contour(X, -Y, Z, 6,
> > colors='k', # negative contours will be dashed by default
> > )
> > clabel(CS, fontsize=9, inline=1)
> > title('Using pcolormesh')
> > show()
> >
> > ---------------------
> >
> >
> > The problem is that we need some of the flexibility of pcolormesh (which
> > is able to map the matrix of value on any deformed mesh), while
> > we would like to use the interpolations available in imshow (which
> > explain why the imshow version is much "smoother" than the pcolormesh
> > one).
> >
> > In fact, what would be needed is not the full flexibility of pcolormesh
> > (which can map the grid to any kind of shape), we "only" have to deal
> > with rectangular grids where x- and y- graduations are irregularly spaced.
> >
> > Is there a drawing function in Matplotlib which would be able to work
> > with such a rectangular non-uniform grid?
> > And if not (and a quick look at the example and the code make me think
> > that indeed the capability is currently not present),
> > what about an extension of imshow which would work as this:
> >
> > im = imshow(Z,x_gridpos=x, y_gridpos=y) #specify the
> > position of the grid's nodes, instead of giving the extend and assuming
> > uniform spacing.
> >
> > Longer term, would a pcolormesh accepting interpolation be possible? The
> > current behavior, averaging the color of the grids node to get a uniform
> > cell color,
> > is quite rough except for a large number of cells...And even then, it
> > soon shows when you zoom in...
> >
> > The best would be to allow the same interpolations as in imshow (or a
> > subset of it), and also allows to use interpolation before colormap
> > lookup (or after),
> > like in Matlab. Indeed, Matlab allows to finely tune interpolation by
> > specifying Gouraud (interpolation after color
> > lookup)/Phong(interpolation before color lookup, i.e. for each pixel).
> > Phong is usually much better but also more CPU intensive. Phong is
> > especially when using discrete colormap, producing banded colors
> > equivalent to countour lines, while Gouraud does not work in those
> > cases.
> >
> > Of course, the performance will be impacted by some of those
> > interpolation options, which would degrade performance in animations for
> > example.... but I think that having the different options available
> > would be very useful, it allows to have the highest map quality, or have
> > a "quick and dirty" map depending on situation (grid spacing, type of
> > map, animation or not, ...).
> >
> > Best regards,
> >
> > Greg.
>
> I have found a method which implement the proposed extension to imshow:
> NonUniformImage...
>
> However, this image instance support only nearest neighbor
> interpolation. Trying to set the interpolation (using the
> set_interpolation method)
> to something akin imshow throw a "NotImplementedError: Only nearest
> neighbor supported" exception....
>
> So basically I am still stuck, it seems that currently there is no way
> in matplotlib to plot interpolated
> colormap on irregular rectangular grid, and even less on arbitrarily
> mapped grid...
>
> Is there any plans to add support for more interpolation in
> NonUniformImage in the future? Or maybe there is another
> drawing function that I did not find yet, with this ability?
>
>
> Best regards,
>
> Greg.
>
>
>
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> Matplotlib-devel mailing list
> Matplotlib-devel@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
Index: src/_image.cpp
===================================================================
--- src/_image.cpp (revision 6030)
+++ src/_image.cpp (working copy)
@@ -1138,21 +1138,190 @@
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.5*(*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.5*(*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.0/(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.0/(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"
;
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 +1329,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");
@@ -1214,14 +1384,16 @@
// Allocate memory for pointer arrays
unsigned int * rowstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int)*rows));
- if (rowstarts == NULL) {
+ float * arows = reinterpret_cast<float *>(PyMem_Malloc(sizeof(float)*rows));
+ if (rowstarts == NULL || arows == NULL ) {
Py_XDECREF(x);
Py_XDECREF(y);
Py_XDECREF(d);
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) {
+ unsigned int * colstarts = reinterpret_cast<unsigned int*>(PyMem_Malloc(sizeof(unsigned int)*cols));
+ float * acols = reinterpret_cast<float*>(PyMem_Malloc(sizeof(unsigned int)*cols));
+ if (colstarts == NULL || acols == NULL) {
Py_XDECREF(x);
Py_XDECREF(y);
Py_XDECREF(d);
@@ -1252,40 +1424,6 @@
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,23 +1435,67 @@
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;
@@ -1328,55 +1510,6 @@
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"
@@ -1477,7 +1610,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);
Index: lib/matplotlib/image.py
===================================================================
--- lib/matplotlib/image.py (revision 6030)
+++ lib/matplotlib/image.py (working copy)
@@ -407,7 +407,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 +452,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):
from matplotlib.pyplot import figure, show,cm
from matplotlib.pylab import meshgrid,bivariate_normal
import numpy as npy
from matplotlib.image import NonUniformImage
#x = (npy.random.rand(10)-0.5)*6
x=npy.arange(-3,3,0.5)
x.sort()#irregular grid -3->3
#y=( npy.random.rand(10)-0.5)*4
y=npy.arange(-2.5,2.5,0.5)
y.sort()#irregular grid -2->2
X, Y = meshgrid(x, y)
Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
# difference of Gaussians
Z = 10.0 * (Z2 - Z1)
cmap=cm.jet
gridon=True
print 'Size %d points' % (len(x) * len(y))
fig = figure(1)
ax = fig.add_subplot(121)
im = NonUniformImage(ax, extent=(-3,3,-3,3),cmap=cmap)
im.set_interpolation('nearest')
print im.get_interpolation()
im.set_data(x, y, Z)
ax.images.append(im)
if gridon:
for xg in x:
ax.plot([xg,xg],[-10,10],color='k', linestyle=':')
for yg in y:
ax.plot([-10,10],[yg,yg],color='k', linestyle=':')
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)
ax = fig.add_subplot(122)
im = NonUniformImage(ax, extent=(-3,3,-3,3),cmap=cmap)
im.set_interpolation('bilinear')
print im.get_interpolation()
im.set_data(x, y, Z)
ax.images.append(im)
if gridon:
for xg in x:
ax.plot([xg,xg],[-10,10],color='k', linestyle=':')
for yg in y:
ax.plot([-10,10],[yg,yg],color='k', linestyle=':')
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)
show()
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
