def Execute(self):

        if self.Image == None:
            self.PrintError('Error: No input image.')

        if self.Image2 == None:
            self.PrintError('Error: No input image2.')

        if self.NegateImage2:
            negateFilter = vtk.vtkImageMathematics()
            negateFilter.SetInputData(self.Image2)
            negateFilter.SetOperationToMultiplyByK()
            negateFilter.SetConstantK(-1.0)
            negateFilter.Update()
            self.Image2 = negateFilter.GetOutput()

        composeFilter = vtk.vtkImageMathematics()
        composeFilter.SetInput1Data(self.Image)
        composeFilter.SetInput2Data(self.Image2)
        if self.Operation == 'min':
            composeFilter.SetOperationToMin()
        elif self.Operation == 'max':
            composeFilter.SetOperationToMax()
        elif self.Operation == 'multiply':
            composeFilter.SetOperationToMultiply()
        elif self.Operation == 'subtract':
            composeFilter.SetOperationToSubtract()
        else:
            self.PrintError('Error: Unsupported operation')
        composeFilter.Update()

        self.Image = composeFilter.GetOutput()

    def IsosurfaceInitialize(self):

        self.PrintLog('Isosurface initialization.')

        if self.Interactive:
            queryString = "Please input isosurface level (\'n\' for none): "
            self.IsoSurfaceValue = self.ThresholdInput(queryString)
        
        imageMathematics = vtk.vtkImageMathematics()
        imageMathematics.SetInput(self.Image)
        imageMathematics.SetConstantK(-1.0)
        imageMathematics.SetOperationToMultiplyByK()
        imageMathematics.Update()

        subtract = vtk.vtkImageMathematics()
        subtract.SetInput(imageMathematics.GetOutput())
        subtract.SetOperationToAddConstant()
        subtract.SetConstantC(self.IsoSurfaceValue)
        subtract.Update()

        self.InitialLevelSets = vtk.vtkImageData()
        self.InitialLevelSets.DeepCopy(subtract.GetOutput())
        self.InitialLevelSets.Update()

        self.IsoSurfaceValue = 0.0

    def ComputeKSpaceOperation(self,kSpace):

        kSpaceMValue = vtk.vtkImageMathematics()
        kSpaceMValue.SetInput(kSpace)
        kSpaceMValue.SetOperationToMultiplyByK()
        kSpaceMValue.SetConstantK(self.MValue)
        kSpaceMValue.Update()

        kSpace = kSpaceMValue.GetOutput()

        if not self.KSpace:
            return kSpace
 
        if self.KSpaceOperation == 'none':
            return kSpace
       
        kSpaceMaths = vtk.vtkImageMathematics()
        if self.KSpaceOperation == 'add':
            kSpaceMaths.SetOperationToAdd()
        elif self.KSpaceOperation == 'subtract':
            kSpaceMaths.SetOperationToSubtract()
        else:
            self.PrintError('KSpaceOperation not supported.')
            return kSpace
            
        kSpaceMaths.SetInput1(self.KSpace)
        kSpaceMaths.SetInput2(kSpace)
        kSpaceMaths.Update()

        return kSpaceMaths.GetOutput()

    def BuildVTKGradientBasedFeatureImage(self):

        cast = vtk.vtkImageCast()
        cast.SetInputData(self.Image)
        cast.SetOutputScalarTypeToFloat()
        cast.Update()

        gradientMagnitude = vtk.vtkImageGradientMagnitude()
        gradientMagnitude.SetInputConnection(cast.GetOutputPort())
        gradientMagnitude.SetDimensionality(self.Dimensionality)
        gradientMagnitude.Update()

        imageAdd = vtk.vtkImageMathematics()
        imageAdd.SetInputConnection(gradientMagnitude.GetOutputPort())
        imageAdd.SetOperationToAddConstant()
        imageAdd.SetConstantC(1.0)
        imageAdd.Update()

        imageInvert = vtk.vtkImageMathematics()
        imageInvert.SetInputConnection(imageAdd.GetOutputPort())
        imageInvert.SetOperationToInvert()
        imageInvert.SetConstantC(1E20)
        imageInvert.DivideByZeroToCOn()
        imageInvert.Update()

        self.FeatureImage = vtk.vtkImageData()
        self.FeatureImage.DeepCopy(imageInvert.GetOutput())

  def sumManualSegmentations(self, manualSegmentationsDirectory, mergedVolume):
    # Get the manual segmentations and create a single summed image
    import glob
    manualSegmentationFilenames = glob.glob(manualSegmentationsDirectory+"/*.mha")

    # Get the first image which each successive image will be added to
    reader = vtk.vtkMetaImageReader()
    reader.SetFileName(manualSegmentationFilenames[0])
    reader.Update()
    summedImage = vtk.vtkImageData()
    summedImage.SetExtent(reader.GetOutput().GetExtent())
    summedImage.AllocateScalars(vtk.VTK_UNSIGNED_CHAR,1)
    summedImage.ShallowCopy(reader.GetOutput())

    # Initialize filter to add images together
    mathFilter = vtk.vtkImageMathematics()

    # Iterate list and add each new image
    for currentFile in manualSegmentationFilenames[1:]:
      # Get new image
      reader.SetFileName(currentFile)
      reader.Update()

      # Add it to existing summation
      mathFilter.SetInput1Data(summedImage)
      mathFilter.SetInput2Data(reader.GetOutput())
      mathFilter.Update()

      # Get new summation
      summedImage.ShallowCopy(mathFilter.GetOutput())

    # Add summed image to slicer scene
    mergedVolume.SetRASToIJKMatrix(self.rasToIjk)
    mergedVolume.SetIJKToRASMatrix(self.ijkToRas)
    mergedVolume.SetAndObserveImageData(summedImage)

    def Execute(self):

        if self.Image == None:
            self.PrintError('Error: No input image.')

        if self.LevelSets == None:
            self.PrintError('Error: No input level sets.')

        if self.NegateLevelSets:
            negateFilter = vtk.vtkImageMathematics()
            negateFilter.SetInputData(self.LevelSets)
            negateFilter.SetOperationToMultiplyByK()
            negateFilter.SetConstantK(-1.0)
            negateFilter.Update()
            self.LevelSets = negateFilter.GetOutput()

        sigmoid = vtkvmtk.vtkvmtkLevelSetSigmoidFilter()
        sigmoid.SetInputData(self.Image)
        sigmoid.SetLevelSetsImage(self.LevelSets)
        sigmoid.SetSigma(self.Sigma)
        sigmoid.SetScaleValue(self.ScaleValue)
        if self.ComputeScaleValueFromInput:
            sigmoid.ComputeScaleValueFromInputOn()
        else:
            sigmoid.ComputeScaleValueFromInputOff()
        sigmoid.Update()
 
        self.Image = sigmoid.GetOutput()

    def subtractionCompare(self):

        imagePoints = self.Image.GetNumberOfPoints()
        referencePoints = self.ReferenceImage.GetNumberOfPoints()

        self.InputInfo('Image Points: ' + str(imagePoints))
        self.InputInfo('Reference Image Points: ' + str(referencePoints))

        if abs(imagePoints - referencePoints) > 0 :
            self.ResultLog = 'Uneven NumberOfPoints'
            return False
            
        imageScalarType = self.Image.GetScalarType()
        referenceScalarType = self.ReferenceImage.GetScalarType()
        if imageScalarType != referenceScalarType:
            self.PrintError('Error: Input image and reference image are not of \
                the same type. Please cast images to the same type.')

        imageMath = vtk.vtkImageMathematics()
        imageMath.SetInput1Data(self.Image) 
        imageMath.SetInput2Data(self.ReferenceImage) 
        imageMath.SetOperationToSubtract()
        imageMath.Update()
        differenceImage = imageMath.GetOutput()
        differenceRange = differenceImage.GetPointData().GetScalars().GetRange()

        self.InputInfo('Difference Range: ' + str(differenceRange))
 
        if max([abs(d) for d in differenceRange]) < self.Tolerance:
            return True
        
        return False

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkImageMathematics(), 'Processing.',
         ('vtkImageData', 'vtkImageData'), ('vtkImageData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

def getIndexedImageData(inputVtkImageData, indexVal):
    """
    @type  inputVtkImageData: C{vtkImageData}
    @param inputVtkImageData: Image data to be processed
    
    @type  indexVal: C{unsigned int}
    @param indexVal: Index value that will be assigned to all non-zero values.
    
    Function replaces all non-zero valued pixels of provided
    C{inputVtkImageData} with C{indexVal} and returns the result.
    
    @rtype: C{vtkImageData}
    @return: UnsignedInt vtkImageData with two values 0 and provided
             C{indexVal}.
    """
    # Convert normalized input image to unsigned int scalars
    cast = vtk.vtkImageCast()
    cast.SetInput(getImageMask(inputVtkImageData))
    cast.SetOutputScalarTypeToUnsignedInt()
    
    # Then multiply the normalized image data by priovided constant
    multip = vtk.vtkImageMathematics()
    multip.SetOperationToMultiplyByK()
    multip.SetConstantK(indexVal)
    multip.SetInput(cast.GetOutput())
    return multip.GetOutput()

    def __init__(self, module_manager):
        # initialise our base class
        ModuleBase.__init__(self, module_manager)


        self._imageMath = vtk.vtkImageMathematics()
        self._imageMath.SetInput1(None)
        self._imageMath.SetInput2(None)
        
        module_utils.setup_vtk_object_progress(self, self._imageMath,
                                           'Performing image math')
        
                                           
        self._config.operation = 'subtract'
        self._config.constantC = 0.0
        self._config.constantK = 1.0


        configList = [
            ('Operation:', 'operation', 'base:str', 'choice',
             'The operation that should be performed.',
             tuple(OPS_DICT.keys())),
            ('Constant C:', 'constantC', 'base:float', 'text',
             'The constant C used in some operations.'),
            ('Constant K:', 'constantK', 'base:float', 'text',
             'The constant C used in some operations.')]

        ScriptedConfigModuleMixin.__init__(
            self, configList,
            {'Module (self)' : self,
             'vtkImageMathematics' : self._imageMath})

        self.sync_module_logic_with_config()

    def test_vtk_pyexception_deadlock(self):
        """Test if VTK has been patched to release the GIL during all
        VTK method calls.
        """

        import vtk
        # this gives floats by default
        s = vtk.vtkImageGridSource()
        c1 = vtk.vtkImageCast()
        c1.SetOutputScalarTypeToShort()
        c1.SetInput(s.GetOutput())
        c2 = vtk.vtkImageCast()
        c2.SetOutputScalarTypeToFloat()
        c2.SetInput(s.GetOutput())

        m = vtk.vtkImageMathematics()

        # make sure we are multi-threaded
        if m.GetNumberOfThreads() < 2:
            m.SetNumberOfThreads(2)
        m.SetInput1(c1.GetOutput())
        m.SetInput2(c2.GetOutput())

        # without the patch, this call will deadlock forever
        try:
            # with the patch this should generate a RuntimeError
            m.Update()
        except RuntimeError:
            pass
        else:
            self.fail(
            'Multi-threaded error vtkImageMathematics did not raise '
            'exception.')

    def _maskSourceExecute(self):
        inputI = self._inputImage
        if inputI:
            inputI.Update()
                
            self._markerSource.Update()
            outputJ = self._markerSource.GetStructuredPointsOutput()
            # we now have an outputJ

            if not inputI.GetScalarPointer() or \
               not outputJ.GetScalarPointer() or \
               not inputI.GetDimensions() > (0,0,0):
                vtk.vtkOutputWindow.GetInstance().DisplayErrorText(
                    'modifyHomotopy: Input is empty.')
                return

            iMath = vtk.vtkImageMathematics()
            iMath.SetOperationToMin()
            iMath.SetInput1(outputJ)
            iMath.SetInput2(inputI)
            iMath.GetOutput().SetUpdateExtentToWholeExtent()
            iMath.Update()

            outputI = self._maskSource.GetStructuredPointsOutput()
            outputI.DeepCopy(iMath.GetOutput())

    def subtractionCompare(self):

        imagePoints = self.Image.GetNumberOfPoints()
        referencePoints = self.ReferenceImage.GetNumberOfPoints()

        self.InputInfo('Image Points: ' + str(imagePoints))
        self.InputInfo('Reference Image Points: ' + str(referencePoints))

        if abs(imagePoints - referencePoints) > 0 :
            self.ResultLog = 'Uneven NumberOfPoints'
            return False
            
        imageScalarType = self.Image.GetScalarType()
        referenceScalarType = self.ReferenceImage.GetScalarType()
        minScalarType = min(imageScalarType,referenceScalarType)
        self.Image.SetScalarType(minScalarType) 
        self.ReferenceImage.SetScalarType(minScalarType) 

        imageMath = vtk.vtkImageMathematics()
        imageMath.SetInput1Data(self.Image) 
        imageMath.SetInput2Data(self.ReferenceImage) 
        imageMath.SetOperationToSubtract()
        imageMath.Update()
        differenceImage = imageMath.GetOutput()
        differenceRange = differenceImage.GetPointData().GetScalars().GetRange()

        self.InputInfo('Difference Range: ' + str(differenceRange))
 
        if max([abs(d) for d in differenceRange]) < self.Tolerance:
            return True
        
        return False

    def MergeLevelSet(self):

        if self.LevelSets == None:
            self.LevelSets = self.LevelSetsOutput
        else:
            minFilter = vtk.vtkImageMathematics()
            minFilter.SetOperationToMin()
            minFilter.SetInput1Data(self.LevelSets)
            minFilter.SetInput2Data(self.LevelSetsOutput)
            minFilter.Update()
            self.LevelSets = minFilter.GetOutput()

    def MergeLevelSets(self):

        if self.MergedInitialLevelSets == None:
            self.MergedInitialLevelSets = vtk.vtkImageData()
            self.MergedInitialLevelSets.DeepCopy(self.InitialLevelSets)
        else:
            minFilter = vtk.vtkImageMathematics()
            minFilter.SetOperationToMin()
            minFilter.SetInput1Data(self.MergedInitialLevelSets)
            minFilter.SetInput2Data(self.InitialLevelSets)
            minFilter.Update()
            self.MergedInitialLevelSets = minFilter.GetOutput()

    def __init__(self, module_manager):
        
        ModuleBase.__init__(self, module_manager)

        self._imageMathSubtractH = vtk.vtkImageMathematics()
        self._imageMathSubtractH.SetOperationToAddConstant()
        
        self._reconstruct = vtkdevide.vtkImageGreyscaleReconstruct3D()
        # second input is marker
        self._reconstruct.SetInput(1, self._imageMathSubtractH.GetOutput())

        self._imageMathSubtractR = vtk.vtkImageMathematics()
        self._imageMathSubtractR.SetOperationToSubtract()

        self._imageMathSubtractR.SetInput(1, self._reconstruct.GetOutput())

        module_utils.setup_vtk_object_progress(self, self._imageMathSubtractH,
                                           'Preparing marker image.')

        module_utils.setup_vtk_object_progress(self, self._reconstruct,
                                           'Performing reconstruction.')

        module_utils.setup_vtk_object_progress(self, self._imageMathSubtractR,
                                           'Subtracting reconstruction.')

        self._config.h = 50

        configList = [
            ('H-dome height:', 'h', 'base:float', 'text',
             'The required difference in brightness between an h-dome and\n'
             'its surroundings.')]

        ScriptedConfigModuleMixin.__init__(
            self, configList,
            {'Module (self)' : self,
             'ImageMath Subtract H' : self._imageMathSubtractH,
             'ImageGreyscaleReconstruct3D' : self._reconstruct,
             'ImageMath Subtract R' : self._imageMathSubtractR})

        self.sync_module_logic_with_config()

    def __init__(self, module_manager):
        # initialise our base class
        ModuleBase.__init__(self, module_manager)
        NoConfigModuleMixin.__init__(self)

        # these will be our markers
        self._inputPoints = None

        # we can't connect the image input directly to the masksource,
        # so we have to keep track of it separately.
        self._inputImage = None
        self._inputImageObserverID = None

        # we need to modify the mask (I) as well.  The problem with a
        # ProgrammableFilter is that you can't request GetOutput() before
        # the input has been set... 
        self._maskSource = vtk.vtkProgrammableSource()
        self._maskSource.SetExecuteMethod(self._maskSourceExecute)
        
        # we'll use this to synthesise a volume according to the seed points
        self._markerSource = vtk.vtkProgrammableSource()
        self._markerSource.SetExecuteMethod(self._markerSourceExecute)
        # second input is J (the marker)

        # we'll use this to change the markerImage into something we can use
        self._imageThreshold = vtk.vtkImageThreshold()
        # everything equal to or above 1.0 will be "on"
        self._imageThreshold.ThresholdByUpper(1.0)
        self._imageThresholdObserverID = self._imageThreshold.AddObserver(
            'EndEvent', self._observerImageThreshold)
        
        self._viewFrame = self._createViewFrame(
            {'Module (self)' : self})

        # we're not going to give imageErode any input... that's going to
        # to happen manually in the execute_module function :)
        self._imageErode = vtk.vtkImageContinuousErode3D()
        self._imageErode.SetKernelSize(3,3,3)

        module_utils.setup_vtk_object_progress(self, self._imageErode,
                                           'Performing greyscale 3D erosion')
        

        self._sup = vtk.vtkImageMathematics()
        self._sup.SetOperationToMax()
        self._sup.SetInput1(self._imageErode.GetOutput())
        self._sup.SetInput2(self._maskSource.GetStructuredPointsOutput())

        # pass the data down to the underlying logic
        self.config_to_logic()
        # and all the way up from logic -> config -> view to make sure
        self.syncViewWithLogic()

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 MergeLevelSet(self):

        if self.LevelSets == None:
            self.LevelSets = self.LevelSetsOutput
        else:
            minFilter = vtk.vtkImageMathematics()
            minFilter.SetOperationToMin()
            minFilter.SetInput1(self.LevelSets)
            minFilter.SetInput2(self.LevelSetsOutput)
            minFilter.Update()
            self.LevelSets = minFilter.GetOutput()

        if self.LevelSets.GetSource():
            self.LevelSets.GetSource().UnRegisterAllOutputs()

 def subImage(self, Images2Sub, timep):
     '''subtract volumes based on indicated postS'''
     sub_preMat = vtk.vtkImageMathematics()
     sub_preMat.SetOperationToSubtract()
     sub_preMat.SetInput1(Images2Sub[timep])
     sub_preMat.SetInput2(Images2Sub[0])
     sub_preMat.Update()
                 
     sub_pre = vtk.vtkImageData()
     sub_pre =  sub_preMat.GetOutput()
     # define image based on subtraction of postS -preS
     subtractedImage = sub_pre
     
     return subtractedImage