def get_histogram(image, maxval=0, minval=0, maxrange=0):
    """
	Return the histogram of the image as a list of floats
	"""
    accu = vtk.vtkImageAccumulate()
    accu.SetInputConnection(image.GetProducerPort())

    if maxval == 0:
        x0, x1 = getImageScalarRange(image)
        # x1 = int(math.floor(x1))
        # x0 = int(math.ceil(x0))
    else:
        # x0, x1 = (int(math.ceil(minval)),int(math.floor(maxval)))
        x0, x1 = (minval, maxval)

        # accu.SetComponentExtent(0, x1 - x0, 0, 0, 0, 0)
        # accu.SetComponentSpacing(1, 0, 0)
    if maxrange:
        accu.SetComponentExtent(0, x1 - x0, 0, 0, 0, 0)
        accu.SetComponentSpacing(1, 0, 0)
    else:
        accu.SetComponentExtent(0, 255, 0, 0, 0, 0)
        accu.SetComponentSpacing((x1 - x0 + 1) / 256.0, 0, 0)

    accu.SetComponentOrigin(x0, 0, 0)
    accu.Update()
    data = accu.GetOutput()

    values = []
    x0, x1, y0, y1, z0, z1 = data.GetWholeExtent()

    for i in range(x0, x1 + 1):
        c = data.GetScalarComponentAsDouble(i, 0, 0, 0)
        values.append(c)
    return values

  def TestSection_1_AddRemoveSegment(self):
    # Add/remove segment from segmentation (check display properties, color table, etc.)
    logging.info('Test section 1: Add/remove segment')

    # Get baseline values
    displayNode = self.inputSegmentationNode.GetDisplayNode()
    self.assertIsNotNone(displayNode)
    colorTableNode = displayNode.GetColorNode()
    self.assertIsNotNone(colorTableNode)
    self.assertEqual(colorTableNode.GetNumberOfColors(), 4)
    # If segments are not found then the returned color is the pre-defined invalid color
    bodyColor = displayNode.GetSegmentColor('Body_Contour')
    self.assertTrue(int(bodyColor[0]*100) == 33 and int(bodyColor[1]*100) == 66 and bodyColor[2] == 0.0)
    tumorColor = displayNode.GetSegmentColor('Tumor_Contour')
    self.assertTrue(tumorColor[0] == 1.0 and tumorColor[1] == 0.0 and tumorColor[2] == 0.0)

    # Create new segment
    sphere = vtk.vtkSphereSource()
    sphere.SetCenter(0,50,0)
    sphere.SetRadius(50)
    sphere.Update()
    spherePolyData = vtk.vtkPolyData()
    spherePolyData.DeepCopy(sphere.GetOutput())
    
    self.sphereSegment = vtkSegmentationCore.vtkSegment()
    self.sphereSegment.SetName(self.sphereSegmentName)
    self.sphereSegment.SetDefaultColor(0.0,0.0,1.0)
    self.sphereSegment.AddRepresentation(self.closedSurfaceReprName, spherePolyData)

    # Add segment to segmentation
    self.inputSegmentationNode.GetSegmentation().AddSegment(self.sphereSegment)
    self.assertEqual(self.inputSegmentationNode.GetSegmentation().GetNumberOfSegments(), 3)
    self.assertEqual(colorTableNode.GetNumberOfColors(), 5)
    sphereColor = displayNode.GetSegmentColor(self.sphereSegmentName)
    self.assertTrue(sphereColor[0] == 0.0 and sphereColor[1] == 0.0 and sphereColor[2] == 1.0)
    
    # Check merged labelmap
    imageStat = vtk.vtkImageAccumulate()
    imageStat.SetInputData(self.inputSegmentationNode.GetImageData())
    imageStat.SetComponentExtent(0,4,0,0,0,0)
    imageStat.SetComponentOrigin(0,0,0)
    imageStat.SetComponentSpacing(1,1,1)
    imageStat.Update()
    self.assertEqual(imageStat.GetVoxelCount(), 1000)
    imageStatResult = imageStat.GetOutput()
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(0,0,0,0), 814)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(1,0,0,0), 0)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(2,0,0,0), 175)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(3,0,0,0), 4)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(4,0,0,0), 7)

    # Remove segment from segmentation
    self.inputSegmentationNode.GetSegmentation().RemoveSegment(self.sphereSegmentName)
    self.assertEqual(self.inputSegmentationNode.GetSegmentation().GetNumberOfSegments(), 2)
    self.assertEqual(colorTableNode.GetNumberOfColors(), 5)
    sphereColorArray = [0]*4
    colorTableNode.GetColor(4,sphereColorArray)
    self.assertTrue(int(sphereColorArray[0]*100) == 50 and int(sphereColorArray[1]*100) == 50 and int(sphereColorArray[2]*100) == 50)
    sphereColor = displayNode.GetSegmentColor(self.sphereSegmentName)
    self.assertTrue(sphereColor[0] == 0.5 and sphereColor[1] == 0.5 and sphereColor[2] == 0.5)

    def __init__(self, binsCount, reader):
        self.Reader = reader
        self.BinsCount = binsCount
        self.BinStep = 0
        self.ScalarRange = [0, 0]
        self.FreqArray = vtk.vtkDataArray.CreateDataArray(vtk.VTK_INT)
        self.FreqArray.SetNumberOfComponents(1)

        extract = vtk.vtkImageExtractComponents()
        extract.SetInputConnection(self.Reader.GetOutputPort())
        extract.SetComponents(0)
        extract.Update()
        extract.GetOutput().GetScalarRange(self.ScalarRange)

        self.BinStep = (self.ScalarRange[1] - self.ScalarRange[0]) / self.BinsCount

        histo = vtk.vtkImageAccumulate()
        histo.SetInputConnection(extract.GetOutputPort())
        histo.SetComponentExtent(0, self.BinsCount, 0, 0, 0, 0)
        histo.SetComponentOrigin(self.ScalarRange[0], 0, 0)
        histo.SetComponentSpacing(self.BinStep,  0, 0)
        histo.SetIgnoreZero(False)
        histo.Update()

        for j in range(0, self.BinsCount):
            compValue = histo.GetOutput().GetPointData().GetScalars().GetValue(j)
            self.FreqArray.InsertNextTuple1(compValue)

        lastBinValue = self.FreqArray.GetComponent(self.BinsCount - 1, 0)
        binValue = histo.GetOutput().GetPointData().GetScalars().GetValue(self.BinsCount)
        self.FreqArray.SetComponent(self.BinsCount - 1, 0, lastBinValue + binValue)

  def splitPerStructureVolumes(masterNode, mergeNode):
    """Make a separate label map node for each non-empty label value in the
    merged label map"""

    colorNode = mergeNode.GetDisplayNode().GetColorNode()

    accum = vtk.vtkImageAccumulate()
    accum.SetInputConnection(mergeNode.GetImageDataConnection())
    accum.Update()
    lo = int(accum.GetMin()[0])
    hi = int(accum.GetMax()[0])

    thresholder = vtk.vtkImageThreshold()
    for index in xrange(lo,hi+1):
      logging.info( "Splitting label %d..."%index )
      thresholder.SetInputConnection( mergeNode.GetImageDataConnection() )
      thresholder.SetInValue( index )
      thresholder.SetOutValue( 0 )
      thresholder.ReplaceInOn()
      thresholder.ReplaceOutOn()
      thresholder.ThresholdBetween( index, index )
      thresholder.SetOutputScalarType( mergeNode.GetImageData().GetScalarType() )
      thresholder.Update()
      if thresholder.GetOutput().GetScalarRange() != (0.0, 0.0):
        structureName = colorNode.GetColorName(index)
        logging.info( "Creating structure volume %s..."%structureName )
        structureVolume = EditUtil.structureVolume( masterNode, structureName )
        if not structureVolume:
          EditUtil.addStructure( masterNode, mergeNode, index )
        structureVolume = EditUtil.structureVolume( masterNode, structureName )
        structureVolume.GetImageData().DeepCopy( thresholder.GetOutput() )
        EditUtil.markVolumeNodeAsModified(structureVolume)

  def computeStatistics(self, segmentID):
    import vtkSegmentationCorePython as vtkSegmentationCore
    requestedKeys = self.getRequestedKeys()

    segmentationNode = slicer.mrmlScene.GetNodeByID(self.getParameterNode().GetParameter("Segmentation"))

    if len(requestedKeys)==0:
      return {}

    containsLabelmapRepresentation = segmentationNode.GetSegmentation().ContainsRepresentation(
      vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName())
    if not containsLabelmapRepresentation:
      return {}

    segment = segmentationNode.GetSegmentation().GetSegment(segmentID)
    segBinaryLabelName = vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName()
    segmentLabelmap = segment.GetRepresentation(segBinaryLabelName)

    if (not segmentLabelmap
      or not segmentLabelmap.GetPointData()
      or not segmentLabelmap.GetPointData().GetScalars()):
      # No input label data
      return {}

    # We need to know exactly the value of the segment voxels, apply threshold to make force the selected label value
    labelValue = 1
    backgroundValue = 0
    thresh = vtk.vtkImageThreshold()
    thresh.SetInputData(segmentLabelmap)
    thresh.ThresholdByLower(0)
    thresh.SetInValue(backgroundValue)
    thresh.SetOutValue(labelValue)
    thresh.SetOutputScalarType(vtk.VTK_UNSIGNED_CHAR)
    thresh.Update()

    #  Use binary labelmap as a stencil
    stencil = vtk.vtkImageToImageStencil()
    stencil.SetInputData(thresh.GetOutput())
    stencil.ThresholdByUpper(labelValue)
    stencil.Update()

    stat = vtk.vtkImageAccumulate()
    stat.SetInputData(thresh.GetOutput())
    stat.SetStencilData(stencil.GetOutput())
    stat.Update()

    # Add data to statistics list
    cubicMMPerVoxel = reduce(lambda x,y: x*y, segmentLabelmap.GetSpacing())
    ccPerCubicMM = 0.001
    stats = {}
    if "voxel_count" in requestedKeys:
      stats["voxel_count"] = stat.GetVoxelCount()
    if "volume_mm3" in requestedKeys:
      stats["volume_mm3"] = stat.GetVoxelCount() * cubicMMPerVoxel
    if "volume_cm3" in requestedKeys:
      stats["volume_cm3"] = stat.GetVoxelCount() * cubicMMPerVoxel * ccPerCubicMM
    return stats

 def getForegroundVoxelCount(self, imageData):
   if imageData is None:
     logging.error('Invalid input image data')
     return False
   imageAccumulate = vtk.vtkImageAccumulate()
   imageAccumulate.SetInputData(imageData)
   imageAccumulate.SetIgnoreZero(1)
   imageAccumulate.Update()
   return imageAccumulate.GetVoxelCount()

  def TestSection_MergeLabelmapWithDifferentGeometries(self):
    # Merge labelmap when segments containing labelmaps with different geometries (both same directions, different directions)
    logging.info('Test section: Merge labelmap with different geometries')

    self.assertIsNotNone(self.sphereSegment)
    self.sphereSegment.RemoveRepresentation(self.binaryLabelmapReprName)
    self.assertIsNone(self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName))

    # Create new segmentation with sphere segment
    self.secondSegmentationNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLSegmentationNode', 'Second')
    self.secondSegmentationNode.GetSegmentation().SetMasterRepresentationName(self.binaryLabelmapReprName)

    self.secondSegmentationNode.GetSegmentation().AddSegment(self.sphereSegment)

    # Check automatically converted labelmap. It is supposed to have the default geometry
    # (which is different than the one in the input segmentation)
    sphereLabelmap = self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName)
    self.assertIsNotNone(sphereLabelmap)
    sphereLabelmapSpacing = sphereLabelmap.GetSpacing()
    self.assertAlmostEqual(sphereLabelmapSpacing[0], 0.629257364931788, 8)
    self.assertAlmostEqual(sphereLabelmapSpacing[1], 0.629257364931788, 8)
    self.assertAlmostEqual(sphereLabelmapSpacing[2], 0.629257364931788, 8)

    # Create binary labelmap in segmentation that will create the merged labelmap from
    # different geometries so that labelmap is not removed from sphere segment when adding
    self.inputSegmentationNode.GetSegmentation().CreateRepresentation(self.binaryLabelmapReprName)

    # Copy segment to input segmentation
    self.inputSegmentationNode.GetSegmentation().CopySegmentFromSegmentation(self.secondSegmentationNode.GetSegmentation(), self.sphereSegmentName)
    self.assertEqual(self.inputSegmentationNode.GetSegmentation().GetNumberOfSegments(), 3)

    # Check merged labelmap
    # Reference geometry has the tiny patient spacing, and it is oversampled to have similar
    # voxel size as the sphere labelmap with the uniform 0.629mm spacing
    mergedLabelmap = vtkSegmentationCore.vtkOrientedImageData()
    self.inputSegmentationNode.GenerateMergedLabelmapForAllSegments(mergedLabelmap, 0)
    mergedLabelmapSpacing = mergedLabelmap.GetSpacing()
    self.assertAlmostEqual(mergedLabelmapSpacing[0], 0.80327868852459, 8)
    self.assertAlmostEqual(mergedLabelmapSpacing[1], 0.80327868852459, 8)
    self.assertAlmostEqual(mergedLabelmapSpacing[2], 0.377049180327869, 8)

    imageStat = vtk.vtkImageAccumulate()
    imageStat.SetInputData(mergedLabelmap)
    imageStat.SetComponentExtent(0,5,0,0,0,0)
    imageStat.SetComponentOrigin(0,0,0)
    imageStat.SetComponentSpacing(1,1,1)
    imageStat.Update()
    imageStatResult = imageStat.GetOutput()
    for i in range(5):
      logging.info("Volume {0}: {1}".format(i, imageStatResult.GetScalarComponentAsDouble(i,0,0,0)))
    self.assertEqual(imageStat.GetVoxelCount(), 226981000)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(0,0,0,0), 178838889)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(1,0,0,0), 39705288)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(2,0,0,0), 890883)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(3,0,0,0), 7545940)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(4,0,0,0), 0)  # Built from color table and color four is removed in previous test section

  def TestSection_2_MergeLabelmapWithDifferentGeometries(self):
    # Merge labelmap when segments containing labelmaps with different geometries (both same directions, different directions)
    logging.info('Test section 2: Merge labelmap with different geometries')

    self.assertIsNotNone(self.sphereSegment)
    self.sphereSegment.RemoveRepresentation(self.binaryLabelmapReprName)
    self.assertIsNone(self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName))

    # Create new segmentation with sphere segment
    self.secondSegmentationNode = slicer.vtkMRMLSegmentationNode()
    self.secondSegmentationNode.SetName('Second')
    self.secondSegmentationNode.GetSegmentation().SetMasterRepresentationName(self.binaryLabelmapReprName)
    slicer.mrmlScene.AddNode(self.secondSegmentationNode)

    self.secondSegmentationNode.GetSegmentation().AddSegment(self.sphereSegment)

    # Check automatically converted labelmap. It is supposed to have the default geometry
    # (which is different than the one in the input segmentation)
    sphereLabelmap = self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName)
    self.assertIsNotNone(sphereLabelmap)
    sphereLabelmapSpacing = sphereLabelmap.GetSpacing()
    self.assertTrue(sphereLabelmapSpacing[0] == 1.0 and sphereLabelmapSpacing[1] == 1.0 and sphereLabelmapSpacing[2] == 1.0)

    # Create binary labelmap in segmentation that will create the merged labelmap from
    # different geometries so that labelmap is not removed from sphere segment when adding
    self.inputSegmentationNode.GetSegmentation().CreateRepresentation(self.binaryLabelmapReprName)

    # Copy segment to input segmentation
    self.inputSegmentationNode.GetSegmentation().CopySegmentFromSegmentation(self.secondSegmentationNode.GetSegmentation(), self.sphereSegmentName)
    self.assertEqual(self.inputSegmentationNode.GetSegmentation().GetNumberOfSegments(), 3)

    # Check merged labelmap
    # Reference geometry has the tiny patient spacing, and it is oversampled to have smimilar
    # voxel size as the sphere labelmap with the uniform 1mm spacing
    mergedLabelmap = vtkSegmentationCore.vtkOrientedImageData()
    self.inputSegmentationNode.GenerateMergedLabelmapForAllSegments(mergedLabelmap, 0)
    mergedLabelmapSpacing = mergedLabelmap.GetSpacing()
    self.assertAlmostEqual(mergedLabelmapSpacing[0], 1.2894736842, 8)
    self.assertAlmostEqual(mergedLabelmapSpacing[1], 1.2894736842, 8)
    self.assertAlmostEqual(mergedLabelmapSpacing[2], 0.6052631578, 8)

    imageStat = vtk.vtkImageAccumulate()
    imageStat.SetInputData(mergedLabelmap)
    imageStat.SetComponentExtent(0,5,0,0,0,0)
    imageStat.SetComponentOrigin(0,0,0)
    imageStat.SetComponentSpacing(1,1,1)
    imageStat.Update()
    self.assertEqual(imageStat.GetVoxelCount(), 54872000)
    imageStatResult = imageStat.GetOutput()
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(0,0,0,0), 43573723)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(1,0,0,0), 10601312)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(2,0,0,0), 274360)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(3,0,0,0), 422605)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(4,0,0,0), 0)  # Built from color table and color four is removed in previous test section

 def imageDataContainsData(self, imageData):
   if imageData is None:
     logging.error('Invalid input image data')
     return False
   acc = vtk.vtkImageAccumulate()
   acc.SetInputData(imageData)
   acc.SetIgnoreZero(1)
   acc.Update()
   if acc.GetVoxelCount() > 0:
     return True
   return False

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self,
         module_manager,
         vtk.vtkImageAccumulate(),
         "Processing.",
         ("vtkImageData", "vtkImageStencilData"),
         ("vtkImageData",),
         replaceDoc=True,
         inputFunctions=None,
         outputFunctions=None,
     )

 def ResetWindowLevel(self):
     if self.IsLocked():
         return
     
     if self.Image:
         return
     
     self.Image.UpdateInformation()
     self.Image.SetUpdateExtent(self.Image.GetWholeExtent())
     self.Image.Update()
     
     if (self.Image.GetScalarType() == vtk.VTK_UNSIGNED_CHAR and 
         (self.Image.GetNumberOfScalarComponents() == 3 or 
          self.Image.GetNumberOfScalarComponents() == 4 )):
         return
     range = self.Image.GetScalarRange()
     histogram = vtk.vtkImageAccumulate()
     histogram.SetInput(self.Image)
     
     histogram.SetComponentExtent(0, 1000, 0, 0, 0, 0)
     histogram.SetComponentSpacing((range[1]-range[0])/1000.0, 0.0, 0.0)
     histogram.SetComponentOrigin(range[0], 0.0, 0.0)
     histogram.Update()
     
     output = histogram.GetOutput()
     ptData = output.GetPointData().GetScalars()
     if ptData:
         raise RuntimeError, "Error: Cannot cast point data to integers."
     numVox = histogram.GetVoxelCount()
     onePercent = numVox/100.0
     
     start = 1
     currentPercent = 0.0
     while (currentPercent<0.1 and start<999):
         ptData.GetTuple(start, tuple)
         currentPercent = currentPercent + tuple/onePercent
         start = start+1
         
     currentPercent = 0.0
     end = 999
     while (currentPercent<0.1 and end>0):
         ptData.GetTuple(start, tuple)
         currentPercent = currentPercent + tuple/onePercent
         end = end-1
     
     window = (end-start)*(range[1]-range[0])/1000.0
     level = 0.5*(start + end)*(range[1]-range[0])/1000.0
     window = (window-self.Shift)/self.Scale
     level = (level-self.Shift)/self.Scale
     
     self.SetWindow(window)
     self.SetLevel(level)
     del histogram

 def CalculateHistogram(self):
     proj = prj.Project()
     image = proj.imagedata
     r = int(image.GetScalarRange()[1] - image.GetScalarRange()[0])
     accumulate = vtk.vtkImageAccumulate()
     accumulate.SetInput(image)
     accumulate.SetComponentExtent(0, r -1, 0, 0, 0, 0)
     accumulate.SetComponentOrigin(image.GetScalarRange()[0], 0, 0)
     accumulate.Update()
     n_image = numpy_support.vtk_to_numpy(accumulate.GetOutput().GetPointData().GetScalars())
     ps.Publisher().sendMessage('Load histogram', (n_image,
                                                  image.GetScalarRange()))

    def test_seedconnect(self):
        """Test whether we can load and run a full network, select a point and 
        do a region growing. This broke with the introduction of vtk 5.6.1 due
        to more strict casting.
        """
                
        # load our little test network #####
        self._ge._load_and_realise_network(
            os.path.join(self._devide_testing.get_networks_dir(), 
                         'seedconnect.dvn'))
        
        # run the network once    
        self._ge._handler_execute_network(None)
        
        self._ge.canvas.redraw()
            
        # now find the slice3dVWR #####
        mm = self._devide_app.get_module_manager()
        svmod = mm.get_instance("svmod")
        # let's show the control frame
        svmod._handlerShowControls(None)
        
        if True:
            # we're doing this the long way to test more code paths
            svmod.sliceDirections.setCurrentCursor([20.0, 20.0, 20.0, 1.0])
            # this handler should result in the whole network being auto-executed
            # but somehow it blocks execution (the vktImageSeedConnect sticks at 0.0)
            svmod.selectedPoints._handlerStoreCursorAsPoint(None)
            
        else:        
            # it seems to block here as well: the whole network is linked up,
            # so it tries to execute when the storeCursor is called, and that
            # blocks everything. WHY?!
            #svmod.selectedPoints._storeCursor((20.0,20.0,20.0,1.0))
            #self.failUnless(len(svmod.selectedPoints._pointsList) == 1)

            # execute the network
            self._ge._handler_execute_network(None)

        # now count the number of voxels in the segmented result
        import vtk
        via = vtk.vtkImageAccumulate()
        scmod = mm.get_instance("scmod")
        via.SetInput(scmod.get_output(0))
        via.Update()
        # get second bin of output histogram: that should be the
        # number of voxels
        s = via.GetOutput().GetPointData().GetScalars()
        print s.GetTuple1(1)
        self.failUnless(s.GetTuple1(1) == 26728)
        via.SetInput(None)
        del via

 def getLabelsFromLabelMap(self, labelMapNode):
   if not labelMapNode:
     return
   accum = vtk.vtkImageAccumulate()
   accum.SetInput(labelMapNode.GetImageData())
   accum.UpdateWholeExtent()
   data = accum.GetOutput()
   data.Update()
   numBins = accum.GetComponentExtent()[1]
   nonZeroLabels = []
   for i in range(0, numBins + 1):
     numVoxels = data.GetScalarComponentAsDouble(i,0,0,0)
     if (numVoxels != 0):
       nonZeroLabels.append(i)
   return nonZeroLabels

  def split(self):
    """split the merge volume into individual structures"""

    self.statusText( "Splitting..." )
    merge = self.merge
    if not merge:
      return
    colorNode = merge.GetDisplayNode().GetColorNode()

    accum = vtk.vtkImageAccumulate()
    if vtk.VTK_MAJOR_VERSION <= 5:
      accum.SetInput(merge.GetImageData())
    else:
      accum.SetInputConnection(merge.GetImageDataConnection())
    accum.Update()
    lo = int(accum.GetMin()[0])
    hi = int(accum.GetMax()[0])

    # TODO: pending resolution of bug 1822, run the thresholding
    # in single threaded mode to avoid data corruption observed on mac release
    # builds
    thresholder = vtk.vtkImageThreshold()
    thresholder.SetNumberOfThreads(1)
    for i in xrange(lo,hi+1):
      self.statusText( "Splitting label %d..."%i )
      if vtk.VTK_MAJOR_VERSION <= 5:
        thresholder.SetInput( merge.GetImageData() )
      else:
        thresholder.SetInputConnection( merge.GetImageDataConnection() )
      thresholder.SetInValue( i )
      thresholder.SetOutValue( 0 )
      thresholder.ReplaceInOn()
      thresholder.ReplaceOutOn()
      thresholder.ThresholdBetween( i, i )
      thresholder.SetOutputScalarType( merge.GetImageData().GetScalarType() )
      thresholder.Update()
      if thresholder.GetOutput().GetScalarRange() != (0.0, 0.0):
        labelName = colorNode.GetColorName(i)
        self.statusText( "Creating structure volume %s..."%labelName )
        structureVolume = self.structureVolume( labelName )
        if not structureVolume:
          self.addStructure( i, "noEdit" )
        structureVolume = self.structureVolume( labelName )
        structureVolume.GetImageData().DeepCopy( thresholder.GetOutput() )
        EditUtil.markVolumeNodeAsModified(structureVolume)

    self.statusText( "Finished splitting." )

  def TestSection_2_MergeLabelmapWithDifferentGeometries(self):
    # Merge labelmap when segments containing labelmaps with different geometries (both same directions, different directions)
    logging.info('Test section 2: Merge labelmap with different geometries')

    self.assertIsNotNone(self.sphereSegment)
    self.sphereSegment.RemoveRepresentation(self.binaryLabelmapReprName)
    self.assertIsNone(self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName))

    # Create new segmentation with sphere segment
    self.secondSegmentationNode = vtkMRMLSegmentationNode()
    self.secondSegmentationNode.SetName('Second')
    self.secondSegmentationNode.GetSegmentation().SetMasterRepresentationName(self.binaryLabelmapReprName)
    slicer.mrmlScene.AddNode(self.secondSegmentationNode)

    self.secondSegmentationNode.GetSegmentation().AddSegment(self.sphereSegment)

    # Check automatically converted labelmap. It is supposed to have the default geometry
    # (which is different than the one in the input segmentation)
    sphereLabelmap = self.sphereSegment.GetRepresentation(self.binaryLabelmapReprName)
    self.assertIsNotNone(sphereLabelmap)
    sphereLabelmapSpacing = sphereLabelmap.GetSpacing()
    self.assertTrue(sphereLabelmapSpacing[0] == 1.0 and sphereLabelmapSpacing[1] == 1.0 and sphereLabelmapSpacing[2] == 1.0)

    # Copy segment to input segmentation
    self.inputSegmentationNode.GetSegmentation().CopySegmentFromSegmentation(self.secondSegmentationNode.GetSegmentation(), self.sphereSegmentName)
    self.assertEqual(self.inputSegmentationNode.GetSegmentation().GetNumberOfSegments(), 3)

    # Check merged labelmap
    mergedLabelmap = self.inputSegmentationNode.GetImageData()
    self.assertIsNotNone(mergedLabelmap)
    mergedLabelmapSpacing = sphereLabelmap.GetSpacing()
    self.assertTrue(mergedLabelmapSpacing[0] == 1.0 and mergedLabelmapSpacing[1] == 1.0 and mergedLabelmapSpacing[2] == 1.0)

    imageStat = vtk.vtkImageAccumulate()
    imageStat.SetInputData(mergedLabelmap)
    imageStat.SetComponentExtent(0,5,0,0,0,0)
    imageStat.SetComponentOrigin(0,0,0)
    imageStat.SetComponentSpacing(1,1,1)
    imageStat.Update()
    self.assertEqual(imageStat.GetVoxelCount(), 54872000)
    imageStatResult = imageStat.GetOutput()
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(0,0,0,0), 46678738)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(1,0,0,0), 0)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(2,0,0,0), 7618805)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(3,0,0,0), 128968)
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(4,0,0,0), 0) # Built from color table and color four is removed in previous test section
    self.assertEqual(imageStatResult.GetScalarComponentAsDouble(5,0,0,0), 445489)

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

        self._imageAccumulate = vtk.vtkImageAccumulate()
        module_utils.setup_vtk_object_progress(self, self._imageAccumulate, "Calculating histogram")

        self._viewFrame = None
        self._createViewFrame()
        self._bindEvents()

        self._config.numberOfBins = 256
        self._config.minValue = 0.0
        self._config.maxValue = 256.0

        self.config_to_logic()
        self.logic_to_config()
        self.config_to_view()

        self.view()

  def compareVolumes(self, volume1, volume2):
    subtractFilter = vtk.vtkImageMathematics()
    subtractFilter.SetOperationToSubtract()
    subtractFilter.SetInput1Data(volume1.GetImageData())
    subtractFilter.SetInput2Data(volume2.GetImageData())
    subtractFilter.Update()

    histogramFilter = vtk.vtkImageAccumulate()
    histogramFilter.SetComponentSpacing([1,0,0])
    histogramFilter.SetInputData(subtractFilter.GetOutput())
    histogramFilter.Update()

    # Both should be zero if generated scanlines == groundtruth
    maxValue = histogramFilter.GetMax()
    minValue = histogramFilter.GetMin()

    if (maxValue[0] == 0 and minValue[0] == 0):
      return True
    else:
      return False

  def addSegmentLabelmapStatistics(self):
    import vtkSegmentationCorePython as vtkSegmentationCore

    containsLabelmapRepresentation = self.segmentationNode.GetSegmentation().ContainsRepresentation(
      vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName())
    if not containsLabelmapRepresentation:
      return

    for segmentID in self.statistics["SegmentIDs"]:
      segment = self.segmentationNode.GetSegmentation().GetSegment(segmentID)
      segmentLabelmap = segment.GetRepresentation(vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationBinaryLabelmapRepresentationName())

      # We need to know exactly the value of the segment voxels, apply threshold to make force the selected label value
      labelValue = 1
      backgroundValue = 0
      thresh = vtk.vtkImageThreshold()
      thresh.SetInputData(segmentLabelmap)
      thresh.ThresholdByLower(0)
      thresh.SetInValue(backgroundValue)
      thresh.SetOutValue(labelValue)
      thresh.SetOutputScalarType(vtk.VTK_UNSIGNED_CHAR)
      thresh.Update()

      #  Use binary labelmap as a stencil
      stencil = vtk.vtkImageToImageStencil()
      stencil.SetInputData(thresh.GetOutput())
      stencil.ThresholdByUpper(labelValue)
      stencil.Update()

      stat = vtk.vtkImageAccumulate()
      stat.SetInputData(thresh.GetOutput())
      stat.SetStencilData(stencil.GetOutput())
      stat.Update()

      # Add data to statistics list
      cubicMMPerVoxel = reduce(lambda x,y: x*y, segmentLabelmap.GetSpacing())
      ccPerCubicMM = 0.001
      self.statistics[segmentID,"LM voxel count"] = stat.GetVoxelCount()
      self.statistics[segmentID,"LM volume mm3"] = stat.GetVoxelCount() * cubicMMPerVoxel
      self.statistics[segmentID,"LM volume cc"] = stat.GetVoxelCount() * cubicMMPerVoxel * ccPerCubicMM

    def execute_module(self):
        # FIXME: if this module ever becomes anything other than an experiment, build
        # this logic into yet another ProgrammableSource
        
        # make sure marker is up to date
        self._markerSource.GetStructuredPointsOutput().Update()
        self._maskSource.GetStructuredPointsOutput().Update()

        tempJ = vtk.vtkStructuredPoints()
        tempJ.DeepCopy(self._markerSource.GetStructuredPointsOutput())

        self._imageErode.SetInput(tempJ)

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

        self._accum = vtk.vtkImageAccumulate()
        self._accum.SetInput(self._diff.GetOutput())

        # now begin our loop
        stable = False
        while not stable:
            # do erosion, get supremum of erosion and mask I
            self._sup.GetOutput().Update()
            # compare this result with tempJ
            self._diff.SetInput1(tempJ)
            self._diff.SetInput2(self._sup.GetOutput())
            self._accum.Update()

            print "%f == %f ?" % (self._accum.GetMin()[0], self._accum.GetMax()[0])
            if abs(self._accum.GetMin()[0] - self._accum.GetMax()[0]) < 0.0001:
                stable = True
            else:
                # not stable yet...
                print "Trying again..."
                tempJ.DeepCopy(self._sup.GetOutput())