# Surface Viewing

By Hal Canary, 2013-04-20 20:00:26 (link)
#computer-science;#computers-code

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()
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):
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)
.......
```

(back)