def __init__(self, pdata):
        """
        Constructor
        @param pdata vtkPolyData instance
        """
        self.polydata = pdata

        # save the points and point data as separate class members
        # so as not to pollute pdata
        self.points = vtk.vtkPoints()
        self.points.DeepCopy(pdata.GetPoints())

        self.pointData = {}
        pd = pdata.GetPointData()
        for i in range(pd.GetNumberOfArrays()):
            arr = pd.GetArray(i)
            name = arr.GetName()
            self.pointData[name] = vtk.vtkDoubleArray()
            self.pointData[name].DeepCopy(arr)

        self.cellData = {}
        cd = pdata.GetCellData()
        for i in range(cd.GetNumberOfArrays()):
            name = cd.GetArray(i).GetName()
            self.cellData[name] = vtk.vtkDoubleArray()
            self.cellData[name].SetName(name)

    def initVTK(self, key):
        print "initializing VTK objects"
        self.__setMinMaxRange(key)
        self.__normalize(key)

        # Create our Data object
        vtk.vtkStructuredGrid()
        self.vtkStructuredGrid = vtk.vtkStructuredGrid()
        self.vtkStructuredGrid.SetDimensions([128]*3)

        # Set up Points
        self.vtkPoints = vtk.vtkPoints()
        self.vtkPoints.Allocate(128*3)

        # Set up Cells
        # self.vtkCells = vtk.vtkCellArray()

        # Setup the velocity vectors
        self.vtkVectors = vtk.vtkDoubleArray()
        self.vtkVectors.SetNumberOfComponents(3)
        self.vtkVectors.SetNumberOfTuples(self.sdfIdent.size)

        # Setup the Scalars
        self.vtkScalars = vtk.vtkDoubleArray()
        self.vtkScalars.SetName(key)

        # Allocate points, cells, scalars, and vector fields
        self.AllocateData(self.sdfIdent.size, key)

        # Now attach newly allocated objects to the grid
        self.vtkStructuredGrid.SetPoints(self.vtkPoints)
        # self.vtkStructuredGrid.SetVerts(self.vtkCells)
        self.vtkStructuredGrid.GetPointData().SetVectors(self.vtkVectors)
        self.vtkStructuredGrid.GetPointData().SetScalars(self.vtkScalars)
        self.vtkStructuredGrid.GetPointData().SetActiveScalars(key)

def SetVtkGrid(x,y,z):
    """Set up the vtk rectilinear grid using x, y, z data.
    
    Parameters:
        x, y, z -- the points in the x, y and z directions respectively
    Returns:
        grid -- a vtkRectilinearGrid object
    """
    grid = vtkRectilinearGrid();
    grid.SetDimensions(len(x),len(y),len(z));
     
    xArray = vtkDoubleArray();
    for xCoord in x: xArray.InsertNextValue(xCoord)
     
    yArray = vtkDoubleArray();
    for yCoord in y: yArray.InsertNextValue(yCoord)
    
    zArray = vtkDoubleArray();
    for zCoord in z: zArray.InsertNextValue(zCoord)
     
    grid.SetXCoordinates(xArray);
    grid.SetYCoordinates(yArray);
    grid.SetZCoordinates(zArray);
    
    print "There are " + str(grid.GetNumberOfPoints()) + " points.";
    print "There are " + str(grid.GetNumberOfCells()) + " cells.";

    return grid

def cell_average(model, bucket):
    """ Calculate a volume fraction estimate at the level of the grid."""

    ugrid = vtk.vtkUnstructuredGrid()
    ugrid.DeepCopy(model)

    locator = vtk.vtkCellLocator()
    locator.SetDataSet(ugrid)
    locator.BuildLocator()

    volfrac = numpy.zeros(ugrid.GetNumberOfCells())
    volume = numpy.zeros(ugrid.GetNumberOfCells())
    temperature = numpy.zeros(ugrid.GetNumberOfCells())
    velocity = numpy.zeros((ugrid.GetNumberOfCells(),3))

    for particle in bucket.particles:
        cell_id = locator.FindCell(particle.pos)
        volume[cell_id] += particle.volume
        velocity[cell_id, :] += particle.volume*particle.vel

    for _ in range(ugrid.GetNumberOfCells()):
        if volume[_] >1.0e-12:
            velocity[_, :] /= volume[_]
        volfrac[_] = volume[_] / get_measure(ugrid.GetCell(_))

    for particle in bucket.particles:
        cell_id = locator.FindCell(particle.pos)
        temperature[cell_id] += particle.volume*distance2(particle.vel,velocity[cell_id, :])

    for _ in range(ugrid.GetNumberOfCells()):
         if volume[_] >1.0e-12:
             temperature[_] /= volume[_]
        
    data = [vtk.vtkDoubleArray()]
    data[0].SetName('SolidVolumeFraction')
    data.append(vtk.vtkDoubleArray())
    data[1].SetName('SolidVolumeVelocity')
    data[1].SetNumberOfComponents(3)
    data.append(vtk.vtkDoubleArray())
    data[2].SetName('GranularTemperature')
#    data.append(vtk.vtkDoubleArray())
#    data[3].SetName('SolidPressure')

    for _ in range(ugrid.GetNumberOfCells()):
        data[0].InsertNextValue(volume[_])
        data[1].InsertNextTuple3(*(velocity[_]))
        data[2].InsertNextValue(temperature[_])
#        data[3].InsertNextValue(solid_pressure[_])

    pdata = vtk.vtkDoubleArray()
    pdata.SetName('Time')
    
    for _ in range(ugrid.GetNumberOfPoints()):
        pdata.InsertNextValue(bucket.time)

    for _ in data:
        ugrid.GetCellData().AddArray(_)
    ugrid.GetPointData().AddArray(pdata)
        
    return ugrid

def initializearray(polydata, arrayname, isscalar=True, ispointdata=True):
    """Initialize a data array. Choose pointdata or celldata and scalar or
    vector. Array values are initialized with 0.0 or [0.0, 0.0, 0.0]."""
    if ispointdata:  # add array to pointdata
        numberofpoints = polydata.GetNumberOfPoints()
        array = vtk.vtkDoubleArray()
        array.SetName(arrayname)
        if isscalar:  # array holds scalars
            array.SetNumberOfValues(numberofpoints)
            array.FillComponent(0, 0.0)
            polydata.GetPointData().AddArray(array)
        else:  # array holds vectors
            array.SetNumberOfComponents(3)
            array.SetNumberOfTuples(numberofpoints)
            for j in range(3):
                array.FillComponent(j, 0.0)
            polydata.GetPointData().AddArray(array)
    else:  # add array to celldata
        numberofcells = polydata.GetNumberOfCells()
        array = vtk.vtkDoubleArray()
        array.SetName(arrayname)
        if isscalar:  # array holds scalars
            array.SetNumberOfValues(numberofcells)
            array.FillComponent(0, 0.0)
            polydata.GetCellData().AddArray(array)
        else:  # array holds vectors
            array.SetNumberOfComponents(3)
            array.SetNumberOfTuples(numberofcells)
            for j in range(3):
                array.FillComponent(j, 0.0)
            polydata.GetCellData().AddArray(array)
    return array

    def __init__(self):

        pypes.pypeScript.__init__(self)

        self.Surface = None
        self.DeformedSurface = None
        self.SourcePoints = vtk.vtkPoints()
        self.TargetPoints = vtk.vtkPoints()
        self.DisplacementNorms = vtk.vtkDoubleArray()
        self.Displacements = vtk.vtkDoubleArray()
        self.Displacements.SetNumberOfComponents(3)
        self.SourceSpheres = vtk.vtkPolyData()
        self.TargetSpheres = vtk.vtkPolyData()
        self.SourceSpheres.SetPoints(self.SourcePoints)
        self.TargetSpheres.SetPoints(self.TargetPoints)
        self.SourceSpheres.GetPointData().SetScalars(self.DisplacementNorms)
        self.SourceSpheres.GetPointData().SetVectors(self.Displacements)
        self.vmtkRenderer = None
        self.OwnRenderer = 0
        self.DisplayDeformed = False
        self.SurfaceMapper = None
        self.Opacity = 1.0
        self.SourceSpheresActor = None
        self.TargetSpheresActor = None

        self.SetScriptName("vmtkthinplatesplinedeformation")
        self.SetInputMembers(
            [
                ["Surface", "i", "vtkPolyData", 1, "", "the input surface", "vmtksurfacereader"],
                ["Opacity", "opacity", "float", 1, "(0.0,1.0)", "object opacities in the scene"],
                ["vmtkRenderer", "renderer", "vmtkRenderer", 1, "", "external renderer"],
            ]
        )
        self.SetOutputMembers([["DeformedSurface", "o", "vtkPolyData", 1, "", "", "vmtksurfacewriter"]])

def make_unstructured_grid(mesh, velocity, pressure, time, outfile=None):
    """Given a mesh (in Gmsh format), velocity and pressure fields, and a
    time level, store the data in a vtkUnstructuredGridFormat."""

    pnts = vtk.vtkPoints()
    pnts.Allocate(len(mesh.nodes))

    node2id = {}
    for k, point in mesh.nodes.items():
        node2id[k] = pnts.InsertNextPoint(point)

    ugrid = vtk.vtkUnstructuredGrid()
    ugrid.SetPoints(pnts)

    for element in mesh.elements.values():
        id_list = vtk.vtkIdList()
        for node in element[2]:
            id_list.InsertNextId(node2id[node])

        ugrid.InsertNextCell(TYPE_DICT[element[0]], id_list)

    data = []
    data.append(vtk.vtkDoubleArray())
    data[0].SetNumberOfComponents(3)
    data[0].Allocate(3*pnts.GetNumberOfPoints())
    data[0].SetName('Velocity')

    data.append(vtk.vtkDoubleArray())
    data[1].Allocate(pnts.GetNumberOfPoints())
    data[1].SetName('Pressure')

    data.append(vtk.vtkDoubleArray())
    data[2].Allocate(pnts.GetNumberOfPoints())
    data[2].SetName('Time')

    for k in range(len(mesh.nodes)):
        if hasattr(velocity, '__call__'):
            data[0].InsertNextTuple3(*velocity(ugrid.GetPoints().GetPoint(k)))
        else:
            data[0].InsertNextTuple3(*velocity[k, :])
        if hasattr(pressure, '__call__'):
            data[1].InsertNextValue(pressure(ugrid.GetPoints().GetPoint(k)))
        else:
            data[1].InsertNextValue(pressure[k])
        data[2].InsertNextValue(time)


    for _ in data:
        ugrid.GetPointData().AddArray(_)

    if outfile:
        write_to_file(ugrid, outfile)

    return ugrid

def Gather(c, arr, root):
    vtkArr = vtk.vtkDoubleArray()
    count = len(arr)
    vtkArr.SetNumberOfTuples(count)
    for i in range(count):
        vtkArr.SetValue(i, arr[i])
    vtkResult = vtk.vtkDoubleArray()
    c.Gather(vtkArr, vtkResult, root)
    result = [vtkResult.GetValue(i) for i in range(vtkResult.GetNumberOfTuples())]
    return [ tuple(result[i : i + count]) \
                for i in range(0, vtkResult.GetNumberOfTuples(), count) ]

def CreateParentArteryPatches(parentCenterlines, clipPoints):
    numberOfDaughterPatches = parentCenterlines.GetNumberOfCells()

    patchedCenterlines = vtk.vtkPolyData()
    patchedCenterlinesPoints = vtk.vtkPoints()
    patchedCenterlinesCellArray = vtk.vtkCellArray()
    patchedRadiusArray = vtk.vtkDoubleArray()

    clipIds, numberOfPatchedCenterlinesPoints = ExtractPatchesIds(parentCenterlines, clipPoints)
    print "Clipping Point Ids ", clipIds

    radiusArray = vtk.vtkDoubleArray()
    radiusArray.SetNumberOfComponents(1)
    radiusArray.SetName(radiusArrayName)
    radiusArray.SetNumberOfTuples(numberOfPatchedCenterlinesPoints)
    radiusArray.FillComponent(0, 0.0)

    numberOfCommonPatch = clipIds[0] + 1
    patchedCenterlinesCellArray.InsertNextCell(numberOfCommonPatch)

    count = 0
    for i in range(0, numberOfCommonPatch):
        patchedCenterlinesPoints.InsertNextPoint(parentCenterlines.GetPoint(i))
        patchedCenterlinesCellArray.InsertCellPoint(i)
        radiusArray.SetTuple1(i, parentCenterlines.GetPointData().GetArray(radiusArrayName).GetTuple1(i))
        count += 1

    for j in range(numberOfDaughterPatches):
        cell = vtk.vtkGenericCell()
        parentCenterlines.GetCell(j, cell)

        numberOfCellPoints = cell.GetNumberOfPoints()
        startId = clipIds[j + 1]
        patchNumberOfPoints = numberOfCellPoints - startId
        patchedCenterlinesCellArray.InsertNextCell(patchNumberOfPoints)

        for i in range(startId, cell.GetNumberOfPoints()):
            point = cell.GetPoints().GetPoint(i)
            patchedCenterlinesPoints.InsertNextPoint(point)
            patchedCenterlinesCellArray.InsertCellPoint(count)
            radiusArray.SetTuple1(
                count, parentCenterlines.GetPointData().GetArray(radiusArrayName).GetTuple1(cell.GetPointId(i))
            )
            count += 1

    patchedCenterlines.SetPoints(patchedCenterlinesPoints)
    patchedCenterlines.SetLines(patchedCenterlinesCellArray)
    patchedCenterlines.GetPointData().AddArray(radiusArray)

    return patchedCenterlines

def vtkbasis(mesh, etob, fname, coeffs):
    ''' Find the directions from a (non-uniform) plane wave basis and output a VTK-compatible file
    
        It's possible that this needs to be updated to work with recent changes to ElementToBasis
    '''
    try:
        import vtk
        
        points = vtk.vtkPoints()
        vectors = vtk.vtkDoubleArray()
        vectors.SetNumberOfComponents(3)
        scalars = vtk.vtkDoubleArray()
        
        nc = 0
        for e in range(mesh.nelements):
            c = paa.origin(mesh, e)
            bs = etob[e]
            cc = np.zeros(3)
            cc[:len(c)] = c
            nondir = 0
            ndir = 0
            for b in bs:
                if hasattr(b, "directions"):
                    for d in b.directions.transpose():
                        dd = np.zeros(3)
                        dd[:len(d)] = d
                        if coeffs is not None: dd *= abs(coeffs[nc])
                        points.InsertNextPoint(*cc)
                        vectors.InsertNextTuple3(*dd)
                        ndir+=1
                        nc+=1
                else:
                    nondir += np.sqrt(np.sum(coeffs[nc:nc+b.n]**2))
                    nc += b.n
            for _ in range(ndir): scalars.InsertNextValue(nondir)
                    
        g = vtk.vtkUnstructuredGrid()
        g.SetPoints(points)
        gpd = g.GetPointData()
        gpd.SetVectors(vectors)
        gpd.SetScalars(scalars)
        writer = vtk.vtkXMLUnstructuredGridWriter()
        writer.SetFileName(fname)
        writer.SetInput(g)
        writer.Write()
    except ImportError as e:
        print "Unable to write basis to file: ",e
                
                

 def ReadTetrIniFile(self,filename,arrayname):
     self.PrintLog('Reading '+filename+'.')
     if (self.Compressed == 1):
         f=gzip.open(filename, 'r')
     else:
         f=open(filename, 'r')
     line = f.readline()
     while (line!='$vel_old1'):
         line = f.readline()
     line = f.readline()
     line = f.readline()
     numberOfTuples = int(line)
     outputArray = vtk.vtkDoubleArray()
     outputArray.SetName(arrayname)
     outputArray.SetNumberOfComponents(3)
     outputArray.SetNumberOfTuples(numberOfTuples)
     for i in range(numberOfTuples):
         line = f.readline()
         splitline = line.split(' ')
         value0 = float(splitline[0])
         value1 = float(splitline[1])
         value2 = float(splitline[2])
         outputArray.SetComponent(i,0,value0)
         outputArray.SetComponent(i,1,value1)
         outputArray.SetComponent(i,2,value2)
     self.Mesh.GetPointData().AddArray(outputArray)

 def ReadTetrVelFile(self,filename,arrayname):
     self.PrintLog('Reading '+filename+'.')
     if self.UnNormalize == 1:
         self.ReadTetrInDimensionalParameters()
     if (self.Compressed == 1):
         f=gzip.open(filename, 'r')
     else:
         f=open(filename, 'r')
     lines = f.readlines()
     line = lines[0]
     lineoffset = 1
     splitline = line.split(' ')
     numberOfTuples = int(splitline[0])
     iteration = float(splitline[1])
     outputArray = vtk.vtkDoubleArray()
     outputArray.SetName(arrayname)
     outputArray.SetNumberOfComponents(3)
     outputArray.SetNumberOfTuples(numberOfTuples)
     velocityUnNormalizationFactor = self.ReD / (2.0 * self.Radius) * self.Viscosity / self.Density 
     for i in range(numberOfTuples):
         line = lines[i+lineoffset]
         splitline = line.split(' ')
         value0 = float(splitline[0])
         value1 = float(splitline[1])
         value2 = float(splitline[2])
         if self.UnNormalize ==1:
             value0 *= velocityUnNormalizationFactor
             value1 *= velocityUnNormalizationFactor
             value2 *= velocityUnNormalizationFactor
         outputArray.SetComponent(i,0,value0)
         outputArray.SetComponent(i,1,value1)
         outputArray.SetComponent(i,2,value2)
     self.Mesh.GetPointData().AddArray(outputArray)

 def addVectorToNode(self,fieldValues,name):
     datas = vtk.vtkDoubleArray()
     datas.SetNumberOfComponents(3)
     datas.SetName(name)
     for value in fieldValues:
         datas.InsertNextTuple3(value[0],value[1],value[2])
     self._hexs.GetPointData().AddArray(datas)

 def convertArray2vtkImage(self, nparray, t_ImagedataVTK, npImagesandMask): 
     """ Takes a numpy.ndarray and converts it to a vtkimageData. require npImagesandMask to pass on image info """
     # Create vtk object
     size_array = npImagesandMask['dims'][0]*npImagesandMask['dims'][1]*npImagesandMask['dims'][2]
     flatim = nparray.transpose(2,1,0)
     flatim = flatim.flatten()
     
     # create vtk image
     vtk_image = vtk.vtkImageData()
     vtk_image.DeepCopy(t_ImagedataVTK)
     vtk_image.SetNumberOfScalarComponents(1)
     vtk_image.SetScalarTypeToDouble()
     vtk_image.AllocateScalars()
     
     # Get scalars from numpy
     image_array = vtk.vtkDoubleArray() 
     image_array.SetNumberOfValues(size_array)
     image_array.SetNumberOfComponents(1) 
     
     # not too efficient convertion of np.array to vtk. Far from ideal
     for k in range(size_array):
         image_array.InsertTuple1(k,flatim[k])
         
     vtk_image.GetPointData().SetScalars(image_array) 
     vtk_image.Update()
       
     return vtk_image   

def buildATPMesh(polydata, filename):
    
    centroidFilter = vtk.vtkCellCenters()
    centroidFilter.VertexCellsOn()
    centroidFilter.SetInputData(polydata)
    
    newPolydata = vtk.vtkPolyData()
    newPolydata = centroidFilter.GetOutput()
    centroidFilter.Update()
    
    ATPValues = vtk.vtkDoubleArray()
    ATPValues.SetName("initialATP")
    
    _, _, yMin, yMax, _, _ = polydata.GetBounds()
    yRange = yMax - yMin
    
    for pointId in range(0, newPolydata.GetNumberOfPoints()):
        _, y, _ = newPolydata.GetPoint(pointId)
        ATPValue = y / (yRange * 1.0)
        ATPValues.InsertNextValue(ATPValue)
        
    newPolydata.GetCellData().SetScalars(ATPValues)
    
    polyDataWriter = vtk.vtkXMLPolyDataWriter()
    polyDataWriter.SetFileName(filename)
    polyDataWriter.SetInputData(newPolydata)    
    polyDataWriter.Write()

def createTrail(ts):
    Points = vtk.vtkPoints()
    id_array = vtk.vtkIntArray()
    #id_array.SetNumberofComponents(1)
    id_array.SetName("haloid")
    
    phi_array = vtk.vtkDoubleArray()
    phi_array.SetName("phi")
    for i in range(0,ts+1):
        px,py,pz,phid,pphi = readParticle(i)
        for j in range(0,len(px)):
            Points.InsertNextPoint(px[j],py[j],pz[j])
            id_array.InsertNextTuple1(phid[j])
            phi_array.InsertNextTuple1(pphi[j])
    
    polydata = vtk.vtkPolyData()
    polydata.SetPoints(Points)
    polydata.GetPointData().AddArray(id_array)
    polydata.GetPointData().AddArray(phi_array)
    
    if vtk.VTK_MAJOR_VERSION <= 5:
        polydata.Update()
        
    
    outputFile = "/home/subhashis/VisData/merger_trees/particleTrail/time" + str(ts) + ".vtp" 
        
    writer = vtk.vtkXMLPolyDataWriter();
    writer.SetFileName(outputFile);
    if vtk.VTK_MAJOR_VERSION <= 5:
        writer.SetInput(polydata)
    else:
        writer.SetInputData(polydata)
    writer.Write()
    print "Done generating output for time %d" %ts

 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)

def make_sphereActor (x, r, rgb, opacity):
    points = vtk.vtkPoints()
    points.InsertNextPoint(x[0], x[1], x[2])

    diameter = vtk.vtkDoubleArray()
    diameter.SetNumberOfComponents(1)
    diameter.InsertNextTuple1(2.0*r)

    pData = vtk.vtkPolyData()
    pData.SetPoints(points)
    pData.GetPointData().SetScalars(diameter)

    pSource = vtk.vtkSphereSource()
    pSource.SetPhiResolution(16)
    pSource.SetThetaResolution(16)

    pGlyph = vtk.vtkGlyph3D()
    pGlyph.SetSource(pSource.GetOutput())
    pGlyph.SetInput(pData)
    pGlyph.ScalingOn()
    pGlyph.SetScaleModeToScaleByScalar()

    pMapper = vtk.vtkPolyDataMapper()
    pMapper.ScalarVisibilityOff()
    pMapper.SetInput(pGlyph.GetOutput())

    pActor = vtk.vtkActor()
    pActor.SetMapper(pMapper)
    pActor.GetProperty().SetColor(rgb[0], rgb[1], rgb[2])
    pActor.GetProperty().SetOpacity(opacity)

    return pActor

    def GridSurface(self):
        #----------------------------------------------------------------------
        # ELECTRODE CONTACT DELAUNAY TRIANGULATION ----------------------------
        deln = vtk.vtkDelaunay3D()
        deln.SetInput(self.electrodePolyData)
        deln.SetTolerance(0.01)
        tmapper = vtk.vtkTextureMapToSphere()
        tmapper.SetInputConnection(deln.GetOutputPort())
        #tmapper.PreventSeamOn()
        mapper = vtk.vtkDataSetMapper()
        mapper.SetInputConnection(tmapper.GetOutputPort())

        # TEST TEXTURE PART
        atext = vtk.vtkOpenGLTexture()

        self.contactScalarData = vtk.vtkDoubleArray()
        self.contactScalarData.SetNumberOfComponents(1)
        for tupleID in range(self.electrodePolyData.GetNumberOfPoints()):
            self.contactScalarData.InsertNextValue(random.uniform(0, 1))

        self.UpdateGridSurface()
        atext.SetInput(self.electrodePolyData)
        self.triangulation = vtk.vtkOpenGLActor()
        self.triangulation.SetMapper(mapper)
        self.triangulation.SetTexture(atext)
        self.triangulation.GetProperty().SetOpacity(1)

def array_to_vtk(array_in, dtype=None):
    """Get vtkFloatArray/vtkDoubleArray from the input numpy array."""
    if dtype == None:
        dtype = _numpy.dtype(array_in.dtype)
    else:
        dtype = _numpy.dtype(dtype)
    if dtype == _numpy.float32:
        float_array = _vtk.vtkFloatArray()
    elif dtype == _numpy.float64:
        float_array = _vtk.vtkDoubleArray()
    elif dtype == _numpy.uint8:
        float_array = _vtk.vtkUnsignedCharArray()
    elif dtype == _numpy.int8:
        float_array = _vtk.vtkCharArray()
    else:
        raise ValueError("Wrong format of input array, must be float32 or float64")
    if len(array_in.shape) != 1 and len(array_in.shape) != 2:
        raise ValueError("Wrong shape: array must be 1D or 2D.")
    #float_array.SetNumberOfComponents(_numpy.product(array_in.shape))
    # if len(array_in.shape) == 2:
    #     float_array.SetNumberOfComponents(array_in.shape[1])
    # elif len(array_in.shape) == 1:
    #     float_array.SetNumberOfComponents(1)
    float_array.SetNumberOfComponents(1)
    array_contiguous = _numpy.ascontiguousarray(array_in, dtype)
    float_array.SetVoidArray(array_contiguous, _numpy.product(array_in.shape), 1)
    float_array._contiguous_array = array_contiguous  # Hack to keep the array of being garbage collected
    # if len(array_in.shape) == 2:
    #     print "set tuple to {0}".format(array_in.shape[1])
    #     #float_array.SetNumberOfTuples(array_in.shape[1])
    #     float_array.Resize(array_in.shape[1])
    #     float_array.Squeeze()
    return float_array