def Frechet_distances_2(input_vtk_polydata_n, input_vtk_polydata_m):
    # compute Frechet distance from an array of fibers to another array

    number_of_lines_m = input_vtk_polydata_m.GetNumberOfLines()
    number_of_lines_n = input_vtk_polydata_n.GetNumberOfLines()
    all_fibers_m = range(0,number_of_lines_m)
    all_fibers_n = range(0,number_of_lines_n)
    distances = numpy.zeros([number_of_lines_m,number_of_lines_n])
    
    #input_vtk_polydata2 = input_vtk_polydata
    
    input_vtk_polydata_m.GetLines().InitTraversal()
    line1_ptids = vtk.vtkIdList()
    
    inpoints1 = input_vtk_polydata_m.GetPoints()
    inpoints2 = input_vtk_polydata_n.GetPoints()
    
    for lidx1 in all_fibers_m:
        input_vtk_polydata_m.GetLines().GetNextCell(line1_ptids)
        
        input_vtk_polydata_n.GetLines().InitTraversal()
        line2_ptids = vtk.vtkIdList()
        
        for lidx2 in all_fibers_n:
            input_vtk_polydata_n.GetLines().GetNextCell(line2_ptids)

            distances[lidx1,lidx2] = _frechet_distance_internal_use(inpoints1,inpoints2,line1_ptids,line2_ptids)

    return distances

  def createGlyph(self):
    self.polyData = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    lines = vtk.vtkCellArray()
    self.polyData.SetPoints( points )
    self.polyData.SetLines( lines )
    prevPoint = None
    firstPoint = None

    for x,y in ((0,0),)*4:
      p = points.InsertNextPoint( x, y, 0 )
      if prevPoint is not None:
        idList = vtk.vtkIdList()
        idList.InsertNextId( prevPoint )
        idList.InsertNextId( p )
        self.polyData.InsertNextCell( vtk.VTK_LINE, idList )
      prevPoint = p
      if firstPoint is None:
        firstPoint = p

    # make the last line in the polydata
    idList = vtk.vtkIdList()
    idList.InsertNextId( p )
    idList.InsertNextId( firstPoint )
    self.polyData.InsertNextCell( vtk.VTK_LINE, idList )

  def copyFirstNLines(self, sourcePolyData, lineCount):
    """make a polydata with only the first N polylines"""

    polyData = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    polyData.SetPoints(points)

    lines = vtk.vtkCellArray()
    polyData.SetLines(lines)

    sourcePoints = sourcePolyData.GetPoints()
    sourceLines = sourcePolyData.GetLines()
    sourceIdList = vtk.vtkIdList()
    sourceLines.InitTraversal()
    while sourceLines.GetNextCell(sourceIdList):
        pointCount = sourceIdList.GetNumberOfIds()
        idList = vtk.vtkIdList()
        for idIndex in range(pointCount):
            sourceId = sourceIdList.GetId(idIndex)
            point = sourcePoints.GetPoint(sourceId)
            id = points.InsertNextPoint(point)
            idList.InsertNextId(id)
        lines.InsertNextCell(idList)
        if lines.GetNumberOfCells() > lineCount:
            break

    return polyData

def getNumberOfPointsSharedByTwoCells( pd, iCell1, iCell2 ):
    '''Compute the number of shared points between iCell1 and iCell2 in the vtkPolyData pd.'''
    cell1_points = vtk.vtkIdList()
    pd.GetCellPoints( iCell1, cell1_points )
    cell2_points = vtk.vtkIdList()
    pd.GetCellPoints( iCell2, cell2_points )
    cell1_points.IntersectWith( cell2_points )
    return cell1_points.GetNumberOfIds()

 def __init__(self):
     self._Surface = None
     self._SeedIds = None
     self._SourceSeedIds = vtk.vtkIdList()
     self._TargetSeedIds = vtk.vtkIdList()
     self.PrintError = None
     self.PrintLog = None
     self.InputText = None
     self.OutputText = None

  def createGlyph(self, polyData):
    """
    create a brush circle of the right radius in XY space
    - assume uniform scaling between XY and RAS which
      is enforced by the view interactors
    """
    sliceNode = self.sliceWidget.sliceLogic().GetSliceNode()
    self.rasToXY.DeepCopy(sliceNode.GetXYToRAS())
    self.rasToXY.Invert()
    maximum, maxIndex = 0,0
    for index in range(3):
      if abs(self.rasToXY.GetElement(0, index)) > maximum:
        maximum = abs(self.rasToXY.GetElement(0, index))
        maxIndex = index
    point = [0, 0, 0, 0]
    point[maxIndex] = self.radius
    xyRadius = self.rasToXY.MultiplyPoint(point)
    import math
    xyRadius = math.sqrt( xyRadius[0]**2 + xyRadius[1]**2 + xyRadius[2]**2 )

    if self.pixelMode:
      xyRadius = 0.01

    # make a circle paint brush
    points = vtk.vtkPoints()
    lines = vtk.vtkCellArray()
    polyData.SetPoints(points)
    polyData.SetLines(lines)
    PI = 3.1415926
    TWOPI = PI * 2
    PIoverSIXTEEN = PI / 16
    prevPoint = -1
    firstPoint = -1
    angle = 0
    while angle <= TWOPI:
      x = xyRadius * math.cos(angle)
      y = xyRadius * math.sin(angle)
      p = points.InsertNextPoint( x, y, 0 )
      if prevPoint != -1:
        idList = vtk.vtkIdList()
        idList.InsertNextId(prevPoint)
        idList.InsertNextId(p)
        polyData.InsertNextCell( vtk.VTK_LINE, idList )
      prevPoint = p
      if firstPoint == -1:
        firstPoint = p
      angle = angle + PIoverSIXTEEN

    # make the last line in the circle
    idList = vtk.vtkIdList()
    idList.InsertNextId(p)
    idList.InsertNextId(firstPoint)
    polyData.InsertNextCell( vtk.VTK_LINE, idList )

def create_cylinder(r=.5, h=1., res=30):
    pts = vtk.vtkPoints()
    pts.SetNumberOfPoints(10*res)

    pd = vtk.vtkPolyData()
    pd.Allocate(3*res, 1)

    ang = 2*math.pi/res

    ind = 0

    for sign in [-1, 1]:
        y = sign*h/2
        for i in range(res):
            c0, s0 = math.cos(i*ang), math.sin(i*ang)
            c1, s1 = math.cos((i+1)*ang), math.sin((i+1)*ang)

            x0, z0 = round(r*c0, 6), round(r*s0, 6)
            x1, z1 = round(r*c1, 6), round(r*s1, 6)

            pts.InsertPoint(ind, 0, y, 0)
            pts.InsertPoint(ind+1, x0, y, z0)
            pts.InsertPoint(ind+2, x1, y, z1)

            tri = vtk.vtkIdList()
            [ tri.InsertNextId(ind+j) for j in range(3) ]

            t = pd.InsertNextCell(vtk.VTK_TRIANGLE, tri)

            ind += 3

            if sign == 1:
                pd.ReverseCell(t)

                pts.InsertPoint(ind, x0, -h/2, z0)
                pts.InsertPoint(ind+1, x0, y, z0)

                pts.InsertPoint(ind+2, x1, y, z1)
                pts.InsertPoint(ind+3, x1, -h/2, z1)

                poly = vtk.vtkIdList()
                [ poly.InsertNextId(ind+j) for j in range(4) ]

                pd.InsertNextCell(vtk.VTK_POLYGON, poly)

                ind += 4

    pd.SetPoints(pts)

    prod = vtk.vtkTrivialProducer()
    prod.SetOutput(pd)

    return prod

    def VtkLoadElemMesh(self,field,defFScale=0.0,eigenMode=None):
        '''Load the element mesh

        :param field: scalar field to be represented
        :param defFScale: factor to apply to current displacement of nodes 
                  so that the display position of each node equals to
                  the initial position plus its displacement multiplied
                  by this factor. In case of modal analysis, the displayed 
                  position of each node equals to the initial position plus
                  its eigenVector multiplied by this factor.
                  (Defaults to 0.0, i.e. display of initial/undeformed shape)
        :param eigenMode: eigenvibration mode if we want to display the deformed 
                 shape associated with it when a modal analysis has been carried out. 
                 Defaults to None: no modal analysis.

        '''
        # Define grid
        self.nodes= vtk.vtkPoints()
        self.gridRecord.uGrid= vtk.vtkUnstructuredGrid()
        self.gridRecord.uGrid.SetPoints(self.nodes)
        eSet= self.gridRecord.xcSet
        eSet.numerate()
        self.gridRecord.uGrid.name= eSet.name+'_grid'
        # Scalar values.
        nodeSet= eSet.getNodes
        if(field):
            arr= field.fillArray(nodeSet)
            field.creaLookUpTable()      
        # Load nodes in vtk
        setNodes= eSet.getNodes
        if eigenMode==None:
            for n in setNodes:
                pos= n.getCurrentPos3d(defFScale)
                self.nodes.InsertPoint(n.getIdx,pos.x,pos.y,pos.z)
        else:
            for n in setNodes:
                pos= n.getEigenPos3d(defFScale,eigenMode)
                self.nodes.InsertPoint(n.getIdx,pos.x,pos.y,pos.z)
         # Load elements in vtk
        setElems= eSet.getElements
        for e in setElems:
            vertices= xc_base.vector_int_to_py_list(e.getIdxNodes)
            vtx= vtk.vtkIdList()
            for vIndex in vertices:
                vtx.InsertNextId(vIndex)
            if(e.getVtkCellType!= vtk.VTK_VERTEX):
                self.gridRecord.uGrid.InsertNextCell(e.getVtkCellType,vtx)
        setConstraints= eSet.getConstraints
        for c in setConstraints:
            vtx= vtk.vtkIdList()
            vtx.InsertNextId(c.getNodeIdx)
            if(c.getVtkCellType!= vtk.VTK_LINE):
                self.gridRecord.uGrid.InsertNextCell(c.getVtkCellType,vtx)

def getDual( pd ):
    '''Get the dual of a vtkPolyData. The finite parts only.'''
    pd.BuildLinks()
    cells = vtk.vtkCellArray()
    points = vtk.vtkPoints()
    for iPt in range(pd.GetNumberOfPoints()):
        neighbor_cellIds = vtk.vtkIdList()
        pd.GetPointCells( iPt, neighbor_cellIds )
        if neighbor_cellIds.GetNumberOfIds() < 3:
            continue
        # sort the neighbor_cellIds into a ring around iPt
        sorted_neighbor_cellIds = [ neighbor_cellIds.GetId( 0 ) ]
        for it in range( neighbor_cellIds.GetNumberOfIds() - 1 ):
            for iicell in range( 1, neighbor_cellIds.GetNumberOfIds() ):
                icell = neighbor_cellIds.GetId( iicell )
                if icell in sorted_neighbor_cellIds:
                    continue
                # does this cell share exactly two vertices with the last one in the list?
                if getNumberOfPointsSharedByTwoCells( pd, sorted_neighbor_cellIds[-1], icell ) == 2:
                    sorted_neighbor_cellIds += [ icell ]
                    break
        if len( sorted_neighbor_cellIds ) < neighbor_cellIds.GetNumberOfIds():
            continue # was a boundary vertex or non-manifold
        if not getNumberOfPointsSharedByTwoCells( pd, sorted_neighbor_cellIds[-1], sorted_neighbor_cellIds[0] ) == 2:
            continue # boundary vertex, in the case where cell id 0 was on the boundary
        # make a face around this vertex: a new point at each centroid of the neighboring cells
        cells.InsertNextCell( neighbor_cellIds.GetNumberOfIds() )
        for id in sorted_neighbor_cellIds:
            # find centroid of this cell
            neighbor_verts = vtk.vtkIdList()
            pd.GetCellPoints( id, neighbor_verts )
            c = (0,0,0)
            for iiv in range(neighbor_verts.GetNumberOfIds()):
                iv = neighbor_verts.GetId(iiv)
                p = pd.GetPoint( iv )
                c = add( c, p )
            c = mul( c, 1.0 / neighbor_verts.GetNumberOfIds() )
            # insert the centroid as a point and as an index into the new face
            cells.InsertCellPoint( points.InsertNextPoint( c ) )
    dual_pd = vtk.vtkPolyData()
    dual_pd.SetPoints( points )
    dual_pd.SetPolys( cells )

    # merge duplicate points
    cleaner = vtk.vtkCleanPolyData()
    if vtk.vtkVersion.GetVTKMajorVersion() >= 6:
        cleaner.SetInputData( dual_pd )
    else:
        cleaner.SetInput( dual_pd )
    cleaner.SetTolerance(0.0001)
    cleaner.Update()
    return cleaner.GetOutput()

 def GetConnectedVertices(self, connectedVerticesIDList, polyData, pointID):
     # Return IDs of all the vertices that compose the first neighbor.
     cellList = vtk.vtkIdList()
     connectedVerticesIDList.InsertUniqueId(pointID)
     # Get cells that vertex 'pointID' belongs to
     polyData.GetPointCells(pointID, cellList)
     numberOfIds = cellList.GetNumberOfIds()
     for i in range(0, numberOfIds):
         # Get points which compose all cells
         pointIdList = vtk.vtkIdList()
         polyData.GetCellPoints(cellList.GetId(i), pointIdList)
         for j in range(0, pointIdList.GetNumberOfIds()):
             connectedVerticesIDList.InsertUniqueId(pointIdList.GetId(j))
     return connectedVerticesIDList

 def __get_common_cells(self,p1,p2,exclude=-1):
     """
     Get all the cells containing p1 and p2
     """
     cell_list_1 = vtk.vtkIdList()
     self.__vtk_model.GetPointCells(p1,cell_list_1)
     nb_ids = cell_list_1.GetNumberOfIds()
     cell_list_1 = [cell_list_1.GetId(j) for j in xrange(nb_ids)]
     cell_list_2 = vtk.vtkIdList()
     self.__vtk_model.GetPointCells(p2,cell_list_2)
     nb_ids = cell_list_2.GetNumberOfIds()
     cell_list_2 = [cell_list_2.GetId(j) for j in xrange(nb_ids)]
     common_cells = [cell for cell in cell_list_1 if cell in cell_list_2 and cell!=exclude]
     return common_cells

 def __init__(self):
     vmtkSeedSelector.__init__(self)
     self.PickedSeedIds = vtk.vtkIdList()
     self.PickedSeeds = vtk.vtkPolyData()
     self.vmtkRenderer = None
     self.OwnRenderer = 0
     self.Script = None

def compute_max_array_along_lines(pd, array_name, output_array_name):
    lines = pd.GetLines()
    point_array = pd.GetPointData().GetArray(array_name)
    
    point_array_max = vtk.vtkFloatArray()
    point_array_lines_list = list()
    point_array_max_list = list()
    lines.InitTraversal()
    for lidx in range(0, pd.GetNumberOfLines()):
        if (lidx % 100) == 0:
            print lidx, '/', pd.GetNumberOfLines()
        pts = vtk.vtkIdList()       
        lines.GetNextCell(pts)
        # compute mean POINT_ARRAY for this line
        if pts.GetNumberOfIds():
            point_array_list = list()
            for pidx in range(0, pts.GetNumberOfIds()):
                point_array_list.append(point_array.GetTuple1(pts.GetId(pidx)))
            point_array_max.InsertNextTuple1(numpy.max(numpy.array(point_array_list)))
            #fa_max.InsertNextTuple1(numpy.median(numpy.array(fa_list)))
        else:
            point_array_max.InsertNextTuple1(0.0)
            
    point_array_max.SetName(output_array_name)

    outpd = pd
    outpd.GetCellData().AddArray(point_array_max)
    outpd.GetCellData().SetActiveScalars(output_array_name)

    return outpd

 def createCurtain( self, **args ):
     trajectory_points = self.getTrajectoryPoints( **args )
     extent =  self.input().GetExtent() 
     spacing =  self.input().GetSpacing() 
     nStrips = extent[5] - extent[4] 
     zmax = spacing[2] * nStrips
     z_inc = zmax / nStrips
     polydata = vtk.vtkPolyData()
     stripArray = vtk.vtkCellArray()
     stripData = [ vtk.vtkIdList() for ix in range(nStrips) ] 
     points = vtk.vtkPoints() 
     for iPt in range( trajectory_points.GetNumberOfPoints() ):
         pt_coords = trajectory_points.GetPoint( iPt )
         z = 0.0
         for iLevel in range( nStrips ):
             vtkId = points.InsertNextPoint( pt_coords[0], pt_coords[1], z )
             sd = stripData[ iLevel ]
             sd.InsertNextId( vtkId )               
             sd.InsertNextId( vtkId+1 )
             z = z + z_inc 
         points.InsertNextPoint( pt_coords[0], pt_coords[1], z )
                    
     for strip in stripData:
         stripArray.InsertNextCell(strip)
         
     polydata.SetPoints( points )
     polydata.SetStrips( stripArray )
     return polydata

def cpnonzerocolor(data):
    CellIdList = vtk.vtkIdList()
    CellIdList.Reset()
    colors = data.GetCellData().GetScalars()
    count = 0
    # Get indices from nonzero colored cells
    for i in range(data.GetNumberOfCells()):
        if ( colors.GetTuple(i) != (0.0, 0.0, 0.0, 0.0) ):
            CellIdList.InsertId(count, i)
            count = count + 1
#            print "Number Of cell:\t",i
#            print "Number of cell added:\t",count
        else:
            pass
    # Extract cells with nonzero color from data
    extractor = vtk.vtkExtractCells()
    extractor.SetInputData(data)
    extractor.SetCellList(CellIdList)
    extractor.Modified()
    extractor.Update()
    print "Number of nodes in clipped subvolume:\t",extractor.GetOutput().GetNumberOfPoints()
    print "Number of egdes in clipped subvolume:\t",extractor.GetOutput().GetNumberOfCells()
    # rearrange the id's of the cells and points - consecutive increasing ids

    print "...extracted zerocolor cells."
    
    return extractor.GetOutput()

def pd_to_array(inpd, dims=225):
    count_vol = numpy.ndarray([dims,dims,dims])
    ptids = vtk.vtkIdList()
    points = inpd.GetPoints()
    data_vol = []    
    # check for cell data
    cell_data = inpd.GetCellData().GetScalars()
    if cell_data:
        data_vol = numpy.ndarray([dims,dims,dims])
    # loop over lines
    inpd.GetLines().InitTraversal()
    print "<filter.py> Input number of points: ",\
        points.GetNumberOfPoints(),\
        "lines:", inpd.GetNumberOfLines() 
    # loop over all lines
    for lidx in range(0, inpd.GetNumberOfLines()):
        # progress
        #if verbose:
        #    if lidx % 1 == 0:
        #        print "<filter.py> Line:", lidx, "/", inpd.GetNumberOfLines()
        inpd.GetLines().GetNextCell(ptids)
        num_points = ptids.GetNumberOfIds()
        for pidx in range(0, num_points):
            point = points.GetPoint(ptids.GetId(pidx))
            # center so that 0,0,0 moves to 100,100,100
            point = numpy.round(numpy.array(point) + 110)         
            count_vol[point[0], point[1], point[2]] += 1
            if cell_data:
                data_vol[point[0], point[1], point[2]] += cell_data.GetTuple(lidx)[0]
    return count_vol, data_vol

def testTwoTrianglesCoplanar():

    "Two triangles"

    # create set of points
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(4)
    points.SetPoint(0, [0., 0., 0.])
    points.SetPoint(1, [1., 0., 0.])
    points.SetPoint(2, [0., 1., 0.])
    points.SetPoint(3, [1., 1., 0.])

    # create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)

    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)

    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 1)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    for order in range(1, 6):
        lslm = PoissonSolver(pdata,
                             max_edge_length=1000.,
                             order=order)
        print('order = ', order)
        print('g matrix: ', lslm.getGreenMatrix())

def testStartingWithTwoCells():

    points = vtk.vtkPoints()
    points.InsertNextPoint((0., 0., 0.))
    points.InsertNextPoint((2., 0., 0.))
    points.InsertNextPoint((2., 1., 0.))
    points.InsertNextPoint((0., 1., 0.))

    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)

    pdata.Allocate(2, 1)
    ptIds = vtk.vtkIdList()
    ptIds.SetNumberOfIds(3)
    ptIds.SetId(0, 0)
    ptIds.SetId(1, 1)
    ptIds.SetId(2, 2)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)
    ptIds.SetId(0, 2)
    ptIds.SetId(1, 3)
    ptIds.SetId(2, 0)
    pdata.InsertNextCell(vtk.VTK_POLYGON, ptIds)

    rs = RefineSurface(pdata)
    rs.refine(max_edge_length=1.1)
    assert(rs.getVtkPolyData().GetNumberOfPolys() == 8)

 def ReadTecplotSurfaceFile(self):
     self.PrintLog('Reading Tecplot surface file.')
     f=open(self.InputFileName, 'r')
     line = f.readline()
     if line.split()[0] == 'TITLE':
         line = f.readline()
     if (line.split()[0] == 'VARIABLES') | (line.split('=')[0] == 'VARIABLES'):
         arrayNames = line.split('=')[1].strip().split(',')
         arrayNames[0:3] = []
         self.PrintLog("ArrayNames" + str(arrayNames))
         line = f.readline()
     if line.split()[0] == 'ZONE':
         lineNid = line.find('N=')
         lineN = line[lineNid : lineNid+line[lineNid:].find(',') ].split('=')[1]
         numberOfNodes = int(lineN)
         lineEid = line.find('E=')
         lineE = line[lineEid : lineEid+line[lineEid:].find(',') ].split('=')[1]
         numberOfElements = int(lineE)
         elementType = 'TRIANGLE'
         if line.find('ET=') != -1:
             if 'TRIANGLE' in line:
                 elementType = 'TRIANGLE'
             elif 'QUADRILATERAL' in line:
                 elementType = 'QUADRILATERAL'
     self.PrintLog("Reading " + str(numberOfNodes)+" nodes.")
     points = vtk.vtkPoints()
     cells = vtk.vtkCellArray()
     points.SetNumberOfPoints(numberOfNodes)
     self.Surface = vtk.vtkPolyData()
     self.Surface.SetPoints(points)
     self.Surface.SetPolys(cells)
     for arrayName in arrayNames:
         array = vtk.vtkDoubleArray()
         array.SetName(arrayName)
         array.SetNumberOfTuples(numberOfNodes)
         self.Surface.GetPointData().AddArray(array)
     self.Surface.Update()
     data = f.read().split()
     dataCounter = 0
     for i in range(numberOfNodes):
         point = [float(data[dataCounter]),float(data[dataCounter+1]),float(data[dataCounter+2])]
         dataCounter += 3
         points.SetPoint(i,point)
         for j in range(len(arrayNames)):
             self.Surface.GetPointData().GetArray(arrayNames[j]).SetComponent(i,0,float(data[dataCounter]))
             dataCounter += 1
     self.PrintLog("Reading " + str(numberOfElements)+" elements.")
     cellIds = vtk.vtkIdList()
     for i in range(numberOfElements):
         cellIds.Initialize()
         cellIds.InsertNextId(int(data[dataCounter])-1)
         dataCounter += 1
         cellIds.InsertNextId(int(data[dataCounter])-1)
         dataCounter += 1
         cellIds.InsertNextId(int(data[dataCounter])-1)
         dataCounter += 1
         if elementType == "QUADRILATERAL":
             cellIds.InsertNextId(int(data[dataCounter])-1)
             dataCounter += 1
         cells.InsertNextCell(cellIds)

 def setSourceFromExpression(self, expression):
     """
     Set the source from expression
     @param expression expression of x, y, and z
     """
     from math import sqrt, pi, sin, cos, tan, log, exp
     
     n = self.pdata.GetNumberOfPolys()
     sourceData = vtk.vtkDoubleArray()
     sourceData.SetNumberOfComponents(1)
     sourceData.SetNumberOfTuples(n)
     sourceData.SetName(self.sourceName)
     midPoint = numpy.zeros((3,), numpy.float64)
     ptIds = vtk.vtkIdList()
     cells = self.pdata.GetPolys()
     cells.InitTraversal()
     for i in range(n):
         cell = cells.GetNextCell(ptIds)
         npts = ptIds.GetNumberOfIds()
         midPoint *= 0 # reset
         for j in range(npts):
             midPoint += self.points.GetPoint(ptIds.GetId(j))
         midPoint /= float(npts)
         x, y, z = midPoint
         v = eval(expression)
         sourceData.SetTuple(i, [v])
     self.pdata.GetCellData().AddArray(sourceData)