Mesh a label image and write cell info to file (pyhton)

This example takes a label image in Meta format and meshes each label. It then prints area and normal of each mesh cell to a file. An example data can be found here Media:labels.tgz

<source lang="python"> import vtk

input='labels.mhd'

1. Prepare to read the file

readerVolume = vtk.vtkMetaImageReader() readerVolume.SetFileName(input) readerVolume.Update()

1. Extract the region of interest

voi = vtk.vtkExtractVOI() voi.SetInput(readerVolume.GetOutput())

1. voi.SetVOI(0,517, 0,228, 0,392)

voi.SetSampleRate(1,1,1)

1. voi.SetSampleRate(3,3,3)

voi.Update()#necessary for GetScalarRange() srange= voi.GetOutput().GetScalarRange()#needs Update() before! print "Range", srange

1. 1. Prepare surface generation
2. contour = vtk.vtkContourFilter()
3. contour = vtk.vtkMarchingCubes()

contour = vtk.vtkDiscreteMarchingCubes() #for label images contour.SetInput( voi.GetOutput() ) contour.ComputeNormalsOn()

1. 1. run through all labels

for index in range(1, int(srange[1]) + 1):

   print "Doing label", index

   contour.SetValue(0, index)
   contour.Update() #needed for GetNumberOfPolys() !!!

   smoother= vtk.vtkWindowedSincPolyDataFilter()
   smoother.SetInput(contour.GetOutput());
   smoother.SetNumberOfIterations(5);
   #smoother.BoundarySmoothingOff();
   #smoother.FeatureEdgeSmoothingOff();
   #smoother.SetFeatureAngle(120.0);
   #smoother.SetPassBand(.001);
   smoother.NonManifoldSmoothingOn();
   smoother.NormalizeCoordinatesOn();
   smoother.Update();

   ##calc cell normal
   triangleCellNormals= vtk.vtkPolyDataNormals()
   triangleCellNormals.SetInput(smoother.GetOutput())
   triangleCellNormals.ComputeCellNormalsOn()
   triangleCellNormals.ComputePointNormalsOff()
   triangleCellNormals.ConsistencyOn()
   triangleCellNormals.AutoOrientNormalsOn()
   triangleCellNormals.Update() #creates vtkPolyData

   ##calc cell area
   triangleCellAN= vtk.vtkMeshQuality()
   triangleCellAN.SetInput(triangleCellNormals.GetOutput())
   triangleCellAN.SetTriangleQualityMeasureToArea()
   triangleCellAN.SaveCellQualityOn() #default
   triangleCellAN.Update() #creates vtkDataSet

   PointNormalArray = triangleCellNormals.GetOutput().GetCellData().GetNormals()
   qualityArray = triangleCellAN.GetOutput().GetCellData().GetArray("Quality")

   if (PointNormalArray.GetNumberOfTuples() != qualityArray.GetNumberOfTuples()):
       print "Error! Sizes of normal array and area array dont equal!"
       exit(1)

   f= open('label_stat' + "_%.4d" % index + ".dat", 'w')
   print >> f, "#cell_index\tarea\tn_x\tn_y\tn_z"

   for i in range(0, PointNormalArray.GetNumberOfTuples()):
       
       pointNormal= PointNormalArray.GetTuple3(i) #this is for 3D data in python
       area= qualityArray.GetValue(i)
       print >> f, i, area, pointNormal[0], pointNormal[1], pointNormal[2]
       
   f.close()