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)

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

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)
         ren1.AddActor( delActor )

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)

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 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)
        ren1.AddActor( delActor )

Thanks Andrew.

I have tried octave with the octive-forge extensions. In fact I have tried two different versions with three different versions of octave-forge. 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