Decimate Tets with VTK

decimatetets.py

import vtk
import os,sys,math

# read unstructured grid
r=vtk.vtkXMLUnstructuredGridReader()
r.SetFileName(sys.argv[1]) # first command line argument after script name
r.Update()

orig=r.GetOutput()
pts=orig.GetPoints()

# calculate center of mass
com=vtk.vtkCenterOfMass()
com.SetInputData(orig)
com.SetUseScalarsAsWeights(False)
com.Update()
center=com.GetCenter()

# A field is needed to determine "nearness" apparently. 
# There's mention of this only in their test case: https://www.vtk.org/gitweb?p=VTK.git;a=blob;f=Filters/Core/Testing/Cxx/TestUnstructuredGridQuadricDecimation.cxx
radius=vtk.vtkDoubleArray()
radius.SetName('radius')
radius.SetNumberOfComponents(1)
radius.SetNumberOfTuples(pts.GetNumberOfPoints())

# calculate radius as each point's distance from the CoM
for i in range(pts.GetNumberOfPoints()):
    pt=pts.GetPoint(i)
    dist=math.sqrt(sum((p-c)**2 for p,c in zip(pt,center)))
    radius.SetTypedTuple(i,(dist,))

orig.GetPointData().SetScalars(radius)

# finally we can decimate the mesh
deci=vtk.vtkUnstructuredGridQuadricDecimation()
deci.SetInputData(orig)
#deci.SetTargetReduction(0.01) # I don't understand what this value means so use the following
deci.SetNumberOfTetsOutput(int(orig.GetNumberOfCells()*0.01)) # reduce to 1%
deci.SetScalarsName('radius')
deci.Update()
decimated=deci.GetOutput()

# save out with modified filename
w=vtk.vtkXMLUnstructuredGridWriter()
w.SetCompressorTypeToNone()
w.SetDataModeToAscii()
w.SetFileName(os.path.splitext(sys.argv[1])[0]+'_decimated.vtu')
w.SetInputData(decimated)

w.Write()