Surface Viewing

In the class I am sitting in on this semester - Visual Solid Shape, we use a piece of software called shapemonger.

I wanted to get a better view at monge patch (a surface defined by (x,y,f(x,y)) so I wrote a little Python-VTK program to generate a vtkPolyData of the surface.

#!/usr/bin/vtkpython
import vtk

function = lambda x,y: x**2 - y**2
outputName = '/tmp/monge.vtp'
resolution = 50
xmin, xmax = -1.0, 1.0
ymin, ymax = -1.0, 1.0

######
polyData = vtk.vtkPolyData()
points = vtk.vtkPoints()
index = 0
cells = vtk.vtkCellArray()
def makeQuad(cells,a,b,c,d):
    cells.InsertNextCell(4)
    cells.InsertCellPoint(a)
    cells.InsertCellPoint(b)
    cells.InsertCellPoint(c)
    cells.InsertCellPoint(d)
for i in xrange(resolution):
    x = xmin + (xmax-xmin) * i / float(resolution - 1)
    for j in xrange(resolution):
        y = ymin + (ymax - ymin) * j / float(resolution - 1)
        points.InsertPoint(index, x,y,function(x,y))
        if (i > 0) and (j > 0):
            makeQuad(
                cells,index,index-resolution,
                index-resolution-1,index-1)
        index += 1
polyData.SetPoints(points)
polyData.SetPolys(cells)
writer = vtk.vtkXMLPolyDataWriter()
writer.SetInput(polyData)
writer.SetFileName(outputName)
writer.Write()
print 'wrote to', outputName

The surface can be viewed in Paraview or vtkviewer.

To make this surface look better, calculate the normals. The (unnormalized) normal direction is (-∂f/∂x,-∂f/∂y,1).

.......
def partial(func,v,i):
    epsilon = 5e-7
    def f(v,i,d):
        return v\[:i\] + \[v\[i\] + d\] + v\[(i+1):\]
    return ( (1.0 / 12.0) \* func(\*f(v,i,-2\*epsilon)) +
             (-2.0 / 3.0) \* func(\*f(v,i,-1\*epsilon)) +
             ( 2.0 / 3.0) \* func(\*f(v,i,+1\*epsilon)) +
             (-1.0 /12.0) \* func(\*f(v,i,+2\*epsilon))
              ) / epsilon
def makenormal(func,x,y):
    n = 045;partial(func,\[x,y\],0),-partial(func,\[x,y\],1),1.0)
    f = 1.0 / math.sqrt(sum(x\*\*2 for x in n))
    return tuple(f \* x for x in n)
.......
normals = vtk.vtkDoubleArray()
normals.SetNumberOfComponents(3)
normals.SetNumberOfTuples(resolution\*\*2)
.......
        .......
        normals.InsertTupleValue(index, makenormal(function,x,y))
        .......
.......
polyData.GetPointData().SetNormals(normals)
.......