From: Curtis Cooper <curtis@hi...>  20040930 23:54:46

Hi, My research is in computational fluid dynamics (specifically, the meteorologies of planetary atmospheres). Working contour and vector plots in matplotlib would make it possible for me to make 2D meteorological maps of atmospheric layers, etc. I noticed for the first time in the goals page that contour plots are being worked on, apparently by STSci. I have been considering implementing these two plot types as sets of line collections, but now that I know contour plots are being worked on, and vector plots are simpler to implement (in 2D), I will work on making vector plots. The mathematics is fairly straightforward. I just need to learn how to use the class library. About contour plots, however, I have a couple of questions. How is it being implemented? I was about to try to write a marching squares contouring routine, although it might have been painfully slow in Python. Does anyone have experience with this? Cheers, Curtis 
From: Perry Greenfield <perry@st...>  20041001 02:17:03

Curtis Cooper writes: > My research is in computational fluid dynamics (specifically, the > meteorologies of planetary atmospheres). Working contour and vector plots > in matplotlib would make it possible for me to make 2D meteorological maps > of atmospheric layers, etc. > > I noticed for the first time in the goals page that contour plots are > being worked on, apparently by STSci. I have been considering > implementing these two plot types as sets of line collections, but now > that I know contour plots are being worked on, and vector plots are > simpler to implement (in 2D), I will work on making vector plots. The > mathematics is fairly straightforward. I just need to learn how to use > the class library. > > About contour plots, however, I have a couple of questions. How is it > being implemented? I was about to try to write a marching squares > contouring routine, although it might have been painfully slow in Python. > Does anyone have experience with this? > We are trying to adapt the C contour program that is used by gist (and can be found in the contour routine used by xplt in scipy). It would be best to look at the source for the precise description of the algorithm it uses (note though that gist apparently uses two different pieces of contour code for its contour tasks. The one we are looking to adapt, mainly because it appears much easier to isolate from the gist environment is the gcntr.c version). I would be amazed if one could find a pure Python algorithm to do contouring that was fast enough. Our current plan is to use these C routines to generate the contour segments, and do the plotting from within Python (as well as any contour labeling). If you have expertise in this area you may be able to do it better and faster than we can. Currently it is being worked on part time so we aren't able to do it as fast as we would like. I'm hoping that we will have at least a basic version (e.g., no labeling) in a couple weeks. If you want me to send or point you to the source code we are using as the basis, let me know. Perry Greenfield 
From: Curtis Cooper <curtis@hi...>  20041001 03:29:35

> We are trying to adapt the C contour program that is used by gist > (and can be found in the contour routine used by xplt in scipy). > It would be best to look at the source for the precise description > of the algorithm it uses (note though that gist apparently uses > two different pieces of contour code for its contour tasks. The > one we are looking to adapt, mainly because it appears much easier > to isolate from the gist environment is the gcntr.c version). > I would be amazed if one could find a pure Python algorithm to do > contouring that was fast enough. Our current plan is to use these > C routines to generate the contour segments, and do the plotting > from within Python (as well as any contour labeling). > > If you have expertise in this area you may be able to do it better > and faster than we can. Currently it is being worked on part time > so we aren't able to do it as fast as we would like. I'm hoping that we > will have at least a basic version (e.g., no labeling) in a couple > weeks. > > If you want me to send or point you to the source code we are > using as the basis, let me know. Hi again, Yes, I had the thought that using C for the algorithm would be easier as well. There are actually some very wellwritten marching squares contouring algorithms in C already out there. I will try to find such an implementation and point you to it or send you the source code. The second half is just the drawing, which should be implemented in matplotlib using the line collections class. Since vector plotting is not that hard, I will try to get that working first. Then, someone can take my source code and adapt it easily to the contouring problem, once an effective and sufficiently highperformance algorithm implementation can be found. Cheers, Curtis 
From: Perry Greenfield <perry@st...>  20041001 03:32:48

> Hi again, > > Yes, I had the thought that using C for the algorithm would be easier as > well. There are actually some very wellwritten marching squares > contouring algorithms in C already out there. I will try to find such an > implementation and point you to it or send you the source code. > Thanks, that would be helpful. In my search I didn't come across many. Keep in mind the license needs to be compatible with that of matplotlib. > The second half is just the drawing, which should be implemented in > matplotlib using the line collections class. Since vector plotting is not Yeah, that's what we have in mind. > that hard, I will try to get that working first. Then, someone can take > my source code and adapt it easily to the contouring problem, once an > effective and sufficiently highperformance algorithm implementation can > be found. > > Cheers, > Curtis > OK, Perry 
From: John Hunter <jdhunter@ac...>  20041001 13:45:09

From: Helge Avlesen <avle@ii...>  20041001 15:44:49

John Hunter <jdhunter@...> writes:  I was concerned by the fact that the lines were not smooth  if you  plot a connected line they line jumps from side to side. But it  does get the contour right, and is implemented in pure numeric, and  so it occurs to me that it might be easier to fix this problem than  start from scratch. Perhaps Helge or one of you has some insight  into how to fix this.   I'm attaching a modified version of the tarfile Helge initially sent  me. I've included a script testkont_mpl.py that calls Helge's lib.  Change the '.' linestyle to '' to see the problem I discussed. Hi, not sure if I have matplotlib 100% correctly installed, but this is what I see using your example script: http://www.ii.uib.no/~avle/mpl/c0.png (and with the current algorithm, more or less what I would expect...) to get straight lines you must plot segments one by one since they are not ordered. if I use gist for this(see the script at the end) I get http://www.ii.uib.no/~avle/mpl/c1.png the first points of the segments are given by the vectors (x1,y1) the second (x2,y2). you can get pretty lines in matplotlib as well, but only by using the scattered line drawing methods of gtk. (something like self.area.window.draw_segments(self.gc, zip( x1,y1,x2,y2)?) if you want do do it "right" in matplotlib, you should implement a contour following algorithm (in C)  with this I mean an routine that returns the linesegments defining each countour in bundles. the current alg. is sort of marching cubes in 2D, a simplified version of CONREC http://astronomy.swin.edu.au/~pbourke/projection/conrec/ but only using 2 triangles per square. doing contour following alg. it is also much easier to implement automatic contour labelling. I suspect python loops are too slow for such algorithms  it may perhaps be possible to do them in Numeric, but it will still be much slower than my simple library. I think you may use the GPL'ed PLPLOT (C) for an example of contour following alg. Helge from matplotlib.matlab import * import hutil delta = 0.05 x = y = arange(3.0, 3.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) Z = Z2Z1 print Z.shape #fsm = ones(Z.shape, Z.typecode()) fsm = ones(Z.shape, 'l') zmax, zmin = hutil.maxmin(Z) depths=linspace(zmin, zmax, 10) x1,y1,x2,y2 = hutil.contour2(Z, fsm, depths ) #imshow(Z, origin='lower', interpolation='nearest') #plot(y2,x2,'') #show() import gist gist.pldefault(dpi=100,style='framed.gs') gist.palette('rainbow.gp') gist.pli(transpose(Z)) gist.pldj(x1,y1,x2,y2) # draw disjoint segments 
From: John Hunter <jdhunter@ac...>  20041001 16:04:34

>>>>> "Helge" == Helge Avlesen <avle@...> writes: Helge> http://www.ii.uib.no/~avle/mpl/c1.png Helge> the first points of the segments are given by the vectors Helge> (x1,y1) the second (x2,y2). you can get pretty lines in Helge> matplotlib as well, but only by using the scattered line Helge> drawing methods of gtk. (something like Helge> self.area.window.draw_segments(self.gc, zip( x1,y1,x2,y2)?) OK, I see. I didn't fully understand that x1,x2,y1,y2 were the verts of unordered line segments. Then one can easily use a LineCollection to draw these efficiently in matplotlib  script below and screenshot http://nitace.bsd.uchicago.edu:8080/files/share/kontour.png. Jeez, I feel bad for sitting on this since February! Helge> if you want do do it "right" in matplotlib, you should Helge> implement a contour following algorithm (in C)  with this Helge> I mean an routine that returns the linesegments defining Helge> each countour in bundles. the current alg. is sort of Helge> marching cubes in 2D, a simplified version of CONREC Helge> http://astronomy.swin.edu.au/~pbourke/projection/conrec/ Helge> but only using 2 triangles per square. Do you have any thoughts on how we might do labels with your code? Helge> doing contour following alg. it is also much easier to Helge> implement automatic contour labelling. I suspect python Helge> loops are too slow for such algorithms  it may perhaps be Helge> possible to do them in Numeric, but it will still be much Helge> slower than my simple library. I think you may use the Helge> GPL'ed PLPLOT (C) for an example of contour following alg. We have a problem in that we cannot use GPL'd code in matplotlib because the GPL does not allow redistribution of closed code, which the matplotlib (and python license) do. If we decide to go with your routines, at least for the time being until we can "do it right", would you be willing to contribute your code to matplotlib under the matplotlib license (PSF inspired, free for commercial and noncommercial reuse)? Thanks! JDH from matplotlib.matlab import * from matplotlib.collections import LineCollection import hutil delta = 0.05 x = y = arange(3.0, 3.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) Z = Z2Z1 print Z.shape fsm = ones(Z.shape, Z.typecode()) zmax, zmin = hutil.maxmin(Z) depths=linspace(zmin, zmax, 10) x1,y1,x2,y2 = hutil.contour2(Z, fsm, depths ) imshow(Z, origin='lower', interpolation='nearest') segments = [ ( (thisy1, thisx1), (thisy2, thisx2) ) for thisx1, thisy1, thisx2, thisy2 in zip( x1,y1,x2,y2)] coll = LineCollection(segments) gca().add_collection(coll) savefig('kontour') show() 
From: Helge Avlesen <avle@ii...>  20041001 16:52:23

John Hunter <jdhunter@...> writes:   Do you have any thoughts on how we might do labels with your code? 2 ways: 1) automatically: define a coarser(user defined coarseness) mesh on the function to be labelled, and add the labels in the vertices of this mesh. the angle of each label can easily be found from the closest contour line segment. this method avoids clusters of labels effectively, but will not be good in areas with high variability. 2) manually: let the user point and click on contours, I implemented this for use with gist that you could take a look at (clabel below). matlab also does this, and I think this is the best way for real publication quality. I have never seen automatic routines that do labelling well. TECPLOT is close but not quite there. I also attached a routine to do this for 2D stretched coordinates (see the contour plots with labels on http://www.ii.uib.no/~avle/python.html for examples) which is called vclabel   If we decide to go with your routines, at least for the time being  until we can "do it right", would you be willing to contribute your  code to matplotlib under the matplotlib license (PSF inspired, free  for commercial and noncommercial reuse)? sure, no problem. Helge def clabel(z,clevels,opa=1,col='black',meth='std',digits=1): """ clabel(z,clevels,opa=1,col='black') At the point where the mouse is clicked, print the contour level from clevels that are closest to the interpolated value in this point. Meant to be useful for labeling contour lines manually... optional arguments: opa=0 : Transparent text. opa=1 : Erase background under label col='white' : text color. meth= 'std' bilinear interpo on a std grid, might give error near coast 'grid' values not given cellcentered but on corners of cells. 'bigrid' values taken from a bilinearly interpol fine mesh from futil.contour 'cell' takes value from cell directly Helge Avlesen <avle@...> """ print """ Insert contour levels by left clicking, middle button display values, right button finishes. """ button=0 while button<>3: mus=gist.mouse() button=mus[9] if meth=='bigrid': i=int(2*mus[0])+1 ; j=int(2*mus[1])+1 x=2*mus[0]i+1 ; y=2*mus[1]j+1 elif meth=='std': i=int( mus[0] ); j=int( mus[1] ) x=mus[0]i ; y=mus[1]j elif meth=='grid': xm=mus[0]+0.5 ; ym=mus[1]+0.5 i=int( xm ); j=int( ym ) x=xmi ; y=ymj elif meth=='cell': print mus[0], mus[1] i=int(round(mus[0])) ; j=int(round(mus[1])) if meth=='cell': val=z[i,j] print val,i,j else: # bilinear interpolation to find value a00=z[i,j] a10=z[i+1,j]a00 a01=z[i,j+1]a00 a11=z[i+1,j+1](a00+a10+a01) val=a00 + a10*x + a01*y + a11*x*y print val,i,j,x,y if button==1: # compare this value to the selected levels diff= abs( clevelsval ) # use the closest label=fpformat.fix( clevels[ Numeric.argmin(diff) ], digits ) gist.plt(label, mus[0], mus[1], opaque=opa, tosys=1, \ height=8, justify="CH", color=col ) def vclabel(z,sx,sy,clevels,opa=1,col='black',digits=1): """ manual(z,clevels,opa=1,col='black') At the point where the mouse is clicked, print the contour level from clevels that are closest to the interpolated value in this point. Meant to be useful for labeling contour lines manually... sx[i,j],sy[i,j] is the x,z coordinate of point z[i,j]. if [:,1] denotes the top layer, [:,kb1] the bottom (common in oceanography) j_is_down will be true. (z is always positive in the upward direction, but the indice j may go downwards) optional arguments: opa=0 : Transparent text. opa=1 : Use background color for text. col='white' : text color. digits: number of decimals in label Helge Avlesen <avle@...> """ print """ Insert contour levels by left clicking, middle button display depth, right button finishes. """ kb=z.shape[1] im=z.shape[0] j_is_down=0 if sy[0,0]>sy[0,1]: j_is_down=1 button=0 while button<>3: mus=gist.mouse() button=mus[9] # bisection search for the indices i=hbisect( sx[:,0], mus[0] ) if i<0 or i>im1: print 'outside:',i continue x=(mus[0]sx[i,0])/(sx[i+1,0]sx[i,0]) if j_is_down: finn=hbisect( (1.x)*sy[i,::1] + x*sy[i+1,::1] , mus[1] ) j=kb2finn if finn<0 or finn>kb1: print 'outside',i,finn continue xa=Numeric.array((sx[i,j+1]+x*(sx[i+1,j+1]sx[i,j+1]),\ sy[i,j+1]+x*(sy[i+1,j+1]sy[i,j+1]))) if mus[1]xa[1]<0: print 'below' continue xb=Numeric.array((sx[i,j]+x*(sx[i+1,j]sx[i,j]),\ sy[i,j]+x*(sy[i+1,j]sy[i,j]))) y=(((mus[0]xa[0])**2 + (mus[1]xa[1])**2 )\ /((xb[0]xa[0])**2 + (xb[1]xa[1])**2 ))**0.5 # bilinear interpolation to find value a1=z[i,j+1] a2=z[i+1,j+1]a1 a3=z[i,j]a1 a4=z[i+1,j](a1+a2+a3) val=a1 + a2*x + a3*y + a4*x*y else: print 'increasing j upwards not yet implemented' continue print 'x,y=',x,y,' i,j=',i,j, 'val=',val if button==1: # compare this value to the selected levels diff= abs( clevelsval ) # use the closest label=fpformat.fix( clevels[ Numeric.argmin(diff) ], 1) if opa==1: label=' '+label+' ' gist.plt(label, mus[0], mus[1], opaque=opa, tosys=1, \ height=8, justify="CH", color=col ) 
From: John Hunter <jdhunter@ac...>  20041001 16:37:29

>>>>> "Perry" == Perry Greenfield <perry@...> writes: >> Hi again, >> >> Yes, I had the thought that using C for the algorithm would be >> easier as well. There are actually some very wellwritten >> marching squares contouring algorithms in C already out there. >> I will try to find such an implementation and point you to it >> or send you the source code. >> Perry> Thanks, that would be helpful. In my search I didn't come Perry> across many. Keep in mind the license needs to be Perry> compatible with that of matplotlib. Of course, in addition to the license, there is the patent issue. I believe marching squares is patented. I know marching cubes is. http://patft.uspto.gov/netacgi/nphParser?Sect1=PTO1&Sect2=HITOFF&d=PALL&p=1&u=/netahtml/srchnum.htm&r=1&f=G&l=50&s1=4,710,876.WKU.&OS=PN/4,710,876&RS=PN/4,710,876 I checked the header of vtkMarchingSquares.cxx which states Program: Visualization Toolkit Module: $RCSfile: vtkMarchingSquares.cxx,v $ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. THIS CLASS IS PATENTED UNDER UNITED STATES PATENT NUMBER 4,710,876 "System and Method for the Display of Surface Structures Contained Within the Interior Region of a Solid Body". Application of this software for commercial purposes requires a license grant from GE. Contact: but the patent number they reference which is linked above begins with A method and apparatus for displaying *three dimensional surface images* includes the utilization of a case table for rapid retrieval of surface approximation information. emphasis mine. So I don't know for sure what the patent status of the 2D algorithm is. JDH 