<gh_stars>0
import os
import unittest
from __main__ import vtk, qt, ctk, slicer
from slicer.ScriptedLoadableModule import *
import logging
import SimpleITK as sitk
import sitkUtils
import math
import numpy
import LabelStatistics

#
# ComputeT2Star
#

class ComputeT2Star(ScriptedLoadableModule):
  """Uses ScriptedLoadableModule base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def __init__(self, parent):
    ScriptedLoadableModule.__init__(self, parent)
    self.parent.title = "ComputeT2Star" # TODO make this more human readable by adding spaces
    self.parent.categories = ["IGT"]
    self.parent.dependencies = []
    self.parent.contributors = ["<NAME> (Brigham and Women's Hospital)"] # replace with "Firstname Lastname (Organization)"
    self.parent.helpText = """
    This is an example of scripted loadable module bundled in an extension.
    It performs a simple thresholding on the input volume and optionally captures a screenshot.
    """
    self.parent.acknowledgementText = """
    This module was developed based on a template created by <NAME>, Kitware Inc.
    and <NAME>, Isomics, Inc. and was partially funded by NIH grant 3P41RR013218-12S1.
""" # replace with organization, grant and thanks.

#
# ComputeT2StarWidget
#

class ComputeT2StarWidget(ScriptedLoadableModuleWidget):
  """Uses ScriptedLoadableModuleWidget base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setup(self):

    ScriptedLoadableModuleWidget.setup(self)

    #--------------------------------------------------
    # For debugging
    #
    # Reload and Test area
    reloadCollapsibleButton = ctk.ctkCollapsibleButton()
    reloadCollapsibleButton.text = "Reload && Test"
    self.layout.addWidget(reloadCollapsibleButton)
    reloadFormLayout = qt.QFormLayout(reloadCollapsibleButton)

    reloadCollapsibleButton.collapsed = True
    
    # reload button
    # (use this during development, but remove it when delivering
    #  your module to users)
    self.reloadButton = qt.QPushButton("Reload")
    self.reloadButton.toolTip = "Reload this module."
    self.reloadButton.name = "ComputeT2Star Reload"
    reloadFormLayout.addWidget(self.reloadButton)
    self.reloadButton.connect('clicked()', self.onReload)
    #
    #--------------------------------------------------

    # Instantiate and connect widgets ...

    #
    # Parameters Area
    #
    parametersCollapsibleButton = ctk.ctkCollapsibleButton()
    parametersCollapsibleButton.text = "Parameters"
    self.layout.addWidget(parametersCollapsibleButton)

    # Layout within the dummy collapsible button
    parametersFormLayout = qt.QFormLayout(parametersCollapsibleButton)

    #
    # input volume selector
    #
    self.inputTE1Selector = slicer.qMRMLNodeComboBox()
    self.inputTE1Selector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.inputTE1Selector.selectNodeUponCreation = True
    self.inputTE1Selector.addEnabled = True
    self.inputTE1Selector.removeEnabled = True
    self.inputTE1Selector.noneEnabled = True
    self.inputTE1Selector.renameEnabled = True
    self.inputTE1Selector.showHidden = False
    self.inputTE1Selector.showChildNodeTypes = False
    self.inputTE1Selector.setMRMLScene( slicer.mrmlScene )
    self.inputTE1Selector.setToolTip( "Pick the first volume" )
    parametersFormLayout.addRow("Input Volume 1: ", self.inputTE1Selector)

    #
    # input volume selector
    #
    self.inputTE2Selector = slicer.qMRMLNodeComboBox()
    self.inputTE2Selector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.inputTE2Selector.selectNodeUponCreation = True
    self.inputTE2Selector.addEnabled = True
    self.inputTE2Selector.removeEnabled = True
    self.inputTE2Selector.noneEnabled = True
    self.inputTE2Selector.renameEnabled = True
    self.inputTE2Selector.showHidden = False
    self.inputTE2Selector.showChildNodeTypes = False
    self.inputTE2Selector.setMRMLScene( slicer.mrmlScene )
    self.inputTE2Selector.setToolTip( "Pick the second volume" )
    parametersFormLayout.addRow("Input Volume 2: ", self.inputTE2Selector)


    #
    # reference ROI selector
    #
    self.referenceROISelector = slicer.qMRMLNodeComboBox()
    self.referenceROISelector.nodeTypes = ( ("vtkMRMLLabelMapVolumeNode"), "" )
    self.referenceROISelector.selectNodeUponCreation = False
    self.referenceROISelector.addEnabled = True
    self.referenceROISelector.removeEnabled = True
    self.referenceROISelector.noneEnabled = True
    self.referenceROISelector.renameEnabled = True
    self.referenceROISelector.showHidden = False
    self.referenceROISelector.showChildNodeTypes = False
    self.referenceROISelector.setMRMLScene( slicer.mrmlScene )
    self.referenceROISelector.setToolTip( "Reference ROI for scaling factor and noise estimation" )
    parametersFormLayout.addRow("Reference ROI: ", self.referenceROISelector)

    #
    # outputT2Star volume selector
    #
    self.outputT2StarSelector = slicer.qMRMLNodeComboBox()
    self.outputT2StarSelector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.outputT2StarSelector.selectNodeUponCreation = True
    self.outputT2StarSelector.addEnabled = True
    self.outputT2StarSelector.removeEnabled = True
    self.outputT2StarSelector.noneEnabled = True
    self.outputT2StarSelector.renameEnabled = True
    self.outputT2StarSelector.showHidden = False
    self.outputT2StarSelector.showChildNodeTypes = False
    self.outputT2StarSelector.setMRMLScene( slicer.mrmlScene )
    self.outputT2StarSelector.setToolTip( "Pick the T2Star output volume." )
    parametersFormLayout.addRow("T2Star Output Volume: ", self.outputT2StarSelector)

    #
    # outputR2Star volume selector
    #
    self.outputR2StarSelector = slicer.qMRMLNodeComboBox()
    self.outputR2StarSelector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.outputR2StarSelector.selectNodeUponCreation = True
    self.outputR2StarSelector.addEnabled = True
    self.outputR2StarSelector.removeEnabled = True
    self.outputR2StarSelector.noneEnabled = True
    self.outputR2StarSelector.renameEnabled = True
    self.outputR2StarSelector.showHidden = False
    self.outputR2StarSelector.showChildNodeTypes = False
    self.outputR2StarSelector.setMRMLScene( slicer.mrmlScene )
    self.outputR2StarSelector.setToolTip( "Pick the R2Star output volume." )
    parametersFormLayout.addRow("R2Star Output Volume: ", self.outputR2StarSelector)

    #
    # First TE
    #
    self.TE1SpinBox = qt.QDoubleSpinBox()
    self.TE1SpinBox.objectName = 'TE1SpinBox'
    self.TE1SpinBox.setMaximum(100.0)
    self.TE1SpinBox.setMinimum(0.0000)
    self.TE1SpinBox.setDecimals(8)
    self.TE1SpinBox.setValue(0.00007)
    self.TE1SpinBox.setToolTip("TE for Input Volume 1")
    parametersFormLayout.addRow("TE1 (s): ", self.TE1SpinBox)

    #
    # Second TE
    #
    self.TE2SpinBox = qt.QDoubleSpinBox()
    self.TE2SpinBox.objectName = 'TE2SpinBox'
    self.TE2SpinBox.setMaximum(100.0)
    self.TE2SpinBox.setMinimum(0.0000)
    self.TE2SpinBox.setDecimals(8)
    self.TE2SpinBox.setValue(0.002)
    self.TE2SpinBox.setToolTip("TE for Input Volume 2")
    parametersFormLayout.addRow("TE2 (s): ", self.TE2SpinBox)

    #
    # Scaling Factor
    #
    self.ScaleSpinBox = qt.QDoubleSpinBox()
    self.ScaleSpinBox.objectName = 'ScaleSpinBox'
    self.ScaleSpinBox.setMaximum(10.000)
    self.ScaleSpinBox.setMinimum(0.000)
    self.ScaleSpinBox.setDecimals(8)
    self.ScaleSpinBox.setValue(1.000)
    self.ScaleSpinBox.setToolTip("Scaling factor to adjust magnitude of volume 2.")
    parametersFormLayout.addRow("Scaling Factor: ", self.ScaleSpinBox)

    #
    # R2* for Scale calibration 
    #
    self.scaleCalibrationR2sSpinBox = qt.QDoubleSpinBox()
    self.scaleCalibrationR2sSpinBox.objectName = 'scaleCalibrationR2sSpinBox'
    self.scaleCalibrationR2sSpinBox.setMaximum(1000.0)
    self.scaleCalibrationR2sSpinBox.setMinimum(0.0)
    self.scaleCalibrationR2sSpinBox.setDecimals(8)
    self.scaleCalibrationR2sSpinBox.setValue(129.565)
    self.scaleCalibrationR2sSpinBox.setToolTip("Scale Calibration")
    parametersFormLayout.addRow("Scale Clibration R2* (s^-1): ", self.scaleCalibrationR2sSpinBox)

    #
    # Echo 1/2 signal lower input threshold
    #
    self.Echo1InputThresholdSpinBox = qt.QDoubleSpinBox()
    self.Echo1InputThresholdSpinBox.objectName = 'Echo1InputThresholdSpinBox'
    self.Echo1InputThresholdSpinBox.setMaximum(65536.000)
    self.Echo1InputThresholdSpinBox.setMinimum(0.000)
    self.Echo1InputThresholdSpinBox.setDecimals(6)
    self.Echo1InputThresholdSpinBox.setValue(0.000)
    self.Echo1InputThresholdSpinBox.setToolTip("Lower input threshold for echo 1.")
    parametersFormLayout.addRow("Lower Input Threshold (Echo 1): ", self.Echo1InputThresholdSpinBox)

    self.Echo2InputThresholdSpinBox = qt.QDoubleSpinBox()
    self.Echo2InputThresholdSpinBox.objectName = 'Echo2InputThresholdSpinBox'
    self.Echo2InputThresholdSpinBox.setMaximum(65536.000)
    self.Echo2InputThresholdSpinBox.setMinimum(0.000)
    self.Echo2InputThresholdSpinBox.setDecimals(6)
    self.Echo2InputThresholdSpinBox.setValue(0.000)
    self.Echo2InputThresholdSpinBox.setToolTip("Lower input threshold for echo 2.")
    parametersFormLayout.addRow("Lower Input Threshold (Echo 2): ", self.Echo2InputThresholdSpinBox)

    self.MinT2sSpinBox = qt.QDoubleSpinBox()
    self.MinT2sSpinBox.objectName = 'MinT2sSpinBox'
    self.MinT2sSpinBox.setMinimum(1.000)
    self.MinT2sSpinBox.setMinimum(0.000)
    self.MinT2sSpinBox.setDecimals(6)
    self.MinT2sSpinBox.setValue(0.00125)
    self.MinT2sSpinBox.setToolTip("Minimum T2* for output (maximum R2* = 1 / (minimum T2*)).")
    parametersFormLayout.addRow("Minimum T2* for output (s): ", self.MinT2sSpinBox)
    
    #
    # Check box to correct noise
    #
    self.useNoiseCorrectionFlagCheckBox = qt.QCheckBox()
    self.useNoiseCorrectionFlagCheckBox.checked = 1
    self.useNoiseCorrectionFlagCheckBox.setToolTip("If checked, correct noise based on the estimated noise level.")
    parametersFormLayout.addRow("Use Noise Correction", self.useNoiseCorrectionFlagCheckBox)

    #
    # Noise Level
    #
    self.Echo1NoiseSpinBox = qt.QDoubleSpinBox()
    self.Echo1NoiseSpinBox.objectName = 'Echo1NoiseSpinBox'
    self.Echo1NoiseSpinBox.setMaximum(500.0)
    self.Echo1NoiseSpinBox.setMinimum(0.0)
    self.Echo1NoiseSpinBox.setDecimals(6)
    self.Echo1NoiseSpinBox.setValue(0.0)
    self.Echo1NoiseSpinBox.setToolTip("Noise level for 1st echo noise correction.")
    parametersFormLayout.addRow("Noise Level (Echo 1): ", self.Echo1NoiseSpinBox)

    self.Echo2NoiseSpinBox = qt.QDoubleSpinBox()
    self.Echo2NoiseSpinBox.objectName = 'Echo2NoiseSpinBox'
    self.Echo2NoiseSpinBox.setMaximum(500.0)
    self.Echo2NoiseSpinBox.setMinimum(0.0)
    self.Echo2NoiseSpinBox.setDecimals(6)
    self.Echo2NoiseSpinBox.setValue(0.0)
    self.Echo2NoiseSpinBox.setToolTip("Noise level for 1st echo noise correction.")
    parametersFormLayout.addRow("Noise Level (Echo 2): ", self.Echo2NoiseSpinBox)

    #
    # check box to use threshold
    #
    self.useOutputThresholdFlagCheckBox = qt.QCheckBox()
    self.useOutputThresholdFlagCheckBox.checked = 1
    self.useOutputThresholdFlagCheckBox.setToolTip("If checked, apply the threshold to limit the pixel value ranges.")
    parametersFormLayout.addRow("Use OutputThreshold", self.useOutputThresholdFlagCheckBox)

    #
    # Upper threshold - We set threshold value to limit the range of intensity 
    #
    self.upperOutputThresholdSpinBox = qt.QDoubleSpinBox()
    self.upperOutputThresholdSpinBox.objectName = 'upperOutputThresholdSpinBox'
    self.upperOutputThresholdSpinBox.setMaximum(1000000.0)
    self.upperOutputThresholdSpinBox.setMinimum(-1000000.0)
    self.upperOutputThresholdSpinBox.setDecimals(6)
    self.upperOutputThresholdSpinBox.setValue(1000000.0)
    self.upperOutputThresholdSpinBox.setToolTip("Upper threshold for the output")
    parametersFormLayout.addRow("Upper OutputThreshold (s): ", self.upperOutputThresholdSpinBox)

    #
    # Lower threshold - We set threshold value to limit the range of intensity 
    #
    self.lowerOutputThresholdSpinBox = qt.QDoubleSpinBox()
    self.lowerOutputThresholdSpinBox.objectName = 'lowerOutputThresholdSpinBox'
    self.lowerOutputThresholdSpinBox.setMaximum(1000000.0)
    self.lowerOutputThresholdSpinBox.setMinimum(-1000000.0)
    self.lowerOutputThresholdSpinBox.setDecimals(6)
    self.lowerOutputThresholdSpinBox.setValue(-1000000.0)
    self.lowerOutputThresholdSpinBox.setToolTip("Lower threshold for the output")
    parametersFormLayout.addRow("Lower OutputThreshold (s): ", self.lowerOutputThresholdSpinBox)

    #
    # Apply Button
    #
    self.applyButton = qt.QPushButton("Apply")
    self.applyButton.toolTip = "Run the algorithm."
    self.applyButton.enabled = False
    parametersFormLayout.addRow(self.applyButton)

    # connections
    self.applyButton.connect('clicked(bool)', self.onApplyButton)
    self.inputTE1Selector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.inputTE2Selector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.referenceROISelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.outputT2StarSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.outputR2StarSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.useOutputThresholdFlagCheckBox.connect('toggled(bool)', self.onUseOutputThreshold)
    self.useNoiseCorrectionFlagCheckBox.connect('toggled(bool)', self.onUseNoiseCorrection)

    # Add vertical spacer
    self.layout.addStretch(1)

    # Refresh Apply button state
    self.onSelect()

  def cleanup(self):
    pass

  def onSelect(self):

    if self.referenceROISelector.currentNode():
      self.ScaleSpinBox.enabled = False
    else:
      self.ScaleSpinBox.enabled = True

    if self.useNoiseCorrectionFlagCheckBox.checked and self.referenceROISelector.currentNode() == None:
      self.Echo1NoiseSpinBox.enabled = True
      self.Echo2NoiseSpinBox.enabled = True
    else:
      self.Echo1NoiseSpinBox.enabled = False
      self.Echo2NoiseSpinBox.enabled = False

    self.applyButton.enabled = self.inputTE1Selector.currentNode() and self.inputTE1Selector.currentNode() and (self.outputT2StarSelector.currentNode() or self.outputR2StarSelector.currentNode())

  def onUseOutputThreshold(self):
    if self.useOutputThresholdFlagCheckBox.checked == True:
      self.lowerOutputThresholdSpinBox.enabled = True;      
      self.upperOutputThresholdSpinBox.enabled = True;      
    else:
      self.lowerOutputThresholdSpinBox.enabled = False;      
      self.upperOutputThresholdSpinBox.enabled = False;

  def onUseNoiseCorrection(self):
    if self.useNoiseCorrectionFlagCheckBox.checked == True and self.referenceROISelector.currentNode() == None:
      self.Echo1NoiseSpinBox.enabled = True;
      self.Echo2NoiseSpinBox.enabled = True;      
    else:
      self.Echo1NoiseSpinBox.enabled = False;
      self.Echo2NoiseSpinBox.enabled = False;      
    

  def onApplyButton(self):
    logic = ComputeT2StarLogic()
    #enableScreenshotsFlag = self.enableScreenshotsFlagCheckBox.checked
    #imageOutputThreshold = self.imageOutputThresholdSliderWidget.value
    t2name = ''
    r2name = ''
    if self.outputT2StarSelector.currentNode():
      t2name = self.outputT2StarSelector.currentNode().GetName()
    if self.outputR2StarSelector.currentNode():
      r2name = self.outputR2StarSelector.currentNode().GetName()
      
    inputThreshold = [self.Echo1InputThresholdSpinBox.value, self.Echo2InputThresholdSpinBox.value]
    
    minT2s = self.MinT2sSpinBox.value
    
    outputThreshold = None
    if self.useOutputThresholdFlagCheckBox.checked == True:
      outputThreshold = [self.lowerOutputThresholdSpinBox.value, self.upperOutputThresholdSpinBox.value]
      
    noiseLevel = None
    if self.useNoiseCorrectionFlagCheckBox.checked == True:
      noiseLevel = [self.Echo1NoiseSpinBox.value, self.Echo2NoiseSpinBox.value]

    scaleFactor = self.ScaleSpinBox.value
    noiseLevel = None

    if self.referenceROISelector.currentNode():
      inputTE1VolumeNode = self.inputTE1Selector.currentNode()
      inputTE2VolumeNode = self.inputTE2Selector.currentNode()
      ROINode = self.referenceROISelector.currentNode()

      imageTE1 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE1VolumeNode.GetID()), sitk.sitkFloat64)
      imageTE2 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE2VolumeNode.GetID()), sitk.sitkFloat64)
      roiImage = sitk.Cast(sitkUtils.PullFromSlicer(ROINode.GetID()), sitk.sitkInt8)
      
      noiseEcho1 = logic.CalcNoise(imageTE1, None, roiImage)
      noiseEcho2 = logic.CalcNoise(imageTE2, None, roiImage)
      noiseLevel = [noiseEcho1, noiseEcho2]
      print "noises = [%f, %f]\n" % (noiseEcho1, noiseEcho2)

      logic.CorrectNoise(imageTE1, noiseEcho1)
      logic.CorrectNoise(imageTE2, noiseEcho2)

      scaleFactor = logic.CalcScalingFactor(imageTE1, imageTE2, roiImage, 
                                            self.TE1SpinBox.value, self.TE2SpinBox.value,
                                            self.scaleCalibrationR2sSpinBox.value)

      print "scale = %f\n" % scaleFactor

      self.ScaleSpinBox.value = scaleFactor
      self.Echo1NoiseSpinBox.value = noiseEcho1
      self.Echo2NoiseSpinBox.value = noiseEcho2
    else:
      if self.useNoiseCorrectionFlagCheckBox.checked:
        noiseLevel = [self.Echo1NoiseSpinBox.value, self.Echo2NoiseSpinBox.value]
    
    logic.run(self.inputTE1Selector.currentNode(), self.inputTE2Selector.currentNode(),
              self.outputT2StarSelector.currentNode(), self.outputR2StarSelector.currentNode(),
              self.TE1SpinBox.value, self.TE2SpinBox.value, scaleFactor,
              noiseLevel, outputThreshold, inputThreshold, minT2s)

    ### Since PushToSlicer() called in logic.run() will delete the original node, obtain the new node and
    ### reset the selector.
    t2Node = slicer.util.getNode(t2name)
    r2Node = slicer.util.getNode(r2name)
    self.outputT2StarSelector.setCurrentNode(t2Node)
    self.outputR2StarSelector.setCurrentNode(r2Node)

  def onReload(self, moduleName="ComputeT2Star"):
    # Generic reload method for any scripted module.
    # ModuleWizard will subsitute correct default moduleName.

    globals()[moduleName] = slicer.util.reloadScriptedModule(moduleName)


#
# ComputeT2StarLogic
#

class ComputeT2StarLogic(ScriptedLoadableModuleLogic):

  def __init__(self):
    ScriptedLoadableModuleLogic.__init__(self)
  
  def isValidInputOutputData(self, inputTE1VolumeNode, inputTE2VolumeNode):
    """Validates if the output is not the same as input
    """
    if not inputTE1VolumeNode:
      logging.debug('isValidInputOutputData failed: no input volume node for TE1 image defined')
      return False
    if not inputTE2VolumeNode:
      logging.debug('isValidInputOutputData failed: no input volume node for TE2 image defined')
      return False
    return True

  
  def CalcNoise(self, image1, image2, roiImage):

    LabelStatistics = sitk.LabelStatisticsImageFilter()
      
    if image2:
      #Aimage1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32)
      #Aimage2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32)
      #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(roiImageNode.GetID()), sitk.sitkInt8)
      subImage = sitk.Subtract(image1, image2)
      absImage = sitk.Abs(subImage)
            
      LabelStatistics.Execute(absImage, roiImage)
      meanAbsDiff = LabelStatistics.GetMean(1)
      return (meanAbsDiff/math.sqrt(math.pi/2.0))

    else:
      LabelStatistics.Execute(image1, roiImage)
      SD = LabelStatistics.GetSigma(1)
      return SD

  def CorrectNoise(self, image, noiseLevel):

    squareImage = sitk.Pow(image, 2)
    subImage = sitk.Subtract(squareImage, noiseLevel*noiseLevel)
    subImagePositive1 = sitk.Threshold(subImage,0.0, float('Inf'), 0.0)
    image = sitk.Sqrt(subImagePositive1)


  def CalcScalingFactor(self, image1, image2, roiImage, TE1, TE2, scaleCalibrationR2s):

    #image1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32)
    #image2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32)
    #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(ROINode.GetID()), sitk.sitkInt8)

    LabelStatistics = sitk.LabelStatisticsImageFilter()
    LabelStatistics.Execute(image1, roiImage)
    echo1 = LabelStatistics.GetMean(1)
    LabelStatistics.Execute(image2, roiImage)
    echo2 = LabelStatistics.GetMean(1)

    scale = echo1 / (echo2 * numpy.exp(scaleCalibrationR2s*(TE2-TE1)))

    return (scale)


  def run(self, inputTE1VolumeNode, inputTE2VolumeNode, outputT2StarVolumeNode, outputR2StarVolumeNode, TE1, TE2, scaleFactor, noiseLevel, outputThreshold, inputThreshold, minT2s):
    """
    Run the actual algorithm
    """

    echo1NoiseLevel = 0.0
    echo2NoiseLevel = 0.0
    upperOutputThreshold = 0.0
    lowerOutputThreshold = 0.0

    if not self.isValidInputOutputData(inputTE1VolumeNode, inputTE2VolumeNode):
      slicer.util.errorDisplay('Input volume is the same as output volume. Choose a different output volume.')
      return False

    logging.info('Processing started')

    imageTE1 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE1VolumeNode.GetID()), sitk.sitkFloat64)
    imageTE2 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE2VolumeNode.GetID()), sitk.sitkFloat64)

    # Noise correction
    # Echo 1
    if noiseLevel != None:
      echo1NoiseLevel = noiseLevel[0]
      echo2NoiseLevel = noiseLevel[1]

      squareImage1 = sitk.Pow(imageTE1, 2)
      subImage1 = sitk.Subtract(squareImage1, echo1NoiseLevel*echo1NoiseLevel)
      subImagePositive1 = sitk.Threshold(subImage1,0.0, float('Inf'), 0.0)
      imageTE1 = sitk.Sqrt(subImagePositive1)

      squareImage2 = sitk.Pow(imageTE2, 2)
      subImage2 = sitk.Subtract(squareImage2, echo2NoiseLevel*echo2NoiseLevel)
      subImagePositive2 = sitk.Threshold(subImage2,0.0, float('Inf'), 0.0)
      imageTE2 = sitk.Sqrt(subImagePositive2)
    else:
      # Simply remove negative values (not needed?)
      imageTE1 = sitk.Threshold(imageTE1,0.0, float('Inf'), 0.0)
      imageTE2 = sitk.Threshold(imageTE2,0.0, float('Inf'), 0.0)

    ## Apply scaling factor to the second echo
    imageTE2 = sitk.Multiply(imageTE2, scaleFactor)

    ## Create mask to exclude invalid pixels
    #  Criteria for validation
    #   1. First or second echo signal is less than the input threshold
    #   2. Second echo signal is greater than the first
    
    mask = None
    imaskFillT2s = None
    imaskFillR2s = None

    if inputThreshold != None:
      mask1 = sitk.BinaryThreshold(imageTE1, inputThreshold[0], float('Inf'), 1, 0) 
      mask2 = sitk.BinaryThreshold(imageTE2, inputThreshold[1], float('Inf'), 1, 0) 
      #mask3 = sitk.Greater(imageTE2, imageTE1, 1, 0)
      mask = sitk.And(mask1, mask2)
      #mask = sitk.And(mask, mask3)

      imask = sitk.Not(mask)
      imaskFloat = sitk.Cast(imask, sitk.sitkFloat64)
      imaskFillT2s = imaskFloat * minT2s
      imaskFillR2s = 0.0

      if minT2s > 0:
        imaskFillR2s = imaskFloat * (1/minT2s)

    if outputThreshold != None:
      lowerOutputThreshold = outputThreshold[0]
      upperOutputThreshold = outputThreshold[1]


    if outputT2StarVolumeNode:
      imageT2Star = sitk.Divide(TE1-TE2, sitk.Log(sitk.Divide(imageTE2, imageTE1)))
      if inputThreshold != None:
        imageT2Star = sitk.Mask(imageT2Star, mask)
        imageT2Star = sitk.Add(imageT2Star, imaskFillT2s)
      if outputThreshold != None:
        imageT2StarThreshold = sitk.Threshold(imageT2Star, lowerOutputThreshold, upperOutputThreshold, 0.0)
        sitkUtils.PushToSlicer(imageT2StarThreshold, outputT2StarVolumeNode.GetName(), 0, True)
      else:
        sitkUtils.PushToSlicer(imageT2Star, outputT2StarVolumeNode.GetName(), 0, True)

    if outputR2StarVolumeNode:
      imageR2Star = sitk.Divide(sitk.Log(sitk.Divide(imageTE2, imageTE1)), TE1-TE2)
      if inputThreshold != None:
        imageR2Star = sitk.Mask(imageR2Star, mask)
        imageR2Star = sitk.Add(imageR2Star, imaskFillR2s)
      if outputThreshold != None:
        imageR2StarThreshold = sitk.Threshold(imageR2Star, lowerOutputThreshold, upperOutputThreshold, 0.0)
        sitkUtils.PushToSlicer(imageR2StarThreshold, outputR2StarVolumeNode.GetName(), 0, True)
      else:
        sitkUtils.PushToSlicer(imageR2Star, outputR2StarVolumeNode.GetName(), 0, True)

    logging.info('Processing completed')

    return True


class ComputeT2StarTest(ScriptedLoadableModuleTest):
  """
  This is the test case for your scripted module.
  Uses ScriptedLoadableModuleTest base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setUp(self):
    """ Do whatever is needed to reset the state - typically a scene clear will be enough.
    """
    slicer.mrmlScene.Clear(0)

  def runTest(self):
    """Run as few or as many tests as needed here.
    """
    self.setUp()
    self.test_ComputeT2Star1()

  def test_ComputeT2Star1(self):
    """ Ideally you should have several levels of tests.  At the lowest level
    tests should exercise the functionality of the logic with different inputs
    (both valid and invalid).  At higher levels your tests should emulate the
    way the user would interact with your code and confirm that it still works
    the way you intended.
    One of the most important features of the tests is that it should alert other
    developers when their changes will have an impact on the behavior of your
    module.  For example, if a developer removes a feature that you depend on,
    your test should break so they know that the feature is needed.
    """

    pass

    #self.delayDisplay("Starting the test")
    ##
    ## first, get some data
    ##
    #import urllib
    #downloads = (
    #    ('http://slicer.kitware.com/midas3/download?items=5767', 'FA.nrrd', slicer.util.loadVolume),
    #    )
    #
    #for url,name,loader in downloads:
    #  filePath = slicer.app.temporaryPath + '/' + name
    #  if not os.path.exists(filePath) or os.stat(filePath).st_size == 0:
    #    logging.info('Requesting download %s from %s...\n' % (name, url))
    #    urllib.urlretrieve(url, filePath)
    #  if loader:
    #    logging.info('Loading %s...' % (name,))
    #    loader(filePath)
    #self.delayDisplay('Finished with download and loading')
    #
    #volumeNode = slicer.util.getNode(pattern="FA")
    #logic = ComputeT2StarLogic()
    #self.assertTrue( logic.hasImageData(volumeNode) )
    #self.delayDisplay('Test passed!')

import os
import unittest
from __main__ import vtk, qt, ctk, slicer
from slicer.ScriptedLoadableModule import *
import logging
import SimpleITK as sitk
import sitkUtils
import math
import numpy
import LabelStatistics

#
# ComputeT2Star
#

class ComputeT2Star(ScriptedLoadableModule):
  """Uses ScriptedLoadableModule base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def __init__(self, parent):
    ScriptedLoadableModule.__init__(self, parent)
    self.parent.title = "ComputeT2Star" # TODO make this more human readable by adding spaces
    self.parent.categories = ["IGT"]
    self.parent.dependencies = []
    self.parent.contributors = ["<NAME> (Brigham and Women's Hospital)"] # replace with "Firstname Lastname (Organization)"
    self.parent.helpText = """
    This is an example of scripted loadable module bundled in an extension.
    It performs a simple thresholding on the input volume and optionally captures a screenshot.
    """
    self.parent.acknowledgementText = """
    This module was developed based on a template created by <NAME>, Kitware Inc.
    and <NAME>, Isomics, Inc. and was partially funded by NIH grant 3P41RR013218-12S1.
""" # replace with organization, grant and thanks.

#
# ComputeT2StarWidget
#

class ComputeT2StarWidget(ScriptedLoadableModuleWidget):
  """Uses ScriptedLoadableModuleWidget base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setup(self):

    ScriptedLoadableModuleWidget.setup(self)

    #--------------------------------------------------
    # For debugging
    #
    # Reload and Test area
    reloadCollapsibleButton = ctk.ctkCollapsibleButton()
    reloadCollapsibleButton.text = "Reload && Test"
    self.layout.addWidget(reloadCollapsibleButton)
    reloadFormLayout = qt.QFormLayout(reloadCollapsibleButton)

    reloadCollapsibleButton.collapsed = True
    
    # reload button
    # (use this during development, but remove it when delivering
    #  your module to users)
    self.reloadButton = qt.QPushButton("Reload")
    self.reloadButton.toolTip = "Reload this module."
    self.reloadButton.name = "ComputeT2Star Reload"
    reloadFormLayout.addWidget(self.reloadButton)
    self.reloadButton.connect('clicked()', self.onReload)
    #
    #--------------------------------------------------

    # Instantiate and connect widgets ...

    #
    # Parameters Area
    #
    parametersCollapsibleButton = ctk.ctkCollapsibleButton()
    parametersCollapsibleButton.text = "Parameters"
    self.layout.addWidget(parametersCollapsibleButton)

    # Layout within the dummy collapsible button
    parametersFormLayout = qt.QFormLayout(parametersCollapsibleButton)

    #
    # input volume selector
    #
    self.inputTE1Selector = slicer.qMRMLNodeComboBox()
    self.inputTE1Selector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.inputTE1Selector.selectNodeUponCreation = True
    self.inputTE1Selector.addEnabled = True
    self.inputTE1Selector.removeEnabled = True
    self.inputTE1Selector.noneEnabled = True
    self.inputTE1Selector.renameEnabled = True
    self.inputTE1Selector.showHidden = False
    self.inputTE1Selector.showChildNodeTypes = False
    self.inputTE1Selector.setMRMLScene( slicer.mrmlScene )
    self.inputTE1Selector.setToolTip( "Pick the first volume" )
    parametersFormLayout.addRow("Input Volume 1: ", self.inputTE1Selector)

    #
    # input volume selector
    #
    self.inputTE2Selector = slicer.qMRMLNodeComboBox()
    self.inputTE2Selector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.inputTE2Selector.selectNodeUponCreation = True
    self.inputTE2Selector.addEnabled = True
    self.inputTE2Selector.removeEnabled = True
    self.inputTE2Selector.noneEnabled = True
    self.inputTE2Selector.renameEnabled = True
    self.inputTE2Selector.showHidden = False
    self.inputTE2Selector.showChildNodeTypes = False
    self.inputTE2Selector.setMRMLScene( slicer.mrmlScene )
    self.inputTE2Selector.setToolTip( "Pick the second volume" )
    parametersFormLayout.addRow("Input Volume 2: ", self.inputTE2Selector)


    #
    # reference ROI selector
    #
    self.referenceROISelector = slicer.qMRMLNodeComboBox()
    self.referenceROISelector.nodeTypes = ( ("vtkMRMLLabelMapVolumeNode"), "" )
    self.referenceROISelector.selectNodeUponCreation = False
    self.referenceROISelector.addEnabled = True
    self.referenceROISelector.removeEnabled = True
    self.referenceROISelector.noneEnabled = True
    self.referenceROISelector.renameEnabled = True
    self.referenceROISelector.showHidden = False
    self.referenceROISelector.showChildNodeTypes = False
    self.referenceROISelector.setMRMLScene( slicer.mrmlScene )
    self.referenceROISelector.setToolTip( "Reference ROI for scaling factor and noise estimation" )
    parametersFormLayout.addRow("Reference ROI: ", self.referenceROISelector)

    #
    # outputT2Star volume selector
    #
    self.outputT2StarSelector = slicer.qMRMLNodeComboBox()
    self.outputT2StarSelector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.outputT2StarSelector.selectNodeUponCreation = True
    self.outputT2StarSelector.addEnabled = True
    self.outputT2StarSelector.removeEnabled = True
    self.outputT2StarSelector.noneEnabled = True
    self.outputT2StarSelector.renameEnabled = True
    self.outputT2StarSelector.showHidden = False
    self.outputT2StarSelector.showChildNodeTypes = False
    self.outputT2StarSelector.setMRMLScene( slicer.mrmlScene )
    self.outputT2StarSelector.setToolTip( "Pick the T2Star output volume." )
    parametersFormLayout.addRow("T2Star Output Volume: ", self.outputT2StarSelector)

    #
    # outputR2Star volume selector
    #
    self.outputR2StarSelector = slicer.qMRMLNodeComboBox()
    self.outputR2StarSelector.nodeTypes = ( ("vtkMRMLScalarVolumeNode"), "" )
    self.outputR2StarSelector.selectNodeUponCreation = True
    self.outputR2StarSelector.addEnabled = True
    self.outputR2StarSelector.removeEnabled = True
    self.outputR2StarSelector.noneEnabled = True
    self.outputR2StarSelector.renameEnabled = True
    self.outputR2StarSelector.showHidden = False
    self.outputR2StarSelector.showChildNodeTypes = False
    self.outputR2StarSelector.setMRMLScene( slicer.mrmlScene )
    self.outputR2StarSelector.setToolTip( "Pick the R2Star output volume." )
    parametersFormLayout.addRow("R2Star Output Volume: ", self.outputR2StarSelector)

    #
    # First TE
    #
    self.TE1SpinBox = qt.QDoubleSpinBox()
    self.TE1SpinBox.objectName = 'TE1SpinBox'
    self.TE1SpinBox.setMaximum(100.0)
    self.TE1SpinBox.setMinimum(0.0000)
    self.TE1SpinBox.setDecimals(8)
    self.TE1SpinBox.setValue(0.00007)
    self.TE1SpinBox.setToolTip("TE for Input Volume 1")
    parametersFormLayout.addRow("TE1 (s): ", self.TE1SpinBox)

    #
    # Second TE
    #
    self.TE2SpinBox = qt.QDoubleSpinBox()
    self.TE2SpinBox.objectName = 'TE2SpinBox'
    self.TE2SpinBox.setMaximum(100.0)
    self.TE2SpinBox.setMinimum(0.0000)
    self.TE2SpinBox.setDecimals(8)
    self.TE2SpinBox.setValue(0.002)
    self.TE2SpinBox.setToolTip("TE for Input Volume 2")
    parametersFormLayout.addRow("TE2 (s): ", self.TE2SpinBox)

    #
    # Scaling Factor
    #
    self.ScaleSpinBox = qt.QDoubleSpinBox()
    self.ScaleSpinBox.objectName = 'ScaleSpinBox'
    self.ScaleSpinBox.setMaximum(10.000)
    self.ScaleSpinBox.setMinimum(0.000)
    self.ScaleSpinBox.setDecimals(8)
    self.ScaleSpinBox.setValue(1.000)
    self.ScaleSpinBox.setToolTip("Scaling factor to adjust magnitude of volume 2.")
    parametersFormLayout.addRow("Scaling Factor: ", self.ScaleSpinBox)

    #
    # R2* for Scale calibration 
    #
    self.scaleCalibrationR2sSpinBox = qt.QDoubleSpinBox()
    self.scaleCalibrationR2sSpinBox.objectName = 'scaleCalibrationR2sSpinBox'
    self.scaleCalibrationR2sSpinBox.setMaximum(1000.0)
    self.scaleCalibrationR2sSpinBox.setMinimum(0.0)
    self.scaleCalibrationR2sSpinBox.setDecimals(8)
    self.scaleCalibrationR2sSpinBox.setValue(129.565)
    self.scaleCalibrationR2sSpinBox.setToolTip("Scale Calibration")
    parametersFormLayout.addRow("Scale Clibration R2* (s^-1): ", self.scaleCalibrationR2sSpinBox)

    #
    # Echo 1/2 signal lower input threshold
    #
    self.Echo1InputThresholdSpinBox = qt.QDoubleSpinBox()
    self.Echo1InputThresholdSpinBox.objectName = 'Echo1InputThresholdSpinBox'
    self.Echo1InputThresholdSpinBox.setMaximum(65536.000)
    self.Echo1InputThresholdSpinBox.setMinimum(0.000)
    self.Echo1InputThresholdSpinBox.setDecimals(6)
    self.Echo1InputThresholdSpinBox.setValue(0.000)
    self.Echo1InputThresholdSpinBox.setToolTip("Lower input threshold for echo 1.")
    parametersFormLayout.addRow("Lower Input Threshold (Echo 1): ", self.Echo1InputThresholdSpinBox)

    self.Echo2InputThresholdSpinBox = qt.QDoubleSpinBox()
    self.Echo2InputThresholdSpinBox.objectName = 'Echo2InputThresholdSpinBox'
    self.Echo2InputThresholdSpinBox.setMaximum(65536.000)
    self.Echo2InputThresholdSpinBox.setMinimum(0.000)
    self.Echo2InputThresholdSpinBox.setDecimals(6)
    self.Echo2InputThresholdSpinBox.setValue(0.000)
    self.Echo2InputThresholdSpinBox.setToolTip("Lower input threshold for echo 2.")
    parametersFormLayout.addRow("Lower Input Threshold (Echo 2): ", self.Echo2InputThresholdSpinBox)

    self.MinT2sSpinBox = qt.QDoubleSpinBox()
    self.MinT2sSpinBox.objectName = 'MinT2sSpinBox'
    self.MinT2sSpinBox.setMinimum(1.000)
    self.MinT2sSpinBox.setMinimum(0.000)
    self.MinT2sSpinBox.setDecimals(6)
    self.MinT2sSpinBox.setValue(0.00125)
    self.MinT2sSpinBox.setToolTip("Minimum T2* for output (maximum R2* = 1 / (minimum T2*)).")
    parametersFormLayout.addRow("Minimum T2* for output (s): ", self.MinT2sSpinBox)
    
    #
    # Check box to correct noise
    #
    self.useNoiseCorrectionFlagCheckBox = qt.QCheckBox()
    self.useNoiseCorrectionFlagCheckBox.checked = 1
    self.useNoiseCorrectionFlagCheckBox.setToolTip("If checked, correct noise based on the estimated noise level.")
    parametersFormLayout.addRow("Use Noise Correction", self.useNoiseCorrectionFlagCheckBox)

    #
    # Noise Level
    #
    self.Echo1NoiseSpinBox = qt.QDoubleSpinBox()
    self.Echo1NoiseSpinBox.objectName = 'Echo1NoiseSpinBox'
    self.Echo1NoiseSpinBox.setMaximum(500.0)
    self.Echo1NoiseSpinBox.setMinimum(0.0)
    self.Echo1NoiseSpinBox.setDecimals(6)
    self.Echo1NoiseSpinBox.setValue(0.0)
    self.Echo1NoiseSpinBox.setToolTip("Noise level for 1st echo noise correction.")
    parametersFormLayout.addRow("Noise Level (Echo 1): ", self.Echo1NoiseSpinBox)

    self.Echo2NoiseSpinBox = qt.QDoubleSpinBox()
    self.Echo2NoiseSpinBox.objectName = 'Echo2NoiseSpinBox'
    self.Echo2NoiseSpinBox.setMaximum(500.0)
    self.Echo2NoiseSpinBox.setMinimum(0.0)
    self.Echo2NoiseSpinBox.setDecimals(6)
    self.Echo2NoiseSpinBox.setValue(0.0)
    self.Echo2NoiseSpinBox.setToolTip("Noise level for 1st echo noise correction.")
    parametersFormLayout.addRow("Noise Level (Echo 2): ", self.Echo2NoiseSpinBox)

    #
    # check box to use threshold
    #
    self.useOutputThresholdFlagCheckBox = qt.QCheckBox()
    self.useOutputThresholdFlagCheckBox.checked = 1
    self.useOutputThresholdFlagCheckBox.setToolTip("If checked, apply the threshold to limit the pixel value ranges.")
    parametersFormLayout.addRow("Use OutputThreshold", self.useOutputThresholdFlagCheckBox)

    #
    # Upper threshold - We set threshold value to limit the range of intensity 
    #
    self.upperOutputThresholdSpinBox = qt.QDoubleSpinBox()
    self.upperOutputThresholdSpinBox.objectName = 'upperOutputThresholdSpinBox'
    self.upperOutputThresholdSpinBox.setMaximum(1000000.0)
    self.upperOutputThresholdSpinBox.setMinimum(-1000000.0)
    self.upperOutputThresholdSpinBox.setDecimals(6)
    self.upperOutputThresholdSpinBox.setValue(1000000.0)
    self.upperOutputThresholdSpinBox.setToolTip("Upper threshold for the output")
    parametersFormLayout.addRow("Upper OutputThreshold (s): ", self.upperOutputThresholdSpinBox)

    #
    # Lower threshold - We set threshold value to limit the range of intensity 
    #
    self.lowerOutputThresholdSpinBox = qt.QDoubleSpinBox()
    self.lowerOutputThresholdSpinBox.objectName = 'lowerOutputThresholdSpinBox'
    self.lowerOutputThresholdSpinBox.setMaximum(1000000.0)
    self.lowerOutputThresholdSpinBox.setMinimum(-1000000.0)
    self.lowerOutputThresholdSpinBox.setDecimals(6)
    self.lowerOutputThresholdSpinBox.setValue(-1000000.0)
    self.lowerOutputThresholdSpinBox.setToolTip("Lower threshold for the output")
    parametersFormLayout.addRow("Lower OutputThreshold (s): ", self.lowerOutputThresholdSpinBox)

    #
    # Apply Button
    #
    self.applyButton = qt.QPushButton("Apply")
    self.applyButton.toolTip = "Run the algorithm."
    self.applyButton.enabled = False
    parametersFormLayout.addRow(self.applyButton)

    # connections
    self.applyButton.connect('clicked(bool)', self.onApplyButton)
    self.inputTE1Selector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.inputTE2Selector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.referenceROISelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.outputT2StarSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.outputR2StarSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
    self.useOutputThresholdFlagCheckBox.connect('toggled(bool)', self.onUseOutputThreshold)
    self.useNoiseCorrectionFlagCheckBox.connect('toggled(bool)', self.onUseNoiseCorrection)

    # Add vertical spacer
    self.layout.addStretch(1)

    # Refresh Apply button state
    self.onSelect()

  def cleanup(self):
    pass

  def onSelect(self):

    if self.referenceROISelector.currentNode():
      self.ScaleSpinBox.enabled = False
    else:
      self.ScaleSpinBox.enabled = True

    if self.useNoiseCorrectionFlagCheckBox.checked and self.referenceROISelector.currentNode() == None:
      self.Echo1NoiseSpinBox.enabled = True
      self.Echo2NoiseSpinBox.enabled = True
    else:
      self.Echo1NoiseSpinBox.enabled = False
      self.Echo2NoiseSpinBox.enabled = False

    self.applyButton.enabled = self.inputTE1Selector.currentNode() and self.inputTE1Selector.currentNode() and (self.outputT2StarSelector.currentNode() or self.outputR2StarSelector.currentNode())

  def onUseOutputThreshold(self):
    if self.useOutputThresholdFlagCheckBox.checked == True:
      self.lowerOutputThresholdSpinBox.enabled = True;      
      self.upperOutputThresholdSpinBox.enabled = True;      
    else:
      self.lowerOutputThresholdSpinBox.enabled = False;      
      self.upperOutputThresholdSpinBox.enabled = False;

  def onUseNoiseCorrection(self):
    if self.useNoiseCorrectionFlagCheckBox.checked == True and self.referenceROISelector.currentNode() == None:
      self.Echo1NoiseSpinBox.enabled = True;
      self.Echo2NoiseSpinBox.enabled = True;      
    else:
      self.Echo1NoiseSpinBox.enabled = False;
      self.Echo2NoiseSpinBox.enabled = False;      
    

  def onApplyButton(self):
    logic = ComputeT2StarLogic()
    #enableScreenshotsFlag = self.enableScreenshotsFlagCheckBox.checked
    #imageOutputThreshold = self.imageOutputThresholdSliderWidget.value
    t2name = ''
    r2name = ''
    if self.outputT2StarSelector.currentNode():
      t2name = self.outputT2StarSelector.currentNode().GetName()
    if self.outputR2StarSelector.currentNode():
      r2name = self.outputR2StarSelector.currentNode().GetName()
      
    inputThreshold = [self.Echo1InputThresholdSpinBox.value, self.Echo2InputThresholdSpinBox.value]
    
    minT2s = self.MinT2sSpinBox.value
    
    outputThreshold = None
    if self.useOutputThresholdFlagCheckBox.checked == True:
      outputThreshold = [self.lowerOutputThresholdSpinBox.value, self.upperOutputThresholdSpinBox.value]
      
    noiseLevel = None
    if self.useNoiseCorrectionFlagCheckBox.checked == True:
      noiseLevel = [self.Echo1NoiseSpinBox.value, self.Echo2NoiseSpinBox.value]

    scaleFactor = self.ScaleSpinBox.value
    noiseLevel = None

    if self.referenceROISelector.currentNode():
      inputTE1VolumeNode = self.inputTE1Selector.currentNode()
      inputTE2VolumeNode = self.inputTE2Selector.currentNode()
      ROINode = self.referenceROISelector.currentNode()

      imageTE1 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE1VolumeNode.GetID()), sitk.sitkFloat64)
      imageTE2 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE2VolumeNode.GetID()), sitk.sitkFloat64)
      roiImage = sitk.Cast(sitkUtils.PullFromSlicer(ROINode.GetID()), sitk.sitkInt8)
      
      noiseEcho1 = logic.CalcNoise(imageTE1, None, roiImage)
      noiseEcho2 = logic.CalcNoise(imageTE2, None, roiImage)
      noiseLevel = [noiseEcho1, noiseEcho2]
      print "noises = [%f, %f]\n" % (noiseEcho1, noiseEcho2)

      logic.CorrectNoise(imageTE1, noiseEcho1)
      logic.CorrectNoise(imageTE2, noiseEcho2)

      scaleFactor = logic.CalcScalingFactor(imageTE1, imageTE2, roiImage, 
                                            self.TE1SpinBox.value, self.TE2SpinBox.value,
                                            self.scaleCalibrationR2sSpinBox.value)

      print "scale = %f\n" % scaleFactor

      self.ScaleSpinBox.value = scaleFactor
      self.Echo1NoiseSpinBox.value = noiseEcho1
      self.Echo2NoiseSpinBox.value = noiseEcho2
    else:
      if self.useNoiseCorrectionFlagCheckBox.checked:
        noiseLevel = [self.Echo1NoiseSpinBox.value, self.Echo2NoiseSpinBox.value]
    
    logic.run(self.inputTE1Selector.currentNode(), self.inputTE2Selector.currentNode(),
              self.outputT2StarSelector.currentNode(), self.outputR2StarSelector.currentNode(),
              self.TE1SpinBox.value, self.TE2SpinBox.value, scaleFactor,
              noiseLevel, outputThreshold, inputThreshold, minT2s)

    ### Since PushToSlicer() called in logic.run() will delete the original node, obtain the new node and
    ### reset the selector.
    t2Node = slicer.util.getNode(t2name)
    r2Node = slicer.util.getNode(r2name)
    self.outputT2StarSelector.setCurrentNode(t2Node)
    self.outputR2StarSelector.setCurrentNode(r2Node)

  def onReload(self, moduleName="ComputeT2Star"):
    # Generic reload method for any scripted module.
    # ModuleWizard will subsitute correct default moduleName.

    globals()[moduleName] = slicer.util.reloadScriptedModule(moduleName)


#
# ComputeT2StarLogic
#

class ComputeT2StarLogic(ScriptedLoadableModuleLogic):

  def __init__(self):
    ScriptedLoadableModuleLogic.__init__(self)
  
  def isValidInputOutputData(self, inputTE1VolumeNode, inputTE2VolumeNode):
    """Validates if the output is not the same as input
    """
    if not inputTE1VolumeNode:
      logging.debug('isValidInputOutputData failed: no input volume node for TE1 image defined')
      return False
    if not inputTE2VolumeNode:
      logging.debug('isValidInputOutputData failed: no input volume node for TE2 image defined')
      return False
    return True

  
  def CalcNoise(self, image1, image2, roiImage):

    LabelStatistics = sitk.LabelStatisticsImageFilter()
      
    if image2:
      #Aimage1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32)
      #Aimage2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32)
      #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(roiImageNode.GetID()), sitk.sitkInt8)
      subImage = sitk.Subtract(image1, image2)
      absImage = sitk.Abs(subImage)
            
      LabelStatistics.Execute(absImage, roiImage)
      meanAbsDiff = LabelStatistics.GetMean(1)
      return (meanAbsDiff/math.sqrt(math.pi/2.0))

    else:
      LabelStatistics.Execute(image1, roiImage)
      SD = LabelStatistics.GetSigma(1)
      return SD

  def CorrectNoise(self, image, noiseLevel):

    squareImage = sitk.Pow(image, 2)
    subImage = sitk.Subtract(squareImage, noiseLevel*noiseLevel)
    subImagePositive1 = sitk.Threshold(subImage,0.0, float('Inf'), 0.0)
    image = sitk.Sqrt(subImagePositive1)


  def CalcScalingFactor(self, image1, image2, roiImage, TE1, TE2, scaleCalibrationR2s):

    #image1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32)
    #image2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32)
    #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(ROINode.GetID()), sitk.sitkInt8)

    LabelStatistics = sitk.LabelStatisticsImageFilter()
    LabelStatistics.Execute(image1, roiImage)
    echo1 = LabelStatistics.GetMean(1)
    LabelStatistics.Execute(image2, roiImage)
    echo2 = LabelStatistics.GetMean(1)

    scale = echo1 / (echo2 * numpy.exp(scaleCalibrationR2s*(TE2-TE1)))

    return (scale)


  def run(self, inputTE1VolumeNode, inputTE2VolumeNode, outputT2StarVolumeNode, outputR2StarVolumeNode, TE1, TE2, scaleFactor, noiseLevel, outputThreshold, inputThreshold, minT2s):
    """
    Run the actual algorithm
    """

    echo1NoiseLevel = 0.0
    echo2NoiseLevel = 0.0
    upperOutputThreshold = 0.0
    lowerOutputThreshold = 0.0

    if not self.isValidInputOutputData(inputTE1VolumeNode, inputTE2VolumeNode):
      slicer.util.errorDisplay('Input volume is the same as output volume. Choose a different output volume.')
      return False

    logging.info('Processing started')

    imageTE1 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE1VolumeNode.GetID()), sitk.sitkFloat64)
    imageTE2 = sitk.Cast(sitkUtils.PullFromSlicer(inputTE2VolumeNode.GetID()), sitk.sitkFloat64)

    # Noise correction
    # Echo 1
    if noiseLevel != None:
      echo1NoiseLevel = noiseLevel[0]
      echo2NoiseLevel = noiseLevel[1]

      squareImage1 = sitk.Pow(imageTE1, 2)
      subImage1 = sitk.Subtract(squareImage1, echo1NoiseLevel*echo1NoiseLevel)
      subImagePositive1 = sitk.Threshold(subImage1,0.0, float('Inf'), 0.0)
      imageTE1 = sitk.Sqrt(subImagePositive1)

      squareImage2 = sitk.Pow(imageTE2, 2)
      subImage2 = sitk.Subtract(squareImage2, echo2NoiseLevel*echo2NoiseLevel)
      subImagePositive2 = sitk.Threshold(subImage2,0.0, float('Inf'), 0.0)
      imageTE2 = sitk.Sqrt(subImagePositive2)
    else:
      # Simply remove negative values (not needed?)
      imageTE1 = sitk.Threshold(imageTE1,0.0, float('Inf'), 0.0)
      imageTE2 = sitk.Threshold(imageTE2,0.0, float('Inf'), 0.0)

    ## Apply scaling factor to the second echo
    imageTE2 = sitk.Multiply(imageTE2, scaleFactor)

    ## Create mask to exclude invalid pixels
    #  Criteria for validation
    #   1. First or second echo signal is less than the input threshold
    #   2. Second echo signal is greater than the first
    
    mask = None
    imaskFillT2s = None
    imaskFillR2s = None

    if inputThreshold != None:
      mask1 = sitk.BinaryThreshold(imageTE1, inputThreshold[0], float('Inf'), 1, 0) 
      mask2 = sitk.BinaryThreshold(imageTE2, inputThreshold[1], float('Inf'), 1, 0) 
      #mask3 = sitk.Greater(imageTE2, imageTE1, 1, 0)
      mask = sitk.And(mask1, mask2)
      #mask = sitk.And(mask, mask3)

      imask = sitk.Not(mask)
      imaskFloat = sitk.Cast(imask, sitk.sitkFloat64)
      imaskFillT2s = imaskFloat * minT2s
      imaskFillR2s = 0.0

      if minT2s > 0:
        imaskFillR2s = imaskFloat * (1/minT2s)

    if outputThreshold != None:
      lowerOutputThreshold = outputThreshold[0]
      upperOutputThreshold = outputThreshold[1]


    if outputT2StarVolumeNode:
      imageT2Star = sitk.Divide(TE1-TE2, sitk.Log(sitk.Divide(imageTE2, imageTE1)))
      if inputThreshold != None:
        imageT2Star = sitk.Mask(imageT2Star, mask)
        imageT2Star = sitk.Add(imageT2Star, imaskFillT2s)
      if outputThreshold != None:
        imageT2StarThreshold = sitk.Threshold(imageT2Star, lowerOutputThreshold, upperOutputThreshold, 0.0)
        sitkUtils.PushToSlicer(imageT2StarThreshold, outputT2StarVolumeNode.GetName(), 0, True)
      else:
        sitkUtils.PushToSlicer(imageT2Star, outputT2StarVolumeNode.GetName(), 0, True)

    if outputR2StarVolumeNode:
      imageR2Star = sitk.Divide(sitk.Log(sitk.Divide(imageTE2, imageTE1)), TE1-TE2)
      if inputThreshold != None:
        imageR2Star = sitk.Mask(imageR2Star, mask)
        imageR2Star = sitk.Add(imageR2Star, imaskFillR2s)
      if outputThreshold != None:
        imageR2StarThreshold = sitk.Threshold(imageR2Star, lowerOutputThreshold, upperOutputThreshold, 0.0)
        sitkUtils.PushToSlicer(imageR2StarThreshold, outputR2StarVolumeNode.GetName(), 0, True)
      else:
        sitkUtils.PushToSlicer(imageR2Star, outputR2StarVolumeNode.GetName(), 0, True)

    logging.info('Processing completed')

    return True


class ComputeT2StarTest(ScriptedLoadableModuleTest):
  """
  This is the test case for your scripted module.
  Uses ScriptedLoadableModuleTest base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setUp(self):
    """ Do whatever is needed to reset the state - typically a scene clear will be enough.
    """
    slicer.mrmlScene.Clear(0)

  def runTest(self):
    """Run as few or as many tests as needed here.
    """
    self.setUp()
    self.test_ComputeT2Star1()

  def test_ComputeT2Star1(self):
    """ Ideally you should have several levels of tests.  At the lowest level
    tests should exercise the functionality of the logic with different inputs
    (both valid and invalid).  At higher levels your tests should emulate the
    way the user would interact with your code and confirm that it still works
    the way you intended.
    One of the most important features of the tests is that it should alert other
    developers when their changes will have an impact on the behavior of your
    module.  For example, if a developer removes a feature that you depend on,
    your test should break so they know that the feature is needed.
    """

    pass

    #self.delayDisplay("Starting the test")
    ##
    ## first, get some data
    ##
    #import urllib
    #downloads = (
    #    ('http://slicer.kitware.com/midas3/download?items=5767', 'FA.nrrd', slicer.util.loadVolume),
    #    )
    #
    #for url,name,loader in downloads:
    #  filePath = slicer.app.temporaryPath + '/' + name
    #  if not os.path.exists(filePath) or os.stat(filePath).st_size == 0:
    #    logging.info('Requesting download %s from %s...\n' % (name, url))
    #    urllib.urlretrieve(url, filePath)
    #  if loader:
    #    logging.info('Loading %s...' % (name,))
    #    loader(filePath)
    #self.delayDisplay('Finished with download and loading')
    #
    #volumeNode = slicer.util.getNode(pattern="FA")
    #logic = ComputeT2StarLogic()
    #self.assertTrue( logic.hasImageData(volumeNode) )
    #self.delayDisplay('Test passed!')
# # ComputeT2Star # Uses ScriptedLoadableModule base class, available at:   https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py # TODO make this more human readable by adding spaces # replace with "Firstname Lastname (Organization)" This is an example of scripted loadable module bundled in an extension.     It performs a simple thresholding on the input volume and optionally captures a screenshot. This module was developed based on a template created by <NAME>, Kitware Inc.     and <NAME>, Isomics, Inc. and was partially funded by NIH grant 3P41RR013218-12S1. # replace with organization, grant and thanks. # # ComputeT2StarWidget # Uses ScriptedLoadableModuleWidget base class, available at:   https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py #-------------------------------------------------- # For debugging # # Reload and Test area # reload button # (use this during development, but remove it when delivering #  your module to users) # #-------------------------------------------------- # Instantiate and connect widgets ... # # Parameters Area # # Layout within the dummy collapsible button # # input volume selector # # # input volume selector # # # reference ROI selector # # # outputT2Star volume selector # # # outputR2Star volume selector # # # First TE # # # Second TE # # # Scaling Factor # # # R2* for Scale calibration # # # Echo 1/2 signal lower input threshold # # # Check box to correct noise # # # Noise Level # # # check box to use threshold # # # Upper threshold - We set threshold value to limit the range of intensity # # # Lower threshold - We set threshold value to limit the range of intensity # # # Apply Button # # connections # Add vertical spacer # Refresh Apply button state #enableScreenshotsFlag = self.enableScreenshotsFlagCheckBox.checked #imageOutputThreshold = self.imageOutputThresholdSliderWidget.value ### Since PushToSlicer() called in logic.run() will delete the original node, obtain the new node and ### reset the selector. # Generic reload method for any scripted module. # ModuleWizard will subsitute correct default moduleName. # # ComputeT2StarLogic # Validates if the output is not the same as input #Aimage1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32) #Aimage2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32) #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(roiImageNode.GetID()), sitk.sitkInt8) #image1 = sitk.Cast(sitkUtils.PullFromSlicer(image1Node.GetID()), sitk.sitkFloat32) #image2 = sitk.Cast(sitkUtils.PullFromSlicer(image2Node.GetID()), sitk.sitkFloat32) #roiImage = sitk.Cast(sitkUtils.PullFromSlicer(ROINode.GetID()), sitk.sitkInt8) Run the actual algorithm # Noise correction # Echo 1 # Simply remove negative values (not needed?) ## Apply scaling factor to the second echo ## Create mask to exclude invalid pixels #  Criteria for validation #   1. First or second echo signal is less than the input threshold #   2. Second echo signal is greater than the first #mask3 = sitk.Greater(imageTE2, imageTE1, 1, 0) #mask = sitk.And(mask, mask3) This is the test case for your scripted module.   Uses ScriptedLoadableModuleTest base class, available at:   https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py Do whatever is needed to reset the state - typically a scene clear will be enough. Run as few or as many tests as needed here. Ideally you should have several levels of tests.  At the lowest level     tests should exercise the functionality of the logic with different inputs     (both valid and invalid).  At higher levels your tests should emulate the     way the user would interact with your code and confirm that it still works     the way you intended.     One of the most important features of the tests is that it should alert other     developers when their changes will have an impact on the behavior of your     module.  For example, if a developer removes a feature that you depend on,     your test should break so they know that the feature is needed. #self.delayDisplay("Starting the test") ## ## first, get some data ## #import urllib #downloads = ( #    ('http://slicer.kitware.com/midas3/download?items=5767', 'FA.nrrd', slicer.util.loadVolume), #    ) # #for url,name,loader in downloads: #  filePath = slicer.app.temporaryPath + '/' + name #  if not os.path.exists(filePath) or os.stat(filePath).st_size == 0: #    logging.info('Requesting download %s from %s...\n' % (name, url)) #    urllib.urlretrieve(url, filePath) #  if loader: #    logging.info('Loading %s...' % (name,)) #    loader(filePath) #self.delayDisplay('Finished with download and loading') # #volumeNode = slicer.util.getNode(pattern="FA") #logic = ComputeT2StarLogic() #self.assertTrue( logic.hasImageData(volumeNode) ) #self.delayDisplay('Test passed!')