def OnExportSurface(self, pubsub_evt):
        filename, filetype = pubsub_evt.data
        if (filetype == const.FILETYPE_STL) or\
           (filetype == const.FILETYPE_VTP) or\
           (filetype == const.FILETYPE_PLY) or\
           (filetype == const.FILETYPE_STL_ASCII):

            # First we identify all surfaces that are selected
            # (if any)
            proj = prj.Project()
            polydata_list = []

            for index in proj.surface_dict:
                surface = proj.surface_dict[index]
                if surface.is_shown:
                    polydata_list.append(surface.polydata)

            if len(polydata_list) == 0:
                utl.debug("oops - no polydata")
                return
            elif len(polydata_list) == 1:
                polydata = polydata_list[0]
            else:
                polydata = pu.Merge(polydata_list)

            # Having a polydata that represents all surfaces
            # selected, we write it, according to filetype
            if filetype == const.FILETYPE_STL:
                writer = vtk.vtkSTLWriter()
                writer.SetFileTypeToBinary()
            elif filetype == const.FILETYPE_STL_ASCII:
                writer = vtk.vtkSTLWriter()
                writer.SetFileTypeToASCII()
            elif filetype == const.FILETYPE_VTP:
                writer = vtk.vtkXMLPolyDataWriter()
            #elif filetype == const.FILETYPE_IV:
            #    writer = vtk.vtkIVWriter()
            elif filetype == const.FILETYPE_PLY:
                writer = vtk.vtkPLYWriter()
                writer.SetFileTypeToASCII()
                writer.SetColorModeToOff()
                #writer.SetDataByteOrderToLittleEndian()
                #writer.SetColorModeToUniformCellColor()
                #writer.SetColor(255, 0, 0)

            if filetype in (const.FILETYPE_STL, const.FILETYPE_PLY):
                # Invert normals
                normals = vtk.vtkPolyDataNormals()
                normals.SetInputData(polydata)
                normals.SetFeatureAngle(80)
                normals.AutoOrientNormalsOn()
                #  normals.GetOutput().ReleaseDataFlagOn()
                normals.UpdateInformation()
                normals.Update()
                polydata = normals.GetOutput()

            filename = filename.encode(wx.GetDefaultPyEncoding())
            writer.SetFileName(filename)
            writer.SetInputData(polydata)
            writer.Write()

def CreateCrossSections(files, outputAreas, outputGeometryFile, outputOutlineFile):
  areaFile = open(outputAreas, 'w')

  planeAppender = vtk.vtkAppendPolyData()
  outlineAppender = vtk.vtkAppendPolyData()
  for idx in range(len(files)):
    print 'Processing contour %d' % idx
    file = files[idx]
    (plane, outline) = CreatePlanarCrossSectionPolyDataFromFile(file)
    areaFile.write(str(ComputePolyDataArea(plane)))
    areaFile.write('\n')
    planeAppender.AddInputData(plane)
    outlineAppender.AddInputData(outline)
    
  planeWriter = vtk.vtkXMLPolyDataWriter()
  planeWriter.SetFileName(outputGeometryFile)
  planeWriter.SetInputConnection(planeAppender.GetOutputPort())
  planeWriter.Update()
  
  outlineWriter = vtk.vtkXMLPolyDataWriter()
  outlineWriter.SetFileName(outputOutlineFile)
  outlineWriter.SetInputConnection(outlineAppender.GetOutputPort())
  outlineWriter.Update()

  areaFile.close()

    def __init__(self, module_manager):

        # call parent constructor
        ModuleBase.__init__(self, module_manager)

        self._writer = vtk.vtkXMLPolyDataWriter()

        module_utils.setup_vtk_object_progress(
            self, self._writer,
            'Writing VTK PolyData')

        

        self._writer.SetDataModeToBinary()

        # ctor for this specific mixin
        FilenameViewModuleMixin.__init__(
            self,
            'Select a filename',
            'VTK PolyData (*.vtp)|*.vtp|All files (*)|*',
            {'vtkXMLPolyDataWriter': self._writer},
            fileOpen=False)
            

        # set up some defaults
        self._config.filename = ''

        self.sync_module_logic_with_config()

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 save_polydata(polydata, file_name, binary=False, color_array_name=None):
    # get file extension (type)
    file_extension = file_name.split(".")[-1].lower()

    # todo better generic load
    # todo test all
    if file_extension == "vtk":
        writer = vtk.vtkPolyDataWriter()
    elif file_extension == "vtp":
        writer = vtk.vtkPolyDataWriter()
    elif file_extension == "fib":
        writer = vtk.vtkPolyDataWriter()
    elif file_extension == "ply":
        writer = vtk.vtkPLYWriter()
    elif file_extension == "stl":
        writer = vtk.vtkSTLWriter()
    elif file_extension == "xml":
        writer = vtk.vtkXMLPolyDataWriter()
    elif file_extension == "obj":
        raise "mni obj or Wavefront obj ?"
    #    writer = set_input(vtk.vtkMNIObjectWriter(), polydata)

    writer.SetFileName(file_name)
    writer = set_input(writer, polydata)
    if color_array_name is not None:
        writer.SetArrayName(color_array_name);
    
    if binary :
        writer.SetFileTypeToBinary()
    writer.Update()
    writer.Write()

def save_lesion(output_file, input_file, field, limits):
    reader = v.vtkXMLUnstructuredGridReader()
    reader.SetFileName(input_file)
    reader.Update()

    threshold = v.vtkThreshold()
    threshold.SetInput(reader.GetOutput())

    if limits[0] is None:
        threshold.ThresholdByLower(limits[1])
    elif limits[1] is None:
        threshold.ThresholdByUpper(limits[0])
    else:
        threshold.ThresholdBetween(*limits)

    threshold.SetInputArrayToProcess(0, 0, 0, v.vtkDataObject.FIELD_ASSOCIATION_CELLS, field)
    threshold.Update()

    extract_surface = v.vtkDataSetSurfaceFilter()
    extract_surface.SetInput(threshold.GetOutput())
    extract_surface.Update()

    writer = v.vtkXMLPolyDataWriter()
    writer.SetFileName(output_file)
    writer.SetInput(extract_surface.GetOutput())
    writer.Write()

def pointToCellData():
    
    print "Using angle: " + str(angle)
    atpFiles = glob.glob(str(angle) + '/*.vtp')
    if not atpFiles:
        exit("No atp files found")
        
    if not os.path.exists(str(angle) + '_new'):
        os.makedirs(str(angle) + '_new')
    
    for inputFile in atpFiles: 

        print 'Reading', inputFile
        atpReader = vtk.vtkXMLPolyDataReader()
        atpReader.SetFileName(inputFile)
        atpReader.Update()
    
        atpDataset = atpReader.GetOutput()
        
        pointToCell = vtk.vtkPointDataToCellData()
        pointToCell.SetInputData(atpDataset)
        pointToCell.PassPointDataOn()
        pointToCell.Update()
        
        newArray = pointToCell.GetOutput()
        
        # Remove the point data arrays, they exist as cell data
        newArray.GetPointData().RemoveArray('ATP')
        newArray.GetPointData().RemoveArray('tau')
        newArray.GetCellData().RemoveArray('p')
    
        atpWriter = vtk.vtkXMLPolyDataWriter()
        atpWriter.SetInputData(newArray)
        atpWriter.SetFileName(str(angle) + '_new/' + inputFile[3:])
        atpWriter.Update()

    def exportModel(self, model, outDir):
        
        modelDisplay = model.GetDisplayNode()
        if not modelDisplay or not modelDisplay.GetVisibility():
            return None
        
        print('  exporting %s' % model.GetName())
        
        writer = vtk.vtkXMLPolyDataWriter()
        if vtk.VTK_MAJOR_VERSION <= 5:
            writer.SetInput(model.GetPolyData())
        else:
            writer.SetInputData(model.GetPolyData())
        fileName = os.path.join(outDir, model.GetName() + '.vtp')
        writer.SetFileName(fileName)
        writer.Write()

        data = {}
        data['name'] = model.GetName()
        data['color'] = modelDisplay.GetColor()
        data['opacity'] = modelDisplay.GetOpacity()
        data['filename'] = fileName
        
        representation = modelDisplay.GetRepresentation()
        if representation == slicer.vtkMRMLDisplayNode.WireframeRepresentation:
            data['geometry_mode'] = 'wireframe'
        elif representation == slicer.vtkMRMLDisplayNode.SurfaceRepresentation \
          and modelDisplay.GetEdgeVisibility():
            data['geometry_mode'] = 'surface_with_edges'
        elif representation == slicer.vtkMRMLDisplayNode.PointsRepresentation:
            data['geometry_mode'] = 'points'

        return data

def convert_mris(input_name, output_name):
  """
  Converts a FreeSurfer surface to VTK file format
  """
  # convert surface to VTK format
  if output_name.endswith('.vtp'): temp_name = output_name[0:-1] + 'k'
  else:                            temp_name = output_name
  if not temp_name.endswith('.vtk'):
    raise RuntimeError('Output file name extension must be either .vtk or .vtp')
  check_call(['mris_convert', input_name, temp_name])
  # get surface RAS translation
  out = check_output(["mris_info", input_name], stderr=STDOUT)
  m = re.search("c_\(ras\)\s:\s\((-?\d+\.\d+),\s(-?\d+\.\d+),\s(-?\d+\.\d+)\)", out)
  if m is None: raise RuntimeError('Could not find c_(ras) coordinates in mris_info output!')
  tx = float(m.group(1))
  ty = float(m.group(2))
  tz = float(m.group(3))
  # transform vertex positions to scanner RAS of orig.mgz
  reader = vtkPolyDataReader()
  reader.SetFileName(temp_name)
  reader.Update()
  surface = reader.GetOutput()
  points = surface.GetPoints()
  for i in range(points.GetNumberOfPoints()):
    x, y, z = points.GetPoint(i)
    points.SetPoint(i, x + tx, y + ty, z + tz)
  surface.SetPoints(points)
  if output_name.endswith('.vtp'): writer = vtkXMLPolyDataWriter()
  else:                            writer = vtkPolyDataWriter()
  writer.SetFileName(output_name)
  writer.SetInput(surface)
  writer.Write()
  if temp_name != output_name:
    remove(temp_name)

	def writePolyData(self, polyData, filename, callback = None):
		"""
		Writes the given vtkPolyData instance to disk
					 as .vtp-file with the given filename
		"""
		writer = vtk.vtkXMLPolyDataWriter()
		writer.SetFileName(filename)
		writer.SetInput(polyData)
		def f(obj, evt):
			if obj and callback:
				callback(obj.GetProgress())
			if scripting.mainWindow:
				scripting.mainWindow.updateProgressBar(obj, evt, obj.GetProgress(),"Writing surface %s"%os.path.basename(filename), 0)
		writer.AddObserver("ProgressEvent", lib.messenger.send)
		lib.messenger.connect(writer, "ProgressEvent", f)
		Logging.info("Writing polydata to file %s"%filename,kw="pipeline")
		try:
			ret = writer.Write()
			
			if ret == 0:
				Logging.error("Failed to write polygonal data",
				"Failed to write vtkPolyData object to file %s" % filename)
				return
		except Exception, ex:
			Logging.error("Failed to write poly data",
			"Failed to write vtkPolyData object to file %s" % filename, ex)
			return

def create_mesh(M,P,T):
	n_naca_pts = 50
	write_test_file(M,P,T,-5.,5.,nsections=5)
	wing=create_wing('current_wing','output')
	n_sections=len(wing)
	n_section_pts = 2*n_naca_pts-1

	# build mesh
	vtk_model = vtk.vtkStructuredGrid()
	vtk_model.SetDimensions(n_section_pts,n_sections,1)
	# build points
	vtk_points = vtk.vtkPoints()

	for j in xrange(n_sections):
	    upper_pts = numpy.array([wing[j][1],wing[j][3]]).T
	    lower_pts = numpy.array([wing[j][2],wing[j][4]]).T
	    section_pts = numpy.concatenate((lower_pts[::-1],upper_pts[1:]))
	    for i in xrange(n_section_pts):
		vtk_points.InsertNextPoint(section_pts[i,0],wing[j][0],section_pts[i,1])
	# set points
	vtk_model.SetPoints(vtk_points)

	# convert to poly data    
	pdata_filter = vtk.vtkGeometryFilter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    pdata_filter.SetInput(vtk_model)
	else:
	    pdata_filter.SetInputData(vtk_model)
	pdata_filter.Update()
	poly_data = pdata_filter.GetOutput()

	# compute normals
	norms = vtk.vtkPolyDataNormals()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    norms.SetInput(poly_data)
	else:
	    norms.SetInputData(poly_data)
	norms.ComputePointNormalsOff()
	norms.ComputeCellNormalsOn()
	norms.ConsistencyOn()
	norms.Update()

	# clean poly data
	clean_poly = vtk.vtkCleanPolyData()
	clean_poly.ToleranceIsAbsoluteOn()
	clean_poly.SetAbsoluteTolerance(1.e-6)
	if vtk.VTK_MAJOR_VERSION <= 5:
	    clean_poly.SetInput(norms.GetOutput())
	else:
	    clean_poly.SetInputData(norms.GetOutput())
	clean_poly.Update()

	# write output mesh
	writer = vtk.vtkXMLPolyDataWriter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    writer.SetInput(clean_poly.GetOutput())
	else:
	    writer.SetInputData(clean_poly.GetOutput())
	writer.SetFileName('output.vtp')
	writer.Write()

def main(argv=None):
    dimension = 2
    numprocs = 2
    level = 0  # startlevel

    no_multigridlevels = 0

    # check how many processors generated aggregation output
    # while check_files_for_next_level(0,numprocs, "aggs_level%LEVEL_proc%PROC.out") == True:
    #  numprocs = numprocs + 1
    # numprocs = numprocs - 1
    # print "Aggregtaion information for " + str(numprocs) + " processors found"

    # process all multigrid levels
    while check_files_for_next_level(level, numprocs, "aggs_level%LEVEL_proc%PROC.out"):
        global_nodecoords = []

        print "Level " + str(level)

        if level == 0:  # read in coordinates (finest level
            global_nodecoords, dimension = read_finelevel_nodecoords_from_file("example.txt")
        else:
            global_nodecoords, dimension = read_nodecoords_from_file("nodes" + str(level) + ".txt")

        # read aggregates
        aggid2nodes, aggid2procs = readin_aggregates("aggs_level%LEVEL_proc%PROC.out", numprocs, level)

        # vtk polygon for output
        aggpolygons = vtk.vtkAppendPolyData()

        # collect all aggregates
        for aggid, agg_nodes in aggid2nodes.iteritems():
            # build an aggregate
            if dimension == 2:
                prepareDelaunayData(
                    dimension, agg_nodes, global_nodecoords, aggpolygons, aggid, aggid2nodes, aggid2procs
                )
            else:
                prepareDelaunayData3d(
                    dimension, agg_nodes, global_nodecoords, aggpolygons, aggid, aggid2nodes, aggid2procs
                )

                # aggpolygons.GetOutput().GetPointData().SetVectors(vtkDisplacementVector)
                # aggpolygons.Update()

        writer = vtk.vtkXMLPolyDataWriter()
        fname = "aggs" + str(level) + ".vtp"
        writer.SetFileName(fname)
        writer.SetInput(aggpolygons.GetOutput())
        writer.Write()

        write_nodes_file("nodes" + str(level + 1) + ".txt", aggid2nodes, global_nodecoords, dimension)

        # increment number of multigrid levels that have been found in the files
        if no_multigridlevels < level:
            no_multigridlevels = level

        print "VTK Export for level " + str(level) + " finished...\r\n"

        level = level + 1

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkXMLPolyDataWriter(), 'Writing vtkXMLPolyData.',
         ('vtkXMLPolyData',), (),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

    def addMeshFromFile(self, name, filename):
        """
        Add a mesh shape from a file.
        Accepted format : .stl or mesh encoded in VTK .vtp format
        """
        if name not in self._ref:

            if os.path.splitext(filename)[-1][1:] == 'stl':
                reader = vtk.vtkSTLReader()
                reader.SetFileName(filename)
                reader.Update()

                with tmpfile() as tmpf:
                    writer=vtk.vtkXMLPolyDataWriter()
                    writer.SetInputData(reader.GetOutput())
                    writer.SetFileName(tmpf[1])
                    writer.Write()

                    shape_data = str_of_file(tmpf[1])

            else:
                assert os.path.splitext(filename)[-1][1:] == 'vtp'
                shape_data = str_of_file(filename)


            self.addMeshShapeFromString(name, shape_data)

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 save_one_file(i):
    pts = vtk.vtkPoints()
    df1 = df[df['T']==unique_times[i]]
    x = df1['COORX']
    y = df1['COORY']
    z = df1['COORZ']
    tp = df1['ITYPE']
    pts.SetNumberOfPoints(df1.shape[0])

    if df1.shape[0]==0:
        empty_file[i]=1
    else:
        for j in range(df1.shape[0]):         
            pts.SetPoint(j, (x.iloc[j],y.iloc[j],z.iloc[j]))

        types = vtk.vtkShortArray()
        types.SetNumberOfComponents(1)
        types.SetNumberOfTuples(pts.GetNumberOfPoints())
        types.SetName('Type')
        
        for j in range(df1.shape[0]):         
            types.SetTuple1(j, tp.iloc[j])

        pd = vtk.vtkPolyData()
        pd.SetPoints(pts)
        pd.GetPointData().AddArray(types)

        wr= vtk.vtkXMLPolyDataWriter()
        wr.SetInputData(pd)
        wr.SetDataModeToBinary()
        filename = "pts_{:010d}.vtp".format(i)
        wr.SetFileName(filename)
        wr.Write()

def write_polydata(polydata, filename):
    """Write polydata as vtkPolyData format, according to extension."""

    if VERBOSE:
        print "Writing ", filename, "..."

    basename, extension = os.path.splitext(filename)

    if   (extension == '.vtk'):
        writer = vtk.vtkPolyDataWriter()
        writer.SetFileTypeToBinary()
    elif (extension == '.vtp'):
        writer = vtk.vtkXMLPolyDataWriter()
        writer.SetDataModeToBinary()
    else:
        print 'Cannot recognize model file format'
        return None

    writer.SetFileName(filename)
    if (vtk.vtkVersion().GetVTKMajorVersion() >= 6.0):
        writer.SetInputData(polydata)
    else:
        writer.SetInput(polydata)
    writer.Update()

    del writer

    if VERBOSE:
        print "Done writing ", filename
        print "Number of lines found:", outpd.GetNumberOfLines()

def rescaleATPToRange():
    # This is where the data is for testing purposes.
    print "Current working directory:", os.getcwd()

    print 'Reading', inputFile
    atpReader = vtk.vtkXMLPolyDataReader()
    atpReader.SetFileName(inputFile)
    atpReader.Update()

    atpDataset = atpReader.GetOutput()
    atp = atpDataset.GetPointData().GetArray('ATP')

    inMin, inMax = atp.GetRange()

    atpR = vtk.vtkFloatArray()
    atpR.SetName('initialATP')    
    
    for i in range(atp.GetNumberOfTuples()):
        val = atp.GetTuple(i)[0]
        atpR.InsertNextValue(rescale(val, inMin, inMax))

    atpDataset.GetCellData().AddArray(atpR)
    
    atpDataset.GetPointData().RemoveArray('ATP')
    atpDataset.GetPointData().RemoveArray('tau_w')
    atpDataset.GetCellData().RemoveArray('initialJPLC')

    print 'Writing', outputFile
    atpWriter = vtk.vtkXMLPolyDataWriter()
    atpWriter.SetInput(atpDataset)
    atpWriter.SetFileName(outputFile)
    atpWriter.Update()

    print 'Reading', inputSurfaceFile
    surfaceReader = vtk.vtkXMLPolyDataReader()
    surfaceReader.SetFileName(inputSurfaceFile)
    surfaceReader.Update()

    surface = surfaceReader.GetOutput()
    surface.GetCellData().AddArray(atpR)

    print 'Writing', outputSurfaceFile
    surfaceWriter = vtk.vtkXMLPolyDataWriter()
    surfaceWriter.SetInput(surface)
    surfaceWriter.SetFileName(outputSurfaceFile)
    surfaceWriter.Update()

def geom_viz(geomData, filename='sample.vtp'):
    """
    DESCRIPTION
    -----------
    geom_viz(geomData, filename='geom.vtp')
    Generate a simple point centered vtp file to visualize grain ID.
    PARAMETERS
    ----------
    geomData: numpy.array
    Grain ID populated in a numpy array representing the extruded
    microstructure.
    filename: str
    Output VTP file name.
    RETURNS
    -------
    NOTES
    ----
    """
    polydata = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    vertices = vtk.vtkCellArray()

    gids = vtk.vtkFloatArray()
    gids.SetNumberOfComponents(1)
    gids.SetName("GrainID")

    # iterating through CPFFT data
    cnt = 0
    x, y, z = geomData.shape
    for i in range(x):
        for j in range(y):
            for k in range(z):
                # set position
                points.InsertNextPoint(i, j, k)
                vertex = vtk.vtkVertex()
                vertex.GetPointIds().SetId(0, cnt)
                cnt += 1
                vertices.InsertNextCell(vertex)
                # set Grain ID
                gids.InsertNextTuple1(geomData[i, j, k])
    # finish the vtp object
    polydata.SetPoints(points)
    polydata.SetVerts(vertices)
    polydata.GetPointData().SetScalars(gids)
    polydata.Modified()

    if vtk.VTK_MAJOR_VERSION <= 5:
        polydata.Update()

    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetFileName(filename)

    if vtk.VTK_MAJOR_VERSION <= 5:
        writer.SetInput(polydata)
    else:
        writer.SetInputData(polydata)

    writer.Write()

def copyScalarsToVtk(nifti_img, vtk_file, output_file, scalar_name='Scalar'):
    img = nib.load(nifti_img)
    affine = np.matrix(img.get_affine())
    print affine
    img_data = img.get_data()

    ph_data = np.zeros(img_data.shape)
    dims = ph_data.shape
    for i,v in enumerate(np.linspace(0,10,num=dims[0])):
        ph_data[i] += v

    coords = []
    vals = []

    aff_inv = affine.I

    from vtkFileIO import vtkToStreamlines
    import vtk

    streams, vtkdata = vtkToStreamlines(vtk_file)
    points = np.concatenate(streams)
    print 'world'
    max = np.amax(points, axis=0)
    print max
    min = np.amin(points, axis=0)
    print min


    points = np.hstack((points,np.ones((points.shape[0],1))))
    #coords = np.hstack((coords,np.ones((coords.shape[0],1))))

    #coords_world = coords * affine
    points_ijk = aff_inv * points.T
    points_ijk = points_ijk.T[:,:3]

    print 'ijk'
    max = np.amax(points_ijk, axis=0)
    print max
    min = np.amin(points_ijk, axis=0)
    print min

    import time
    s = time.time()
    scalars = trilinear_interpolator_speedup(img_data, points_ijk)

    vtkScalars = vtk.vtkFloatArray()
    vtkScalars.SetName(scalar_name)
    vtkScalars.SetNumberOfComponents(1)
    #vtkScalars.SetNumberOfTuples(len(scalars))
    for i,v in enumerate(scalars):
        vtkScalars.InsertNextTuple1(v)
    vtkdata.GetPointData().SetScalars(vtkScalars)

    writer = vtk.vtkXMLPolyDataWriter()
    #writer = vtk.vtkPolyDataWriter()
    writer.SetInput(vtkdata)
    writer.SetFileName(output_file)
    writer.Write()