 # pcolor in action...

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)

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

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 nice-looking results. They cite a reference in either the help file or the .m file; one can presumably read the reference and re-implement 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)

Cheers!
Andrew

PGP.sig (155 Bytes)

···

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)

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 nice-looking results. They cite a reference in either the help file or the .m file; one can presumably read the reference and re-implement 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)