def jaccard3D_pd(pd1,pd2):
    """
    computes the jaccard distance error for two polydata objects

    returns (error) float
    """
    union = vtk.vtkBooleanOperationPolyDataFilter()
    union.SetOperationToUnion()
    union.SetInputData(0,pd1)
    union.SetInputData(1,pd2)
    union.Update()
    u = union.GetOutput()
    massUnion = vtk.vtkMassProperties()
    massUnion.SetInputData(u)

    intersection = vtk.vtkBooleanOperationPolyDataFilter()
    intersection.SetOperationToIntersection()
    intersection.SetInputData(0,pd1)
    intersection.SetInputData(1,pd2)
    intersection.Update()
    i = intersection.GetOutput()
    massIntersection = vtk.vtkMassProperties()
    massIntersection.SetInputData(i)

    return 1 - massIntersection.GetVolume()/massUnion.GetVolume()

def obj(x,sF,devF,R,sv,mv,tv,material,spatial,rc,sc):
    nF = np.dot(R,x[0]*sF)+np.dot(R,x[1]*devF)
    #Make a copy of material configuration and deform this with nF
    nm = vtk.vtkPolyData()
    nm.DeepCopy(material)
    pcoords = vtk.vtkFloatArray()
    pcoords.SetNumberOfComponents(3)
    pcoords.SetNumberOfTuples(nm.GetNumberOfPoints())
    for i in xrange(nm.GetNumberOfPoints()):
        p = [0.,0.,0.]
        nm.GetPoint(i,p)
        p = np.dot(nF,p-rc)
        p.flatten()
        pcoords.SetTuple3(i,p[0]+sc[0],p[1]+sc[1],p[2]+sc[2])

    points = vtk.vtkPoints()
    points.SetData(pcoords)
    nm.SetPoints(points)
    nm.GetPoints().Modified()

    #calculate both the intersection and the union of the two polydata
    intersect = vtk.vtkBooleanOperationPolyDataFilter()
    intersect.SetOperation(1)
    intersect.SetInputData(0,nm)
    intersect.SetInputData(1,spatial)
    intersect.Update()

    union = vtk.vtkBooleanOperationPolyDataFilter()
    union.SetOperation(0)
    union.SetInputData(0,nm)
    union.SetInputData(1,spatial)
    union.Update()

    unionmass = vtk.vtkMassProperties()
    unionmass.SetInputConnection(union.GetOutputPort())
    vol_union = unionmass.GetVolume()

    if intersect.GetOutput().GetNumberOfPoints() > 0:
        intmass = vtk.vtkMassProperties()
        intmass.SetInputConnection(intersect.GetOutputPort())
        vol_int = intmass.GetVolume()
    else:
        #penalize with distance between centroids 
        w = sc.T-np.dot(nF,rc.T)
        c = np.linalg.norm(w)
        vol_int = c*vol_union 

    diff = max([abs(sv-vol_int),abs(sv-vol_union)])
    return diff

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

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

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

    containsClosedSurfaceRepresentation = segmentationNode.GetSegmentation().ContainsRepresentation(
      vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationClosedSurfaceRepresentationName())
    if not containsClosedSurfaceRepresentation:
      return {}

    segment = segmentationNode.GetSegmentation().GetSegment(segmentID)
    closedSurfaceName = vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationClosedSurfaceRepresentationName()
    segmentClosedSurface = segment.GetRepresentation(closedSurfaceName)

    # Compute statistics
    massProperties = vtk.vtkMassProperties()
    massProperties.SetInputData(segmentClosedSurface)

    # Add data to statistics list
    ccPerCubicMM = 0.001
    stats = {}
    if "surface_mm2" in requestedKeys:
      stats["surface_mm2"] = massProperties.GetSurfaceArea()
    if "volume_mm3" in requestedKeys:
      stats["volume_mm3"] = massProperties.GetVolume()
    if "volume_cm3" in requestedKeys:
      stats["volume_cm3"] = massProperties.GetVolume() * ccPerCubicMM
    return stats

        def MakeText(primitive):

            tf.update({primitive: vtk.vtkTriangleFilter()})
            tf[primitive].SetInputConnection(primitive.GetOutputPort())

            mp.update({primitive: vtk.vtkMassProperties()})
            mp[primitive].SetInputConnection(tf[primitive].GetOutputPort())

            # here we capture stdout and write it to a variable for processing.
            summary = StringIO.StringIO()
            # save the original stdout
            old_stdout = sys.stdout
            sys.stdout = summary

            print mp[primitive]
            summary = summary.getvalue()

            startSum = summary.find("  VolumeX")
            endSum = len(summary)
            print summary[startSum:]
            # Restore stdout
            sys.stdout = old_stdout

            vt.update({primitive: vtk.vtkVectorText()})
            vt[primitive].SetText(summary[startSum:])

            pdm.update({primitive: vtk.vtkPolyDataMapper()})
            pdm[primitive].SetInputConnection(vt[primitive].GetOutputPort())

            ta.update({primitive: vtk.vtkActor()})
            ta[primitive].SetMapper(pdm[primitive])
            ta[primitive].SetScale(0.2, 0.2, 0.2)
            return ta[primitive]

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

    def _readstls(self):
        """
        Reads in all STL files contained in directories indicated by **ref_dir** and **def_dir**.
        Also calls **_make3Dmesh()** to create 3-D tetrahedral meshes.

        Returns
        -------
        rsurfs, dsurfs
        """
        for fname in sorted(os.listdir(self._ref_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(
                    str(os.path.normpath(self._ref_dir + os.sep + fname)))
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.rsurfs.append(triangles.GetOutput())
                massProps = vtk.vtkMassProperties()
                massProps.SetInputData(self.rsurfs[-1])
                massProps.Update()
                print(("Generating tetrahedral mesh from {:s}".format(fname)))
                self._make3Dmesh(
                    str(os.path.normpath(self._ref_dir + os.sep + fname)),
                    'MATERIAL', massProps.GetVolume())

        for fname in sorted(os.listdir(self._def_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(
                    str(os.path.normpath(self._def_dir + os.sep + fname)))
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.dsurfs.append(triangles.GetOutput())
                massProps = vtk.vtkMassProperties()
                massProps.SetInputData(self.dsurfs[-1])
                massProps.Update()
                print(("Generating tetrahedral mesh from {:s}".format(fname)))
                self._make3Dmesh(
                    str(os.path.normpath(self._def_dir + os.sep + fname)),
                    'SPATIAL', massProps.GetVolume())

def createMassProperties(
        pdata,
        verbose=1):

    myVTK.myPrint(verbose, "*** createMassProperties ***")

    mass_properties = vtk.vtkMassProperties()
    mass_properties.SetInputData(pdata)

    return mass_properties

def getMassProperties(
        pdata,
        verbose=0):

    mypy.my_print(verbose, "*** getMassProperties ***")

    mass_properties = vtk.vtkMassProperties()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        mass_properties.SetInputData(pdata)
    else:
        mass_properties.SetInput(pdata)

    return mass_properties

def createMassProperties(
        pdata,
        verbose=1):

    myVTK.myPrint(verbose, "*** createMassProperties ***")

    mass_properties = vtk.vtkMassProperties()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        mass_properties.SetInputData(pdata)
    else:
        mass_properties.SetInput(pdata)

    return mass_properties

    def surface_area_vtk(self):
        if self.VTK_installed is False:
            raise VTK_Exception('VTK must be installed to access the surface_area_vtk property')
        if self._surface_area_vtk is None:
            tri_converter = vtk.vtkTriangleFilter()
            tri_converter.SetInputDataObject(self.vtp_mesh)
            tri_converter.Update()
            tri_mesh = tri_converter.GetOutput()
            mass_props = vtk.vtkMassProperties()
            mass_props.SetInputDataObject(tri_mesh)
            self._surface_area_vtk = mass_props.GetSurfaceArea()

            print 'Calculated mesh surface area using VTK Python bindings'

        return self._surface_area_vtk

    def volume_vtk(self):
        if self.VTK_installed is False:
            raise VTK_Exception('VTK must be installed to access the volume_vtk property')

        if self._volume_vtk is None:
            tri_converter = vtk.vtkTriangleFilter()
            tri_converter.SetInputDataObject(self.vtp_mesh)
            tri_converter.Update()
            tri_mesh = tri_converter.GetOutput()
            mass_props = vtk.vtkMassProperties()
            mass_props.SetInputDataObject(tri_mesh)
            self._volume_vtk = mass_props.GetVolume()

            print 'Calculated mesh volume using VTK library'

        return self._volume_vtk

def calculateVolume(polyData):
    """
    Calculates the volume of a polydata cell.

    Input:
        polyData - vtkPolyData object - Contains a cell defined by triangular polygons and
            needs to be waterthight.

    Returns:
        float - The volume of the cell in the respecitve units.
    """
    # Mass prop
    mp = vtk.vtkMassProperties()
    # Set the input
    mp.SetInputData(polyData)
    return float(mp.GetVolume())

def _test():
    # import the stl surfaces
    bonePolydata = importStlSurf('oks001_MRC_FMB_LVTIT_10.stl')
    # cartPolydata = importStlSurf('oks001_MRC_FMC_LVTIT_10.stl')
    #
    # boolFilter = vtk.vtkBooleanOperationPolyDataFilter()
    # boolFilter.SetOperationToIntersection() # Define the type of boolean filter
    # boolFilter.SetInputData(0, bonePolydata)
    # boolFilter.SetInputData(1, cartPolydata)
    # # boolFilter.Update()
    # overlapPolydata = boolFilter.GetOutput()

    volFilter = vtk.vtkMassProperties()
    volFilter.SetInputData(bonePolydata)
    volFilter.Update()
    print volFilter.GetVolume()

    return

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

  appender = vtk.vtkAppendPolyData()
  mp = vtk.vtkMassProperties()
  for idx in range(len(files)):
    file = files[idx]
    polydata = CreateCrossSectionPolyDataFromFile(file)
    mp.SetInputData(polydata)
    mp.Update()
    areaFile.write(str(mp.GetSurfaceArea()))
    areaFile.write('\n')
    appender.AddInputData(polydata)
    
  writer = vtk.vtkXMLPolyDataWriter()
  writer.SetFileName(outputGeometryFile)
  writer.SetInputConnection(appender.GetOutputPort())
  writer.Update()

  areaFile.close()

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: No input surface.')

        cleaner = vtk.vtkCleanPolyData()
        cleaner.SetInput(self.Surface)
        cleaner.Update()

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInput(cleaner.GetOutput())
        triangleFilter.Update()

        massProps = vtk.vtkMassProperties()
        massProps.SetInput(triangleFilter.GetOutput())
        massProps.Update()

        self.SurfaceArea = massProps.GetSurfaceArea()
        self.Volume = massProps.GetVolume()
        self.ShapeIndex = massProps.GetNormalizedShapeIndex()

  def addSegmentClosedSurfaceStatistics(self):
    import vtkSegmentationCorePython as vtkSegmentationCore

    containsClosedSurfaceRepresentation = self.segmentationNode.GetSegmentation().ContainsRepresentation(
      vtkSegmentationCore.vtkSegmentationConverter.GetSegmentationClosedSurfaceRepresentationName())
    if not containsClosedSurfaceRepresentation:
      return

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

      # Compute statistics
      massProperties = vtk.vtkMassProperties()
      massProperties.SetInputData(segmentClosedSurface)

      # Add data to statistics list
      ccPerCubicMM = 0.001
      self.statistics[segmentID,"CS surface mm2"] = massProperties.GetSurfaceArea()
      self.statistics[segmentID,"CS volume mm3"] = massProperties.GetVolume()
      self.statistics[segmentID,"CS volume cc"] = massProperties.GetVolume() * ccPerCubicMM

    def UpdateMathValues(self):

        connection = self._app_states[self._current_image_index].GetFactory().GetInputConnection()
        object = connection.GetProducer().GetOutput()

        # determine number of connected components
        _f = vtk.vtkPolyDataConnectivityFilter()
        _f.SetInputConnection(connection)
        _f.Update()
        nregions = _f.GetNumberOfExtractedRegions()

        self._app_states[self._current_image_index]._poly_count = object.GetNumberOfPolys()
        self._app_states[self._current_image_index]._volume, self._app_states[self._current_image_index]._area = (
            0.0,
            0.0,
        )
        self._app_states[self._current_image_index]._region_count = nregions

        if self._app_states[self._current_image_index]._poly_count > 0:
            math0 = vtk.vtkMassProperties()
            math0.SetInputConnection(connection)
            math0.Update()
            self._app_states[self._current_image_index]._volume, self._app_states[self._current_image_index]._area = (
                math0.GetVolume(),
                math0.GetSurfaceArea(),
            )

        self.UpdateStatisticsText()

        # log some results
        logger = logging.getLogger("results")
        logger.info(
            "Isosurface plugin (%s):" % (component.getUtility(ICurrentViewportManager).GetPageState().GetTitle())
        )
        logger.info(u"\tArea: %0.3f mm\u00B2" % self._app_states[self._current_image_index]._area)
        logger.info(u"\tVolume: %0.3f mm\u00B3" % (self._app_states[self._current_image_index]._volume))
        logger.info("\tPolygon Count: %d" % (self._app_states[self._current_image_index]._poly_count))
        logger.info("\tRegion Count: %d" % (self._app_states[self._current_image_index]._region_count))
        logger.info("")

def volume(surface):
    """Compute volume of a closed triangulated surface mesh."""
    properties = vtk.vtkMassProperties()
    properties.SetInput(surface)
    properties.Update()
    return properties.GetVolume()

#.........这里部分代码省略.........
			labelMap = labelStatistics.GetOutput()
			numberOfLabels = labelMap.GetNumberOfLabelObjects()
		else:
			labelShape = itk.LabelShapeImageFilter[labelITK].New()
			labelShape.SetInput(labelITK)
			data = labelShape.GetOutput()
			data.Update()
			numberOfLabels = labelShape.GetNumberOfLabels()
		
			if self.parameters["AvgInt"]:
				avgintCalc = itk.LabelStatisticsImageFilter[origITK,labelITK].New()
				avgintCalc.SetInput(origITK)
				avgintCalc.SetLabelInput(labelITK)
				avgintCalc.Update()

		self.progressObj.setProgress(0.2)
		self.updateProgress(None, "ProgressEvent")

		# Area calculation pipeline
		if self.parameters["Area"]:
			voxelArea = x*y*2 + x*z*2 + y*z*2
			largestSize = labelITK.GetLargestPossibleRegion().GetSize()
			# if 2D image, calculate area using volume
			if largestSize.GetSizeDimension() > 2 and largestSize.GetElement(2) > 1:
				areaSpacing = labelVTK.GetSpacing()
				objectThreshold = vtk.vtkImageThreshold()
				objectThreshold.SetInput(labelVTK)
				objectThreshold.SetOutputScalarTypeToUnsignedChar()
				objectThreshold.SetInValue(255)
				objectThreshold.SetOutValue(0)
				marchingCubes = vtk.vtkMarchingCubes()
				#marchingCubes.SetInput(labelVTK)
				marchingCubes.SetInput(objectThreshold.GetOutput())
				massProperties = vtk.vtkMassProperties()
				massProperties.SetInput(marchingCubes.GetOutput())
				areaDiv = (areaSpacing[0] / x)**2
				
				if self.parameters["Smoothness"]:
					smoothDecimate = vtk.vtkDecimatePro()
					smoothProperties = vtk.vtkMassProperties()
					smoothDecimate.SetTargetReduction(0.9)
					smoothDecimate.PreserveTopologyOff()
					smoothDecimate.SetInput(marchingCubes.GetOutput())
					smoothProperties.SetInput(smoothDecimate.GetOutput())

		# Filter needed for axes calculations
		if self.parameters["Axes"] and newITKStatistics:
			labelGeometry = itk.LabelGeometryImageFilter[labelITK,labelITK].New()
			labelGeometry.SetCalculateOrientedBoundingBox(1)
			labelGeometry.SetInput(labelITK)
			labelGeometry.Update()

		# Get results and do some calculations for each object
		tott = 0
		voxelSize = voxelSizes[0] * voxelSizes[1] * voxelSizes[2]
		for i in range(startIntensity, numberOfLabels+1):
			areaInUm = 0.0
			avgInt = 0
			avgIntStdErr = 0.0
			roundness = 0.0
			objIntSum = 0.0
			minorLength = 0.0
			majorLength = 0.0
			elongation = 0.0
			angleMinX = 0.0
			angleMinY = 0.0

#.........这里部分代码省略.........
        normals.SetFeatureAngle(80)
        normals.AutoOrientNormalsOn()
        normals.GetOutput().ReleaseDataFlagOn()
        normals.Update()
        del polydata
        polydata = normals.GetOutput()
        #polydata.Register(None)
        polydata.SetSource(None)
        del normals

        # Improve performance
        stripper = vtk.vtkStripper()
        stripper.ReleaseDataFlagOn()
        stripper_ref = weakref.ref(stripper)
        stripper_ref().AddObserver("ProgressEvent", lambda obj,evt:
                        UpdateProgress(stripper_ref(), _("Creating 3D surface...")))
        stripper.SetInput(polydata)
        stripper.PassThroughCellIdsOn()
        stripper.PassThroughPointIdsOn()
        stripper.GetOutput().ReleaseDataFlagOn()
        stripper.Update()
        del polydata
        polydata = stripper.GetOutput()
        #polydata.Register(None)
        polydata.SetSource(None)
        del stripper

        # Map polygonal data (vtkPolyData) to graphics primitives.
        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInput(polydata)
        mapper.ScalarVisibilityOff()
        mapper.ReleaseDataFlagOn()
        mapper.ImmediateModeRenderingOn() # improve performance

        # Represent an object (geometry & properties) in the rendered scene
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        del mapper
        #Create Surface instance
        if overwrite:
            surface = Surface(index = self.last_surface_index)
        else:
            surface = Surface(name=surface_name)
        surface.colour = colour
        surface.polydata = polydata
        del polydata

        # Set actor colour and transparency
        actor.GetProperty().SetColor(colour)
        actor.GetProperty().SetOpacity(1-surface.transparency)

        prop = actor.GetProperty()

        interpolation = int(ses.Session().surface_interpolation)

        prop.SetInterpolation(interpolation)

        proj = prj.Project()
        if overwrite:
            proj.ChangeSurface(surface)
        else:
            index = proj.AddSurface(surface)
            surface.index = index
            self.last_surface_index = index

        session = ses.Session()
        session.ChangeProject()

        # The following lines have to be here, otherwise all volumes disappear
        measured_polydata = vtk.vtkMassProperties()
        measured_polydata.ReleaseDataFlagOn()
        measured_polydata.SetInput(to_measure)
        volume =  float(measured_polydata.GetVolume())
        surface.volume = volume
        self.last_surface_index = surface.index
        del measured_polydata
        del to_measure

        Publisher.sendMessage('Load surface actor into viewer', actor)

        # Send actor by pubsub to viewer's render
        if overwrite and self.actors_dict.keys():
            old_actor = self.actors_dict[self.last_surface_index]
            Publisher.sendMessage('Remove surface actor from viewer', old_actor)

        # Save actor for future management tasks
        self.actors_dict[surface.index] = actor

        Publisher.sendMessage('Update surface info in GUI',
                                    (surface.index, surface.name,
                                    surface.colour, surface.volume,
                                    surface.transparency))
        
        #When you finalize the progress. The bar is cleaned.
        UpdateProgress = vu.ShowProgress(1)
        UpdateProgress(0, _("Ready"))
        Publisher.sendMessage('Update status text in GUI', _("Ready"))
        
        Publisher.sendMessage('End busy cursor')
        del actor