def GetRawDICOMData(filenames,fileID):
  print filenames,fileID
  vtkRealDcmReader = vtk.vtkDICOMImageReader()
  vtkRealDcmReader.SetFileName("%s/%s"%(rootdir,filenames[0]) )
  vtkRealDcmReader.Update()
  vtkRealData = vtk.vtkImageCast()
  vtkRealData.SetOutputScalarTypeToFloat()
  vtkRealData.SetInput( vtkRealDcmReader.GetOutput() )
  vtkRealData.Update( )
  real_image = vtkRealData.GetOutput().GetPointData() 
  real_array = vtkNumPy.vtk_to_numpy(real_image.GetArray(0)) 

  vtkImagDcmReader = vtk.vtkDICOMImageReader()
  vtkImagDcmReader.SetFileName("%s/%s"%(rootdir,filenames[1]) )
  vtkImagDcmReader.Update()
  vtkImagData = vtk.vtkImageCast()
  vtkImagData.SetOutputScalarTypeToFloat()
  vtkImagData.SetInput( vtkImagDcmReader.GetOutput() )
  vtkImagData.Update( )
  imag_image = vtkImagData.GetOutput().GetPointData() 
  imag_array = vtkNumPy.vtk_to_numpy(imag_image.GetArray(0)) 

  vtkAppend = vtk.vtkImageAppendComponents()
  vtkAppend.SetInput( 0,vtkRealDcmReader.GetOutput() )
  vtkAppend.SetInput( 1,vtkImagDcmReader.GetOutput() )
  vtkAppend.Update( )

  # write raw data
  vtkRawData = vtkAppend.GetOutput()
  vtkRawData.SetSpacing( spacing )
  vtkDcmWriter = vtk.vtkDataSetWriter()
  vtkDcmWriter.SetFileTypeToBinary()
  vtkDcmWriter.SetFileName("s%d/rawdata.%04d.vtk" % (dirID,fileID) )
  vtkDcmWriter.SetInput( vtkRawData )
  vtkDcmWriter.Update()

  # write raw phase data
  vtkPhase = vtk.vtkImageMathematics()
  vtkPhase.SetInput1( vtkRealData.GetOutput() )
  vtkPhase.SetInput2( vtkImagData.GetOutput() )
  vtkPhase.SetOperationToATAN2( )
  vtkPhase.Update( )
  vtkPhaseData = vtkPhase.GetOutput()
  vtkPhaseData.SetSpacing( spacing )
  vtkDcmWriter = vtk.vtkDataSetWriter()
  vtkDcmWriter.SetFileTypeToBinary()
  vtkDcmWriter.SetFileName("s%d/phase.%04d.vtk" % (dirID,fileID) )
  vtkDcmWriter.SetInput( vtkPhaseData)
  vtkDcmWriter.Update()

  return (real_array,imag_array)

    def write(self, mesh_points, filename):
        """
        Writes a vtk file, called filename, copying all the
        structures from self.filename but the coordinates.
        `mesh_points` is a matrix that contains the new coordinates
        to write in the vtk file.

        :param numpy.ndarray mesh_points: it is a `n_points`-by-3
            matrix containing the coordinates of the points of the
            mesh
        :param string filename: name of the output file.
        """
        self._check_filename_type(filename)
        self._check_extension(filename)
        self._check_infile_instantiation()

        self.outfile = filename

        reader = vtk.vtkDataSetReader()
        reader.SetFileName(self.infile)
        reader.ReadAllVectorsOn()
        reader.ReadAllScalarsOn()
        reader.Update()
        data = reader.GetOutput()

        points = vtk.vtkPoints()
        for i in range(data.GetNumberOfPoints()):
            points.InsertNextPoint(mesh_points[i, :])

        data.SetPoints(points)

        writer = vtk.vtkDataSetWriter()
        writer.SetFileName(self.outfile)
        writer.SetInputData(data)
        writer.Write()

def GetRawDICOMData(filenames,fileID):
  print filenames,fileID
  vtkRealDcmReader = vtk.vtkDICOMImageReader()
  vtkRealDcmReader.SetFileName("%s/%s"%(rootdir,filenames[0]) )
  vtkRealDcmReader.Update()
  vtkRealData = vtk.vtkImageCast()
  vtkRealData.SetOutputScalarTypeToFloat()
  vtkRealData.SetInput( vtkRealDcmReader.GetOutput() )
  vtkRealData.Update( )
  real_image = vtkRealData.GetOutput().GetPointData() 
  real_array = vtkNumPy.vtk_to_numpy(real_image.GetArray(0)) 

  vtkImagDcmReader = vtk.vtkDICOMImageReader()
  vtkImagDcmReader.SetFileName("%s/%s"%(rootdir,filenames[1]) )
  vtkImagDcmReader.Update()
  vtkImagData = vtk.vtkImageCast()
  vtkImagData.SetOutputScalarTypeToFloat()
  vtkImagData.SetInput( vtkImagDcmReader.GetOutput() )
  vtkImagData.Update( )
  imag_image = vtkImagData.GetOutput().GetPointData() 
  imag_array = vtkNumPy.vtk_to_numpy(imag_image.GetArray(0)) 

  vtkAppend = vtk.vtkImageAppendComponents()
  vtkAppend.SetInput( 0,vtkRealDcmReader.GetOutput() )
  vtkAppend.SetInput( 1,vtkImagDcmReader.GetOutput() )
  vtkAppend.Update( )

  vtkDcmWriter = vtk.vtkDataSetWriter()
  vtkDcmWriter.SetFileName("rawdata.%04d.vtk" % fileID )
  vtkDcmWriter.SetInput(vtkAppend.GetOutput())
  vtkDcmWriter.Update()

  return (real_array,imag_array)

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self,
         module_manager,
         vtk.vtkDataSetWriter(),
         "Writing vtkDataSet.",
         ("vtkDataSet",),
         (),
         replaceDoc=True,
         inputFunctions=None,
         outputFunctions=None,
     )

def dump2VTK(obj,fnm=None):
  global dumps
  if fnm is None:
    fnm="foo%.3i.vtk" % dumps
    dumps+=1
  dsw = vtk.vtkDataSetWriter()
  dsw.SetFileName(fnm)
  try:
    dsw.SetInputData(obj)
  except:
    dsw.SetInputConnection(obj.GetOutputPort())

  dsw.Write()

 def write(self):
     writer=vtk.vtkDataSetWriter()
     writer.SetInput(self.reslice.GetOutput())
     
     #proboje wczytac workspace z pliku, jesli sie nie uda to otwiera folder w ktorym sie znajduje plik
     try:
         dir=ReadFile().read_variable('output_folder:')
         print dir
     except:
         dir=""
     self.filename=asksaveasfilename(initialdir=dir,filetypes=[("allfiles","*"),("VTKfiles","*.vtk")])
     writer.SetFileName(self.filename)
     writer.Write() 

	def write(self, output_values, filename, output_name=None, write_bin=False):
		"""
		Writes a mat file, called filename. output_values is a matrix that contains the new values of the output 
		to write in the mat file.

		:param numpy.ndarray output_values: it is a `n_points`-by-1 matrix containing the values of the chosen output.
		:param string filename: name of the output file.
		:param string output_name: name of the output of interest inside the mat file. 
			If it is not passed, it is equal to self.output_name.
		:param bool write_bin: flag to write in the binary format. Default is False.
		"""
		
		self._check_filename_type(filename)
		self._check_extension(filename)
		self._check_infile_instantiation(self.infile)
		
		if output_name is None:
			output_name = self.output_name
		else:
			self._check_filename_type(output_name)
		
		reader = vtk.vtkDataSetReader()
		reader.SetFileName(self.infile)
		reader.ReadAllVectorsOn()
		reader.ReadAllScalarsOn()
		reader.Update()
		data = reader.GetOutput()
	
		output_array = ns.numpy_to_vtk(num_array=output_values,array_type=vtk.VTK_DOUBLE)
		output_array.SetName(output_name)
		
		if self.cell_data is True:
			data.GetCellData().AddArray(output_array)
		else:
			data.GetPointData().AddArray(output_array)	
	
		writer = vtk.vtkDataSetWriter()
		
		if write_bin:
			writer.SetFileTypeToBinary()
		
		writer.SetFileName(filename)
		
		if vtk.VTK_MAJOR_VERSION <= 5:
			writer.SetInput(data)
		else:
			writer.SetInputData(data)
	
		writer.Write()

def TestXdmfConversion(dataInput, fileName):
  global CleanUpGood, timer
  fileName = OutputDir + fileName
  xdmfFile = fileName + ".xmf"
  hdf5File = fileName + ".h5"
  vtkFile = fileName + ".vtk"

  xWriter = vtk.vtkXdmf3Writer()
  xWriter.SetLightDataLimit(LightDataLimit)
  xWriter.WriteAllTimeStepsOn()
  xWriter.SetFileName(xdmfFile)
  xWriter.SetInputData(dataInput)
  timer.StartTimer()
  xWriter.Write()
  timer.StopTimer()
  print "vtkXdmf3Writer took", timer.GetElapsedTime(), "seconds to write",\
    xdmfFile

  ds = vtk.vtkDataSet.SafeDownCast(dataInput)
  if ds:
    dsw = vtk.vtkDataSetWriter()
    dsw.SetFileName(vtkFile)
    dsw.SetInputData(ds)
    dsw.Write()

  if not DoFilesExist(xdmfFile, None, None, False):
    message = "Writer did not create " + xdmfFile
    raiseErrorAndExit(message)

  xReader = vtk.vtkXdmf3Reader()
  xReader.SetFileName(xdmfFile)
  timer.StartTimer()
  xReader.Update()
  timer.StopTimer()
  print "vtkXdmf3Reader took", timer.GetElapsedTime(), "seconds to read",\
    xdmfFile

  rOutput = xReader.GetOutputDataObject(0)

  fail = DoDataObjectsDiffer(dataInput, rOutput)

  if fail:
    raiseErrorAndExit("Xdmf conversion test failed")
  else:
    if ds:
      DoFilesExist(xdmfFile, hdf5File, vtkFile, CleanUpGood)
    else:
      DoFilesExist(xdmfFile, hdf5File, None, CleanUpGood)

def WriteVTKPoints(vtkpoints,OutputFileName):
   # loop over points an store in vtk data structure
   #vtkpoints = vtk.vtkPoints()
   vertices= vtk.vtkCellArray()
   for idpoint in range(vtkpoints.GetNumberOfPoints()):
       #vertices.InsertNextCell( 1 ); vertices.InsertCellPoint( vtkpoints.InsertNextPoint(point) )
       vertices.InsertNextCell( 1 ); vertices.InsertCellPoint( idpoint )

   # set polydata
   polydata = vtk.vtkPolyData()
   polydata.SetPoints(vtkpoints)
   polydata.SetVerts( vertices )

   # write to file
   polydatawriter = vtk.vtkDataSetWriter()
   polydatawriter.SetFileName(OutputFileName)
   polydatawriter.SetInput(polydata)
   polydatawriter.Update()

def dump2VTK(obj,fnm=None):
  global dumps
  if fnm is None:
    fnm="foo.vtk" % dumps
  if fnm[:-4].lower()!=".vtk":
    fnm+=".vtk"
  if fnm in dumps:
    dumps[fnm]+=1
    fnm=fnm[:-4]+"%.3i.vtk" % dumps[fnm]
  else:
    dumps[fnm]=0
  dsw = vtk.vtkDataSetWriter()
  dsw.SetFileName(fnm)
  try:
    dsw.SetInputData(obj)
  except:
    dsw.SetInputConnection(obj.GetOutputPort())

  dsw.Write()

    def _save_polydata(self, data, write_bin=False):
        """
        This private method saves into `filename` the `data`. `data` is a
        vtkPolydata. It is possible to specify format for `filename`: if
        `write_bin` is True, file is written in binary format, otherwise in
        ASCII format. This method save cached polydata to reduce number of IO
        operations.

        :param vtkPolyData data: polydatat to save.
        :param bool write_bin: for binary format file.
        """
        self._cached_data = data

        writer = vtk.vtkDataSetWriter()

        if write_bin:
            writer.SetFileTypeToBinary()

        writer.SetFileName(self._filename)
        writer.SetInputData(data)
        writer.Write()

def WriteVTKTemplateImage( TemplateFilename ):
  import vtk
  import vtk.util.numpy_support as vtkNumPy 
  import numpy
  # set Image Template Dimensions
  femBounds  = (0.0,0.04,-0.04, 0.04, -0.03,0.06)
  origin = (femBounds[0], femBounds[2], femBounds[4])
  spacing = ( (femBounds[1]-femBounds[0])/ imageDimensions[0] ,
              (femBounds[3]-femBounds[2])/ imageDimensions[1] ,
              (femBounds[5]-femBounds[4])/ imageDimensions[2]  
            )
  print femBounds, origin, spacing
  # imports raw data and stores it.
  dataImporter = vtk.vtkImageImport()
  # array is converted to a string of chars and imported.
  # numpy array stored as ROW MAJOR
  # MUST write out in COLUMN MAJOR format to be the same as VTK
  data_string = numpy.zeros(imageDimensions,dtype=numpy.float,order='F').tostring(order='F')
  dataImporter.CopyImportVoidPointer(data_string, len(data_string))
  # The type of the newly imported data is set to unsigned char (uint8)
  dataImporter.SetDataScalarTypeToDouble()
  # Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer
  # must be told this is the case.
  dataImporter.SetNumberOfScalarComponents(1)
  # The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this
  # simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case.
  # I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although
  # VTK complains if not both are used.
  dataImporter.SetDataExtent( 0, imageDimensions[0]-1, 0, imageDimensions[1]-1, 0, imageDimensions[2]-1)
  dataImporter.SetWholeExtent(0, imageDimensions[0]-1, 0, imageDimensions[1]-1, 0, imageDimensions[2]-1)
  dataImporter.SetDataSpacing( spacing )
  dataImporter.SetDataOrigin(  origin  )
  dataImporter.SetScalarArrayName(  "scalars" )
  dataImporter.Update()
  vtkTemplateWriter = vtk.vtkDataSetWriter()
  vtkTemplateWriter.SetFileName( TemplateFilename )
  vtkTemplateWriter.SetInput( dataImporter.GetOutput() )
  vtkTemplateWriter.Update()
  return 

def WriteVTKPoints(self,vtkpoints,OutputFileName):
   # loop over points an store in vtk data structure
   # write in meters
   MillimeterMeterConversion = .001;
   scalevtkPoints = vtk.vtkPoints()
   vertices= vtk.vtkCellArray()
   for idpoint in range(vtkpoints.GetNumberOfPoints()):
       point = MillimeterMeterConversion * numpy.array(vtkpoints.GetPoint(idpoint))
       vertices.InsertNextCell( 1 ); vertices.InsertCellPoint( scalevtkPoints.InsertNextPoint(point) )
       #vertices.InsertNextCell( 1 ); vertices.InsertCellPoint( idpoint )

   # set polydata
   polydata = vtk.vtkPolyData()
   polydata.SetPoints(scalevtkPoints )
   polydata.SetVerts( vertices )

   # write to file
   print "WriteVTKPoints: writing",OutputFileName
   polydatawriter = vtk.vtkDataSetWriter()
   polydatawriter.SetFileName(OutputFileName)
   polydatawriter.SetInput(polydata)
   polydatawriter.Update()

LandmarkTransform = vtk.vtkLandmarkTransform()
LandmarkTransform.SetSourceLandmarks(SourceLMReader.GetOutput().GetPoints() )
LandmarkTransform.SetTargetLandmarks(TargetLMReader.GetOutput().GetPoints() )
LandmarkTransform.SetModeToRigidBody()
LandmarkTransform.Update()
print LandmarkTransform.GetMatrix()

# apply transform
transformFilter = vtk.vtkTransformFilter()
#transformFilter.SetInput(vtkCylinder.GetOutput() ) 
transformFilter.SetInput(trianglefilter.GetOutput() ) 
transformFilter.SetTransform( LandmarkTransform) 
transformFilter.Update()

# write model
modelWriter = vtk.vtkDataSetWriter()
modelWriter.SetInput(transformFilter.GetOutput())
modelWriter.SetFileName("needle.vtk")
modelWriter.SetFileTypeToBinary()
modelWriter.Update()

# read image/ROI
ImageReader = vtk.vtkDataSetReader()
ImageReader.SetFileName("newimage.vtk")
ImageReader.Update()

# resample needle to image to create mask
vtkResample = vtk.vtkCompositeDataProbeFilter()
vtkResample.SetInput( ImageReader.GetOutput() )
vtkResample.SetSource( transformFilter.GetOutput() ) 
vtkResample.Update()

#.........这里部分代码省略.........
  # vtkVOIExtract.SetVOI( VOI )
  # vtkVOIExtract.Update()
  # mrti_point_data= vtkVOIExtract.GetOutput().GetPointData()
  # mrti_array = vtkNumPy.vtk_to_numpy(mrti_point_data.GetArray('scalars'))
  # #print mrti_array
  # #print type(mrti_array)

  # Interpolate FEM onto imaging data structures
  vtkExodusIIReader = vtk.vtkExodusIIReader()
  vtkExodusIIReader.SetFileName( fem_mesh_file )
  vtkExodusIIReader.SetPointResultArrayStatus("u0",1)
  vtkExodusIIReader.SetPointResultArrayStatus("u0*",1)
  vtkExodusIIReader.SetPointResultArrayStatus("u1",1)

  matsize = int(dimensions[1])#get matrix size
  #preallocate size of arrays
  u0_array_1 = scipy.zeros((matsize,matsize))
  u0_array_2 = scipy.zeros((matsize,matsize,ntime*nsubstep))
  u1_array_1 = scipy.zeros((matsize,matsize))
  u1_array_2 = scipy.zeros((matsize,matsize,ntime*nsubstep))
  u0star_array_1 = scipy.zeros((matsize,matsize))
  u0star_array_2 = scipy.zeros((matsize,matsize,ntime*nsubstep))

  #for timeID in range(1,2):
  for timeID in range(1,ntime*nsubstep):
    vtkExodusIIReader.SetTimeStep(timeID-1)
    vtkExodusIIReader.Update()
    
    # apply the transform
    TransformedFEMMesh = None
    if vtkExodusIIReader.GetOutput().IsA("vtkMultiBlockDataSet"):
      AppendBlocks = vtk.vtkAppendFilter()
      iter = vtkExodusIIReader.GetOutput().NewIterator()
      iter.UnRegister(None)
      iter.InitTraversal()
      # loop over blocks...
      while not iter.IsDoneWithTraversal():
          curInput = iter.GetCurrentDataObject()
          vtkTransformFEMMesh = vtk.vtkTransformFilter()
          vtkTransformFEMMesh.SetTransform( AffineTransform )
          vtkTransformFEMMesh.SetInput( curInput )
          vtkTransformFEMMesh.Update()
          AppendBlocks.AddInput( vtkTransformFEMMesh.GetOutput() )
          AppendBlocks.Update( )
          iter.GoToNextItem();
      TransformedFEMMesh = AppendBlocks.GetOutput()
    else:
      vtkTransformFEMMesh = vtk.vtkTransformFilter()
      vtkTransformFEMMesh.SetTransform( AffineTransform )
      vtkTransformFEMMesh.SetInput( vtkExodusIIReader.GetOutput() )
      vtkTransformFEMMesh.Update()
      TransformedFEMMesh = vtkTransformFEMMesh.GetOutput()

    # reflect
    #vtkReflectX = vtk.vtkReflectionFilter()
    #vtkReflectX.SetPlaneToXMin()
    #vtkReflectX.SetInput( TransformedFEMMesh )
    #vtkReflectX.Update()

    # reflect
    #vtkReflectZ = vtk.vtkReflectionFilter()
    #vtkReflectZ.SetPlaneToZMax()
    #vtkReflectZ.SetInput( vtkReflectX.GetOutput() )
    #vtkReflectZ.Update()

    # reuse ShiftScale Geometry
    vtkResample = vtk.vtkCompositeDataProbeFilter()
    vtkResample.SetInput( templateImage )
    vtkResample.SetSource( TransformedFEMMesh )
    #vtkResample.SetSource( vtkReflectZ.GetOutput() )
    vtkResample.Update()
    fem_point_data= vtkResample.GetOutput().GetPointData()
    u0_array = vtkNumPy.vtk_to_numpy(fem_point_data.GetArray('u0'))
    u0star_array = vtkNumPy.vtk_to_numpy(fem_point_data.GetArray('u0*'))
    u1_array = vtkNumPy.vtk_to_numpy(fem_point_data.GetArray('u1'))
 
    
    #go from 1x256^2 array to 256X256 array for each timestep
    for nn in range(0,matsize):
     u0_array_1[nn,:]=u0_array[nn*matsize:(nn+1)*matsize]
     u0star_array_1[nn,:]=u0star_array[nn*matsize:(nn+1)*matsize]
     u1_array_1[nn,:]=u1_array[nn*matsize:(nn+1)*matsize]
   
    #apply proper rotations/reflections and combine 2D arrays into 3D array 
    u0_array_2[:,:,timeID-1]=numpy.fliplr(numpy.rot90(u0_array_1,k=3))
    u0star_array_2[:,:,timeID-1]=numpy.fliplr(numpy.rot90(u0star_array_1,k=3))
    u1_array_2[:,:,timeID-1]=numpy.fliplr(numpy.rot90(u1_array_1,k=3))
    
    # write numpy to disk in matlab
    #scipyio.savemat("MS795.%04d.mat" % (timeID), {'u0':u1_array,'MRTI0':MRTI0_array })

    # write output
    print "writing ", timeID
    vtkStatsWriter = vtk.vtkDataSetWriter()
    vtkStatsWriter.SetFileTypeToBinary()
    vtkStatsWriter.SetFileName("test.%04d.vtk" % timeID )
    vtkStatsWriter.SetInput(vtkResample.GetOutput())
    vtkStatsWriter.Update()

  scipyio.savemat("S695.mat",{'ModelFluence':u1_array_2,'MRTI':u0star_array_2,'ModelTemp':u0_array_2})

#.........这里部分代码省略.........
  simplicity but there are a number of disadvantages and caveats as well as
  described in this section.
  
  The simple array files are made up of a) a header followed by b) the data
  itself. This layout of data and header is illustrated below. The header has a
  single 32-bit integer to indicate the number of dimensions of the dataset
  followed by one integer for each dimension to indicate the length of that
  dimension. The data follows immediately after the header. The data is stored
  such that the first dimension is the fastest moving dimension, second dimension
  is second fastest, etc. The header contains no information about the size of
  each individual data element and consequently the user needs to know what type
  of data is contained in the array. In general, the Gadgetron uses 3 different
  types of data and the convention is to use the file extension to indicate the
  data type in the file:
  
     16-bit unsigned short. File extension: *.short
     32-bit float. File extension: *.real
     32-bit complex float. Two 32-bit floating point values per data element. File extension: *.cplx
  
  The Gadgetron framework provides function for reading these files in C++. The
  functions are located in toolboxes/ndarray/hoNDArray_fileio.h in the Gadgetron
  source code distribution.
  
  It is also trivial to read the files into Matlab. Below is a function which
  detects the data type based on the file extension and reads the file into
  Matlab.
  
  """
    import vtk.util.numpy_support as vtkNumPy
    import numpy
    import scipy.io as scipyio

    # read the header
    imagedimension = numpy.fromfile(input_filename, dtype=numpy.int32, count=1, sep="")[0]
    fileheader = numpy.fromfile(input_filename, dtype=numpy.int32, count=1 + imagedimension, sep="")
    print imagedimension, fileheader
    if imagedimension == 1:
        dims = [fileheader[1], 1, 1, 1]
    elif imagedimension == 2:
        dims = [fileheader[1], fileheader[2], 1, 1]
    elif imagedimension == 3:
        dims = [fileheader[1], fileheader[2], fileheader[3], 1]
    else:
        raise RuntimeError("unknown dimension %d " % imagedimension)

    # the extension is the datatype
    extension = input_filename.split(".").pop()
    dataImporter = vtk.vtkImageImport()
    WriteImageData = True
    if extension == "short":
        datatype = numpy.short
        dataImporter.SetDataScalarTypeToShort()
    elif extension == "real":
        datatype = numpy.float32
        dataImporter.SetDataScalarTypeToFloat()
    elif extension == "cplx":
        datatype = numpy.float32
        dataImporter.SetDataScalarTypeToFloat()
        dims[3] = 2
        WriteImageData = False
    else:
        raise RuntimeError("unknown data type %s " % extension)

    # offset the data read by the header
    datafile = open(input_filename, "rb")  # reopen the file
    headeroffset = len(fileheader) * 4  # 4 bytes per integer
    datafile.seek(headeroffset, os.SEEK_SET)  # seek
    numpy_data = numpy.fromfile(datafile, dtype=datatype, sep="")

    # error check
    ExpectedImageSize = 1
    for pixelsize in dims:
        ExpectedImageSize = ExpectedImageSize * pixelsize
    if ExpectedImageSize != numpy_data.size:
        raise RuntimeError("file read error: expected size %d, size found %d " % (ExpectedImageSize, numpy_data.size))

    # FIXME Hack for trajectory Points
    if WriteImageData:
        # convert to vtk
        spacing = [1.0, 1.0, 1.0]
        dataImporter.SetNumberOfScalarComponents(dims[3])
        dataImporter.SetDataExtent(0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1)
        dataImporter.SetWholeExtent(0, dims[0] - 1, 0, dims[1] - 1, 0, dims[2] - 1)
        dataImporter.SetDataSpacing(spacing[0], spacing[1], spacing[2])
        # numpy_data= numpy_data.reshape(dims[0],dims[1],dims[2])
        # numpy_data= numpy_data.transpose(1,0,2)
        data_string = numpy_data.tostring()
        dataImporter.CopyImportVoidPointer(data_string, len(data_string))
        dataImporter.SetScalarArrayName(input_filename.split(".").pop(0))

        # write vtk file
        print "writing ", output_filename
        vtkImageDataWriter = vtk.vtkDataSetWriter()
        vtkImageDataWriter.SetFileTypeToBinary()
        vtkImageDataWriter.SetFileName(output_filename)
        vtkImageDataWriter.SetInput(dataImporter.GetOutput())
        vtkImageDataWriter.Update()
    else:
        # TODO need to correct spacing
        WriteVTKPoints(numpy_data.reshape(dims[0], dims[3]), output_filename)

    Vertices.InsertNextCell(1)
    Vertices.InsertCellPoint(pointID)

  Polydata.SetPoints(Points)
  Polydata.SetVerts(Vertices)
  Polydata.Modified()
  if vtk.VTK_MAJOR_VERSION <= 5: Polydata.Update()
 
# ------------------------------------------ output result ---------------------------------------  

  if name:
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetCompressorTypeToZLib()
    writer.SetDataModeToBinary()
    writer.SetFileName(os.path.join(os.path.split(name)[0],
                                    os.path.splitext(os.path.split(name)[1])[0] +
                                    '.' + writer.GetDefaultFileExtension()))
  else:
    writer = vtk.vtkDataSetWriter()
    writer.SetHeader('# powered by '+scriptID)
    writer.WriteToOutputStringOn()
  
  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
  else:                          writer.SetInputData(Polydata)

  writer.Write()

  if name is None: sys.stdout.write(writer.GetOutputString()[:writer.GetOutputStringLength()])      # limiting of outputString is fix for vtk <7.0

  table.close()

   femImaging.ProjectImagingToFEMMesh("Background",0.0,v2,eqnSystems)  
   femImaging.ProjectImagingToFEMMesh("ImageMask" ,0.0,v3,eqnSystems)  
   bgSystem.SystemSolve( ) 
   exodusII_IO.WriteTimeStep(MeshOutputFile,eqnSystems, timeID+1, timeID )  
   # FIXME: The FEM ordering does not match the imaging
   # FIXME:  Need to use exodus file for now
   # write numpy to disk in matlab
   fem_orig_array = origSystem.GetSolutionVector( )[...]
   fem_drop_array = bgSystem.GetSolutionVector( )[...]
   fem_mask_array = maskSystem.GetSolutionVector( )[...]
   scipyio.savemat("Processed/background.%04d.mat"%(timeID), {'pixelsize':npixelROI, 'original':fem_orig_array , 'drop':fem_drop_array , 'mask':fem_mask_array } )

   # write numpy to vtk 
   ROIOrigin = (xbounds[0],ybounds[0],zbounds[0])
   vtkOrigImage = ConvertNumpyVTKImage(fem_orig_array, "orig" , npixelROI, ROIOrigin )
   vtkOrigWriter = vtk.vtkDataSetWriter()
   vtkOrigWriter.SetFileName("Processed/original.%04d.vtk" % timeID )
   vtkOrigWriter.SetInput( vtkOrigImage )
   vtkOrigWriter.Update()
   
   vtkDropImage = ConvertNumpyVTKImage(fem_drop_array, "drop" , npixelROI, ROIOrigin )
   vtkDropWriter = vtk.vtkDataSetWriter()
   vtkDropWriter.SetFileName("Processed/drop.%04d.vtk" % timeID )
   vtkDropWriter.SetInput( vtkDropImage )
   vtkDropWriter.Update()

   ## FIXME: bug putting data back into imaging data structures
   ##  # get libMesh Background Solution as numpy data structure
   ##  maxwell_data  = phase_curr.copy()
   ##  # reshape from colume major Fortran-like storage
   ##  maxwell_data[ROI[0][0]+0:ROI[0][1]+0,

def write(data,filename):
    writer = vtk.vtkDataSetWriter()
    writer.SetFileName(filename + ".vtk")
    writer.SetInputData(data)
    writer.Write()
    print "Data written to file:\t " + filename +".vtk"

  def __init__(self, SEMDataDirectory,variableDictionary  ):
     
     self.DataDictionary = {}
     self.DebugObjective = True
     self.DebugObjective = False
     self.ctx   = cl.create_some_context()
     self.queue = cl.CommandQueue(self.ctx)
     self.prg   = cl.Program(self.ctx, """
              __kernel void diff_sq(__global const float *a,
              __global const float *b, __global float *c)
              {
                int gid = get_global_id(0);
                c[gid] = (a[gid] - b[gid]) * (a[gid] - b[gid]);
              }
              """).build()

     # FIXME  should this be different ?  
     self.SEMDataDirectory = SEMDataDirectory 
  
     # FIXME vtk needs to be loaded AFTER kernel is built
     import vtk
     import vtk.util.numpy_support as vtkNumPy 
     print "using vtk version", vtk.vtkVersion.GetVTKVersion()
     print "read SEM data"

     start = time.clock()
     vtkSEMReader = vtk.vtkXMLUnstructuredGridReader()
     vtufileName = "%s/%d.vtu" % (self.SEMDataDirectory,0)
     vtkSEMReader.SetFileName( vtufileName )
     vtkSEMReader.SetPointArrayStatus("Temperature",1)
     vtkSEMReader.Update()
     elapsed = (time.clock() - start)
     print "read SEM data", elapsed
  
     # get registration parameters
  
     # register the SEM data to MRTI
     AffineTransform = vtk.vtkTransform()
     AffineTransform.Translate([ 
       float(variableDictionary['x_displace']),
       float(variableDictionary['y_displace']),
       float(variableDictionary['z_displace'])
                               ])
     # FIXME  notice that order of operations is IMPORTANT
     # FIXME   translation followed by rotation will give different results
     # FIXME   than rotation followed by translation
     # FIXME  Translate -> RotateZ -> RotateY -> RotateX -> Scale seems to be the order of paraview
     AffineTransform.RotateZ( float(variableDictionary['z_rotate'  ] ) ) 
     AffineTransform.RotateY( float(variableDictionary['y_rotate'  ] ) )
     AffineTransform.RotateX( float(variableDictionary['x_rotate'  ] ) )
     AffineTransform.Scale([1.e0,1.e0,1.e0])
     self.SEMRegister = vtk.vtkTransformFilter()
     self.SEMRegister.SetInput(vtkSEMReader.GetOutput())
     self.SEMRegister.SetTransform(AffineTransform)
     self.SEMRegister.Update()
  
     print "write transform output"
     if ( self.DebugObjective ):
        vtkSEMWriter = vtk.vtkDataSetWriter()
        vtkSEMWriter.SetFileTypeToBinary()
        semfileName = "%s/semtransform%04d.vtk" % (self.SEMDataDirectory,0)
        print "writing ", semfileName 
        vtkSEMWriter.SetFileName( semfileName )
        vtkSEMWriter.SetInput(self.SEMRegister.GetOutput())
        vtkSEMWriter.Update()