def addPlot(self,_plotName,_style="Lines"):   # called directly from Steppable; add a (possibly more than one) plot to a plot window

        self.plotWindowInterfaceMutex.lock()
#        self.plotWindowMutex.lock()

#        return
#        print MODULENAME,'   addPlot():  _plotName= ',_plotName
#        import pdb; pdb.set_trace()
        
#        self.plotData[_plotName] = [array([],dtype=double),array([],dtype=double),False]  # 'array': from PyQt4.Qwt5.anynumpy import *
        
        self.chart = vtk.vtkChartXY()
#        self.chart.GetAxis(vtk.vtkAxis.LEFT).SetLogScale(True)
#        self.chart.GetAxis(vtk.vtkAxis.BOTTOM).SetLogScale(True)
#        self.numCharts += 1
        self.plotData[_plotName] = [self.chart]

        self.view = vtk.vtkContextView()
        self.ren = self.view.GetRenderer()
#        self.renWin = self.qvtkWidget.GetRenderWindow()
        self.renWin = self.pW.GetRenderWindow()
        self.renWin.AddRenderer(self.ren)

        # Create a table with some points in it
        self.table = vtk.vtkTable()

        self.arrX = vtk.vtkFloatArray()
        self.arrX.SetName("xarray")

        self.arrC = vtk.vtkFloatArray()
        self.arrC.SetName("yarray")

        numPoints = 5
        numPoints = 15
        inc = 7.5 / (numPoints - 1)

#        for i in range(0,numPoints):
#            self.arrX.InsertNextValue(i*inc)
#            self.arrC.InsertNextValue(math.cos(i * inc) + 0.0)

#        self.arrX.InsertNextValue(0.0)
#        self.arrC.InsertNextValue(0.0)
#        self.arrX.InsertNextValue(0.1)
#        self.arrC.InsertNextValue(0.1)

        self.table.AddColumn(self.arrX)
        self.table.AddColumn(self.arrC)

        # Now add the line plots with appropriate colors
        self.line = self.chart.AddPlot(0)
        self.line.SetInput(self.table,0,1)
        self.line.SetColor(0,0,255,255)
        self.line.SetWidth(1.0)


        self.view.GetRenderer().SetBackground([0.6,0.6,0.1])
        self.view.GetRenderer().SetBackground([1.0,1.0,1.0])
        self.view.GetScene().AddItem(self.chart)

        self.plotWindowInterfaceMutex.unlock()

    def removeOldGeometry(self, fileName):
        if fileName is None:
            self.grid = vtk.vtkUnstructuredGrid()
            self.gridResult = vtk.vtkFloatArray()
            #self.emptyResult = vtk.vtkFloatArray()
            #self.vectorResult = vtk.vtkFloatArray()
            self.grid2 = vtk.vtkUnstructuredGrid()
            self.scalarBar.VisibilityOff()
            skipReading = True
        else:
            self.TurnTextOff()
            self.grid.Reset()
            self.grid2.Reset()
            self.gridResult = vtk.vtkFloatArray()
            self.gridResult.Reset()
            self.gridResult.Modified()
            self.eidMap = {}
            self.nidMap = {}

            self.resultCases = {}
            self.nCases = 0
            try:
                del self.caseKeys
                del self.iCase
                del self.iSubcaseNameMap
            except:
                print "cant delete geo"
                pass
            ###
            #print dir(self)
            skipReading = False
        self.scalarBar.VisibilityOff()
        self.scalarBar.Modified()
        return skipReading

def ndarray_to_vtkarray(colors, radius, nbat):
    # define the colors
    color_scalars = vtk.vtkFloatArray()
    for c in colors:
            color_scalars.InsertNextValue(c)
    color_scalars.SetName("colors")
    
    # some radii
    radius_scalars = vtk.vtkFloatArray()
    for r in radius:
        radius_scalars.InsertNextValue(r)
    radius_scalars.SetName("radius")
    
    # the original index
    index_scalars = vtk.vtkIntArray()
    for i in range(nbat):
        index_scalars.InsertNextValue(i)
    radius_scalars.SetName("index")
    
    scalars = vtk.vtkFloatArray()
    scalars.SetNumberOfComponents(3)
    scalars.SetNumberOfTuples(radius_scalars.GetNumberOfTuples())
    scalars.CopyComponent(0, radius_scalars ,0 )
    scalars.CopyComponent(1, color_scalars ,0 )
    scalars.CopyComponent(2, index_scalars ,0 )
    scalars.SetName("scalars")
    return scalars 

def BlenderToPolyData(me, uvlayer=None):

	pcoords = vtk.vtkFloatArray()
	pcoords.SetNumberOfComponents(3)
	pcoords.SetNumberOfTuples(len(me.verts))
	for i in range(len(me.verts)):
		p0 = me.verts[i].co[0]
		p1 = me.verts[i].co[1]
		p2 = me.verts[i].co[2]
		pcoords.SetTuple3(i, p0, p1, p2)

	points = vtk.vtkPoints()
	points.SetData(pcoords)

	polys = vtk.vtkCellArray()
	lines = vtk.vtkCellArray()
	for face in me.faces:
		if len(face.v) == 4:
			polys.InsertNextCell(4)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
			polys.InsertCellPoint(face.v[3].index)
		elif len(face.v) == 3:
			polys.InsertNextCell(3)
			polys.InsertCellPoint(face.v[0].index)
			polys.InsertCellPoint(face.v[1].index)
			polys.InsertCellPoint(face.v[2].index)
		elif len(face.v) == 2:
			lines.InsertNextCell(2)
			lines.InsertCellPoint(face.v[0].index)
			lines.InsertCellPoint(face.v[1].index)

	for edge in me.edges:
		lines.InsertNextCell(2)
		lines.InsertCellPoint(edge.v1.index)
		lines.InsertCellPoint(edge.v2.index)

	pdata =vtk.vtkPolyData()
	pdata.SetPoints(points)
	pdata.SetPolys(polys)
	pdata.SetLines(lines)

	if me.faceUV:
		if uvlayer:
			uvnames = me.getUVLayerNames()
			if uvlayer in uvnames:
				me.activeUVLayer = uvlayer
		tcoords = vtk.vtkFloatArray()
		tcoords.SetNumberOfComponents(2)
		tcoords.SetNumberOfTuples(len(me.verts))
		for face in me.faces:
			for i in range(len(face.verts)):
				uv = face.uv[i]
				tcoords.SetTuple2(face.v[i].index, uv[0], uv[1])
		pdata.GetPointData().SetTCoords(tcoords);

	pdata.Update()

	return pdata

    def testLinePlot(self):
        "Test if line plots can be built with python"

        # Set up a 2D scene, add an XY chart to it
        view = vtk.vtkContextView()
        view.GetRenderer().SetBackground(1.0,1.0,1.0)
        view.GetRenderWindow().SetSize(400,300)
        chart = vtk.vtkChartXY()
        view.GetScene().AddItem(chart)

        # Create a table with some points in it
        table = vtk.vtkTable()

        arrX = vtk.vtkFloatArray()
        arrX.SetName("X Axis")

        arrC = vtk.vtkFloatArray()
        arrC.SetName("Cosine")

        arrS = vtk.vtkFloatArray()
        arrS.SetName("Sine")

        arrS2 = vtk.vtkFloatArray()
        arrS2.SetName("Sine2")

        numPoints = 69
        inc = 7.5 / (numPoints - 1)

        for i in range(0,numPoints):
            arrX.InsertNextValue(i*inc)
            arrC.InsertNextValue(math.cos(i * inc) + 0.0)
            arrS.InsertNextValue(math.sin(i * inc) + 0.0)
            arrS2.InsertNextValue(math.sin(i * inc) + 0.5)

        table.AddColumn(arrX)
        table.AddColumn(arrC)
        table.AddColumn(arrS)
        table.AddColumn(arrS2)

        # Now add the line plots with appropriate colors
        line = chart.AddPlot(0)
        line.SetInput(table,0,1)
        line.SetColor(0,255,0,255)
        line.SetWidth(1.0)

        line = chart.AddPlot(0)
        line.SetInput(table,0,2)
        line.SetColor(255,0,0,255);
        line.SetWidth(5.0)

        line = chart.AddPlot(0)
        line.SetInput(table,0,3)
        line.SetColor(0,0,255,255);
        line.SetWidth(4.0)

        view.GetRenderWindow().SetMultiSamples(0)

        img_file = "TestLinePlot.png"
        vtk.test.Testing.compareImage(view.GetRenderWindow(),vtk.test.Testing.getAbsImagePath(img_file),threshold=25)
        vtk.test.Testing.interact()

def ndarray_to_vtkarray(colors, radius, nbat):
    # define the colours
    color_scalars = vtk.vtkFloatArray()
    color_scalars.SetNumberOfValues(colors.shape[0])
    for i,c in enumerate(colors):
            color_scalars.SetValue(i,c)
    color_scalars.SetName("colors")
    
    # some radii
    radius_scalars = vtk.vtkFloatArray()
    radius_scalars.SetNumberOfValues(radius.shape[0])
    for i,r in enumerate(radius):
        radius_scalars.SetValue(i,r)
    radius_scalars.SetName("radius")
    
    # the original index
    index_scalars = vtk.vtkIntArray()
    index_scalars.SetNumberOfValues(nbat)
    for i in range(nbat):
        index_scalars.SetValue(i,i)
    index_scalars.SetName("index")
    
    scalars = vtk.vtkFloatArray()
    scalars.SetNumberOfComponents(3)
    scalars.SetNumberOfTuples(radius_scalars.GetNumberOfTuples())
    scalars.CopyComponent(0, radius_scalars ,0 )
    scalars.CopyComponent(1, color_scalars ,0 )
    scalars.CopyComponent(2, index_scalars ,0 )
    scalars.SetName("scalars")
    return scalars 

    def testBufferShared(self):
        """Test the special buffer_shared() check that VTK provides."""
        a = bytearray(b'hello')
        self.assertEqual(vtk.buffer_shared(a, a), True)
        b = bytearray(b'hello')
        self.assertEqual(vtk.buffer_shared(a, b), False)

        a = vtk.vtkFloatArray()
        a.SetNumberOfComponents(3)
        a.InsertNextTuple((10, 7, 4))
        a.InsertNextTuple((85, 8, 2))

        b = vtk.vtkFloatArray()
        b.SetVoidArray(a, 6, True)
        self.assertEqual(vtk.buffer_shared(a, b), True)

        c = vtk.vtkFloatArray()
        c.DeepCopy(a)
        self.assertEqual(vtk.buffer_shared(a, c), False)

        if sys.hexversion >= 0x02070000:
            m = memoryview(a)
            self.assertEqual(vtk.buffer_shared(a, m), True)

        if sys.hexversion < 0x03000000:
            m = buffer(a)
            self.assertEqual(vtk.buffer_shared(a, m), True)

 def get_vtk_data(self,geo):
     """Returns dictionary of VTK data arrays from rock types.  The geometry object geo must be passed in."""
     from vtk import vtkIntArray,vtkFloatArray,vtkCharArray
     arrays={'Block':{'Rock type index':vtkIntArray(),'Porosity':vtkFloatArray(),
                      'Permeability':vtkFloatArray(),'Name':vtkCharArray()},'Node':{}}
     vector_properties=['Permeability']
     string_properties=['Name']
     string_length=5
     nele=geo.num_underground_blocks
     array_length={'Block':nele,'Node':0}
     for array_type,array_dict in arrays.items():
         for name,array in array_dict.items():
             array.SetName(name)
             if name in vector_properties:
                 array.SetNumberOfComponents(3)
                 array.SetNumberOfTuples(array_length[array_type])
             elif name in string_properties:
                 array.SetNumberOfComponents(string_length)
                 array.SetNumberOfTuples(array_length[array_type])
             else: 
                 array.SetNumberOfComponents(1)
                 array.SetNumberOfValues(array_length[array_type])
     natm=geo.num_atmosphere_blocks
     rindex=self.rocktype_indices[natm:]
     for i,ri in enumerate(rindex):
         arrays['Block']['Rock type index'].SetValue(i,ri)
         rt=self.rocktypelist[ri]
         arrays['Block']['Porosity'].SetValue(i,rt.porosity)
         k=rt.permeability
         arrays['Block']['Permeability'].SetTuple3(i,k[0],k[1],k[2])
     for i,blk in enumerate(self.blocklist[natm:]):
         arrays['Block']['Name'].SetTupleValue(i,blk.name)
     return arrays

  def section_CreateTable(self):
    self.delayDisplay("Create table",self.delayMs)

    tableNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode", self.tableName)
    self.assertIsNotNone(tableNode)
    table = tableNode.GetTable()
    self.assertIsNotNone(table)

    # Create X, Y1, and Y2 series

    arrX = vtk.vtkFloatArray()
    arrX.SetName(self.xColumnName)
    table.AddColumn(arrX)

    arrY1 = vtk.vtkFloatArray()
    arrY1.SetName(self.y1ColumnName)
    table.AddColumn(arrY1)

    arrY2 = vtk.vtkFloatArray()
    arrY2.SetName(self.y2ColumnName)
    table.AddColumn(arrY2)

    # Fill in the table with some example values
    import math
    numPoints = 69
    inc = 7.5 / (numPoints - 1)
    table.SetNumberOfRows(numPoints)
    for i in range(numPoints):
      table.SetValue(i, 0, i * inc )
      table.SetValue(i, 1, math.cos(i * inc))
      table.SetValue(i, 2, math.sin(i * inc))

    def findLocalCsys(self, filename):
        poly = vtk.vtkGeometryFilter()
        poly.SetInputData(self.vtkMeshes[filename])
        poly.Update()

        distanceFilter = vtk.vtkImplicitPolyDataDistance()
        distanceFilter.SetInput(poly.GetOutput())

        gradients = vtk.vtkFloatArray()
        gradients.SetNumberOfComponents(3)
        gradients.SetName("LevelSet Normals")

        distances = vtk.vtkFloatArray()
        distances.SetNumberOfComponents(1)
        distances.SetName("Signed Distances")

        N = self.vtkMeshes[filename].GetCellData().GetArray(
            "Centroids").GetNumberOfTuples()
        for i in range(N):
            g = np.zeros(3, np.float32)
            p = self.vtkMeshes[filename].GetCellData().GetArray(
                "Centroids").GetTuple3(i)
            d = distanceFilter.EvaluateFunction(p)
            distanceFilter.EvaluateGradient(p, g)
            g = old_div(np.array(g), np.linalg.norm(np.array(g)))
            gradients.InsertNextTuple(g)
            distances.InsertNextValue(d)
        self.vtkMeshes[filename].GetCellData().AddArray(gradients)
        self.vtkMeshes[filename].GetCellData().AddArray(distances)

def transferTCoords():
	input = paf.GetInputDataObject(0,0)
	refin = paf.GetInputList().GetItem(0)
	output = paf.GetPolyDataOutput()
	# output.CopyStructure(input)
	# output.CopyAttributes(input)
	
	TCorig = refin.GetPointData().GetTCoords()
	Norig = refin.GetPointData().GetNormals()
	
	TC = vtk.vtkFloatArray()
	TC.SetNumberOfComponents(TCorig.GetNumberOfComponents())
	TC.SetNumberOfTuples(TCorig.GetNumberOfTuples())
	TC.SetName('Texture Coordinates')
	for ii in range(TCorig.GetNumberOfTuples()):
		ff = TCorig.GetTuple2(ii)
		TC.SetTuple2(ii,ff[0],ff[1])
	
	NN = vtk.vtkFloatArray()
	NN.SetNumberOfComponents(Norig.GetNumberOfComponents())
	NN.SetNumberOfTuples(Norig.GetNumberOfTuples())
	NN.SetName('normals')
	for ii in range(Norig.GetNumberOfTuples()):
		ff = Norig.GetTuple3(ii)
		NN.SetTuple3(ii,ff[0],ff[1],ff[2])
		
	output.GetPointData().AddArray(TC)
	output.GetPointData().SetActiveTCoords('Texture Coordinates')
	output.GetPointData().AddArray(NN)
	output.GetPointData().SetActiveNormals('normals')

    def generateMesh(self, x, y, z, in_x, in_y, in_z):
        x_coord = vtk.vtkFloatArray()
        x_coord.InsertNextValue(x)

        y_coord = vtk.vtkFloatArray()
        y_coord.InsertNextValue(y)

        z_coord = vtk.vtkFloatArray()
        z_coord.InsertNextValue(z)

        grid = vtk.vtkRectilinearGrid()
        grid.SetDimensions(self.nx if in_x else 1, self.ny if in_y else 1, self.nz if in_z else 1)
        grid.SetXCoordinates(self.x_coords if in_x else x_coord);
        grid.SetYCoordinates(self.y_coords if in_y else y_coord);
        grid.SetZCoordinates(self.z_coords if in_z else z_coord);

        # We're going to generate all of the IDs of the cells in that rectilinear grid so we can extract them as an UnstructuredGrid
        # Why would we do such a thing?  Because there is some bug associated with RecitilinearGrids and clipping / rendering
        # So, instead I'm going to use RectilinearGrid for generating the cells and then "copy" it to an UnstructuredGrid using ExtractCells
        num_cells = grid.GetNumberOfCells()
        id_list = vtk.vtkIdList()
        for i in xrange(num_cells):
            id_list.InsertNextId(i)

        extract = vtk.vtkExtractCells()
        if vtk.VTK_MAJOR_VERSION <= 5:
            extract.SetInput(grid)
        else:
            extract.SetInputData(grid)
        extract.SetCellList(id_list)

        return extract.GetOutput()

def compute_scalar_measures(pd):
    tensors = pd.GetPointData().GetTensors()
    lines = pd.GetLines()

    fa = vtk.vtkFloatArray()
    lambda_parallel = vtk.vtkFloatArray()
    lambda_perp = vtk.vtkFloatArray()
    #fa_array = list()
    
    for tidx in range(tensors.GetNumberOfTuples()):
        if (tidx % 5000) == 0:
            print tidx, '/', tensors.GetNumberOfTuples()
        D = tensors.GetTuple9(tidx)
        D = numpy.array(D).reshape(3,3)
        w, v = numpy.linalg.eig(D)
        fa_tidx = fractional_anisotropy(w)
        lambda_parallel_tidx = w[2]
        lambda_perp_tidx = (w[0] + w[1])/2.0
        fa.InsertNextTuple1(fa_tidx)
        lambda_parallel.InsertNextTuple1(lambda_parallel_tidx)
        lambda_perp.InsertNextTuple1(lambda_perp_tidx)

    fa_avg = vtk.vtkFloatArray()
    fa_lines_list = list()
    fa_avg_list = list()
    lines.InitTraversal()
    for lidx in range(0, pd.GetNumberOfLines()):
        if (lidx % 100) == 0:
            print lidx, '/', pd.GetNumberOfLines()
        pts = vtk.vtkIdList()       
        #lines.GetCell(lidx, pts)
        lines.GetNextCell(pts)
        # compute average FA for this line
        if pts.GetNumberOfIds():
            fa_list = list()
            for pidx in range(0, pts.GetNumberOfIds()):
                fa_list.append(fa.GetTuple1(pts.GetId(pidx)))
            #fa_avg.InsertNextTuple1(numpy.mean(numpy.array(fa_list)))
            fa_avg.InsertNextTuple1(numpy.median(numpy.array(fa_list)))
            fa_lines_list.append(fa_list)
            fa_avg_list.append(numpy.median(numpy.array(fa_list)))
        else:
            fa_avg.InsertNextTuple1(0.0)
            fa_lines_list.append([0.0])
            
    fa_avg.SetName('mean_FA')
    fa.SetName('FA')
    lambda_parallel.SetName('parallel_diffusivity')
    lambda_perp.SetName('perpendicular_diffusivity')

    outpd = pd
    outpd.GetCellData().AddArray(fa_avg)
    outpd.GetCellData().SetActiveScalars('mean_FA')
    outpd.GetPointData().AddArray(fa)
    outpd.GetPointData().AddArray(lambda_parallel)
    outpd.GetPointData().AddArray(lambda_perp)
    outpd.GetPointData().SetActiveScalars('FA')

    return outpd, fa_lines_list, fa_avg_list

def writeVTR(vtr_name, scalar_fields, vector_fields, vtkX, vtkY, vtkZ, localZrange):
    """Writes a single VTR file per Python processor/variable

    Parameters:
     vtr_name - name of the VTR file
     scalar_fields - dictionary with scalar field arrays ordered [x, y, z], e.g. {'p': array[nx,ny,nz], 'rho0': array[nx,ny,nz]}
     vector_fields - dictionary with vector fields ordered [3, x, y, z], e.g. {'J': array[3,nx,ny,nz], 'B': array[3,nx,ny,nz]}
     vtkX, vtkY, vtkZ - VTR coordinates, see createVtrCoordinates()
     localZrange - local range for Z indices

    """

    Nx = vtkX.GetNumberOfTuples() - 1
    Ny = vtkY.GetNumberOfTuples() - 1
    Nz=0 #2D
    numpoints = (Nx+1)*(Ny+1)*(Nz+1)
    rtg = vtk.vtkRectilinearGrid()
    rtg.SetExtent(0, Nx, 0, Ny, localZrange[0], localZrange[1]) 
    rtg.SetXCoordinates(vtkX)
    rtg.SetYCoordinates(vtkY)
    vtk_data = []
    array_list = []
    for f in scalar_fields:
        vtk_data.append(vtk.vtkFloatArray())
        vtk_data[-1].SetNumberOfTuples(numpoints)
        vtk_data[-1].SetNumberOfComponents(1)
        array_list.append(scalar_fields[f].swapaxes(0,2).flatten().astype("float32"))
        vtk_data[-1].SetVoidArray(array_list[-1], numpoints, 1)
        vtk_data[-1].SetName(f)
        if f == scalar_fields.keys()[0]:
            rtg.GetPointData().SetScalars(vtk_data[-1])
        else:
            rtg.GetPointData().AddArray(vtk_data[-1])
    #end for
    for f in vector_fields:
        vtk_data.append(vtk.vtkFloatArray())
        vtk_data[-1].SetNumberOfTuples(numpoints*3)
        vtk_data[-1].SetNumberOfComponents(3)
        array_list.append(vector_fields[f].swapaxes(0,3).swapaxes(1,2).flatten().astype("float32"))
        vtk_data[-1].SetVoidArray(array_list[-1], numpoints*3, 1)
        vtk_data[-1].SetName(f)
        if f == vector_fields.keys()[0]:
            rtg.GetPointData().SetVectors(vtk_data[-1])
        else:
            rtg.GetPointData().AddArray(vtk_data[-1])
    #end for
    try:
        writer = vtk.vtkXMLRectilinearGridWriter()                                                                                              
        writer.SetFileName(vtr_name)
        writer.SetInput(rtg)
        writer.Write()
    except:
        print 'Error writing VTR file: ', vtr_name

    def __init__(self, volume, level=None):
        self._surface_algorithm = None
        self._renderer = None
        self._actor = None
        self._mapper = None
        self._volume_array = None

        self._float_array = _vtk.vtkFloatArray()
        self._image_data = _vtk.vtkImageData()
        self._image_data.GetPointData().SetScalars(self._float_array)
        self._setup_data(_numpy.float32(volume))

        self._surface_algorithm = _vtk.vtkMarchingCubes()
        self._surface_algorithm.SetInputData(self._image_data)
        self._surface_algorithm.ComputeNormalsOn()

        if level is not None:
            try:
                self.set_multiple_levels(iter(level))
            except TypeError:
                self.set_level(0, level)

        self._mapper = _vtk.vtkPolyDataMapper()
        self._mapper.SetInputConnection(self._surface_algorithm.GetOutputPort())
        self._mapper.ScalarVisibilityOn() # new
        self._actor = _vtk.vtkActor()
        self._actor.SetMapper(self._mapper)

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 array2vtk(z):    
    """Converts a numpy Array to a VTK array object directly. The
    resulting array copies the data in the passed  array.  The
    array can therefore be deleted safely.  This works for real arrays.
    """ 
    arr_vtk = {'c':vtkConstants.VTK_UNSIGNED_CHAR,
               'b':vtkConstants.VTK_UNSIGNED_CHAR,
               '1':vtkConstants.VTK_CHAR,
               's':vtkConstants.VTK_SHORT,
               'i':vtkConstants.VTK_INT,
               'l':vtkConstants.VTK_LONG,
               'f':vtkConstants.VTK_FLOAT,
               'd':vtkConstants.VTK_DOUBLE,
               'F':vtkConstants.VTK_FLOAT,
               'D':vtkConstants.VTK_DOUBLE }

    # A dummy array used to create others.
    f = vtk.vtkFloatArray()
    # First create an array of the right type by using the typecode.
    tmp = f.CreateDataArray(arr_vtk[z.dtype.char])
    tmp.SetReferenceCount(2) # Prevents memory leak.
    zf = N.ravel(z)
    tmp.SetNumberOfTuples(len(zf))
    tmp.SetNumberOfComponents(1)
    tmp.SetVoidArray(zf, len(zf), 1)
    # Now create a new array that is a DeepCopy of tmp.  This is
    # required because tmp does not copy the data from the NumPy array
    # and will point to garbage if the NumPy array is deleted.
    arr = f.CreateDataArray(arr_vtk[z.dtype.char])
    arr.SetReferenceCount(2) # Prevents memory leak.
    arr.DeepCopy(tmp)
    return arr

def make_spreadsheet(column_names, table):
    # column_names is a list of strings
    # table is a 2D numpy.ndarray
    # returns a vtkTable object that stores the table content

    # Create a vtkTable to store the output.
    rows = table.shape[0]

    if (table.shape[1] != len(column_names)):
        print('Warning: table number of columns differs from number of '
              'column names')
        return

    from vtk import vtkTable, vtkFloatArray
    vtk_table = vtkTable()
    for (column, name) in enumerate(column_names):
        array = vtkFloatArray()
        array.SetName(name)
        array.SetNumberOfComponents(1)
        array.SetNumberOfTuples(rows)
        vtk_table.AddColumn(array)

        for row in range(0, rows):
            array.InsertValue(row, table[row, column])

    return vtk_table

    def draw(self, renderer):
        ### FIXME this should be made faster (C++ Module? How to deal with C++ linkage problems?)
        ### Sticking the vtkPoints objects in a cache would help somewhat but not on the first view.
        ### - Jack
        if not self.drawn:
            vtk_points = vtkPoints()
            points = self.visualiser.getQuantityPoints(self.quantityName, dynamic=self.dynamic)
            nPoints = len(points)
            vtk_points.SetNumberOfPoints(nPoints)
            setPoint = vtkPoints.SetPoint
            for i in xrange(nPoints):
                z = points[i] * self.zScale + self.offset
                setPoint(vtk_points, i, self.visualiser.xPoints[i], self.visualiser.yPoints[i], z)

            polyData = vtkPolyData()
            polyData.SetPoints(vtk_points)
            polyData.SetPolys(self.visualiser.vtk_cells)
            mapper = vtkPolyDataMapper()
            mapper.SetInput(polyData)
            setValue = vtkFloatArray.SetValue
            if hasattr(self.colour[0], '__call__'):
                scalars = self.colour[0](self.visualiser.getQuantityDict())
                nScalars = len(scalars)
                vtk_scalars = vtkFloatArray()
                vtk_scalars.SetNumberOfValues(nScalars)
                for i in xrange(nScalars):
                    setValue(vtk_scalars, i, scalars[i])
                polyData.GetPointData().SetScalars(vtk_scalars)
                mapper.SetScalarRange(self.colour[1:3])
            mapper.Update()
            self.actor.SetMapper(mapper)
        Feature.draw(self, renderer)

 def pca(self):
   """Performs Principle Component Analysis on the alignedGrids. Also calculates the mean shape.
   """
   
   logging.info("running pca")
   
   size = len(self.alignedGrids)
   
   self.filterPCA.SetNumberOfInputs(size)
   
   for id, grid in enumerate(self.alignedGrids):    
     self.filterPCA.SetInput(id, grid)
   
   self.filterPCA.Update()
      
   #Get the eigenvalues
   evals = self.filterPCA.GetEvals()
   
   #Now let's get mean ^^
   b = vtk.vtkFloatArray()
   b.SetNumberOfComponents(0)
   b.SetNumberOfTuples(0)
   mean = vtk.vtkUnstructuredGrid()
   mean.DeepCopy(self.alignedGrids[0])
   #Get the mean shape:
   self.filterPCA.GetParameterisedShape(b, mean)
   self.meanShape.append(mean)   
   
   #get the meanpositions
   for pos in range(self.meanShape[0].GetNumberOfCells()):
     bounds = self.meanShape[0].GetCell(pos).GetBounds()
     self.meanPositions.append((bounds[0],bounds[2]))
     
   logging.info("done")