From: Peter Groszkowski <pgroszko@ge...>  20040414 18:45:35

Hello Everyone: Before I reinvent the wheel, I'd thought I'd ask: I have a set of x,y points in a plane, and to each a corresponding z. They are spread out in circular concentric orbits. I need to do a flat surface plot (via pcolor) of these z values. A simple algorithm would be: #get data for x, y, z.. this would be some fixed points x,y inside #a box 4<=x<=4 and 4<=y<=4 with a z value for each. #Create a grid xi=linspace(4,4,200) yi=linspace(4,4,200) [Xi,Yi]=meshgrid(xi,yi) # This is the key step; find a z value for each grid point. Zi=griddata(x,y,z,Xi,Yi) pcolor(Xi, Yi, Zi, shading='flat') show() The key here is this griddata function. It interpolates the surface at the grid points based on the few data points I provide it with (x,y,z). It uses Delaunay triangulation, and is part of the standard MATLAB distribution (more about it here: http://www.mathworks.com/access/helpdesk/help/techdoc/ref/griddata.html#1007265). My question is whether anyone has come across its implementation in python. I could not find any information in any of the "standard" math packages. How do people do this? Thanks, Peter 
From: Andrew Straw <strawman@as...>  20040414 20:44:32
Attachments:
PGP.sig

On Apr 14, 2004, at 11:37 AM, Peter Groszkowski wrote: > Hello Everyone: > > Before I reinvent the wheel, I'd thought I'd ask: > > I have a set of x,y points in a plane, and to each a corresponding z. > They are spread out in circular concentric orbits. > I need to do a flat surface plot (via pcolor) of these z values. A > simple algorithm would be: > > #get data for x, y, z.. this would be some fixed points x,y inside > #a box 4<=x<=4 and 4<=y<=4 with a z value for each. > > #Create a grid > xi=linspace(4,4,200) > yi=linspace(4,4,200) > [Xi,Yi]=meshgrid(xi,yi) > > # This is the key step; find a z value for each grid point. > Zi=griddata(x,y,z,Xi,Yi) > > pcolor(Xi, Yi, Zi, shading='flat') > > show() > > The key here is this griddata function. It interpolates the surface at > the grid points based on the few data points I provide it with > (x,y,z). It uses Delaunay triangulation, and is part of the standard > MATLAB distribution (more about it here: > http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ > griddata.html#1007265). My question is whether anyone has come across > its implementation in python. I could not find any information in any > of the "standard" math packages. How do people do this? VTK, including the Python API, does Delaunay triangulation, and I've used it to do what you're talking about. IIRC, there are other Delaunay triangulation codes with Python bindings around. However, the Matlab function is a little more complicated; they use interpolation in addition to a straight Delaunay triangulation, and it often produces nicelooking results. They cite a reference in either the help file or the .m file; one can presumably read the reference and reimplement it for use in Python. FYI, here's some relevant Python VTK code: profile = vtkPolyData() profile.SetPoints(points) delny = vtk.vtkDelaunay2D() delny.SetInput(profile) delny.SetTolerance(0.01) delny.BoundingTriangulationOff() normals = vtk.vtkPolyDataNormals() normals.SetInput(delny.GetOutput()) lo,hi = zrange elevation = vtk.vtkElevationFilter() elevation.SetInput(delny.GetOutput()) elevation.SetLowPoint(0, 0, lo) elevation.SetHighPoint(0, 0, hi) elevation.SetScalarRange(lo, hi) elevation.ReleaseDataFlagOn() delMesh = vtk.vtkPolyDataMapper() delMesh.SetInput(elevation.GetOutput()) delMesh.SetLookupTable(lutLand) # doesn't work delActor = vtk.vtkActor() delActor.SetMapper(delMesh) ren1.AddActor( delActor ) Cheers! Andrew 
From: Peter Groszkowski <pgroszko@ge...>  20040415 21:54:58

> > VTK, including the Python API, does Delaunay triangulation, and I've > used it to do what you're talking about. IIRC, there are other > Delaunay triangulation codes with Python bindings around. However, > the Matlab function is a little more complicated; they use > interpolation in addition to a straight Delaunay triangulation, and > it often produces nicelooking results. They cite a reference in > either the help file or the .m file; one can presumably read the > reference and reimplement it for use in Python. > > FYI, here's some relevant Python VTK code: > > profile = vtkPolyData() > profile.SetPoints(points) > > delny = vtk.vtkDelaunay2D() > delny.SetInput(profile) > delny.SetTolerance(0.01) > delny.BoundingTriangulationOff() > > normals = vtk.vtkPolyDataNormals() > normals.SetInput(delny.GetOutput()) > > lo,hi = zrange > elevation = vtk.vtkElevationFilter() > elevation.SetInput(delny.GetOutput()) > elevation.SetLowPoint(0, 0, lo) > elevation.SetHighPoint(0, 0, hi) > elevation.SetScalarRange(lo, hi) > elevation.ReleaseDataFlagOn() > > delMesh = vtk.vtkPolyDataMapper() > delMesh.SetInput(elevation.GetOutput()) > delMesh.SetLookupTable(lutLand) # doesn't work > delActor = vtk.vtkActor() > delActor.SetMapper(delMesh) > ren1.AddActor( delActor ) Thanks Andrew. I have tried octave with the octiveforge extensions. In fact I have tried two different versions with three different versions of octaveforge. All seem to be broken in one way or another. I get random segmentation faults here and there on data that other times gets processed properly. Their implementation of griddata gives me some areas (triangles) which are empty even though MATLAB produces nice plots  perhaps some convergence issue. In any case I tried using octave's griddata and plot the result using matplotlib's pcolor. If not for the empty triangles (due to griddata) the plots are precisely what I need. I will look into VTK and see what it can do. Peter 